半?yún)?shù)模型在GPS季節(jié)性信號(hào)分離中的應(yīng)用

盧辰龍 趙衛(wèi)林 匡翠林 劉鐵柱

1 鄭州市市政工程勘測(cè)設(shè)計(jì)研究院，鄭州市民生路1號(hào)，450052

2 中南大學(xué)地球科學(xué)與信息物理學(xué)院，長(zhǎng)沙市麓山南路932號(hào)，410083

GPS坐標(biāo)序列不僅受構(gòu)造運(yùn)動(dòng)的影響，還受非構(gòu)造運(yùn)動(dòng)的影響，因此GPS坐標(biāo)序列中含有速率及季節(jié)性信號(hào)，尤其是在高程方向上。研究表明，GPS時(shí)間序列不僅存在異常高頻周期性信號(hào)［1］，而且周期信號(hào)的振幅也是時(shí)變的［2－3］。因此，簡(jiǎn)單地使用Nikolaidis［4］提出的模型（下文稱其為傳統(tǒng)模型）去除趨勢(shì)及季節(jié)性信號(hào)，可能會(huì)忽略一些有用的周期信號(hào)，引入模型誤差，從而造成速率的偏差［3，5］。此外，周期信號(hào)還會(huì)影響ITRF實(shí)現(xiàn)的精度［6］。目前，世界上有成千上萬(wàn)個(gè)連續(xù)運(yùn)行測(cè)站，積累了大量的GPS數(shù)據(jù)，因此在實(shí)際應(yīng)用中，有必要采用新的模型對(duì)其進(jìn)行處理。

半?yún)?shù)模型是一種既包含參數(shù)又包含非參數(shù)分量的統(tǒng)計(jì)回歸模型，且具有較強(qiáng)的解釋能力。Nikolaidis提出的模型分離季節(jié)性信號(hào)會(huì)造成模型殘差，而半?yún)?shù)模型則可以彌補(bǔ)其系統(tǒng)偏差，因此，半?yún)?shù)模型是一種理想的分離GPS季節(jié)性信號(hào)的模型。解算半?yún)?shù)模型有正則化矩陣R選取和平滑因子α確定兩個(gè)關(guān)鍵問(wèn)題。由于GPS時(shí)間序列是連續(xù)的，且相鄰時(shí)刻的模型誤差相差不大，可采用文獻(xiàn)［7］的方法確定正則化矩陣R。確定平滑因子α的方法已有多種，例如廣義交叉證認(rèn)法［8］、L 曲線法［9］等，但這些方法都存在一些弊端，例如廣義交叉證認(rèn)法需事先給定平滑因子α的取值范圍及其步長(zhǎng)，對(duì)于單個(gè)時(shí)間序列并無(wú)太大影響，但若處理大量數(shù)據(jù)，其效率則至關(guān)重要；L曲線法也需事先選擇一些平滑因子α，且以α為自變量的加權(quán)范數(shù)應(yīng)該包含L 曲線的特征點(diǎn)，因此也不適用于大量數(shù)據(jù)處理?；诖?，本文基于噪聲時(shí)間序列Hurst指數(shù)的特性，提出一種二分搜索法確定平滑因子α。通過(guò)對(duì)模擬數(shù)據(jù)及實(shí)際的GPS坐標(biāo)序列進(jìn)行分析，驗(yàn)證該方法的有效性，尤其適用于大量GPS站時(shí)間序列季節(jié)性信號(hào)的分離。

1 基于Hurst指數(shù)確定半?yún)?shù)模型的平滑因子

1.1 半?yún)?shù)方法

半?yún)?shù)模型為：

其中，L為n維觀測(cè)向量，X為u維參數(shù)向量，A為系數(shù)矩陣，Δ為噪聲，S為信號(hào)，是非隨機(jī)未知量。

方程（1）的誤差方程形式為：

其對(duì)應(yīng)的估計(jì)準(zhǔn)則為：

其中，R為正則化矩陣，描述了非參數(shù)的光滑性；α為平滑因子，用于調(diào)節(jié)擬合部分VTPV與光滑部分的 平 衡。

在補(bǔ)償最小二乘原則下，根據(jù)拉格朗日極值法求解，可分別得到非參數(shù)分量S與參數(shù)分量X的估值：

