采用改進(jìn)神經(jīng)網(wǎng)絡(luò)的CORS站高程混沌時間序列分析

張思慧，黃聲和

（1.福建船政交通職業(yè)學(xué)院 土木工程學(xué)院，福建 福州 350007；2.福建省特種設(shè)備檢驗(yàn)研究院，福建 福州 350008）

隨著我國IGS（國際衛(wèi)星全球定位服務(wù)）網(wǎng)絡(luò)不斷發(fā)展，CORS站（連續(xù)運(yùn)行參考站）的觀測數(shù)據(jù)逐步應(yīng)用到各個領(lǐng)域.例如，數(shù)字中國、地殼運(yùn)動監(jiān)測和導(dǎo)航等.但是由于CORS基準(zhǔn)站所處的觀測環(huán)境復(fù)雜多變，因此觀測到的高程時間序列含有“噪聲”.含有“噪聲”的高程時間序列在未經(jīng)處理的情況下，無法滿足應(yīng)用要求.因此對高程時間序列的數(shù)據(jù)處理顯得尤為重要.

目前對于高程時間序列常用的處理方法有：LFS法(最小二乘法)、Spline插值（三次樣條插值）和功率譜估計(jì)等[1-6].采用以上數(shù)據(jù)處理方法的前提是所處理的時間序列必須滿足Gauss-Markov 條件或正態(tài)性(平穩(wěn)性)條件.然而北京房山觀測站所處的環(huán)境復(fù)雜多變，其得到的觀測數(shù)據(jù)可能無法滿足Gauss-Markov 條件或正態(tài)性.若不滿足上述條件，直接對數(shù)據(jù)進(jìn)行擬合、濾波和譜估計(jì)等處理，其結(jié)果將嚴(yán)重失真.因此，在進(jìn)行擬合、濾波和譜分析等數(shù)據(jù)處理之前，必須對數(shù)據(jù)進(jìn)行穩(wěn)定性檢驗(yàn)、非線性檢驗(yàn)和系統(tǒng)確定性檢驗(yàn).由文獻(xiàn)[11]得出北京房山站高程時間序列為部分確定性機(jī)制低階混沌動力系統(tǒng).

1 數(shù)據(jù)來源及數(shù)據(jù)處理流程

采用GAMIT/GLOBK10.35 軟件對我國29 個CORS 基準(zhǔn)站與國際IGS 基準(zhǔn)站進(jìn)行聯(lián)合解算，得到北京房山站1999—2009 年單天解算時間序列，作為本文的數(shù)據(jù)來源.對高程時間序列進(jìn)行正則化RBF（徑向基）神經(jīng)網(wǎng)絡(luò)相空間重構(gòu)，使其滿足Gauss-Markov 條件，然后采用小波神經(jīng)網(wǎng)絡(luò)對高程時間序列進(jìn)行濾波分析，對消除“噪聲”的數(shù)據(jù)進(jìn)行加窗譜估計(jì)，然后采用BP 神經(jīng)網(wǎng)絡(luò)對北京房山站高程時間序列進(jìn)行預(yù)測，其數(shù)據(jù)處理流程圖如圖1所示.

圖1 數(shù)據(jù)處理流程圖

2 正則化RBF網(wǎng)絡(luò)和小波神經(jīng)網(wǎng)絡(luò)濾波

2.1 正則化RBF神經(jīng)網(wǎng)絡(luò)時間序列重構(gòu)

依據(jù)文獻(xiàn)[11]的研究結(jié)果可知，在對北京房山站所觀測到的數(shù)據(jù)進(jìn)行進(jìn)一步利用前，需要對其進(jìn)行相空間重構(gòu).結(jié)合北京房山站高程時間序列的數(shù)據(jù)特點(diǎn)，本文采用正則化RBF神經(jīng)網(wǎng)絡(luò)進(jìn)行相空間重構(gòu)[10-15].數(shù)據(jù)模型選用54個輸入神經(jīng)元、2個隱含層和1個輸出神經(jīng)元.

2.2 小波神經(jīng)網(wǎng)絡(luò)濾波

基于房山站數(shù)據(jù)量和隱節(jié)點(diǎn)數(shù)，設(shè)計(jì)正則化RBF 網(wǎng)絡(luò)的輸入層神經(jīng)元個數(shù)為30個，隱層為2個，輸出神經(jīng)元為1個.隱層選取Morlet函數(shù)作為小波神經(jīng)元的激勵函數(shù).

