基于植被生產力的黃土高原地區生態脆弱性及其控制因子分析

楊艷 ，周德成，宮兆寧 ，劉子源，張良俠*

1. 南京信息工程大學應用氣象學院，江蘇 南京 210044；2. 首都師范大學資源環境與旅游學院，北京 100048；3. 資源環境與地理信息系統北京市重點實驗室，北京 100048

生態環境是人類賴以生存的基本條件。近年來，隨著全球氣候變化和人類活動加劇，生態系統的自我調節能力快速下降，呈現出脆弱化的趨勢，嚴重威脅著區域生態安全和經濟社會的可持續發展（Hu et al.，2021；徐廣才等，2009；徐興良等，2022）。依據IPCC第五次評估報告，生態脆弱性被定義為生態系統容易受到不利影響的傾向或者習性，它包括生態系統對不利影響的敏感性以及對新的外界條件的適應性兩個方面（IPCC，2014；蘇勝亮等，2022）。科學合理的對區域生態脆弱性進行量化評估，對理解生態系統的功能和過程以及因地制宜的進行生態系統修復和保護具有重要意義。

生態脆弱性評價的關鍵在于構建科學的評價指標體系。當前相關研究采用的評價指標體系大致可以分為綜合指標體系和單一指標體系兩種（楊飛等，2019）。其中，綜合指標體系大多基于不同的評估模型框架，從自然風險、人為干擾和社會經濟等方面選取多個指標構建而成（張學玲等，2018）。同時，在進行脆弱性評估過程中，往往需要依據不同的方法（如主成分分析、層次分析和專家打分法）對選取的多個指標進行權重分配（鄢繼堯等，2020）。例如，張學淵等（2021）基于“壓力-狀態-響應”評估模型框架，從氣候、土壤、植被、地形地貌和社會經濟等方面選取了 13個指標，利用空間主成分分析法對 2000—2018年西北干旱區生態脆弱性進行了評估；朱琪等（2021）基于“成因及結果表現”評估模型框架，從地形、地表、氣象和社會經濟等方面選取 13個指標，結合空間主成分分析法評估了 2005—2015年東北森林帶生態脆弱性狀況；金麗娟等（2022）基于“敏感-恢復-壓力”評估模型框架，選取了 17個涵蓋地形、地表、氣象、結構、功能和壓力等多方面的評價因子，并基于層次分析法對 2005—2018年四川省進行生態脆弱性評價。然而，由于考慮范圍廣和指標數量多，綜合指標體系存在指標因子權重分配具有主觀性以及因子之間相關性高等問題。相較而言，采用能夠綜合反映區域生態脆弱性的主要因子建立單一指標評價體系具有較強的可靠性和可操作性（Zhang et al.，2017；肖桐等，2010）。

生態系統生產力是生態系統功能的重要指標和生態系統承載力的物質基礎，能夠有效反映出生態系統的健康狀態，對外界環境變化具有高度敏感性，是衡量生態脆弱性的重要依據（于貴瑞等，2020）。國內外學者基于生態系統凈初級生產力（net primary productivity，NPP）和總初級生產力（gross primary productivity，GPP）已成功地對不同地區的生態脆弱性進行了評估。如以NPP作為評價指標，分別對1989—2007年的三江源地區、1961—2000年長江中下游地區、2000—2010年岷江上游地區進行了生態脆弱性評估（Zhang et al.，2017；肖桐等，2010；於琍等，2012）；以 GPP為評價指標，分別對2000—2015年西南地區以及2000—2014年天山—塔里木河綠洲地區生態脆弱性進行了評估（何敏等，2019；王鶴松等，2021）。這些研究工作有助于提高對生態脆弱性問題的理解，但在探究生態脆弱性控制因子時仍存在一定局限性。如通過對各控制因子（如溫度、降水、高程、坡度）作分區統計，來分析其對生態脆弱性的影響（肖桐等，2010；王鶴松等，2021）；分別基于沿各控制因子的變異梯度隨機選取的40個采樣點數據和基于整個研究區的所有數據利用線性回歸分析方法，分析每個控制因子與生態脆弱性的相關關系（Zhang et al.，2017；何敏等，2019）。受研究方法的限制，這些研究僅初步分析了區域尺度上各控制因子對生態脆弱性的單一影響，難以剔除各控制因子間的相互作用，更不能定量表達各控制因子對生態脆弱性的貢獻率。地理探測器是探測空間分異性并揭示其背后驅動力的一組統計學方法，該方法不需要對影響因子預設限定，且無線性假設，計算得出的q值具有明確的物理含義，能夠客觀反映出自變量對因變量的解釋力，因此能夠較全面的探究生態脆弱性時空分布格局的控制因子（王勁峰等，2017）。