可以看出，解算半?yún)?shù)模型的關(guān)鍵是選取合適的正則化矩陣R與平滑因子α。在實(shí)際應(yīng)用中，若L是一個(gè)觀測(cè)序列，且相鄰時(shí)刻的模型誤差相差不大時(shí)，可令正則化矩陣?。?/p>

其中，

此時(shí)，由于正則化矩陣R秩虧（rank（R）＝n－1＜n），為確保估計(jì)準(zhǔn)則（3）有唯一解，可增加一個(gè)約束條件［7］：

1.2 基于Hurst指數(shù)確定平滑因子

Hurst指數(shù)是一種判別時(shí)間序列是否對(duì)于時(shí)間有依賴的參數(shù)。一般來(lái)說(shuō)，當(dāng)時(shí)間序列的Hurst指數(shù)為0.5時(shí)，表示該時(shí)間序列是布朗運(yùn)動(dòng)，也即白噪聲；當(dāng)Hurst指數(shù)＜0.5 時(shí)，則表示該時(shí)間序列是高斯隨機(jī)變量過(guò)程；當(dāng)Hurst指數(shù)＞0.5時(shí)，表明該時(shí)間序列是長(zhǎng)期記憶相關(guān)的。在實(shí)際應(yīng)用中，一般將Hurst指數(shù)≤0.5的時(shí)間序列數(shù)據(jù)當(dāng)作噪聲［10］，若半?yún)?shù)模型將信號(hào)完全提取出來(lái)，那么其殘差序列的Hurst指數(shù)必定等于0.5。因此，可以利用Hurst指數(shù)來(lái)獲取最佳的平滑因子α。Hurst指數(shù)的計(jì)算方法有多種，本文選取受數(shù)據(jù)長(zhǎng)度及趨勢(shì)或周期影響較小的聚合方差法，其基本原理如下：

1）對(duì)于給定的時(shí)間序列Xi，i＝1，2，…，N，根據(jù)長(zhǎng)度m，將時(shí)間序列Xi分割成［N／m］＝k個(gè)子區(qū)間；

2）計(jì)算各子區(qū)間的平均值Xm（k）；

3）計(jì)算各子區(qū)間平均值的方差var（Xm）；

4）采用最小二乘擬合雙對(duì)數(shù)圖上（m，var（Xm））的斜率p；

5）Hurst指數(shù)為H＝（p／2）＋1。

二分搜索法的原理如圖1所示。

2 模擬實(shí)驗(yàn)分析

2.1 模擬數(shù)據(jù)

為驗(yàn)證本文提出的基于Hurst指數(shù)的二分搜索法確定半?yún)?shù)模型平滑因子的可靠性以及半?yún)?shù)模型分離GPS季節(jié)性信號(hào)的有效性，通過(guò)下面公式模擬了3 000個(gè)數(shù)據(jù)：

其中，y（ti）是時(shí)變季節(jié)性信號(hào)的模擬數(shù)據(jù)；a、b、d、e是常數(shù)；ti為GPS年積日；ε（ti）為隨機(jī)噪聲，服從正態(tài)分布N（0，0.25）；c（ti）是振幅變化因子，其形式如下：

圖2為模擬數(shù)據(jù)，其中圖2（a）、（b）分別為固定振幅的年周期、半年周期信號(hào)與加入振幅變化因子后的年周期、半年周期信號(hào)，圖2（c）為周期信號(hào)疊加得到的時(shí)變季節(jié)性信號(hào)，圖2（d）為加入噪聲后的時(shí)變季節(jié)性信號(hào)。

圖1 二分法確定平滑因子α的流程Fig.1 The flowchart of dichotomy to determine smoothing parameterα

圖2 季節(jié)性信號(hào)坐標(biāo)序列模擬數(shù)據(jù)Fig.2 Simulated coordination time series with seasonal signals

2.2 基于Hurst指數(shù)的二分搜索法的可行性檢驗(yàn)

