丘陵小流域氮磷流失通量監(jiān)測(cè)頻次的優(yōu)化研究*

蔡燕秋 唐翔宇 張建強(qiáng) 程 軍 唐家良 章熙峰

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

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

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

本研究以四川省鹽亭縣萬安丘陵小流域截流、大興及萬安3個(gè)河流斷面為對(duì)象，比較不同監(jiān)測(cè)時(shí)間間隔條件下總氮(TN)、總磷(TP)日流失通量估算值相比實(shí)測(cè)值的相對(duì)誤差，以確定既經(jīng)濟(jì)可行又結(jié)果可靠的最低監(jiān)測(cè)頻次，并揭示小流域?qū)傩缘目臻g尺度差異所產(chǎn)生的影響。

1 材料與方法

1.1 研究區(qū)

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

圖1 萬安小流域概貌與監(jiān)測(cè)點(diǎn)位置Fig.1 Map of Wan’an catchment and monitoring locations

1.2 LOADEST模型

LOADEST模型參數(shù)的估計(jì)方法有極大似然估計(jì)(MLE)法、校正極大似然估計(jì)(AMLE)法和最小絕對(duì)偏差(LAD)法3種。一般采用AMLE法，但當(dāng)殘差不符合正態(tài)假設(shè)時(shí)，則考慮采用LAD法。當(dāng)源文件中將模型選擇設(shè)置為0時(shí)，LOADEST模型運(yùn)用Akaike信息準(zhǔn)則(AIC)[16-18]和Schwarz后驗(yàn)概率準(zhǔn)則(SPPC)[19]，自動(dòng)從9個(gè)回歸方程(見式(1)至式(9))中取得最小AIC值和SPPC值的方程為最優(yōu)方程，并進(jìn)行自變量中心化消除多重共線性[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為回歸系數(shù)；Q為徑流量；td為分?jǐn)?shù)化的時(shí)間與中心化時(shí)間之差；參數(shù)單位均視具體情況而定。

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

1.3 方法比較與優(yōu)化

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

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

(10)

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

2 結(jié)果與討論

2.1 污染物流失通量的時(shí)空變化特征

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

2.2 回歸方程的有效性

表1給出了不同監(jiān)測(cè)頻次下3個(gè)河流斷面的TN和TP流失通量回歸方程的R2。對(duì)于水質(zhì)監(jiān)測(cè)時(shí)間間隔相同的不同數(shù)據(jù)分組應(yīng)選擇不同的最優(yōu)回歸方程?？梢钥闯?，大部分方程的擬合程度均比較好，除大興斷面外，R2在0.75以上。SCR為0.01～0.57，表明部分模型擬合受到共線性影響，且受影響的模型為水質(zhì)監(jiān)測(cè)時(shí)間間隔為30～40 d的數(shù)據(jù)分組。PPCC均在0.9以上，表明各污染物流失通量回歸方程的參數(shù)估計(jì)均可采用AMLE法。

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

2.3 不同監(jiān)測(cè)頻次污染物總流失量估算值的相對(duì)誤差分析

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

圖3 監(jiān)測(cè)時(shí)間間隔為10 d條件下污染物流失通量LOADEST模型估算值與實(shí)測(cè)值的時(shí)間動(dòng)態(tài)變化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

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

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

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

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

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

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

表2 各河流斷面模型估算結(jié)果在25%允許誤差范圍內(nèi)的占比

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

2.4 模型外推可行性分析

利用2015—2016年污染物濃度數(shù)據(jù)與2015—2017年徑流量數(shù)據(jù)外推得出2017年氮磷流失通量，與用2017年水質(zhì)指標(biāo)和徑流量實(shí)測(cè)數(shù)據(jù)計(jì)算出的結(jié)果相比較，驗(yàn)證外推的可行性。由表3可以看出，外推估算值與實(shí)測(cè)值差異較大，相對(duì)誤差最高可達(dá)182.70%。因此，外推估算法不可行，對(duì)于該丘陵小流域，當(dāng)年徑流量與水質(zhì)實(shí)測(cè)數(shù)據(jù)都是準(zhǔn)確估算當(dāng)年TN、TP總流失量的必要數(shù)據(jù)，氣象水文情況的年際變化會(huì)對(duì)估算結(jié)果產(chǎn)生不可忽視的影響。

表3 污染物年總流失量的LOADEST模型外推估算值和實(shí)測(cè)值

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

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

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

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

3 結(jié) 論

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