鄭博文，杜向峰，詹松輝，吳希文，王 華

（1. 廣東工業(yè)大學(xué) 土木與交通工程學(xué)院, 廣東 廣州 510006；2. 廣東工貿(mào)職業(yè)技術(shù)學(xué)院 測(cè)繪遙感信息學(xué)院, 廣東 廣州 510510）

目前，利用GPS位置時(shí)間序列估計(jì)速度已廣泛用于研究地表形變、板塊運(yùn)動(dòng)等[1-2]，然而GPS速度值的精度如何評(píng)定尚需探討[3-5]。由于GPS位置時(shí)間序列中不僅含有隨機(jī)噪聲(白噪聲)，還含有與時(shí)間相關(guān)的噪聲。如果只使用純白噪聲模型，速度的精度會(huì)被嚴(yán)重高估。大量研究表明，僅用純白噪聲模型進(jìn)行精度估計(jì)，速度的精度會(huì)被高估3～38倍[6-13]。因此，必須采用正確的噪聲模型來估計(jì)GPS速度的精度。本文以“中國地殼運(yùn)動(dòng)觀測(cè)網(wǎng)絡(luò)”和“中國大陸構(gòu)造環(huán)境監(jiān)測(cè)網(wǎng)絡(luò)”(CNONOC-I/II)所有GPS連續(xù)站的坐標(biāo)時(shí)間序列為基礎(chǔ)，利用極大似然估計(jì)法和頻譜分析法計(jì)算其光譜指數(shù)，同時(shí)采用白噪聲+冪律噪聲(White Noise+Power Law Noise, WN+PL)、白噪聲+閃爍噪聲(White Noise+Flicker Noise, WN+FN)和白噪聲+閃爍噪聲+隨機(jī)游走噪聲(White Noise+Flicker Noise+Random Walk Noise, WN+FN+RW)等模型計(jì)算了速度的不確定性，并繪制了陸態(tài)網(wǎng)絡(luò)連續(xù)站的速度場(chǎng)。

1 數(shù)據(jù)來源與GPS時(shí)間序列解算

本文使用了CMONOC-I/II共264個(gè)GPS連續(xù)站的數(shù)據(jù)。其中CMONOC-I是從1999年觀測(cè)至今的27個(gè)連續(xù)站點(diǎn)，CMONOC-II于2009年開始觀測(cè)，共有237個(gè)連續(xù)站點(diǎn)[1-2]；時(shí)間跨度在4～9 a的有59個(gè)測(cè)站，10 a及10 a以上的有205個(gè)測(cè)站，如圖1所示。本文使用武漢大學(xué)的PANDA軟件進(jìn)行解算[14]，采用精密單點(diǎn)定位模式獲取所有測(cè)站的單日解。在解算中，先后進(jìn)行了接收機(jī)差分碼偏差改正、地球自轉(zhuǎn)改正和相對(duì)天線相位中心模型改正，具體解算的參數(shù)參見文獻(xiàn)[2]的表1。

圖1 GPS連續(xù)站分布圖Fig.1 Map of GPS continuous station distribution

2 基本原理與方法

基于GPS坐標(biāo)時(shí)間序列，利用加權(quán)最小二乘法和極大似然估計(jì)方法計(jì)算測(cè)站的速度及其不確定性。先對(duì)原始觀測(cè)值進(jìn)行粗差及異常階躍的改正，獲取干凈的時(shí)間序列。隨后通過單位權(quán)陣獲取觀測(cè)值的改正值，根據(jù)改正值和選取的噪聲模型采用極大似然估計(jì)獲取最適合描述觀測(cè)值的協(xié)方差矩陣。最后通過確定的協(xié)方差矩陣計(jì)算出正確的測(cè)站速度及其不確定性。

2.1 粗差及異常階躍檢測(cè)

GPS觀測(cè)值常常由于野外測(cè)量中出現(xiàn)的突發(fā)狀況(如地震、接收機(jī)故障等)出現(xiàn)異常。因此本文使用式(1)所示的中位數(shù)和四分位范圍(Interquartile Range, IQR)統(tǒng)計(jì)[15]來確定粗差，再利用地震目錄中的地震時(shí)刻檢測(cè)時(shí)間序列中由地震引起的階躍，其余未知的階躍使用改進(jìn)的啟發(fā)式分割算法進(jìn)行階躍探測(cè)[16]。圖2為YANC測(cè)站的時(shí)間序列圖，其中綠色為正常點(diǎn)，黑色為粗差點(diǎn)，橙色為探測(cè)出的階躍。

