阿什河流域氮磷輸出負荷及其與下墊面特征的關系研究

高曉曦，左德鵬，馬廣文，徐宗學，李佩君

(1.北京師范大學水科學研究院 城市水循環與海綿城市技術北京市重點實驗室，北京 100875;2.中國環境監測總站 國家環境保護環境監測質量控制重點實驗室，北京100012)

水是人類賴以生存的自然資源，水質好壞對人類生產生活至關重要。隨著工業和農業的快速發展，對水資源的需求日益增加，而工業廢水、生活污水亂排亂放，農藥化肥肆意使用，使不同地區水資源及水生態環境均遭受到不同程度的污染和破壞，環境問題已經成為全社會關注的焦點，其中水環境問題尤其突出。

為解決這一問題，世界上許多國家最初都在點源污染防治上投入了大量人力和物力。但把水環境安全僅僅寄托在點源污染防治上，并不能從根本上解決問題[1,2]。越來越多國家逐漸開始關心非點源污染在水環境污染方面所扮演的角色。美國環保局2003年調查結果顯示，農業非點源污染是美國河流和湖泊污染的最大污染源，是造成40%河流和湖泊水體水質變差、地下水體污染和濕地退化的主要原因[3]；歐洲國家農業非點源污染也是造成水體，尤其是地下水污染和地表水磷富集的主要原因，其中由農業非點源排放所造成污染占地表水污染總負荷的24%～71%[4]；我國大量研究結果同樣表明[5,6]，農業非點源污染在非點源污染中占決定性作用。土壤中的氮磷等營養物質在降水沖刷和泥沙攜帶作用下進入河道，造成水生態環境惡化[7,8]。農業非點源污染形成過程隨機性大、影響因子復雜多樣，其分布范圍廣、潛伏周期長等特點導致研究和控制非點源污染難度較大[5,9]。

阿什河是松花江干流南岸支流，位于黑龍江省南部。近些年來，隨著該地區經濟快速增長和糧食產量不斷增加，流域水體污染愈發嚴重，已經被《松花江流域水污染防治“十二五”規劃》列為“哈爾濱市污染最為嚴重的河流”[10]。隨著阿什河流域點源污染治理的穩步推進，非點源污染所引發的問題逐漸顯現。阿什河流域是我國重要的糧食生產基地之一，化肥農藥用量高，施用方式不合理，造成河水污染負荷較高。流域面臨的農業生產壓力和非點源污染風險，使得開展非點源污染研究迫在眉睫。

基于以上背景，本研究通過在阿什河流域構建SWAT模型，在水文相應單元尺度上分析了該地區的總磷、總氮的空間分布特征，重點分析了不同土地利用類型、不同土壤類型、不同坡度等下墊面特征上氮磷輸出負荷的變化規律，以期為阿什河流域非點源污染控制和土地管理提供科學依據。

1 研究區概況

阿什河位于黑龍江省哈爾濱市轄區內(126°40′E～127°43′E，45°05′N～45°50′N)，屬于松花江上游右岸一級支流。阿什河發源于尚志蜜蜂山，流經尚志、五常、阿城和哈爾濱等區縣，于哈爾濱市東郊注入松花江[11]。流域面積3 545 km2，河流全長257 km。流域內地形起伏不大，上游以山地為主，中下游主要是平原和丘陵。阿什河流域屬大陸性季風氣候，具體表現為春季干燥多風，夏季降水集中，秋季降溫迅速，冬季干旱漫長[12]。流域年內降水量分布不均，多年平均降水量565 mm，但6-9月降水占全年降水量的77%左右[13]。阿什河屬山溪性河流，主要以降水補給為主，融雪和地下水補給為輔。流域內徑流流量與地勢高程呈正相關，上游地勢高，降水量量大，徑流量相應較大；而下游地勢較低，降雨較少，使得徑流量也相對較少[14]。阿什河流域是我國重要的糧食生產基地，流域內主要農作物有水稻、玉米和大豆，而農業生產大量施用化肥、農藥所造成的非點源污染，導致阿什河阿城以下河段水質污染嚴重[15]。研究區地理位置如圖1所示。

圖1 阿什河流域位置示意圖Fig.1 The location of the Ashi River basin

2 數據和方法

2.1 數據來源

