利用形变时序分形特征识别高山冰川区滑坡

摘要

【目的】探究利用时序形变分形特征识别高山冰川区滑坡的方法并分析其适用性。

【方法】基于查莫利滑坡及其相邻冰川的形变时间序列描述其斜率（平均形变速率）及分形特征差异，利用聚类分析区分滑坡区域与冰川并进行影响因素分析。

【结果】与冰川相比，滑坡的形变时序具有较高的分形维数和较低的分形拟合优度。虽然滑坡与冰川在形变时序的斜率（平均形变速率）上也存在较大差异，但仅使用形变速率难以对滑坡进行聚类识别，准确率仅为61.70%；而使用形变时序的分形指标（包括分形维数和分形拟合优度）可将聚类分析的准确率显著提升至近84.00%。基于形变时序分形特征进行高山冰川区滑坡识别的适用性，根本原因在于滑坡和冰川在物质组成、影响因素和发展演化等方面存在差异。相较冰川，滑坡物质组成更复杂、更易受多种因素影响、形变时序的波动性更强。

【结论】利用形变时序分形特征能够成功识别高山冰川区滑坡，在全球变暖背景下，该方法预期可为高山冰川区的滑坡识别、进而为高山冰川区的防灾减灾提供一定的支撑。

高山冰川区滑坡识别形变时间序列分形维数分形拟合优度平均形变速率

引言

在全球变暖加剧背景下，冰川出现大幅退缩和消融等失常状态，导致冰川景观和冰川地貌快速改变和淡水资源减少等高山冰川区的地景改变，虽然大多源于冰川的流动消融，但雪崩和滑坡等灾害事件也是重要的肇因之一。相对冰川形变，高山冰川区的崩塌滑坡具有突发性强、破坏力大、移动范围广等特点，可对人民生命财产造成更大影响，应受更多重视。因此，在高山冰川区中识别出潜在的崩塌滑坡危险区，对评估和减轻灾害风险至关重要。

大型滑坡在发生完全破坏之前存在若干先兆，包括明显的地表形变。因此，地表形变测量成为识别潜在滑坡的常用途径，且这种识别主要是通过发现方向一致的地表形变来实现的。当前，常用的地表形变遥感反演方法为像素偏移跟踪（Pixel Offset Tracking，POT）和合成孔径雷达干涉测量（Interferometric Synthetic Aperture Radar，InSAR）两种。光学POT在形变量超过厘米级别时表现更好。

高山冰川区滑坡识别的挑战

冰川与滑坡形变相似

冰川与潜在滑坡均可导致方向一致的地表形变；特别是在冰川与滑坡具有相近运动速率的情况下，这对通过地表形变进行滑坡识别造成挑战。

环境复杂多变

高山冰川区环境恶劣，气候条件复杂，地形陡峭，增加了监测和识别的难度。

数据获取困难

高山冰川区的地理位置偏远，实地调查困难，主要依赖遥感技术获取数据，但遥感数据的质量和时效性受多种因素影响。

分形几何在地理学中的应用

分形几何是一种用于研究具有复杂、不规则和自相似性现象的非线性数学方法。众多自然现象既表现出空间形态与结构的复杂性，也表现出较为显著的自相似性，因而成为分形研究的重要研究对象。

地貌分析

在云南永胜古滑坡堰塞湖沉积物的研究中，不同的分形维数对应着不同的堰塞湖沉积条件。

景观生态

对上海市景观破碎化的定量分析中，分形拟合优度与分形维数的有机组合能够定量地刻画景观破碎化中的粒度变化特征。

滑坡研究

近年来，分形也被应用于分析滑坡形变时间序列。相关研究表明，滑坡形变时序的分形特征可反映滑坡形变在过程与机理上的差异。

研究目标与创新点

本文拟探讨基于形变时序分形特征的高山冰川区滑坡识别。主要研究内容包括：

