趙貴章,王淑麗,李志萍*,龔建師,王赫生
(1.華北水利水電大學地球科學與工程學院,河南 鄭州 450011;2.中國地質調查局南京地質調查中心,江蘇 南京 210016)
由于水環境系統具有高度復雜性,且水質指標數據具有不確定性、非線性等特征,因此不斷進行水質變化規律分析和預測方法改進是必要的[1]。渦河作為淮河流域第二大支流,水質的變化嚴重影響著淮河水質及沿岸地區國民經濟的進一步發展。通過分析渦河流域長期水質監測數據,掌握流域水質的時間格局特征,進而對水質進行有效預測,對深入了解渦河水資源循環規律,保護渦河流域生態環境至關重要。
小波分析具有時頻局部化特征,它可以揭示水文序列的周期變化特征和多尺度變化規律[2],BP神經網絡能大規模學習和存儲輸入輸出模式的映射關系,且能夠較好處理非線性系統,符合水質預測的要求[3]。因此,采用小波多分辨分析能力和神經網絡的自學習能力特性相結合的方式,可以更好地探究水質變化規律。國內外學者也將小波分析應用于水質時間序列分析的研究中。一部分學者通過小波分析對水質序列原始數據進行降噪處理,剔除干擾數據,發現能夠提高預測結果的準確度;另一部分學者研究了小波分析在水質變化規律中的應用,探討了水質指標時間序列的周期模式以及時間格局特征[4-6];還有一部分學者將小波分析與神經網絡相結合對水質指標進行預測,發現該種組合模式預測精度更高[7-8]。
小波分析具有多分辨分析的特點,可以對信號進行不同尺度上的分解,通過小波分解得到的小波系數、小波方差等信息,可以用于識別各水質指標時間序列的周期變化規律及突變點,該方法對水質時間特征的演化研究更精細。神經網絡與小波分析相結合是國際上研究的熱點之一。小波神經網絡既有神經網絡的自學習能力特性,又有小波的局部特性,彌補了以往神經網絡預測時不分主次的不足,提高了網絡的穩定性。本研究以渦河流域水質指標為研究對象,以期為小波分析在水質序列變化規律及預測方面的應用提供參考。
渦河是淮河第二大支流,也是淮北平原區主要河道,起源于河南省開封市賈魯河,流經開封、通許、扶溝、太康、鹿邑和安徽省亳州等地,最終注入淮河,具體位置見圖1[9]。渦河流域地勢平坦,河道坡降較小,受季風氣候影響,區域年降水量分配不均,主要集中在6—8月,多年平均年降水量765 mm,多年平均徑流量1.5億m3。20世紀以來,隨著工農業生產的快速發展,污水大量排入河流,導致渦河入境水質不佳,渦河水體受到嚴重污染。進入21世紀,伴隨著河長制的推行,以及渦河生態建設示范帶的推進,渦河的流域水質明顯好轉,目前主要是Ⅱ、Ⅲ類水。

圖1 渦河流域位置
小波分析可以對時間和頻率進行局部化分析,最終實現高頻的時間細分和低頻的頻率細分。它是繼傅里葉變換以來科學方法的重大突破。在實際場景中,有效信號頻率較低且相對穩定,而噪聲頻率一般較高。小波變換就是找到一組低通和高通濾波器組,對原始信號進行連續向下分解,得到一系列具有不同分辨率的低頻概貌信號caN和高頻細節信號cd1,…,cdN,從而大大增強了局部信號。對這些低頻和高頻信號進行重構,得到去噪后的信號,提高了模型的抗噪性。
小波是一種特殊的波形,通常用小波母函數ψ(τ)來表示,對ψ(τ) 進行一定的時間上的平移和尺度上的伸縮,就可以得到小波函數基:
(1)
(2)
式中b——平移因子;a——尺度因子。
小波變換的含義是尺度參數a和平移參數b的函數,是一個時間-尺度分析。連續小波變換公式:
(3)

