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

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

0 引言

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

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

1 極大似然估計法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖1 GSLZ站時間序列

圖2 CSTS站時間序列

3 噪聲模型的建立及分析

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

由MLE的原理可知，不同噪聲模型下GPS時間序列的MLE值不同，且MLE值的大小反映了所得結(jié)果的可靠性.一般情況下，MLE值越大，噪聲模型越有效[16].對于不同的噪聲模型，蒙卡特羅實驗表明：95%的顯著水平下，當兩種噪聲模型的MLE之差大于3.0時，兩種模型具有可區(qū)分性[17].針對甘肅省境內(nèi)的所有CMONOC站點，剔除時間序列中的奇異值，利用CATS軟件計算并統(tǒng)計WN+FN、WN+RWN、WN+FN+RWN、WN+GM、WN+RWN+GM、WN+PL六種噪聲模型的平均MLE值與WN模型的平均MLE值之間的差值，結(jié)果均大于3.0，表明甘肅CMONOC站時間序列中同時含有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站坐標時間序列的噪聲特性.

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

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

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

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

表2 兩種噪聲模型下坐標時序參數(shù)估計比較

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

5 顧及有色噪聲的甘肅省水平速度場分析

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

將ITRF14框架下的水平速度場轉(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)定歐亞框架下的速度場，如圖4所示.結(jié)果與中國地震局GNSS數(shù)據(jù)產(chǎn)品服務平臺（http://www.cgps.ac.cn）結(jié)果基本一致.表4給出了歐亞框架下甘肅省境內(nèi)CMONOC站的三維速度估值和中誤差統(tǒng)計結(jié)果，可以看出，基準站E方向速度平均值為6.37mm/a，N方向速度平均值為1.21mm/a；水平方向運動的平均速率為6.49mm/a，優(yōu)勢方向為NEE 79.23°.

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

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

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

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

6 結(jié)束語

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

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

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

3）甘肅省CMONOC基準站在ITRF14框架下水平方向運動的平均速率為34.54mm/a，運動方向為SEE 98.07°；相對于歐亞板塊的水平方向運動的平均速率為6.49mm/a，運動方向為NEE 79.23°.

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