基于InVEST模型的疏勒河上游產水量時空變化特征分析

魏培潔,吳明輝,賈映蘭,高雅月,徐浩杰,劉章文,陳生云1,,6,*

1 中國科學院西北生態環境資源研究院 冰凍圈科學國家重點實驗室 疏勒河源冰凍圈與生態環境綜合監測研究站,蘭州 730000 2 中國科學院大學,北京 100049 3 祁連山國家公園國家長期科研基地,西寧 810000 4 蘭州大學草地農業科技學院草地農業生態系統國家重點實驗室,蘭州 730020 5 中國科學院西北生態環境資源研究院 黑河上游生態水文試驗研究站,蘭州 730000 6 青海師范大學高原科學與可持續發展研究院 地理科學學院,西寧 810008

生態系統服務是指人類從各種生態系統中獲取的所有惠益,由生態系統的支持功能、供給功能、文化功能和調節功能及其相互作用形成[1—2]。生態系統功能是生態系統服務的基礎和物質保障,其對人類生存和社會發展至關重要,特別是與水資源相關的功能[3]。隨著全球氣候變化、地表水資源短缺和水環境問題加劇,空間定量化、可視化評估區域產水量已成為當前研究的熱點議題之一[4—8]。產水量是流域特別是干旱地區重要的調節功能,它對工業、農業、漁業、人類消費和娛樂活動至關重要[9—10]。一方面,產水量是人類生存和區域經濟可持續發展的重要物質基礎；另一方面,產水量與區域自然地理條件和人類活動密切相關[10]。因此,研究區域產水量變化的驅動因素變得十分必要,其可為區域水資源的有效管理和可持續發展提供科學依據[10—11]。

自1970年以來,遙感技術和水文模型發展迅速,越來越多的學者嘗試通過模型模擬方法來定量化、可視化和精細化地評估和分析區域產水量,如MIKESHE(MIKE System Hydrological European)模型、TOPMODEL模型(TOPgraphy based hydrological Model)、SWAT模型(Soil and Water Assessment Tool)和InVEST模型等[12—13]。InVEST是一款由美國斯坦福大學、世界自然基金會、大自然保護協會和其他相關機構聯合開發用于生態系統服務功能評估的模型方法[14]。InVEST模型產水量模塊基于水量平衡原理,充分考慮不同土地利用類型下土壤滲透性的空間差異及地形、地表粗糙程度等因素對徑流影響,以柵格為單元定量估算不同土地利用類型的產水量[12—13]。其模型由于參數調整靈活、評估結果空間表達性強,且能夠實現時空間和多目標權衡[15],被國外許多國家和地區廣泛應用。如：Leh等[16]、Marquèsa等[17]和Redhead等[18]利用InVEST模型分別對西非科特迪瓦、西班牙東北部弗蘭科利河流域和英國多條河流的產水量進行了評估,并取得了很好的應用效果。此外,中國學者利用InVEST模型的產水量模塊分別對北京山區[19]和貴水流域[20]、陜西省商洛市五大水系流域[21]、三江并流區(云南、四川和西藏的交界處)[22]以及高原湖泊[23]等區域的產水量及其時空分布特征進行了科學評價,模型在各研究區均取得了很好的應用效果。同時,InVEST模型已成功用于青藏高原地區。其中,潘濤等[7]和呂樂婷等[24]分別評估了三江源區1981—2010、1980—2016年產水量的時間變化特征與空間分布狀況。王雅等[25]和趙亞茹等[10]分別對位于祁連山區中部的黑河流域、東部的石羊河流域的產水量空間分異特征進行了分析,這些研究對西北地區水資源規劃和管理都具有重要意義。

