試驗海區(qū)水文參數(shù)時間序列的多周期重構(gòu)預測方法

艾銳峰 程 杰 歐陽軍 孫云鵬

（解放軍63850部隊 白城 137001）

1 引言

在進行海上裝備試驗時，需要提前掌握海區(qū)水文要素（水溫、鹽度、海面風場、海流、波浪等）的變化規(guī)律，從而為試驗方案的制定提供決策依據(jù)。動力學的氣象海洋預報可以提供一定程度的參考［1］，但對部分水文參數(shù)（如海水溫度）能力有限。目前的處理手段是通過歷史數(shù)據(jù)的統(tǒng)計處理，提供海區(qū)水文參數(shù)不同時間周期上（年、月、日等）的平均值、方差、最大最小值等統(tǒng)計量作為未來值的估計。由于水文參數(shù)影響因素眾多、變化機理復雜、空變時變特性，此類估計量作為預測值相對粗糙。

更進一步的分析方法有回歸分析、馬爾可夫概型分析和時間序列分析等［2］?；貧w分析通過變量和因變量樣本數(shù)據(jù)的擬合，建立回歸模型，利用模型對未來時刻值作預報估計?；貧w分析需要變量、因變量之間確實存在物理上的關(guān)聯(lián)關(guān)系，且需要二者的同步樣本值，若只有本要素的時間測量序列，則無法進行回歸分析［3］。馬爾可夫概型分析以時間序列內(nèi)部概率分布結(jié)構(gòu)為基礎(chǔ)，建立馬爾可夫鏈，通過待分析參數(shù)狀態(tài)的逐步轉(zhuǎn)移，實現(xiàn)對參數(shù)時間演變過程的分析與預測［4］。它要求參數(shù)序列滿足無后效特點，適用場合有限［5］。基于ARIMA模型的時間序列分析方法是一種適用性更好的方法，它只基于待分析參數(shù)的歷史時間序列數(shù)據(jù)，建立參數(shù)化的模型，實現(xiàn)對未來時刻值的遞推預測［6～7］。

一種水文要素受到多種不同效力因素的影響，其最終呈現(xiàn)的時間序列吸收了所有影響因素的作用效應(yīng)，從而導致其可能是非平穩(wěn)的時間序列。當利用ARIMA模型分析法進行建模與預測時，雖然通過差分方法可以進行序列的平穩(wěn)化預處理，但當平穩(wěn)化程度有限時預測效果下降［8～9］。分析時間序列的影響因子，當單一因子作用時，產(chǎn)生的影響相對確定，且不同的影響因子可以認為作用于不同的尺度上。鑒于此，本文運用小波變換與ARIMA模型分析法，通過時間序列的多周期拆分、建模預測、反向組合的方式，設(shè)計了一種基于多周期重構(gòu)的預測方法。利用海洋浮標測量數(shù)據(jù)對該方法與傳統(tǒng)的ARIMA方法進行對照分析。

2 信號的多尺度分解與ARIMA建模

2.1 基于小波變換的信號多尺度分解原理

對水文參數(shù)時間序列的多周期拆分基于小波分析［10］。小波變換通過逐漸精細的取樣步長聚焦到信號的細節(jié)，從而在時域和頻域同時具有良好的局部化性質(zhì)。它相當于通過不同尺度參數(shù)的選擇將信號分解到對應(yīng)的尺度空間上去［11］。

設(shè) f(t)為空間L2(R)上的信號函數(shù)，則其小波變換為

式 （1） 中 φ(t) 為 母 小 波 函 數(shù) ，φa,b(t)= ||a-1/2φ((t-b)/a)是 φ(t)經(jīng)過伸縮平移變換得到的一系列基函數(shù)。其中a為尺度參數(shù)（或者伸縮參數(shù)），b為平移參數(shù)。通過小波逆變換可以重構(gòu)原信號函數(shù)，如式（2）所示：

基于小波變換可以實現(xiàn)對信號的多尺度分解。設(shè){φ (t)j,k}j,k∈Z是正交小波基函數(shù)，對應(yīng)于L2(R)空間的分解。其中，Z為整數(shù)。則可以將 f(t)展開為

式（3）中 cj,k=＜f(t),φj,k(t)＞ 為小波級數(shù)。令，則 fj(t)為信號 f(t)在子空間Vj上的投影。于是可以如圖1所示，將信號分解到不同的尺度空間上去。

圖1 信號的小波分解示意圖

2.2 時間序列的ARIMA建模方法

