張恒璟，崔東東，程鵬飛

1. 遼寧工程技術(shù)大學(xué)測繪與地理科學(xué)學(xué)院，遼寧 阜新 123000； 2. 中國測繪科學(xué)研究院，北京 100830

CORS站坐標(biāo)運動特征預(yù)測分析是維持2000中國大地坐標(biāo)系統(tǒng)(China geodetic coordinate system 2000，CGCS2000)框架準(zhǔn)確性和現(xiàn)勢性的重要基礎(chǔ)[1]。通過對國內(nèi)多個CORS站高程坐標(biāo)分量序列特征分析發(fā)現(xiàn)，CORS站在垂直方向表現(xiàn)出明顯的非線性和周期性特征，并且不同的CORS站對應(yīng)周期性變化模式和運動規(guī)律也不盡相同[2]。影響全球定位系統(tǒng)(global positioning system，GPS)高程分量周期性變化的因素很多，主要包括地殼運動、冰雪、土壤水分、大氣壓負(fù)載以及CORS站測量技術(shù)的高低等[3-4]。在分析CORS站高程坐標(biāo)分量運動特征時，只考慮線性運動，直接在線性最小二乘擬合模型中對周期項賦予固定數(shù)值[5]，而忽略其非線性運動特征，導(dǎo)致建立的CORS站高程線性速度場模型對于垂向速率估計和誤差估計存在很大的偏差[6-8]。同時CORS站周圍存在物理、地理因素及測量儀器等產(chǎn)生的噪聲，也會影響CORS站高程運動的規(guī)律。在對我國CORS站高程時間序列數(shù)據(jù)噪聲模型研究時，文獻[9]利用功率譜分析CORS高程序列的噪聲性質(zhì)，并采用極大似然法定量地估計噪聲序列中有色噪聲的分量。文獻[10]利用不同組合的噪聲模型對CORS站坐標(biāo)時間序列的噪聲進行分析，并計算了積雪深度、大氣壓負(fù)載、土壤濕度負(fù)載等對CORS站位移的影響，得到中國區(qū)域CORS站坐標(biāo)序列的噪聲特性主要表現(xiàn)為白噪聲+閃爍噪聲和白噪聲+帶通冪律噪聲。文獻[11]利用AR模型對CORS站高程噪聲序列建模，并與CORS站高程非線性擬合模型結(jié)合進行高程運動的預(yù)報。

研究結(jié)果表明，CORS站高程分量建模的殘差序列表現(xiàn)為有色噪聲，從統(tǒng)計意義上是否為非平穩(wěn)的隨機過程，需要進一步研究。統(tǒng)計學(xué)采用異方差法分析時間序列的平穩(wěn)性特性，如果存在異方差特性，則時間序列是非平穩(wěn)的隨機過程。統(tǒng)計學(xué)中對時間序列是否存在異方差特性進行了研究，大多數(shù)研究者對截面數(shù)據(jù)產(chǎn)生的異方差給予了足夠多的關(guān)注而放松了對時間序列數(shù)據(jù)也會產(chǎn)生異方差的警惕[12-13]。文獻[14]于1982年提出了在時間序列背景下也有可能出現(xiàn)異方差性的觀點，并從理論上提出了一種觀測時間序列方差變動的方法，自回歸條件異方差(autoregressive conditional heteroskedasticity，ARCH)檢驗。ARCH檢驗認(rèn)為時間序列存在的異方差性為 ARCH過程，通過檢驗該過程是否成立去判斷時間序列是否存在異方差[15]。除了ARCH方法檢驗異方差外，還有殘差平方圖分析法、帕克檢驗法、等級相關(guān)系數(shù)檢驗法、White檢驗法等[16]。

