內(nèi)蒙古陸態(tài)網(wǎng)絡坐標時間序列噪聲分析及速度場修正

李鑫，姚德新

（ 1. 蘭州交通大學 測繪與地理信息學院, 蘭州 730070；2. 地理國情監(jiān)測技術(shù)應用國家地方聯(lián)合工程研究中心, 蘭州 730070；3. 甘肅省地理國情監(jiān)測工程實驗室, 蘭州 730070 ）

0 引 言

中國大陸構(gòu)造環(huán)境監(jiān)測網(wǎng)絡(CMONOC)主要用于監(jiān)測中國大陸地殼運動、維護坐標參考框架[1]，對高精度坐標基準的建立和維持等科學研究提供基礎資料和產(chǎn)品. 截至目前，陸態(tài)網(wǎng)絡基準站已經(jīng)建立了260個連續(xù)觀測站[2]，其中內(nèi)蒙古已經(jīng)建立15座陸態(tài)網(wǎng)絡基準站，這些站點構(gòu)成了內(nèi)蒙古地區(qū)的三維(3D)坐標框架，對于內(nèi)蒙古地區(qū)板塊運動研究和高等級控制網(wǎng)的維護、更新起到了決定性的作用. 自陸態(tài)網(wǎng)絡建立以來，已經(jīng)積攢了大量的觀測數(shù)據(jù)，可用作區(qū)域速度場的維護和地殼運動的監(jiān)測[3]，由文獻[4]可知，利用這些觀測數(shù)據(jù)對某一區(qū)域的速度場進行可靠性估計，需要分析該區(qū)域陸態(tài)網(wǎng)絡基準站至少連續(xù)5年的時間序列觀測數(shù)據(jù). 本文的研究區(qū)域中，只有NMBT站連續(xù)觀測時間為6年，其他站點連續(xù)觀測時間近10年，時間序列觀測數(shù)據(jù)完全滿足要求. 張錫越等[2]對北京地區(qū)13個連續(xù)運行參考站系統(tǒng)(CORS)的時間序列觀測數(shù)據(jù)進行處理，得到各站點三方向的最優(yōu)噪聲模型組合白噪聲+閃爍噪聲+隨機漫步噪聲(WN+FN+RWN)，最終得到北京地區(qū)速度場估計. 分析內(nèi)蒙古站點時間序列觀測數(shù)據(jù)中，發(fā)現(xiàn)主要存在四種噪聲，分別是：FN、冪律噪聲(PL)、WN和RWN.對上述四種噪聲在北(N)、東(E)、天頂(U)三方向選取最佳噪聲模型后，可對內(nèi)蒙古地區(qū)站點進行速度場估計. 內(nèi)蒙古區(qū)域遼闊，東西跨度大，地質(zhì)變化跨度也大，因此需要在N、E、U三個方向做完整的速度場估計.

1 基準站坐標時間序列及其精度

本文選取整個內(nèi)蒙古地區(qū)的陸態(tài)網(wǎng)絡連續(xù)站，包括NMAG、NMAL、NMAY、NMAZ、NMBT、NMDW、NMEJ、NMEL、NMER、NMTK、NMWH、NMWJ、NMWL、NMWT、NMZL共15個基準站，觀測數(shù)據(jù)時間段為：2010-06—2021-06，各站點具體時間段如表1所示，圖1為站點分布情況. 在ITRF14框架下，使用GAMIT/GLOBK軟件解算[5-7]，得到內(nèi)蒙古自治區(qū)15個陸態(tài)網(wǎng)絡連續(xù)站近10年的坐標時間序列觀測數(shù)據(jù)，具體解算流程參考文獻[8].

圖1 內(nèi)蒙古連續(xù)運行的15個陸態(tài)網(wǎng)絡基準站分布

表1 內(nèi)蒙古陸態(tài)網(wǎng)絡各站點連續(xù)觀測時間段

加權(quán)均方根誤差(WRMS)是評價GNSS時間序列精度常用指標，由文獻[9]可知，對應站點N、E、U方向的WRMS值越小說明時間序列離散程度越小，重復性越好，精度越高. 計算公式為

式中：nexp是預測值；ncale是實際觀測值； σexp是預測值中誤差. 根據(jù)文獻[10]可知，理想情況下，N、E方向重復性為1～2 mm，U方向重復性為3～10 mm. 內(nèi)蒙地區(qū)各個站點的N、E、U方向的WRMS值如表2所示. 由表2可知，僅NMWT的N、E方向WRMS值較大，其余站點原始時間序列均滿足理想情況.

表2 內(nèi)蒙地區(qū)各站點時間序列的N、E、U方向的WRMS值

2 坐標時間序列噪聲分析

2.1 坐標時間序列預處理

原始坐標時間序列可表示為[11]

式中：y(ti) 為 原始坐標時間序列；ti是某時刻的歷元，單位為a；參數(shù)a表示各陸態(tài)網(wǎng)絡站點初始坐標；b表示線性變化速率；c、d，e、f分別表示年周期和半年周期運動參數(shù)；表示外界環(huán)境因素造成的階遷偏移，其中H(t) 表示階躍函數(shù)、參數(shù)Tgi指代同震偏移；vi是初始觀測值與預期結(jié)果不符所產(chǎn)生的偏差，即殘差[12].

