城市場次降雨徑流污染負(fù)荷快速估算方法*

何勝男，陳文學(xué) ，廖定佳，周 瑾，穆祥鵬

(1：中國水利水電科學(xué)研究院，北京 100038) (2：中國水利水電科學(xué)研究院流域水循環(huán)模擬與調(diào)控國家重點實驗室，北京 100038) (3：深圳市水文水質(zhì)中心，深圳518055) (4：深圳市東部水源管理中心，深圳 518172)

根據(jù)《2018年中國生態(tài)環(huán)境狀況公報》，2018年全國地表水中Ⅳ類和Ⅴ類占比達(dá)22.3%,劣Ⅴ類占比6.7%,全國地下水水質(zhì)監(jiān)測點中Ⅳ類占70.7%，Ⅴ類占15.5%[1]，我國地表及地下水面臨著不同程度的污染. 水污染控制和水環(huán)境保護(hù)的關(guān)鍵是控制污染物總量的排放，而其基礎(chǔ)和關(guān)鍵是定量水體污染負(fù)荷總量[2]. 對于城市水環(huán)境而言，城市降雨徑流污染是城市水質(zhì)惡化的重要原因之一[3-4]，定量計算城市降雨徑流污染負(fù)荷，可為城市水環(huán)境治理和污染控制提供科學(xué)依據(jù).

定量非點源污染負(fù)荷計算最早是美國在1960s-1970s開展的[5]，中國始于1980s[6]. 目前，定量城市降雨徑流污染負(fù)荷的計算方法較多，如濃度法[7-9]、統(tǒng)計方法[10]、推算方法[11]、反算方法[12-13]，這些方法大多用于估算年污染負(fù)荷，通常需要監(jiān)測多場場次徑流污染負(fù)荷，監(jiān)測難度和工作量大、費用高，且估算精度與監(jiān)測的降雨場次有關(guān). 目前，關(guān)于估算場次降雨徑流污染負(fù)荷研究相對較少. 隨著計算機(jī)的高速發(fā)展，數(shù)學(xué)模型的優(yōu)勢逐漸顯現(xiàn). 目前，通用的計算模型較多，如SWMM[14]、STORM[15]、DR3M-QUAL[16]、SLAMM[17]、HydroWorks[18]、HSPF[19]、MOUSE[20-21]等. 數(shù)學(xué)模型是定量分析非點源污染負(fù)荷的重要手段之一，通常利用有限場次降雨的監(jiān)測數(shù)據(jù)便可估算研究區(qū)降雨徑流污染負(fù)荷. 數(shù)學(xué)模型存在的主要問題是建模和參數(shù)率定時間較長，并且在預(yù)測分析平原城市的污染負(fù)荷時存在子流域劃分困難、河道流向難以確定等問題[22]，導(dǎo)致徑流量計算偏差較大，增加了污染物負(fù)荷計算的難度和預(yù)測精度.

為此，本文提出一種適用地區(qū)較廣的快速估算場次降雨徑流污染負(fù)荷的數(shù)學(xué)模型，即以污染物累積-沖刷理論為基礎(chǔ)，提出了“特征面積”的概念和計算公式，構(gòu)建了場次降雨徑流污染負(fù)荷數(shù)學(xué)模型，并結(jié)合多個案例，對數(shù)學(xué)模型的有效性、預(yù)測精度和適用性等方面進(jìn)行評價，以期為場次降雨徑流污染負(fù)荷的計算提供一種快速估算方法.

1 場次降雨徑流污染負(fù)荷估算方法

地表徑流污染負(fù)荷是指由降雨引起地表徑流排放的污染物總量，其中由一場降雨引起的稱為場次降雨污染負(fù)荷，由一年中多場降雨引起的稱為年污染負(fù)荷[23]. 對于徑流污染負(fù)荷，有學(xué)者采用徑流量與對應(yīng)污染物濃度的乘積進(jìn)行計算，但數(shù)據(jù)需要現(xiàn)場監(jiān)測，難度較大；另有學(xué)者通過多年監(jiān)測數(shù)據(jù)統(tǒng)計出的估算模型進(jìn)行計算，所需資料少，應(yīng)用簡便，如Schueler提出了城市開發(fā)區(qū)年地表徑流污染物計算模型[10]，該模型中污染負(fù)荷與地表徑流系數(shù)、匯水區(qū)面積、降雨量和污染物徑流量加權(quán)平均濃度呈正比. 地表徑流系數(shù)和污染物徑流量加權(quán)平均系數(shù)對模型的精度影響甚大. Thomson等[24]研究表明，污染物徑流量加權(quán)平均濃度的準(zhǔn)確估算至少需要15～20場降雨徑流觀測數(shù)據(jù).