目前,氣候變化和土地利用變化是影響產水量的關鍵因素[26]。其中,氣候變化在宏觀尺度上可以通過改變流域的降水量和蒸散量來影響產水量[27—28],土地利用可通過改變流域下墊面間接影響產水量[29—30]。一些研究已探究了氣候變化和土地變化對流域產水量的影響。例如：Pessacg等[9]利用降水數據集評估了楚布特河流域的產水量,表明降水量與產水量之間存在正相關關系；Stone等[31]重點研究了氣候變化對密蘇里河流域產水量的影響,證實氣候變化對產水量影響顯著；Zhang等[32]和Li等[33]通過分析土地利用變化對產水量影響,得出不同土地利用方式的改變對產水量的影響不同。目前,全球氣候變化對中國水資源的影響日益加劇[34]。青藏高原作為對全球環境變化特別是氣候變化響應最敏感的地區之一,近年來受氣候變暖影響,50年來氣溫顯著升高0.3 ℃/(10a),部分區域降水波動較大[35],并且氣候暖濕化促進了植被生長,NDVI顯著增加[36]。然而,青藏高原地區氣候變化和土地利用變化對產水量的相關研究仍然缺乏,尤其植被根系深度和NDVI對產水量的影響更是鮮有研究。

水在維持流域生態系統平衡和生態承載力可持續性方面扮演著極其重要的角色[37]。干旱區水資源空間分布不均衡,降水多集中于上游山區,上游的出山口徑流量基本代表著流域的總水資源[38]。疏勒河流域地處青藏高原東北緣祁連山西段,是我國河西走廊內陸干旱區三大內陸河流域之一,也是河西農牧民群眾的“生命線”和“天然水塔區”[39]。然而,疏勒河流域產水量相關的研究報道甚少。開展上游產水量研究,對整個流域水資源的動態評估、生態平衡的維系以及氣候調節具有十分重要的意義[10,15]。鑒于此,基于InVEST模型產水量模塊,本研究對疏勒河上游高寒山區產水量時空變化特征進行了分析,并進一步探討其對環境要素變化的響應,旨在為流域水資源的有效管理和可持續發展提供科學參考。

1 研究區概況

疏勒河流域地處我國西北干旱區腹地,是河西走廊僅次于黑河的第二大內陸河[40—41],其水源補給主要源于冰川和凍土融水以及山區降水[37,42]。疏勒河上游(96°37′12″—98°59′24″ E,38°13′12″—39°52′12″ N)是指以昌馬堡水文站(96°51′ E,39°49′ N)為出水控制點的上部區域,為疏勒河干流的形成區[42],流域面積約為10973.9 km2,海拔2031—5763 m(圖1),多年平均徑流量為13.81×108m3。地形上主要由疏勒南山、托來南山和疏勒河谷地組成,山地地形陡峭,河谷地形相對低緩[38]。氣候類型為典型的大陸性氣候,按陸地表層干濕指數區劃屬于半干旱區—干旱區邊緣地帶,呈現干寒氣候條件下的草原—荒漠景觀[38—40]。土壤類型主要為高山寒漠土、高山草甸草原土、栗鈣土、淡栗鈣土和山地灰鈣土[40],土壤質地主要以壤土為主,砂土和粘土占比均較小。疏勒河上游是整個流域的水源涵養區和產流區,為中下游山前平原和綠洲地區賴以生存和發展提供了水源保障。

圖1 研究區位置和水文站點分布圖Fig.1 Location of study area and hydrological station

2 研究方法與數據處理

2.1 產水量模型

InVEST產水量模塊基于水量平衡原理,用各柵格單元的降水量減去實際蒸散后的水量,進而得到該柵格產水量。該模型假設每個柵格單元所有的產水量通過地下徑流或地表徑流的方式到達流域出水口[14]。模型主要算法如下[7—8,43]：

(1)

式中,Yxj為第j土地利用類型柵格x的產水量(mm),AETxj為第j土地利用類型柵格x的年實際蒸散發量(mm),Px為柵格x的多年平均降水量(mm)。

(2)

式中,Rxj為第j土地利用類型柵格x的Budyko干燥指數,ωx表示植被有效含水量與年均降水量的比值：

(3)

(4)

AWCx=Min(DS,Dr)×PAWCx

(5)

式中,AWCx為柵格x的植被有效含水量(mm)(由有效土壤深度和土壤質地決定)；Z為季節參數(即Zhang系數)用于表征降水的季節性特征,kxj為第j土地利用類型柵格x的植被蒸散系數,ET0x表示柵格x的參考蒸散量(mm);Min為minimum的縮寫,DS為土壤深度(mm),Dr為根系深度(mm),PAWCx為柵格單元x的植被可利用水分(mm)。

2.2 參數來源及分析

2.2.1模型參數

