基于三維熒光時間序列雙閾值的飲用水污染事件檢測方法研究

薛方家,喻 潔,尹 航,夏戚宇,施杰根,侯迪波,黃平捷,張光新

浙江大學控制科學與工程學院,工業控制技術國家重點實驗室,浙江 杭州 310058

引 言

供水安全是國家安全中極為重要的一環,它是社會正常運作的重要保障之一[1]。近幾十年來,雖然隨著公民素質的提高,生活廢水對水體的污染事件明顯減少,但是偷排工業廢水仍數見不鮮,突發性污染事故更是頻頻發生。根據對國內水污染事件的統計,污染環節主要發生在水源污染、管網污染以及二次供水污染,且污染原因以化學性污染為主,其中有機物污染尤為嚴重[2]。因此為了減少突發性污染產生的損失,迫切需要一種更加快捷實時的水質檢測方法。

飲用水污染檢測的常用方法有化學分析法[3]、氣相色譜法[4]、質譜法[4]、光譜法[5]等方法。相較于其他方法,光譜法無需復雜且耗時的樣品富集、化學試劑添加和其他預處理操作,具有儀器成本更低、檢測操作更簡單、檢測結果獲取更快速、無二次污染等優點。此外,三維熒光靈敏度高、信息豐富,能夠更有效地檢測低濃度有機物污染,并對污染物進行識別[6-8]。

近年來,三維熒光光譜越來越多地被應用于飲用水污染檢測、城市供水突發污染檢測與預警中。Peiris[9]等利用主成分分析法(PCA)對三維熒光原始數據進行特征提取,根據提取的特征值對超濾和納濾飲用水的系統性能進行檢測,實現了對飲用水處理過程中膜污染事件的預警;Heibati[10]等采用平行因子法(PARAFAC)對飲用水的三維熒光數據進行特征提取,發現溶解有機物(DOM)的熒光組分能夠反映飲用水是否受到污染。PCA和PARAFAC能夠有效提取三維熒光光譜圖中的信息,使得三維熒光檢測方法能夠檢測水體中的污染物,但也存在一定的局限性:由于在特征提取過程中,對三維熒光光譜數據進行線性化處理,會損失部分特征信息因此在檢測低濃度污染物時,難以有效區分低強度熒光峰與外界干擾[11]。而飲用水中出現有機污染物質的濃度大多較低(μg·L-1級別),因此常規的基于PCA和PARAFAC的水質異常檢測方法較難實現飲用水突發性污染的有效檢測。

為了解決PCA和PARAFAC對低濃度污染物檢出率低的問題,許多學者提出了其他熒光分析方法來彌補這一缺陷。Huang[12]等使用2-D Gabor小波從三維熒光光譜中提取特征,并將其與SVM結合以識別水中存在的污染物;Yu[13]等采用了對PARAFAC方法進行改進的三線性分解(ATLD)方法來分析正常飲用水樣品的特征,通過殘差矩陣、閾值方法使研究人員能夠確定飲用水是否受到污染;陳方[14]等通過設計了不同環境下的實驗,對不同水質環境得到的水樣本進行分析,對比了基于殘差平方,奇異值分解和三線性分解模型三種特征提取方法,通過對比不同方法在不同水環境中的效果,提高了三維熒光水質檢測在不同污染環境下的準確率。這些方法雖然有效提高了飲用水中低濃度污染物的檢出率,但是對飲用水背景波動缺乏適應性。飲用水三維熒光光譜(背景光譜)會隨著時間、外界環境變化等各種干擾,不斷發生變化。在檢測低濃度污染時,傳統特征提取結合分類器的方法由于未考慮三維熒光光譜在時間維度上的相關性,易受到背景波動和外界干擾,進而造成誤報和漏報。

飲用水背景變化造成的外界干擾,可以分為周期性干擾與隨機性干擾。在時間序列上,干擾與真正污染發生時的三維熒光特征值變化情況存在明顯區別:周期性干擾發生時,特征值周期性超出閾值范圍;隨機性干擾發生時,特征值表現出突然上升突然下降的特點,同時突發性污染事件特征值常存在持續上升,趨于平穩,持續下降的過程?；谕饨绺蓴_與實際污染在時間維度上不同的表現,可以對時間序列三維熒光數據進行時間序列分析,實現污染事件的檢測。此前,于紹慧[15]等利用時間維和光譜維的內在特性,采用時間維補償矩陣和相鄰時段三維熒光矩陣的累計相似度,對按時間順序排列的三維熒光光譜數據中的異常值進行了檢測。采用時間序列累積的方法,雖然能放大異常值的三維熒光信息,但同樣會放大外界干擾,難以有效區分外界干擾與實際污染事件。

