張承水源涵養區土地利用演變及產水量的響應

張 薇，王鳳春，3，4，賈 悅，萬虹麟，5，牛天浩云，杜雨森

（1.河北水利電力學院水利工程系，河北滄州 061001；2.滄州市遙感與智慧水利技術創新中心，河北滄州 061001；3.河北省高校水利自動化與信息化應用技術研發中心，河北滄州 061001；4.河北省數據中心相變熱管理技術創新中心，河北滄州 061001；5.河北省巖土工程安全與變形控制重點實驗室，河北滄州 061001）

0 引 言

近年來隨著水資源短缺問題的日益突出和生態系統服務功能研究的不斷深入，水源供給作為生態系統服務的重要功能之一，被廣泛關注［1］。生態系統的水源供給功能概念較廣，其內涵隨著研究的不斷深入而不斷完善，主要表現為生態系統攔蓄降水、調節徑流、影響降雨、凈化水質等，研究表明不同的土地利用/覆被（LULC）類型通過影響流域水循環過程中的蒸散、下滲和持水，進而改變流域水源供給服務功能［2］。目前，眾多學者在水源供給服務功能方面做了大量研究，取得了系列成果，為相關研究提供了有益參考［3-6］。特別是InVEST 模型，相較于其他水資源評估模型以其原理簡單，易于被不同學科背景的使用者理解，且需要的數據量相對較少、易于獲取，近年來被廣泛應用于水源供給服務功能評價相關研究中［7-9］。

密云水庫上游流域張承地區作為密云水庫的水源涵養區，擔負著向北京市提供飲用水源的重要任務，其水源供給服務功能質量將直接關系首都的水生態安全。近年來，隨著經濟的發展和人類活動干擾的加劇，流域生態環境脆弱、庫區水生態質量下降。其中涉及水質評價、污染源及污染評價、流域水生態改善機制等方面的研究較多，但分析產水量時空變化及其與土地利用變化的響應關系方面的研究相對較少，且不同學者研究的流域范圍也不盡相同［10-12］。從土地利用結構驅動水源涵養生態服務功能的角度出發，定量評估區域產水量并分析其時空變化特征，揭示其變化與土地利用類型變化的內在關系，深入分析哪種土地利用方式有利于提高區域水源供給能力，有助于實施精準靶向生態補償，提高補償效率和效果。

基于此，本文借助InVEST 模型產水模塊，定量評估張承水源涵養區2000-2019年的產水量并分析其時空變化特征，進一步探討不同土地利用類型變化與產水量變化的響應關系，以期為張承水源涵養區的國土規劃和京冀生態協同發展提供科學依據。

1 數據來源與研究方法

1.1 研究區概況

密云水庫有白河和潮河兩大入庫河流。其中，白河起源于張家口沽源縣，經赤城、延慶、懷柔入庫；潮河源于承德豐寧縣，經灤平、古北入庫。本文研究區為密云水庫上游流域的河北省赤城、豐寧和灤平三縣，位于40°30′～41°37′N，115°30′～117°32′E，總面積約111.30 萬hm2（圖1）。該區域地處山地與平原的過渡地區，地勢北高南低。該區域屬于中溫帶向暖溫帶過渡、半干旱向半濕潤過渡的中緯度大陸性季風氣候，年平均氣溫在7 ℃左右，且隨海拔的上升有所下降。全年降雨量約為490 mm，且大多數集中在7-9月份，降雨所形成的地表徑流是河流的主要補給形式。人工林和天然次生林是區域內的主要植被類型，同時，灌叢有較廣的分布。

圖1 研究區范圍Fig.1 Scope of study

1.2 研究方法

（1）土地利用轉移矩陣。引進土地利用轉移矩陣來定量描述土地利用類型轉移情況。轉移矩陣是將各地類變化的轉移面積按矩陣形式排列，定量分析不同時期土地利用的變化結構和方向，直觀反映一定時期內各地類間的轉移方向和數量［13］，其數學模型為：

式中：U為土地利用轉移矩陣；Uij表示轉移前的i地類轉換成轉以后的j地類的面積。

（2）InVEST 模型產水模塊。InVEST 模型產水模塊基于Budyko 水熱耦合平衡假設（1974）和年平均降水量數據，年產水量等于年降雨總量與年蒸散總量的差值［14］。該模型運行基于柵格地圖，確定研究區每個柵格單元的年產水量［15］，公式如下：