基于分形维数和分形拟合优度等关键参数，提出一套基于形变时序分形特征的高山冰川区滑坡识别方案；
以查莫利高山冰川区滑坡为例，对基于形变时序分形特征的高山冰川区滑坡识别方案进行应用，验证其实用性；
初步探讨利用形变时序分形特征进行高山冰川区滑坡识别的适用原因，并分析进一步应用的前景和可行性。

本研究预期可为高山冰川区的高危崩塌滑坡识别提供新思路和新方法。

研究方法

技术流程

本文提出的利用时序分形特征识别高山冰川区滑坡的方法包括3个主要步骤：

图1 基于形变时序分形特征的高山冰川区滑坡识别技术流程

形变时序数据采集与处理

采用光学像素跟踪偏移方法实现地表形变反演。具体步骤如下：

输入不同时期的2幅图像组合作为图像对，前一幅图像为参考，后一幅图像为目标图像
输出3个图层，分别为NS（北/南）方向位移、EW（东/西）方向位移和信噪比
基于参考图像，使用滑动窗口检测目标图像的变化
通过稳定区校准初始NS向和EW向地表变形，以消除图像配准错误
通过合成NS向和EW向变形得出总变形
利用形变时序的线性拟合p值和单调性等指标去掉噪声，提取可靠的地表形变结果

形变时序分形特征计算方法

分形维数和分形拟合优度是分形几何中两个相关的概念。分形维数是一种用于量化分形结构复杂度的度量，较高的分形维数代表测量对象结构更加复杂且变动性更强。分形拟合优度则是指在将某一数据集或现象拟合为一个分形模型时，评估该拟合程度的好坏。

分形维数计算

形变时序各数据点可连接成为时间-形变量二维空间的一条曲线，其分形维数采用粗糙曲线的圆规曲线法计算，所得出的分形维数又称为Mandelbrot-维数。

利用固定长度为r₀的量尺沿形变时序曲线测量其长度，假如使用尺子的次数为N₀，其长度则为N₀×r₀。当r₀取不同值时，N₀也会有所不同。对于具有分形特征的曲线，N₀与r₀呈幂次关系：

N₀ = ar₀^-d

其中，a为常数，d即为该曲线的分形维数。两边取对数，则有：

lnN₀ = -d lnr₀ + lna = -d lnr₀ + c

利用所得的若干组(r₀，N₀)值，通过一元线性回归即可算出分形维数d值。

分形拟合优度计算

分形拟合优度使用线性回归模型的决定系数（Coefficient of Determination, R²）量化：

R² = 1 - R_ss/T_ss

其中，R_SS为残差平方和：

R_ss = Σ(y_i - ŷ_i)²

式中：y_i是实际观测值；ŷ_i是通过分形模型预测的值；T_ss = Σ(y_i - ȳ_i)²是总平方和，ȳ_i是观测值的均值。

计算步骤

形变曲线准备：将研究区各观测点的形变时序数据转化为时间-形变量二维空间的曲线。
测量尺度确定：根据曲线节点之间的距离确定量尺的长度，以对数等间距确定量尺长度。
曲线长度测量：采用以起始点为球心、以量尺长度为半径作球的方法，计算所需量尺的数量。
分形维数和分形拟合优度的计算：基于不同量尺大小及其对应的形变时序曲线测量计数，进行对数变换和线性回归。

基于形变时序分形特征的滑坡识别方案

根据高山冰川区滑坡和冰川形变时序在几何上的差异，本研究提出3个或可用以区分滑坡和冰川的指标，分别是形变时序曲线的斜率、分形维数和分形拟合优度。

斜率（平均形变速率）

时间-形变量数据的线性拟合线的斜率，可反映观测时段内的平均形变速率，指示了形变的快慢。

分形维数

描述了形变时序的复杂性，分形维数越大，表示形变过程越复杂，可能与形变对象的不稳定性质和复杂影响因素有关。