自回歸積分滑動平均模型（Autoregressive Integrated Moving Average Model，ARIMA）分析法是一種利用參數(shù)模型對隨機分布信號序列進行處理，從而估計出模型參數(shù)、建立遞推模型、進行數(shù)據(jù)序列的擬合與預測的方法［6］。其一般步驟：通過差分將非平穩(wěn)時間序列轉(zhuǎn)換為平穩(wěn)序列，再經(jīng)過模型類型識別、模型定階、模型參數(shù)估計建立ARMA模型，然后利用ARMA模型進行歷史數(shù)據(jù)序列的擬合或者對未來值進行遞推預測。

設(shè)某時間序列為 fn，n=1,2,…,N。若其非平穩(wěn)，則對其通過差分實現(xiàn)平穩(wěn)化處理，即令

f′n=fn-fn-1。其 ARMA(P,Q)模型為式（4）所示：

其中，模型參數(shù) λp、θq和白噪聲序列 αn的方差為待估計項。通過對序列進行相關(guān)性分析［7，8］，判斷適用的模型為AR、MA還是ARMA。模型階次(P,Q)的確定采用AIC準則［8］。模型參數(shù) λp、θq可通過最小二乘法解算得到［9］。

利用建立的ARMA(P,Q)模型可以實現(xiàn)對未來數(shù)據(jù)的遞推預測。設(shè)(P,Q)=(2,3)，利用預測。以k=1步預測為例，如式（5）所示，進行遞推求解：

3 時間序列的多周期重構(gòu)預測方法

3.1 時間序列的多周期重構(gòu)預測模型

以海面氣溫或水溫為例，設(shè) fn=tn，溫度受季節(jié)因素影響，溫度時間序列tn呈現(xiàn)季節(jié)變化的特征。一日當中晝夜太陽輻射變化，溫度時間序列tn呈現(xiàn)日變化的特征。同時溫度還受到其他各種因素影響［12］。所有這些因素共同作用導致所測量的溫度時間序列tn呈現(xiàn)隨機性和規(guī)律性并存的局面。ARMA模型分析方法可以實現(xiàn)對時間序列的參數(shù)化建模，并對未來值進行遞推預測。在理論上，當序列是平穩(wěn)隨機序列時，ARMA方法的性能較好；當序列非平穩(wěn)時，雖然可以通過差分，趨勢項剔除的方法對原序列進行平穩(wěn)化處理，但是平穩(wěn)化處理后的序列也并非完全平穩(wěn)。從物理機理上理解，假設(shè)tn只受季節(jié)變化的影響，則測量值tn為平穩(wěn)的周期序列，變化周期與季節(jié)變化同步；假設(shè)tn只受晝夜變化的影響，則tn為平穩(wěn)的周期序列，變化周期與晝夜變化同步。以此類推，可以假設(shè)tn受到m=1,2,…,M個因素影響，單個因素產(chǎn)生的序列ζm,n可以認為是平穩(wěn)序列，于是所觀測的序列tn可以用式（6）表示：

其中εn為剩余其他因素產(chǎn)生的序列。

一般而言，M個影響因素擁有不同的作用時間尺度?；诖思俣?，設(shè)有時間序列xn，n=1,2,…,N，構(gòu)建如下的時間序列多周期重構(gòu)預測模型，以實現(xiàn)利用序列xn預測xn+k。

首先令 xn=xm,n，其中m=1,2,…,M 為分解的次數(shù)。利用小波變換對xm,n進行分解，得到序列xm,n的細節(jié)部分[vm,1,n,vm,2,n,…vm,L,n]和概貌部分sm,n；剔除細節(jié)部分，取sm,n作周期為Tm的滑動平均，得到 ym,n；基于ARIMA方法，利用 ym,n建立模型 ARIMA(p,I,q)，利用 ARIMA(p,I,q)遞推預測y?m,n+k。同時，令 zm,n=xm,n-ym,n，對差序列 zm,n做下一步的分解，即令m=m+1，循環(huán)操作。最后將分解后的預測序列組合起來重構(gòu)為最終的預測序列 x?n+k。

圖2 時間序列多周期重構(gòu)預測模型

3.2 具體的預測方法

下面為具體的預測方法：

令m=1,2,…,M，即共進行M重分解預測操作。

1）取歷史數(shù)據(jù)序列 xn，令 xn=xm,n，m=1，作第一個周期的分解預測操作。首先，運用2.1節(jié)的方法，對xm,n進行小波分解，小波分解層數(shù)設(shè)為L。由分解得到的小波系數(shù)[cD1,cD2,…cDL,cA]逆變換得到序列 xm,n的細節(jié)部分[vm,1,n,vm,2,n,…vm,L,n]和概貌部分sm,n。

