基于InVEST模型的大理州永平縣水生態系統服務功能時空變化分析*

杜佳衡，王錦

(西南林業大學，云南 昆明 650224)

水生態系統作為自然生態系統中的一個重要類型，為人類的生存繁衍提供了多種物質產品及資源，具有多種服務功能，如集水、產水、水質凈化，同時兼具一定的美學和娛樂游憩價值[1-3]。其中，產水和水質凈化作為水生態系統主要的服務功能，具有提高民生福祉及改善城鎮居民飲用水安全的重要作用[4-5]。然而，在當前快速城市化進程中，人類社會生產活動對水生態系統的干擾逐漸增強，導致其服務功能日益衰退[6]。為明確區域水生態系統的健康狀況，開展產水和水質凈化服務功能的空間可視化定量分析，揭示兩者的時空變化分布特征，對于實現區域水資源管理和可持續發展意義深遠。

生態系統服務自20世紀70年代提出以來，成為生態學、生態經濟學等學科領域中的重要研究方向。在SCEP(study of critical environmental problems)[7]、Costanza等[8]、Daily[9]的研究基礎上，生態系統服務的概念、分類和評估方法逐漸被明確，水生態系統服務功能的相關研究也越來越廣泛[10]。從研究對象和尺度上看，以往學者主要以流域[11-12]、濕地[13]、河流[14]等類型來進行水生態系統服務功能評估研究。從研究方法來看，目前針對水生態系統產水服務的研究主要基于水量平衡法、蓄水能力法等，運用實測數據開展產水量的評估[15-16]。水質凈化服務往往通過單因子法、灰色綜合評價法等方法分析評估區域水質狀況[17]。以上方法受水文站點實測數據影響較大，僅適用于小尺度范圍內的研究，且多采用價值量評價法進行估算，結果未能實現空間化。近年來，在生態系統服務物質量評價法的基礎上，隨著ArcGIS、RS等技術的迅速發展，對于生態系統服務“模型化、空間化”的研究成為新趨勢，各類水文研究模型隨之得到廣泛應用，包括InVEST模型、SWAT模型等[18]。其中，InVEST模型以簡便及快速的空間表達手段成為當下最熱門的研究工具，越來越多的學者借助該模型的產水模塊(water yield model)和水質凈化模塊(nutrient delivery ratio model)進行區域水生態系統服務功能的研究。如白楊等[19]、陳龍等[20]利用InVEST模型分別對美國肯塔基州和瀾滄江流域的水源涵養、水土保持、水質凈化功能進行了評估；吳瑞等[21]、韓會慶等[17]分別對官廳水庫流域和珠江流域的產水和水質凈化服務進行了模擬，并對其時空分布差異性進行驅動因素的探究。綜上所述，InVEST模型因易操作、數據易獲取等特點已被國內外學者廣泛應用于不同尺度的研究當中，但目前針對國內西南部山區縣域尺度的相關研究仍然較少。

永平縣地處云南大理州西部山區，近年來，當地發展面臨著諸多壓力和挑戰，特別是開發活動的不斷擴張造成了區域內資源的快速消耗甚至是生態破壞。為緩解當地工農業生產與自然過程之間的矛盾，有必要對其水生態系統產水服務和水質凈化服務功能展開分析與評估，而目前，該區域內對于這兩種服務功能的研究尚未形成一套較為完整的研究思路及方法。因此，本文基于InVEST模型產水量和水質凈化量兩個子模塊模擬2000年、2010年、2020年大理州永平縣的年均產水量及年均氮磷輸出量，分析這兩個水生態系統主要服務功能的時空變化，并劃分重要性等級，以期為當地經濟健康發展、水源水質保護提供科學依據。

1 研究地概況

永平縣(99°17′～99°56′E、25°03′～25°45′N)位于云南省西部，總面積約為2 790 km2，屬大理州下轄縣?？h內與云臺、博南兩山相連，地勢狹長且高差明顯，自西北向東南降低，海拔最高2 974 m，最低1 105 m。銀江河分布于永平縣中部成為主要干流，瀾滄江、順濞河經縣域東西兩側穿越而過，形成“兩山三河”的天然山水格局和山、谷、壩錯綜復雜的地形地貌[22]。水資源豐富，平均水資源總量超過10×108m3，人均水資源量約7 000 m3，有開發利用價值的河流27條(圖1)。