分形拟合优度

用于评估分形模型的可靠性和准确性，高拟合优度表示分形模型能够较好地描述形变过程的复杂性。

基于以上3个几何指标，本研究建立了3个维度的滑坡识别方案：

一维方案：仅采用单一指标（斜率、分形维数、分形拟合优度）对形变时序数据进行聚类分析。
二维方案：采用两个指标联合对数据进行聚类分析，包括斜率与分形维数组合、斜率与分形拟合优度组合、分形维数与分形拟合优度组合。
三维方案：将利用斜率、分形维数和分形拟合优度全部3个指标结合进行聚类分析。

聚类采用K-Means方法，拟划定的簇个数为2，也即滑坡和冰川2类。

实验区概况及数据来源

实验区概况

本研究选择印度北阿坎德邦的查莫利地区为研究区。其靠近喜马拉雅山西段的南达昆提峰，气候主要受印度夏季风控制，部分受西风带影响。年均降水量约为1770 mm，其中夏季（6—8月）的降水量占全年总量的一半以上。年平均气温约为2.1℃，夏季（6—8月）和冬季（12—次年2月）的平均气温分别为10.2℃和-6.3℃，山体冰川有明显活动迹象。

查莫利滑坡源区位于Ronti Peak（朗提峰）北侧的陡峭斜坡，海拔约5600 m，岩体表面常年被冰雪覆盖。滑坡坍塌体的体积约为27×10⁶ m³，坡度在35°~45°之间。从滑坡源头到山谷，海拔下降约1800 m，滑坡体在下落过程中积聚了巨大动量，并在山谷底部卷起大量冰碛。

图2 查莫利滑坡研究区位

地质背景

从区域构造和地层来看，查莫利滑坡位于喜马拉雅造山带次级构造单元内，地质构造复杂，山体走向与岩层走向基本一致，主要为NW-SE和NE-SW 2组方向；区内新构造活动频繁，岩石剪切作用显著，导致易发生山体滑坡。滑坡区出露的岩性主要为混辉岩片麻岩、夹带大理岩的变质岩。

形变时序数据来源

本文以2013年5月—2020年10月查莫利滑坡及其紧邻冰川的Landsat 8影像为基础，利用光学POT方法反演得出研究区的地表形变时间序列。

图3 本研究所使用Landsat8影像时间线

从google earth引擎下载了Landsat 8影像数据，选择了35张未被雾霾或云层覆盖的全色15 m分辨率的影像用于反演地面变形，时间跨度2704天，最短时间间隔16天，最长时间间隔192天。

结果及分析

形变时序特征

图4 查莫利滑坡及其相邻冰川形变时序的研究指标空间分布

研究指标	类别	最小值	最大值	中位数	平均值	标准偏差	变异系数
斜率(平均形变速率)/(cm/day)	滑坡	0.676	3.231	1.892	1.900	0.681	0.359
斜率(平均形变速率)/(cm/day)	冰川	0.427	6.066	3.134	3.194	1.316	0.412
分形维数	滑坡	1.015	1.146	1.043	1.053	0.029	0.028
分形维数	冰川	1.000	1.165	1.006	1.011	0.016	0.016
分形拟合优度	滑坡	0.997	1.000	0.999	0.999	0.001	0.001
分形拟合优度	冰川	0.997	1.000	1.000	1.000	0.000	0.000

斜率（平均形变速率）特征

从形变时序的斜率（平均形变速率）来看，查莫利滑坡及其相邻冰川都呈现由中心向边缘递减的同心圆分布。无论对于值域范围、中位数还是平均值，冰川的平均形变速率都显著大于查莫利滑坡。查莫利滑坡形变时序的斜率在东侧更大（1.86~3.68 cm/day）、西侧更小（0.43~2.29 cm/day），表明东侧变形速率较大，是滑坡的主要物源区。相邻冰川形变时序的最大斜率超过5.05 cm/day，在上部边缘区域也可达2.42 cm/day以上。

