基于ICEEMDAN与环境负载的GNSS坐标时序非线性形变去除

王勇，曹慧鹏，李锁，闫勇，杨军

( 1. 天津城建大学 地质与测绘学院, 天津 300384；2. 中国科学院精密测量科学与技术创新研究院 大地测量与地球动力学国家重点实验室, 武汉430077；3. 天津市地质工程勘测设计院有限公司, 天津 300191 )

0 引 言

全球卫星导航系统(GNSS)连续运行基准站坐标时序在导航与位置服务、精密卫星定轨、地质灾害监测等工程和科学应用中具有重要的作用[1-2]. GNSS坐标时间序列可反映测站趋势变化外，也能反映出测站非线性变化[3]. 在参考框架构建时，非线性因素会对参考框架的精度产生一定的影响，研究基准站的非线性形变去除方法对于修正误差模型，提高GNSS基准站位置精度，为研究地球物理、地球环境动态变化的检测和研究提供依据具有重要的意义[4-6]. 姜卫平等[7]研究发现中国区域基准站位移呈现显著的区域性特征，环境负载改正能够削弱中国区域大多数GPS基准站(约70%)垂向(U)以及东(E)方向位移时间序列的非线性形变，其计算环境负载改正由负载改正和GPS坐标时间序列相加得到；陈岸等[8]研究了华北地区的GNSS基站站高程时间序列的非线性变化，研究表明GNSS时间序列存在地区差异性，且在高程方向非线性形变存在显著的年周期和半年周期，水文负载对其的改正效果最好，其仅研究了高程方向的时间序列变化；尹涛涛等[9]使用了经验模态分解(EMD)的方法先将测站的坐标时间序列进行分解，提取出相应的周期部分，然后再进行环境负载的改正，但EMD方法进行坐标时间序列分解时会产生模态混叠现象，而且噪声层中也会含有大量的有效信息，该方法无法对噪声层中的有效信息进行提取，这会造成一定的信息损失[10]. 因此有必要探寻一种新的方法进行GNSS坐标时序三维方向的非构造形变去除，减少时序分解模态混叠现象，提高GNSS时序精度.

本文拟采用改进的自适应噪声总体集合经验模态分解(ICEEMDAN)方法和环境负载改正相结合的方法开展GNSS测站非线性形变去除研究. 首先使用GMIS软件将GNSS坐标时序补充完整并去除粗差，然后使用ICEEMDAN方法对GNSS坐标时序进行分解，使用排列熵算法选取包含噪声和非线性形变的高频分量，最后使用环境负载对高频分量进行去除，与EMD方法和环境负载结合的方法进行去除效果对比.

1 研究数据与研究方法

1.1 GNSS坐标时序