產水量模塊所需輸入參數有氣象、土地利用、土壤深度、植被可利用水、生物物理參數表、流域以及Zhang系數。另外,相關分析涉及NDVI。本研究所有數據的空間分辨率均重采樣到500 m,且坐標系均使用世界大地坐標系統(WGS1984)。

(1)降水量和參考蒸散發

利用中國氣象科學數據中心(http://data.cma.cn/user/toLogin.html)收集周邊氣象站2001—2019年的降水量和氣溫數據。用ANUSPLIN模型進行空間插值,得到降水量和氣溫的空間柵格數據。本文選用Modified—Hargreaves法計算參考蒸散發,與其它方法相比,該方法只需最高氣溫、最低氣溫和太陽大氣頂層輻射便可計算[45]。公式如下：

(6)

式中,Tmax和Tmin分別為月最高氣溫和月最低氣溫(℃),Ra是月太陽大氣頂層輻射(MJ m-2d-1),C、E和T是常數,其值參考胡慶芳等[46]推薦青藏高原C、E和T值。

(2)土地利用

土地利用數據源于美國宇航局地球科學數據系統(https://search.earthdata.nasa.gov/search)中的MCD12Q1數據產品。它是是根據一年的Terra和Aqua觀測所得的數據經過處理,描述土地覆蓋的類型。MCD12Q1采用了五種不同的土地覆蓋分類方案,信息提取的主要技術是監督決策樹分類。本文采用國際地圈—生物圈計劃(IGBP)的分類方法。研究區土地利用類型主要包括草地、永久性濕地、農田、建設用地、農作物—自然植被交錯帶、永久性冰雪和裸地。

(3)土壤數據和PAWC

基于蘭州市寒區旱區科學數據中心下載的第二次土壤調查數據庫生成土壤深度和土壤質地柵格圖。PAWC是指土壤土層中為植物生長提供的水量所占比例,為0—1的小數。周文佐等[47]提出了一種根據物理和化學性質計算PAWC的方法,具體算法如下;

PAWC=54.509-0.132Sand-0.003(Sand)2-0.055Silt-0.006(Silt)2-0.738Clay+0.007(Clay)2-2.668C+

0.501(C)2

(7)

式中,Sand為土壤砂粒含量(%),Silt為土壤粉粒含量(%),Clay為土壤粘粒含量(%)；C為土壤有機質含量(%)。

(4)流域和生物物理參數表

數字高程模型(DEM)由資源環境科學與數據中心(http://www.resdc.cn/)提供的雷達地形測繪SRTM(90 m)數據。基于數字高程模型數據通過ArcGIS 10.4中的水文分析工具生成的子流域。生物物理表主要用來反映土壤覆蓋度和LULC的屬性,包括LULC編碼、植物蒸散系數(Kc)和根系深度等,生物物理表中每種LULC類型的系數來源于文獻[45,48]和《InVEST用戶指南》[49]。

(5)Zhang系數

Zhang系數是反映當地降水格局和水文地質特征的氣候季節因子,其值在1—30之間變化。根據昌馬鋪水文站多年平均徑流量觀測資料,采用水量平衡法對Zhang系數進行校驗[7],當Zhang系數為3.33時,InVEST產水模塊的模擬效果最好。

(6)NDVI與植被根系深度

NDVI可以準確反映地表植被覆蓋狀況,是各種植被指數中最常用的表征地表植被狀況的指標[36]。本研究中所涉及的2001—2019年NDVI數據源于資源環境科學與數據中心(http://www.resdc.cn/)。該數據集是基于SPOT/VEGETATION NDVI衛星遙感數據,在月數據的基礎上采用最大值合成法生成年度NDVI數據集。植被根系深度柵格數據是InVEST模型基于生物物理參數表生成的。

2.2.2時空變化趨勢分析

本研究主要通過OriginPro 9.1中的線性擬合方法分析研究目標(氣溫、降水量和實際蒸散量、NDVI和產水量)的時間變化趨勢。利用Matlab—R2015B對研究目標的空間變化趨勢進行線性回歸分析,該方法廣泛應用于時序數據動態檢測,線性回歸的斜率是量化動態變化趨勢的最佳指標。具體算法如下[50]：

(8)

式中,n為研究年限(本研究n=19),i為年份序號,Mi為第i年的研究目標值,Slope>0表示19年來研究目標變化趨勢為正,Slope<0則表示研究目標變化趨勢為負。

利用線性回歸的相關系數對研究目標動態趨勢進行顯著性檢驗,相關系數的計算公式如下：

(9)

式中,i為研究期間的年份序號,Mi為第i年的研究目標值,cov為協方差函數,var為方差函數。

2.2.3產水功能重要性分級

合理分配、利用和保護水資源的前提是需對區域產水功能重要性進行分級。本研究基于InVEST模型對研究區產水量進行了定量評估,并依據《國家生態保護紅線—生態功能紅線劃定技術指南(試行)》中涵養水源重要生態功能區的分級方法[51],利用ArcGIS 10.4中的Quantile(分位數)功能對2001—2019年的年均產水量進行分類。按產水量值將其重要性由低到高依次劃分為4個重要性級別,即一般重要[0.00—98.38)、中等重要[98.38—227.33)、高度重要[227.33—344.76)和極重要[344.76—488.72)。

2.2.4統計分析

利用MODIS Reprojection Tools與ArcGIS對所獲取的柵格數據進行拼接、提取、裁剪等操作,以得到本研究所需的數據集。其次,將2001—2019年昌馬堡水文站實測地表年徑流量與模擬年總產水量進行線性擬合,并計算調整R2、P以校驗模型結果；擬合直線斜率與1∶1直線斜率間的差異通過SPSS 22.0中的一般線性模型進行分析。另外,利用R v4.0.4對環境要素與產水量之間相關性進行Pearson相關矩陣分析。

3 結果與分析

3.1 模型校驗

將模型預測結果與實測數據進行比較,通過不斷調整參數,得到最適于研究區的預測模型則為模型校驗[14]。本研究利用昌馬堡水文站地表徑流量實測數據對InVEST模型產水量的預測結果進行了校驗。結果顯示(圖2)：模型模擬的多年平均總產水量為13.66×108m3,水文站實測多年平均總徑流量為13.81×108m3,其相對誤差為1.11%。同時,模型模擬的總產水量與實測地表總徑流量之間存在較強的線性關系(R2= 0.986P<0.01)。雖然,模擬總產水量相較實測地表總徑流量略低估了0.37×108m3/a,但兩者所擬直線斜率與1沒有顯著差異(P>0.05)。可見,InVEST模型在本研究區適用性很好。

圖2 InVEST 模擬總產水量與實測地表總徑流量對比 Fig.2 InVEST modelled total water yield vs measured total surface runoff紅線代表斜率為0,截距為1的直線

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

2001—2019年疏勒河上游產水量顯著增加,增加速率為0.26×108m3/a(圖3)。2003年總產水量最低(9.76×108m3),2017年總產水量最高(17.25×108m3)。從空間分布格局來看,山地高海拔區產水量較高,河谷低海拔區產水量較低,尤其是東南部低海拔區與西北部東側最低(圖4)。2001—2019年,所有柵格的產水變化量空間趨勢如圖4所示。產水量顯著增加區域約占整個研究區的32.01%,主要分布在山地高海拔地區；僅有8.39%的區域顯著下降且集中分布在研究區的西北部。

圖3 2001—2019年產水量時間變化Fig.3 Temporal variation of water yield in 2001—2019

圖4 產水量空間分布及變化Fig.4 Spatial distribution and its variation of water yield

3.3 產水功能重要性分級

研究區產水功能重要性分級總體呈現出由河谷低海拔區向山地高海拔區遞增的階梯分布趨勢(圖5)。其中,一般重要區域所占面積最大,比例為47.26%,主要分布在東南部低海拔區與西北部東側；中等重要區域所占比例為24.48%,主要分布在西北部；高度重要和極重要區域所占比例分別為14.51%和13.74%,集中分布于山地高海拔區。

圖5 產水功能重要性空間分級Fig.5 Importance spatial classification of water yield function

3.4 影響流域產水量的因素分析

疏勒河上游的平均海拔高于3945 m,導致該地區氣候寒冷干燥,年均氣溫為-6.06℃,年降水量325.96 mm。2001—2019年研究區年均氣溫、年降水量和年實際蒸散發的時間變化趨勢如圖6所示。年均氣溫與年降水量分別以0.02 ℃/a (P>0.05)和5.92 mm/a(P<0.05)的速率增加,研究區氣候呈現暖濕化趨勢。在此影響下,年實際蒸散量和NDVI分別以3.00 mm/a (P>0.05)和0.003/a (P<0.05)的速率遞增。

圖6 年均氣溫、年降水量、年實際蒸散量和年NDVI的時間變化Fig.6 Temporal variations of the mean annual air temperature, annual precipitation, annual actual evapotranspiration and annual NDVI

研究區氣溫、降水量、實際蒸散量和NDVI的空間分布均具有明顯的空間異質性。具體而言,氣溫的空間變化范圍為-16.92—5.33℃,自西北向東南逐漸下降且河谷低海拔區高于山地高海拔區,2001—2019年研究區大部分區域(76.08%)的氣溫呈明顯上升趨勢,但未達到統計學顯著水平(P>0.05)。空間上,降水量在127.42—527.13 mm之間變化,總體呈西北低、東南高且河谷低海拔區低于山地高海拔區的分布格局(圖7),與氣溫的分布情況恰好相反。研究期間,降水量整體顯著增加(P<0.05)且東南局部地區(24.41%)極顯著增加(P<0.01)。實際蒸散量的空間變化范圍為0.00—357.54 mm, 其高值區主要位于東南部低海拔區與西北部東側；2001—2019年實際蒸散量顯著增加區域約為33.39%(P<0.05),主要分布河谷低海拔區,而顯著下降區域約為4.01%(P<0.05),散落分布在整個研究區。NDVI空間分布格局與實際蒸散量的分布情況一致,研究期間,NDVI顯著增加區域約為56.92%(P<0.05),分布在河谷低海拔局部地區,極顯著下降區域約占11.27%(P<0.05),主要分布在山地高海拔區。此外,植被根系深度的空間分布與NDVI一致。

圖7 環境要素空間分布及其變化Fig.7 Spatial distributions and the variations of environmental elements

Pearson相關分析表明(圖8),產水量與降水量顯著正相關,卻與氣溫、實際蒸散量、NDVI和根系深度顯著負相關(P<0.001)。不同環境因素之間的相關性表明,實際蒸散量與根系深度、NDVI和氣溫顯著正相關(P<0.001)。氣溫升高或NDVI、根系深度增加均可能會促進實際蒸散量增加,從而對產水量產生負向影響。

圖8 產水量與影響因子的相關矩陣Fig.8 Correlation matrix of water yield and driving factors***, P<0.001

4 討論

產水量是一項重要的生態功能指標,直接影響著社會經濟和生態系統的可持續發展[34]。本研究著重對疏勒河上游產水量的時空變化特征以及影響因素進行了探討分析。結果表明,2001—2019年疏勒河上游產水量顯著上升,這同呂樂婷等[24]對1980—2016年三江源區的研究結果類似,卻同王雅等[25]對2000—2014年祁連山黑河中游的研究結果相反。這可能與各研究區的氣候條件以及研究的時間尺度不同有關。產水量是區域水循環過程中綜合考慮收(降水量)支(實際蒸散量)平衡的產物[10],因此,區域降水量和實際蒸散量以及兩者之間平衡關系是影響產水量和產水能力的主導要素,這是許多產水量研究中的一致結論[7,52—53]。潘韜等[7]和呂樂婷等[24]對三江源產水量的研究得出,降水和實際蒸散發是決定生態系統產水量的兩個關鍵環節,其中降水量與產水量顯著正相關,實際蒸散量與其顯著負相關。劉軍會等[52]和劉婭等[53]分別對青藏高原和京津冀地區的產水量進行了研究,結果表明降水量和實際蒸散量綜合影響區域產水變化量。然而,氣溫、NDVI和植被根系深度會通過影響實際蒸散量來間接影響產水量。其中,NDVI和植被根系深度主要通過下墊面植被變化來影響實際蒸散量,進而影響產水量。Lu等[35]對青藏高原拉薩河流域產水量的研究發現,降水量和NDVI變化是影響區域產水量的關鍵要素,特別是通過植被變化來推動產水量變化。19年來,疏勒河上游實際蒸散增加趨勢(3.00 mm/a)低于降水增加趨勢(5.92 mm/a),降水量對產水量的補給遠超過實際蒸散對產水量的消耗,因此,疏勒河上游產水量整體呈顯著增加趨勢。

疏勒河上游發育大面積多年凍土,其中山地高海拔區的土層存在凍結面,這使得降水無法下滲而在凍結面匯流形成徑流。與此同時,山地高海拔區的降水量高,氣溫和實際蒸散量低且土地利用類型主要為冰川和裸地。冰川屬于水域,最易融化形成徑流；裸地的植被覆蓋度低且根系淺,截留能力弱,更多降水也直接形成徑流[24,54—55]。相比之下,河谷低海拔區由于凍土分布面積和深度受限,降水或冰水可以滲入土壤深層,徑流形成的幾率大大降低。另外,河谷低海拔區的主要土地利用類型為草地,其植被覆蓋度和實際蒸散量都很高,并且植被根系深而密集,可以有效的截留降水,使得產水量偏低。2001—2019年山地高海拔區產水量顯著增加,主要是因為降水量補給顯著增加而蒸散量無顯著變化,致使生態產水以增加為主。與其相反,西北地區的裸地被草地和農田大面積取代,相比裸地,草地與農田的根系較深且植被生長消耗的水量也比較大,因此產水量顯著下降。

疏勒河上游產水重要性分級與產水量空間分布格局基本一致,這與孫艷偉等[14]和王玉純等[51]研究結論一致。具體而言,產水極重要和高度重要區降水量較多、蒸散量較小,并且多分布冰川和裸地,是疏勒河流域水資源集成區,而產水一般重要和中等重要地區為相對干旱的堆積平原和盆地,為水資源消耗區。為使研究區實現健康可持續發展,極重要和高度重要區需作為重點和優先保護的水資源地而納入政策和管理決策之中[56]。

本研究對疏勒河上游產水量的模擬結果較為理想,基本符合區域的實際情況。但因模型自身設定、參數算法以及數據精度等原因,使研究仍存在一些不確定性。首先,InVEST模型簡化了匯流過程,沒有區分地表徑流、土壤中徑流和基流,并假定所有柵格的產水量都到達同一出水口[14],并且產水量常與社會經濟和人類活動密切相關,而模型的輸入數據只有自然數據。特別是,疏勒河上游是一個典型的高寒山區,InVEST模型未考慮冰川和凍土融水,這會增加水文模擬的不確定性。因此,InVEST模型在這些方面都有待進一步改進。其次,降水量和參考蒸散量是精確模擬產水量的重要參數,而西北地區水文站點較少且氣候變化具有復雜性,利用ANUSPLIN模型空間插值生成的降水數據以及通過Modified—Hargreaves法計算的參考蒸散發的準確性都會受到氣象數據精度的影響[57],從而影響模型的模擬精度。盡管存在這些不確定性,但研究結果仍可以有效反映產水量的時間變化和空間分布,揭示產水量與環境要素之間的關系,為疏勒河流域水資源的有效管理、合理利用和保護提供科學參考。

5 結論

基于InVEST模型,結合環境要素對疏勒河上游產水量的時空分布格局及影響因素進行分析,得出：

(1) 模型中的Zhang系數為3.3時,研究區產水量模擬效果最佳；研究區年均總產水量為13.66×108m3,2001—2019年研究區產水量以0.26×108m3/a速率顯著增加。從空間分布來看,山地高海拔區產水量較高,河谷低海拔區產水量較低,尤其東南部低海拔區和西北部東側最低。研究期間,產水量顯著增加區約占總面積的32.01%,主要分布在山地高海拔區；僅有8.39%的研究區顯著下降且集中分布于西北部。

(2) 產水量與降水量顯著正相關,而與氣溫、實際蒸散量、NDVI和根系深度顯著負相關；實際蒸散量與氣溫、NDVI和植被根系深度顯著正相關。

(3) 在對研究區產水量進行定量評估的基礎上,進行了產水功能重要性空間分級。得出：一般重要區域所占面積最大,約占整個研究區的47.26%,主要位于東南部低海拔區與西北部東側；高度重要和極重要區域共占總面積的28.25%,集中分布于山地高海拔區。