黃土高原處于干旱、半干旱和半濕潤區的過渡帶，植被退化和水土流失問題嚴重，相對落后的經濟、不合理的土地利用以及煤炭資源的開發更加劇了該地區生態環境的退化，屬于中國典型的生態脆弱區（賀鵬等，2022；張家政等，2022）。本研究以黃土高原為研究區，選擇NPP作為評價指標，依據IPCC中有關生態脆弱性的定義首先計算敏感性和適應性，進而對 2001—2020年黃土高原生態脆弱性的空間分布格局進行了評估，并結合地理探測器定量分析了生態脆弱性的控制因子及其交互作用。研究結果有助于深化對該地區生態系統對氣候變化和人類活動響應的認識，并為當地采取針對性措施以減少外部脅迫對生態系統產生的不良影響、促進可持續發展和規劃提供科學參考和理論支持。

1 材料和方法

1.1 研究區概況

黃 土 高 原 （ 33°43′— 41°16′N ， 100°54′—114°33′E）地處中國中部偏北，位于黃河中上游和海河上游地區，總面積約6.40×105km2（劉靜等，2020），覆蓋了青海、甘肅、寧夏、內蒙古、陜西、山西和河南等7個省區（圖1）。該區具有冬季寒冷干燥、夏季濕潤溫暖的暖溫帶大陸性季風氣候特征，年均溫度4—12 ℃；降水分布不均勻，自西北向東南逐漸增多，年均降水量400—600 mm（陳劍南等，2022）。區域內以黃土覆蓋為主，土體疏松，抗侵蝕能力較弱，從而形成了溝壑縱橫支離破碎的典型黃土地貌。土地利用類型以草地、耕地和林地為主，面積占比分別為41%、32%和15%。植被由東南向西北可劃分為森林帶、森林草原帶、典型草原帶、荒漠草原和草原化荒漠帶（修麗娜等，2019）。能礦資源豐富，是全國能源安全保障基地和生態安全屏障的重要組成部分（王佳寧等，2020）。

圖1 2010年黃土高原地區土地利用情況Figure 1 Land use in the Loess Plateau in 2010

1.2 數據來源與預處理

本文使用的數據主要包括遙感、氣象、高程、土地利用/覆被、人均GDP、人口密度和坡度數據。遙感數據包括NPP和歸一化植被指數（normalized difference vegetation index，NDVI）兩種，時間跨度為 2001—2020年，均從美國地質調查局網站（https://lpdaac.usgs.gov）下載獲取。其中，NPP來自MODIS數據集中的MOD17A3H產品，時間分辨率為年，空間分辨率為500 m；NDVI來自MODIS數據集中的MOD13A2產品，時間分辨率為16 d，空間分辨率為 1 km。本研究利用最大值合成法將NDVI數據由16 d尺度轉換為年尺度。氣象數據為覆蓋研究區及其周圍148個氣象站點的2001—2020年的年均氣溫和年均降水量，來源于中國氣象科學數據共享網（http://data.cma.cn）。本研究基于氣象觀測數據，利用AUSPLIN插值軟件結合數字高程模型（digital elevation model，DEM）進行空間插值，生成1 km×1 km空間分辨率的氣象要素柵格圖像。DEM、2010年土地利用情況、人均GDP和人口密度數據均來自于中國科學院資源環境科學數據中心(http://www.resdc.cn)。其中，DEM數據空間分辨率為90 m，土地利用、人均GDP和人口密度數據空間分辨率為1 km。坡度數據由DEM計算得到。為了便于分析不同指標因子與生態脆弱性之間的關系，本研究將所有數據統一重采樣為Albers投影下1 km×1 km的柵格數據。

