劉雙童，徐鑒民，高志鈺，陳晨，楊樹文

(1.甘肅農(nóng)業(yè)職業(yè)技術(shù)學(xué)院，蘭州 730020；2.中國(guó)地震局地質(zhì)研究所地震動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100029；3.中煤航測(cè)遙感集團(tuán)有限公司，西安 710100；4.蘭州交通大學(xué)測(cè)繪與地理信息學(xué)院，蘭州 730070)

0 引言

近年來，隨著GPS坐標(biāo)時(shí)間序列數(shù)據(jù)的不斷積累，對(duì)于GPS時(shí)間序列中噪聲特性的研究越來越深.許多前人的研究結(jié)果發(fā)現(xiàn)，GPS站坐標(biāo)時(shí)間序列中在存在白噪聲（WN）的同時(shí)，也存在其他有色噪聲（CN），例如閃爍噪聲（FN）、隨機(jī)漫步噪聲（RWN）、冪律噪聲（PL）、帶通噪聲（BP）、一階高斯馬爾科夫噪聲（FOGM）等，這些噪聲對(duì)GPS時(shí)間序列分析具有重要影響[1-6].Zhang等[1]將噪聲分析引入到GPS坐標(biāo)時(shí)間序列分析中，發(fā)現(xiàn)具有“WN＋FN”的噪聲特征.Langbein[2]分析了加利福尼亞南部及內(nèi)華達(dá)南部地區(qū)連續(xù)GPS站的噪聲模型，認(rèn)為該地區(qū)GPS站的最佳噪聲模型可以表示為FN或RWN.黃立人[3]對(duì)中國(guó)大陸的3個(gè)GPS基準(zhǔn)站時(shí)間序列噪聲特性進(jìn)行了分析，結(jié)果表明所有9個(gè)位置分量（N、E、U）的噪聲時(shí)間序列都具有“WN＋FN”的特點(diǎn).蔣志浩等[4]認(rèn)為我國(guó)國(guó)家連續(xù)運(yùn)行參考站（CORS）網(wǎng)基準(zhǔn)站坐標(biāo)時(shí)間序列噪聲表現(xiàn)為WN、FN及RWN.李昭等[5]對(duì)中國(guó)區(qū)域11個(gè)國(guó)際GNSS服務(wù)（IGS）基準(zhǔn)站坐標(biāo)時(shí)間序列進(jìn)行地表質(zhì)量負(fù)荷改正后，得到的最優(yōu)噪聲模型，主要以WN＋FN為主，同時(shí)也存在PL、BP及GM.陳晨等[6]對(duì)香港衛(wèi)星參考站網(wǎng)坐標(biāo)時(shí)間序列進(jìn)行分析，認(rèn)為WN＋FN是主要噪聲模型，僅考慮WN而忽視CN的做法會(huì)大大低估參數(shù)估計(jì)的不確定度.由此可見，不同地區(qū)GPS坐標(biāo)時(shí)間序列噪聲特性存在著明顯的差異.此外，諸多學(xué)者利用GPS技術(shù)對(duì)全球或區(qū)域地殼運(yùn)動(dòng)進(jìn)行了不斷研究與探索，取得了一些顯著成果[7-12].曾波等[10]、袁鵬等[11]、徐杰等[12]分別建立了山西省、安徽省、山東省三維速度場(chǎng)，這對(duì)于研究該區(qū)域地殼形變具有重要的意義.而針對(duì)甘肅省地區(qū)的速度場(chǎng)方面的研究較少.

針對(duì)上述存在的問題，本文選取甘肅境內(nèi)的19個(gè)CMONOC基準(zhǔn)站坐標(biāo)時(shí)間序列，結(jié)合CATS軟件[13]采用極大似然估計(jì)法（MLE）對(duì)噪聲特性進(jìn)行分析，并首次建立甘肅省地殼水平運(yùn)動(dòng)速度場(chǎng)，為研究甘肅省的地殼運(yùn)動(dòng)規(guī)律提供了相應(yīng)的參考，具有一定的實(shí)際意義.

1 極大似然估計(jì)法

MLE分析坐標(biāo)時(shí)序的噪聲特征時(shí)，對(duì)坐標(biāo)分量時(shí)間序列建立如下參數(shù)模型：

式中：t為GPS單日坐標(biāo)解歷元；a、b分別為測(cè)站初始位置和運(yùn)動(dòng)速度；參數(shù)c、d、e、f分別為年周期和半年周期特征的測(cè)站運(yùn)動(dòng)振幅系數(shù)；H為Heaviside階梯函數(shù)；為跳變改正項(xiàng)，表示發(fā)生在歷元Tgj處振幅為gi的nj個(gè)跳變量；vti表示殘差噪聲序列.

