基于CATS軟件的GPS時間序列分析

葉錦強

（1.惠州市國土資源局 大亞灣經濟技術開發區分局，廣東 惠州 516081）

GPS已被廣泛運用到地殼形變的研究上，而將GPS資料制成時間序列后，可以分析GPS的連續觀測資料得到測站點與時間及空間的關系，進而推估出各種構造與地殼形變之間的關聯性。

GPS的觀測數據中存在噪聲。數據噪聲可以簡單地分為2種：與時間無關的白噪聲和與時間相關的有色噪聲。白噪聲是一種隨機的信號，在功率頻譜上呈現水平的功率譜密度。白噪聲在光功率頻譜內位于可見光波段，當眼睛接受光的三原色分量相等時為白光。白噪聲之外的噪聲統稱為有色噪聲。

在多數的GPS時間序列的研究中，僅將觀測資料的誤差視為與觀測時間無關的白噪聲，因為白噪聲的數值模型和計算比較容易。但大量研究表明，如果GPS時間序列分析中不考慮與時間相關的有色噪聲，將會低估地殼形變速度的誤差[1-3]。尤其是在一些地殼活動速度比較低的地區測得的GPS位移量很小，若低估誤差可能導致地殼變形的大小及方向上的錯誤，這時候考慮有色噪聲的影響就格外重要。本文中將有色噪聲分為閃爍噪聲和隨機游走噪聲。前者可能來自GPS系統中未被模型化的大氣效應、衛星軌道及星歷誤差、多路徑效應、衛星海拔高度及方位角所造成的效應[3-5]；后者來自于地表潛移或降雨使儀器支柱產生布朗運動所造成的效應[6]。

本文利用時間序列分析軟件CATS[7]，基于兩種不同的方法來進行資料的噪聲分析，分別是時間域上的最大似然估計法和頻率域上的頻譜分析法。

1 GPS資料處理與噪聲分析方法

對于GPS的時間序列的分析，通常考慮長周期運動、同震變形以及震后變形。GPS單站、單分量坐標序列一般用下列周期性模型來表示[8]：

式中，ti以年為單位；系數a、b分別代表地殼位置和線性變換率；系數c、d、e和f描述了周年和半周年運動振幅；系數gj代表地震造成的同震偏移（offset）；H（t）為階躍函數（heaviside step function）；vi為殘差值，代表觀測值與預測值間的差異，也就是本文需要分析的噪聲信號。

由于研究目的是與地殼形變無關的噪聲信息，所以針對GPS時間序列都必須按照式（1）事先扣除已知的信號，包括變形速率、各種原因造成的偏移量、季節性的周期變化（在此為年周期和半年周期），研究區站點分布見圖1。

圖1 研究地區站點分布圖

圖2 SUIY測站時間序列處理事例

圖2a中從上至下分別是南北分量、東西分量、垂直分量的原始時間序列（紅色線）及其擬合線（綠色線），圖2b中從上至下分別是南北分量、東西分量、垂直分量的殘差時間序列，也就是本文所要研究的噪聲時間序列。

1.1 頻譜分析

GPS原始時間序列扣除擬合值后的噪聲時間序列，可以用頻譜定性分析噪聲的類型。頻譜分析是一種在頻率域上分析信號的方法。噪聲時間序列在頻譜域上的功率譜可以使用Power Law的形式表示[6]：

式中，α為頻譜指數；P0為常數。α越大，表示噪聲序列的時間相關性越高。

將式（2）取對數進行曲線擬合，運用最小二乘法可以得到P0及α值。α通常是介于-1～3之間的任意實數，其中整數α代表一些特殊的噪聲類型：α=0時，是標準的白噪聲；α=1時，是標準的閃爍噪聲；α=2時，則是標準的隨機游走噪聲。另外，除了標準的白噪聲以外，其余的部分統稱為有色噪聲。

1.2 最大似然估計法

最大似然法是一種非線性的最小二乘法，目的在于找到與時間序列最相近的模型參數。通過調整協方差矩陣使得似然函數取得最大值，即可得到與該時間序列最相近的噪聲模型，這樣就可以通過協方差矩陣估算出時間序列中的噪聲振幅大小。與頻譜分析相比，最大似然估計法可以定量地計算噪聲大小。

最大似然估計法的方程式如下[1]：

或

式中，ln是自然對數；lik（x，C）是概略值；n是資料點數；det是行列式；C是資料的協方差矩陣。確定協方差矩陣的方法可以參考文獻 [1]、[2]、[9]、[10]。

