寒冻风化控制的祁连山风化碎屑的空间分布

洪 洋， 耿豪鹏，2， 潘保田，2

（1.兰州大学资源环境学院西部环境教育部重点实验室，甘肃兰州 730000；2.甘肃省石羊河流域野外科学观测研究站，甘肃兰州 730000）

0 引言

地貌是由岩体的抬升、变形和风化，以及地表的侵蚀、搬运和堆积等过程所塑造的［1］，其中风化作用是地貌侵蚀过程的先导［2］。寒冻风化是冰缘地貌区及高寒山区的主要物理风化过程［3-5］。20世纪80年代起，崔之久等通过野外定位观测探讨了寒冻风化的主要特点及影响因素［6-7］，发现寒冻风化提供的风化碎屑物质是冰缘地貌区及高寒山区崩塌、泥石流等灾害性地貌过程［8］及倒石堆的重要物质来源［9］，寒冻风化可以通过影响倒石堆的运动速率，进而破坏公路等基础设施［10］，由此开启了我国研究寒冻风化过程的序幕。

寒冻风化的作用机制存在两种不同的假说。第一种是冻融作用（freeze-thaw）相关的体积膨胀理论，即岩石孔隙中的水从液态变为固态时体积膨胀，对岩壁施加压力使岩体破碎，易于崩解［11］。冻融作用要求存在一个饱和且密闭的空间，使体积膨胀所产生的压力可以全部施加于岩壁［12］。但自然界中基本不存在饱和且密闭的孔隙，不饱和孔隙中体积膨胀所产生的压力仅为饱和孔隙的10%左右［13］，因此冻融作用对真实地貌的有效性一直存在争议［14］。为此学者提出分凝冰理论（segregation ice growth）来解释真实环境中的寒冻风化机制。分凝冰理论认为高寒环境下岩石孔隙中广泛存在分凝冰，由于界面融化和曲率效应，过冷水在冰-岩界面间形成薄膜［15］，使液态水从温度较高的区域流向分凝冰所在孔隙，促进分凝冰不断生长［16］，冰-岩界面间的范德华力及静电力使岩体破碎［17］。分凝冰理论在实验室得到了广泛的验证［18-20］，以分凝冰理论为基础构建的寒冻风化模型也广泛应用于计算不同区域的寒冻风化强度，进而解释高寒山区侵蚀速率的时空变化及现代地貌单元的空间分布特征［21-23］。因此本文中的寒冻风化指基于分凝冰理论的岩体破碎过程。

广泛分布的风化碎屑是高寒山区及冰缘地貌区寒冻风化的主要产物［24］，极易受多种因素诱发转化为滑坡、泥石流等地质灾害［25］。因此探明风化碎屑空间分布的主要控制因素对预防地质灾害的发生具有重要意义。由于分凝冰生长速率及液态水从高温向低温的流动过程均与温度相关［12，17］，因此温度是寒冻风化的主要控制因素。同时寒冻风化对岩体破碎的影响是缓慢而持续的累积过程，因此在探讨寒冻风化对现代地貌单元空间分布特征的影响及其在长期地貌演化中的作用时考虑过去温度变化（包含一个或多个冰期—间冰期旋回）的影响十分必要［26-28］。同时前人研究表明，末次冰盛期以来的温度变化计算的寒冻风化强度可以很好地约束风化碎屑空间分布的下界［29-30］。但上述研究均在新西兰南阿尔卑斯山、意大利阿尔卑斯山等较为湿润的山区开展，干旱半干旱高寒山区温度变化时间尺度的选取尚需进一步探究。

前期野外考察结合遥感影像发现祁连山北部风化碎屑物质广泛分布，且存在由风化碎屑转化产生的滑坡、泥石流等地质灾害。因此本文拟以祁连山北部为研究区，探讨干旱半干旱高寒山区风化碎屑空间分布的主要控制因素。基于以上目标，我们首先利用哨兵二号遥感影像，采取NDVI（归一化植被指数）阈值、海拔筛选结合目视解译的方法提取祁连山北部风化碎屑的空间分布；其次利用寒冻风化模型计算不同时间尺度温度变化情景下的寒冻风化强度的空间分布特征；最后对比分析风化碎屑及寒冻风化强度的空间相关性，探讨干旱半干旱高寒山区风化碎屑空间分布的主要控制因素。

1 研究区概况

