简介一套覆盖水准网间接平差案例的Python、C与MATLAB代码实现资源面向测绘工程、地理信息系统及测量数据处理学习者通过三种语言完整演示基于最小二乘原理的水准网间接平差流程包括高差观测值组织、未知点高程参数化、误差方程建立与法方程求解等关键步骤。压缩包共3个文件分别对应Python脚本.py、MATLAB脚本.m与C程序.cpp代码均从原始数据出发涵盖数据预处理、最小二乘优化、残差计算与平差后高程输出等环节结构清晰便于直接运行和逐段比对。资源整体仅3KB轻量聚焦适合具备一定编程基础、希望快速上手水准网平差实现的在校生或工程技术人员。目前已有1044人学习下载既可用于课程实验和毕业设计验证也可作为调用scipy.optimize、Eigen库和MATLAB优化工具箱的拓展范例直观感受三种语言在矩阵运算、内存管理和数值计算效率上的差异深化对间接平差核心原理的理解。1. 水准网间接平差从高差观测到最小二乘求解水准测量外业结束后最磨人的往往不是算高差而是闭合差超限后把每条路线在手簿上反复校核。水准网间接平差就是把“已知点高程 观测高差 待定点”整张网丢进一个线性模型用最小二乘同时求出全部待定点高程的最优估值并根据协因数阵给出每个点的高程精度。这个模型覆盖三等、四等加密工程沉降监测甚至高速铁路 CPIII 高程网的概算。本文用同一组案例数据给出 Python、C、Matlab 三套完整实现覆盖误差方程、法方程求解、定权和精度评定也会把初学者最容易出错的点号索引和权矩阵定义逐一挑明。2. 误差方程与法方程水准网间接平差的数学模型水准测量拿到的高差是直接观测值但最终需要的是点位高程。间接平差把每个待定点的高程设为参数把每一条观测高差都写成这些参数的线性函数然后一次性解算。与条件平差要为每个闭合环单独列条件方程不同间接平差的误差方程列立规则全局统一换一个网形、加一条观测边程序逻辑完全不用动。这也是市面上绝大部分平差软件内部采用间接平差的原因。2.1 从观测高差到误差方程设某条观测边的起点为 s终点为 e观测高差为 h(s,e)起终点高程平差值为 X_e、X_s理论上存在关系h(s,e) v X_e - X_s其中 v 是该观测值的改正数也就是残差。X 本身是未知数直接列方程不便求解因此先给每个待定点一个近似高程 X0代入 X X0 δx整理得到如下误差方程v (δx_e - δx_s) - l这里常数项l h(s,e) - (X0_e - X0_s)。若起点或终点是已知点对应的 δx 直接取 0那一项就消失。把所有 n 条观测边竖向堆叠得到矩阵形式V B * δx - lB 是一个 n 行、t 列的矩阵t 为未知点个数。由于水准方程是终点高程减起点高程B 的每一行只会出现 0、1、-1 三个值。这里有一个值得注意的性质近似高程 X0 怎么选不改变最终平差结果。哪怕某条推算路线选得粗糙l 向量会自动吸收它与观测值之间的差异。实际工程中直接从已知点沿任意可达路线推算近似高程即可程序对此非常宽容。符号含义维度n观测高差条数标量t未知点个数必要观测数标量r n - t多余观测数标量B误差方程系数矩阵n × tl误差方程常数项向量n × 1P观测权矩阵n × nδx近似高程改正数向量t × 1V高差改正数残差向量n × 12.2 最小二乘原理与法方程间接平差采用的准则是加权残差平方和最小即目标函数V^T * P * V min对 δx 求导并令结果为零得到法方程N * δx W其中N B^T * P * BW B^T * P * lN 是 t × t 对称正定矩阵。只要每个未知点都有观测边连入网中N 就非奇异解唯一。几十个点的工程网直接用稠密矩阵求解即可上千点的控制网则应改用稀疏 Cholesky 分解。顺带说明必要观测的概念t 就是水准网的必要观测数一个未知点至少要有一条观测边把它与已知部分连接。多余观测数 r n - t闭合环和附合路线带来冗余r 越大平差结果的检验能力越强单位权中误差也越有统计意义。若 r 0各边只是逐段推算不存在平差问题。2.3 权的确定与精度评定水准测量最常用的定权方式是按路线长度P_i c / s_i其中 s_i 是第 i 段水准路线长度单位 km也可按测站数定权P_i c / n_i。用什么权取决于外业记录里哪个量更可信。需要澄清的是常数 c 取任何正数都不改变 δx 的解但它会整体缩放 P进而改变 V^T P V 与 σ0 的数值。因此成果报告在写“单位权中误差”时必须说清单位权定义例如“以 1 km 路线观测高差为单位权”否则该精度数字没有可比性。整体解算步骤如下确定未知点近似高程 X0对每条边计算常数项 l。按选定的权模型生成对角权矩阵 P。组成法方程并求解 δx。计算 X X0 δx回代得 V再计算 σ0 与各点高程精度。精度评定公式为σ0 sqrt(V^T * P * V / r)第 j 个未知点的高程中误差为σ_j σ0 * sqrt(Q_jj)其中Q N^-1为协因数阵Q_jj 是该点对应的对角元素。最终成果表里每个点高程后面的正负号数字就是这个公式算出来的。3. Python 实现间接平差用 NumPy 搭出完整解算流程Python 实现平差的最大好处是能把教科书公式逐一对应到 NumPy 调用上。环境要求很简单Python 3.10 及以上版本安装 NumPy 即可。下面用一个 5 条观测边、2 个未知点的水准网把完整流程走一遍这个网闭合环多但规模小每一行代码都能对应到上一章的公式。3.1 一个三口闭合环的水准网案例数据已知点 A 高程 10.000 mB 高程 12.000 m待定点为 C、D。5 条观测边构成两个闭合环A-C-D-A 和 A-C-B不边 3 与边 5 在 C、D、B 处形成交叉条件。数据如下边号起点终点高差观测值 h (m)路线长度 s (km)1AC1.2351.22CD2.4550.83DB-1.6901.04AD3.6881.55CB0.7661.1近似高程用边 1 推得 C0 11.235 m用边 1 边 2 推得 D0 13.690 m。边 4 对应闭合路径 A-C-D-A闭合差为 1.235 2.455 - 3.688 0.002 m边 5 与边 2、3 构成另一组闭合关系。两个闭合环都在毫米量级适合做平差演示。数据为示例值但采用的定权方式与解算流程与生产环境一致。3.2 构造误差方程与解算法方程的核心代码下面是可直接运行的 Python 脚本。点号直接用 0、1、2、3 对应 A、B、C、D未知点 C、D 在参数列表中的逻辑位置是 2、3。import numpy as np # 已知点与近似高程A0, B1, C2, D3 H0 np.array([10.0, 12.0, 11.235, 13.690]) known {0, 1} # 已知点集合不参与参数解算 # 观测边: (起点, 终点, 高差观测值, 路线长度 km) edges [ (0, 2, 1.235, 1.2), (2, 3, 2.455, 0.8), (3, 1, -1.690, 1.0), (0, 3, 3.688, 1.5), (2, 1, 0.766, 1.1), ] n len(edges) t len(H0) - len(known) # 未知参数个数C、D 共 2 个 B np.zeros((n, t)) # 误差方程系数矩阵 l np.zeros(n) # 常数项 P np.eye(n) # 观测权矩阵先初始化为等权 for i, (s, e, h, dist) in enumerate(edges): l[i] h - (H0[e] - H0[s]) P[i, i] 1.0 / dist # 按距离定权与第 2 章公式对应 # 起点对参数列贡献 -1终点贡献 1已知点对应的 δx 为 0 if s not in known: B[i, s - 2] -1.0 # 点号 2 - 参数列 0点号 3 - 参数列 1 if e not in known: B[i, e - 2] 1.0 # 组成法方程并求解 N B.T P B W B.T P l dx np.linalg.solve(N, W) # 平差值、残差与精度 X H0.copy() X[2:] dx V B dx - l sigma0 np.sqrt(V P V / (n - t)) Q np.linalg.inv(N) for idx, name in enumerate([C, D]): sigma sigma0 * np.sqrt(Q[idx, idx]) print(f{name} {X[2 idx]:.4f} m, sigma {sigma:.4f} m) print(fsigma0 {sigma0:.4f} m)B 矩阵的构造是本段代码的核心逻辑起点是未知点时该行列填 -1终点是未知点时填 1已知点那端因改正数为 0 不参与。l 的计算符号与第 2 章推导一致即观测高差减近似高差。权矩阵取对角线元素 1/s等于 1 km 观测高差被定义为单位权。np.linalg.solve直接解 N、W比先求逆再相乘更稳定平差结果中 V 反映每条边被“调整”的量可用于后续粗差筛查。按这套数据运行输出 C11.2344 m、D13.6893 mσ0≈0.83 mm。3.3 等权与按距离定权改一个参数看结果变化把上面代码中的P[i, i] 1.0 / dist改成P[i, i] 1.0就退化为普通最小二乘。等权解算会得到 C11.2344、D13.6893数值与定权结果相差不到 0.1 mm但 σ0 的含义完全不同一个表示“1 km 路线观测高差”的中误差一个表示“这条观测边”的中误差。这个例子说明权的选择直接影响“最优”的含义也影响精度指标的读数。实际水准网规范通常明确要求按距离或测站数定权。如果需要从 CSV 读取外业数据只需把上面的 edges 初始化替换为如下片段即可列顺序按起点、终点、高差、距离import csv edges [] with open(leveling.csv) as f: for row in csv.reader(f): s, e, h, dist map(float, row) edges.append((int(s), int(e), h, dist))4. C 实现间接平差矩阵手动实现与编译运行C 没有 NumPy 这类默认矩阵库程序化水准网平差需要自己处理矩阵运算。小规模工程网建议用自写高斯消元几十个点时完全够用项目里需要更高性能或更大网形再上 Eigen。下面的完整程序不依赖任何第三方库结构上与 Python 版本保持一致方便逐行对照。4.1 数据结构与高斯消元求解点号同样约定 0A、1B、2C、3D。程序输出 X3、X4即 C、D 两点的高程解算结果。#include iostream #include vector #include cmath using namespace std; struct Edge { int s, e; double h, dist; }; // 列主元高斯消元求解 A * x b vectordouble gauss(vectorvectordouble A, vectordouble b) { int n A.size(); for (int col 0; col n; col) { int p col; for (int r col 1; r n; r) if (fabs(A[r][col]) fabs(A[p][col])) p r; swap(A[col], A[p]); swap(b[col], b[p]); for (int r col 1; r n; r) { double f A[r][col] / A[col][col]; for (int c col; c n; c) A[r][c] - f * A[col][c]; b[r] - f * b[col]; } } vectordouble x(n); for (int r n - 1; r 0; --r) { x[r] b[r]; for (int c r 1; c n; c) x[r] - A[r][c] * x[c]; x[r] / A[r][r]; } return x; } int main() { const int n_known 2; // 已知点 A、B vectordouble H0 {10.0, 12.0, 11.235, 13.690}; vectorEdge edges { {0, 2, 1.235, 1.2}, {2, 3, 2.455, 0.8}, {3, 1, -1.690, 1.0}, {0, 3, 3.688, 1.5}, {2, 1, 0.766, 1.1} }; int n edges.size(); int t H0.size() - n_known; // 未知点数量 vectorvectordouble B(n, vectordouble(t, 0.0)); vectordouble l(n, 0.0), P(n, 1.0); vectordouble W(t, 0.0); vectorvectordouble N(t, vectordouble(t, 0.0)); for (int i 0; i n; i) { const Edge e edges[i]; l[i] e.h - (H0[e.e] - H0[e.s]); P[i] 1.0 / e.dist; // 按距离定权 if (e.s n_known) B[i][e.s - n_known] -1.0; if (e.e n_known) B[i][e.e - n_known] 1.0; // 直接累加组成法方程 N、W for (int a 0; a t; a) { W[a] B[i][a] * P[i] * l[i]; for (int b 0; b t; b) N[a][b] B[i][a] * P[i] * B[i][b]; } } vectordouble dx gauss(N, W); vectordouble X H0; for (int k 0; k t; k) X[k n_known] dx[k]; double VtPV 0.0; for (int i 0; i n; i) { double v 0.0; for (int a 0; a t; a) v B[i][a] * dx[a]; v - l[i]; VtPV v * P[i] * v; } double sigma0 sqrt(VtPV / (n - t)); const char* names[] {C, D}; for (int k 0; k t; k) cout names[k] X[k n_known] m\n; cout sigma0 sigma0 m\n; return 0; }代码中 Edge 结构体直接使用 H0 的下标省去字符串点号解析。组法方程时采用双重循环累加而不是先构造 B^T P 再乘 B避免临时对象代码也更直观。P 在这里是一维数组本质上就是对角权矩阵的紧凑表示。注意if (e.s n_known)这个条件它等价于判断起点是否未知点列索引用点号减 2与已知点的数量绑定。如果你的已知点编号不是从 0、1 开始这里要换成点号到参数列的映射表。编译运行方式如下g -stdc17 adjust.cpp -o adjust ./adjustWindows 下可以在 VS Code 配置 C/C 环境选择 g 作为编译器后在终端直接编译若使用 MSVC 工具链则注意安装对应架构的 Visual C 运行库。程序输出结果与 Python 版本一致C、D 高程分别约为 11.2344 m 与 13.6893 m。4.2 用 Eigen/LDLT 替换自写求解器当网形扩大自写高斯消元的平方级增长会变得明显。Eigen 是 C 生态中最常用的线性代数库面对对称正定法方程LDLT 分解比通用高斯消元更快、更省内存。改造时只需将 N、W 填入 Eigen 矩阵与向量然后调用一次求解#include Eigen/Dense Eigen::MatrixXd NE(t, t); Eigen::VectorXd WE(t); // 填充 NE、WE内容与上面的 N、W 累加循环一致 Eigen::VectorXd dx NE.ldlt().solve(WE);LDLT 不显式求逆数值稳定性好适合法方程这类对称正定结构。实际生产环境中若观测边达到几万条、未知点上千个还应该配合稀疏矩阵存储但本标题的示例网用不到这一步。5. Matlab 实现间接平差矩阵语法的直接映射Matlab 的优势在于矩阵运算被当作一等公民平差脚本几乎能把数学公式逐行翻译成语句。R2018a 及之后各版本运行以下脚本均无兼容问题。与 Python、C 版本相比Matlab 版代码量最少适合课程设计、算法原型和快速验证。5.1 完整脚本与运行结果Matlab 下标从 1 开始因此点号换算需要多一步。脚本仍以 0A、1B、2C、3D 定义边数据转成 Matlab 下标时统一加 1% 水准网间接平差C、D 为待定点 H0 [10.0; 12.0; 11.235; 13.690]; % A、B 已知C、D 近似 edges [0 2 1.235 1.2; ... % 每行: 起点 终点 h s 2 3 2.455 0.8; ... 3 1 -1.690 1.0; ... 0 3 3.688 1.5; ... 2 1 0.766 1.1]; n size(edges, 1); t length(H0) - 2; % C、D 两个未知点 B zeros(n, t); l zeros(n, 1); P eye(n); for i 1:n s edges(i, 1) 1; % 转成 Matlab 1 基下标 e edges(i, 2) 1; l(i) edges(i, 3) - (H0(e) - H0(s)); if s 2, B(i, s-2) -1; end % 起点为未知点 - -1 if e 2, B(i, e-2) 1; end % 终点为未知点 - 1 P(i, i) 1 / edges(i, 4); % 按距离定权 end N B * P * B; W B * P * l; dx N \ W; % 左除等价 N^-1 * W X H0; X(3:end) X(3:end) dx; V B * dx - l; sigma0 sqrt(V * P * V / (n - t)); Q inv(N); for k 1:t fprintf(X%d %.4f m, sigma %.4f m\n, k2, X(k2), sigma0*sqrt(Q(k,k))); endMatlab 版与前面两种语言有三个容易踩的差异点。其一B(i, s-2) 里的索引换算依赖未知点编号恰好从点号 2 开始若换网形建议单独建立点号到列的字典。其二N \ W是左除等效于 inv(N) * W但实际走 LU 或 Cholesky 分解速度与稳定性更好。其三inv(N) 在这个规模下没有问题但未知点很多时建议换成Q N \ eye(t)既算协因数阵又避免显式求逆。运行结果与 Python、C 一致输出 X311.2344、X413.6893单位权中误差约 0.83 mm。X3 即 C 点、X4 即 D 点。5.2 三种语言实现对照与选型建议对比项PythonCMatlab矩阵基础NumPy自写或 Eigen内置矩阵类型代码量短最长最短部署环境需 Python 运行时可编译为独立可执行文件需 MATLAB 环境适用阶段数据处理流程、沉降监测自动化采集软件集成、生产系统教学、算法原型、可视化选型不只看语言本身。Python 适合把平差写进自动化流程例如自动读取 CSV、批量计算并生成报告C 适合嵌入采集程序或无人值守服务不依赖解释器启动快且可控Matlab 适合验证新算法、快速画残差图和分析精度。三种实现对同一份数据输出结果一致本身就是一次很好的交叉验证。6. 平差结果的校验与粗差定位技巧平差算完不能直接写进成果表数值正确性要靠闭合差和残差双重检核。闭合差反映整张网的内部一致性残差则能把问题定位到具体某一条观测边。6.1 闭合差、残差与 3σ0 筛选闭合环闭合差在外业阶段就应该初步检查平差后还要再看一遍。以 A-C-D-A 环为例edges[0][2] edges[1][2] - edges[3][2]就是它的闭合差本例约 2 mm。若一条边在多个闭合环中反复出现超限这条边往往就是问题所在。平差后的残差 V 是更细的定量指标逐边计算 v/σ0 并排序能直接给出可疑观测列表limit 3 * sigma0 for i, v in enumerate(V): ratio abs(v) / sigma0 flag -- suspicious if ratio 3 else print(fedge {i1}: v {v*1000:6.2f} mm, |v|/sigma0 {ratio:5.2f}{flag})正常情况下残差应当在 ±2σ0 范围内随机分布。出现 |v| 3σ0 时先查外业手簿与电子记录确认仪器、尺垫、转点设置是否有误再决定是否复测而不是直接删观测值。6.2 可疑观测的处理原则与重算对比剔除观测边会改变网形必须保证剩余观测仍然把所有未知点连接到已知部分。删边后重新平差如果 σ0 显著下降且某个未知点高程变化超过 2 倍该点中误差说明这条边对整个网形的拉力很强处理时需要格外谨慎。对真正确认有误的观测复测值参与替换后再平差流程与最初完全一致只需改一条数据重新运行一次脚本。把每次重算得到的高程与中误差、σ0 保存在一张表里观察变化趋势比单次结果更有说服力。上面那段筛选脚本输出的第一行就是下一个外业工作日的复测点名。本文还有配套的精品资源点击获取