東北陸態(tài)網(wǎng)坐標時序有色噪聲及環(huán)境負載分析

黃裕，隋哲民，金澤林

(1.遼寧工程技術(shù)大學 測繪與地理科學學院,遼寧 阜新 123000；2.蘭州交通大學 測繪與地理信息學院,蘭州 730070；3.中國能源建設(shè)集團遼寧電力勘測設(shè)計院有限公司,沈陽110000)

0 引言

中國大陸環(huán)境監(jiān)測網(wǎng)絡(luò)(CMONOC)在全國范圍內(nèi)已陸續(xù)建成260 個連續(xù)運行參考站(CORS) 目前已經(jīng)積累了數(shù)量可觀且連續(xù)的觀測數(shù)據(jù)[1].通過對觀測數(shù)據(jù)進行處理，可以分析相應(yīng)地區(qū)的速度場，進而為地殼運動形變規(guī)律研究以及地震分析預(yù)報等提供數(shù)據(jù)基礎(chǔ)[1].謝方等[2]通過最小二乘配置移動擬合推估法以及克里金插值法于2015 年初步建立東北地塊的速度場.于吉鵬等[3]通過對東北地區(qū)2011—2013 年的GAMIT 結(jié)果的解算初步獲得了該地區(qū)在全球框架下的速度場.朱李忠等[4]選取2012—2017 年中國東北及俄羅斯遠東地區(qū)東南部的GPS 觀測數(shù)據(jù)初步得到了該地區(qū)基于ITRF2014 下的速度場.但是，由于東北地區(qū)陸態(tài)網(wǎng)絡(luò)建立時間較晚，上述研究并未使用較長時間段即5 a 以上的時間序列，或者未對其有色噪聲(CN)及環(huán)境負載進行分析研究.倘若忽視坐標時間序列數(shù)據(jù)中存在的CN，將無法有效分離CN 與觀測的時間序列數(shù)據(jù)，進而會對速度場的解算產(chǎn)生很大的影響.管雅慧等[5]的研究結(jié)果也證實了在分析速度場時，需要積累超過5 a 的時間序列觀測數(shù)據(jù)才能獲得較為可靠且穩(wěn)定的時間序列最優(yōu)噪聲模型.同時，姜衛(wèi)平等[6]也指出在研究測站時間序列及其速度場時環(huán)境負載不可忽視，尤其是在沿海地區(qū)需要重視非潮汐海洋負載.其研究成果也證明，采用環(huán)境負載數(shù)據(jù)對坐標時間序列數(shù)據(jù)進行改正后，得到的非潮汐海洋負載位移對沿海區(qū)域測站的 GPS 會進一步優(yōu)化[6].

為此，本文選取東北地區(qū)16 個陸態(tài)網(wǎng)絡(luò)連續(xù)站2010—2021 年間 10 a 的時間序列觀測數(shù)據(jù)及同一時段的環(huán)境負載數(shù)據(jù)進行研究，確定了三個方向分量的最優(yōu)噪聲模型類型，并剔除環(huán)境負載的影響，進而得到了東北地區(qū)陸態(tài)網(wǎng)絡(luò)連續(xù)觀測站經(jīng)改正后基于國際地球參考框架(ITRF)下的速度場.為東北地區(qū)高精度坐標框架的研究提供了數(shù)據(jù)基礎(chǔ)，為我國東北地區(qū)板塊運動、地殼形變監(jiān)測、礦區(qū)生產(chǎn)建設(shè)等提供了重要的參考資料.

1 時間序列獲取

1.1 坐標時間序列的構(gòu)建

坐標時間序列函數(shù)模型可表示為[7-8]

式中：時間ti的單位為年，y(ti) 為測站ti時刻下的位移；參數(shù)a為截距，b為時間序列中線性變化的速率；第j個諧波對應(yīng)的振幅、角頻率以及相位分別表示為Aj、ωj和 φj；n為諧波的數(shù)量；gk為階躍項；Tk為產(chǎn)生階躍的時刻；H為階躍函數(shù)；參數(shù)gi為同震偏移(offset)；ε(ti)為隨機過程[9].

本文選取GAMIT/GLOBK 軟件[10]解算得到的基于ITRF14 框架下的東北地區(qū)16 個陸態(tài)網(wǎng)絡(luò)連續(xù)站近10 a 的原始坐標時間序列數(shù)據(jù)，具體的數(shù)據(jù)解算流程可參考文獻[11-12]，陸態(tài)網(wǎng)絡(luò)連續(xù)觀測站分布如圖1 所示.

圖1 東北地區(qū)CMONOC 測站分布圖

