岷江上游生態脆弱性時空分異及驅動因子交互效應分析
——以雜谷腦河流域為例

肖成志，計揚，李建忠，張志，巴仁基，曹亞廷

1. 中國地質調查局軍民融合地質調查中心，四川 成都 610036；2. 中國地質大學（武漢）地球物理與空間信息學院，湖北 武漢 430074；3. 中國地質大學（武漢）工程學院，湖北 武漢 430074

岷江上游雜谷腦河流域是典型的山地型生態脆弱區（劉鐵冬，2014），是中國重要的生態功能區，生態屏障戰略地位十分突出，對維護長江流域生態安全，促進長江流域經濟社會可持續發展具有舉足輕重的作用（師嘉祺等，2021；汪璇等，2022）。該區域地形起伏大，易發生崩塌、滑坡泥石流等地質災害（余恩旭等，2023），特別是干旱河谷區地質條件、自然條件非常復雜，生態十分脆弱，經濟發展非常落后（陳小平等，2014），加上2008 年5?12 汶川地震影響，使得該區域生態脆弱性更加典型。因此了解岷江上游地震前后生態環境狀況對該區域生態修復與重建至關重要。

生態系統脆弱性是一個生態系統在時間和空間上承受壓力的能力（Williams et al.，2000）。早期的研究主要基于單一風險因子，以單維度、單要素的評價為主，隨著脆弱性評價一步研究，逐漸發展為多維度、多要素的綜合評價（李平星等，2014），其過程主要是選取生態脆弱性指標，在此基礎上建立區域生態脆弱性的評估指標體系，然后進行全面評估。評價方法主要有專家打分法、層次分析法、空間主成分分析法、熵權法、物元法、人工神經網絡法等（Niu et al.，2020；Qiao et al.，2020；Rezaeisabzevar et al.，2020）。生態系統作為一個動態系統，具有一定的自我調節機制，即適應能力（肖桐等，2010），不同特征生態系統受到的自然過程（如降水、極端事件等）和人類活動影響（如森林砍伐、綠化、土地使用變化等）不同（Wanyama et al.，2021）。如Polsky et al.（2007）認為人類和自然環境不是獨立的系統是緊密耦合的，并且對威脅的暴露、敏感和適應能力不同，“脆弱性”的認知需要對人類與環境相互作用的生物物理、認知和社會維度進行調查，并在此基礎上提出了暴露度-敏感性-適應能力（Vulnerability Scoping Diagram，VSD）模型。Wang et al.（2008）將生態脆弱性指標分為“土地資源條件”、“水熱氣象條件”和“地形條件”三大類。Kang et al.（2018）將生態脆弱性指標分為“自然環境”和“社會經濟”兩個代表性類別。以上研究大多以自然環境和社會經濟方面來研究生態脆弱性，但對地質條件在生態脆弱性研究中起到的作用研究較少，未單獨進行討論。因此本文基于生態地質學理論和雜谷腦河流域生態本底條件，從“自然環境-地質地貌-地質災害-社會環境”4 個維度選取評價指標，以探究復雜地質條件對研究區生態脆弱性的作用，為該區域綠色發展提供支撐。

1 研究區及數據

1.1 研究區概況

雜谷腦河流域（圖1）位于青藏高原東緣，阿壩藏族羌族自治州東南部，岷江干流西部。地理位置31°26′—31°33′N，103°9—103°36′E，東西長約94 km，南北寬約56 km，流域面積為4 629 km2。區內山巒綿延、層嶺疊嶂、山高谷深、地勢西北高、東南低，區內坡度較陡多在20°—60°之間，流域內地貌類型為：侵蝕、深切河谷；龍門山剝蝕、侵蝕高山、中山；岷山、邛崍山剝蝕、侵蝕高山；岷山、邛崍山冰蝕、冰緣極高山、高山。氣候相對復雜，垂向氣候特征顯著，水熱環境具有明顯的區域性和垂向變異特征，特別是干旱河谷區降水量較少。植被類型主要有針葉林、闊葉林、灌草，及櫻桃樹、李子樹等。土壤類型主要有暗棕壤、棕壤、褐土、高山草甸土、亞高山草甸土、高山寒漠土等。

圖1 研究區概況示意圖Figure 1 Overview of the study area

1.2 數據源及處理

