朱凱航，陳磊*，王怡雯，劉國(guó)王辰，顏小曼，郭晨茜，張亮，沈珍瑤

（1.北京師范大學(xué)環(huán)境學(xué)院，北京 100875；2.中國(guó)科學(xué)院精密測(cè)量科學(xué)與技術(shù)創(chuàng)新研究院，武漢 430077）

水土流失和面源污染是兩個(gè)密切相關(guān)的全球環(huán)境問題[1-2]。在我國(guó)，人為擾動(dòng)和不當(dāng)?shù)耐恋乩梅绞郊铀倭怂亮魇н^程，2021 年全國(guó)水土流失面積為267.42 萬km2，占國(guó)土面積的27.96%。水土流失過程中不僅土壤顆粒被沖刷進(jìn)入河道，引起泥沙淤積，而且其伴生的面源污染還會(huì)引起水質(zhì)惡化，導(dǎo)致湖泊富營(yíng)養(yǎng)化等環(huán)境問題[3-4]。水土流失、面源污染二者互相影響、互相制約，共同導(dǎo)致了水土流失型面源污染。

水土流失型面源污染涉及水、土、污染物等多介質(zhì)耦合作用，主要包含降雨徑流、土壤侵蝕和泥沙/污染物運(yùn)輸3 個(gè)過程[5]，涉及農(nóng)田、溝渠、水塘、濕地、河濱等多個(gè)環(huán)節(jié)和田塊、灌排單元、集水區(qū)、流域等多個(gè)尺度[6]。水土流失型面源污染的產(chǎn)生和傳輸具有不確定性和復(fù)雜性，對(duì)其關(guān)鍵過程的正確模擬與評(píng)估對(duì)污染入水體量有重要影響[7]。目前多采用面源污染模型對(duì)水土流失和面源污染進(jìn)行模擬，如SWAT（Soil and Water Assessment Tool）、HSPF（Hydrological Simu?lation Program：Fortran）、AGNPS（Agricultural Nonpoint Pollution Source Model）等[8-9]。但是此類模型集總式或半分布式的建模思路忽略了集水區(qū)內(nèi)部地形特征、徑流路徑、景觀格局等因素的差異對(duì)水土流失和污染物傳輸造成的影響，并且無法考慮水土流失和污染物在傳輸中的損失[7]，不適用于我國(guó)中西部景觀破碎度較高、空間異質(zhì)性強(qiáng)的丘陵地區(qū)。因此，充分考慮水土流失型面源污染的微觀機(jī)理，精準(zhǔn)模擬水土流失型面源污染的陸面?zhèn)鬏斶^程，識(shí)別水土流失型面源污染優(yōu)先控制區(qū)，是準(zhǔn)確評(píng)估水土流失型面源污染風(fēng)險(xiǎn)和形成防控技術(shù)的基礎(chǔ)[10-11]。

本文針對(duì)現(xiàn)有半分布式模型的局限性，構(gòu)建了耦合田塊尺度富集/流失、流域尺度傳輸及河道遷移過程的分布式水土流失型面源污染模型，開發(fā)出高效柵格匯流算法，破解了對(duì)丘陵山區(qū)水土流失型面源污染精準(zhǔn)評(píng)估和防控的技術(shù)瓶頸問題。進(jìn)而本研究以三峽庫(kù)區(qū)菁林溪流域?yàn)榘咐齾^(qū)，對(duì)新模型進(jìn)行驗(yàn)證，基于新模型識(shí)別了三峽庫(kù)區(qū)典型流域的水土流失型面源污染時(shí)空分布規(guī)律，評(píng)價(jià)了典型措施的防控效果，以期為水土流失型面源污染防控提供機(jī)理和技術(shù)支持。

1 材料與方法

1.1 研究區(qū)概況

菁林溪流域位于三峽庫(kù)區(qū)澎溪河上游（31.19°~31.27°N，108.34°~108.38°E，圖1），流域面積9.51 km2，高程位于665~1 117 m之間，流域上游高程差較大，中下游相對(duì)平坦，農(nóng)業(yè)種植較發(fā)達(dá)[12]。該地多年平均氣溫22.1 ℃，平均降水1 293 mm，5—9 月的降水量占全年的66%以上。本研究日氣象數(shù)據(jù)來源于重慶市統(tǒng)計(jì)年鑒。DEM 數(shù)據(jù)來源于中國(guó)科學(xué)院資源環(huán)境科學(xué)與數(shù)據(jù)中心（http：//www.resdc.cn/）。土壤類型數(shù)據(jù)參考呂明權(quán)[13]對(duì)研究區(qū)土壤樣品進(jìn)行實(shí)測(cè)得到的數(shù)據(jù)，主要有三類土壤：粗骨土、水稻土和紫色土。土地利用類型來源于對(duì)高精度遙感影像的人工劃分，主要由林地、水田、坡耕地、建設(shè)用地、河道及池塘等水體組成。作物管理數(shù)據(jù)（作物種類、種植時(shí)間、化肥施用情況）通過實(shí)地調(diào)研獲取。本研究開展了多年的流域水量水質(zhì)觀測(cè)，使用自動(dòng)水位檢測(cè)儀HOBO UL20獲取每日的流量數(shù)據(jù)，通過人工采集后在實(shí)驗(yàn)室測(cè)定獲取水質(zhì)數(shù)據(jù)，采樣點(diǎn)位于流域主干道，于每月初、月中以及降雨日進(jìn)行采樣，使用0.45μm濾膜過濾水樣，以鉬酸銨分光光度法測(cè)定水樣中的總磷與溶解態(tài)磷濃度，吸附態(tài)磷濃度為二者之差，相關(guān)數(shù)據(jù)用于后續(xù)模型構(gòu)建與驗(yàn)證。

