基于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站點坐標時序.