简介面向测绘、GIS与工程测量方向的开发者和学生这份资源包提供了水准网间接平差的Python、C、MATLAB三种实现代码。包内共3个文件分别对应.py、.m、.cpp各一个整体仅3KB体积小巧便于快速阅读和移植。已有1044人学习下载适合需要理解间接平差建模与最小二乘求解的入门及进阶者。代码覆盖了水准网平差的核心流程先对观测高差进行异常值检测等预处理再依据水准测量原理构建线性观测方程随后调用相应优化算法完成最小二乘参数求解最后输出平差高程并评估精度。Python版依托numpy和scipy实现简洁科学计算C版基于Eigen库并使用Levenberg-Marquardt迭代注重性能MATLAB版借助lsqnonlin等内置函数代码直观且便于可视化验证。通过对比三种实现既能加深对间接平差数学模型的理解也能快速将代码改造迁移至实际工程或毕业设计项目中。1. 水准网间接平差的 Python、C、MATLAB 三语言实现不只是演示水准网间接平差是测量平差里最容易“看似会、上手算错”的模型。它把重合路线上互差几毫米的高差观测值放进同一个最小二乘模型最后得到一组唯一的未知点高程和一整套精度指标。用 Python、C、MATLAB 各写一遍不是为了重复造轮子而是因为生产环境里三套代码要互相验证Python 适合出图和快速建模C 适合嵌入采集程序MATLAB 适合做教学和报告。下面按同一个水准网案例把建模、法方程、结果校验完整走一遍换成自己的路线数据时只需要改观测表和已知点高程。2. 水准网间接平差的观测方程与法方程符号定对三套代码才能复现间接平差又叫参数平差关键动作只有一个把每条边的高差观测值表示成未知点高程的线性函数。水准网不像 GPS 网需要先给近似坐标它是一个线性模型观测方程列对符号后面三套代码结果自然一致。2.1 高差观测怎样变成设计矩阵 A设已知点 A 的高程是 100.000 m未知点 P1、P2、P3 的高程分别是 X1、X2、X3。以路线起点到终点的高差 h 为准平差后的关系是h v H_to - H_from对 A→P1 这条边方程是v X1 - (100 h)对 P1→P2 这条边方程是v X2 - X1 - h。把每个未知点设为一列上面两条边的设计矩阵行就是[1, 0, 0]和[-1, 1, 0]。序号起点终点高差 h (m)距离 (km)1AP11.2561.02P1P2-0.5820.83P2P31.7641.24AP20.6721.55P1P31.1832.06AP32.4362.57P3A-2.4411.8这张表故意包含了两个闭合环和 4 个多余观测方便验证协因数阵和单位权中误差。构造设计矩阵时最核心的代码只有一段分支# 只列行系数和常数项完整程序见第 3 章 if f in known and t not in known: row[t - 1] 1.0 const known[f] h elif f not in known and t in known: row[f - 1] -1.0 const h - known[t] else: row[f - 1] -1.0 row[t - 1] 1.0 const h这里的常数项最容易写错。已知点推向未知点时未知点高程等于已知高程加实测高差所以用known[f] h当常数反方向时常数是h - known[t]两端都是未知点时常数就是h。代码里的row·X - const是残差表示同一时刻的观测改正数。2.2 权阵采用路线长度倒数水准测量误差近似随路线长度累积常用做法是取P diag(1 / L)L 以 km 为单位。短边权重高长边权重低。现场可获得的先验信息权路线长度 L (km)w 1 / L测站数 nw 1 / n先验高差中误差 σ_hw 1 / σ_h²如果观测记录里只有测站数把1 / dist换成1 / n_station即可如果有仪器标称精度也只要改权这一行。已知点之间的高差观测不要放进法方程它没有未知参数参与平差只会干扰自由度统计应单独做检核。2.3 法方程、单位权方差和协因数阵把残差写成一列V A X - d目标是最小化0.5 * Vᵀ P V。对 X 求导并令其为零得到法方程N Aᵀ P Arhs Aᵀ P dN X rhs量公式作用法方程系数N Aᵀ P A阶数等于未知点数常数向量rhs Aᵀ P d由高差和已知高程共同决定未知点高程X N⁻¹ rhs平差结果残差V A X - d观测改正数单位权方差σ₀² Vᵀ P V / (m - u)整体精度协因数阵Qxx N⁻¹用于计算中误差这里u是未知点个数不是观测数。本案例有 7 条观测边、3 个未知点自由度是7 - 3 4单位权中误差才有实际意义。若网形没有多余观测或没有固定已知点法方程会奇异或无法估计精度。3. Python、C、MATLAB 三套水准网间接平差代码实现三种语言用同一张观测表、同一个权重约定。下面代码不是伪码是能直接跑的最小实现。3.1 Python NumPyimport numpy as np # (from, to, dh, dist_km) obs [ (0, 1, 1.256, 1.0), (1, 2, -0.582, 0.8), (2, 3, 1.764, 1.2), (0, 2, 0.672, 1.5), (1, 3, 1.183, 2.0), (0, 3, 2.436, 2.5), (3, 0, -2.441, 1.8), ] known {0: 100.000} # 节点 0 为已知点 A u 3 # 未知点 P1, P2, P3 A, d, w [], [], [] for f, t, h, dist in obs: row np.zeros(u) if f in known and t not in known: row[t - 1] 1.0 const known[f] h elif f not in known and t in known: row[f - 1] -1.0 const h - known[t] else: row[f - 1] -1.0 row[t - 1] 1.0 const h A.append(row) d.append(const) w.append(1.0 / dist) A np.array(A) P np.diag(w) d np.array(d) N A.T P A rhs A.T P d X np.linalg.solve(N, rhs) V A X - d sigma0_2 V P V / (len(obs) - u) Qxx np.linalg.inv(N) print(X:, X) print(V:, V) print(sigma0:, np.sqrt(sigma0_2)) print(std:, np.sqrt(np.diag(sigma0_2 * Qxx)))Python 的关键是用np.linalg.solve而不是inv(N) rhs后者在法方程接近病态时会放大误差只有需要协因数阵时才显式取逆。P在这里是u known之外的对角矩阵水准网观测通常只有几十到几百条直接构造稠密对角阵没有问题。换路线时只需改obs列表known字典支持多个已知点。3.2 C 标准库实现C 实现不需要第三方库核心是自写高斯消元。小网的法方程是稠密u × u标准库足够。#include vector #include algorithm #include cmath #include iostream #include iomanip using Vec std::vectordouble; using Mat std::vectorVec; Vec gauss(Mat A, Vec b) { int n (int)b.size(); for (int col 0; col n; col) { int r0 col; for (int i col 1; i n; i) if (std::abs(A[i][col]) std::abs(A[r0][col])) r0 i; std::swap(A[col], A[r0]); std::swap(b[col], b[r0]); for (int i col 1; i n; i) { double t A[i][col] / A[col][col]; for (int j col; j n; j) A[i][j] - t * A[col][j]; b[i] - t * b[col]; } } Vec x(n); for (int i n - 1; i 0; --i) { x[i] b[i]; for (int j i 1; j n; j) x[i] - A[i][j] * x[j]; x[i] / A[i][i]; } return x; } struct Obs { int from, to; double dh, dist; }; int main() { const int m 7, n 3; std::vectorObs obs { {0, 1, 1.256, 1.0}, {1, 2, -0.582, 0.8}, {2, 3, 1.764, 1.2}, {0, 2, 0.672, 1.5}, {1, 3, 1.183, 2.0}, {0, 3, 2.436, 2.5}, {3, 0, -2.441, 1.8} }; Mat A(m, Vec(n, 0.0)); Vec d(m), w(m), X(n), V(m); double H0 100.0; auto is_known [](int k) { return k 0; }; for (int k 0; k m; k) { int f obs[k].from, t obs[k].to; double h obs[k].dh; if (is_known(f) !is_known(t)) { A[k][t - 1] 1.0; d[k] H0 h; } else if (!is_known(f) is_known(t)) { A[k][f - 1] -1.0; d[k] h - H0; } else { A[k][f - 1] -1.0; A[k][t - 1] 1.0; d[k] h; } w[k] 1.0 / obs[k].dist; } Mat N(n, Vec(n, 0.0)); Vec rhs(n, 0.0); for (int k 0; k m; k) { for (int i 0; i n; i) { rhs[i] A[k][i] * w[k] * d[k]; for (int j i; j n; j) { N[i][j] A[k][i] * w[k] * A[k][j]; N[j][i] N[i][j]; } } } X gauss(N, rhs); for (int k 0; k m; k) { V[k] 0.0; for (int i 0; i n; i) V[k] A[k][i] * X[i]; V[k] - d[k]; } double s2 0.0; for (int k 0; k m; k) s2 w[k] * V[k] * V[k]; s2 / (m - n); std::cout std::fixed std::setprecision(5); std::cout X: X[0] X[1] X[2] \n; std::cout sigma0: std::sqrt(s2) \n; return 0; }C 代码的 N 组装是累加写法N[i][j]和N[j][i]同时赋值避免浮点累加顺序造成轻微不对称。工程上点数超过几百时应该改用稀疏矩阵存储但在小网场景下这个稠密版本更直观。若需要协因数阵再对 N 解 n 个单位向量Mat Qxx(n, Vec(n, 0.0)); for (int j 0; j n; j) { Vec e(n, 0.0); e[j] 1.0; Vec col gauss(N, e); for (int i 0; i n; i) Qxx[i][j] col[i]; } // 第 i 个未知点的中误差 sqrt(s2 * Qxx[i][i])3.3 MATLABMATLAB 的矩阵表达最接近教科书适合边算边核对中间量。obs [ 0 1 1.256 1.0; 1 2 -0.582 0.8; 2 3 1.764 1.2; 0 2 0.672 1.5; 1 3 1.183 2.0; 0 3 2.436 2.5; 3 0 -2.441 1.8 ]; known 0; H0 100.0; m size(obs, 1); u 3; A zeros(m, u); d zeros(m, 1); P zeros(m, m); for k 1:m f obs(k, 1); t obs(k, 2); h obs(k, 3); if f known A(k, t) 1; d(k) H0 h; elseif t known A(k, f) -1; d(k) h - H0; else A(k, f) -1; A(k, t) 1; d(k) h; end P(k, k) 1 / obs(k, 4); end N A * P * A; rhs A * P * d; X N \ rhs; V A * X - d; sigma0_2 V * P * V / (m - u); Qxx inv(N); fprintf(X: %.5f %.5f %.5f\n, X); fprintf(sigma0: %.5f\n, sqrt(sigma0_2)); disp(sqrt(diag(sigma0_2 * Qxx)));MATLAB 里数组下标从 1 开始但观测表中的节点编号 0 只是数据值。进入f known分支时t 一定是未知点A(k, t)不会越界进入elseif t known分支时f 一定是未知点A(k, f)同样安全。更稳妥的做法是把未知点从 1 开始编号已知点用单独变量而不是编进索引表。语言依赖求解入口数据结构PythonNumPynp.linalg.solvendarrayC标准库自写高斯消元vectorMATLAB基础环境N \ rhsmatrix三份代码的输入顺序都是from, to, dh, dist权重都取1 / dist_km。现场如果按测站数定权只需要把三个文件里的1.0 / dist全部替换掉。4. 三份水准网间接平差结果互相校验同一组高程和方差同一个网三种实现的差异只在语法不应该出现在结果里。校验时不要只看 X还要看 V、σ₀ 和协因数阵。4.1 应该看到的输出按上面的输入跑完三份代码应得到同一组数程序X1 (m)X2 (m)X3 (m)σ₀ (m)Python101.25562100.67343102.438310.00145C101.25562100.67343102.438310.00145MATLAB101.25562100.67343102.438310.00145残差向量 V 的符号是“平差值减观测值”单位是 m换成 mm 后如下观测序号残差 (mm)1-0.382-0.1830.8741.435-0.3162.3172.69第 6 和第 7 条边都连接 A 和 P3但两条路线的矛盾一正一负P3 的最终高程被两组观测按权拉到了 102.43831 m。4.2 统计量如何手工核对自由度是7 - 3 4加权残差平方和约等于8.40e-6于是σ₀² 8.40e-6 / 4 2.10e-6σ₀ ≈ 0.00145 m。协因数阵由Qxx N⁻¹得到Qxx ≈ [[0.555, 0.325, 0.239], [0.325, 0.599, 0.289], [0.239, 0.289, 0.594]]最终协方差阵是Cx σ₀² * Qxx对角线开方得到三个未知点的中误差约为 1.08 mm、1.12 mm、1.12 mm。只用 X 做平差报告是不够的没有中误差就无法回答“这个点高程差了多少毫米”这种问题。4.3 数值一致性判据手动对比时可以用一个小检查脚本expected np.array([101.25562, 100.67343, 102.43831]) np.testing.assert_allclose(X, expected, atol1e-6)如果三个程序的结果差得超过 1e-6优先检查两处第一处是常数d的三种分支第二处是权重公式。X 差在毫米级通常是h - known[t]和known[t] - h写反X 差在厘米级以上通常是高差符号方向不统一或漏了某条观测边。5. 把水准网间接平差最小实现封装成可复用的平差函数三份代码直接演示有余直接进业务还差一层封装。输入不应该是硬编码观测表而应该是一个(起点, 终点, 高差, 权重)的列表输出统一成(X, V, sigma0, Cx)四个量。调用方只关心前两个精度部分留给报告生成。后验残差一定要做粗差检核。多余观测只有 4 个时不能指望 3σ 法则特别灵敏但它仍然能暴露最明显的错误dof len(obs) - u sigma0 np.sqrt(V P V / dof) normalized np.abs(V) / sigma0 for k in np.where(normalized 3)[0]: print(check obs, k 1, obs[k], normalized[k])最后补三个工程上常用的处理点第一外业反测记录的h必须在导入处统一成路线方向起点和终点不要随意交换第二多个已知点会让网形处于过约束状态如果 N 的条件数很大先查已知点高程是否存在系统差而不是直接改权第三观测边超过一千条时C 换用 Eigen 的稀疏 CholeskyPython 换用scipy.sparse.linalgMATLAB 把A转成sparse再解算原理不变。按同一张观测表跑完三个程序如果 X、V、σ₀ 能对上小数点后 5 位剩下的工作就是把观测数据从代码常量换成数据文件。本文还有配套的精品资源点击获取