基于奇異譜分析的SPBO-ANFIS月徑流組合預(yù)測模型

張亞杰，崔東文

(1.云南省玉溪市易門縣水利局，云南 玉溪 651100；2.云南省文山州水務(wù)局，云南 文山 663000)

中長期徑流預(yù)報的準(zhǔn)確性和預(yù)見性對于區(qū)域水資源規(guī)劃與管理具有重要意義。由于徑流受氣候變化、地環(huán)境和人類活動的多重影響，呈現(xiàn)出較強(qiáng)的非線性、非平穩(wěn)性和多尺度等特征，使得預(yù)測難度大大增加，傳統(tǒng)單一模型或方法難以獲得滿意的預(yù)測效果[1-2]。研究表明，提高徑流時間序列預(yù)測精度的關(guān)鍵在于原始序列的有效分解和預(yù)測器的合理選取。時序數(shù)據(jù)分解方法有經(jīng)驗?zāi)B(tài)分解(EMD)、總體經(jīng)驗?zāi)B(tài)分解(EEMD)、小波分解(WD)、奇異譜分析(SSA)、小波包分解(WPD)等，均在時間序列預(yù)測研究中得到廣泛應(yīng)用；預(yù)測器有自回歸滑動平均模型[3]、人工神經(jīng)網(wǎng)絡(luò)模型[4-5]、支持向量機(jī)(SVM)[6]、相關(guān)向量機(jī)(RVM)[7]、LSTM神經(jīng)網(wǎng)絡(luò)[8-9]等。

奇異譜分析(Singular Spectrum Analysis，SSA)是近年來興起的一種研究非線性時間序列數(shù)據(jù)的有效方法，它根據(jù)所觀測到的時間序列構(gòu)造出軌跡矩陣，并對軌跡矩陣進(jìn)行分解、重構(gòu)，從而獲得更具規(guī)律的子序列分量，目前已在月徑流時間序列預(yù)報中得到應(yīng)用[6]。月徑流時間序列預(yù)測器主要有BP、SVM、LSTM等，BP神經(jīng)網(wǎng)絡(luò)存在設(shè)置參數(shù)多、易陷入局部最優(yōu)等缺點；SVM模型存在對參數(shù)敏感、大容量樣本預(yù)報中表現(xiàn)不佳等不足；LSTM模型預(yù)測性能較好，但存在內(nèi)存資源消耗大、運行時間長等缺陷。本文預(yù)測器選用自適應(yīng)神經(jīng)模糊推理系統(tǒng)(Adaptive Network based Fuzzy Inference System，ANFIS)。ANFIS具有較強(qiáng)的學(xué)習(xí)能力和表達(dá)能力，但其條件參數(shù)、結(jié)論參數(shù)的選取是制約ANFIS性能提高的關(guān)鍵因素，采用一種新型群體智能算法——學(xué)生心理學(xué)優(yōu)化(Student Psychology Based Optimization，SPBO)算法優(yōu)化ANFIS條件參數(shù)、結(jié)論參數(shù)，以期提高ANFIS預(yù)測器性能和智能化水平?；诖?，為有效提高徑流時間序列預(yù)測精度，克服單一模型的不足，基于SSA方法、SPBO算法、ANFIS系統(tǒng)，提出SSA-SPBO-ANFIS月徑流時間序列組合預(yù)測模型。主要內(nèi)容包括：①以云南省某水文站1960年1月至2012年12月共636個月月徑流預(yù)測為例，利用SSA對原始月徑流時間序列數(shù)據(jù)特征進(jìn)行充分提取，將原序列分解為若干獨立子序列；②簡要介紹SPBO算法，對驗證SPBO算法尋優(yōu)能力，選取8個典型單峰和多峰測試函數(shù)在不同維度條件下對SPBO算法的尋優(yōu)能力進(jìn)行仿真測試，并與當(dāng)前尋優(yōu)效果較好的教學(xué)優(yōu)化(LTBO)算法的尋優(yōu)結(jié)果進(jìn)行比較；③選取自適應(yīng)神經(jīng)模糊推理系統(tǒng)(ANFIS)作為月徑流時間序列預(yù)測模型，利用SPBO算法優(yōu)化ANFIS條件參數(shù)和結(jié)論參數(shù)，建立SSA-SPBO-ANFIS組合模型對分解后的各子序列進(jìn)行預(yù)測，疊加后作為最終月徑流預(yù)測結(jié)果，并構(gòu)建基于集合經(jīng)驗?zāi)B(tài)分解(EEMD)的EEMD-SPBO-ANFIS模型和未經(jīng)分解的SPBO-ANFIS模型作對比分析模型。

