謝智峰，吳潤秀，呂 莉

(1.江西科技學院 信息工程學院，江西 南昌 330098；2.南昌工程學院 信息工程學院，江西 南昌 330099)

徑流預測是根據(jù)徑流形成規(guī)律，借助已知信息對未來某一段時間內(nèi)的徑流量做出定性或定量的預測，其在防汛抗旱、水資源配置和水庫調(diào)度等方面發(fā)揮著重要作用[1]。徑流形成受太陽活動、大氣環(huán)流和人類活動等諸多因素影響，因此徑流預測屬于受多因素影響的非線性復雜問題。目前，徑流預測的方法主要有小波分析(Wavelet Analysis，WA)[2]、時間序列分析(Time Series Analysis，TSA)[3]、人工神經(jīng)網(wǎng)絡(Artificial Neural Networks，ANN)[4]和支持向量回歸(Support Vector for Regression，SVR)[5-7]等。其中，支持向量回歸因在解決非線性復雜問題方面表現(xiàn)突出，已成為水文預測領域的常用方法之一。

群智能算法(Swarm Intelligence Algorithm，SIA)是一種新興的仿生進化算法，具有模型簡單、實現(xiàn)容易、魯棒性強等特點。典型的群智能算法包括粒子群算法(Particle Swarm Optimization，PSO)[8-9]、人工蜂群算法(Artificial Bee Colony，ABC)[10-11]和螢火蟲算法(Firefly Algorithm，F(xiàn)A)[12-14]等。由于這類算法不受優(yōu)化問題的連續(xù)或可微的限制，已廣泛應用于無線傳感器優(yōu)化[15-17]、圖像壓縮[18-19]和徑流預測[20-22]等領域。

為獲取最佳的SVR核參數(shù)組合，很多學者將群智能算法應用于SVR的參數(shù)尋優(yōu)，不僅提高了參數(shù)尋優(yōu)效率，而且獲得了較好的泛化性能。王宏偉[23]等將遺傳算法(Genetic Algorithm，GA)嵌入SVR，構建了GA-SVR預報模型，取得了較好的預報效果；暢明琦[24]等設計了基于不同預報因子的兩組SVR模型，實驗結果表明，多預報因子的SVR模型預測效果更佳；蔣林利[25]等提出基于粒子群和模擬退火的混合算法，實現(xiàn)SVR參數(shù)優(yōu)化選取，實驗驗證了混合算法模型具備較好的泛化性能，預測精度高且穩(wěn)定。

Yang[26]于2008年提出螢火蟲算法。該算法思想來源于自然界中發(fā)光螢火蟲向比自身更亮的螢火蟲飛行的生物習性。螢火蟲離散分布于一定的活動范圍內(nèi)，每只螢火蟲因體內(nèi)不同量的熒光素而發(fā)出強弱不一的光，隨著時間推移，該范圍內(nèi)的螢火蟲會逐漸聚集于最亮螢火蟲附近，形成若干個相似亮度的聚集中心，類似于尋優(yōu)過程。

為兼顧種群全局探索和局部開采的平衡，提高算法求解精度，提出多策略融合學習螢火蟲算法(Multi-Strategy Fusion Learning Firefly Algorithm，MSFLFA)。首先，算法為最優(yōu)螢火蟲隨機選擇兩只不同螢火蟲作為學習對象，制定新更新方式，進行定量次數(shù)的進化，提高算法探索全局最優(yōu)解的能力；然后，在非最優(yōu)螢火蟲的既有更新策略上添加新的擾動項，動態(tài)識別搜索域，增強算法的局部尋優(yōu)能力；最后，不同螢火蟲交替迭代進化，完成尋優(yōu)。為驗證MSFLFA性能，設計基準測試函數(shù)的仿真實驗，結果表明，MSFLFA尋優(yōu)精度高且綜合性能優(yōu)于文中涉及的比較算法。

1 支持向量回歸和螢火蟲算法

1.1 支持向量回歸