圖1 菁林溪流域地理位置與土地利用狀況Figure 1 Location and land use of Jinglinxi watershed

1.2 模型構(gòu)建

面源污染物在降雨、地形等因素驅(qū)動(dòng)下離開產(chǎn)污單元，并伴隨泥沙和徑流進(jìn)入河道。本研究的柵格模型考慮了土壤原位剝離、溶解態(tài)和吸附態(tài)面源污染物流失及其微觀傳輸規(guī)律，重點(diǎn)包括田塊富集/流失模塊、流域匯流/攔截模塊、河道傳輸模塊。

1.2.1 田塊富集/流失模塊

（1）降雨產(chǎn)流模擬

由于大部分流域都沒有降雨-徑流的歷史數(shù)據(jù)，因此選取SCS-CN 法，該方法因?qū)?shù)據(jù)要求較低、對(duì)徑流模擬效果較好而得到廣泛應(yīng)用[14]。在SCS-CN方法中，首先對(duì)研究區(qū)的降雨量數(shù)據(jù)進(jìn)行空間插值，計(jì)算月降雨的空間分布情況，再結(jié)合相關(guān)研究，對(duì)研究區(qū)曲線數(shù)（CN）和初損量進(jìn)行修正，計(jì)算地表徑流，計(jì)算公式為：

式中：Qsurf為地表徑流，mm；Rday為降雨量，mm；Ia為初損量，包括產(chǎn)流前的地面填洼量、植物截留量和下滲量，mm；S為最大蓄水量，mm。

式中：CN表示某天的曲線數(shù)。Ia通常近似為0.2S，因此公式（1）變?yōu)椋?/p>

（2）土壤侵蝕模擬

模型采用RUSLE 方程對(duì)泥沙的原位負(fù)荷量進(jìn)行模擬[15]：

式中：sed為某柵格的泥沙日負(fù)荷量，t·d-1；R為降雨侵蝕力因子，MJ·mm·hm-2·h-1·d-1；K為土壤可蝕性因子，t·hm2·h·MJ-1·mm-1·hm-2；L為坡長(zhǎng)因子，無量綱；S為坡度因子，無量綱；C為植被覆蓋與管理因子，無量綱；P為水土保持措施因子，無量綱；cellsize為柵格的面積，30 m×30 m的柵格面積為0.09 hm2。

降雨侵蝕力因子（R）為降雨作用對(duì)土壤剝蝕的能量因子，本研究選擇基于逐日降水?dāng)?shù)據(jù)的R因子計(jì)算方式，計(jì)算公式為：

式中：Pi為當(dāng)日降水量，降雨量不足12 mm 時(shí)取0；Pd12和Py12分別為降水量大于12 mm的日平均降水量和年平均降水量，mm；α和β為計(jì)算參數(shù)，與降雨強(qiáng)度和降雨頻率相關(guān)。

土壤可蝕性因子（K）定量描述了土壤對(duì)降雨侵蝕的抵抗能力，其主要受土壤質(zhì)地、土壤顆粒組成和土壤結(jié)構(gòu)等物理因素的影響。由于研究區(qū)位于三峽庫(kù)區(qū)，本研究采用吳昌廣等[16]基于重慶市和湖北省各土種的理化性質(zhì)數(shù)據(jù)庫(kù)提出的三峽庫(kù)區(qū)土壤可蝕性。

土壤侵蝕量隨著坡長(zhǎng)和坡度的增加而增加，坡長(zhǎng)因子（L）和坡度因子（S）全面反映了地形對(duì)土壤侵蝕造成的影響，計(jì)算公式為：

式中：λ為坡長(zhǎng)的水平投影，m；m為坡長(zhǎng)指數(shù)，無量綱；θ為坡度，（°）。

植被覆蓋與管理因子（C）反映了植被對(duì)降雨徑流和水土流失的攔截作用。本研究基于歸一化植被指數(shù)推算C因子值，遙感數(shù)據(jù)來自美國(guó)國(guó)家航空航天局（NASA）的 EOS/MODIS 數(shù)據(jù)產(chǎn)品（https：//ladsweb.modaps.eosdis.nasa.gov），計(jì)算公式為：