1 SSA-SPBO-ANFIS預(yù)測模型

1.1 奇異譜分析(SSA)

SSA是一種非線性時序數(shù)據(jù)處理方法，對于研究周期振蕩行為具有獨到的優(yōu)勢，其原理基于時間序列的動力重構(gòu)，是經(jīng)驗函數(shù)正交分解的一種特殊應(yīng)用，能夠有效地從包含噪聲的有限尺度時間序列中提取特征信息。SSA包含分解與重構(gòu)2個階段[10-13]。

a)分解。選取窗口長度L構(gòu)造軌跡矩陣X如下：

(1)

式中M——時間序列長度；L——窗口長度，即嵌入維數(shù)。

對矩陣X進(jìn)行SVD分解：

(2)

式中λi——第i個特征值；Ui——與第i個特征值對應(yīng)的特征向量；Vi——第i個主成分；d=rank(X)=max{i∶λi>0}。

(3)

式中 向量數(shù)量K=M-L+1；L*=min(L,K)；K*=max(L,K)。

1.2 學(xué)生心理學(xué)優(yōu)化(SPBO)算法

1.2.1SPBO算法簡述

學(xué)生心理學(xué)優(yōu)化(SPBO)算法是Bikash Das等人基于班級學(xué)生通過努力提高學(xué)習(xí)成績、爭當(dāng)優(yōu)等生心理而提出的一種新型群體智能優(yōu)化算法。該算法根據(jù)班級學(xué)生表現(xiàn)將學(xué)生分為最優(yōu)學(xué)生、好學(xué)生、普通學(xué)生和嘗試隨機(jī)提高學(xué)生4類，通過模擬此4類學(xué)生提高學(xué)習(xí)成績來達(dá)到求解優(yōu)化問題的目的。SPBO算法基于以下假設(shè)：①學(xué)生成績通過考試得分衡量，在考試中獲得最高分?jǐn)?shù)的學(xué)生被認(rèn)為是班級最優(yōu)學(xué)生；②班級中其他學(xué)生會通過努力來提高考試成績，并成為班級最優(yōu)學(xué)生，考試成績提高途徑因?qū)W生不同而各異；③要成為班級最優(yōu)學(xué)生，需要在所學(xué)科目上付出更多努力，任何一個學(xué)生對任何一門學(xué)科的努力都取決于學(xué)生的能力、效率以及對該學(xué)科的興趣；④除優(yōu)秀學(xué)生外，其余3類學(xué)生(好學(xué)生、普通學(xué)生和嘗試隨機(jī)提高學(xué)生)的學(xué)科選擇視為一個隨機(jī)過程[14]。

SPBO算法數(shù)學(xué)描述簡述如下[14]。

a)最優(yōu)學(xué)生。最優(yōu)學(xué)生始終通過努力獲得全班最高分來維持自己的位置。為獲得最高分并維持自己的位置，最優(yōu)學(xué)生需要對每個學(xué)科付出比其他學(xué)生更多的努力。最優(yōu)學(xué)生考試成績描述如下：

Xbestnew=Xbest+(-1)k×rand×(Xbest-Xj)

(4)

式中Xbestnew——當(dāng)前最優(yōu)學(xué)生某一學(xué)科考試成績；Xbest——最優(yōu)學(xué)生某一學(xué)科考試成績；Xj——第j名學(xué)生某一學(xué)科考試成績；rand——[0,1]范圍內(nèi)隨機(jī)數(shù)；k——隨機(jī)選擇為1或2的參數(shù)。