為驗(yàn)證基于Hurst指數(shù)的二分搜索法確定半?yún)?shù)模型平滑因子的可行性，首先采用不同的平滑因子α代入半?yún)?shù)模型，計(jì)算殘差序列的Hurst指數(shù)。圖3是不同的平滑因子α對(duì)應(yīng)的半?yún)?shù)模型殘差序列的Hurst指數(shù)。

圖3 不同平滑因子對(duì)應(yīng)的半?yún)?shù)模型殘差序列的Hurst指數(shù)Fig.3 The Hurst index of semi－parametric model residual series for different smoothing parameter

從圖3 可以看出，平滑因子α在接近0 時(shí)，Hurst指數(shù)出現(xiàn)遞減，這是由于平滑因子α極小，造成信號(hào)被過(guò)度提?。▓D4），從而導(dǎo)致其殘差序列極度平穩(wěn)，使得Hurst指數(shù)在估計(jì)時(shí)呈現(xiàn)下降，但該變化并不影響本文提出的二分搜索法。當(dāng)平滑因子α超過(guò)某個(gè)值之后，其對(duì)應(yīng)的Hurst指數(shù)呈現(xiàn)逐漸遞增，這說(shuō)明采用二分搜索法是合理的，能夠確定半?yún)?shù)模型的平滑因子α。上述實(shí)驗(yàn)分析證明二分搜索法是可行的，但并未能反映出二分搜索法確定的α是較優(yōu)的。為此，本文又分別計(jì)算了不同平滑因子α對(duì)應(yīng)模擬數(shù)據(jù)分離的SRMS與NRMS，其定義見文獻(xiàn)［11］。

圖4 不同平滑因子對(duì)應(yīng)的均方根誤差Fig.4 Different smooth factors corresponding to the root mean square error

從圖4看出，當(dāng)平滑因子α極小時(shí)，其信號(hào)被過(guò)度提取，從而使得噪聲也被分離到系統(tǒng)誤差中，使得SRMS值增大，NRMS值噪聲變?。划?dāng)平滑因子α極大時(shí)，系統(tǒng)誤差并沒有被充分提取，從而使得SRMS值與NRMS值都增大。從各分圖的局部放大圖可以看出，當(dāng)平滑因子α在100左右時(shí)，其Hurst指數(shù)為0.35，對(duì)應(yīng)的SRMS值最小，約為0.076，其NRMS值為0.485，與模擬噪聲僅差0.015，說(shuō)明此時(shí)分離得到的信號(hào)與模擬信號(hào)最接近，一般認(rèn)為該平滑因子α的分離效果最好。但SRMS值是在信號(hào)已知情況下得到的，而在實(shí)際應(yīng)用中，信號(hào)一般是未知的，因此，無(wú)法通過(guò)SRMS值確定平滑因子α。二分搜索法確定的平滑因子為171.56，盡管與模擬最優(yōu)平滑因子α相差較大，但從圖4的局部放大圖中可以看出，當(dāng)平滑因子α為171.56時(shí)，其對(duì)應(yīng)的SRMS值為0.084，NRMS 值為0.493，與模擬最優(yōu)平滑因子α對(duì)應(yīng)的SRMS 值與NRMS 值相差微乎其微，與模擬信號(hào)及噪聲也十分接近。

圖5 模擬數(shù)據(jù)的擬合殘差Fig.5 The fitting residuals of simulated data

2.3 季節(jié)性信號(hào)提取比較分析

為比較半?yún)?shù)模型相對(duì)于傳統(tǒng)模型的優(yōu)越性，分別采用兩種模型提取季節(jié)性信號(hào)，其結(jié)果如圖5所示，其中，圖5（a）表示兩種模型提取的信號(hào)與模擬信號(hào)的殘差序列，圖5（b）表示提取的信號(hào)與加噪模擬信號(hào)的殘差序列。由圖5 可以看出，傳統(tǒng)模型最小二乘擬合提取季節(jié)性信號(hào)有較大的系統(tǒng)偏差，尤其是在振幅變化處；由基于Hurst指數(shù)二分搜索法確定光滑因子α的半?yún)?shù)模型則基本完整地提取出其季節(jié)性信號(hào)，避免了系統(tǒng)偏差殘留，且模擬季節(jié)性信號(hào)的殘差波動(dòng)平穩(wěn)，基本為0。這充分說(shuō)明，在振幅變化的季節(jié)性信號(hào)中，再使用傳統(tǒng)模型最小二乘擬合已經(jīng)無(wú)法滿足其實(shí)際需要，而半?yún)?shù)模型則能很好地解決系統(tǒng)偏差問(wèn)題，但其平衡因子α的選取是關(guān)鍵。而本文提出的基于Hurst指數(shù)的二分搜索法則彌補(bǔ)了該缺陷，能夠極大地提高其效率，從而可以推廣到GPS時(shí)間序列分析應(yīng)用中。