由于連續小波變換的過程會產生許多冗余數據,在實際使用中連續小波必須加以離散化。因此將尺度參數a和平移參數b離散化后得到離散小波變換。通常情況下,離散化公式分別取作:
(4)
式中,a0>0,j∈Z。對應的離散小波基可表示為:
(5)
離散小波變換公式:
(6)

最后,將所有低頻信號和高頻信號相加實現重構,還原公式如下:
(7)
式中S(t)——還原之后的數據;ψj,k(t)——低通濾波器;h(ψj,k(t))——高通濾波器。
本研究采用常用的Morlet小波,通過小波分解,得到小波系數、小波方差等信息,識別渦河流域2005—2018年的水質指標多時間尺度變化特征及其突變特征。Daubechies(dbN)正交小波基,由于其良好的時頻分析性能,已在眾多領域應用[10]。基于Daubechies(dbN)小波,對主要影響指標進行小波多分辨率分析,為小波神經網絡預測模型的建立進行預處理[7,11-13]。
BP神經網絡能夠較好處理模糊的非線性關系,因而可以很好適用于水質預測等復雜問題的研究[14]。本研究進一步探討,基于小波分析的BP神經網絡模型在水質預測方面的適用性。所采用的是反向傳播人工神經網絡BPNN(Back-Propagation Artificial Neural Network),其典型結構為3層網絡模型,即輸入層-隱含層-輸出層,具體原理見圖2。
輸入層、隱含層與輸出層之間的神經元連接可描述為:
a1=f(IW·Pi+b1)
(8)
a2=f(LW·Pi+b2)
(9)
式中a1、a2——隱含層和輸出層的計算值;Pi——輸入量;IW、LW——隱含層連接輸入層和輸出層之間的權重;b1、b2——隱含層與輸入層和輸出層之間的偏差值。
其中,隱含層的傳遞函數采用雙曲正切S型函數:
(10)
輸出層的傳遞函數為純線性函數:
f(x)=x
(11)
本文通過小波分析對渦河流域水質指標進行多尺度變化規律分析,并結合BPNN神經網絡對其進行預測。首先,使用主成分分析方法確定關鍵水質參數,這作為水質預測的目標值,提高預測效率。再通過相關性分析方法將與目標參數相關性較強的水質參數篩選出來并作為預測目標值的輸入值,提高預測精度。最后通過Morlet小波分析渦河流域水質指標多尺度變化規律,為水質預測周期提供參考。
在模擬實測水質數據時,需要預處理工具去噪來提高模型的精度,而小波分解能將原始數據中不同頻段的信息進行分解,極大地降低數據復雜度,再分別對這些數據進行預測從而提高預測精度。因此,通過比較幾個Daubechies小波,選用最為合適的小波基和分解尺度對目標參數的時間序列數據進行降噪處理,得到低頻概貌分量(caN)和多個高頻細節分量(cdN)。
通過BPNN神經網絡,分別對多個分量進行預測,并將其預測結果相加即為目標值的最終預測結果,BPNN模型均采用3層網絡結構,主要采用自動完成的“試錯法”來訓練模型并獲取最佳的模擬效果,可調整的主要參數包括隱含層神經元數目、迭代次數、動量因子和學習速率等。具體的技術路線見圖3。

圖3 技術路線
本研究數據來源于渦河流域中下游的4個監測斷面2005—2018年(共168個月)的水質指標每月一次監測值,其中玄武閘斷面位于河南省鹿邑縣,大寺閘斷面、渦陽閘斷面和蒙城閘斷面位于安徽省內,具體斷面位置見圖1。應用的水質指標為溶解氧(DO)、高錳酸鹽指數(KMnO4)、氨氮(NH4-N)、化學需氧量(COD)和總磷(TP)共5項,其監測數據完整且連續。
為提高預測效率,對溶解氧、高錳酸鹽指數、氨氮、化學需氧量和總磷5個水質指標進行主成分分析選取預測參數。具體主成分分析結果見表1。由表1可知,第一主成分方差貢獻率為47.56%,其反映了原變量的大部分信息,其中化學需氧量、氨氮和高錳酸鹽指數得分較高,且化學需氧量得分最高;第一和第二主成分方差累計貢獻率為 69.71%,其中溶解氧得分最高為0.81;第三主成分中,總磷得分最高達到了1.07。

