北京大学学报(自然科学版) 第60卷 第2期 2024年3月

Acta Scientiarum Naturalium Universitatis Pekinensis, Vol. 60, No. 2 (Mar. 2024)

doi: 10.13209/j.0479-8023.2024.003

国家重点研发计划(2017YFB0502703)资助

收稿日期: 2023–04–11;

修回日期: 2023–05–23

利用结合土地覆盖类型的自适应DS-InSAR方法监测矿区地表形变

张志亮1,2 曾琪明1,2,† 杨立功1,2

1.北京大学遥感与地理信息系统研究所, 北京大学地球与空间科学学院, 北京 100871; 2.空间信息集成与3S工程应用北京市重点实验室, 北京 100871; †通信作者, E-mail: qmzeng@pku.edu.cn

摘要 针对传统时序 InSAR 技术在测量点选取过程中多采用全局单一阈值, 在受时空失相干影响严重和短时间内形变梯度较大的矿区容易出现测量点稀疏且空间采样不足, 因而无法获取矿区地面沉降完整信息的问题, 以覆盖大同矿区的 22 景 ALOS-1 数据为例, 利用结合土地覆盖类型的自适应 DS-InSAR 方法(ADSI-CLC)获取该区域长时间序列地表形变信息。ADSI-CLC 方法形变测量结果的时空分布模式与 StaMPS-SBAS 方法具有相似性, 且与高分辨率光学遥感影像上的采矿设施保持良好的空间相关性。ADSI-CLC 方法可以显著地提高测量点的数量和空间分布密度, 在研究区识别出的 DS 点数量约是 StaMPS-SBAS 方法的 4 倍。两种方法的测量结果在无形变和小形变区域吻合程度较高; 在大形变区域, StaMPS-SBAS 方法无法有效地获取形变结果, 而 ADSI-CLC 方法由于测量点数量的增加, 反演得到的形变测量结果基本上符合由采矿活动造成的漏斗状的空间分布特征, 间接地验证了该方法在矿区地表形变监测中的可靠性和有效性, 证明 ADSI-CLC 方法能够提供更详细的时空形变细节信息, 更好地服务于矿区地表稳定性的监测和预警。

关键词 ADSI-CLC 方法; 分布式散射体; 矿区形变监测; 大同矿区

长期高强度的资源开采, 容易破坏地层原有的应力状态, 导致矿区生态环境遭到破坏, 引发地表塌陷、地下水倒灌及山体滑坡等各类地质灾害, 进而造成房屋倒塌等现象, 对经济社会的可持续发展造成极大的影响[1]。因此, 开展矿区地面沉降监测具有非常重要的现实意义。目前, 传统的监测技术(如水准测量、全站仪和 GNSS 等)已在矿山安全监测领域发挥重要作用, 虽然监测精度高, 但存在耗时长、效率低和时空分辨率低等诸多问题[2], 难以全面地掌握矿区形变特征。

近年来, 由 D-InSAR (differential interferometric synthetic aperture radar)技术发展而来的时序 InSAR测量技术(multi-temporal InSAR, MT-InSAR)凭借其高精度、大范围、高时空分辨率和长时间序列等优势, 广泛地应用于矿区地表形变监测中[3–9]。MT-InSAR 技术基于长时间序列 SAR 影像生成的干涉图组合, 仅对受时空失相干影响较小的测量点进行解算, 根据不同相位成分的时空变化特征建模, 估计并分离形变相位和各类误差相位。然而, 矿区多位于自然地表区域, 人工建筑物稀疏, 地表通常被耕地和植被等地物覆盖, 在长时间序列上容易受时空失相干影响[4–6]

传统的时序 InSAR 技术[10–12]在测量点选取过程中多采用全局单一的选取阈值, 容易造成矿区测量点分布不均匀, 出现大片区域测量点稀疏或无测量点的极端情况, 进而导致解算结果不可靠, 甚至无解, 难以对矿区沉降进行有效的测量和分析。降低 DS(distributed scatterer, 分布式散射体)点的选取阈值, 虽然可在一定程度上增加其数量, 但增加的点不一定分布在感兴趣区域(region of interest, ROI); 若要保证感兴趣区域测量点的数量, 则可能引入大量低质量像元, 影响时序 InSAR 解算结果的精度, 并增加解算时间。

