張幸，魏冠軍

（ 1. 蘭州交通大學 測繪與地理信息學院, 蘭州 730070；2. 地理國情監(jiān)測技術應用國家地方聯(lián)合工程研究中心, 蘭州 730070；3. 甘肅省地理國情監(jiān)測工程實驗室, 蘭州 730070 ）

0 引 言

目前，中國大陸構(gòu)造環(huán)境監(jiān)測網(wǎng)絡(CMONOC)積累了大量連續(xù)運行參考站(CORS)的原始觀測數(shù)據(jù)，為研究青海地區(qū)的板塊運動、地殼形變、速度場模型等提供了十分重要的基礎數(shù)據(jù)[1]. 由于觀測數(shù)據(jù)中存在的有色噪聲(CN)降低了測站運動參數(shù)的估計精度，且不同地區(qū)的噪聲特性差異明顯[2-7]，因此，選擇合適的噪聲模型來研究連續(xù)站的穩(wěn)定性及地表規(guī)律具有重要意義. ZHANG等[5]在加利福尼亞南部地區(qū)時間序列研究中引入噪聲分析，發(fā)現(xiàn)時間序列噪聲具有白噪聲+閃爍噪聲(WN+FN)的特征. 蔣志浩等[6]的研究表明，CORS站坐標時間序列具有WN、FN和隨機漫步噪聲(RW)的噪聲特征. 管雅慧等[7]對云南地區(qū)CMONOC連續(xù)站坐標時間序列進行分析，認為時間序列觀測數(shù)據(jù)超過5年時，噪聲模型才能逐漸趨于穩(wěn)定. 王偉等[8]基于CMONOC觀測數(shù)據(jù)獲得了青藏高原的水平運動速度場及應變場，認為該地區(qū)應變場的整體分布特征與該地區(qū)活動構(gòu)造的空間分布及地震活動十分一致. 綜合上述研究表明，應變場分析應當選擇較長時間段的觀測數(shù)據(jù)并顧及CN的影響.

針對目前青海地區(qū)應變場研究中未顧及CN的影響且未對較長時間段(10年或以上)坐標時間序列噪聲進行分析，本文通過分析青海地區(qū)14個CMONOC連續(xù)站2010—2020年的時間序列數(shù)據(jù)，分別確定了3個方向的最優(yōu)噪聲模型，得到國際地球參考框架(ITRF)下經(jīng)最優(yōu)噪聲模型改正后的青海速度場，并基于修正后的速度場建立整體旋轉(zhuǎn)與線性應變模型來分析青海地區(qū)的應變特征.

1 時間序列噪聲分析方法

目前大部分測繪工作者主要使用時間序列分析軟件CATS和HECToR等軟件對不同區(qū)域的GPS觀測數(shù)據(jù)進行時間序列噪聲分析. 其中，CATS軟件使用最大似然估計法(MLE)且具有較為精確的算法和模型. 但是，在處理時間跨度較大并且噪聲類型較為復雜的時間序列數(shù)據(jù)時，CATS軟件的處理速度十分緩慢. 為更高效處理含有大量數(shù)據(jù)的時間序列，BOS等[9]研發(fā)了HECToR軟件，該軟件通過算法改進提高了數(shù)據(jù)處理的速度，很好地克服了CATS軟件的弊端. 改進后的MLE即為貝葉斯信息準則(BIC).

1.1 MLE

目前坐標時間序列分析普遍使用MLE，該方法可以得到GPS坐標時間序列中所涉及到的絕大部分參數(shù). 坐標時間序列表示為

式中：ti為歷元，以a為單位；a為初始位置；b為線性速率；c、d、e、f分別為周期項的系數(shù)；H為海維西特階梯函數(shù)；gj為同震偏移；vi為殘差.

1.2 BIC

Hector軟件使用BIC來評價擬合結(jié)果的優(yōu)良性，該軟件能夠?qū)r間序列的線性趨勢項、高階多項式、季節(jié)項、其他周期項信號以及多種噪聲模型的組合進行估計[10]. 使用BIC分析噪聲模型的公式為

