新疆陸態(tài)網(wǎng)絡(luò)坐標(biāo)時(shí)間序列噪聲模型分析

隋哲民，李建章，佟雪佳，高志鈺

（1.蘭州交通大學(xué) 測(cè)繪與地理信息學(xué)院/地理國(guó)情監(jiān)測(cè)技術(shù)應(yīng)用國(guó)家地方聯(lián)合工程研究中心/甘肅省地理國(guó)情監(jiān)測(cè)工程實(shí)驗(yàn)室，蘭州 730070；2.中國(guó)地震局地質(zhì)研究所 地震動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100029）

0 引言

迄今為止，我國(guó)大陸構(gòu)造環(huán)境監(jiān)測(cè)網(wǎng)絡(luò)即陸態(tài)網(wǎng)絡(luò)（crustal movement observation network of China,CMONOC）的連續(xù)運(yùn)行參考站（continuously operating reference stations,CORS）積累了數(shù)量可觀且連續(xù)的原始觀測(cè)數(shù)據(jù)，通過(guò)分析目標(biāo)測(cè)站坐標(biāo)時(shí)間序列的原始數(shù)據(jù)可以求出相應(yīng)地區(qū)板塊的速度場(chǎng)，進(jìn)而為地殼形變研究以及地震分析預(yù)報(bào)等提供堅(jiān)實(shí)的數(shù)據(jù)基礎(chǔ)[1]。自新疆地區(qū)陸態(tài)網(wǎng)絡(luò)建立完成后，許多國(guó)內(nèi)外專(zhuān)家學(xué)者對(duì)該地區(qū)的坐標(biāo)時(shí)間序列數(shù)據(jù)做了大量的相關(guān)研究，并建立了各自不同的速度場(chǎng)模型。文獻(xiàn)[2]通過(guò)對(duì)新疆加密帕米爾東北側(cè)地區(qū)的全球定位系統(tǒng)（global positioning system,GPS）監(jiān)測(cè)網(wǎng)的解算與復(fù)測(cè)，得到了該地區(qū)地殼形變速率圖及GPS基準(zhǔn)站的時(shí)間序列[2]。文獻(xiàn)[3]通過(guò)PODAP 軟件解算2011—2014 年間3 a的新疆陸態(tài)網(wǎng)絡(luò)連續(xù)站時(shí)間序列，得到了相應(yīng)的速度場(chǎng)，發(fā)現(xiàn)新疆西南部至東北部，基準(zhǔn)站的南北分量總體呈現(xiàn)遞減的趨勢(shì)，東西分量呈現(xiàn)出先遞增后遞減再次遞增的拱形變化趨勢(shì)[3]。文獻(xiàn)[4]采用最小二乘法來(lái)配置并估計(jì)隨機(jī)信號(hào)，在解算數(shù)據(jù)時(shí)引入赫爾默特（Helmert）方差分量估計(jì)來(lái)建立精度更高的速度場(chǎng)模型[4]。但是，由于新疆地區(qū)陸態(tài)網(wǎng)絡(luò)建立時(shí)間較晚，導(dǎo)致上述研究并未對(duì)較長(zhǎng)時(shí)間段即5 a 以上的時(shí)間序列及其有色噪聲進(jìn)行分析研究。倘若忽視所得坐標(biāo)時(shí)間序列數(shù)據(jù)中存在的噪聲，將無(wú)法有效分離有色噪聲與觀測(cè)得到的時(shí)間序列數(shù)據(jù)，進(jìn)而會(huì)對(duì)速度場(chǎng)的解算產(chǎn)生很大的影響。文獻(xiàn)[5]研究結(jié)果證實(shí)，要想獲得較為可靠且穩(wěn)定的時(shí)間序列最優(yōu)噪聲模型，至少需要積累超過(guò)5 a的時(shí)間序列觀測(cè)數(shù)據(jù)[5]。為此，本文選取新疆31 個(gè)陸態(tài)網(wǎng)絡(luò)連續(xù)站2010—2020 年間10 a的時(shí)間序列觀測(cè)數(shù)據(jù)進(jìn)行研究，確定3 個(gè)方向分量的最優(yōu)噪聲模型類(lèi)型及其各自所占比重，進(jìn)而得到新疆地區(qū)經(jīng)最優(yōu)噪聲模型改正后基于國(guó)際地球參考框架（international terrestrial reference frame,ITRF）下的速度場(chǎng)。新疆時(shí)間序列最優(yōu)噪聲模型的確定可為新疆高精度坐標(biāo)框架的研究提供參考，并為板塊運(yùn)動(dòng)、地殼形變監(jiān)測(cè)研究，以及新疆地區(qū)生產(chǎn)建設(shè)等提供思路。

