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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

其中，

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

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

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

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

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

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

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

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

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

2 模擬實驗分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 GPS實測坐標(biāo)序列分析

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

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

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

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

4 結(jié) 語

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

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

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

［1］田云鋒.GPS坐標(biāo)時間序列中的異常高頻周期性噪聲［J］.測繪科學(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é)性信號［J］.大地測量與地球動力學(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］.大地測量與地球動力學(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é)報：信息科學(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］.測繪工程，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］王振杰，歐吉坤，曲國慶，等.用L－曲線法確定半?yún)?shù)模型中的平滑因子［J］.武漢大學(xué)學(xué)報：信息科學(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é)報：信息科學(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é)信號提取［J］.大地測量與地球動力學(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）