袁壯志，楊貴軍，劉明星

（1.遼寧工程技術(shù)大學(xué) 測繪與地理科學(xué)學(xué)院，遼寧 阜新 123000；2.北京農(nóng)業(yè)信息技術(shù)研究中心，北京 100097）

0 引言

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

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

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

1 CORS 站高程序列建模預(yù)測方法

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

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

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

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

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

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

1.2 噪聲項AR 建模方法

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

式中： Wt為噪聲序列的實(shí)際值； εt為誤差；α 、β 為對應(yīng)的系數(shù)項；p、q 為利用自相關(guān)函數(shù)和偏自相關(guān)函數(shù)確定，殘差序列的自相關(guān)函數(shù)定義為

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

偏自相關(guān)函數(shù)的數(shù)學(xué)表達(dá)式定義為

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

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

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

式中：n 為預(yù)測數(shù)據(jù)個數(shù)； dΔ 為預(yù)測值與真實(shí)值之間的偏差程度。

2 數(shù)據(jù)預(yù)處理

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

圖1 粗差剔除結(jié)果

3 實(shí)驗(yàn)結(jié)果與分析

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

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

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

表1 非線性建模結(jié)果值

圖2 為建模結(jié)果，其中縱軸代表CORS 站高程坐標(biāo)，單位為mm；橫軸代表與各個CORS 站高程坐標(biāo)所相應(yīng)的時間；灰色散點(diǎn)代表CORS 站原始高程值，黑色粗線代表模型擬合得到的高程值對應(yīng)曲線，4 個CORS 站的高程數(shù)據(jù)建模精度在6 mm 左右。

圖2 非線性建模結(jié)果值

3.2 噪聲項AR 建模預(yù)測

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

表2 AR 建模精度值

圖3 殘差序列AR 建模結(jié)果值

3.3 高程數(shù)據(jù)預(yù)測

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

圖4 4 個CORS 站高程數(shù)據(jù)預(yù)測結(jié)果

表3 預(yù)測結(jié)果精度值（RMSE）

4 結(jié)束語

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

1）未考慮噪聲項的非線性最小二乘周期擬合模型，對CORS 站高程數(shù)據(jù)建模是不完善的，因?yàn)楦叱虜?shù)據(jù)中一定會含有噪聲，建立包含噪聲項、周期項、趨勢項數(shù)學(xué)模型才能進(jìn)行高程數(shù)據(jù)的準(zhǔn)確預(yù)測。

2）從4 個CORS 站預(yù)測結(jié)果的精度值統(tǒng)計可知，提出的NL/AR 模型對CORS 站高程數(shù)據(jù)的預(yù)測較NL 模型預(yù)測精度高。