假設(shè)vti由振幅分別為m和n的WN及冪律譜噪聲組成，則有

觀測(cè)值協(xié)方差陣可表示為

式中：I為單位陣；Jα對(duì)應(yīng)譜指數(shù)為α的冪律譜噪聲協(xié)方差陣，表示為

轉(zhuǎn)換矩陣T的表達(dá)式為

式中：φn=[?α/2(1?α/2)···(n?1?α/2)]/n!，n>0；φ0=1，n=0，?T為采樣間隔.當(dāng)數(shù)據(jù)采樣不等間隔時(shí)，對(duì)式（5）中矩陣各列分別乘上相應(yīng)系數(shù)α=?1時(shí)表示FN，α=?2時(shí)表示RWN.

式（1）中，a、b、c、d、e、f、g和噪聲分量振幅m、n為待求參數(shù).按照MLE準(zhǔn)則同時(shí)確定a、b、c、d、e、f、g和噪聲分量振幅m、n，使得坐標(biāo)序列的殘差與其協(xié)方差的聯(lián)合概率密度值達(dá)到最大：

即聯(lián)合概率函數(shù)值的對(duì)數(shù)達(dá)到最大：

式中：N為時(shí)間序列長(zhǎng)度；C為協(xié)方差陣；det表示矩陣C的行列式；lik為似然值.

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

中國(guó)大陸構(gòu)造環(huán)境監(jiān)測(cè)網(wǎng)絡(luò)（簡(jiǎn)稱陸態(tài)網(wǎng)，CMONOC）是中國(guó)地殼運(yùn)動(dòng)觀測(cè)網(wǎng)（簡(jiǎn)稱網(wǎng)絡(luò)工程）的延續(xù)，在中國(guó)境內(nèi)總共建設(shè)有260個(gè)全球衛(wèi)星導(dǎo)航系統(tǒng)（GNSS）基準(zhǔn)站，包含27個(gè)網(wǎng)絡(luò)工程基準(zhǔn)站.在甘肅省境內(nèi)，總共包含19個(gè)站點(diǎn)，站點(diǎn)信息如表1所示.選取GAM IT/GLOBK軟件[14]解算得到的在ITRF14框架下的原始坐標(biāo)時(shí)間序列為基礎(chǔ)數(shù)據(jù)，數(shù)據(jù)處理流程參考文獻(xiàn)[7].同一測(cè)站的最優(yōu)噪聲模型和速度場(chǎng)差與所選時(shí)段有較大的關(guān)系，因此在研究分析測(cè)站最優(yōu)噪聲模型、速度場(chǎng)時(shí)需要指出獲取該時(shí)間序列的時(shí)段[15].本文所選站點(diǎn)對(duì)應(yīng)的時(shí)間段如表1所示.本文所使用的數(shù)據(jù)來源于中國(guó)地震局GNSS數(shù)據(jù)產(chǎn)品服務(wù)平臺(tái)（http://www.cgps.ac.cn）.

表1 甘肅省境內(nèi)CMONOC站點(diǎn)概況

由于篇幅有限，僅繪制GSLZ、GSTS站時(shí)間序列曲線，如圖1、2所示.由圖1、2可以看出，N、E方向分量時(shí)間序列線性趨勢(shì)變化明顯，N方向分量上隨時(shí)間的變化而遞減且變化幅度較小，E方向分量上隨時(shí)間的變化而遞增且變化幅度較大；U方向分量具有一定的周期性變化趨勢(shì)，且以一年周期表征最為明顯；時(shí)間序列中存在明顯的奇異值（或稱為外野值），進(jìn)行數(shù)據(jù)分析前必須剔除.

圖1 GSLZ站時(shí)間序列

圖2 CSTS站時(shí)間序列

3 噪聲模型的建立及分析

GPS坐標(biāo)時(shí)間序列中不僅存在WN，而且存在CN.WN可分為簡(jiǎn)單白噪聲與可變白噪聲（VW）.而VW不能作為參考站最優(yōu)噪聲模型模型，它僅能反映測(cè)站分量的質(zhì)量好壞[5].有色噪聲包含的種類較多，主要選取FN、RWN、PL、GM四種.因此本文選取以下七種噪聲模型，分別為：WN、WN+FN、WN+RWN、WN+FN+RWN、WN+GM、WN+RWN+GM、WN+PL.

