施文飛， 陳俊鴻， 趙振鑫， 翟 婷

(桂林理工大學 環(huán)境科學與工程學院， 廣西環(huán)境污染控制理論與技術重點實驗室， 巖溶地區(qū)水污染控制與用水安全保障協(xié)同創(chuàng)新中心， 廣西 桂林 541004)

水污染是中國面臨的主要環(huán)境問題。人口持續(xù)增加和經(jīng)濟快速發(fā)展帶來的水資源短缺，生活污水、工業(yè)廢水的大量排放以及農(nóng)業(yè)面源污染等問題，導致地表水和地下水污染嚴重、生態(tài)環(huán)境惡化，已成為制約中國經(jīng)濟社會持續(xù)發(fā)展的重要因素。因此，探究河流污染負荷對于流域水環(huán)境承載力精準解析具有重要意義[1-3]。漓江是中國著名的風景名勝區(qū)，漓江風景名勝區(qū)的主體部分，其獨特秀麗的山水享譽海內外，有“世界十大水上奇跡”的美譽，同時肩負著流域內工農(nóng)業(yè)用水、居民生活生產(chǎn)用水、生態(tài)用水等重任。漓江整體水質狀況良好，其干流水質常年達到國家地表水Ⅱ類水質標準，但依然存在區(qū)間入河污染，其中以生活污水的排放為主，主要污染物為磷和氨氮。且由于枯水期水量減少，水體的自凈能力下降，漓江流域水環(huán)境污染事件時有發(fā)生。近年來不少學者針對漓江流域點面源污染及重金屬污染展開了一系列研究[4-5]，但漓江流域主要污染物通量估算及水質分析等鮮有研究。

污染物通量估算是掌握研究斷面河流污染負荷總量，表征流域污染物產(chǎn)輸特征的方法，但通常需要高頻率的水量水質監(jiān)測數(shù)據(jù)支撐。目前，國內河流流量監(jiān)測系統(tǒng)比較完善，但常規(guī)水質監(jiān)測受環(huán)境變化、監(jiān)測成本等限制，通常僅能獲得每月1次的離散數(shù)據(jù)，難以反映出河流水質的連續(xù)變化[6-7]，且無法反映出河流污染物通量逐日變化情況。LOADEST模型是由美國地質勘測局(United States Geological Survey，USGS)研發(fā)的一款用于估算河流污染物通量的模型。在實際應用中，利用離散、有限的水質數(shù)據(jù)及河流連續(xù)日徑流量數(shù)據(jù)，建立對應污染物通量與流量相關的回歸方程，從而估算河流在不同時空尺度下的污染通量。目前該模型已應用于美國北卡羅來納州紐斯河[8]、美國印第安納州小巴克溪[9]、美國印第安納州、威斯康星州和密歇根州流域[10]、寶象河[11]、淮河[12]、深圳河[13]等國內外多個流域，是估算低頻率水質監(jiān)測河流污染物通量的有效途徑之一[14]。然而，受數(shù)據(jù)觀測誤差的影響，模型擬合效果不好，其估算結果與實測值之間存在一定誤差。

因此，本文在已有模型研究的基礎上，提出基于“預處理-回歸”過程鏈的河流污染物通量估算方法。首先采用箱線圖法對水質數(shù)據(jù)進行預處理；隨后根據(jù)桂林及陽朔斷面逐日流量數(shù)據(jù)與逐月水質數(shù)據(jù)，基于LOADEST模型建立兩水文站污染物通量回歸方程，在驗證方程有效性的基礎上估算NH3-N(氨氮)、CODMn(高錳酸鹽指數(shù))及TP(總磷)通量，為漓江流域污染物治理提供科學依據(jù)。

1 材料與方法

1.1 研究區(qū)概況

