基于LOADEST模型和小波變換的河流氮磷污染動態(tài)分析

金亞楠，張柏發(fā)，郝 韻，鄔建紅，呂 軍

(浙江大學 環(huán)境與資源學院，浙江 杭州310058)

非點源污染致使水體富營養(yǎng)化，是自然水體水質退化的主要原因[1]。盡管流域氮、磷等營養(yǎng)物質的非點源污染過程十分復雜，但最終都會反映在河流的水質變化上。水中氮、磷污染物濃度變化的直接監(jiān)測和污染物通量估算，是進行河流水質評價和實施水污染控制的基礎工作。目前，我國大多數(shù)河流都已建立了連續(xù)的流量監(jiān)測系統(tǒng)，但水質的自動化監(jiān)測尚未普及。一些已有的在線自動監(jiān)測裝置，也大多存在著受環(huán)境變化影響較大、監(jiān)測閾值限制，以及安裝、使用和維護的成本較高等問題[2]。傳統(tǒng)的人工采樣和水質測定方法雖然仍在水質監(jiān)測中廣泛應用，但由于耗時費力，很難獲得連續(xù)的或高頻的水質監(jiān)測資料；離散的每月1～2次的水質采樣分析結果，卻又難以反映河流水質的連續(xù)變化。LOADEST模型是美國地質調查局(USGS)為此專門開發(fā)的，用于根據(jù)河流離散的水文水質監(jiān)測資料建立營養(yǎng)物負荷量與水流通量之間關系的多元統(tǒng)計模型[3-5]。應用該模型就可以在每日流量資料的基礎上估算氮、磷負荷量在河流中的連續(xù)變化。這一方法已被國內外廣泛采用[6-8]。

LOADEST模型是優(yōu)化的統(tǒng)計回歸模型，而統(tǒng)計方法的基本假設是樣本數(shù)據(jù)在時空上有獨立性。但實際上，驅動河流水質變化的主要因素(水文氣象因素和氮、磷的投、排放變化等)具有明顯的時間上的周期性和延續(xù)性，以及空間上的傳承性和加和性。小波變換是在序列數(shù)據(jù)的時間-頻率分析領域迅速發(fā)展的一種新技術[9]，它可以通過解構序列數(shù)據(jù)，由粗到細地描述各種復雜信號的具體結構，從而實現(xiàn)對整個數(shù)據(jù)序列的局部化特征和多層次變化規(guī)律的分析。本研究在對浙江東部典型農業(yè)流域長樂江出口斷面進行每月一次水質監(jiān)測的基礎上，分析其2003—2016年間河流水文水質的變化特征，建立LOADEST模型來核定河流總氮(TN)和總磷(TP)的每日負荷量，并采用小波分析方法，進一步探索河流水質和營養(yǎng)物負荷通量的動態(tài)變化，以期揭示我國東部地區(qū)農業(yè)流域河流非點源污染變化規(guī)律及其時間格局特征，為農業(yè)流域河流水質控制提供科學依據(jù)。

1 材料與方法

1.1 研究區(qū)概況

長樂江流域位于浙江省紹興市嵊州市(120°35′56″～120°49′03″E，29°27′98″～29°35′12″N)，是我國東部地區(qū)典型的低山河谷平原農業(yè)流域，流域總面積約為653 km2，多年平均降水量1 307.4 mm，但全年約75%的降水集中在4—9月，多年平均氣溫17.48 ℃。土地利用類型主要有林地(46.14%)、水田(21.04%)、園地(18.30%)、旱地(5.87%)、人居地(5.71%)、水域(2.31%)和其他(0.62%)。流域內的主干河流——長樂江是曹娥江上游，全長70.5 km，平均河寬55 m，河道比降為3.6‰，多年平均流量為16.39 m3·s-1，輸沙量和侵蝕模數(shù)分別為10.9萬t·a-1和127 t·km-2·a-1。

1.2 水質監(jiān)測方法與水文氣象數(shù)據(jù)來源