式中：Yx，j為j類土地利用類型柵格x的產水量；Px為柵格x的年降水量；AETx，j為j類土地利用類型x的年實際蒸散量［16］；Rx，j是土地利用類型j柵格x的布德科干燥度指數，是潛在蒸散與降水的比值；Wx是一個非物理參數；AWCx是植被可利用的體積含水量，mm；Z為Zhang 系數，為一個經驗值，其值在1～10，根據實際情況調整［17］。

式中：ET0為柵格x內的潛在蒸散量；RA為太陽大氣頂層輻射；Tagv是日最高溫平均值與日最低溫平均值的均值；TD是最高溫平均值與日最低溫平均值之差；P為柵格x的年降水量［18］；AWC植被可利用的體積含水量；SAN%土壤沙粒含量百分比；SIL%土壤粉粒含量百分比；CLA%土壤黏粒含量百分比；C%土壤有機碳含量百分比。

（3）情景分析法。InVEST模型基于水量平衡原理模擬產水量，降水和實際蒸散量是影響模型模擬結果的主要因素，而蒸散量主要受植被和土地利用類型影響。為進一步探究產水量與土地利用/覆被變化的關系，這里設計了模擬情景：在2000年模型數據的基礎上，只改變土地利用覆被數據，分別輸入2010年和2019年土地利用數據和對應的作物生物物理參數，可研究2000-2010 和2000-2019 土地利用覆被變化對產水量的影響［19］。

式中：Rl土地利用覆被變化對研究區產水量變化的貢獻率；l0基準產水量；li模擬產水量。

1.3 數據來源及預處理

選取2000年、2010年、2019年3 期數據為研究時序?；谏鲜龇椒ǎ占搜芯繀^的流域界限、行政區劃、土地利用、氣象、土壤屬性、DEM 等數據用于模型計算，同時還整理了研究區歷年水資源量數據對模型結果進行驗證。其中，土地利用數據主要通過遙感解譯獲得，遙感數據由地理空間數據云（http：//www.gscloud.cn/sources/）下載；氣象數據來源于國家氣象信息中心中國氣象數據網（http：//data.cma.cn/），包括降雨、氣溫、日照等數據；土壤數據來源于世界土壤屬性數據庫（HWSD）；DEM數據也來源于地理空間數據云；水資源量來源于《河北省水資源公報》（http：//slt.hebei.gov.cn/）。

土地利用數據由Landsat OLI/ TM 遙感影像，利用ENVI5.3通過監督分類結合人工目視解譯的方法獲取，參照中國科學院資源環境數據中心土地分類系統，將土地利用類型分為6 個一級類：耕地、林地、草地、建設用地、水域、裸地，其綜合精度達87%以上；降雨量空間數據由河北省13 個站點實測數據，經ArcGIS10.2空間插值、裁剪所得；小流域劃分空間數據理論上應該由DEM 經ArcGIS10.2 投影變換、裁剪、水文分析流域提取所得，但這里由于研究區的不是一個完整的流域，且為了后期分析統計方便，小流域劃分空間數據由行政區劃矢量數據代替；潛在蒸發量空間數據是由氣象站點實測數據，經公式計算后進行空間插值、裁剪得到；植物有效含水率基于土壤屬性數據（沙粒、粉粒、黏粒、有機碳含量）經ArcGIS10.2 屬性字段計算得到；土壤深度數據由全國土壤深度數據經ArcGIS10.2 裁剪得到（見圖2）。

圖2 基礎處理數據（2010年6月）Fig.2 Processed basic data（2010，Jun）

2 結果與討論

2.1 土地利用結構時空變化分析

由遙感解譯得到2000年、2010年和2019年3 期的土地利用分布圖（見圖3）。分析可知：林地的分布最為廣泛，其次是草地、耕地和建設用地，耕地分布在河流沿岸地勢平坦地區，草地則沿河流走勢分布在農業生產用地的外圍，林地則主要位于離河較遠的丘陵、山地區域。

圖3 2000-2019年研究區土地利用空間分布Fig.3 Spatial distribution of Land use of the study area in 2000-2019