國內(nèi)外CORS站已經(jīng)積累了20余年的連續(xù)觀測數(shù)據(jù)，為研究地殼運動、地震活動等提供了重要的信息。本文采用SOPAC(scripps orbit and permanent array center)全球數(shù)據(jù)中心6個國內(nèi)外CORS站20余年的高程時間序列數(shù)據(jù)，設(shè)計非線性周期模型迭代算法，建立CORS站高程坐標(biāo)分量非線性最小二乘擬合模型，并以均方根誤差(root mean square error，RMSE)和擬合優(yōu)度(goodness of fit，目前沒有固定的簡寫形式，文中記為RNEW)為精度指標(biāo)衡量模型的優(yōu)劣。建立CORS站高程非線性擬合模型后，以殘差平方序列為研究對象[17]，利用ARCH檢驗法、殘差平方序列圖分析法判斷殘差平方序列的異方差特性，對具有異方差特性的殘差平方序列建立廣義自回歸條件異方差(generalized autoregressive conditional heteroskedasticity，GARCH)波動模型，描述其隨著時間變化的非平穩(wěn)波動情況[18]，驗證GARCH模型對非平穩(wěn)殘差平方序列建模的可行性，為CORS站高程殘差序列有色噪聲建模提供一種途徑。

1 CORS站高程時間序列非線性建模方法

y=a+b×t+A1sin(2π×t)+A2sin(4π×t)+A3sin(π×t)

(1)

非線性最小二乘模型將各個周期項作為未知值，加入了每個未知周期項對應(yīng)的初相和頻率未知數(shù)，例如式(2)中，M=[abA1f1φ1A2f2

φ2A3f3φ3]T是模型的未知參數(shù)。

y=a+b×t+A1sin(2×π×f1×t+φ1)+A2sin(2×

π×f2×t+φ2)+A3sin(2×π×f3×t+φ3)

(2)

線性模型(1)的解作為非線性模型(2)的初值，記為M0。將模型(2)在初值M0處按泰勒級數(shù)展開至一階項，得到線性化后的誤差方程式為V=Bm-l。

l=[l1l2…li]

LS求解非線性模型時，給定未知參數(shù)初值的精度不同，一次平差結(jié)果可能不滿足參數(shù)的精度要求，采用高斯-牛頓迭代算法，計算流程如圖1所示，其中|a|=|an-an-1|。

圖1 CORS站高程非線性建模流程Fig.1 The nonlinear modeling process of CORS stations height

為了評價CORS站高程時間序列非線性建模的效果，采用均方根誤差RMSE、可決系數(shù)(coefficient of determination，記為R2)、擬合優(yōu)度RNEW 3個精度評價指標(biāo)。

均方根誤差RMSE反映了測量值與真實值之間的偏離程度[19]

(3)

式中，n為測量次數(shù)；Δd為測量值與真實值之間的偏差。

可決系數(shù)R2是線性回歸方程對觀測值的擬合程度[20]

(4)

擬合優(yōu)度在統(tǒng)計學(xué)中用于判定非線性回歸方程對觀測值的擬合程度[21]，在對回歸方程擬合程度的解釋上，RNEW與R2意義相同，當(dāng)RNEW越接近于1，非線性模型擬合程度越高

(5)

式中，Q=∑r2；Y=∑y2；r是擬合殘差；y是序列觀測值。

2 CORS站高程非線性建模試驗

從SOPAC全球數(shù)據(jù)中心獲得國內(nèi)外6個CORS站20余年干凈的帶有趨勢項的高程坐標(biāo)時間序列數(shù)據(jù)，對于缺值數(shù)據(jù)，使用三次樣條法進行插值處理[22]。采用線性模型(1)和非線性模型(2)分別建立CORS站高程分量擬合模型。

表1是線性模型(1)的未知參數(shù)估計值。作為CORS站高程非線性高斯-牛頓迭代的初值，其中YAR2和AUCK兩站只給出了一年和半年周期項的初值，這兩個站的兩年周期項振幅估計值均小于0.1 mm，兩年周期項在高程分量運動模型中所占的比重可以忽略不計。

表1 非線性高斯-牛頓迭代法的初值

