小流域非點源污染模擬與生態修復影響評價

吳敬東, 楊勝天, 葉芝菡, 胡曉靜, 張耀方

(1.北京市水科學技術研究院， 北京100048； 2.北京師范大學 地理學與遙感科學學院， 北京 100875)

非點源污染或面源污染(non-point source pollution)是指在降雨徑流(包括灌溉)的淋溶和沖刷作用下，大氣、地面和土壤中的污染物進入江河、湖泊、水庫和海洋等水體而造成的水污染[1]。目前，大量研究證明，非點源污染是導致水環境惡化的主要原因之一。據美國、日本等國的科學研究報道，即使點源污染得到全面控制，江河的水質達標率也僅為65%，湖泊的水質達標率為42%，海域水質達標率為78%。在我國，云南省滇池、天津市于橋水庫、安徽省巢湖、江蘇省太湖和遼寧省大伙房水庫等水域，其非點源污染比例均超過點源污染，非點源污染已上升為威脅人類飲用水的主要原因[2]。

當前，對密云庫上游非點源污染產生、排放和遷移轉換機理研究不是很深入，對流域不同區域種植業、養殖業和農村生活污水產生、排放系數及源強等重要參數還不掌握，對河湖水庫的水污染變化規律也不是很清楚，對各種非點源污染控制措施的效果效益研究也較少。因此，在點源控制技術相對成熟之后，對非點源污染機理及防控技術和措施的研究，成為當前非點源污染研究的熱點。密云水庫是北京城市供水的唯一大型地表水源，也是南水北調工程通水后余水存蓄庫，保護好密云水庫水源對北京城市供水系統的保障、調節和儲備安全有重大意義。由于人為活動的影響，密云水庫流域的生態環境發生了急劇的變化，入庫徑流污染物濃度有增加趨勢，水庫水質安全壓力大增。本文通過開展水源地非點源污染研究，摸清非點源污染主要來源、傳輸途徑及污染貢獻，旨在為保護密云水庫水源地、開展流域山水林田湖一體化修復和生態環境建設提供參考。

1 研究區概況

選擇北京密云水庫上游蛇魚川小流域作為研究區。蛇魚川小流域位于北京市密云水庫西北側，東經116°47′—116°57′，北緯40°35′—40°39′，屬于白河水系蛇魚川河流域，蛇魚川河直接注入密云水庫，小流域地跨密云水庫一級(西灣子村)和二級(黃峪口村)水源保護區。蛇魚川小流域面積25.66 km2。蛇魚川小流域屬燕山山脈，地勢北高南低，海拔高度為145～940 m。全流域坡度變化于0°～79°，坡度大于35°的面積占到全流域面積的56%，陡坡在流域內普遍存在。土壤以褐土為主，包括淋溶褐土和普通褐土。土壤質地較粗，多為砂壤土，顆粒松散，黏粒含量低，保水性較差，容易產生土壤侵蝕。流域屬暖溫帶大陸性半濕潤半干旱氣候，四季分明，冬季干燥寒冷，夏季炎熱，雨量集中且多暴雨，春季干旱多風。流域多年平均降水量為652 mm，75%集中在6—9月份。1月份平均氣溫-6.6 ℃，7月份平均氣溫24.8 ℃。無霜期為176 d，年日照時數為2 762 h，最大凍土深為85 cm。小流域植被類型有人工的針葉林、落葉闊葉林、灌叢、灌草從、草叢、水生植物以及經濟林、農田等植被。

2 研究方法

本研究從流域土地利用類型出發，利用EcoHAT-LCM模型和傳統大尺度非點源污染模型相結合的模式，分別計算各流域場次降雨徑流過程和非點源污染負荷，利用EcoHAT-LCM模型計算流域水文過程，再利用大尺度非點源污染模型，分別模擬流域多年平均、5 a一遇和10 a一遇暴雨產生的溶解態非點源污染負荷和吸附態非點源污染負荷，設置情景分析，評價小流域生態建設方案對密云水庫流域水文過程和污染防控效果影響。

2.1 EcoHAT-LCM模型

EcoHAT-LCM模型是劉昌明等[3-8]在陜西、甘肅、四川、青海、西藏和新疆等地區做了大量坡地降雨實驗，通過實驗途徑探討了暴雨徑流的超滲產流理論，并在這一理論的基礎上提出的暴雨徑流模型。