研究所采用的基礎數據包括地理空間數據云提供的90 m DEM高程數據、中國科學院資源環境科學數據中心提供的1∶10萬土地利用柵格數據、1∶100萬土壤矢量數據集以及1∶25萬河網圖，采用的這些基礎地理信息數據已經過嚴格的質量檢查，數據質量較為可靠，且這些基礎數據已被廣泛應用于水文模型建模中[16-18]，具有一定的代表性。采用SWAT模型中國大氣同化驅動數據集(CMADS V1.0)以及中國氣象數據共享網氣象資料作為模型的氣象驅動數據。水文資料及水質監測數據分別來源于水利部阿城水文站和哈爾濱市環境監測中心，時間序列長度為2010-2015年。研究區農業用地的化肥、農藥施用量是根據哈爾濱統計年鑒以及根據實際情況折算得到。由于阿什河產沙量較低，該流域沒有進行泥沙監測；選取的水質指標總氮總磷只有8個月的例行監測，監測時間為1月、2月及5-10月。

2.2 研究方法

2.2.1 SWAT模型介紹

SWAT模型是由美國農業研究部農業研究中心開發、具有很強物理機制的分布式水文模型。模型根據流域高程將流域劃分為若干子流域，每一個子流域根據土地利用、土壤、坡度的不同進一步劃分為不同的水文響應單元(HRU)。模型可用于模擬地表水和地下水水質水量、預測土地管理措施對不同土壤類型、土地利用方式和管理條件的大面積復雜流域的水文、泥沙和農業污染物的影響[19]。

SWAT模型的水文模擬主要是基于水量平衡，其水文循環基于的水量平衡方程為[20]：

(1)

式中:SWt為土壤最終含水量；SW0為土壤初始含水量；t為時間；Rd為第i天降水量；Qsurf為第i天地表徑流；Ea為第i天蒸發量；Wseep為第i天土壤剖面底層的滲透量和測流量；Qgw為第i天的地下水出流量。

SWAT模型采用修正的MULSE方程來模擬土壤侵蝕過程。修正的MULSE方程用徑流因子代替降水動能，可清楚地分離和輸送泥沙的能量，并可用于單次暴雨事件。具體公式[20]為：

Y=11.8(Qpr)0.56KUSLECUSLEPUSLELSUSLECFRG

(2)

式中：Y為土壤侵蝕量，t；Q為地表徑流，mm；pr為洪峰徑流，m3/s；KUSLE為土壤侵蝕因子；CUSLE為植被覆蓋和作物管理因子；PUSLE為保持措施因子；LSUSLE為地形因子；CFRG為粗碎快土壤成分計算因子。

土壤中的氮磷等營養物隨地表徑流和泥沙輸移進入河道，而附著在泥沙上的營養物隨泥沙進入地表徑流的量為[21]：

(3)

式中：N0為隨泥沙進入地表徑流的營養物量；C0為地表10 mm土層中有機物濃度；Y為模擬步長內的泥沙產量；An為HRU面積；εN營養物富集系數。

而隨地表徑流進入主河道的營養物量表示為[21]：

(4)

2.2.2 適用性評價

SWAT涉及的目標函數類型一共有9種，本文采用目前國內外常用的兩個評價指標—決定系數(R2)和納什效率系數[22](ENS)作為參數優化的評價標準，其表達式如下所示：

(5)

(6)

其中R2表征模擬值和實測值變化趨勢的一致性，R2越接近于1，表示兩者的擬合程度越高；ENS表示實測值和模擬值的偏離程度，ENS越接近于1表示偏離程度越小。當ENS≤0.36時，認為模擬效果不好；當0.36≤ENS<0.75時表示模擬效果令人滿意；當ENS≥0.75時，認為模擬效果好[23]。

3 結果與分析

3.1 模型構建

阿什河流域上游地形起伏較大，主要為山區，下游為適合農業耕種及居住的平原地帶。由于平原區地形起伏不大，直接采用D8算法容易使得柵格水流流向一致而生成不合理的平行偽河道[24]。為了減少這種誤差，采用Burn in算法導入流域實際數字水系圖加以修正。基于阿什河流域90 m分辨率DEM高程數據，通過設置集水區閾值2 500 hm2，將流域劃分為89個子流域。通過對研究區土地利用類型進行重分類和編碼，將研究區劃分為包括森林(FRSD)、灌叢(RNGB)、草地(PAST)、水田(RICE)等10種土地利用類型，詳細的研究區土地利用見圖5(a)。基于中國土壤數據庫，查詢和計算研究區內各類土壤的相關屬性參數，研究區具體的土壤類型可見圖5(b)。根據劃分得到各子流域內土地利用、土壤類型及坡度范圍，最終將研究區劃分為1 255個水文響應單元(HRUs)。

