李情，唐代生，孟岑，李裕元，吳金水

（1.中南林業(yè)科技大學(xué)林學(xué)院，長(zhǎng)沙 410004；2.中國(guó)科學(xué)院亞熱帶農(nóng)業(yè)生態(tài)研究所亞熱帶農(nóng)業(yè)生態(tài)過(guò)程重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙410125；3.中國(guó)科學(xué)院亞熱帶農(nóng)業(yè)生態(tài)研究所長(zhǎng)沙農(nóng)業(yè)環(huán)境觀測(cè)研究站，長(zhǎng)沙 410125）

磷作為水生生態(tài)系統(tǒng)富營(yíng)養(yǎng)化的限制性營(yíng)養(yǎng)元素，其主要來(lái)源包括點(diǎn)源和面源，相較于點(diǎn)源，面源磷污染的控制和管理仍存在巨大的挑戰(zhàn)[1]，主要因?yàn)槠渚哂锌臻g分布的廣泛性和不確定性等特點(diǎn)[2]。面源磷負(fù)荷的精準(zhǔn)核算，關(guān)鍵源區(qū)和空間分布規(guī)律的明晰，是實(shí)現(xiàn)基于流域面源污染高效控制的重要前提。

面源磷污染的隨機(jī)性、遷移和轉(zhuǎn)化的復(fù)雜性和滯后性決定了其負(fù)荷估算和控制難度大，因此通過(guò)模型模擬是量化污染負(fù)荷、識(shí)別關(guān)鍵源區(qū)的重要手段[3]。近年來(lái)我國(guó)面源磷負(fù)荷估算模型在構(gòu)建與應(yīng)用方面取得了長(zhǎng)足的進(jìn)展。根據(jù)參數(shù)特征，模型可分為集總式和分布式。分布式模型如AGNPS、APEX、HSPF、SWAT 等應(yīng)用廣泛，分布式模型考慮了污染物遷移過(guò)程，涉及多介質(zhì)、多尺度的污染物產(chǎn)匯模擬和農(nóng)業(yè)管理等措施評(píng)估[4]。這些模型多建立在大量基礎(chǔ)資料和物理機(jī)制的基礎(chǔ)上，對(duì)于特定區(qū)域環(huán)境能夠?qū)崿F(xiàn)精確模擬，但卻限制了其在區(qū)域外的應(yīng)用和推廣[5]。通用土壤流失模型作為國(guó)內(nèi)外廣泛應(yīng)用的土壤侵蝕模型，具有結(jié)構(gòu)簡(jiǎn)潔、參數(shù)物理意義明確、計(jì)算簡(jiǎn)單等優(yōu)勢(shì)[6]，但其僅能用于吸附態(tài)磷的模擬估算。近年來(lái)有研究將其與輸出系數(shù)等集總式模型進(jìn)行耦合模擬流域總磷負(fù)荷量及空間分布[7-8]，但多數(shù)研究未對(duì)模型模擬效率進(jìn)行驗(yàn)證，且模型缺乏參數(shù)本地化校準(zhǔn)。

我國(guó)長(zhǎng)江中下游亞熱帶丘陵區(qū)是我國(guó)重要的農(nóng)業(yè)糧食產(chǎn)區(qū)，農(nóng)業(yè)集約化程度高，水系分布密集，面源磷污染給當(dāng)?shù)厣鷳B(tài)環(huán)境帶來(lái)了嚴(yán)重的負(fù)面影響[9]。同時(shí)該區(qū)域具有相似的景觀格局[10]，相比中大尺度區(qū)域，選取典型農(nóng)林流域?yàn)閷?duì)象構(gòu)建磷負(fù)荷估算模型具有較好的可操作性和推廣性。本研究選取長(zhǎng)沙縣金井河流域?yàn)檠芯繀^(qū)域，基于改進(jìn)的輸出系數(shù)和RULSE 模型構(gòu)建流域磷負(fù)荷模型，對(duì)面源磷負(fù)荷進(jìn)行模擬估算，解析不同形態(tài)磷負(fù)荷的空間分布特征，以期為亞熱帶農(nóng)業(yè)流域面源磷負(fù)荷的核算及控制管理提供有效手段和科學(xué)依據(jù)。

