测量平差入门:水准网闭合差处理与最小二乘精度评定 咱们搞测量的最怕的不是仪器不好用而是数据测回来之后自己都说不清楚问题出在哪。你辛辛苦苦测完一段水准闭合差超了限监理问你怎么处理你说“重测”理论上没问题但成本谁出如果每次都要靠重测解决问题这活就没法干了。误差理论与测量平差基础这门课解决的就是测完之后“怎么把矛盾合理消化掉并且给出让各方都信服的精度结论”。这篇笔记八正好是整门课从概念走向计算的转折点前面几篇把偶然误差特性、方差传播律、权与定权讲完之后从这里开始真正进入解算环节。如果你正在学这门课或者工作中需要经常处理控制网数据这篇内容可以帮你把“条件平差”“间接平差”“最小二乘”“精度评定”这几块东西串起来。我尽量不堆公式而是用一套完整的水准网数据带你从头算到尾中间该解释的“为什么”我都会解释。1. 笔记八的位置这两节课为什么让很多人开始怀疑人生1.1 误差理论到底在解决什么问题先说一个容易被忽略的事实测量学教你怎么测误差理论教你怎么面对“测不准”这件事。任何观测都会有误差这不是仪器不够好而是客观上必然存在。误差按性质可以分成系统误差、偶然误差、粗差三类其中偶然误差是平差理论主要处理的对象。偶然误差满足四个特性有界性、单峰性、对称性、抵偿性。这四个特性你可以这样理解误差不会无限大小误差出现概率高正负误差机会均等观测次数足够多时误差之和趋近于零。这四个特性是所有平差方法的统计学基础。因为误差存在测出来的数据就会出现相互矛盾的情况。比如一条附合水准路线你测了五段高差按道理从A点传到B点累计高差必须等于B点高程减A点高程但实际测量结果往往不等于这个理论值多出来的那部分就是闭合差。闭合差怎么处理就是平差的问题。好多人学到这里开始犯晕是因为把误差理论当成了一门纯数学课天天背公式却不理解这些公式在解决什么现实问题。实际上这门课的核心逻辑非常朴素先承认观测有误差再用统计学方法估计出最可靠的结果同时给出这个结果有多可靠。1.2 为什么平差不是“凑数”而是“加权分配矛盾”平差这个词听起来像“把数据抹匀”实际上完全不是这个意思。平差的本质是在满足几何条件的前提下找到一组改正数加到观测值上让改正后的观测值满足所有理论约束并且让这组改正数在某种原则下最优。这个“最优原则”通常是最小二乘法也就是让改正数的加权平方和最小。为什么要加权因为不同观测值的精度不一样同样一段路线距离短的水准测量误差小距离长的误差大如果你把闭合差平均分配到每一段上相当于把短距离的高精度观测和长距离的低精度观测一视同仁这显然不合理。加权之后精度高的观测分到的改正数小精度低的观测分到的改正数大这才符合误差分布规律。从数学上理解平差问题是一个带约束的优化问题。观测值有n个必要观测数是t个多余观测数是r n - t个。r大于零时观测值之间会产生多种组合理论上都能得到唯一解但不同组合的结果互不一致这就是矛盾来源。平差就是在这个矛盾空间里按加权最小二乘原则找一个最优点。笔记八通常讲到这里就是要开始把概念模型变成可计算的矩阵模型了。接下来我拿一个具体的水准网把整个流程走一遍。2. 核心模型实例一个水准网如何被拆成数学问题2.1 案例设计与观测数据我设计一个最简单的附合水准网。已知点A的高程是10.000米已知点B的高程是12.000米中间有三个待定点P1、P2、P3共观测了五段高差。观测数据如下测段起点终点观测高差(米)路线长(千米)h1AP10.8021.0h2P1P21.5050.8h3P2P3-0.3931.2h4P3B0.0971.5h5P2B-0.3020.6注意看数据里的矛盾。从A经过P1、P2、P3到B累计高差是0.802加1.505减0.393加0.097等于2.011米。但B点高程减A点高程是2.000米闭合差11毫米。另一条路A到P1到P2再到B累计高差是0.802加1.505减0.302等于2.005米闭合差5毫米。两条路线推出来的P2高程也不完全一样。这些矛盾就是我们要通过平差消除的。这是一个典型的多余观测案例。必要观测数t等于3观测值数n等于5多余观测数r等于2。r大于0说明这个网具备平差条件也有能力对观测质量做检验。2.2 误差方程怎么列间接平差的第一步是选参数。这里三个待定点的高程就是最自然的参数记为x1、x2、x3。观测高差可以表示成这些参数的线性函数。逐条写出来第一条高差是从A到P1理论上h1应该等于P1高程减A点高程也就是x1减10.000。观测值加改正数等于理论值误差方程写成v1 x1 - 10.000 - 0.802 x1 - 10.802第二条高差是从P1到P2理论上等于x2减x1所以v2 x2 - x1 - 1.505第三条高差是从P2到P3理论上等于x3减x2v3 x3 - x2 - (-0.393) x3 - x2 0.393第四段是从P3到B理论上等于B点高程减x3也就是12.000减x3v4 12.000 - x3 - 0.097 11.903 - x3这里我习惯把v4写成x3的函数形式整理成统一格式时它会变成v4 -x3 11.903后面矩阵化的时候要注意符号。第五段是从P2到B理论上等于12.000减x2v5 12.000 - x2 - (-0.302) 12.302 - x2写成矩阵形式就是 V Bx - l。这里的x是参数平差值l是常数项。2.3 最小二乘解算的完整流程间接平差的准则只有一个V^T P V min。其中P是权阵因为各段观测相互独立P是主对角阵对角元素是每段观测的权。水准测量中权通常取路线长度的倒数也就是P_i 1 / s_i。路线越长观测误差越大权越小。这里用C1路线长度除以千米五段权的值分别是1、1.25、0.8333、0.6667、1.6667。推导过程在很多教材里有我这里只说结论。令V Bx - l代入V^T P V对x求导并令其等于零得到法方程B^T P B x B^T P l令N B^T P BW B^T P l解出来就是x N^{-1} W这个x就是三个待定点高程的最小二乘估计值。整个过程看起来很简单但每一步矩阵运算都对应一个实际意义。N矩阵反映的是整个网的几何结构强度W向量反映的是观测数据与近似值的差异。参数估计完成后再把x代回V Bx - l就能得到每段观测的改正数。改正数的大小可以直接反映观测质量。3. 手算与Python复现把课堂公式变成能跑的结果3.1 中间矩阵的计算过程按照上面的误差方程写出B矩阵和l向量。B矩阵是5行3列行对应观测值列对应参数B [1, 0, 0] [-1, 1, 0] [0, -1, 1] [0, 0, -1] [0, -1, 0]注意第四行我之前把v4写成-x3 11.903所以B矩阵第四行对应x3的系数是-1。l向量是l [10.802] [1.505] [-0.393] [11.903] [12.302]权阵P是对角阵对角元是[1.0, 1.25, 0.8333, 0.6667, 1.6667]接下来计算N B^T P B。这一步手算容易出错但结果非常有规律。算出来N [2.25, -1.25, 0] [-1.25, 3.75, -0.8333] [0, -0.8333, 1.5]W B^T P l算出来W [8.92075] [22.712] [7.6078]解方程组得到x [10.7986] [12.302] [11.906]把x代回误差方程得到各段改正数v [-0.0034] [-0.0016] [-0.003] [0.003] [0]这个结果很合理。五段残差都在几毫米量级说明观测数据质量可以平差后的高程值也符合已知点约束。3.2 Python代码实现手算一遍是为了理解原理实际项目里肯定用程序。这里用Python的numpy库写一个完整实现代码可以直接复制运行。import numpy as np # 已知点高程 HA 10.000 HB 12.000 # 观测高差和路线长度 h np.array([0.802, 1.505, -0.393, 0.097, -0.302]) s np.array([1.0, 0.8, 1.2, 1.5, 0.6]) # 权与路线长度成反比 P np.diag(1.0 / s) # 误差方程 V Bx - l B np.array([ [1, 0, 0], [-1, 1, 0], [0, -1, 1], [0, 0, -1], [0, -1, 0] ]) l np.array([ HA h[0], h[1], h[2], HB - h[3], HB - h[4] ]) # 法方程 N B.T P B W B.T P l # 解算参数 x np.linalg.solve(N, W) print(平差高程, x) # 改正数 V B x - l print(改正数, V) # 单位权中误差 r len(h) - len(x) sigma0 np.sqrt(V.T P V / r) print(单位权中误差, sigma0) # 协因数阵 Qxx np.linalg.inv(N) print(协因数阵, Qxx) # 高程中误差 sigma_x sigma0 * np.sqrt(np.diag(Qxx)) print(高程中误差, sigma_x)运行结果和我手算的一致。这里要提醒一句numpy的矩阵乘法符号是Python 3.5之后才有的如果你的环境比较老改成np.dot也是可以的。3.3 判断结果是否合理算完不是结束还要判断结果到底能不能用。我最常看两个指标。第一个是改正数V如果某一段的改正数明显比其他段大很多比如五段残差都是毫米级突然有一段是厘米级那这一段很可能存在粗差需要重点检查。第二个是单位权中误差σ0它是平差后单位权观测值的中误差。如果σ0和仪器标称精度在同一量级说明观测质量正常如果σ0明显偏大说明整网观测质量有问题或者定权不合理或者存在未发现的系统误差。在实际项目中我还会把平差结果带回原始观测条件验证。比如算出来的P2高程是12.302米那h2平差后的值是12.302减10.7986等于1.5034米与原观测值1.505差1.6毫米这个残差在容许范围内。这种“还原验证”看着简单但能防止很多低级别错误。4. 精度评定平差值什么时候能交付使用4.1 单位权中误差平差算出的坐标或高程只是第一步甲方要的不是一个数字而是这个数字有多可靠。精度评定就是回答这个问题。单位权中误差σ0是整网精度的一个综合估计公式是σ0 sqrt(V^T P V / r)这里的r是多余观测数也就是自由度。自由度越大精度估计越可靠。如果r等于0观测值个数正好等于必要观测数这时V等于0σ0算不出来平差也就失去了意义。这从另一个角度解释了为什么多余观测是必要的。计算V^T P V时要注意P是权阵不是单位阵。如果所有权都乘以同一个常数V不变但V^T P V同比例变化σ0也会变。这看起来是个问题实际上不是。因为权的绝对值不影响参数估计结果只影响精度估计的基准。在同一个工程项目里只要定权方式一致σ0的对比就有意义。4.2 未知数协因数阵与中误差参数x的精度不是均匀的每个点的精度都不一样。精度信息藏在协因数阵里Qxx N^{-1}协因数阵对角线的平方根乘以单位权中误差就是对应参数的中误差σ_xi σ0 * sqrt(Q_ii)拿刚才的例子Qxx算出来大致是[[0.5635, 0.2143, 0.1190], [0.2143, 0.3857, 0.2143], [0.1190, 0.2143, 0.7857]]σ0是0.0038米所以三个点的高程中误差分别是0.0029米、0.0024米、0.0034米。P3的位置离已知点A最远精度最低这个规律符合直觉。平差理论的价值就在这里你不用等测完再做实验光从网形和观测精度就能提前预判哪些点位精度高哪些点位精度低这正好是网形设计阶段需要的功能。4.3 平面网中的点位误差与误差椭圆上面讲的是水准网是一维问题精度用一个数就能描述。但平面控制网是两个方向X方向和Y方向的误差合起来会形成一个二维分布用误差椭圆表示更完整。误差椭圆的两个半轴由协方差矩阵的两个特征值决定长半轴方向对应误差最大的方向短半轴方向对应误差最小的方向。刚接触这个概念的测量员容易犯一个错误只关注点位中误差M_p sqrt(σx^2 σy^2)完全不看误差椭圆的方向。实际工作中如果甲方需要控制横向误差比如桥梁墩台定位误差椭圆长轴方向沿桥轴线还是垂直桥轴线结果完全不一样。点位中误差把方向信息抹掉了误差椭圆保留了方向信息。做高精度工程控制网时误差椭圆必须画出来看不能只看合成中误差。5. 常见问题与排查技巧实录5.1 误差方程符号经常搞反怎么办我做学生的时候符号问题是犯错最多的地方。V Bx - l使用频率最高但这个式子里的l既不是观测值本身也不是观测值取负而是“用近似参数计算得到的观测值近似值”。不同教材对l的定义略有差异有些教材写成V Bx l有些写成V l - Bx。关键是选定一种后必须坚持到底不能混合使用。我的习惯是先写出物理含义明确的观测方程比如“h2的理论值 x2 - x1”然后写成“v2 x2 - x1 - h2”。这样可以避免只记公式不记含义带来的混乱。遇到符号不确定时就回到物理定义去推最多几分钟比死记硬背可靠得多。5.2 法方程病态如何识别法方程病态在测量控制网里不算罕见尤其是点位分布差、观测结构差的网。病态的表现是N矩阵的行列式接近零解出来的参数对观测值微小变化极其敏感。一个很实用的检查方法是看N的条件数条件数可以用numpy的np.linalg.cond计算。条件数超过1e6基本可以判定病态。处理病态的方法有几个一是重新选择参数尽量选独立性强、相关性低的量作为未知数二是增加观测加强网形结构三是改变参数化方式比如用基线向量代替绝对坐标。这里最忌讳的是不管病态不病态直接硬建模算出来的数字可能很漂亮但实际一点参考价值没有。5.3 定权标准不一致导致精度失真定权是平差里最容易被忽视但又最影响结果的一步。水准测量通常用路线长度定权导线测量通常用测站数或边长定权GNSS网通常用基线解算精度定权。如果你一会儿用距离定权一会儿用测站数定权算出来的σ0就失真了。项目上我常做的检查是平差完成后把σ0和仪器标称精度做对比。比如一台标称每公里偶然中误差为1毫米的水准仪一段1公里路线高程观测中误差应该大约是1毫米如果σ0算出来是3毫米要么是观测确实有问题要么是权定得不对。这个对比能帮你快速判断是否需要重新定权。5.4 残差最大的观测不一定就是要重测的观测残差大说明这个观测值和平差后整体结果不一致但这不意味着它一定是粗差。在闭合差内部一个观测值有粗差会把误差分摊到好几个观测值的残差上有时候最大残差反而出现在没有粗差的观测段上。比较稳妥的办法是看标准化残差也就是把残差除以对应观测值的中误差再做判断。标准化残差超过2到3倍时才需要重点怀疑该观测含有粗差。另外粗差检测最好在平差完成后系统做不要看一个残差大就马上重测先把数据分布、网形结构、定权方式全部检查一遍。我整理了一个常见问题速查表平时遇到问题可以直接对照现象可能原因排查方向法方程解不出来N矩阵奇异检查多余观测数是否大于0参数是否独立改正数全部偏大定权不合理或观测值有粗差重新检查权定义逐项检查观测记录某点精度异常低该点周围观测数量不足或网形差增加该点相关观测调整网形σ0远超仪器标称存在系统误差或粗差未剔除检查仪器、观测条件做粗差检测不同平差软件结果不一致参数选择、定权方式、约束方式不同统一输入文件和定权标准后重算已知点参与平差后位移过大已知点本身有误差或未做稳定性检查先检测已知点之间的兼容性6. 最后再分享几点实际操作中的习惯我做了这么多年测量数据处理有一个体会特别深平差不是把数据交给软件就完事的流程而是需要你不断根据结果做判断的决策过程。初学的时候我建议大家手算一遍完整流程再做一遍编程实现最后再上商业软件。手算能让你理解每个矩阵元的来源编程能让你理解软件底层在做什么之后再用商业软件才不会变成一个只会点按钮的人。还有一个小习惯每次平差前先把观测值单位统一。角度用度分秒还是弧度高差用米还是毫米距离用公里还是米这些单位不统一轻则数值差几个量级重则法方程病态。我见过太多人在这上面栽跟头数据算到一半发现结果离谱回头一查原来是某一路段的距离忘了从米换算成公里。最后一个建议一定要保留平差前的原始观测文件和全部中间计算记录。工程验收或者数据复核时别人不光要你的最终成果还要你整个处理过程的可追溯性。误差理论这门课真正带给你的不是让你会背公式而是让你养成一种习惯每一个数字都说得清来源每一项成果都经得起推敲。