基于GNFA-SVR污水出水總氮軟測量研究

郭利進(jìn)，李博侖

（天津工業(yè)大學(xué)控制科學(xué)與工程學(xué)院，天津 300387）

氮是微生物繁殖的必要元素，但也是引發(fā)水體富營養(yǎng)化的元素之一〔1〕。出水總氮（TN）是污水排放的重要指標(biāo)之一，當(dāng)出水TN檢測結(jié)果不達(dá)標(biāo)時(shí)，污水廠會(huì)采取應(yīng)急措施，通過調(diào)整某個(gè)工藝段的工藝參數(shù)來解決問題。因此，根據(jù)實(shí)時(shí)進(jìn)水水質(zhì)情況判斷未來出水TN是否存在超標(biāo)可能并提前對工藝參數(shù)進(jìn)行調(diào)整，對污水處理系統(tǒng)的優(yōu)化具有重要意義。

在污水水質(zhì)預(yù)測創(chuàng)建模型方面，眾多學(xué)者從污水處理過程機(jī)理和數(shù)據(jù)驅(qū)動(dòng)方式去構(gòu)建模型。O.S.DJANDJA等〔2〕提出一種根據(jù)生化反應(yīng)溫度與原料中元素含量關(guān)系的前饋神經(jīng)網(wǎng)絡(luò)氮含量預(yù)測模型。喬俊飛等〔3〕采用一種基于遞歸徑向基（RBF）神經(jīng)網(wǎng)絡(luò)算法的氨氮軟測量方法。盡管神經(jīng)網(wǎng)絡(luò)算法在預(yù)測領(lǐng)域有著不錯(cuò)的性能表現(xiàn)，但也存在一定局限性，如收斂速度慢、易陷入局部極值、過度擬合等〔4〕。

與神經(jīng)網(wǎng)絡(luò)算法相比，支持向量回歸機(jī)（SVR）具有訓(xùn)練時(shí)間短、泛化能力強(qiáng)、預(yù)測精度高的優(yōu)點(diǎn)，并且克服了神經(jīng)網(wǎng)絡(luò)算法在小樣本、非線性問題中的重復(fù)性差、過度擬合等不足，更適用于預(yù)測研究〔5-8〕。然而，SVR模型構(gòu)建的關(guān)鍵是解決核參數(shù)σ和懲罰因子c的最優(yōu)取值問題，初期算法的參數(shù)選取多依賴于主觀經(jīng)驗(yàn)，缺乏嚴(yán)謹(jǐn)?shù)臄?shù)學(xué)依據(jù)。本研究通過引入基于高斯函數(shù)遞減步長的螢火蟲算法（GNFA）對SVR參數(shù)進(jìn)行優(yōu)化，可有效提高模型的預(yù)測準(zhǔn)確度。

1 基本原理

1.1 支持向量回歸機(jī)

支持向量回歸機(jī)（SVR）是一種監(jiān)督性機(jī)器學(xué)習(xí)算法，是在支持向量機(jī)（SVM）理論基礎(chǔ)上延伸出的一種對回歸分析問題的應(yīng)用模型。SVR通過引入不敏感損失系數(shù)ε讓所有的樣本點(diǎn)與超平面距離和最小，得到一個(gè)最佳條狀區(qū)域（2ε寬度），以達(dá)到期望風(fēng)險(xiǎn)最小的目的〔9〕。

給定n個(gè)訓(xùn) 練 樣 本D={(xi，yi)|i=1，2，…，n}，且xi∈Rn，yi∈Rn，xi為預(yù)測模型的相關(guān)輸入向量，yi為對應(yīng)的輸出向量，構(gòu)造回歸決策函數(shù)，見式（1）。其中，ω是權(quán)向量；φ(x)是原始樣本數(shù)據(jù)的非線性映射函數(shù)，可將數(shù)據(jù)映射至特征空間；b為偏置值。

引入不敏感損失系數(shù)ε與懲罰因子c，yi為實(shí)際值，f(xi)為模型的輸出預(yù)測值，當(dāng)yi與f(xi)差的絕對值大于ε時(shí)計(jì)算損失；反之則認(rèn)為預(yù)測正確，不計(jì)算損失〔10〕。不敏感損失系數(shù)ε對模型影響較小，選取0.15即可。根據(jù)最小結(jié)構(gòu)風(fēng)險(xiǎn)原則，建立約束條件〔式（2）〕與目標(biāo)函數(shù)〔式（3）〕。