污染物累積過程和降雨沖刷過程是決定降雨徑流污染負(fù)荷的兩個重要過程，描述累積過程的模型主要有線性、指數(shù)、對數(shù)等形式，其中以指數(shù)及一些變形形式相對比較成熟，而基于反映動力學(xué)假設(shè)的沖刷過程模型被人們廣泛使用[25]. 應(yīng)用廣泛的SWMM模型中采用的是污染物線性累積模型和指數(shù)沖刷模型[26]，見公式(1)和(2). 該模型能較好地反映污染物負(fù)荷機(jī)理.

(1)

式中，C1為最大累積量，kg/m2；C2為半飽和常數(shù)(達(dá)到最大累積物一半時所用的時間),d；t為晴天時間，d.

W=c1·qc2·B

(2)

式中，c1為沖刷系數(shù)；c2為沖刷指數(shù)；q為單位面積徑流速率，mm/h；B為污染物累積總量，kg/mm.

本文以污染物累積-沖刷理論為基礎(chǔ)，提出“特征面積”概念，以表征污染物在各類土地上的污染負(fù)荷特性即污染物累積特性和沖刷特性. 污染物累積特性用權(quán)重系數(shù)表示，即以污染物在各類土地上的最大累積量作為各類土地面積的權(quán)重，用以表征污染物累積強(qiáng)度，并以歸一化的方式計算特征面積的權(quán)重系數(shù)，即各類土地面積的權(quán)重系數(shù)等于污染物在各類土地上的最大累積量與在研究區(qū)所有土地上最大累積量的最大值的比值. 污染物沖刷特性用影響系數(shù)表示，以表征不同土地類型下降雨徑流沖刷特性的差異. 影響降雨沖刷特性的因素較多，主要可以分為兩類：一是降雨的基本參數(shù)和雨型，二是降雨區(qū)域的相關(guān)特征，包括城市下墊面類型、功能區(qū)類型和排水體制類型等，其中，下墊面條件對徑流量和徑流污染負(fù)荷有非常重要的影響. 本文將匯水區(qū)的土地利用類型分為透水區(qū)和不透水區(qū)兩大類，分別設(shè)置透水區(qū)影響系數(shù)和不透水區(qū)影響系數(shù). 特征面積與污染物最大累積量、土地利用性質(zhì)、匯水區(qū)面積、權(quán)重系數(shù)和影響系數(shù)有關(guān)，計算公式為：

Sa=wim·Sim+wm·Sm

(3)

根據(jù)污染物累積-沖刷模型，降雨徑流污染負(fù)荷與降雨量、匯水面積、污染物累積量、地表徑流系數(shù)等因素呈正比，而特征面積表征了污染物累積特性和不同土地利用類型的沖刷特性，因此，降雨徑流污染負(fù)荷與降雨量和特征面積呈正比，降雨徑流污染負(fù)荷數(shù)學(xué)模型可表示為：

M=k·Sa·P+b=k1·Sim·P+k2·Sm·P+b

(4)

式中，M為污染物負(fù)荷，kg；Sa為特征面積，hm2；Sim為不透水區(qū)特征面積，hm2；Sm為透水區(qū)特征面積，hm2；P為降雨總量，mm；k為斜率，kg/(hm2·mm)；b為截距，kg.k1=k·wim，k2=k·wm.

圖1 研究區(qū)各匯水區(qū)分布Fig.1 Distribution of catchments