經統計得2000-2019年各地類面積及變化情況，如表1 所示。張承水源涵養區土地利用結構中林地和草地是主要的用地類型，2019年兩種用地類型面積分別為793 461.70 和138 779.25 hm2，各占總面積的71.3%和12.5%，作為密云水庫上游流域涵養區，林地、草地有利于水量的儲存和凈化。2000-2019年的土地利用結構轉型中，耕地、裸地和林地總體減少；草地、建設用地和水域總體增加，其中建設用地面積增長較快，19年間增加了4 268.32 hm2，年均漲幅2.1%，表明隨著城鎮化水平的不斷提高，研究區建設用地需求量增加。

表1 2000-2019年各地類面積及變化 hm2Tab.1 2000-2019 Area and change of each category

進一步分析各土地利用類型相互轉化的方向和數量，在圖3 的3 期土地利用分類圖的基礎上，利用ArcGIS 的空間分析功能對不同時期的土地利用圖進行空間分析，獲得研究區2000-2010年和2010-2019年的土地利用功能類型的轉移矩陣（表2、3）。結果總體表現為建設用地面、草地面積增加，林地、耕地面積減少。2000-2010年，耕地、林地、建設用地面積均有所增加，增加的面積大多都是由水域、草地轉化而來，其中草地分別向耕地、林地、建設用地轉化了1 115.59、5 773.18 和92.08 hm2，裸地有86.3%開墾成了耕地。2010-2019年，林草地面積之和增加，耕地和裸地面積均有減少，且減少的耕地主要轉化為林草地、水域和建設用地，分別為34 655.8、384.68 和4 231.44 hm2；建設用地持續增加，林草地、耕地、裸地、水域均不同程度向其轉化，轉化率分別為：0.8%、2.4%、2.2%和1.5%；水域面積有所增加，增加的面積由林地、耕地、建設用地轉化而來。研究區作為密云水庫的水源涵養區，長期以來始終以保護生態環境作為地區發展的首要條件，充分利用地區林地、草地面積較多的優勢，優化區域水源涵養生態服務功能，按照《京津冀協同規劃綱要》的相關規定，多措并舉推進生態清潔小流域建設。特別是隨著《河北省張承地區生態保護和修復實施方案》等一系列政策的落地實施，對研究區土地利用結構變化起到一定的促進作用。

表2 2000-2010土地利用變化矩陣 hm2Tab.2 2000-2010 Land use change matrix

表3 2010-2019土地利用變化矩陣 hm2Tab.3 2010-2019 Land use change matrix

2.2 產水量時空分布及變化特征

InVEST 模型在估算年產水量時除了數據預處理中所列的數據外，還需要輸入Zhang系數和生物物理參數表。Zhang系數是InVEST 模型需要輸入的重要參數，取值范圍1～10，可以由研究區年降雨發生次數進行初步估算，再結合模型檢校數據進行反復調試。生物物理參數表至少包含各土地利用/覆被類型的植被蒸散系數Kc和最大根系深度兩項參數，其中植被蒸散系數Kc是作物實際蒸散量與參考蒸散量的比值，經多次試驗，發現研究區模型產水量的預測結果對該指標敏感性較高。而以往的研究中Kc的確定多是直接采用聯合國糧農組織作物參考值［20］，由于植被蒸散系數受植被的生長狀況、氣候及水分條件等諸多因素的影響，最好結合當地實際情況進行校正。目前基于遙感估算作物蒸散系數Kc的研究相對較多［21，22］，在韓文霆、李霞［23，24］等人研究成果基礎上，結合楊潔、吳瑞［25，26］等人的研究確定本文3 期生物物理參數如表4 所示。以河北省水資源公報數據進行模型檢校，經過反復調試當3 期Zhang 系數分別為1.1、3.5、3.0 時，模型模擬產水量與實際產水量最接近，相對誤差均控在20%以內（表5）。

電裝公司是汽車零部件及系統的全球知名供應商，于1949年在日本成立，目前在全球30多個國家和地區,設有180多家關聯公司，擁有員工約14萬名。電裝公司在環保、安全、舒適和便利等廣泛領域中，推進技術開發，并提供眾多產品。2003年，電裝（中國）投資有限公司正式成立，并在長春、天津、上海、廣州、武漢、濟南和重慶設有分公司。

表4 InVEST模型的各土地利用類型生物物理參數Tab.4 Biophysical parameters of land use description in InVEST model

表5 InVEST年產水模型結果校驗Tab.5 Verification of result about InVEST water yield model