3 GPS實(shí)測(cè)坐標(biāo)序列分析

本實(shí)驗(yàn)采用IGS 站ZIMM 的GPS 坐標(biāo)序列，該站坐標(biāo)時(shí)間序列不僅跨度較長(zhǎng)，而且季節(jié)性信號(hào)呈現(xiàn)明顯的振幅變化。圖6 為ZIMM 站U方向的坐標(biāo)序列半?yún)?shù)模型與傳統(tǒng)參數(shù)模型的對(duì)比。從圖6（b）可以看出，半?yún)?shù)模型充分地分離出GPS坐標(biāo)序列中的趨勢(shì)信號(hào)以及季節(jié)性信號(hào)，克服了傳統(tǒng)參數(shù)模型描述實(shí)際模型不準(zhǔn)確的缺陷。圖6（b）的殘差序列整體平穩(wěn)，相對(duì)于圖6（c）沒有顯著的模型殘差，這也再次說(shuō)明本文提出的基于Hurst指數(shù)的二分搜索法確定半?yún)?shù)模型的平滑因子α是有效的。傳統(tǒng)參數(shù)模型由于沒有顧及到季節(jié)性信號(hào)的時(shí)變特性，其擬合的殘差序列出現(xiàn)了明顯的系統(tǒng)偏差。同時(shí)，GPS坐標(biāo)時(shí)間序列不僅包含年周期、半年周期信號(hào)，可能還存在其他頻率的周期信號(hào)，而傳統(tǒng)參數(shù)模型并未建模，從而可能會(huì)造成較大的系統(tǒng)偏差，半?yún)?shù)模型的非參數(shù)部分則恰好可以彌補(bǔ)該部分系統(tǒng)偏差。因此，采用半?yún)?shù)模型建模相對(duì)于傳統(tǒng)參數(shù)模型具有顯著的優(yōu)勢(shì)。

圖6 ZIMM 站U 方向信號(hào)擬合及殘差Fig.6 The signal fitting and residuals in U direction at ZIMM station

為進(jìn)一步說(shuō)明半?yún)?shù)模型的優(yōu)勢(shì)，又采用Lomb－Scargle功率譜分析法，對(duì)原始坐標(biāo)序列及兩種模型擬合殘差序列進(jìn)行分析，結(jié)果如圖7所示。從圖7中的豎線可以看出，半?yún)?shù)模型與傳統(tǒng)參數(shù)模型都減少了年周期信號(hào)與半年周期信號(hào)的功率，但半?yún)?shù)模型相對(duì)于傳統(tǒng)參數(shù)模型在低頻帶0.1～1.5cpy處有顯著的差異，這是因?yàn)榘雲(yún)?shù)模型將低頻信號(hào)基本分離出來(lái)，使其殘差序列減小，信噪比增大，從而使其功率譜密度減小。

圖7 ZIMM 站U 方向的功率譜分析Fig.7 Power spectral density in U direction at ZIMM station

4 結(jié) 語(yǔ)

1）半?yún)?shù)模型提取周期信號(hào)的效果要優(yōu)于傳統(tǒng)參數(shù)模型。在時(shí)變周期信號(hào)的分離中，半?yún)?shù)模型克服了傳統(tǒng)參數(shù)模型描述實(shí)際模型不準(zhǔn)確的缺陷，尤其是GPS坐標(biāo)序列中包含年周期信號(hào)與半年周期信號(hào)時(shí)。

2）驗(yàn)證了基于Hurst指數(shù)的二分搜索法確定半?yún)?shù)模型平滑因子的有效性與可靠性，明顯提高了半?yún)?shù)模型數(shù)據(jù)處理的效率，為大量、自動(dòng)化GPS時(shí)間序列分析提供了一種新的選擇。

