基于KPCA-PSO-ELM算法的地表水化學需氧量紫外-可見吸收光譜檢測研究

鄭培超, 周椿棪, 王金梅, 尹義同, 張 莉, 呂 強, 曾金銳, 何雨欣

重慶郵電大學光電工程學院, 光電信息感測與傳輸技術重慶重點實驗室, 重慶 400065

引 言

隨著社會和工業(yè)發(fā)展, 大量生產(chǎn)生活廢水直接排入地表水域中, 造成嚴重污染。 由水污染帶來的問題已經(jīng)嚴重影響社會發(fā)展。 化學需氧量(chemical oxygen demand, COD)是水質檢測的重要指標, 用于表征水體有機物含量。 化學需氧量越大, 說明水體污染越嚴重, 對生物危害越大。 目前, 化學需氧量檢測方法主要是重鉻酸鉀耗氧量(CODcr)和高錳酸鉀耗氧量(CODMn)等實驗室化學方法。 傳統(tǒng)化學方法存在分析時間長、 二次污染等缺點, 難以滿足實際場景下的檢測要求。 為了防止水體污染繼續(xù)加重, 亟需發(fā)展一種水質COD的快速、 實時檢測方法。

近年來, 由于光譜法快速、 無污染等優(yōu)點, 被廣泛應用到水質檢測中。 針對水質污染物常用的光譜檢測技術有熒光光譜、 近紅外光譜及紫外吸收光譜等[1-4]。 其中, 紫外-可見吸收光譜(UV-Vis)作為一種高效、 實時、 精確的光譜方法。 通過測量水體吸光度實現(xiàn)水質COD的定量分析, 是當前水質檢測中運用最為廣泛的光譜技術之一[5-7]。

目前, 使用紫外吸收光譜進行COD的檢測多采用基于單波長和多波長的方法, 采用COD在其某些特定波長下的吸收峰實現(xiàn)其定量分析。 Mai等[8]采用吸收光譜在254 nm處的吸光度結合偏最小二乘回歸和主成分回歸對印染廢水COD進行檢測, 將結果相對誤差控制在5%以內。 方坷昊等[9]采用COD在240～300 nm波段內的特征波長和560 nm處的特征波長建立非線性預測模型, 標準溶液預測結果誤差控制在3.5%內。 Zhang等[10]采用400和600 nm處斜率對紫外-吸收光譜的一階導數(shù)光譜進行補償, 將PLS預測模型R-squared(R2)提升至0.99。 然而, 實際水體組分復雜, 存在濁度、 色度等干擾, 單波長或多波長模型難以滿足預測的需求[11-13]。

基于全光譜的化學需氧量檢測方法, 逐漸成為當前研究熱點。 張崢等[14]使用主成分分析對全光譜數(shù)據(jù)進行壓縮, 并聯(lián)合粒子群優(yōu)化極限學習機(particle swarm optimization extreme learning machine, PSO-ELM)算法對COD進行檢測, 較ELM模型運行效率提升了一個量級。 Li等[15]使用改進的bagging算法對全光譜數(shù)據(jù)建模, 測試集R2達到了0.931 7, RMSEP降低到5.39 mg·L-1。 對于實際水體的檢測, 全光譜預測模型較單波長和多波長預測模型具備良好的抗干擾能力和較高的預測精度。 基于全光譜數(shù)據(jù)模型存在信息冗余、 高特征維度、 干擾復雜等缺點, 會導致模型無法收斂或過擬合。 對光譜信息進行預處理和數(shù)據(jù)壓縮方法直接影響全光譜吸收光譜法對COD的預測精度。

由于實際水體成分復雜, 紫外-吸收光譜存在非線性變化。 內核主成分分析作為一種非線性數(shù)據(jù)壓縮方法, 可以有效處理線性不可分的數(shù)據(jù)集。 本工作提出了一種基于內核主成分分析的全光譜水體化學需氧量檢測方法。 使用內核主成分分析對全光譜數(shù)據(jù)進行壓縮, 建立基于粒子群優(yōu)化的極限學習機回歸模型, 以實現(xiàn)對COD快速、 實時檢測。