b)好學(xué)生。如果學(xué)生對某學(xué)科感興趣，那么她/他會在該學(xué)科上付出更多努力，從而使她/他的整體成績得到提高。這類學(xué)生的選擇是一個隨機(jī)過程，為通過在考試中獲得最高分而成為最優(yōu)學(xué)生，該類學(xué)生會嘗試做出與最優(yōu)學(xué)生相似或更多的努力。好學(xué)生考試成績描述如下：

(5)

式中Xnewi——當(dāng)前隨機(jī)選擇學(xué)生某一學(xué)科考試成績；Xi——第i名學(xué)生某一學(xué)科考試成績；Xmean——某一學(xué)科平均考試成績；其他參數(shù)意義同上。

c)普通學(xué)生。學(xué)生努力程度取決于學(xué)生對所學(xué)科目的興趣，如果該類學(xué)生對某學(xué)科不太感興趣，那么該類學(xué)生將為該科目付出平均努力。在平均努力學(xué)習(xí)該科目的同時，該類學(xué)生將嘗試在其他學(xué)科上付出更多努力來提高他們的整體成績。根據(jù)學(xué)生的不同心理，普通學(xué)生的選拔也是一個隨機(jī)過程。普通學(xué)生考試成績描述如下：

Xnewi=Xi+[rand×(Xmean-Xi)]

(6)

d)嘗試隨機(jī)提高學(xué)生。除上述3類學(xué)生外，有些學(xué)生試圖自己提高學(xué)習(xí)成績。該類學(xué)生嘗試隨機(jī)對科目進(jìn)行努力，以提高考試的整體成績。該類學(xué)生考試成績描述如下：

Xnewi=Xmin+[rand×(Xmax-Xmin)]

(7)

式中Xmax、Xmin——某一學(xué)科考試最高分、最低分。

1.2.2SPBO算法仿真驗證

選取Sphere等8個典型測試函數(shù)在不同維度條件下對SPBO算法進(jìn)行仿真驗證，并與當(dāng)前具有較好尋優(yōu)效果的TLBO算法的仿真結(jié)果進(jìn)行比較，見表1。設(shè)置SPBO、TLBO算法最大迭代次數(shù)T=500，群體規(guī)模N=50；其中TLBO算法參數(shù)TF為1～10之間隨機(jī)整數(shù)。

表1 函數(shù)優(yōu)化對比結(jié)果

對于單峰函數(shù)，SPBO算法尋優(yōu)精度較TLBO算法提高20個數(shù)量級以上。對于多峰函數(shù)Griewank，SPBO算法在低維(10維)條件尋優(yōu)效果遠(yuǎn)優(yōu)于TLBO算法，其他維度條件下與TLBO算法相同。對于Rastrigin函數(shù)，SPBO算法在不同維度條件下20次尋優(yōu)均獲得理論最優(yōu)值0，具有理想的尋優(yōu)效果。對于Penalized函數(shù)，SPBO算法在不同維度條件下尋優(yōu)精度均在4.71E-32以上，在低維(10維)條件下，SPBO算法尋優(yōu)精度略優(yōu)于TLBO算法；在高維條件下，SPBO算法尋優(yōu)精度遠(yuǎn)優(yōu)于TLBO算法。對于多峰函數(shù)Levy，SPBO算法在不同維度條件下尋優(yōu)精度均為1.50E-32，在低維(10維)條件下，SPBO算法尋優(yōu)精度略優(yōu)于TLBO算法；在高維條件下，SPBO算法尋優(yōu)精度遠(yuǎn)優(yōu)于TLBO算法。

可見，在上述條件下，SPBO算法尋優(yōu)效果優(yōu)于TLBO算法，具有較好的尋優(yōu)性能。

1.3 自適應(yīng)神經(jīng)模糊推理系統(tǒng)(ANFIS)

ANFIS兼顧了神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)機(jī)制和模糊系統(tǒng)語言推理能力的優(yōu)點，已在多個領(lǐng)域得到了成功應(yīng)用。ANFIS結(jié)構(gòu)表示見式(8)[15-19]：