支持向量回歸常用的核函數(shù)[27]主要包括線性核函數(shù)、sigmoid 核函數(shù)、多項式核函數(shù)和徑向基核函數(shù)。線性核函數(shù)適合線性可分問題；sigmoid 核函數(shù)類似于多層神經(jīng)網(wǎng)絡；多項式核函數(shù)和徑向基核函數(shù)均可將低維非線性問題轉(zhuǎn)化為高維線性問題。對于受多因子影響的非線性徑流預報問題，線性核函數(shù)明顯不適用；基于sigmoid 核函數(shù)的SVR回歸預報結果不理想；徑向基核函數(shù)相較于多項式核函數(shù)涉及的參數(shù)少，計算復雜度低?；趶碗s度和易操作性原則，本文選取徑向基核函數(shù)作為SVR核函數(shù)。

SVR模型可抽象為如下形式：

y=f(x|C，ε，σ)，

(1)

其中C為懲罰系數(shù)；ε為不敏感損失系數(shù)；σ為徑向基核參數(shù)。

1.2 螢火蟲算法

根據(jù)螢火蟲算法尋優(yōu)原理，建立算法優(yōu)化模型，其相關數(shù)學表達式如下：

定義1螢火蟲的光亮強度I：

I=I0e-γrij.

(2)

定義2螢火蟲的吸引度β：

β=β0e-γrij，

(3)

其中I0和β0分別表示螢火蟲的固有光亮強度和吸引度；γ表示光吸收系數(shù)；rij表示螢火蟲i和螢火蟲j之間的歐氏距離，見定義3。

定義3螢火蟲i和螢火蟲j的歐氏距離rij：

(4)

其中xi和xj分別表示螢火蟲i和螢火蟲j的搜索空間位置；xid(t)和xjd(t)分別為螢火蟲i和j的第t代第d維位置；D表示搜索空間維度。

定義4螢火蟲i向螢火蟲j的學習后的位置xid(t+1)：

xid(t+1)=xid(t)+β(xjd(t)-xid(t))+αi(t)ε，

(5)

其中αi(t)表示螢火蟲i的第t代步長因子；ε服從均勻分布，取值范圍為[-0.5，0.5]。

2 多策略融合學習螢火蟲算法

在標準螢火蟲算法中，每只螢火蟲均向其余螢火蟲依次學習，采用這種全吸引模型更新方式，既損耗更新資源，又易使種群出現(xiàn)振蕩現(xiàn)象，導致種群過早收斂，陷入局部最優(yōu)[28]。為有效利用更新資源，避免種群“早熟”，本文為不同螢火蟲制定了2種學習策略。文獻[29]為廣義中心粒子設計了深度學習策略，在優(yōu)勝劣汰的評估機制下，不斷挖掘更為優(yōu)異的廣義中心粒子，引導種群進化。為此，本文將深度學習策略應用于最優(yōu)螢火蟲，提高最優(yōu)螢火蟲的質(zhì)量，加快完成全局挖掘最優(yōu)解任務；非最優(yōu)螢火蟲采用新添加的領域搜索項更新模式，在隨機吸引模型[28]的指導下，完成局部開采種群優(yōu)勢信息任務。最優(yōu)螢火蟲和非最優(yōu)螢火蟲分別采用不同的學習策略，構成一個進化單元，迭代尋優(yōu)。

2.1 最優(yōu)螢火蟲學習策略

最優(yōu)螢火蟲隨機選擇兩只不同的螢火蟲作為學習對象，根據(jù)式(6)進行Deep_count次深度學習，使其最大限度挖掘全局最優(yōu)信息。

(6)

種群進化前期，種群處于發(fā)散狀態(tài)，螢火蟲之間位置相距較遠，較大的搜索步長，可有效提高種群全局探索最優(yōu)解能力；種群進化后期，種群處于收斂狀態(tài)，螢火蟲位置相距較近，較小的搜索步長，可有效提高種群局部開采最優(yōu)解能力。為此，通過柯西函數(shù)動態(tài)調(diào)整步長，并采用貪婪更新策略評估學習狀態(tài)。評估結果為優(yōu)，則更新狀態(tài)；反之，保留原有狀態(tài)。最優(yōu)螢火蟲學習策略定義如下：

(1)記錄最優(yōu)螢火蟲當前位置XBest和適應值f(XBest)；

(3)按照式(6)更新最優(yōu)螢火蟲當前位置；