本研究主要分析長樂江流域出口斷面監(jiān)測點(圖1中S7)的氮、磷營養(yǎng)物變化。水樣采集和監(jiān)測分析時間為2003年8月—2016年12月，采樣監(jiān)測頻率為每月1次。在監(jiān)測斷面河流主流道距水面30 cm左右位置采集水樣，裝入聚乙烯瓶，用稀硫酸酸化后，當天運回實驗室就TN、TP濃度等水質指標進行分析測定。TN濃度參照GB 11894—89，采用堿性過硫酸鉀消化-紫外分光光度法測定；TP濃度參照GB 11893—89，采用過硫酸鉀消化-鉬銻抗比色法測定。

圖1 長樂江流域和水質監(jiān)測點的位置示意圖Fig.1 Location of Changle River watershed and water quality monitoring points

河流日流量數(shù)據(jù)由浙江省水文局提供；氣象數(shù)據(jù)來自于中國氣象科學數(shù)據(jù)共享服務網(wǎng)(http://www.escience.gov.cn/metdatapageindex.html)。

1.3 LOADEST模型

USGS專門開發(fā)并發(fā)布了基于河流流量變化與水質變化關系的河流污染物負荷量估算和優(yōu)化模型，即LOADEST模型和相應的軟件[3]。LOADEST模型的主要思路是河流營養(yǎng)物的負荷量是河流流量的函數(shù)，并在不同的條件下伴隨著不同的時間周期性變化。LOADEST模型包含11個河流營養(yǎng)物負荷量與流量關系的多元回歸方程，基于水質監(jiān)測得到的營養(yǎng)物負荷量與流量數(shù)據(jù)進行各方程的擬合，通過赤池信息量準則(AIC)[10-11]和Schwarz后驗概率準則(SPPC)[12]比較各方程的擬合效果，最后采用具有最小AIC值和SPPC值的方程作為最優(yōu)的描述河流污染物負荷量與流量關系的回歸方程。

LOADEST模型軟件提供了3種方程參數(shù)的擬合估計方法，即最小方差無偏估計(MVUE)[13]、漸進極大似然估計(AMLE)[14]和最小絕對偏差估計(LAD)[15]，可對非正態(tài)分布數(shù)據(jù)資料、刪失型數(shù)據(jù)資料(具有不完整的或異常的序列數(shù)據(jù)資料)和變量共線性問題進行規(guī)范、方便的處理。

1.4 小波分析方法

小波分析的基本思想是，一個由變化的序列數(shù)據(jù)組成的函數(shù)(母函數(shù))，可以被分解為多組不同變頻的子函數(shù)來表示或逼近，從而可以通過對子函數(shù)的分析來表征母函數(shù)變化的細節(jié)。比如河流水體中的營養(yǎng)物負荷量的變化(母函數(shù))，可以通過正交小波變換分解為不同變頻的子函數(shù)組合，從而通過對子函數(shù)變化的特殊紋理和突變信息的分析，揭示負荷量變化的細節(jié)和規(guī)律。因此，小波分析又被稱為數(shù)字信號的“數(shù)學放大鏡”。小波分析的基本原理是通過增加或減小伸縮尺度(a)來得到信號的低頻或高頻信息，然后分析信號的概貌或細節(jié)，從而實現(xiàn)對信號不同時間尺度和空間局部特征的分析。在實際研究中，最主要的就是要由小波變換方程得到小波系數(shù)，然后通過這些系數(shù)來分析時間序列的時頻變化特征[16]。通過小波變換分析得到的小波系數(shù)，可以描述不同時期的信號做出周期性變化(信號振蕩)所需要的時間跨度(時間尺度)。如果要更深入地了解各時期信號振蕩周期的整體性規(guī)律，就需要借助方差分析的方法[17]。將小波系數(shù)的平方差在b域上進行積分，就可得到小波方差。小波方差隨尺度a的變化過程，稱為小波方差圖。小波方差圖可用來描述信號在不同時間尺度擾動的相對強度，或者是不同振蕩強度所對應的主要時間跨度(尺度)，即信號變化(振蕩)的主次周期。

2 結果與分析

2.1 河流水文水質特征

2003—2016年間，長樂江出口斷面的多年平均流量為16.39 m3·s-1，河流年平均流量呈顯著增加趨勢。河流水質監(jiān)測結果表明，實測水體中的總氮和總磷濃度平均值分別為(3.74±0.91)mg·L-1和(0.136±0.086)mg·L-1，遠高于該河流總氮和總磷濃度的參照標準(即未受人為污染時的河流水質指標參考標準：總氮1.70 mg·L-1，總磷0.055 mg·L-1)[18]。從年內的分布來看，月平均流量的變異較大，流量最大的6月份的多年平均流量為枯水期(10月—次年2月)月平均流量的2～3倍(圖2-a)。

圖2-b是研究期內河流監(jiān)測斷面上TN濃度的月平均分布狀況，可以看出，枯水季節(jié)河流TN濃度具有高于豐水季節(jié)的趨勢，這與我國南方其他地區(qū)的研究結果一致[19]，意味著多雨季節(jié)徑流對河流的稀釋作用強于對土壤的沖刷和對養(yǎng)分的淋溶作用。但在圖2-c中，河流中TP濃度的月平均變化在枯水期和豐水期之間未呈現(xiàn)出顯著差異，推測是在本研究區(qū)良好的植被和以水稻為主的種植制度下[20]，多雨季節(jié)徑流對土壤沖刷造成的河流中磷的增加與稀釋作用的效果基本相當。

圖2 2003—2016年長樂江出口斷面流量(a)、TN濃度(b)、TP濃度(c)各月份實測平均值Fig.2 River flow (a)， TN concentration (b) and TP concentration (c) measured at outlet of Changle River in each month from 2003 to 2016

河流中的養(yǎng)分通量主要取決于2個因素——河流流量與養(yǎng)分濃度。但相關分析表明，河流中TN和TP的負荷量都只與河流的流量呈極顯著(P<0.01)相關，而與它們的濃度相關性很弱(圖3)。實測河流TN負荷量與河流流量呈極顯著線性相關(R2=0.732 9，P<0.01)，而與TN濃度相關性較弱(R2=0.247 5，P>0.05)；實測河流TP負荷量與河流流量也呈極顯著線性相關(R2=0.676 1，P<0.01)，而與TP濃度相關性較弱(R2=0.191 4，P>0.05)?？梢?，雖然河流養(yǎng)分的負荷量來自于河流流量與養(yǎng)分濃度的乘積，但在研究期間長樂江TN和TP的負荷量變化都主要取決于河流流量的變化，養(yǎng)分濃度變化的影響比較弱。本研究顯示，河流養(yǎng)分負荷量與流量之間具有顯著的相關關系，這與許多前人的研究結果一致[1，21-22]；而污染物負荷量與濃度的相關關系則取決于河流的狀態(tài)，通常流量較穩(wěn)定的河流，其河流污染物負荷量與濃度的相關性會高于流量漲落起伏較大的源頭地區(qū)的河流。

圖3 長樂江出口斷面實測TN、TP負荷量與河流流量，以及TN、TP濃度的相關分析Fig.3 Correlation of measured TN， TP load with river flow and TN， TP concentration at outlet section of Changle River

2.2 河流氮、磷負荷量估算

為了更好地對河流流量和養(yǎng)分濃度的每日連續(xù)性變化進行表征，采用LOADEST模型，先建立河流營養(yǎng)物負荷量與河流流量之間的關系，再利用每日的流量資料來估算河流中每日的營養(yǎng)物負荷量。

首先，把2003—2016年的水質監(jiān)測數(shù)據(jù)和水文資料分成奇數(shù)年和偶數(shù)年2組，基于奇數(shù)年各月份的水質監(jiān)測資料與監(jiān)測當日的流量資料，通過LOADEST模型來建立和優(yōu)化長樂江出口斷面上的TN、TP負荷量與河流流量之間的回歸方程。LOADEST模型運行結果表明，通過AIC和SPPC優(yōu)選，2003—2016年間TN、TP負荷量與出口斷面流量之間的最優(yōu)模型均為7參數(shù)模型(即LOADEST系列模型中的第9個方程)，模型各參數(shù)的優(yōu)化標定結果如表1所示。然后，用偶數(shù)年水質監(jiān)測當日的平均流量來模擬估算所有偶數(shù)年水質監(jiān)測當日的TN和TP負荷量，并將計算結果與實測負荷量進行比較，來驗證回歸方程模擬結果的準確性和可靠性(圖4)。結果表明，在參數(shù)標定期(奇數(shù)年)，TN擬合方程的決定系數(shù)R2=0.89，納氏模擬效率系數(shù)Nse=0.86；驗證期(偶數(shù)年)模擬驗證結果的決定系數(shù)R2=0.78，納氏模擬效率系數(shù)Nse=0.76，達到了“非常好”的水平(R2和Nse均大于0.75)[23]。TP的方程擬合和驗證效果與TN模型相比稍差，但也達到了“好”(R2和Nse均大于0.60)的水平：標定期(奇數(shù)年)模型擬合的決定系數(shù)R2=0.76，納氏模擬效率系數(shù)Nse=0.69；驗證期(偶數(shù)年)模型擬合的決定系數(shù)R2=0.67，納氏模擬效率系數(shù)Nse=0.66。

圖4 長樂江出口斷面TN、TP負荷量的LOADEST模型擬合效果和驗證結果Fig.4 Calibrations and validations of LOADEST model for TN and TP load in outlet section of Changle River

表1 長樂江2003—2016奇數(shù)年TN和TP負荷量的LOADEST模型率定結果Table 1 Calibration results for LOADEST model of TN and TP loads in Changle River using 2003-2016 odd years’ water quality mornitoring data

驗證結果表明，所優(yōu)化的LOADEST模型可以直接用于所有研究期間長樂江每日TN和TP負荷量的模擬。為了進一步提高模擬的精度和可靠性，采用這14 a的全部實測資料，重新擬合優(yōu)化LOADEST模型，所得到的TN擬合模型的R2=0.82，Nse=0.80；TP擬合模型的R2=0.76，Nse=0.67。模型具體形式如下：

lnLTN=7.483 4+1.141 8lnQ-0.004 1lnQ2+0.035 1sin(2πT)+0.221 4cos(2πT)-0.004 1T+0.003 4T2；

(1)

lnLTP=4.300 2+1.053 2lnQ-0.048 4lnQ2-0.048 5sin(2πT)+0.056 8cos(2πT)+0.038 9T-0.005 1T2。

(2)

式(1)、(2)中，LTN和LTP分別為通過監(jiān)測斷面的TN和TP的日負荷量(kg·d-1)，Q為日平均流量(m3·s-1)，T為分數(shù)形式的日期(如2016年2月5日記為36/366)?；?1)式和(2)式，應用2003—2016年的每日流量資料，估算長樂江出口斷面研究期間的每日TN和TP負荷量，進而獲得各個時間尺度(月、年等)的TN、TP通量(圖5)。結果表明，TN的年負荷量為989.03～3 578.46 t·a-1，平均年負荷量為(2 207.71±862.48)t·a-1。山區(qū)河谷地區(qū)河流流量(特別是洪峰流量)變化巨大，導致TN負荷量的日平均值變化很大，為(6.04±11.81)t·d-1。TP的年負荷量為25.80～122.21 t·a-1，平均年負荷量為(71.43±31.56)t·d-1，平均日負荷量為(0.19±0.26)t·d-1。長樂江TN和TP的負荷量在2003—2016年間總體呈上升趨勢。從年內變化來看，多年平均的TN月負荷量在12月份最低(1 395.94 t)，而最高值出現(xiàn)在6月份(4 316.65 t)，與流量最大值的出現(xiàn)時間一致；多年平均的TP月負荷量最低值出現(xiàn)在1月份(49.18 t)，最高值與TN一樣出現(xiàn)在6月份，為143.53 t，變化趨勢與在其他河流上的研究結果類似[21-22]。需要注意的是，在特別高的流量條件下(如洪峰期間)，TN、TP負荷量估算可能存在低估的現(xiàn)象[1]。