式中：NDVI為歸一化植被指數(shù)；NDVImin和NDVImax為一定置信區(qū)間內(nèi)的NDVI最小值和最大值（3%~97%）。

水土保持措施因子（P）為采取特定水土保持措施下的水土流失量與未采取任何保持措施情況下的水土流失量之比，林地、旱地、水田、水域和建設(shè)用地的P因子分別為0.90、0.35、0.25、0和1.0[17]。

（3）土壤磷富集、流失模擬

模型的磷循環(huán)模塊借鑒SWAT 模型磷庫(kù)的思路[18]，將土壤中的無機(jī)磷劃分為溶解態(tài)無機(jī)磷庫(kù)、活性態(tài)無機(jī)磷庫(kù)和穩(wěn)定態(tài)無機(jī)磷庫(kù)，將有機(jī)磷劃分為活性態(tài)有機(jī)磷庫(kù)、穩(wěn)定態(tài)有機(jī)磷庫(kù)和新生有機(jī)磷庫(kù)，磷循環(huán)模塊在柵格的尺度上進(jìn)行，磷循環(huán)過程如圖2所示。

圖2 土壤磷循環(huán)流程Figure 2 Circulation of soil phosphorus

土壤中的可溶性磷主要通過擴(kuò)散作用遷移，由于地表徑流僅與表層10 mm 土層中的部分可溶性磷相互作用，因此某柵格某月的溶解態(tài)磷負(fù)荷量計(jì)算公式為：

式中：Psol為某柵格的溶解態(tài)磷月負(fù)荷量，kg；Psol，surf，i為某月第i天表層10 mm 土層中的可溶性磷量，kg·hm-2；Qsurf，i為某月第i天的地表徑流量，mm；ρb為表層10 mm 土層的容重，Mg·m-3；depthsurf為土壤表層的深度，10 mm；kd，surf為磷的土壤分配系數(shù)，是表層 10 mm土層中的可溶性磷濃度與地表徑流中的可溶性磷濃度的比值，m3·Mg-1；cellsize為柵格的面積，0.09 hm2；m為某月的天數(shù)，d。

吸附態(tài)磷以泥沙的形式進(jìn)入河道，某柵格某月的吸附態(tài)磷負(fù)荷量為：

式中：Psed為某柵格的吸附態(tài)磷月負(fù)荷量，kg；concsedP，i為某月第i天表層10 mm 土層中吸附在泥沙上的磷量，g·t-1；sedi為某月第i天的泥沙月負(fù)荷量，t；εP：sed為磷的富集比。

表層10 mm 土層中吸附在泥沙上的磷量（conc?sedP）計(jì)算公式為：

式中：minPact，surf為表土層中活性態(tài)無機(jī)磷庫(kù)中的磷含量，kg·hm-2；minPsta，surf為表土層中穩(wěn)定態(tài)無機(jī)磷庫(kù)中的磷含量，kg·hm-2；orgPhum，surf為表土層中腐殖質(zhì)有機(jī)磷庫(kù)中的磷含量，kg·hm-2；orgPfrsh，surf為表土層中新生有機(jī)磷庫(kù)中的磷含量，kg·hm-2；ρb為表層 10 mm 土層的容重，Mg·m-3；depthsurf為土壤表層的深度，10 mm。

富集比（εP：sed）定義為隨泥沙遷移的磷含量與土壤表層中磷含量的比值，計(jì)算公式為：

式中：concsed，surf為地表徑流中的含沙量，Mg·m-3；sed為某日的產(chǎn)沙量，t；cellsize為柵格的面積，0.09 hm2；Qsurf為某日的地表徑流量，mm。

1.2.2 流域匯流/攔截模塊

傳統(tǒng)的柵格匯流算法是基于源頭到出口的正向思路實(shí)現(xiàn)的[19]，但是正向計(jì)算的思路存在計(jì)算效率低的問題，因此本研究提出了結(jié)合反向遞歸及哈希表的優(yōu)化算法。具體步驟如下：

（1）基于D8算法獲得柵格之間的上下游關(guān)系，構(gòu)建一個(gè)字典保存柵格之間的上下游關(guān)系，字典的值保存了某一柵格所有上游柵格的集合，鍵則為下游柵格（圖3b）。字典基于哈希算法實(shí)現(xiàn)，鍵名是唯一、不能重復(fù)的，因此只需要以接近常量的時(shí)間便可以查詢某一柵格的所有上游柵格，無需遍歷字典中的所有元素。在查詢某一柵格的上游柵格時(shí)，查詢時(shí)間并不會(huì)隨著流域面積（字典長(zhǎng)度）的增大而增大。

