王 冶，薛忠財(cái)*，武智勇，靳 松，姜白楊

(1.河北民族師范學(xué)院資源與環(huán)境科學(xué)學(xué)院，河北承德 067000；2.河北省山區(qū)地質(zhì)環(huán)境重點(diǎn)實(shí)驗(yàn)室，河北承德 067000；3.承德市古代建筑保護(hù)研究所，河北承德 067000)

產(chǎn)水服務(wù)與氮磷凈化能力對流域的水質(zhì)及水量的變化有直接的影響，是流域生態(tài)系統(tǒng)研究中非常重要的兩項(xiàng)服務(wù)[1]。其中，產(chǎn)水服務(wù)常被作為區(qū)域生態(tài)系統(tǒng)狀況的重要指示器，而氮磷凈化能力則是水質(zhì)凈化功能的重要體現(xiàn)[2]。產(chǎn)水服務(wù)研究大多是基于大量實(shí)測數(shù)據(jù)進(jìn)行估算，主要包括水量平衡法、降水存儲量法等[3]，但上述方法難以獲取較高精度數(shù)據(jù)，從而進(jìn)行大區(qū)域估算時(shí)會產(chǎn)生較大誤差。水質(zhì)評價(jià)方法有單因子法、水質(zhì)指數(shù)法、多變量分析法等[4]，這些方法多用來評價(jià)水質(zhì)狀況，無法針對水質(zhì)凈化功能進(jìn)行評估[5]。近年來，隨著計(jì)算機(jī)科學(xué)的不斷進(jìn)步，諸如水文評估模型Swat[6]、半分布式水文水質(zhì)模型HSPF[7]、水文生態(tài)模型InVEST[8]等已被應(yīng)用于產(chǎn)水服務(wù)與氮磷凈化能力分布及動態(tài)特征研究。其中，由斯坦福大學(xué)開發(fā)的InVEST模型在國內(nèi)外得到廣泛應(yīng)用[9]。Leh等[10]成功利用InVEST模型對西非國家加納和科特迪瓦的產(chǎn)水量進(jìn)行了評估，Goldstein等[11]利用InVEST模型分析了不同情境下的夏威夷島嶼水質(zhì)凈化功能空間分異情況；我國學(xué)者利用InVEST模型的產(chǎn)水模塊分析了黃河流域[12]、太湖流域[13]、黑河流域[14]、三江源地區(qū)[15]、石羊河流域[16]等地區(qū)的產(chǎn)水量及水源涵養(yǎng)時(shí)空分布特征，并且利用InVEST模型的NDR模塊(水質(zhì)凈化模塊)對珠江流域[5]、海河流域[17]、官廳水庫流域[18]、海南島[19]等流域的氮磷凈化功能開展研究，并取得較好的評估效果。

武烈河流域是承德地區(qū)主要供水水源地，也是京津冀水源涵養(yǎng)功能區(qū)重要組成部分，其水量與水質(zhì)有重要的生態(tài)價(jià)值與研究意義，但是截至目前仍缺乏針對該區(qū)域產(chǎn)水服務(wù)和水質(zhì)凈化服務(wù)的研究。當(dāng)前針對武烈河流域的研究主要集中于徑流模擬[20]、水體污染特征[21]、水質(zhì)污染源解析[22]等，但從生態(tài)系統(tǒng)服務(wù)角度對水量與水質(zhì)的研究鮮見報(bào)道。筆者以武烈河流域?yàn)檠芯繀^(qū)域，基于InVEST模型產(chǎn)水模塊與NDR模塊，分析了1990—2018年產(chǎn)水服務(wù)及氮磷凈化能力的時(shí)空分布特征，并利用情景模擬法量化了氣候要素與土地利用/覆被變化對其影響的程度，以期為京津冀水源涵養(yǎng)功能區(qū)的合理規(guī)劃提供科學(xué)支撐。

1 資料與方法