针对时序 InSAR 技术在测量点选取过程中存在的问题, 本研究组在前期工作中提出一种结合土地覆盖类型的自适应 DS-InSAR 方法(adaptive distribu-ted scatterer InSAR combined with land cover, ADSI-CLC)[13], 实现 DS 自适应选取, 并将其用于四川省丹巴县甲居滑坡长时间序列的地表形变监测, 有针对性地提高了滑坡区域测量点的数量和空间分布密度。与 StaMPS-SBAS 方法和传统全局单一阈值法相比, 该方法可以提供滑坡区域更丰富的时空形变细节信息, 但是否可以解决矿区形变监测中存在的测量点稀疏问题有待研究。此外, 在估计时序地表形变时, 传统的解算方法均基于单个测量点的多次观测来构建解算方程, 通过最小二乘法, 直接求解时序地表形变, 未考虑空间维度上测量点的质量差异对模型解算的影响。实际处理过程中, 研究区测量点质量存在差异, 尤其是采用 DS 自适应选取方法后, 势必为了提高测量点数量而引入一些质量较差的点, 这时就需要考虑形变解算过程中由测量点质量差异带来的影响[14]

鉴于上述背景, 本文选取大同矿区为研究对象, 以 2007—2011 年的 22 景 L 波段 ALOS-1 PALSAR-1影像为数据源, 利用 ADSI-CLC 方法获取该区域研究时段内的视线(line of sight, LOS)向地表形变, 并选取典型形变区域开展分级优化解算, 在此基础上分析矿区在监测时段内的地表形变时空分布特征, 以期为采空区塌陷灾害预警及治理提供技术和数据支撑。

1 技术方法

1.1 ADSI-CLC方法

本研究采用的 ADSI-CLC 方法核心步骤如下: 第一步, 结合土地覆盖类型的 DS 自适应选取; 第二步, 基于相干性加权的空间自适应滤波算法的干涉相位优化[15]。其中, DS 自适应选取是 ADSI-CLC方法最重要的步骤。该方法的技术流程如图 1 所示, 技术细节参阅文献[13,16]。

DS 自适应选取算法包括同质像元(statistically homogeneous pixels, SHP)识别(即 DSC(distributed scatterer candidates)选取)和 DS 选取两部分, 与经典算法(如 SqueeSARTM[10]和 ADSI[17])的主要差别如下。

width=453.6,height=391.2

第一行图中五角星代表中心像元, 曲线代表土地覆盖类型的边界, 黑色圆点代表窗口内与中心像元同质的像元

图1 DS自适应选取流程

Fig. 1 DS adaptive selection process

1)在同质像元识别部分, 传统方法仅利用统计检验(如 K-S 检验[10]和 A-D 检验[18]等)来识别同质像元, 不仅容易导致非同类像元参与检验, 增加计算负担, 而且会增加非同类像元被误判为同质像元的概率。ADSI-CLC 方法是在传统方法的基础上, 引入具有地理学属性的土地覆盖类型图作为同质像元识别的约束条件, 当且仅当二者同类, 且土地覆盖类型为非水体时, 才会开展后续统计检验, 即统计检验仅在相同土地覆盖类型的边界内进行, 有效地降低非同类像元被误选为同质像元的概率, 提高同质像元的识别效率, 纯化同质像元集合和 DSC 选取集合。通过引入土地覆盖类型图, 实现同质像元识别过程中散射机理、统计检验和地表覆盖类型三者的协调统一。

2)在 DS 选取部分, ADSI-CLC 方法引入土地覆盖类型图作为研究区分区依据, 在不同的土地覆盖类型区域, 自适应地确定不同的相干性选取阈值(如式(1)所示), 有针对性地提高 DS 的数量和空间分布密度, 是一种考虑像元自然属性的自适应选取方法。

width=158.65,height=21.65(1)

其中, corthreshold(I)表示第 I 种土地类型下的 DS 相干性选取阈值, min(·)表示取括号内两个数值中较小者,width=29.25,height=17.55代表研究区所有 DSC 时序平均相干性的空间均值,width=41.55,height=17.55表示各土地覆盖类型区域所有DSC 时序平均相干性的空间均值。

1.2 考虑DS质量差异以及空间相关性的解算方法

