江 頌，蒙吉軍?，陳奕云

(1.北京大學城市與環(huán)境學院，地表過程分析與模擬教育部重點實驗室，100871，北京；2.武漢大學資源與環(huán)境科學學院，430079，武漢)

隨著社會經(jīng)濟的發(fā)展，水資源已成為全球多地區(qū)發(fā)展的限制因子之一[1]。作為生態(tài)系統(tǒng)的重要服務功能，區(qū)域水資源的產(chǎn)生與供給對確保灌溉農(nóng)業(yè)、工業(yè)及旅游業(yè)的發(fā)展，促進居民生活水平的提高具有重要意義[2]；因此，流域生態(tài)系統(tǒng)產(chǎn)水服務及其影響機制的研究逐漸成為水文學、流域規(guī)劃與管理等領域的研究熱點[3]。已有研究表明：在長時間尺度上，氣候是降水、產(chǎn)水等水資源變化的主導因素；但在較短時間尺度上，土地利用/覆被變化(land use and land cover change，LUCC)才是水文效應的重要影響因素[4]。揭示和預測LUCC的水文效應，對指導區(qū)域土地利用規(guī)劃，優(yōu)化區(qū)域水土資源配置，以及流域的可持續(xù)發(fā)展至關重要[5]。近年來，關于土地利用與水文過程關系的研究正逐漸增加，特征參數(shù)[6]、衛(wèi)星遙感[7]、水文模型[8]等方法技術被綜合應用，植被變化、農(nóng)田增減、城鎮(zhèn)化等用地變化被認為可能影響水文過程[9-11]。目前，大部分研究都關注于土地利用比例對水文過程的影響，而對景觀格局的水文效應研究較少[12]；而且，在探討土地利用水文效應時，較多研究采用相關分析和多元線性回歸分析的方法，沒有考慮到空間因子，忽略了空間依賴性和空間非平穩(wěn)性對水文效應的影響。

自20世紀60年代起，隨著計算機技術的發(fā)展，水文效應的研究方法由傳統(tǒng)的統(tǒng)計分析向水文模型轉(zhuǎn)變。InVEST模型便是目前應用廣泛的一種。InVEST(the integrated valuation of ecosystem services and tradeoffs tool)模型，全稱“生態(tài)系統(tǒng)服務功能綜合估價和權(quán)衡得失評估模型”，由美國斯坦福大學、大自然保護協(xié)會(TNC)和世界自然基金會(WWF)聯(lián)合開發(fā)而成[8]。使用該模型能夠快捷地對包括產(chǎn)水功能在內(nèi)的生態(tài)系統(tǒng)多種服務功能進行量化，從而很好地為政府及相關部門提供決策依據(jù)。

基于此，本文選擇水資源供需矛盾突出的黑河中游為研究區(qū)，建立InVEST模型計算流域產(chǎn)水：一方面，將空間因子考慮在內(nèi)，通過普通最小二乘回歸、空間回歸以及地理加權(quán)回歸等多種方法，分析土地利用類型與流域產(chǎn)水的關系；另一方面，結(jié)合主成分分析篩選景觀指數(shù)，進一步揭示景觀格局與區(qū)域產(chǎn)水的關系，從而為區(qū)域土地利用規(guī)劃與決策提供參考，促進水土資源的持續(xù)利用。

1 研究區(qū)概況

黑河是中國第2大內(nèi)陸河，發(fā)源于祁連山中段，流經(jīng)青海、甘肅、內(nèi)蒙古3省自治區(qū)[13]。黑河鶯落峽至正義峽之間為中游，東靠武威和金昌，西至嘉峪關和酒泉，南與青海省接壤，北和內(nèi)蒙古毗鄰，具體范圍介于E 97°20′～102°12′，N 37°28′～39°57′。中游干流全長185 km，總面積約1.96萬km2。黑河中游地處青藏高原向內(nèi)蒙古高原的過渡地帶，海拔在1 234～4 886 m之間，地勢東南高西北低，其北部山地和走廊平原屬溫帶大陸性干旱氣候，南部祁連山區(qū)為高寒半干旱氣候。年降水量僅54.9～436.2 mm，且年內(nèi)分布不均，多集中在6—9月；但全年蒸發(fā)量高達1 700 mm[14]。黑河中游是典型的灌溉農(nóng)業(yè)綠洲，也是整個流域人類活動最為密集、綠洲最為集中、經(jīng)濟最為發(fā)達的地區(qū)[15]。流域人均可利用水資源量為1 250 m3(僅為全國平均水平的57%)，是典型的資源型缺水地區(qū)。