祁连山位于青藏高原东北缘，晚新生代起开始快速抬升［31］。本研究所关注的祁连山北部海拔较高（平均海拔约3 700 m），坡度较陡（超三分之一区域坡度超过20°），年平均气温约2.3℃，年降水量在100～600 mm·a-1之间变化，且呈现明显的东高西低的分布趋势［32］。植被表现出明显的垂直地带性：2 000～2 300 m主要为荒漠草原带，2 300～2 600 m主要为温性草原带，2 600～3 200 m主要为森林草原带，3 200～3 700 m主要为灌丛草原带，3 700～4 000 m主要为高寒草甸带，4 000 m以上主要为冰雪带［33］；同时植被还表现出明显的水平差异，植被垂直带谱由东向西趋于简化［34］。祁连山北部是疏勒河、黑河及石羊河三大内陆河流域的发源地，其中共有58条河流的多年平均出山径流量大于0.1 m3·s-1［32］。本文选取祁连山北部乌鞘岭以西为研究区（图1）。

图1 研究区数字高程模型（Shuttle Radar Topography Mission DEM，http：//srtm.csi.cgiar.org/）及气象站点位置（黄色圆点）Fig.1 The digital elevation model（Shuttle Radar Topography Mission DEM，http：//srtm.csi.cgiar.org/）of the study area with locations of meteorological stations（yellow circles）

2 数据和方法

2.1 风化碎屑边界的提取

我们采取NDVI（归一化植被指数）阈值、海拔筛选结合遥感影像目视解译的方式进行风化碎屑边界的提取。风化碎屑与植被由于反射率不同在遥感影像上表现出的光谱特征差异是我们采用NDVI阈值提取与目视解译的基础。我们使用哨兵二号10 m分辨率的遥感影像（European Space Agency，https：//sentinel.esa.int/web/sentinel/home）用于计算研究区的NDVI，所用影像选取原则为2021年6—9月以最大限度去除积雪的影响，且云覆盖量低于10%。由于植被垂直带谱在东西方向存在差异，因此NDVI阈值的选取随植被特征的变化而变化，范围为0.09～0.15。谷歌地球0.6 m分辨率的遥感影像主要用于目视解译，对NDVI阈值提取结果进行修正，去除水体、云雾等因素干扰。同时为了剔除低海拔区域其他成因裸岩斑块的影响，我们依据研究区植被的垂直地带性特征结合野外考察、遥感影像目视解译及前人研究成果［35］，以风化碎屑萎缩或扩张变化的海拔下限3 800 m为阈值对提取结果进行筛选。我们利用30 m分辨率的数字高程模型（Shuttle Radar Topography Mission DEM，http：//srtm.csi.cgiar.org/）进行研究区海拔属性的提取。

2.2 寒冻风化模型

Hales等［36］以分凝冰理论为基础设计了寒冻风化模型，以寒冻风化强度（frost cracking intensity，Fci）作为评价温度对岩体破坏程度的指标。依据Hales等［36］的模型，我们先用一维热传导公式模拟温度随时间（t）和深度（z）的变化（图2）：

式中：T是不同时间（t，单位为天）不同深度（z，单位为cm）下的温度（℃）；TMA是年平均气温（℃）；Ta是年气温变幅的一半（℃）；α是热扩散系数（设定为1 mm2·s-1）；Py指气温变化的周期，在本研究中为1年。寒冻风化模型的构建存在三个假设前提：（1）当温度在-3～-8℃时（即寒冻风化窗）分凝冰生长速率最快［12］；（2）当温度高于0℃时存在液态水；（3）液态水从高温向低温单向流动。依据Hales等［36］的定义，寒冻风化强度（Fci，°C·cm-1）是寒冻风化窗内温度梯度随时间的积分［37］：

式中：∆T是温度梯度；n是符合三个假设前提的天数。同时，由于研究区属于干旱半干旱区，还需考虑液态水在流动过程中因无法及时补充造成的可能损耗，因此我们引入Anderson等［27］提出的距离函数（e-D/D*），以Fci与距离函数的乘积表征适用于干旱半干旱区的寒冻风化强度［37］。通过对修改后的寒冻风化强度进行深度积分，得到不同年平均气温条件下的总寒冻风化强度。依据气温和海拔之间的线性关系可得到不同年平均气温所对应的海拔，即可得到不同海拔条件下的总寒冻风化强度［37］：

式中：D是0℃所在位置到-3℃所在位置的距离；D*是长度指数，依据Anderson等［27］设置为50 cm。但研究发现某些日期温度随深度的变化并非单调的，而是先增加后减小（图2紫色实线），因此存在一个温度转折点，导致0℃同时出现在两个不同的深度。因此在本研究中，将D设定为转折点所在位置到-3℃所在位置的距离以保证温度随深度的单调变化。