为了提高感兴趣区域 DS 点的数量, ADSI-CLC方法势必引入一些低质量的点, 从而影响矿区地表形变监测。为削弱这种影响, 本文利用考虑 DS 质量差异和空间相关性的解算方法, 对典型形变区域进行分级优化解算, 对比优化解算前后形变速率的差异。具体地, 通过引入地表应力应变模型来建立邻近点间的空间相关性约束[19–20]; 利用时间相干性作为质量分级指标, 将 DS 点划分为高质量一级点和低质量二级点, 一级点的初始解算结果可靠, 无需进行优化, 可以作为二级点优化解算的初值和约束条件, 二级点的初始解算结果可靠性较差, 亟需优化。在分级优化解算的求解过程中, 采用附加约束条件的间接平差策略。具体流程如下。

1)初始解算。

2)基于地表应力应变模型的 InSAR 干涉相位优化。

3)一级分区: 像空间规则分块。

4)测量点分级: 分为一级点和二级点。分级过程中, 需要兼顾测量点质量以及空间分布的均匀性和连通性。

5)局部Delaunay三角网的构建。

6)优化解算: ①以弧段上的干涉相位差为观测量, 求解各弧段上的累积地表形变差; ②以空间各弧段上的累积地表形变差为观测量, 以一级点的初始解算结果为约束条件, 利用附加约束条件的间接平差求解各二级点的约束解算结果; ③将二级点初始解算结果与约束解算结果进行加权融合, 得到各二级点的优化解算结果。

根据一级点初始解算结果, 可以建立如下约束方程:

width=67.3,height=15.8(2)

其中, 矩阵 C 为系数矩阵, 标记(Q–1)个一级点的位置, 大小为(Q–1)×(H–1); 矩阵 X为待求矩阵, 大小为(H–1)×1; Mx为(Q–1)个一级点的初始解算结果, 大小为(Q–1)×1, 为已知量。利用附加约束条件的间接平差策略, 可以得到待求量 X 的平差结果:

width=204.9,height=15.8(3)

其中,

width=81.95,height=15.8

width=71.4,height=15.8

width=71.4,height=15.8

其中, Warc 是大小为 G×G 的对角权重矩阵, 其中每个元素值刻画一级子区域中各条弧段的质量, 利用测量点的相干性对不同弧段进行定权; Lref 为弧段上的累积地表形变差, 为观测量; Uref 为系数矩阵, 每一行代表一条弧段, 每一列代表一个测量点。

2 研究区概况及数据

大同市位于山西省最北端, 与河北、内蒙接壤, 是中国最大的煤炭能源基地之一, 多年的煤炭开采引发地面沉降、生态破坏以及空气污染等诸多问题[21], 由地面沉降引起的滑坡、地表塌陷、地裂缝等地质灾害严重地影响着矿区周围村庄道路及地表建筑物的安全。

如图 2(a)所示, 研究区位于大同市西南方向, 地形起伏较大, 地势西北高东南低, 大多为高山和起伏较大的丘陵地貌, 存在大面积的塌陷区。研究区地理坐标为东经 112.89°—113.25°, 北纬 39.96°—40.24°。该区气候干燥寒冷, 属大陆性气候, 早晚温度低, 降雨集中在每年的 7—9 月。如图 2(b)所示, 研究区地表覆盖以草地和耕地为主, 存在部分林地和居民用地。

本研究利用的 SAR 数据为 2007 年 2 月至 2011年 1 月覆盖山西省大同矿区的 22 景 ALOS-1 PA-LSAR-1 数据, 具体的参数信息和获取时间如表 1 所示。外部 DEM 数据为 30m 分辨率的 SRTM1 DEM数据, 土地覆盖类型图为 Gong 等[22]提供的 10m 分辨率的土地覆盖类型数据集。

采用时空基线阈值方法[23]选取差分干涉对, 将相干性阈值(γ)、时间基线阈值(Bt)和空间基线阈值(B)分别设置为 0.1, 300 天和 2500m。共选中 49 幅干涉图, 对应的时空基线如图 3 所示。

3 结果与讨论

3.1 ADSI-CLC方法的形变测量结果

