MATLAB——测量平差数据处理

野性酷女 2022-09-15 13:56 57阅读 0赞

\*\*

## 分别使用间接平差、秩亏自由网平差、拟稳平差求解最基本的水准网平差，仅做理解展示，代码较为简单。 ##

\*\*

## 第一题 ##

题目如下：  
![在这里插入图片描述][0cfde11838e744dc882c20090f4eca9a.png]  
![在这里插入图片描述][9d463def9da54d2cbc2cd60ebc3bae43.png]

![在这里插入图片描述][b5120b3aa7dd4b97bf3af1fb97f663dd.png]  
![在这里插入图片描述][efc926f1e9954ae69b1746ea094d44df.png]  
![在这里插入图片描述][bc8302f03d3147268b86e3c934a77ca2.png]  
（1）间接平差

format short
    clc 
    clear all
    %已知数据：
    h1=12.345;
    h2=3.478;
    h3=-15.817;
    %近似值：各点近似高程为：
    H10=0;
    H20=12.345;
    H30=15.823;
    
    %间接平差
    n=3;t=1;r=2;
    %设x1点高程已知，为0m，x2,x3的平差值为xc2,xc3
    x1=0;
    %计算x2,x3的近似值x02,x03
    x02=x1+h1;
    x03=x1-h3;
    %列观测方程
    %hc1=Xc2-x1;
    %hc2=-Xc2+Xc3;
    %hc3=-Xc3+x1;
    
    %设s路线长的观测误差中误差为单位权重中误差
    %则权阵P
    P=diag([1 1 1]);
    
    B=[1 0;-1 1;0 -1]; %系数矩阵
    l=[0 6 0];%常数阵，单位为mm
    l=l';
    N=(B')*P*B;%法方程系数阵
    W=B'*P*l;%法方程常数阵
    
    xc=inv(N)*W %参数计算
    
    xc1=xc.*0.001;
    Qxx=inv(N);%协因数阵
    Qxx=sym(Qxx)%分数形式显示
    
    M=l'*P*l-W'*xc;
    sita=sqrt(M/r);%单位权中误差
    
    xc2=xc1(1);
    xc3=xc1(2);
    
    %列误差方程
    v1=xc2;
    v2=-xc2+xc3-0.006;
    v3=-xc3;
    
    V=[v1 v2 v3];%改正数矩阵
    V=V'
    
    
    X0=[x02 x03];X0=X0';
    Xc=X0+xc1 %参数平差值
    
    h=[h1 h2 h3];
    h=h';
    H=h+V %高程平差值

（2）秩亏自由网平差

format short
    clc 
    clear all
    %已知数据：
    h1=12.345;
    h2=3.478;
    h3=-15.817;
    %近似值：各点近似高程为：
    H01=0;
    H02=12.345;
    H03=15.823;
    
    %秩亏自由网平差
    n=3;t=1;r=2;
    
    %列观测方程
    %hc1=-Xc1+Xc2;;
    %hc2=-Xc2+Xc3;
    %hc3=Xc1-Xc3；
    
    
    %设s路线长的观测误差中误差为单位权重中误差
    %则权阵P
    P=diag([1 1 1]);
    
    
    B=[-1 1 0;0 -1 1;1 0 -1]; %系数矩阵
    l=[0 0 6];%常数阵，单位为mm
    l=l';
    d=t-rank(B);
    N=(B')*P*B;%法方程系数阵
    W=B'*P*l;%法方程常数阵
    Nm=N*pinv(N*N);%求伪逆
    xc=Nm*W %参数计算
    
    
    xcl=xc.*0.001;
    
    
    Qxx=Nm*Nm*N;%协因数阵
    Qxx=sym(Qxx)%分数形式显示
    M=l'*P*l-W'*xc;
    sita=sqrt(M/r);%单位权中误差
    
    
    xc1=xcl(1);
    xc2=xcl(2);
    xc3=xcl(3);
    
    %列误差方程
    v1=-xc1+xc2;
    v2=-xc2+xc3;
    v3=xc1-xc3-0.006;
    
    
    V=[v1 v2 v3];%改正数矩阵
    V=V';
    
    
    X0=[H01 H02 H03];X0=X0';
    Xc=X0+xcl %高程参数平差值
    
    h=[h1 h2 h3];
    h=h';
    H=h+V %高差平差值

## 第二题 ##

题目如下：  
![在这里插入图片描述][80ae7a8340a642bc85975a781e752b0c.png]![在这里插入图片描述][watermark_type_ZHJvaWRzYW5zZmFsbGJhY2s_shadow_50_text_Q1NETiBAbm9jaGVuZ3pp_size_11_color_FFFFFF_t_70_g_se_x_16]![在这里插入图片描述][watermark_type_ZHJvaWRzYW5zZmFsbGJhY2s_shadow_50_text_Q1NETiBAbm9jaGVuZ3pp_size_20_color_FFFFFF_t_70_g_se_x_16]![在这里插入图片描述][351e6f1c4ba04391813daefb636f7978.png]  
![在这里插入图片描述][watermark_type_ZHJvaWRzYW5zZmFsbGJhY2s_shadow_50_text_Q1NETiBAbm9jaGVuZ3pp_size_20_color_FFFFFF_t_70_g_se_x_16 1]  
拟稳平差

format short
    clc 
    clear all
    %已知据 ：
    h=[86.809 25.714 31.225 -71.952 -61.084 -44.178 10.847 16.350 11.409];
    s=[204.1 188.7 344.8 149.2  142.9 250.0 128.2 98.0 196.1];
    P=diag([0.049 0.053 0.029 0.067 0.070 0.040 0.078 0.102 0.051]);
    %近似值：各点近似高程为
    X0=[-0.004 86.806 14.846 25.706 31.216 42.626];
    sita=1.85;
    
    %经典自由网平差
    %{
        
    n=9;t=5;r=4;
    %以6号点为固定点
    x6=0;
    %设x1,x2，x3，x4，x5的平差值为xc1，xc2，xc3，xc4，xc5
    %计算x1,x2，x3，x4，x5的的近似值x01,x02,x03,x04,x05
    x01=x6-h(6)-h(1);
    x02=x6-h(6);
    x03=x02+h(4);
    x04=x02+h(5);
    x05=x6-h(9);
    
    h01=-x01+x02;
    h02=-x01+x04;
    h03=-x01+x05;
    h04=-x02+x03;
    h05=-x02+x04;
    h06=-x02+x6;
    h07=-x03+x04;
    h08=-x03+x05;
    h09=-x05+x6;
    
    h0=[h01 h02 h03 h04 h05 h06 h07 h08 h09];
    
    %{
        
    列观测方程
    hc1=-xc1+xc2;
    hc2=-xc1+xc4;
    hc3=-xc1+xc5;
    hc4=-xc2+xc3;
    hc5=-xc2+xc4;
    hc6=-xc2+x6;
    hc7=-xc3+xc4;
    hc8=-xc3+xc5;
    hc9=-xc5+x6;
    %}
    
    
    
    %计算常数阵
    l=[0 0 0 0 0 0 0 0 0];
    for i=1:9
        l(i)=h(i)-h0(i);
    end
    l=l';
    l=l.*1000; %常数项扩大1000倍
    
    B=[-1 1 0 0 0  
       -1 0 0 1 0 
       -1 0 0 0 1 
        0 -1 1 0 0 
        0 -1 0 1 0 
        0 -1 0 0 0
        0 0 -1 1 0 
        0 0 -1 0 1 
        0 0 0 0 -1 ];
    
    N=(B')*P*B;%法方程系数阵
    W=B'*P*l;%法方程常数阵
    
    
    xc=inv(N)*W %参数计算
    xcl=xc.*0.001
    Qxx=inv(N);%协因数阵
    Qxx=sym(Qxx)%分数形式显示
    M=l'*P*l-W'*xc;
    sita=sqrt(M/r);%单位权中误差
    
    xc1=xcl(1);
    xc2=xcl(2);
    xc3=xcl(3);
    xc4=xcl(4);
    xc5=xcl(5);
    
    l2=l./1000;%常数项缩小1000倍
    %列误差方程
    v1=-xc1+xc2-l2(1);
    v2=-xc1+xc4-l2(2);
    v3=-xc1+xc5-l2(3);
    v4=-xc2+xc3-l2(4);
    v5=-xc2+xc4-l2(5);
    v6=-xc2+x6-l2(6);
    v7=-xc3+xc4-l2(7);
    v8=-xc3+xc5-l2(8);
    v9=-xc5+x6-l2(9);
    
    V=[v1 v2 v3 v4 v5 v6 v7 v8 v9 ];%改正数矩阵
    V=V'
    
    
    X00=[ x01 x02 x03 x04  x05];X00=X00';
    Xc=X00+xcl %参数平差值
    
    h=h';
    H=h+V %高程平差值
    
    %}
    
    
    %秩亏自由网平差
    %{
        
    n=9;t=5;r=4;
    %水准网中基准条件d=1
    S=[1 1 1 1 1 1];
    S=S';
    Px=diag([1 1 1 1 1 1 ]);
    %设x1,x2，x3，x4，x5,x6的平差值为xc1，xc2，xc3，xc4，xc5,xc6
    %x1,x2，x3，x4，x5,x6的近似值x01,x02,x03,x04,x05,x06
    x0=[-0.004 86.806 14.846 25.706 31.216 42.626];
    
    %{
        
    %列观测方程
    hc1=-xc1+xc2;
    hc2=-xc1+xc4;
    hc3=-xc1+xc5;
    hc4=-xc2+xc3;
    hc5=-xc2+xc4;
    hc6=-xc2+x6;
    hc7=-xc3+xc4;
    hc8=-xc3+xc5;
    hc9=-xc5+x6;
    %}
    
    h01=-x0(1)+x0(2);
    h02=-x0(1)+x0(4);
    h03=-x0(1)+x0(5);
    h04=-x0(2)+x0(3);
    h05=-x0(2)+x0(4);
    h06=-x0(2)+x0(6);
    h07=-x0(3)+x0(4);
    h08=-x0(3)+x0(5);
    h09=-x0(5)+x0(6);
    
    h0=[h01 h02 h03 h04 h05 h06 h07 h08 h09];
    
    %计算常数阵
    l=[0 0 0 0 0 0 0 0 0];
    for i=1:9
        l(i)=h(i)-h0(i);
    end
    l=l';
    l=l.*1000; %常数项扩大1000倍
    
    B=[-1 1 0 0 0 0
       -1 0 0 1 0 0
       -1 0 0 0 1 0
        0 -1 1 0 0 0
        0 -1 0 1 0 0
        0 -1 0 0 0 1
        0 0 -1 1 0 0
        0 0 -1 0 1 0
        0 0 0 0 -1 1 ];
    
    N=(B')*P*B;%法方程系数阵
    W=B'*P*l;%法方程常数阵
    
    Z=Px*S*S'*Px;
    xc=inv(N+Z)*W;
    xcl=xc.*0.001
    Qxx=inv(N+Z)*N*inv(N+Z);%协因数阵
    Qxx=sym(Qxx)%分数形式显示
    M=l'*P*l-W'*xc;
    sita=sqrt(M/r)%单位权中误差
    
    
    
    xc1=xcl(1);
    xc2=xcl(2);
    xc3=xcl(3);
    xc4=xcl(4);
    xc5=xcl(5);
    xc6=xcl(6);
    
    l2=l./1000;%常数项缩小1000倍
    %列误差方程
    v1=-xc1+xc2-l2(1);
    v2=-xc1+xc4-l2(2);
    v3=-xc1+xc5-l2(3);
    v4=-xc2+xc3-l2(4);
    v5=-xc2+xc4-l2(5);
    v6=-xc2+xc6-l2(6);
    v7=-xc3+xc4-l2(7);
    v8=-xc3+xc5-l2(8);
    v9=-xc5+xc6-l2(9);
    
    V=[v1 v2 v3 v4 v5 v6 v7 v8 v9 ];%改正数矩阵
    V=V'
    
    
    x0=x0';
    Xc=x0+xcl %参数平差值
    
    
    h=h';
    H=h+V %高程平差值
    %}
    
    
    %拟稳平差
    n=9;t=5;r=4;
    %水准网中基准条件d=1
    S=[1 1 1 1 1 1];
    S=S';
    Px=diag([1 1 0 0 1 1 ]);
    %设x1,x2，x3，x4，x5,x6的平差值为xc1，xc2，xc3，xc4，xc5,xc6
    %x1,x2，x3，x4，x5,x6的近似值x01,x02,x03,x04,x05,x06
    x0=[-0.004 86.806 14.846 25.706 31.216 42.626];
    x0=x0';
    
    %{
        
    %列观测方程
    hc1=-xc1+xc2;
    hc2=-xc1+xc4;
    hc3=-xc1+xc5;
    hc4=-xc2+xc3;
    hc5=-xc2+xc4;
    hc6=-xc2+x6;
    hc7=-xc3+xc4;
    hc8=-xc3+xc5;
    hc9=-xc5+x6;
    %}
    
    h01=-x0(1)+x0(2);
    h02=-x0(1)+x0(4);
    h03=-x0(1)+x0(5);
    h04=-x0(2)+x0(3);
    h05=-x0(2)+x0(4);
    h06=-x0(2)+x0(6);
    h07=-x0(3)+x0(4);
    h08=-x0(3)+x0(5);
    h09=-x0(5)+x0(6);
    
    h0=[h01 h02 h03 h04 h05 h06 h07 h08 h09];
    
    %计算常数阵
    l=[0 0 0 0 0 0 0 0 0];
    for i=1:9
        l(i)=h(i)-h0(i);
    end
    l=l';
    l=l.*1000; %常数项扩大1000倍
    
    B=[-1 1 0 0 0 0
       -1 0 0 1 0 0
       -1 0 0 0 1 0
        0 -1 1 0 0 0
        0 -1 0 1 0 0
        0 -1 0 0 0 1
        0 0 -1 1 0 0
        0 0 -1 0 1 0
        0 0 0 0 -1 1 ];
    
    N=(B')*P*B;%法方程系数阵
    W=B'*P*l;%法方程常数阵
    
    Z=Px*S*S'*Px;
    xc=inv(N+Z)*W;
    xcl=xc.*0.001
    Qxx=inv(N+Z)*N*inv(N+Z); %协因数阵
    Qxx=sym(Qxx)%分数形式显示
    M=l'*P*l-W'*xc;
    sita=sqrt(M/r);%单位权中误差
    
    xc1=xcl(1);
    xc2=xcl(2);
    xc3=xcl(3);
    xc4=xcl(4);
    xc5=xcl(5);
    xc6=xcl(6);
    
    l2=l./1000;%常数项缩小1000倍
    
    %列误差方程
    v1=-xc1+xc2-l2(1);
    v2=-xc1+xc4-l2(2);
    v3=-xc1+xc5-l2(3);
    v4=-xc2+xc3-l2(4);
    v5=-xc2+xc4-l2(5);
    v6=-xc2+xc6-l2(6);
    v7=-xc3+xc4-l2(7);
    v8=-xc3+xc5-l2(8);
    v9=-xc5+xc6-l2(9);
    
    V=[v1 v2 v3 v4 v5 v6 v7 v8 v9 ];%改正数矩阵
    V=V'
    
    
    Xc=x0+xcl %参数平差值
    
    
    h=h';
    H=h+V %高程平差值

[0cfde11838e744dc882c20090f4eca9a.png]: /images/20220828/38e26ebce19d48c99553a4af567a2c88.png
[9d463def9da54d2cbc2cd60ebc3bae43.png]: /images/20220828/69799aaf146641fab046b29fd00706a1.png
[b5120b3aa7dd4b97bf3af1fb97f663dd.png]: /images/20220828/5f5ea30d5d9743bba125130b6561f262.png
[efc926f1e9954ae69b1746ea094d44df.png]: /images/20220828/fbfd673398a24566b431b82f13c38424.png
[bc8302f03d3147268b86e3c934a77ca2.png]: /images/20220828/ddfa73c958b54b7798942db5c4ac216b.png
[80ae7a8340a642bc85975a781e752b0c.png]: https://img-blog.csdnimg.cn/80ae7a8340a642bc85975a781e752b0c.png
[watermark_type_ZHJvaWRzYW5zZmFsbGJhY2s_shadow_50_text_Q1NETiBAbm9jaGVuZ3pp_size_11_color_FFFFFF_t_70_g_se_x_16]: https://img-blog.csdnimg.cn/086dce5455ca4262a7001efb9ae116e0.png?x-oss-process=image/watermark,type_ZHJvaWRzYW5zZmFsbGJhY2s,shadow_50,text_Q1NETiBAbm9jaGVuZ3pp,size_11,color_FFFFFF,t_70,g_se,x_16
[watermark_type_ZHJvaWRzYW5zZmFsbGJhY2s_shadow_50_text_Q1NETiBAbm9jaGVuZ3pp_size_20_color_FFFFFF_t_70_g_se_x_16]: https://img-blog.csdnimg.cn/b26d57109bb74efe937d55ebdf73d93b.png?x-oss-process=image/watermark,type_ZHJvaWRzYW5zZmFsbGJhY2s,shadow_50,text_Q1NETiBAbm9jaGVuZ3pp,size_20,color_FFFFFF,t_70,g_se,x_16
[351e6f1c4ba04391813daefb636f7978.png]: /images/20220828/e0507bf0a212464f8b6e400daf23fc5a.png
[watermark_type_ZHJvaWRzYW5zZmFsbGJhY2s_shadow_50_text_Q1NETiBAbm9jaGVuZ3pp_size_20_color_FFFFFF_t_70_g_se_x_16 1]: /images/20220828/2fe0ff5dcd394b56aaccbc7be811aab3.png