条件平差观测方程A*L_hatA_00%误差方程A*VW0法方程Naa*KW0间接平差观测方程L_hatB*X_hatd%误差方程VB*x_hat-l法方程Nbb*x_hat-B*P*l0clc; clear; %% 1. 基本观测数据输入 % 已知点高程 HA 5.016; HB 6.016; % 观测高差单位m L [ 1.359; 2.009; 0.363; 1.012; 0.657; 0.238; -0.595 ]; % 路线长度单位km S [ 1.1; 1.7; 2.3; 2.7; 2.4; 1.4; 2.5 ]; n length(L); % 观测数 %% 2. 系数矩阵 B 和常数项 d对应图中 hat{h}_i BX d B [ 1 0 0; % h1 X1 - HA 0 1 0; % h2 X2 - HA 1 0 0; % h3 X1 - HB 0 1 0; % h4 X2 - HB -1 1 0; % h5 -X1 X2 -1 0 1; % h6 -X1 X3 0 0 -1 % h7 -X3 HB ]; % 保留图1中的常数项注意单位一致 d [ -HA; % h1 -HA; % h2 -HB; % h3 -HB; % h4 0; % h5 0; % h6 HB % h7 ]; %% 3. 权阵 P diag(1/S) P diag(1 ./ S); % 对角阵S 是路线长度权 ∝ 1/距离 %% 4. 初始高程估值 X0根据图中定义 X0 [ HA L(1); % X1^0 HA h1 HA L(2); % X2^0 HA h2 HA L(1) L(6) % X3^0 HA h1 h6 ]; %% 5. 构造法方程右边常数项l L - (B*X0 d) l L - (B * X0 d); % 观测值 - 模型值 l l * 1000; % 单位转为 mm增强数值精度 %% 6. 解法方程求改正数 x_hat并更新平差值 Nbb B * P * B; x_hat inv(Nbb) * B * P * l; X_hat X0 x_hat / 1000; % 记得转回 m %% 7. 改正数 V 和观测值平差值 L_hat V B * x_hat - l; % 改正数单位 mm L_hat B * X_hat d; % 平差观测值单位 m %% 8. 输出结果 fprintf( 间接平差结果 \n); fprintf(未知点高程平差值 X_hat:\n); fprintf(H_P1 %.4f m\n, X_hat(1)); fprintf(H_P2 %.4f m\n, X_hat(2)); fprintf(H_P3 %.4f m\n, X_hat(3)); fprintf(\n观测值改正数 V单位mm\n); disp(V); fprintf(平差后观测值 L_hat单位m\n); disp(L_hat); %% 9. 计算单位权中误差单位 mm r n - length(X0); % 自由度 观测数 - 参数数 sigma0_hat sqrt((V * P * V) / r); % 单位权中误差mm %% 10. 协方差阵 参数中误差单位 m Qxx inv(Nbb); % 协方差矩阵单位 mm^2 m_X sqrt(diag(Qxx)) * sigma0_hat / 1000; % 转成 m %% 11. 观测值中误差单位 mm Qvv B * Qxx * B; % 改正数协方差阵 m_L sqrt(diag(Qvv)) * sigma0_hat; %% 12. 输出精度信息 fprintf(\n 精度评定 \n); fprintf(单位权中误差 sigma0_hat %.4f mm\n, sigma0_hat); fprintf(\n未知点高程中误差单位m\n); fprintf(m_H_P1 %.4f m\n, m_X(1)); fprintf(m_H_P2 %.4f m\n, m_X(2)); fprintf(m_H_P3 %.4f m\n, m_X(3)); fprintf(\n观测值中误差单位mm\n); disp(m_L);