隋哲民,李建章,佟雪佳,高志鈺
(1.蘭州交通大學(xué) 測繪與地理信息學(xué)院/地理國情監(jiān)測技術(shù)應(yīng)用國家地方聯(lián)合工程研究中心/甘肅省地理國情監(jiān)測工程實驗室,蘭州 730070;2.中國地震局地質(zhì)研究所 地震動力學(xué)國家重點實驗室,北京 100029)
迄今為止,我國大陸構(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è)等提供思路。
目前,測繪工作者大多采用卡茨(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ù)值分析法。
以往在進行時間序列分析時普遍使用極大似然估計法,利用此類估計方法可以求解時間序列數(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]。
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)的組合噪聲模型。
實驗選取中國新疆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ù)站的速度場進行修正。
中國大陸構(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 站原始時間序列
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ù)處理后的時間序列
通過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”。
繪制未考慮最優(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。
繪制經(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é)果的可靠性。
本文選取中國地震局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)模型是十分必要的。