基于小波分析的水質(zhì)變化及預(yù)測研究
——以渦河為例

趙貴章，王淑麗，李志萍*，龔建師，王赫生

(1.華北水利水電大學(xué)地球科學(xué)與工程學(xué)院，河南 鄭州 450011;2.中國地質(zhì)調(diào)查局南京地質(zhì)調(diào)查中心，江蘇 南京 210016)

由于水環(huán)境系統(tǒng)具有高度復(fù)雜性，且水質(zhì)指標(biāo)數(shù)據(jù)具有不確定性、非線性等特征，因此不斷進(jìn)行水質(zhì)變化規(guī)律分析和預(yù)測方法改進(jìn)是必要的[1]。渦河作為淮河流域第二大支流，水質(zhì)的變化嚴(yán)重影響著淮河水質(zhì)及沿岸地區(qū)國民經(jīng)濟(jì)的進(jìn)一步發(fā)展。通過分析渦河流域長期水質(zhì)監(jiān)測數(shù)據(jù)，掌握流域水質(zhì)的時間格局特征，進(jìn)而對水質(zhì)進(jìn)行有效預(yù)測，對深入了解渦河水資源循環(huán)規(guī)律，保護(hù)渦河流域生態(tài)環(huán)境至關(guān)重要。

小波分析具有時頻局部化特征，它可以揭示水文序列的周期變化特征和多尺度變化規(guī)律[2]，BP神經(jīng)網(wǎng)絡(luò)能大規(guī)模學(xué)習(xí)和存儲輸入輸出模式的映射關(guān)系，且能夠較好處理非線性系統(tǒng)，符合水質(zhì)預(yù)測的要求[3]。因此，采用小波多分辨分析能力和神經(jīng)網(wǎng)絡(luò)的自學(xué)習(xí)能力特性相結(jié)合的方式，可以更好地探究水質(zhì)變化規(guī)律。國內(nèi)外學(xué)者也將小波分析應(yīng)用于水質(zhì)時間序列分析的研究中。一部分學(xué)者通過小波分析對水質(zhì)序列原始數(shù)據(jù)進(jìn)行降噪處理，剔除干擾數(shù)據(jù)，發(fā)現(xiàn)能夠提高預(yù)測結(jié)果的準(zhǔn)確度；另一部分學(xué)者研究了小波分析在水質(zhì)變化規(guī)律中的應(yīng)用，探討了水質(zhì)指標(biāo)時間序列的周期模式以及時間格局特征[4-6]；還有一部分學(xué)者將小波分析與神經(jīng)網(wǎng)絡(luò)相結(jié)合對水質(zhì)指標(biāo)進(jìn)行預(yù)測，發(fā)現(xiàn)該種組合模式預(yù)測精度更高[7-8]。

小波分析具有多分辨分析的特點，可以對信號進(jìn)行不同尺度上的分解，通過小波分解得到的小波系數(shù)、小波方差等信息，可以用于識別各水質(zhì)指標(biāo)時間序列的周期變化規(guī)律及突變點，該方法對水質(zhì)時間特征的演化研究更精細(xì)。神經(jīng)網(wǎng)絡(luò)與小波分析相結(jié)合是國際上研究的熱點之一。小波神經(jīng)網(wǎng)絡(luò)既有神經(jīng)網(wǎng)絡(luò)的自學(xué)習(xí)能力特性，又有小波的局部特性，彌補(bǔ)了以往神經(jīng)網(wǎng)絡(luò)預(yù)測時不分主次的不足，提高了網(wǎng)絡(luò)的穩(wěn)定性。本研究以渦河流域水質(zhì)指標(biāo)為研究對象，以期為小波分析在水質(zhì)序列變化規(guī)律及預(yù)測方面的應(yīng)用提供參考。

1 研究區(qū)概況