2）剔除細節(jié)部分，取sm,n作周期為Tm的滑動平均，得到 ym,n，即

做差序列zm,n=xm,n-ym,n，從而將原序列xm,n分解為主序列 ym,n和差序列zm,n。

3）令 m=m+1，對差序列 zm,n重新進行1）、2）的操作，進行進一步的分解。經(jīng)過M重分解后，可以將原序列 xn=x1,n分解為 y1,n,y2,n,…,yM,n,zM,n。即式（8）、（9）所示：

4）運用2.2節(jié)的方法，對 y1,n,y2,n,…,yM,n,zM,n分別進行ARIMA建模，利用建立的模型ARIMA(p,I,q)分別進行k步遞推預測 y?m,n+k。

5）將預測值組合起來，實現(xiàn)對xn+k的預測，即

4 數(shù)據(jù)分析

利用實際數(shù)據(jù)序列，將所構(gòu)建的時間序列多周期重構(gòu)預測方法（下面簡稱為MTC方法）與直接的ARIMA方法的預測性能進行比較。

按照3.2節(jié)的處理步驟，對某海洋浮標站的海面氣溫測量數(shù)據(jù)進行分析。浮標站所處位置為（123°E，38°N）。數(shù)據(jù)采集時間區(qū)間為2010年7月28日至2012年7月27日，采樣間隔為每小時一次，共計12100個數(shù)據(jù)。如圖3所示，將溫度時間序列分為兩部分，前20%作為樣本進行模型構(gòu)建，后80%作為模型遞推預測效果的比對。

圖3 海洋浮標溫度測量數(shù)據(jù)

差分操作可以對數(shù)據(jù)序列進行平穩(wěn)化處理，圖3（b）為溫度數(shù)據(jù)序列一次差分后的結(jié)果。根據(jù)MTC方法的處理流程，對溫度時間序列進行三重分解即 M=3，滑動平均周期分別為T1=720（月）、T2=168（周）、T1=24（日）。原始序列被分解為X1=Y1+Y2+Y3+Z3。

圖4為溫度時間序列第一重分解。其中X1為原始數(shù)據(jù)序列，S1為X1經(jīng)小波變換剔除細節(jié)部分后的序列，將之以周期T1=720滑動平均后的序列為Y1，Z1為X1減去Y1后的差序列。

圖4 溫度時間序列第一重分解

圖5 為將第一重分解后的差序列Z1（即X2）繼續(xù)分解后的處理結(jié)果。

圖5 溫度時間序列第二重分解

圖6 為將第二重分解后的差序列Z2（即X3）繼續(xù)分解后的處理結(jié)果。

根據(jù)上面的多周期拆分、再逐一建模遞推預測、再反向組合，得到圖7（a）的一步遞推預測結(jié)果，將之與直接ARIMA方法預測結(jié)果進行比較。

圖6 溫度時間序列第三重分解

圖7 一步遞推預測結(jié)果

圖8 為兩種方法的預測誤差比較，直接ARIMA方法預測誤差的標準差為1.6℃，多周期重構(gòu)方法預測誤差的標準差為0.72℃，性能得到提升。

圖8 一步遞推預測結(jié)果的誤差

5 結(jié)語

本文通過時間序列影響因子的分析，假定最終的時間序列由一個個單一因子作用下的序列組合而成，且各因子作用于不同的時間尺度?；诖思僭O(shè)，運用小波變換和滑動平均計算對原始序列進行若干周期的層層分解；再運用ARIMA建模方法對分解后的序列分別進行ARIMA模型構(gòu)建與遞推預測，通過預測結(jié)果的反向重構(gòu)實現(xiàn)對原始序列未來值的預測，從而構(gòu)建了一種時間序列預測的新方法。

實際時間序列數(shù)據(jù)的建模與預測結(jié)果表明，此方法相比直接的ARIMA方法，預測準確度有一定程度的提升?？蓪υ囼灪^(qū)水文參數(shù)序列變換采樣率，利用此模型對歷史數(shù)據(jù)進行分析，獲取海區(qū)水文參數(shù)的規(guī)律及未來時刻的預測值，為試驗決策提供參考。另外，在程序計算中發(fā)現(xiàn)，分解周期及次數(shù)的選擇、小波變換分層數(shù)的選擇均對該方法的預測效果有影響，需要在后續(xù)工作中進行深入的分析，研究其與實際水文參數(shù)影響因素及作用機理的關(guān)聯(lián)關(guān)系，形成相匹配的預測模型；其次，多步預測的效果也有待研究。