1 實驗部分

1.1 裝置

實驗裝置如圖1所示, 主要包括光源、 光譜儀、 樣品槽、 光纖等。 光源采用氘鹵燈(愛萬提斯, Avalight-Hal-Cal-Mini), 光源的出射光經(jīng)光纖耦合到10 mm×10 mm×30 mm石英比色皿。 產(chǎn)生的透射光經(jīng)光纖光譜儀(復享光學, NOVA)接收。 采集到的光譜數(shù)據(jù)傳輸?shù)接嬎銠C中進行數(shù)據(jù)處理與建模。

圖1 實驗裝置示意圖

1.2 樣品

采集了重慶主城內長江、 嘉陵江、 盤龍溪、 白溪、 白云水庫、 豐收水庫所采集到的6類地表水共217個樣本。 所有樣本經(jīng)實驗室化學法(CODMn)測定后, 作為光譜法建模結果的參考。 采用紫外-可見吸收光譜對水體有機物定量分析的理論基礎是朗伯-比爾定律, 不同的有機物有不同吸收峰, 不同濃度的污染物吸收強度也不同。 圖2是以去離子水為參比, 待測樣品所采集到的紫外-可見吸收部分光譜, 可以看出吸收段主要集中在400 nm之前, 吸收峰在200 nm附近。

圖2 地表水紫外-可見吸收光譜

1.3 數(shù)據(jù)處理

光譜數(shù)據(jù)處理流程如圖3所示。 將數(shù)據(jù)按照7∶3的比例隨機劃分成訓練集和測試集, 然后對訓練集和測試集的光譜數(shù)據(jù)進行去噪、 增強、 特征提取等處理, 對訓練集光譜數(shù)據(jù)進行回歸模型訓練, 得到光譜-COD濃度的反演模型; 再將測試集的光譜數(shù)據(jù)輸入訓練集的反演模型, 進行光譜數(shù)據(jù)與COD濃度反演, 獲得測試集的COD濃度。

圖3 光譜數(shù)據(jù)處理流程

1.3.1 光譜數(shù)據(jù)預處理

光譜數(shù)據(jù)中既有信息, 同時又包含了噪聲。 采用Savitzky-Golay(SG)濾波算法對光譜進行濾波處理, 可以有效抑制隨機噪聲, 提高信噪比。 原始輸入數(shù)據(jù)用給定的多項式的擬合值代替, 對每個數(shù)據(jù)反復此操作, 直到得到所有數(shù)據(jù)的平滑值[16]。 采用窗口大小為21、 多項式階數(shù)為3的SG濾波, 濾波算法也可以看作是一種加權平均的過程。 為了降低由于信號波動帶來的影響, 再將濾波后每個波長下的吸光度用區(qū)間積分值代替, 即每一個點的值都用該點前后10個點的積分值代替。 不同樣本光譜數(shù)據(jù)積分后, 會產(chǎn)生數(shù)據(jù)量綱帶來誤差, 還需要將所有光譜數(shù)據(jù)進行歸一化。 經(jīng)過濾波、 增強和歸一化等預處理后的數(shù)據(jù)較未處理的光譜更為平滑, 此時不同樣本在此吸收段內的區(qū)別更明顯, 如圖4所示。

圖4 預處理后紫外-可見吸收光譜圖

1.3.2 核主成分分析算法

內核主成分分析(kernel principal component analysis, KPCA)作為主成分分析(principal component analysis, PCA)改進算法, 是一種針對非線性數(shù)據(jù)的特征壓縮方法。 KPCA借助映射函數(shù)將數(shù)據(jù)從原始空間映射到高維特征空間H中, 使得數(shù)據(jù)在高維空間中可分[17],H中映射函數(shù)數(shù)據(jù)的協(xié)方差矩陣表示為式(1)

