丘陵小流域氮磷流失通量監測頻次的優化研究*

蔡燕秋 唐翔宇 張建強 程 軍 唐家良 章熙峰

(1.西南交通大學地球科學與環境工程學院，四川 成都 611756；2.中國科學院、水利部成都山地災害與環境研究所，四川 成都 610041)

河流是流域內污染物遷移進入湖泊、水庫等水體的主要途徑。河流的水質指標和徑流量是污染物流失通量計算的重要依據。由于高頻次水質監測的成本較高，我國各地區水質監測頻率一般為1～4次/月，在污染物流失通量計算時需要通過線性插值法或非線性插值法進行估算。在人類活動強度高的區域內，水質與水文情況復雜，線性插值法嚴重偏離實際情況[1]。已有多種模型如SWAT、雨水管理模型(SWMM)等[2-6]被應用于較大流域污染物流失通量的估算，但此類模型并不適用于小流域。

河流污染物流失通量與徑流量間一般存在著較好的相關關系，美國地質調查局(USGS)據此開發了LOADEST模型。該模型能利用連續的徑流量數據和離散的水質數據，建立污染物流失通量與徑流量之間的非線性回歸方程，因其符合污染物波動式流失的常態，并且操作簡便、數據需求量小，已被廣泛應用[7-10]。宋方方[11]通過LOADEST模型對趙家溪污染物流失通量進行估算，為趙家溪污染控制方案的制定提供數據支撐；DWIVEDI等[12]研究發現，LOADEST模型用于估算低水平范圍的大腸桿菌負荷有較好的效果；JHA等[13]用LOADEST模型對美國紐斯河的氮磷流失通量進行估算，發現估算值與實測值相比總體偏高；PARK等[14]指出回歸方程不僅與水質參數有關，而且與暴雨事件樣本所占比例有關；高信娟[15]在九龍江的研究結果表明，水質數據缺失越少，污染物流失通量的估算值越接近實測值。

本研究以四川省鹽亭縣萬安丘陵小流域截流、大興及萬安3個河流斷面為對象，比較不同監測時間間隔條件下總氮(TN)、總磷(TP)日流失通量估算值相比實測值的相對誤差，以確定既經濟可行又結果可靠的最低監測頻次，并揭示小流域屬性的空間尺度差異所產生的影響。

1 材料與方法

1.1 研究區

研究區位于四川省鹽亭縣萬安小流域(105°27′E，31°16′N)，海拔為389～426 m，屬中亞熱帶季風氣候。2015年、2016年和2017年的年降雨量分別為1 029.3、886.7、628.7 mm，其中5—10月降雨量占全年總降雨量的比例分別為87%、83%和76%。該小流域內截流、大興、萬安集水面積分別為0.35、4.80、12.36 km2。研究區概況見圖1。

圖1 萬安小流域概貌與監測點位置Fig.1 Map of Wan’an catchment and monitoring locations

1.2 LOADEST模型

LOADEST模型參數的估計方法有極大似然估計(MLE)法、校正極大似然估計(AMLE)法和最小絕對偏差(LAD)法3種。一般采用AMLE法，但當殘差不符合正態假設時，則考慮采用LAD法。當源文件中將模型選擇設置為0時，LOADEST模型運用Akaike信息準則(AIC)[16-18]和Schwarz后驗概率準則(SPPC)[19]，自動從9個回歸方程(見式(1)至式(9))中取得最小AIC值和SPPC值的方程為最優方程，并進行自變量中心化消除多重共線性[20]。

lnL=a0+a1lnQ

(1)

lnL=a0+a1lnQ+a2lnQ2

(2)

lnL=a0+a1lnQ+a2td

(3)

lnL=a0+a1lnQ+a2sin(2πtd)+a3cos(2πtd)

(4)

lnL=a0+a1lnQ+a2lnQ2+a3td

(5)

lnL=a0+a1lnQ+a2lnQ2+a3sin(2πtd)+
a4cos(2πtd)

(6)

lnL=a0+a1lnQ+a2sin(2πtd)+a3cos(2πtd)+a4td

(7)

lnL=a0+a1lnQ+a2lnQ2+a3sin(2πtd)+
a4cos(2πtd)+a5td

(8)

(9)

式中：L為污染物流失通量；a0～a6為回歸系數；Q為徑流量；td為分數化的時間與中心化時間之差；參數單位均視具體情況而定。

