基于地理探測器的青藏高原多年凍土分布影響因子分析

肖 瑤， 趙 林， 鄒德富， 劉世博， 馬 露， 應 雪， 劉藝闐

（1.中國科學院西北生態環境資源研究院冰凍圈科學國家重點實驗室藏北高原冰凍圈特殊環境與災害國家野外科學觀測研究站，甘肅蘭州730000；2.中國科學院大學，北京100049）

0 引言

從凍土物理學角度出發，多年凍土是巖石圈-土壤-大氣圈熱質交換過程中相互作用形成和發展的［1］，自然界的許多因素參與了這一過程。在全球尺度上，氣候的地帶性規律控制著多年凍土的變化和分布差異；在區域尺度和局地尺度上，植被、土壤、積雪以及微地形等的差異對多年凍土分布的影響在一定條件下會超過大的氣候背景成為驅動多年凍土空間分布的主要因素，造成多年凍土在空間上的非氣候帶分布規律［2-6］。青藏高原發育著全球中低緯度地區海拔最高、面積最大的多年凍土區［7］，其多年凍土面積達1.06×106km2［8］，在區域氣候和地質條件等因素的綜合影響下，青藏高原多年凍土分布特征具有顯著的空間差異。因此，研究多年凍土分布的影響因素是區域凍土研究的重要內容之一，建立多年凍土空間分布與環境因子之間的關系也是多年凍土預測、建模必不可少的基礎。

溫度是影響多年凍土分布的主要因素，是驅動土壤熱狀態的主要因子［9-10］。目前，對青藏高原多年凍土分布的研究主要通過模型建立溫度與多年凍土分布之間的關系來實現。如凍土制圖模型中，傳統的經驗統計模型（高程模型、年平均地溫模型）以及簡單半經驗-半物理機制模型［凍結數模型、TTOP（Top Temperature Of Permafrost）模型］都將海拔、氣溫或地溫作為驅動，采用多年凍土分布下界或者設定閾值來判定多年凍土的分布范圍［11-14］，但是，這些模型模擬過程中未充分考慮影響多年凍土分布的其他局地因素。例如，程國棟［3］發現坡度、坡向等局地因素可以通過影響地溫進而影響多年凍土的空間分布。Keller等［15］提出的山地多年凍土分布模型（PERMAKART），考慮到坡度和坡向等地形因子的影響，模擬了瑞士地區山地多年凍土的分布。Frauenfelder等［16］在山地多年凍土分布模型基礎上考慮到積雪、植被等因素的影響，模擬了阿爾卑斯山地區多年凍土分布。以上研究結果都是建立在足夠的調查資料基礎上，但實際上大部分多年凍土分布于氣候惡劣、可達性差的地區，調查資料稀少，而隨著遙感技術的應用和發展，科學家們可通過衛星遙感數據快速的獲取不同區域的下墊面信息，如K??b［5］利用數字地形模型（DTM）、可見光遙感以及微波遙感等數據源，綜合考慮與多年凍土分布密切相關的因素，對多年凍土災害進行評估和管理。Etzelmüller［17］考慮影響多年凍土分布相關的因子來確定蒙古國庫蘇古爾地區多年凍土存在的可能性大小。總的來說，局地因素通過影響地氣能量交換和到達地表的太陽輻射，從而造成多年凍土空間分布上的差異［18-20］，而且不同區域環境下，其主導影響因子也有較大差異。同理，相同因子作用在不同地區也會造成對多年凍土分布的差異。目前的研究多針對各因子如何影響多年凍土分布，但是如何定量描述驅動因子對多年凍土分布影響程度，尤其是地表局地因素對多年凍土分布影響的定量化研究還存在不足。地理探測器能夠探測空間分異性，以及揭示其背后驅動力的一種新的統計學方法，其基于如果某個自變量對某個因變量有重要影響，那么自變量和因變量的空間分布應該具有相似性［21］，其不僅可以定量分析影響多年凍土空間分布各因子的驅動力相對大小，還能探測各驅動因子的交互作用。目前，地理探測器模型能夠應用到自然科學、環境科學、社會科學等各個方面［22-23］。