3.2 模型適用性分析

基于SWAT模型，根據阿什河流域地理信息數據和屬性數據進行研究區模型構建，使用SUFI-2算法對模型水量和水質相關參數進行率定。阿什河含沙量較少，該流域并沒有進行泥沙方面的監測，考慮到一部分磷的流失主要是與泥沙的輸移有關，因此在實際總磷參數的率定過程中，將泥沙參數和總磷參數統一率定。率定的總體原則是先進行水量參數的率定，然后再分別對總氮和總磷進行率定。本研究選取2008-2009年作為模型預熱期，以降低初始條件的影響，尤其是初始土壤水的影響，但其模擬結果不參與模型評價，將2010-2013年作為率定期，2014、2015作為驗證期。阿城站月徑流、總氮、總磷在率定期和驗證期的模擬效果如圖2所示。模型適用性評價指標結果如表1所示。

從圖2可知，在率定期和驗證期，月徑流模擬效果總體較總氮、總磷要好，月徑流納什效率系數(ENS)在率定期和驗證期分別為0.70和0.67，確定系數(R2)則分別為0.62和0.60，擬合效果較好。總磷的模擬值和實測值變化趨勢較為吻合，其中在2012年份出現一定的偏差，但月總磷納什效率(ENS)在率定期和驗證期分別為0.49和0.46，確定系數(R2)分別為0.47和0.45，擬合效果符合模擬的精度要求，因此可以根據校準好的模型進行進一步的分析。

圖2 阿什河流域徑流、總氮、總磷模擬值與實測值擬合圖Fig.2 Fitting results of simulated and measured values of runoff, TN and TP in the Ashi River basin

表1 阿什河流域月徑流模擬效果評價指標Tab.1 Evaluation index of monthly runoff, TN and TP simulation in the Ashi River basin

3.3 總氮總磷流失負荷空間分布

基于模型輸出結果，在水文響應單元尺度上統計每一個水文響應單元在2010-2015年總氮總磷輸出負荷大小，最終得到的阿什河流域多年總氮總磷年平均負荷空間分布如圖3和圖4所示。根據阿什河流域水體重污染現狀，從空間分布上針對性的分析總氮總磷的輸出規律。

圖3 總氮負荷空間分布圖Fig.3 The spatial distribution of TN load

由圖3和圖4可知，阿什河流域總氮總磷負荷空間分布都主要集中在下游的平原區以及上游的河岸兩側，總氮負荷最大可達到180.91 kg/hm2，總磷負荷最大的為4.59 kg/hm2。氮磷輸出負荷較小的區域則主要分布在上游的林區，非點源污染空間變化差異較大。氮磷的空間分布特征與流域地形及植被覆蓋類型有很大關系，在上游地形起伏較大、地表覆蓋主要為林地的山區，總氮總磷負荷均比較小，只有在河流附近的農田區域氮磷負荷相對較高；在流域下游、適宜居住和耕種的平原區，總氮總磷負荷相對上游普遍偏高。而在下游，氮磷輸出負荷高的區域集中在地表覆蓋類型為水田的區域，尤其對于總磷；不同的是，土壤類型為草甸黑土的區域，總氮的污染負荷比總磷高。

阿什河流域作為我國重要的糧食生產基地，種植面積大、耕種強度高，使得該流域很容易產生農業非點源污染[25]。化肥施用強度高，而農作物實際吸收較少，多余的營養物則在降水的溶解和沖刷下進入河道，造成河流水質下降和富營養化。考慮到研究區不同區域所受的人類活動影響不同以及研究區本身所固有的自然特征變化，因此分析下墊面特征與氮磷輸出負荷變化規律顯得十分必要。

3.4 氮磷流失負荷與下墊面特征的關系