產流模型是LCM暴雨徑流模型的核心內容，其通過LCM入滲公式計算單位時段內入滲量，并耦合水量平衡方程，以參數優化選擇方法率定模型不確定性參數，再求算各個水文分量。對于地表及地表以下的產流過程，LCM模型采用出流系數法計算壤中流和地下徑流。

LCM模型中產流模型將整個過程分為地表層、土壤層和地下水層3步進行：首先計算土壤的下滲能力，采用經驗公式計算：

f=R·Pr

(1)

R=0.878 1×ln(r)+1.342 2

(2)

式中：f——下滲量/mm；P——降雨量/mm；R，r——經驗系數，與土地利用/覆蓋和土壤水含量有關，是計算產流過程的關鍵。

地表徑流量的計算基于水量平衡方程，即地表徑流量等于降水量減去下滲量：

Qd=P-f=P-R·Pr

(3)

計算壤中流量時，考慮壤中流與土壤濕度和降水入滲量成正比，因此，經驗公式為：

Ql=La·(Ws/Wsm)·f

(4)

式中：La——壤中流系數；Ws——非飽和土壤含水量(mm)；Wsm——土壤最大蓄水量(mm)；Ql——壤中流量(mm)。

地下水補給量與基流的計算同樣基于水量平衡方程，計算公式為：

REC=Rc·(Ws/Wsm)·(f-Ql)

(5)

Qb=Kb·(GWs+REC)

(6)

式中：REC——地下水補給量(mm)；Rc——地下水補給系數；Kb——基流系數； GWs——地下水儲量(mm)。

坡面匯流采用等流時線法。基于DEM數據，進行子流域劃分，結合流域等流時線分布圖，劃分子流域內等流時線空間分布，以等流時線區間為單元，計算流域內坡面匯流過程。子流域等流時線計算公式為：

Tsub_i=Ti-Tmin_i

(7)

式中：Tsub_i——調整后第i個子流域內某點匯流到子流域出口時長；Ti——第i個子流域內某點匯流到全流域出口時間；Tmin_i——第i個子流域內匯流到全流域出口最短時間。

河道匯流采用馬斯京根法。基于ArcGIS子流域屬性計算結果，結合子流域上下游拓撲關系，應用馬斯京根方程對子流域間河道匯流過程進行模擬計算，計算公式為：

Qout,2=C1·Qin,2+C2·Qin,1+C3·Qout,1

(8)

式中：Qin,1，Qin,2——河道時段初和時段末入流量(m3/s)；Qout,1，Qout,2——河道時段初和時段末出流量(m3/s)；C1，C2，C3——調節系數。

2.2 大尺度非點源污染模型

大尺度非點源污染模型是郝芳華、楊勝天、程紅光等[9-10]在充分借鑒統計性經驗模型和機理性過程模型優勢的基礎上，結合我國在非點源污染調查工作中的實際情況構建的，將非點源污染負荷產生和運移過程分成溶解態污染負荷和吸附態污染負荷分別進行：在溶解態污染模型中，按城市徑流、農村生活、畜禽養殖和農田徑流4種類型進行模擬；在吸附態污染模型中，先調用土壤侵蝕模型，估算出研究區內的土壤侵蝕量，然后借助構建的表層土壤氮磷含量數據庫，完成水蝕環境下吸附態氮磷的負荷量估算。

2.2.1 模型原理 從影響非點源污染產生要素來看，污染形成主要取決于兩個方面[11-15]： ①地表污染物的量級大小及存在形態，即源強要素； ②降水沖刷產流過程的影響程度，即產流要素。前者是流域社會經濟特征的綜合反映，包括居民的生活方式(如污水、垃圾排放)、城市發展水平(如排水管道鋪設、垃圾處理)、農業活動(如化肥、農藥的施用)和畜牧業發展情況等；后者則體現出區域的自然環境特征，包括氣候條件(如降水頻次、雨強、雨量)、土壤類型、地形特征和植被蓋度等，這兩方面因素是模型構建過程中主要考慮的因素。其次，從污染物產生、輸移和轉化規律看，非點源污染物有兩類：一類為溶解態污染物，另一類為吸附態污染物。溶解態污染物即表明該物質具有水溶性，能隨地表徑流一起發生遷移，整個過程受水循環控制，主要來自城市徑流、農田徑流、農村生活和畜禽養殖等人類活動；吸附態污染物即污染物通過附著于顆粒體(泥沙)實現遷移運動，其過程同水土流失密切相關。