图2 年平均气温为-3℃及Ta=15℃时温度随时间和深度的变化［灰色阴影为寒冻风化窗所在范围（温度在-3～-8℃时），垂直的黑色虚线为0℃所在位置；紫色实线（t=182）代表温度随深度的非单调变化］Fig.2 Schematic representation of the temperature variation as a function of depth using TMA of-3°C and Ta of 15°C［The grey shadow represents the frost cracking window，between-3 and-8°C；The black dashed line represents the location of 0°C；The purple line（t=182）is an example of a depth profile with two 0°C locations and an inflection point］

为了计算平均寒冻风化强度来评估不同时间尺度温度变化对寒冻风化的累积影响，我们从Vostok冰芯记录中以500年为时间间隔提取不同年代的年平均气温（TMA）并假定Ta不变，得到不同年代的总寒冻风化强度，取其平均值即为平均寒冻风化强度［37］：

式中：i是时间间隔的数量；Fci-ti是第i个时间间隔相对应年代的总寒冻风化强度。

2.3 温度数据

年平均气温可以通过11个气象站点获得（图1中的黄色圆点，表1），数据来源为国家气象信息中心（https：//data.cma.cn/），记录了1981—2010年的日平均气温的平均值。我们利用气象站记录的气温与站点海拔之间的线性回归关系来计算气温直减率及年平均气温0℃所在的海拔高度（图3）。我们利用中国区域地面气象要素驱动数据集（1979—2018）空间分辨率为0.1°，时间分辨率为3小时的气温数据［38-39］计算年气温变幅（日均最高气温与最低气温的差值，其数值的一半即Ta，图4）。由于古里雅冰芯及敦德冰芯的年代存在争议［40］，南极Vostok冰芯［41］所记录的末次冰盛期（Last Glacial Maximum，LGM）的降温幅度与青藏高原6～9℃的降温幅度较为一致［42］，其用于计算LGM以来平均寒冻风化强度的有效性在祁连山东段得到了验证［37］，且其温度记录的时间范围为420 ka至今，满足本研究对不同时间尺度温度变化的需求，因此我们选择南极Vostok冰芯作为长期温度变化的参考（图5）。

表1 11个国家气象信息中心气象站基本信息表Table 1 Information of 11 national meteorological stations

图3 11个气象站点海拔与年平均气温散点图Fig.3 Scatter plots of elevation against TMA

图4 祁连山Ta空间分布图Fig.4 The spatial distribution of Ta in the Qilian Mountains

3 结果与分析

3.1 风化碎屑边界的空间分布特征

风化碎屑的范围表现为自东向西逐渐增大的趋势（图6）。为了进一步验证风化碎屑边界提取的准确性，我们在祁连山东中西段各选择一个典型区域对比提取的风化碎屑边界与哨兵二号遥感影像［图6（b）～6（d）］。结果显示祁连山东段［图6（b）］与中段［图6（c）］的风化碎屑边界与遥感影像所显示的风化碎屑边界具有较好的一致性。而祁连山西段年降水量较低，植被覆盖较差（NDVI均值0.05），且多为斑块状分布，导致利用NDVI阈值提取风化碎屑边界较为困难，因此风化碎屑边界与遥感影像所显示的风化碎屑边界存在部分区域的不匹配［图6（d）］。

图6 祁连山风化碎屑边界分布图［（b）～（d）是黑色实线框区域内哨兵二号遥感影像及风化碎屑边界］Fig.6 The spatial distribution of the boundary of scree in the Qilian Mountains：（b）～（d）are Sentinel-2 imageries and the boundary of scree in the range of black rectangle

3.2 寒冻风化强度的时空变化特征

25～16 ka（图7蓝色实线）及16 ka以来（图7紫色实线）平均寒冻风化强度随海拔的变化曲线表现为单峰特征。由于从16 ka开始温度快速升高且在一段时间内维持在相对高温状态（图5虚线框），寒冻风化窗向高海拔区域移动，导致16 ka以来平均寒冻风化强度的峰值海拔高于25～16 ka，且两者寒冻风化峰值海拔差距较大，因此25 ka以来（包含一个冰期-间冰期旋回）平均寒冻风化强度随海拔的变化曲线表现为双峰特征（图7橙色实线）。十万年尺度平均寒冻风化强度随海拔的变化表现为单峰特征（图7），最大值在8～9℃·cm-1之间，且均在海拔3 900～4 700 m之间达到最大值（图7灰色阴影），不受随机选取的时间范围的影响。其与25 ka平均寒冻风化强度随海拔变化的双峰特征存在较大差异，主要表现为海拔5 000 m以上，25 ka平均寒冻风化强度依然较高。而在海拔5 000 m以下，25 ka平均寒冻风化强度随海拔的变化趋势与十万年尺度平均寒冻风化较为一致。