(1)

計算協(xié)方差矩陣CH的特征值和特征向量, 其特征值為λ, 特征矢量β。 此時定義矩陣M=[Mi,j]n×n,Mi,j=(φi)·(φj), 可通過核函數(shù)來確定。βk是β的第k個特征矢量, 對其進行歸一化處理, 即βkβk=1則可得到原始空間中任意一樣本x的映射數(shù)據(jù)φ(x)在特征矢量βk上的投影為式(2)

(2)

將映射數(shù)據(jù)中心化得到數(shù)據(jù)第k維的非線性主成分為式(3)

(3)

1.3.3 粒子群優(yōu)化算法

粒子群優(yōu)化算法模擬鳥群覓食行為, 每個粒子都代表種群中的個體。 每個粒子初始位置都是隨機, 通過不斷迭代從潛在解中找到當前最優(yōu)解[18-19]。 計算每次迭代后適應度函數(shù)以更新個體和群體極值。 粒子群優(yōu)化算法中粒子速度和位置的更新公式如式(4)和式(5)

Vid(t+1)=wvid(t)+c1r1[Pid-Xid(t)]+

c2r2[gid-Xid(t)]

(4)

χid(t+1)=χid(t)+vid(t+1)

(5)

式(4)和式(5)中,w為慣性權重,c1和c2分別為個體和群體的學習因子,r1和r2為[0, 1]之間隨機數(shù)。 通過不斷迭代確定粒子的最優(yōu)解。 權重w設置為0.6, 學習因子c1,c2設置為2。

1.3.4 極限學習機算法

假定存在N個不同樣本(xi,yi), 其中特征xi∈Rn, 標簽yi∈Rm, 具有K個隱含層節(jié)點前饋神經(jīng)網(wǎng)絡輸出可以表示為

(6)

式(6)中,G(x)為網(wǎng)絡的激勵函數(shù),αi是第i個連接輸入層到隱含層節(jié)點的權值,bi則是第i隱含層節(jié)點的偏差,βi是第i個連接隱含層節(jié)點到輸出的權值。 若存在K個隱含層節(jié)點的前饋神經(jīng)網(wǎng)絡能以零誤差逼近全部個樣本, 則存在αi,bi,βi使得式(7)和式(8)成立

(7)

Wβ=Y

(8)

(9)

式(9)中,W+為隱含層輸出矩陣W的Moore-penrose廣義逆。

1.3.5 KPCA-PSO-ELM回歸模型

ELM的初始輸入權值和隱含層偏差是隨機給定的, 而輸出權值矩陣則是由輸入權值矩陣和隱含層偏差計算得到, 通常會導致部分隱含層節(jié)點失效; 在實際應用中ELM可能需要設置一定數(shù)量的神經(jīng)元才能達到理想的精度。 因此采用粒子群優(yōu)化算法對ELM的輸入權值矩陣和隱含層偏差進行優(yōu)化, 可以得到最優(yōu)的神經(jīng)網(wǎng)絡。 對于模型中隱含層神經(jīng)元, 一般問題設置為20～40, 此處設為30, 將RMSE作為適應度函數(shù), 即模型的評價指標。 使用KPCA對預處理過的光譜數(shù)據(jù)進行壓縮, 再將數(shù)據(jù)輸入到ELM模型中, 最后使用PSO對ELM模型進行優(yōu)化, KPCA-PSO-ELM流程圖如圖5所示。

圖5 KPCA-PSO-ELM流程圖

2 結果與討論

