顧及噪聲的CORS 站高程非線性建模方法

袁壯志，楊貴軍，劉明星

（1.遼寧工程技術大學 測繪與地理科學學院，遼寧 阜新 123000；2.北京農業信息技術研究中心，北京 100097）

0 引言

連續運行參考站(continuously operating reference stations, CORS）的水平與高程坐標分量運動特性預測分析，是建立和維持2000 國家大地坐標系（China geodetic coordinate system 2000,CGCS2000）準確系統框架和現勢性的基礎[1]。現有的研究結果表明：CORS 站高程序列分量的變化同時具有線性和非線性特征[2]；不僅包含白噪聲，也包含有色噪聲[3]。因此，建立1 個完整、準確的CORS 站高程分量變化規律模型，需要包含噪聲項。

現有文獻研究主要集中在對CORS 站高程分量數據建模預測、噪聲模型分析等方面。文獻[4]引入了奇異譜分析的技術進行高程數據非線性建模，并指出該技術相比地球物理效應分析建模方法具有很明顯的優勢；文獻[5]通過功率譜分析的方法，探測得到CORS 站垂直方向主要存在半年周期運動和年周期運動，在此基礎上建立了CORS 站高程坐標分量序列的非線性速度場模型；文獻[6]首先利用給定固定周期項的方法，建立線性最小二乘模型，然后將其解作為非線性最小二乘模型的迭代初值，實現了CORS 站高程時間序列數據的非線性周期擬合建模；文獻[7]采用功率譜分析高程殘差序列的噪聲性質，并利用極大似然法定量地估計殘差序列中有色噪聲的分量；文獻[8]首先利用非線性擬合的方法對高程數據建模，然后對未來時間的高程數據進行預測，該預測方法存在的問題是，對高程數據建模時只考慮了周期項和趨勢項，忽視了噪聲項的影響，使得預測的模型僅適用于短時間內的預測，而在較長時間的預測數據，往往與實際數據會存在較大的誤差；文獻[9]利用不同組合的噪聲模型，分析了CORS 站高程坐標時間序列的噪聲，并且考慮了積雪深度、土壤濕度、大氣壓負載等外界因素對CORS 站位移的影響，其結果證明了中國區域內CORS 站坐標序列的噪聲特性主要有2 類，即白噪聲與閃爍噪聲，白噪聲與帶通冪律噪聲。

截至目前，國內外眾多CORS 站已經記錄了長達20 余年的連續觀測數據，該數據為全球研究地殼板塊運動與地質災害等的學者提供了重要的信息[10]，本論文從斯克里普斯軌道和常駐陣列中心(Scripps Orbit and Permanent Array Center,SOPAC）獲取CORS 站高程數據，分別對CORS 站高程分量數據建立非線性最小二乘周期擬合模型（nonlinear least square periodic fitting model, NL）、噪聲項自回歸（auto-regressive, AR）模型。然后將NL 模型和AR 模型進行組合，對CORS 站的高程坐標序列進行預測[11]。

1 CORS 站高程序列建模預測方法

目前CORS 站高程分量建模方法，包括固定周期項的線性最小二乘建模和將周期值作為未知項的非線性最小二乘建模，這2 種方法主要對高程數據的趨勢項和周期項進行建模，未考慮噪聲項的影響。本文首先利用非線性最小二乘周期擬合方法對高程數據建模，其次對得到的殘差序列（噪聲項）建立AR 模型，最后將NL 模型和AR 模型進行組合，對CORS 站的高程坐標序列進行預測，后文中將組合模型記為NL/AR 模型。

1.1 非線性最小二乘周期擬合建模方法

非線性最小二乘周期擬合模型是將各個周期項作為未知值，加入了每個未知周期項對應的初相和頻率未知數，得到非線性最小二乘擬合模型為

式中：a 為常數項；b 為線性速度值；t 為CORS 站高 程 值 對 應 的 單 日 歷 元 時 間； A1、 A2、 A3分 別 對 應了各自周期項的振幅值； f1、 f2、 f3分別對應了各自周期項的頻率值； φ1、 φ2、 φ3分別對應了各自周期項初相值。

將線性最小二乘模型的解作為非線性最小二乘周期擬合模型式（1）的初始值，記為H0，并將非線性最小二乘模型式（1）在初始值H0處，按照泰勒級數展開至1 階項，即得到線性化之后的誤差方程式為

然后根據最小二乘原理TV PV =min，求解出未知參數。

1.2 噪聲項AR 建模方法

對利用非線性最小二乘周期擬合模型計算得到得殘差序列（噪聲項），建立AR 模型，具體表達式為

式中： Wt為噪聲序列的實際值； εt為誤差；α 、β 為對應的系數項；p、q 為利用自相關函數和偏自相關函數確定，殘差序列的自相關函數定義為

式中：ACF 表示殘差平方的自相關系數值； r ( s, t)為序列自協方差； D( X (t )) 和D( X ( s)) 表示不同時刻方差。其中序列自協方差、不同時刻方差定義為：

偏自相關函數的數學表達式定義為

式中：kkφ 表示滯后數為k 的偏自相關函數值；D 為系數行列式；kD 為將D 的第k 列換為常數項，ρi為樣本的自相關系數。

本文采用均方根誤差（root mean squared error，RMSE）作為指標，來評價NL/AR 模型的高程預測效果。

RMSE 反映了預測值與真實值之間的偏差程度[12]，其計算方法為

式中：n 為預測數據個數； dΔ 為預測值與真實值之間的偏差程度。

2 數據預處理

本文從SOPAC 網站上下載了近20 余年的高程分量數據。發現每個CORS 站高程數據均有缺值現象，利用3 次樣條法[13]對缺省的數據進行插值處理。插值情況如下：BJFS 站插值110 個，WUHN 站插值210 個，BARH 站插值150 個，IENG 站插值189 個。對插值后的數據進行粗差剔除，利用中誤差公式剔除了BJFS 站的12 個數據、WUHN 站的25 個數據、BARH 站的18 個數據、IENG 站的14 個數據。剔除后的數據作為后續建模預測的研究對象，圖1 給出了剔除數據的結果，其中黑色小圓圈代表超出3 倍中誤差的高程數據（粗差）。

圖1 粗差剔除結果

3 實驗結果與分析

本文以進行插值和粗差剔除后的CORS 站高程數據為研究對象，分別做了BJFS 站、WUHN 站、BARH 站、IENG 站等4 個CORS 站的高程數據非線性最小二乘周期擬合建模、噪聲項AR 建模、NL/AR 模型的高程數據預測。最后對比了NL/AR模型與未考慮噪聲項的非線性最小二乘周期擬合模型預測高程數據的效果。

3.1 非線性最小二乘周期擬合建模分析

BJFS 站建模區間（單日歷元時間）為1999.804～2018.741（單位為 d），WUHN 站為 1999.067～2013.999，BRAH 站為1998.752～2018.673，IENG 站為2004.001～2018.640。表1 為4 個CORS 站高程數據非性最小二乘周期擬合建模得到的均方根誤差和未知參數值。

表1 非線性建模結果值

圖2 為建模結果，其中縱軸代表CORS 站高程坐標，單位為mm；橫軸代表與各個CORS 站高程坐標所相應的時間；灰色散點代表CORS 站原始高程值，黑色粗線代表模型擬合得到的高程值對應曲線，4 個CORS 站的高程數據建模精度在6 mm 左右。

圖2 非線性建模結果值

3.2 噪聲項AR 建模預測

3.1 節建立的非線性最小二乘周期擬合模型實際上是對CORS 站高程數據中的周期項、趨勢項進行的建模分析，未對噪聲項建模，故模型的殘差項（真實高程值與擬合高程值之差）即為噪聲項。利用AR 對4 個CORS 站高程數據的噪聲序列進行建模預測：表2 給出了建模預測精度值；圖3 給出了預測時間分別為0.5、1 a 的殘差序列預測結果，縱軸為殘差真值和預測值，單位是 m，橫軸為預測時間，其中殘差真值為非線性最小二乘周期擬合建模得到的殘差值，殘差預測值表示利用AR模型對殘差序列建模得到的預測值。

表2 AR 建模精度值

圖3 殘差序列AR 建模結果值

3.3 高程數據預測

本節對CORS 站高程數據建立了NL/AR 預測模型，并與非線性最小二乘周期擬合模型預測進行對比，預測時長為0.5 和1 a。表3 為0.5 a 和1 a 預測精度結果精度值，評價指標為RMSE，其中NL表示非線性最小二乘周期擬合預測，NL/AR 表示組合模型預測。圖4 給出了非線性最小二乘周期擬合和NL/AR 預測模型的結果，縱軸表示高程真實值和2 種不同方法的預測高程值，單位是mm，橫軸表示預測時間。由表3 的預測結果精度值統計可得：對于2 種預測方法，預測區間越長，相應的預測結果精度越低，反之，預測時間短，則預測精度高；采用NL/AR 模型的預測方法相對于NL模型的預測提高了精度（RMSE 越小，精度越高），說明文章提出的顧及噪聲的CORS 站高程非線性建模方法是可行的，由于加入了噪聲項的AR 模型，使得高程預測模型更加完整（包含了周期項、趨勢項、噪聲項）。

圖4 4 個CORS 站高程數據預測結果

表3 預測結果精度值（RMSE）

4 結束語

本文以國內外4 個CORS 站近20 余年的高程時間數據為研究對象，提出了1 種顧及噪聲的CORS站高程非線性建模方法，分別建立了非線性最小二乘周期擬合模型和AR 模型，最后將2 種模型組合進行CORS 站高程數據的預測，得出了如下結論：

1）未考慮噪聲項的非線性最小二乘周期擬合模型，對CORS 站高程數據建模是不完善的，因為高程數據中一定會含有噪聲，建立包含噪聲項、周期項、趨勢項數學模型才能進行高程數據的準確預測。

2）從4 個CORS 站預測結果的精度值統計可知，提出的NL/AR 模型對CORS 站高程數據的預測較NL 模型預測精度高。