盡管基于Hurst指數(shù)的二分搜索法確定半?yún)?shù)模型平滑因子具有上述優(yōu)點(diǎn)，但Hurst指數(shù)的計(jì)算精度受其估計(jì)方法以及時(shí)間序列跨度的影響，應(yīng)予重視。

［1］田云鋒.GPS坐標(biāo)時(shí)間序列中的異常高頻周期性噪聲［J］.測(cè)繪科學(xué)，2011，36（1）：26－28（Tian Yunfeng.Anomalous High Frequency Seasonal Noises in GPS Positions Time Series［J］.Science of Surveying and Mapping，2011，36（1）：26－28）

［2］Freymueller J T.Seasonal Position Variations and Regional Reference Frame Realization［A］／／Geodetic Reference Frames［M］.Springer Berlin Heidelberg，2009

［3］盧辰龍，匡翠林，戴吾蛟，等.采用變系數(shù)回歸模型提取GPS 坐標(biāo)序列季節(jié)性信號(hào)［J］.大地測(cè)量與地球動(dòng)力學(xué)，2014，34（5）：94－100（Lu Chenlong，Kuang Cuilin，Dai Wujiao，et al.Extracting Seasonal Signals from Continuous GPS Time Series Based on Varying－Coefficient Regression Models［J］.Journal of Geodesy and Geodynamics，2014，34（5）：94－100）

［4］Nikolaidis R.Observation of Geodetic and Seismic Deformation with the Global Positioning System［D］.San Diego：University of California，2002

［5］Blewitt G，Lavallée D.Effect of Annual Signals on Geodetic Velocity［J］.Journal of Geophysical Research，2002，107（B7）：2 145

［6］鄒蓉，丁開華，楊少敏，等.利用GRACE 改正地球參考框架的季節(jié)性變化［J］.大地測(cè)量與地球動(dòng)力學(xué)，2014，34（6）：50－54（Zou Rong，Ding Kaihua，Yang Shaomin，et al.Incorporating Seasonal Variations of Global and Regional ITRF Solutions Based on GRACE Measurements［J］.Journal of Geodesy and Geodynamics，2014，34（6）：50－54）

［7］孫海燕，吳云.半?yún)?shù)回歸與模型精化［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2002，27（2）：172－174（Sun Haiyan，Wu Yun.Semiparametric Regression and Mode Refining［J］.Geomatics and Information Science of Wuhan University，2002，27（2）：172－174）

［8］胡宏昌，孫海燕.正規(guī)矩陣R及平滑因子α的選取［J］.測(cè)繪工程，2004，12（4）：5－8（Hu Hongchang，Sun Haiyan.Choice of the Regular MatrixRand Smoothing parameterα［J］.Engineering of Surveying and Mapping，2004，12（4）：5－8）

［9］王振杰，歐吉坤，曲國(guó)慶，等.用L－曲線法確定半?yún)?shù)模型中的平滑因子［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2004，29（7）：651－653（Wang Zhenjie，Ou Jikun，Qu Guoqing，et al.Determining the Smoothing Parameter in Semi－parametric Model Using L－Curve Method［J］.Geomatics and Information Science of Wuhan University，2004，29（7）：651－653）

［10］盧辰龍，匡翠林，易重海，等.奇異譜分析濾波法在消除GPS多路徑中的應(yīng)用［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2015，40（7）：924－931（Lu Chenlong，Kuang Cuilin，Yi Zhonghai，et al.Singular Spectrum Analysis Filter Method for Mitigation of GPS Multipath Error［J］.Geomatics and Information Science of Wuhan University，2015，40（7）：924－931）

［11］羅勇，匡翠林，盧辰龍，等.基于SSA 的GPS坐標(biāo)序列去噪及季節(jié)信號(hào)提取［J］.大地測(cè)量與地球動(dòng)力學(xué)，2015，35（3）：391－395（Luo Yong，Kuang Cuilin，Lu Chenlong，et al.GPS Coordinate Series Denoising and Seasonal Signal Extraction Based on SSA［J］.Journal of Geodesy and Geodynamics，2015，35（3）：391－395）