（2）基于遞歸算法，從流域出口反向索引柵格之間的上下游關(guān)系，對(duì)全流域柵格進(jìn)行匯流計(jì)算。遞歸算法是在函數(shù)中直接或調(diào)用自身的算法。在匯流計(jì)算時(shí)，對(duì)流域出口柵格進(jìn)行反向遞歸索引：首先在字典中查詢?cè)摉鸥?，如果不存在則表明該柵格不存在上游柵格，因此可直接基于該柵格的源圖層進(jìn)行計(jì)算；如果該柵格存在于字典中，說明該柵格存在上游柵格，則需對(duì)其所有上游柵格再次進(jìn)行反向遞歸索引，重復(fù)以上判斷，待所有上游節(jié)點(diǎn)完成模擬并返回結(jié)果時(shí)，再基于上游柵格的模擬結(jié)果及該柵格的源圖層進(jìn)行模擬。以圖3a 所示的傳輸路徑為例，其模擬步驟如下：

圖3 模型構(gòu)建流程圖Figure 3 Flow chart of model construction

①?gòu)某隹跂鸥? 號(hào)柵格開始索引，在傳輸字典中索引到了它的上游柵格2 號(hào)柵格，于是對(duì)2 號(hào)柵格進(jìn)行索引，待2號(hào)柵格的模擬結(jié)果返回后，再對(duì)4號(hào)柵格進(jìn)行模擬。

②對(duì)2 號(hào)柵格進(jìn)行索引，在傳輸字典中索引到了它的上游柵格1 號(hào)柵格和3 號(hào)柵格，于是先對(duì)1 號(hào)柵格和3 號(hào)柵格進(jìn)行模擬，待1 號(hào)柵格和3 號(hào)柵格的模擬結(jié)果返回后，再對(duì)2號(hào)柵格進(jìn)行模擬。

③對(duì)1 號(hào)柵格進(jìn)行索引，發(fā)現(xiàn)1 號(hào)柵格無上游柵格，則直接進(jìn)行該柵格的模擬，并將模擬結(jié)果返回給2號(hào)柵格；對(duì)3號(hào)柵格進(jìn)行索引，索引到了它的上游柵格5號(hào)柵格和6號(hào)柵格，于是先對(duì)5號(hào)柵格和6號(hào)柵格進(jìn)行模擬，基于5 號(hào)柵格和6 號(hào)柵格的入流情況對(duì)3號(hào)柵格進(jìn)行模擬。

④對(duì)5 號(hào)柵格和6 號(hào)柵格進(jìn)行索引，發(fā)現(xiàn)二者無上游柵格，則直接進(jìn)行模擬，并將模擬結(jié)果返回給3號(hào)柵格。

假設(shè)某流域一共有n個(gè)柵格，從流域出口反向遞歸索引保證一次模擬只需要進(jìn)行n次計(jì)算，而正向計(jì)算需要不斷地判斷上下游關(guān)系，計(jì)算次數(shù)可能會(huì)達(dá)到n2次甚至n3次。經(jīng)過遞歸算法和哈希表的優(yōu)化，所構(gòu)建分布式模型的運(yùn)行時(shí)間只會(huì)隨流域面積的增大及內(nèi)部方程的增多線性增長(zhǎng)而不是指數(shù)增長(zhǎng)，這為后續(xù)的高效計(jì)算和率定奠定了基礎(chǔ)。

在傳輸過程中，主要考慮水塘與植被過濾帶對(duì)水、沙、污染物的攔截作用。

（1）水塘攔截

水塘的攔沙效率計(jì)算公式為：

式中：trap為攔沙效率，%；Vset為泥沙沉降速率，m·d-1；Vovfl為出流速率，m·d-1；Q0表示水塘的出流量，m3·s-1；SAres表示水塘的水面面積，hm2。

（2）植被過濾帶

本模型基于鄧娜等[20]提出的植被過濾帶凈化效果多元非線性回歸方程，模擬傳輸路徑上植被對(duì)泥沙、徑流和吸附態(tài)磷的攔截效果，計(jì)算公式為：

式中：Itr為植被過濾帶對(duì)泥沙、徑流和吸附態(tài)磷的攔截比例，%；C為植被覆蓋與管理因子，無量綱；Qsurf為地表徑流，mm；conc為入流的泥沙、吸附態(tài)磷的濃度，mg·L-1；θ為坡度，（°）；WB為植被過濾帶寬度，m。

1.2.3 河道傳輸模塊

（1）泥沙沉降/再懸浮模擬

本模型采用MOLINAS 等[21]提出的基于通用河流功率的泥沙搬運(yùn)方程計(jì)算河道搬運(yùn)的最大含沙量：

式中：concsed，ch，mx為河道能搬運(yùn)的最大含沙量，t·m-3；Sg為河道中固體的相對(duì)密度，2.65；Cw為河道含沙量，kg·m-3；ψ為通用河道功率；depth表示河道深度，m；D50表示泥沙的中值粒徑，mm。

大于河道最大含沙量部分的泥沙儲(chǔ)存在當(dāng)前的河道柵格中，直到下次河道含沙量小于最大含沙量時(shí)，發(fā)生再懸浮。

（2）河道磷循環(huán)模擬

在水體中，藻類死亡后藻體磷轉(zhuǎn)化為有機(jī)磷，有機(jī)磷可通過礦化作用轉(zhuǎn)化為能被藻類吸收的可溶性磷，也可通過沉降作用從水流中去除。本模型借鑒SWAT 模型中河道磷循環(huán)過程的相關(guān)方程，某日有機(jī)磷濃度的變化量為：