斜率k1和k2表征降雨事件中不透水區(qū)和透水區(qū)降雨徑流污染沖刷特性，其值與不透水率或徑流系數(shù)有關(guān). 特征面積與降雨量的乘積具有體積單位，可定義為“特征體積”，因此，斜率k1和k2可理解為單位體積的降雨徑流污染負(fù)荷. 理論上講，當(dāng)降雨量為零時，污染量負(fù)荷為零，截距b也應(yīng)該為零. 但是，受模型簡化和測量誤差等因素的影響，截距可能并不等于零，因此，公式(4)中增加了截距項，以提高模型的預(yù)測精度.

2 場次降雨徑流污染負(fù)荷數(shù)學(xué)模型評價

2.1 材料與數(shù)據(jù)來源

2.1.1 研究區(qū)概況 以安徽省某縣城為例，研究區(qū)的總面積為43.6 km2，地勢平緩，地面高程平均海拔26.5～33.5 m，地面自然坡降為1/9000. 研究區(qū)分為9個匯水片區(qū)，分別為ST1、ST2、ST3、ST4、ST5、ST6、ST7、ST8和ST9，各匯水片區(qū)范圍見圖1. 研究區(qū)的土地類型分為路面、屋面和綠地，各匯水片區(qū)的面積及其包含各土地類型的占比和面積見表1.

表1 各匯水區(qū)土地利用占比情況

2.1.2 非點源污染模型 在建立場次降雨徑流污染負(fù)荷數(shù)學(xué)模型時，需要多場污染負(fù)荷結(jié)果. 由于SWMM非點源污染模型易于生成不同降雨徑流條件下的污染負(fù)荷樣本，因此本文采用SWMM非點源污染模型的模擬結(jié)果驗證模型. 以SWMM為平臺，對研究區(qū)和排水系統(tǒng)進(jìn)行概化，子匯水區(qū)690個，總面積為42.135 km2，管道876根，節(jié)點876個. 各子匯水區(qū)的不透水率和坡度利用地形資料和遙感影像借助ArcGIS計算得到. 管網(wǎng)計算選擇動力波，下滲模型選擇Horton模型，前期干旱時間選擇10 d，清掃街道去除率為60%，匯水區(qū)寬度系數(shù)、不透水區(qū)糙率、透水區(qū)糙率、不透水區(qū)洼蓄量、透水區(qū)洼蓄量、最大入滲率、最小入滲率、衰減常數(shù)、晴天時間分別取0.8、0.013、0.17、1 mm、3 mm、76.2 mm/h、3.81 mm/h、2 h-1、7 d[27]. 設(shè)計暴雨雨型選擇芝加哥雨型[28]，雨峰系數(shù)取0.4，設(shè)計重現(xiàn)期分別選取1、2、5、10、12、15和20 a(分別記為P1、P2、P5、P10、P12、P15、P20)，根據(jù)當(dāng)?shù)氐谋┯陱?qiáng)度計算公式，各降雨重現(xiàn)期下的降雨量分別為51、76、108、132、138、146和156 mm.

選取COD、TP和TN 3種污染物作為參考研究區(qū)降雨徑流污染負(fù)荷情況的指標(biāo). 根據(jù)SWMM手冊[26]和相關(guān)文獻(xiàn)[29-30]得到污染物累積和沖刷的參數(shù)，分別見表2和表3.

表2 污染物累積參數(shù)

表3 污染物沖刷參數(shù)