圖6給出了利用KPCA-PSO-ELM模型的參數(shù)優(yōu)化結果。 由圖6(a)可以看出, 采用KPCA對光譜數(shù)據(jù)進行壓縮, 方差貢獻率隨著主成分數(shù)量增加快速上升, 在第5個核主成分時候方差累計貢獻率到達0.999, 內核主成分分析對光譜數(shù)據(jù)有明顯的降維效果。 圖6(b)顯示粒子群優(yōu)化算法經(jīng)過前100次迭代后適應度快速下降到一個穩(wěn)定值, 200次迭代后不再出現(xiàn)明顯下降, 經(jīng)過500次迭代后模型訓練集的RMSE下降至0.363 0 mg·L-1。

圖6 KPCA-PSO-ELM模型參數(shù)優(yōu)化

以實驗數(shù)據(jù)建立的KPCA-PSO-ELM模型對地表水COD預測結果, 如圖7所示。 可以看出測試集樣本包含在訓練集中, 訓練集擬合優(yōu)度R2為0.930 2[見圖7(a)], 測試集擬合優(yōu)度R2為0.932 0[見圖7(b)], 訓練集和測試集R2均在0.93以上, 僅出現(xiàn)了少數(shù)幾個離群點, 說明該模型性能滿足實際檢測要求, 泛化性較好, 針對地表水有好的預測結果。

圖7 KCA-PSO-ELM模型回歸結果

此外, 本文還對比了PCA、 套索回歸(least absolute shrinkage and selection operator, LASSO)、 連續(xù)投影算法(successive projection algorithm, SPA)光譜數(shù)據(jù)特征方法在PSO-ELM模型上的表現(xiàn)。 表1為數(shù)據(jù)經(jīng)過PCA壓縮后的結果, 可以看出第一主成分占據(jù)絕大部貢獻率, 前5個主成分累計貢獻率已經(jīng)達到99.94%。 選取方差累計貢獻率99.9%以上主成分建模。 Lasso返回非零權重的特征, 由于核參數(shù)lambda對結果影響極大, 對該參數(shù)進行了優(yōu)化, 結果如圖8所示。 SPA保留交叉驗證下10個投影最大的特征。

表1 PCA主成分貢獻率

圖8 Lambda交叉驗證曲線

以上四種光譜特征方法在PSO-ELM模型測試集RMSE分別為0.400 7、 0.715 1、 0.473 7和0.412 6 mg·L-1。 由圖9可以看出, 基于KPCA特征提取模型的相對誤差相比于其他三種模型, 誤差波動最小。

表2列出了測試集樣本的預測值和參考值以及相對誤差。 以上四種特征提取方法模型的相對誤差絕對值的均值為9.26%、 17.18%、 11.73%、 11.74%, 其中KPCA僅有一個樣本誤差超過40%, 證明在此模型中KPCA優(yōu)于其他三種方法。

表2 測試集樣本COD濃度ELM回歸預測值

表3為不同特征方法模型的參數(shù)比較, KPCA-PSO-ELM模型訓練集的RMSE為0.363 0 mg·L-1, 擬合優(yōu)度R2為0.930 2, 測試集的RMSE為0.400 7 mg·L-1, 擬合優(yōu)度R2為0.931 9。 KPCA-PSO-ELM模型的訓練集和測試集的R2均在0.93以上, 說明所提出的模型優(yōu)于其他模型。

表3 不同特征提取模型評價參數(shù)比較

3 結 論

以地表水COD作為研究對象, 建立了基于KPCA-PSO-ELM結合紫外-可見吸收光譜模型。 對比了ELM模型在PCA、 LASSO、 SPA等特征提取方法下, 模型的表現(xiàn)情況。 基于KPCA-PSO-ELM預測模的訓練集和測試集模型R2均在0.93以上, 測試集RMSE為0.400 7 mg·L-1, 相比于其他模型的RMSE分別降低了78.46%、 18.22%、 2.97%。 結果表明, KPCA-PSO-ELM作為一種結合了非線性特征提取算法的預測模型, 能夠實現(xiàn)對地表水COD快速、 實時的檢測。