張幸，魏冠軍

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

0 引 言

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

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

1 時(shí)間序列噪聲分析方法

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

1.1 MLE

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

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

1.2 BIC

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

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

其中， σ 為殘差序列的標(biāo)準(zhǔn)差，其計(jì)算公式為

detcA=cNdetA

由 可得

于是可得

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

2 時(shí)間序列分析

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

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

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

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

圖2 測站原始時(shí)間序列圖

2.2 原始時(shí)間序列處理

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

圖3 處理后測站殘差時(shí)間序列圖

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

對(duì)14個(gè)CMONOC連續(xù)站E、N、U方向的時(shí)間序列數(shù)據(jù)進(jìn)行解算，并將WN與FN、PL、GGM及FN+RW進(jìn)行組合，利用這四種組合模型計(jì)算每個(gè)測站所對(duì)應(yīng)E、N、U方向的BIC值，對(duì)比各個(gè)模型所得BIC值即可獲得所有站點(diǎn)各方向的最優(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ù)站的水平運(yùn)動(dòng)速度場. 由圖4可知，青海地區(qū)的整體運(yùn)動(dòng)方向近E方向，東部地區(qū)的運(yùn)動(dòng)方向近E方向，其余地區(qū)的運(yùn)動(dòng)方向?yàn)镋、N方向.未經(jīng)最優(yōu)噪聲模型修正的青海CMONOC連續(xù)站在ITRF14框架下水平運(yùn)動(dòng)方向?yàn)?8°55′10″NEE，平均速度為39.58 mm/a. E方向的平均速度為39.57 mm/a，標(biāo)準(zhǔn)差為5.09 mm；N方向的平均速度為0.75 mm/a，標(biāo)準(zhǔn)差為3.76 mm.

圖4 修正前的水平速度場

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

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

4.1 應(yīng)變分析模型

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

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

4.2 應(yīng)變場分析

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

圖6 主應(yīng)變圖

圖7 最大剪應(yīng)變圖

圖8 面膨脹圖

5 結(jié) 論

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

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

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

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