1.1 研究區(qū)概況武烈河流域位于河北省承德市境內(nèi)，屬灤河一級支流，古有熱河之稱，是承德市最重要的供水水源地之一，干流全長110 km，流域面積2 589 km2，主要由武烈河干流、玉帶河、鸚鵡河、興隆河這4條水系構(gòu)成(圖1a)，界于117°45′～118°26′E、40°54′～41°42′N；屬于暖溫帶和寒溫帶過渡地帶，四季分明，氣候?yàn)榈湫痛箨懷嗌缴降貧夂颍昃鶜鉁?.9 ℃，年降雨量為537.2 mm，年均徑流量2.60×108m3，多年平均流量6.93 m3/s，降水的年內(nèi)分配十分不均勻，降水期與汛期主要集中在6—9月，占全年降雨量的77.9%。流域內(nèi)主要土地覆被類型為林地、草地，占總面積的 75.1%，其次為耕地與建筑用地，分別占總面積的19.3%、4.7%(圖1b)，主要土壤類型為棕壤、褐土和草甸土。

圖1 研究區(qū)地理位置(a)和土地利用類型(b)Fig.1 Geographical location(a)and land use type(b)of the study area

1.2 研究方法

1.2.1產(chǎn)水服務(wù)計(jì)算原理。InVEST模型產(chǎn)水量計(jì)算是根據(jù)水循環(huán)的原理，通過降水、植物蒸騰、地表蒸散發(fā)、根系深度、土壤深度和植被可利用水等參數(shù)經(jīng)過模型計(jì)算獲得。計(jì)算公式如下：

(1)

(2)

(3)

(4)

AWCxj=min(MSDxj,RDxj)×PAWCxj

(5)

式中，Yxj表示土地利用/覆被類型j中柵格單元x上的年產(chǎn)水量；AETxj為土地利用/覆被類型j上柵格單元x的年平均蒸散量；Pxj為柵格單元x土地利用/覆被類型j的年降雨量；Rxj表示土地利用/覆被類型j上柵格單元x的干燥指數(shù)；Kcxj為柵格單元x上土地利用/覆被類型j的作物系數(shù)；LAI表示植被葉面積指數(shù)；ET0表示潛在蒸散發(fā)量；Z為Zhang系數(shù)，是用于表征研究區(qū)降水特征的常數(shù)，取值為1～10[23]；ωxj是氣候土壤非物理參數(shù)；AWCxj表示柵格單元x上土地利用/覆被類型j的有效可利用水；MSDxj為土層厚度最大值(m)；RDxj為根系深度(m)；PAWCxj表示柵格單元x上土地利用/覆被類型j的植被可利用水，計(jì)算方法采用周文佐等[24]提出的模型進(jìn)行估算，公式如下：

PAWCxj=54.509-0.132SSAN-0.003(SSAN)2-0.055SSIL-0.006(SSIL)2-0.738CCLA+0.007(CCLA)2-0.688OM+0.501(OM)2

(6)

式中,SSAN為土壤砂粒含量(%)；SSIL表示土壤粉粒含量(%)；CCLA為土壤黏粒含量(%)；OM為土壤有機(jī)質(zhì)含量(%)。

1.2.2氮磷凈化能力計(jì)算原理。氮磷凈化能力采用InVEST模型中的NDR模塊(水質(zhì)凈化模塊)進(jìn)行計(jì)算，主要運(yùn)用植被和土壤可以通過儲存和轉(zhuǎn)換等方式移除或減少徑流中的營養(yǎng)鹽污染物以達(dá)到凈化水質(zhì)的作用的機(jī)理，該模型只考慮非點(diǎn)源污染中的總氮(TN)、總磷(TP)，二者輸出量越高，表明氮磷凈化能力越低，用于評估生態(tài)系統(tǒng)中植被和土壤的水質(zhì)凈化服務(wù)，其主要算法如下：

ALVx=HSSx×polx

(7)

式中，ALVx是柵格x的調(diào)整后輸出量，HSSx為柵格x的水文敏感性得分，polx為柵格x的輸出系數(shù)。計(jì)算得到TN、TP輸出量之后，再結(jié)合各土地利用/覆被類型對污染物的移除效率來計(jì)算得出營養(yǎng)物保持(持留)量。

1.2.3情景模擬法。基于情景模擬的方法，研究氣候要素與土地利用/覆被變化對武烈河流域產(chǎn)水服務(wù)功能、氮磷凈化能力的貢獻(xiàn)程度。為排除異常氣候數(shù)據(jù)的影響，采用1990—1999、2000—2009、2010—2018年3個(gè)時(shí)期的平均模擬產(chǎn)水量、TN輸出量、TP輸出量進(jìn)行研究。設(shè)計(jì)情景如表1所示。