2.3 最小差距法

最小差距法能滿足擬合誤差平方和最小原則，并且在用于擬合不滿足Gauss-Markov 條件的線性模型時，仍具有無偏性、最小方差性等良好的統(tǒng)計(jì)特性.結(jié)合北京房山高程時間序列的特性，本文采用最小差距法對數(shù)據(jù)進(jìn)行5階擬合.擬合所得函數(shù)為：

3 譜估計(jì)、周期擬合和結(jié)果分析

3.1 譜估計(jì)

通過正則化RBF 和小波神經(jīng)網(wǎng)絡(luò)對北京房山站高程時間序列進(jìn)行剔除趨勢項(xiàng)處理，但其含有的周期項(xiàng)依然存在.為了進(jìn)一步分析高程時間序列的數(shù)據(jù)特性，本文結(jié)合上述的數(shù)據(jù)處理結(jié)果，采用哈明窗、巴特萊特窗和Parzen窗譜估計(jì)進(jìn)行周期譜估計(jì).時間序列周期譜估計(jì)如圖2所示.

圖2 時間序列加窗譜估計(jì)

根據(jù)MATLAB12.0的計(jì)算結(jié)果可知，北京房山站高程時間序列的平均譜功率為7.684×103W.2000—2003 年的年周期為1.003～1.034 a，半年周期為0.438～0.526 a，季周期為0.253～0.347 a；2004—2007 年的年周期為0.994～1.109 a，半年周期為0.478～0.617 a，季周期為0.247～0.378 a.由此可知，北京房山站高程時間序列的年周期性明顯，半年周期性和季周期性不是很明顯且波動比較大.這是因?yàn)槭艿酱箨懓鍓K運(yùn)動、地下水周期性變化、地殼巖漿運(yùn)動和太陽黑子周期性活動等因素的影響.

3.2 周期擬合

采用快速傅里葉變換方法對北京房山高程時間序列進(jìn)行周期擬合，結(jié)果如下：

3.3 模型計(jì)算結(jié)果分析

根據(jù)上述2.1 和2.2 構(gòu)建的數(shù)學(xué)模型，采用MATLAB12.0 軟件對北京房山站2000—2006 年觀測到的高程時間序列進(jìn)行濾波和擬合計(jì)算，其計(jì)算結(jié)果如圖3所示.

圖3 時間序列濾波和擬合

擬合殘差和去除趨勢FFT周期擬合函數(shù)結(jié)果如圖4所示.

圖4 5階最小差距擬合殘差和去除趨勢FFT周期擬合

由圖3 可知，2000—2001 年北京房山站觀測到的數(shù)據(jù)傾斜趨勢波動較大，但是2002—2003 年北京房山站高程時間序列呈現(xiàn)緩慢增長趨勢，2004 年甚至出現(xiàn)了負(fù)傾斜趨勢，2005 年北京房山站高程時間序列急劇增長，這表明這段時間北京房山站所處位置的地殼運(yùn)動較為劇烈.北京房山站觀測數(shù)據(jù)未能呈現(xiàn)單一的變化趨勢，其原因是高程時間序列不是受單一因素的影響，可能受到大陸板塊運(yùn)動、地下水周期性變化和地殼巖漿運(yùn)動等因素中的兩個或多個因素的影響.圖4 表明2000—2006 年北京房山站高程時間序列的年周期為0.997～1.037 a，半年周期為0.451～0.629 a，且存在較大波動，季周期為0.232～0.269 a，且有較大波動.由此可知，北京房山站高程時間序列存在明顯的年周期性、半年周期性和不明顯的季節(jié)周期性.

4 BP神經(jīng)網(wǎng)絡(luò)時間序列預(yù)測

BP網(wǎng)絡(luò)預(yù)測模型只能處理長期記憶性數(shù)據(jù)，無法精確預(yù)測短期記憶性(非平穩(wěn))數(shù)據(jù).但是經(jīng)過上述正則化RBF、小波神經(jīng)網(wǎng)絡(luò)處理后，北京房山站高程時間序列轉(zhuǎn)化為穩(wěn)態(tài)數(shù)據(jù)，進(jìn)而可以進(jìn)行BP網(wǎng)絡(luò)預(yù)測.