漓江位于廣西東北部(109°45′～110°40′ E，24°18′～25°41′ N)[15]，是珠江水系支流桂江上游河段的通稱，起點為桂江源頭越城嶺貓兒山，終點為平樂三江口，從北向南依次流經(jīng)興安縣、靈川縣、桂林市區(qū)、陽朔縣及平樂縣，全長214 km，流域總面積5 857.11 km2。受季風氣候影響，年內降水分配極不均勻，具有典型的雨季和旱季交替特征，雨季降水占全年降水的72.26%，而1月月平均流量不足20 m3/s，由此導致豐水期與枯水期水量差異明顯。漓江干流水環(huán)境總體狀況良好，全年期滿足Ⅰ～Ⅲ類水質河長超過總長的99%，但支流水質仍有待改善。城區(qū)生活污水及城郊工業(yè)廢水、農(nóng)業(yè)灌溉用水等排入漓江及其支流，以及城區(qū)汽車、游輪等排放的含氮尾氣溶于水中，導致氮磷污染時有發(fā)生。桂林至陽朔段是漓江通航主干道，也是漓江風景最優(yōu)美的河段，保證漓江“黃金河道”的水質達標顯得格外重要。因此，本文以漓江桂林-陽朔河段為研究區(qū)域(圖1)，收集了2010—2019年漓江流域桂林及陽朔水文站逐月的水質資料以及逐日的流量資料，數(shù)據(jù)由桂林市水文中心提供。

圖1 研究區(qū)域示意圖

1.2 箱線圖法異常檢測機制

箱線圖[16]是一種用作顯示一組數(shù)據(jù)分散情況資料的統(tǒng)計圖，其常用于數(shù)據(jù)異常值處理。圖2展示了箱線圖位置點與標準正態(tài)分布之間的對應關系，Q1為1/4分位數(shù)點，Q3為3/4分位數(shù)點，Median為中位數(shù)，IQR為Q3與Q1值之差，可以看出，服從正態(tài)分布的樣本數(shù)據(jù)位于箱形范圍以外的概率僅為0.7%，因此可認為將處于Q1-1.5IQR及Q3+1.5IQR范圍外的點視為異常值[17]。在進行分析數(shù)據(jù)時，利用箱線圖可以判斷數(shù)據(jù)集的偏向及離散程度(觀察盒子的長度及上下隔間的形狀)；同時能直觀地識別數(shù)據(jù)集中的異常值(查看離群點)。

圖2 箱線圖異常值判定

1.3 LOADEST模型原理

污染物通量一般指在一定時段內污染物通過研究斷面的總量，可表示為

(1)

式中：Lτ為τ時段內斷面污染物通量，kg/d；Q為流量，m3/s;C為污染物濃度隨時間變化的函數(shù),mg/L。但在實際研究中，由于水質數(shù)據(jù)監(jiān)測頻率遠低于流量數(shù)據(jù)，無法直接使用式(1)計算河流斷面污染物通量，所以考慮時段平均流量與瞬時污染物濃度之積作為估算河流污染物通量的方法。在此基礎上，LOADEST模型利用多元線性回歸的方法對河流斷面污染物通量進行估算[18]。其通量估算公式通常表達為

(2)

(3)

表1 LOADEST模型提供的河流污染物通量回歸方程

LOADEST模型提供了3種方程參數(shù)的擬合估算方法，分別是最大似然估計法(MLE)[21]、漸進極大似然估計法(AMLE)[22]、最小絕對偏差法(LAD)[23]。通量估算中使用的主要估算方法是AMLE法，當數(shù)據(jù)(流量、水質)包含刪失數(shù)據(jù)時，則選用ALME法，當數(shù)據(jù)為非刪失型數(shù)據(jù)時，則采用MLE法。AMLE法和MLE法應用取決于模型殘差服從正態(tài)分布的假設，如果殘差不符合正態(tài)分布假設，則使用AMLE和MLE法估算結果可能不是最優(yōu)的，應考慮用LAD估算法來代替AMLE法。LAD通量估算法不依賴正態(tài)性假設。