圖2 YANC測(cè)站原始坐標(biāo)序列(綠色為正常點(diǎn)，黑色為粗差點(diǎn)，橙色豎線為探測(cè)出的階躍)Fig.2 Continuous GPS time series of YANC station (The green is the normal point, the black is gross error, and the orange vertical line is the detected step)

式中：v?i為擬合后的殘差，w為窗口大小，n是用于設(shè)置拒絕水平的整數(shù)因子，一般設(shè)置為3，median為中位數(shù)函數(shù)，IQR為四分位范圍統(tǒng)計(jì)函數(shù)。

2.2 加權(quán)最小二乘估計(jì)速度及其精度

本文采用加權(quán)最小二乘進(jìn)行GPS測(cè)站的速度及速度不確定性的計(jì)算，參數(shù)方程如式(2)所示。

式中：ti表示GPS坐標(biāo)時(shí)間序列的單日解歷元，以年為單位；a為測(cè)站起始年份第一天的測(cè)站位置(橫軸截距)；b為線性速度；c和d為年周期振幅；e和f為半年周期振幅；jk為由于各種原因引起的階躍式突變；tjk為發(fā)生突變的歷元；H為赫維賽德階躍函數(shù)(Heaviside step function)，當(dāng)ti小于tjk時(shí)，H值為0，當(dāng)ti大于或等于tjk時(shí)，H值為1；vi為測(cè)量誤差。

假設(shè)參數(shù)矩陣x=[a b c d e f jk]T，GPS時(shí)間序列的觀測(cè)方程為

式中：y為單日解，B為系數(shù)矩陣，v為隨機(jī)誤差vi的矢量。利用加權(quán)最小二乘法(式(4))計(jì)算出每個(gè)測(cè)站的速度、年周期、半年周期及階躍。

2.3 極大似然估計(jì)

采用極大似然估計(jì)方法選取最優(yōu)的噪聲模型協(xié)方差矩陣[17]，見式(5)。

式中：ln為自然對(duì)數(shù)函數(shù)，lik為極大似然函數(shù)，det表示矩陣的行列式，是式(4)中所估計(jì)的最小二乘方殘差，C是所選噪聲模型的協(xié)方差矩陣，N為測(cè)站時(shí)間序列的歷元個(gè)數(shù)。其中C由多個(gè)噪聲模型組合而成，如式(6)所示。

式中：aw為白噪聲振幅大小，I為單位陣，n為噪聲組合個(gè)數(shù)，pi為其余噪聲類型的振幅，Ji是不同噪聲模型的協(xié)方差矩陣，所有的矩陣都為N×N的矩陣。其中，aw、pi為待估計(jì)參數(shù)。閃爍噪聲的近似噪聲模型J-1為

式中：n為行數(shù)；m為列數(shù)。

隨機(jī)游走噪聲的模型J-2為

式中：Δti=ti-t0。

冪律噪聲的噪聲模型J-3[18]可以表示為

其中T的轉(zhuǎn)化矩陣為

式中：ΔT為時(shí)間間隔，φn可以表示為

式中：k為譜指數(shù)，將作為未知參數(shù)與一同估計(jì)。

2.4 頻譜分析

本文利用Lomb-Scargle周期圖法獲得每個(gè)測(cè)站、每個(gè)方向殘差噪聲的功率譜圖，該方法的優(yōu)點(diǎn)是只對(duì)測(cè)量時(shí)間的時(shí)間序列數(shù)據(jù)進(jìn)行評(píng)估，不需要連續(xù)的時(shí)間序列即可計(jì)算出功率譜。如式(12)所示。

式中：f為時(shí)間頻率，t為式(2)中的單日解歷元，yi為式(3)中的單日解，為其算術(shù)平均值，σ2為其方差，N為歷元個(gè)數(shù)，常數(shù)τ 對(duì)應(yīng)的是使周期圖以任意的常數(shù)量獨(dú)立于移位ti的時(shí)滯。由于本文中每個(gè)測(cè)站的時(shí)間序列都足夠長，因此光譜指數(shù)的估計(jì)可以通過簡(jiǎn)單最小二乘擬合來實(shí)現(xiàn)。根據(jù)功率譜P和頻率f，利用式(13)計(jì)算出譜指數(shù)K。

2.5 最優(yōu)噪聲模型評(píng)價(jià)原則