圖1 研究地概況

2 研究方法

2.1 數據來源及預處理

InVEST產水(water yield model)和水質凈化(nutrient delivery ratio model)模型運行參數主要包括以下內容。

(1)土地利用數據 數據來自資源環境數據云，選取2000年、2010年、2020年3期柵格數據，統一分辨率為30 m。根據劉紀遠等[23]提出的體系，按照二級分類標準將研究區土地利用劃分為5大類12小類(圖2)。

圖2 大理州永平縣土地利用類型

(2)氣象數據 包括降水、氣溫等日值數據，均來自中國氣象數據網，包括研究區及其周邊28個氣象站點的數據，為使數據具備代表性，分別選取1996—2000年、2006—2010年、2016—2020年3個時段的數據滑動平均值，表征2000年、2010年、2020年的氣象數據。其中，年均降水量在ArcGIS中通過反距離權重法(IDW)處理插值后獲得柵格數據；潛在蒸散發量在參考FAO推薦的Penman-Monteith公式[24-25]計算后，通過IDW插值為柵格數據。

(3)DEM 數據來自地理空間數據云，經投影、填洼處理后借助水文分析工具提取流域及子流域。

(4)土壤數據 數據來自中科院南京土壤所，土壤深度數據在ArcGIS中經投影變換、裁剪處理獲取柵格數據，植被可利用含水率數據參考周文佐等[26]擬合的經驗公式計算提取。

(5)生物物理系數表 參考前人研究成果[27-28]及模型用戶使用手冊[29]賦值設定(表1)。

表1 InVEST模型運行參數生物物理系數表

(6)經驗參數 產水模塊Z系數通過《大理白族自治州水資源公報》調整設定。由于數據獲取限制，僅收集到研究區2010年水資源統計成果?；诖藢Ξa水模塊輸出結果進行校驗，發現當Z系數賦值為30時，模型的產水量更接近實測結果。其他經驗參數如水質凈化模塊水文連通性K系數等使用模型默認值設定。

2.2 產水服務模擬

InVEST 產水模型(water yield model)是基于Budyko水熱耦合平衡原理[30]提出的一種水量平衡估算方法，即以各柵格單元降水量與實際蒸散量之差作為該柵格單元產水量，模型不考慮地下水的補給。產水量計算公式如下：

式中：Y(xj)為地類j在柵格單元上x的年均產水量(mm)，AET(xj)為柵格單元x上地類j的年均實際蒸散量(mm)，P(x)為柵格單元x上的年均降水量(mm)。其中，AET(xj)/P(x)為Zhang等[31]在參考Budyko曲線后提出的方程，計算公式如下。

PET(xj)=Kc(lc)×ET0(x)

AWC(c)=Min(Rest.layar.depth,root.depth)×PAWC

AET(x)=Min〔Kc(lc)×ET0(x),P(x)〕

式中：PET(xj)為地類j在柵格單元x上的年均潛在蒸散發量(mm)，ω(x)為自然氣候—土壤屬性的非物理參數；ET0(x)為參考植被在柵格單元x上的蒸散量(mm)，Kc(lx)為地類x的植被蒸散系數；AWC(x)為柵格單元x的土壤有效含水率(mm)，PAWC為植被可利用水；Z即Zhang系數，季節常數，取值在1～30之間。

2.3 水質凈化服務模擬

水質凈化服務通過InVEST NDR Model進行計算，模型基于土地利用/覆被類型在年尺度上計算柵格單元氮、磷營養鹽的保持量和輸出量，利用氮、磷含量表征水質狀況，TN、TP輸出量越大表明柵格單元水質凈化能力越弱。模型僅考慮植被和土壤對氮、磷營養鹽非點源污染的截留及過濾功能，不涉及化學和生物之間的相互作用。水質凈化服務計算公式如下：

Xexpj=loadsurf,j×NDRsurf,j+loadsubs,j×NDRsubs,j

式中：xexptot為流域年均氮磷營養鹽輸出總量(kg/a)；xexpj為柵格單元j的氮磷營養鹽輸出量(kg/a)，loadsurf,j和loadsubs,j為柵格單元j地表和地下氮磷營養鹽負荷量(kg/a)，NDRsurf,j和NDRsubs,j為柵格單元j地表和地下氮磷營養鹽輸送效率。其中，氮磷營養鹽負荷量計算公式如下：