本研究數據包括數字高程模型（Digital Elevation Model，DEM）、Landsat 數據、土地利用數據、潛在蒸發數據、降雨數據、人口密度數據、地層巖性數據、地貌分區數據、地質災害點數據、干旱河谷區邊界數據、道路矢量數據和土壤類型數據。其中，DEM 數據為美國航空航天局（National Aeronautics and Space Administration，NASA）和日本宇宙航空研究開發機構（Japan Aerospace Exploration Agency，JAXA）共同開發的數字高程模型，可在地理空間數據云（https://www.gscloud.cn/）下載。Landsat 數據來源于美國地質勘探局（https://earthexplorer.usgs.gov/），考慮到數據質量、植被生長旺盛時期、1998 年洪水及2008 年5?12 地震影響，震前選擇了2001 年6 月13 日數據，震后選擇了2009 年6 月9 日數據，現狀選擇了2020 年7 月19 日數據，道路矢量數據由Google Earth 軟件解譯得到，其他數據來源如表1 所示。所有數據空間參考均統一為WGS_1984_UTM_48N，其空間分辨率統一至30 m×30 m 大小計算。

表1 數據源Table 1 Data Sources

2 研究方法

2.1 評價指標

根據指標體系構建原則，本文從“自然環境-地質地貌-地質災害-社會環境”4 個維度選取了14 個生態脆弱性指標，指標計算和分級方法如表2 所示，生態脆弱性指標量綱不同，為統一量綱本文將所有指標分為5 級（圖2、3）。有參考標準的參考相關標準和前人研究成果，無參考標準的根據現場實際調查情況或自然間斷點法分級如表3 所示。

表2 雜谷腦河流域生態脆弱性評價指標體系Table 2 Ecological vulnerability assessment index system in the Zagunao River basin

表3 雜谷腦河流域生態脆弱性評價指標分級標準Table 3 Classification criteria of ecological vulnerability indexes in the Zagunao River basin

圖3 穩定指標空間分布特征Figure 3 Spatial distribution characteristics of stability indicators

2.2 空間主成分分析

空間主成分分析可以通過對原始空間軸的旋轉，將相關空間變量轉換為少數幾個互不相關的綜合指標，從而實現對數據結構簡化的目的（Abson et al.，2012；李志學等，2021），因此該方法可有效解決生態脆弱性評價指標之間相關性過大和指標冗余問題（孫宇晴等，2021）。本文選取累積貢獻率在85%以上的主要成份，以其為替代因子，來評估生態脆弱程度，其計算公式為：

式中：

Wj——第j個替代指標（主成分）；

Yi——第i個初始指標的標準化值；

Zij——第i個初始指標在第j個主成分所對應的特征向量；

n——初始指標個數。

主成分分析需要對相關指標進行中心化和歸一化處理。一般可將指標分為正面指標和負面指標兩類。對于正向指標，指標值越大，脆弱性就越強烈，而對于負向指標，指標值越大，脆弱性就越低。最后生態脆弱性指數計算公式如下：

式中：

E——生態脆弱性指數；

Wj——第j個主成分；

Qj——第j個主成分的貢獻率；

m——累計貢獻率大于85%的前幾個主成分個數。生態指數指數值的大小表示生態脆弱性程度，其值越大生態脆弱性越高生態環境越差，其值越小生態脆弱性越低，生態環境越好。為進行多年生態脆弱性的比較分析，對EI 進行歸一化處理，歸一化公式如下：

式中：

V——標準化值；

Ei——生態脆弱性指數第i個區域內的值；

Emin和Emax——生態脆弱性指數在研究區的最小值和最大值。

參考已有分級標準采用等間距法將生態脆弱性指數分為5 級：潛在脆弱（0—0.2）、微度脆弱（0.2—0.4）、輕度脆弱（0.4—0.6）、中度脆弱（0.6—0.8）、重度脆弱（0.8—1）。生態脆弱性值越大表示生態脆弱性越強，生態環境質量越差，生態脆弱性值越小表示生態脆弱性越低，生態環境質量越好。

2.3 地理探測器模型

地理探測器是一種基于數據挖掘和統計學方法的空間分析工具，旨在分析和解釋不同因素對空間分布現象的影響，包括因素探測、相互作用探測、風險探測和生態探測4 個模塊（Ren et al.，2022）。本文基于因子探測模塊和交互探測模塊對雜谷腦河流域的驅動因子進行探討。因子檢測的基本目的是識別地理差異的原因，即自變量X對生態脆弱性變化屬性Y的解釋力，并用下式進行計算：

式中：