同质像元识别过程中, 设置同质像元检测窗口的大小为 25×9 (方位向×距离向), 同质像元检验的置信度水平为 0.5, DSC 阈值选取为 25, 得到的同质像元识别结果和 DSC 选取结果如图 4(a)所示。之后, 根据土地覆盖类型图, 将研究区划分为不同的土地覆盖类型区域。不同土地覆盖类型区域自适应地进行计算, 并采取不同的 DS 阈值, 实现 DS 的自适应选取, 各 DS 的时序平均相干性计算结果如图4(b)所示。

width=476.25,height=223.9

图2 研究区地形(a)及土地覆盖类型(b)

Fig. 2 Topographic map (a) and land cover types (b) of the study area

表1 ALOS-1 PALSAR-1数据的详细参数信息

Table 1 Detail acquisition information of ALOS-1 PALSAR-1 data

参数项 参数值 轨道方向升轨 轨道号454 波长23.6 cm 波段L波段 重访周期46天 距离向分辨率9.4 m 方位向分辨率3.1 m 极化方式HH/HV 成像模式FBS/FBD 影像数量22 时间跨度2007年2月至2011年1月

基于上述 DS 选取和干涉相位优化结果, 根据时序 InSAR 处理流程, 获取研究区 LOS 向地表形变速率, 如图 5(a1)所示, 图5(a2)为对应的直方图。为了对比, 使用国际流行的开源软件 StaMPS 中的SBAS 方法[23]进行地表形变测量, 将得到的结果作为参考标准, 用于验证 ADSI-CLC 方法在矿区地表形变监测中的可靠性。为了客观地比较 ADSI-CLC和 SBAS 两种方法, 在实际处理过程中, 除 DS 的选取和干涉相位的优化两个环节外, 其余环节均采用完全相同的参数和误差去除策略。StaMPS-SBAS方法的形变测量结果如图 5(b1)所示, 图 5(b2)为对应的直方图。可见, 两种方法反演得到的地表形变的空间分布模式和量级基本上一致, 均反映出该区域由于采煤导致的地面沉降, 最大形变速率达到−180mm/a。

width=221.15,height=155.85

每一个蓝色圆点代表一景 SAR 影像, 每一条黑色线代表一幅干涉图

图3 干涉图的时空基线网络

Fig. 3 Spatial-temporal baseline network of interferometric pairs

width=470.5,height=161.5

图4 DSC (a)和DS (b)的分布

Fig. 4 Distribution of DSC (a) and DS (b)

width=413.85,height=260.75

(a1)和(a2) ADSI-CLC方法; (b1)和(b2) StaMPS-SBAS方法。红色五角星代表空间参考点, 黑色空心圆圈代表交叉比对点, L1~L5为形变区域编号

图5 大同矿区LOS向地表形变速率及其直方图

Fig. 5 LOS surface deformation velocity and histograms in the Datong mining area

对研究区的 DS 数量及其空间分布密度进行统计, 结果如表 2 所示。与经典的 StaMPS-SBAS 方法相比, ADSI-CLC 方法可以保留更多的测量点。对本研究区而言, ADSI-CLC 方法选取得到的 DS 数量和空间分布密度大约是 StaMPS-SBAS 方法的 4 倍, 数量提高 3 倍, 而时序平均相干性均值仅损失11.7%, 以较小的相干性损失换来 DS 数量的大幅度提高。特别是在形变区域 L1~L4, DS 数量的提高尤为明显。这是一种有针对性地提高 DS 数量的策略, 有助于对矿区地面沉降区域展开分析。由此可见, ADSI-CLC 方法可以提供更丰富的时空形变细节信息, 更好地服务于地面沉降区域的识别和定位。我们基于 ADSI-CLC 的测量结果, 在研究区识别出至少 20 个明显的沉降区域。

表2 DS的数量及其空间分布密度

Table 2 Number and spatial distribution density of DS

方法 时序平均相干性均值DS数量DS空间分布密度/(DS·km−2) StaMPS-SBAS0.4958788307825 ADSI-CLC0.437831486493296

结合图 2(b)可以看出, ADSI-CLC 方法主要提高了林地、耕地和草地区域的 DS 分布数量, 分别是StaMPS-SBAS 方法的 5.6, 4.7 和 3.9 倍, 林地区域数量的提高尤为明显。特别地, 水体在 SAR 影像上表现出的相干性很差, 属于无效的测量点, 如果这些测量点参与解算, 会在增加解算时间的同时降低解算结果的精度。ADSI-CLC 方法以土地覆盖类型图作为约束条件, 直接剔除土地覆盖类型为水体的像元, 避免了无效的测量点被选中。

