全球导航卫星系统(global navigation satellite system，GNSS)具有速度快、全天候、自动化程度高等优点，已广泛应用于大型工程的安全施工与变形监测中^[1-2]。但是，由于变形监测工作环境的复杂性和不确定性，GNSS信号容易受到强烈干扰，造成变形监测时间序列不可避免地出现粗差^[3]，而且粗差呈现出变化波动范围较大的特征^[4]。这些大规模粗差往往会导致无法准确反映变形体的实际变形，对其稳定性分析及评价甚至灾害预警应急都会带来不可忽视的影响。

GNSS时间序列的粗差探测一直都是国内外测量数据处理的热点问题之一，可分为统计量粗差探测法和非统计量粗差探测法两大类型。统计量粗差探测法是基于观测数据的误差分布特征，构造与之相适应的粗差判别统计量，采用假设检验对粗差进行探测^[5-6]。该类方法应用较为广泛，最常用的是3σ和四分位间距(inter-quartile range，IQR)两种统计量。3σ法作为一种经典的粗差探测方法，结合最小二乘法得到观测数据的残差序列，并对其进行粗差探测。该方法在先验模型准确、粗差数据较少情况下的效果会比较好。但最小二乘法自身抗御粗差的能力有限，造成经其计算得到的残差容易受到粗差的污染，最终粗差探测效果往往不佳^[6]；同时，最小二乘法的拟合效果在很大程度上依赖于先验模型的准确度，当实际模型与先验模型存在偏差时，3σ法探测效果也必然会受到影响^[7-8]。IQR法是一种稳健的统计分析方法，其中位值和标准四分位间距不易受粗差影响，用它们分别代替统计方法中的平均值和标准偏差来进行总体估计，再通过稳健Z比分数统计量来探测数据中的极端异常值（粗差）。IQR法对于偏离假定模型的观测数据具有良好的适应性，但对于离散度较大的观测数据，其四分位间距随之变大，使得稳健Z比分数统计量对于偏离度较小的粗差判定不够敏感^[9]。非统计量粗差探测法主要以小波分析为主，它从小波变换后小波系数具有的模量极大值特性出发，寻找出观测数据的突变点来实现粗差探测^[10-11]，其在数据统计规律不明确的情况下，对粗差含量少、偏离度较大的数据序列具有一定的探测效果。但小波变换信号的奇异性分为两种情况：一是信号幅值发生突变，这是信号中的粗差；二是信号虽然光滑，但其一阶微分不连续也会引起小波变换后信号发生突变，因此难以判断序列奇异点是否是粗差点^[12]，从而导致该方法在实际应用中具有一定人为主观性，难以适应高可靠性的GNSS变形监测实际需求。

本文以3σ统计量粗差探测法为基础，提出了一种基于小波分析的改进型3σ粗差探测方法，既解决了GNSS变形监测时间序列在非理想先验模型和大规模粗差影响下无法准确提取变形趋势的难题，也改善了3σ法因残差受到粗差污染影响而无法准确预估中误差所带来的影响，提高了3σ法对GNSS变形监测时间序列粗差探测的适应性。

1 常规统计量粗差探测方法

1.1   3σ粗差探测法

3σ法不仅要求观测数据所含误差满足高斯分布，而且也要求重复测量次数达到一定数量级。根据随机误差高斯分布，观测数据误差分布在±3σ以外被认为是粗差。在GNSS变形监测时间序列粗差探测时，一般结合最小二乘法得到去除趋势项的残差序列R_n，利用贝塞尔公式计算其标准差σ。如果在测量数据中有大于3σ的残差，即：

则被认为是粗差，置信水平为99.7%。

1.2   IQR粗差探测法

对于一组离散的观测数据序列，将其从小到大排列后，可以采用中位值和标准化IQR对数据的集中度和离散度进行量化。中位值是从小到大排列后的观测数据序列的中间值，IQR则是低四分位数和高四分位数的差值（低四分位数和高四分位数分别是观测数据序列由小到大排列后第25%和75%的数值），表示为：

式中，Q₁和Q₃分别为低四分位数和高四分位数。IQR法下的稳健Z比分数统计量为：

式中，X_i是观测数据；median(X_i)表示观测数据的中位数；s是标准化IQR，它是观测数据变异性的度量，s=0.741 3IQR。当统计量|Z|>3，则认为对应的观测值为粗差，置信水平为99%。

2 改进型3σ粗差探测法

2.1   小波分析原理

小波分析是近年来数学分析中发展起来的信号处理方法，在时域和频域都具有表征信号局部特征的能力，能够对原始信号的各频率成分进行分解，从而实现有用信息和干扰信息的有效分离^[13]。

对于任意函数$f\left( x \right) \in {L^2}\left( R \right)$，连续小波变换为：