q值——影響因素解釋力；

h=1, 2, 3...,L——影響因子x分層，即分類或分區；

Gh和G——h層和整個區域的樣本數；

W和T——層內方差和全區的方差。q的范圍是0—1，q的值越大，因變量的解釋力就越大。在極端情況下，q值為1 表示自變量的空間完全分異，而q值為0 則表示自變量隨機分布。

交互檢測可以定量表示兩個驅動因素對生態脆弱性的分布格局的相互作用，表4 中概述的交互檢測類型可用于識別各驅動因素之間如何交互作用。

表4 交互探測類型Table 4 Types of interactive detection

2.4 交互效應變異比

驅動因子相互作用效應可以看作是間接效應，是指兩個驅動因子相互作用時對因變量的綜合作用。交互效應變異比（Interactive Effect Variation Ratio，IEVR）是一種量化驅動因子的個體效應（直接影響）與交互效應（間接影響）之間差異的模型（Deng et al.，2021）。驅動因子IEVR 值越大驅動因子對因變量間接影響越強，IEVR 值越小驅動因子對因變量直接影響越強。其原理是計算某影響因子獨立影響和與其他因子交互影響值的差值，然后計算該因子獨立影響值和交互影響值所有差值的均值，表達式如下：

式中：

Rxi——變量xi與其他自變量的交互效應相對于其獨立影響力的平均變異率；

n——因子個數；

P(xi)——xi的獨立影響力；

P(xi∩xj)——變量xi和xj的交互影響力。

3 結果分析

3.1 生態脆弱性時空分布及演變

本文基于14 個生態脆弱性評價指標進行空間主成分分析，選取累計貢獻率大于85%的主成分如表5 所示，基于公式（2）、公式（3）計算EVI，按照生態脆弱性分級標準得到生態脆弱性等級圖（圖4），由圖4 可知研究區生態脆弱性變化活躍的區域主要為域干旱河谷區以及西部、北部高山極高山區域，生態脆弱性分布特征為西部和北部高于中部和東部，生態脆弱性隨海拔升高呈現出先減弱后增強趨勢。計算2001、2009 和2020 年3 個時期生態脆弱性面積占比（圖5），3 個時期研究區均為潛在脆弱區面積占比最大面積占比最大分別為46.4%、41.6%和43.5%，重度脆弱區面積占比最小分別為7.1%、7.5%和5.3%，微度脆弱、輕度脆弱和中度脆弱區面積占比接近在13%—21.8%之間。圖6 為各時期生態脆弱性轉移?；鶊D，顯示出2001—2009 年生態脆弱性主要變化類型為潛在脆弱轉微度脆弱、微度脆弱轉輕度脆弱、輕度脆弱轉中度脆弱和中度脆弱轉重度脆弱。2009—2020 生態脆弱性主要變化類型為，潛在脆弱轉微度脆弱、微度脆弱轉潛在脆弱、輕度脆弱轉微度脆弱和中度脆弱轉輕度脆弱。

表5 各主成分的特征值、貢獻率及累計貢獻率表Table 5 Eigenvalue, contribution rate and cumulative contribution rate of principal component

圖4 2001、2009 和2020 年雜谷腦河流域生態脆弱性分級圖Figure 4 Ecological vulnerability classification diagram of the Zagunao River Basin in 2001, 2009, 2020

圖6 不同時期生態脆弱性轉移?；鶊DFigure 6 Sankey diagram of ecological vulnerability transfer in different periods

圖7a 為干旱河谷（a1），侵蝕、深切河谷（a2），龍門山剝蝕、侵蝕高山、中山（a3），岷山、邛崍山剝蝕、侵蝕高山（a4）和岷山、邛崍山冰蝕、冰緣極高山、高山（a5）在3 個時期的生態脆弱性均值，結果顯示3 個時期差距較小但不同區域之間生態脆弱性具有明顯差異。其中干旱河谷區和岷山、邛崍山冰蝕、冰緣極高山、高山生態脆弱性較強（生態脆弱性多年均值分別為2.27 和3.29），侵蝕深切河谷、龍門山、剝蝕、侵蝕高山、中山和岷山、邛崍山剝蝕、侵蝕高山生態脆弱性弱（生態脆弱性多年均值分別為1.28、1.15 和1.19）。

圖7 干旱河谷區、不同地貌類型區生態脆弱性均值及面積比例Figure 7 Mean value and area proportion of ecological vulnerability in arid valley and different geomorphic types

