肖 瑶， 赵 林， 邹德富， 刘世博， 马 露， 应 雪， 刘艺阗

（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）探测不同区域尺度下各因子对多年冻土分布的影响程度，随着尺度的增加，地表温度的影响程度逐渐增大，反映了在大尺度上温度是影响多年冻土分布的主要因素，而坡度和坡向的影响程度逐渐减少，表明在小尺度范围内尤其是山地多年冻土分布研究时，坡度和坡向等局域因素对多年冻土分布的影响不可忽视。

本文目的是使用地理探测器模型探究不同因子对冻土分布的影响程度，只选取了部分影响因子，如果需要进行进一步深入探讨影响多年冻土分布的局域因素，需要加入更多局域因子如土壤含水量、积雪厚度等，因变量可以选取活动层厚度等与冻土密切相关的参数，灵活利用地理探测器模型和其他方法，综合考虑尺度问题，研究局地因素对多年冻土分布的影响强弱，为区域冻土模拟和制图提供更好的决策基础。