基于EEMD和螢火蟲算法優化SVM的溶解氧預測

劉 晨，李 莎，叢孫麗，朱正偉

(1. 無錫太湖學院江蘇省物聯網應用技術重點建設實驗室，江蘇 無錫 214064；2. 常州大學信息科學與工程學院， 江蘇 常州213164)

1 引言

溶解氧(DO)濃度是漁業養殖水質安全的重要指標，也是影響養殖水產品品質的關鍵因素，對其進行實時預測具有重要的意義[1-2]。而影響水中溶解氧含量的各因素關系比較復雜，很難進行數學建模，水中溶解氧預測模型是一個復雜的非線性系統[3-4]。針對養殖池塘中溶解預測問題，國內外學者將研究都集中在溶解氧含量預測上，取得了很大的成就。

目前溶解氧預測方法主要有模糊方法、貝葉斯模型、人工神經網絡和支持向量機(SVM)等。神經網絡法雖然在非線性預測領域應用較廣，但是基于樣本無限大的假設和經驗風險最小化原則為基礎，存在易于局部最優值、過學習、不適于高維數、小樣本預測等缺陷嚴重影響了其應用效果[5]。支持向量機具有計算復雜度低、魯棒性強、收斂精度較高、非線性擬合能力較好等特點[6]。文獻[7]采用支持向量機作為水質預警方法，對漁業養殖池塘中水質參數進行預測。文獻[8]提出了基于支持向量機的花粉濃度預報模型。雖然支持向量機有較好地預測效果，但它對核函數參數的選取十分敏感，不合適的SVM核函數參數將會使得預測性能大大降低[9]。如何獲得最優的參數組合，目前比較常見的包括遺傳算法、粒子群算法、蛙跳算法等群智能算法[10]。這些算法雖然在不同程度上滿足SVM核函數參數選取的要求，但算法精度偏低、局部性能較差，很容易陷入局部極值，出現早熟現象，對于全局與局部性能的平衡，以及魯棒性、簡易性和計算精度的兼顧上有所欠缺。而螢火蟲算法以其極強的局部與全局優化性能、較高的魯棒性等天然優勢受到了學者的廣泛關注，目前以及成功應用于多模態函數優化、傳感器的噪聲測試、計算機模擬人工機器人等人工智能領域。本文利用螢火蟲算法進行 SVM 核函數參數優化，通過改進螢火蟲位置更新公式大幅提升收斂速度，并且在移動過程中引入亮度特征，避免算法后期出現的收斂不穩定現象，顯著地提升了算法的求解精度，從而確定了最佳的 SVM 核函數參數，提高SVM的預測性能。

由于溶解氧序列具有非線性和非平穩性的特點，如果上述方法直接在原始時序數據上建模，難以充分挖掘和利用各時頻特征信息，會影響預測性能。為了對溶解氧序列進行降噪處理，文獻[11]利用經驗模態分解(empirical mode decomposition，EMD)處理時序數據，EMD算法使用簡單方便，分解速度快，能較好地捕捉不同尺度下的特征信息。但經驗模態分解方法會帶來模態混疊的問題，集合經驗模態分解(ensemble empirical mode decomposition，EEMD)方法，利用高斯白噪聲的頻率均勻分布的統計特性彌補模態經驗分解的缺陷[12]。

為了提高預測精度和泛化性能，本文提出一個基于EEMD和螢火蟲算法優化SVM的組合模型來預測溶解氧。利用EEMD將DO時間序列分解成一組相對穩定的子序列，并利用相空間重構PSR重建分解子序列，在相空間中用螢火蟲優化SVM對各子序列進行建模預測。驗證結果表明本文提出的EEMD-FA-SVM組合預測模型能夠滿足現代漁業養殖水質精細化管理的高需求。

2 相關基本原理簡介

2.1 集合經驗模態分解

EEMD算法目的在于能夠根據信號自身的特點，自適應地將非線性、非平穩性的多模態信號分解為若干個平穩單一模態的固有模態函數(Intrinsic Mode Function，IMF)分量和一個余項[13]。傳統的EMD方法中IMF分量的不連續造成相鄰波形模態混疊現象。為了彌補這一缺陷，文獻[14]利用高斯白噪聲具有頻率均勻分布的統計特性，保證模態分解的準確性。分解出IMF需要滿足兩個條件；①在整個信號序列中，極值點的個數與過零點次數必須相等或者最多相差一點；②在任意時間點上，由信號局部極大值和極小值確定的上下包絡線的均值為零[15]。