本研究所需數(shù)據(jù)包括基礎地理信息數(shù)據(jù)、土地利用數(shù)據(jù)、氣象數(shù)據(jù)和土壤數(shù)據(jù)等。基礎地理信息數(shù)據(jù)來自中國科學院寒區(qū)旱區(qū)科學數(shù)據(jù)中心的黑河計劃數(shù)據(jù)管理中心(www.heihedata.org)。土地利用數(shù)據(jù)包括黑河中游2000、2010和2014年土地利用柵格數(shù)據(jù)。其中，2000年土地利用數(shù)據(jù)來源于中國科學院“八五”重大應用項目[16]，2010、2014年土地利用數(shù)據(jù)是通過黑河中游Landsat影像進行人機交互式目視解譯得到，解譯精度均>88%。研究中，將土地利用類型劃分為耕地、草地、林地、水域、建筑用地和未利用地6類，數(shù)據(jù)的空間分辨率為30 m。

氣象數(shù)據(jù)來自中國氣象數(shù)據(jù)網(wǎng)(cdc.nmic.cn/home.do)。對研究區(qū)及其周邊的11個氣象站點求得多年平均降水量，再通過普通克里金插值得到黑河中游多年平均降水量。使用“改進的哈格里夫斯”(modified Hargreaves)公式計算并插值得到多年平均參考蒸散量[17]。土壤深度、植物可利用含水量數(shù)據(jù)均來自世界土壤數(shù)據(jù)庫(HWSD)，空間分辨率為1 km。

DEM數(shù)據(jù)來自美國馬里蘭大學全球土地覆被數(shù)據(jù)庫，分辨率為30 m。基于DEM數(shù)據(jù)，在ArcGIS中使用Hydrology工具將研究區(qū)劃分為20個子流域(圖1)。

圖1 黑河中游子流域劃分結(jié)果Fig.1 Divided subwatersheds in the middle reaches of Heihe River

2 研究方法

2.1 基于InVEST模型的產(chǎn)水量模擬

InVEST模型產(chǎn)水模塊的計算原理，主要是通過降水量減去蒸發(fā)蒸騰水量來得到產(chǎn)水量。該模型忽略地表水與地下水的交互作用，不區(qū)分地表徑流、地下徑流以及基流。

產(chǎn)水模塊是基于Budyko曲線和年平均降水量進行計算的。年產(chǎn)水量的計算公式為

(1)

式中：Y(x)為x類景觀像素的年產(chǎn)水量，mm；AET(x)為年實際蒸散發(fā)量，mm；P(x)為年降水量，mm。

2.2 回歸分析模型

1)普通最小二乘回歸(ordinary least squares，簡稱OLS)是最簡單、最基本且應用最為廣泛的模型之一。模型的一般形式為

(2)

式中：Y為因變量，mm；Xi為第i個自變量，mm；βi為回歸系數(shù)；ε為隨機誤差項，mm。

2)空間誤差模型(spatial error model，簡稱SEM)通過空間誤差項來引入空間依賴性。其表達式為

Y=Xβ+ε,

(3)

ε=λWε+μ。

(4)

式中：λ為自回歸參數(shù)，衡量樣本觀測值中的空間依賴作用；μ為空間殘差項，mm；W為空間權(quán)重矩陣。