由式(2)可知，將原始坐標時間序列作為觀測值，基于最小二乘法進行平差計算，求解各項參數(shù)，得到坐標殘差時間序列.利用HECTOR軟件(研發(fā)的一款基于C++語言的時間序列分析軟件[8])剔除時間序列奇異值，首先要去除線性項和周期項計算殘差，然后利用四分位數(shù)粗差探測方法探測原始觀測數(shù)據(jù)所包含的粗差，并剔除其中的奇異值. 以NMEL站為例，圖2(a)、(b)、(c)分別為此站點在N、E、U三方向上處理前的時間序列圖像，圖3(a)、(b)、(c)分別為此站點在N、E、U三方向上處理后的時間序列圖像.

圖2 NMEL站N、E、U方向原始時間序列

圖3 NMEL站N、E、U方向處理后的時間序列

2.2 坐標時間序列噪聲模型分析

HECTOR時間序列分析軟件基于貝葉斯數(shù)值分析法，該方法是最大似然估計法(MLE)的優(yōu)化[13]. 此軟件相較目前使用的CATS等時間序列分析軟件，在處理大量數(shù)據(jù)時，解算速度有大幅度提升，且此軟件可分析WN、PL和分形自回歸聚合滑動模型(ARFIMA)噪聲等，還可以將單個噪聲模型任意組合形成新的組合噪聲模型[13]. 使用貝葉斯數(shù)值分析法對不同噪聲模型或者不同組合噪聲模型分析的函數(shù)式為

式中：N為實際觀測數(shù)據(jù)的數(shù)量(不包括缺失的數(shù)據(jù))；r為殘差矩陣. 協(xié)方差陣C可分解為

式中：C為各類噪聲的總和；σ是殘差序列的中誤差.

因detCA=CNdetA，結(jié)合式(4)、(5)，式(3)可寫為

于是，BIC值的具體計算為

式中，k代表模型所需估計的參數(shù)個數(shù)[13].

通過上式計算得到BIC值，如果BIC值越小，則說明噪聲模型選擇越好；反之，BIC值越大，則噪聲模型選擇越差.

通過基于BIC數(shù)值分析法的HECTOR時間序列分析軟件對全部連續(xù)站點的三個方向進行解算，然后利用四種模型組合[14-15](WN+高斯化-廣義匹配(GGM)、WN+PL、WN+FN、WN+FN+RWN)分別計算出各連續(xù)站點在N、E、U三方向上的BIC值，圖4～6為各測站三方向四種噪聲模型解算得到的BIC差值結(jié)果.

圖4 各站點N方向差值

圖5 各站點E方向差值

圖6 各站點U方向差值

由上述各站點三個方向(N、E、U)的不同噪聲模型組合得到的BIC值，可以找到每個測站各個方向所對應的最小BIC值，最終確定內(nèi)蒙古自治區(qū)坐標時間序列在E方向的最優(yōu)模型為WN+FN，在N和U方向的最優(yōu)噪聲模型為WN+PL.

3 內(nèi)蒙古地區(qū)速度場估計

結(jié)合坐標時間序列噪聲模型分析結(jié)果，可得到內(nèi)蒙古地區(qū)陸態(tài)網(wǎng)絡連續(xù)站的速度場估計模型. 根據(jù)BIC數(shù)值分析得到的內(nèi)蒙古自治區(qū)坐標時間序列在E方向的最優(yōu)模型為WN+FN，在N和U方向的最優(yōu)噪聲模型WN+PL，對內(nèi)蒙古地區(qū)15個陸態(tài)網(wǎng)絡連續(xù)站的N、E、U方向的速度場進行估計，給出修正后的水平速度和垂直速度以及不確定度的統(tǒng)計結(jié)果如表3所示.

表3 N、E、U方向速度估值和不確定度統(tǒng)計結(jié)果

內(nèi)蒙古地區(qū)在水平方向運動的平均速率是30.635 mm/a，水平運動方向為26°57′51″SEE，圖7為水平方向的速度場；內(nèi)蒙古地區(qū)在豎直方向運動的平均速率是0.792 mm/a，圖8為垂直方向的速度場.

圖7 內(nèi)蒙古地區(qū)修正后的水平方向速度場

圖8 內(nèi)蒙古地區(qū)修正后的垂直方向速度場

統(tǒng)計內(nèi)蒙古地區(qū)15個站點的原始速度場數(shù)據(jù)，求得原始N方向標準差2.668 mm，E方向標準差2.125 mm，U方向標準差1.669 mm. 考慮有色噪聲后，內(nèi)蒙古地區(qū)陸態(tài)網(wǎng)連續(xù)站在E方向的最優(yōu)噪聲模型為WN+FN，在N、U方向采用最優(yōu)噪聲模型為WN+PL. 統(tǒng)計改正內(nèi)蒙古地區(qū)15個站點N、E、U三方向的速度場估計后，得出N方向標準差為2.528 mm，E方向標準差為2.091 mm，U方向標準差為1.632 mm，與原始速度場對比精度有了一定的提高，具體解算結(jié)果如表4所示.

表4 各站點在N、E、U方向修正后的速度場估計 mm·a-1

4 結(jié) 語

本文主要使用內(nèi)蒙古地區(qū)15個陸態(tài)網(wǎng)絡基準站近10年的時間序列觀測數(shù)據(jù)，通過原始觀測數(shù)據(jù)的處理，結(jié)合最優(yōu)噪聲模型，對內(nèi)蒙古速度場進行重新估計，最終確定內(nèi)蒙古自治區(qū)坐標時間序列在E方向的最優(yōu)模型為WN+FN，在N和U方向的最優(yōu)噪聲模型為WN+PL；在ITRF14框架下，最終得到內(nèi)蒙古境內(nèi)陸態(tài)網(wǎng)絡水平方向平均運動速率為30.653 mm/a，水平運動方向為26°57′51″SEE，垂直方向上平均運動速率為0.792 mm/a.