渦河是淮河第二大支流，也是淮北平原區(qū)主要河道，起源于河南省開封市賈魯河，流經(jīng)開封、通許、扶溝、太康、鹿邑和安徽省亳州等地，最終注入淮河，具體位置見圖1[9]。渦河流域地勢平坦，河道坡降較小，受季風(fēng)氣候影響，區(qū)域年降水量分配不均，主要集中在6—8月，多年平均年降水量765 mm，多年平均徑流量1.5億m3。20世紀(jì)以來，隨著工農(nóng)業(yè)生產(chǎn)的快速發(fā)展，污水大量排入河流，導(dǎo)致渦河入境水質(zhì)不佳，渦河水體受到嚴(yán)重污染。進(jìn)入21世紀(jì)，伴隨著河長制的推行，以及渦河生態(tài)建設(shè)示范帶的推進(jìn)，渦河的流域水質(zhì)明顯好轉(zhuǎn)，目前主要是Ⅱ、Ⅲ類水。

圖1 渦河流域位置

2 研究方法

2.1 小波分析

小波分析可以對時間和頻率進(jìn)行局部化分析，最終實現(xiàn)高頻的時間細(xì)分和低頻的頻率細(xì)分。它是繼傅里葉變換以來科學(xué)方法的重大突破。在實際場景中，有效信號頻率較低且相對穩(wěn)定，而噪聲頻率一般較高。小波變換就是找到一組低通和高通濾波器組，對原始信號進(jìn)行連續(xù)向下分解，得到一系列具有不同分辨率的低頻概貌信號caN和高頻細(xì)節(jié)信號cd1,…,cdN，從而大大增強(qiáng)了局部信號。對這些低頻和高頻信號進(jìn)行重構(gòu)，得到去噪后的信號，提高了模型的抗噪性。

小波是一種特殊的波形，通常用小波母函數(shù)ψ(τ)來表示，對ψ(τ) 進(jìn)行一定的時間上的平移和尺度上的伸縮，就可以得到小波函數(shù)基：

(1)

(2)

式中b——平移因子；a——尺度因子。

小波變換的含義是尺度參數(shù)a和平移參數(shù)b的函數(shù)，是一個時間-尺度分析。連續(xù)小波變換公式:

(3)

由于連續(xù)小波變換的過程會產(chǎn)生許多冗余數(shù)據(jù)，在實際使用中連續(xù)小波必須加以離散化。因此將尺度參數(shù)a和平移參數(shù)b離散化后得到離散小波變換。通常情況下，離散化公式分別取作：

(4)

式中，a0>0,j∈Z。對應(yīng)的離散小波基可表示為：

(5)

離散小波變換公式：

(6)

最后,將所有低頻信號和高頻信號相加實現(xiàn)重構(gòu)，還原公式如下:

(7)

式中S(t)——還原之后的數(shù)據(jù)；ψj,k(t)——低通濾波器；h(ψj,k(t))——高通濾波器。

本研究采用常用的Morlet小波，通過小波分解，得到小波系數(shù)、小波方差等信息，識別渦河流域2005—2018年的水質(zhì)指標(biāo)多時間尺度變化特征及其突變特征。Daubechies(dbN)正交小波基，由于其良好的時頻分析性能，已在眾多領(lǐng)域應(yīng)用[10]?；贒aubechies(dbN)小波，對主要影響指標(biāo)進(jìn)行小波多分辨率分析，為小波神經(jīng)網(wǎng)絡(luò)預(yù)測模型的建立進(jìn)行預(yù)處理[7,11-13]。

2.2 BP神經(jīng)網(wǎng)絡(luò)

BP神經(jīng)網(wǎng)絡(luò)能夠較好處理模糊的非線性關(guān)系，因而可以很好適用于水質(zhì)預(yù)測等復(fù)雜問題的研究[14]。本研究進(jìn)一步探討，基于小波分析的BP神經(jīng)網(wǎng)絡(luò)模型在水質(zhì)預(yù)測方面的適用性。所采用的是反向傳播人工神經(jīng)網(wǎng)絡(luò)BPNN(Back-Propagation Artificial Neural Network)，其典型結(jié)構(gòu)為3層網(wǎng)絡(luò)模型，即輸入層-隱含層-輸出層，具體原理見圖2。

