張新華，李紅霞，肖玉成，陳 奕

（四川大學(xué) 水力學(xué)與山區(qū)河流開發(fā)保護(hù)國家重點(diǎn)實(shí)驗室，四川 成都 610065）

1 研究背景

流域水污染防控中，點(diǎn)源污染的防治得到了充分的重視，并取得了一定的治理效果。但是，非點(diǎn)源污染引起的水體污染等問題日益突出，如湖庫水體富營養(yǎng)化、水華和水藻爆發(fā)都與流域內(nèi)的非點(diǎn)源污染有密切關(guān)系。由于非點(diǎn)源污染物的成分復(fù)雜、類型多樣、排放分散、量大面廣，定量分析較為困難［1-2］。如何綜合考慮點(diǎn)源和非點(diǎn)源污染負(fù)荷進(jìn)行總量控制與分配，是國內(nèi)外研究的熱點(diǎn)。

在污染總量控制與分配方面，負(fù)荷歷時曲線（LDC）由于其簡單易用，對資料要求較低，并且可以對點(diǎn)源和非點(diǎn)源污染同時分析，得到廣泛應(yīng)用，尤其在資料缺乏地區(qū)。如程艷等［3］以洱海彌苴河流域作為典型研究區(qū)域利用LDC進(jìn)行水質(zhì)分析，計算出超標(biāo)負(fù)荷的總削減率，丁京濤［4］利用LDC進(jìn)行了大寧河巫溪段水體總磷估算與分配研究。但LDC方法研究大都應(yīng)用在整個流域，沒有細(xì)致計算流域內(nèi)部單元的污染物分配情況。

贛江袁河流域流量、水質(zhì)監(jiān)測資料都較少，進(jìn)行污染模擬及總量控制與分配較為困難。本文以該流域為例，在進(jìn)行流域控制單元劃分的基礎(chǔ)上，利用SWAT模型的模擬獲得各控制單元的流量，然后利用LDC方法進(jìn)行資料缺乏流域的污染總量控制與分配方案研究。本文通過對污染負(fù)荷在不同水文條件下、不同控制單元和不同水質(zhì)目標(biāo)下的點(diǎn)源和非點(diǎn)源之間的分配的研究，構(gòu)建袁河流域分區(qū)（控制單元）、分期（豐水期、平水期、枯水期）、分類（點(diǎn)源和非點(diǎn)源）和分級（不同水質(zhì)目標(biāo)）的容量總量方案。

2 方法簡介

2.1 SWAT模型SWAT（Soil and Water Assessment Tool）模型，由美國農(nóng)業(yè)部農(nóng)業(yè)研究中心（US?DA-ARS）開發(fā)，是在SWRRB模型基礎(chǔ)上發(fā)展起來的一個大、中尺度的連續(xù)的流域環(huán)境模擬模型。該模型具有很強(qiáng)的水文物理機(jī)制，可以對不同的土壤條件、土地利用類型、氣候狀況和在人類活動的干擾下做出有效的徑流、水質(zhì)模擬預(yù)測［5］，SWAT模型在很多流域得到了成功應(yīng)用［6-10］。

SWAT模型將研究流域劃分為若干個水文響應(yīng)單元（Hydrologic Response Units，HRU），以HRU為最小的水文模擬單元進(jìn)行模擬［11］。每一個HRU內(nèi)的水量平衡是基于降水、蒸散、滲透、地表徑流、壤中流、地下水回流和河道運(yùn)移損失來計算。模型可以模擬流域內(nèi)部的多種水文循環(huán)物理過程：水的運(yùn)動、泥沙的輸移、植物的生長以及營養(yǎng)物質(zhì)的遷移轉(zhuǎn)化等。模型的整個模擬過程分為兩大部分：子流域模塊（產(chǎn)流和坡面匯流部分）和匯流演算模塊（河道匯流部分）。前者控制著每個子流域內(nèi)主河道的水、沙、營養(yǎng)物質(zhì)和化學(xué)物質(zhì)的輸入量，后者決定了物質(zhì)從河網(wǎng)向流域出口的輸移運(yùn)動及負(fù)荷的演算匯總過程。

