張思慧，黃聲和

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

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

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

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

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

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

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

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

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

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

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

2.3 最小差距法

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

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

3.1 譜估計(jì)

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

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

根據(jù)MATLAB12.0的計(jì)算結(jié)果可知，北京房山站高程時(shí)間序列的平均譜功率為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.由此可知，北京房山站高程時(shí)間序列的年周期性明顯，半年周期性和季周期性不是很明顯且波動(dòng)比較大.這是因?yàn)槭艿酱箨懓鍓K運(yùn)動(dòng)、地下水周期性變化、地殼巖漿運(yùn)動(dòng)和太陽黑子周期性活動(dòng)等因素的影響.

3.2 周期擬合

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

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

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

圖3 時(shí)間序列濾波和擬合

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

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

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

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

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

本文的BP 網(wǎng)絡(luò)選用54 個(gè)輸入神經(jīng)元和1 個(gè)輸出神經(jīng)元，第一層(隱層)傳遞函數(shù)選用tansig()，第二層(輸出層)是單個(gè)神經(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時(shí)，隱層節(jié)點(diǎn)的輸出為

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

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

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

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

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

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

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

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

閾值修正為

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

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

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

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

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

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

5 結(jié)論

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

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