1 材料與方法

1.1 研究區(qū)域概況

金井河流域位于湖南省長(zhǎng)沙縣金井鎮(zhèn)（27°55′~28°40′N(xiāo)、112°56′~113°30′E），湘江一級(jí)支流撈刀河上游區(qū)域，流域總面積為134.92 km2。流域內(nèi)以低山丘陵為主，中北部地勢(shì)高，南部地勢(shì)低，海拔在56~440 m 之間，平均坡度為18°，坡度大于25°的面積約占流域面積的1/3。該地為典型南方紅壤丘陵地貌，土壤主要類(lèi)型為紅壤和水稻土，其氣候?qū)儆诘湫蛠啛釒駶?rùn)季風(fēng)氣候，年平均降水量為1 200~1 500 mm。流域水利設(shè)施以水庫(kù)為主，其中有1 座Ⅰ中型水庫(kù)（金井水庫(kù)，位于流域東南部），2 座?、裥退畮?kù)，15 座小Ⅱ型水庫(kù)，2 422 口山塘，總蓄水容量為2 150.38 萬(wàn)m3。土地利用方式以林地、農(nóng)田、居民地為主，其面積占比分別為65.1%、30.2%和1.4%。林地主要分布于丘陵崗地等海拔坡度較高的區(qū)域，主要樹(shù)種有馬尾松（Pinus massoniana）、杉木（Cunninghamia lanceolata）和油茶（Camellia oleifera）等。農(nóng)田分布在南部河流沖擊平原及丘陵溝谷之間，以水稻田為主，平均施磷量為 6 500 kg·km-2·a-1，并輔以有機(jī)肥。流域內(nèi)無(wú)中大型城鎮(zhèn)和工廠等大型點(diǎn)源污染。金井河流域研究區(qū)域如圖1所示。

1.2 模型框架構(gòu)建

1.2.1 模型數(shù)據(jù)庫(kù)

流域內(nèi)設(shè)有1個(gè)標(biāo)準(zhǔn)氣象觀測(cè)站和3個(gè)小型氣象觀測(cè)站。流域出口設(shè)有長(zhǎng)期定位水文監(jiān)測(cè)點(diǎn)（2014—2020年），選取2014—2019年觀測(cè)的溶解態(tài)磷負(fù)荷用于模型參數(shù)優(yōu)化。利用ArcGIS 軟件的水文分析模塊和D8 流向法將金井河流域劃分為20 個(gè)集水區(qū)，選取其中部分集水區(qū)建立水文水質(zhì)監(jiān)測(cè)點(diǎn)（圖1），并將2020 年觀測(cè)得到的溶解態(tài)磷負(fù)荷與總磷負(fù)荷用于模型效果驗(yàn)證。本研究所需其他基礎(chǔ)資料包括遙感影像數(shù)據(jù)、DEM 數(shù)據(jù)（分辨率5 m×5 m）、土地利用數(shù)據(jù)、土壤數(shù)據(jù)、降雨數(shù)據(jù)等的來(lái)源見(jiàn)表1。

表1 金井河流域面源磷污染負(fù)荷模型數(shù)據(jù)來(lái)源Table 1 Data sources of nonpoint source phosphorus pollution load model in Jinjing watershed

圖1 金井河流域研究范圍Figure 1 Study area of Jinjing watershed

1.2.2 模型構(gòu)建

面源磷負(fù)荷可分為溶解態(tài)磷負(fù)荷和吸附態(tài)磷負(fù)荷。由于土壤侵蝕是驅(qū)動(dòng)吸附態(tài)磷輸出的主要因素，因此基于RUSLE 實(shí)現(xiàn)模型對(duì)吸附態(tài)磷的模擬。由于磷在河道遷移過(guò)程中存在衰減，因此將傳統(tǒng)輸出系數(shù)模型與一階河道衰減系數(shù)耦合構(gòu)建溶解態(tài)磷模擬模塊。模型框架示意見(jiàn)圖2。

圖2 金井河流域面源磷負(fù)荷模型框架Figure 2 Modeling framework of nonpoint source phosphorus load model in Jinjing watershed

1.3 溶解態(tài)磷負(fù)荷模型構(gòu)建