分形维数特征

从形变时序的分形维数看，滑坡分形维数总体大于冰川，统计表明其中位数和平均值（1.043和1.053）均高于冰川（1.006和1.011）。同时，滑坡形变时序分形维数的标准偏差和变异系数（0.029和0.028）也高于冰川（0.016和0.016），表明滑坡形变时序分形维数的数值更高、分布更为分散、变化幅度更大。滑坡分散分布的高分形维数值表明，其形变时序具有更高的复杂性和长期依赖性，且滑坡内部差异较大。

分形拟合优度特征

从形变时序的分形拟合优度来看，滑坡和冰川的最小值和最大值都非常接近（0.997和1.000），值域范围都是0.003，表明两组数据的数值都分布在一个狭窄的、接近于1的范围内，分形模型都具有较高的拟合度。滑坡形变时序分形拟合优度的峰度为4.115、偏度为-1.271，而冰川的峰度为64.610、偏度为-6.303，说明两者的数据分布都向左偏斜（即更多的数据点集中在平均值的右侧）。冰川的峰度和偏度绝对值都要远远大于滑坡，表明冰川数据分布非常尖锐，大多数数据点集中在平均值附近，数据分布更加不对称。

聚类识别结果与分析

图6 聚类识别方案结果的空间分布

一维聚类识别方案

当仅使用斜率（平均形变速率）作为单一变量进行聚类时，滑坡的识别区域覆盖率是最大的；并且，在真实滑坡区域中被误识为冰川的概率很小，仅在中心区域存在小部分。但是，在此方案的结果中，存在大面积的冰川被误识为滑坡；也即，在冰川外围区域，滑坡的面积被大幅度放大，整体识别准确率不高。

在仅使用分形维数以及分形拟合优度的情况下，冰川的绝大部分区域均能很好的识别；尤其是在使用分形拟合优度的方案中，几乎可以将整个冰川区完整的识别出来。但是，在此方案的结果中，滑坡的东侧物源区会与冰川产生一部分混淆现象，这可能与查莫利滑坡的发育演化过程有一定的联系。

二维聚类识别方案

斜率（平均形变速率）与分形拟合优度以及分形维数与分形拟合优度两种指标组合的聚类准确度要明显优于斜率（平均形变速率）与分形维数组合。斜率与分形维数组合识别方案，相较仅使用斜率（平均形变速率）的方案，识别出了更加清晰的滑坡冰川界限；但是，冰川区域误识为滑坡的区域集中到了冰川的西南部，也就是滑坡区域分成了两个部分，这对滑坡的识别是非常不利的。

斜率（平均形变速率）与分形拟合优度以及分形维数与分形拟合优度2个方案，识别效果以及空间分布差别不大；与仅使用分形维数以及仅使用分形拟合优度的结果相比，效果有小幅度提升，可作为识别滑坡区域的最优方案。

三维聚类识别方案

同时使用斜率（平均形变速率）、分形维数以及分形拟合优度的方案，准确度与仅使用分形维数与分形拟合优度的方案相比几乎没有提升。

适用性验证

维度	指标	准确率
一维	斜率（平均形变速率）	61.70%
	分形维数	83.76%
	分形拟合优度	82.52%
二维	斜率+分形维数	78.34%
	斜率+分形拟合优度	83.92%
	分形维数+分形拟合优度	83.98%
三维	斜率+分形维数+分形拟合优度	84.01%

从结果可以看出，对于一维聚类识别方案，当仅使用斜率（平均形变速率）作为单一变量进行聚类时，准确率仅为61.70%，相对较低。这是因为，滑坡与冰川的形变速率在数值分布上有较大重叠，单一使用形变速率变量难以有效区分滑坡和冰川。