图5 420 ka以来Vostok冰芯温度记录（修改自Petit等［41］，紫色虚线框为25 ka以来温度变化）Fig.5 420 ka continuous temperature record based on the Vostok ice core，Antarctica：modified from Petit et al［41］.Purple dashed box is 25 ka continuous temperature record

图7 不同时间尺度平均寒冻风化强度随海拔的变化（灰色阴影代表十万年尺度平均寒冻风化强度峰值的海拔范围）Fig.7 The variation of time averaged Fci with elevation for different temporal scale（The grey shadows represent the peak elevation range of time averaged frost cracking intensity）

3.3 年气温变幅对寒冻风化强度空间变化的影响

为了探讨Ta对寒冻风化强度空间变化的影响，我们计算了不同Ta条件下总寒冻风化强度随海拔的变化［图8（a）］，并以25 ka平均寒冻风化强度为例，计算了不同Ta条件下平均寒冻风化强度随海拔的变化［图8（b）］。结果显示，Ta的增大可以使相同海拔条件下总寒冻风化强度增加，同时总寒冻风化强度峰值所在的海拔范围增大，进而导致平均寒冻风化强度增加，同时其峰值海拔范围增大。因此本文采取的差异化Ta参数是准确计算研究区寒冻风化强度的关键。

图8 不同Ta条件下总寒冻风化强度（a）及25 ka平均寒冻风化强度（b）随海拔的变化（灰色阴影代表Ta=15℃时寒冻风化强度峰值的海拔范围）Fig.8 The variation of total Fci（a）and time averaged Fci（b）with elevation with different Ta（The grey shadows represent the peak elevation range of Fci with Ta=15℃）

3.4 风化碎屑边界与寒冻风化强度的空间关系

由于万年尺度与十万年尺度平均寒冻风化强度的空间分布特征存在差异（图7），因此我们分别制作了风化碎屑边界与万年尺度［25 ka，图9（a）］和十万年尺度［100 ka，图9（b）］平均寒冻风化强度的空间分布图，并以风化碎屑与寒冻风化强度高值重叠面积占风化碎屑面积的百分比（表2）来定量探讨风化碎屑边界与不同时间尺度温度变化控制的平均寒冻风化强度的空间关系。结果表明，25 ka和100 ka平均寒冻风化强度的空间分布与风化碎屑边界的空间分布均存在较好的相关性，重叠面积占比大于65%，其中祁连山西段的重叠面积可达90%以上。由于Ta呈自东向西增大的趋势（图4），平均寒冻风化强度的峰值海拔范围增大［图8（b）］，导致祁连山西段寒冻风化强度高值的覆盖面积占比（大于65%）大于祁连山东段（小于50%）。但祁连山西段100 ka平均寒冻风化强度高值覆盖面积占比（76.49%）较25 ka平均寒冻风化强度高值覆盖面积占比（65.89%）的显著增加（超过10%）并没有导致重叠面积占比的大幅增加（仅增加3.4%），表明100 ka平均寒冻风化强度高值面积的增大并没有很好地提高其对风化碎屑边界的约束，而25 ka平均寒冻风化强度的高值以更小的覆盖面积约束了90%的风化碎屑范围，即25 ka平均寒冻风化强度的高值与风化碎屑边界的空间相关性优于100 ka。可能原因为祁连山西段海拔5 000 m以上区域所占比重较大（约1.5%），而25 ka平均寒冻风化强度在海拔5 000 m以上依然较高（图7橙色实线），可以较好地约束风化碎屑的边界。而对于祁连山中东段，100 ka平均寒冻风化强度高值覆盖面积占比（47.21%）较25 ka平均寒冻风化强度高值覆盖面积占比（27.78%）的显著增加（约20%）导致重叠面积占比的大幅增加（约20%），表明25 ka和100 ka平均寒冻风化强度的高值均可以很好地约束风化碎屑边界。可能原因为祁连山中段及东段海拔较低（99.8%的区域海拔低于5 000 m，图1），在海拔5 000 m以下，25 ka和100 ka平均寒冻风化强度随海拔的变化趋势较为一致（图7）。

表2 风化碎屑及不同时间尺度平均寒冻风化强度面积统计Table 2 Statistics of area of scree and different temporal of time averaged Fci

