丹江口庫區產水量時空動態與情景模擬*

張 建 雷 剛 漆良華,

(1.國際竹藤中心 國家林業和草原局/北京市共建竹藤科學與技術重點實驗室 北京 100102； 2.國際竹藤中心安徽太平試驗中心 太平 245700)

產水量是流域重要的生態系統服務功能，可表征一定時空范圍內生態系統水分保持的過程與能力，既影響區域水資源的整體水平，又反映區域自然環境和人為活動的關系(陳姍姍, 2016； Langetal., 2018)。因此，深入了解區域產水量時空動態變化特征及其主要驅動因素具有重要意義。

隨著3S技術與水文生態模型的發展，越來越多的學者通過水文模型分析和評估區域生態系統產水量(Lehetal., 2013)。其中InVEST模型可以較好地把握總體格局，將量化的生態系統服務功能和情景以地圖的形式表達，可設置不同情景來模擬未來土地利用條件下生態系統服務功能和多種服務之間的權衡關系(劉智方等, 2016； 唐堯等, 2015)。該模型由美國斯坦福大學、大自然保護協會、世界自然資金會共同開發，輸出結果的空間化表達使生態系統服務的重要區域易于識別(Rimaletal., 2019)。InVEST模型產水量模塊被Leh等(2013)、Neson等(2009)應用于西非加納和科特迪瓦、美國俄勒岡州的威拉米特河流域； 顧晉飴等( 2018)利用InVEST模型發現太湖流域水源涵養空間分布差異性和不均衡性不斷增大，降水量對水源涵養能力的影響大于土地利用/覆被類型的影響； 程一凡(2019)利用模型評估了三江源國家公園近十年水源涵養量變化，得到不同時期不同區域水源涵養量變化特性，但對驅動因素的探討依然空缺； 邱問心等(2018)以浙江臨安區水濤莊水庫集水區為例，基于水量平衡法、綜合蓄水法和InVEST模型法評價研究區水源涵養功能，結果表明模型結果精度可達83.36%，InVEST模型測算水源涵養功能是可行的且具有較高推廣價值。

丹江口市為南水北調中線工程水源區，近年來城鄉發展、人口聚集和土地利用/覆被類型變化等影響了庫區水源涵養服務能力，因此為學者們所關注。丁霞等( 2019)對庫區馬尾松(Pinusmassoniana)林水源涵養價值進行了研究，徐慧( 2017)利用遙感技術探究了南水北調中線工程的初始運行對水源地帶來的影響，這些研究多局限于某一種生態系統類型的單一生態服務價值(李亦秋等, 2011； 柳晶輝等, 2008)，缺乏時空動態變化和驅動因素的研究。因此，本研究以南水北調中線工程水源區丹江口市為研究區域，采用InVEST模型研究庫區2003、2013和2018年產水量時空動態，并利用情景模擬方法探討年降水量和土地利用/覆被類型變化對產水量的影響，以期為庫區水源涵養能力提升與調控提供科學依據。

1 研究區概況

丹江口市地處湖北省西北部，十堰市東部，漢江中上游(110°47′53″—110°34′47″E，32°14′10″—32°58′10″N)，是南水北調中線工程重要水源地，東西最大橫距73 km，南北最大縱距81 km。南水北調中線工程于2003年開工，同年丹江口市開展移民搬遷工程； 2013年南水北調中線水源地工程完工，開展調水試驗，同年移民搬遷工程圓滿完成； 截止2018年底丹江口庫區累計供水2.22億m3。丹江口市屬亞熱帶半濕潤季風氣候區，氣候較溫和，四季分明，光照時間長，雨熱同季且無霜期較長，年均氣溫15.9 ℃，年均降水量約800 mm，年蒸發量1 979.1 mm。土壤以黃棕壤和黃壤為主，成土母質由石灰巖、片麻巖等發育而成，質地疏松。馬尾松、栓皮櫟(Quercusvariabilis)、柏木(Cupressusfunebris)、柑橘(Citrusreticulata)等在庫區均有分布。

2 研究方法

2.1 數據來源及處理

InVEST模型所需數據包括研究區土地利用/覆被類型圖、年均降水量、年均潛在蒸散量、土層厚度、植物可利用水含量、集水區范圍、季節性因子Z值等。

1) 土地利用/覆被類型圖 通過遙感影像圖進行解譯分類得到。原始遙感影像數據源于美國地質調查局(http:∥www.usgs.gov/)，選擇丹江口市2003年Landsat 7 ETM+、2013年和2018年Landsat 8 OLI_TIRS 影像，分辨率為30 m，云量均3%以下，影像質量良好，在ENVI軟件中通過幾何校正、融合鑲嵌、裁剪處理，采用監督分類方法進行解譯，將研究區土地利用類型分為林地、耕地、園地、建設用地、水域和未利用地6類(圖1)。