(8)

式中x、y——輸入；A1、B1、A2、B2——模糊語言；a1、b1、c1——規(guī)則1的結(jié)論參數(shù)；a2、b2、c2——規(guī)則2的結(jié)論參數(shù)；f1、f2——規(guī)則輸出。

ANFIS推理系統(tǒng)通常由5層數(shù)學(xué)模型組成，在給定條件參數(shù)后，ANFIS輸出可以表示成結(jié)論參數(shù)的線性組合：

(9)

1.4 SSA-SPBO-ANFIS建模流程

步驟一利用SSA將原非平穩(wěn)月徑流時間序列分解為若干獨立子序列IMF；通過偏自相關(guān)函數(shù)法(PAFM)、自相關(guān)函數(shù)法(AFM)綜合確定各IMF輸入向量，合理劃分訓(xùn)練樣本和預(yù)測樣本。流程見圖1。

圖1 月徑流時間序列預(yù)測流程

步驟二利用訓(xùn)練樣本均方誤差構(gòu)建目標(biāo)函數(shù)：

(10)

步驟三初始化SPBO算法學(xué)生規(guī)模N，最大迭代次數(shù)Tmax和算法終止條件。

步驟四初評估班級每個學(xué)生的初始成績，根據(jù)初始成績將學(xué)生分為4類：最優(yōu)學(xué)生、好學(xué)生、普通學(xué)生和嘗試隨機(jī)提高學(xué)生。找到并保存當(dāng)前最優(yōu)學(xué)生，并令當(dāng)前迭代次數(shù)t=1。

步驟五判斷學(xué)生所屬類型：利用式(4)更新最優(yōu)學(xué)生成績；利用式(5)更新好學(xué)生成績；利用(6)更新普通學(xué)生成績；利用式(7)更新嘗試隨機(jī)提高學(xué)生成績。

步驟六評估班級每個學(xué)生的當(dāng)前成績，找到并保存當(dāng)前最優(yōu)學(xué)生。

步驟七令t=t+1。判斷終止條件，輸出全局最優(yōu)學(xué)生，算法結(jié)束；否則返回步驟五。

步驟八輸出最優(yōu)學(xué)生，即ANFIS最佳條件參數(shù)和結(jié)論參數(shù)。

步驟九將最優(yōu)ANFIS最佳條件參數(shù)和結(jié)論參數(shù)代入SPBO-ANFIS模型對各子序列分量進(jìn)行預(yù)測，將各分量預(yù)測結(jié)果加和重構(gòu)后即為月徑流最終預(yù)測結(jié)果。

步驟十采用平均絕對百分比誤差(MAPE)、平均絕對誤差(MAE)、納什系數(shù)(NSE)、合格率(PR)對模型有效性進(jìn)行評估，見式(11)。

(11)

2 實例應(yīng)用

2.1 數(shù)據(jù)來源

云南省某水文站控制徑流面積203 km2，屬國家基本站，觀測項目有水位、流量、降雨。本文研究數(shù)據(jù)來源于該站1960年1月至2012年12月共636個實測月徑流時間序列。從圖2可以看出，該站月徑流時間序列表現(xiàn)出強(qiáng)烈的非線性、非平穩(wěn)性和多層次性特征。

圖2 某水文站1960年1月至2012年12月徑流變化曲線

2.2 月徑流時間序列數(shù)據(jù)多尺度分解

2.2.1SSA分解

SSA是一種處理非線性時間序列數(shù)據(jù)的常用方法。研究表明，SSA分解確定子序列數(shù)量至關(guān)重要，若子序列數(shù)量過少，則不足以將原始序列中蘊含的不同成分提取出來；若子序列數(shù)量過多，則增加了模型復(fù)雜程度和建模工作量[19]。本文通過傅里葉變換(Fourier Transform，F(xiàn)T)方法，并結(jié)合月徑流顯著的年周特性，設(shè)置SSA窗口長度L=12，依據(jù)SSA原理對實例月徑流時間序列數(shù)據(jù)進(jìn)行多尺度分解，即將原月徑流時間序列分解為12個更具規(guī)律的子序列IMF，利用SPBO-ANFIS模型分別對12個IMF進(jìn)行預(yù)測，將預(yù)測結(jié)果疊加即得到最終月徑流預(yù)測結(jié)果。月徑流時間序列SSA分解結(jié)果見圖3。