圖5 2003—2016年長樂江出口斷面TN和TP年負荷總量和平均月負荷量分布Fig.5 Annual loadsand mean monthly loads of TN and TP at toutlet section of Changle River from 2003 to 2016

2.3 河流氮磷負荷量變化小波分析

采用Morlet小波對2003—2016年長樂江出口斷面的月TN和TP負荷進行連續(xù)小波變換，獲得它們在研究期間的小波系數(shù)圖譜(圖6)。圖6-a展示了長樂江出口斷面以月為單位的TN負荷量小波系數(shù)變化(振蕩)，大約以12～24個月的時間跨度(時間尺度)為周期，發(fā)生明顯的高低循環(huán)變化；特別是自2009年起變化加劇，反映了TN負荷量振蕩強度加劇，振蕩強度的正負值差加大。此外，在時間跨度為3～12個月的尺度上，也有一些較明顯的周期性變化。TP負荷量小波系數(shù)的變化周期與TN負荷量的變化類似，但變化強度相對較小，并且變化強度在研究期間一直呈增強趨勢，最大振蕩的周期出現(xiàn)在2014年以后(圖6-b)。

分析研究期間監(jiān)測斷面上的河流流量(圖6-c)、TN濃度(圖6-d)和TP濃度(圖6-e)的小波系數(shù)變化頻譜，可以發(fā)現(xiàn)，河流流量的小波系數(shù)變化與TN負荷量小波系數(shù)的變化最為相似，包括高低值振蕩的時間跨度(時間尺度)和振蕩發(fā)生的日期等；TN濃度小波系數(shù)的變化在時間尺度上與TN負荷量的變化類似，但變化強度大約僅為河流流量變化的1/200。TP濃度小波系數(shù)的變化，在振蕩的主要時間尺度和振蕩強度上均與TP負荷量小波系數(shù)的變化圖譜相差較大。因此，雖然TN和TP的負荷量是河流流量與TN濃度的乘積，但小波系數(shù)頻譜分析表明，由于流量的變化遠大于TN和TP濃度的變化，所以長樂江TN和TP負荷量的變化主要取決于流量的變化，而TN和TP濃度變化的影響相對較小。由于河流非點源污染過程主要受水文循環(huán)驅動，因此盡管氣候和物候條件并不相同，不同地區(qū)河流氮、磷負荷量的年度變化周期性卻既有類似的規(guī)律，也有區(qū)域性的特征。例如在松花江哈爾濱段上的類似分析表明，污染物負荷量在1 a(12個月)和2 a(24個月)時間尺度上的2個變化主周期基本與本研究一致。這是因為，無論是在浙江還是在黑龍江，河流污染物負荷量變化都與地表徑流過程的變化基本一致；但同時它們也有各自其他時間變化尺度上的特點[24]。