3)地理加權(quán)回歸(geographical weighted regression，簡稱GWR)是一種對空間變化關系建模的線性回歸的局部形式，是對普通回歸模型的擴展。該模型假定了其擬合系數(shù)是變量空間位置的函數(shù)，將變量的地理位置引入到回歸參數(shù)中，再利用局部加權(quán)最小二乘法進行逐點的參數(shù)估計[18]。其中，權(quán)是回歸點所在空間位置到其他各采樣點空間位置的距離函數(shù)。其表達式一般為

(5)

式中：Yi為第i個采樣點的因變量，即待估測值，mm；Xik為第i個采樣點上的第k個解釋變量，mm；(ui,vi)為第i個采樣點的坐標；βk(ui,vi)為第i個采樣點上的第k個回歸參數(shù)，是地理位置的函數(shù)。

在回歸分析中，分別以土地利用類型和景觀指數(shù)指標作為自變量。其中，土地利用類型指標以每個子流域的各類用地面積占該子流域總面積的比值來表示，計算公式如下：

(6)

式中：Xij為第i個子流域中第j類用地面積指標；Sij為第i個子流域內(nèi)第j類用地的面積，m2；Si為第i個子流域的總面積，m2。

3 結(jié)果分析

3.1 黑河中游產(chǎn)水量時空格局

2000、2010和2014年的產(chǎn)水量分別為1.73億、1.63億和1.38億m3，呈現(xiàn)出逐年減少趨勢，尤其是2010年至2014年，產(chǎn)水量減少幅度更大。

圖2中，虛線是整個黑河中游當年的平均產(chǎn)水深度，對比可知子流域間產(chǎn)水能力差異極大。三期中，產(chǎn)水集中于15～20號這6個子流域，15、16、18和19號子流域的產(chǎn)水深度遠高于平均產(chǎn)水深度，而1～14號子流域產(chǎn)水深度則遠低于平均產(chǎn)水深度，部分子流域甚至為0。

圖2 黑河中游2000(a)、2010(b)和2014(c)年各子流域產(chǎn)水深度Fig.2 Water yield depths of subwatersheds in the middle reaches of Heihe River in 2000 (a), 2010 (b), and 2014 (c)

由圖3可見，三期產(chǎn)水量均呈現(xiàn)出東南多、西北少的空間分布。對比產(chǎn)水量與降水量分布格局，發(fā)現(xiàn)盡管產(chǎn)水量和降水量都有著東南多、西北少的大格局特征，但在小格局上的分布規(guī)律并不完全一致。說明產(chǎn)水量受到降水量的控制，但子流域間的局部差異還受到下墊面的影響。

圖3 黑河中游2000(a)、2010(b)和2014(c)年各子流域產(chǎn)水深度空間分布Fig.3 Spatial distribution of water yield depths of subwatersheds in the middle reaches of Heihe River in 2000 (a), 2010 (b), and 2014 (c)

3.2 黑河中游土地利用的水文效應

3.2.1 土地利用水文效應的普通最小二乘回歸 在土地利用水文效應的回歸建模中，6類指標間具有線性關系，若同時將其作為自變量，會造成變量冗余；因此，先去掉草地(16.29%)、未利用地(57.38%)這2類用地(兩者在研究區(qū)所占面積最大，會對其它類型用地面積指標值影響較大)，只將耕地、林地、水體、建設用地4類用地指標作為自變量，將極差標準化后的產(chǎn)水指標作為因變量，進行普通最小二乘回歸(表1)。

表1 黑河中游4類用地對產(chǎn)水深度的普通最小二乘回歸結(jié)果