1.3 脆弱性計算方法

選取NPP作為生態脆弱性評價指標，根據IPCC第五次評估報告中有關脆弱性的概念進行生態脆弱性的計算（IPCC，2014）。生態脆弱性被定義為生態系統容易受到外界破壞或傷害的程度，可由敏感性和適應性的差值計算得到：

式中：

V——生態脆弱性；

S——敏感性，表征了生態系統受氣候變化或其他外界擾動的影響；

A——適應性，表征了生態系統受自身調節和恢復能力的制約。

S由NPP的年際波動情況來表示，計算公式如下：

式中：

Fi——該地區第i年的NPP值（n=20）；

——黃土高原NPP的平均值。

A由NPP年際變率線性擬合趨勢線的斜率來表示，計算公式如下：

式中：

x——自然數 1、2、3…，對應于 2001—2020年的時間序列；

y——NPP年際變率的絕對變化量。

值得注意的是，由公式（2）和（3）計算得到的S和A數值不在同一量綱，難以直接用于數據分析和生態脆弱性計算。因此，本研究基于極差標準化方法（Hu et al.，2021）對S和A數據進行了標準化處理。

為了進一步揭示黃土高原地區生態脆弱性空間變異特征，需要對生態脆弱性指數進行分級。常用的數據分級方法包括自定義間隔、等距、標準差和自然斷點等，其中自然斷點分級法是一種根據數值統計分布規律分級和分類的統計方法，可使得各級內部方差之和最小且差異最大化（Liu et al.，2017）。本研究采用自然斷點分級法確定生態脆弱性指數分級閾值，將其劃分為微度脆弱、輕度脆弱、中度脆弱、重度脆弱和極度脆弱5級。

1.4 控制因子分析

為進一步探究黃土高原生態脆弱性的控制因子，選取年均溫、年均降水量、海拔、坡度、人均GDP、人口密度和NDVI 7個指標（表1），利用地理探測器模型定量分析了探測器包括因子探測、交互作用探測、風險區探測和生態各因子對生態脆弱性空間分布的影響。利用地理探測器模型定量分析了各因子對生態脆弱性空間分布的影響。地理探測器包括因子探測、交互作用探測、風險區探測和生態探測四個部分（王勁峰等，2017）。其中因子探測利用因子解釋力衡量自變量因子對因變量變化產生的貢獻力大小，可以反映出黃土高原生態脆弱性的主要控制因子；交互作用探測通過比較兩個指標的單因子解釋力之和及其交互作用的解釋力，可用于探測自變量因子兩兩之間的交互后對生態脆弱性的影響程度。地理探測器模型主要用q值度量解釋力結果，其值越大說明自變量對于因變量空間分異的解釋力越強，反之則越弱。首先，將生態脆弱性指數V作為因變量Y，將選取的7個控制因子指標作為自變量X，借助ArcGIS 10.7軟件在15 km×15 km格網尺度統計脆弱性和各自變量因子平均值；其次，利用自然斷點法將各自變量因子重分類為5級，由數值量轉變為類型量，并獲取脆弱性數據與因子分類數據的空間關聯屬性表；最后將屬性表導入地理探測器運行得到q值，對影響黃土高原生態脆弱性空間分布的控制因子進行探測分析。

表1 影響因子Table 1 Impact Factor

2 結果

2.1 黃土高原地區生態脆弱性空間分布格局

2001—2020年黃土高原生態脆弱性整體較高，中度及以上脆弱等級地區面積占比為61%（表2），且呈現出西北高而東南低的空間分布格局（圖2c）。對比而言，適應性空間分布與脆弱性剛好相反，呈現西北低而東南高的格局，而敏感性呈現中部高而東部及南部低的格局（圖2a、b）。具體而言，微度和輕度脆弱區所占面積分別為 1.05×105km2和1.46×105km2，主要分布在黃土高原東南部，如南部的六盤山和子午嶺地區，以及東部的呂梁太行山區，具有敏感性低且適應性極高的特點；重度和極度脆弱區所占面積則分別為 1.49×105km2和 6.80×104km2，主要分布在黃土高原中部的丘陵溝壑區及西南部的高塬溝壑區，敏感程度高且適應性較低。此外，不同土地利用類型間的脆弱性差異也較大（圖3）。林地整體脆弱性較低，以微度、輕度和中度脆弱為主，面積占比超過 85%；其次為建設用地、耕地和草地，中度及以下脆弱區面積占比分別為76%、66%和 59%；未利用地的脆弱性程度最高，以中度、重度和極度脆弱為主，面積占比超過86%。