LOADEST模型主要通過以下3種方式檢驗回歸方程的有效性：(1)判定系數(R2)越接近1，證明方程的擬合程度越好；(2)殘差序列相關系數(SCR)越接近0，證明方程的各殘差之間越獨立；(3)概率曲線相關系數(PPCC)越趨近于1，表明殘差越接近正態分布。

1.3 方法比較與優化

于2015—2017年，在上述3個河流斷面進行徑流量的連續監測，并對河水TN和TP兩種污染物指標開展采樣監測，頻次為雨季(5—10月)每周1次、非雨季(11月至次年4月)每10天1次。

采用LOADEST模型以不同時間間隔條件下的水質數據和連續的徑流量數據(日數據)分別對3個河流斷面的TN和TP日流失通量進行估算。為避免監測日期選擇的偶然性，將用于流失通量估算的污染物濃度數據起算日期錯開得到時間間隔相同的不同數據分組。采用式(10)計算2015—2016年較大與較小時間間隔條件下污染物總流失量估算值的相對誤差：

(10)

式中：RE為相對誤差，%；Li和L10為用LOADEST模型估算得到的時間間隔為i及10 d時某斷面的污染物兩年總流失量，kg；i為時間間隔，d。

2 結果與討論

2.1 污染物流失通量的時空變化特征

由圖2可見，在時間分布特征上，雨季各斷面的TN、TP日流失通量的平均值明顯高于非雨季，且在雨季的波動范圍比非雨季大。從空間尺度來看，截流斷面的TN、TP日流失通量的平均值明顯低于大興和萬安斷面；大興斷面的TN日流失通量略高于萬安斷面，其TP流失通量平均值與萬安斷面在非雨季基本持平，而在雨季略高于萬安斷面；截流斷面在雨季與非雨季的TN、TP日流失通量波動最大，動態變化特征最為突出。

2.2 回歸方程的有效性

表1給出了不同監測頻次下3個河流斷面的TN和TP流失通量回歸方程的R2。對于水質監測時間間隔相同的不同數據分組應選擇不同的最優回歸方程。可以看出，大部分方程的擬合程度均比較好，除大興斷面外，R2在0.75以上。SCR為0.01～0.57，表明部分模型擬合受到共線性影響，且受影響的模型為水質監測時間間隔為30～40 d的數據分組。PPCC均在0.9以上，表明各污染物流失通量回歸方程的參數估計均可采用AMLE法。

表1 截流、大興及萬安斷面污染物流失通量回歸方程的檢驗結果1)

2.3 不同監測頻次污染物總流失量估算值的相對誤差分析

為驗證水質監測時間間隔為10 d條件下小流域污染物流失通量LOADEST模型估算值的可靠性，將其與流失通量實測值(污染物濃度實測值與日徑流量實測值的乘積)的時間變化特征對比，結果表明：污染物日流失通量的估算值與實測值吻合較好(估算值相對誤差在-25%～25%的數據占49.59%)，動態變化特征基本一致(見圖3)。因此，可采用水質監測時間間隔為10 d的污染物流失通量估算值作為參比值，驗證較大的監測時間間隔條件下污染物流失通量估算值的相對誤差。

圖3 監測時間間隔為10 d條件下污染物流失通量LOADEST模型估算值與實測值的時間動態變化Fig.3 Dynamics of the estimated pollutant loss flux in the river using LOADEST with the measured pollutant loss flux at a time interval of 10 d

較大的監測時間間隔條件下污染物兩年總流失量估算值(大興斷面部分數據由于無法收斂而缺乏參數估算結果)的相對誤差見圖4。結果表明，對于同一個河流斷面，在校準文件中輸入不同時間間隔條件下的污染物濃度及日徑流量實測數據時，所得到的污染物總流失量估算值存在差異。相對誤差絕對值最小為0.07%，最大為148.08%，說明不同時間間隔條件下獲得的兩年總流失量估算值之間存在極大差異。從相同時間間隔的不同數據分組來看，在40 d間隔條件下，截流斷面TN總流失量估算值的相對誤差變化范圍很大，但萬安斷面TN總流失量估算值卻較為穩定。盡管截流斷面在不同監測時間間隔條件下獲得的TN流失通量最優回歸方程R2均在0.94以上，但其兩年總流失量的模型估算值之間卻存在較大差異，這表明影響截流斷面TN流失的相關因素及過程隨時間推移呈現隨機、非單調性變化。

圖4 不同監測時間間隔條件下兩年總流失量的LOADEST模型估算結果及相對誤差Fig.4 The estimated 2-year total exports of pollutants based on monitoring data obtained at different time intervals and their relative errors

