山西陸態(tài)網(wǎng)絡(luò)坐標(biāo)時(shí)間序列噪聲及速度場分析

隋哲民 李建章 王思凡 高志鈺

1 蘭州交通大學(xué)測繪與地理信息學(xué)院，蘭州市安寧西路88號(hào)，7300702 地理國情監(jiān)測技術(shù)應(yīng)用國家地方聯(lián)合工程研究中心，蘭州市安寧西路88號(hào)，730070 3 甘肅省地理國情監(jiān)測工程實(shí)驗(yàn)室，蘭州市安寧西路88號(hào)，7300704 中國地震局地質(zhì)研究所地震動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，北京市華嚴(yán)里甲1號(hào)，100029

中國大陸構(gòu)造環(huán)境監(jiān)測網(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ù)的觀測資料，通過分析這些數(shù)據(jù)可以得到相應(yīng)地區(qū)的速度場，進(jìn)而為地殼運(yùn)動(dòng)形變規(guī)律研究及地震分析預(yù)報(bào)等提供數(shù)據(jù)基礎(chǔ)[1]。山西省陸態(tài)網(wǎng)絡(luò)初步建成于2009年，共有10個(gè)連續(xù)觀測站，其中SXCZ站由于滑坡被搬遷重建并更名為SXCH站。曾波等[2]利用2009～2011年的數(shù)據(jù)初步確定了山西省陸態(tài)網(wǎng)絡(luò)的速度場；劉瑞春等[3]通過分析2010～2014年山西斷陷帶陸態(tài)網(wǎng)絡(luò)的時(shí)間序列數(shù)據(jù)，并結(jié)合GPS基線結(jié)果，從整體上研究了山西斷陷帶南北2個(gè)區(qū)域的變形特征；成誠等[4]采用GAMIT/GLOBK軟件計(jì)算時(shí)間序列，得到山西省GNSS基準(zhǔn)站在水平方向呈東偏南向運(yùn)動(dòng)，垂直方向呈周期運(yùn)動(dòng)的結(jié)論。但由于山西地區(qū)陸態(tài)網(wǎng)絡(luò)建立的時(shí)間較晚，上述研究使用的數(shù)據(jù)積累時(shí)長均不超過5 a，且未對(duì)有色噪聲進(jìn)行分析研究。若忽視坐標(biāo)時(shí)間序列數(shù)據(jù)中存在的噪聲，將無法有效分離有色噪聲與觀測數(shù)據(jù)，進(jìn)而對(duì)速度場的解算產(chǎn)生很大影響。管雅慧等[5]的研究結(jié)果也證實(shí)，要想獲得較為可靠且穩(wěn)定的坐標(biāo)時(shí)間序列最優(yōu)噪聲模型，至少需積累5 a以上的觀測數(shù)據(jù)。

為此，本文選取山西地區(qū)10個(gè)陸態(tài)網(wǎng)絡(luò)連續(xù)站2010-11～2021-01近 10 a的觀測數(shù)據(jù)進(jìn)行研究，確定N、E、U分量的最優(yōu)噪聲模型類型，進(jìn)而得到山西地區(qū)經(jīng)最優(yōu)噪聲模型改正后基于國際地球參考框架(international terrestrial reference frame, ITRF)下的速度場，為山西地區(qū)高精度坐標(biāo)框架的研究提供支撐，也為板塊運(yùn)動(dòng)、地殼形變監(jiān)測研究及山西礦區(qū)生產(chǎn)建設(shè)等提供重要的參考資料。

1 噪聲分析方法

目前，測繪工作者大多采用CATS軟件來分析坐標(biāo)時(shí)間序列的噪聲，該軟件解算時(shí)間序列時(shí)采用的算法及相應(yīng)的模型都較為準(zhǔn)確，但在處理較大量的數(shù)據(jù)時(shí)，解算速度十分緩慢，無法給實(shí)際工作提供便利[6]。為此，本文使用Hector時(shí)間序列分析軟件，該軟件使用貝葉斯數(shù)值分析法(Bayesian information criterion,BIC)[7]，數(shù)據(jù)處理速率得到大幅提高。

2 時(shí)間序列獲取

山西省CMONOC測站具體分布如圖1所示。本文以GAMIT/GLOBK軟件[8]解算得到的基于ITRF2014框架下的山西地區(qū)10個(gè)陸態(tài)網(wǎng)絡(luò)連續(xù)站近10 a的原始坐標(biāo)時(shí)間序列數(shù)據(jù)為基礎(chǔ)，具體數(shù)據(jù)解算流程可參考文獻(xiàn)[9]。

