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

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

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

0 引言

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

1 噪聲分析方法

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

1.1 極大似然估計法

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

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

1.2 BIC 數(shù)值分析法

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

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

式中：σ為殘差序列的標準差；為各類噪聲的總和。σ的計算方法為

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

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

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

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

2 時間序列分析

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

2.1 時間序列獲取

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

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

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

圖2 XJBY 站原始時間序列

圖3 XJKE 站原始時間序列

圖4 XJYC 站原始時間序列

2.2 原始時間序列處理

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖10 修正前的垂向速度場

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

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

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

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

3.2 修正后的速度場模型

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

圖11 修正后的水平速度場

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

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

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

圖12 修正后的垂向速度場

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

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

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

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

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

4 結(jié)束語

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

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

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

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