鑒于相同測站在不同時段所求得最優(yōu)噪聲模型特性不同，進行白噪聲(WN)分析時需選取10 a 以上的時間序列數(shù)據(jù)才能獲得較為穩(wěn)定的噪聲模型[13].本文選取的連續(xù)站時間序列數(shù)據(jù)均來源于中國地震局GNSS 數(shù)據(jù)產(chǎn)品服務(wù)平臺(http://www.cgps.ac.cn)，所選測站對應(yīng)的時間段如表1 所示.

表1 東北地區(qū)CMONOC 站點概況

在評價GNSS 時間序列數(shù)據(jù)精度時，通常使用加權(quán)均方根誤差(WRMS)來評價，WRMS 的值越小，表明對應(yīng)時間序列數(shù)據(jù)的重復性更佳[14].對應(yīng)站點北(N)、東(E)、天頂(U)方向的WRMS 值越小說明時間序列離散程度越小，重復性越好，精度越高.計算公式為

式中：nexp為預(yù)測值；ncale為實際觀測值；σexp為預(yù)測值中誤差.

其中，水平方向上的理想WRMS 值約為1～2 mm，垂直方向上的理想WRMS 值約為3～10 mm[14].東北地區(qū)陸態(tài)網(wǎng)絡(luò)連續(xù)站時間序列在N、E、U 方向上的WRMS 值統(tǒng)計結(jié)果如表2 所示.由表2 可知，東北地區(qū)大部分站點原始時間序列均滿足理想情況，但HLFY、HLMH 以及SUIY 站在N 方向上WRMS 值偏大，CHUN、HLMH、HRBN、JLCB、JLYJ 以及SUIY站在E 方向上WRMS 值偏大.

表2 東北連續(xù)站時間序列各分量WRMS 值 mm

1.2 殘差時間序列獲取

以LNDD 站為例，圖2 為該站的殘差時間序列圖像.可以看出，N、E 方向分量觀測數(shù)據(jù)均有明顯的線性變化趨勢，N 方向隨時間變化呈現(xiàn)向南運動趨勢且斜率較小，E 方向隨時間變化呈現(xiàn)向東運動趨勢且斜率較大；U 方向的變化趨勢呈現(xiàn)出明顯的周年項波動特性；原始時間序列觀測數(shù)據(jù)中存在明顯的奇異值(或稱為外野值)，必須在數(shù)據(jù)處理前對其進行剔除.

圖2 LNDD 站殘差時間序列圖像

Hector 軟件使用最小二乘法估計原始坐標時間序列的線性趨勢與周期性趨勢來剔除時間序列奇異值，進而通過去除線性趨勢和周期性趨勢計算殘差.在去趨勢項后利用四分位數(shù)粗差探測方法剔除殘差.其具體步驟為：1)基于最小二乘法，利用式(1)獲取殘差時間序列；2)基于殘差時間序列，利用四分位距探測原始數(shù)據(jù)所包含的粗差，并剔除其中的奇異值；3)重新擬合新獲取的殘差時間序列并再次進行步驟2，直至完全剔除粗差[11].

2 噪聲分析

目前，速度場研究大多采用CATS 軟件來分析坐標時間序列噪聲，但是在處理較大量的數(shù)據(jù)時，解算速度將十分緩慢，無法為數(shù)據(jù)解算提供便利[6].為了克服此難題，BOS 等[15]研究人員于2012 年開發(fā)了Hector 時間序列分析軟件，該軟件通過優(yōu)化算法大幅提高了數(shù)據(jù)處理速率.這種經(jīng)過算法優(yōu)化后的極大似然估計法(MLE)即為貝葉斯信息準則(BIC)數(shù)值分析法.

2.1 BIC 數(shù)值分析法

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

式中：N為實際觀測數(shù)據(jù)的數(shù)量(不包括缺失的數(shù)據(jù))；r為殘差矩陣；協(xié)方差陣C可分解為

式中：C為各類噪聲的總和；σ 為殘差序列的標準差

由 d etCA=CNdetA，式(3)～(5)可得

通過上式可以粗略知曉所解算模型的復雜程度以及該模型擬合數(shù)據(jù)的準確程度[15]，具體可定義為

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

2.2 噪聲模型分析

通過Hector 軟件對16 個連續(xù)站三個方向分量的數(shù)據(jù)進行解算，分別使用這4 種組合模型計算出每個站所對應(yīng)的N、E、U 三個方向的BIC 值，然后在同一分量中找到每個測站所對應(yīng)的最小BIC 值.全部測站E 方向上的噪聲模型解算所得BIC 差值，如圖3所示.

圖3 各測站E 方向BIC 差值

通過對不同組合噪聲模型對應(yīng)的BIC 值比較可知三個方向時間序列的不同組合噪聲模型間的關(guān)系，且N、E、U 分量具有不同的噪聲特性：在N 方向上，最優(yōu)噪聲模型為WN+FN，占比為56.25%；在E 方向，最優(yōu)噪聲模型為WN+PL，占比為56.25%；U 方向上，最優(yōu)噪聲模型為WN+FN，占比為87.5%.

3 顧及CN 的東北地區(qū)速度場分析

繪制經(jīng)最優(yōu)噪聲模型改正后的ITRF14 框架下東北地區(qū)CMONOC 基準站的水平運動速度場，如圖4 所示.東北地區(qū)整體呈現(xiàn)近SE 向運動，只考慮CN 修正后東北陸態(tài)網(wǎng)絡(luò)基準站在ITRF14 框架下N 方向上運動的平均速率為-12.902 mm/a，E 方向上運動的平均速率為26.852 mm/a，U 方向上運動的平均速率為0.312 mm/a.