式中：α2為藻體磷占藻類生物質(zhì)的比例，無量綱；ρa(bǔ)為當(dāng)?shù)卦孱惖暮粑鼜?qiáng)度或死亡率，d-1；algea為某日初始時(shí)刻藻類的生物質(zhì)濃度，mg·L-1；βP，4為有機(jī)磷的成礦速率常數(shù)，d-1；orgPstr為某日初始時(shí)刻有機(jī)磷的濃度，mg·L-1；σ5為有機(jī)磷的沉降速率系數(shù)，d-1；TT為河段的水流運(yùn)動(dòng)時(shí)間，d。

某日可溶性磷的濃度變化量為：

式中：σ2為沉積物提供可溶性磷的速率，mg·m-2·d-1；depth為河道水深，m；μa為當(dāng)?shù)卦孱惖纳L(zhǎng)率，d-1。與藻類相關(guān)的參數(shù)均為待率定參數(shù)。

1.3 模型應(yīng)用

本模型代碼實(shí)現(xiàn)主要基于Python 3.9 編程語(yǔ)言，使用pandas（csv 文件的讀取與操作）和gdal（柵格文件與csv 文件的互相轉(zhuǎn)換）等基礎(chǔ)包實(shí)現(xiàn)計(jì)算。模擬時(shí)將土地利用、土壤類型、徑流流向、坡度、降雨量、產(chǎn)流量、水土流失量等基礎(chǔ)數(shù)據(jù)轉(zhuǎn)換為csv文件，將單元格視作柵格，在Python 腳本中模擬計(jì)算水土流失、徑流與污染物的源-流-匯過程，模型的運(yùn)行步長(zhǎng)為月。

本模型共涉及35 個(gè)參數(shù)，以納什效率系數(shù)（NSE）、決定系數(shù)（R2）和百分比偏差（PBIAS）為優(yōu)化目標(biāo)進(jìn)行模型參數(shù)率定。

式中：Oi為每月觀測(cè)值；Pi為每月模擬值為每月觀測(cè)值的平均值為每月模擬值的平均值。

由于涉及35 個(gè)待率定參數(shù)，以及3 個(gè)率定指標(biāo)，因此模型率定屬于典型NP-HARD 問題，在解決該問題的方法中以帶有精英保留策略的NSGA-Ⅱ算法運(yùn)用最為廣泛[22]。但在實(shí)際應(yīng)用中NSGA-Ⅱ算法的性能會(huì)隨著優(yōu)化目標(biāo)的擴(kuò)展和優(yōu)化問題的復(fù)雜化而受到顯著影響。因此本研究使用2014 年DEBK 提出的NSGA-Ⅲ算法對(duì)模型的參數(shù)進(jìn)行率定。通過引入廣泛分布參考點(diǎn)維持種群多樣性并對(duì)擁擠排序進(jìn)行改編，NSGA-Ⅲ算法具有良好的搜索高維空間中解的能力，在參數(shù)率定中得到了廣泛的應(yīng)用[23]。各參數(shù)設(shè)置如表1所示。

表1 NSGA-Ⅲ算法參數(shù)取值Table 1 Parameters of NSAG-Ⅲalgorithm

分布式模型提供的高精度的污染源空間分布為防控措施的模擬提供了基礎(chǔ)?；诓煌臇鸥癯叨仍磮D層，如水土流失圖層、土壤磷含量圖層、土地利用圖層等，本研究分別評(píng)估了退耕還林、化肥減量、坡改梯、濱岸緩沖帶等措施對(duì)面源污染的削減效率，為水土流失型面源污染的防控提供依據(jù)，提出了5 種模擬情景。情景1：基準(zhǔn)情景，不做任何的修改；情景2：基于流域多年水土流失情況，對(duì)水土流失強(qiáng)度前30%柵格上的坡耕地進(jìn)行退耕還林；情景3：基于土壤磷含量圖層，對(duì)土壤磷含量前30%的坡耕地/水田化肥減量50%；情景4：基于坡度圖層，對(duì)坡度>10°的坡耕地實(shí)施坡改梯措施；情景5：基于河道圖層與污染傳輸圖層，識(shí)別在總磷入河量前30%的濱岸柵格，并設(shè)置30 m濱岸緩沖帶。

2 結(jié)果與討論

2.1 模型可靠性

利用2015年5月至2017年9月的徑流、總磷監(jiān)測(cè)數(shù)據(jù)對(duì)模型進(jìn)行率定，利用2019年12月至2021年12月的徑流、總磷監(jiān)測(cè)數(shù)據(jù)以及2021 年9 月至12 月的吸附態(tài)磷、溶解態(tài)磷監(jiān)測(cè)數(shù)據(jù)對(duì)模型進(jìn)行驗(yàn)證，率定結(jié)果如圖4所示。結(jié)果表明，率定期和驗(yàn)證期內(nèi)徑流的NSE和R2均大于0.5，PBIAS均在±25%范圍內(nèi)，說明模擬效果較好。模擬效率方面，研究區(qū)共劃分為9 635個(gè)30 m×30 m的柵格，模擬一次平均耗時(shí)4.196 s，率定共耗時(shí)6.99 h，模型效率較高。

