山西省GNSS基準站時間序列噪聲特征分析

成 誠，穆慧敏，平 旗，宋志英，李 云，楊林森

(1.山西省地震局臨汾地震監測中心站，山西 臨汾 041000；2.山西省地震局，山西 太原 030021；3.太原大陸裂谷動力學國家野外科學觀測研究站，山西 太原 030025；4.山西省懷仁市地震局，山西 懷仁 038300)

0 引言

“中國大陸構造環境監測網絡”(以下簡稱陸態網絡)項目是由中國地震局牽頭，多部委共同建設的國家重大科技基礎設施項目。2012年3月通過國家驗收至今，陸態網絡運行穩定，其建設經驗已推廣到海外[1-2]。陸態網絡為在中國大陸及其周邊地區開展大范圍定量化地殼形變監測提供第一手數據，也為我國大地測量基準體系的建立、維持和自主衛星導航定位系統的建立奠定基礎。此外，其還在探究地震孕育機理及地震預報研究中發揮出作用[3]。

山西地區共建設9個陸態網絡連續基準站，積累多年觀測數據，使用專業處理軟件計算得到山西基準站高精度時間序列。時間序列中存在的噪聲，分為與時間無關的白色噪聲和與時間相關的有色噪聲。白色噪聲是一種隨機信號，在功率譜上呈水平的功率譜密度。白色噪聲之外的噪聲統稱有色噪聲。Nilolaidis認為時間序列的噪聲更接近于白色噪聲和閃爍噪聲的疊加。因忽略了這一點，速率估計的結果會存在一定偏差，尤其在地質條件和板塊運動較為穩定的地區。由于位移量小，若低估觀測噪聲的誤差影響，會導致地殼形變大小甚至方向上出現錯誤[4]。該文使用Simon D.P.Williams研發時間序列分析軟件CATS，對山西省GNSS基準站時間序列進行頻譜計算和最大似然估計值計算，并分析其噪聲特征。

1 山西GNSS網時間序列與噪聲分析方法

1.1 時間序列函數模型

時間序列中包含構造和非構造形變信息，非構造形變周期性成分不會對提取大尺度的構造運動產生影響,會造成參考框架存在周期性的變化[5]。該文使用中國地震臺網中心發布的山西省GNSS基準站時間序列數據結果。

文獻[6]的研究結果表明，從山西部分GNSS基準站時間序列圖發現，N、E、U三個方向變化趨勢相同，N、E方向時間序列呈線性變化，N方向向南運動，E方向向東運動，U方向時間序列呈周期性變化。對于GNSS時間序列的分析，通常要考慮長周期、同震變形及震后變形，一般用下列周期性模型表示[7]：

(1)

式中：ti以年為單位；系數a、b分別表示地殼位置和線性變化率；系數c、d、e和f描述周年和半周年運動振幅；系數gj表示地震造成的同震偏移；Ht為階躍函數；vi為殘差值即下文所指的噪聲信號。應用該參考模型可有效估計基準站時間序列的季節性變化，同時，也能反映連續基準站的位置和速度信息。

由于研究與地殼形變無關的噪聲信息，故對基準站時間序列按式(1)扣除線性項、周期項及各種原因造成的偏移量。采用Williams研發的時間序列分析軟件CATS擬合式(1)中各系數，扣除時間序列的已知項，得到殘差序列，對噪聲特性進行分析。

1.2 頻譜分析法

原始時間序列扣除已知項的噪聲時間序列，可用頻譜定性分析噪聲類型。頻譜分析是一種在頻率域上分析信號的方法，可視為一種冪次法則。噪聲時間序列在頻譜域上的功率譜可使用Power Law的形式表示[8]：

(2)

式中：f為時間頻率；P0和f0是正則化常數；k為譜指數。k值越大，表示噪聲序列的時間相關性越高。

將式(2)取對數進行曲線擬合，運用最小二乘法可得到P0及k值。k值通常為-3

1.3 最大似然估計法

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

(3)

式中：detC為矩陣的行列式；C為假定噪聲的協方差陣；N為時間序列的長度；v為線性擬合后的殘差，由C采用加權最小二乘法求得。協方差矩陣C可表達若干隨機噪聲過程，如，白色噪聲(WH)、閃爍噪聲(FN)、隨機游走噪聲(RWN)等及其各類組合[10-11]。

2 噪聲計算結果與特征分析

CATS軟件是由Simon D.P.Williams研發計算時間序列的程序，能得到時間序列中包含線性部分和非線性部分的處理結果。其中，線性部分使用最小二乘法計算截距、斜率、偏移及周期信號的振幅，非線性部分采用線性計算后的殘差求特定噪聲模型的參數和振幅大小[12]。

2.1 頻譜指數分析結果

在自然界噪聲中，不同的譜指數對應不同的噪聲特性。應用CATS軟件頻譜分析模塊計算山西各站扣除已知項時間序列分量的譜指數(見表1)。

表1 山西GNSS基準站功率譜計算結果

2.2 計算各基準站最大似然估計值

由于山西各基準站所處環境不同，產生噪聲的來源可能會有差異，噪聲特性可能也不相同。根據前人研究成果，選取WH、WH+FN、WH+RWN、WH+FN+RW 4種噪聲模型，采用CATS軟件計算9個參考站各類噪聲模型組合的最大似然估計值(見第40頁表2)。