選取土地利用類型、土壤類型和坡度3個反映下墊面特征的要素來分別分析其與氮磷輸出負荷間的關系。從圖5可以看出，阿什河流域下墊面特征變化較為復雜。從土地利用類型上來看[圖5(a)]，阿什河流域土地利用類型主要以耕地和林地為主。耕地主要分布在中下游，耕地面積占整個流域面積的46.42%，其中水田和旱地分別占全流域的6.58%和38.05%。林地次之，主要分布在流域上游，其面積約占全流域的44.59%，落葉闊葉林和灌叢則分別占流域面積的38.15%和1.02%。而草地、水體、聚落、荒漠所占面積較小，分別占流域面積的1.99%、0.98%、3.28%、0.17%。從空間分布上來看，水田主要以帶狀分布在中下游河岸兩側，在流域中游和河口附近分布有少量面積的建設用地。從土壤類型上來看[圖5(b)]，流域的土壤類型主要有3種，分別是暗棕壤、黑土、草甸土，各占流域面積的37.51%、29.79%、21.87%，其余類型土壤則分布較少，從空間分布上看，暗棕壤主要分布在流域上游，黑土分布在流域下游，而草甸土則沿河道兩側分布，從上游延伸至下游。從土壤坡度上來看[圖5(c)]，研究區主要坡度在9.42°以下，9.42°～ 17.98°(緩坡)以及大于17.98°(陡坡)主要分布在流域上游。流域整體上呈現上游陡峭，中下游平緩。

圖5 研究區土地利用、土壤類型、坡度分布Fig.5 The distribution of land use, soil type and slope in Ashi River Basin

3.4.1 氮磷流失負荷與土地利用類型的關系

對研究區主要土地利用上的氮磷單位輸出負荷及輸出占比進行統計，統計結果如圖6所示。從不同土地利用氮磷輸出負荷方面來看，農業用地的總氮輸出負荷最大，為50.383 kg/hm2，其次為林地、草地、灌叢和水體，分別為11.69、4.87、2.47和0.005 kg/hm2；對于磷來說，輸出負荷最高的土地利用是農業用地和草地，分別為0.951和0.920 kg/hm2，其次為灌叢、濕地和森林，分別占0.483、0.370、0.265 kg/hm2。由此可見，耕地是造成非點源污染的主要土地類型，農業生產所施用的化肥、農藥、殺蟲劑等污染物經遷移進入河道，造成河流水質的下降。不同的是，草地總磷的輸出負荷接近于農業用地，可能是畜牧業的發展有關，動物牲畜的糞便以及牧草施肥都會對草地的磷負荷產生影響。從不同土地利用類型氮磷輸出總量占比來看，不同土地利用總氮總磷輸出占比次序相似，輸出占比最高的是農業用地、其次為林地、草地和濕地。林地雖然輸出負荷比較低，但其所占流域面積很大，枯枝敗葉以及動物尸體等經微生物分解的有機物都會成為氮磷營養物的來源。

圖6 不同土地利用總氮總磷輸出負荷及總量占比Fig.6 The yield of TN and TP and proportion of different land use

3.4.2 氮磷流失負荷與土壤類型的關系

對研究區的不同土壤類型上的氮磷輸出負荷及占比進行統計，統計結果如圖7所示。從不同土壤類型氮磷輸出負荷來看，總氮輸出負荷最大的為城區，達131.25 kg/hm2，其后為水稻土、黑土和草甸土，分別為65.22、39.81、39.18 kg/hm2，比較小的土類為沼澤土、暗棕壤、白漿土，依次是20.81、18.06和14.45 kg/hm2。對于總磷，不同土壤類型的輸出負荷排序與總氮相似，最大的為城區，其后分別是水稻土、黑土、草甸土、白漿土、沼澤土以及暗棕壤，各輸出負荷分別為1.34、0.87、0.83、0.76、0.50、0.46、0.37 kg/hm2。城區污染負荷大的原因可能與降水對城市地表的沖刷有關。城區地面含有許多污染物質，諸如城市垃圾、動物糞便、施工場地堆積物、空氣沉降物等，同時城區地表相對于自然地表，不透水率增加而糙率減少，極易在降水的作用下引發污染物的沖刷和遷移而進入地表水體，最終對水體水質產生很大的影響[26,27]。不同土壤的氮磷輸出負荷不同，除人工施加化肥等因素影響外，土壤的物理化學性質也不可忽略。研究區水稻土質地為黏土，表層土壤成碎塊狀結構，緊實，透水不良，有機質分解慢，養分不易釋放，再加上其本身養分就比較充足，使得土壤固持能力有限，所以損失較為嚴重。研究區黑土土壤表層為黏土，雖呈團塊狀結構，但土質較松，相比水稻土，淋溶作用更明顯。而草甸土與二者相比，腐蝕質層較厚，表層雖為黏土，但土質疏松，多根系，淋溶量以及土壤的吸附和吸收作用更強。不同土地利用磷負荷差異不同于氮負荷的原因可能與磷本身的物理性質有關。磷肥進入土壤后，較難溶解，且溶解的磷容易被土壤顆粒吸附而形成比較穩定的物質，不易被釋放和移動[28]。從不同土壤類型氮磷輸出占比來看，最高的依次均為黑土、草甸土和暗棕壤，三者累計占氮磷總輸出分別為89.8%、90.3%。