黃色和綠色表示接近小波系數(shù)的平均值，紅色和藍色代表小波系數(shù)偏高和偏低的狀態(tài)，線條為偏高(正值)或偏低(負值)相對量的等值線。Yellow and green represented data close to the average of wavelet coefficients. Red and blue represented high and low wavelet coefficients. The line represented the isoline of relative quantity of high (positive) or low (negative).圖6 長樂江出口斷面2003—2016年TN(a)和TP(b)負荷量、河流流量(c)、TN(d)和TP(e)濃度小波系數(shù)頻譜等值線圖Fig.6 Contour of wavelet coefficient for TN (a) and TP (b) loads， river flow (c)， TN (d) and TP (e) concentrations at outlet section of Changle River from 2003 to 2016

利用小波方差可以更明確地定量判別不同時間尺度下所發(fā)生的TN和TP負荷量振蕩的相對強弱和各種振蕩發(fā)生的時間周期。如圖7所示，TN和TP的小波方差的峰值分別出現(xiàn)在6、12和18個月的時間尺度上(圖7-a和7-b)，意味著它們負荷量從低值到高值的輪回變化平均分別有6、12和18個月3個周期，其中18個月時間尺度上的小波方差最大，是變化最明顯的主周期。TN和TP負荷的這種變化規(guī)律首先是由河流流量的變化周期決定的(圖7-c)，而TN濃度(主周期為19個月)和TP濃度(24個月)的變化影響較小，因為它們的變化強度遠遠小于流量的變化強度(圖6)。