表1 氣候變化和土地利用/覆被變化情景的各時(shí)段數(shù)據(jù)輸入

假設(shè)各時(shí)期的氣候要素按前一時(shí)期輸入，土地利用/覆被按實(shí)際情況輸入，計(jì)算得到氣候變化情境下的模擬值，該值與真實(shí)情景模擬值的差即為氣候變化的影響量；同理，假設(shè)各時(shí)期的土地利用/覆被按前一時(shí)期輸入，氣候要素按實(shí)際情況輸入，土地利用/覆被變化情景下的模擬值與真實(shí)情景模擬值的差即為土地利用/覆被變化的影響量?；诖耍瑲夂蛞嘏c土地利用/覆被影響的貢獻(xiàn)比值可以用以下公式進(jìn)行量化：

(8)

(9)

(10)

式中，RWY、RTN、RTP分別表示氣候要素與土地利用/覆被對產(chǎn)水量、TN輸出量、TP輸出量影響的貢獻(xiàn)比值，該值越接近1表明氣候要素與土地利用/覆被的影響程度相近，越大表明氣候要素影響程度大于土地利用/覆被，越小表明土地利用/覆被影響程度大于氣候要素；WYO、TNO、TPO分別表示真實(shí)情景下的各時(shí)期年均產(chǎn)水量(m3)、TN輸出量(t)、TP輸出量(t)；WYC、TNC、TPC分別表示氣候變化情境下的各時(shí)期年均模擬產(chǎn)水量(m3)、TN輸出量(t)、TP輸出量(T)；WYL、TNL、TPL分別表示土地利用/覆被變化情景下的各時(shí)期年均模擬產(chǎn)水量(m3)、TN輸出量(t)、TP輸出量(t)。

1.3 數(shù)據(jù)來源與處理

1.3.1數(shù)據(jù)來源。該研究收集了承德市武烈河流域氣象、土地利用/覆被、土壤屬性、流域矢量邊界、數(shù)字高程模型(DEM)、生物物理屬性以及氮磷輸出負(fù)荷數(shù)據(jù)等用于模型計(jì)算，同時(shí)還整理了武烈河水文站徑流監(jiān)測數(shù)據(jù)與河流斷面總氮、總磷監(jiān)測數(shù)據(jù)對模型結(jié)果進(jìn)行驗(yàn)證。其中，DEM 數(shù)據(jù)來源于地理空間數(shù)據(jù)云平臺(http：//www.gscloud.cn)，7期土地利用/覆被數(shù)據(jù)(分別為1990、1995、2000、2005、2010、2015、 2018年)來源于資源環(huán)境數(shù)據(jù)云平臺(http：//www.resdc.cn)；土壤屬性數(shù)據(jù)來源于世界土壤屬性數(shù)據(jù)庫(HWSD)；氣象數(shù)據(jù)(包括降水、氣溫、氣壓、風(fēng)速等)來源于國家氣象信息中心(http：//data.cma.cn)與國家青藏高原科學(xué)數(shù)據(jù)中心(http：//data.tpdc.ac.cn)；生物物理屬性以及氮磷輸出負(fù)荷數(shù)據(jù)參考InVEST 用戶指南和查閱研究區(qū)附近區(qū)域的相關(guān)文獻(xiàn)[25-26]后估算獲得，相關(guān)參數(shù)見表2。

表2 各土地利用/覆被類型生物物理屬性及氮磷輸出主要參數(shù)

1.3.2數(shù)據(jù)處理。降水?dāng)?shù)據(jù)通過使用ArcGIS對國家青藏高原科學(xué)數(shù)據(jù)中心下載的全國月均降水量柵格數(shù)據(jù)進(jìn)行疊加、裁剪獲取。土地利用與土壤質(zhì)地經(jīng)ArcGIS對世界土壤屬性數(shù)據(jù)庫進(jìn)行投影變換、裁剪獲取。集水區(qū)利用ArcGIS水文分析模塊的工具，根據(jù)流域與子流域的邊界將研究區(qū) DEM數(shù)據(jù)分割成小流域，共產(chǎn)生4個(gè)子流域。植被可利用水是基于世界土壤屬性數(shù)據(jù)庫(HWSD)中研究區(qū)內(nèi)土壤黏粒含量、土壤粉粒含量、土壤砂粒含量、土壤有機(jī)質(zhì)含量由經(jīng)驗(yàn)公式(6)計(jì)算得到。根據(jù)史曉亮等[20]對武烈河徑流的模擬和產(chǎn)水量的估算，通過公式(2)、(3)、(4)對Zhang系數(shù)反復(fù)驗(yàn)證校準(zhǔn)，最終確定Zhang系數(shù)為3.7。潛在蒸散發(fā)量采用經(jīng)過系數(shù)校正后的Hargreaves公式進(jìn)行計(jì)算獲取[27]。