本研究基于不同土地利用類(lèi)型（林地、農(nóng)田、居民地）輸出系數(shù)及一階河道衰減系數(shù)構(gòu)建溶解態(tài)磷負(fù)荷模型，公式如下：

式中：LSRP為溶解態(tài)磷負(fù)荷，kg·a-1；ei、mi分別為第i種土地利用溶解態(tài)磷輸出系數(shù)和面積，kg·km-2·a-1和km2；k為溶解態(tài)磷河流衰減系數(shù)，d?1；t為溶解態(tài)磷從集水區(qū)河道到流域出水口的輸移時(shí)間，由河道長(zhǎng)度與平均流速計(jì)算得到，d；Lj為第j個(gè)支流或河道上游入口的溶解態(tài)磷負(fù)荷，kg·a-1；tj為第j個(gè)支流或上游流入的溶解態(tài)磷從對(duì)應(yīng)河道起點(diǎn)到集水區(qū)出口的輸移時(shí)間，d。

此外對(duì)于活塞流遠(yuǎn)大于分散流的中小流域，一般假設(shè)流域面源污染物沿集水區(qū)內(nèi)的給定河段進(jìn)入，平均損失超過(guò)河流長(zhǎng)度的1/2（相當(dāng)于t/2）[11]。本研究?jī)H對(duì)瞬時(shí)流速進(jìn)行監(jiān)測(cè)，因此基于已有研究中瞬時(shí)流速與河道平均流速的經(jīng)驗(yàn)公式計(jì)算得到平均流速[12]，計(jì)算公式為：

式中：d為流程，m；v為平均流速，m·s-1；Q為集水區(qū)出口河流檢測(cè)瞬時(shí)流量，m3·s-1；Q1為河流年均流量，m3·s-1；s為集水區(qū)面積，m2；s1為無(wú)量綱的積水區(qū)面積；Q2為無(wú)量綱的流量；g為重力加速度，9.8 m·s-2。

采用貝葉斯-蒙特卡洛方法（Markov Chain Monte Carlo，MCMC）對(duì)公式（1）中不同土地利用磷輸出系數(shù)和河道磷衰減系數(shù)后驗(yàn)分布進(jìn)行校準(zhǔn)計(jì)算[13]。具體方法如下：假設(shè)Y為響應(yīng)變量，即流域溶解態(tài)磷負(fù)荷的自然對(duì)數(shù)。設(shè)Yh（h=1，…，l）服從正態(tài)分布，其均值為μh（對(duì)應(yīng)每一次采樣流域日溶解態(tài)磷負(fù)荷Yh）、方差為σ2（對(duì)應(yīng)所有Yh）：

溶解態(tài)磷瞬時(shí)負(fù)荷的對(duì)數(shù)服從正態(tài)分布。均值μh是一個(gè)空間過(guò)程，通過(guò)模型得到；方差σ2是一個(gè)常數(shù)，采用觀測(cè)數(shù)據(jù)得到。

通過(guò)前期不同土地利用溶解態(tài)磷負(fù)荷和磷的河道衰減系數(shù)獲取模型參數(shù)先驗(yàn)分布，設(shè)模型參數(shù)為相互獨(dú)立的先驗(yàn)分布：

由于模型方差沒(méi)有先驗(yàn)信息，故對(duì)所有模型參數(shù)使用模糊先驗(yàn)分布，即對(duì)輸出系數(shù)（ei）先驗(yàn)分布取自然對(duì)數(shù)，對(duì)河道衰減系數(shù)（k）先驗(yàn)分布開(kāi)方。為模型參數(shù)指定正態(tài)分布，并為參數(shù)指定逆gamma（IG）[14]：

采用 Gibbs sampling 算法實(shí)現(xiàn) MCMC 計(jì)算[15]，利用已知數(shù)據(jù)（LSRP、mi、t）計(jì)算獲得溶解態(tài)磷負(fù)荷模型中e1（林地）、e2（農(nóng)田）、e3（居民地）和k的后驗(yàn)分布，并進(jìn)行校正后輸出模型參數(shù)結(jié)果[13，16]。

1.4 吸附態(tài)磷負(fù)荷模型構(gòu)建