(4)評估最優(yōu)螢火蟲，結果為優(yōu)則記錄新位置new_XBest和新適應值new_f(XBest)，反之保留更新之前的位置及適應值，同時記錄評估次數(shù)FEs=FEs+1；

(5)判斷最優(yōu)螢火蟲是否完成Deep_count次深度學習，完成跳轉(zhuǎn)至(6)，反之跳轉(zhuǎn)至步驟(1)；

(6)記錄最優(yōu)螢火蟲最終位置XBest和適應值f(XBest)，最優(yōu)螢火蟲更新完成。

2.2 非最優(yōu)螢火蟲學習策略

位置更新公式(5)由3項組成，第1項表示螢火蟲Xi當前位置，第2項表示螢火蟲Xi向螢火蟲Xj移動的距離，第3項表示螢火蟲Xi完成學習任務后，以自身為中心，搜索附近最優(yōu)信息移動的步長。在整個種群進化過程中，前兩項組合保證了螢火蟲算法特有的學習方式，而第3項直至種群尋優(yōu)任務結束，螢火蟲始終做輕微擾動，無法根據(jù)當前搜索狀態(tài)，有效地適應不同搜索域，易陷入局部最優(yōu)，導致種群出現(xiàn)停滯。

為有效控制螢火蟲在更新過程中陷入局部最優(yōu)，同時適應不同的搜索空間，本文將式(7)作為新位置更新方式，以提高螢火蟲逃離局部最優(yōu)能力。

xid(t+1)=xid(t)+β(xjd(t)-xid(t)) +αi(t)ε(xjd(t)-xid(t)).

(7)

式(7)第3項在原有基礎之上乘以兩只螢火蟲的實時位置差，未額外添加參數(shù)。在進化前期，螢火蟲i和螢火蟲j之間的距離較大，螢火蟲j對螢火蟲i吸引力很弱，第2項使得螢火蟲i向螢火蟲j移動距離較小，學習程度不夠，為此，第3項結合螢火蟲i和螢火蟲j的實際距離，增大學習步長，增強學習程度，探索能力顯著加強；在進化后期，兩螢火蟲之間相距很近，第2項和第3項值均相對較小，有效保證螢火蟲進行局部開采；因此，式(7)動態(tài)地兼顧了種群在整個搜索過程的探索和開采能力。同時，該更新策略適用于不同搜索域，有效地提高了算法的泛化能力。

最優(yōu)螢火蟲完成任務之后，非最優(yōu)螢火蟲啟動更新進程。為了充分利用更新資源，非最優(yōu)螢火蟲使用隨機吸引模型替代全吸引模型，且每只螢火蟲使用式(7)完成位置的搜索。至此，由最優(yōu)和非最優(yōu)螢火蟲學習策略組合構成種群進化單元，不同螢火蟲按照進化單元既定規(guī)則協(xié)作迭代，完成種群尋優(yōu)任務。非最優(yōu)螢火蟲學習策略定義如下：

(1)記錄除最優(yōu)螢火蟲外N-1個螢火蟲當前位置Xi和適應值f(Xi)；

(2)設置隨機吸引模型作為N-1個螢火蟲位置更新方式；

(3)按照式(7)更新N-1個螢火蟲當前位置；

(4)評估最優(yōu)螢火蟲，記錄新位置new_Xi和新適應值new_f(Xi)，同時記錄評估次數(shù)FEs=FEs+1；

(5)N-1個螢火蟲完成更新任務并記錄最終位置Xi和適應值f(Xi)，非最優(yōu)螢火蟲更新完成。

2.3 算法流程

本文將評估次數(shù)作為算法終止條件。其中，種群大小為N，總評估次數(shù)為MAX_FEs，當前評估次數(shù)為FEs。MSFLFA流程如下：

Step1：初始化種群，記錄FEs=N；

Step2：執(zhí)行最優(yōu)螢火蟲學習策略設定的進化流程；

Step3：執(zhí)行非最優(yōu)螢火蟲學習策略設定的進化流程；

Step4：若FEs>MAX_FEs，則跳轉(zhuǎn)至Step5，反之，跳轉(zhuǎn)至Step2；

Step5：輸出種群最優(yōu)值和最優(yōu)解。