圖1為SVR的示意圖，黑色實(shí)心的數(shù)據(jù)點(diǎn)表示模型預(yù)測結(jié)果是有效的。

圖1 支持向量回歸示意Fig.1 Support vector regression diagram

懲罰因子c的取值影響模型的結(jié)構(gòu)風(fēng)險(xiǎn)，取值過大，結(jié)構(gòu)風(fēng)險(xiǎn)越大，模型容易出現(xiàn)過擬合；取值過小，模型會(huì)過于簡單化，易出現(xiàn)欠擬合問題〔11〕。

引入松弛變量ξi與ξ*i，建立新的目標(biāo)函數(shù)〔式（4）〕。

引入拉格朗日乘子α、α*，將式（4）變換為對偶形式〔式（5）〕。式中，K(xi，xj)=φ(xi)φ(xj)為核函數(shù)。

根據(jù)式（5）的函數(shù)關(guān)系式可求出α與α*，最終可得到回歸函數(shù)，表達(dá)式見式（6）。

映射時(shí)發(fā)現(xiàn)，當(dāng)變量增多時(shí)，映射到高維空間的維度呈指數(shù)增長，導(dǎo)致計(jì)算難度驟增；而引入核函數(shù)可使低維計(jì)算結(jié)果等效于高維。核函數(shù)根據(jù)需要選取，研究中選用徑向核函數(shù)，表達(dá)式見式（7）。σ為徑向核函數(shù)的核參數(shù)，其取值大小也會(huì)影響模型效果，取值過小會(huì)導(dǎo)致模型的泛化能力減弱，取值過大會(huì)令模型出現(xiàn)過擬合現(xiàn)象。

綜上，SVR模型構(gòu)建的關(guān)鍵是解決核參數(shù)σ和懲罰因子c的最優(yōu)取值問題。

1.2 改進(jìn)螢火蟲算法

螢火蟲算法（FA）是一種啟發(fā)式算法，是根據(jù)自然界中螢火蟲的發(fā)光行為提出的〔12〕。在FA中，空間里的每一個(gè)解就好比螢火蟲種群里的每一只螢火蟲，空間里的初始解可以理解為螢火蟲種群的初始位置，螢火蟲通過個(gè)體之間的吸引進(jìn)行移動(dòng)，完成位置的更新，即完成解的更新〔13〕。

搜索過程涉及到螢火蟲的發(fā)光亮度和相互吸引度〔14〕，這2個(gè)參數(shù)都與螢火蟲之間的距離成反比，這與自然現(xiàn)象中光在空間傳播時(shí)被傳播介質(zhì)吸收而逐漸衰減的特性一致。

FA算法中螢火蟲發(fā)光強(qiáng)度的計(jì)算公式見式（8）。

式中，Ir表示發(fā)光強(qiáng)度；rij表示螢火蟲i與螢火蟲j之間的距離，rij=||Xi-Xj||；吸收系數(shù)γ的改變可以理解為傳播介質(zhì)的改變，可設(shè)置為常數(shù)；當(dāng)螢火蟲個(gè)體與光源的距離逐漸趨近于0時(shí)，螢火蟲的發(fā)光亮度最為明亮，即為I0〔15〕。

螢火蟲的吸引度β見式（9），與螢火蟲發(fā)光強(qiáng)度同理，當(dāng)螢火蟲個(gè)體與光源的距離逐漸趨近于0時(shí)，吸引度最大，即為β0〔16〕。

螢火蟲i被吸引向螢火蟲j移動(dòng)的位置更新由式（10）決定。

式中，X′i為下一次迭代的位置；α(rand-0.5)為擾動(dòng)項(xiàng)；α為步長因子，是介于0～1之間的常數(shù)；rand是［0，1］上服從均勻分布的隨機(jī)因子。擾動(dòng)項(xiàng)的加入是為了擴(kuò)大尋優(yōu)區(qū)域〔17〕。

在文獻(xiàn)〔18〕中，優(yōu)化粒子群算法引入了線性遞減權(quán)重。受該思路啟發(fā)，不妨將螢火蟲的步長進(jìn)行線性遞減優(yōu)化，定義該種改進(jìn)算法為線性遞減步長螢火蟲算法（LFA），具體步長變化表達(dá)式見式（11）。