表1 主成分分析結果
目前影響渦河流域水質的主要因子為以化學需氧量為代表的,包括氨氮、高錳酸鹽指數在內的第一主成分因子。它反映出的信息一方面為有機物指標,如化學需氧量和高錳酸鹽指數,另一方面為水體營養鹽的指標,如氨氮。說明目前導致渦河流域污染的主要原因為有機污染,因此利用主要污染因子,即用化學需氧量的變化來反映渦河流域污染的變化是可行的。通過分析化學需氧量與各指標間的相關性,為化學需氧量的預測選擇合適的輸入參數,具體相關性分析結果見表2。由表2可知,化學需氧量與高錳酸鹽指數、氨氮的相關性較強,相關系數分別達到0.604、0.650,因此將其作為化學需氧量預測模型的輸入變量。

表2 各指標相關關系
基于主成分分析法選取了渦河流域水質的主要影響因子,即化學需氧量。通過對比分析選擇合適的分解尺度對化學需氧量原序列進行分解,具體對比結果見圖4。

a)小波基(s)
由圖4a—4e可以看出,隨著分解尺度N的增加,曲線變得更平滑,并趨近原序列,但同時會丟失許多細節成分,當N=5時,曲線太過平滑,無法反映原始序列的峰值。綜合比較表明,db2既能顯示化學需氧量原序列s的大致概貌,細節成分保留也較完整,而且能很好地反映峰值。由于國內外缺乏最佳小波基的規范性選取方法,但不同小波基的預測結果差別不大,因此,綜合比較選擇db2小波基對化學需氧量序列進行小波分解。
由圖4f—4i可以看出,低頻概貌部分體現了原始信號的基本形狀,可以看成是原始序列的一個近似。隨著分解尺度N的增大,分解序列變得越來越光滑,然而當N=3時,曲線太光滑,無法表現原序列的總體概貌。綜合比較當N取2時,最適合于化學需氧量時間序列的小波分析。此時化學需氧量原序列s分解如下:
s=ca2+cd2+cd1
(12)
式中 ca2——趨勢成分;cd2——周期成分;cd1——隨趨成分。
構建化學需氧量小波BPNN模型時,其中ca2分量為原序列的趨勢項,應考慮與化學需氧量變化有關的相關變量,因此選取與主要影響因子相關系數較大的指標,即高錳酸鹽指數和氨氮作為預測ca2的輸入變量,ca2為輸出變量;cd2為原序列的周期項,cd1為原序列的隨趨成分,與各水質參數相關性不大,基本只與原始序列自身的變化有關,因此選擇化學需氧量本身作為其輸入,cd2和cd1分別作為輸出,然后將所有分變量預測結果疊加作為化學需氧量最終預測結果。經過多次試驗測試與參數調整,在保證網絡快速收斂的同時又具有較高的預測精度,使模型達到最優結果。訓練擬合中db2-2相關參數的最終配置結果見表3。

表3 模型關鍵參數
基于Morlet小波,對渦河流域河流2005—2018年,共168個月的溶解氧、高錳酸鹽指數、化學需氧量、氨氮和總磷進行多尺度變化分析。通過小波實部等值線圖和小波方差來識別不同水質指標的主周期及變化特征,具體結果見圖5。

a)溶解氧等值線