在極大似然估計(jì)中，不同的噪聲模型會(huì)計(jì)算出不同的極大似然估計(jì)值，數(shù)值越大往往越可靠。但由于噪聲模型參數(shù)的數(shù)量也會(huì)對(duì)模型的可靠性造成影響，因此，為了保證結(jié)果的可靠性，本文采用貝葉斯信息(Bayesian Information Criterion, BIC)準(zhǔn)則[19]對(duì)不同模型的BIC值進(jìn)行計(jì)算(見式(14))，依據(jù)BIC值最小的原則選取最優(yōu)模型。

式中：H為參數(shù)個(gè)數(shù)，n為樣本個(gè)數(shù)。

3 結(jié)果

3.1 光譜指數(shù)

本文在剔除GPS坐標(biāo)時(shí)間序列中的粗差、趨勢(shì)項(xiàng)、周期項(xiàng)、半周期項(xiàng)和階躍后，得到GPS噪聲時(shí)間序列。通過利用頻譜分析法，估計(jì)了所有陸態(tài)連續(xù)站ENU 3個(gè)方向的功率譜。圖3為KMIN站的功率譜圖，由圖3可以看出時(shí)間間隔大于14 d(即頻率低于1/14)的噪聲功率是恒定的，此時(shí)主要以白噪聲為主[6]，隨著時(shí)間間隔的增加，噪聲信號(hào)的相對(duì)功率明顯上升。因此本文選取時(shí)間間隔高于14 d的數(shù)據(jù)進(jìn)行簡(jiǎn)單的最小二乘擬合，圖中以紅色直線標(biāo)出，直線的斜率為光譜指數(shù)。圖4(a)為光譜指數(shù)估計(jì)的分布直方圖。由圖可知，東方向光譜指數(shù)的范圍為-1.31～-0.48，北方向?yàn)?1.23～-0.55，垂直方向?yàn)?0.31～-1.244，不同方向上噪聲光譜的特性沒有明顯差異。經(jīng)計(jì)算，東、北、垂直方向的的加權(quán)平均數(shù)分別為-0.82，-0.84，-0.63，所有光譜指數(shù)的加權(quán)平均值為-0.76。已有研究表明[4]，光譜指數(shù)是一個(gè)可以描述噪聲源的指標(biāo)，光譜指數(shù)的范圍為-1＜K＜3，在-1＜K＜1中，K=0時(shí)，噪聲源為白噪聲，K=-1時(shí)為閃爍噪聲，K=-2時(shí)為隨機(jī)游走噪聲。因而中國大陸GPS位置時(shí)間序列主要以白噪聲和閃爍噪聲為主。同時(shí)，從直方圖4(a)可以看出，東、北方向頻次峰值為-0.75，而垂直方向頻次峰值為-0.5。本文認(rèn)為造成該現(xiàn)象的原因是垂直方向的白噪聲過大，掩蓋了閃爍噪聲的特性。圖4(b)為光譜指數(shù)分布直方圖。其中，東方向光譜指數(shù)的范圍為-1.56～-0.58；北方向?yàn)?1.43～-0.70；垂直方向?yàn)?1.40～-0.48。它們的加權(quán)平均數(shù)為-1.02，-1.01，-0.83。從圖4上看，各個(gè)方向沒有明顯的差異，垂直方向的光譜指數(shù)相對(duì)水平方向偏低，但光譜指數(shù)普遍在-1附近。上述兩種結(jié)果表明：中國陸態(tài)連續(xù)站普遍具有閃爍噪聲。同時(shí)，利用極大似然估計(jì)得到的光譜指數(shù)普遍比頻譜分析的要小。在頻譜分析中，光譜指數(shù)可能會(huì)被高估[20]。在本文中也驗(yàn)證了該結(jié)論。

圖3 KMIN站功率譜圖(紅色為利用頻率和功率譜值擬合的一條直線，它的斜率為譜指數(shù))Fig.3 Power spectrum of the KMIN station (Red is a straight line fitted with frequency and power spectrum values, and its slope is the spectral index)

圖4 不同方法所得出的光譜指數(shù)對(duì)比圖Fig.4 Comparison chart of spectral indices obtained by different methods

3.2 速度不確定性分析

