在GPS变形监测中，可实时获取监测点的三维坐标。在短基线差分解算中可以很好地消除电离层延迟、对流层延迟、轨道误差等空间相关性强的误差。但通过差分解算对基站多路径等空间相关性不强的误差影响却不能很好消除，因此变形监测的动态定位精度一般只达到cm级^[1]。经验模态分解^[2](empircal mode decomposition, EMD)是一种自适应时域分析方法，可以将原始信号从高频到低频逐级分解，形成一系列本征模态(intrinsic mode function，IMF)，适用于分析非平稳信号^[2]。戴吾蛟等^[4]使用EMD滤波去噪削弱多路径影响，并提出了一种分解层数的选择标准，但分解过程中如果遇到数据中某些频率尺度的分量是不连续的，会影响不同尺度信号的有效分离，即出现模态混叠的现象。通过在信号中加入白噪声的总体经验模态分解(ensemble empirical mode decomposition，EEMD)可以减小模态混叠，但EEMD存在的噪声残留及模态数量不等等问题。Torres等^[4]提出了自适应噪声的经验模态分解(complete ensemble empirical mode decomposition，CEEMD), 在分解的每一阶段添加一个特定白噪声，并计算一个唯一残差以得到每个IMF，以此完成精确的原始信号重构。本文采用CEEMD对坐标残差序列进行分解，采用主成分分析结合KS检验判断有效信号层，并将处理结果与传统方法进行对比。

1 坐标域的多路径信号CEEMD分解方法

CEEMD分解信号后将获得一系列从高频到低频排列的细节信号成分，即本征模态IMF，其获得第一个IMF的方式与EEMD相同，即向原始分解信号(此处信号为坐标残差序列X)中加入I次均值为零的单位方差白噪声，产生I组不同的信号$\mathit{\boldsymbol{\tilde X}}$，之后对信号$\mathit{\boldsymbol{\tilde X}}$进行经验模态分解，取第一层模态${\widetilde {{\rm{IMF}}}_{\rm{1}}}$ (I组)的平均值，得IMF₁，表达式为^[5]：

$ {\rm{IM}}{{\rm{F}}_1} = \frac{1}{I}\sum\limits_{i = 1}^I {\widetilde {{\rm{IMF}}}_1^i} $

(1)

所得余项 r₁= X － IMF₁。定义算子E_j(·)表示对输入信号的EMD分解过程，并获得第j个本征模态。对 r₁作式(2)处理：

$ {\mathit{\boldsymbol{r}}_{k + 1}} = {\mathit{\boldsymbol{r}}_k} + {\beta _k} \cdot {E_k}(\mathit{\boldsymbol{\omega }}{\;^j}) $

(2)

式中，k=1…N，j=1…I，β为信噪比，ω为零均值单位方差白噪声。对加入的I组噪声分别进行EMD分解并获取其第k个本征模态，分别加入 r_k, 当k=1时，即得I组r₂。

再对 r₂进行EMD分解，取第一层模态的平均，得CEEMD分解的第二个本征模态 IMF₂, 表达式为：

$ {\rm{IM}}{{\rm{F}}_2} = \frac{1}{I}\sum\limits_{i = 1}^I {{E_1}(\mathit{\boldsymbol{r}}_2^i)} $

(3)

需要注意的是，将高斯白噪声经EMD分解后的第k层模态加入第k阶余项后，进行EMD分解并取第一层模态。以此类推，直至解析信号余项不能被分解。此时原始分解信号表示为：

$ \mathit{\boldsymbol{X}} = \sum\limits_{k = 1}^K {{\rm{IM}}{{\rm{F}}_k} + {\mathit{\boldsymbol{r}}_{k + 1}}\;} $

(4)

式中，K表示分解层数。通过不同IMF信号层的组合可以组成高通、带通、低通滤波器提取关注信号。多路径信号表现为低频特征，存在于高阶模态中，本文采用主成分分析结合KS检验方法确定多路径信号层l, 取l层以后的信号层进行组合，以此实现从原始坐标残差中提取多路径信号。

2 多路径信号提取

在主成分分析中，为得到信号的主成分，利用奇异值分解^[6]：

$ [\mathit{\boldsymbol{V}}{\;^{\rm{T}}}, \mathit{\boldsymbol{S}}, \mathit{\boldsymbol{Q}}] = {\rm{SVD}}\left( \mathit{\boldsymbol{Y}} \right)\; $

(5)

$ {\mathit{\boldsymbol{Y}}_{m \times n}} \approx {\mathit{\boldsymbol{Q}}_{m \times n}}{\mathit{\boldsymbol{S}}_{n \times n}}{\mathit{\boldsymbol{V}}_{n \times n}} $

(6)

式中，Y代表原始信号，此处是IMF所构成的矩阵；Q为特征向量矩阵，每一列代表一个特征向量；V的每一行正交，代表与IMF相联系的主成分分量；S是一个对角阵，S =diag(s₁, s₂, s₃, s₄, …, s_n)。如果Y为可逆阵，则 S对角线上为Y的特征值，且按照从大到小排列。

$ \frac{{\sum\limits_{i = 1}^l {{s_i}} }}{{\sum\limits_{i = 1}^n {{s_i}} }}\; > D $

(7)

式中，D表示噪声贡献率，按式(7)找出级数l，l级之前的本征模态当作噪声，l之后的叠加起来，当作信号，即多路径误差。检验D值是否选取合适，采用Kolmogorov-Smirnov(KS)检验进行验证。白噪声分布符合高斯分布，对提取的高频信号求功率谱，反算其高斯分布参数σ和u，并与符合相同参数的高斯分布数据频率谱进行KS检验^[7]。如果KS检验结果接受，说明残余量(即l级之前的量之和)是符合高斯分布的，即所选取的D值是可接受的。实验时，由于数据中的噪声含量未知，将D值设置较大，通过KS检验时，获取的l层之前的高频信号不可避免地包含部分多路径信号，那么KS检验将不通过, 然后以一定步长减小D值进行迭代，直至KS检验通过。

