高程数据粗差探测全解析：3σ、MAD、Dixon检验与稳健估计实战

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σ准则的核心问题是均值与标准差被粗差污染，那自然的思路就是改用抗污染的统计量。中位数就是最典型的稳健统计量。

中位数对粗差的抵抗力极强，哪怕一组数据里混入几个极端值，中位数也几乎不受影响。类似地，MAD（Median 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检验在处理小样本（n=3到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是阈值系数，通常取c1=2.5、c2=4.0，可以根据粗差比例适当调整。

这个方法的优势在于，它不只是“找出粗差”，而是把粗差处理融入了平差流程本身。尤其是处理控制网数据时，高程异常点会被自动降权或剔权，最终平差结果更稳健。缺点是计算流程比前面几种方法复杂，需要迭代运算，而且阈值参数需要经验积累。

3. 高程数据粗差探测的系统设计与工具实现

3.1 工具的总体架构与处理流程

基于上面的方法梳理，我设计了一个高程粗差探测工具。这个工具不是只跑一个算法，而是串起一条完整的数据处理链路，从数据读入到最终输出报告，中间每一步都有明确的逻辑和可追溯的记录。

整体流程分为五个环节：

1. 数据读入与规范化处理：接收点号、高程值、可能的平面坐标等字段，统一格式，处理缺失值。
2. 粗差初筛与可视化：计算基本统计量，生成分布直方图和序序图，让数据形态直观可见。
3. 多方法联合探测：依次使用中位数+MAD、3σ准则、Dixon检验，对可疑点进行标记并输出各方法的判定结果。
4. 人工复核与决策：结合点位空间分布和项目背景，确认最终剔除名单。
5. 输出处理结果与报告：生成剔点前后对比表、统计指标变化、XY方向分布图。

其中第二步的可视化是我特别想强调的。很多人在粗差探测时只看数字结果，忽略了图形直觉。高程数据在空间上的分布是有规律的，一个点在高程异常序列中可能不明显，但在平面位置和高程值的三维关系里往往一目了然。比如某个GNSS水准控制点距离其余点很远，偶然高程误差又大，这个点在统计检验中可能不达标，但实际可能是点位本身代表性不足，而不是观测粗差，两者的处理方式完全不同。

3.2 核心代码实现

下面给出我用Python实现的核心算法代码，包括多方法联合探测模块。注意，这里的代码是工程简化版，重点在于演示处理思路，实际项目请根据数据特征和精度要求做适配。

import numpy as np import pandas as pd from scipy import stats def mad_outlier_detect(values, threshold=3.0): """ 基于中位数和MAD的粗差探测 values: array-like, 高程观测值 threshold: Z分数阈值, 默认3.0 return: 布尔数组, True表示正常, False表示粗差 """ arr = np.asarray(values, dtype=float) 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, dtype=float) mean_val = np.mean(arr) sigma_val = np.std(arr, ddof=1) return np.abs(arr - mean_val) < 3 * sigma_val def dixon_outlier_detect(values, alpha=0.05): """ Dixon检验法(Q检验), 一次检测一个极值 values: array-like alpha: 显著性水平 return: 布尔数组 """ arr = np.asarray(values, dtype=float) 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临界值表(简化版, 实际使用请查完整表) # 这里给出n=5~10时alpha=0.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, dtype=bool) if Q > crit_val: normal_idx[sorted_idx[-1]] = False return normal_idx def ensemble_detect(df, value_col, methods=None): """ 多方法联合探测, 汇总每个方法的判定结果 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(index=df.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), dtype=bool) 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(axis=1) 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，高程值如下（单位：米）：

先肉眼扫一眼，P12的120.350明显异常，但P07的53.240有没有问题，以及是否有更隐蔽的小幅粗差，肉眼很难判断。

用上面的代码跑一遍联合探测，得到结果：

P12毫无疑问被所有方法标记为粗差，这种情况直接剔除，没有争议。

P07的情况比较复杂。中位数和3σ都认为异常，但Dixon检验认为正常。为什么会这样？因为Dixon检验是基于排序位置关系的，P07虽然是最大值，但和次大值P04（47.250）之间的差距，相比于整体极差（120.350-42.778）并不算大。也就是说，P12的存在掩盖了P07的异常程度。当P12被剔除后，Dixon检验的Q值会发生明显变化，P07就可能被判定为异常了。

这个案例很好地说明了为什么粗差探测要“迭代进行”而不是“一锤定音”。剔除P12后，需要在剩余11个点上重新做一轮探测。

第二轮结果：

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. 从粗差探测到数据质量闭环管理的一点体会

做测绘内业这些年，我越来越觉得粗差探测不是一道独立的工序，而是数据质量管理闭环中的一环。

这个闭环包含几个环节：外业观测规范的执行、数据入库时的完整性检查、粗差探测与剔除、敏感性分析、平差质量评估、成果精度评定。每一环都有人为误差的可能性，粗差探测做的就是拦截其中一部分问题，但后面还有平差检验、成果复核等多道防线。

我见过很多同行在粗差探测上过度依赖统计方法，把工具跑一遍就完事。实际上，工具给出的是概率意义上的异常标记，它的判断建立在一系列统计假设之上。数据是否满足这些假设，样本量是否足够支撑方法的灵敏度，异常点是否有非粗差原因，这些都需要人的判断。

所以最终我的工作习惯是：统计方法负责“圈定可疑范围”，人工负责“最终定性”，两者互为补充。圈定可疑范围时，宁宽勿严，避免漏网；人工定性时，宁问到底，不留隐患。

这个工具后续还可以扩展，比如加入空间聚类算法来识别区域性异常、用机器学习模型学习历史数据的异常模式、对接平差软件自动生成剔除标记和精度报告。但核心的东西不会变：粗差探测的本质是对数据质量的敬畏，是对后续平差模型可靠性的保护。

如果你正在做高程数据处理，建议先把这几种基础方法吃透，用自己的数据多跑几遍，感受不同参数下的输出差异。统计方法不是玄学，它的每个阈值背后都有数学逻辑，理解了逻辑再去用工具，才不会在异常情况面前手足无措。