3 實驗仿真及分析

3.1 基準測試函數(shù)及參數(shù)設置

在仿真實驗中，采用了經(jīng)典基準測試函數(shù)[30]對MSFLFA進行測試。測試集函數(shù)包含12個，f1～f7是單峰函數(shù)，f8～f12是多峰函數(shù)，所有測試函數(shù)維度D設置為30。

分別選取同類算法FA，MFA，CFA，WSSFA，VSSFA，RaFA，ApFA算法與本文提出的MSFLFA進行實驗比較。所有算法的最大評估次數(shù)均設置MAX_FEs為5e5。MSFLFA的螢火蟲種群大小N為20，步長因子α為0.2，初始吸引度β0為1.0，光吸收系數(shù)γ為1.0，深度學習次數(shù)Deep_count為2 400(該參數(shù)設置依據(jù)參考文獻[29])，其它比較算法的參數(shù)設置及仿真實驗數(shù)據(jù)均來自參考文獻[30]。

3.2 實驗結果分析

表1給出了8種算法在基準測試函數(shù)上獨立運行20次的平均適應值。分析結果可知，MSFLFA在12個測試函數(shù)上的精度均優(yōu)于FA，VSSFA和WSSFA算法；對于單峰函數(shù)，MSFLFA在5個測試函數(shù)上的精度優(yōu)于所有算法；對于多峰函數(shù)，MSFLFA在4個測試函數(shù)上的精度優(yōu)于所有算法；對于單峰函數(shù)f6，MSFLFA，MFA，CFA，RaFA和ApFA均搜索到了理論最優(yōu)值；MSFLFA在多峰函數(shù)f12上的精度顯著優(yōu)于所有算法。

表2分類給出了8種算法在各類型函數(shù)集上的Friedman值，結果表明MSFLFA在所有函數(shù)、單峰和多峰函數(shù)的綜合性能均擁有最小Friedman值，算法綜合優(yōu)化性能最佳。

表1 8種算法在基準試函數(shù)上的平均適應值

表2 8種算法在基準測試函數(shù)上的Friedman值

為了更加直觀觀察算法的收斂速度和求解精度，圖1顯示了基準測試函數(shù)的所有算法收斂曲線圖。從圖1(a)，(b)，(c)，(d)，(f)和(i)中可以看出，MSFLFA收斂曲線在整個尋優(yōu)進程中一直呈線性遞減的趨勢，表明該算法不斷地搜索全局最優(yōu)值，未出現(xiàn)陷入局部最優(yōu)的情況。雖然圖1(e)，(h)，(j)和(l)中的MSFLFA收斂曲線在搜索進程后期趨于平緩，但其都處在所有算法的收斂曲線下方，表明MSFLFA搜索精度均優(yōu)于其它算法。圖1(f)和(i)中的MSFLFA收斂曲線在尋優(yōu)進程前期已收斂完成，表明該算法提前完成尋優(yōu)任務并已搜索到全局最優(yōu)值。

圖1 8種算法在基準測試函數(shù)的收斂曲線圖

4 MSFLFA在年徑流預測中的應用

4.1 MSFLFA與SVR預測模型融合

選擇以實測值和預測值的均方誤差(Mean Squared Error，MSE)作為MSFLFA的目標函數(shù)，其定義如下：

(8)

確定MSFLFA目標函數(shù)之后，對{C，ε，σ}參數(shù)進行尋優(yōu)的描述如下：

(1)確定模型的訓練和測試樣本，設置{C，ε，σ}參數(shù)的搜索范圍；

(2)初始化種群數(shù)量N、步長因子α、初始吸引度β0、光吸收系數(shù)γ；

(3)設定搜索范圍，初始化位置，將{C，ε，σ}參數(shù)與螢火蟲的位置關聯(lián)，即xi=(xi1，xi2，xi3)=(Ci1，εi2，σi3)；

(4)執(zhí)行最優(yōu)和非最優(yōu)螢火蟲學習策略設置的進化流程。若MSFLFA達到終止條件，則跳轉(zhuǎn)至步驟(5)，否則，重復執(zhí)行步驟(4)；

(5)輸出最優(yōu)螢火蟲位置信息，即SVR模型的{C，ε，σ}參數(shù)值，SVR利用參數(shù)值進行模型預測。