圖4 模型率定結(jié)果Figure 4 Calibration results of the model

2.2 水土流失型面源污染時(shí)空分布規(guī)律

2.2.1 空間分布

如圖5a 所示，受地形地貌、人類活動(dòng)等的影響，年均水土流失、吸附態(tài)磷和溶解態(tài)磷負(fù)荷量在空間上分布差異明顯。研究區(qū)水土流失量在0~149 t·hm-2·a-1之間，平均流失量為 17.23 t·hm-2·a-1。根據(jù)水利部發(fā)布的《土壤侵蝕分類分級(jí)標(biāo)準(zhǔn)》（SL 190—2007）規(guī)定，丘陵紅壤區(qū)容許的土壤流失量為5 t·hm-2·a-1，研究區(qū)超過此限值的區(qū)域占流域面積的40.2%。吸附態(tài)磷負(fù)荷量在0~1.72 kg·hm-2·a-1之間，平均值為1.22 kg·hm-2·a-1，溶解態(tài)磷負(fù)荷量在 0~1.61 kg·hm-2·a-1之間，平均值為0.56 kg·hm-2·a-1。水土流失、吸附態(tài)磷和溶解態(tài)磷的入河量較負(fù)荷量有所下降，分別為9 032.2 t·a-1、601.3 kg·a-1和 541.7 kg·a-1，分別占負(fù)荷量的54.32%、50.87%和72.99%。水土流失和吸附態(tài)磷的入河系數(shù)隨著入河距離的增大而減小，因此水土流失和吸附態(tài)磷入河量的熱點(diǎn)區(qū)域主要分布于河道兩側(cè)。由于水塘、植被過濾帶對(duì)溶解態(tài)磷的攔截作用有限，隨入河距離的增大溶解態(tài)磷的入河系數(shù)并未顯著減小。

此外，將柵格尺度的數(shù)據(jù)進(jìn)行加和平均，升尺度至集水區(qū)尺度，在集水區(qū)的空間尺度上計(jì)算了水土流失、吸附態(tài)磷、溶解態(tài)磷的負(fù)荷量、入河量和入河系數(shù)（圖5b）。

圖5 水土流失、吸附態(tài)磷、溶解態(tài)磷的產(chǎn)生和入河空間分布情況Figure 5 Spatial distribution of production and inflow of and soil erosion，adsorbed phosphorus and dissolved phosphorus

由于林地和坡耕地多集中于中上游流域，因此水土流失和吸附態(tài)磷的關(guān)鍵源區(qū)多位于中上游的集水區(qū)。溶解態(tài)磷則主要由水田密集的下游集水區(qū)貢獻(xiàn)。中上游流域河網(wǎng)稀疏，入河距離較長(zhǎng)，入河系數(shù)較河網(wǎng)密集的下游流域小。

基于分布式的空間分布可以在柵格尺度上識(shí)別負(fù)荷量與入河量的熱點(diǎn)區(qū)域，進(jìn)而基于柵格尺度的水土流失、面源污染入河量、土壤磷含量等圖層，精準(zhǔn)識(shí)別防控措施的最佳實(shí)施位置，為水土流失型面源污染的防控提供依據(jù)。盡管基于半分布式的建模方法能大致得出水土流失、吸附態(tài)磷、溶解態(tài)磷的分布規(guī)律，并在集水區(qū)尺度上識(shí)別面源污染關(guān)鍵源區(qū)[24-26]，但集水區(qū)尺度的結(jié)果難以進(jìn)一步空間降尺度，因此對(duì)防控措施指導(dǎo)意義較差。對(duì)比分布式與半分布式入河系數(shù)的空間分布，分布式模型揭示出柵格入河系數(shù)與其入河距離呈負(fù)相關(guān)關(guān)系，但半分布式模型無法得出這一規(guī)律，難以為精細(xì)化的關(guān)鍵源區(qū)識(shí)別及防控措施實(shí)施提供依據(jù)。

表2 為不同土地利用類型下吸附態(tài)磷、溶解態(tài)磷和水土流失的負(fù)荷情況。林地、水田、坡耕地的年均水土流失負(fù)荷分別為30.204、9.529 t·hm-2和 15.111 t·hm-2，呈現(xiàn)林地>坡耕地>水田的趨勢(shì)。林地面積占36.37%，水土流失量占比卻高達(dá)62.83%，主要是因?yàn)榱值卮蠖嗉杏诟叱滩钶^大的上游流域，此處坡度因子較大，更容易發(fā)生水土流失。而面積占比為23.37%和28.26%的水田和坡耕地的水土流失量則分別占比為12.74%和24.42%。春夏季水田和坡耕地施肥使土壤含磷量較林地高，因此水田和坡耕地吸附態(tài)磷的占比較其水土流失的占比有所上升。溶解態(tài)磷負(fù)荷主要受土壤含磷量的影響，因此呈現(xiàn)坡耕地>水田>林地的趨勢(shì)。