式中，a是用于控制伸缩的变量，称为尺度因子；b是用于控制平移的变量，称为平移因子；$\psi \left( x \right)$为小波基函数；$\bar \psi \left[ {\left( {x - b} \right)/a} \right]$为$\psi \left[ {\left( {x - b} \right)/a} \right]$的复共轭函数。

若尺度函数$\varphi \left( t \right) \in {V_0}$是标准正交函数，设有信号f(t)，其最高频率为w，采用小波分析常用的Mallat算法对其分解^[14]，得到信号的不同频率成分：

式中，等号右边第1项为信号频率不超过2^-jw的低频部分；第2项为信号频率在2^-jw和2^-j+1w之间的高频部分；j为小波分解的层数；c_j，k为最大尺度上小波变换的低频系数；d_j，k为各层分解的高频系数。根据小波分析多分辨率的特点，从而实现信号中有效信息和误差干扰的分离。

2.2   基于小波分析的改进型3σ粗差探测法

在对粗差进行探测时，GNSS变形监测时间序列可以看成是一个被噪声和粗差污染的趋势信号，时效(趋势)变化部分表现为低频率(长周期)的变化^{[7, 15]}。故通过变形监测数据在频域方面的差异分离变形趋势和误差，比传统的时域分析不易受到先验模型和粗差数据的影响。同时，由于测量数据误差信息主要集中在较高频，高频系数中包含了测量条件和测量仪器引起的误差，尤其当GNSS监测数据在某些时段受到的干扰较大时，其误差主要集中在第1层高频系数中^[16]。

基于上述小波分析特性，采用小波分析方法改进传统的3σ法粗差探测方法，其主要思路是通过小波分解提取低频系数，对其进行重构得到变形趋势，较传统3σ法得到的变形趋势更为准确；利用小波分解出的高频系数计算得到观测数据的中误差估值σ，其不易受粗差大小和分布的影响。在此基础上，利用传统的3σ法进行粗差探测，具体改进策略如下：

1) 偏离度较大的粗差初步探测及插补。采用3σ法直接对原始GNSS变形监测时间序列进行粗差探测，对于探测到的粗差点X_i，选用平均值法进行插补，插补方法为：

2) GNSS变形监测时间序列的小波分解。需要考虑小波基函数在正则性、紧支撑度和消失矩等方面的要求，选择合适的小波基函数，再结合变形监测时间序列的复杂程度确定分解层次$j$（$j$一般取3~8），对变形监测时间序列X_n进行小波分解^[17]。

由于在小波基函数选择方面至今尚无统一标准，故本文在选择小波基函数时，依据经验选取变形分析中常用的sym7小波基函数，因其集紧支撑、正交性和对称性于一体，具有良好的时频特性，能够有效分离信号中的随机误差和有用变形信息，并在重构信号时能获得更好的平滑效果^[17-18]。在确定最佳分解层次时，采用均方根误差(root mean square error，RMSE)作为判断指标。

3) GNSS变形监测时间序列中误差σ的估计。利用小波分解得到第1、2层高频系数，通过计算其绝对值中值来估计中误差σ：

式中，${{w_{d1}}}$、${{w_{d2}}}$分别为小波分解第1、2层小波高频系数；$p$为中值参数，一般取0.674 5；${\alpha _1} = \sqrt {{E_1}} /\left( {\sqrt {{E_1}} + \sqrt {{E_2}} } \right)$；${\alpha _2} = \sqrt {{E_2}} /\left( {\sqrt {{E_1}} + \sqrt {{E_2}} } \right)$；E₁和E₂分别为第1、2层高频系数对应能量。

4) GNSS变形监测时间序列变形趋势${{X'}_n}$的提取。利用小波分解的低频系数进行计算重构，得到变形监测时间序列的变形趋势${{X'}_n}$。