loadsurf,j=(1-proportion_subsurfacej)×modified.load_nj

loadsubs,j=proportion_subsurfacej×modified.load_nj

modified.loadxj=loadxj×RPIxj

式中：modified.loadxj為修正后柵格單元j的氮磷營養鹽負荷量(kg/a)，RPIxj為徑流代理指數，RPj為營養鹽徑流代理，RPav為平均徑流代理。

氮磷營養鹽輸送效率計算公式如下：

式中：NDRsurf,j為地表營養鹽輸送效率；IC0和k是校準參數；ICj為地形指數；NDR0,j為柵格單元j和徑流之間的養分最大持留效率；effLULCj為地類j的最大持留效率；eff’downj為柵格單元j的下游有效持留效率；sj為步長因子；ljdown為柵格單元j到下游相鄰柵格單元流動的路徑長度；lLULCj為柵格單元土地利用類型j的持留距離；Dup和Ddn為上下坡貢獻面積的平均坡度；NDRsubs,j為地下營養鹽輸送效率；effsubs為通過地下徑流可達到的最大營養鹽持留效率；lsubs為地下營養鹽流動持留效率；li表示從柵格單元到徑流的距離。

2.4 冷熱點分析

3 結果與分析

3.1 產水服務時空變化模擬

基于InVEST模型產水(Water Yield)模塊運行結果，大理州永平縣2000年、2010年、2020年產水量空間分布格局如圖3所示。

圖3 大理州永平縣產水服務時空分布特征

20年間，研究區產水量總體趨于下降趨勢，不同年份間，產水量的空間分布格局基本保持穩定，大體呈現出西南高東北低的態勢。3個時期流域的平均產水深度分別為642.37、348.95、275.59 mm，高值區集中于龍門鄉東北部、博南鎮中部、北斗鄉西南部；較高值區集中于杉陽鎮西部、龍街鎮西部、廠街鄉中部、水泄鄉中北部；低值區集中連片分布于北斗鄉中部、龍街鎮東部、博南鎮西北部和南部。

同一年份中，單位面積產水深度在空間上的分布具有異質性。按照一級分類標準，對土地利用重新分類后進行分區統計，結果如圖4a所示。從各地類上看，20年來建設用地產水量始終占據高值，且始終高于流域年均產水深度，這可能與快速城市化進程改變的不透水下墊面有關，這一過程改變了水量的產匯流，使得地表產流增加而入滲減少，產水能力相對較高；耕地產水能力次之，這可能與人類引水灌溉等農耕干擾行為存在一定的關聯性；林地和草地產水能力中等，這可能是由于植被及其枯落物對降水和地表徑流的攔截作用所導致的，同時受到水量入滲土壤的影響，削減了洪峰流量，因此產水能力相對較弱；而水域的產水量始終處于低值，這是由于水域的蒸散發量較高，大于植被的蒸散發量，因此產水能力最弱。

從鄉鎮尺度上看，3個時期產水深度如圖4b所示，產水深度最高的鄉鎮分別為龍門鄉、杉陽鎮、水泄鄉；產水深度低值區則集中在北斗鄉、龍街鎮、博南鎮?？紤]到土地利用類型的影響，這種分布空間差異可能與各時期鄉鎮間的用地類型比例有關，可以看出杉陽鎮、水泄鄉、龍門鄉等用地類型以耕地為主導的鄉鎮產水量最高，而北斗鄉、龍街鎮、博南鎮等以林草地為主要用地類型的鄉鎮產水量則相對較少。同時，各時期氣候變化導致的降水量、蒸散發量的大小也會對產水深度產生一定影響，從而影響產水量的時空分布格局特征。

圖4 大理州永平縣各地類、各鄉鎮產水量

3.2 水質凈化服務時空變化模擬