針對上述問題,本文提出了一種將三維熒光與時間序列異常檢測相結合的方法,該方法不僅能檢測飲用水體中的低濃度有機物污染,更加能夠解決低濃度污染事件檢測過程中飲用水背景波動的問題。本方法先采用主成分分析法提取三維熒光光譜的前兩個主元特征值,訓練線性自回歸(AR)模型并對未來時段主元特征值進行預測,與實測主元特征值作差得到特征值差值,用于排除周期性干擾的影響;并對實測特征值變化率進行計算,用于排除隨機性干擾的影響。最終設置特征值差值-特征值變化率雙閾值確定污染的起始點與結束點,從而確定整個污染事件。并且設計實驗,選取苯酚作為模擬污染物,分別模擬了高濃度有機物突發性污染、低濃度有機物突發性污染以及飲用水背景大幅度波動下的低濃度有機物突發性污染,并將時間序列雙閾值方法與傳統的三維熒光檢測方法在檢測準確率上進行了對比,證明了方法的有效性。

1 實驗部分

1.1 實驗設計

為了證實本文提出的時間序列雙閾值法對飲用水特征污染物進行定性判別檢測的有效性,本文搭建了飲用水污染模擬事件檢測系統,如圖1所示,其包含了污染模擬系統和光譜檢測系統。污染模擬系統中,飲用水和污染物分別由A、B泵抽取混合,模擬突發性污染事件,并可以通過改變雙泵流量,實時控制污染物濃度;光譜檢測系統中,熒光測量儀器采用Horiba公司的Aqualog為檢測設備,激發波長設置范圍為240～800 nm,發射波長設置范圍為243.544～823.84 nm,積分時間0.1 s,每90 s就可以實現一次三維熒光光譜采樣,可以得到時間序列三維熒光光譜數據。

圖1 實驗裝置示意圖Fig.1 Schematic diagram of experimental device

在污染物選取方面,苯酚作為飲用水中常見的污染物,已經被我國《生活飲用水衛生標準》(GB5749—2022)列為生活飲用水水質的擴展指標[16],因此實驗選取苯酚作為模擬污染物,依據《生活飲用水衛生標準》(GB5749—2022)中苯酚限值,設計了三組不同環境下模擬突發性污染實驗進行驗證。

第一組實驗用于模擬突發性污染發生在檢測點附近的高濃度污染事件,實驗中苯酚濃度梯度設置為40～200 μg·L-1,共5個污染事件,實驗共進行120 min,獲得三維熒光光譜數據35個。

第二組實驗用于模擬突發性污染發生點距離檢測點較遠,由擴散導致的低濃度污染事件,實驗中苯酚濃度梯度設置為5～20 μg·L-1,共3個污染事件,獲得三維熒光光譜數據85個。

第三組實驗選用sipper實驗裝置直接從管道取水,由于sipper裝置內水流量不斷發生變化,外界干擾也明顯大于之前兩組實驗,與真實管道相類似,用于模擬水質背景波動下的低濃度污染事件,實驗中苯酚濃度設置為2 μg·L-1,共1個污染事件,獲得三維熒光光譜數據137個,其中前40個數據點用于訓練AR模型,后97個數據點用于測試算法有效性。

1.2 方法

算法的流程如圖2所示,預處理實驗數據得到的三維熒光在線數據,使用主成分分析法(PCA)進行降維得到主元特征值數據,之后通過線性自回歸(AR)方法進行特征值預測,得到特征值差值與特征值變化率,根據設定的閾值確定污染事件起始和結束點,從而確定整個污染事件,并將檢測的準確率與支持向量機檢測結果,以及基于支持向量機(SVM)的時間序列修正結果進行對比,對算法進行評估。

圖2 方法流程圖Fig.2 Algorithm framework diagram

1.2.1 數據預處理