3.2 ADSI-CLC测量结果可靠性分析

由于未能获取该区域研究时间段内的 GPS 或水准数据, 本文以 StaMPS-SBAS 方法的测量结果为参考, 进行 ADSI-CLC 方法的可靠性验证。

选取 18810 个地理坐标完全一致的同名测量点, 对 StaMPS-SBAS 方法和 ADSI-CLC 方法得到的形变速率进行对比, 计算得到形变速率误差的均值、平均绝对误差、标准差和均方根误差值分别为 3.3, 10.5, 15.6 和 15.9mm/a, 皮尔逊相关系数为 0.5。在无形变和小形变区域, 两种方法的测量结果吻合程度较高; 在大形变区域, StaMPS-SBAS 方法无法获取有效的形变结果, 而 ADSI-CLC 方法得到的形变测量结果符合由采矿活动造成的漏斗状的空间分布特征[3], 如图 6 所示。

在区域 L1~L4 中的无形变、小形变和大形变区域, 随机选取 12 个交叉比对点(具体位置如图 5 所示), 比较两种方法获取的 LOS 向累积地表形变, 结果如图 7 所示。可以看出, ADSI-CLC 和 StaMPS-SBAS 方法获取的长时间序列地表形变结果在所有交叉比对点上均表现出完全一致的形变趋势; 无形变和小形变的交叉比对点上形变量级一致, 而大形变的交叉比对点上形变量级差异较大。我们推测是由于 StaMPS-SBAS 方法在大形变区域无法获取足够数量的测量点而发生解缠错误, 导致解算结果不准确。

上述实验结果间接地验证了 ADSI-CLC 方法在矿区地表形变监测中的可靠性。接下来, 我们对其在典型形变区域(4 个明显的沉降区域 L1~L4 和一个无形变区域 L5)中的表现进行分析。

3.3 典型形变区域监测结果分析

对比区域 L1~L5 的形变速率场(图 8 中第 4~5行)可知, ADSI-CLC 方法显著地提高了这些区域DS点的数量, 因此可以提供更详细的形变细节信息。区域 L1 在研究时段地表形变范围最大, 且形变程度最高; 区域 L5 为大同市城区, 未发生明显的地表形变。

由图 8 中区域 L1 的原始缠绕干涉图可知, 该区域的南部和东部在研究时段均出现明显的形变条纹, 利用 StaMPS-SBAS 方法在该区域识别出的 DS点较少, 空间采样严重不足, 易发生解缠错误, 因此无法准确地探测沉降发生的位置及沉降分布的范围; 相反, ADSI-CLC 方法在该区域识别出的 DS 点数量显著增加, 既提高了相位解缠过程中的空间采样, 也提供了更丰富的形变细节信息, 进而可以有效地确定采矿引起地表形变的影响范围。此外, 从ADSI-CLC 方法的测量结果可以看出, 该区域有两处明显的沉降中心, 分布于南侧和东侧, 与原始干涉图中干涉条纹出现的位置吻合, 区域内受影响较严重的村庄是夏庄村和红墙村。

对区域 L2 而言, 在研究时段, 利用 StaMPS-SBAS 方法仅能获取少量漏斗状沉降盆地边缘的形变信息, 稀疏且分布不均匀的 DS 点难以准确地反映地表形变的真实情况; ADSI-CLC 方法则可以有效地探测出地表沉降发生的大体位置和影响范围, 与该区域出现明显干涉条纹的位置吻合, 地表沉降中心区域的形变速率约为−150mm/a, 边缘区域的形变速率约为−60mm/a。针对区域 L3, 两种方法均未能有效地探测出地面沉降区域的大体位置和影响范围, 但 ADSI-CLC 方法在边缘区域保留了更多的测量点, 可以提供局部的地表形变信息。对于区域L4, ADSI-CLC 方法可以识别出分布在其中部、南部和东南部的 3 处明显的沉降中心, 这些沉降中心的位置与原始缠绕干涉图中干涉条纹的位置基本上吻合; StaMPS-SBAS 方法由于测量点稀疏, 未能准确地识别区域 L4 发生地表形变的范围。

3.4 联合光学遥感影像的地表形变分析

