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

（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ù)，通過分析目標(biāo)測站坐標(biāo)時間序列的原始數(shù)據(jù)可以求出相應(yīng)地區(qū)板塊的速度場，進而為地殼形變研究以及地震分析預(yù)報等提供堅實的數(shù)據(jù)基礎(chǔ)[1]。自新疆地區(qū)陸態(tài)網(wǎng)絡(luò)建立完成后，許多國內(nèi)外專家學(xué)者對該地區(qū)的坐標(biāo)時間序列數(shù)據(jù)做了大量的相關(guān)研究，并建立了各自不同的速度場模型。文獻[2]通過對新疆加密帕米爾東北側(cè)地區(qū)的全球定位系統(tǒng)（global positioning system,GPS）監(jiān)測網(wǎng)的解算與復(fù)測，得到了該地區(qū)地殼形變速率圖及GPS基準(zhǔn)站的時間序列[2]。文獻[3]通過PODAP 軟件解算2011—2014 年間3 a的新疆陸態(tài)網(wǎng)絡(luò)連續(xù)站時間序列，得到了相應(yīng)的速度場，發(fā)現(xiàn)新疆西南部至東北部，基準(zhǔn)站的南北分量總體呈現(xiàn)遞減的趨勢，東西分量呈現(xiàn)出先遞增后遞減再次遞增的拱形變化趨勢[3]。文獻[4]采用最小二乘法來配置并估計隨機信號，在解算數(shù)據(jù)時引入赫爾默特（Helmert）方差分量估計來建立精度更高的速度場模型[4]。但是，由于新疆地區(qū)陸態(tài)網(wǎng)絡(luò)建立時間較晚，導(dǎo)致上述研究并未對較長時間段即5 a 以上的時間序列及其有色噪聲進行分析研究。倘若忽視所得坐標(biāo)時間序列數(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)噪聲模型的確定可為新疆高精度坐標(biāo)框架的研究提供參考，并為板塊運動、地殼形變監(jiān)測研究，以及新疆地區(qū)生產(chǎn)建設(shè)等提供思路。

1 噪聲分析方法

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

1.1 極大似然估計法

以往在進行時間序列分析時普遍使用極大似然估計法，利用此類估計方法可以求解時間序列數(shù)據(jù)中所包含的絕大部分參數(shù)，例如截距、偏移、斜率、振幅等。故GPS 坐標(biāo)時間序列一般可以表示[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 信息判別準(zhǔn)則[7]，該準(zhǔn)則是評價統(tǒng)計模型擬合結(jié)果優(yōu)良性的一種標(biāo)準(zhǔn)；采用對數(shù)似然函數(shù)來解算數(shù)據(jù)，并在解算過程中添加參數(shù)k來避免過度擬合。此對數(shù)函數(shù)的定義為

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

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

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

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

式中：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ù)貝葉斯準(zhǔn)則確定最優(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 框架下的原始坐標(biāo)時間序列數(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 方向隨時間的變化而遞增且斜率較?。籙（天頂）方向的變化趨勢存在一定的周期性，且以1 a 周期來看最為顯著。原始時間序列觀測數(shù)據(jù)中存在明顯的奇異值（或稱為外野值），必須在數(shù)據(jù)處理前對其進行剔除。

圖2 XJBY 站原始時間序列

圖3 XJKE 站原始時間序列

圖4 XJYC 站原始時間序列

2.2 原始時間序列處理

Hector 軟件剔除時間序列奇異值是首先使用最小二乘法估計原始坐標(biāo)時間序列的線性趨勢、周期性趨勢，并去除線性趨勢和周期性趨勢以計算殘差，然后利用四分位數(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基準(zhǔn)站的水平運動速度場，如圖9 所示。新疆東部運動方向近E 向，而西南部地區(qū)為EN 向運動。未修正的新疆陸態(tài)網(wǎng)絡(luò)基準(zhǔn)站在ITRF14 框架下水平方向運動的平均速率為30.966 mm/a，運動方向為NEE73°26′06″。

表2 給出了經(jīng)修正后ITRF14 框架下新疆境內(nèi)CMONOC基準(zhǔn)站的水平方向速度估值和中誤差統(tǒng)計結(jié)果。另外，31 個連續(xù)站的統(tǒng)計結(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 方向即整體水平方向運動的平均速率為30.966 mm/a，優(yōu)勢方向為NEE73°26′06″。

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

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

圖10 修正前的垂向速度場

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

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

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

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

3.2 修正后的速度場模型

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

圖11 修正后的水平速度場

表4 給出了經(jīng)最優(yōu)噪聲模型修正后ITRF14 框架下新疆境內(nèi)CMONOC基準(zhǔn)站的水平方向速度估值和中誤差統(tǒng)計結(jié)果。統(tǒng)計表明，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)噪聲模型修正的原始速度場模型。表中可以看出，基準(zhǔn)站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ò)基準(zhǔn)站在ITRF14 框架下垂向運動的平均速率為0.153 mm/a。

圖12 修正后的垂向速度場

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

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

比較修正前后的標(biāo)準(zhǔn)差：本文獲取的新疆CMONOC基準(zhǔn)站的垂向速度精度較高，U 方向速度的標(biāo)準(zhǔn)差為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)?？傮w向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的坐標(biāo)時間序列觀測數(shù)據(jù)，采用貝葉斯信息準(zhǔn)則獲取了新疆地區(qū)各坐標(biāo)分量對應(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)模型是十分必要的。