圖2 2001—2020年黃土高原地區敏感性（a）、適應性（b）和生態脆弱性（c）空間分布特征Figure 2 Spatial distributions of sensitivity, adaptability and ecological vulnerability in the Loess Plateau in 2001-2020

圖3 2001—2020年黃土高原不同土地利用類型下各生態脆弱性等級面積占比Figure 3 Area proportion of each ecological vulnerability under different land use types in the Loess Plateau in 2001-2020

表2 黃土高原不同脆弱性等級所占面積及占比Table 2 Area and proportion of different vulnerability levels in the Loess Plateau

2.2 黃土高原地區生態脆弱性空間格局的控制因子分析

7個影響因子的q值大小排序為NDVI＞降水＞坡度＞溫度＞海拔＞人均GDP＞人口密度（表3）。其中NDVI解釋力最大，且遠高于其他指標，q值為0.59；其次為降水，q值為0.48；坡度次之，q值為0.18；海拔、人均GDP和人口密度的影響最小，q值均為0.09。這表明植被覆蓋度和年降水量變化對2001—2020年黃土高原生態脆弱性變化的解釋程度相對較大，其他因子對生態脆弱性變化的解釋程度相對較低。

表3 黃土高原地區各影響因子對生態脆弱性的解釋力Table 3 Explanatory power of different drivers on ecological vulnerability in the Loess Plateau in 2001-2020

生態脆弱性因子交互作用值均大于單個因子的最大值（表 4），這說明各指標對于黃土高原生態脆弱性的影響并不是相互獨立，而是相互作用的結果。海拔與降水、海拔與人均GDP、海拔與人口密度、坡度與降水和坡度與溫度之間表現為非線性增強，其余都表現為雙因子增強。雖然單因子探測結果顯示海拔對黃土高原地區生態脆弱性影響較小，但是當它與NDVI交互作用時，對生態脆弱性具有最強的解釋力（q=0.66）。綜上，NDVI和降水分別與其他因子交互作用后產生的解釋力都遠高于其他因子之間的相互作用（q＞0.51）。

表4 黃土高原地區生態脆弱性各影響因子間的交互作用Table 4 Interactions among drivers of ecological vulnerability in the Loess Plateau

3 討論

本研究表明，2001—2020年黃土高原地區生態脆弱性水平整體較高，且呈現出明顯的西北高而東南低的分布格局，該結果與張良俠等（2022）基于“暴露-敏感性-適應力”評估模型框架，結合層次分析和空間主成分分析方法得出的黃土高原生態脆弱性評估結果相似，且與劉會軍等（2015）在全國尺度上界定的生態脆弱區相一致。極度脆弱區主要集中在黃土高原西南部的高塬溝壑區，該地區氣候干旱且植被覆蓋稀疏，水土流失嚴重，生態脆弱性程度最高，與陳楓等（2018）和蘇勝亮等（2022）的研究結果相一致。重度和部分極度脆弱區主要分布在陜北、晉北和鄂爾多斯東部地區，其中包含一些重要的煤炭基地。該區位于干旱和半干旱地帶，煤炭能源的開采不僅破壞了地表原本稀疏的植被且進一步加劇了水資源的短缺，從而導致生態環境高度脆弱，這與王佳寧等（2020）的研究結果相一致。中度脆弱區分布在黃土高原中北部的毛烏素沙地和庫布齊沙漠地區，該地區沙化草原面積占比大，生態系統對氣候變化敏感且抵御外界干擾能力差（Sun et al.，2021）。然而，本研究結果與部分前人研究存在一些不同，例如楊雯娜等（2021）發現黃河流域中的極度脆弱區集中在寧夏北部以及楊晴青等（2019）發現陜西佳縣的生態脆弱性為中等脆弱，產生差異的原因可能是受研究尺度、指標選取和評價方法等因素的影響，生態脆弱性評價結果是相對性的，區域之間不具有可比性，這也說明了在區域尺度上建立系統客觀的評價體系的重要性。