2.2.2 溶解態模型 溶解態污染負荷采用二元結構模型，認為影響非點源產生和分布的因子主要是兩個方面，即區域的自然環境特征和社會經濟活動。這兩個方面的影響又分別歸納為土壤類型、土地利用、地形、降雨以及開發強度等。在上述因素中，地形、降雨和土壤結構屬于自然因子，而土地利用和開發強度屬于社會因子。綜合考慮社會、自然二元因子的影響，按農田、農村居民點、畜禽養殖(大牲畜、小牲畜)幾種類型分別計算。

(9)

式中：i——非點源污染類型，分別是城市徑流、農村居民點、畜禽養殖和農田徑流4個類型；n——非點源污染類型數；C——單位面積非點源污染負荷(t/km2)；Qi——單位面積非點源污染源強(t/km2)；ε——徑流系數；ε0——標準徑流系數，反映不透水硬化地面情況；k——地面沖刷系數；R——標準雨強(mm/d)；t——降雨歷時(d)；Ni——自然因子修正系數，表示自然因子修正，如在非點源污染產生過程中需要進一步修正的坡度、植被覆蓋和土壤等因子；Si——社會因子修正系數，表示社會發展程度對非點源污染源強的削弱程度。

2.2.3 吸附態模型 土壤侵蝕與非點源污染是一對密不可分的共生現象，在農業非點源污染中表現的尤其明顯。從本質上來說，土壤侵蝕產生的泥沙本身就是一種特殊的非點源污染物，會增加水體濁度，破壞水生生物的棲息生存環境，而且泥沙中往往攜帶了大量的有機物、磷酸鹽、銨離子等，從而造成受納水體污染。非點源污染類型泥沙的產生受到眾多外界因素的影響，污染源與受納水體間的空間距離、地形坡度與土壤糙度、地表植物的滯留作用、積水和低洼區等均會影響進入水體的輸沙量。吸附態非點源污染負荷受土壤侵蝕狀況的控制，其計算模型為：

Ca=X·Qa·η

(10)

式中：Ca——吸附態非點源污染負荷(t/km2)；X——土壤侵蝕量(t/km2)；Qa——吸附態非點源污染源強(%)；η——吸附態非點源污染富集系數。

2.2.4 土壤侵蝕模型 本研究采用修正的MUSLE方程計算土壤侵蝕量，計算公式為：

Sed=11.8·(RS·qpeak·Apixel)0.56·

KUSLE·CUSLE·PUSLE·LSUSLE·CFRG

(11)

式中：Sed——土壤流失量(t)；RS——地表徑流量(mm)；qpeak——洪峰徑流(m3/s)；Apixel——柵格單元面積(hm2)；KUSLE——土壤可蝕性因子；CUSLE——植被覆蓋和作物管理因子；PUSLE——保持措施因子；LSUSLE——地形因子； CFRG——粗碎屑因子。

洪峰徑流計算公式為：

(12)

式中：αtc——匯流時段內日降雨比例，對于日尺度模型取值1；tconc——柵格匯流時間，對于日尺度模型取值24； 3.6——單位轉換因子。

土地利用信息提取采用2017年GF-1衛星多光譜數據，數據空間分辨率為10 m；同時應用無人機機載多光譜傳感器獲取研究區內敏感區域下墊面數據，數據空間分辨率為1 m，再結合地面測量，按照3級分類標準，進行土地利用信息提取，獲取整個研究區的土地利用與土地覆蓋信息。

3 結果與分析

3.1 小流域土地利用特征

根據土地利用特征分析結果進行蛇魚川小流域土地利用類型專題圖制作，如附圖3和表1所示，小流域內林地面積最大，占流域面積的93%，其中天然林地2 056.49 hm2，占流域面積的80%，以板栗樹和核桃樹為主的經濟林(屬其他林地)占13%。林地面積與旱地和草地面積合計，流域植被覆蓋面積達到2 451.07 hm2，占全流域的96%。

表1 蛇魚川小流域土地利用類型面積統計

3.2 小流域EcoHAT-LCM模型計算

應用EcoHat-LCM模型分別模擬蛇魚川小流域多年平均、5 a一遇和10 a一遇的洪水徑流過程線。由于流域內無雨量站，無法獲取歷史雨量資料，選取離蛇魚川小流域最近的張家墳雨量站來代表當地的雨量狀況。

根據《北京水文手冊》，獲得張家墳雨量站多年平均、5 a一遇、10 a一遇的24 h暴雨量，不同頻率場次暴雨量詳見表2，同時雨量累積曲線采用《北京水文手冊》中確定的北京山區24 h降雨雨型。