輸入層、隱含層與輸出層之間的神經(jīng)元連接可描述為：

a1=f(IW·Pi+b1)

(8)

a2=f(LW·Pi+b2)

(9)

式中a1、a2——隱含層和輸出層的計算值；Pi——輸入量；IW、LW——隱含層連接輸入層和輸出層之間的權(quán)重；b1、b2——隱含層與輸入層和輸出層之間的偏差值。

其中，隱含層的傳遞函數(shù)采用雙曲正切S型函數(shù)：

(10)

輸出層的傳遞函數(shù)為純線性函數(shù)：

f(x)=x

(11)

2.3 技術(shù)方法

本文通過小波分析對渦河流域水質(zhì)指標(biāo)進(jìn)行多尺度變化規(guī)律分析，并結(jié)合BPNN神經(jīng)網(wǎng)絡(luò)對其進(jìn)行預(yù)測。首先，使用主成分分析方法確定關(guān)鍵水質(zhì)參數(shù)，這作為水質(zhì)預(yù)測的目標(biāo)值，提高預(yù)測效率。再通過相關(guān)性分析方法將與目標(biāo)參數(shù)相關(guān)性較強(qiáng)的水質(zhì)參數(shù)篩選出來并作為預(yù)測目標(biāo)值的輸入值，提高預(yù)測精度。最后通過Morlet小波分析渦河流域水質(zhì)指標(biāo)多尺度變化規(guī)律，為水質(zhì)預(yù)測周期提供參考。

在模擬實測水質(zhì)數(shù)據(jù)時，需要預(yù)處理工具去噪來提高模型的精度，而小波分解能將原始數(shù)據(jù)中不同頻段的信息進(jìn)行分解，極大地降低數(shù)據(jù)復(fù)雜度，再分別對這些數(shù)據(jù)進(jìn)行預(yù)測從而提高預(yù)測精度。因此，通過比較幾個Daubechies小波，選用最為合適的小波基和分解尺度對目標(biāo)參數(shù)的時間序列數(shù)據(jù)進(jìn)行降噪處理，得到低頻概貌分量(caN)和多個高頻細(xì)節(jié)分量(cdN)。

通過BPNN神經(jīng)網(wǎng)絡(luò)，分別對多個分量進(jìn)行預(yù)測，并將其預(yù)測結(jié)果相加即為目標(biāo)值的最終預(yù)測結(jié)果，BPNN模型均采用3層網(wǎng)絡(luò)結(jié)構(gòu)，主要采用自動完成的“試錯法”來訓(xùn)練模型并獲取最佳的模擬效果，可調(diào)整的主要參數(shù)包括隱含層神經(jīng)元數(shù)目、迭代次數(shù)、動量因子和學(xué)習(xí)速率等。具體的技術(shù)路線見圖3。

圖3 技術(shù)路線

3 數(shù)據(jù)處理及分析

3.1 數(shù)據(jù)來源

本研究數(shù)據(jù)來源于渦河流域中下游的4個監(jiān)測斷面2005—2018年(共168個月)的水質(zhì)指標(biāo)每月一次監(jiān)測值，其中玄武閘斷面位于河南省鹿邑縣，大寺閘斷面、渦陽閘斷面和蒙城閘斷面位于安徽省內(nèi)，具體斷面位置見圖1。應(yīng)用的水質(zhì)指標(biāo)為溶解氧(DO)、高錳酸鹽指數(shù)(KMnO4)、氨氮(NH4-N)、化學(xué)需氧量(COD)和總磷(TP)共5項，其監(jiān)測數(shù)據(jù)完整且連續(xù)。

3.2 主成分分析

