平差期末作业

测绘专用 MATLAB

7-6-51

Xa = 4728.008;Ya = 227.880;
Xb = 4604.993;Yb = 362.996;
Xc = 4750.191;Yc = 503.152;
Xp = 4881.27;Yp = 346.86;
%输入数据
%将角度（度、分、秒）转换为弧度
L1 = [94 29 32];L2 = [44 20 36.3];
L3 = [47 19 43.3];L4 = [85 59 51.6];
L10 = dms2degrees(L1);
L20 = dms2degrees(L2);
L30 = dms2degrees(L3);
L40= dms2degrees(L4);
L10=deg2rad(L10);
L20=deg2rad(L20);
L30=deg2rad(L30);
L40=deg2rad(L40);
p=206265;
%计算边长
S1=sqrt((Xa-Xp)^2+(Ya-Yp)^2);
S2=sqrt((Xb-Xp)^2+(Yb-Yp)^2);
S3=sqrt((Xc-Xp)^2+(Yc-Yp)^2);
%计算角度观测近似值
L10 = rad2deg(atan2(Yb-Ya,Xb-Xa)) - rad2deg(atan2(Yp-Ya,Xp-Xa)) ;
L10 = degrees2dms(L10);
L20 = (360 - rad2deg(atan2(Ya-Yb,Xa-Xb))) + rad2deg(atan2(Yp-Yb,Xp-Xb)) ;
L20 = L20-360;
L20 = degrees2dms(L20);
L30 = rad2deg(atan2(Yc-Yb,Xc-Xb)) - rad2deg(atan2(Yp-Yb,Xp-Xb)) ;
L30 = degrees2dms(L30);
L40 = rad2deg(atan2(Yp-Yc,Xp-Xc)) - rad2deg(atan2(Yb-Yc,Xb-Xc)) ;
L40 = degrees2dms(L40);
%计算误差方程系数B
B = [
    -206265*(Ya-Yp)/(S1^2*100) 206265*(Xa-Xp)/(S1^2*100);
    -206265*(Yb-Yp)/(S2^2*100) 206265*(Xb-Xp)/(S2^2*100);
    -206265*(Yb-Yp)/(S2^2*100) 206265*(Xb-Xp)/(S2^2*100);
    -206265*(Yc-Yp)/(S3^2*100) 206265*(Xc-Xp)/(S3^2*100)
        ];
    %计算误差方程系数l
    l1 = L1-L10;
    l2 = L2-L20;
    l3 = L3-L30;
    l4 = L4-L40;
    l=[
        l1(3);
        l2(3);
        l3(3);
        l4(3)
        ];
    P = [
        1 0 0 0;
        0 1 0 0;
        0 0 1 0;
        0 0 0 1
        ];
   %计算Nbb、W
   Nbb = (B')*P*B;
   W = (B')*P*l;
   %解法方程
   x=(Nbb^-1)*W;
      Xp0 = Xp+x(1)/10; 
      Yp0 = Yp+x(2)/10;
      sprintf('%.3f',Xp0)
      sprintf('%.3f',Yp0)
      %计算协因数
      Q = vpa(Nbb^-1,4);

7-7-58

Xa = 586.843; Ya = 488.027;
Xb = 776.407;Yb = 568.693;
Xc = 795.565; Yc = 191.581;
Xp1 = 880.267; Yp1 = 367.025;
Xp2 = 585.832; Yp2 = 238.972;
S1 = 249.115;
S2 = 380.913;
S3 = 317.406;
S4 = 226.930;
S5 = 321.154;
S6 = 215.109;
S7 = 194.845;
%输入数据
%计算边长近似值
S10 = sqrt((Xp2-Xa)^2+(Yp2-Ya)^2);
S20 = sqrt((Xp2-Xb)^2+(Yp2-Yb)^2);
S30 = sqrt((Xp1-Xa)^2+(Yp1-Ya)^2);
S40 = sqrt((Xp1-Xb)^2+(Yp1-Yb)^2);
S50 = sqrt((Xp2-Xp1)^2+(Yp2-Yp1)^2);
S60 = sqrt((Xp2-Xc)^2+(Yp2-Yc)^2);
S70 = sqrt((Xp1-Xc)^2+(Yp1-Yc)^2);
%计算B矩阵
B = [0 0 -(Xa-Xp2)/S10 -(Ya-Yp2)/S10;
        0 0 -(Xb-Xp2)/S10 -(Yb-Yp2)/S20;
        -(Xa-Xp1)/S30 -(Ya-Yp1)/S30 0 0;
        -(Xb-Xp1)/S40 -(Yb-Yp1)/S40 0 0;
        (Xp1-Xp2)/S50 (Yp1-Yp2)/S50 -(Xp1-Xp2)/S50 -(Yp1-Yp2)/S50;
        0 0 (Xp2-Xc)/S60 (Yp2-Yc)/S60;
        (Xp1-Xc)/S70 (Yp1-Yc)/S70 0 0];
%计算常数项
l = [S1-S10;
        S2-S20;
        S3-S30;
        S4-S40;
        S5-S50;
        S6-S60;
        S7-S70];
    %计算权阵
  P = [100/S1 0 0 0 0 0 0;
        0 100/S2 0 0 0 0 0;
        0 0 100/S3 0 0 0 0;
        0 0 0 100/S4 0 0 0;
        0 0 0 0 100/S5 0 0;
        0 0 0 0 0 100/S6 0;
        0 0 0 0 0 0 100/S7];
    %计算法方程系数Nbb、W
    Nbb = (B')*P*B;
    W = (B')*P*l;
    %解法方程
    x = (Nbb^-1)*W;
    %求观测值改正数
    V = B*x-l;
    %求观测值平差值
    v = [
        S1+V(1);
        S2+V(2);
        S3+V(3);
        S4+V(4);
        S5+V(5);
        S6+V(6);
        S7+V(7);
        ];
    %求坐标平差值
    Xp10 = Xp1+x(1);
    Yp10 = Yp1+x(2);
    Xp20 = Xp2+x(3);
    Yp20 = Yp2+x(4);
%单位权中误差
e0 = sqrt((V')*P*V/(7-4));
 P0 = Nbb^-1;
 exp1 = e0*sqrt(P0(1,1));
 eyp1 = e0*sqrt(P0(2,2));
 ep1 = sqrt(exp1^2+eyp1^2);
 exp2 = e0*sqrt(P0(3,3));
 eyp2 = e0*sqrt(P0(4,4));
 ep2 = sqrt(exp2^2+eyp2^2);
 PP = [exp1 eyp1 ep1;
            exp2 eyp2 ep2];   

7-8-61