確定多年凍土的主要影響因子及貢獻率是構建多年凍土空間分布和預測模型的基礎。本文利用可見光和熱紅外遙感技術獲取與多年凍土分布有關的地表溫度、植被、積雪和地表反照率等數據，利用地理探測器模型探測不同因子對各區域多年凍土分布的影響程度，以及對影響多年凍土空間分布因子的強弱進行定量分析，從而為不同環境條件下多年凍土分布規律和區域多年凍土建模制圖提供決策。

1 數據與方法

1.1 研究區概況

青藏高原凍土分布圖采用Zou等［8］利用TTOP模型的模擬結果。在此基礎上選取五個研究區，分別為東部溫泉地區、西部西昆侖地區、南部改則地區、北部阿爾金地區和中部不凍泉地區。這些區域均分布于多年凍土與季節凍土的過渡地帶（圖1），多為不連續多年凍土。

圖1 研究區分布Fig.1 Spatial distribution of the research areas on the Qinghai-Tibet Plateau

各研究區內地形和氣候特征存在較大差異，溫泉地區海拔介于3 430～5 300 m之間，地形以山地和丘陵為主，區域內年平均氣溫在1.3℃以下，年降水量約400～600 mm［24-25］。不凍泉地區地形呈現明顯起伏，青藏鐵路斜穿而過，區內海拔多在4 000 m以上，氣象站資料顯示年平均氣溫為-6.2℃，年平均降水量為270 mm［26］。阿爾金地區以丘陵和山間盆地為主，地勢從西南向東北呈現降低趨勢，海拔在3 400～5 400 m之間，年平均氣溫低于3.6℃，年降水稀少，極為干旱［27］。西昆侖地區地勢總體表現為東北低，中部多為高山，南部較緩，區內海拔除去圖中東北區域外，多集中在5 000 m以上，為研究區內平均海拔最高的區域，冰川分布廣泛，年平均氣溫在-4.5℃以下，年平均降水量在60～120 mm之間［28］。改則地區地勢大致為西北高，東南低，無明顯起伏，海拔在4 400～6 300 m之間，區內地形以山地和山間盆地為主，區內年平均氣溫為0℃左右，年平均降水量為150 mm［29］。

1.2 數據

地表溫度作為地表能量平衡中的主要參數，是影響多年凍土分布的主要因子，而其他局域因子也可以通過影響地表能量和輻射平衡，從而影響多年凍土的空間分布。本研究以遙感反演產品為基礎，選取對多年凍土分布影響較大的六個因子，包括地表溫度、積雪、植被、地表反照率、坡度和坡向。青藏高原上，尤其在多年凍土區，植被以高寒草地為主，冬季植被較少，且冬季容易形成穩定的積雪，因此采用7—9月NDVI和冬半年積雪日數為植被和積雪的指標。數據指標和數據來源見表1，其中地表溫度、積雪、植被和反照率均采用MODIS數據產品，該產品來自于美國國家航空航天局（NASA），數據下載網址為https：//search.earthdata.nasa.gov/，時間范圍為2003—2012年。

表1 數據來源Table 1 Sources of the used data

1.3 方法與模型

本研究利用地理探測器模型，確定不同因子對各研究區多年凍土分布的影響程度。地理探測器模型［21，33］能探測同一區域內變量的相似性、不同區域間變量的差異性。假設自變量X和因變量Y的空間分布趨于一致，那么它們之間存在統計關聯性，進而可以揭示其因果關系。兩個變量空間分布的一致性、空間分異性可以用地理探測器中的q值度量，而q值可直接用地理探測器模型計算（http：//www.geodetector.org/）。本文所使用地理探測器功能包括［34］：（1）因子探測；（2）交互作用探測，地理探測器模型的公式表達如下：