1 噪聲分析方法

目前，測(cè)繪工作者大多采用卡茨（CATS）軟件以及赫克托（Hector）軟件來(lái)分析GPS 坐標(biāo)時(shí)間序列噪聲，其中CATS 軟件主要采用極大似然估計(jì)法（maximum likelihood estimation,MLE）來(lái)分析GPS 坐標(biāo)時(shí)間序列的噪聲特征。由于該款軟件解算時(shí)間序列時(shí)所采用的算法以及相應(yīng)的模型都較為準(zhǔn)確，故能滿足測(cè)繪工作者在解算時(shí)間序列數(shù)據(jù)時(shí)的精度需求，但是此軟件在處理較大量的數(shù)據(jù)時(shí)，解算數(shù)據(jù)的速度十分緩慢，無(wú)法給實(shí)際工作提供便利[6]。為了克服此類(lèi)難題，能夠既好又快地處理大量的時(shí)間序列數(shù)據(jù)，2012 年，葡萄牙研究人員Bos 等人開(kāi)發(fā)了Hector 時(shí)間序列分析軟件，該款軟件通過(guò)改進(jìn)算法使得數(shù)據(jù)處理速率得到了大幅提升[7]。這種經(jīng)過(guò)算法優(yōu)化后的極大似然估計(jì)法即為貝葉斯信息準(zhǔn)則（Bayesian information criterion,BIC）數(shù)值分析法。

1.1 極大似然估計(jì)法

以往在進(jìn)行時(shí)間序列分析時(shí)普遍使用極大似然估計(jì)法，利用此類(lèi)估計(jì)方法可以求解時(shí)間序列數(shù)據(jù)中所包含的絕大部分參數(shù)，例如截距、偏移、斜率、振幅等。故GPS 坐標(biāo)時(shí)間序列一般可以表示[8-10]為

式中：ti表示儒略日（modified Julian date,MJD）時(shí)間，單位為年；a為站點(diǎn)起始地殼位置參數(shù)；b為時(shí)間序列中線性變化的速率；年周期和半年周期的振幅分別使用參數(shù)c、d、e、f來(lái)表示；H（t）表示階躍函數(shù)（heaviside step function）；T指代發(fā)生階躍的時(shí)刻；參數(shù)gi指代同震偏移（offset）；vi表示觀測(cè)結(jié)果的殘差，即當(dāng)原始觀測(cè)值與預(yù)期結(jié)果不同時(shí)所產(chǎn)生的偏差[6]。

1.2 BIC 數(shù)值分析法

Hector 軟件可以用來(lái)估計(jì)線性趨勢(shì)項(xiàng)、高階多項(xiàng)式、季節(jié)項(xiàng)、其他周期項(xiàng)信號(hào)以及多種噪聲模型的組合。在解算最優(yōu)的噪聲模型時(shí)，Hector 軟件采用BIC 信息判別準(zhǔn)則[7]，該準(zhǔn)則是評(píng)價(jià)統(tǒng)計(jì)模型擬合結(jié)果優(yōu)良性的一種標(biāo)準(zhǔn)；采用對(duì)數(shù)似然函數(shù)來(lái)解算數(shù)據(jù)，并在解算過(guò)程中添加參數(shù)k來(lái)避免過(guò)度擬合。此對(duì)數(shù)函數(shù)的定義為

式中：N為實(shí)際觀測(cè)數(shù)據(jù)的數(shù)量（不包括缺失的數(shù)據(jù)）；r為由觀察殘差和遺漏殘差所組成的殘差矩陣；協(xié)方差矩陣C可分解為

式中：σ為殘差序列的標(biāo)準(zhǔn)差；為各類(lèi)噪聲的總和。σ的計(jì)算方法為