表2是非線性模型(2)的未知參數(shù)估計值及精度評價指標(biāo)。為了便于對比線性和非線性模型擬合結(jié)果的差異，將線性模型的建模精度指標(biāo)結(jié)果也放在了表2，其中(nl)rmse和(l)rmse分別表示非線性建模與線性建模的均方根誤差指標(biāo)。周期方面：6個CORS站高程非線性建模的周期值(頻率f的倒數(shù))均不是整數(shù)，其周期性變化并不是經(jīng)典的CORS站高程時間序列模型所采用的固定整年、整半年或兩年；近似年周期的估計誤差，除了AUCK站達到12%外，其他5個站均在1%；近似半年周期的估計誤差不超過18%；近似兩年周期的估計誤差不超過6%。振幅方面：近似年周期項對應(yīng)的振幅值介于5～8 mm，近似半年周期項的振幅介于1～2 mm，近似兩年周期項對應(yīng)的振幅值介于0.1～1.3 mm，說明6個CORS站的高程方向以近似年和半年周期性變化為主，近似兩年周期性變化基本忽略不計。精度評價指標(biāo)方面：6個CORS站的非線性模型RMSE精度在5 mm水平，低于線性模型1～2 mm，線性模型R2明顯低于非線性模型RNEW指標(biāo)0.1左右，非線性建模精度整體優(yōu)于線性模型。

圖2是6個CORS站高程分量線性和非線性建模結(jié)果，縱軸表示CORS站高程坐標(biāo)值，橫軸表示高程坐標(biāo)對應(yīng)的時間值，縱軸單位是mm，橫軸單位是年；綠色散點代表CORS站原始高程值，黑色粗實線代表線性模型擬合曲線，紅色粗實線代表非線性模型擬合曲線。從圖2可以看出，各個CORS站建立的非線性模型擬合效果優(yōu)于線性模型。

3 ARCH效應(yīng)檢驗與GARCH(p, q)波動模型

為了克服AR模型只能對平穩(wěn)序列數(shù)據(jù)建模的不足，引入經(jīng)濟學(xué)的序列波動分析方法對非平穩(wěn)序列數(shù)據(jù)建模，驗證序列波動GARCH分析方法對CORS站高程非平穩(wěn)殘差平方序列建模的可行性，為CORS站高程殘差序列建模和非線性速度場重構(gòu)奠定基礎(chǔ)。

3.1 ARCH效應(yīng)檢驗方法

異方差性是相對于同方差而言，若一個殘差平方序列隨時間的變化其值波動幅度較大，則稱該殘差平方序列具有異方差性。對CORS站高程非線性建模后的殘差平方序列檢驗其異方差性是否存在，即殘差平方序列是否存在ARCH效應(yīng)[23]。圖3是CORS站高程殘差平方序列ARCH效應(yīng)檢驗流程，分為兩步：首先檢驗殘差平方序列自相關(guān)特性，利用Q檢驗法及殘差平方序列自相關(guān)系數(shù)(autocorrelation coefficient，ACF)檢驗如圖4所示；接著檢驗殘差平方序列的異方差性，構(gòu)造t統(tǒng)計量檢驗流程如圖5所示。若CORS站高程非線性建模后的殘差平方序列同時存在自相關(guān)性和異方差性，則其具有ARCH效應(yīng)，即CORS站高程非線性建模后的殘差平方序列具有非平穩(wěn)性[24]。

表2 非線性高斯-牛頓迭代法求解出的模型參數(shù)值及精度評價值

圖2 6個CORS站高程數(shù)據(jù)建模結(jié)果Fig.2 The height data modeling results of 6 CORS stations

圖3 CORS站高程殘差平方序列非平穩(wěn)性檢驗流程Fig.3 The nonstationarity test procedure of CORS stations height residual square sequence

圖4 自相關(guān)性Q檢驗流程Fig.4 The autocorrelation of Q test procedure