2.1.3 數(shù)據(jù)來源 本文中使用的兩個主要數(shù)據(jù)集：(1)研究區(qū)遙感數(shù)據(jù)，是從歐洲航天局(ESA，European Space Agency)的哨兵系列衛(wèi)星科研數(shù)據(jù)中心(Sentinels Scientific Data Hub)中以10 m的空間分辨率檢索土地利用數(shù)據(jù)(https://scihub.copernicus.eu/dhus/#/home)；(2)研究區(qū)的地形數(shù)據(jù)，是當(dāng)?shù)厝嗣裾峁┑?∶1000 城市地形規(guī)劃圖.

2.2 數(shù)學(xué)模型的有效性

對7種重現(xiàn)期分別進(jìn)行模擬，模擬時間為24 h(含退水時間). 假定有i個匯水片區(qū)，各匯水片區(qū)包含j個排水口，根據(jù)SWMM模型模擬結(jié)果即各排水口地表徑流量qi,j(τ)、對應(yīng)的地表徑流污染物含量ci,j(τ)和模擬的總時間，積分計算出一個排水口的污染負(fù)荷，對j個排水口的污染負(fù)荷進(jìn)行累積，即可得到各匯水片區(qū)污染負(fù)荷MSTi，其計算公式見公式(5). 研究區(qū)的污染負(fù)荷即所有匯水片區(qū)污染負(fù)荷的累加. 以污染物COD為例，各匯水片區(qū)在不同降雨條件下污染負(fù)荷見表4. 各匯水片區(qū)在不同降雨重現(xiàn)期下污染負(fù)荷與“特征體積”之間的關(guān)系分別見圖2.

(5)

式中，T為模擬的總時間, s；qi,j(τ)為第i個匯水片區(qū)包含的第j個排水口對應(yīng)匯水區(qū)徑流量, m3/s；ci,j(τ)為第i個匯水片區(qū)包含的第j個排水口對應(yīng)匯水區(qū)污染物濃度含量,mg/L.

圖2 不同降雨條件下COD(a)、TP(b)和TN(c)累積總量與特征體積的關(guān)系曲線Fig.2 The relationship between cumulated COD (a), TP (b) and TN (c) and characteristic volume under different rainfall conditions

從圖2中可見，降雨徑流污染負(fù)荷與特征體積具有較好的線性關(guān)系，其中，TN和TP的相關(guān)系數(shù)為0.97，COD的相關(guān)系數(shù)為0.96，說明降雨徑流污染負(fù)荷與“特征體積”即特征面積和降雨量的乘積呈正比，不同場次降雨下污染負(fù)荷可以用同一個公式進(jìn)行描述，也從側(cè)面說明了本文提出的估算模型是有效的. 由于各相關(guān)系數(shù)并不等于1，所以本文提出的數(shù)學(xué)模型的預(yù)測精度還需要做進(jìn)一步分析.

2.3 數(shù)學(xué)模型的預(yù)測精度

研究區(qū)路面和屋面為不透水區(qū)，綠地為透水區(qū)，根據(jù)各匯水片區(qū)中各類土地面積(表1)、污染物在各類土地上的最大累積量(表2)，結(jié)合公式(3)，計算各匯水片區(qū)對應(yīng)的不透水區(qū)特征面積和透水區(qū)特征面積；根據(jù)各降雨量和各匯水片區(qū)污染負(fù)荷(表4)，結(jié)合公式(4)，率定場次降雨徑流污染負(fù)荷數(shù)學(xué)模型. 本文設(shè)置3種情景，即分別選擇1場降雨(P5)、2場降雨(P5和P10)和3場降雨(P2、P5和P10)下的污染負(fù)荷，率定場次降雨徑流污染負(fù)荷數(shù)學(xué)模型，預(yù)測其他場次降雨徑流污染負(fù)荷，并與SWMM模型計算出的實際污染負(fù)荷進(jìn)行對比，分析預(yù)測精度.

表4 各匯水片區(qū)和研究區(qū)在不同降雨條件下的COD負(fù)荷

2.3.1 各匯水片區(qū)降雨徑流污染負(fù)荷預(yù)測精度 以COD為例，3種情景下的率定結(jié)果見表5，各匯水片區(qū)的污染負(fù)荷預(yù)測誤差如圖3所示.

圖3 各匯水片區(qū)污染物COD負(fù)荷預(yù)測誤差(a：情景一；b：情景二；c：情景三)Fig.3 Prediction error of COD pollutant load under different storm frequencies in each catchment (a: scenario 1; b: scenario 2; c: scenario 3)

表5 模型系數(shù)和影響系數(shù)計算結(jié)果

從表5中可以看出，不透水區(qū)影響系數(shù)幾乎是透水區(qū)影響系數(shù)的兩倍. 透水區(qū)不僅能減少地表徑流，也可控制降雨徑流污染，低影響開發(fā)技術(shù)即是利用透水路面、植物透水磚、生物滯留池等技術(shù)增加透水區(qū)面積，降低降雨徑流污染負(fù)荷. 影響系數(shù)的差異也表明本文提出的“特征面積”的合理性，反映了匯水區(qū)上污染物累積特性和不透水率差異引起的地表徑流沖刷特性的差異.

從圖3可以看出，3種情景下，P1和P2的降雨徑流污染負(fù)荷的預(yù)測誤差較大，P2以上的預(yù)測誤差相對較小，預(yù)測精度較高.

圖4 3種情景下各污染物在不同重現(xiàn)期下的預(yù)測誤差Fig.4 Prediction error of each pollutant load under the three scenarios and different storm frequencies

2.3.2 研究區(qū)降雨徑流污染負(fù)荷預(yù)測精度 以整個研究區(qū)為對象，對3種情景下降雨徑流污染負(fù)荷進(jìn)行預(yù)測精度分析，計算誤差如圖4所示.

從圖4中可以看出，對于COD，情景一時，P1的預(yù)測誤差為-28.85%,P2及以上的預(yù)測誤差在-7.3%～9.8%之間；情景二時，P2及以下的誤差在-18.6%～30.0%之間,P12及以上的預(yù)測誤差在-1.4%～7.0%之間；情景三時，P2及以下的誤差在-29.2%以內(nèi)，P12及以上的預(yù)測誤差在-5.0%～0.78%之間. 對于TP，情景一時，P1的預(yù)測誤差為37.62%,P2的預(yù)測誤差為14.02%，P10及以上的預(yù)測誤差在-10.10%～-1.82%之間；情景二時，P1的預(yù)測誤差為36.02%,P2的預(yù)測誤差為13.7%，P12及以上的預(yù)測誤差在-9.6%～-9.2%之間；情景三時，P1的預(yù)測誤差為26.61%,P12及以上的預(yù)測誤差約為-11%. 對于TN，預(yù)測誤差均在±11%以內(nèi).

綜上，本文提出場次降雨徑流污染負(fù)荷數(shù)學(xué)模型，對于P2以上的降雨徑流污染負(fù)荷預(yù)測精度較高，其相對誤差小于±11%，此外，降雨場次對數(shù)學(xué)模型的預(yù)測精度影響不大. 模型中有3個系數(shù)需要率定，因此，利用3場以上場次污染負(fù)荷結(jié)果便可預(yù)測其他場次的污染負(fù)荷，從而大大減小獲取污染負(fù)荷數(shù)據(jù)的工作量. 對于2 a一遇以下的降雨情況，模型預(yù)測精度相對較低，其主要原因是：小降雨情況下，下墊面的洼蓄量和雨水下滲量在總降雨量中占比相對較大，采用大降雨事件的污染負(fù)荷結(jié)果率定模型參數(shù)與小降雨事件下有一定的差異，使得計算誤差相對較大. 因此，對于小降雨工況而言，可以采用降雨量相近的數(shù)據(jù)進(jìn)行率定，以提高模型的預(yù)測精度.

2.4 數(shù)學(xué)模型的適用性

選擇武漢市和蘇州市兩個地區(qū)的實際案例進(jìn)行分析和說明.

1)案例一