4.2 MSFLFA-SVR預測模型應用

實驗數(shù)據(jù)引自參數(shù)文獻[23]，即以陜西府谷縣黃甫川水文站1979—2008年的實測徑流值為預測對象，同時以年降水量和蒸發(fā)量為年預測因子。

實測徑流樣本數(shù)量共30 a，將前25 a實測徑流值作為MSFLFA-SVR模型的訓練樣本，后5 a實測徑流值作為MSFLFA-SVR模型的預測樣本，同時以實測徑流值和預測徑流值的MSE作為MSFLFA-SVR的評價函數(shù)。設螢火蟲種群大小N為20，迭代次數(shù)ITERs為5000，步長因子α為0.2，初始吸引度β0為1.0，光吸收系數(shù)γ為1.0，深度學習次數(shù)Deep_count為24，懲罰系數(shù)C、不敏感損失系數(shù)ε和核參數(shù)σ的取值區(qū)間分別為[0.01，500]、[0，1]和[0.1，100]。

4.3 MSFLFA-SVR預測模型性能評價

通過MSFLFA對SVR模型進行參數(shù)優(yōu)化，得到一組最優(yōu)的參數(shù)組合{C，ε，σ}={37.87，0.02，4.14}，使得MSFLFA-SVR模型預測效果最佳。獲得最佳模型之后，將30 a的預報因子輸入到模型中，得到最佳的年預測徑流值。同時，將年預測徑流值與FA-SVR、WSSFA-SVR、MFA-SVR、RaFA-SVR及文獻[24]中BP-ANN-SVR和PPR-SVR模型結果進行比較分析，預測結果見圖2。

圖2 1979—2008年各算法模型年徑流值訓練和預測結果

為對MSFLFA-SVR模型預測性能進行評價，選取平均絕對誤差R1、平均相對誤差R2%和確定性系數(shù)DC作為模型評價指標[31]。經(jīng)計算，獲得上述7種模型性能評價指標值，計算結果見表3。

表3 7種預報模型的擬合和預報值

分析表3，通過比較3個指標，對于訓練樣本，F(xiàn)A-SVR模型的3個指標均優(yōu)于其它6種模型，說明該模型對數(shù)據(jù)擬合結果較好，MFA-SVR、WSSFA-SVR、RaFA-SVR和MSFLFA-SVR模型的3個指標值近似相等，與FA-SVR模型相比較，擬合效果差于FA-SVR模型，PPR-SVR和BP-ANN-SVR模型的3個指標基本一致，相較于其它5種訓練結果，該2種模型擬合效果較差。對于測試樣本，F(xiàn)A-SVR模型的3個指標發(fā)生了較大變動，R1和R2%分別增大了0.44和0.10，DC減小了0.80，表明該模型預測性能極差，其它6種模型的3個評價指標波動范圍相似，相較于訓練擬合結果，預測性能均有所下降。

為綜合評估7種模型的性能，對上述評價指標進行秩均值檢驗(訓練和測試樣本各3個)，秩均值越小，表明模型綜合預報性能好，測試結果見表4。通過秩均值排序。由表4數(shù)據(jù)可知，MSFLFA-SVR模型的秩均值最小，其綜合預報性能均優(yōu)于其它6種模型。

表4 7種模型評價指標的秩均值

5 結束語

MSFLFA通過為最優(yōu)螢火蟲和非最優(yōu)螢火蟲分配不同的學習策略，實現(xiàn)了種群的全局探索和局部開采，且在測試函數(shù)上尋優(yōu)性能表現(xiàn)優(yōu)異。為驗證MSFLFA-SVR模型的性能，將其應用于陜西府谷縣黃甫川年徑流量的預測，與其它6種預測模型對比，MSFLFA-SVR模型取得了更好的預測效果，為徑流預報提供了一種有效的方法。MSFLFA僅從更新策略改善算法性能，未對算法固有參數(shù)做優(yōu)化，同時，算法應用方面僅進行了年徑流預報，預報粒度較大，無法滿足不同預報需求。因此，算法固有參數(shù)設置、月徑流和日徑流預報是下一步研究工作重點。