因?yàn)镕A算法在實(shí)際問題的尋優(yōu)過程中是復(fù)雜非線性的，所以步長也呈動(dòng)態(tài)非線性變化預(yù)計(jì)會(huì)出現(xiàn)較好的尋優(yōu)結(jié)果。結(jié)合高斯函數(shù)分布的特性，猜想當(dāng)步長的隨機(jī)變化呈正態(tài)分布時(shí)，能夠提高算法的尋優(yōu)精度。因此，根據(jù)高斯函數(shù)分布特性，把相應(yīng)的步長參數(shù)帶入步長變化表達(dá)式，實(shí)現(xiàn)對步長α的非線性映射，得到最終的高斯遞減函數(shù)調(diào)整步長策略，見式（12）。

式中，常數(shù)k可以影響曲線的凹凸程度，即曲線的變化率可通過k進(jìn)行調(diào)整。

以步長取值范圍為0～1、最大迭代次數(shù)為100時(shí)為例，圖2給出了不同k值對應(yīng)步長的變化曲線。

圖2 k取不同值時(shí)步長的變化Fig.2 Change of step with different k

由圖2可看出，當(dāng)k≤0.1時(shí)，步長過早收斂于0，導(dǎo)致尋優(yōu)效果過于局限；當(dāng)k≥0.3時(shí)，步長在整個(gè)迭代過程中的變化趨勢較為平緩，不利于算法收斂。所以本研究中k取0.2，高斯函數(shù)遞減步長的螢火蟲算法（GNFA）位置更新表達(dá)式見式（13）。

2 GNFA算法性能測試

為驗(yàn)證改進(jìn)算法的性能，分別將FA、LFA與提出的GNFA進(jìn)行比較。選擇Sphere函數(shù)和Ackley函數(shù)作為測試算法性能的函數(shù)，Sphere函數(shù)多用來檢驗(yàn)算法的收斂速度和精度，Ackley函數(shù)可以測試算法的局部尋優(yōu)能力。這2個(gè)函數(shù)是最常見的測試函數(shù)，被廣泛應(yīng)用在算法測試性能實(shí)驗(yàn)中〔19〕。

在仿真環(huán)境中，螢火蟲步長改進(jìn)優(yōu)化算法的具體參數(shù)設(shè)置為：螢火蟲數(shù)目M=20，光強(qiáng)吸收系數(shù)γ=1.0，初始吸引度β0=0.8，步長因子αmax、αmin分別設(shè)置為0.8、0.1，最大迭代次數(shù)T=500，α的變化停止閾值Th設(shè)置為10-3。而FA算法步長設(shè)置為0.4，其他參數(shù)和改進(jìn)算法設(shè)置一致。測試函數(shù)的搜索范圍、期望目標(biāo)值、搜索空間維數(shù)及函數(shù)類型等見表1。

(2)仔豬出生后7～l5日齡開始飼喂優(yōu)質(zhì)全價(jià)飼料，換料應(yīng)逐步更換，并適量補(bǔ)充礦物質(zhì)，多喂富含維生素的青綠飼料，供給充足、清潔、新鮮飲水。

表1 測試函數(shù)特征Table 1 Test function characteristics

測試函數(shù)的優(yōu)化結(jié)果見表2，平均最優(yōu)值隨迭代次數(shù)增加的變化曲線見圖3。

表2 測試函數(shù)的優(yōu)化結(jié)果Table 2 The optimization results of the test functions

從表2的優(yōu)化結(jié)果可看出，相比LFA與標(biāo)準(zhǔn)FA算法，GNFA算法綜合表現(xiàn)更佳，說明當(dāng)采用高斯非線性遞減方式改進(jìn)步長策略時(shí)，螢火蟲個(gè)體尋優(yōu)性能強(qiáng)，且算法執(zhí)行時(shí)間也得到優(yōu)化。

從圖3可看出，3種算法對單峰值Sphere函數(shù)的優(yōu)化結(jié)果在尋優(yōu)精度上相差不大，但GNFA收斂速度快于其他2種算法。3種算法對多峰值A(chǔ)ckley函數(shù)的優(yōu)化結(jié)果說明，在迭代后期，標(biāo)準(zhǔn)FA算法定值步長容易導(dǎo)致算法陷入局部最優(yōu)；而LFA算法的性能表現(xiàn)介于其他2種算法。雖然線性遞減改進(jìn)步長策略并非最好的優(yōu)化方式，但也說明了步長α是影響算法性能的一個(gè)關(guān)鍵參數(shù)。總之，與傳統(tǒng)FA算法相比，GNFA算法在局部尋優(yōu)能力、全局搜索精度及收斂速度方面都得到改善。