式中：h為影響多年凍土分布因子的分類，h∈[1，L]；Nh為層h的單元數，N為研究區全部單元數和σ2分別是變量某一層h方差和全區方差。q∈[0，1]，q值越大說明分層異質性越明顯；由于多年凍土分布是由于各個因素共同作用的結果，q值越大表示影響因子對多年凍土分區的解釋力越強，相反則越弱；q=1時表明影響因子完全控制多年凍土的空間分布；q=0時表明多年凍土的空間分布不受該因子影響。本研究中Y為多年凍土，X為影響因子，基于地理探測器模型，依據不同區域多年凍土分布的結果，通過計算并比較各因子q值大小，分析其對多年凍土的影響程度。

交互探測分別計算因子X1和X2的q值q（X1）和q（X2），然后計算因子X1和X2交互作用的q值q（X1∩X2）；根據q（X1）、q（X2）和q（X1∩X2）的關系，交互作用可分為非線性減弱、單因子非線性減弱、雙因子增強、獨立、非線性增強（表2）。本文基于交互作用的探測結果，評估在雙因子共同作用下影響多年凍土空間分布的能力是否增強［21，35］。

表2 交互作用關系Table 2 Redefined interaction relationships

2 結果與分析

2.1 典型多年凍土區影響因子對比與分析

對溫泉地區、西昆侖地區、改則地區、阿爾金地區和不凍泉地區五個研究區的多年凍土分布以及研究區內2003—2012年平均地表溫度（LST）、冬半年積雪日數（ASW）、7—9月植被指數（NDVI）、地表短波反照率（Albedo）、坡度（Slope）和坡向（Aspect）進行分類并可視化（圖2）。

山區和低溫區域容易形成積雪，其地表反照率較高［36］。高海拔地區其植被覆蓋較低［37］，其地表反照率較高。由于冰川和湖泊具有高反照率、低NDVI等特點，因此去除冰川和湖泊的影響，對研究區內各因子的空間分布（圖2）和箱線圖（圖3）進行對比分析：西昆侖地區高度差最大，同時其地表溫度、冬半年積雪日數和地表反照率相差也最大，區內年平均地表溫度為研究區內最低（-0.68℃），冬半年平均積雪日數最大（46 d），年平均地表反照率最高（0.27）。不凍泉地區年平均地表溫度為-0.38℃，僅高于西昆侖地區；溫泉地區雖然年平均地表溫度僅低于改則地區，但區內多山地，其海拔相差也較大，最高坡度和平均坡度為研究區內最高，從而形成溫泉地區和不凍泉地區冬半年平均積雪日數（37 d）僅次于西昆侖地區，但是由于區內降水量高于其他研究區（在青藏高原，降水量與植被成正相關［38］），而造成這兩個地區夏季年平均NDVI遠高于其他研究區，分別為0.42和0.32，且集中分布在0.2以上，使其區內地表反照率低于其他研究區；改則地區和阿爾金地區年地表反照率特征相似，多集中在0.2～0.3之間；改則地區年平均地表溫度為研究區內最高（1.22℃），其冬半年平均積雪日數（16 d）為研究區內最低，其平均海拔為研究區最高（4 997 m）；阿爾金地區海拔雖然多分布在5 000 m以下，但其降水量為研究區內最低，造成NDVI為研究區最低，且集中在0.1以下。

圖2 影響因子的空間分布Fig.2 Spatial distribution of influencing factors

圖3 不同區域環境因子箱線圖Fig.3 The box-plot of environment factors in the different areas

2.2 單因子對典型區多年凍土分布的影響