圖4和圖5中，滯后數(shù)Lags值由殘差平方序列自相關(guān)系數(shù)值對應(yīng)的最大滯后數(shù)確定，試驗選擇的BJFS和BOGO兩個CORS站時間序列長度為21年，其殘差平方序列的自相關(guān)系數(shù)值對應(yīng)的最大滯后數(shù)均為21，即Lags=21。自相關(guān)性檢驗的原假設(shè)為殘差平方序列不存在自相關(guān)性，異方差性檢驗的原假設(shè)為殘差平方序列不存在異方差性，顯著性水平alpha=0.05，表示殘差平方序列相關(guān)性和異方差性的檢驗置信度為95%，h為檢驗殘差平方序列自相關(guān)性和異方差性程序中設(shè)置的邏輯值，若h=1表示拒絕原假設(shè)，h=0表示接受原假設(shè)。

自相關(guān)性檢驗的Q統(tǒng)計量可表示為

(6)

式中，ρj是殘差平方序列的j階自相關(guān)系數(shù)；T是殘差平方序列的總長度；p是設(shè)定的滯后階數(shù)。

殘差平方序列的自相關(guān)系數(shù)定義為

ACF=r(s,t)/[(DX(t)·DX(s))0.5]

(7)

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

r(s,t)=E[(X(s)-E(X(s)))·(X(t)-

E(X(t)))]

(8)

(9)

殘差平方序列的偏自相關(guān)系數(shù)(partial autocorrelation coefficient，PACF)定義如式(10)，φkk表示滯后數(shù)為k的偏自相關(guān)系數(shù)，由克萊姆法則求得，D為系數(shù)行列式，ρi為樣本的自相關(guān)系數(shù)，Dk為將D的第k列換為常數(shù)項。

圖5 異方差性ARCH檢驗流程Fig.5 The ARCH test procedure of heteroscedasticity

(10)

圖5異方差性t檢驗時，回歸統(tǒng)計矩陣X由殘差平方序列值對應(yīng)的時間構(gòu)成，反應(yīng)變量y為殘差平方序列值，預(yù)測響應(yīng)yhat為殘差平方預(yù)測值。t統(tǒng)計量stat定義為

(11)

3.1.1 殘差平方序列自相關(guān)性檢驗

自相關(guān)性指序列之間不是完全相互獨立的，而是存在某種相互關(guān)系[25]。BJFS和BOGO兩站的殘差平方序列如圖6所示，縱軸表示高程時間序列數(shù)據(jù)非線性建模后的殘差平方值，橫軸表示時間。為了檢驗CORS站高程殘差平方序列的自相關(guān)性，繪制兩個CORS站的殘差平方序列自相關(guān)如圖7所示，縱坐標(biāo)ACF表示CORS站殘差平方序列的自相關(guān)系數(shù)，橫坐標(biāo)Lags表示自相關(guān)系數(shù)對應(yīng)的滯后階數(shù)值。自相關(guān)系數(shù)ACF隨滯后階數(shù)Lags的變化并不是快速趨近于0，而是緩慢地變化為0，即殘差平方序列在不同時刻存在某種相互關(guān)系，并不是相互獨立的，通過自相關(guān)圖定性地判斷出兩站殘差平方序列具有自相關(guān)性。

圖6 兩個CORS站高程的殘差平方序列Fig.6 The residual squared sequence of two CORS stations height

圖7 兩個CORS站高程殘差平方序列的自相關(guān)圖Fig.7 The autocorrelation graph of two CORS stations height residual square sequences

為了定量地判斷兩站殘差平方序列的自相關(guān)性，采用Q檢驗法對兩站的殘差平方序列自相關(guān)性進行檢驗，利用式(7)計算殘差平方序列的自相關(guān)系數(shù)ACF，以最大滯后數(shù)(Lags=21)計算式(6)的Q統(tǒng)計量，得到BJFS、BOGO兩站的Q統(tǒng)計量分別為4.12和4.23，均大于顯著性水平為0.05時的Q臨界值3.48，故拒絕原假設(shè)，即邏輯量h=1，兩站的殘差平方序列具有自相關(guān)性。

3.1.2 殘差平方序列異方差性檢驗