2.2 LDC模型負(fù)荷歷時曲線法是Searcy在1959年首次提出，近年來在美國污染總量控制中得到廣泛應(yīng)用［12-13］，污染負(fù)荷歷時曲線法只考慮水體的稀釋過程，無需水力學(xué)數(shù)據(jù)（流速、斷面形狀、坡度、糙率、離散系數(shù)等），而只需要水文數(shù)據(jù)（流量或水位）。對于單一控制斷面（有水文水質(zhì)同步信息的斷面）的容量分析，其容量分布特征可以直接從負(fù)荷歷時曲線考察得出。將流量系列乘以特定污染物濃度標(biāo)準(zhǔn)得到負(fù)荷通量作為縱坐標(biāo)，建立的歷時曲線即為該污染物對應(yīng)流量系列的LDC。其步驟如下。

（1）繪制流域內(nèi)控制站點(diǎn)的日流量歷時曲線。根據(jù)流量資料（無或缺資料地區(qū)可以采用水文學(xué)上的面積比或相似流域法獲得），將流量Q從大到小排序就可以得到流量歷時曲線。

（2）將流量曲線轉(zhuǎn)變?yōu)槿肇?fù)荷歷時曲線。確定控制站所在區(qū)域的水功能區(qū)劃及其對應(yīng)的水質(zhì)要求，再將各種需要控制的污染物的容許濃度與流量歷時曲線各縱坐標(biāo)相乘得到允許負(fù)荷量，將允許負(fù)荷量與相應(yīng)的保證率建立一一對應(yīng)的關(guān)系，繪制允許負(fù)荷頻率曲線圖，即得到負(fù)荷歷時曲線。

（3）總量控制與分配。將實(shí)測（或預(yù)測的規(guī)劃水平年）的水質(zhì)數(shù)據(jù)轉(zhuǎn)化為日負(fù)荷并點(diǎn)繪制到日負(fù)荷歷時曲線圖上，比較實(shí)測負(fù)荷量與允許負(fù)荷量，即可知不同季節(jié)情況下是否超標(biāo)以及相應(yīng)的削減量。

3 研究區(qū)域及數(shù)據(jù)

袁河屬贛江水系一級支流，發(fā)源于江西省萍鄉(xiāng)縣羅霄山脈武功山七盤嶺金頂峰北麓，自西南向東北流經(jīng)萍鄉(xiāng)、宜春、分宜、新余以及樟樹部分地區(qū)，至樟樹鎮(zhèn)金鳳洲注入贛江（圖1）。該河全長279km，河床比降0.585‰，集水面積6 262km2。多年平均降水量為1 602.6mm，多年平均氣溫為17.3℃，多年平均蒸發(fā)量為1 282.9mm。流域內(nèi)植被良好，土壤大多為紅壤。

所收集到的數(shù)據(jù)包括DEM、土壤類型、土地利用類型、降雨、氣象、徑流、泥沙、水質(zhì)等。本文中袁河DEM為30m×30m；土壤類型圖為1∶20萬；土地利用類型圖為1∶100萬；降雨數(shù)據(jù)（1999—2008年）來自袁河流域上16個雨量站（圖1），氣象數(shù)據(jù)（1999—2008年）來源于袁河流域上4個氣象站（圖1）。

徑流數(shù)據(jù)來自袁河上中游的蘆溪站和茅洲2個水文站，其中蘆溪徑流數(shù)據(jù)從1978—2008年，而茅洲數(shù)據(jù)為1978—1994年?？紤]到數(shù)據(jù)長度的一致性，本文中模型率定階段為1999—2008年，因此需要將茅洲數(shù)據(jù)進(jìn)行延長。本文通過建立茅洲月平均徑流和蘆溪月平均徑流的線性相關(guān)方程（R2=0.91），將茅洲徑流延長到2008年。