對于P>0.05的自變量，其系數(shù)不具有顯著性，尚不能認為有統(tǒng)計學意義，在此不予討論[19]。由表1可知，在2000年的回歸模型中，林地和水體的系數(shù)為正，表明林地和水體促進產(chǎn)水。在InVEST產(chǎn)水模型中，產(chǎn)水量與產(chǎn)流量的意義有所不同，模型計算的產(chǎn)水量包含在之后可能轉(zhuǎn)化為的土壤含水量、凋落物持水量以及冠層截留量；因此，林地因其有樹蔭遮蔽，整體溫度較低，蒸發(fā)速度慢，且其葉片、根系等具有持水功能，整體上表現(xiàn)為促進產(chǎn)水。而水體的蒸發(fā)形式是水面蒸發(fā)，主要受控于氣象條件，由于水體比熱容較大，其溫度會低于陸面，使蒸發(fā)量往往更小，從而促進產(chǎn)水[20]。2010年，耕地的回歸系數(shù)為正，表明耕地也有利于產(chǎn)水。建設用地的系數(shù)為負，表明建設用地對產(chǎn)水有抑制作用，這可能由于城鎮(zhèn)的建筑用地具有熱島效應，致使水分的蒸散發(fā)量較大，從而不利于產(chǎn)水。2014年的回歸結(jié)果與前2期相似，耕地和林地促進產(chǎn)水，建設用地抑制產(chǎn)水。

同樣地，將未利用地和草地這2類用地指標作為自變量，將極差標準化后的產(chǎn)水指標作為因變量，進行普通最小二乘回歸(表2)。

由表2可見，草地的系數(shù)均不顯著，2000、2010和2014年中未利用地的系數(shù)均為負，表明未利用地抑制產(chǎn)水。黑河中游的未利用地主要分布在北部地區(qū)，其地面多表現(xiàn)為荒漠，土壤類型主要為黏重鹽沙，蓄水能力差，蒸散量大，故不容易產(chǎn)水。

3.2.2 土地利用水文效應的空間回歸 空間依賴性檢驗能對模型選擇做出判斷。通過查看各普通最小二乘回歸模型的空間依賴性檢驗結(jié)果，發(fā)現(xiàn)2010、2014年的回歸模型中拉格朗日乘數(shù)均不顯著(P>0.05)，表明對其保留OLS回歸結(jié)果即可，而以2000年草地和未利用地為自變量的回歸模型中，LM-lag不顯著(P=0.28)而LM-error顯著(P=0.03)，表明有必要使用空間誤差回歸模型進行改進。

于是，以2000年草地和未利用地指標為自變量，以產(chǎn)水深度為因變量，建立空間誤差模型(表3)。

由表3可知，空間回歸與普通最小二乘回歸結(jié)果相同，未利用地系數(shù)為負，抑制產(chǎn)水。且空間誤差模型的R2為0.71，而普通最小二乘回歸的R2為0.59，可見模型性能得到的了較大提升。為檢驗空間回歸模型效果，分別計算OLS和SEM的回歸殘差的全局Moran’s I，從而檢驗回歸模型的殘差是否具有空間隨機性。經(jīng)過計算得到，OLS殘差的Moran’s I為-0.43(P=0.03)，表明OLS的殘差具空間負相關性。而SEM殘差的Moran’s I為0.12(P=0.11)，與OLS相比，指數(shù)的絕對值明顯減小，結(jié)合指數(shù)對應的P值，可以認為空間誤差模型消除了殘差的空間負相關性。

表2 黑河中游2類用地對產(chǎn)水深度的普通最小二乘回歸結(jié)果

表3 2000年黑河中游2類用地對產(chǎn)水深度的空間誤差回歸模型

注：λ：空間誤差系數(shù)；z值：z檢驗統(tǒng)計量。Notes:λ: Spatial error coefficient.z-value:ztest statistics.

3.2.3 土地利用水文效應的地理加權(quán)回歸 根據(jù)地理學第一定律，OLS中的平穩(wěn)性假設并非恒成立，而地理加權(quán)回歸能夠容許一些不平穩(wěn)的數(shù)據(jù)被模擬。經(jīng)計算，2000年黑河中游產(chǎn)水深度的Moran’s I為0.32(P=0.02)，表明因變量確實具有空間正相關性，適合采用地理加權(quán)回歸模型進行擬合。

以2000年草地和未利用地指標為自變量，以產(chǎn)水深度為因變量，建立地理加權(quán)回歸模型。得到的R2為0.90，遠高于OLS的R2，證明地理加權(quán)回歸明顯地提升模型性能。地理加權(quán)回歸模型對局部進行了回歸，因此每個子流域都有1個回歸系數(shù)，未利用地在各子流域的回歸系數(shù)如圖4所示。