由MLE的原理可知，不同噪聲模型下GPS時(shí)間序列的MLE值不同，且MLE值的大小反映了所得結(jié)果的可靠性.一般情況下，MLE值越大，噪聲模型越有效[16].對(duì)于不同的噪聲模型，蒙卡特羅實(shí)驗(yàn)表明：95%的顯著水平下，當(dāng)兩種噪聲模型的MLE之差大于3.0時(shí)，兩種模型具有可區(qū)分性[17].針對(duì)甘肅省境內(nèi)的所有CMONOC站點(diǎn)，剔除時(shí)間序列中的奇異值，利用CATS軟件計(jì)算并統(tǒng)計(jì)WN+FN、WN+RWN、WN+FN+RWN、WN+GM、WN+RWN+GM、WN+PL六種噪聲模型的平均MLE值與WN模型的平均MLE值之間的差值，結(jié)果均大于3.0，表明甘肅CMONOC站時(shí)間序列中同時(shí)含有WN和CN.WN+RWN和WN+GM兩種噪聲模型與WN模型的MLE差值小于其它4種噪聲模型與WN模型的MLE差值，表明WN+FN、WN+FN+RWN、WN+RWN+GM和WN+PL噪聲模型相比WN+RWN和WN+GM噪聲模型更適合表征甘肅CMONOC站坐標(biāo)時(shí)間序列的噪聲特性.

由于噪聲模型中包含的未知參數(shù)會(huì)直接影響到MLE的估值.因此為了確保結(jié)果的可靠性，不能只以MLE值的大小來判定最優(yōu)噪聲模型.采用Langbein保守估計(jì)準(zhǔn)則[17]來選擇19個(gè)基準(zhǔn)站各分量方向上的最優(yōu)噪聲模型，相應(yīng)的實(shí)驗(yàn)結(jié)果表明：在57個(gè)分量方向上，作為最優(yōu)噪聲模型的WN+FN組合模型所占比例最大，為68.42%。因此，WN+FN組合模型可作為甘肅省CMONOC基準(zhǔn)站最優(yōu)噪聲模型.

4 顧及有色噪聲對(duì)參數(shù)估計(jì)的影響

通過CATS軟件估計(jì)出的年周期項(xiàng)系數(shù)a、b及其相應(yīng)的誤差項(xiàng)σa、σb，根據(jù)公式（8）～（11）可計(jì)算振幅A、相位φ及其相應(yīng)的誤差項(xiàng)σA和σφ.

為探討CN對(duì)參數(shù)估計(jì)的影響，本文僅列舉GSAX、GSLZ、GSML、GSTS四個(gè)基準(zhǔn)站W(wǎng)N+FN模型與WN模型估計(jì)的線性速率、年周期振幅和相位三個(gè)參數(shù)比較結(jié)果，如表2所示.

表2 兩種噪聲模型下坐標(biāo)時(shí)序參數(shù)估計(jì)比較

由表2可以看出，4個(gè)站點(diǎn)WN與WN+FN兩種噪聲模型所估計(jì)的參數(shù)差異明顯.WN與WN+FN兩種噪聲模型所估計(jì)得到的線性速度最大差異為GSTS站的U方向（0.56 mm/a）；年周期振幅最大差異的為GSTS站的N方向（0.15mm）；相位最大差異的為GSLZ站的N方向（2.95 rad）.WN+FN模型所得到的線性速率、年周期振幅和相位三個(gè)參數(shù)估計(jì)的最大不確定度是WN模型的14.5倍、6.8倍和4.7倍.結(jié)合其余15個(gè)站兩種模型所估計(jì)參數(shù)結(jié)果可以看出，考慮CN獲得的速度不確定度大約是僅考慮WN的4～15倍，僅考慮WN所獲得的速度估計(jì)量的中誤差不能反映實(shí)際精度.因此，在進(jìn)行坐標(biāo)時(shí)間序列分析時(shí)，CN對(duì)參數(shù)估計(jì)的影響是不容忽視的.

5 顧及有色噪聲的甘肅省水平速度場(chǎng)分析

