资讯详情

资讯详情

误差椭圆详解:从协方差阵到点位精度分析

1. 为什么笔记十二要单独写误差椭圆误差理论与测量平差基础这门课大家最熟悉的肯定是协方差传播、权、条件平差、间接平差这些大块头。等这些基础过了之后随之而来的一个非常实际的问题就是平差算出的坐标点到底有多可靠我们说“点位精度”的时候其实不是在说某个简单数字而是一组方向上的精度差异。笔记写到第十二篇我决定把误差椭圆单独拉出来写。原因很简单这是从“会算平差”到“会看平差结果”的分水岭也是后续做控制网优化设计和施工测量放样时绕不开的精度分析工具。这个话题适合测绘工程、地理信息、土木工程相关专业的学生也适合刚入行的测量从业者。如果你还在死记公式不理解为什么平差之后要算误差椭圆或总是把误差椭圆和误差曲线搞混那这篇笔记就是给你整理的。我会尽量把推导、计算、绘图和工程应用串在一起把教材里分散的内容收拢成一个可以“抄作业”的流程也会把我自己踩过的坑一并写出来。1.1 点位误差不是一根直线能描述的很多同学学完协方差传播后习惯性认为点位的精度就是一个方差或中误差比如“这个点中误差正负3毫米”。严格说这只是在平面坐标两个方向上误差相同且不相关时的简化说法。实际情况中一个点的X和Y坐标误差往往不相等而且两者之间存在相关性比如边长观测误差主要影响一个方向角度观测误差会影响另一个方向。这个时候用一个统一的“圆半径”来描述点位精度会丢掉很多信息。换句话说点位误差在不同方向上是变化的。如果把这个变化规律画出来会得到一个以真值为中心、概率密度等值线为椭圆形状的图这就是误差椭圆。误差椭圆最直接的价值就是告诉你这个点在哪个方向上最不可靠哪个方向上最可靠以及可靠和不可靠的方向具体偏向哪里。1.2 误差椭圆在课程体系中的位置从课程顺序上看误差椭圆一般放在各种平差方法之后、平差程序设计或工程案例之前。它本质上是平差成果精度评定的一部分。就像你算完了一个工程控制网得到每个点的平差坐标如果只把坐标表格交出去别人其实是没法判断这个网到底可靠的。必须要配合点位中误差、误差椭圆、边长相对中误差等指标一起交付才是一份完整的成果。误差椭圆解决的是“点位在哪几个方向上精度如何”的问题和边角精度指标互相补充。而且误差椭圆的概念也不是孤立的。它的底座是协方差阵而协方差阵又是从平差的法方程逆矩阵里取出来的。所以如果前面的条件平差或间接平差没打通误差椭圆算起来就会很吃力。反过来一旦你理解了误差椭圆很多原来抽象的矩阵知识会突然变得具体起来。2. 误差椭圆的参数从哪来协方差阵与特征值计算既然误差椭圆要描述方向相关的精度我们需要从平差结果里拿到点位坐标的协方差阵。这个协方差阵在间接平差里最容易获得因为未知数就是坐标平差值。2.1 从平差结果到协方差阵先回顾一下间接平差的整体流程因为协方差阵不是凭空冒出来的。以间接平差为例误差方程是$$V Bx - l$$法方程为$$Nx W$$其中 $N B^T P B$$W B^T P l$。未知数 $x$ 的协方差阵为$$D_x \sigma_0^2 N^{-1}$$这里 $\sigma_0^2$ 是单位权方差一般在平差后用 $V^T P V / r$ 来估计$r$ 是多余观测数。$N^{-1}$ 是法方程系数阵的逆阵。如果未知数排列为 $x [X_1, Y_1, X_2, Y_2, \ldots]^T$那么每个点对应的 2×2 协方差子阵就是$$ D_i \begin{bmatrix} \sigma_{X}^2 \sigma_{XY} \ \sigma_{YX} \sigma_{Y}^2 \end{bmatrix} $$这才是描述该点误差分布的东西。注意$\sigma_{XY}$ 就是协方差不要因为数值小就忽略它。即使协方差不大它也会影响椭圆的长轴方向尤其在边角网中忽略协方差会导致误差椭圆方向完全扭曲。2.2 特征值与特征向量椭圆长短半轴和方位角的来源误差椭圆的长半轴和短半轴其实是这个 2×2 协方差阵的两个特征值的平方根只是需要乘上单位权方差之后开方。具体说取 $D_i$ 的特征值为 $\lambda_1$ 和 $\lambda_2$$\lambda_1 \ge \lambda_2$则长半轴为$$a \sigma_0 \sqrt{\lambda_1}$$短半轴为$$b \sigma_0 \sqrt{\lambda_2}$$特征向量方向就是长轴和短轴的方向。如果把坐标系的X轴指向北方向、Y轴指向东方向长半轴与X轴的夹角就是误差椭圆的长轴方位角。为什么非得用特征值因为协方差阵本质上描述了误差分布在坐标系中的一个“拉伸加旋转”状态。对角线元素是沿坐标轴方向的方差非对角线元素体现的是旋转。而特征值分解正好把一个任意方向的椭圆变换到它自己的主轴坐标系里。可以用一个生活类比来理解把一块圆形橡皮泥在某个方向上拉伸一下它的主轴就是拉伸方向和垂直于拉伸的方向拉伸程度就是半轴长度拉伸方向就是特征向量方向。2.3 参数计算示例一个模拟点位的完整手算过程我们不用真实项目直接模拟一个点位的协方差子阵方便你对着算。假设某控制点经间接平差后坐标协方差子阵如下单位mm²$$ D \begin{bmatrix} 4.0 1.2 \ 1.2 2.0 \end{bmatrix} $$从矩阵可以看出X方向的方差是4Y方向方差是2两者协方差是1.2说明误差在两个方向存在正相关。求解特征值用特征方程$$|D - \lambda I| (4-\lambda)(2-\lambda) - 1.2^2 0$$展开$$\lambda^2 - 6\lambda 8 - 1.44 0$$ $$\lambda^2 - 6\lambda 6.56 0$$解出$$\lambda_1 \approx 4.82\lambda_2 \approx 1.18$$对应特征向量把 $\lambda_1$ 代回 $(D-\lambda I)v0$得到 $v_1$ 方向约为 $[0.85, 0.53]$与X轴的夹角约为32°。假设单位权中误差 $\sigma_0$ 已经算出为1.5 mm那么长半轴 $a 1.5 \times \sqrt{4.82} \approx 3.29$ mm短半轴 $b 1.5 \times \sqrt{1.18} \approx 1.63$ mm看起来X方向方差比较大但长轴并不完全指向X轴而是偏向32°说明因为协方差的存在最不可靠的方向并不是单纯的X方向。这个例子虽然简单却说明一个问题只要忽略了协方差你就永远找不准误差椭圆的长轴方向。3. 误差椭圆怎么画工具、步骤、精度检查参数算出来是一回事真正画到图上看才是工程习惯。误差椭圆通常叠加在控制网图上用来快速判断哪些点精度差以及误差方向是否合理。3.1 用Excel一步步画误差椭圆Excel不是专业的绘图工具但胜在随处可用适合学习阶段检查和验收。做法是这样的先准备一列角度 $t$从0°到360°步长建议取2°。然后根据参数方程计算椭圆边线在长轴坐标系里的坐标$$u a \cos t$$$$v b \sin t$$再把长轴坐标系转换到原坐标系$$X X_0 u \cos\theta - v \sin\theta$$$$Y Y_0 u \sin\theta v \cos\theta$$其中 $X_0$、$Y_0$ 是待画点的平差坐标$\theta$ 是长轴方位角。在Excel里逐行填公式最后用XY散点图画出坐标对就可以了。这里有个关键点Excel散点图默认的X轴和Y轴比例可能不一样这会导致椭圆被“压扁”或“拉长”干扰你对形状的判断。最简单的解决办法是手动把横纵坐标显示范围设成一致或者直接看数值而不是只看形状。如果你只是验证椭圆方向和大小保持纵横尺度一致还是很重要的。3.2 用Python脚本批量绘制与检查实际项目里点位多动辄几十上百个点一个个在Excel里画不仅枯燥还容易因为人工填写坐标而抄错。我建议直接用Python脚本批量处理。下面这个脚本把平差软件导出的坐标和协方差阵读取出来按点位绘制误差椭圆。代码很基础但作用直接适合作为学习阶段的脚手架也能直接套用到小规模的工程图上。如果你还没有现成数据可以在脚本里把示例的协方差阵替换成2.3节的算例验证一下输出结果是否和手算一致。import numpy as np import matplotlib.pyplot as plt def plot_error_ellipse(ax, x0, y0, cov, scale1.0, **kwargs): # 特征值分解 eigvals, eigvecs np.linalg.eigh(cov) # 协方差阵已经是平差后的成果协方差阵直接开方得到半轴 a np.sqrt(eigvals[1]) * scale b np.sqrt(eigvals[0]) * scale theta np.arctan2(eigvecs[1, 1], eigvecs[0, 1]) # 角度序列 t np.linspace(0, 2 * np.pi, 100) u a * np.cos(t) v b * np.sin(t) ell_x x0 u * np.cos(theta) - v * np.sin(theta) ell_y y0 u * np.sin(theta) v * np.cos(theta) ax.plot(ell_x, ell_y, **kwargs) ax.plot(x0, y0, marker, linestylenone) # 示例三个点位 coves [ (100.0, 200.0, np.array([[4.0, 1.2], [1.2, 2.0]])), (300.0, 180.0, np.array([[9.0, -0.8], [-0.8, 1.0]])), (500.0, 220.0, np.array([[1.0, 0.5], [0.5, 2.5]])), ] fig, ax plt.subplots(figsize(6, 6)) for x0, y0, cov in coves: plot_error_ellipse(ax, x0, y0, cov, scale1.0, labelerror ellipse) ax.set_aspect(equal) ax.set_xlabel(X (mm)) ax.set_ylabel(Y (mm)) plt.show()这段代码里有个细节需要说清楚特征值本质上是方差在主方向上的度量开平方才是对应方向的标准差。如果你的协方差阵已经是“成果协方差阵”也就是已经乘过单位权方差那么 $\sqrt{\lambda}$ 直接就是半轴长度如果协方差阵还只是先验单位权下的协因数阵那还得再多乘一个 $\sigma_0$。很多同学在这里把单位权方差重复乘了一次导致椭圆大得离谱。我的习惯是在代码注释里明确写清这一步。3.3 绘图前必须做的单位与坐标系核对画误差椭圆前最常踩的坑有三个。第一是单位不一致。平差软件里协方差阵的单位可能是毫米平方而坐标常以米为单位。画图时如果坐标用米、误差椭圆用毫米椭圆就会缩成一个点。解决方法很简单把坐标统一换算成毫米或者把椭圆半轴除以1000后再画。我见过不少同学卡了半天最后发现只是单位没统一这种低级错误最可惜。第二是坐标轴方向。测量坐标系中通常X方向是北方向Y方向是东方向与数学坐标系中的X、Y习惯不同。绘图时如果不保持轴比例一致椭圆就会被视觉拉伸。尤其是用Excel或某些默认绘图库时横纵坐标比例不一致很容易让人误判长轴方向。建议设置ax.set_aspect(equal)这是Python绘图里很关键的一步。第三是特征值排序不统一。有些库函数默认返回升序有些是降序。如果你不管顺序就直接取sqrt(eigvals[0])很可能把短半轴当成画长半轴导致椭圆整体旋转90度。稳妥的做法永远是先比较两个特征值的大小再确定哪个对应长轴而不是假设索引0就是长轴。4. 误差椭圆在工程中的应用思路误差椭圆不是课堂上的概念玩具在控制测量设计和结果分析中非常实用。4.1 控制网精度评估看椭圆方向对不对拿到一个GPS网或边角网的平差结果后把每个点的误差椭圆画到图上马上就能发现一些规律通常离已知点较远的点误差椭圆会大一些在网的边缘椭圆长轴往往指向远离已知点的一侧。如果某个点的误差椭圆方向出现明显反常比如短边密集区域的点反而长轴很长那很可能是观测值之间有粗差或设计网形欠佳。这时候不要急着调网先对误差椭圆方向做一个整体判断再决定是否需要补测。我在实际项目里见过一种情况某个区域所有点的误差椭圆长轴都指向同一个方向刚开始以为是北方向系统误差后来检查才发现是已知点约束时某个数据存在尺度异常。把异常点剔除后椭圆方向的整体偏向就消失了。误差椭圆在这时起到了“全局体检”的作用。4.2 隧道贯通测量中的误差椭圆隧道贯通中的横向贯通误差是施工测量最关心的问题之一。隧道两端分别从两个洞口向中间掘进两端控制点误差都会传递到贯通面。把两端控制点和导线网进行联合平差后贯通面待定点的误差椭圆长轴方向如果偏向横向说明横向误差风险高如果偏向纵向则横向误差较小。这里的处理手法是在贯通面中心设置一个模拟点把从两端推算到该点的误差椭圆分别算出来再看长轴方向与隧道轴线的关系。如果长轴明显垂直于隧道轴线那就必须加强横向方向的观测比如增设短边或采用更高精度的测角仪器。这个方法在选线阶段就能为洞外控制网设计提供量化依据比单纯看点位中误差要直观得多。4.3 用误差椭圆指导放样与点位布设在施工放样中放样点位通常由控制点通过极坐标法或前方交会得到。这时放样点位的误差椭圆是控制点误差、放样测角误差、测边误差三者的综合结果。我们可以提前计算不同放样方案的误差椭圆选择放样误差在工程允许范围内的方案。比如桥梁桥墩中心线放样假设桥墩容许横向偏差是5毫米我们就需要对比从两岸控制点直接放样与通过后方交会放样两种方案的误差椭圆。如果直接放样方案的椭圆长轴刚好落在横向那就该换方案。这种分析在平差软件里已经做得比较成熟但核心思想仍然是看椭圆方向而不是只看椭圆大小。5. 误差椭圆相关的高频坑点与排查方法这一节是我自己一路踩坑总结出来的按出现频率从高到低排希望能帮你节省不少排查时间。5.1 协方差阵非正定怎么办理论上协方差阵必须是正定的但实际平差中偶尔会算出半正定甚至负定的子阵。原因可能是网形存在秩亏或者约束数量不足也可能是浮点计算精度问题。排查方法是先看整个法方程矩阵的秩。如果秩亏说明无约束或约束不足需要先加最小约束再平差如果只是某个点的子阵异常看看这个点是不是属于多余观测数太少的“悬空点”比如由两条短边交会形成的细长三角。解决思路是增加观测值或重测而不是强行修正协方差阵。另一个容易忽略的原因是在提取子阵时坐标顺序搞反导致协方差阵不对称出现虚特征值。检查一下提取的行列号对应点位是否一致。5.2 长半轴方向角算错的原因长半轴方向角计算是手算最容易错的地方常见原因有两个。第一个是把特征向量的两个分量搞反或者少写了负号导致角度符号错误。处理办法是画出坐标轴把特征向量画在图上验证如果特征向量指向第一象限但你的角度却是负的那肯定是符号问题。第二个是没有把角度换算到测量方位角。数学上的arctan返回的是-90°到90°而测量方位角是0°到360°需要根据X、Y分量的符号判断实际象限。比如特征向量X分量为正、Y分量为负时真实方位角应该在第四象限可以按360°减去绝对角度来处理。建议写一个自动转换函数而不是手动判断。5.3 误差椭圆与误差曲线的区别我一直觉得教材里最容易被忽略的坑点是误差曲线。误差曲线不是误差椭圆它是以点位到曲线的向量长度表示对应方向中误差的图形形状通常不是椭圆而误差椭圆则是固定概率下的等值线。两者经常被搞混的原因是图示上都有“椭圆”的趋势但误差曲线的径向长度直接代表那个方向的标准差误差椭圆的径向长度则需要乘以一个置信系数才能反映对应概率。在成果报告里如果你画的是误差曲线就必须标注方向向量对应的是“σ”如果画的是误差椭圆一般要标注置信概率比如95%或一倍中误差。把两者混着标注是很不专业的做法。5.4 相对误差椭圆与点位误差分解在桥梁、滑坡监测等场合真正关心的是两个点之间的相对精度而不是单点绝对精度。两个点的相对误差椭圆可以由两点的协方差子阵叠加得到$$D_{rel} D_1 D_2 - 2D_{12}$$如果平差软件没有输出两点间的互协方差 $D_{12}$那就不能用简单相加。有些初学者直接把两个点的协方差阵相加忽略了负的互协方差项结果相对椭圆偏大很多。正确做法是检查平差软件是否提供了完整协方差阵如果只提供点位子阵就需要从法方程逆阵中重新提取包含两个点的完整子块。这一步是工程成果精度分析中最容易出问题的地方一定不要偷懒。6. 笔记十二之外从误差椭圆到现代平差笔记写到这里误差椭圆似乎已经讲完了但我觉得还是应该把视野稍微拉开一点。误差椭圆的思想在后面的课程和实践里还会延伸出很多内容。6.1 稳健估计与粗差定位中的椭圆意义在做粗差检测时我们常看到“数据探测法”中对残差做标准化检验。标准化残差其实和误差椭圆的思路同源因为残差向量也带有协方差阵。如果某个观测值的残差相对其理论误差明显异常就可以通过马氏距离来判断。马氏距离的定义就是点到分布中心的距离在协方差椭球尺度下的归一化数值。理解了误差椭圆你对马氏距离就不陌生了。在GNSS基线解算中三维点位协方差阵对应的是误差椭球。二维误差椭圆是三维椭球在平面上的投影。很多GNSS后处理软件直接给出平面坐标和协方差子阵其实就是在平面方向做了投影后的结果。如果三维椭球的长轴方向基本垂直于地面那么平面误差椭圆会很小这时我们就要警惕高程方向的精度是否合理。这也是为什么看GNSS结果时不能只盯着平面椭圆的原因。6.2 误差椭球与三维平差的联系三维平差中点的协方差阵是3×3特征值分解后得到三个半轴长度构成误差椭球。实际工程中往往取其投影到二维平面后的椭圆来使用。投影的方法是把三维协方差阵投影到局部切平面上而不是简单丢掉高程分量。如果高程分量和平面分量相关简单去掉一维会丢信息。用来进行高精度施工控制网设计时最好保持三维计算再投影到施工坐标系得到的误差椭圆才是真实的二维精度。6.3 可以继续做的延伸练习如果你看到这里我建议你找一套现成的边角网数据手动做一次全流程间接平差解算坐标提取协方差阵计算每个控制点的误差椭圆参数再用Python画在网图上。然后试着改变其中一个观测值的权观察误差椭圆会发生什么变化。这个练习会帮你真正理解权与精度的关系比背十遍公式都管用。我自己是在一次模拟项目X的网形优化中通过反复调整观测权和误差椭圆参数才彻底明白为什么增加一个对角线观测会大幅压缩某个方向的长半轴。这种体会很难直接从书本上获得但如果你动手做一遍就会发现误差椭圆其实就是为这种“方向敏感性”而生的。这篇笔记写到这里我突然想起第一次画误差椭圆时卡了两天也没搞明白的旧事。当时我把协方差阵提取错了导致所有椭圆都旋转了45°后来逐项核对数据才找到问题。所以现在每次出图前我都会先挑一个精度已知的简单控制网做验证确认椭圆的尺度、方向和理论预期一致再开始处理正式数据。这个习惯替我挡掉了很多返工。如果你刚接触误差椭圆建议也能从一个小网练起把坐标、协方差阵、画图脚本全部打通再逐步扩展到复杂网形。
觉得有用,分享给同行:

为您的企业打造数字门面

稳重轻奢商务风格,端正雅致视觉,长效耐看不易过时。

立即咨询 →