圖4 修正后的水平速度場

繪制經(jīng)最優(yōu)噪聲模型改正后的ITRF14 框架下東北地區(qū)水平運動速度場，如圖4 所示.

修正后東北陸態(tài)網(wǎng)絡(luò)基準站的水平速度精度較高，E 方向速度的標準差為2.462 mm，N 方向標準差為1.369 mm，且整體精度優(yōu)于未經(jīng)修正的原始速度場模型(原始E 方向速度標準差為2.664 mm，N 方向速度標準差為1.511 mm).

繪制經(jīng)最優(yōu)噪聲模型改正后的ITRF14 框架下東北地區(qū)垂向運動速度場，如圖5 所示.

圖5 修正后的垂向速度場

本文獲取的東北CMONOC 基準站的垂向速度精度較高，U 方向速度標準差為0.892 mm，且整體精度優(yōu)于未經(jīng)修正的原始速度場模型(原始U 方向速度標準差為0.981 mm).

對比修正前后的速度場模型可知：經(jīng)最優(yōu)噪聲模型修正后的速度場精度明顯優(yōu)于未修正的速度場且二者在水平方向上的運動速率最大差值為0.564 mm/a；垂直方向上的運動速率最大差值為1.285 mm/a.因而，在對時間序列數(shù)據(jù)進行處理以及解算速度場之前需要考慮CN 的最優(yōu)模型.

4 顧及CN 及環(huán)境負載的東北地區(qū)速度場估計

在計算環(huán)境負載時采用GFZ (德國地學中心)的數(shù)據(jù)：非潮汐大氣負載時間分辨率為3 h，空間分辨率為0.5°×0.5°；非潮汐海洋負載時間分辨率為3 h，空間分辨率為0.5°×0.5°；水文負載時間分辨率為24 h，空間分辨率為0.5°×0.5°[16].在解算時選取與各站點一致的觀測時段，這樣可以避免環(huán)境負載中非線性趨勢對原始時間序列數(shù)據(jù)的影響，進而避免在測站速度估計時產(chǎn)生偏差.圖6 為各測站不同方向上時間序列環(huán)境負載的平均值.

圖6 各測站環(huán)境負載時間序列圖像

結(jié)合GFZ 的負載數(shù)據(jù)分析，由圖6 可知U 方向上環(huán)境負載對東北地區(qū)時間序列引起的負荷形變尤為嚴重.

由表3 可知，在考慮環(huán)境負載后，部分測站的噪聲會發(fā)生變化：在N 方向上，最優(yōu)噪聲模型WN+FN的占比提高至75%；在E 方向，最優(yōu)噪聲模型WN+PL的占比提高至68.75%；U 方向上最優(yōu)噪聲模型WN+FN的占比降低至81.25%.結(jié)合環(huán)境負載影響對東北地區(qū)16 個陸態(tài)網(wǎng)絡(luò)連續(xù)觀測站的N、E、U 方向速度場進行估計，結(jié)果如表4 所示.

表3 環(huán)境負載改正后各方向最優(yōu)噪聲模型變化

表3 (續(xù))

表4 東北陸態(tài)網(wǎng)絡(luò)連續(xù)站區(qū)域速度場估計 mm/a

東北地區(qū)陸態(tài)網(wǎng)絡(luò)連續(xù)站總體上向E 方向沿順時針運動這與文獻[3]的研究結(jié)果一致.在垂直方向上，可以看出北部地區(qū)垂直運動速率相對較低且垂直向運動速率差值較大，進而可以推斷西南部地區(qū)沉降速率大于東部地區(qū)[4]，進一步證明了文獻[4]研究結(jié)果的可靠性.

5 結(jié)束語

本文采用東北地區(qū)16 個陸態(tài)網(wǎng)連續(xù)站10 a的坐標時間序列觀測數(shù)據(jù)，利用BIC 信息準則確定了東北地區(qū)各坐標分量對應(yīng)的最優(yōu)噪聲模型，并在考慮環(huán)境負載的基礎(chǔ)上確定了修正后的東北地區(qū)陸態(tài)網(wǎng)絡(luò)速度場.最后得到如下結(jié)論：

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

2)考慮CN 及環(huán)境負載后，東北地區(qū)陸態(tài)網(wǎng)絡(luò)基于ITRF14 框架下在N 方向上運動的平均速率為-13.003 mm/a，E 方向上運動的平均速率為27.020 mm/a，U 方向上運動的平均速率為0.528 mm/a.

3)考慮CN 及環(huán)境負載的速度場精度明顯優(yōu)于未修正的速度場，二者在水平方向上的運動速率最大差值為0.802 mm/a；垂直方向上的運動速率最大差值為1.377 6 mm/a.因此，在處理東北地區(qū)時間序列數(shù)據(jù)以及解算速度場時考慮CN 及環(huán)境負載影響是十分必要的.