a)IMF1

b)IMF2

c)IMF3

d)IMF4

e)IMF5

f)IMF6

g)IMF7

h)IMF8

i)IMF9

j)IMF10

k)IMF11

l)IMF12

2.2.2EEMD分解

利用EEMD方法對實例636個連續(xù)月徑流數(shù)據(jù)進(jìn)行分解。經(jīng)分解，原始月徑流數(shù)據(jù)分解為9個IMF分量和一個趨勢項rse。趨勢項rse大致反映原月徑流時間序列的變化趨勢；IMF1—IMF9反映了原徑流時間序列的波動情況。

2.3 確定輸入向量

合理確定各IMF輸入、輸出向量是決定月徑流時間序列預(yù)測精度高低的關(guān)鍵。目前普遍采用相空間重構(gòu)的方法確定時間序列輸入、輸出向量。相空間重構(gòu)目的是確定延遲時間和嵌入維數(shù)，目前延遲時間的確定方法有自相關(guān)函數(shù)法(AFM)等，嵌入維數(shù)確定方法有G-P法、C-C法等。研究表明，對于同一水文時間序列，不同方法確定的延遲時間和嵌入維數(shù)可能不同，甚至出差較大，這給確定時間序列輸入、輸出向量帶來困難。

經(jīng)過反復(fù)調(diào)試，采用偏自相關(guān)函數(shù)法(PAFM)或自相關(guān)函數(shù)法(AFM)綜合確定各IMF輸入向量，即將各IMF分量自相關(guān)系數(shù)或偏自相關(guān)系數(shù)最大時所對應(yīng)的滯后數(shù)H視為最優(yōu)嵌入維數(shù)，即將預(yù)測月前H個徑流數(shù)據(jù)作為輸入向量，預(yù)測月作為輸出向量。各IMF輸入、輸出向量確定結(jié)果見表2。本文將實例后120個月作為預(yù)測樣本。

EEMD分解得到的9個IMF分量和1個趨勢項rse的輸入、輸出向量采用同樣的方法確定，限于篇幅從略。

2.4 參數(shù)設(shè)置及預(yù)測分析

2.4.1參數(shù)設(shè)置

SSA-SPBO-ANFIS、EEMD-SPBO-ANFIS、SPBO-ANFIS模型：SPBO算法學(xué)生規(guī)模N=50，最大迭代次數(shù)T=100；ANFIS模型條件參數(shù)和結(jié)論參數(shù)搜索范圍為[-25,25]，初始聚類數(shù)目5，分類矩陣指數(shù)5，最大迭代次數(shù)100，目標(biāo)誤差0.001。

2.4.2預(yù)測分析

利用SSA-SPBO-ANFIS、EEMD-SPBO-ANFIS模型對實例各子序列進(jìn)行預(yù)測，將預(yù)測結(jié)果加和重構(gòu)后得到實例最終月徑流預(yù)測結(jié)果；并利用SPBO-ANFIS模型對未經(jīng)分解的原月徑流時間序列進(jìn)行訓(xùn)練及預(yù)測。同時，采用平均絕對百分比誤差MAPE(%)、平均絕對誤差MAE(m3/s)、納什系數(shù)NSE和合格率PR(%)對各模型性能進(jìn)行評估，結(jié)果見表3，預(yù)測效果見圖4。

表2 各IMF相關(guān)系數(shù)、嵌入維數(shù)及序列長度

表3 實例月徑流時間序列預(yù)測結(jié)果對比

a)徑流量

b)相對誤差

c)絕對誤差

從表3、圖4可以得出以下結(jié)論。