由殘差平方序列圖6定性的判斷殘差平方序列的異方差性[26]：兩站的殘差平方序列隨時間變化，并不是一個常數(shù)，并且每隔一年左右都會出現(xiàn)集聚現(xiàn)象(尖峰厚尾)，即殘差平方序列隨時間波動幅度較大，定性的判斷出了兩站殘差平方序列具有異方差特性。

采用ARCH檢驗方法定量地判定殘差平方序列的異方差特性。ARCH過程如式(12)，i為ARCH過程的階數(shù)；α0為常數(shù)項；αi為ARCH過程階數(shù)對應(yīng)殘差平方值的系數(shù)；υt為隨機誤差；α0>0,αi≥0(i=1,2,…,p)

選擇2017年4月—2018年4月在本院就診的膝關(guān)節(jié)置換術(shù)患者330例作為研究對象。根據(jù)護理方式的不同，將其分為兩組，對照組和觀察組。對照組患者中，男82例，女83例，年齡為57～78歲，平均年齡為（72.2±4.6）歲，病程為3～12年，平均病程為（6.3±1.3）；觀察組組患者中，男85例，女80例，年齡為59～83歲，平均年齡為（73.2±5.1）歲，病程為2～11年，平均病程為（6.4±1.4）年。兩組患者在一般臨床資料比較，差異不具有統(tǒng)計學(xué)意義（P＞0.05），具有可比性。

(12)

由圖5計算兩站殘差平方序列的t統(tǒng)計量stat：BJFS站stat=2.98、BOGO站stat=2.81，均大于顯著性水平為0.05時的t臨界值2.10，故拒絕原假設(shè)，即邏輯值h=1，殘差平方序列存在異方差性，證明了兩站高程殘差平方序列是非平穩(wěn)隨機過程。

3.2 殘差平方序列GARCH(p, q)波動模型

為了定量描述殘差平方序列的波動，引入GARCH(p,q)模型對CORS站高程非平穩(wěn)殘差平方序列建模，GARCH模型可以有效地擬合具有長期記憶性的異方差函數(shù)，而ARCH模型是GARCH模型的一個特例[27]，即p=0的GARCH(p,q)模型。GARCH (p,q)模型定義如式(13)所示，i=1,2,…,p，j=1,2,…,q，其中p是模型拖尾值，它的閾值由序列自相關(guān)系數(shù)值基本趨于穩(wěn)定時對應(yīng)的滯后數(shù)Lags確定；q是模型截尾值，它的閾值由序列偏自相關(guān)系數(shù)值基本趨于穩(wěn)定時對應(yīng)的滯后數(shù)Lags確定

(13)

根據(jù)殘差平方序列自相關(guān)圖7判斷拖尾值p：BJFS站ACF在滯后階數(shù)Lags=6后趨于穩(wěn)定，模型拖尾值p的閾值為6，即p≤6；BOGO站ACF在滯后階數(shù)Lags=10后趨于穩(wěn)定，模型拖尾值p的閾值為10，即p≤10。

圖8是殘差平方序列偏自相關(guān)圖：縱坐標(biāo)PACF表示CORS站殘差平方序列的偏自相關(guān)系數(shù)值，橫坐標(biāo)Lags表示偏自相關(guān)系數(shù)值對應(yīng)的滯后階數(shù)值。BJFS站PACF在滯后階數(shù)Lags=5后趨于穩(wěn)定，模型截尾值q的閾值為5，即q≤5。BOGO站PACF在滯后階數(shù)Lags=4后趨于穩(wěn)定，模型截尾值q的閾值為4，即q≤4。

圖8 兩個CORS站高程殘差平方序列的偏相關(guān)圖Fig.8 The partial correlogram graph of two CORS stations height residual square sequences

由殘差平方序列自相關(guān)圖和偏自相關(guān)圖綜合判斷出了兩站GARCH模型的拖尾值p和截尾值q的閾值：BJFS站(p≤6，q≤5)，BOGO站(p≤10，q≤4)。