袁河上的非點(diǎn)源污染來源主要考慮農(nóng)業(yè)化肥、農(nóng)村生活污水以及禽畜養(yǎng)殖。點(diǎn)源污染中主要考慮城市生活污水和工業(yè)污水，袁河沿干支流分布有62個排污口。袁河水質(zhì)監(jiān)測結(jié)果表明，袁河主要以Ⅲ類水為主，污染較為嚴(yán)重的區(qū)域主要在新余工業(yè)區(qū)。

4 SWAT模擬結(jié)果分析

根據(jù)研究區(qū)自然情況以及研究需要，本文將袁河流域共劃分為38個子流域（見圖1）。采用試錯法對SWAT模型進(jìn)行參數(shù)率定，即通過比較實(shí)測與模擬的徑流擬合情況來調(diào)整相關(guān)模型參數(shù)。由于袁河數(shù)據(jù)較少，本文將所獲得數(shù)據(jù)都用來進(jìn)行參數(shù)率定。通過對參數(shù)的敏感性分析以及參照其他文獻(xiàn)確定需要率定的參數(shù)，主要徑流參數(shù)包括：CN2（SCS曲線數(shù)）、ESCO（土壤蒸發(fā)補(bǔ)償系數(shù)）、AWC（土壤可供水量）、SLOP（子流域坡度）等。首先利用蘆溪和茅洲（延長后的）的徑流數(shù)據(jù)進(jìn)行徑流參數(shù)率定，然后輸出各控制單元的流量。

本文采用兩個評價指標(biāo)：線性相關(guān)系數(shù)R2和納西效率系數(shù)NSE，對模擬結(jié)果進(jìn)行分析。R2和NSE的計算公式如下：

式中：Xobs和Xsim分別為實(shí)測值和模擬值；分別為實(shí)測值和模擬值的算術(shù)平均值。

R2和NSE越接近于1表示實(shí)測與模擬值擬合越好，模擬精度越高。

圖2和圖3為蘆溪站和茅洲站的徑流的模擬結(jié)果。其中蘆溪徑流模擬的R2=0.92，NSE=0.78，茅洲徑流模擬的R2=0.89，NSE=0.73（表1）。一般認(rèn)為，當(dāng)R2＞0.6和NSE＞0.5時，模型模擬結(jié)果可以接受［12］。因此可以看出，SWAT對兩站的徑流模擬效果較好。此外，蘆溪站徑流模擬值比實(shí)測值略低，茅洲站徑流模擬值比實(shí)測值略高。但二者的模擬精度都遠(yuǎn)遠(yuǎn)超過了月徑流模擬的標(biāo)準(zhǔn)。

5 總量分配結(jié)果

5.1 控制單元劃分為了更有針對性地制定水環(huán)境保護(hù)方案，本文首先對研究流域進(jìn)行分區(qū)，然后進(jìn)行容量總量分配。分區(qū)原則主要是基于流域水環(huán)境生態(tài)系統(tǒng)的特征差異，即結(jié)合水文過程的集水范圍、行政分區(qū)、水功能分區(qū)以及劃分后單元便于管理統(tǒng)計等，體現(xiàn)流域水文特性、集水特性以及水功能保護(hù)。本文最終將袁河流域范圍劃分為4個一級控制單元（蘆溪、宜春、新余、樟樹），每個一級控制單元下又劃分為2～3個二級控制單元，即整個袁河流域共分為10個控制單元，從上游到下游依次為蘆溪保護(hù)區(qū)、蘆溪工業(yè)區(qū)、蘆溪保流區(qū)、宜春保留區(qū)1、宜春工業(yè)區(qū)、宜春保留區(qū)2、新余保留區(qū)1、新余工業(yè)區(qū)、新余保留區(qū)2和樟樹保留區(qū)，見圖4。

5.2 LDC總量分配結(jié)果根據(jù)水功能區(qū)劃可知各控制單元的水質(zhì)目標(biāo)，其中，保護(hù)區(qū)水質(zhì)目標(biāo)為Ⅱ類，保留區(qū)水質(zhì)目標(biāo)為Ⅲ類，工業(yè)區(qū)水質(zhì)目標(biāo)為Ⅳ類。通過比較污染負(fù)荷現(xiàn)狀值和水質(zhì)目標(biāo)值即判斷污染物的超標(biāo)情況。