2.3 各基準站最佳噪聲模型分析

譜指數結果可直接揭示山西GNSS基準站時間序列的噪聲特性，同時，還可采用不同噪聲模型下最大似然估計值的比較來研究噪聲特性[13]。山西GNSS基準站坐標時間序列的譜指數結果如表1所示。表中站點為GNSS基準站站名，N、E、U分別為各基準站坐標時間序列的三個方向分量譜指數值。可以看出，山西9個GNSS基準站時間序列各分量的噪聲譜指數都為負值(-2～0)，表明各分量不僅有白色噪聲，也包含有色噪聲。

為分析噪聲類型或噪聲組合，根據式(3)計算得到4種不同噪聲模型組合的最大似然估計值。表2為山西GNSS基準站時間序列各分量在不同噪聲模型的最大似然估計值，數值越大，結果越可靠。僅靠這一指標確定最優模型存在較大片面性，Langbein提出的保守估計方法評價，是以某種噪聲模型的最大似然值為參考值，另外幾種噪聲模型最大似然估計值與參考值作差，差值最大為最優噪聲模型[13]。研究中以WN模型下最大似然估計值為參考值，分別將WN+FN、WN+RWN、WN+FN+RWN 3種模型的最大似然值與之作差，得到的結果如第41頁圖1至圖3所示。

圖1 N方向最大似然值差值

圖2 E方向最大似然值差值

圖3 U方向最大似然值差值

表2 山西GNSS基準站最大似然估計值

可以看出，山西省GNSS基準站時間序列各分量的WH+FN、WH+RWN、WH+FN+RWN與WH噪聲模型的最大似然估計值之差值大于0，表明各分量中存在白色噪聲，同時也存在有色噪聲；WH+FN、WH+FN+RWN的噪聲模型最大似然估計值差值大于WH+RWN，且這兩種模型組合差值彼此極為接近，表明WH+FN、WH+FN+RWN兩種噪聲模型組合優于WH+RWN。

計算結果顯示，山西省GNSS基準站時間序列各分量中，存在WH+FN或WH+FN+RWN噪聲模型。根據蒙特卡羅模擬實驗證明，WH+FN與WH+FN+RWN兩種噪聲模型無可分性[14]。因此，根據前人研究結果，以WH+FN+RWN為最佳噪聲模型組合，驗證噪聲模型中各種噪聲相對應的噪聲分量估值，確定隨機漫步噪聲是否存在[15]。

第42頁表3為山西基準站時間序列WH+FN+RWN噪聲分量數值，可以得出，山西省GNSS基準站時間序列各分量具有不同噪聲特性。其中，SXCZ、SXCH、SXKL、SXLF、SXLQ站時間序列N方向存在WH+FN類型的噪聲,SXDT、SXGX、SXTY、SXXX站時間序列在N方向存在WH+FN+RWN噪聲；山西各基準站時間序列E方向存在WH+FN+RWN噪聲;SXDT的U方向時間序列存在WH+FN+RWN噪聲，其余基準站U方向時間序列包含WH+FN噪聲[16-17]。選用WH+FN+RWN為山西 GNSS基準站的最佳模型并不科學。因此，山西省GNSS基準站N方向的時間序列噪聲模型為WH+FN或WH+FN+RWN，E方向的為WH+FN+RWN，U方向的為WH+FN。

表3 基準站坐標分量時間序列噪聲分量估值

3 結論與討論

運用CATS軟件里譜指數方法、最大似然估計方法分析山西GNSS基準站時間序列中的噪聲特征，得到以下結論：

(1)應用CATS軟件計算山西省基準站坐標時間序列的譜指數，結果表明，基準站N、E、U方向時間序列的噪聲分量均具有明顯的負譜指數，譜指數值介于-2～0之間，基準站坐標、時間序列不僅包括白色噪聲，也包括有色噪聲。

(2)采用CATS軟件WH+FN、WH+RWN、WH+FN+RWN、FN+RWN四種噪聲模型，對山西省基準站時間序列進行最大似然估計值計算，與WN模型的最大似然估計值作差比較，初步判定WN+FN+RWN為最佳噪聲模型。

(3)山西省GNSS基準站時間序列在N、E、U方向上具有不同的噪聲特性。基準站點在N方向上有5個站點包含WN+FN+RWN噪聲，4站點包含WH+FN噪聲；E方向上有9個站點包含WH+FN+RWN噪聲；U方向上大部分站點只包含WH+FN+RWN噪聲，有1個站點包含WH+FN噪聲。所以，山西省GNSS基準站時間序列采用WH+FN+RWN模型作為N、E、U方向的最佳噪聲模型并不合理，其最佳噪聲模型應為：N方向上采用WH+FN+RWN或WH+FN模型，E方向上采用WH+FN模型，U方向上采用WH+FN+RWN模型。可見，三個方向的最佳噪聲模型可以不一致。

感謝中國地震臺網中心高級工程師李瑜對該文的指導，文章所用數據源于中國地震臺網中心和中國地震局GNSS數據產品服務平臺。