以形变范围最大的区域 L1 为例, 结合光学遥感影像, 分析地表形变情况。大同矿区多采用井下开采模式, 开采工作面多分布于采矿设施附近。图9(a)显示, 区域 L1 存在多处采矿设施。煤炭资源开采会造成采矿设施附近区域发生地面沉降, 严重时会形成塌陷坑。如图 9(b)所示, InSAR 测量结果识别的地面沉降区域均位于采矿设施附近, 这一现象与由井下资源开采导致的可能发生地面沉降的区域相吻合。截至光学遥感影像的获取时间(2010 年 9月), 采矿设施附近地表尚未形成明显的塌陷坑, 如果仅依据光学遥感影像, 无法准确地识别地面沉降区域。为详细地了解区域 L1 的形变信息, 掌握形变的整体发展情况, 图 10 给出该区域 2007—2011年的累积地表形变。可以看出, 在研究时段早期, 沉降中心位于该区域南侧; 随着时间推移, 在东侧出现新的沉降中心, 且影响范围比南侧大。

width=476.25,height=167.25

图6 漏斗状地面沉降中心

Fig. 6 Funnel-shaped ground subsidence centers

width=470.5,height=277.8

第1~4列依次为区域L1~L4, 第1~3行依次为无形变点、小形变点和大形变点

图7 交叉比对点处的累积地表形变

Fig. 7 Cumulative ground deformations at cross-comparison points

width=416.6,height=368.5

第 1~5 列依次为区域 L1~L5; 第 1~3 行依次为干涉对 20070213 和 20070701, 20081003 和 20090103 以及 20100709和20110109; 第 4 行为 StaMPS-SBAS 方法得到的形变速率, 第 5 行为 ADSI-CLC 方法得到的形变速率

图8 区域L1~L5缠绕相位及形变速率

Fig. 8 Wrapped phase and deformation velocity of regions L1–L5

width=430.9,height=181.4

(a)光学遥感影像, 红色曲线圈出采矿设施分布范围; (b) LOS向地表形变速率

图9 矿区地表形变分析结果

Fig. 9 Results of ground deformations analysis in mining area

width=476.25,height=317.5

图10 区域L1的长时间序列LOS向形变

Fig. 10 Time series LOS deformation maps of region L1

width=467.75,height=246.6

第1~5列依次为区域L1~L5

图11 区域L1~L5地表形变速率优化解算结果(第一行)、优化解算前后形变速率差异(第二行)和相干性(第三行)

Fig. 11 Optimization results of surface deformation rate in regions L1–L5 (first row), deformation rate difference before and after optimization (second row) and coherence (third row)

3.5 典型形变区域优化解算结果对比

典型形变区域 L1~L5 地表形变速率的优化解算结果如图 11 所示。整体来看, 该研究区中典型形变区域优化解算前后 LOS 向地表形变速率的空间分布模式和量级基本一致。具体地, 相比于初始解算结果, 优化后的 LOS 向形变速率在高质量 DS 点处未出现显著差异, 在低质量DS点处差异较大, 最大达到 40mm/a。

在区域 L1~L5, 与 ADSI-CLC 方法的整体解算策略相比, 优化解算方法虽然可以对部分低相干测量点的解算结果进行优化, 且相干性越差, 优化效果越明显, 但未能改变 ADSI-CLC 方法识别的沉降区域范围。因此, 为了提高 DS 点数量而引入的一些低质量点不会对 ADSI-CLC 方法的形变解算结果产生太大的影响, 其识别的矿区地面沉降范围基本上是准确的。在实际地表形变监测过程中, 仅采用ADSI-CLC 算法就可以实现矿区地面沉降范围的准确识别和定位。若需要某些点位的形变数值或需要基于 InSAR 测量结果开展地层参数反演, 则可以在ADSI-CLC 测量结果的基础上, 针对所关注的区域采用分级优化解算方法开展进一步的处理。

4 结论与展望

本研究以 2007—2011 年覆盖大同矿区的 22 景ALOS-1 PALSAR-1 影像为例, 利用 ADSI-CLC 方法获取该区域研究时段内的 LOS 向地表形变, 并与StaMPS-SBAS 方法的结果进行比较, 得到如下主要结论。