流域內(nèi)控制單元內(nèi)均設(shè)有水質(zhì)監(jiān)測站，相應(yīng)站點(diǎn)的流量由SWAT模型模擬獲得，由此可建立其相應(yīng)的負(fù)荷歷時曲線。LDC計算結(jié)果如圖5所示（由于受篇幅限制只列出所有控制單元的氨氮及超標(biāo)控制單元的總氮圖）。

從各控制單元的負(fù)荷歷時曲線上可以看出，明顯超標(biāo)的有新余工業(yè)區(qū)和新余保留區(qū)2的總氮以及新余工業(yè)區(qū)的氨氮。為了判斷最關(guān)鍵水質(zhì)條件，采用上述代表值分別估算不同流量與實(shí)測水質(zhì)樣本的日負(fù)荷通量系列第90百分位數(shù)的濃度值相乘，與標(biāo)準(zhǔn)水質(zhì)對比得到削減率。

此外，LDC法可反映污染時的水文條件，其中，特豐條件：0～10%（歷時范圍）；豐水條件：10%～40%；平水條件：40%～60%；枯水條件：60%～90%；特枯條件：90%～100%。從圖中查看負(fù)荷超標(biāo)發(fā)生的季節(jié)，即可判斷污染類型：若發(fā)生在豐水期，則表示主要為非點(diǎn)源污染；若水質(zhì)超標(biāo)發(fā)生在枯水期，則表示主要為點(diǎn)源污染；若發(fā)生在平水期，則為點(diǎn)源與非點(diǎn)源的綜合。

由已計算出的LDC圖，按照上述算法算得各控制區(qū)各項指標(biāo)削減率并分析其污染類型，見表1。從圖5和表1可以看出，新余工業(yè)區(qū)內(nèi)枯水期氨氮負(fù)荷超標(biāo)率約為76%，豐水期超標(biāo)率約為47%。由于枯水期基本是點(diǎn)源污染，豐水期以非點(diǎn)源為主，同時包括點(diǎn)源污染。因此說明新余工業(yè)區(qū)氨氮超標(biāo)以點(diǎn)源為主；同樣，新余保留區(qū)枯水期總氮超標(biāo)約68%，豐水期超標(biāo)約53%，較新余工業(yè)區(qū)的非點(diǎn)源污染比重有所上升，但仍以點(diǎn)源為主。

在進(jìn)行削減分配計算中，視枯水期污染超標(biāo)全部為點(diǎn)源污染造成，豐水期污染為點(diǎn)源污染與非點(diǎn)源污染的綜合，即豐水期污染超標(biāo)總量扣除枯水期污染超標(biāo)量即為非點(diǎn)源污染分配量。同時應(yīng)考慮到，工業(yè)區(qū)位于保留區(qū)上游，當(dāng)工業(yè)區(qū)點(diǎn)源削減后，保留區(qū)點(diǎn)源削減量應(yīng)當(dāng)扣除工業(yè)區(qū)已經(jīng)削減的數(shù)值，若扣除后削減值小于等于零，則點(diǎn)源污染經(jīng)由上游削減后，下游可不用再削減；若扣除后大于零，則點(diǎn)源污染還需削減。此外，同一個控制分區(qū)內(nèi)，如有不同污染負(fù)荷削減，取其中最大值作為此控制分區(qū)的最終削減值，具體計算結(jié)果見表2。