表2 蛇魚川小流域不同頻率場次暴雨量

基于ASTER-GDEM數據進行研究區河網提取，子流域劃分和等流時線計算[16]，并對各子流域內部的等流時線進行調整，將等流時線劃分到各子流域中，實現對子流域內部坡面單元的劃分，為坡面匯流過程提供空間數據支持。

由于北京地區近年來降水持續偏少，氣候偏旱，因此模型土壤入滲系數選擇主要參考張亦弛[17]的黃河多沙粗沙區LCM模型分布式構建與應用碩士論文結果來分析蛇魚川小流域不同頻次設計暴雨的洪峰流量，入滲參數選擇如表3所示。

表3 偏旱狀況下的土壤入滲系數查找表

通過對設計暴雨數據進行IDW(inverse distance weight)空間插值處理，獲取小時尺度面雨量數據，結合DEM流域匯流屬性數據[18]，分布式構建確定LCM模型對降雨空間異質性的響應機制，基于土地利用、土壤質地等柵格數據對LCM模型下滲系數實現空間離散化處理，最后應用分布式LCM模型對蛇魚川小流域多年平均、5 a和10 a這3種設計頻率暴雨分別進行徑流模擬。偏旱狀況下蛇魚川小流域在多年平均、5 a一遇和10 a一遇場次暴雨下的洪峰流量，分別為34.07，76.45和131.80 m3/s，與吳敬東[19]用WMS模型的模擬結果近似，分別為39.82，81.43和132.89 m3/s。

3.3 小流域非點源污染負荷計算

應用大尺度非點源污染模型結合LCM模型的模擬結果，分別模擬蛇魚川小流域多年平均、5 a一遇和10 a一遇暴雨產生的溶解態非點源污染負荷和吸附態非點源污染負荷。首先模擬得到流域土壤侵蝕狀況，不同頻率場次暴雨下研究區土壤侵蝕的空間分布如圖1所示。多年平均暴雨下，蛇魚川小流域大部分坡面區域土壤侵蝕屬微度侵蝕，輕度侵蝕主要分布在山區旱地。5 a一遇暴雨下，蛇魚川小流域大部分坡面區域土壤侵蝕仍屬微度侵蝕，但在北部和西部部分山區，出現了輕度侵蝕。10 a一遇暴雨下，蛇魚川小流域大部分坡面區域以輕度侵蝕為主，其次是微度侵蝕，山區旱地為中度侵蝕或強度侵蝕。

對研究區在不同頻率場次暴雨下的土壤侵蝕量進行統計(詳見表4)。多年平均暴雨下全流域平均侵蝕模數為5.99 t/km2，微度侵蝕占全流域面積的98.3%，輕度侵蝕占1.7%；5 a一遇暴雨下全流域平均侵蝕模數為13.86 t/km2，微度侵蝕占80.7%，輕度侵蝕占19.3%，中度侵蝕占0.1%；10 a一遇暴雨下全流域平均侵蝕模數為24.96 t/km2，微度侵蝕占42%，輕度侵蝕占57.5%，中度侵蝕占0.4%，強度侵蝕占0.1%。

總體上，蛇魚川小流域由于植被覆蓋率高，占全流域的96%，土壤侵蝕以微度侵蝕為主，但是隨著暴雨量的增加，侵蝕級別向輕度侵蝕方向移動，當暴雨量為10 a一遇時，將出現強度侵蝕。不同頻率場次暴雨產生的溶解態非點源污染負荷和吸附態非點源污染負荷的空間分布分別如圖2—4所示。流域內溶解態N，P負荷最高的為經濟林地和居民點，其次為旱地，這3種土地利用類型受人類活動影響，施用大量化肥或產生大量垃圾、生活污水，導致溶解態非點源N，P負荷高。吸附態N，P負荷與土壤流失量密切相關，由于流域內土壤侵蝕主要發生在林地，通過林地水土流失帶走的吸附態N，P負荷相對較高。

表4 不同頻率場次暴雨下蛇魚川小流域土壤侵蝕統計情況

圖2 多年平均暴雨下產生的溶解態非點源污染負荷