本文采用ICEEMDAN和环境负载相结合的方法开展GNSS测站非线性形变去除，涉及GNSS坐标时序，采用中国大陆构造环境监测网络(CMONOC)提供的北京4个站点坐标时序开展相关研究. BJFS(北京房山)和BJSH(北京十三陵)建立于1995年与1998年，建站较早、观测时间长，BJGB(北京古北口)、BJYQ(北京延庆) 2011年起存档历史观测数据，本文研究GNSS站点坐标时序时间跨度为2011-01-01—2021-05-04. 坐标时序由中国地震局第一监测中心提供(https://www.fmac.ac.cn).

1.2 环境负载数据

极潮、非潮汐海洋负载、大气负载、水文负载等环境负载是造成测站非线性形变的主要原因[11-13]. 本次实验采用大气压潮负载、海潮负载、海洋极潮、以及地球极移4种负载改正产品，由地球潮汐负荷影响与形变监测计算系统ETideLoad4.0计算得到.

大气压潮负载数据采用欧洲中期气候预报中心(https://www.ecmwf.int/) ECMWF-DCDA2006的大气压周日S1、半周日S2、半年Ssa和年周期Sa风潮的0.5°×0.5°全球调和常数格网生成；海潮负载采用AVISO+(https://www.aviso.altimetry.fr/en/home.html)

的FES2014b潮高调和常数模型、经ETideLosd4.0(分潮球谐分析与负荷球谐系数模型构建)生成. 海洋极潮采用IERS2010协议(https://hpiers.obspm.fr/)标准自适应海洋潮汐质量均衡算法生成；地球极移采用IERS2010协议(https://hpiers.obspm.fr/)给出的地球定向参数产品EOPC04生成.

1.3 研究方法

由于EMD方法进行坐标时间序列分解时会产生模态混叠现象，本文采用ICEEMDAN方法进行GNSS坐标时序三维方向的非构造形变去除.ICEEMDAN对测站的坐标时间序列进行分解，重新提取本征模函数(IMF)，在增加高斯白噪声的基础上也会添加相应特殊噪声，也就是该白噪声(WN)能够被EMD分解之后获取的k层IMF，随后借助于获取其中存在着特殊残差(1个)，使得该IMF所对应的定义就是残差信号和有关局部均值所对应的差值[14]. 结果显示，IMF所对应的剩余噪声问题得到很好地改善，且基于EMD生成的IMF中的均值问题也能相应地解决. 经过ICEEMDAN方法分解得到的IMF分层能够更好地体现原有坐标时间序列，避免了有效信息的损失，并且解决了模态混叠的问题.

1.4 评价指标

采用的评价指标有两种，第一个是均方根(RMS)作为分析环境负载对区域GNSS基准站坐标时间序列影响的一个指标，其定义为[15]

式中：n表示观测个数；m(j) 表示环境负载造成的位移.

第二个评价指标采用加权均方根(WRMS)来进行评价，使用负载改正前后GNSS时间序列的加权均方根差值来评估环境负载对GNSS坐标时间序列的影响，其定义为[16]

式中：w(j)=1/sig(j)2；G(j)、sig(j) 分别表示基准站在j时刻的位移(坐标分量或位置相对于参考值的差值)及其精度(不确定度)；n表示观测个数.

计算基准站 W RMS 变 化，WRMS变化=WRMS原始-WRMS改正， W RMS变化大于0表示环境负载改正后GNSS坐标时间序列 W RMS 减小，即非线性形变变化幅度减弱.

2 GNSS坐标时序非线性形变去除

2.1 GNSS坐标时序预处理

GNSS永久参考站自动运行，因为通信链路终端或电源故障，GNSS坐标时间序列中数据会经常丢失，丢失的数据是随机的，可能是3～5天，也可能是3～5个月，这将对数据分析造成很大的干扰，因此实验先使用GMIS[17](GNSS Missing Data Interpolation)软件对数据进行修复，采用KKF(Kriged Kalman Filter)[18]的方法. 插值法又因为存在多路径、轨道异常以及测站自身等因素会导致坐标时间序列产生粗差，因此在使用KKF插值方法对数据进行补充后，再使用3σ准则对粗差进行识别并剔除，以此来获取一个相对“干净且完整”的坐标序列. 完善并剔除粗差后的测站坐标时间序列结果如图1所示(以BJFS站和BJSH站为例).

如图1所示，黑线为原始的时间序列，可以看到存在缺失的部分，而且存在异常值，这些异常值就是由于各种原因造成的粗差，红线为使用KKF法和3σ准则进行插值与去粗差后的结果，能明显的看出缺失部分得到了补充，粗差得到了处理，而且补充后的数据与原始数据的趋势基本一样.

图1 完整去噪后站点时间序列

2.2 基于ICEEMDAN方法的GNSS坐标时序高频信号提取

使用ICEEMDAN方法和EMD方法对测站的坐标时间序列进行分解，可获得由高频到低频的IMF分量，为了更好地区分IMF分量中噪声与非线性分量，使用了排列熵[19](PE)来衡量时间序列的复杂程度. PE具有运算效率高、抗干扰能力强、输出结果简单的优点，用于分析非线性、不规则的信号，对于区分IMF分量的频率具有良好的适用性. 时间序列越规则，其值越小；时间序列越复杂，即包含的噪声和非线性形变越多，其值越大. 本实验参照参考文献[20]将PE值的模糊隶属度设定为0.6. 当时间序列IMF分量大于0.6时，可认为是包含噪声和非线性形变的高频信号；PE值小于0.6时，可看作干净信号.计算每一个IMF分量的PE值，选出高频信号然后和环境负载相结合来减少非线性形变. 得到的每个站点的PE值分布如图2所示.

由图2可知，每个测站的时间序列经ICEEMDAN方法分解后在U方向的IMF分量都是11层，在水平N方向都是8层， E方向都是7层，而且3个方向计算得到的PE值大于0.6的分层都在前4层；经EMD方法分解后在U方向的IMF分量基本都是10层，只有BJGB站的坐标时间序列分解后IMF分量是9层，在N方向都是7层， E方向都是6层，而且3个方向计算得到的PE值大于0.6的分层都是在前3层；然后将两种方法提取出来的高频IMF分量进行非线性形变去除.

图2 ICEEMDAN与EMD分解IMF分量PE值比较

2.3 环境负载对GNSS坐标时间序列RMS影响

将4个测站坐标时间序列使用EMD和ICEEMDAN方法进行分解后，再使用排列熵筛选出包含噪声和非线性形变的高频层，然后分别用ETideLoad4.0计算得到的海潮负载、大气压潮负载、地球极移与海洋极潮等环境负载进行改正，得到的RMS变化情况如表1所示.

由表1可知，环境负载改正对于坐标时间序列的非线性形变具有一定的改正效果，在U方向变化最为明显，E方向变化次之，N方向的变化量最小. EMD方法和环境负载结合改正后U方向的RMS变化范围在5.037 0 mm (BJGB，海洋极潮)～7.663 1 mm (BJSH，海潮负载)，E方向的RMS变化次之，变化范围在0.733 6 mm (BJSH，海洋极潮)～1.890 2 mm (BJGB，海潮负载)，N方向的RMS变化量最小，其范围是0.869 9 mm (BJYQ，地球极移)～1.326 7 mm (BJFS，海潮负载)；ICEEMDAN方法和环境负载结合改正后U方向的RMS变化量在5.269 6 mm (BJSH，海洋极潮)～6.715 0 mm (BJSH，地球极移)，E方向的RMS变化量在0.817 3 mm (BJGB，海洋极潮)～1.697 2 mm(BJGB，海潮负载)，N方向的RMS变化量在0.816 2 mm(BJGB，地球极移)～1.513 0 mm (BJSH，海潮负载).

表1 测站时间序列经环境负载改正后RMS变化 mm

通过分别对比不同环境负载改正对坐标时间序列造成的RMS变化量可以看出： 1)变化最大的是海潮负载改正，EMD方法和海潮负载改正结合中N、E、U 3个方向RMS变化都是最大的，ICEEMDAN方法和海潮负载改正结合中N、E、U 3个方向中RMS变化最大的有两个方向分别是N(BJSH站)和E(BJGB站)方向；2)变化最小的是海洋极潮改正，EMD方法和海洋极潮改正结合中N、E、U 3个方向中RMS变化最小的有两个方向分别是U(BJGB站)和E(BJSH站)方向，ICEEMDAN方法和海潮负载改正结合中N、E、U 3个方向中RMS变化最小的有两个分别是U(BJSH站)和E(BJGB站)站；剩下的两个环境负载改正中地球极移改正造成的RMS变化相较于大气压潮改正造成的RMS变化在U、E两个方向的变化量较大，在N方向的变化量较小. EMD方法和环境负载改正结合造成的RMS变化量比ICEEMDAN方法和环境负载结合造成的RMS变化量大，进一步说明了EMD方法分解得到的IMF分量中包含的噪声和非线性形变更多，因此环境负载改正造成的RMS变化量更大，ICEEMDAN方法的分解效果更好，IMF分量能更好的保留原始信息.

2.4 环境负载对GNSS坐标时间序列WRMS变化影响

对北京4个测站的坐标时间序列使用EMD、ICEEMDAN 2种方法进行分解并使用环境负载进行改正后，得到的WRMS变化如表2所示.

由表2可知，北京区域4个测站坐标时间序列环境负载改正的效果不尽相同，EMD方法和环境负载改正结合后在E方向的改正效果最好，WRMS变化全部为正，变化范围是0.037 1 mm (BJSH，海潮负载)～0.472 5 mm (BJYQ，海洋极潮)，说明在E方向经过环境负载改正后非线性形变得到了减弱，U方向和N方向经过环境负载改正后其中WRMS变化为正的都是9个，约占全部改正的56.25%，U方向的WRMS变化为正的范围是0.203 6 mm (BJYQ，海洋极潮)～1.535 6 mm (BJSH，海洋极潮)，WRMS变化为负的范围是-0.048 1 mm (BJYQ，地球极移)～ -1.279 5 mm(BJGB，地球极移)，可以看出改正后非线性形变增加的幅度相较于非线性形变减小的幅度较小；N方向的WRMS变化为正的范围是0.000 3 mm (BJGB，大气压潮)～0.085 1 mm (BJFS，大气压潮)，WRMS变化为负的范围是-0.001 0 mm (BJGB，地球极移)～ -0.128 3 mm(BJYQ，海潮负载)，N方向的整体环境负载改正WRMS变化最小，E方向整体环境负载改正WRMS变化次之，U方向的环境负载改正WRMS变化幅度最大.

表2 测站时间序列经环境负载改正后WRMS变化 mm

ICEEMDAN方法和环境负载改正结合后N方向的改正效果最好，WRMS变化全为正，变化范围是0.082 3 mm (BJGB，海潮负载)～0.230 4 mm (BJFS，大气压潮)，N方向的非线性形变全部得到了减弱，E方向的WRMS变化为正的有12个，约占全部改正75%，其中WRMS变化为正的范围是0.022 5 mm(BJFS，海洋极潮)～1.198 5 mm (BJSH，大气压潮)，WRMS变化为负的范围是-0.001 6 mm (BJGB，地球极移)～ -0.138 2 mm (BJFS，海潮负载)，非线性形变减弱的幅度大于非线性形变增加的幅度；U方向WRMS变化中为正的有10个，约占全部改正62.5%，其中WRMS变化为正的范围是0.078 8 mm (BJGB，海洋极潮)～1.960 6 mm (BJFS，大气压潮)，WRMS变化为负的范围是-0.087 4 mm (BJSH，海潮负载)～ -1.076 1 mm(BJSH，地球极移)；N方向的整体环境负载改正WRMS变化幅度最大，E方向的整体环境负载改正WRMS变化幅度次之，U方向的整体环境负载改正变化幅度最小. 为了更加明显的看出坐标时间序列经环境负载改正前后的变化，将坐标时间序列用图3表示出来，因为环境负载改正变化较小，若坐标轴时间较长则在图中看不出明显的变化，因此在2011—2021年的坐标时间序列变化中随机选取出2015年的坐标时间序列来对比，以BJYQ站为例，如图3所示.

由图3可知，在经过环境负载改正后坐标时间序列的峰值有了一定的变化，大多有减少的趋势，说明在环境负载改正后，坐标时间序列中的非线性形变得到了一定的减少.

图3 BJYQ站环境负载改正前后对比图

通过分别对比不同环境负载改正对坐标时间序列造成的WRMS变化量可以看出，海潮负载改正后WRMS变化为负的有8个，大气压潮改正后WRMS变化为负的有4个，地球极移改正后WRMS变化为负的有9个，海洋极潮改正后WRMS变化为负的有3个，说明海洋极潮改正对于去除非线性形变具有更好的效果，接下来是大气压潮改正，地球极移改正效果最不明显.

相较于EMD方法和环境负载改正结合后的WRMS变化，ICEEMDAN方法和环境负载改正结合后的WRMS变化中变化为正的数量明显较多，且变化幅度变大，说明ICEEMDAN方法和环境负载改正结合能够更加有效地减少测站时间序列中非线性形变的影响.

3 结 论

本文选取了北京区域的GNSS测站坐标时间序列，计算了测站对应的环境负载位移，分别分析了EMD和ICEEMDAN两种方法和环境负载改正结合后对测站坐标时间序列的影响，获得以下结论：

1)使用EMD和ICEEMDAN两种方法和环境负载改正结合后，各测站RMS变化各有区别，总体而言U方向变化最为明显，E方向变化次之，N方向变化最小；就4种环境负载改正造成的RMS变化来看，变化最大的是海潮负载，变化最小的是海洋极潮改正；而且EMD方法和环境负载改正结合造成的RMS变化量比ICEEMDAN方法和环境负载结合造成的RMS变化量大，说明EMD方法分解得到的IMF高频分量中包含的噪声和非线性形变更多，因此环境负载改正造成的RMS变化量更大，ICEEMDAN方法的分解效果更好，IMF分量能更好地保留原始信息.

2)环境负载改正对测站的3个方向的非线性形变大多具有一定的改善效果. EMD方法和环境负载改正结合后在E方向WRMS变化全部为正，U方向和N方向经过环境负载改正后其中WRMS变化为正的有9个，约占全部改正的56.25%；ICEEMDAN方法和环境负载改正结合后N方向WRMS变化全为正，E方向的WRMS变化为正的有12个，约占全部改正的75%，U方向的WRMS变化中为正的有10个，约占全部改正的62.5%. 对比4种环境负载改正造成的WRMS变化量可以看出海洋极潮改正对于去除非线性形变具有更好的效果，接下来是大气压潮改正，地球极移改正效果最不明显；ICEEMDAN方法和环境负载改正结合后的WRMS变化对比与EMD方法和环境负载改正结合变化为正的数量明显较多，且变化幅度变大，说明ICEEMDAN方法和环境负载改正结合能够更加有效的减少测站时间序列中非线性形变的含量.

致谢：感谢中国地震局第一监测中心提供论文研究的北京GNSS站点坐标时序.