2000-2019年，研究區產水深度在7.2～65.1 mm。2000、2010 和2019年的總產水量分別為2.61、3.46 和2.71 億m3，說明研究區產水量先增加后減少（圖4）。不同時期的產水量表現出空間分異格局，整體變化趨勢為產水量由西高東低逐漸變為東高西低。2000年產水主要集中在赤城縣和豐寧縣西北部；2010年赤城中東部和豐寧西部地區產水量減少，豐寧東部和灤平縣地區產水量增加；2019年研究區產水量呈東高西低的趨勢更加明顯。這種格局與研究區年均降水量和土地利用類型分布有直接關系，降水量多且植被蒸散量低的區域產水能力強，反之產水能力較弱。

圖4 2000-2019產水深度空間分布Fig.4 Spatial distribution of the water yield depth in 2000-2019

由圖5 可以看出，2000-2019年，研究區產水量先增加后減少，總產水量略有增加。2000-2010年增加區域主要集中在灤平縣和豐寧縣東北部，以及赤城縣的西北部，增加的面積占總面積的76.5%；2010-2019年產水量總體減少，其中赤城縣、豐寧縣中部和灤平減少最多。2000-2019年總產水量略有增加，其中以巴克什營鎮、澇洼鄉和兩間房鄉產水深度增加最多，分別增加了28.5、31.5 和29.4 mm；產水量減少的區域主要集中在赤城縣和豐寧縣西部地區。2000-2019年研究區的產水量變化與不同時期的降雨和土地利用結構有直接的關系。

圖5 2000-2019產水深度變化的空間分布Fig.5 Spatial distribution of the water yield depth change in 2000-2019

為了定量描述3 個區縣2000-2019年產水量的變化情況，分別對赤城、豐寧、灤平3個縣域三期的總產水情況進行了統計分析（表6）。2000-2019年3 個縣域內產水整體呈現先增加后減少趨勢，跟整個研究區三期產水趨勢吻合。其中，赤城縣2000-2019年產水總量總體減少0.177 億m3，豐寧縣和灤平縣2000-2019年產水總量總體分別增加了0.005和0.277 億m3。

表6 2000-2019各縣模擬產水量及變化 億m3Tab.6 Water yield and change in each county of 2000-2019

2.3 不同土地利用類型的產水深度和產水量

InVEST 模型的產水模塊是一種基于水量平衡法的估算方法，某柵格單元的降水量減去實際蒸發量即為該柵格的產水量，即總體的產水能力與降雨量成正比，與蒸散發量成反比，與坡度高程等地形因子無關。由于不同土地利用類型的蒸散發能力、土壤含水量、凋落物持水能力及冠層截留量存在差異，導致同一時期不同土地利用類型的產水深度不同。2000-2019年研究區的各地類平均產水深度如圖6所示，分析可知：不同土地利用類型平均產水深度從大到小依次為：建設用地、裸地、水域、耕地、林地、草地，各地類單位面積的平均產水深度分別為：30.786、27.638、27.564、27.304、26.572 和26.547 mm。其中建設用地植被覆蓋較少，相對蒸散發量也就較少，再加之不透水面的增加，使得降水滲入地下較少，導致該地類產水量最高；草地和林地地表調落物攔截地表徑流，延遲降雨匯流時間，增加土壤入滲量，同時植被冠層蒸散量相對較強，導致林地和草地產水量相對較少。

各地類的產水能力和面積占比是影響該地類總產水量的主要因素，通過圖6 分析除建設用地外，其他5 種地類平均產水能力相差不是很大，這樣研究區各地類總產水量主要取決于各地類面積占比。研究區主要的用地類型是林地，2019年，林地占整個研究區的71.3%，產水量占產水總量的72.0%。

圖6 不同土地覆被類型平均產水深度的時間變化Fig.6 Temporal variation of average water yield depth of different land cover types

2.4 產水量與土地利用相關性

