岷江上游土地利用與生態系統服務價值的動態變化

左玲麗,彭文甫,*,陶 帥,祝 聰,徐新良

1 四川師范大學地理與資源科學學院, 成都 610068 2 四川師范大學西南土地資源評價與監測教育部重點實驗室, 成都 610068 3 新疆大學資源與環境科學學院, 烏魯木齊 830046 4 中國科學院資源環境科學數據中心, 北京 100101

生態系統服務是指生物體從維持地球生命支持的生態系統中獲得的利益[1],其可持續供給是經濟社會可持續發展的基礎[2]。生態系統服務價值(Ecosystem Services Values，ESV)作為評價區域生態系統服務強弱的一種方法[3],開展ESV的科學評估并對其時空演化特征進行定量描述,對區域制定合理的生態保護政策具有重要的意義[4]。土地利用/覆被變化是人類活動與自然環境相互作用最直接的表現形式[5],是引起全球環境變化的重要因素[6]。土地利用變化影響生態系統格局與過程,改變著生態系統提供產品與服務的能力[7],從而影響著ESV的變化。近年來,隨著3S技術的發展及運用,土地利用變化與生態系統服務價值的影響研究成為熱點。

自Costanza等[8]實現了全球ESV的量化后,掀起了ESV估算熱潮[9]。國內,謝高地等[10- 11]基于Costanza等[8]研究建立了中國生態系統單位面積服務價值當量表,并被國內學者廣泛使用。由于生態系統存在異質性和復雜性,直接引用當量因子法或者只考慮ESV的靜態估算價值,將難以完全解釋區域ESV的動態變化[12]。因此,修正當量表與研究區當量表之間的差異成為了研究熱點。目前,有學者在對謝高地等制定的當量因子進行修正的基礎上開展了區域ESV的研究。如姜棟棟等[13]利用植被凈初級生產力(Net Primary Productivity,NPP)進行區域修正；王永琪等[14]利用謝高地等制作的生物量因子表進行修正；劉倩等[15]采用區域差異系數、社會發展系數進行調整；雷軍成等[16]利用植被覆蓋度指數進行修正；謝高地等[17]、潘洪義等[12]利用凈初級生產力、降水量、土壤保持量數據計算其調整系數。但在修正過程中,同時考慮自然和社會經濟因素對ESV的影響仍然關注較少。隨著ESV研究的不斷深入,針對不同的空間尺度(行政區域、流域等)、生態系統(森林、河流、草地等)以及功能都有了大量研究[18],但鮮有專門針對山地等典型地區的研究,尤其是我國西南生態環境脆弱山區ESV的時空分布特征及演變趨勢仍難以完全被表征。

岷江上游地處我國地形階梯中第一級階梯與第二級階梯過渡地帶[19],是重要的生態屏障區、生物多樣性區和典型生態環境脆弱區。除此之外,岷江上游也是集資源供給、生態服務、環境調節等多種功能與經濟社會發展于一體的復雜系統,對整個岷江流域乃至長江上游地區的生態平衡具有重要意義[20]。鑒于上述情況,利用空間異質系數、社會經濟調整系數、資源稀缺系數對當量進行修正,構建了更能反映研究區實際情況的ESV模型。且基于遙感解譯的2000—2015年土地利用數據,利用CA-Markov模型預測了2035年岷江上游土地利用格局,進一步分析岷江上游土地利用與生態系統服務價值的時空演變趨勢,為區域可持續發展提供科學依據。

1 研究區及數據來源

1.1 研究區概況

岷江上游介于 102°32′—104°15′E,30°45′—33°09′N 之間,位于四川西北部,地處我國青藏高原東麓、橫斷山脈北端與川西北高山峽谷東側的結合部(圖1)[20]。地形結構復雜,以高山峽谷地貌為主,河谷深切,山高坡陡。區域內多地震、滑坡、泥石流等自然災害,生態環境極其脆弱。該區域是長江上游生態保護屏障的重要組成部分,也是成都平原生態保護的生命線[20]。岷江上游流經的行政區域大部分在阿壩藏族羌族自治州,主要包括茂縣、松潘、汶川、理縣和黑水五縣,面積 24741.05 km2。根據阿壩州2016年年鑒,2015年岷江上游總人口約為39.01×104人,其中農業人口約28.77×104人,非農業人口約10.24×104人。區域內地廣人稀,羌、藏、回等少數民族聚居,是民族遷移的“走廊”和交匯地帶[21]。地區經濟發展整體落后,是中國傳統農耕文化區向游牧文化區的交匯面[21]。

