1. 粗差探测到底在解决什么问题搞大地测量的人几乎每天都要跟高程数据打交道。GNSS水准测量、似大地水准面精化、沉降监测、工程控制网复测这些项目最终的成果都落在一堆高程值上。数据从外业采集回来内业处理的第一步往往不是平差而是先问一句这批数据里有没有“坏点”。所谓坏点在测量领域正式叫法是粗差。粗差不是偶然误差偶然误差还有个正态分布的规律可循粗差是观测过程中因为人为失误、仪器故障、环境突变、记录笔误等原因产生的异常值。比如水准测量读数时读错一位GNSS拟合时天线高量错甚至数据传输过程中某个点位被污染。粗差的特点是量级大、偏离明显、分布没有规律而且对最终结果的污染能力极强。一个很典型的现象是在GNSS高程拟合中如果控制点里有1到2个点的高程存在粗差拟合后的模型会被这些点“带偏”导致整个区域的高程转换参数失真。更麻烦的是粗差对最小二乘平差的影响不是按比例分摊的一个大的粗差能把多个正常点“拉下水”这就是为什么必须在外业数据进入正式平差之前专门做一遍粗差探测与剔除。很多人觉得粗差探测不就是把明显偏离的点去掉吗没那么简单。实际处理中粗差的大小、数量、分布位置都是未知的而且粗差和正常观测值的边界往往模糊。尤其是小范围、点位少的工程数据一个粗差点可能会被误判为正常高值或者一个正常的极值点被误杀。这里面的取舍艺术恰恰是内业处理经验的核心。这篇文章我会从原理讲到实操重点讲清楚几种常用的粗差探测方法各自适用什么场景、参数怎么取、代码怎么写以及我这些年处理高程数据时踩过的坑。2. 粗差探测的常用方法及其原理2.1 3σ准则经典但必须谨慎使用3σ准则拉依达准则是几乎所有测量数据处理教材都会提到的方法。它的逻辑很简单假设观测值服从正态分布那么观测值与均值的偏差超过3倍标准差σ的概率只有约0.27%。根据小概率事件原理发生了这种“不太可能”的事件就认为该观测值是粗差。公式表达为|vi - v_mean| 3σ → 判定为粗差其中v_mean是所有观测值的算术平均值σ是标准差按贝塞尔公式计算σ sqrt( Σ(vi - v_mean)² / (n-1) )用起来确实快一跑就出结果。但我必须提醒3σ准则存在两个致命短板。第一个短板是掩盖效应。当粗差值非常大时它会显著拉高σ的计算结果导致阈值3σ变得很大反而让粗差点“藏在”阈值之内。比如一组数据正常值在100米附近混入一个实测值为150米的粗差均值被拉高σ也被拉大这时候150米这个点可能反而满足“偏差小于3σ”的条件被当作正常值放行。这个现象叫粗差对均值与标准差的“污染”。第二个短板是样本量限制。3σ准则严格来说是建立在n足够大的前提下的当样本量小于10时σ本身的计算就不稳定阈值意义大打折扣。工程控制网中高程点往往就几个到十几个直接用3σ很容易误判。所以我的建议是3σ准则可以用于快速预览判断数据是否存在明显的异常量级但不要作为唯一的判定依据。尤其在点位少的高程拟合数据里更不要裸用。2.2 中位数和MAD抗差性更强的替代方案既然3σ准则的核心问题是均值与标准差被粗差污染那自然的思路就是改用抗污染的统计量。中位数就是最典型的稳健统计量。中位数对粗差的抵抗力极强哪怕一组数据里混入几个极端值中位数也几乎不受影响。类似地MADMedian Absolute Deviation中位数绝对偏差是标准差的一种稳健替代MAD median( |vi - median(v)| )在此基础上可以构造稳健的Z分数Z_i 0.6745 * (vi - median(v)) / MAD这里的0.6745是校正系数目的是让MAD在正态分布下等于σ。当Z_i的绝对值大于某个阈值通常取2.5或3时判定该点为粗差。这个方案的优点是抗污染能力强不会因为个别粗差导致整体阈值失真。缺点是对真值的分布形态有一定依赖如果数据本身不是近似对称分布中位数的表现也会打折扣。我在实际项目中的做法是先用中位数MAD做第一遍粗差探查把明显异常的点标记出来再用3σ准则对剔除后的干净数据做复核。两道工序配合比单一方法可靠得多。2.3 Dixon检验法适合小样本数据的粗差判定Dixon检验法Q检验法专门针对小样本数据的异常值检测在测量数据处理中应用非常广泛。它的核心思想是计算极值与邻近值的差距占整个数据极差的比值如果这个比值超过临界值则判定为离群值。对于升序排列的样本 x1 ≤ x2 ≤ ... ≤ xn检测最大值时用Q (xn - x(n-1)) / (xn - x1)检测最小值时用Q (x2 - x1) / (xn - x1)根据样本量n和置信水平查Dixon临界值表如果计算出的Q大于临界值判定该极端值为粗差。Dixon检验在处理小样本n3到30时比3σ准则稳定得多因为它不依赖均值和方差只用排序后的位置关系。这也是为什么它在测绘领域特别受青睐。但要注意Dixon检验一次只能检测一个极值检测完需要剔除后再对剩余数据重新检验不能一次把所有可疑点全部判定。2.4 稳健估计迭代法最接近工程实用前面几种方法本质上是统计检验还有一个思路是直接改造估计过程本身让粗差在估计阶段的权重自动降为0。这个方法在测绘平差领域叫稳健估计Robust Estimation常用的有丹麦法、IGG法周江文方案等。以IGG法为例它的核心思想是先做一次普通最小二乘估计得到残差然后根据残差大小给每个观测值重新定权残差小正常的权重不变残差中等可疑的权重降低残差大粗差的权重设为0再用新权重重新估计如此迭代直到权重变化稳定。我用一个简化的流程来说明1. 初始权重全部设为1做一次最小二乘平差 2. 计算各观测值残差 vi 和单位权中误差 σ0 3. 对每个观测值重新定权 - |vi| ≤ c1·σ0 → 权重不变 - c1·σ0 |vi| ≤ c2·σ0 → 按函数降权 - |vi| c2·σ0 → 权重置零 4. 用新权重重新平差重复2~4步 5. 当迭代前后估值变化小于阈值时收敛这里的c1、c2是阈值系数通常取c12.5、c24.0可以根据粗差比例适当调整。这个方法的优势在于它不只是“找出粗差”而是把粗差处理融入了平差流程本身。尤其是处理控制网数据时高程异常点会被自动降权或剔权最终平差结果更稳健。缺点是计算流程比前面几种方法复杂需要迭代运算而且阈值参数需要经验积累。3. 高程数据粗差探测的系统设计与工具实现3.1 工具的总体架构与处理流程基于上面的方法梳理我设计了一个高程粗差探测工具。这个工具不是只跑一个算法而是串起一条完整的数据处理链路从数据读入到最终输出报告中间每一步都有明确的逻辑和可追溯的记录。整体流程分为五个环节数据读入与规范化处理接收点号、高程值、可能的平面坐标等字段统一格式处理缺失值。粗差初筛与可视化计算基本统计量生成分布直方图和序序图让数据形态直观可见。多方法联合探测依次使用中位数MAD、3σ准则、Dixon检验对可疑点进行标记并输出各方法的判定结果。人工复核与决策结合点位空间分布和项目背景确认最终剔除名单。输出处理结果与报告生成剔点前后对比表、统计指标变化、XY方向分布图。其中第二步的可视化是我特别想强调的。很多人在粗差探测时只看数字结果忽略了图形直觉。高程数据在空间上的分布是有规律的一个点在高程异常序列中可能不明显但在平面位置和高程值的三维关系里往往一目了然。比如某个GNSS水准控制点距离其余点很远偶然高程误差又大这个点在统计检验中可能不达标但实际可能是点位本身代表性不足而不是观测粗差两者的处理方式完全不同。3.2 核心代码实现下面给出我用Python实现的核心算法代码包括多方法联合探测模块。注意这里的代码是工程简化版重点在于演示处理思路实际项目请根据数据特征和精度要求做适配。import numpy as np import pandas as pd from scipy import stats def mad_outlier_detect(values, threshold3.0): 基于中位数和MAD的粗差探测 values: array-like, 高程观测值 threshold: Z分数阈值, 默认3.0 return: 布尔数组, True表示正常, False表示粗差 arr np.asarray(values, dtypefloat) median_val np.median(arr) mad_val np.median(np.abs(arr - median_val)) if mad_val 1e-12: # 若MAD接近0, 说明数据分布过于集中, 改用标准差 mad_val np.std(arr) # Z分数: 0.6745是正态分布下MAD与sigma之间的转换系数 z_scores 0.6745 * (arr - median_val) / mad_val return np.abs(z_scores) threshold def sigma3_outlier_detect(values): 3σ准则粗差探测 arr np.asarray(values, dtypefloat) mean_val np.mean(arr) sigma_val np.std(arr, ddof1) return np.abs(arr - mean_val) 3 * sigma_val def dixon_outlier_detect(values, alpha0.05): Dixon检验法(Q检验), 一次检测一个极值 values: array-like alpha: 显著性水平 return: 布尔数组 arr np.asarray(values, dtypefloat) n len(arr) if n 3 or n 30: raise ValueError(Dixon检验适用于样本量3到30的数据) # 升序排列 sorted_idx np.argsort(arr) sorted_vals arr[sorted_idx] # Q值计算, 这里以检测最大值为例 Q (sorted_vals[-1] - sorted_vals[-2]) / (sorted_vals[-1] - sorted_vals[0]) # Dixon临界值表(简化版, 实际使用请查完整表) # 这里给出n5~10时alpha0.05的临界值 crit_table {3: 0.970, 4: 0.829, 5: 0.710, 6: 0.628, 7: 0.569, 8: 0.608, 9: 0.717, 10: 0.487} if n in crit_table: crit_val crit_table[n] else: # 超出简化表范围时, 回退到3σ准则 print(样本量超出简化临界表, 回退到3σ准则) return sigma3_outlier_detect(arr) # 标记最大值是否异常 normal_idx np.ones(n, dtypebool) if Q crit_val: normal_idx[sorted_idx[-1]] False return normal_idx def ensemble_detect(df, value_col, methodsNone): 多方法联合探测, 汇总每个方法的判定结果 df: DataFrame, 包含点号和观测值 value_col: str, 高程值所在列名 methods: list, 需要执行的方法名列表 return: 带标记列的DataFrame if methods is None: methods [mad, sigma3, dixon] arr df[value_col].values result pd.DataFrame(indexdf.index) result[点号] df[点号].values result[高程] arr result[正常_mad] mad_outlier_detect(arr) result[正常_sigma3] sigma3_outlier_detect(arr) # Dixon方法每次迭代检测一个极值, 这里循环检测直到无异常 dixon_normal np.ones(len(arr), dtypebool) if dixon in methods: working arr.copy() working_idx np.arange(len(arr)) while True: current_normal dixon_outlier_detect( working[working_idx]) if not np.all(current_normal): current_anom working_idx[~current_normal] dixon_normal[current_anom] False working_idx working_idx[current_normal] if len(working_idx) 3: break else: break result[正常_dixon] dixon_normal # 综合判定: 至少两种方法异常则标记为需复核粗差 methods_col [c for c in result.columns if c.startswith(正常_)] abnormal_count result[methods_col].shape[1] \ - result[methods_col].sum(axis1) result[异常方法数] abnormal_count result[综合标记] abnormal_count 2 return result代码的逻辑不复杂几个地方需要解释。MAD方法中的0.6745系数很多人第一次看到会问为什么是这么个常数。因为在标准正态分布下MAD的期望值是0.6745σ所以用0.6745乘以实际MAD就得到了σ的估计值这样Z分数就能和常规的σ度量对齐。这也是稳健统计中的标准做法。Dixon检验的临界值表我只写了一个简化版本实际使用中必须根据样本量查完整的Dixon临界值表。Dixon检验的临界表有严格的理论依据不同样本量、不同显著性水平下临界值差异很大直接用一个固定阈值会导致误判。联合探测的逻辑是不单一依赖任何一种方法而是把多种方法的判定结果汇总当至少两种方法同时判为异常时才标记为“需复核粗差”。这样做的原因是不同方法的适用条件和失效模式不同3σ 在大粗差情况下容易漏判MAD对于中位数的依赖在数据本身分布不规整时可能出问题Dixon仅适合小样本。联合判定可以有效减少单方法的误判风险。3.3 阈值参数该怎么取参数设置是粗差探测里最容易纠结的地方。阈值取得太大粗差漏网取得太小正常极值被误杀。没有一个通用万能的值必须结合数据和精度要求来定。中位数MAD法的阈值我日常使用3.0作为默认值。如果项目精度要求高、外业观测质量一般我会放宽到2.5尽量不放过可疑点宁可多复核几个。如果数据量大、外业质量稳定就收紧到3.5避免误删。3σ准则的3倍系数是约定俗成的一般不建议调。但要注意n小于10时这个系数的理论意义会弱化最好只用它做参考而不是最终裁决。Dixon检验的显著性水平alpha通常取0.05这是统计学的常规选择。如果数据量特别小比如5个点以内取0.10也可以接受提高检测灵敏度。还有一个项目上的参数是“最大剔除比例”。比如控制点总共12个不能无限制地剔除粗差否则控制网结构就破坏了。我会设定一个阈值比如不超过总点位的30%。一旦超过说明不是个别粗差问题而是数据整体质量问题甚至观测方案有问题需要回到外业重新观测而不是在内业继续“筛数据”。4. 实际案例一组模拟高程数据的完整处理过程光说原理和代码可能还是觉得隔了一层。我构造一组模拟数据完整走一遍流程展示每一步的输入、输出和处理逻辑。假设某工程项目的GNSS高程拟合需要12个控制点点号从P01到P12高程值如下单位米点号高程值P0145.632P0244.789P0346.101P0447.250P0543.987P0645.120P0753.240P0846.005P0942.778P1044.356P1145.881P12120.350先肉眼扫一眼P12的120.350明显异常但P07的53.240有没有问题以及是否有更隐蔽的小幅粗差肉眼很难判断。用上面的代码跑一遍联合探测得到结果点号高程MAD判定3σ判定Dixon判定异常方法数综合标记P01-P06, P08-P11正常区间正常正常正常0否P0753.240异常异常正常2需复核P12120.350异常异常异常3需复核P12毫无疑问被所有方法标记为粗差这种情况直接剔除没有争议。P07的情况比较复杂。中位数和3σ都认为异常但Dixon检验认为正常。为什么会这样因为Dixon检验是基于排序位置关系的P07虽然是最大值但和次大值P0447.250之间的差距相比于整体极差120.350-42.778并不算大。也就是说P12的存在掩盖了P07的异常程度。当P12被剔除后Dixon检验的Q值会发生明显变化P07就可能被判定为异常了。这个案例很好地说明了为什么粗差探测要“迭代进行”而不是“一锤定音”。剔除P12后需要在剩余11个点上重新做一轮探测。第二轮结果点号高程MAD判定3σ判定Dixon判定异常方法数综合标记P0753.240异常异常异常3需复核P07在第二轮中被确认是粗差剔除后剩余10个点的数据才是可以进入平差计算的有效值。再算一下剔除前后的统计量变化剔除前: 均值 52.716, 标准差 21.584 剔除后: 均值 45.372, 标准差 1.349可以很清楚地看到两个粗差点把均值拉高了约7米标准差放大了16倍。如果用未剔除的数据直接做高程拟合结果完全不可用。这就是粗差探测存在的根本价值——它守护的是整个数据链条的入口质量。5. 粗差探测与常见的系统性误差如何区分这一节想讲一个容易被忽略的问题粗差不等于异常异常也不一定是粗差。工程实践中高程数据的异常来源大致分三类第一类是粗差也就是观测过程中的人工疏忽或仪器故障比如读数错误、记录错误、天线高量错。这类误差的特点是孤立、随机、没有空间和时间规律处理方式就是剔除。第二类是系统性偏差比如水准尺零点误差、GNSS高程基准不一致、潮汐改正不足。这类误差的特点是成片出现所有点的高程值整体偏移在统计检验中不会表现为离群值单纯做粗差探测根本查不出来需要通过差分、付里叶分析或外部基准比对来发现。第三类是点位代表性问题这个最容易和粗差混淆。比如某个GNSS控制点建在楼顶楼顶到地面之间有很厚的混凝土层GNSS测的是楼顶大地的相位中心水准联测的是楼下地面点二者之间的高程差被忽略了这个点的高程就会存在一个几十厘米到几米的偏差。从统计上看这个点和整体数据明显不和谐容易被判定为粗差。但它不是观测粗差而是控制点选址或联测方案的问题。错误的做法是直接把点剔除然后当无事发生正确的做法是先排查这个点的来源确认偏差原因后再决定是重测、改用其他点还是在平差中降低权重。我遇到过这样一个项目某次高程拟合数据整体质量不错但有一个控制点反复被各种方法判为粗差。一开始也以为是观测问题后来调出外业记录才发现这个点所在的建筑物在观测期间临时施工影响了卫星信号接收。这个施工干扰导致的系统性误差不是单点观测错误而是外部环境影响。如果直接把这个点当作粗差剔除就会丢失一个有效控制点削弱控制网强度。最终的处理方式是保留该点但降低权重同时在拟合模型中加入该时段误差的补偿项。结果整体精度反而提升了。这个案例给我的教训是粗差探测工具给出的是“异常清单”不是“剔除清单”。工具标记了异常之后必须有人去查清原因。自动剔除 人工复核永远是测绘内业最可靠的工作模式。6. 常见问题与排查技巧实录6.1 大粗差掩盖小粗差怎么处理这是我最常被问到的问题。一组数据里既有明显的巨大粗差也有小幅度的可疑点由于大粗差把均值和方差都污染了小可疑点往往无法被检测出来。处理策略是分步迭代先剔除最极端的点然后重新计算统计量再做下一轮探测。每轮只剔除一个点不做批量剔除。这样虽然效率低一些但能避免“误伤”——因为批量剔除时统计量变化后原本正常的点可能会在新一轮统计中被误判为异常。经验法则一轮探测只处理一个最可疑点剔除后重新探测直到没有新的异常点出现。6.2 样本量太小少于15个用什么方法控制网的高程点一般不多10个点左右的常见情况。此时3σ准则的可靠性会下降因为样本量少σ本身的估计就不稳定。我的做法是把中位数MAD当作主要判定方法Dixon检验当作复核最后用3σ做参考。不要只看单一方法的结论交叉验证比任何单一方法都可靠。6.3 正常区域内高点位误差被误判为粗差怎么办这种情况在起伏地形比较多地表本身的高差变化很大某些天然高点并非观测粗差但在统计上确实偏离整体均值显著。关键是引入空间信息作为辅助判断。如果某个点的高程值虽然偏离整体均值但与该点周边的点位高程趋势一致比如同一条剖面线的其他点也在抬升那么它更可能是地形因素而不是粗差。反之如果该点与周边趋势明显矛盾比如周边都是平坦地形突然冒出一个孤立的高点这才需要怀疑粗差。所以一个完善的工具不能只看高程值的一维分布最好结合平面坐标或剖面信息。这也是为什么我在工具架构里加入可视化环节而不是纯统计输出的原因。6.4 粗差探测结果平差后仍不满意怎么办有一个容易被忽视的场景粗差探测是通过统计方法找出离群值但有些粗差并不表现为离群值。比如系统性的高程基准错误所有点的高程都差一个常数统计上完全正常但对最终绝对高程精度的影响很大。所以粗差探测只是数据质量控制的中间环节不是终点。数据进入平差后还要看平差结果的各种检验指标单位权中误差是否异常、残差是否有系统性趋势、网点间的闭合差是否超限。如果这些指标异常即使粗差探测通过了也要回头重新审视数据考虑是否存在系统性误差。7. 从粗差探测到数据质量闭环管理的一点体会做测绘内业这些年我越来越觉得粗差探测不是一道独立的工序而是数据质量管理闭环中的一环。这个闭环包含几个环节外业观测规范的执行、数据入库时的完整性检查、粗差探测与剔除、敏感性分析、平差质量评估、成果精度评定。每一环都有人为误差的可能性粗差探测做的就是拦截其中一部分问题但后面还有平差检验、成果复核等多道防线。我见过很多同行在粗差探测上过度依赖统计方法把工具跑一遍就完事。实际上工具给出的是概率意义上的异常标记它的判断建立在一系列统计假设之上。数据是否满足这些假设样本量是否足够支撑方法的灵敏度异常点是否有非粗差原因这些都需要人的判断。所以最终我的工作习惯是统计方法负责“圈定可疑范围”人工负责“最终定性”两者互为补充。圈定可疑范围时宁宽勿严避免漏网人工定性时宁问到底不留隐患。这个工具后续还可以扩展比如加入空间聚类算法来识别区域性异常、用机器学习模型学习历史数据的异常模式、对接平差软件自动生成剔除标记和精度报告。但核心的东西不会变粗差探测的本质是对数据质量的敬畏是对后续平差模型可靠性的保护。如果你正在做高程数据处理建议先把这几种基础方法吃透用自己的数据多跑几遍感受不同参数下的输出差异。统计方法不是玄学它的每个阈值背后都有数学逻辑理解了逻辑再去用工具才不会在异常情况面前手足无措。