式中：N為實際觀測數(shù)，r為殘差矩陣；協(xié)方差矩陣C可分解為

其中， σ 為殘差序列的標準差，其計算公式為

detcA=cNdetA

由 可得

于是可得

似然函數(shù)L越小，BIC值越小，噪聲模型相對越好；反之，似然函數(shù)L越大，BIC值越大，噪聲模型相對越差.

2 時間序列分析

2.1 數(shù)據(jù)來源

CMONOC是中國地殼運動觀測網(wǎng)的延續(xù)，在中國境內(nèi)共設置了260個連續(xù)站，青海境內(nèi)共有14個連續(xù)站，選取GAMIT/GLOBK 軟件解算得到的青海地區(qū)2010—2020年的14個CMONOC連續(xù)站為基礎數(shù)據(jù)，測站分布如圖1所示. 本文數(shù)據(jù)來源于中國地震局GNSS數(shù)據(jù)產(chǎn)品服務平臺[11].

圖1 青海地區(qū)CMONOC站點分布圖

此處僅展示QHGE、QHMY和QHTR站的原始坐標時間序列. 如圖2所示，原始坐標時間序列的東(E)、北(N)、天頂(U)方向都存在異常值，必須予以剔除. E方向表現(xiàn)為線性遞增趨勢，斜率較大；N方向表現(xiàn)為線性趨勢，斜率較?。籙方向表現(xiàn)為周期變化趨勢.

圖2 測站原始時間序列圖

2.2 原始時間序列處理

本文利用最小二乘法估計并剔除坐標時間序列中的線性趨勢和周期性趨勢來計算殘差，并使用四分位數(shù)粗差探測方法剔除異常值. 具體步驟包括：1) 采用最小二乘法剔除原始時間序列中的線性趨勢和周期性趨勢，得到時間序列中的殘差；2) 用四分位距探測殘差時間序列包含的異常值，然后剔除；3) 擬合剔除異常值之后的殘差時間序列并重復步驟2)，直到完全剔除[9]. 圖3為處理后的殘差時間序列.

圖3 處理后測站殘差時間序列圖

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

對14個CMONOC連續(xù)站E、N、U方向的時間序列數(shù)據(jù)進行解算，并將WN與FN、PL、GGM及FN+RW進行組合，利用這四種組合模型計算每個測站所對應E、N、U方向的BIC值，對比各個模型所得BIC值即可獲得所有站點各方向的最優(yōu)噪聲模型.表1為各方向噪聲模型所占百分比，由表1可知，各方向具有不同噪聲特征. 在E、N、U方向上，最優(yōu)噪聲模型分別為 WN+PL、WN+GGM和WN+FN.

表1 E、N、U方向噪聲模型百分比 %

3 顧及有色噪聲的青海地區(qū)水平速度場分析

圖4為青海境內(nèi)未經(jīng)最優(yōu)噪聲模型修正的CMONOC連續(xù)站的水平運動速度場. 由圖4可知，青海地區(qū)的整體運動方向近E方向，東部地區(qū)的運動方向近E方向，其余地區(qū)的運動方向為E、N方向.未經(jīng)最優(yōu)噪聲模型修正的青海CMONOC連續(xù)站在ITRF14框架下水平運動方向為88°55′10″NEE，平均速度為39.58 mm/a. E方向的平均速度為39.57 mm/a，標準差為5.09 mm；N方向的平均速度為0.75 mm/a，標準差為3.76 mm.

圖4 修正前的水平速度場

圖5為經(jīng)最優(yōu)噪聲模型改正后的青海CMONOC連續(xù)站的水平運動速度場，修正后的青海CMONOC連續(xù)站在ITRF14框架下水平運動方向為88°57′58″NEE，平均速率為39.45 mm/a. 并且，E方向的平均速度為39.44 mm/a，標準差為5.01 mm；N方向的平均速度為0.71 mm/a，標準差為3.76 mm. 將最優(yōu)噪聲模型修正前后的速度場進行對比分析可知，修正后速度場的精度優(yōu)于修正前速度場的精度，其平均水平運動速率相差0.13 mm/a；水平運動方向相差2′48″NEE.