對選定的DO時間序列x(t)進行經驗模態分解，分解步驟如下：

步驟1：添加一個白色噪聲序列，該序列服從(αε)2)的正態分布；

步驟2：將原始DO數據分解為一組IMF分量和一個余量Res。

步驟3：將步驟1重復n 次，每次添加一個不同振幅的白色噪聲序列。

步驟4：經過 n 次分解后，將不同IMF分量的平均值作為最終結果。

2.2 相空間重構

輸入向量的構造方法影響著SVM預測模型的預測精度和預測效率。同時，在利用SVM實現多步驟預測時，常常將預測值代替實際值。合理的輸入向量結構是減小多步預測累積誤差的有效措施， 可以提高多步預測的精度。

相空間重構是時間序列混沌分析和處理的關鍵步驟，能夠從時間序列中提取更多有用的信息，根據Takens嵌入定理建立可靠的數學基礎，可以從一維混沌時間序列中重構一個多維的時間序列矩陣[16]。混沌時間序列的分析和預測都在這個重構的相空間中進行的。文獻[11]通過對實測一維時間序列X(k)進行相空間重構，得到的多維數序列如式(1)和式(2)所示。

X(k)=[x(k)，x(k+τ)，…，x(k+(D-1)τ)]

(1)

Y(k)=x(k+D)

(2)

其中，D為嵌入維數，τ為延遲參數，X(k)是D維輸入向量，且t=(D-1)τ+1，(D-1)τ，…，N.

2.3 螢火蟲算法

螢火蟲算法(Firefly Algorithm，FA)是一種基于群智能的隨機搜索技術，其思想源于螢火蟲向比更亮的螢火中移動這一生物學特性[17]。在搜索空間，螢火蟲的位置表示優化問題的解，亮度對應優化問題的適應值。螢火蟲不斷向更高的螢火蟲移動，直至達到預設的算法終止條件，完成尋優任務。算法的數學描述如下：

設螢火蟲數量N，維度D，第i只螢火蟲的位置表示為xi=(xi1，xi2，…，xiD)，i=1，2，…，N，和第j只螢火蟲的位置表示為xj=(xj1，xj2，…，xjD)，j=1，2，…，N。螢火蟲i和螢火蟲j的距離rij計算如式(3)所示

(3)

其中，xid和xjd分別表示第i只和第j只螢火蟲的第d維位置。螢火蟲的亮度及吸引度的計算如式(4)和式(5)所示

(4)

(5)

其中，I0表示螢火蟲初始亮度，β0表示螢火蟲初始吸引度，γ為光強吸收系數，熒光會隨著距離的增加和傳播媒介的吸收而減弱。

在螢火蟲算法中，每只螢火蟲的自身亮度與其所在位置目標函數的適應度值有關。適應度值越大熒光素濃度越大，自身亮度越大，位置越優。吸引度與亮度有關，亮度越大，吸引力越強。而在移動過程中亮度和吸引度又都與個體之間的距離成反比。FA算法每次迭代都可以看成是由熒光素濃度更新階段和移動位置更新階段組成。

螢火蟲的熒光素濃度更新如式(6)所示。

li(t)=(1-ρ)li(t-1)+γJ(xi(t))

(6)

其中，li(t)表示第t代第i個螢火蟲的熒光素濃度值；ρ∈(0，1)為熒光素揮發因子，主要控制熒光素值；γ為熒光素更新率，J(xi(t))為目標適應函數，對熒光素值有決定性作用。

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

(7)

其中Xid(t)表示螢火蟲i第t代第d維位置，αi(t)表示螢火蟲i的第t代步長因子，ε服從均勻分布，取值范圍為[-0.5，0.5]。此外，動態決策域相當于螢火蟲的感知范圍，在此范圍內，螢火蟲數量過多就縮小決策域；反之擴大決策域，螢火蟲i的動態決策域半徑的更新如式(8)所示

(8)

其中，rs為感知半徑，ni是鄰居數的閾值，|Ni(t)|為鄰域范圍內螢火蟲的個數，β為動態決策更新率。

2.4 支持向量機

支持向量機核心思想就是尋求某種歸納原則實現結構風險最小化，本質上屬于一種最優化問題，將線性不可分的待測數據空間的樣本通過核函數投影到高維特征空間，使原本不可分的問題在特征空間上實現線性可分[18]。

對于訓練樣本集：(xi，yi)，i=1，2，…N，其中xi是輸入向量，yi是相關的輸出向量。由此產生的SVM優化問題可以通過以下形式表述

(9)