圖7 長樂江出口斷面2003—2016年TN(a)和TP(b)負荷量、流量(c)、TN(d)和TP(e)濃度的小波方差Fig.7 Wavelet variance of TN (a) and TP (b) loads， flow (c)， TN (d) and TP (E) concentrations at outlet section of Changle River from 2003 to 2016

在一些大型流域，地下徑流的水齡可長達數(shù)十年[25]，而通過地下徑流產(chǎn)生的河流非點源污染也具有數(shù)月至數(shù)年的滯后作用[26]。因此，長樂江TN、TP負荷量的變化，除了受河川徑流的影響而形成以18個月為主周期的重要變化規(guī)律以外，還可能與當?shù)剞r業(yè)種植制度、TN和TP的投/排放，以及上游水庫的蓄/放水操作等有關，形成了其他次要的多尺度周期性變化。這些多尺度的變化機制更為復雜，尚待進一步研究。

3 結論

(1)2003—2016年間長樂江出口斷面水質監(jiān)測表明，該河流水質明顯受到氮、磷的污染，水體中的TN和TP濃度平均值分別為(3.74±0.91)mg·L-1和(0.136±0.086)mg·L-1，均高于該河流TN和TP濃度的參照標準。河流中TN和TP的實測負荷量都與河流的流量呈極顯著(P<0.01)相關，決定系數(shù)R2分別為0.732 9和0.676 1；但與TN和TP濃度的相關性較弱，決定系數(shù)R2分別為0.247 5和0.191 4。