3 GPS实验数据结果分析 3.1 实验数据

数据采集用中海达接收机，基站距离200 m, 并将接收机安置在稳定的混凝土桩中进行连续3 d的静态数据采集，卫星高度角为10°，采样间隔为1 s。由于基线很短，坐标残差序列主要由多路径误差和高频随机噪声组成。建筑物的结构振动基本周期在[0.05~0.1]N(N为楼层)之间变化^[8]，得到接收机的基本振动频率为2~4 Hz，远高于50 m内反射物造成的多路径效应典型频率(约0.026 Hz)，对提取多路径效应不造成影响^[9]。由于能量在传播过程中会衰减，对50 m外反射物多路径信号不予与考虑^[10]。

3.2 CEEMD实验结果分析

进行两组实验，连续观测2 d(DOY 254~255、255~256)，采用动态双差逐历元模块解算得到监测点X、Y、Z三维坐标，提取其中UTC 16:00:01~18:46:38时段，共10 000个历元数据并剔除3倍中误差以外的粗差作为实验分析数据。在坐标域内对前一天的坐标双差残差进行去噪，以消除高频噪声，从而提取多路径改正模型。利用多路径重复性的特点对第二天的坐标序列进行改正。图 1(由上而下DOY分别为254和255)为监测点X、Y、Z方向连续2 d的原始坐标残差序列(单历元坐标解算值-坐标真值)，可以看出坐标残差序列存在一定的重复性，且第二天的坐标残差序列比第一天提前了约4 min。去除高频噪声前，用最大相关系数法求取平移时间为241 s、244 s、244 s。以滤波前后坐标序列的RMS值来衡量滤波质量，RMS值越小说明多路径改正越明显。以去噪前后坐标序列的相关系数来衡量去噪效果，相关系数越大去噪效果越明显。

以X方向(DOY 254~255)为例，对2 d的数据进行模态分解，得到12个本征模态IMF及剩余项, 如图 3所示, 由上至下依次表示原始观测信号、12个本征模态(按高频到低频依次排列)及剩余项。对本征模态构成的矩阵进行奇异值分解，选取第7层以后的分量进行重构。频谱分析如图 4所示，由上至下分别对应图 3中各IMF的频谱，可见7层以后频率与多路径频率相近。图 2为通过CEEMD提取的低频信号，前后2 d三个方向相似度均大于0.8，具有较高重复性，认为此信号为多路径。以X方向为例，其相关系数为0.847，相对于原始残差序列的相关系数0.510得到了较大的提高，说明残差序列中含有较多随机噪声。

利用第一天的多路径改正模型对第二天的坐标残差序列进行改正, 比较残差序列改正前后的RMS值，如表 1、2所示。可以看出，剔除第一天多路径改正模型后的残差序列三个方向RMS值均小于剔除前的残差序列RMS值。由图 5知，X、Y、Z三个方向的残差序列稳定性有明显提高，残差序列更加平滑，综上说明多路径误差在一定程度上被削弱。图 6为DOY 255减去自身提取的多路径改正模型后的结果，可见具有很强的随机性，根据KS检验，其符合高斯分布，为随机噪声。

3.3 CEEMD与传统方法对比

将CEEMD方法分别与EMD和移动平均法(91历元移动窗口)进行实验对比。图 7为DOY 254~255X方向三种方法实验结果，从上至下分别为第一天多路径信号、第二天多路径信号、经改正后的第二天坐标序列、减去自身多路径的第二天坐标序列，可以看出三种方法提取的多路径信息较为相似。比较表 4、6可知，三种方法提取多路径信号前后2 d的相关系数都大于0.8，保证了模型较高的相似程度。CEEMD去噪效果优于EMD，移动平均去噪效果最佳。利用多路径重复性，分别经三种方法改正后的坐标序列更加平滑。CEEMD和EMD处理结果相近，随机噪声呈现高斯分布，移动平均的随机噪声带有一定的趋势项。由表 3、5知，CEEMD分解结果在X、Y、Z三个方向上RMS都小于EMD结果，说明利用CEEMD的有效性且处理效果优于EMD，其原因是CEEMD通过加入特定辅助白噪声到原始信号，具有均匀频谱特性的白噪声改善了间断数据分量的时间尺度特性^[11]，从而减小模态混叠，提供精确的原始信号重构。较相关系数更大的移动平均方法，CEEMD和EMD的处理效果更优，其原因是移动平均在消噪过程中会将如图 1(a)中1 500历元处跳跃现象数据平均分配到窗口内其他历元上，造成偏离真实多路径信号值。而经验模态可以有效克服随机跳跃现象，剔除瞬间强噪声，在分离过程中异常值只会在低阶模态中反映出来，而多路径信号存在于高阶模态中。

4 结语

本文采用CEEMD对原始坐标残差序列进行分解，采用主成分分析结合KS检验进行有效信号层的判断。实验结果表明，利用CEEMD可以有效从残差序列中分离高频噪声和低频信息，利用提取的当天的多路径改正模型进行多路径误差改正，从而有效削弱多路径误差的影响，提高坐标序列的稳定性。X、Y、Z三个方向的精度提高约20%；另外此方法可以有效地减少经验模态分解EMD中出现的模态混叠几率，且较移动平均法具有更好的处理效果。