為提高預(yù)測效率，對溶解氧、高錳酸鹽指數(shù)、氨氮、化學(xué)需氧量和總磷5個水質(zhì)指標(biāo)進(jìn)行主成分分析選取預(yù)測參數(shù)。具體主成分分析結(jié)果見表1。由表1可知，第一主成分方差貢獻(xiàn)率為47.56%，其反映了原變量的大部分信息，其中化學(xué)需氧量、氨氮和高錳酸鹽指數(shù)得分較高，且化學(xué)需氧量得分最高；第一和第二主成分方差累計貢獻(xiàn)率為 69.71%，其中溶解氧得分最高為0.81；第三主成分中，總磷得分最高達(dá)到了1.07。

表1 主成分分析結(jié)果

目前影響渦河流域水質(zhì)的主要因子為以化學(xué)需氧量為代表的，包括氨氮、高錳酸鹽指數(shù)在內(nèi)的第一主成分因子。它反映出的信息一方面為有機(jī)物指標(biāo)，如化學(xué)需氧量和高錳酸鹽指數(shù)，另一方面為水體營養(yǎng)鹽的指標(biāo)，如氨氮。說明目前導(dǎo)致渦河流域污染的主要原因為有機(jī)污染，因此利用主要污染因子，即用化學(xué)需氧量的變化來反映渦河流域污染的變化是可行的。通過分析化學(xué)需氧量與各指標(biāo)間的相關(guān)性，為化學(xué)需氧量的預(yù)測選擇合適的輸入?yún)?shù)，具體相關(guān)性分析結(jié)果見表2。由表2可知，化學(xué)需氧量與高錳酸鹽指數(shù)、氨氮的相關(guān)性較強(qiáng)，相關(guān)系數(shù)分別達(dá)到0.604、0.650，因此將其作為化學(xué)需氧量預(yù)測模型的輸入變量。

表2 各指標(biāo)相關(guān)關(guān)系

3.3 小波分析

基于主成分分析法選取了渦河流域水質(zhì)的主要影響因子，即化學(xué)需氧量。通過對比分析選擇合適的分解尺度對化學(xué)需氧量原序列進(jìn)行分解，具體對比結(jié)果見圖4。

a)小波基(s)

由圖4a—4e可以看出，隨著分解尺度N的增加，曲線變得更平滑，并趨近原序列，但同時會丟失許多細(xì)節(jié)成分，當(dāng)N=5時，曲線太過平滑，無法反映原始序列的峰值。綜合比較表明，db2既能顯示化學(xué)需氧量原序列s的大致概貌，細(xì)節(jié)成分保留也較完整，而且能很好地反映峰值。由于國內(nèi)外缺乏最佳小波基的規(guī)范性選取方法，但不同小波基的預(yù)測結(jié)果差別不大，因此，綜合比較選擇db2小波基對化學(xué)需氧量序列進(jìn)行小波分解。

由圖4f—4i可以看出，低頻概貌部分體現(xiàn)了原始信號的基本形狀,可以看成是原始序列的一個近似。隨著分解尺度N的增大，分解序列變得越來越光滑，然而當(dāng)N=3時，曲線太光滑，無法表現(xiàn)原序列的總體概貌。綜合比較當(dāng)N取2時，最適合于化學(xué)需氧量時間序列的小波分析。此時化學(xué)需氧量原序列s分解如下:

s=ca2+cd2+cd1

(12)

式中 ca2——趨勢成分；cd2——周期成分；cd1——隨趨成分。

3.4 神經(jīng)網(wǎng)絡(luò)分析