潘璐[31]研究了武漢市典型校區(qū)的降雨徑流污染負(fù)荷. 研究區(qū)坡度在0.13%～0.50%之間. 經(jīng)現(xiàn)場采樣，TN和CODMn為主要污染物. 地面分為路面、屋面、庭院/廣場、操場和綠地，CODMn在對應(yīng)土地類別上的最大累積量分別為80、60、45、45和30 kg/hm2；TN在對應(yīng)土地類別上的最大累積量分別為6、4、6、6和10 kg/hm2. 采用芝加哥雨型，雨峰系數(shù)0.4，降雨歷時120 min，干旱時間10 d. 研究區(qū)降雨總量和各降雨條件下的污染負(fù)荷等基礎(chǔ)信息如表6所示.

表6 武漢市典型校區(qū)降雨徑流污染負(fù)荷[31]

圖5 案例一中CODMn和TN污染負(fù)荷總量預(yù)測誤差Fig.5 Prediction error of total CODMn and TN pollution load of case 1

由于數(shù)據(jù)數(shù)量相對較少，因此，根據(jù)P1、P3、P5 3場降雨條件下的污染負(fù)荷結(jié)果率定模型，預(yù)測P0.5和P10下的污染負(fù)荷，預(yù)測精度均在98.5%以內(nèi)，精度較高，計算誤差見圖5.