從各區域生態脆弱性百分比來看（圖7b），重度脆弱區只分布在岷山、邛崍山冰蝕、冰緣極高山、高山且輕度和重度脆弱區分布最廣面積百分比為25.8%和36.3%。該區域內植被覆蓋度低，山頂主要為雪被和裸地雪線以下主要為高山草甸，海拔高自然環境惡劣環境，生態環境極其脆弱，生態受損后不易修復，因此該區域需要將強保護，防止區域環境繼續惡化，在保護的同時適地適樹植樹造林可能起到生態修復輔助作用。干旱河谷微度脆弱和輕度脆弱區分布最廣面積百分比為46.6%和33.9%，主要位于汶川縣城威州鎮至理縣縣城雜谷腦鎮。該區域人口密集，人類工程活動強烈，降雨較少，降雨時空分布不均，植被覆蓋度較低，植被類型主要為灌草以及果樹，巖性較為軟弱易發生崩塌滑坡泥石流等地質災害，多種因素綜合導致該區域生態脆弱性變化強烈，該區域在農業活動、工程建造、資源開發時應更加注重生態環境保護。侵蝕深切河谷、龍門山剝蝕侵蝕高山中山和岷山邛崍山剝蝕侵蝕高山主潛在脆弱區分布最廣面積百分比分別為74.5%、87.6 和83.8%生態環境較好。

3.2 生態脆弱性關鍵驅動因子分析

表6為因子探測結果表，結果顯示高程（0.862）、生物豐度指數（0.77）、植被覆蓋度（0.739）、水源涵養指數（0.724）和干旱指數（0.655）解釋力均處于較高水平是生態脆弱性關鍵驅動因子，其余因子解釋力均較低對研究區生態脆弱性空間分布影響較小。圖8 為因子交互探測結果，各驅動因子交互大多為雙因子增強和非線性增強不存在非線性減弱或者獨立的關系，雙因子交互q值均大于單因子q值，說明兩因子交互作用解釋力均大于單一因子解釋力，其中植被覆蓋度、干旱指數、水源涵養指數、生物豐度指數和高程與其他因子交互解釋力最高，進一步說明植被覆蓋度、干旱指數、水源涵養指數、生物豐度指數和高程對研究區生態脆弱性空間分布影響較大。

表6 生態脆弱性驅動因子探測表Table 6 Detection table of driving factors of ecological vulnerability

圖8 驅動因子交互探測熱力圖Figure 8 Thermodynamic diagram of driving factor interaction detection

3.3 關鍵驅動因子交互效應分析

圖9 為3 個時期生態脆弱性主要驅動因子交互效應變異比的變化，顯示出2001、2009 和2020 年植被覆蓋度IEVR 值別為2.41%、5.29%和18.73%，呈現逐漸增大趨勢，一方面是2001—2020 年區域內中等植被覆蓋度和較高植被覆蓋度區域面積不斷減小，另一方面是區域內土壤輕度侵蝕區域不斷增加，地質災害環境影響力不斷加強。干旱指數在2001 年和2020 年IEVR 值相差不大分別為15.2%和14.43%，但在2009 年呈現出較低水準（8.09%），這是因為2009年相較于2001 年和2020 年氣候變化更大，氣候較為干燥，特別是干旱河谷區降雨較少。水源涵養指數IEVR 值分別為8.05%、12.82%和29.98%，IEVR 值逐漸增大，這是因為區域內中等植被覆蓋度和較高植被覆蓋度區域面積不斷減小，而具有高質量生態系統的地區，如森林、灌木、草原和濕地等區域，地表徑流較少，保水能力相對較高（Duan et al.，2022）。生物豐度指數IEVR 值在2001 年和2020 年較小分別為4.71%和3.87%，但在2009 年顯著提升（8.17%），這是由于5?12 地震影響2009 年生物豐度和大大降低，地質災害等因素影響力大大加強，而生物多樣性決定了生態系統的抗變化能力，生物多樣性越單一生態就越脆弱。高程IEVR 值在3 個時期均最低分別為1.48%、1.96%和2.12%，這是因為山區生態環境質量與地質環境息息相關（Liu et al.，2023），研究區垂直分異明顯土壤類型和厚度、植被類型、降雨量、溫度均隨海拔的變化而呈現出不同特征。

圖9 不同時期主要驅動因子交互效應比變化圖Figure 9 Variation diagram of main driver IEVR in different periods

4 討論與結論

