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

李鑫，姚德新

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

0 引 言

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

1 基準(zhǔn)站坐標(biāo)時(shí)間序列及其精度

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

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

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

加權(quán)均方根誤差(WRMS)是評(píng)價(jià)GNSS時(shí)間序列精度常用指標(biāo)，由文獻(xiàn)[9]可知，對(duì)應(yīng)站點(diǎn)N、E、U方向的WRMS值越小說明時(shí)間序列離散程度越小，重復(fù)性越好，精度越高. 計(jì)算公式為

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

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

2 坐標(biāo)時(shí)間序列噪聲分析

2.1 坐標(biāo)時(shí)間序列預(yù)處理

原始坐標(biāo)時(shí)間序列可表示為[11]

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

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

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

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

2.2 坐標(biāo)時(shí)間序列噪聲模型分析

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

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

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

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

于是，BIC值的具體計(jì)算為

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

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

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

圖4 各站點(diǎn)N方向差值

圖5 各站點(diǎn)E方向差值

圖6 各站點(diǎn)U方向差值

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

3 內(nèi)蒙古地區(qū)速度場(chǎng)估計(jì)

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

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

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

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

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

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

表4 各站點(diǎn)在N、E、U方向修正后的速度場(chǎng)估計(jì) mm·a-1

4 結(jié) 語

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