本文的BP 網(wǎng)絡(luò)選用54 個輸入神經(jīng)元和1 個輸出神經(jīng)元，第一層(隱層)傳遞函數(shù)選用tansig()，第二層(輸出層)是單個神經(jīng)元選用purelin()函數(shù)，訓(xùn)練函數(shù)選用traingd().BP 網(wǎng)絡(luò)的學(xué)習(xí)規(guī)則選用負(fù)梯度下降函數(shù)，其負(fù)梯度下降函數(shù)[16-17]為

式中：xk、xk+1為節(jié)點(diǎn)前后的閾值矩陣權(quán)值；ηk為學(xué)習(xí)速率；gk為當(dāng)前函數(shù)梯度.

1）BP算法模型節(jié)點(diǎn)假設(shè)三層BP網(wǎng)絡(luò)，輸入節(jié)點(diǎn)為xi，隱層節(jié)點(diǎn)為yi，輸出節(jié)點(diǎn)為zi.輸入節(jié)點(diǎn)與隱層節(jié)點(diǎn)間的網(wǎng)絡(luò)權(quán)值為wji，隱層節(jié)點(diǎn)與輸出節(jié)點(diǎn)間的網(wǎng)絡(luò)權(quán)值為vji，θj、θl為閾值.當(dāng)輸出節(jié)點(diǎn)的期望值為tl時，隱層節(jié)點(diǎn)的輸出為

輸出節(jié)點(diǎn)為

2）誤差函數(shù)對輸出節(jié)點(diǎn)求導(dǎo)：

Ek是多個zk的函數(shù)，但只有一個zl與vlj有關(guān)，各zk間相互獨(dú)立，則有

3）誤差函數(shù)對隱層節(jié)點(diǎn)求導(dǎo)：

El是多個zl的函數(shù)，針對某一個yj，它與所有zl有關(guān)，則有

由于權(quán)值的修正Δvlj、Δwji正比于誤差函數(shù)，沿梯度下降，則有

式中：δl為節(jié)點(diǎn)誤差；η為學(xué)習(xí)速率.

4）閾值的修正.閾值θ是變化值，誤差函數(shù)對輸出節(jié)點(diǎn)閾值求導(dǎo)為

閾值修正為

5）S型傳遞函數(shù)為

6）AR（自回歸模型）預(yù)測.BP神經(jīng)網(wǎng)絡(luò)AR預(yù)測數(shù)學(xué)結(jié)構(gòu)為

式中：y(t)～y(t-na)為時間t的輸出；na為節(jié)點(diǎn)個數(shù)；an為待定參數(shù)；e(t)為白噪聲輸入數(shù)據(jù).

將經(jīng)過模型處理后的高程時間序列代入建好的BP神經(jīng)網(wǎng)絡(luò)模型進(jìn)行預(yù)測.2003—2007年高程時間序列預(yù)測結(jié)果的殘差如圖5所示.

圖5 BP網(wǎng)絡(luò)預(yù)測殘差圖

由圖5 可知，構(gòu)建的BP 網(wǎng)絡(luò)模型預(yù)測殘差在0～0.035 m 之間，并且通過計(jì)算其相對誤差為-0.028 6%～0.040 1%.通過加窗譜估計(jì)對BP網(wǎng)絡(luò)預(yù)測的2003—2007年的高程時間序列進(jìn)行計(jì)算.計(jì)算結(jié)果顯示，其存在著年周期性、半年周期性和季周期性，年周期為1.005～1.097 a，半年周期為0.561～0.613 a，季周期為0.389～0.371 a，這與2003—2007年的原始時間序列函數(shù)關(guān)系基本一致.由此可知，該模型預(yù)測精度非常高，BP網(wǎng)絡(luò)可以很好地對時間系列進(jìn)行預(yù)測.

5 結(jié)論

1）基于正則化RBF 對北京高程時間序列進(jìn)行相空間重構(gòu)，然后對相空間重構(gòu)的時間序列進(jìn)行小波濾波，對濾波后的時間序列進(jìn)行加窗譜估計(jì).由降噪后的功率譜圖可知，北京房山站高程時間序列具有明顯的周期性，且具有諧波性.

2）采用BP網(wǎng)絡(luò)對2003—2007年高程時間序列進(jìn)行預(yù)測.結(jié)果顯示，其相對誤差為-0.028 6%～0.040 1%；采用BP預(yù)測2003—2007年的周期性與2003—2007年的原始時間序列函數(shù)關(guān)系基本一致.預(yù)測精度非常高，預(yù)測效果良好.