圖3 Sphere函數(shù)（a）和Ackley函數(shù)（b）收斂曲線Fig.3 Convergence curves of Sphere function（a）and Ackley function（b）

3 基于GNFA-SVR的TN軟測量建模

3.1 模型數(shù)據(jù)選取與處理

如果進(jìn)水水質(zhì)組分全部被選作特征向量，其中部分特征變量存在弱相關(guān)性的問題往往會(huì)被忽視，導(dǎo)致維度過高，模型的計(jì)算時(shí)間和計(jì)算難度增加，不利于模型預(yù)測結(jié)果分析。主成分分析（PCA）是一種對高維度特征數(shù)據(jù)降維的統(tǒng)計(jì)分析法，有利于確定特征變量維數(shù)。圖4為使用PCA對數(shù)據(jù)特征變量進(jìn)行降維分析的相關(guān)度，當(dāng)選取6個(gè)特征變量時(shí)，總體的累計(jì)方差貢獻(xiàn)率已經(jīng)達(dá)到了90%。故本研究選取6個(gè)特征變量：進(jìn)水COD、TN、NH4+-N、溫度（T）、濁度（SS）、DO，并以出水總氮（TN）作為模型輸出。

圖4 累計(jì)方差貢獻(xiàn)率Fig.4 Cumulative variance contribution rate

從某污水廠2021年6月至9月水質(zhì)信息庫共獲取657組樣本數(shù)據(jù)，部分采集數(shù)據(jù)見表3。其中200組數(shù)據(jù)作為測試數(shù)據(jù)用于檢驗(yàn)?zāi)Ｐ托阅埽?57組數(shù)據(jù)作為訓(xùn)練集用于預(yù)測模型訓(xùn)練。

表3 部分采集數(shù)據(jù)Table 3 Part of collected data

樣本數(shù)據(jù)通常不直接使用，因?yàn)椴煌愋蛿?shù)據(jù)的量綱不同；且樣本數(shù)值波動(dòng)較大時(shí)，不利于模型的學(xué)習(xí)訓(xùn)練。為消除特征向量間的量綱影響、保證模型運(yùn)行速度以及提高預(yù)測精度，需對訓(xùn)練樣本數(shù)值進(jìn)行歸一化處理。歸一化公式見式（14）。

式中，x是任意一組樣本的某一維數(shù)據(jù)，最大值用max(x)表示，最小值用min(x)表示，x′為歸一化后的數(shù)據(jù)。

3.2 算法流程

步驟1：將經(jīng)過預(yù)處理后的n組水質(zhì)樣本數(shù)據(jù)(x1，y1)，(x2，y2)，…，(xn，yn)歸一化，并將其分為訓(xùn)練集和測試集2部分。

步驟2：對GNFA算法參數(shù)進(jìn)行初始化設(shè)置，包括螢火蟲數(shù)量M、初始吸引度β0、光強(qiáng)吸收系數(shù)γ、步長因子上下限αmax和αmin、最大迭代次數(shù)T。

步驟3：[c，σ]作為SVR的初始參數(shù)，也代表螢火蟲的1個(gè)位置；將訓(xùn)練集樣本作為SVR模型的輸入，并輸出預(yù)測集

步驟4：將預(yù)測值與實(shí)際值的均方根誤差作為GNFA-SVR算法的適應(yīng)度。

步驟5：判斷終止條件，即當(dāng)前迭代次數(shù)是否大于設(shè)置的最大迭代次數(shù)；滿足終止條件執(zhí)行步驟7，不滿足則執(zhí)行下一步驟。

步驟6：根據(jù)螢火蟲的相對亮度決定下一個(gè)位置的移動(dòng)方向并調(diào)整步長α；若目標(biāo)位置劣于之前的位置，螢火蟲位置保持不變，否則螢火蟲進(jìn)行位置更新，產(chǎn)生新的[c，σ]，返回步驟3。