圖1 岷江上游地理位置Fig.1 The geographical position of the upper reaches of the Minjiang River

1.2 數據來源與處理

岷江上游2000、2005、2010和2015年的四期遙感影像來源于地理空間數據云網站(http://www.gscloud.cn/)。選取了云量最少的5—9月份的Landsat5/7和Landsat 8遙感影像數據為數據源,軌道號為130/37、130/38和130/39。參考中國科學院的土地利用和覆被變化分類方法,結合岷江上游實際情況,分為耕地、林地、草地、水域、建設用地和未利用地。利用ENVI軟件進行預處理,具體步驟包括：幾何糾正、圖像拼接和裁剪等。再通過目視解譯,得到四期岷江上游土地利用分類圖,經過隨機抽樣檢查Kappa檢驗精度達到96%以上。

生態系統服務價值估算所用到的耕地面積、糧食作物產量、CPI和修正系數計算所需要的社會經濟數據來源于全國和地方統計年鑒以及國民經濟和社會發展公報(http://www.stats.gov.cn/tjsj/)。2015年糧食平均收購價來源于《中國糧食年鑒》,NPP計算所需要的全國和地方的降水和氣溫資料來源于中國國家氣象信息中心和中國環境狀況公報以及水資源公報。驅動因子數據來源于資源環境科學與數據中心(http://www.resdc.cn/)、地理空間數據云(http://www.gscloud.cn/)和OSM(https://www.openstreetmap.org)。

2 研究方法

2.1 土地利用變化特征分析

土地利用變化幅度指各土地利用類型的面積在數量上的變化,有助于了解研究區土地利用面積的變化趨勢[22]。計算公式如下：

(1)

式中,J為研究期某種土地利用類型的年變化量(hm2)；Ua和Ub分別為研究期初和研究期末某種土地利用類型的面積(hm2)；T為年數。

土地利用轉移圖譜,利用ArcGIS空間分析模塊中地圖代數疊加運算,可以反映不同時期土地利用類型轉換情況。計算公式如下：

Y=10×P+Q

(2)

式中,Y為土地利用圖譜單元類型；P為前期的土地利用類型；Q為后期的土地利用類型。從而生成兩個年份間的土地利用轉移圖譜。

2.2 生態系統服務價值評估

依據Costanza等[8]的研究結果,生態系統服務可分為9種[10]。本文參考謝高地等2007年修訂的“中國陸地生態系統單位面積生態系統服務價值當量表”[10]。由于謝高地等人的研究是基于全國生態系統的平均水平,未能有效反映生態系統的空間異質性、生態系統資源和服務功能的區域稀缺性與需求性[23]、區域社會經濟發展水平的不均衡性。因此,本文參考已有研究[24- 25],并結合研究區域的實際情況,通過空間異質系數、資源稀缺系數和社會經濟調整系數進行修正。

2.2.1單位面積生態系統服務的經濟價值確定

參考謝高地等[26]的研究,單位面積農田食物生產生態系統服務價值(1個標準當量因子)相當于全國平均糧食單產市場價值的1/7。計算公式如下:

(3)

式中,E為單位面積農田生態系統提供食物生產服務功能的經濟價值(元/hm2)；M為研究區每年的糧食總價值(元)；H為研究區耕地面積(hm2)。

為了減少不同年份物價差異帶來的影響,以2015年糧食平均收購價為基準,運用CPI指數對不同年份價值進行修正。

VC=Xij×E

(4)

式中,VC為生態系統服務價值系數(元/hm2)；Xij為i類生態系統第j類生態系統服務的當量因子。耕地、林地、草地、水域、未利用地的價值系數參考謝高地的當量因子表計算,耕地對應農田、林地對應森林、草地對應草地、水域對應河流/湖泊、未利用地對應荒漠。建設用地參考文獻取值[27]。

2.2.2生態系統服務價值系數修正

(1)空間異質性修正

生態系統服務的強弱直接影響生態系統服務價值,而生態系統服務與生物量有著密切聯系[26]。一般情況,生物量越大,生態系統服務也就越強[28]。因此,本文用NPP反映生物量并作為空間異質系數進行修正[15]。生產力的計算要考慮到數據獲取難易程度且又能說明變化,本文選用Thornthwaite Memorial模型[29]。計算公式如下：

(5)

NPP=3000×[1-e-0.0009695×(V-20)]

(6)

(7)

L=300+25×T+0.05×T3

(8)

式中,N表示空間異質系數；NPP表示植被凈初級生產力(kg hm-2a-1)；NPP′表示研究區植被凈初級生產力；NPP″表示全國植被凈初級生產力；V表示實際蒸散量(mm)；r為年降水量(mm);L為平均蒸發量(mm)；T為年平均氣溫(℃)。

(2)社會經濟因素修正

各地區的自然和社會經濟條件存在差異,導致社會經濟發展水平和生態系統服務以及人們對生態系統服務價值的認知具有差異性[30]。因此,本文從支付意愿和支付能力兩個方面對生態系統服務價值系數進行修正。計算公式如下：

D=R×W

(9)

式中,D為社會經濟調整系數；R為ESV的支付意愿；W為居民對ESV的支付能力。

居民的支付意愿會隨著社會經濟發展水平發生變化。利用皮爾生長曲線模型可以表示社會發展水平[31]。計算公式如下：

(10)

(11)

E′n=Ea×Pa+Eb×Pb

(12)

支付能力是個人或者國家經濟能力的表現[32]。利用研究區和全國人均國內生產總值的比值進行ESV支付能力修正。計算公式如下：

(13)

式中,GDP′表示研究區人均國內生產總值； GDP″表示全國人均國內生產總值。

(3)資源稀缺性修正

區域的人口越多,人口密度越大,人們對于資源的需求也就會越大。因此,利用人口密度反映系數修正[33]。計算公式如下：

(14)

式中,F表示岷江上游資源稀缺系數；g為岷江上游的人口密度；G為全國人口密度。

綜上,2000—2015年的生態系統服務價值調整系數見表1。

表1 生態系統服務價值調整系數

修正后的生態系統服務價值計算公式如下：

ESV=∑(A×VC×N×D×F)

(15)

式中,ESV生態系統服務價值(元)；A為各類土地利用類型的面積(hm2)；VC年份生態系統服務價值系數(元/hm2)；N空間異質系數；D社會經濟調整系數；F資源稀缺系數。

2.3 生態系統服務價值的敏感性分析

對敏感性進行分析的目的是確定時間變化后價值系數對生態系統服務價值的影響程度。計算公式如下：

(16)

式中,CS是敏感性指數；VCnk和VCmk分別為調整前后k地類的生態系統服務價值系數(元/hm2)；ESVnk和ESVmk分別為調整前后k地類的生態系統服務價值(元)。

若CS<1,說明ESV對VC不敏感,缺乏彈性,即結合區域特征修正后的ESV是合理[34]。

2.4 土地利用變化模擬

元胞自動機(Cellular Automata,CA)是一種基于不連續的時空動力學模型,其特點是時間、空間和狀態都是離散的,具有強大的空間運算能力,可以有效地模擬系統的空間變化[35-36]。其模型可表示為：

S(t+1)=f[S(t),U]

(17)

式中,S為元胞有限、離散的狀態集合；U為元胞的鄰域；t和t+1表示不同的時刻；f為局部空間元胞狀態的轉化規則。

馬爾科夫(Markov)模型,是根據馬爾科夫隨機過程理論形成的方法,研究和分析隨機事件的變化規律,并預測未來變化[37]。其模型可表示為：

S(t+1)=St×Zij

(18)

式中,St和St+1為該時刻土地利用系統的狀態；Zij為狀態轉移矩陣。

CA-Markov模型綜合了CA模型所具有的空間動態演化優勢及Markov模型的長期預測優勢,可以有效地表現土地利用時間和空間的變化情況[38]。

3 結果分析

3.1 土地利用變化

3.1.1模擬驗證

本文利用CA-Markov模型,對耕地、建設用地進行限制,選用地貌、氣象、距離因素作為影響土地利用變化的驅動因素,包括海拔、坡度、降水、氣溫、GDP、人口、河流距離、道路距離和中心城鎮距離,制作適宜性圖集并對土地利用格局進行預測。為了驗證CA-Markov的準確性,以2010年為基年,預測了岷江上游2015年土地利用格局,并與實際解譯的2015年土地利用狀態進行比較,結果見圖2。檢驗得出,總體Kappa指數為0.9462,隨機指數為0.9578、位置指數為0.9543、分層區位指數為0.9543,說明預測結果可信度較高。通過檢驗后,以2015年為基年,迭代次數為20,選擇5×5連續性濾波器,模擬2035年岷江上游土地利用格局(圖3)。

圖2 2015年模擬和解譯土地利用圖Fig.2 Simulated and interpreted land use in 2015

3.1.2土地利用類型面積變化

由表2和圖3可知,岷江上游土地利用類型以林地和草地為主,建設用地和未利用地面積最少。2000—2015年林地面積持續減少,年均下降438.738 hm2,2015—2035年林地面積年均增加8 593.691 hm2。建設用地和耕地面積持續增加,尤其是耕地在2015—2035年年均增加了2456.901 hm2。草地、水域和未利用地面積呈現波動變化,2000—2015年草地面積年均增加33.906 hm2,但在2015—2035年大幅度下降,年均下降11066.008 hm2。水域和未利用地在2010—2015年均分別下降59.706 hm2和40.554 hm2,在2015—2035年下降81.729 hm2和0.810hm2。

圖3 岷江上游2000—2035土地利用變化Fig.3 Land use change in the upper reaches of the Minjiang River from 2000 to 2035

3.1.3土地利用轉移圖譜分析

由圖4可知,林地和草地在與其它土地類型轉換中較為劇烈。2000—2015年林地轉出到草地、耕地面積較大。建設用地的增加主要來源于林地和耕地的轉入,耕地的增加主要來源于林地和草地的轉入。因此,這一階段主要是林地面積在減少。未利用地的增加主要來源于林地和草地的轉入,說明區域的林地和草地在一定程度上出現退化的現象。從空間上看,林地轉耕地和耕地轉建設用地主要分布在松潘縣和汶川縣東部,林地轉草地主要分布在松潘和黑水縣,林地轉建設用地和草地轉耕地主要分布在松潘和茂縣。草地和林地轉未利用地主要分布在松潘縣。

表2 岷江上游2000—2035年土地利用變化/hm2

2015—2035年,建設用地的增加主要來源于草地和耕地的轉入,草地轉建設用地主要分布在松潘縣和汶川縣,耕地轉建設用地主要分布在松潘和理縣。耕地的增加主要來源于草地和水域的轉入,大致沿著岷江上游河谷低地分布,主要分布在黑水、茂縣、松潘縣以及理縣和汶川縣東部。林地增加,從轉移來看得益于草地和水域的轉入,主要分布在松潘縣境內。水域主要轉入林地、耕地和建設用地。

總體來說,無論是2000—2015年還是預測的2035年,建設用地和耕地的持續增加與林地、草地、水域的轉出關系最為密切,呈現出正相關。岷江上游地形復雜,低地和人口多集中在河谷地區,這些區域成為了耕地、建設用地主要的侵占地,也說明隨著岷江上游經濟的發展,城鎮化和工業化的快速推進,人類活動占用生態空間,導致植被破壞,生態退化。因此,隨著人口的增加,協調人類需求與生態保護會變得愈加重要。

圖4 岷江上游2000—2035年土地利用轉移Fig.4 Land use transfer in the upper reaches of the Minjiang River from 2000 to 2035

3.2 生態系統服務價值的變化

3.2.1生態系統服務價值時間變化

由表3可知,所有土地利用類型產生的ESV中,林地和草地貢獻最大,對岷江上游ESV有著重要影響。2000—2015年岷江上游的ESV總體增加了254756.926×104元,一方面雖然這一階段林地面積減少,但轉出主要是向價值較高的生態系統,如草地。另一個方面社會經濟調整系數和空間異質系數的變化對其造成影響。2010—2015年ESV下降了166743.189×104元,這與林地、水域減少,建設用地、耕地增加帶來的生態負效應以及單位面積農田食物生產的經濟價值差異有關。2015—2035年,ESV增加了56494.557×104元,這與林地面積增加帶來的生態正效應密切相關。

由表4可知,各單項ESV由高到低依次為：維持生物多樣性>保持土壤>氣體調節>水文調節>氣候調節>原材料生產> 廢物處理>提供美學景觀 >食物生產。說明岷江上游生態系統在維持生物多樣性、保持土壤、氣體調節、水文調節和氣候調節這五個方面發揮著重要作用,在未來仍然是區域生態系統提供的重要功能以及保護的對象。總體上,2000—2015年和2015—2035年各單項ESV都呈現增長的趨勢。但是,2015—2035年各單項ESV增長的幅度明顯低于2000—2015年,說明未來研究區土地生態系統提供各種服務的能力在降低。

表3 2000—2035岷江上游各土地用類型的生態系統服務價值/104元

3.2.2生態系統服務價值空間變化

由圖5可知,以村界為基礎,根據自然斷點法從低到高對研究區2000—2035年ESV劃分了5個等級。ESV較高的地區集中在岷江上游西部,尤其是松潘縣的西北部,黑水縣、理縣和汶川縣西部。這些區域海拔高,林地和草地覆蓋面積較大,人口及經濟活動少。岷江上游的ESV低值區主要沿著河谷地帶延伸,集中在中東部地區。由于岷江上游地形地質條件復雜,河谷地帶地勢低平,光熱條件較好,人口集中,人為擾動較大,土地利用變更也較西部更快。因此,這些區域ESV較低。

由圖6可知,對2000—2035年岷江上游ESV變化做冷熱點分析。2000—2015年岷江上游ESV增加熱點區域主要分布在西部地區,損失冷點區域集中分布在中東部河谷地帶。熱點地區主要分布在人類活動較少的山區,冷點地區主要分布在人類活動較為頻繁的河谷地帶,說明人類活動對ESV的變化影響較大。2015—2035年ESV增加熱點區域集中在西北部,損失冷點區域仍然分布在東部,相比2000—2015年冷點區域向理縣東部和汶川東部集中,黑水東部的冷點區域有一定規模的縮小。

綜上, 2035年ESV的高低值分布大致與2000—2015年范圍相同,說明2035年的預測數據能在一定程度上反映岷江上游ESV的情況。并且岷江上游中東部地區和河谷地帶是人口集中,人類活動最為頻繁的地區,也是低值ESV和損失冷點的集中分布區,尤其是黑水、汶川和茂縣的東部地區。因此,這些區域的生態修復是岷江上游地區生態建設的重點。

表4 2000—2035岷江上游單項生態系統服務價值表/104元

圖5 2000—2035年岷江上游生態系統服務價值分布Fig.5 Distribution of the value of ecological services in the upper reaches of the Minjiang River from 2000 to 2035

3.3 敏感性分析

本文對各類土地利用類型的生態系統服務價值系數分別上下調整 50%,根據敏感性指數公式,分別求出 2000年、2005 年、2010 年、2015年和2035年的敏感性指數。由表5可知,其中林地和草地在各時期的敏感性指數均較高,說明兩者的變化對ESV影響很大。主要原因是林地和草地面積在所有地類中占比最大。不同土地利用類型的敏感性指數,在同一時期具有較明顯的差異。但是每一種土地利用類型在不同時期,敏感性指數變化都比較小,且都小于1,說明價值系數的調整對生態系統服務價值影響不大,估算出的生態系統服務價值結果適用于岷江上游區域。

表5 2000—2035年岷江上游敏感性指數

4 結論與討論

4.1 討論

岷江上游作為高原和盆地的過渡地帶,地質結構復雜,自然資源豐富,生態環境脆弱,其地理環境在川西高原極具代表性。在氣候變化和人類活動的雙重作用下,山地水土要素趨向失調、失衡,對區域生態服務產生影響[39]。因此,加強過渡帶研究不僅具有科學內涵[40],對于長江上游生態屏障建設具有重要戰略意義。目前,國內少數研究者對岷江上游的ESV進行了研究[41- 45]且與本文部分結果呈現出一致性。比如林地和草地對ESV的貢獻率最大,林地ESV變化最大,ESV的空間分布呈現西高東低的趨勢等。本文在當量因子法的基礎上進行了系數修正并且從時空尺度對土地利用格局和ESV動態趨勢進行研究,發現土地利用的變化會導致ESV的協同變化,但修正系數的差異也會對ESV產生影響。比如,2000—2005年和2000—2015年林地面積減少幅度較大,但ESV增加幅度也較大,造成這種現象的主要原因是修正系數的變化。因此,合理確定ESV的修正系數,對真實反應研究區ESV的變化具有重要的作用[46]。

林地、草地是岷江上游高ESV的代表[41- 43],且單項功能價值的變化與草地、林地面積變化的相關性大,為了確保岷江上游生態系統發揮生態平衡的作用,建議對這兩種地類進行保護或修復,引導形成合理且利于提高ESV的土地利用格局。以耕作和城市化等人類活動為主導的土地利用方式往往導致區域ESV呈現下降或者減少趨勢[47]。岷江上游河谷低地是耕地和建設用地增加的主要區域也是ESV低值分布的主要區域。因此,建議通過制定合理的國土空間規劃[46],加強生態保護紅線的管控,開展生態系統修復工作,對ESV損失較為嚴重區域的生態系統進行修復并加以保護,從而提高整個區域生態系統提供服務的能力。

就研究方法而言,本文在使用當量法的過程中,為了減少當量表與研究區現狀的脫節,結合區域實際情況進行了系數修正。雖然綜合考慮了影響ESV的自然和社會經濟因素,但岷江上游是青藏高原東緣高山峽谷區,自然地理環境復雜,影響其ESV的因素有很多,如,地貌、自然災害、政策因素等[48]。因此,還需要進一步的探索,以便構建更為精準的ESV估算模型。基于預測的土利用數據模擬ESV變化的過程中,2035年的修正系數沿用了2015年。雖然修正系數會對ESV產生影響,比如2000—2015年ESV的變化。但由于經過敏感性分析,生態系統服務價值對系數缺乏彈性。因此,其影響在可接受的范圍內,仍然能對2035年ESV的變化做一定程度的反映。

4.2 結論

本文利用自然和社會經濟修正因子,構建了符合岷江上游實際情況的ESV估算模型。并且利用CA-Markov模型,模擬了2035年岷江上游土地利用空間格局。從時空尺度分析了西南生態脆弱山區—岷江上游土地利用及ESV的動態變化。具體結論如下：

(1) 2000—2035年,岷江上游土地利用類型以林地和草地為主。建設用地和耕地面積的持續增加與林地、草地和水域的轉出呈現正相關,草地、水域和未利用地面積呈現波動變化。從空間上看,未來,岷江上游人類活動的主要區域仍然在河谷地區。

(2)岷江上游生態系統服務中,林地和草地產生的ESV最高。2000—2015年岷江上游的ESV增加了254756.926×104元,2015—2035年ESV增加了56494.557×104元,其變化與地類的變化密切相關,尤其是林地、草地和耕地的變化。各單項生態系統服務中,維持生物多樣性、保持土壤、氣體調節、水文調節和氣候調節這五項服務在岷江上游生態系統中起主要的作用,是岷江上游生態系統保護應該關注的服務類型。

(3)從空間上看,2000—2035年ESV較高區域和變化熱點區域集中在岷江上游西部,岷江上游的ESV低值區和變化冷點區域主要沿著河谷地帶延伸,集中在東部地區。未來,中東部及河谷地區的生態保護應該加以重視。

(4)通過敏感性分析,各土地利用類型的敏感性指數均小于1,且每一類土地利用類型的敏感性指數變幅較小,說明生態系統服務價值對價值系數不敏感、缺乏彈性,研究結果可信,并且林地和草地是生態系統服務價值變化的敏感因子。