其中，C為正則化常數，ei表示所需輸出和實際輸出之間差異的松弛變量。w是特征空間中分類超平面的系數向量，b是分類面的閾值，決定了超平面的位置。φ(.)是輸入空間與輸出空間之間的映射函數。

(10)

從KKT(Karush-Kuhn-Tucher)條件可知，式(10)必須滿足?L/?w=0， ?L/?b=0，?L/?e=0，?L/?a=0。通過以矩陣形式表示的線性方程組(11)求解。

(11)

其中，a=[a1…al]，Y=[y1…yl]和I表示單位矩陣，定義ψkj=K(xk，xj)，k，j=1，2，…，l，滿足Mercer定理。SVM回歸模型如式(12)所示

(12)

3 預測模型的構建

如上所述，溶解氧時間序列具有非線性和非平穩的特征，導致時間序列預測中的不可克服的困難。為此，本文根據分解和集成原理，采用了EEMD-FA-SVM混合模型，基于EEMD和螢火蟲算法優化SVM的溶解氧預測流程圖如圖1所示。

圖1 基于EEMD-FA-SVM的溶解氧預測流程圖

3.1 數據獲取

以江蘇省溧陽市埭頭鎮黃家蕩特種水產養殖場為試驗場所，養殖場占地153hm2。養殖池塘中有一個池塘循環水系統，配備 pH 傳感器、溶解氧傳感器、水溫傳感器等，增氧泵、無線監控系統等現代化漁業裝置[20]。養殖環境數據來源于水產養殖遠程無線監控系統，每30 min采樣1次，采樣周期為 2018 年 11月1日～2018年11 月30日，共計1440個數據，取前1392個數據為訓練數據，后48個數據為測試數據，DO時間序列如圖2所示。

圖2 DO時間序列

(13)

3.2 數據分解

在使用EEMD將原始DO時間序列分解為9個IMF分量和一個余量Res時，首先確定 n 和α的值。根據文獻[13]， n 和α的值分別為100和0.25。溶解氧時間序列由EEMD分解后生成的 IMF 分量較多，分別為 IMF1、IMF2、IMF3、IMF4、IMF5、IMF6、IMF7、IMF8、IMF9和 Res。分解的結果如圖3所示。

圖3 EEMD分解溶解氧序列結果

3.3 螢火蟲算法優化SVM

將正則化參數C和核函數σ設為螢火蟲的位置坐標，螢火蟲亮度由支持向量機訓練返回誤差決定，然后利用FA算法找出使返回誤差最小的螢火蟲位置，其坐標所對應的C和σ即為尋找的最優參數[21]。以均方誤差函數(MSE)作為目標函數F，其表達式如式(14)所示。

(14)

利用FA算法優化RBF核函數參數流程如下

1)初始化每個螢火蟲的熒光素、動態決策域等參數，設置最大迭代次數為iter_max和螢火蟲初始位置，初始位置隨機分布在預先設定好的參數選擇范圍內。位置的x軸范圍是Cmin～Cmax，y軸范圍是σmin～σmax。

2)按照式(6)進行亮度的更新，計算出每個螢火蟲的適應度值J(xi(t))；

4)計算螢火蟲i向螢火蟲j移動的概率pij(t)；

5)朝pij(t)值最大的方向，根據式7)對位置進行更新，若更新后值優于原先位置，則替換；

6)按照式(8)更新動態決策半徑；

7)回到步驟2進行循環，若已經到達最大循環次數，則終止循環，輸出結果，程序結束。

3.4 多步預測

對每個分量的預測結果進行疊加，可獲得下一個時刻的DO預測值。因為SVM只能實現單步預測，即一次只能獲得一個預測值，運用單點迭代法實現多步預測，即預測未來24h的DO值，其原理是利用預測值代替真實值作為預測輸入預測向量的元素來獲取下一個的預測值。

3.5 誤差分析

為了測試本文提出的EEMD-FA-SVM模型的預測能力，選擇了三種不同的統計指標， 即誤差均方根 (RMSE)、平均相對誤差 (MAPE)、平均絕對誤差 (MAE)來評估預測精度。誤差均方根、平均相對誤差和平均絕對誤差如式(15)、式(16)和式(17)所示。

(15)

(16)

(17)

4 結果與討論

采用本文建模方法得到的DO各個分量的嵌入維數D、延遲時間τ、模型參數C和σ如表1所示。

為了評估EEMD-FA-SVM預測模型的性能，將FA-SVM、EEMD-FA-BP和EEMD-PSO-SVM模型的預測結果與 EEMD-FA-SVM 模型的預測結果進行比較。幾種模型采用了相同的訓練數據集并直接預測。EEMD-FA-SVM的預測結果及其誤差如圖4所示，四種模型的預測結果如圖5所示。