對研究區在不同頻率場次暴雨下產生的非點源N，P污染負荷進行統計，并根據等流時線計算不同子流域污染負荷進入河道的總量(詳見表5)。蛇魚川小流域由于植被覆蓋率高，土壤侵蝕程度較輕，大部分坡面區域土壤侵蝕屬輕微侵蝕，導致流域總體吸附態污染負荷較其他區域偏小。以板栗樹為主的經濟林，對氮肥需求量大，導致流域溶解態N負荷較其他區域偏高。此外，不同頻次暴雨對吸附態N，P負荷的影響大于對溶解態N，P負荷的影響，如表6所示，隨著場次雨量的增加，吸附態N、P負荷占污染負荷總量的比重在逐漸增加，相應地溶解態N，P負荷占污染負荷總量的比重在逐漸減少。

圖3 5 a一遇暴雨下產生的溶解態非點源污染負荷

圖4 10 a一遇暴雨下產生的溶解態非點源污染負荷

表5 不同頻率場次暴雨下產生的非點源N，P污染負荷統計情況

表6不同頻率場次暴雨下溶解態和吸附態N、P污染負荷所占比例%

暴雨頻率溶解態N吸附態N溶解態P吸附態P多年平均435736645 a一遇3367287210 a一遇26742179

3.4 結果驗證

夏立忠等[20]對密云水庫流域的氮素流失進行了監測，密云水庫流域氮素流失通量為12.7 kg/hm2。吳敬東[19]對水土流失、農業生產、畜禽飼養和村莊生活等污染源負荷的計算，得到蛇魚川小流域年總氮流失量為91.21 t, 總磷31.28 t，折全流域平均氮素流失通量17 kg/hm2，磷素流失通量5.20 kg/hm2，農業生產的溶解態氮流失通量6.20 kg/hm2，生活污水的溶解態氮流失通量和溶解態磷流失通量分別為25.49和2.61 kg/hm2。

由于缺乏場次暴雨產污負荷的實測數據，統計研究區在不同頻率場次暴雨下單位面積產生的非點源N，P污染負荷，按照降水量占多年平均降水量的比重對場次暴雨的產污負荷進行折算，與夏立忠[20]和吳敬東[19]等人的研究結果進行對比，結論較為一致(表7)。

表7 蛇魚川場次暴雨產生的非點源污染負荷驗證

3.5 情景分析

情景設置是情景分析中的核心內容，直接決定情景分析的科學性以及由此提出建議、措施的合理性等。根據非點源污染負荷估算結果可知，農村居民點和農業生產活動為蛇魚川小流域非點源污染的主要來源。在此基礎上，本研究以研究區社會、經濟現狀為基礎，分別針對農村居民點和農業生產活動進行情景設置。其中基準情景(S0)采用研究區在現狀條件下3種頻率場次暴雨產生的非點源污染負荷的模擬結果。

(1) 情景1。農村居民點設置0.2的垃圾處理率和0.1的垃圾入網率，其他條件不變下，設置農村居民點具有小城鎮級別的垃圾處理率0.2和垃圾入網率0.1，通過模擬計算，農村生活的N，P污染負荷能夠削減28%(表3—7)。

(2) 情景2。經濟林施肥量削減為原施肥水平的50%，其他條件不變下，將經濟林的施肥量變為原施肥水平的1/2，通過模擬計算，農業生產活動的N，P污染負荷能夠削減43.9%和48.9%(表8)。

表8 不同情景條件下蛇魚川場次暴雨產生的非點源污染負荷模擬結果

注：3種情景設置的TN，TP負荷單位為t。

4 結 論

(1) 偏旱狀況下蛇魚川小流域在多年平均、5 a一遇和10 a一遇場次暴雨下的洪峰流量分別為34.07，76.45，131.80 m3/s。

(2) 蛇魚川小流域植被覆蓋率較高，土壤侵蝕以微度侵蝕為主，但是隨著暴雨量的增加，侵蝕級別向輕度侵蝕方向移動，當暴雨量為10 a一遇時，將出現強度侵蝕。蛇魚川流域由于以板栗樹為主的經濟林，對氮肥需求量大，導致流域溶解態N負荷較高。不同頻次暴雨對吸附態N，P負荷的影響大于對溶解態N，P負荷的影響，隨著場次雨量的增加，吸附態N，P負荷占污染負荷總量的比重在逐漸增加，相應地溶解態N，P負荷占污染負荷總量的比重在逐漸減少。

(3) 小流域內居民點具有小城鎮級別的垃圾處理率和垃圾入網率時，農村生活N，P污染負荷可以削減28%；減少經濟林的現有施肥量，蛇魚川小流域農業生產活動的N，P污染負荷能夠削減43.9%和48.9%。