經過InVEST模型(NDR)模塊計算，大理州永平縣2000年、2010年、2020年水質凈化能力如圖5所示，TN、TP輸出量越大則水質凈化能力越差。結果顯示，20年來研究區TN、TP總輸出量先增大后減小，總體處于增加態勢，表明研究區水質凈化能力總體趨于減弱趨勢。其中，3個時期流域尺度TN年均輸出量最大值分別為5.90、6.10、6.40 kg/(hm2·a)，總輸出量分別為840.75、880.49、877.66 t；TP年均輸出量最大值分別為0.14、0.15、0.16 kg/(hm2·a)，總輸出量分別為78.46、79.46、79.34 t。從空間上看，3個時期氮磷營養鹽輸出空間分布格局基本保持一致，高值區主要集中于龍門鄉中南部、博南鎮中部和西南部、杉陽鎮中東部、廠街鄉西部且多呈面狀分布，其余高值區以點狀分布于水泄鄉東北部、龍街鎮中南部、北斗鄉西北部。

圖5 大理州永平縣水質凈化服務時空分布特征

從土地利用類型上看，3個時期TN年均輸出量高低值依次為耕地>建設用地>草地>林地>水域；TP輸出量依次為耕地>建設用地>林地>草地>水域。從各時期來看，耕地始終為大理州永平縣氮磷營養鹽輸出總量的主要貢獻角色，且輸出量呈現出逐年增加的趨勢，這是由于農耕地所種植的作物生長長期依賴于農藥和化肥的使用，而農耕地多呈面狀分布于河谷地帶，因此成為氮、磷營養鹽污染物的主要來源；建設用地TN、TP輸出量次之，這可能是由于研究區城鎮建設及農村養殖畜牧業發展所導致的；林地、草地氮磷營養鹽輸出量相對較少，這是因為植被具有過濾和截留營養鹽的功能；而水域作為人類生存發展所需的重要資源，則不具備氮磷營養鹽輸出的能力(圖6)。

圖6 大理州永平縣各地類TN、TP輸出量

從不同鄉鎮間來看，杉陽鎮、廠街鄉、博南鎮氮磷營養鹽輸出量始終高于其余4個鄉鎮，從土地利用類型的分布格局進行分析，杉陽鎮、廠街鄉、博南鎮的地類以耕地和建設用地為主導，因此三者的營養鹽輸出量較高；北斗鄉、龍門鄉、龍街鎮的氮磷輸出量相對較少，這是因為這些區域的植被覆蓋較為完好，從而對營養鹽的輸出起到一定的過濾及凈化效用(圖7)。

圖7 大理州永平縣各鄉鎮TN、TP輸出量

3.3 產水和水質凈化服務冷熱點分析

在進行冷熱點分析之前，首先運用ArcGIS空間自相關分析判斷大理州永平縣產水量和水質凈化量在空間上是否存在聚集特征，計算結果如表2所示。其中，2000年、2010年、2020年研究區的產水量和氮磷輸出量莫蘭指數均大于0，表明兩項服務功能均具有空間相關性，同時，Z值均大于2.58，概率P值為0，均小于0.01，表明研究區產水量和氮磷輸出量的空間分布為非隨機狀態，具有聚集特征?；诖耍_展大理州永平縣產水和水質凈化服務冷熱點分析，結果如圖8、圖9所示。

表2 大理州永平縣產水和水質凈化服務功能空間自相關指數

圖8 大理州永平縣產水服務冷熱點空間分布

圖9 大理州永平縣水質凈化服務冷熱點空間分布

3個時期，研究區兩項服務的冷熱點空間分布格局特征具有相似性。其中產水服務冷熱點分析結果顯示，熱點區域主要分布于杉陽鎮、龍門鄉、博南鎮、水泄鄉，這些區域地類多以耕地、建設用地為主，冷點區域主要集中在北斗鄉，該區以林草地為主要用地類型，植被覆蓋度高，人為干擾活動較少，有利于抑制地表徑流；水質凈化服務熱點區域則主要分布于杉陽鎮、博南鎮、廠街鄉等鄉鎮間，這些區域用地多以耕地為主，氮磷輸出量較高，北斗鄉則為研究區的主要冷點區域。總體上看，20年來大理州永平縣產水服務的熱點區域從空間上看主要集中分布于西部地區，并逐漸從北部地區向南部地區發生轉移，冷點區域主要集中分布于中部和東北部地區；水質凈化服務熱點區域分布較為穩定，主要集中在西北部和東南部地區，呈現出與耕地和建設用地一致性的空間分布格局特征。

3.4 大理州永平縣水生態系統服務功能重要性等級劃分