2 結(jié)果與分析

2.1 武烈河流域產(chǎn)水服務(wù)功能與氮磷凈化能力的時(shí)空分布特征經(jīng)分析，武烈河流域多年平均產(chǎn)水深度為104.22 mm，平均產(chǎn)水總量為2.69×108m3；多年平均TN、TP輸出總量分別為291.52、28.22 t，單位面積多年平均TN、TP輸出總量分別為1.17、0.11 kg/(hm2·a)。

為精確量化武烈河流域的產(chǎn)水服務(wù)功能與氮磷凈化能力，將數(shù)據(jù)分為1990—1999、2000—2009、2010—2018年3個(gè)時(shí)期進(jìn)行分析。從時(shí)間上看，產(chǎn)水服務(wù)功能方面(圖2a)，1990—1999年產(chǎn)水量呈減小趨勢，2000—2018年產(chǎn)水量呈增加趨勢，總體呈微弱下降趨勢，多年平均降幅為0.04×108m3/a；最大值出現(xiàn)在2016年，為6.67×108m3，最小值出現(xiàn)在2009年，為0.58×108m3。氮磷凈化能力方面(圖2b)，TN、TP輸出量在1990—2009年變化趨勢不明顯，2010—2018年顯著增加，總體呈增加趨勢，增幅分別為0.02、0.06 t/a；TN輸出量在2017年最大，為293.86 t，在2009年最小，為290.04 t；TP輸出量在2017年最大，為29.66 t，在1999年最小，為27.74 t。

圖2 1990—2018年武烈河流域年均產(chǎn)水量(a)與氮、磷輸出量(b)變化Fig.2 Variation in annual average water yield(a)and TN and TP output(b)in Wulie River Basin in 1990-2008

從空間上看，產(chǎn)水服務(wù)功能方面(圖3a)，產(chǎn)水深度呈東高西低、南高北低的空間分布特征，這主要是由降水、蒸散發(fā)、土地利用/覆被等因素的空間分異導(dǎo)致的；平均產(chǎn)水深度1990—1999年為139.65 mm，2000—2009年為58.19 mm，2010—2018年為104.38 mm。氮磷凈化能力方面(圖3b、c)，TN、TP輸出量高的地區(qū)均分布在地勢平坦和人口分布密集區(qū)，武烈河流域以種植農(nóng)作物和果園為主，農(nóng)藥、化肥的使用是氮、磷的主要來源，有部分農(nóng)業(yè)區(qū)距離河道較近，因此瀕臨河道處應(yīng)擴(kuò)大林地面積，進(jìn)一步增強(qiáng)該區(qū)水質(zhì)凈化能力；年均TN、TP輸出量1990—1999年分別為1.16、0.10 kg/(hm2·a)，2000—2009分別為1.15、0.09 kg/(hm2·a)，2010—2018年分別為1.18、0.12 kg/(hm2·a)。

圖3 1990—2018年武烈河流域各時(shí)段年均產(chǎn)水量(a)與TN輸出量(b)、TP輸出量(c)時(shí)空分布變化Fig.3 Temporal and spatial distribution changes of the average annual water production(a)，nitrogen output(b)，and phosphorus output(c)in the Wulie River Basin from 1990 to 2018