構(gòu)建化學(xué)需氧量小波BPNN模型時，其中ca2分量為原序列的趨勢項，應(yīng)考慮與化學(xué)需氧量變化有關(guān)的相關(guān)變量，因此選取與主要影響因子相關(guān)系數(shù)較大的指標(biāo)，即高錳酸鹽指數(shù)和氨氮作為預(yù)測ca2的輸入變量，ca2為輸出變量；cd2為原序列的周期項，cd1為原序列的隨趨成分，與各水質(zhì)參數(shù)相關(guān)性不大，基本只與原始序列自身的變化有關(guān)，因此選擇化學(xué)需氧量本身作為其輸入，cd2和cd1分別作為輸出，然后將所有分變量預(yù)測結(jié)果疊加作為化學(xué)需氧量最終預(yù)測結(jié)果。經(jīng)過多次試驗測試與參數(shù)調(diào)整,在保證網(wǎng)絡(luò)快速收斂的同時又具有較高的預(yù)測精度,使模型達(dá)到最優(yōu)結(jié)果。訓(xùn)練擬合中db2-2相關(guān)參數(shù)的最終配置結(jié)果見表3。

表3 模型關(guān)鍵參數(shù)

4 結(jié)果與討論

4.1 水質(zhì)指標(biāo)周期變化規(guī)律分析

基于Morlet小波，對渦河流域河流2005—2018年，共168個月的溶解氧、高錳酸鹽指數(shù)、化學(xué)需氧量、氨氮和總磷進(jìn)行多尺度變化分析。通過小波實部等值線圖和小波方差來識別不同水質(zhì)指標(biāo)的主周期及變化特征，具體結(jié)果見圖5。

a)溶解氧等值線

e)化學(xué)需氧量等值線

由圖5小波系數(shù)等值線可以看出，不同水質(zhì)指標(biāo)序列存在以不同尺度隨時間交替變化的特性，圖5小波方差可看出，各指標(biāo)的小波方差圖的主要峰值所對應(yīng)的時間尺度與小波實部等值線圖所顯示的波動變化規(guī)律相吻合。從圖4可知：渦河流域溶解氧8、18、33、64個月尺度交替變化表現(xiàn)較清晰，波動極值點分布規(guī)律明顯，以18個月尺度較為顯著，說明渦河流域溶解氧系列存在以8、18、33、64個月變化的主要周期；高錳酸鹽指數(shù)存在以9、20、31、50個月變化主要周期，以20個月尺度較為顯著；化學(xué)需氧量存在以8、20、31、43個月尺度變化主要周期，以31個月尺度較為顯著；氨氮存在以8、20、31、57個月尺度變化主要周期，以57個月尺度較為顯著；總磷存在以8、12、23個月變化主要周期，以23個月尺度較為顯著。綜合來看，各水質(zhì)指標(biāo)具有多尺度振蕩的特點，且主要存在以8、20、31個月變化的主周期，各水質(zhì)指標(biāo)主周期變化主要發(fā)生在汛期。

為探究主周期形成的主要原因，基于各指標(biāo)多年平均值來分析各指標(biāo)年內(nèi)變化趨勢，具體結(jié)果見圖6。由圖6可以看出，各指標(biāo)在汛期(7—9月)含量明顯低于年內(nèi)其他月份，這可能是由于降雨量的增加，河段水流量也隨之增加，會對污染物產(chǎn)生稀釋作用，各指標(biāo)含量降低，因而各指標(biāo)呈現(xiàn)出汛期周期變化規(guī)律。

圖6 各指標(biāo)年內(nèi)變化趨勢

4.2 水質(zhì)預(yù)測研究

基于上文的各個小波基和分解尺度對比分析結(jié)果，最終選擇db2小波基對化學(xué)需氧量原序列s進(jìn)行2層尺度的多分辨率分析，得到低頻概貌分量ca2和高頻細(xì)節(jié)分量cd2及cd1，見圖7。

a)s

由圖7可知，化學(xué)需氧量1年中的變化情況大致為夏季高、冬季低，春秋相當(dāng)，但年際情況不盡相同。ca2與原序列最為接近，體現(xiàn)了化學(xué)需氧量的總體下降趨勢,且在2005—2009年(第1—50月)波動較大，但具體細(xì)節(jié)變化未得到很好體現(xiàn)。cd2盡管為細(xì)節(jié)部分，但對原序列的一些峰值如第40個月仍有著較好的表達(dá)。cd1在峰值處呈現(xiàn)較大的波動，但更多地含有隨機(jī)成分。