5) GNSS变形监测时间序列粗差的判定。根据求得的中误差估值σ和提取的变形趋势${{X'}_n}$，利用3σ法进行粗差探测，判断式如下：

当${q_i}$ > 3时，则认为${{X_i}}$是粗差点。

3 模拟实验结果与分析

3.1   实验数据与处理

为验证改进算法粗差探测效果的有效性，采用GNSS变形分析中常用的坐标时间序列模型^[8]，模拟生成不含粗差的原始数据，函数表达式如下：

式中，$i$为坐标历元时刻标识；$x\left( {{t_i}} \right)$为GNSS测站某一分量${{t_i}}$时刻的坐标；$b$为截距；${v_0}$为线性速度；m₀为谐波个数；${{a_m}}$和${{b_m}}$是频率为${{f_m}}$时周期项的振幅；${r_{{t_i}}}$为随机噪声，即${r_{{t_i}}} \in N\left( {0, {\sigma ^2}} \right)$。由于变形监测工程项目现场的实际情况受到多种因素影响，凭经验建立的先验理论模型与实际情况可能存在一定的偏差^[10]。为了更加充分地验证改进型3σ法的粗差探测效果，分别基于理想先验模型和非理想先验模型生成两套不含粗差的原始数据。其中，理想先验模型生成的模拟数据的参数为:$b = 5, {v_0} = 2, {m_0} = 2, {a_1} = 5, {b_1} = 5, {a_2} = 3, {b_2} = 3$，f₁和f₂分别为周期为1 a和半年对应的频率。在非理想先验模型生成的模拟数据中，另外加入高频周期项，其参数为a₃=6，对应频率f₃=7。

在模拟粗差数据时，考虑到变形监测工作环境下粗差数量多、大小分布广等特点，实验首先采用标准差为10σ（σ为随机噪声标准差）的正态分布模拟得到一组随机误差序列，然后将大于3σ的误差序列随机加入到不含粗差的原始数据，最终得到粗差污染的GNSS变形监测时间序列。实验数据由Matlab平台生成，历元间隔4 h，共2 191个历元。模拟的粗差总数为145，占总观测历元的6.62%，实验数据如图 1所示。

图  1  两种实验数据

Figure  1.  Two Kinds of Experimental Data

下载: 全尺寸图片 幻灯片

利用改进型3σ法对两种实验数据进行处理，通过筛选比较确定分解层次分别为7和6。表 1表示两种模型情况下，常规3σ法和改进型3σ法的趋势拟合度和中误差估值偏离度比较。其中，趋势拟合度用决定系数表示，决定系数越接近1，说明趋势拟合越准确。通过表 1可以看出，两种模型下，改进型3σ法在趋势拟合精度和中误差估值精度两个方面均优于常规3σ法，尤其在模型由理想转变为非理想时，3σ法效果明显变差，而改进型3σ法则表现出较强的适应性。

表  1  两种模型下趋势提取拟合度和中误差估值偏离度

Table  1.  Fitting Degrees of Trend Extraction and Deviation of Median Error Estimation of Two Models

先验模型	方法	趋势拟合度	中误差估值偏离度/%
理想	3σ法	0.867	38.6
	改进型3σ法	0.946	9.2
非理想	3σ法	0.711	84.4
	改进型3σ法	0.929	14.0

下载: 导出CSV 

| 显示表格

3.2   结果与讨论

利用3种方法（3σ法、IQR法、改进型3σ法）对理想先验模型生成的实验数据进行处理，其结果如图 2(a)~2(c)所示。可以看到，3σ法可以探测到GNSS变形监测时间序列的大部分粗差，IQR法只能探测到少部分粗差，而改进型3σ法则几乎可以探测到所有粗差点。

图  2  两种模型生成实验数据的粗差探测结果

Figure  2.  Gross Error Detection Results of Experimental Data Generated by Two Models

下载: 全尺寸图片 幻灯片

3种粗差探测方法对非理想先验模型生成的实验数据处理结果如图 2(d)~2(f)所示。可以明显看出，较理想先验模型生成的实验数据处理结果，3σ法粗差探测效果明显变差，IQR法探测效果变化不大，只能探测到少部分粗差，改进型3σ法依然保持优异的探测效果，几乎可以探测到所有粗差点。

为简洁描述，以下仅针对非理想先验模型生成数据的处理结果进行详细分析。

1) 统计分析

为定量、直观地比较3种粗差探测方法的探测效果，对3种方法探测粗差的大小和数量进行统计分析，结果如表 2所示。可以看出，3σ法、IQR法和改进型3σ法的粗差总探测率分别为59.3%、32.4%以及90.3%。在粗差探测的总体效果方面，3σ法在7σ以上粗差探测率为100%，在5σ~7σ粗差探测率为50%，而在3σ~5σ粗差探测率仅为15.9%。由于模型为非理想模型，故最小二乘法提取变形监测时间序列的变形趋势较理想先验模型更加不准确，造成残差序列仍含有异常趋势项，再加上粗差的干扰，使得求取的中误差估值更大，从而导致对偏离度较小粗差的探测不够敏感；IQR法随着粗差偏离度变小，其探测效果逐渐变差，甚至在3σ~5σ粗差探测率仅为0.7%，这主要是由于IQR法是从观测数据总体离散度出发去探究其极端离异值，无法从数据发展特性来分析内在规律，导致对于偏离度较小的粗差探测效果不理想；而改进型3σ法在5σ以上粗差探测率为100%，并且在3σ~5σ粗差探测率上同样能达到68.2%的显著效果，这主要是因为改进型3σ法能够利用小波分析得到的GNSS变形监测时间序列的低频信息更准确地提取变形趋势，同时由于小波分析得到的高频信息求取的中误差估值受粗差干扰相对较小，能够更好地趋近真实观测精度，从而提高粗差探测的准确度和敏感度。

表  2  粗差探测统计表/%

Table  2.  Statistics of Gross Error Detection/%

探测方法	按粗差分区的探测率					总探测率
	3σ~5σ	5σ~7σ	7σ~9σ	9σ~11σ	> 11σ
3σ法	15.9	50	100	100	100	59.3
IQR法	0.7	11.4	43.5	86.7	94.7	32.4
改进型3σ法	68.2	100	100	100	100	90.3

下载: 导出CSV 

| 显示表格

2) 粗差探测前后参数估值比较