圖1 2003—2018年丹江口市土地利用/覆被類型Fig. 1 Land use/cover types of the Danjiangkou City in 2003—2018

2) 年均降水量 年降水數據下載于中國氣象數據網(http:∥data.cma.cn/)，為規避單一年份數據代表性低，選取3個時期(2002—2004、2012—2014、2017—2019 年)氣象數據的均值，整理研究區周邊氣象站點及監測設備降水量數據，通過克里格插值得到年均降水量柵格圖。

3) 年均潛在蒸散量 利用中國地面氣象觀測報告，收集研究區周邊4個氣象站數據。利用Hargreaves方程(Droogersetal., 2002)計算年均潛在蒸散量ETo，通過克里格插值得到年均潛在蒸散量的柵格圖。

ETo=0.001 3×0.408×RA×(T+17)×

(TD-0.0123P)0.76。

式中： RA為太陽大氣頂層輻射量；T為年最高溫均值和最低溫均值的平均值；TD為年最高氣溫均值和年最低氣溫均值之差。

4) 土層厚度 根據森林資源二類調查數據結果中的屬性數據，采用克里格插值，得到模型所需圖層厚度的柵格圖。

5) 植物可利用水含量 表征土壤為植被生長所儲蓄的總水量，參考周文佐等(2003)提出的計算公式，將計算結果進行克里格插值得到模型所需柵格圖。

PAWC=(54.509-0.132SAN-0.003SAN2-0.055SIL-0.006SIL2-0.738CLA+0.007CLA2-2.688C+

0.501C2)×100%。

式中： PAWC為植物可利用水含量；SAN為土壤的沙粒含量(%)； SIL為土壤粉粒含量(%)； CLA為土壤黏粒含量(%)；C為土壤有機質含量(%)。

6) 集水區范圍 指被分水嶺所包圍的集水流域，內部包含若干個互不嵌套的子流域。從地理空間數據云網站(http:∥www.gscloud.cn)下載研究區30 m分辨率高程數據，經鑲嵌、裁剪、填洼后生成DEM數據； 基于DEM使用ArcGIS水文分析工具和矢量數據生成流域和子流域，每一級流域賦予唯一數值。

7) 季節性因子Z值 是降水特征的常數，通常用Zhang系數表征(Zhangetal., 2001)。對于降水總量相等的區域，降水次數越多該值越大。本研究經多次調節Z值，確定季節性因子取值為11時，模型評估結果與實測值最為接近。

2.2 年產水量模型原理

InVEST模型產水量模塊是基于Budyko水熱耦合平衡假設和年均降水的水量平衡估算方法建立的(潘韜等, 2013)。模型不區分地表水、地下水和基流，某柵格單元的年降水量減去年蒸發量即年產水量，參考InVEST模型手冊(Redheadetal., 2019)，其算法為：

式中:Yxj為柵格x中土地利用/覆被類型j的年產水量； AETxj為柵格x中土地覆被類型j的年蒸散量；Px為柵格x的年降水量。

式中：Rxj為為第j類土地利用/覆被類型在柵格單元x的干燥度指數，為年蒸散量與年降水量的比值，無量綱；wx為植物年需水量和年降水量的比值，無量綱。

式中：Kxj為柵格單元x中土地利用/覆被類型j的植被蒸散量的系數； ETox為柵格單元x中的年均潛在蒸散量。

式中： AWCx為土壤有效含水量，由土壤質地和有效土壤層厚度決定；Z值為季節性因子，代表季節性降水分布和降水深度，不同研究區Z值存在差異。

2.3 情景模擬法

由于InVEST模型產水量模塊是基于水量平衡法原理開發的，年降水量和年蒸散量是影響結果的主要因素，而年蒸散量主要受植被和土地利用/覆被類型影響(潘韜等, 2013； Terradoetal.， 2014)。為了進一步探究年降水量和土地利用/覆被類型變化對產水量的影響，本研究設計2種情景，情景1： 模型均輸入2003年降水數據，土地利用/覆被數據為2013年和2018年實際數據，與2003年結果相比，即年降水量數據保持不變而土地利用/覆被類型變化，可研究2013年和2018年土地利用/覆被對產水量的影響； 情景2： 模型均輸入2003年土地利用/覆被數據， 2013年和2018年降水量數據為相應的年降水量，與2003年結果相比，可探究年降水量變化對產水量的影響。引入實際產水量總變化量WT(趙亞茹等， 2019； 楊潔等， 2020)，其正負值表示年降水量和土地利用/覆被類型變化使產水量增加或減少。

ΔWC=WO-WC;

ΔWL=WO-WL;

ΔWT=ΔWL-ΔWC。