表2 不同土地利用類型下水土流失、吸附態(tài)磷、溶解態(tài)磷負(fù)荷情況Table 2 Distribution of sediment，adsorbed phosphorus and dissolved phosphorus in different land uses

表3 為不同坡度下吸附態(tài)磷、溶解態(tài)磷和水土流失的負(fù)荷情況。結(jié)果表明，面積占比為27.71%的陡坡地（坡度>20°）水土流失占比達(dá)到了43.80%；相比之下，面積占比為49.13%、坡度<10°的土地水土流失占比卻只有26.04%。由于水田和坡耕地大多分布在坡度<10°的土地上，土壤磷含量較高，因此坡度<10°的土地吸附態(tài)磷占比達(dá)29.88%，溶解態(tài)磷占比達(dá)54.91%。

表3 不同坡度泥沙、吸附態(tài)磷、溶解態(tài)磷負(fù)荷情況Table 3 Distribution of sediment，adsorbed phosphorus and dissolved phosphorus in different slopes

2.2.2 時(shí)間分布

圖6為吸附態(tài)磷和溶解態(tài)磷濃度（圖6a）、入河系數(shù)（圖6b）及其與降雨（圖6c）、侵蝕性降雨（>12 mm，圖6d）相關(guān)性的月際時(shí)間變化。結(jié)果表明，河道吸附態(tài)磷、溶解態(tài)磷及其入河系數(shù)春夏季節(jié)較秋冬季節(jié)高，且與侵蝕性降雨相關(guān)性較強(qiáng)，其中吸附態(tài)磷與侵蝕性降雨的相關(guān)性要略高于溶解態(tài)磷，在4 月至8 月間，吸附態(tài)磷與侵蝕性降雨的相關(guān)性在0.4以上，溶解態(tài)磷與侵蝕性降雨的相關(guān)性在0.3 以上。2 月和3 月水土流失量、流量與侵蝕性降雨的相關(guān)性較高，但是吸附態(tài)磷、溶解態(tài)磷與侵蝕性降雨的相關(guān)性很低，主要是因?yàn)檩剂窒饔虻母N、施肥從4月開始，2月和3 月土壤磷含量較低，導(dǎo)致了磷與侵蝕性降雨之間的低相關(guān)性。此外，溶解態(tài)磷與吸附態(tài)磷入河系數(shù)的變化趨勢(shì)與侵蝕性降雨變化趨勢(shì)基本一致，相關(guān)性在0.3~0.8之間。溶解態(tài)磷的入河系數(shù)在0.6~0.9之間波動(dòng)，而吸附態(tài)磷的入河系數(shù)在0.2~0.6 之間波動(dòng)，這是因?yàn)樗梁椭脖坏臄r截能力有限，在面源污染量較少（即侵蝕性降雨量低）時(shí)能攔截絕大部分的傳輸量，而當(dāng)面源污染較多（即侵蝕性降雨量高）時(shí)只能攔截部分傳輸量，因此侵蝕性降雨量低時(shí)吸附態(tài)磷和溶解態(tài)磷的入河系數(shù)較低，而侵蝕性降雨量高時(shí)二者的入河系數(shù)較高。吸附態(tài)磷、溶解態(tài)磷及其入河系數(shù)與降雨量的相關(guān)性較侵蝕性降雨量的相關(guān)性低，這主要是因?yàn)榻涤炅枯^低時(shí)，幾乎不產(chǎn)生地表徑流與水土流失，此時(shí)水土流失與非點(diǎn)源磷和降雨量之間無顯著關(guān)系。

圖6 吸附態(tài)磷、溶解態(tài)磷濃度和入河系數(shù)的變化情況以及各指標(biāo)與非侵蝕性和侵蝕性降雨的相關(guān)性分析Figure 6 Concentration and introduced ration changes of adsorbed phosphorus and dissolved phosphorus and correlation analysis between each index and non-erosive/erosive rainfall

2.3 防控措施效果評(píng)估