(2)低頻率的人工采樣水質監(jiān)測不能反映河流水體中負荷量的每日連續(xù)變化。利用LOADEST模型和2003—2016年間每月1次的水質監(jiān)測資料，建立了長樂江出口斷面上TN和TP的負荷量估算模型。模型驗證結果較好，其中TN負荷量的模擬結果與實測值對比的決定系數(shù)和納氏模擬效率系數(shù)分別為0.78和0.76；TP負荷量的模擬結果與實測值對比的決定系數(shù)和納氏模擬效率系數(shù)分別為0.67和0.66。

(3)利用LOADEST模型核定了長樂江TN和TP的每日負荷量，進而獲得各個不同時間尺度上TN和TP的負荷量。結果表明，在2003—2016年間，長樂江出口斷面的TN平均年負荷量為2 207.71 t·a-1，TP平均年負荷量為71.43 t·a-1；它們在研究期間均總體呈上升趨勢。多年平均的TN月負荷量在12月份最低(1 395.94 t)，在6月份最高(4 316.65 t)。

(4)2003—2016年以月為單位的TN和TP負荷量小波系數(shù)分析結果表明，長樂江TN負荷量以18個月的時間跨度為主周期發(fā)生著周期性的變化；TP負荷量的變化主周期與TN負荷量的變化類似，但變化強度相對較小。TN和TP負荷量變化的小波方差分析還發(fā)現(xiàn)，長樂江TN和TP負荷量的變化主要取決于河流流量的變化，而受TN和TP濃度變化的影響相對較小，因為TN和TP濃度的變化強度遠小于流量的變化強度。