步驟7：將輸出的最優(yōu)解［c*，σ*］代入SVR預(yù)測模型，利用測試集檢驗(yàn)算法的準(zhǔn)確性和可靠性，評估預(yù)測模型性能。

3.3 污水廠出水TN預(yù)測

設(shè)置模型參數(shù)：螢火蟲數(shù)目M=40，光強(qiáng)吸收系數(shù)γ=1.0，初始吸引度β0=0.8，步長因子αmax=0.8、αmin=0.1，最大迭代次數(shù)T=600，α變化停止閾值Th=10-3，懲罰因子c取值范圍為［0.01，1 000］，σ取值范圍為［0.001，10］。在MATLAB R2018b環(huán)境下借助支持向量回歸機(jī)工具箱運(yùn)行GNFA-SVR模型?；贕NFA-SVR模型對測試集的出水TN預(yù)測結(jié)果見圖5，GNFA-SVR可以較好地預(yù)測出水TN的實(shí)際值，預(yù)測誤差較小。

圖5 出水TN預(yù)測結(jié)果Fig.5 Prediction results of total nitrogen in effluent

為進(jìn)一步驗(yàn)證3種算法在實(shí)際問題上對SVR算法的優(yōu)化效果，現(xiàn)用GNFA-SVR、LFA-SVR和FA-SVR這3種預(yù)測算法對出水TN進(jìn)行預(yù)測。部分測試樣本的預(yù)測結(jié)果見圖6，其中GNFA-SVR與LFA-SVR的模型參數(shù)設(shè)置同上；FA-SVR所需參數(shù)步長α=0.4，其他參數(shù)和GNFA-SVR算法設(shè)置一致。

圖6 出水TN預(yù)測結(jié)果對比Fig.6 Comparison of predicted results of TN in effluent

從3種模型預(yù)測對比圖（圖6）可看出，3種模型的預(yù)測值與實(shí)測值的變化趨勢基本一致，其中GNFA-SVR模型預(yù)測值與實(shí)測值擬合效果最佳，預(yù)測誤差最小。

對3種模型的性能進(jìn)行具體評估，采用平均絕對誤差（MAE）、最大誤差（Maxer）和均方根誤差（RMSE）對預(yù)測結(jié)果進(jìn)行評定，表達(dá)式見式（15）～（17）。

式中，N為樣本總個(gè)數(shù)；為預(yù)測值；yi為實(shí)際值。

3種預(yù)測模型的評價(jià)指標(biāo)見表4。

表4 不同預(yù)測模型的評價(jià)指標(biāo)Table 4 Evaluation indicators of different prediction models

由表4可知，GNFA-SVR算法預(yù)測結(jié)果的各項(xiàng)指標(biāo)均優(yōu)于其他2種算法，其中MAE為0.352 1，預(yù)測誤差控制在±0.5 mg/L以內(nèi)可滿足污水廠實(shí)際生產(chǎn)需求；RMSE較LFA-SVR算法與FA-SVR算法分別降低了29.7%、63.2%。說明GNFA-SVR組合模型在改進(jìn)策略的優(yōu)化下預(yù)測性能有所提升，可實(shí)現(xiàn)對出水TN的精準(zhǔn)預(yù)測。

4 結(jié)論

針對污水處理廠出水水質(zhì)檢測存在大滯后性的問題，本研究提出一種基于GNFA-SVR算法的預(yù)測模型。在SVR算法基礎(chǔ)上，引入螢火蟲算法對其進(jìn)行參數(shù)尋優(yōu)。針對標(biāo)準(zhǔn)FA算法的步長不隨迭代次數(shù)變化、易出現(xiàn)陷入局部最優(yōu)的現(xiàn)象，提出基于高斯函數(shù)非線性遞減步長策略，有效提高了算法的性能。通過對污水處理廠的實(shí)際數(shù)據(jù)進(jìn)行仿真驗(yàn)證，GNFA-SVR算法在平均絕對誤差、最大誤差和均方根誤差3個(gè)評價(jià)指標(biāo)方面均明顯優(yōu)于LFA-SVR與FA-SVR算法，預(yù)測精度最高，誤差控制在±0.5 mg/L以內(nèi)，可實(shí)現(xiàn)對出水TN的精準(zhǔn)預(yù)測，在實(shí)際水質(zhì)預(yù)測方面具有良好的應(yīng)用前景。