由det 可求已知常數(shù)C及N階矩陣A的行列式，又detA=CNdetA，故聯(lián)立式（2）、式（3）、式（4）可得

通過(guò)式（5）可以大致了解所估計(jì)的模型的復(fù)雜程度以及該模型所擬合數(shù)據(jù)的準(zhǔn)確性，具體定義為

式中：B為所求BIC 值；k為模型參數(shù)的個(gè)數(shù)。

似然函數(shù)L越大，對(duì)應(yīng)的模型越優(yōu)，即具有（相對(duì)）最小BIC 值的模型為最優(yōu)模型。在對(duì)不同組合模型進(jìn)行解算后，比對(duì)這些模型的BIC 值，進(jìn)而可以得到最優(yōu)的組合噪聲模型。

2 時(shí)間序列分析

實(shí)驗(yàn)選取中國(guó)新疆31 個(gè)CMONOC 連續(xù)站2010—2020 年的原始時(shí)間序列數(shù)據(jù)。采用Hector軟件對(duì)原始時(shí)間序列數(shù)據(jù)剔除奇異值后進(jìn)行最優(yōu)噪聲模型分析，并求出各假設(shè)模型相應(yīng)的BIC 值。最后根據(jù)貝葉斯準(zhǔn)則確定最優(yōu)噪聲模型，在得到最優(yōu)模型的基礎(chǔ)上對(duì)新疆陸態(tài)網(wǎng)絡(luò)連續(xù)站的速度場(chǎng)進(jìn)行修正。

2.1 時(shí)間序列獲取

中國(guó)大陸構(gòu)造環(huán)境監(jiān)測(cè)網(wǎng)是中國(guó)地殼運(yùn)動(dòng)觀測(cè)網(wǎng)的延續(xù)，在新疆境內(nèi)，總共設(shè)置了31 個(gè)站點(diǎn)，測(cè)站具體分布如圖1 所示。選取加米特（GAMIT）/格洛布克（GLOBK）軟件[11]解算得到的基于ITRF14 框架下的原始坐標(biāo)時(shí)間序列數(shù)據(jù)進(jìn)行解算，具體流程可參考文獻(xiàn)[12]。由于同一測(cè)站在不同長(zhǎng)度時(shí)段內(nèi)所求最優(yōu)噪聲模型、速度場(chǎng)差別較大，在研究分析測(cè)站最優(yōu)噪聲模型、速度場(chǎng)前需要指定該時(shí)間序列的時(shí)段[13]。本文所使用的數(shù)據(jù)來(lái)源于中國(guó)地震局全球衛(wèi)星導(dǎo)航系統(tǒng)（global navigation satellite system,GNSS）數(shù)據(jù)產(chǎn)品服務(wù)平臺(tái)[14]。

圖1 新疆維吾爾自治區(qū)境內(nèi)CMONOC 測(cè)站分布