圖4 EEMD-FA-SVM的預測結果及其誤差

圖5 四種模型的預測結果圖

從圖4可以看出，本文所提模型EEMD-FA-SVM的預測結果接近真實結果，誤差曲線圍繞零點上下波動。其中，6：00～12：00時間段預測結果曲線與實測值曲線趨勢基本一致，波動較小。從圖5可以看出，與其它3種模型進行對比，EEMD-FA-SVM的預測結果曲線是最接近實測值曲線的。

為了更準確地評價模型的精度，本文還計算出四種預測模型中的均方根誤差(RMSE)[11，22]、平均相對誤差 (MAPE)、平均絕對誤差 (MAE)，預測未來6h、12h、18h和24h的預測誤差結果如圖6所示。

圖6 預測未來6h、12h、18h和24h誤差指標

從圖6可以看出，在相似的前提條件下，預測24h DO值時，本文所提模型EEMD-FA-SVM的均方根誤差(RMSE)比FA-SVM模型降低了16.81%，平均相對誤差(MAPE)降低了2.94%，平均絕對誤差(MAE)降低了15.87%，這表明EEMD-FA-SVM模型比FA-SVM模型有更好的預測效果，即集合經驗分解對預測模型建立的重要性。此外，本文所提模型的均方根誤差比EEMD-FA-BP模型降低了8.42%，平均相對誤差降低了2.38%，平均絕對誤差降低了10.36%，可以看出EEMD-FA-SVM模型比EEMD-FA-BP模型的預測精度高。

EEMD-FA-SVM模型的3種誤差指標比EEMD-PSO-SVM分別降低了7.69%，1.35%，9.07%，從中可以發現螢火蟲算法優化支持向量機比粒子群算法優化支持向量機具有更好的預測結果。

本文方法在預測不同時長溶解氧的均方根誤差、平均相對誤差、平均絕對誤差隨預測時長的增加，大體呈現遞增趨勢。即在不同的預測時間下，預測精度有明顯變化，預測時長越短，預測值和實測值吻合度越好，預測精度越高；當預測時長為24h時，預測結果和實測結果偏差較大，這是由于多步預測采用的單點迭代法，起初預測值再輸入向量中占比不大，越往后預測值占比就越來越多，導致預測精度下降，其它三種模型也存在著相同的問題。但是，本文提出的模型在預測12h時的平均相對誤差比預測6h增加了0.1%，預測18h時比預測6h 增加0.17%，預測24h時比預測6h增加0.26%，而FA-SVM在預測12、18和24h 時的平均相對誤差比預測6h分別增加 0.35%，0.84%和0.92%。EEMD-FA-BP在預測12、18和24h 時的平均相對誤差比預測6h分別增加0.29%，0.92%和1.28%，EEMD-PSO-SVM在在預測12、18和24h 時的平均相對誤差比預測6h分別增加0.27%，0.36%，0.59%。在長時間的預測中，本文提出的模型誤差遞增速度比其它三種模型慢，所以本文所提出的模型更具有穩定性，比其它三種模型更加適合長時間的預測。此外，預測未來6h、12h、18h和24h溶解氧的均方根誤差分別為 0.1444，0.1301，0.1745 和 0.1860，預測時長在12h的預測值與實際值之間的偏差小于6h，可能由于本時段其它偶然因素，比如水的溫度，水生生物的活動情況等。也可能由于建立的模型在該時刻的預測值與實際值偏差過大所致。

各項誤差指標證明了EEMD-FA-SVM的組合預測模型的可行性，預測精度要高于常用預測模型。EEMD-FA-SVM模型可以對DO值進行有效預測，在溶解氧預測系統中應用顯著，為水產養殖調控提供科學決策。

5 結論

本文提出了一種“分解-預測-重建”的組合模型，該模型結合了集成經驗模式分解和螢火蟲算法優化的支持向量機。先利用EEMD對DO時間序列進行分解，EEMD可以有效地降低溶解氧時序的非平穩特性，再對分解出來的各個分量進行空間重構，構造輸入輸出向量，利用螢火蟲算法對SVM的超參數進行優化后，再將優化后的SVM用于溶解氧的預測。與傳統方法相比，螢火蟲算法優化后的SVM網絡收斂速度快，預測精度高，泛化能力強。在今后的研究中，將改進螢火算法對支持向量機參數進行優化，進一步提高預測精度和準確性。