2.2 氣候因素與武烈河流域產(chǎn)水功能、氮磷凈化能力的關(guān)系基于線性回歸方法，分析1990—2018年武烈河流域各年均產(chǎn)水量以及TN、TP輸出量與對應(yīng)年份的降水量、潛在蒸散發(fā)量、溫度之間的關(guān)系。結(jié)果表明(圖4)，武烈河流域降水量與產(chǎn)水量呈顯著正相關(guān)(P<0.05)，與TN、TP輸出量均呈正相關(guān)，但不顯著(P>0.05)；潛在蒸散發(fā)量與產(chǎn)水量呈顯著負(fù)相關(guān)(P<0.05)，與TN、TP輸出量均呈負(fù)相關(guān)，但不顯著(P>0.05)；氣溫與產(chǎn)水量呈負(fù)相關(guān)，與TN、TP輸出量均呈正相關(guān)，但相關(guān)性均不顯著(P>0.05)。由此可見，武烈河流域的產(chǎn)水量與降水量(R2=0.987 8)、潛在蒸散發(fā)量(R2=0.539 5)具有較強(qiáng)的相關(guān)性，與氣溫存在一定相關(guān)性但不顯著；流域的TN、TP輸出量與各氣候因素不顯著。因此，各類氣候因素中，降水量、潛在蒸散發(fā)量是影響武烈河流域產(chǎn)水量的主要因素，而對流域的TN、TP輸出量影響相對較小。

圖4 武烈河流域各氣象因素與產(chǎn)水深度、TN輸出量、TP輸出量的關(guān)系Fig.4 Correlation analysis of each climatic factor and water yield depth，TN output and TP output in Wulie River Basin

2.3 土地利用/覆被變化對武烈河流域產(chǎn)水功能與氮磷凈化能力的影響1990—2018年武烈河流域各類土地利用/覆被對產(chǎn)水量、TN輸出量、TP輸出量的多年平均貢獻(xiàn)率與單位面積平均值表明(表3)，林地與草地的產(chǎn)水量貢獻(xiàn)率最高，占比達(dá)73.33%，其次為耕地、建筑用地、水域；建筑用地的單位面積產(chǎn)水能力最強(qiáng)，為1 471.99 m3/(hm2·a)，耕地次之，其次為林地、草地、水域。耕地的TN、TP輸出量貢獻(xiàn)率最高，占比分別為83.15%、63.38%，其次為草地、林地、建筑用地、水域；耕地的單位面積TN、TP輸出量也是最高的，分別為4.65、0.34 kg/(hm2·a)，其次為建筑用地、草地、水域、林地。由此可知，林地、草地雖然產(chǎn)水能力低于耕地和建筑用地，但因其面積占比高，故武烈河流域產(chǎn)水量主要來源于林地、草地；TN、TP輸出能力最強(qiáng)的為耕地、建筑用地，原因是耕地化肥農(nóng)藥的使用增加了氮磷的輸出量，同時(shí)建筑用地城鎮(zhèn)化建設(shè)明顯，人口和人類活動的增加也導(dǎo)致污染物排放遞增，從而導(dǎo)致氮磷凈化能力較弱，而林地、草地對污染物的截留效用較強(qiáng)，導(dǎo)致氮磷凈化能力強(qiáng)。土地利用/覆被的變化會導(dǎo)致局部水汽循環(huán)、地表徑流及污染物截留效用的變化，進(jìn)而導(dǎo)致產(chǎn)水服務(wù)功能、氮磷凈化能力的變化。隨著武烈河流域城市化的加劇，大面積的林地、草地和耕地轉(zhuǎn)化為建筑用地(表4)，導(dǎo)致其攔截降水與污染物截留效用的減弱，同時(shí)，這些區(qū)域更加頻繁的人類活動也導(dǎo)致污染物排放的增加，從而導(dǎo)致產(chǎn)水服務(wù)功能增強(qiáng)及氮磷凈化能力的減弱。

表3 武烈河流域各類土地利用/覆被對產(chǎn)水量、TN輸出量與TP輸出量的貢獻(xiàn)率及單位面積平均值