為劃定大理州永平縣水生態系統服務功能重要性等級，本文僅以2020年產水和水質凈化服務作為劃分對象。研究運用ArcGIS將產水服務和水質凈化服務計算結果進行歸一化處理，并進行空間疊加分析，將研究區水生態系統服務功能重分類為極重要、較重要、一般重要、不重要四個等級，分別歸為極重要區、較重要區、一般重要區和不重要區[33-35]，面積分別為103 089.87、108 449.37、43 441.38、14 644.8 hm2，占總面積的比例為38.24%、40.22%、16.11%、5.43%，結果如圖10所示。

圖10 大理州永平縣2020年水生態系統服務功能重要性等級分布

由圖中各區空間分布特征可知，水生態系統服務功能極重要區主要分布在研究區中部和東部區域，用地以林草地為主，該區植被覆蓋度高，具備良好的水源涵養功能，人為干擾活動少，是繼續推進實施退耕還林還草的重要區域，同時應加強該區內水源地的保護工作，杜絕毀林開荒、毀林采礦等行為；較重要區和一般重要區主要分布在用地以灌草和旱地為主的西北部和南部區域，這些區域應加強退化草地的恢復治理力度，控制好環境污染，切實保護水源水質，減少對該區的干擾；不重要區分布在研究區用地類型以耕地和建設用地為主的區域，該區植被覆蓋度低，人為干擾活動強烈，是重點實施退耕還林還草和污染防治的區域，以提升水生態系統的水源涵養能力。

4 討論與結論

4.1 討論

本文通過運用InVEST模型產水模塊和水質凈化模塊對大理州永平縣產水量和水質凈化量時空分布變化特征進行模擬，并利用冷熱點分析工具識別兩項服務的空間分布聚集情況。研究結論與已有的研究相類似[36-39]，結果均呈現出產水服務以建設用地和耕地占據高值，而林草地占據低值的特征，這是由于快速城市化進程所改變的不透水下墊面導致的[40]。因此，當地應加強林草地種植及恢復力度，提高水生態系統水源涵養能力；水質凈化服務中，耕地的氮磷輸出量始終最高，該結果強調了農藥和化肥的大量使用是造成區域水質凈化服務功能衰退的主要因素，實施退耕還林還草工程是提高當地水質凈化能力的有效措施[41]。同時，在氮磷輸出量高值區，應多種植喜氮、喜磷植物，充分發揮植物的攔截與凈化效用，合理使用農藥和化肥，制定排污標準，降低氮磷進入水體的概率，確保當地居民飲用水安全。此外，通過模擬大理州永平縣兩項服務功能的時空變化，本文將水生態系統服務功能劃分為4個重要區，確定優先重點保護的區域，以此為當地開展水源水質保護提供科學指導。其中，極重要區植被覆蓋度高，具有良好的調節地表徑流、水源涵養能力，是當地優先考慮實施保護的區域，應加強區域內退耕還林還草工程的建設工作，以此促進水源地的健康與穩定。

4.2 結論

(1)2000年、2010年、2020年大理州永平縣的產水服務年均產水深度分別為642.37、348.95、275.59 mm，呈下降趨勢。建設用地產水量始終占據高值，其次為耕地，水域產水量最低；用地類型以耕地為主導的鄉鎮產水量最高。

(2)2000年、2010年、2020年大理州永平縣的水質凈化服務TN年均輸出量最大值分別為5.90、6.10、6.40 kg/(hm2·a)，TP年均輸出量最大值分別為0.14、0.15、0.16 kg/(hm2·a)，均呈上升趨勢。耕地氮磷輸出量最高，其次為建設用地；以耕地和建設用地為主要用地類型的鄉鎮氮磷輸出量較高。

(3)大理州永平縣產水服務熱點區域主要集中在西部地區，并逐漸從北部地區向南部地區發生轉移，冷點區域主要集中分布在中部和東北部地區；水質凈化服務熱點區域呈現出與耕地和建設用地一致性的空間分布規律，冷點區域則集中連片分布于植被覆蓋度高，較少人為活動干擾的區域。

(4)大理州永平縣水生態系統服務功能極重要區占研究區總面積的38.24%，主要分布在用地以林草地為主的中部和東部區域；較重要區和一般重要區主要分布在西北部和南部用地以灌草和旱地為主的區域，占總面積的40.22%、16.11%；不重要區占總面積的5.43%，分布在用地以耕地和建設用地為主的區域。