由InVEST 模型產水模塊的估算原理可知，某柵格產水量主要受降雨和蒸散發兩個因素的影響，由于不同土地利用類型的蒸散發能力不同，且不同地類的凋落物持水能力及冠層截留量也存在一定的差異，導致不同土地利用類型將直接導致產水量的變化。本文采用情景模擬法，在2000年模型數據的基礎上，分別輸入2010年和2019年土地利用數據和對應的作物生物物理參數，進行情景模擬得2010年和2019年研究區總產水量，并分析兩期各地類的平均產水深度，以探究研究區產水量與土地利用的相關性。結果顯示2010年和2019年研究區總產水量相較2000年實際產水總量2.61 億m3分別減少1.42 億m3和增加0.15 億m3，地表覆被變化對產水量變化的貢獻率分別為-54.4%和5.7%。對比分析對應的生物物理參數，2019年各地類作物系數與2000年的相同，而2010年各地類作物系數較2000年略有增加，作物系數代表各類地物的蒸散發能力，進一步說明該產水模型對各地類蒸散系數較敏感，且產水量與作物系數成反比。通過對比分析2019年各地類模擬產水深度，如圖7 所示，發現與2000年相比，2019年各地類平均產水深度均有所增加，其中水域平均產水深度增幅最大由2000年的22.5 mm增加到2019年的27.8 mm，林地平均產水深度增幅最小由2000年的23.5 mm 增加到2019年的24.7 mm。通過情景模擬分析說明土地利用類型之間的相互轉移可能會導致產水量的增加或減少，例如在其他條件不變的情況下，建設用地增加會增加產水量，而林地增加會減少產水量。

圖7 模擬情景下各地類平均產水深度變化Fig.7 Change of average water yield depth in different regions in simulated situations

圖9 2000-2019產水量變化分布Fig.9 Spatial distribution of water yield and change in 2000-2019

為了進一步探討區域土地利用與產水量之間的時空分布關系，將2019年模擬產水量與2000年實際產水量變化以及兩期的土地利用變化做時空對比分析，如圖8、9 所示。通過對比分析發現：赤城縣的西南部、豐寧縣的中部以及灤平縣的南部產水量明顯增多，該地區主要土地利用類型變化為耕地或草地轉為建設用地，林地或草地轉為耕地，轉換后的土地利用類型相較耕地產水能力強，因此該地區整體產水量增加；產水量減少的4 個鄉鎮：赤城縣的白草鎮、豐寧縣的胡麻營鄉、灤平縣的鄧廠滿族鄉和付家店滿族，主要土地利用類型變化為耕地、建設用地、林地轉為草地，轉化后的土地利用類型較轉化前的地類產水能力較弱，導致該區域產水量減少。通過土地利用與產水量時空變化關系的分析，再次證明了在6種土地利用類型中，建設用地產水能力最強，其次是耕地，林地和草地產水能力最弱，蓄水能力較強。

圖8 2000-2019土地利用變化分布Fig.8 Spatial distribution of the land use change in 2000-2019

3 結 語

通過遙感解譯張承水源涵養區2000年、2010年和2019年三期土地利用數據，用InVEST 模型產水模塊從空間上量化評估了2000-2019年研究區的產水量，同時模擬不同土地利用變化對研究區產水功能的影響，最終借助GIS 空間分析功能分析了研究區產水功能的時空分布特征及其變化，并進一步探討了產水功能變化與土地利用類型變化的內在關系，通過研究分析得出以下結論。

（1）2000-2019年，研究區土地利用整體呈現耕地面積減少，林地、草地、建設用地不斷增加的特點。將其劃分為2000-2010年和2010-2019年兩個階段，前期水域和草地減少，轉化為耕地、林地和建設用地；后期耕地和林地減少，轉化為草地、水域和建設用地。

（2）2000年、2010年和2019年研究區的總產水量分別為2.61、3.46 和2.71 億m3。不同時期的產水量表現出空間異質性，整體變化趨勢為由西高東低逐漸變為東高西低。2000-2010年產水量增加區域占總面積的76.5%，主要集中在灤平縣和豐寧縣東北部，以及赤城縣的西北部；2010-2019年產水量總體減少，其中赤城縣、豐寧縣中部和灤平減少最多。

（3）研究區不同土地利用類型平均產水深度從大到小依次為：建設用地、裸地、水域、耕地、林地、草地，通過情景模擬分析發現產水模型對各地類蒸散系數較敏感，在蒸散系數不變的情況下，土地利用由產水能力低的地類轉為產水能力高的地類會造成區域產水量的增加，反之減少。

（4）采用了土地利用的6個一級分類進行了產水量分析，但由于水源涵養功能的涵義較廣，未來有必要將土地利用分類進一步細化，進一步考慮區域土壤的滲透性、地表徑流系數等因素，以便更好地分析土地利用和人類活動對張承地區水源涵養功能的影響。

（5）未來應進一步校驗模型所需的各類參數，精確研究區域的本地化參數，提高模型精度。