圖4 黑河中游2000年未利用地對產(chǎn)水深度的地理加權(quán)回歸系數(shù)Fig.4 Geographical weighted regression coefficient of bare land on water yield depth in the middle reaches of Heihe River in 2000

由圖4可見：與普通最小二乘回歸和空間回歸結(jié)果相一致，未利用地的地理加權(quán)回歸系數(shù)整體上為負，但在各子流域中存在差異。在最南端，回歸系數(shù)為正，表明未利用地在此處促進產(chǎn)水；在北部地區(qū)，系數(shù)為負，表明抑制產(chǎn)水；而中偏南部地區(qū)，回歸系數(shù)達到最小值-1.17，表明未利用地在此處對產(chǎn)水的抑制作用最強。這種分布可能與黑河中游土壤質(zhì)地及局部地形有關，南部土壤松軟，容易蓄水持水，蒸發(fā)量相對減少，因此容易產(chǎn)水；由南向北土壤變得黏重，蓄水能力差，再加上中部主要為平原地形，接受太陽輻射多，蒸發(fā)量相對更大，因此不利于產(chǎn)水；而北部地區(qū)雖然土壤也屬于黏重型，但因山地地形的存在，使當?shù)赜猩降仃幱?、坡地小氣候等對降水、蒸發(fā)產(chǎn)生影響，從而使北部未利用地對產(chǎn)水的抑制作用不及中部偏南地區(qū)強烈。

結(jié)合研究區(qū)2000—2014年土地利用變化趨勢分析，總體來看，黑河中游產(chǎn)水能力的下降受到了未利用地面積增加以及林地、水體、耕地面積減少的影響。

3.3 黑河中游景觀格局的水文效應

3.3.1 景觀格局指數(shù)主成分分析 參考已有研究[21-23]，選擇目前常用、且具有代表性的景觀指數(shù)共9個，包括最大斑塊占景觀面積比例(LPI，largest patch index)、斑塊面積算術平均(AREA_MN，mean of patch area)、斑塊面積標準差(AREA_SD，standard deviation of patch area)、形狀指數(shù)加權(quán)平均(SHAPE_AM，area-weighted mean of shape index)、分形維數(shù)加權(quán)平均(FRAC_AM，area-weighted mean of fractal index)、斑塊多度(PR，patch richness)、斑塊多度密度(PRD，patch richness density)、修正Simpson多樣性指數(shù)(MSIDI，modified Simpson’s diversity index)、修正Simpson均勻度指數(shù)(MSIEI，modified Simpson’s evenness index)?；?000、2010及2014年土地利用數(shù)據(jù)，使用Fragstats 4.2[24]計算出黑河中游的9種景觀指數(shù)。

對黑河中游3期景觀指數(shù)進行主成分分析，結(jié)果如表4所示。

表4 黑河中游2000、2010和2014年景觀指數(shù)主成分分析結(jié)果

注：加粗數(shù)據(jù)表示該景觀指數(shù)在該主成分中載荷較大。Notes: The bold data indicates that the landscape index has a greater load in the principal component. LPI: Largest patch index. AREA_MN: Mean of patch area. AREA_SD: Standard deviation of patch area. SHAPE_AM: Area-weighted mean of shape index. FRAC_AM: Area-weighted mean of fractal index. PR: Patch richness. PRD: Patch richness density. MSIDI: Modified Simpson’s diversity index. MSIEI: Modified Simpson’s evenness index.

由表4可見，3期數(shù)據(jù)的前3個特征根的累計貢獻率均大于85%，即3個年份的景觀格局均可使用前3個主成分描述。對比不同年份的主成分載荷可以發(fā)現(xiàn)，第1主成分主要由LPI、FRAC_AM、MSIDI、MSIEI來解釋，主要代表景觀豐富度、均勻度；第2主成分主要由AREA_SD解釋，代表斑塊面積特征；第3主成分主要由PRD解釋，代表景觀斑塊相互間影響情況[23,25]。