圖1 山西境內(nèi)CMONOC測站分布Fig.1 Distribution of CMONOC stations in Shanxi

由于同一測站在不同長度的時(shí)段內(nèi)求得的最優(yōu)噪聲模型和速度場差別較大，在研究分析測站最優(yōu)噪聲模型和速度場前需要指定時(shí)間序列的時(shí)段長度[10]。本文數(shù)據(jù)來源于中國地震局GNSS數(shù)據(jù)產(chǎn)品服務(wù)平臺(tái)(http:∥www.cgps.ac.cn.)，所選測站對(duì)應(yīng)的時(shí)間段如表1所示。

表1 山西省CMONOC站點(diǎn)概況

通常使用加權(quán)均方根誤差(WRMS)來評(píng)價(jià)GNSS坐標(biāo)時(shí)間序列的精度，WRMS值越小，表明對(duì)應(yīng)數(shù)據(jù)的重復(fù)性越佳[11]。其中，水平方向上的WRMS理想值為1～2 mm，垂直方向上的WRMS理想值為3～10 mm。山西省陸態(tài)網(wǎng)絡(luò)連續(xù)站坐標(biāo)時(shí)間序列在N、E、U方向上的WRMS值統(tǒng)計(jì)結(jié)果如表2(單位mm)所示，可以看出，SXCZ、SXGX、SXTY站在N方向上的WRMS值大于理想值，SXCZ、SXTY站在E方向上的WRMS值大于理想值，其余各站在水平方向和垂直方向上的重復(fù)性均相對(duì)較好。

表2 山西省連續(xù)站坐標(biāo)時(shí)間序列各分量WRMS值

以SXDT站為例，圖2為該站的原始坐標(biāo)時(shí)間序列。可以看出，N、E方向的觀測數(shù)據(jù)具有明顯的線性變化趨勢(shì)，N方向隨時(shí)間遞增且斜率較小，E方向隨時(shí)間遞增且斜率較大；U方向的變化趨勢(shì)存在一定的周期性，且以1 a周期最為顯著；原始坐標(biāo)時(shí)間序列觀測數(shù)據(jù)中存在明顯的奇異值(或稱為外野值)，必須在數(shù)據(jù)處理前進(jìn)行剔除。

圖2 SXDT站原始坐標(biāo)時(shí)間序列Fig.2 Original coordinate time series of SXDT station

2.1 殘差時(shí)間序列獲取

Hector軟件剔除時(shí)間序列奇異值的具體步驟為：1)基于最小二乘法獲取殘差時(shí)間序列；2)基于殘差時(shí)間序列，利用四分位距探測原始數(shù)據(jù)包含的粗差，并剔除其中的奇異值；3)重新擬合新獲取的殘差時(shí)間序列并重復(fù)步驟2)，直至粗差被完全剔除[7]。圖3為SXDT站坐標(biāo)殘差時(shí)間序列。

圖3 SXDT站坐標(biāo)殘差時(shí)間序列Fig.3 Coordinate residual time series of SXDT station

2.2 噪聲模型分析

使用Hector軟件對(duì)10個(gè)連續(xù)站3個(gè)方向分量的觀測數(shù)據(jù)進(jìn)行解算，分別使用4種模型計(jì)算每個(gè)測站對(duì)應(yīng)N、E、U方向的BIC值，然后在同一分量中找到每個(gè)測站對(duì)應(yīng)的最小BIC值，進(jìn)而得出不同方向上的最優(yōu)噪聲模型。全部測站各方向上的噪聲模型分布如圖4所示。

圖4 各方向噪聲模型分布Fig.4 Noise model distribution in eachdirection

由圖4可知，N方向、E方向和U方向具有不同的噪聲特性：在N方向上，最優(yōu)噪聲模型為WN+PL；在E和U方向上，最優(yōu)噪聲模型為WN+FN。

3 顧及有色噪聲的山西速度場估計(jì)

繪制經(jīng)最優(yōu)噪聲模型改正后ITRF2014框架下的山西省CMONOC基準(zhǔn)站水平方向運(yùn)動(dòng)速度場，結(jié)果如圖5所示。可以看出，修正后山西陸態(tài)網(wǎng)絡(luò)基準(zhǔn)站水平方向的運(yùn)動(dòng)速率平均為33.542 mm/a，運(yùn)動(dòng)方向?yàn)镾EE 26°22′12″。

圖5 修正后的水平方向速度場Fig.5 Corrected horizontal velocity field