氣候背景直接影響到地表溫度［39］，各種地質地理因素對多年凍土的影響主要表現為影響到達地面的熱流大小［1］，其相互作用下導致地氣系統能量變化在空間上的不同，造成多年凍土在空間上的分異［40］。在不同區域，各局地因素對多年凍土分布的影響程度也有較大差異。在選取的因子中，探測結果（表3）顯示，地表溫度都為影響各區域多年凍土分布的主要因素，并結合相關分析得到研究區各因子與多年凍土分布之間的相關性（表4），分析討論得出地表溫度越低的區域多年凍土分布越廣泛。根據探測的各因子q值，以0.3為標準進行劃分，探討各局地因素對多年凍土分布的影響強弱，表明：積雪日數、NDVI和地表反照率對溫泉地區和西昆侖地區多年凍土分布影響較大，改則地區和不凍泉地區影響多年凍土分布較強的因子為積雪日數，阿爾金地區多年凍土分布除了受地表溫度影響較強，其他局地因素都表現為弱影響。

表3 驅動因子探測Table 3 The q-statistic of geographical factors

表4 相關分析Table 4 The r coefficient of correlation analysis

植被是影響多年凍土熱狀況的主要因素之一，植被能減少到達地表的太陽輻射，降低地面溫度差，減少了進入土壤中的熱量［1］。NDVI對溫泉地區和西昆侖地區的多年凍土分布影響較強，而對改則地區、阿爾金地區和不凍泉地區的多年凍土分布影響較弱。分析其原因為：植被對地表層熱狀況影響取決于所在地區植被分布特征［41］，在溫泉地區和西昆侖地區，植被與地表溫度的空間分布存在相似性，即地表溫度高的地區，植被發育越好，與兩地主要以山地地形為主有很大的聯系，而地表溫度特征與多年凍土的分布具有明顯的一致性。而在改則地區、阿爾金地區和不凍泉地區，山地地形分布不集中，NDVI與地表溫度在空間上的分布無規律。在溫泉和西昆侖地區，NDVI和多年凍土的分布關系表現為NDVI越小的區域，多年凍土分布越廣泛；而在其他地區表現為NDVI越大的區域，多年凍土分布越廣泛。

積雪的作用主要有：阻礙地表與大氣之間熱交換起隔熱作用、使一部分熱量不能達到地表起冷卻作用、融化時吸收大量的熱量起冷卻作用，高原上不同區域雪蓋對淺層地溫的影響存在較大差異，冷季積雪的綜合作用是削弱地表的冷卻作用，對淺層地溫起到了保溫作用［1，41］。積雪對研究區多年凍土分布的影響表現為：溫泉地區和西昆侖地區影響最強，不凍泉地區和改則地區次之，阿爾金地區影響最弱。分析其原因為：阿爾金地區平地較多且分布集中，溫泉地區和西昆侖地區多為山地地形，海拔相差大；不凍泉地區和改則地區平地分布雖然多，但不集中，積雪日數特征和多年凍土分布在空間上表現為積雪日數越多的區域多年凍土分布越廣泛。

地表反照率是一種綜合指標，受到NDVI和積雪的影響，積雪和NDVI在溫泉地區和西昆侖地區對多年凍土分布的影響較強，造成地表反照率對溫泉地區和西昆侖地區的多年凍土分布影響也較強，對改則地區、阿爾金地區和不凍泉地區的影響較弱。

本研究中使用的各遙感數據產品和多年凍土分布數據分辨率都為1 km，以及各個研究區選取的尺度范圍內，坡度和坡向對研究區多年凍土分布的影響都表現為弱影響。因此，尺度問題在本文研究中不可忽略。

2.3 不同尺度下因子對多年凍土分布的影響

以溫泉地區、阿爾金地區、改則地區和西昆侖地區不同半徑內的范圍作為不同區域尺度（圖4、表5），探究不同尺度下各局地因子對多年凍土分布的影響程度（圖5）。

圖4 研究地區不同尺度范圍Fig.4 The different scales in the research areas