e)化學需氧量等值線
由圖5小波系數等值線可以看出,不同水質指標序列存在以不同尺度隨時間交替變化的特性,圖5小波方差可看出,各指標的小波方差圖的主要峰值所對應的時間尺度與小波實部等值線圖所顯示的波動變化規律相吻合。從圖4可知:渦河流域溶解氧8、18、33、64個月尺度交替變化表現較清晰,波動極值點分布規律明顯,以18個月尺度較為顯著,說明渦河流域溶解氧系列存在以8、18、33、64個月變化的主要周期;高錳酸鹽指數存在以9、20、31、50個月變化主要周期,以20個月尺度較為顯著;化學需氧量存在以8、20、31、43個月尺度變化主要周期,以31個月尺度較為顯著;氨氮存在以8、20、31、57個月尺度變化主要周期,以57個月尺度較為顯著;總磷存在以8、12、23個月變化主要周期,以23個月尺度較為顯著。綜合來看,各水質指標具有多尺度振蕩的特點,且主要存在以8、20、31個月變化的主周期,各水質指標主周期變化主要發生在汛期。
為探究主周期形成的主要原因,基于各指標多年平均值來分析各指標年內變化趨勢,具體結果見圖6。由圖6可以看出,各指標在汛期(7—9月)含量明顯低于年內其他月份,這可能是由于降雨量的增加,河段水流量也隨之增加,會對污染物產生稀釋作用,各指標含量降低,因而各指標呈現出汛期周期變化規律。

圖6 各指標年內變化趨勢
基于上文的各個小波基和分解尺度對比分析結果,最終選擇db2小波基對化學需氧量原序列s進行2層尺度的多分辨率分析,得到低頻概貌分量ca2和高頻細節分量cd2及cd1,見圖7。

a)s
由圖7可知,化學需氧量1年中的變化情況大致為夏季高、冬季低,春秋相當,但年際情況不盡相同。ca2與原序列最為接近,體現了化學需氧量的總體下降趨勢,且在2005—2009年(第1—50月)波動較大,但具體細節變化未得到很好體現。cd2盡管為細節部分,但對原序列的一些峰值如第40個月仍有著較好的表達。cd1在峰值處呈現較大的波動,但更多地含有隨機成分。
考慮到化學需氧量時間序列的長度,以及汛期周期變化的規律,最終選擇以一整年12個月作為預測周期。即以2005—2017年(第1—156月)的月均化學需氧量數據作為訓練對象進行訓練,以2018年各月化學需氧量含量(第157—168月)作為最終預測對象。預測結果的誤差分析見表4,預測值與實測值的擬合結果見圖8。

表4 化學需氧量第157—168月預測結果誤差分析

圖8 化學需氧量預測值與實測值對照結果
由表4可知,在預測期內平均相對誤差控制在8.4%,均方根誤差為1.5。2018年1—6月(第157—162月)誤差較小,7—10月(第163—166月)誤差相對較大,7月份誤差最大,可能與降雨量大小有關,因此在提高水質預測精度時,可以考慮加入季節性因素。圖8可知,預測值與實測值的曲線變化趨勢一致,尤其在短時預測時段擬合更優,總體來看,構建的小波神經網絡模型預測化學需氧量精度較高,且較為穩定,能夠適用于渦河流域的短期水質預測。
基于渦河流域2005—2018年的水質指標月監測數據,即溶解氧、氨氮、高錳酸鹽指數和化學需氧量,探究小波神經網絡在流域水質污染中的應用。通過Morlet小波分析了渦河流域水質指標多尺度變化規律,并基于Daubechies小波和BPNN神經網絡建立了小波神經網絡預測模型,該研究結果證實了小波分析在水質時間序列變化規律及預測方面的可行性。
a)渦河流域水質指標的周期變化規律表明,總磷有3個主要變化周期,其他指標均具有4個變化主周期,各水質指標變化主周期存在一定差異,但主要存在以8、20、30個月左右變化的主周期,各水質指標主周期變化主要發生在汛期;目前渦河流域水質的主要影響因子是以化學需氧量為代表的,包括氨氮、高錳酸鹽指數在內的第一主成分因子。
b)選擇db2小波基對化學需氧量原序列進行了2層尺度的多分辨率分解;通過小波神經網絡得到的化學需氧量預測值與實測值的曲線擬合度較高,平均相對誤差為8.4%,均方根誤差為1.5,該方法較適用于渦河流域的短期水質預測,對于水質預測具有較好的參考價值。
c)由于國內外缺乏最佳小波基的規范性選取方法,且研究發現季節性因素對水質序列預測誤差影響較大,因而在今后提高小波BP神經網絡預測精度時可以考慮從小波基的優化選擇和季節性因素的輸入出發。