由于篇幅有限，僅繪制XJBY、XJKE、XJYC 3 站的原始時(shí)間序列圖像，如圖2～圖4 所示。從圖2～圖4 可以看出，N(北）、E(東）方向分量觀測(cè)數(shù)據(jù)具有明顯線性變化趨勢(shì)：N 方向隨時(shí)間的變化而遞增且斜率較大；E 方向隨時(shí)間的變化而遞增且斜率較?。籙（天頂）方向的變化趨勢(shì)存在一定的周期性，且以1 a 周期來(lái)看最為顯著。原始時(shí)間序列觀測(cè)數(shù)據(jù)中存在明顯的奇異值（或稱(chēng)為外野值），必須在數(shù)據(jù)處理前對(duì)其進(jìn)行剔除。

圖2 XJBY 站原始時(shí)間序列

圖3 XJKE 站原始時(shí)間序列

圖4 XJYC 站原始時(shí)間序列

2.2 原始時(shí)間序列處理

Hector 軟件剔除時(shí)間序列奇異值是首先使用最小二乘法估計(jì)原始坐標(biāo)時(shí)間序列的線性趨勢(shì)、周期性趨勢(shì)，并去除線性趨勢(shì)和周期性趨勢(shì)以計(jì)算殘差，然后利用四分位數(shù)粗差探測(cè)方法對(duì)殘差進(jìn)行剔除。其具體步驟為：①基于最小二乘法，利用式（1）獲取殘差時(shí)間序列；②基于殘差時(shí)間序列，利用四分位距探測(cè)原始數(shù)據(jù)所包含的粗差，并剔除其中的奇異值；③重新擬合新獲取的殘差時(shí)間序列并重復(fù)步驟②，直至粗差完全剔除[7]。圖5～圖7 為XJBY、XJKE、XJYC 站數(shù)據(jù)處理后的時(shí)間序列圖像。

圖5 XJBY 站數(shù)據(jù)處理后的時(shí)間序列

圖6 XJKE 站數(shù)據(jù)處理后的時(shí)間序列

圖7 XJYC 站數(shù)據(jù)處理后的時(shí)間序列

2.3 最優(yōu)噪聲模型的確立

通過(guò)Hector 軟件對(duì)全部31 個(gè)連續(xù)站3 個(gè)方向分量的數(shù)據(jù)進(jìn)行解算，使用白色噪聲（WN）分別與閃爍噪聲（FN）、隨機(jī)漫步噪聲（RW）、冪律噪聲（PN）和一階高斯馬爾科夫噪聲（GGM）進(jìn)行組合。進(jìn)而使用這4 種組合模型計(jì)算出每個(gè)站所對(duì)應(yīng)的N、E、U 3 個(gè)方向的BIC 值，然后在同一分量中找到每個(gè)測(cè)站所對(duì)應(yīng)的最小BIC 值。其中，XJBC、XJBY、XJKC、XJKE、XJYC、XJYN 站在N 方向上噪聲模型解算得到的BIC 差值如圖8 所示。

圖8 不同測(cè)站N 方向上各噪聲模型的BIC 差值

通過(guò)比對(duì)不同模型的BIC 值，可獲悉各測(cè)站每個(gè)方向分量的最優(yōu)噪聲模型。全部測(cè)站不同方向的最優(yōu)噪聲模型所占百分比如表1 所示。

表1 各方向上最優(yōu)噪聲模型所占百分比 %

由表可知3 個(gè)方向分量時(shí)間序列與最優(yōu)噪聲模型之間的關(guān)系，且N 分量和E 分量、U 分量具有不同的噪聲特性：在N 方向上，最優(yōu)噪聲模型為組合模型“WN/GGM”；在E 方向和U 方向上，最優(yōu)噪聲模型為組合模型“WN/FN”。

3 顧及最優(yōu)噪聲模型的新疆速度場(chǎng)分析

3.1 未經(jīng)修正的速度場(chǎng)模型

繪制未考慮最優(yōu)噪聲模型改正的ITRF14 框架下新疆境內(nèi)CMONOC基準(zhǔn)站的水平運(yùn)動(dòng)速度場(chǎng)，如圖9 所示。新疆東部運(yùn)動(dòng)方向近E 向，而西南部地區(qū)為EN 向運(yùn)動(dòng)。未修正的新疆陸態(tài)網(wǎng)絡(luò)基準(zhǔn)站在ITRF14 框架下水平方向運(yùn)動(dòng)的平均速率為30.966 mm/a，運(yùn)動(dòng)方向?yàn)镹EE73°26′06″。

表2 給出了經(jīng)修正后ITRF14 框架下新疆境內(nèi)CMONOC基準(zhǔn)站的水平方向速度估值和中誤差統(tǒng)計(jì)結(jié)果。另外，31 個(gè)連續(xù)站的統(tǒng)計(jì)結(jié)果表明，E 方向速度的標(biāo)準(zhǔn)差為2.199 mm，N 方向標(biāo)準(zhǔn)差為6.384 mm，本文獲取的新疆CMONOC基準(zhǔn)站的水平速度精度較高。修正后基準(zhǔn)站E 方向速度平均值為 29.681 mm/a，N 方向速度平均值為8.828 mm/a；N 和E 方向即整體水平方向運(yùn)動(dòng)的平均速率為30.966 mm/a，優(yōu)勢(shì)方向?yàn)镹EE73°26′06″。

表2 水平速度估值和中誤差統(tǒng)計(jì) mm/a

繪制未考慮最優(yōu)噪聲模型改正的ITRF14 框架下新疆CMONOC基準(zhǔn)站的垂向運(yùn)動(dòng)速度場(chǎng)，結(jié)果如圖10 所示。