為了求出GARCH模型最佳的(p，q)值，本文參考文獻[28]中引用的赤池信息量準(zhǔn)則(Akaike information criterion，AIC)選取最佳模型拖尾和截尾值

AIC=2K-2ln(L)

(14)

式中，K是參數(shù)的數(shù)量(GARCH(p，q)模型的組合個數(shù))；L是似然函數(shù)，即殘差平方序列之和除以GARCH(p，q)模型的組合個數(shù)。

根據(jù)AIC準(zhǔn)則求出了最佳(p,q)值，BJFS站(p=1,q=1)，BOGO站(p=1,q=1)，并對兩站的殘差平方序列建立GARCH(1，1)波動模型

(15)

由極大似然估計準(zhǔn)則求解BJFS和BOGO站波動模型GARCH(1，1)的未知參數(shù)[29-30]分別見表3和表4，constant代表未知參數(shù)M中的常數(shù)項α0，GARCH(1)代表β，ARCH(1)代表α，offset代表u，由于4個值為常數(shù)或某項的系數(shù)，故無單位。

將表3和表4殘差平方序列GARCH(1,1)模型參數(shù)值分別代入式(15)，解算兩站的條件方差序列值，如圖9所示，縱軸表示殘差平方值與條件方差值，橫軸表示時間，黑色實線代表殘差平方序列，綠色點線代表條件方差序列。條件方差是在殘差平方基礎(chǔ)上描述的殘差平方序列波動情況，其數(shù)值比殘差平方值小，為了直觀地反映條件方差與殘差平方隨時間的變化規(guī)律，圖9將條件方差擴大到100倍，與殘差平方值在同一個數(shù)量等級。從圖9可知，條件方差序列隨時間的變化趨勢會在某一段時間急劇升高或降低，與殘差平方序列隨時間的變化趨勢一致，利用GARCH(1,1)模型可以對CORS站高程非平穩(wěn)的殘差平方序列波動情況進行描述。

圖9 殘差平方序列與條件方差序列趨勢變化比較Fig.9 Comparison of the trend between the residual squared sequence and the conditional variance sequence

表3 BJFS站殘差平方序列GARCH(1, 1)模型參數(shù)

Tab.3 GARCH (1, 1) model parameters of BJFS station residual squared sequence

parametervaluestandarderrort statisticconstant0.000000200.000006150.13291100GARCH(1)0.688887000.375217001.83597000ARCH(1)0.167274000.906637000.18450000offset0.000048960.000138280.35401500

表4 BOGO站殘差平方序列GARCH(1, 1)模型參數(shù)

Tab.4 GARCH (1, 1) model parameters of BOGO station residual squared sequence

parametervaluestandarderrort statisticconstant0.000000200.000001700.11776700GARCH(1)0.521826000.648599000.80454200ARCH(1)0.214017000.0120118017.8172000offset0.000026290.000118600.22171300

4 結(jié) 論

(1) 建立了CORS站高程分量非線性速度場周期擬合模型。CORS站高程分量周期性運動以近似年周期項為主，近似半年周期項次之，近似兩年變化的比重幾乎可以忽略不計。周期項并不是嚴(yán)格的年、半年或兩年，固定周期項的線性建模估計結(jié)果存在偏差，采用非線性模型對CORS站高程數(shù)據(jù)建模精度明顯優(yōu)于傳統(tǒng)固定周期項的線性模型。

(2) 建立了描述CORS站高程殘差平方序列非平穩(wěn)波動的GARCH(1,1)模型。非線性建模后的CORS站高程殘差平方序列自相關(guān)性和異方差性檢驗結(jié)果表明，殘差平方序列是非平穩(wěn)的隨機過程，即非線性建模后的CORS站高程殘差序列具有非平穩(wěn)性。GARCH(p，q)模型克服了AR模型只能對平穩(wěn)序列建模的不足，為下一步基于GARCH(p，q)模型對CORS站高程非平穩(wěn)殘差序列建模和非線性速度場重構(gòu)提供了思路。