表3給出修正后ITRF2014框架下的山西省CMONOC基準(zhǔn)站水平方向速度估值和不確定度統(tǒng)計(jì)結(jié)果。由表可知，E方向速度標(biāo)準(zhǔn)差為0.751 mm/a，N方向速度標(biāo)準(zhǔn)差為0.552 mm/a，修正后山西省CMONOC基準(zhǔn)站水平方向的速度精度較高，且整體精度優(yōu)于未經(jīng)修正的原始速度場模型(原始模型中，E方向速度標(biāo)準(zhǔn)差為1.715 mm/a，N方向速度標(biāo)準(zhǔn)差為1.723 mm/a)。

表3 水平速度估值和不確定度統(tǒng)計(jì)

繪制經(jīng)最優(yōu)噪聲模型改正后ITRF2014框架下的山西垂直方向運(yùn)動(dòng)速度場，如圖6所示。

圖6 修正后的垂直方向速度場Fig.6 Corrected vertical velocity field

表4給出經(jīng)最優(yōu)噪聲模型修正后ITRF2014框架下的山西省CMONOC基準(zhǔn)站垂直方向的速度估值和不確定度統(tǒng)計(jì)結(jié)果。由表可知，本文獲取的山西省CMONOC基準(zhǔn)站垂直方向的速度精度較高，速度標(biāo)準(zhǔn)差為0.906 mm/a，且整體精度優(yōu)于未經(jīng)修正的原始速度場模型(原始模型垂直方向的速度標(biāo)準(zhǔn)差為1.161 mm/a)。

表4 垂向速度估值和不確定度統(tǒng)計(jì)

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

采用N方向的最優(yōu)噪聲模型WN+PL與E、U方向的最佳噪聲模型WN+FN，結(jié)合BIC數(shù)值分析法，對(duì)山西省10個(gè)陸態(tài)網(wǎng)絡(luò)連續(xù)觀測站在三方向的速度場進(jìn)行估計(jì)，結(jié)果如表5所示?？梢钥闯?，山西地區(qū)陸態(tài)網(wǎng)絡(luò)連續(xù)站總體呈向E沿順時(shí)針運(yùn)動(dòng)，這與成誠等[4]的研究結(jié)果一致；在垂直方向上，南部地區(qū)的運(yùn)動(dòng)速率相對(duì)較低且運(yùn)動(dòng)速率差值較大，可以推斷，南部地區(qū)沉降速率大于北部地區(qū)，這與曾波等[2]的研究結(jié)果一致。

表5 山西省陸態(tài)網(wǎng)絡(luò)連續(xù)站區(qū)域速度場估計(jì)

4 結(jié) 語

本文選取中國地震局GNSS數(shù)據(jù)產(chǎn)品服務(wù)平臺(tái)(http:∥www.cgps.ac.cn)提供的山西省10個(gè)CMONOC連續(xù)站10 a的坐標(biāo)時(shí)間序列觀測數(shù)據(jù)，利用BIC數(shù)值分析法獲取了山西地區(qū)各坐標(biāo)分量對(duì)應(yīng)的最優(yōu)噪聲模型，并確定了修正后的山西陸態(tài)網(wǎng)絡(luò)速度場，得到如下結(jié)論：

1)通過近10 a的觀測數(shù)據(jù)得到了較為穩(wěn)定且可靠的三方向坐標(biāo)時(shí)間序列與最優(yōu)噪聲模型的關(guān)系，且N方向和E、U方向具有明顯不同的噪聲特性：在N方向上，最優(yōu)噪聲模型為WN+PL；在E、U方向上，最優(yōu)噪聲模型均為WN+FN。

2)山西陸態(tài)網(wǎng)絡(luò)基于ITRF2014框架下的速度場，在水平方向的平均運(yùn)動(dòng)速率為33.542 mm/a，運(yùn)動(dòng)方向?yàn)镾EE 26°22′12″；垂直方向的平均運(yùn)動(dòng)速率為2.214 mm/a。

3)經(jīng)最優(yōu)噪聲模型修正后的速度場精度明顯優(yōu)于未修正的速度場，二者在水平方向上的運(yùn)動(dòng)速率最大差值為4.025 mm/a；水平運(yùn)動(dòng)方向角的最大差值可達(dá)6°50′15″；垂直方向上的運(yùn)動(dòng)速率最大差值為2.627 mm/a。因此，在處理山西地區(qū)坐標(biāo)時(shí)間序列數(shù)據(jù)及解算速度場時(shí)考慮有色噪聲的影響十分必要。