3.3.2 景觀格局水文效應的普通最小二乘回歸 基于主成分分析，將極差標準化后的主成分得分作為自變量，將產(chǎn)水深度作為因變量，進行普通最小二乘回歸(表5)。

由表5可見，3期數(shù)據(jù)中第1主成分的系數(shù)均顯著，且均為正，表明景觀豐富度、均勻度與產(chǎn)水成正相關。2000年中的第2主成分系數(shù)也顯著，其系數(shù)為負，第2主成分主要由斑塊面積標準差AREA_SD解釋，表明斑塊間面積差距越大即景觀均勻程度越低，越不利于產(chǎn)水。綜合推斷，在黑河中游，景觀類型越豐富多樣，景觀斑塊分布越均勻，則越有利于產(chǎn)水。

表5 景觀指數(shù)主成分對產(chǎn)水深度的普通最小二乘回歸結(jié)果

注：pc1：第1主成分；pc2：第2主成分；pc3：第3主成分。Notes: pc1: First principal component. pc2: Second principal component. pc3: Third principal component.

總體來看，2000—2014年期間研究區(qū)景觀多樣性和景觀均勻度的下降，也是黑河中游產(chǎn)水能力下降的原因之一。

4 結(jié)論與討論

1)黑河中游產(chǎn)水能力在2000—2014年期間呈現(xiàn)出減小趨勢，與未利用地面積增加和林地面積減少、以及景觀多樣性和均勻度下降等景觀格局變化有直接關系。2000、2010和2014年產(chǎn)水量分別為1.73億、1.63億和1.38億m3。產(chǎn)水量空間分布不均勻，東南多、西北少，子流域間產(chǎn)水能力差異大，這與王蓓等[26]的研究結(jié)論基本一致。

2)黑河中游水文效應的尺度特征明顯。在土地利用尺度上，林地、水體、耕地等有助于促進產(chǎn)水，而建設用地和未利用地則抑制產(chǎn)水；在景觀尺度上，景觀類型越豐富、斑塊分布越均勻，越有利于產(chǎn)水，這與前人的研究結(jié)論相近[20,27]。因此，結(jié)合黑河中游生態(tài)安全建設[28]，對于北部地區(qū)，加強未利用地的封育保護，增加荒漠植被群落的穩(wěn)定性，防止沙漠化的進一步擴展；對于水土條件較為優(yōu)越的中部偏南地區(qū)，在維持地類多樣性的同時，控制農(nóng)業(yè)綠洲無序擴展，協(xié)調(diào)各類用地比例，提升景觀斑塊的均勻度，合理配置生態(tài)用水；對于產(chǎn)水量較高的南部地區(qū)，加快林草植被的保護與恢復，保障流域源頭的森林生態(tài)屏障。

3)普通最小二乘回歸、空間回歸以及地理加權(quán)回歸結(jié)果均能表明黑河中游未利用地整體上抑制產(chǎn)水。但與普通最小二乘回歸模型相比，空間回歸模型消除了回歸殘差的空間負相關性，使模型性能得到提升。地理加權(quán)回歸使模型性能得到進一步提升，并揭示出未利用地回歸系數(shù)的空間差異。

4)筆者從土地利用和景觀格局2個尺度揭示黑河中游水文效應的影響因素與影響機制，為區(qū)域生態(tài)系統(tǒng)管理提供參考。本研究關注揭示土地利用活動對水文效應的影響，沒有考慮氣候變化對產(chǎn)水的影響。今后在更長時間序列的研究中，可結(jié)合氣候變化來分析對流域產(chǎn)水的綜合影響。另外，研究中所選的景觀格局指數(shù)雖依據(jù)主成分分析和已有研究，但仍難避免主觀因素的影響。未來研究可結(jié)合研究區(qū)實際情況改進景觀格局指數(shù)篩選方法，深入探討景觀格局的水文效應。