a)SSA-SPBO-ANFIS模型對實例月徑流時間序列預(yù)測的MAPE、MAE、NSE和PR分別為5.57%、0.20 m3/s、0.994 0、96.7%，MAPE分別較EEMD-SPBO-ANFIS、SPBO-ANFIS模型降低68.7%和90.1%，MAE分別降低62.3%和87.1%，NSE分別提高2.1%和40.0%，PR分別提高38.1%和223.3%，擬合、預(yù)測效果優(yōu)于EEMD-SPBO-ANFIS、SPBO-ANFIS模型，具有更小的預(yù)測誤差和更高的預(yù)測精度；SSA能科學(xué)將原徑流時間序列數(shù)據(jù)分解成多個更具規(guī)律的子序列，顯著提高水文時間序列月徑流的預(yù)測精度。

b)從表3來看，在長達(dá)120個月的月徑流預(yù)測中，SSA-SPBO-ANFIS、EEMD-SPBO-ANFIS模型預(yù)測的納什系數(shù)NSE分別為0.994 0和0.973 7，高于SPBO-ANFIS模型預(yù)測的納什系數(shù)0.709 8，表明SSA-SPBO-ANFIS、EEMD-SPBO-ANFIS模型預(yù)測結(jié)果可靠。其中，尤以SSA-SPBO-ANFIS模型的預(yù)測效果最好，預(yù)測可信度最高。

c)從圖4來看，SSA-SPBO-ANFIS模型的預(yù)測效果最好，EEMD-SPBO-ANFIS模型次之，SPBO-ANFIS模型最差。SSA-SPBO-ANFIS模型預(yù)測結(jié)果能夠很好地逼近實測徑流，對峰值捕捉效果好，具有較小的預(yù)測誤差和較好的預(yù)測效果，將其用于月徑流時間序列預(yù)測是可行的。

3 結(jié)論

基于奇異譜分解(SSA)、學(xué)生心理學(xué)優(yōu)化(SPBO)算法和自適應(yīng)神經(jīng)模糊推理系統(tǒng)(ANFIS)構(gòu)建水文時間序列月徑流組合預(yù)測模型，利用云南省某水文站月徑流時間序列預(yù)測實例對SSA-SPBO-ANFIS、EEMD-SPBO-ANFIS、SPBO-ANFIS模型進(jìn)行檢驗；并通過對SPBO算法仿真測試和比較MAPE、MAE、NSE、PR 4個指標(biāo)，可以得出以下結(jié)論。

a)介紹了學(xué)生心理學(xué)優(yōu)化(SPBO)算法。通過8個單峰、多峰測試函數(shù)在不同維度條件下對SPBO算法進(jìn)行仿真驗證。結(jié)果顯示SPBO算法尋優(yōu)性能優(yōu)于TLBO算法，將SPBO算法用于ANFIS條件參數(shù)和結(jié)論參數(shù)尋優(yōu)具有可靠性。

b)SSA-SPBO-ANFIS模型對實例月徑流時間序列擬合、預(yù)測效果優(yōu)于EEMD-SPBO-ANFIS模型，遠(yuǎn)優(yōu)于SPBO-ANFIS模型。模型有效地提高了月徑流時間序列的預(yù)測精度，且具有較好的通用性與穩(wěn)定性，為月徑流預(yù)測提供了一種新方法。

c)實例驗證表明，SSA能有效將復(fù)雜非線性和非平穩(wěn)性的月徑流時序數(shù)據(jù)分解為若干更具規(guī)律的子序列，抽取出月徑流時間序列的整體趨勢和不同周期上的波動情況，同時挖掘月徑流時序數(shù)據(jù)的物理特征和結(jié)構(gòu)信息，提高了模型的預(yù)測精度。

d)利用SPBO算法優(yōu)化ANFIS條件參數(shù)和結(jié)論參數(shù)，不但克服ANFIS模型隨機(jī)選取參數(shù)的不足，而且有效提高了模型的預(yù)測精度和智能化水平，可為ANFIS相關(guān)預(yù)測研究提供參考。