吸附態(tài)磷負(fù)荷主要來(lái)源于土壤侵蝕流失，通用土壤流失方程（USLE）和RUSLE 是估算土壤流失的最經(jīng)典模型，本研究以RUSLE 與污染物負(fù)荷方程構(gòu)建吸附態(tài)磷負(fù)荷模型[17]，計(jì)算公式如下：

式中：LXFP為吸附態(tài)磷負(fù)荷，kg·a-1；λ為泥沙輸移比；XP為土壤侵蝕模數(shù)，t·km-2·a-1；CP為磷污染物背景含量，g·kg-1；η為磷污染物富集比；mi為第i種土地利用面積，km2；R為降雨侵蝕力因子，MJ·mm·hm-2·h-1·a-1；K為土壤可蝕性因子，t·hm2·h·hm-2·MJ-1·mm-1；LS為坡長(zhǎng)坡度因子，無(wú)量綱；C為植被覆蓋因子，無(wú)量綱；P為水土保持措施因子，無(wú)量綱。

1.4.1 泥沙輸移比（λ）

泥沙輸移比反映泥沙從土壤侵蝕到進(jìn)入河道某攔截面處的減少程度，定義為流域出口某斷面輸沙量與斷面以上土壤侵蝕產(chǎn)沙量之比，與流域面積、地形、土壤性質(zhì)及降雨徑流密切相關(guān)。計(jì)算方法包括直接計(jì)算和模型估算，參照長(zhǎng)江流域和洞庭湖流域已有的研究結(jié)果確定本研究流域內(nèi)的泥沙輸移比為0.4[18-19]。

1.4.2 磷污染物背景含量（CP）

磷污染物背景含量指土壤磷污染物含量，即土壤中總磷的含量。根據(jù)流域內(nèi)土壤采樣分析，得到流域內(nèi)土壤磷含量均值為0.51 g·kg-1。

1.4.3 磷污染物富集比（η）

磷污染物富集比是土壤侵蝕產(chǎn)生的泥沙磷含量與雨前土壤表層磷含量之比。以往研究發(fā)現(xiàn)泥沙中磷含量常高于土壤中的磷含量，富集比在1~4 之間，且多在2 附近變動(dòng)，故本研究的磷污染物富集比定為2[20-22]。

1.4.4 降雨侵蝕力因子（R）

降雨侵蝕力因子是土壤侵蝕的動(dòng)力因子，表征降雨引起的潛在動(dòng)能，與雨量、雨強(qiáng)、降雨歷時(shí)有關(guān)，本研究選用WISCHEMIER等[23]的研究中的公式計(jì)算：

式中：P為多年平均降雨量，mm；Pi為第i月多年平均降雨量，mm。通過(guò)流域內(nèi)4個(gè)氣象觀測(cè)站2014—2020 年的降雨數(shù)據(jù)得到平均月降雨量和年降雨量，利用Arc?GIS軟件采用泰森多邊形進(jìn)行空間插值繪制R因子空間分布（圖3）。

1.4.5 土壤可蝕性因子（K）

土壤可蝕性因子表示降雨侵蝕力作用下土壤被沖離的難易程度，用土壤侵蝕速率來(lái)衡量各種土壤的侵蝕敏感度，該因子與土壤性質(zhì)密切相關(guān)。本研究選用WILLIAMS等[24]提出的EPIC模型計(jì)算：

式中：SAN為土壤砂粒含量，%；SIL為土壤粉粒含量，%；CLA為土壤黏粒含量，%；C為土壤有機(jī)碳含量，%；SOM為土壤有機(jī)質(zhì)含量，%；SN為土壤粉粒和黏粒含量總和，%。通過(guò)采樣得到本研究區(qū)域土壤砂粒、粉粒、黏粒、有機(jī)質(zhì)含量，并利用ArcGIS 軟件繪制K因子空間分布（圖3）。

1.4.6 坡長(zhǎng)坡度因子（LS）

坡長(zhǎng)坡度因子可根據(jù)地形地貌直接反映坡長(zhǎng)因子與坡度因子對(duì)土壤侵蝕的影響。本研究利用Arc?GIS水文分析功能，根據(jù)坡長(zhǎng)坡度公式計(jì)算得到[25]：

