成 誠(chéng)，穆慧敏，平 旗，宋志英，李 云，楊林森

(1.山西省地震局臨汾地震監(jiān)測(cè)中心站，山西 臨汾 041000；2.山西省地震局，山西 太原 030021；3.太原大陸裂谷動(dòng)力學(xué)國(guó)家野外科學(xué)觀測(cè)研究站，山西 太原 030025；4.山西省懷仁市地震局，山西 懷仁 038300)

0 引言

“中國(guó)大陸構(gòu)造環(huán)境監(jiān)測(cè)網(wǎng)絡(luò)”(以下簡(jiǎn)稱陸態(tài)網(wǎng)絡(luò))項(xiàng)目是由中國(guó)地震局牽頭，多部委共同建設(shè)的國(guó)家重大科技基礎(chǔ)設(shè)施項(xiàng)目。2012年3月通過(guò)國(guó)家驗(yàn)收至今，陸態(tài)網(wǎng)絡(luò)運(yùn)行穩(wěn)定，其建設(shè)經(jīng)驗(yàn)已推廣到海外[1-2]。陸態(tài)網(wǎng)絡(luò)為在中國(guó)大陸及其周邊地區(qū)開展大范圍定量化地殼形變監(jiān)測(cè)提供第一手?jǐn)?shù)據(jù)，也為我國(guó)大地測(cè)量基準(zhǔn)體系的建立、維持和自主衛(wèi)星導(dǎo)航定位系統(tǒng)的建立奠定基礎(chǔ)。此外，其還在探究地震孕育機(jī)理及地震預(yù)報(bào)研究中發(fā)揮出作用[3]。

山西地區(qū)共建設(shè)9個(gè)陸態(tài)網(wǎng)絡(luò)連續(xù)基準(zhǔn)站，積累多年觀測(cè)數(shù)據(jù)，使用專業(yè)處理軟件計(jì)算得到山西基準(zhǔn)站高精度時(shí)間序列。時(shí)間序列中存在的噪聲，分為與時(shí)間無(wú)關(guān)的白色噪聲和與時(shí)間相關(guān)的有色噪聲。白色噪聲是一種隨機(jī)信號(hào)，在功率譜上呈水平的功率譜密度。白色噪聲之外的噪聲統(tǒng)稱有色噪聲。Nilolaidis認(rèn)為時(shí)間序列的噪聲更接近于白色噪聲和閃爍噪聲的疊加。因忽略了這一點(diǎn)，速率估計(jì)的結(jié)果會(huì)存在一定偏差，尤其在地質(zhì)條件和板塊運(yùn)動(dòng)較為穩(wěn)定的地區(qū)。由于位移量小，若低估觀測(cè)噪聲的誤差影響，會(huì)導(dǎo)致地殼形變大小甚至方向上出現(xiàn)錯(cuò)誤[4]。該文使用Simon D.P.Williams研發(fā)時(shí)間序列分析軟件CATS，對(duì)山西省GNSS基準(zhǔn)站時(shí)間序列進(jìn)行頻譜計(jì)算和最大似然估計(jì)值計(jì)算，并分析其噪聲特征。

1 山西GNSS網(wǎng)時(shí)間序列與噪聲分析方法

1.1 時(shí)間序列函數(shù)模型

時(shí)間序列中包含構(gòu)造和非構(gòu)造形變信息，非構(gòu)造形變周期性成分不會(huì)對(duì)提取大尺度的構(gòu)造運(yùn)動(dòng)產(chǎn)生影響,會(huì)造成參考框架存在周期性的變化[5]。該文使用中國(guó)地震臺(tái)網(wǎng)中心發(fā)布的山西省GNSS基準(zhǔn)站時(shí)間序列數(shù)據(jù)結(jié)果。

文獻(xiàn)[6]的研究結(jié)果表明，從山西部分GNSS基準(zhǔn)站時(shí)間序列圖發(fā)現(xiàn)，N、E、U三個(gè)方向變化趨勢(shì)相同，N、E方向時(shí)間序列呈線性變化，N方向向南運(yùn)動(dòng)，E方向向東運(yùn)動(dòng)，U方向時(shí)間序列呈周期性變化。對(duì)于GNSS時(shí)間序列的分析，通常要考慮長(zhǎng)周期、同震變形及震后變形，一般用下列周期性模型表示[7]：

(1)

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

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

1.2 頻譜分析法

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

(2)

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

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

1.3 最大似然估計(jì)法

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

(3)

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

2 噪聲計(jì)算結(jié)果與特征分析

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

2.1 頻譜指數(shù)分析結(jié)果

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

表1 山西GNSS基準(zhǔn)站功率譜計(jì)算結(jié)果

2.2 計(jì)算各基準(zhǔn)站最大似然估計(jì)值

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

2.3 各基準(zhǔn)站最佳噪聲模型分析

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

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

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

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

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

表2 山西GNSS基準(zhǔn)站最大似然估計(jì)值

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

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

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

表3 基準(zhǔn)站坐標(biāo)分量時(shí)間序列噪聲分量估值

3 結(jié)論與討論

運(yùn)用CATS軟件里譜指數(shù)方法、最大似然估計(jì)方法分析山西GNSS基準(zhǔn)站時(shí)間序列中的噪聲特征，得到以下結(jié)論：

(1)應(yīng)用CATS軟件計(jì)算山西省基準(zhǔn)站坐標(biāo)時(shí)間序列的譜指數(shù)，結(jié)果表明，基準(zhǔn)站N、E、U方向時(shí)間序列的噪聲分量均具有明顯的負(fù)譜指數(shù)，譜指數(shù)值介于-2～0之間，基準(zhǔn)站坐標(biāo)、時(shí)間序列不僅包括白色噪聲，也包括有色噪聲。

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

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

感謝中國(guó)地震臺(tái)網(wǎng)中心高級(jí)工程師李瑜對(duì)該文的指導(dǎo)，文章所用數(shù)據(jù)源于中國(guó)地震臺(tái)網(wǎng)中心和中國(guó)地震局GNSS數(shù)據(jù)產(chǎn)品服務(wù)平臺(tái)。