考慮到化學(xué)需氧量時間序列的長度，以及汛期周期變化的規(guī)律，最終選擇以一整年12個月作為預(yù)測周期。即以2005—2017年(第1—156月)的月均化學(xué)需氧量數(shù)據(jù)作為訓(xùn)練對象進(jìn)行訓(xùn)練，以2018年各月化學(xué)需氧量含量(第157—168月)作為最終預(yù)測對象。預(yù)測結(jié)果的誤差分析見表4，預(yù)測值與實測值的擬合結(jié)果見圖8。

表4 化學(xué)需氧量第157—168月預(yù)測結(jié)果誤差分析

圖8 化學(xué)需氧量預(yù)測值與實測值對照結(jié)果

由表4可知，在預(yù)測期內(nèi)平均相對誤差控制在8.4%，均方根誤差為1.5。2018年1—6月(第157—162月)誤差較小，7—10月(第163—166月)誤差相對較大，7月份誤差最大，可能與降雨量大小有關(guān)，因此在提高水質(zhì)預(yù)測精度時，可以考慮加入季節(jié)性因素。圖8可知，預(yù)測值與實測值的曲線變化趨勢一致，尤其在短時預(yù)測時段擬合更優(yōu)，總體來看，構(gòu)建的小波神經(jīng)網(wǎng)絡(luò)模型預(yù)測化學(xué)需氧量精度較高，且較為穩(wěn)定，能夠適用于渦河流域的短期水質(zhì)預(yù)測。

5 結(jié)論

基于渦河流域2005—2018年的水質(zhì)指標(biāo)月監(jiān)測數(shù)據(jù)，即溶解氧、氨氮、高錳酸鹽指數(shù)和化學(xué)需氧量，探究小波神經(jīng)網(wǎng)絡(luò)在流域水質(zhì)污染中的應(yīng)用。通過Morlet小波分析了渦河流域水質(zhì)指標(biāo)多尺度變化規(guī)律，并基于Daubechies小波和BPNN神經(jīng)網(wǎng)絡(luò)建立了小波神經(jīng)網(wǎng)絡(luò)預(yù)測模型，該研究結(jié)果證實了小波分析在水質(zhì)時間序列變化規(guī)律及預(yù)測方面的可行性。

a)渦河流域水質(zhì)指標(biāo)的周期變化規(guī)律表明，總磷有3個主要變化周期，其他指標(biāo)均具有4個變化主周期，各水質(zhì)指標(biāo)變化主周期存在一定差異，但主要存在以8、20、30個月左右變化的主周期，各水質(zhì)指標(biāo)主周期變化主要發(fā)生在汛期；目前渦河流域水質(zhì)的主要影響因子是以化學(xué)需氧量為代表的，包括氨氮、高錳酸鹽指數(shù)在內(nèi)的第一主成分因子。

b)選擇db2小波基對化學(xué)需氧量原序列進(jìn)行了2層尺度的多分辨率分解；通過小波神經(jīng)網(wǎng)絡(luò)得到的化學(xué)需氧量預(yù)測值與實測值的曲線擬合度較高，平均相對誤差為8.4%，均方根誤差為1.5，該方法較適用于渦河流域的短期水質(zhì)預(yù)測，對于水質(zhì)預(yù)測具有較好的參考價值。

c)由于國內(nèi)外缺乏最佳小波基的規(guī)范性選取方法，且研究發(fā)現(xiàn)季節(jié)性因素對水質(zhì)序列預(yù)測誤差影響較大，因而在今后提高小波BP神經(jīng)網(wǎng)絡(luò)預(yù)測精度時可以考慮從小波基的優(yōu)化選擇和季節(jié)性因素的輸入出發(fā)。