優(yōu)選出的回歸方程利用殘差序列相關系數(shù)(SCR)檢驗殘差是否存在序列相關性，SCR值越小，表明殘差之間相互獨立。并通過相關系數(shù)R2來判定方程總體的擬合效果，R2越接近于1，表明擬合效果越好。對于非刪失型數(shù)據(jù)，適用概率曲線相關系數(shù)(PPCC)檢驗其殘差正態(tài)分布，相關系數(shù)越接近1，殘差越接近正態(tài)分布。

2 結果與分析

2.1 河流水文水質特征

由圖3(a)、(b)可知，2010—2019年，桂林及陽朔水文站斷面年平均流量分別為236.22、139.82 m3/s,河流年平均流量整體呈現(xiàn)階段性高低交錯的狀態(tài)，年內月平均流量變化均呈現(xiàn)“峰狀”。陽朔水文站月平均流量高于桂林水文站，且流量峰值相差大，每年豐水期流量峰值近于桂林站的2倍。由于桂林雨季長，降水充足，地下水資源豐富，河流自上游到下游，河流支流補給，導致陽朔站平均流量高于桂林水文站。

圖3 桂林及陽朔水文站斷面年平均流量、多年月平均流量、NH3-N濃度、CODMn濃度、TP濃度

圖3(c)、(d)、(e)表征研究期間桂林及陽朔水文站斷面NH3-N、CODMn、TP月平均濃度變化情況。兩水文站NH3-N、CODMn、TP年平均濃度分別為0.26(桂林)、0.12(陽朔)、1.63(桂林)、1.56(陽朔)、0.056(桂林)、0.056(陽朔) mg/L?？傮w水質情況較好，均能達到Ⅲ類水以上標準。桂林水文站NH3-N含量普遍高于陽朔站，漓江流經(jīng)城區(qū)時，城市生活污水排放城郊農(nóng)業(yè)面源污染、汽車及游船含氮尾氣排放溶于水形成氨氮，均導致桂林水文站NH3-N含量的升高[3]。由圖3(d)、(e)可以看出，兩水文站TP、CODMn含量在枯水期高于豐水期，說明雨季徑流對河流的稀釋作用強于對土壤的沖刷及養(yǎng)分的淋溶作用，這與中國其他南方地區(qū)的研究結果一致[4]。

河流中的污染物通量主要取決于流量及污染物濃度兩個因素。以陽朔水文站為例進行相關分析，分析結果表明：河流中污染物通量與污染物濃度相關性弱，而與河流流量呈極顯著相關(P<0.01)。由圖4可知，河流NH3-N、CODMn、TP通量與河流流量的相關系數(shù)R2分別為0.648 9、0.763 1、0.820 5；與NH3-N、CODMn、TP濃度的相關系數(shù)R2分別為0.316 7、0.360 8、0.296 7?？梢钥闯觯廴疚锿恐饕Q于流量的變化，與污染物濃度變化的影響弱，這與其他學者的研究結果一致[7,11]。

圖4 陽朔水文站斷面實測通量與河流流量及污染物濃度的相關分析

2.2 水質數(shù)據(jù)預處理

由于人工水質監(jiān)測過程中的人為誤差，導致水質數(shù)據(jù)存在異常值，對河流污染物通量估算精度產(chǎn)生了重要影響[24]。因此，本文在進行污染物通量估算之前，首先利用箱線圖法對水質數(shù)據(jù)進行異常值預處理，結果如圖5所示。將桂林及陽朔水文站2010—2019年CODMn、TP、NH3-N污染物濃度按月份進行歸類，兩水文站3種污染物分別有120個濃度數(shù)據(jù)，利用箱線圖法對水質異常數(shù)據(jù)進行有效范圍判斷，得到桂林水文站CODMn、TP、NH3-N污染物濃度分別有異常值7、0、4個，陽朔水文站分別有6、4、9個，其中數(shù)據(jù)偏大的異常值占70%。常用的異常值處理方法有：①用平均值修正；②將異常值直接刪除；③將異常值視為刪失值處理。分別用這3種方法對異常值進行修正，以減小異常值對污染物通量估算的影響，比較預處理與未處理擬合效果，取最優(yōu)修正方法用LOADEST模型進行污染物通量估算。