圖10 修正前的垂向速度場(chǎng)

從圖可以看出，新疆東北部運(yùn)動(dòng)方向近似向上運(yùn)動(dòng)，而西南部地區(qū)為向下運(yùn)動(dòng)。未修正新疆陸態(tài)網(wǎng)絡(luò)基準(zhǔn)站在ITRF14 框架下垂向運(yùn)動(dòng)的平均速率為0.207 mm/a。

表3 給出了未修正的ITRF14 框架下新疆境內(nèi)CMONOC基準(zhǔn)站的豎直方向速度估值和中誤差統(tǒng)計(jì)結(jié)果。

表3 垂向速度估值和中誤差統(tǒng)計(jì) mm/a

統(tǒng)計(jì)表明，U 方向速度的標(biāo)準(zhǔn)差為1.415 mm，可知本文獲取的新疆CMONOC基準(zhǔn)站的垂向速度精度較高。由表可見(jiàn)，基準(zhǔn)站U 方向速度平均值為0.207 mm/a。

3.2 修正后的速度場(chǎng)模型

繪制經(jīng)最優(yōu)噪聲模型改正后的ITRF14 框架下新疆CMONOC基準(zhǔn)站的水平運(yùn)動(dòng)速度場(chǎng)，如圖11所示。結(jié)果與中國(guó)地震局GNSS 數(shù)據(jù)產(chǎn)品服務(wù)平臺(tái)[14]結(jié)果基本一致，即圖9 所示。新疆東部運(yùn)動(dòng)方向?yàn)榻麰 向，而西南部地區(qū)為E—N 向運(yùn)動(dòng)。修正后新疆陸態(tài)網(wǎng)絡(luò)基準(zhǔn)站在ITRF14 框架下水平方向運(yùn)動(dòng)的平均速率為 30.976 mm/a，運(yùn)動(dòng)方向?yàn)镹EE73°23′27″。

圖11 修正后的水平速度場(chǎng)

表4 給出了經(jīng)最優(yōu)噪聲模型修正后ITRF14 框架下新疆境內(nèi)CMONOC基準(zhǔn)站的水平方向速度估值和中誤差統(tǒng)計(jì)結(jié)果。統(tǒng)計(jì)表明，E 方向速度的標(biāo)準(zhǔn)差為2.167 mm，N 方向標(biāo)準(zhǔn)差為6.368 mm，修正后新疆CMONOC基準(zhǔn)站的水平速度精度較高，且整體精度優(yōu)于未經(jīng)最優(yōu)噪聲模型修正的原始速度場(chǎng)模型。表中可以看出，基準(zhǔn)站E 方向速度平均值為29.684 mm/a，N 方向速度平均值為8.854 mm/a。

表4 水平速度估值和中誤差統(tǒng)計(jì)表 mm/a

繪制經(jīng)最優(yōu)噪聲模型改正后的ITRF14 框架下新疆垂向運(yùn)動(dòng)速度場(chǎng)，如圖12 所示。結(jié)果與中國(guó)地震局 GNSS 數(shù)據(jù)產(chǎn)品服務(wù)平臺(tái)[14]結(jié)果基本一致，即如圖10 所示。修正后新疆陸態(tài)網(wǎng)絡(luò)基準(zhǔn)站在ITRF14 框架下垂向運(yùn)動(dòng)的平均速率為0.153 mm/a。

圖12 修正后的垂向速度場(chǎng)

表5 給出了經(jīng)最優(yōu)噪聲模型修正后ITRF14 框架下新疆境內(nèi)CMONOC基準(zhǔn)站的豎直方向速度估值和中誤差統(tǒng)計(jì)結(jié)果。

表5 垂向速度估值和中誤差統(tǒng)計(jì) mm/a

比較修正前后的標(biāo)準(zhǔn)差：本文獲取的新疆CMONOC基準(zhǔn)站的垂向速度精度較高，U 方向速度的標(biāo)準(zhǔn)差為1.458 mm；且整體精度優(yōu)于未經(jīng)最優(yōu)噪聲模型修正的原始速度場(chǎng)模型。