2)案例二

祁繼英[32]研究了蘇州市南園水系排水區(qū)的降雨徑流污染負(fù)荷. 主要污染物為COD、TN和TP，地面分為屋面、綠地和交通道路. 各污染物在對應(yīng)土地類別上的最大累積量見表2. 采用芝加哥雨型，雨峰系數(shù)為0.367，降雨歷時60 min，干旱時間10 d. 研究區(qū)降雨總量和各降雨條件下的污染負(fù)荷總量等基礎(chǔ)信息見表7.

表7 蘇州南園水系降雨徑流污染負(fù)荷[32]

圖6 案例二中COD、TN和TP污染負(fù)荷總量預(yù)測誤差Fig.6 Prediction error of total COD and TN pollution load in case 2

根據(jù)P0.5、P1、P2 3場降雨條件下的污染負(fù)荷結(jié)果率定模型，預(yù)測P0.25、P3和P5下的污染負(fù)荷，預(yù)測精度均在99.5%以內(nèi)，精度較高，計算誤差見圖6.

2.5 數(shù)學(xué)模型的局限性

本文構(gòu)建了場次降雨徑流污染負(fù)荷數(shù)學(xué)模型，預(yù)測精度較好，但其也有一定的局限性：

1)對于小降雨事件，模型預(yù)測精度比其他降雨條件相對較低. 為提高模型預(yù)測精度，建議采用降雨量相近的觀測結(jié)果率定模型參數(shù).

2)降雨間隔影響污染物在地面上的最大累積量，而不同土地上的污染物最大累積量是影響場次降雨徑流污染數(shù)學(xué)模型預(yù)測精度的關(guān)鍵參數(shù)，該參數(shù)與土地利用性質(zhì)和狀況、綠化條件、交通狀況以及土地裸露程度直接相關(guān)[4]，通?？稍诟黝愅恋厣线x取有代表性的1 m2區(qū)域定時監(jiān)測分析獲得[33]，也可參考相關(guān)文獻(xiàn)，或根據(jù)經(jīng)驗確定(如查閱SWMM手冊等). 有分析表明[34]，一場典型降雨48 h后地表污染物才能恢復(fù)原狀. 為此，對于連續(xù)降雨或者降雨間隔時間較短導(dǎo)致地面污染物在各類土地上的最大累積量發(fā)生變化，當(dāng)全域最大累積量與局部最大累積量線性變化時，權(quán)重系數(shù)并不會發(fā)生較大變化，不會影響模型預(yù)測精度. 倘若發(fā)生較大變化且不是線性變化時，為了提高模型的計算精度，模型還需重新率定.

3 結(jié)論

本文以污染物累積-沖刷理論為基礎(chǔ)，提出了反映污染負(fù)荷特性的“特征面積”，建立了場次降雨徑流污染負(fù)荷數(shù)學(xué)模型，并結(jié)合多個案例，分析了數(shù)學(xué)模型的有效性、預(yù)測精度、適用性和局限性. 主要結(jié)論如下：

1)特征面積較好地反映了污染物在各類土地上的污染負(fù)荷特性即污染物累積特性和沖刷特性，經(jīng)實例驗證場次降雨徑流污染負(fù)荷與特征面積和降雨量的乘積呈正比.

2)利用3場及以上降雨徑流污染負(fù)荷結(jié)果，可較好地率定場次降雨徑流污染負(fù)荷數(shù)學(xué)模型，從而可快速且較準(zhǔn)確地估算單場次降雨徑流污染負(fù)荷. 該方法簡單實用，獲取數(shù)據(jù)工作量小，適用地區(qū)廣.

3)對于小降雨事件，建議采用降雨量相近的觀測結(jié)果對模型進(jìn)行率定，以提高模型的預(yù)測精度.