基于ICEEMDAN與環(huán)境負(fù)載的GNSS坐標(biāo)時序非線性形變?nèi)コ?/h1>
2022-05-19 05:18:42王勇曹慧鵬李鎖閆勇楊軍
全球定位系統(tǒng) 2022年2期
王勇，曹慧鵬，李鎖，閆勇，楊軍
( 1. 天津城建大學(xué) 地質(zhì)與測繪學(xué)院, 天津 300384；2. 中國科學(xué)院精密測量科學(xué)與技術(shù)創(chuàng)新研究院 大地測量與地球動力學(xué)國家重點(diǎn)實(shí)驗室, 武漢430077；3. 天津市地質(zhì)工程勘測設(shè)計院有限公司, 天津 300191 )
0 引 言
全球衛(wèi)星導(dǎo)航系統(tǒng)(GNSS)連續(xù)運(yùn)行基準(zhǔn)站坐標(biāo)時序在導(dǎo)航與位置服務(wù)、精密衛(wèi)星定軌、地質(zhì)災(zāi)害監(jiān)測等工程和科學(xué)應(yīng)用中具有重要的作用[1-2]. GNSS坐標(biāo)時間序列可反映測站趨勢變化外，也能反映出測站非線性變化[3]. 在參考框架構(gòu)建時，非線性因素會對參考框架的精度產(chǎn)生一定的影響，研究基準(zhǔn)站的非線性形變?nèi)コ椒▽τ谛拚`差模型，提高GNSS基準(zhǔn)站位置精度，為研究地球物理、地球環(huán)境動態(tài)變化的檢測和研究提供依據(jù)具有重要的意義[4-6]. 姜衛(wèi)平等[7]研究發(fā)現(xiàn)中國區(qū)域基準(zhǔn)站位移呈現(xiàn)顯著的區(qū)域性特征，環(huán)境負(fù)載改正能夠削弱中國區(qū)域大多數(shù)GPS基準(zhǔn)站(約70%)垂向(U)以及東(E)方向位移時間序列的非線性形變，其計算環(huán)境負(fù)載改正由負(fù)載改正和GPS坐標(biāo)時間序列相加得到；陳岸等[8]研究了華北地區(qū)的GNSS基站站高程時間序列的非線性變化，研究表明GNSS時間序列存在地區(qū)差異性，且在高程方向非線性形變存在顯著的年周期和半年周期，水文負(fù)載對其的改正效果最好，其僅研究了高程方向的時間序列變化；尹濤濤等[9]使用了經(jīng)驗?zāi)B(tài)分解(EMD)的方法先將測站的坐標(biāo)時間序列進(jìn)行分解，提取出相應(yīng)的周期部分，然后再進(jìn)行環(huán)境負(fù)載的改正，但EMD方法進(jìn)行坐標(biāo)時間序列分解時會產(chǎn)生模態(tài)混疊現(xiàn)象，而且噪聲層中也會含有大量的有效信息，該方法無法對噪聲層中的有效信息進(jìn)行提取，這會造成一定的信息損失[10]. 因此有必要探尋一種新的方法進(jìn)行GNSS坐標(biāo)時序三維方向的非構(gòu)造形變?nèi)コ?，減少時序分解模態(tài)混疊現(xiàn)象，提高GNSS時序精度.
本文擬采用改進(jìn)的自適應(yīng)噪聲總體集合經(jīng)驗?zāi)B(tài)分解(ICEEMDAN)方法和環(huán)境負(fù)載改正相結(jié)合的方法開展GNSS測站非線性形變?nèi)コ芯? 首先使用GMIS軟件將GNSS坐標(biāo)時序補(bǔ)充完整并去除粗差，然后使用ICEEMDAN方法對GNSS坐標(biāo)時序進(jìn)行分解，使用排列熵算法選取包含噪聲和非線性形變的高頻分量，最后使用環(huán)境負(fù)載對高頻分量進(jìn)行去除，與EMD方法和環(huán)境負(fù)載結(jié)合的方法進(jìn)行去除效果對比.
1 研究數(shù)據(jù)與研究方法
1.1 GNSS坐標(biāo)時序
本文采用ICEEMDAN和環(huán)境負(fù)載相結(jié)合的方法開展GNSS測站非線性形變?nèi)コ婕癎NSS坐標(biāo)時序，采用中國大陸構(gòu)造環(huán)境監(jiān)測網(wǎng)絡(luò)(CMONOC)提供的北京4個站點(diǎn)坐標(biāo)時序開展相關(guān)研究. BJFS(北京房山)和BJSH(北京十三陵)建立于1995年與1998年，建站較早、觀測時間長，BJGB(北京古北口)、BJYQ(北京延慶) 2011年起存檔歷史觀測數(shù)據(jù)，本文研究GNSS站點(diǎn)坐標(biāo)時序時間跨度為2011-01-01—2021-05-04. 坐標(biāo)時序由中國地震局第一監(jiān)測中心提供(https://www.fmac.ac.cn).
1.2 環(huán)境負(fù)載數(shù)據(jù)
極潮、非潮汐海洋負(fù)載、大氣負(fù)載、水文負(fù)載等環(huán)境負(fù)載是造成測站非線性形變的主要原因[11-13]. 本次實(shí)驗采用大氣壓潮負(fù)載、海潮負(fù)載、海洋極潮、以及地球極移4種負(fù)載改正產(chǎn)品，由地球潮汐負(fù)荷影響與形變監(jiān)測計算系統(tǒng)ETideLoad4.0計算得到.
大氣壓潮負(fù)載數(shù)據(jù)采用歐洲中期氣候預(yù)報中心(https://www.ecmwf.int/) ECMWF-DCDA2006的大氣壓周日S1、半周日S2、半年Ssa和年周期Sa風(fēng)潮的0.5°×0.5°全球調(diào)和常數(shù)格網(wǎng)生成；海潮負(fù)載采用AVISO+(https://www.aviso.altimetry.fr/en/home.html)
的FES2014b潮高調(diào)和常數(shù)模型、經(jīng)ETideLosd4.0(分潮球諧分析與負(fù)荷球諧系數(shù)模型構(gòu)建)生成. 海洋極潮采用IERS2010協(xié)議(https://hpiers.obspm.fr/)標(biāo)準(zhǔn)自適應(yīng)海洋潮汐質(zhì)量均衡算法生成；地球極移采用IERS2010協(xié)議(https://hpiers.obspm.fr/)給出的地球定向參數(shù)產(chǎn)品EOPC04生成.
1.3 研究方法
由于EMD方法進(jìn)行坐標(biāo)時間序列分解時會產(chǎn)生模態(tài)混疊現(xiàn)象，本文采用ICEEMDAN方法進(jìn)行GNSS坐標(biāo)時序三維方向的非構(gòu)造形變?nèi)コ?ICEEMDAN對測站的坐標(biāo)時間序列進(jìn)行分解，重新提取本征模函數(shù)(IMF)，在增加高斯白噪聲的基礎(chǔ)上也會添加相應(yīng)特殊噪聲，也就是該白噪聲(WN)能夠被EMD分解之后獲取的k層IMF，隨后借助于獲取其中存在著特殊殘差(1個)，使得該IMF所對應(yīng)的定義就是殘差信號和有關(guān)局部均值所對應(yīng)的差值[14]. 結(jié)果顯示，IMF所對應(yīng)的剩余噪聲問題得到很好地改善，且基于EMD生成的IMF中的均值問題也能相應(yīng)地解決. 經(jīng)過ICEEMDAN方法分解得到的IMF分層能夠更好地體現(xiàn)原有坐標(biāo)時間序列，避免了有效信息的損失，并且解決了模態(tài)混疊的問題.
1.4 評價指標(biāo)
采用的評價指標(biāo)有兩種，第一個是均方根(RMS)作為分析環(huán)境負(fù)載對區(qū)域GNSS基準(zhǔn)站坐標(biāo)時間序列影響的一個指標(biāo)，其定義為[15]
式中：n表示觀測個數(shù)；m(j) 表示環(huán)境負(fù)載造成的位移.
第二個評價指標(biāo)采用加權(quán)均方根(WRMS)來進(jìn)行評價，使用負(fù)載改正前后GNSS時間序列的加權(quán)均方根差值來評估環(huán)境負(fù)載對GNSS坐標(biāo)時間序列的影響，其定義為[16]
式中：w(j)=1/sig(j)2；G(j)、sig(j) 分別表示基準(zhǔn)站在j時刻的位移(坐標(biāo)分量或位置相對于參考值的差值)及其精度(不確定度)；n表示觀測個數(shù).
計算基準(zhǔn)站 W RMS 變 化，WRMS變化=WRMS原始-WRMS改正， W RMS變化大于0表示環(huán)境負(fù)載改正后GNSS坐標(biāo)時間序列 W RMS 減小，即非線性形變變化幅度減弱.
2 GNSS坐標(biāo)時序非線性形變?nèi)コ?/h2>2.1 GNSS坐標(biāo)時序預(yù)處理
GNSS永久參考站自動運(yùn)行，因為通信鏈路終端或電源故障，GNSS坐標(biāo)時間序列中數(shù)據(jù)會經(jīng)常丟失，丟失的數(shù)據(jù)是隨機(jī)的，可能是3～5天，也可能是3～5個月，這將對數(shù)據(jù)分析造成很大的干擾，因此實(shí)驗先使用GMIS[17](GNSS Missing Data Interpolation)軟件對數(shù)據(jù)進(jìn)行修復(fù)，采用KKF(Kriged Kalman Filter)[18]的方法. 插值法又因為存在多路徑、軌道異常以及測站自身等因素會導(dǎo)致坐標(biāo)時間序列產(chǎn)生粗差，因此在使用KKF插值方法對數(shù)據(jù)進(jìn)行補(bǔ)充后，再使用3σ準(zhǔn)則對粗差進(jìn)行識別并剔除，以此來獲取一個相對“干凈且完整”的坐標(biāo)序列. 完善并剔除粗差后的測站坐標(biāo)時間序列結(jié)果如圖1所示(以BJFS站和BJSH站為例).
如圖1所示，黑線為原始的時間序列，可以看到存在缺失的部分，而且存在異常值，這些異常值就是由于各種原因造成的粗差，紅線為使用KKF法和3σ準(zhǔn)則進(jìn)行插值與去粗差后的結(jié)果，能明顯的看出缺失部分得到了補(bǔ)充，粗差得到了處理，而且補(bǔ)充后的數(shù)據(jù)與原始數(shù)據(jù)的趨勢基本一樣.
圖1 完整去噪后站點(diǎn)時間序列
2.2 基于ICEEMDAN方法的GNSS坐標(biāo)時序高頻信號提取
使用ICEEMDAN方法和EMD方法對測站的坐標(biāo)時間序列進(jìn)行分解，可獲得由高頻到低頻的IMF分量，為了更好地區(qū)分IMF分量中噪聲與非線性分量，使用了排列熵[19](PE)來衡量時間序列的復(fù)雜程度. PE具有運(yùn)算效率高、抗干擾能力強(qiáng)、輸出結(jié)果簡單的優(yōu)點(diǎn)，用于分析非線性、不規(guī)則的信號，對于區(qū)分IMF分量的頻率具有良好的適用性. 時間序列越規(guī)則，其值越??；時間序列越復(fù)雜，即包含的噪聲和非線性形變越多，其值越大. 本實(shí)驗參照參考文獻(xiàn)[20]將PE值的模糊隸屬度設(shè)定為0.6. 當(dāng)時間序列IMF分量大于0.6時，可認(rèn)為是包含噪聲和非線性形變的高頻信號；PE值小于0.6時，可看作干凈信號.計算每一個IMF分量的PE值，選出高頻信號然后和環(huán)境負(fù)載相結(jié)合來減少非線性形變. 得到的每個站點(diǎn)的PE值分布如圖2所示.
由圖2可知，每個測站的時間序列經(jīng)ICEEMDAN方法分解后在U方向的IMF分量都是11層，在水平N方向都是8層， E方向都是7層，而且3個方向計算得到的PE值大于0.6的分層都在前4層；經(jīng)EMD方法分解后在U方向的IMF分量基本都是10層，只有BJGB站的坐標(biāo)時間序列分解后IMF分量是9層，在N方向都是7層， E方向都是6層，而且3個方向計算得到的PE值大于0.6的分層都是在前3層；然后將兩種方法提取出來的高頻IMF分量進(jìn)行非線性形變?nèi)コ?
圖2 ICEEMDAN與EMD分解IMF分量PE值比較
2.3 環(huán)境負(fù)載對GNSS坐標(biāo)時間序列RMS影響
將4個測站坐標(biāo)時間序列使用EMD和ICEEMDAN方法進(jìn)行分解后，再使用排列熵篩選出包含噪聲和非線性形變的高頻層，然后分別用ETideLoad4.0計算得到的海潮負(fù)載、大氣壓潮負(fù)載、地球極移與海洋極潮等環(huán)境負(fù)載進(jìn)行改正，得到的RMS變化情況如表1所示.
由表1可知，環(huán)境負(fù)載改正對于坐標(biāo)時間序列的非線性形變具有一定的改正效果，在U方向變化最為明顯，E方向變化次之，N方向的變化量最小. EMD方法和環(huán)境負(fù)載結(jié)合改正后U方向的RMS變化范圍在5.037 0 mm (BJGB，海洋極潮)～7.663 1 mm (BJSH，海潮負(fù)載)，E方向的RMS變化次之，變化范圍在0.733 6 mm (BJSH，海洋極潮)～1.890 2 mm (BJGB，海潮負(fù)載)，N方向的RMS變化量最小，其范圍是0.869 9 mm (BJYQ，地球極移)～1.326 7 mm (BJFS，海潮負(fù)載)；ICEEMDAN方法和環(huán)境負(fù)載結(jié)合改正后U方向的RMS變化量在5.269 6 mm (BJSH，海洋極潮)～6.715 0 mm (BJSH，地球極移)，E方向的RMS變化量在0.817 3 mm (BJGB，海洋極潮)～1.697 2 mm(BJGB，海潮負(fù)載)，N方向的RMS變化量在0.816 2 mm(BJGB，地球極移)～1.513 0 mm (BJSH，海潮負(fù)載).
表1 測站時間序列經(jīng)環(huán)境負(fù)載改正后RMS變化 mm
通過分別對比不同環(huán)境負(fù)載改正對坐標(biāo)時間序列造成的RMS變化量可以看出： 1)變化最大的是海潮負(fù)載改正，EMD方法和海潮負(fù)載改正結(jié)合中N、E、U 3個方向RMS變化都是最大的，ICEEMDAN方法和海潮負(fù)載改正結(jié)合中N、E、U 3個方向中RMS變化最大的有兩個方向分別是N(BJSH站)和E(BJGB站)方向；2)變化最小的是海洋極潮改正，EMD方法和海洋極潮改正結(jié)合中N、E、U 3個方向中RMS變化最小的有兩個方向分別是U(BJGB站)和E(BJSH站)方向，ICEEMDAN方法和海潮負(fù)載改正結(jié)合中N、E、U 3個方向中RMS變化最小的有兩個分別是U(BJSH站)和E(BJGB站)站；剩下的兩個環(huán)境負(fù)載改正中地球極移改正造成的RMS變化相較于大氣壓潮改正造成的RMS變化在U、E兩個方向的變化量較大，在N方向的變化量較小. EMD方法和環(huán)境負(fù)載改正結(jié)合造成的RMS變化量比ICEEMDAN方法和環(huán)境負(fù)載結(jié)合造成的RMS變化量大，進(jìn)一步說明了EMD方法分解得到的IMF分量中包含的噪聲和非線性形變更多，因此環(huán)境負(fù)載改正造成的RMS變化量更大，ICEEMDAN方法的分解效果更好，IMF分量能更好的保留原始信息.
2.4 環(huán)境負(fù)載對GNSS坐標(biāo)時間序列WRMS變化影響
對北京4個測站的坐標(biāo)時間序列使用EMD、ICEEMDAN 2種方法進(jìn)行分解并使用環(huán)境負(fù)載進(jìn)行改正后，得到的WRMS變化如表2所示.
由表2可知，北京區(qū)域4個測站坐標(biāo)時間序列環(huán)境負(fù)載改正的效果不盡相同，EMD方法和環(huán)境負(fù)載改正結(jié)合后在E方向的改正效果最好，WRMS變化全部為正，變化范圍是0.037 1 mm (BJSH，海潮負(fù)載)～0.472 5 mm (BJYQ，海洋極潮)，說明在E方向經(jīng)過環(huán)境負(fù)載改正后非線性形變得到了減弱，U方向和N方向經(jīng)過環(huán)境負(fù)載改正后其中WRMS變化為正的都是9個，約占全部改正的56.25%，U方向的WRMS變化為正的范圍是0.203 6 mm (BJYQ，海洋極潮)～1.535 6 mm (BJSH，海洋極潮)，WRMS變化為負(fù)的范圍是-0.048 1 mm (BJYQ，地球極移)～ -1.279 5 mm(BJGB，地球極移)，可以看出改正后非線性形變增加的幅度相較于非線性形變減小的幅度較??；N方向的WRMS變化為正的范圍是0.000 3 mm (BJGB，大氣壓潮)～0.085 1 mm (BJFS，大氣壓潮)，WRMS變化為負(fù)的范圍是-0.001 0 mm (BJGB，地球極移)～ -0.128 3 mm(BJYQ，海潮負(fù)載)，N方向的整體環(huán)境負(fù)載改正WRMS變化最小，E方向整體環(huán)境負(fù)載改正WRMS變化次之，U方向的環(huán)境負(fù)載改正WRMS變化幅度最大.
表2 測站時間序列經(jīng)環(huán)境負(fù)載改正后WRMS變化 mm
ICEEMDAN方法和環(huán)境負(fù)載改正結(jié)合后N方向的改正效果最好，WRMS變化全為正，變化范圍是0.082 3 mm (BJGB，海潮負(fù)載)～0.230 4 mm (BJFS，大氣壓潮)，N方向的非線性形變?nèi)康玫搅藴p弱，E方向的WRMS變化為正的有12個，約占全部改正75%，其中WRMS變化為正的范圍是0.022 5 mm(BJFS，海洋極潮)～1.198 5 mm (BJSH，大氣壓潮)，WRMS變化為負(fù)的范圍是-0.001 6 mm (BJGB，地球極移)～ -0.138 2 mm (BJFS，海潮負(fù)載)，非線性形變減弱的幅度大于非線性形變增加的幅度；U方向WRMS變化中為正的有10個，約占全部改正62.5%，其中WRMS變化為正的范圍是0.078 8 mm (BJGB，海洋極潮)～1.960 6 mm (BJFS，大氣壓潮)，WRMS變化為負(fù)的范圍是-0.087 4 mm (BJSH，海潮負(fù)載)～ -1.076 1 mm(BJSH，地球極移)；N方向的整體環(huán)境負(fù)載改正WRMS變化幅度最大，E方向的整體環(huán)境負(fù)載改正WRMS變化幅度次之，U方向的整體環(huán)境負(fù)載改正變化幅度最小. 為了更加明顯的看出坐標(biāo)時間序列經(jīng)環(huán)境負(fù)載改正前后的變化，將坐標(biāo)時間序列用圖3表示出來，因為環(huán)境負(fù)載改正變化較小，若坐標(biāo)軸時間較長則在圖中看不出明顯的變化，因此在2011—2021年的坐標(biāo)時間序列變化中隨機(jī)選取出2015年的坐標(biāo)時間序列來對比，以BJYQ站為例，如圖3所示.
由圖3可知，在經(jīng)過環(huán)境負(fù)載改正后坐標(biāo)時間序列的峰值有了一定的變化，大多有減少的趨勢，說明在環(huán)境負(fù)載改正后，坐標(biāo)時間序列中的非線性形變得到了一定的減少.
圖3 BJYQ站環(huán)境負(fù)載改正前后對比圖
通過分別對比不同環(huán)境負(fù)載改正對坐標(biāo)時間序列造成的WRMS變化量可以看出，海潮負(fù)載改正后WRMS變化為負(fù)的有8個，大氣壓潮改正后WRMS變化為負(fù)的有4個，地球極移改正后WRMS變化為負(fù)的有9個，海洋極潮改正后WRMS變化為負(fù)的有3個，說明海洋極潮改正對于去除非線性形變具有更好的效果，接下來是大氣壓潮改正，地球極移改正效果最不明顯.
相較于EMD方法和環(huán)境負(fù)載改正結(jié)合后的WRMS變化，ICEEMDAN方法和環(huán)境負(fù)載改正結(jié)合后的WRMS變化中變化為正的數(shù)量明顯較多，且變化幅度變大，說明ICEEMDAN方法和環(huán)境負(fù)載改正結(jié)合能夠更加有效地減少測站時間序列中非線性形變的影響.
3 結(jié) 論
本文選取了北京區(qū)域的GNSS測站坐標(biāo)時間序列，計算了測站對應(yīng)的環(huán)境負(fù)載位移，分別分析了EMD和ICEEMDAN兩種方法和環(huán)境負(fù)載改正結(jié)合后對測站坐標(biāo)時間序列的影響，獲得以下結(jié)論：
1)使用EMD和ICEEMDAN兩種方法和環(huán)境負(fù)載改正結(jié)合后，各測站RMS變化各有區(qū)別，總體而言U方向變化最為明顯，E方向變化次之，N方向變化最??；就4種環(huán)境負(fù)載改正造成的RMS變化來看，變化最大的是海潮負(fù)載，變化最小的是海洋極潮改正；而且EMD方法和環(huán)境負(fù)載改正結(jié)合造成的RMS變化量比ICEEMDAN方法和環(huán)境負(fù)載結(jié)合造成的RMS變化量大，說明EMD方法分解得到的IMF高頻分量中包含的噪聲和非線性形變更多，因此環(huán)境負(fù)載改正造成的RMS變化量更大，ICEEMDAN方法的分解效果更好，IMF分量能更好地保留原始信息.
2)環(huán)境負(fù)載改正對測站的3個方向的非線性形變大多具有一定的改善效果. EMD方法和環(huán)境負(fù)載改正結(jié)合后在E方向WRMS變化全部為正，U方向和N方向經(jīng)過環(huán)境負(fù)載改正后其中WRMS變化為正的有9個，約占全部改正的56.25%；ICEEMDAN方法和環(huán)境負(fù)載改正結(jié)合后N方向WRMS變化全為正，E方向的WRMS變化為正的有12個，約占全部改正的75%，U方向的WRMS變化中為正的有10個，約占全部改正的62.5%. 對比4種環(huán)境負(fù)載改正造成的WRMS變化量可以看出海洋極潮改正對于去除非線性形變具有更好的效果，接下來是大氣壓潮改正，地球極移改正效果最不明顯；ICEEMDAN方法和環(huán)境負(fù)載改正結(jié)合后的WRMS變化對比與EMD方法和環(huán)境負(fù)載改正結(jié)合變化為正的數(shù)量明顯較多，且變化幅度變大，說明ICEEMDAN方法和環(huán)境負(fù)載改正結(jié)合能夠更加有效的減少測站時間序列中非線性形變的含量.
致謝：感謝中國地震局第一監(jiān)測中心提供論文研究的北京GNSS站點(diǎn)坐標(biāo)時序.
猜你喜歡
時序坐標(biāo)《學(xué)習(xí)方法報》歷史中考版(2024年8期)2024-12-31 00:00:00 GNSS鐘差估計中的兩種測站選取策略分析地理空間信息(2022年11期)2022-11-26 05:47:02 柳永·望海潮韶關(guān)學(xué)院學(xué)報(2022年4期)2022-04-11 11:12:16 基于Sentinel-2時序NDVI的麥冬識別研究中國農(nóng)業(yè)信息(2021年3期)2021-11-22 06:44:48 全球GPS測站垂向周年變化統(tǒng)計改正模型的建立測繪學(xué)報(2018年10期)2018-10-26 06:12:16 測站分布對GPS解算ERP的影響分析城市勘測(2018年1期)2018-03-15 03:37:02 望海潮·八里湖老友(2017年7期)2017-08-22 02:36:50 一種毫米波放大器時序直流電源的設(shè)計電子制作(2016年15期)2017-01-15 13:39:08 望海潮·楓林谷新農(nóng)業(yè)(2016年23期)2016-08-16 03:38:17 望海潮·慶嫦三落月太空探索(2014年2期)2014-07-12 09:52:15
王勇，曹慧鵬，李鎖，閆勇，楊軍
( 1. 天津城建大學(xué) 地質(zhì)與測繪學(xué)院, 天津 300384；2. 中國科學(xué)院精密測量科學(xué)與技術(shù)創(chuàng)新研究院 大地測量與地球動力學(xué)國家重點(diǎn)實(shí)驗室, 武漢430077；3. 天津市地質(zhì)工程勘測設(shè)計院有限公司, 天津 300191 )
0 引 言
全球衛(wèi)星導(dǎo)航系統(tǒng)(GNSS)連續(xù)運(yùn)行基準(zhǔn)站坐標(biāo)時序在導(dǎo)航與位置服務(wù)、精密衛(wèi)星定軌、地質(zhì)災(zāi)害監(jiān)測等工程和科學(xué)應(yīng)用中具有重要的作用[1-2]. GNSS坐標(biāo)時間序列可反映測站趨勢變化外，也能反映出測站非線性變化[3]. 在參考框架構(gòu)建時，非線性因素會對參考框架的精度產(chǎn)生一定的影響，研究基準(zhǔn)站的非線性形變?nèi)コ椒▽τ谛拚`差模型，提高GNSS基準(zhǔn)站位置精度，為研究地球物理、地球環(huán)境動態(tài)變化的檢測和研究提供依據(jù)具有重要的意義[4-6]. 姜衛(wèi)平等[7]研究發(fā)現(xiàn)中國區(qū)域基準(zhǔn)站位移呈現(xiàn)顯著的區(qū)域性特征，環(huán)境負(fù)載改正能夠削弱中國區(qū)域大多數(shù)GPS基準(zhǔn)站(約70%)垂向(U)以及東(E)方向位移時間序列的非線性形變，其計算環(huán)境負(fù)載改正由負(fù)載改正和GPS坐標(biāo)時間序列相加得到；陳岸等[8]研究了華北地區(qū)的GNSS基站站高程時間序列的非線性變化，研究表明GNSS時間序列存在地區(qū)差異性，且在高程方向非線性形變存在顯著的年周期和半年周期，水文負(fù)載對其的改正效果最好，其僅研究了高程方向的時間序列變化；尹濤濤等[9]使用了經(jīng)驗?zāi)B(tài)分解(EMD)的方法先將測站的坐標(biāo)時間序列進(jìn)行分解，提取出相應(yīng)的周期部分，然后再進(jìn)行環(huán)境負(fù)載的改正，但EMD方法進(jìn)行坐標(biāo)時間序列分解時會產(chǎn)生模態(tài)混疊現(xiàn)象，而且噪聲層中也會含有大量的有效信息，該方法無法對噪聲層中的有效信息進(jìn)行提取，這會造成一定的信息損失[10]. 因此有必要探尋一種新的方法進(jìn)行GNSS坐標(biāo)時序三維方向的非構(gòu)造形變?nèi)コ?，減少時序分解模態(tài)混疊現(xiàn)象，提高GNSS時序精度.
本文擬采用改進(jìn)的自適應(yīng)噪聲總體集合經(jīng)驗?zāi)B(tài)分解(ICEEMDAN)方法和環(huán)境負(fù)載改正相結(jié)合的方法開展GNSS測站非線性形變?nèi)コ芯? 首先使用GMIS軟件將GNSS坐標(biāo)時序補(bǔ)充完整并去除粗差，然后使用ICEEMDAN方法對GNSS坐標(biāo)時序進(jìn)行分解，使用排列熵算法選取包含噪聲和非線性形變的高頻分量，最后使用環(huán)境負(fù)載對高頻分量進(jìn)行去除，與EMD方法和環(huán)境負(fù)載結(jié)合的方法進(jìn)行去除效果對比.
1 研究數(shù)據(jù)與研究方法
1.1 GNSS坐標(biāo)時序
本文采用ICEEMDAN和環(huán)境負(fù)載相結(jié)合的方法開展GNSS測站非線性形變?nèi)コ婕癎NSS坐標(biāo)時序，采用中國大陸構(gòu)造環(huán)境監(jiān)測網(wǎng)絡(luò)(CMONOC)提供的北京4個站點(diǎn)坐標(biāo)時序開展相關(guān)研究. BJFS(北京房山)和BJSH(北京十三陵)建立于1995年與1998年，建站較早、觀測時間長，BJGB(北京古北口)、BJYQ(北京延慶) 2011年起存檔歷史觀測數(shù)據(jù)，本文研究GNSS站點(diǎn)坐標(biāo)時序時間跨度為2011-01-01—2021-05-04. 坐標(biāo)時序由中國地震局第一監(jiān)測中心提供(https://www.fmac.ac.cn).
1.2 環(huán)境負(fù)載數(shù)據(jù)
極潮、非潮汐海洋負(fù)載、大氣負(fù)載、水文負(fù)載等環(huán)境負(fù)載是造成測站非線性形變的主要原因[11-13]. 本次實(shí)驗采用大氣壓潮負(fù)載、海潮負(fù)載、海洋極潮、以及地球極移4種負(fù)載改正產(chǎn)品，由地球潮汐負(fù)荷影響與形變監(jiān)測計算系統(tǒng)ETideLoad4.0計算得到.
大氣壓潮負(fù)載數(shù)據(jù)采用歐洲中期氣候預(yù)報中心(https://www.ecmwf.int/) ECMWF-DCDA2006的大氣壓周日S1、半周日S2、半年Ssa和年周期Sa風(fēng)潮的0.5°×0.5°全球調(diào)和常數(shù)格網(wǎng)生成；海潮負(fù)載采用AVISO+(https://www.aviso.altimetry.fr/en/home.html)
的FES2014b潮高調(diào)和常數(shù)模型、經(jīng)ETideLosd4.0(分潮球諧分析與負(fù)荷球諧系數(shù)模型構(gòu)建)生成. 海洋極潮采用IERS2010協(xié)議(https://hpiers.obspm.fr/)標(biāo)準(zhǔn)自適應(yīng)海洋潮汐質(zhì)量均衡算法生成；地球極移采用IERS2010協(xié)議(https://hpiers.obspm.fr/)給出的地球定向參數(shù)產(chǎn)品EOPC04生成.
1.3 研究方法
由于EMD方法進(jìn)行坐標(biāo)時間序列分解時會產(chǎn)生模態(tài)混疊現(xiàn)象，本文采用ICEEMDAN方法進(jìn)行GNSS坐標(biāo)時序三維方向的非構(gòu)造形變?nèi)コ?ICEEMDAN對測站的坐標(biāo)時間序列進(jìn)行分解，重新提取本征模函數(shù)(IMF)，在增加高斯白噪聲的基礎(chǔ)上也會添加相應(yīng)特殊噪聲，也就是該白噪聲(WN)能夠被EMD分解之后獲取的k層IMF，隨后借助于獲取其中存在著特殊殘差(1個)，使得該IMF所對應(yīng)的定義就是殘差信號和有關(guān)局部均值所對應(yīng)的差值[14]. 結(jié)果顯示，IMF所對應(yīng)的剩余噪聲問題得到很好地改善，且基于EMD生成的IMF中的均值問題也能相應(yīng)地解決. 經(jīng)過ICEEMDAN方法分解得到的IMF分層能夠更好地體現(xiàn)原有坐標(biāo)時間序列，避免了有效信息的損失，并且解決了模態(tài)混疊的問題.
1.4 評價指標(biāo)
采用的評價指標(biāo)有兩種，第一個是均方根(RMS)作為分析環(huán)境負(fù)載對區(qū)域GNSS基準(zhǔn)站坐標(biāo)時間序列影響的一個指標(biāo)，其定義為[15]
式中：n表示觀測個數(shù)；m(j) 表示環(huán)境負(fù)載造成的位移.
第二個評價指標(biāo)采用加權(quán)均方根(WRMS)來進(jìn)行評價，使用負(fù)載改正前后GNSS時間序列的加權(quán)均方根差值來評估環(huán)境負(fù)載對GNSS坐標(biāo)時間序列的影響，其定義為[16]
式中：w(j)=1/sig(j)2；G(j)、sig(j) 分別表示基準(zhǔn)站在j時刻的位移(坐標(biāo)分量或位置相對于參考值的差值)及其精度(不確定度)；n表示觀測個數(shù).
計算基準(zhǔn)站 W RMS 變 化，WRMS變化=WRMS原始-WRMS改正， W RMS變化大于0表示環(huán)境負(fù)載改正后GNSS坐標(biāo)時間序列 W RMS 減小，即非線性形變變化幅度減弱.
2 GNSS坐標(biāo)時序非線性形變?nèi)コ?/h2>2.1 GNSS坐標(biāo)時序預(yù)處理
GNSS永久參考站自動運(yùn)行，因為通信鏈路終端或電源故障，GNSS坐標(biāo)時間序列中數(shù)據(jù)會經(jīng)常丟失，丟失的數(shù)據(jù)是隨機(jī)的，可能是3～5天，也可能是3～5個月，這將對數(shù)據(jù)分析造成很大的干擾，因此實(shí)驗先使用GMIS[17](GNSS Missing Data Interpolation)軟件對數(shù)據(jù)進(jìn)行修復(fù)，采用KKF(Kriged Kalman Filter)[18]的方法. 插值法又因為存在多路徑、軌道異常以及測站自身等因素會導(dǎo)致坐標(biāo)時間序列產(chǎn)生粗差，因此在使用KKF插值方法對數(shù)據(jù)進(jìn)行補(bǔ)充后，再使用3σ準(zhǔn)則對粗差進(jìn)行識別并剔除，以此來獲取一個相對“干凈且完整”的坐標(biāo)序列. 完善并剔除粗差后的測站坐標(biāo)時間序列結(jié)果如圖1所示(以BJFS站和BJSH站為例).
如圖1所示，黑線為原始的時間序列，可以看到存在缺失的部分，而且存在異常值，這些異常值就是由于各種原因造成的粗差，紅線為使用KKF法和3σ準(zhǔn)則進(jìn)行插值與去粗差后的結(jié)果，能明顯的看出缺失部分得到了補(bǔ)充，粗差得到了處理，而且補(bǔ)充后的數(shù)據(jù)與原始數(shù)據(jù)的趨勢基本一樣.
圖1 完整去噪后站點(diǎn)時間序列
2.2 基于ICEEMDAN方法的GNSS坐標(biāo)時序高頻信號提取
使用ICEEMDAN方法和EMD方法對測站的坐標(biāo)時間序列進(jìn)行分解，可獲得由高頻到低頻的IMF分量，為了更好地區(qū)分IMF分量中噪聲與非線性分量，使用了排列熵[19](PE)來衡量時間序列的復(fù)雜程度. PE具有運(yùn)算效率高、抗干擾能力強(qiáng)、輸出結(jié)果簡單的優(yōu)點(diǎn)，用于分析非線性、不規(guī)則的信號，對于區(qū)分IMF分量的頻率具有良好的適用性. 時間序列越規(guī)則，其值越??；時間序列越復(fù)雜，即包含的噪聲和非線性形變越多，其值越大. 本實(shí)驗參照參考文獻(xiàn)[20]將PE值的模糊隸屬度設(shè)定為0.6. 當(dāng)時間序列IMF分量大于0.6時，可認(rèn)為是包含噪聲和非線性形變的高頻信號；PE值小于0.6時，可看作干凈信號.計算每一個IMF分量的PE值，選出高頻信號然后和環(huán)境負(fù)載相結(jié)合來減少非線性形變. 得到的每個站點(diǎn)的PE值分布如圖2所示.
由圖2可知，每個測站的時間序列經(jīng)ICEEMDAN方法分解后在U方向的IMF分量都是11層，在水平N方向都是8層， E方向都是7層，而且3個方向計算得到的PE值大于0.6的分層都在前4層；經(jīng)EMD方法分解后在U方向的IMF分量基本都是10層，只有BJGB站的坐標(biāo)時間序列分解后IMF分量是9層，在N方向都是7層， E方向都是6層，而且3個方向計算得到的PE值大于0.6的分層都是在前3層；然后將兩種方法提取出來的高頻IMF分量進(jìn)行非線性形變?nèi)コ?
圖2 ICEEMDAN與EMD分解IMF分量PE值比較
2.3 環(huán)境負(fù)載對GNSS坐標(biāo)時間序列RMS影響
將4個測站坐標(biāo)時間序列使用EMD和ICEEMDAN方法進(jìn)行分解后，再使用排列熵篩選出包含噪聲和非線性形變的高頻層，然后分別用ETideLoad4.0計算得到的海潮負(fù)載、大氣壓潮負(fù)載、地球極移與海洋極潮等環(huán)境負(fù)載進(jìn)行改正，得到的RMS變化情況如表1所示.
由表1可知，環(huán)境負(fù)載改正對于坐標(biāo)時間序列的非線性形變具有一定的改正效果，在U方向變化最為明顯，E方向變化次之，N方向的變化量最小. EMD方法和環(huán)境負(fù)載結(jié)合改正后U方向的RMS變化范圍在5.037 0 mm (BJGB，海洋極潮)～7.663 1 mm (BJSH，海潮負(fù)載)，E方向的RMS變化次之，變化范圍在0.733 6 mm (BJSH，海洋極潮)～1.890 2 mm (BJGB，海潮負(fù)載)，N方向的RMS變化量最小，其范圍是0.869 9 mm (BJYQ，地球極移)～1.326 7 mm (BJFS，海潮負(fù)載)；ICEEMDAN方法和環(huán)境負(fù)載結(jié)合改正后U方向的RMS變化量在5.269 6 mm (BJSH，海洋極潮)～6.715 0 mm (BJSH，地球極移)，E方向的RMS變化量在0.817 3 mm (BJGB，海洋極潮)～1.697 2 mm(BJGB，海潮負(fù)載)，N方向的RMS變化量在0.816 2 mm(BJGB，地球極移)～1.513 0 mm (BJSH，海潮負(fù)載).
表1 測站時間序列經(jīng)環(huán)境負(fù)載改正后RMS變化 mm
通過分別對比不同環(huán)境負(fù)載改正對坐標(biāo)時間序列造成的RMS變化量可以看出： 1)變化最大的是海潮負(fù)載改正，EMD方法和海潮負(fù)載改正結(jié)合中N、E、U 3個方向RMS變化都是最大的，ICEEMDAN方法和海潮負(fù)載改正結(jié)合中N、E、U 3個方向中RMS變化最大的有兩個方向分別是N(BJSH站)和E(BJGB站)方向；2)變化最小的是海洋極潮改正，EMD方法和海洋極潮改正結(jié)合中N、E、U 3個方向中RMS變化最小的有兩個方向分別是U(BJGB站)和E(BJSH站)方向，ICEEMDAN方法和海潮負(fù)載改正結(jié)合中N、E、U 3個方向中RMS變化最小的有兩個分別是U(BJSH站)和E(BJGB站)站；剩下的兩個環(huán)境負(fù)載改正中地球極移改正造成的RMS變化相較于大氣壓潮改正造成的RMS變化在U、E兩個方向的變化量較大，在N方向的變化量較小. EMD方法和環(huán)境負(fù)載改正結(jié)合造成的RMS變化量比ICEEMDAN方法和環(huán)境負(fù)載結(jié)合造成的RMS變化量大，進(jìn)一步說明了EMD方法分解得到的IMF分量中包含的噪聲和非線性形變更多，因此環(huán)境負(fù)載改正造成的RMS變化量更大，ICEEMDAN方法的分解效果更好，IMF分量能更好的保留原始信息.
2.4 環(huán)境負(fù)載對GNSS坐標(biāo)時間序列WRMS變化影響
對北京4個測站的坐標(biāo)時間序列使用EMD、ICEEMDAN 2種方法進(jìn)行分解并使用環(huán)境負(fù)載進(jìn)行改正后，得到的WRMS變化如表2所示.
由表2可知，北京區(qū)域4個測站坐標(biāo)時間序列環(huán)境負(fù)載改正的效果不盡相同，EMD方法和環(huán)境負(fù)載改正結(jié)合后在E方向的改正效果最好，WRMS變化全部為正，變化范圍是0.037 1 mm (BJSH，海潮負(fù)載)～0.472 5 mm (BJYQ，海洋極潮)，說明在E方向經(jīng)過環(huán)境負(fù)載改正后非線性形變得到了減弱，U方向和N方向經(jīng)過環(huán)境負(fù)載改正后其中WRMS變化為正的都是9個，約占全部改正的56.25%，U方向的WRMS變化為正的范圍是0.203 6 mm (BJYQ，海洋極潮)～1.535 6 mm (BJSH，海洋極潮)，WRMS變化為負(fù)的范圍是-0.048 1 mm (BJYQ，地球極移)～ -1.279 5 mm(BJGB，地球極移)，可以看出改正后非線性形變增加的幅度相較于非線性形變減小的幅度較??；N方向的WRMS變化為正的范圍是0.000 3 mm (BJGB，大氣壓潮)～0.085 1 mm (BJFS，大氣壓潮)，WRMS變化為負(fù)的范圍是-0.001 0 mm (BJGB，地球極移)～ -0.128 3 mm(BJYQ，海潮負(fù)載)，N方向的整體環(huán)境負(fù)載改正WRMS變化最小，E方向整體環(huán)境負(fù)載改正WRMS變化次之，U方向的環(huán)境負(fù)載改正WRMS變化幅度最大.
表2 測站時間序列經(jīng)環(huán)境負(fù)載改正后WRMS變化 mm
ICEEMDAN方法和環(huán)境負(fù)載改正結(jié)合后N方向的改正效果最好，WRMS變化全為正，變化范圍是0.082 3 mm (BJGB，海潮負(fù)載)～0.230 4 mm (BJFS，大氣壓潮)，N方向的非線性形變?nèi)康玫搅藴p弱，E方向的WRMS變化為正的有12個，約占全部改正75%，其中WRMS變化為正的范圍是0.022 5 mm(BJFS，海洋極潮)～1.198 5 mm (BJSH，大氣壓潮)，WRMS變化為負(fù)的范圍是-0.001 6 mm (BJGB，地球極移)～ -0.138 2 mm (BJFS，海潮負(fù)載)，非線性形變減弱的幅度大于非線性形變增加的幅度；U方向WRMS變化中為正的有10個，約占全部改正62.5%，其中WRMS變化為正的范圍是0.078 8 mm (BJGB，海洋極潮)～1.960 6 mm (BJFS，大氣壓潮)，WRMS變化為負(fù)的范圍是-0.087 4 mm (BJSH，海潮負(fù)載)～ -1.076 1 mm(BJSH，地球極移)；N方向的整體環(huán)境負(fù)載改正WRMS變化幅度最大，E方向的整體環(huán)境負(fù)載改正WRMS變化幅度次之，U方向的整體環(huán)境負(fù)載改正變化幅度最小. 為了更加明顯的看出坐標(biāo)時間序列經(jīng)環(huán)境負(fù)載改正前后的變化，將坐標(biāo)時間序列用圖3表示出來，因為環(huán)境負(fù)載改正變化較小，若坐標(biāo)軸時間較長則在圖中看不出明顯的變化，因此在2011—2021年的坐標(biāo)時間序列變化中隨機(jī)選取出2015年的坐標(biāo)時間序列來對比，以BJYQ站為例，如圖3所示.
由圖3可知，在經(jīng)過環(huán)境負(fù)載改正后坐標(biāo)時間序列的峰值有了一定的變化，大多有減少的趨勢，說明在環(huán)境負(fù)載改正后，坐標(biāo)時間序列中的非線性形變得到了一定的減少.
圖3 BJYQ站環(huán)境負(fù)載改正前后對比圖
通過分別對比不同環(huán)境負(fù)載改正對坐標(biāo)時間序列造成的WRMS變化量可以看出，海潮負(fù)載改正后WRMS變化為負(fù)的有8個，大氣壓潮改正后WRMS變化為負(fù)的有4個，地球極移改正后WRMS變化為負(fù)的有9個，海洋極潮改正后WRMS變化為負(fù)的有3個，說明海洋極潮改正對于去除非線性形變具有更好的效果，接下來是大氣壓潮改正，地球極移改正效果最不明顯.
相較于EMD方法和環(huán)境負(fù)載改正結(jié)合后的WRMS變化，ICEEMDAN方法和環(huán)境負(fù)載改正結(jié)合后的WRMS變化中變化為正的數(shù)量明顯較多，且變化幅度變大，說明ICEEMDAN方法和環(huán)境負(fù)載改正結(jié)合能夠更加有效地減少測站時間序列中非線性形變的影響.
3 結(jié) 論
本文選取了北京區(qū)域的GNSS測站坐標(biāo)時間序列，計算了測站對應(yīng)的環(huán)境負(fù)載位移，分別分析了EMD和ICEEMDAN兩種方法和環(huán)境負(fù)載改正結(jié)合后對測站坐標(biāo)時間序列的影響，獲得以下結(jié)論：
1)使用EMD和ICEEMDAN兩種方法和環(huán)境負(fù)載改正結(jié)合后，各測站RMS變化各有區(qū)別，總體而言U方向變化最為明顯，E方向變化次之，N方向變化最??；就4種環(huán)境負(fù)載改正造成的RMS變化來看，變化最大的是海潮負(fù)載，變化最小的是海洋極潮改正；而且EMD方法和環(huán)境負(fù)載改正結(jié)合造成的RMS變化量比ICEEMDAN方法和環(huán)境負(fù)載結(jié)合造成的RMS變化量大，說明EMD方法分解得到的IMF高頻分量中包含的噪聲和非線性形變更多，因此環(huán)境負(fù)載改正造成的RMS變化量更大，ICEEMDAN方法的分解效果更好，IMF分量能更好地保留原始信息.
2)環(huán)境負(fù)載改正對測站的3個方向的非線性形變大多具有一定的改善效果. EMD方法和環(huán)境負(fù)載改正結(jié)合后在E方向WRMS變化全部為正，U方向和N方向經(jīng)過環(huán)境負(fù)載改正后其中WRMS變化為正的有9個，約占全部改正的56.25%；ICEEMDAN方法和環(huán)境負(fù)載改正結(jié)合后N方向WRMS變化全為正，E方向的WRMS變化為正的有12個，約占全部改正的75%，U方向的WRMS變化中為正的有10個，約占全部改正的62.5%. 對比4種環(huán)境負(fù)載改正造成的WRMS變化量可以看出海洋極潮改正對于去除非線性形變具有更好的效果，接下來是大氣壓潮改正，地球極移改正效果最不明顯；ICEEMDAN方法和環(huán)境負(fù)載改正結(jié)合后的WRMS變化對比與EMD方法和環(huán)境負(fù)載改正結(jié)合變化為正的數(shù)量明顯較多，且變化幅度變大，說明ICEEMDAN方法和環(huán)境負(fù)載改正結(jié)合能夠更加有效的減少測站時間序列中非線性形變的含量.
致謝：感謝中國地震局第一監(jiān)測中心提供論文研究的北京GNSS站點(diǎn)坐標(biāo)時序.
猜你喜歡
《學(xué)習(xí)方法報》歷史中考版(2024年8期)2024-12-31 00:00:00
地理空間信息(2022年11期)2022-11-26 05:47:02
韶關(guān)學(xué)院學(xué)報(2022年4期)2022-04-11 11:12:16
中國農(nóng)業(yè)信息(2021年3期)2021-11-22 06:44:48
測繪學(xué)報(2018年10期)2018-10-26 06:12:16
城市勘測(2018年1期)2018-03-15 03:37:02
老友(2017年7期)2017-08-22 02:36:50
電子制作(2016年15期)2017-01-15 13:39:08
新農(nóng)業(yè)(2016年23期)2016-08-16 03:38:17
太空探索(2014年2期)2014-07-12 09:52:15