式中：λ是由 DEM 提取出的水平投影坡長(zhǎng)，m；θ是由DEM 提取出的坡度，（°）。本研究區(qū)域LS因子空間分布如圖3所示。

1.4.7 植被覆蓋因子（C）

植被覆蓋因子用于反映植被覆蓋度和管理水平對(duì)土壤流失量的作用關(guān)系。本研究采用蔡崇法等[26]提出的模型方法計(jì)算：

式中：NDVI為歸一化的植被指數(shù)，通過(guò)遙感影像波譜計(jì)算得到。利用ArcGIS軟件繪制的C因子空間分布見(jiàn)圖3。

圖3 金井河流域R、K、LS和C因子空間分布Figure 3 Spatial distribution of R，K，LS and C factor in Jinjing watershed

1.4.8 水土保持措施因子（P）

水土保持措施因子是采用水土保持措施下與直接順坡耕種下土壤流失量之比。本研究采用賦值法，綜合參照已有研究成果[27-28]，將林地、農(nóng)田、居民地、水域的P因子分別賦值為1、0.35、0.01、0。

2 結(jié)果與討論

2.1 流域磷負(fù)荷驗(yàn)證

基于改進(jìn)的輸出系數(shù)模型和RUSLE 方程分別模擬得到流域不同土地利用溶解態(tài)磷輸出負(fù)荷以及空間尺度上的吸附態(tài)磷負(fù)荷，并進(jìn)一步擬合得到流域面源磷輸出負(fù)荷。為了對(duì)模型的模擬效果進(jìn)行驗(yàn)證，分別將總磷、溶解態(tài)磷、吸附態(tài)磷負(fù)荷模型模擬值與2020 年觀測(cè)值進(jìn)行對(duì)比，結(jié)果表明總磷模擬值與觀測(cè)值的誤差絕對(duì)值為3.3%~27.1%（表2）。溶解態(tài)磷負(fù)荷決定系數(shù)（R2）和納什效率系數(shù)（NSE）分別為0.84和0.82，誤差絕對(duì)值為3.9%~24.7%。顆粒態(tài)磷負(fù)荷R2和NSE分別為 0.68 和 0.65，誤差絕對(duì) 值為 8.9%~30.0%（表2，圖4）。模型模擬結(jié)果具有一定精度，可用于亞熱帶典型農(nóng)林流域面源磷污染的空間估算，且溶解態(tài)磷模擬效果優(yōu)于吸附態(tài)磷。

圖4 溶解態(tài)磷和吸附態(tài)磷負(fù)荷觀測(cè)值與模擬值比較Figure 4 Observed and the modeled of dissolved and adsorbed phosphorus load

表2 金井河流域面源磷負(fù)荷估算與驗(yàn)證Table 2 Estimation and validation of nonpoint source phosphorus load in Jinjing watershed

2.2 流域磷負(fù)荷空間分布

模擬結(jié)果表明，研究區(qū)溶解態(tài)磷和吸附態(tài)磷負(fù)荷的空間分布存在較大的差異（圖5），其中溶解態(tài)磷負(fù)荷平均為29.6 kg·km-2·a-1（12.5~37.4 kg·km-2·a-1），高值區(qū)多分布在中部和西南部的農(nóng)田和居民地聚集區(qū)，而東南部和北部區(qū)域相對(duì)較低。吸附態(tài)磷負(fù)荷平均為34.7 kg·km-2·a-1（20.5~59.4 kg·km-2·a-1），高值區(qū)集中在北部多林地區(qū)域。研究表明不同形態(tài)面源磷負(fù)荷污染來(lái)源和輸移機(jī)制存在差異，溶解態(tài)磷主要來(lái)源于化肥、生活污水和畜禽糞污，而吸附態(tài)磷主要源自高流量水文過(guò)程下的土壤侵蝕。

圖5 金井河流域溶解態(tài)磷和吸附態(tài)磷負(fù)荷集水區(qū)空間分布Figure 5 Sub-watershed spatial distribution of dissolved and adsorbed phosphorus load in Jinjing watershed

2.3 流域磷負(fù)荷源解析