在最小二乘估计中，粗差的大规模存在会导致参数估计严重偏离其真实值。基于此，利用粗差探测前后参数估值的变化，从另一个角度评估3种方法的粗差探测质量。其中，对于探测到的粗差点采用平均值法进行插补。粗差探测前后参数估值的变化见图 3。可以发现，含粗差的原始数据的参数估值相对于真实值偏离度较大，尤其是参数v₀的估值偏离严重。3种探测方法均使参数估值在一定程度上趋近于理想参数值，但效果各有差别。其中IQR法效果最差，改进型3σ法效果最好，其总体参数估值相对误差最小。

图  3  粗差剔除前后参数估值的相对误差

Figure  3.  Relative Error of Parameter Estimation Before and After Gross Error Elimination

下载: 全尺寸图片 幻灯片

4 工程案例结果与分析

4.1   工程数据与处理

本文基于陕西省渭南市金堆城露天矿GPS变形监测预警系统数据开展工程案例研究。该系统于2014年1月建成，主要针对采场边帮、排土场和机车轨道边坡的变形进行监测。由于矿区周边植被茂盛，现场高压电线以及信号基站密布，对变形监测信号接收终端造成了一定干扰，监测数据序列含有一些无规律的粗差。此次工程数据选用矿山北部边坡B12 GPS监测点，时间跨度为2014-03-22 12:00—2014-07-20 12:00的Z方向变形数据进行分析，其数据获取方式采用时段解，采样间隔为4 h，共计721历元数，如图 4（a）所示。

图  4  原始变形监测数据及其粗差探测结果

Figure  4.  Original Deformation Monitoring Data and Its Results of Gross Error Detection

下载: 全尺寸图片 幻灯片

4.2   结果与讨论

从图 4（a）可以看出，原始变形监测数据中的粗差在时间域和大小上分布范围较广。采用3种方法分别对其进行粗差探测，其结果如图 4（b）~4（d）所示，其中改进型3σ法确定的分解层次为8。可以看出，3种方法均可以探测出一定数量的粗差数据，尤其是偏离度较大的粗差数据。但在偏离度较小的粗差数据探测效果方面，改进型3σ法明显优于3σ法和IQR法。

由于仅仅根据识别粗差的数量来评价方法的有效性存在一定的缺陷，因而本文通过粗差剔除前后的原始变形监测数据的RMSE来评价其粗差探测效果，结果如表 3所示。为客观评价3种方法的探测效果，利用奇异谱分析提取的变形趋势作为参照^[19]。可以看出，三者都具有一定的粗差探测能力，但效果依然是改进型3σ法明显优于其他两种方法。所以无论从粗差探测的数量还是质量方面分析，改进型3σ法的效果更优。

表  3  3种方法的粗差探测数量及剔除前后的RMSE

Table  3.  Quantities of Gross Error Detection and RMSE of Before and After Gross Error Elimination of Three Methods

方法	探测数量	RMSE/mm
原始数据	—	9.941
3σ法	39	4.201
IQR法	24	5.472
改进型3σ法	55	3.499

下载: 导出CSV 

| 显示表格

5 结语

本文针对GNSS变形监测应用中监测环境的复杂性和不确定性，考虑到变形监测时间序列中粗差数量多、大小分布广的特点，提出了一种改进的粗差探测方法，一方面解决了GNSS变形监测时间序列在非理想先验模型和大规模粗差影响下无法准确提取变形趋势的难题，另一方面也改善了3σ法因残差受到粗差污染而无法准确预估中误差所带来的影响，在GNSS变形监测时间序列粗差探测效果上表现优异。

从粗差探测结果统计定量分析和探测前后参数估值比较结果来看，改进型3σ法相比常规方法，在保证对复杂情况下GNSS变形监测时间序列探测更好的适应性的同时，对3σ~5σ粗差探测效果具有明显优势，能够更好地适应GNSS变形监测项目的实际需求。

本文利用小波对GNSS变形监测时间序列进行分析时，在小波基函数选取和分解层次上依赖于前人的相关经验和重复的实验筛选，这是改进型3σ法探测粗差的关键环节。今后需要结合GNSS时间序列组成和小波分析，形成系统的小波分析方法，从而达到更高效、更有效的粗差探测效果。