情景2 至情景4 措施實(shí)施的空間位置如圖7 所示，5 種情景對(duì)吸附態(tài)磷、溶解態(tài)磷和水土流失量的削減情況如表4 所示。在情景2 中，將水土流失強(qiáng)度前30% 的坡耕地進(jìn)行退耕還林（占流域面積的8.48%），由于坡耕地水土流失占比僅為24.42%，退耕還林措施僅削減了4.78%的水土流失量。而退耕還林措施對(duì)面源污染的治理效果較好，分別削減了9.32%的吸附態(tài)磷和6.01%的溶解態(tài)磷。在情景3中，對(duì)土壤磷含量前30%的坡耕地/水田（占流域面積的15.49%）化肥減量50%，分別有17.97%的吸附態(tài)磷和20.81%的溶解態(tài)磷被削減。在情景4 中，對(duì)坡度大于10°的坡耕地（占流域面積的13.61%）進(jìn)行坡改梯，梯田坡度較坡耕地小，對(duì)水土流失、徑流、污染物的截留效果更好，分別削減了8.08%的水土流失量、9.03%的吸附態(tài)磷和2.37%的溶解態(tài)磷。情景5 中，在總磷入河量前30%的濱岸柵格設(shè)置了30 m 濱岸緩沖帶，經(jīng)過上游節(jié)點(diǎn)的層層傳輸，濱岸入河?xùn)鸥竦乃亮魇Ъ懊嬖次廴镜牧孔畲螅诖颂幉荚O(shè)緩沖帶對(duì)水土流失及面源污染的削減效果最好。結(jié)果表明，情景5 能分別削減26.22%、22.97%、25.01%的吸附態(tài)磷、溶解態(tài)磷及水土流失。

圖7 措施實(shí)施位置識(shí)別Figure 7 Identification of implementation position of measures

表4 措施削減效果Table 4 Reducing effect of measures

分布式模擬的思路解決了傳統(tǒng)水保措施效果評(píng)估方法的諸多不足。LI 等[27]基于SWAT 模型，發(fā)現(xiàn)在坡度>25°的坡耕地上實(shí)施退耕還林能降低59.53%的總磷負(fù)荷，但是基于半分布式模型進(jìn)行效果評(píng)估僅能模擬每個(gè)集水區(qū)的削減效果，并且忽略了水土流失和污染物在傳輸過程中的損失。水土流失和污染物在傳輸過程中由于自然沉降和植被攔截會(huì)遺留在田塊上，在之后的降雨事件中會(huì)被二次沖刷[28]，因此半分布式模型可能會(huì)對(duì)退耕還林的減污效果有所高估，并且難以精準(zhǔn)定位退耕還林的最佳實(shí)施位置。分布式模型則能模擬水塘、植被緩沖帶等地表景觀對(duì)水土流失、污染物的截留，以最終入河量削減為評(píng)估目標(biāo)，更加真實(shí)可靠。另外，分布式模型能在柵格尺度上模擬水土流失和污染物的負(fù)荷量、入河量、土壤磷含量等信息，為識(shí)別退耕還林、化肥減量、坡改梯、濱岸緩沖帶等措施的最佳實(shí)施位置提供依據(jù)。例如在設(shè)置濱岸緩沖帶時(shí)，半分布式模型無法識(shí)別面源污染入河量最高的濱岸柵格，在所有濱岸柵格設(shè)置濱岸緩沖帶則成本過高，隨機(jī)設(shè)置濱岸緩沖帶則效果不如設(shè)置在關(guān)鍵濱岸柵格，因此分布式模型比半分布式模型對(duì)管理措施的實(shí)施更具指導(dǎo)意義。

3 結(jié)論

（1）本研究針對(duì)半分布式模型的局限性，考慮了降雨徑流、土壤侵蝕、田塊流失、水塘攔截、植被過濾、泥沙沉降、河道磷循環(huán)等分布式水土流失型面源污染關(guān)鍵過程，構(gòu)建了分布式水土流失型面源污染模型，模型模擬精度較高；并基于反向遞歸與哈希表提出了高效匯流算法，破解了分布式模擬面臨的計(jì)算瓶頸問題。

（2）三峽庫(kù)區(qū)菁林溪流域內(nèi)水土流失、吸附態(tài)磷和溶解態(tài)磷的負(fù)荷強(qiáng)度分別為17.23 t·hm-2、1.22 kg·hm-2和 0.56 kg·hm-2，其中分別有 54.32%、50.87% 和72.99%的水土流失、吸附態(tài)磷和溶解態(tài)磷最終進(jìn)入河道。

（3）丘陵山區(qū)水土流失面源防控應(yīng)重點(diǎn)關(guān)注坡耕地影響，林地、水田、坡耕地的吸附態(tài)磷流失強(qiáng)度分別為 1.911、1.179、2.192 kg·hm-2，溶解態(tài)磷流失強(qiáng)度分別為0.604、1.084、1.086 kg·hm-2，面積占比28.26%的坡耕地即貢獻(xiàn)了全流域38.91%的溶解態(tài)磷和39.45%的吸附態(tài)磷。

（4）溶解態(tài)磷的入河系數(shù)在0.6~0.9 之間波動(dòng)，吸附態(tài)磷的入河系數(shù)在0.2~0.6 之間波動(dòng)，并且二者與降雨量之間呈現(xiàn)較強(qiáng)的相關(guān)性。

（5）相對(duì)于化肥減量等經(jīng)典農(nóng)藝措施，在入河負(fù)荷最大的濱岸柵格上設(shè)置濱岸緩沖帶對(duì)丘陵山區(qū)面源防控效果更佳，其能分別減少26.22%的吸附態(tài)磷和22.97%的溶解態(tài)磷的入河量。