從不同的監測頻次來看，時間間隔為20 d條件下，不同河流斷面的污染物總流失量估算值的平均相對誤差為10.70%；時間間隔為30 d條件下，平均相對誤差為15.02%；時間間隔為40 d條件下，平均相對誤差為31.76%。這說明水質監測時間間隔越長，污染物總流失量的估算值與實測值的偏差可能越大。

注：STN、STP分別為單位集水面積TN、TP日流失通量，g/(d·hm2)。圖2 2015—2016年河流斷面的非雨季與雨季單位集水面積實測污染物日流失通量Fig.2 Daily pollutant loss flux per unit drainage area in non-rainy and rainy season during 2015-2016

從不同污染物來看，TN的平均相對誤差為28.93%，TP的平均相對誤差為20.90%。這說明LOADEST模型對TP流失通量的估算更為可靠，但對TN流失通量估算的適用性相對較差，此結果與KIM等[21]得出的結果相似。

以25%作為污染物總流失量LOADEST模型估算的允許誤差，由表2可以看出，截流斷面TN和TP流失通量監測最大時間間隔為20 d；大興斷面TN監測最大時間間隔可設定為30 d，而TP流失通量則不能用LOADEST模型估算；萬安斷面TN和TP監測最大時間間隔為40 d。

表2 各河流斷面模型估算結果在25%允許誤差范圍內的占比

丘陵小流域水文過程一般存在顯著的空間尺度變化規律，溪溝河流的徑流系數隨集水面積增加先增大而后趨于恒定，降雨徑流量和污染物流失通量隨時間的變化也逐漸趨于平緩。本研究的數據表明，當集水面積過小，污染物流失通量的估算值相對于實測值的誤差不穩定，因此需要加大監測頻率才能增加模型輸出結果的可靠性。

2.4 模型外推可行性分析

利用2015—2016年污染物濃度數據與2015—2017年徑流量數據外推得出2017年氮磷流失通量，與用2017年水質指標和徑流量實測數據計算出的結果相比較，驗證外推的可行性。由表3可以看出，外推估算值與實測值差異較大，相對誤差最高可達182.70%。因此，外推估算法不可行，對于該丘陵小流域，當年徑流量與水質實測數據都是準確估算當年TN、TP總流失量的必要數據，氣象水文情況的年際變化會對估算結果產生不可忽視的影響。

表3 污染物年總流失量的LOADEST模型外推估算值和實測值

2.5 模型估算誤差的來源分析

2015—2016年各河流斷面污染物濃度的變異系數分析表明，大興斷面TN、TP濃度的變異系數均為最高，截流斷面次之，萬安斷面的變異系數最小(見表4)。濃度分布越離散，就需要越多的數據才能實現模型的良好擬合。因此，在相同監測頻次下LOADEST模型在大興斷面表現較差，而在萬安斷面則表現良好。

表4 河流各斷面污染物濃度的變異程度及集水區基本屬性

結合集水區的基本屬性，可辨識污染物濃度的變異原因。一般而言：集水區面積越大，緩沖性能越好，濃度變異系數就越小；坡降越大，越有利于污染物遷移，其濃度變異系數就越大；農地和居民地占比越大，施肥、人為排放的污染物越多，污染物濃度越易受到人類活動的影響。由表4可見，對于面積小、坡降大的集水區(截流斷面)污染物徑流流失過程(特別是河道徑流)漲退變化快、峰值多；而對于面積大、坡降小的集水區(萬安斷面)，河道徑流的漲退變化較慢、峰值少。因此，前者需要較高的監測頻次才能準確刻畫污染物流失通量隨時間的變化規律。大興斷面則有所不同，處于河谷平壩區，受居民地和農田兩方面的影響均較大，故而氮磷流失通量的變異系數較高。

3 結 論

為探究LOADEST模型在丘陵小流域的適用性，優化污染物流失通量估算方法，本研究以四川省鹽亭萬安丘陵小流域為對象，估算基流條件下的氮磷流失通量及其年際變化，比較不同監測時間間隔的TN、TP流失通量估算值與實測值的差異，并分析其成因，結果表明：對于集水面積大、坡降小的河流斷面，允許的最大監測時間間隔為20 d；而對于集水面積大、坡降小的河流斷面，允許的最大監測時間間隔為40 d；LOADEST模型在TP流失通量估算中的表現優于TN；外推估算法不適用于該小流域的污染物流失通量估算。