相比之下，单一使用分形维数的聚类准确率为83.76%，表明分形维数是一个更强的分类指标，能更好地区分滑坡和冰川。单一使用分形拟合优度的聚类准确率为82.52%，与分形维数类似，也表现出较强的分类能力。

在二维和三维聚类方案中，分形维数与分形拟合优度组合的准确率达到83.98%，三维方案的准确率为84.01%，相较于一维方案有所提升，但提升幅度有限。

讨论

滑坡与冰川形变时序特征差异的成因

物质组成差异

滑坡由岩石、土壤等多种物质组成，结构复杂多变；而冰川主要由冰组成，物质相对单一。这导致滑坡形变时序的分形维数更高，表现出更复杂的形变模式。

影响因素差异

滑坡受地质结构、降雨、地震等多种因素影响；冰川主要受温度变化影响，影响因子较为单一。多种因素的综合作用使滑坡形变时序更为复杂。

发展演化差异

滑坡形变通常表现为非线性、阶段性变化过程；冰川流动则相对稳定、连续。这导致滑坡形变时序的分形特征更为显著。

滑坡形变时序分形维数高于冰川，表明滑坡形变过程更为复杂，这与滑坡的物质组成、影响因素和发展演化特性密切相关。滑坡由多种地质材料组成，受多种因素影响，形变过程呈现非线性特征；而冰川物质相对单一，主要受温度影响，形变过程较为稳定。

滑坡形变时序分形拟合优度略低于冰川，说明滑坡形变时序的分形特性不如冰川那样"纯粹"，这可能是因为滑坡受多种因素影响，形变过程更为复杂多变，导致其分形特性受到一定程度的"污染"。

方法的适用性与局限性

适用性

该方法在冰川与滑坡形变速率相近的情况下，仍能有效区分两者，解决了传统方法的局限性
分形特征能够反映形变过程的内在机制差异，具有物理意义
该方法不依赖于特定的形变速率阈值，具有较好的通用性
可与现有的基于形变速率的方法结合，提高识别准确率

局限性

需要较长时间序列的形变数据，对数据的时间跨度和质量要求较高
计算分形特征的过程相对复杂，增加了处理的复杂度
在某些特殊地质条件下，滑坡与冰川的分形特征可能存在重叠，影响识别效果
目前仅在查莫利滑坡案例中验证，需要更多案例来检验方法的普适性

结论

本研究探讨了基于形变时序分形特征识别高山冰川区滑坡的方法，通过对查莫利滑坡及其相邻冰川的形变时序特征分析，得出以下主要结论：

1

滑坡与冰川在形变时序的分形特征上存在显著差异。与冰川相比，滑坡的形变时序具有更高的分形维数和稍低的分形拟合优度，反映了滑坡形变过程的复杂性和多样性。

2

基于形变时序分形特征的聚类分析能够有效识别高山冰川区的滑坡。单独使用分形维数或分形拟合优度进行聚类的准确率均超过82%，远高于仅使用形变速率的准确率（61.70%）。

3

分形维数与分形拟合优度组合的二维聚类方案效果最佳，准确率达到83.98%，几乎与三维聚类方案（84.01%）相当，表明分形特征是区分滑坡与冰川的有效指标。

4

滑坡与冰川形变时序特征差异的根本原因在于两者在物质组成、影响因素和发展演化等方面的不同。相较于冰川，滑坡物质组成更为复杂，受多种因素影响，形变过程更为动态多变。

本研究提出的基于形变时序分形特征识别高山冰川区滑坡的方法，为高山冰川区滑坡识别提供了新的思路和技术手段。在全球变暖背景下，该方法有望为高山冰川区的滑坡早期识别和防灾减灾工作提供支持。

未来研究可以进一步探索分形特征与滑坡发展阶段的关系，建立更加完善的滑坡预警模型；同时，可以将该方法应用于更多高山冰川区案例，验证其普适性，并结合机器学习等先进技术，进一步提高识别精度。