4.1 討論

1）本文選取14 個生態脆弱性評價指標，利用主成分分析法對雜谷腦河流域2001、2009 和2020年3 個時期生態脆弱性進行評價，并對3 個時期生態脆弱性變化進行分析。研究結果表明，受2008 年5?12 地震影響該區域植被大面積損毀地質災害多發，加劇了該區域水土流失損害了原有山地生態系統功能（吳寧等，2008），相對于震前震后區域內生態脆弱性增強。自5?12 地震以來政府部門實施了許多生態保護修復措施，從而促進了雜谷腦河流域植被恢復（姚昆等，2020），現今區域內生態脆弱性降低。研究結果與楊渺等（2020）對5?12 地震極重災區生態恢復進行評估結果類似。其中干旱河谷生態脆弱性較高和高寒地區生態脆弱性最高印證了前人干旱河谷區和高寒地區易形成生態脆弱去的結論。但本文生態脆弱性評價也存在局限性，一是指標分級存在主觀性，如高程分級是按照現場調查情況和植被垂直分異劃分，該標準只適用于雜谷腦河流域。二是本文評價指標較多，本文僅在評價方法上去除了評價指標的相關性，未對評價指標的相關性具體分析。因此，如何構建適用于不同區域的指標分級標準，深入分析各指標之間復雜關系避免評價指標冗余是下一步工作的重點方向。

2）研究區除主要影響因子外土壤侵蝕因素具有一定的影響力且不斷加強，這與生態環境因素中植被覆蓋度影響力逐漸減小相關，植被可以防止水土流失和土壤侵蝕，保持土壤濕度和養分促進水循環，如果土壤易于流失和侵蝕，水資源的凈化能力下降，生態系統服務能力和生產力都會受到影響，生態脆弱性就會增加。地質災害密度的影響力較小2009 和2020 年較2001 年其影響力稍有提高，這是因為地質災害點主要區域位于雜谷腦河沿線特別是汶川縣至理縣縣城地質災害點較多，但對于整個流域生態脆弱性而言影響較弱。人口密度和離道路距離在所有時期的影響力均處低水平且變化不大，主要因為研究區人口分布密度較大的區域主要為雜谷腦河河谷區域，特別是威州鎮至雜谷腦鎮區域人口最為集中，區域內城鎮建設和農業活動基本位于雜谷腦河河流階地以及緩坡區域，對于研究區其他區域影響力較小。地質災害與社會環境數據不夠豐富，在今后的研究中需進進一步補充和分析。

4.2 結論

基于雜谷腦河流域環境背景和人文特征選取指標，利用空間主成分分析法得出5?12 地震震前、震后以及現狀的生態脆弱性，利用地理探測器和交互效應變異比，得出了主要驅動因子及其交互效應。主要有以下結論：

1）雜谷腦河流域總體生態脆弱性較低以潛在脆弱區為主，但局部區域生態脆弱性較高。震后（2009 年）相對于震前（2001 年）雜谷腦河流域潛在和輕度脆弱區面積減少，微度、中度和重度脆弱區面積增加，生態脆弱性增強生態環境呈變差趨勢，現今（2020 年）相對于震后（2009 年）潛在脆弱、微度脆弱和輕度脆弱區域面積增加，中度脆弱、重度脆弱區域面積有不同程度的減少，生態脆弱性降低生態環境得到改善。生態脆弱性變化強列區域主要位于研究區東部干旱河谷和西部北部及高山區域。

2）雜谷腦河流域不同區域生態脆弱性分布特征明顯，東部干旱河谷區生態脆弱性相對較高，主要為微度脆弱和輕度脆弱區。重度脆弱區只分布在西部北部岷山、邛崍山冰蝕、冰緣極高山、高山，同時該區域輕度和中度脆弱區面積占比最大，生態脆弱性最高。其他區域潛在脆弱區面積占比最大生態脆弱性較低。

3）雜谷腦河生態脆弱性空間分異主要驅動因子為植被覆蓋度、干旱指數、水源涵養指數、生物豐度指數和高程。對比主要驅動因子各時期IEVR值，結果表明高程是研究區生態脆弱性最直接影響因子，植被覆蓋度和水源涵養指數對研究區生態脆弱性間接影響力不斷加強，干旱指數在氣候干燥時期對研究區生態脆弱性直接影響力增強，生物豐度指數對研究區生態脆弱性主要為直接影響，在震后對研究區生態脆弱性間接影響增強。