近10年漓江干流水量水質(zhì)分析及污染物通量估算

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

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

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

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

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

1 材料與方法

1.1 研究區(qū)概況

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

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

1.2 箱線圖法異常檢測(cè)機(jī)制

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

圖2 箱線圖異常值判定

1.3 LOADEST模型原理

污染物通量一般指在一定時(shí)段內(nèi)污染物通過(guò)研究斷面的總量，可表示為

(1)

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

(2)

(3)

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

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

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

2 結(jié)果與分析

2.1 河流水文水質(zhì)特征

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

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

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

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

圖4 陽(yáng)朔水文站斷面實(shí)測(cè)通量與河流流量及污染物濃度的相關(guān)分析

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

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

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

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

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

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

表2 桂林及陽(yáng)朔水文站未處理與預(yù)處理結(jié)果對(duì)比

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

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

圖6 基于LOADEST模型的桂林及陽(yáng)朔水文站污染物通量擬合效果及驗(yàn)證結(jié)果

表5 桂林及陽(yáng)朔站污染物通量回歸方程

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

表6 桂林及陽(yáng)朔水文站參數(shù)檢驗(yàn)結(jié)果

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

3 結(jié)論

1)2010—2019年漓江流域桂林站斷面水質(zhì)分析表明，桂林站斷面NH3-N濃度大于陽(yáng)朔站斷面，而陽(yáng)朔站斷面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水質(zhì)情況較好，均能達(dá)到Ⅲ類水以上水平。且污染物通量與河流流量的相關(guān)系數(shù)大于污染物濃度，說(shuō)明污染物通量主要取決于流量的變化。

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

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

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

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