4 顧及有色噪聲的青海地區(qū)應變特征分析

4.1 應變分析模型

以往研究表明，地殼板塊是彈塑性體或黏彈性體，在周圍板塊的作用下會產(chǎn)生整體平移、旋轉(zhuǎn)和內(nèi)部形變[12-13]. 塊體內(nèi)部應變的實質(zhì)是板塊內(nèi)部質(zhì)點的運動，并且這些質(zhì)點運動的應變參數(shù)往往是與位置相關的線性函數(shù)[14-15]. 顧及板塊的整體運動與線性形變可得塊體的整體旋轉(zhuǎn)與線性應變模型，其形式為[16]

式中：A0～A2、B0～B2、C0～C2為板塊的應變參數(shù)； ωx、ωy、ωz為板塊的歐拉矢量的3個分量；r為地球半徑；λ為質(zhì)點的大地經(jīng)度；φ為質(zhì)點的大地緯度.

4.2 應變場分析

本文利用經(jīng)最優(yōu)噪聲模型改正后的速度場建立青海地區(qū)地殼水平運動的整體旋轉(zhuǎn)與線性應變模型.通過分析計算得到青海地區(qū)內(nèi)部形變的空間分布特征并繪制該區(qū)域的主應變圖、面膨脹圖以及最大剪應變圖. 由圖6可知，青海地區(qū)東北部表現(xiàn)為明顯的擠壓應變特征，而該地區(qū)的西南部主要表現(xiàn)為拉張?zhí)卣?，最大主壓應變和最大主張應變分別出現(xiàn)在青海地區(qū)的東北角與西南角，其值分別可達-4.31×10-8和4.56×10-8. 由圖7可知，青海地區(qū)的東北、西南部為最大剪應變高值區(qū)，對應構(gòu)造活動較強烈. 面膨脹率反映擠壓與拉張作用的相對大小，正、負值分別表示以張性活動或壓性活動為主. 由圖8可知，從東北向西南，擠壓應變逐漸減小，拉張應變逐漸增大，總體表現(xiàn)為擠壓應變.

圖6 主應變圖

圖7 最大剪應變圖

圖8 面膨脹圖

5 結(jié) 論

本文基于最優(yōu)噪聲模型，對青海境內(nèi)14個CMONOC連續(xù)站2010—2020年期間的觀測數(shù)據(jù)進行分析. 首先，利用最優(yōu)噪聲模型得到修正后的青海速度場. 其次，使用整體旋轉(zhuǎn)與線性應變模型對速度場進行插值，并求出應變參數(shù). 最后，根據(jù)所得結(jié)果對青海地區(qū)的應變特征進行分析，結(jié)論如下：

1)青海地區(qū)CMONOC連續(xù)站坐標時間序列在E、N、U方向上具有不同的噪聲特性，其最優(yōu)噪聲模型分別為組合模型PL+WN、GGM+WN和 FN+WN；

2)考慮CN影響后的速度場精度明顯優(yōu)于未經(jīng)修正的速度場精度，修正前后的水平運動方向和平均水平運動速率差值可達2′48″NEE和0.13 mm/a，修正后青海地區(qū)CMONOC連續(xù)站基于ITRF2014框架下水平方向運動的平均速率為39.45 mm/a，運動方向為88°57′58″NEE；

3)青海地區(qū)的東北、西南部構(gòu)造活動較強烈，東北部表現(xiàn)為明顯的擠壓應變特征，而西南部主要表現(xiàn)為拉張?zhí)卣? 從東北往西南，擠壓應變逐漸減小，拉張應變逐漸增大，總體表現(xiàn)為擠壓應變.