1～12為桂林站逐月水質數(shù)據(jù)；13～24為陽朔站逐月水質數(shù)據(jù)。圖5 水質數(shù)據(jù)箱線圖

利用箱線圖法對水質數(shù)據(jù)進行預處理，將異常值分別用①平均值修正；②將異常值直接刪除；③將異常值視為刪失值處理進行修正，然后進行LOADEST模型模擬。未處理和3種預處理結果對比見表2，最優(yōu)預處理方法為①平均值修正。將未處理水質數(shù)據(jù)進行LOADEST模型模擬，得出桂林站的NH3-N、CODMn、TP相關系數(shù)R2依次為0.696 6、0.925 5、0.875 6；陽朔站的NH3-N、CODMn、TP相關系數(shù)R2依次為0.733 1、0.916 4、0.871 1。將模擬結果與未處理計算結果相比較，桂林水文站NH3-N、CODMn、TP擬合效果分別提高了2.37%、2.82%、0%；陽朔水文站NH3-N、CODMn、TP擬合效果分別提高了2.89%、1.72%、3.75%。結果表明，在已有模型研究的基礎上，提出基于“預處理-回歸”過程鏈的河流污染物通量估算方法是行之有效的。

2.3 漓江流域污染物通量估算

利用離散、有限的水質數(shù)據(jù)及河流連續(xù)日徑流量數(shù)據(jù)，通過LOADEST模型，建立污染物通量與日均流量的回歸方程，根據(jù)河流流量數(shù)據(jù)估算出污染物逐日通量。首先，將2010—2019年的水文水質監(jiān)測數(shù)據(jù)分成兩組，2010—2017年作為標定數(shù)據(jù)，2018—2019年作為驗證數(shù)據(jù)?；?010—2017年逐日流量數(shù)據(jù)與逐月水質數(shù)據(jù)，通過LOADEST模型來建立桂林及陽朔水文站斷面NH3-N、TP及CODMn污染通量與河流流量的回歸方程?；贚OADEST模型的桂林及陽朔水文站污染物通量回歸方程參數(shù)率定結果見表3、表4。然后，利用2018—2019年桂林及陽朔水文站的實測逐月水質及流量數(shù)據(jù)，將實測通量值與估測值進行對比，檢驗回歸方程的可靠性和準確性(圖6)。驗證結果表明，在標定期，桂林及陽朔水文站CODMn(0.92，0.91)及TP(0.88，0.86)擬合方程的相關系數(shù)R2與驗證期模擬驗證結果(0.86，0.70；0.86，0.83)的相關系數(shù)R2均大于0.70，表明回歸方程擬合程度非常好。相較于CODMn與TP，NH3-N的擬合效果稍差，但也達到了較好水平(R2>0.60)，桂林及陽朔水文站標定期擬合方程相關系數(shù)R2分別為0.64，0.75；驗證期相關系數(shù)R2分別為0.67，0.60。參數(shù)估計與驗證結果表明，利用LOADEST模型優(yōu)選出的污染物通量回歸方程可以用于估算兩水文站NH3-N、CODMn及TP通量，其結果合理、可信。利用2010—2019年桂林及陽朔水文站實測預處理(用平均值修正異常值)水質及流量數(shù)據(jù)，基于LOADEST模型優(yōu)選出斷面NH3-N、CODMn及TP通量回歸方程見表5。

表2 桂林及陽朔水文站未處理與預處理結果對比

表3 桂林水文站污染物通量回歸方程參數(shù)率定結果

表4 陽朔水文站污染物通量回歸方程參數(shù)率定結果

