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

肖 瑤， 趙 林， 鄒德富， 劉世博， 馬 露， 應(yīng) 雪， 劉藝闐

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

0 引言

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

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

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

1 數(shù)據(jù)與方法

1.1 研究區(qū)概況

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

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

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

1.2 數(shù)據(jù)

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

表1 數(shù)據(jù)來源Table 1 Sources of the used data

1.3 方法與模型

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

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

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

表2 交互作用關(guān)系Table 2 Redefined interaction relationships

2 結(jié)果與分析

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

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

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

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

圖3 不同區(qū)域環(huán)境因子箱線圖Fig.3 The box-plot of environment factors in the different areas

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

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

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

表4 相關(guān)分析Table 4 The r coefficient of correlation analysis

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

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

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

本研究中使用的各遙感數(shù)據(jù)產(chǎn)品和多年凍土分布數(shù)據(jù)分辨率都為1 km，以及各個研究區(qū)選取的尺度范圍內(nèi)，坡度和坡向?qū)ρ芯繀^(qū)多年凍土分布的影響都表現(xiàn)為弱影響。因此，尺度問題在本文研究中不可忽略。

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

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

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

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

表5 不同尺度劃分依據(jù)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 雙因子對典型區(qū)多年凍土分布的影響

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

3 結(jié)論

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

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

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

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

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

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

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