對(duì)比修正前后的速度場(chǎng)模型可知：經(jīng)最優(yōu)噪聲模型修正后的速度場(chǎng)精度明顯優(yōu)于未修正的速度場(chǎng)，二者在水平方向上的運(yùn)動(dòng)速率最大差值為1.394 mm/a；水平運(yùn)動(dòng)方向角最大差值可達(dá)1°50′59″；垂直方向上的運(yùn)動(dòng)速率最大差值為1.730 mm/a。因此，在處理時(shí)間序列數(shù)據(jù)以及解算速度場(chǎng)時(shí)，有必要考慮有色噪聲的最優(yōu)模型。

在新疆陸態(tài)網(wǎng)絡(luò)31 個(gè)測(cè)站中，南北分量以塔什庫(kù)爾干站（TASH）的速率值最大，為22.699 mm/a，芨芨臺(tái)子站（XJJJ）的為最小，速率為0.467 mm/a；東西分量以烏拉斯臺(tái)站（XJWL）為最大，速率為32.188 mm/a，布倫口站（XJBL）為最小，速率為24.283 mm/a；垂向分量以克拉瑪依站（XJKL）為最大，速率為6.008 mm/a，木壘站（XJML）為最小，速率為0.004 mm/a?？梢钥闯?，東西分量運(yùn)動(dòng)速率的最值之間差距為7.905 mm/a，差距并不大。其原因可能是新疆西南地區(qū)受印度板塊的擠壓，使其向N 方向運(yùn)動(dòng)速率增大，而阿勒泰及天山東部地區(qū)沿N 方向運(yùn)動(dòng)速率較小，也就是說(shuō)從新疆西南部至東北部，在N 方向上的運(yùn)動(dòng)速率依次遞減，且過(guò)渡平穩(wěn)?？傮w向E 方向沿順時(shí)針運(yùn)動(dòng)，這與文獻(xiàn)[4]的研究結(jié)果一致。在垂直方向上，可以看出盆山結(jié)合部垂直運(yùn)動(dòng)速率相對(duì)較高且垂直向運(yùn)動(dòng)速率差值較大，進(jìn)而可以推斷天山及其周邊區(qū)域整體呈現(xiàn)隆起的趨勢(shì)[3]，進(jìn)一步證明了文獻(xiàn)[3]研究結(jié)果的可靠性。

4 結(jié)束語(yǔ)

本文選取中國(guó)地震局GNSS 數(shù)據(jù)產(chǎn)品服務(wù)平臺(tái)提供的新疆境內(nèi)31 個(gè)CMONOC 連續(xù)站10 a的坐標(biāo)時(shí)間序列觀測(cè)數(shù)據(jù)，采用貝葉斯信息準(zhǔn)則獲取了新疆地區(qū)各坐標(biāo)分量對(duì)應(yīng)的最優(yōu)噪聲模型，并確定了修正后的新疆陸態(tài)網(wǎng)絡(luò)速度場(chǎng)。最后得到如下結(jié)論：

1）通過(guò)10 a的觀測(cè)數(shù)據(jù)，得到了較為穩(wěn)定且可靠的3 個(gè)方向分量時(shí)間序列與最優(yōu)噪聲模型之間的關(guān)系，且N 分量和E、U 分量具有明顯不同的噪聲特性：在N 方向上，最優(yōu)噪聲模型為組合模型“WN/GGM”；在E 方向和U 方向上，最優(yōu)噪聲模型均為組合模型“WN/FN”。

2）量化了3 個(gè)方向分量上最優(yōu)噪聲模型所占百分比：在N 方向上，組合模型“WN/GGM”噪聲模型占比70.97%；E 方向上，組合模型“WN/FN”噪聲模型占比51.61%；在U 方向上，組合模型“WN/FN”噪聲模型占比80.65%。

3）經(jīng)最優(yōu)噪聲模型修正后的速度場(chǎng)精度明顯優(yōu)于未修正的速度場(chǎng)，二者在水平方向上的運(yùn)動(dòng)速率最大差值為1.394 mm/a；水平運(yùn)動(dòng)方向角最大差值可達(dá)1°50′59″；垂直方向上的運(yùn)動(dòng)速率最大差值為1.730 mm/a。因此，在處理新疆時(shí)間序列數(shù)據(jù)以及解算速度場(chǎng)時(shí)考慮有色噪聲的最優(yōu)模型是十分必要的。