图9 祁连山25 ka（a）及100 ka（b）平均寒冻风化强度及风化碎屑边界空间分布图Fig.9 The spatial distribution of time averaged Fci since 25 ka（a）and 100 ka（b）and the boundary of scree in the Qilian Mountains

4 讨论

万年尺度平均寒冻风化强度随海拔的变化趋势在5 000 m以下与十万年尺度平均寒冻风化一致，表明寒冻风化在万年及以上时间尺度对该区域的作用是持续且稳定的。这种持续且稳定的寒冻风化为滑坡、泥石流等灾害性地貌过程提供充足的风化碎屑物质。因此在探讨寒冻风化对长时间（万年及以上）尺度地貌演化的累积影响时，时间尺度的选取对结果的影响可以忽略不计。我们推测的可能原因为：在冰期-间冰期旋回中，冰期持续时间长，间冰期持续时间短（图5），因此当时间尺度足够长时（例如万年及以上），相对短暂的间冰期在累积寒冻风化中的影响相对于冰期可以忽略不计。因此，在探讨寒冻风化在长期地貌演化中的累积影响时，如果研究区温度记录无法覆盖研究所需的较长时间范围（如百万年），可以采用万年或十万年尺度温度记录计算的平均寒冻风化强度作为替代。

25 ka平均寒冻风化强度与祁连山西段风化碎屑边界的空间相关性优于100 ka平均寒冻风化强度，表明当探讨寒冻风化对现代地貌单元空间分布特征的影响时，间冰期的快速升温过程显得尤为重要。温度快速升高且在一段时间内（千年—万年）维持在相对高温状态（图5虚线框），导致寒冻风化窗向高海拔区域移动，形成万年尺度平均寒冻风化强度随海拔变化曲线中的高海拔峰值（图7）。但由于快速升温为非连续性过程，寒冻风化对该峰值海拔范围的影响是阶段性的。在冰期大规模冰进的背景下，前一个间冰期所积累的碎屑物质无法保存，即寒冻风化的累积影响近乎清零。因此在探讨寒冻风化对现代地貌单元空间分布特征的影响时应采用最近一次冰期-间冰期旋回的温度变化数据（即末次冰盛期）。

在全球变暖的背景下，前人对由冻土退化、冰川退缩、冰湖扩张等引发的灾害性过程的研究已较为成熟［44-46］，但对受寒冻风化影响的灾害性地貌过程缺乏重视。本研究表明，温度升高使寒冻风化的作用区向高海拔区域移动，其提供的风化碎屑物质是滑坡、泥石流等灾害性地貌过程重要物质来源（图10），进而产生新的高风险灾害区域。寒冻风化不仅可以通过控制滑坡、泥石流等灾害性地貌过程的发生，成为滑坡-堰塞-溃决洪水这一灾害链［47］的源头，还可以加速与冰川退缩相关的灾害性地貌过程的发生［48］。本研究强调了全球变暖背景下寒冻风化对灾害性地貌过程的调控，为全球变暖背景下灾害性地貌过程的预测提供了新的思路，可以成为防灾减灾决策的重要参考。本研究是对崔之久先生寒冻风化如何影响生产建设并服务于社会经济建设的继承与发展。

图10 祁连山典型风化碎屑积累实体例证Fig.10 The accumulation of scree in the Qilian Mountains

5 结论

本文通过NDVI阈值、海拔筛选结合遥感影像目视解译，提取了祁连山北部的风化碎屑边界，通过对比分析风化碎屑边界及不同时间尺度平均寒冻风化强度的空间相关性，探讨干旱半干旱高寒山区风化碎屑空间分布的主要控制因素。获得主要结论如下：

（1）万年尺度平均寒冻风化强度随海拔的变化曲线表现为双峰特征，十万年尺度平均寒冻风化强度随海拔的变化曲线表现为单峰特征。海拔5 000 m以下万年尺度与十万年尺度平均寒冻风化强度随海拔的变化趋势一致。

（2）在祁连山北部万年尺度及十万年平均寒冻风化强度的高值与风化碎屑边界均存在较好的空间一致性，祁连山西段万年尺度平均寒冻风化强度的高值与风化碎屑边界的空间相关性优于十万年尺度。

（3）探讨寒冻风化对长时间（万年及以上）尺度地貌演化的累积影响时，时间尺度的选取对结果的影响可以忽略不计；探讨寒冻风化对现代地貌单元空间分布特征的累积影响时，应采用最近一次冰期-间冰期旋回的温度变化数据。

（4）在全球变暖背景下，应重视寒冻风化对灾害性地貌过程的调控。