1)在研究区, ADSI-CLC 方法识别的 DS 数量以及空间分布密度是 StaMPS-SBAS 方法的约 4 倍。特别地, 针对复杂形变区域, DS 数量同样可以得到显著提高。因此, ADSI-CLC 方法解决了 StaMPS-SBAS 方法在矿区地表形变监测中存在的测量点稀疏或分布不均匀的问题, 可以提供更详细的矿区时空形变细节信息。

2)ADSI-CLC 方法与 StaMPS-SBAS 方法获取的地表形变时空分布模式具有相似性, 二者都识别出研究区内存在多处由采矿引起的地面沉降区。两种方法获取的无形变和小形变区域吻合程度较高; 大形变区域有明显的差异, StaMPS-SBAS 方法无法有效地获取形变信息, 而 ADSI-CLC 方法获取的形变结果基本上符合由采矿活动造成的漏斗状沉降区的空间分布特征, 最大形变速率达到−180mm/a, 由此验证了 ADSI-CLC 方法在矿区地表形变监测中的可靠性。StaMPS-SBAS 方法可以反映出开采沉降的影响范围和沉降区的空间分布特征, 为矿区大范围、精细化的形变调查和分析提供技术支撑。

3)ADSI-CLC 方法对矿区最大沉降量的获取能力有限, 需联合其他方法进行形变测量, 以期获取更准确的矿区地表形变信息。

参考文献

[1] 刘童谣. 基于SBAS-DInSAR技术的煤矿地表沉降监测研究[D]. 太原: 太原理工大学, 2019: 1–2

[2] 董玉森, Ge Linlin, Chang Hsingchun, 等. 基于差分雷达干涉测量的矿区地面沉降监测研究. 武汉大学学报(信息科学版), 2007, 32(10): 888–891

[3] 朱建军, 邢学敏, 胡俊, 等. 利用InSAR技术监测矿区地表形变. 中国有色金属学报, 2011, 21(10): 2564–2576

[4] 王波, 谭志祥, 邓喀中. 基于DS-InSAR的西部矿区地表时序沉降监测与分析. 金属矿山, 2022(5): 160–169

[5] 赵立峰, 范洪冬, 渠俊峰, 等. 基于DS-InSAR的 张双楼煤矿长时序地表形变监测方法. 金属矿山, 2021(8): 142–149

[6] 李毅, 蒋金雄, 杜玉玲, 等. 融合分布式目标的矿区采动地表时序InSAR监测. 中国矿业大学学报, 2020, 49(6): 1199–1206

[7] Yang Chengsheng, Lu Zhong, Zhang Qin, et al. Ground deformation and fissure activity in Datong basin, China 2007–2010 revealed by multi-track InSAR. Geomatics Natural Hazards & Risk, 2019, 10(1): 465–482

[8] Xi Ning, Mei Gang, Liu Ziyang, et al. Automa- tic identification of mining-induced subsidence using deep convolutional networks based on time-series In-SAR data: a case study of Huodong mining area in Shanxi Province, China. Bulletin of Engineering Geo-logy and the Environment, 2023, 82(3): 78

[9] Yao Sheng, He Yi, Zhang Lifeng, et al. A ConvLSTM neural network model for spatiotemporal prediction of mining area surface deformation based on SBAS-InSAR monitoring data. IEEE Transactions on Geo-science and Remote Sensing, 2023, 61: 1–22

[10] Ferretti A, Fumagalli A, Novali F, et al. A new algorithm for processing interferometric data-stacks: SqueeSAR. IEEE Transactions on Geoscience and Re-mote Sensing, 2011, 49(9): 3460–3470

[11] Ferretti A, Prati C, Rocca F. Permanent scatterers in SAR interferometry. IEEE Transactions on Geoscience and Remote Sensing, 2001, 39(1): 8–20

[12] Berardino P, Fornaro G, Lanari R, et al. A new algo-rithm for surface deformation monitoring based on small baseline differential SAR interferograms. IEEE Transactions on Geoscience and Remote Sensing, 2002, 40(11): 2375–2383

[13] Zhang Zhiliang, Zeng Qiming, Jiao Jian. Deformations monitoring in complicated-surface areas by adaptive distributed scatterer InSAR combined with land cover: taking the Jiaju landslide in Danba, China as an exam-ple. ISPRS Journal of Photogrammetry and Remote Sen-sing, 2022, 186: 102–122