式中：WO為實際產水量；WC為年降水量不變而土地利用/覆被類型變化的產水量；WL為土地利用/覆被類型不變而年降水量變化的產水量；WT為各模擬的總變化量，即在土地利用/覆被類型變化和年降水量變化雙重影響下產水量的變化量； ΔWC為實際情況與年降水量不變情景差值； ΔWL為實際情況與土地利用/覆被類型不變情景的差值。

2.4 模型檢驗與校正

土地利用/覆被圖是模型運行的數據基礎，其精度直接影響結果的輸出，采用Kappa分析方法分別對3期土地利用/覆被圖結果進行精度驗證。

在模型其他參數確定的情況下，季節性因子Z值是影響模型準確性的重要參數。將產水量模塊運行結果與實際監測數據進行擬合，通過調整季節性因子Z值大小，確定最佳值。選擇《十堰市水資源公報(2013年)》中丹江口市2013年產水量為標準數據，與模型計算的產水量進行檢驗校正。

3 結果與分析

3.1 模型校驗

在對土地利用/覆被類型的校驗中，利用谷歌地圖獲取3個時期的高清歷史影像圖，結合野外調查數據和歷史影像圖建立樣本區，對監督分類解譯的數據修正，經Kappa系數檢驗，3期土地利用/覆被圖精度達到0.84、0.86、0.85，解譯結果符合精度要求。

將所需數據按要求輸入模型后，在其他數據已確定情況下，其運行結果受季節性因子Z值的影響。將年產水量模型運行結果與監測數據進行擬合，結果發現隨著季節性因子Z值增加，產水量呈減少趨勢(圖2)，確定季節性因子Z值為11時，模型結果與監測值最為接近，誤差為0.72%。

圖2 模型校驗結果Fig. 2 Model verification result

3.2 產水量時空動態

整體上看(圖3)，丹江口市近15年產水量呈先減弱后增強的趨勢， 2003年單位面積產水量最深可達646.96 mm，遠高于2013年(451.25 mm)和2018年(527.53 mm)。2003年丹江口市產水深度310.09 mm，產水量為13.6億m3，產流系數為0.324； 2013年產水能力最弱，產水深度146.67 mm，產水量為5.25億m3，產流系數0.281； 2018年產水量相較2013年增加了27.88%，產水量為7.28億m3，產水深度為166.06 mm，產流系數0.296。

圖3 2003、2013和2018年丹江口市產水深度時空分布Fig. 3 Spatial and temporal distribution of water yield depth in Danjiangkou in 2003, 2013， and 2018

同一年內，丹江口市產水能力在空間上分布不均勻，北部丘陵山區產水深度高出南部39.84%～57.79%； 時間序列下， 2003—2013年，全市面積的85.10%為減弱區，13.50%的區域表現為增加趨勢(主要分布在林地建設區)； 2013—2018年，丹江口市89.60%的區域產水深度有所增加，而10.20%的區域表現出減少的趨勢(圖4)。

圖4 2003—2018 年丹江口市產水深度變化Fig. 4 Variation of water yield depth in Danjiangkou from 2003 to 2018

下墊面是影響產水量變化的重要因素，利用分區統計功能得到不同土地利用/覆被類型下的產水深度變化情況(圖5)。丹江口市4種土地利用/覆被類型中，產水深最高的是建設用地(530.35 mm)，其次分別為耕地(292.85 mm)、園地(284.07 mm)和林地(273.76 mm)； 相較于其他土地利用/覆被類型，建設用地無植被截留降水，蒸散量小，因此產水深度最好。從變化趨勢上可看出， 2003—2018年各土地利用/覆被類型的產水深度呈先減小后緩升的趨勢； 總體而言，耕地保持增長而其他土地利用/覆被類型則降低。

3.3 情景模擬分析

總體上看(圖6)，在年降水量不變而土地利用/覆被類型變化情景下，情景1中2013和2018年產水量分別為10.5億m3和8.6億m3，表明土地利用/覆被類型變化是導致產水量減少的重要因素； 在土地利用/覆被類型不變而年降水量變化情景下，情景2中2013年和2018年產水量變幅較小，分別為5.68億m3和9.7億m3，說明年降水量是影響產水量變化的重要因素。

圖5 不同土地利用/覆被類型產水情況Fig. 5 Water yield of different land use types

表1表明了年降水量變化和土地利用/覆被類型變化對產水量的影響程度。2013年，兩種情景模擬結果顯示，土地利用/覆被類型變化影響下實際產水量減少了0.43億m3，年降水量變化使實際產水量減少了5.25億m3，二者的貢獻率分別為7.57%和92.43%，說明2003—2013年降水量與土地利用/覆被類型變化均使得產水量減少，且年降水量是影響產水量的主要因素。而2018年，土地利用/覆被類型變化導致實際產水量減少了2.42億m3，貢獻率為64.71%，年降水量變化使實際產水量減少了1.32億m3，說明相較于年降水量，土地利用/覆被類型變化對2013—2018年產水量的影響更大。