圖6 基于LOADEST模型的桂林及陽朔水文站污染物通量擬合效果及驗證結果

表5 桂林及陽朔站污染物通量回歸方程

桂林及陽朔水文站的參數(shù)檢驗結果見表6，結果表明：兩水文站斷面CODMn及TP通量擬合效果非常好(R2>0.90)，NH3-N擬合效果良好(R2>0.70)，且殘差不存在序列相關性、服從正態(tài)分布(SCR值接近于0，PPCC值接近于1)。

表6 桂林及陽朔水文站參數(shù)檢驗結果

用LOADEST模型基于上述回歸方程估算出兩水文站的污染物通量，結果如圖7所示。由圖7可知：桂林及陽朔水文站斷面NH3-N、CODMn、TP年平均通量分別為987.302、1 254.781 t，6 784.407、14 964.549 t，265.945、756.055 t。從年際變化來看，兩水文站污染物通量變化趨勢大體一致，陽朔站斷面CODMn及TP通量均高于桂林站斷面，且總量是桂林站斷面的2.21和2.84倍。NH3-N總通量兩水文站差距較小，在枯水年2011及2018年，桂林站斷面NH3-N總通量高于陽朔站斷面。相較于陽朔站，桂林站斷面污染物通量波動幅度較小。從年內變化來看，污染物通量逐月變化趨勢呈現(xiàn)單峰狀，兩水文站NH3-N通量“峰型”一致，豐水期降雨量大，約占全年降雨量的72.5%，污染物通過降雨-徑流沖刷形式進入河流。由圖7(d)、圖7(f)可知，陽朔站斷面CODMn及TP通量單峰遠高于桂林站斷面，各污染物通量最大值出現(xiàn)在5、6月份，最小值一般出現(xiàn)在12—2月份，這與河流豐水、枯水期變化一致。

3 結論

1)2010—2019年漓江流域桂林站斷面水質分析表明，桂林站斷面NH3-N濃度大于陽朔站斷面，而陽朔站斷面CODMn及TP濃度更高。兩水文站NH3-N、CODMn、TP年平均濃度分別為0.26、0.12 mg/L，1.63、1.56 mg/L，0.056、0.056 mg/L?？傮w水質情況較好，均能達到Ⅲ類水以上水平。且污染物通量與河流流量的相關系數(shù)大于污染物濃度，說明污染物通量主要取決于流量的變化。

2)根據(jù)箱線圖法異常值判定原理，對水質數(shù)據(jù)進行預處理，將異常值分別用平均值修正、直接刪除、視為刪失值處理進行修正，然后進行LOADEST模型模擬，取最優(yōu)擬合結果。得到用平均值修正處理效果最佳，與未處理結果相比較，桂林水文站NH3-N、CODMn、TP擬合效果分別提高了2.37%、2.82%、0%；陽朔水文站NH3-N、CODMn、TP擬合效果分別提高了2.89%、1.72%、3.75%，表明基于“預處理-回歸”過程的河流污染物通量估算方法提高了擬合精度。

3)利用2010—2019年漓江流域桂林及陽朔水文站的逐日流量數(shù)據(jù)及逐月水質數(shù)據(jù)，基于LOADEST模型分別建立了桂林及陽朔斷面的污染物通量回歸方程，以2010—2018年作為標定期，2018—2019年為驗證期，對回歸方程參數(shù)進行率定。結果表明，標定期和驗證期各污染物回歸方程相關系數(shù)R2均大于0.60，方程擬合效果較好。

圖7 桂林及陽朔水文站污染物逐年及多年月平均總通量

4)桂林及陽朔水文站斷面NH3-N、CODMn、TP年平均通量分別為987.302、1 254.781 t，6 784.407、14 964.549 t，265.945、756.055 t。河流污染物通量估算是科學保護漓江、維持優(yōu)良生態(tài)環(huán)境的前提，估算通量結果可為漓江流域污染物治理提供科學依據(jù)。