從三維熒光儀器中得到的數據,需要經過去散射和歸一化預處理:

(1)去散射:在三維熒光光譜的測量中,由于復雜的能級躍遷,瑞利散射和拉曼散射同時發生,散射光的干擾嚴重影響了光譜的靈敏度及光譜數據分析[17],會對污染物檢測造成影響,而通過去散射可以使污染物的特征峰更加明顯,增強污染檢測的精度。本文采用插值法對原始光譜數據進行預處理來減少瑞利散射對檢測的影響,并采用扣除空白溶劑背景的方法消除拉曼散射的影響。

(2)歸一化:三維熒光儀測得的三維熒光光譜數據會出現小于0的情況,三維熒光強度并不會小于0,因此屬于誤報,需要將三維熒光數據小于0的值改成0。而對于不同時間獲得的水樣,容易受到溫度,濕度以及其他外界環境干擾造成數據不統一,因此在完成每組實驗之前先檢測純水水拉曼值,每組數據除以當日的水拉曼值后進行歸一化處理。在水拉曼處理后,再進行標準化處理此處采用min-max標準化,其數學表達式如式(1)所示

(1)

式(1)中,X′為歸一化處理后的數據,X為原數據,data_max為三維熒光數據中的最大值,data_min為三維熒光數據中的最小值。

1.2.2 主成分分析法

主成分分析法(PCA)[18]是一種數據降維算法,通過正交變換將一組存在相關性的變量轉換為一組線性不相關的變量,轉換后的這組變量稱為主成分,主成分按其方差大小排序,第一成分含有原始數據最多的信息,第二成分含原始數據次多信息,依次類推,保留低階主成分就等同于保留了原始數據主要特征信息。假設有m個n維向量xi,它們之間存在一定相關性,要提取它們主要特征信息,PCA算法實現如下:

1.2.3 特征值差值-特征值變化率雙閾值方法

特征值差值-特征值變化率雙閾值算法基于干擾與真實污染事件不同的特點進行設計,在時間全維度上,外界干擾發生時特征值表現出突然上升突然下降的特點,而突發性污染事件發生時特征值常有持續上升,趨于平穩,持續下降的過程。其具體判斷流程如下:

(2)

Qt=(Yt-Yt-1)/Δt

Qt+1=(Yt+1-Yt)/Δt

(3)

式(3)中,Δt為三維熒光光譜采樣間隔時間。

(4)

式(4)中,P,Q1,Q2分別為特征值差值以及特征值變化率的閾值,本文中,P=0.5,Q1=0.5 min-1,Q2=-0.5 min-1。

依照以上算法,可以對污染發生和結束點進行準確判斷,從而確定整個污染事件,算法對污染起始點的判斷邏輯框架圖如圖3所示。

圖3 算法邏輯框架圖Fig.3 Algorithm logic framework diagram

2 結果與討論

2.1 高濃度污染檢測結果

高濃度污染檢測結果如圖4所示,可以看出,由于污染物濃度較高,各方法都能很好地檢測出突發性污染事件。該實驗用于模擬突發性污染發生在污染監測點附近的情況,苯酚濃度梯度設置為40～200 μg·L-1,共設置5個污染事件,各個事件苯酚濃度分別為200、160、120、80和40 μg·L-1。

圖4 高濃度污染檢測結果Fig.4 High-level pollution test results

2.2 低濃度污染檢測結果

低濃度污染檢測結果如圖5所示,可以看出,傳統SVM分類由于低濃度苯酚的三維熒光光譜特征峰較低,容易與自來水中受外界干擾出現的干擾峰相混淆,因此出現大量低濃度苯酚污染被誤報的情況。該實驗用于模擬突發性污染發生點距離檢測點較遠,由擴散導致的低濃度污染事件,實驗中苯酚濃度梯度設置為5～20 μg·L-1,共3個污染事件。

圖5 低濃度污染檢測結果Fig.5 Low-level pollution test results

對于低濃度的突發性污染中傳統SVM分類結果進行分析后可以發現,很多誤報漏報點都以孤立的形式出現,如圖6,在低濃度污染檢測結果中,A、B、C點就是明顯的誤報漏報點?？梢愿鶕r間序列上SVM對前后幾個時間點是否發生污染的分類結果,對該點是否存在污染進行重新判斷,如式(5)所示