流域面源磷負(fù)荷年輸出總量為7 310.3 kg·a-1，其中溶解態(tài)磷和吸附態(tài)磷負(fù)荷分別占39.0%和61.0%（圖6）。不同土地利用輸出負(fù)荷總量為林地（5 004.6 kg·a-1）>農(nóng)田（2 170.8 kg·a-1）>居民地（134.9 kg·a-1）。林地溶解態(tài)磷和吸附態(tài)磷負(fù)荷總量分別為1 186.3 kg·a-1和 3 818.3 kg·a-1，農(nóng)田溶解態(tài)磷和吸附態(tài)磷負(fù)荷總量分別為1 544.2、626.6 kg·a-1，居民地溶解態(tài)磷和吸附態(tài)磷負(fù)荷總量分別為157.8、17.2 kg·a-1。

圖6 金井河流域不同土地利用類(lèi)型面源磷輸出情況Figure 6 Nonpoint source phosphorus exports of different land use types in Jinjing watershed

不同土地利用面源磷輸出負(fù)荷分別為林地56.5 kg·km-2·a-1、農(nóng) 田 56.9 kg·km-2·a-1、居民 地 35.6 kg·km-2·a-1。林地吸附態(tài)磷負(fù)荷（43.1 kg·km-2·a-1）顯著高于溶解態(tài)磷負(fù)荷（13.4 kg·km-2·a-1，P<0.01）。農(nóng)田溶解態(tài)磷負(fù)荷（40.5 kg·km-2·a-1）顯著高于吸附態(tài)磷負(fù)荷（16.5 kg·km-2·a-1）。居民地溶解態(tài)磷負(fù)荷（31.1 kg·km-2·a-1）顯著高于吸附態(tài)磷負(fù)荷（4.5 kg·km-2·a-1）。研究流域作為典型的亞熱帶丘陵區(qū)的典型農(nóng)林流域，其農(nóng)業(yè)生產(chǎn)活動(dòng)強(qiáng)度大、施肥量高是造成農(nóng)田溶解態(tài)磷負(fù)荷高的主要原因。此外居民地普遍缺乏污水收集和處理設(shè)施，地表硬化率較高，溶解態(tài)磷易隨地表徑流進(jìn)入河道。林地溶解態(tài)磷負(fù)荷雖然較低，但其吸附態(tài)磷負(fù)荷高，主要由于林地多在流域內(nèi)的山崗高地，坡度大、土壤侵蝕量大。

3 結(jié)論

（1）基于相關(guān)評(píng)價(jià)指標(biāo)所構(gòu)建的面源磷負(fù)荷估算模型模擬具有一定精度?；诟倪M(jìn)的輸出系數(shù)模型構(gòu)建溶解態(tài)磷負(fù)荷模型，并采用貝葉斯方法求解模型參數(shù)能夠有效解決模型參數(shù)本地化問(wèn)題；基于RULSE 模型構(gòu)建的吸附態(tài)磷負(fù)荷模型能夠較為準(zhǔn)確地模擬亞熱帶丘陵區(qū)土壤侵蝕狀況，但受降雨監(jiān)測(cè)頻率及模型參數(shù)限制，時(shí)間精度有待進(jìn)一步提升。

（2）金井河流域面源總磷負(fù)荷平均為64.3 kg·km-2·a-1，其中溶解態(tài)磷負(fù)荷為29.6 kg·km-2·a-1，吸附態(tài)磷負(fù)荷為34.7 kg·km-2·a-1。流域溶解態(tài)磷負(fù)荷高值區(qū)集中在西南部平原區(qū)，吸附態(tài)磷負(fù)荷高值區(qū)主要分布在西北部林坡地區(qū)域。

（3）農(nóng)田、居民地溶解態(tài)磷負(fù)荷顯著高于林地，平均輸出負(fù)荷分別為 40.5、31.1、13.4 kg·km-2·a-1；林地吸附態(tài)磷負(fù)荷顯著高于農(nóng)田和居民地，負(fù)荷均值分別為 43.1、16.5、4.5 kg·km-2·a-1。亞熱帶丘陵區(qū)面源磷負(fù)荷的治理需要根據(jù)不同形態(tài)磷負(fù)荷輸出空間分布差異有針對(duì)性地開(kāi)展管控措施。