2 時間序列分析軟件CATS

CATS軟件包含2個程序：線性部分使用最小二乘，利用式（1）計算截距、斜率、偏移以及周期信號的振幅；非線性部分則利用線性計算后的殘差求得特定噪聲模型的參數及振幅大小。

假設模型中只有一種噪聲，則模型的協方差矩陣C=σ2J（σ表示噪聲振幅大小），則式（4）可以改寫為：

此時，可以由式（5）得到噪聲振幅的關系式，此關系式可以令式（5）達到最大值：

當模型中不只一種噪聲模型時，CATS用下面的方法來計算最大似然估計值。設2個噪聲振幅分別σ1和σ2，將振幅以角度φ和一個比例常數r表示：

而協方差矩陣可以表示成：

代入式（4），就可以簡單地運用角度φ計算最大似然估計法的數值。

3 最佳噪聲模型分析

本文使用3個噪聲模型來分析噪聲時間序列：白噪聲（WN）、白噪聲+閃爍噪聲（FL）、白噪聲+隨機游走噪聲（RW），利用CATS軟件，使用最大似然估計法在時間域上求最佳噪聲模型，并使用頻譜分析法從頻率域上找出最佳模型。

表1中針對3種模型計算出各測站噪聲三分量的MLE值。從表中可以統計出，16個GPS測站的48個分量中有90%以上的FL有較大的MLE值（相比RW和WN的MLE值）。只有極少分量（JLYJ的E分量和SUIY的N、E分量）的RW有較大MLE值，這可能是因為觀測或解算中存在誤差，也不排除與該站的基座有關。雖然存在上述異常情況，但整體上可以說明，白噪聲+閃爍噪聲是描述研究地區GPS觀測網噪聲時間序列的較好模型。

表1 GPS觀測網中各測站的MLE值

另外需要注意的是，“陸態網絡”的GPS儀器基座為水泥柱，屬于深基及淺基支柱型的較穩定的儀器基座，相對來說不容易產生隨機游走噪聲，這與計算的結果也是一致的。

利用CATS軟件計算出各測站三分量的有色噪聲頻譜指數，將頻譜分析的模型設定為噪聲加上一個Power Law的有色噪聲，計算結果如表2。從表中可以看出，研究地區的GPS資料的頻譜指數基本介于0～-1之間，與前人的研究結果相同[2,3]，與最大似然估計法的計算結果也相符。

表2 觀測網中各測站的頻譜指數

4 結 語

本文對GPS觀測數據進行噪聲分析。基于最大似然估計法計算結果表明，白噪聲+閃爍噪聲是較近似的噪聲模型。少數幾個噪聲接近白噪聲+隨機游走噪聲模型的測站，可能是由于儀器基座不穩定或者時間序列長度較長所產生的效應。而頻譜分析所得到的結果與最大似然估計法一致。不同的噪聲模型將對地殼形變速度的誤差造成很大的影響，而研究更精確的模型對于反演地球自身的運動狀態尤為重要。

[1]Langbein J,Johnson H.Correlated Errors in Geodetic Timeseries:Implications for Time-dependent Deformation[J].Geophys Res,1997,102(1):591-603

[2]Zhang J,Bock Y,Johnson H,et al.Southern California Permanent GPS Geodetic Array:Error Analysis of Daily Position Estimates and Site Velocities[J].Geophys Res,1997,102(8): 35-55

[3]Mao A,Harrison C G A,Dixon T H.Noise in GPS Coordinate Time Series[J].Geophys Res, 1999,104(1):2 797-2 816

[4]Williams S,Bock Y,Fang P,et al.Error Analysis of Continuous GPS Position Time Series[J].Geophys Res,2004,109(3):412-421

[5]Langbein J.Noise in Two-color Electronic Distance Meter Measurements Revisited[J].Geophys Res,2004,109(4):1-16

[6]Agnew D C.The Time-domain Behavior of Power Law Noises[J].Geophysical Research Letters,1992,19(4):333-336

[7]Williams S D P.The Effect of Coloured Noise on the Uncertainties of Rates Estimated from Geodetic Time Series[J].Geodesy, 2003,76(10):483-494

[8]Nikolaidis R.Observation of Geodetic and Deismic Deformation with the Global Positioning System[D].University of California,2002

[9]黃立人.GPS基準站坐標分量時間序列的噪聲特性分析[J].大地測量與地球動力學,2006,26(2):31-34

[10]Gardner M.White and Brown Music, Fractal Curves and Oneover-fluctuations[J].Scientific American, 1978,238(4):16-32