圖6 模擬不同情景下產水量變化Fig. 6 Simulate the change of water yield under different scenarios

表1 年降水量與土地利用/覆被類型變化對產水量的影響Tab.1 Annual precipitation and land use/cover type change influence water yield

在僅考慮年降水量不變而土地利用/覆被類型變化情況下， 2013—2018年水域和林地面積持續增加，特別是林地面積增加了10.03%，耕地和園地面積不斷減少，建設用地先減少后增加，未利用地變化不明顯(表2)。由于建設用地基本為不透水面，耕地和園地的植物密度、根系深度及人為勞作帶來的地表變化等都有助于產水量的增加，而林地通過冠層截留降水、枯落物層吸收降水、土壤層蓄滲降水以及其較深的根系可有效攔截降水，且有強大蒸騰作用，從而使得產水量較少。

表2 丹江口市2003-2018年不同土地利用/覆被類型面積變化Tab.2 Area changes of different land use/cover types in Danjiangkou from 2003 to 2018

4 討論

Lü等(2013)研究了中國東南部南北樣帶1981—2000年產水量和氣候及土地利用/覆被類型變化之間的動態響應，其中產水量的時空變化受年降水量影響最顯著，與本研究結果相似； 王耕等( 2018)利用年產水量數據進一步對大凌河上游匯水區水源涵養功能定量評價，發現產水量與水源涵養均有明顯的空間異質性； Wang等(2020)評估了時空變化下丹江口庫區景觀格局和土地利用/覆被對水資源的影響，對未來土地利用管理和政策的實施具有重要意義。當前產水量研究方法中均存在一定的優勢和局限性，而InVEST模型考慮了蒸散發、降水量和多年徑流系數等數據，研究結果更為準確。

本研究利用InVEST模型對丹江口市2003—2018年產水量開展了研究，結果表明近15年產水量時空變化較大，產水量分布不均，且各土地利用/覆被類型的產水深度存在差異。研究結果與太湖流域(顧晉飴等， 2018； 陳駿宇等， 2016)、三江源流域(程一凡， 2019； 呂樂婷等， 2020)大凌河流域(王耕等， 2018)、石羊河流域(趙亞茹等， 2019)、大連市(呂樂婷等， 2019)等流域或地區的結果相似，其中楊潔等(2020)研究了黃河流域不同土地利用/覆被類型的產水深度，發現永久性冰川雪地、裸巖石質地、城鎮建設用地的產水深最高，與本研究結果中建設用地產水量最高相同。

產水量大小受該區域年降水量、蒸散量以及二者之間平衡關系的影響，而實際蒸散量除了受氣象因素影響外，還直接受制于土地利用/覆被類型的影響。Lang等(2017)量化了氣候變化和土地利用/覆被對產水量的影響，發現年降水量對年產水量變化的貢獻率高達97.44%，與本研究情景1中的模擬結果相似。土地利用/覆被類型變化通過改變下墊面結構和類型直接或間接影響產水量，其次是土壤孔隙度、土壤質地、結構均對產水量帶來間接影響(Jatinetal., 2018； Pamukcuetal., 2016)。因此產水量變化是多因子耦合驅動的過程(Lüetal., 2017； Suetal., 2013)，集合了包括氣候變化與土地利用/覆被類型變化在內的自然、社會、經濟等多方面因素的共同作用。

本研究忽略了上游對該區域的水源補給，可能會導致InVEST模型擬合結果偏大； 季節性因子Z值最終設定為11，但由于小流域氣候等因素也會存在一定誤差； 模型結構的本土化不足以及未能獲得長期野外觀測數據，會對研究結果的準確性產生影響。在今后研究中，應加強對野外數據的監測收集、調整模型參數使其本土化，確保評估結果的可靠性。

5 結論

本研究基于InVEST模型，以丹江口市域范圍為研究區，定量研究了2003、2013和2018年產水量的動態變化特征，并運用情景模擬法探討了年降水量和土地利用/覆被類型變化對產水量的影響。結果表明： 1) 丹江口市2003—2018年產水量表現出先減弱后增強的趨勢；2) 單位面積產水量在空間上分布不均，平均而言北部丘陵山區高出南部39.84%～57.79%；3) 各土地利用/覆被類型產水深度均有差異，其中產水深度最高的是建設用地(530.35 mm)；4) 情景模擬表明年降水量和土地利用/覆被類型是影響產水量的主要因素， 2013年產水量中年降水量貢獻率為92.43%，土地利用/覆被類型變化對2018年產水量貢獻率為64.71%。