就水平速度場(chǎng)而言，顧及有色噪聲獲得的速度不確定度大約是僅考慮WN的4～11倍.繪制估計(jì)有色噪聲的ITRF14框架下甘肅省水平運(yùn)動(dòng)速度場(chǎng)，如圖3所示.結(jié)果與中國(guó)地震局GNSS數(shù)據(jù)產(chǎn)品服務(wù)平臺(tái)（http://www.cgps.ac.cn）結(jié)果基本一致.甘肅省東部運(yùn)動(dòng)方向近E向，而西部地區(qū)為E、N向運(yùn)動(dòng).表3給出了ITRF14框架下甘肅省境內(nèi)CMONOC站的三維速度估值和中誤差統(tǒng)計(jì)結(jié)果.由表3可知，本文獲取的甘肅省境內(nèi)CMONOC基準(zhǔn)站在ITRF2014框架下的水平速度精度較高，E方向速度的標(biāo)準(zhǔn)差優(yōu)于0.39 mm，N方向標(biāo)準(zhǔn)差優(yōu)于0.24mm.也可以看出，基準(zhǔn)站E方向速度平均值為34.20mm/a，N方向速度平均值為?4.85 mm/a；水平方向運(yùn)動(dòng)的平均速率為34.54mm/a，優(yōu)勢(shì)方向?yàn)镾EE 98.07°.

將ITRF14框架下的水平速度場(chǎng)轉(zhuǎn)換為ITRF08框架下，然后根據(jù)Wang等[9]給出的歐拉矢量（wx=?0.087 mas/a，wy=?0.514 mas/a和wz=0.741 mas/a），將ITRF08框架下的速度轉(zhuǎn)換為穩(wěn)定歐亞框架下的速度場(chǎng)，如圖4所示.結(jié)果與中國(guó)地震局GNSS數(shù)據(jù)產(chǎn)品服務(wù)平臺(tái)（http://www.cgps.ac.cn）結(jié)果基本一致.表4給出了歐亞框架下甘肅省境內(nèi)CMONOC站的三維速度估值和中誤差統(tǒng)計(jì)結(jié)果，可以看出，基準(zhǔn)站E方向速度平均值為6.37mm/a，N方向速度平均值為1.21mm/a；水平方向運(yùn)動(dòng)的平均速率為6.49mm/a，優(yōu)勢(shì)方向?yàn)镹EE 79.23°.

圖3 ITRF14框架下甘肅省境內(nèi)CMONOC站點(diǎn)水平速度場(chǎng)

圖4 歐亞框架下甘肅省境內(nèi)CMONOC站點(diǎn)水平速度場(chǎng)

表3 ITRF14框架下甘肅省境內(nèi)CMONOC站的水平速度估值和中誤差統(tǒng)計(jì)表(mm·a?1)

表4 歐亞框架下甘肅省境內(nèi)CMONOC站的水平速度估值和中誤差統(tǒng)計(jì)表(mm·a?1)

6 結(jié)束語

本文選取中國(guó)地震局GNSS數(shù)據(jù)產(chǎn)品服務(wù)平臺(tái)（http://www.cgps.ac.cn）提供的由GAM IT/GLOBK軟件結(jié)算得到的甘肅省境內(nèi)19個(gè)CMONOC基準(zhǔn)站坐標(biāo)時(shí)間序列結(jié)果，采用MLE探討了各基準(zhǔn)站坐標(biāo)分量對(duì)應(yīng)的最優(yōu)噪聲模型，分析了顧及有色噪聲對(duì)各參數(shù)估計(jì)的影響，并確定了各基準(zhǔn)站的速度場(chǎng)，分析了甘肅省現(xiàn)今地殼運(yùn)動(dòng)情況.得出以下結(jié)論：

1）甘肅省境內(nèi)CMONOC基準(zhǔn)站各坐標(biāo)分量噪聲特性存在較大的差異，其中“WN+FN”能夠更好地描述基準(zhǔn)站坐標(biāo)時(shí)間序列3分量上的噪聲特性；

2）估計(jì)CN估算得到的速度不確定度是僅考慮WN時(shí)的4～15倍，僅考慮WN所獲得的中誤差不能反映實(shí)際精度；

3）甘肅省CMONOC基準(zhǔn)站在ITRF14框架下水平方向運(yùn)動(dòng)的平均速率為34.54mm/a，運(yùn)動(dòng)方向?yàn)镾EE 98.07°；相對(duì)于歐亞板塊的水平方向運(yùn)動(dòng)的平均速率為6.49mm/a，運(yùn)動(dòng)方向?yàn)镹EE 79.23°.

致謝：感謝中國(guó)地震局GNSS數(shù)據(jù)產(chǎn)品服務(wù)平臺(tái)（http://www.cgps.ac.cn）提供數(shù)據(jù)支撐.部分圖是使用Generic Mapping Tools（GMT）軟件[18]（Wessel et al.,2019）繪制的.