圖7 不同土壤類型總氮總磷輸出負荷及總量占比Fig.7 The yield of TN and TP and proportion of different soil type

3.4.3 氮磷失負荷與不同坡度的關系

采用自然斷點法，將研究區坡度分為4個等級，并對研究區不同坡度等級的氮磷輸出負荷進行統計，統計結果如圖8所示。從不同坡度等級氮磷輸出負荷大小來看，氮磷輸出負荷隨坡度的增加大致呈遞減趨勢，對于總氮，4個坡度等級總氮輸出負荷分別為42.86、22.86、17.95、15.40 kg/hm2，對于總磷，分別為0.91、0.47、0.27、0.28 kg/hm2。原因是研究區耕地主要分布于坡度較低的平原區，耕地因農業耕作化肥施用量高，故單位負荷輸出高。當坡度增加到一定程度，不宜種植業的發展，土壤中氮磷將由地表覆蓋和自身土壤條件決定。坡度的增加使得地表所受的雨水條件將發生變化，坡度越高，流速越快，侵蝕作用也相應增加。因為磷易被土壤固定，而氮素容易被淋洗而被徑流攜帶損失，所以氮素主要以徑流流失為主，而磷容易吸附到泥沙上，隨泥沙遷移而流失。張夢等[29]研究表明坡度較小時，磷素流失途徑以徑流流失為主，隨著坡度的增加，磷素的流失途徑以泥沙流失為主。這可能是坡度在18.74～33.07和大于33.07磷元素負荷量表現不同于氮負荷的主要原因。從不同坡度等級氮磷輸出占比來看，占比最高的主要分布在坡度范圍是7.16～18.76和0～7.16區域內，0～18.76坡度范圍氮磷輸出占總量的百分比依次是72.83%、78.15%，要占到全流域總流失量的絕大部分。

圖8 不同坡度范圍總氮總磷輸出負荷及總量占比Fig.8 The yield of TN and TP and proportion of different slope

4 結 語

(1)本研究使用SWAT模型分析了阿什河流域氮磷輸出負荷與下墊面特征的關系，結果表明，率定期徑流的決定系數和納什效率系數分別為0.62、0.70,驗證期分別為0.60、0.67；總氮的率定期決定系數和納什效率系數分別為0.50、0.67，驗證期為0.49、0.75；總磷決定系數和納什效率系數分別為0.46、0.49，驗證期分別為0.45、0.47。除總磷部分年份模擬有一定偏差外，總體模擬結果較好，模型適用于流域非點源污染模擬。

(2)在空間尺度上，流域內總氮總磷輸出負荷空間分布較為相似。總氮總磷單位輸出負荷較高的區域分布在流域中下游主河道附近。在進行流域非點源污染的治理時，可首先對這部分敏感區域進行集中治理，可在一定程度上降低治理難度和提高治理成效。

(3)流域內不同土地利用氮輸出負荷最高的是農業用地、其次為森林、草地等，而對于磷輸出負荷，最大也為農業用地，之后則依次為草地、灌叢、濕地等；不同土壤類型單位輸出負荷最高的是城區，然后為水稻土、黑土、草甸土等，不同土壤類型總磷輸出負荷排序與總氮類似；研究區不同坡度范圍的氮磷輸出負荷隨坡度的增加而呈下降趨勢。

□