GNSS网平差-间接平差法（附Python和Matlab代码）

一、实验目的与要求

1. 掌握GNSS网无约束平差的原理；

2. 理解间接平差原理；

3. 掌握GNSS网无约束平差计算步骤。

二、实验内容、步骤及结果

1. 根据GNSS网，确定平差方法（约束平差还是无约束平差，或与地面数据联合平差），因只给一个已知点坐标，采用无约束平差；
2. 确定参数t和间接平差方程数；
3. 列出间接平差方程；
4. 计算P阵；
5. 写出间接平差计算过程和中间结果
6. 计算出参数（未知点坐标）。

Python代码

import numpy as np
from scipy.linalg import block_diag


# load covariance and calculate P
D1=np.mat([[2.32099900000000e-07,-5.09700800000000e-07,-4.37140100000000e-07],[-5.09700800000000e-07,1.33993100000000e-06,1.10935600000000e-06],[-4.37140100000000e-07,1.10935600000000e-06,1.00859200000000e-06]])
D2=np.mat([[1.04489400000000e-06,-2.39653300000000e-06,-2.31968300000000e-06],[-2.39653300000000e-06,6.34129100000000e-06,5.902876E-06],[-2.31968300000000e-06,5.902876E-06,6.03557700000000e-06]])
D3=np.mat([[5.85006400000000e-07,-1.32962000000000e-06,-1.25237400000000e-06],[-1.32962000000000e-06,3.36254800000000e-06,3.06982000000000e-06],[-1.25237400000000e-06,3.06982000000000e-06,3.01923300000000e-06]])
D4=np.mat([[1.20531900000000e-06,-2.63670200000000e-06,-2.17410600000000e-06],[-2.63670200000000e-06,6.85858500000000e-06,5.48074500000000e-06],[-2.17410600000000e-06,5.48074500000000e-06,4.82012500000000e-06]])
D5=np.mat([[9.66265700000000e-06,-2.17547600000000e-05,-1.97146800000000e-05],[-2.17547600000000e-05,5.19477700000000e-05,4.63356500000000e-05],[-1.97146800000000e-05,4.63356500000000e-05,4.32411000000000e-05]])
D=block_diag(D1,D2,D3,D4,D5)
P=np.linalg.inv(D/(0.00298**2))

# matrix of coefficients
_ones=-np.diag([1,1,1])
ones=np.diag([1,1,1])
zeros=-np.diag([0,0,0])
B=np.r_[
np.c_[_ones,zeros,zeros],
np.c_[zeros,zeros,-ones],
np.c_[ones,zeros,-ones],
np.c_[ones,-ones,zeros],
np.c_[zeros,ones,-ones]
]
l=np.mat([-0.00100000000000000,0.000700000000000000,0.00150000000000000,-0.00700000000000000,0.0140000000000000,0.001150000000000000,-0.0110000000000000,0.0243000000000000,0.0250000000000000,0.00400000000000000,-0.00970000000000000,-0.0120000000000000,0,0,0],
).transpose()

# Solving the normal equation
Nbb=np.matmul(np.matmul((B.transpose()),P),B)
W=np.matmul(np.matmul((B.transpose()),P),l)
# W=(B.transpose())*P*l
x=np.matmul(np.linalg.inv(Nbb),W)
# x=np.linalg.inv(Nbb)*W

# Calculate the result
X0Y0Z0=np.mat([-1973420.17400000,4591054.04670000,3951407.20500000, -1974825.70100000,4591232.19400000,3950235.81300000, -1974909.19800000,4590518.04100000,3951265.01200000]).transpose()
X0Y0Z0_=X0Y0Z0+x
print("The final result of adusting")
print(X0Y0Z0_)

Matlab代码

clc;clear;

%% load covariance and calculate P
D1=[2.32099900000000e-07,-5.09700800000000e-07,-4.37140100000000e-07;-5.09700800000000e-07,1.33993100000000e-06,1.10935600000000e-06;-4.37140100000000e-07,1.10935600000000e-06,1.00859200000000e-06];
D2=[1.04489400000000e-06,-2.39653300000000e-06,-2.31968300000000e-06;-2.39653300000000e-06,6.34129100000000e-06,5.902876E-06;-2.31968300000000e-06,5.902876E-06,6.03557700000000e-06];
D3=[5.85006400000000e-07,-1.32962000000000e-06,-1.25237400000000e-06;-1.32962000000000e-06,3.36254800000000e-06,3.06982000000000e-06;-1.25237400000000e-06,3.06982000000000e-06,3.01923300000000e-06];
D4=[1.20531900000000e-06,-2.63670200000000e-06,-2.17410600000000e-06;-2.63670200000000e-06,6.85858500000000e-06,5.48074500000000e-06;-2.17410600000000e-06,5.48074500000000e-06,4.82012500000000e-06];
D5=[9.66265700000000e-06,-2.17547600000000e-05,-1.97146800000000e-05;-2.17547600000000e-05,5.19477700000000e-05,4.63356500000000e-05;-1.97146800000000e-05,4.63356500000000e-05,4.32411000000000e-05];
D = blkdiag(D1,D2,D3,D4,D5);
P=(D/(0.00298^2))^(-1);
P1=(D1/(0.00298^2))^(-1);
%% matrix of coefficients
one3=ones(1,3);
B=[-diag(one3) zeros(3) zeros(3);
    zeros(3) zeros(3) -diag(one3);
    diag(one3) zeros(3)  -diag(one3);
    diag(one3) -diag(one3) zeros(3);
    zeros(3) diag(one3) -diag(one3)];
l=[-0.00100000000000000;0.000700000000000000;0.00150000000000000;-0.00700000000000000;0.0140000000000000;0.001150000000000000;-0.0110000000000000;0.0243000000000000;0.0250000000000000;0.00400000000000000;-0.00970000000000000;-0.0120000000000000;0;0;0];

%% Solving the normal equation
Nbb=B'*P*B;
W=B'*P*l;
x=Nbb^(-1)*W;

%% accuracy assessment
V=B*x-l;
sigma0=sqrt((V'*P*V)/(15-9));
Q=inv(P);
disp("The accuracy");
sigma=sigma0*(diag(Q)).^0.5

%% Calculate the result
X0Y0Z0=[-1973420.17400000,4591054.04670000,3951407.20500000 -1974825.70100000,4591232.19400000,3950235.81300000 -1974909.19800000,4590518.04100000,3951265.01200000]';
XYZ=[-1218.56100000000,-1039.22700000000,1737.72000000000 270.457000000000,-503.208000000000,1879.92300000000 1489.01300000000,536.030000000000,142.218000000000 1405.53100000000,-178.157000000000,1171.38000000000 83.4970000000000,714.153000000000,-1029.19900000000]';
XYZ_=XYZ+V;
disp("The result of adjusting");
X0Y0Z0_=X0Y0Z0+x