2.4 氣候和土地利用/覆被變化的影響程度量化分析氣候要素與土地利用/覆被變化是導(dǎo)致武烈河流域產(chǎn)水服務(wù)功能、氮磷凈化能力變化的主要原因，兩者的影響程度存在差異，基于情景模擬法計(jì)算結(jié)果表明(表5)，1990—1999年，氣候要素在產(chǎn)水服務(wù)功能、氮磷凈化能力方面分別是土地利用/覆被影響程度的58.17、0.78、0.64倍；2000—2009年，氣候要素在產(chǎn)水服務(wù)功能、氮磷凈化能力方面分別是土地利用/覆被影響程度的26.17、0.49、0.11倍；2010—2018年，氣候要素在產(chǎn)水服務(wù)功能、氮磷凈化能力方面分別是土地利用/覆被影響程度的38.68、0.62、0.32倍。由此可見，1990—2018年在產(chǎn)水服務(wù)功能、氮磷凈化能力方面，氣候要素與土地利用/覆被影響程度比的平均值分別為41.01、0.63、0.36倍。氣候要素變化對產(chǎn)水服務(wù)功能的影響程度高于土地利用/覆被變化，而土地利用/覆被變化對氮磷凈化能力的影響程度高于氣候要素變化。該結(jié)論與前人研究結(jié)果相似，如王亞慧等[28]研究表明橫斷山區(qū)的產(chǎn)水量與降水量存在極顯著正相關(guān)；趙亞茹等[23]研究表明氣候要素是石羊河上游產(chǎn)水量的主要影響因素；竇攀烽等[29]通過量化氣候與土地利用/覆被對寧波地區(qū)產(chǎn)水量的影響，確定氣候變化對產(chǎn)水量的貢獻(xiàn)率在97%以上；張建等[2]研究表明丹江口市的氮磷輸出量與土地利用/覆被的景觀形狀指數(shù)呈顯著正相關(guān)；吳瑞等[18]研究表明官廳水庫區(qū)域的氮磷凈化能力的強(qiáng)弱與土地利用/覆被類型呈顯著相關(guān)。

表4 武烈河流域各土地利用/覆被類型面積占比轉(zhuǎn)移矩陣Table 4 Area proportion transfer matrix of various land use/land cover types in Wulie River Basin %

表5 氣候要素與土地利用/覆被變化對武烈河流域產(chǎn)水服務(wù)功能、氮磷凈化能力的影響程度

3 結(jié)論

該研究以武烈河流域?yàn)檠芯繉ο?，運(yùn)用InVEST模型對武烈河流域1990—2018年產(chǎn)水服務(wù)功能、氮磷凈化能力的時(shí)空變化進(jìn)行計(jì)算，分析了與不同氣候因素的相關(guān)性，并量化了各類土地利用/覆被對產(chǎn)水量、TN輸出量與TP輸出量的貢獻(xiàn)率，采用情景模擬法對影響武烈河流域產(chǎn)水功能與氮磷凈化能力的因素進(jìn)行量化。具體結(jié)果如下：

(1)1990—2018年，武烈河流域產(chǎn)水量呈現(xiàn)先減小后增加，總體呈微弱減小的變化趨勢，多年平均產(chǎn)水量為2.69×108m3；TN、TP輸出量呈增加趨勢，多年平均TN、TP輸出量分別為291.52、28.22 t。近29年來，武烈河流域產(chǎn)水服務(wù)功能與氮磷凈化能力均減弱。

(2)武烈河流域產(chǎn)水量與降水量呈顯著正相關(guān)，與潛在蒸散發(fā)量呈顯著負(fù)相關(guān)，與氣溫呈負(fù)相關(guān)但不顯著；TN、TP輸出量與各氣候要素呈一定相關(guān)關(guān)系，但均不顯著。氣候要素與產(chǎn)水量的相關(guān)性較強(qiáng)，與TN、TP輸出量的相關(guān)性較弱。

(3)武烈河流域產(chǎn)水量主要來源于林地和草地，占比達(dá)73.33%，其次為耕地、建筑用地、水域；TN、TP輸出量主要來源于耕地，占比分別達(dá)83.15%、63.38%，其次為草地、林地、建筑用地、水域。建筑用地產(chǎn)水能力最強(qiáng)，單位面積產(chǎn)水量為1 471.99 m3/(hm2·a)；耕地TN、TP輸出能力最強(qiáng)，單位面積輸出量分別為4.65、0.34 kg/(hm2·a)。

(4)氣候要素變化對產(chǎn)水服務(wù)功能的影響更為顯著，而土地利用/覆被變化對氮磷凈化能力的影響更為顯著。1990—2018年產(chǎn)水服務(wù)功能方面，氣候要素是土地利用/覆被影響程度的41.01倍；氮磷凈化能力方面，氣候要素分別是土地利用/覆被影響程度的0.63、0.36倍。