圖6 低濃度污染檢測誤報點Fig.6 False alarm point of low-level pollution test

rt=0.4ct+0.2ct-1+0.2ct-2+0.2ct+1

(5)

式(5)中,ct為原分類結果,rt為修正值,1代表存在污染,0代表無污染。

時間序列修正的過程如圖7所示,經過簡單的時間序列修正,如表1所示,低濃度污染檢測結果準確率明顯上升。

表1 修正準確率對比Table 1 Correction accuracy comparison

圖7 時間序列修正過程圖Fig.7 Time-series correction process

2.3 水質背景波動下的低濃度污染檢測結果

水質背景波動選取模擬真實管道作為實驗對象,在水流量不斷變化、高外界干擾的情況下,實驗中飲用水背景三維熒光光譜圖在外界干擾下不斷發生變化,如圖8所示,四張圖均為無污染發生時的正常飲用水背景三維熒光光譜圖,但是其特征峰明顯存在不同,因此PCA提取得到的主元特征值也隨著時間不斷變化,如圖10(a)所示,無法采用SVM對污染事件和背景波動正常值進行有效區分。

圖8 水質背景三維熒光光譜波動圖(a)—(d)為不同時間段的飲用水光譜圖Fig.8 Three dimensional fluorescence spectrum fluctuation of water quality background(a)—(d) are the EEM of drinking water in different time periods

在這種環境中,各算法檢測結果如圖9所示,可以看出,在高干擾和污染物濃度很低的情況下,傳統SVM方法出現大量誤報和漏報,難以有效檢測污染事件。

圖9 水質背景波動下的低濃度污染檢測結果Fig.9 Low concentration pollution detection results under fluctuating water quality background

因此需要對高干擾環境使用時間序列雙閾值方法進一步改善準確性。以主元貢獻率最高的兩個特征值為對象,建立線性自回歸模型,綜合考慮主元特征值閾值與變化率進行異常檢測。此處將前40個未發生污染的數據選作訓練集,選取其PCA貢獻度最高的兩個主元作為輸入,經過最小二乘法,確定AR模型如式(6)所示,其中貢獻度最高的主元特征值預測結果與實測值對比圖如圖10(b)所示,主元特征值變化率如圖10(c)所示

圖10 主元特征值圖像(a):主元特征值變化曲線;(b):正常水背景主元特征值預測結果;(c):主元特征值變化率Fig.10 The image of principal component(a):The change curve of principal component;(b):The prediction result of principal component in the normal water background;(c):The change rate of principal component

(6)

從預測與實測值曲線可以看出,在未發生污染事件時,預測結果與實測值擬合效果較好,污染發生后,兩者出現了明顯的偏差,先采用閾值法對苯酚污染進行檢測,綜合考慮第一主元與第二主元,預測值與實測值差值計算公式如式(7)所示

(7)

表2 雙閾值檢測準確率對比Table 2 Comparison of double threshold detection accuracy

如表3所示,高干擾環境采用時間序列雙閾值檢測方法后,檢測準確率相較簡單的時間序列修正上升了11.4%,尤其是對非邊緣的誤報區有了很大程度的改善,能夠有效地將污染事件與水樣本身波動區分開來,檢測背景波動較大環境下的低濃度污染事件。

表3 不同算法在不同環境下準確率對比表Table 3 Comparison table of accuracy of different algorithms in different environments

3 結 論

針對三維熒光方法檢測低濃度污染物質時,傳統判別方法容易受到水樣本身波動以及檢測設備噪聲等影響,進而產生污染事件誤報和漏報等不足,提出了一種基于時間序列雙閾值法的三維熒光水質異常事件檢測模型,通過挖掘三維熒光數據時間維度信息,設置特征值差值、特征值變化率閾值,有效解決了三維熒光在污染檢測過程中受水質背景波動影響等問題。實驗結果表明,該方法不僅在高濃度污染事件中檢測準確,在低濃度污染、水質背景波動下的低濃度污染檢測中,準確率分別達到了98.8%和99.0%,相較于傳統的判別方法,檢測準確率明顯上升,說明本方法可以有效檢測水質背景波動下的低濃度污染事件,為三維熒光在線監測飲用水的有機物突發性污染提供了幫助。