從表2可以看出，對于新余工業(yè)區(qū)，豐水期和枯水期氨氮超標(biāo)，其中點(diǎn)源削減75.66%，非點(diǎn)源削減0%；同時豐水期和枯水期TN超標(biāo)，其中點(diǎn)源削減45.63%，非點(diǎn)源削減6.75%；取上述氨氮和TN削減比例的最大值，則新余工業(yè)區(qū)最終點(diǎn)源削減75.66%，非點(diǎn)源削減量為6.75%。對于新余保留區(qū)2，除去上游新余工業(yè)區(qū)的削減量后，新余保留區(qū)2的TN還應(yīng)削減：點(diǎn)源削減35.22%，非點(diǎn)源削減19.46%。總體來看，袁河新余區(qū)點(diǎn)源超標(biāo)嚴(yán)重，非點(diǎn)源略有超標(biāo)。

表1 LDC結(jié)果分析

表2 點(diǎn)源和非點(diǎn)源分配結(jié)果

以上分析表明，袁河流域污染物超標(biāo)河段主要在新余工業(yè)區(qū)及其下游的新余保留區(qū)，超標(biāo)污染物主要是點(diǎn)源污染，因此袁河流域應(yīng)嚴(yán)格控制流域內(nèi)點(diǎn)源的排放，同時對非點(diǎn)源污染進(jìn)行一定的防治。

6 結(jié)論

本文利用分布式SWAT模型和污染負(fù)荷歷時曲線法，構(gòu)建了資料缺乏的贛江袁河流域非點(diǎn)源污染模擬與總量控制與分配方案，實(shí)現(xiàn)了對贛江袁河流域的主要污染物的模擬以及基于“分區(qū)、分期、分類、分級”思想的水污染物總量分配。研究結(jié)果表明該方法在無資料流域應(yīng)用的可行性和適用性。本文所用技術(shù)和方法對其它資料缺乏流域的非點(diǎn)源污染模擬及總量分配也具有一定的借鑒意義和參考價值。

［1］郝芳華，程紅光，楊勝天.非點(diǎn)源污染模型——理論方法與應(yīng)用［M］.北京：中國環(huán)境科學(xué)出版社，2006.

［2］李海杰.吉林省雙陽水庫匯水區(qū)農(nóng)業(yè)非點(diǎn)源污染研究［D］.長春：吉林大學(xué)，2007.

［3］程艷，李炳花，和壽芳，等.水污染LDC在洱海彌苴河流域水質(zhì)分析中的應(yīng)用［J］.水電能源科學(xué)，2008，8：27-30.

［4］丁京濤.大寧河巫溪段水體總磷估算及分配研究［D］.北京：北京化工大學(xué)，2009.

［5］Arnold J G，Allen P M，Muttiah R，et al.Automated base-flow separation and recession analysis techniques［J］.Ground Water，1995，33（6）：1010-1018.

［6］Bosch N.The influence of impoundments on riverine nutrient transport：An evaluation using the Soil and Water Assessment Tool［J］.Journal of Hydrology，2008，355：131-147.

［7］Grizzetti B，Bouraoui F，Granlund K，et al.Modelling diffuse emission and retention of nutrients in the Vantaan?joki watershed（Finland）using the SWAT model［J］.Ecological Modelling，2003，169（1）：25-38.

［8］Santhi C，Arnold J G，William J R，et al.Validation of the SWAT model on a large river basin with point and nonpoint sources［J］.Journal of the American Water Resources Association，2001，37（5）：1169-1188.

［9］Saleh A，Arnold J G，Gassman P W，et al.Application of SWAT for the Upper North Bosque River Watershed［J］.Transactions of the ASAE，2000，43（5）：1077-1087.

［10］Neitsch S L，Arnold J G，Kiniry J R，et al.Soil and Water Assessment Tool Theoretical Documentation：Version 2000［Z］.Agricultural Research Service，Temple，Texas，2002.

［11］龐靖鵬.非點(diǎn)源污染分布式模擬--以密云水庫水源地保護(hù)為例［D］.北京：北京師范大學(xué)，2007.

［12］An Approach for Using Load Duration Curves in the Development of TMDLs［Z］.EPA 841-B-07-006，2007.

［13］Singh K P，Malik A，et al.Multivariate statistical techniques for the evaluation of spatial and temporal variations in water quality of Gomti River（India）——a case study［J］.Water Research，2004，38：3980-3992．