本研究發現黃土高原地區生態脆弱性主要的兩個控制因子為植被覆蓋度和降水，這與張良俠等（2022）和王麗霞等（2021）在黃土高原區的研究結果具有一致性。降水是影響干旱和半干旱區植被生長發育的重要因素，植被覆蓋度高意味著其生態功能較強，二者的增加客觀上都有利于生態脆弱性的降低（Li et al.，2021；陳玉蘭等，2022）。例如，在黃土高原東南部，較高的年降水量適宜植被生長，生態脆弱性等級整體較低。而在中部半干旱區，植被覆蓋度較低，生態系統較脆弱。研究時段內黃土高原地區年降水量整體呈增加趨勢（劉洋洋等，2019），同時大規模實施的生態恢復工程顯著提高了植被覆蓋度（Xiao，2014；楊丹等，2022），這兩者均會降低黃土高原地區的生態脆弱性。然而，何敏等（2019）發現西南地區生態脆弱性主導因子為溫度和海拔，李路等（2021）發現氣溫、地形和植被覆蓋度是新疆喀什地區生態脆弱性的主要控制因子，張學淵等（2021）發現西北干旱區生態脆弱性受土壤和地形影響最大，產生分歧的原因主要在于所選研究區不同，自然條件存在差異。相較于西南地區或者新疆喀什地區，黃土高原所在區域整體熱量條件的變化以及生態系統受地形的制約程度相對較低，因而溫度和海拔對于黃土高原生態脆弱性的影響較小。本研究表明，海拔與植被覆蓋度相互作用對黃土高原生態脆弱性具有最大的解釋力。

相較于以往的研究（張學淵等，2021；朱琪等，2021；金麗娟等，2022），本研究采用的將NPP作為單一指標的生態脆弱性評價方法，回避了綜合指標體系中的指標選取具有重疊性和權重分配主觀性高等問題，能夠有效評估黃土高原地區生態脆弱性的空間分布格局。基于該方法得到的結果與前人研究結果相一致，證明本研究方法是科學可行的。地理探測器的引入也為定量識別生態脆弱性的控制因子提供了評估方法。然而，本研究仍存在一些不足。首先，本研究基于2001—2020年NPP的動態變化特征來計算生態脆弱性，結果只能表征研究區生態脆弱性的空間變異特征，無法反映其時間變異特征。其次，受數據限制，本研究時段相對較短，且所選生態脆弱性驅動因子沒有考慮其他人為活動（如生態工程、耕作活動和牲畜數量），退耕還林和退牧還草等生態政策對黃土高原脆弱性的影響還有待更深入的研究。未來還需進一步完善生態脆弱性的計算方法，利用長時序高空間分辨率數據，并結合研究區的社會經濟發展狀況，選取系統綜合的指標進行生態脆弱性的主控因子分析，以期為黃土高原生態保護和建設提供具體的科學參考和理論支持。

4 結論

基于2000—2020年黃土高原NPP數據，依據IPCC中有關生態脆弱性的定義，對黃土高原生態脆弱性的空間分布格局進行評估，并結合地理探測器定量分析了生態脆弱性的控制因子及其交互作用，研究主要得到以下結論：

（1）2001—2020年，黃土高原地區生態脆弱性整體較高，中度及以上等級脆弱區所占面積比例約為61%，且呈現西北高東南低的空間分布格局。

（2）不同土地利用類型間的生態脆弱性差異較大，林地脆弱性程度較低，草地、耕地和未利用土地的脆弱性程度相對較高。

（3）黃土高原生態脆弱性的主要驅動因子是植被覆蓋度和降水，且所有指標因子間存在較強的交互作用，植被覆蓋度和海拔的交互作用對生態脆弱性具有最強的解釋力。