[14] 黄俊松, 曾琪明, 高胜, 等. 适用于自然地表形变反演的小基线集方法. 地球信息科学学报, 2018, 20(4): 440–451

[15] 黄俊松, 曾琪明, 焦健, 等. 一种适用于时序InSAR的同质像元加权干涉相位滤波方法. 北京大学学报(自然科学版), 2018, 54(6): 1242–1250

[16] 曾琪明, 张志亮, 焦健. 基于土地类型的InSAR时序分析中DS自适应选取方法[P]. 北京: CN 1121 30148B, 2021–12–28

[17] Jia Hongguo, Zhang Hao, Liu Luyao, et al. Landslide deformation monitoring by adaptive distributed scatte-rer interferometric synthetic aperture radar. Remote Sensing, 2019, 11(19): 2273

[18] Goel K, Adam N. A distributed scatterer interferometry approach for precision monitoring of known surface deformation phenomena. IEEE Transactions on Geo-science and Remote Sensing, 2014, 52(9): 5454–5468

[19] Liu Jihong, Hu Jun, Bürgmann R, et al. A strain-model based InSAR time series method and its application to the Geysers Geothermal Field, California. Journal of Geophysical Research-Solid Earth, 2021, 126(8): e2021 JB021939

[20] Liu Jihong, Hu Jun, Li Zhiwei, et al. A method for measuring 3-D surface deformations with InSAR based on strain model and variance component estimation. IEEE Transactions on Geoscience and Remote Sen-sing, 2018, 56(1): 239–250

[21] 李学军. InSAR技术在大同矿区地面沉降监测中的应用研究[D]. 太原: 太原理工大学, 2007: 26–30

[22] Gong Peng, Liu Han, Zhang Meinan, et al. Stable classification with limited sample: transferring a 30-m resolution sample set collected in 2015 to mapping 10-m resolution global land cover in 2017. Science Bulletin, 2019, 64(6): 370–373

[23] Hooper A. A multi-temporal InSAR method incur-porating both persistent scatterer and small baseline approaches. Geophysical Research Letters, 2008, 35 (16): 96–106

Investigation of Ground Deformations in Mining Areas Using the Adaptive DS-InSAR Method Combined with Land Cover

ZHANG Zhiliang1,2, ZENG Qiming1,2,†, YANG Ligong1,2

1. Institute of Remote Sensing and Geographic Information System, School of Earth and Space Sciences, Peking University, Beijing 100871; 2. Beijing Key Lab of Spatial Information Integration & Its Applications, Beijing 100871; † Corresponding author, E-mail: qmzeng@pku.edu.cn

Abstract Traditional multi-temporal InSAR (MT-InSAR) technology commonly employs a global single threshold in selecting valid pixels. In mining areas where there are significant temporal and spatial decorrelation effects or where the deformation gradient is large in a short period, the issue of sparse measuring points and inadequate spatial sampling frequently arises, further causing the incomplete acquisition of ground subsidence information. Here, we take 22 scenes of ALOS-1 images covering the Datong mining area as an example and use the adaptive distributed scatterer InSAR combined with land cover (ADSI-CLC) method to measure the long term time series of ground deformations in this area. The results indicate that the spatial-temporal distribution pattern of the deformation measurements obtained by ADSI-CLC method is similar to that of the StaMPS-SBAS method, and has a good corre-lation with mining facilities on high-resolution optical remote sensing images. However, the ADSI-CLC method can significantly increase the number and spatial distribution density of measuring points, identifying approximately four times as many distributed scatterers (DS) points as the StaMPS-SBAS method in this study area. Specifically, the consistency between the two methods is high in areas with zero deformation or minor deformation. In regions with large deformation, the StaMPS-SBAS method cannot effectively obtain deformation measurements, while the ADSI-CLC method can derive deformation measurement results that are consistent with the funnel-shaped spatial distribu-tion characteristics caused by mining activities, because of the increased number of measuring points. These results indirectly verify the reliability and effectiveness of the ADSI-CLC method in the ground deformations monitoring of the mining areas. Overall, the ADSI-CLC method can provide more detailed temporal and spatial deformation information, and can be used to monitor and warn the surface stability in mining areas.

Key words ADSI-CLC method; distributed scatterer; deformations monitoring of mining areas; Datong mining area