在選擇的四個實驗區域內，可以發現隨著尺度范圍的增大，地表溫度對多年凍土分布的影響程度都有明顯的增強趨勢，側面反映了在大尺度上，多年凍土的分布主要受溫度影響。局地因素在不同區域，隨著尺度的變化，對多年凍土分布的影響有著不同的規律。而對于坡度和坡向，隨著尺度范圍的增大，其對多年凍土分布的影響程度都呈現減弱趨勢。坡向和坡度直接影響到地面接收太陽輻射的強度。而青藏高原的太陽輻射強烈，致使山坡方位作用增強，因此在小范圍尤其是山地多年凍土分布時，各局地因素與坡度和坡向的共同作用不可忽略。

表5 不同尺度劃分依據Table 5 The principle of different scales

圖5 不同尺度下的q值探測Fig.5 The q-statistic of geographical factors in different scales in the research areas

2.4 雙因子對典型區多年凍土分布的影響

使用交互探測分析兩個因子共同作用下對不同區域多年凍土分布的影響程度，并對單因子探測的結果進行對比，選取各個區域的交互作用最強的前三位（表6）和每個區域內與該因子最大的交互值與交互類型（表6，表7），可以得出：雙因子交互作用下的影響程度都表現為雙因子增強和非線性增強。交互作用主要揭示影響多年凍土分布的因子在共同作用下的影響程度要大于任何一個單一因子的影響程度。在研究區內，各交互作用最強的都是與地表溫度的交互，在西昆侖地區、改則地區和不凍泉地區，交互作用最強的為地表溫度和NDVI的交互，在溫泉地區和阿爾金地區，交互作用最強的為地表溫度和積雪日數。說明各個局地因子共同作用下，通過直接或間接的影響地表淺層熱狀況，從而影響所在區域內地表溫度進而造成多年凍土分布的差異。

3 結論

多年凍土的空間分布是多個因素共同作用的結果，大的氣候背景是多年凍土分布的主要控制因素，但是在局域尺度上，植被、積雪、地形等的差異也會對凍土分布的影響起到重要作用。本研究基于地理探測器模型，以地表溫度、NDVI、積雪日數、地表反照率、坡度和坡向等影響因子為自變量，青藏高原典型區多年凍土為因變量，分析各因子以及雙因子間交互作用對多年凍土分布的影響程度，結論如下：

表6 各地區交互探測和類型（最強的前三位）Table 6 The interaction relationships in research areas（the top three strongest）

表7 各地區交互探測和類型（與該因子最大的交互值）Table 7 The interaction relationships in research areas（maximum interaction value with this factor）

（1）因子探測結果中，五個典型多年凍土研究區內，地表溫度為影響多年凍土分布的最主要因子，其q值均大于0.4；交互探測結果表明兩個變量相互作用下對多年凍土的分布影響程度要大于一種變量的獨立影響，增強類型為雙因子增強和非線性增強。

（2）本文選取的各局域因子對研究區內多年凍土分布的影響不同。溫泉地區和西昆侖地區積雪日數、NDVI和地表反照率的影響最強；改則地區和不凍泉地區積雪日數的影響最強；阿爾金地區各局地因子的影響都比較弱。此外，坡度和坡向在研究區為弱影響因子，q值均小于0.1。同一因子對不同研究區多年凍土分布的影響也存在較大差異，其差異與各區域該因子分布特征有關。

（3）探測不同區域尺度下各因子對多年凍土分布的影響程度，隨著尺度的增加，地表溫度的影響程度逐漸增大，反映了在大尺度上溫度是影響多年凍土分布的主要因素，而坡度和坡向的影響程度逐漸減少，表明在小尺度范圍內尤其是山地多年凍土分布研究時，坡度和坡向等局域因素對多年凍土分布的影響不可忽視。

本文目的是使用地理探測器模型探究不同因子對凍土分布的影響程度，只選取了部分影響因子，如果需要進行進一步深入探討影響多年凍土分布的局域因素，需要加入更多局域因子如土壤含水量、積雪厚度等，因變量可以選取活動層厚度等與凍土密切相關的參數，靈活利用地理探測器模型和其他方法，綜合考慮尺度問題，研究局地因素對多年凍土分布的影響強弱，為區域凍土模擬和制圖提供更好的決策基礎。