本文計(jì)算了792個(gè)GPS坐標(biāo)時(shí)間序列，發(fā)現(xiàn)如果僅使用白噪聲模型進(jìn)行速度不確定性的估計(jì)，速度不確定性會(huì)被嚴(yán)重低估。圖5為其中10個(gè)測(cè)站使用不同的噪聲類型所估計(jì)的速度不確定性。從圖上可以看出利用白噪聲模型估計(jì)速度的不確定性均低于0.05 mm/a，在白噪聲模型中加入閃爍噪聲模型、隨機(jī)游走噪聲模型以及冪律噪聲模型后，速度不確定性會(huì)擴(kuò)大10倍以上。因此不同的噪聲模型會(huì)對(duì)測(cè)站速度以及速度不確定性造成影響，選擇合適的噪聲模型尤為重要。本文選擇了貝葉斯信息準(zhǔn)則即式(13)計(jì)算每個(gè)模型的極大似然估計(jì)值，以此選取每個(gè)測(cè)站的最優(yōu)噪聲模型。圖6給出根據(jù)貝葉斯信息準(zhǔn)則選取的264個(gè)連續(xù)站共792個(gè)分量的最優(yōu)噪聲模型分布情況。整體而言，所有的噪聲分量以WN+FN為主。在水平方向上，東方向和北方向WN+FN噪聲模型分別高達(dá)87%和90%。而在垂直方向上，以WN+PL噪聲模型為主，占比為59%。因此，CMONOC連續(xù)站主要以WN+FN噪聲模型為主。這與田云峰等[7]、李昭等[8]認(rèn)為中國大陸CMONOC連續(xù)站的噪聲模型為白噪聲加閃爍噪聲的結(jié)果一致。本文利用計(jì)算出來的速度及其不確定性，繪制了中國大陸連續(xù)站的速度場(chǎng)(見圖7)，箭頭表示速度，橢圓表示速度不確定性。由圖7可得，部分點(diǎn)(TJWQ(117.130°E, 39.375°N,HETS(118.295°E, 39.736°N等)的不確定性遠(yuǎn)大于其他點(diǎn)，這是由于通過貝葉斯信息準(zhǔn)則判斷出這些點(diǎn)更加符合WN+FN+RW噪聲模型，不確定性會(huì)被擴(kuò)大。當(dāng)測(cè)站的GPS時(shí)間序列具有隨機(jī)游走噪聲特性時(shí)，測(cè)站石墩的穩(wěn)定性可能較弱[20]。因此根據(jù)速度不確定性的解算結(jié)果，這些測(cè)站石墩可能存在不穩(wěn)定的情況。

圖5 不同噪聲模型所得出的速度不確定性Fig.5 Velocity accuracy derived from different noise models

圖 6 噪聲模型不同分量占比Fig.6 Proportion of different components of noise model

圖7 中國大陸連續(xù)站速度場(chǎng)及其速度不確定性Fig.7 Velocity field and velocity accuracy of continents in Chinese mainland

3.3 白噪聲振幅大小與緯度關(guān)系

文獻(xiàn)[20]利用23個(gè)全球連續(xù)站發(fā)現(xiàn)白噪聲振幅在赤道上最大。因此，本文對(duì)空間范圍緯度為15°～50°，經(jīng)度為75°～135°區(qū)域的白噪聲振幅進(jìn)行空間分析。圖8顯示了東、北和垂直方向上白噪聲振幅與站點(diǎn)緯度的關(guān)系圖。圖中各個(gè)方向的關(guān)系圖存在個(gè)別離散點(diǎn)離擬合直線較遠(yuǎn)，這是由于部分GPS測(cè)站受外界因素影響(天線、接收機(jī)故障等)而白噪聲增大。在圖中能夠看出垂直方向的白噪聲振幅是水平方向上的2～3倍。另外3個(gè)方向的白噪聲振幅都會(huì)隨著緯度的下降而隨之增大，這表明白噪聲具有緯度依賴性。

圖8 各個(gè)方向白噪聲振幅與緯度關(guān)系圖Fig.8 The relationship between white noise amplitude and latitude in all directions

4 討論

本文對(duì)利用純白噪聲模型與最佳噪聲模型所計(jì)算的速度不確定性進(jìn)行了對(duì)比，最優(yōu)噪聲模型計(jì)算的速度不確定性是只使用白噪聲模型的6～43倍。具體如表1所示。同時(shí)雖然利用最大似然估計(jì)的方法計(jì)算的速度不確定性是準(zhǔn)確的，但計(jì)算效率較低。因此本文與在GLOBK軟件中利用FOGMEx模型計(jì)算速度不確定性的方法也進(jìn)行了比較。該軟件通過縮放協(xié)方差矩陣的方式來近似估計(jì)速度不確定性[21]，無需利用極大似然估計(jì)方法，具有速度快的特點(diǎn)。圖9是白噪聲模型、選取的最佳噪聲模型和GLOBK軟件利用FOGMEx模型所計(jì)算的速度及其不確定性對(duì)比圖。在圖9可看出東、北方向的白噪聲模型所計(jì)算的速度不確定性均低于0.1 mm/a,垂直方向的速度不確定性少數(shù)在0.1～0.2 mm/a之間，利用最佳噪聲模型所計(jì)算的速度不確定性比利用FOGMEx模型所計(jì)算的略大而垂直方向上具有明顯差異。對(duì)于測(cè)站速度，在東、北方向上兩種方法所計(jì)算的速度差異不大，94%的測(cè)站速度差值的絕對(duì)值小于0.2 mm/a，影響不大。在垂直方向中，可以看出利用本文的最佳噪聲模型所計(jì)算的速度不確定性普遍高于利用FOGMEx模型所計(jì)算的速度不確定性，具有明顯差異，并且22%的測(cè)站會(huì)存在速度差值過大的情況，具體差異如表1所示。根據(jù)表1可以看出在東、北方向上，速度及其速度不確定性的差異平均值都為0.1 mm/a左右，不具有明顯差異。而垂直方向上，速度的不確定性差異平均值為(0.53±1.09) mm/a，具有明顯差異。本文也對(duì)兩者的計(jì)算結(jié)果統(tǒng)計(jì)了倍數(shù)關(guān)系，如表1所示，可以發(fā)現(xiàn)兩者并不具有明顯的倍數(shù)關(guān)系，不能利用簡(jiǎn)單的相乘系數(shù)來獲取精確的速度不確定性。由此可以得出：對(duì)于毫米級(jí)地表形變和板塊運(yùn)動(dòng)研究來說，在東、北方向上，兩種方式所計(jì)算的速度不確定性不具有明顯的倍數(shù)關(guān)系來獲取準(zhǔn)確的速度不確定性，但是兩者方法的結(jié)果不具有明顯差異，可以利用GLOBK軟件中的FOGMEx模型較為快速地獲取測(cè)站的速度不確定性。而垂直方向上兩者存在差異，利用GLOBK軟件無法獲得準(zhǔn)確的速度不確定性。

圖9 WN、所選最佳噪聲模型和FOGMEx模型速度以及速度不確定性對(duì)比Fig.9 Comparison of WN, Best Optimum Model and FOGMEx model velocities and velocity accuracy

表1 白噪聲、最佳噪聲模型與FOGMEx模型對(duì)比統(tǒng)計(jì)Table 1 The Optimum noise model is compared with the FOGMEx model

5 結(jié)論

本文利用頻譜分析和極大似然估計(jì)的方法對(duì)CMONOC連續(xù)站的792個(gè)GPS坐標(biāo)時(shí)間序列噪聲特性進(jìn)行了頻譜分析。并且利用極大似然估計(jì)的方法估計(jì)并且判斷出每個(gè)測(cè)站的噪聲類型和速度不確定性。通過實(shí)驗(yàn)分析可以得出以下結(jié)論：

(1) 在頻譜分析和極大似然估計(jì)方法中，頻譜分析所計(jì)算出的光譜指數(shù)普遍要比極大似然估計(jì)方法得出來的光譜指數(shù)大。但兩種結(jié)果均表明CMONOC連續(xù)站的噪聲主要為白噪聲和閃爍噪聲。

(2) CMONOC連續(xù)站存在噪聲多樣性，并且不同方向上時(shí)間序列，其噪聲類型也有所不用。其中WN+FN噪聲模型占所有GPS坐標(biāo)時(shí)間序列的73%。7%的GPS時(shí)間序列屬于WN+FN+RW噪聲模型，20%的GPS時(shí)間序列屬于白噪聲加冪律噪聲模型。

(3) 通過噪聲分析，發(fā)現(xiàn)有17個(gè)測(cè)站的隨機(jī)游走噪聲振幅偏大，極有可能是觀測(cè)石墩不穩(wěn)定造成的。

(4) CMONOC連續(xù)站GPS坐標(biāo)時(shí)間序列的白噪聲具有緯度依賴性，隨著緯度的降低，白噪聲振幅增大。

(5) 在東、北、垂直方向上，最優(yōu)噪聲模型計(jì)算的速度不確定性分別是只使用白噪聲模型的(11.5±2.1)倍、(12.9±2.9)倍和(14.8±3.7)倍。利用極大似然估計(jì)方法與FOGMEx方法所計(jì)算的速度差距不大，利用極大似然估計(jì)方法所計(jì)算的速度不確定性分別是FOGMEx方法的(2.8±1.5)倍、(1.5±0.7)倍和(3.5±2.8)倍。