融合DS-InSAR与PIM的煤矿开采影响范围确定方法

刘金霖，谭志祥，范洪冬

(中国矿业大学 环境与测绘学院，江苏 徐州 221116)

煤矿开采会导致地面沉降，对沉陷区内的房屋、铁路、公路等工建(构)筑物造成一定程度的损害。我国部分煤矿资源位于村庄下方，此类煤矿开采时一定程度上会导致村庄房屋损毁；无开采作业时，很多房屋因质量等原因也会造成房屋出现裂缝等问题，因此煤矿企业仅希望承担因煤炭开采引起的损毁责任，然而损害原因不清经常造成矿村纠纷，影响矿区和谐发展，因此开采损害技术鉴定是一项非常必要的工作[1]。

开采损害技术鉴定包括确定煤矿的开采影响边界和影响程度两个方面的内容。其中，开采损害影响程度可根据相关规范要求进行判定，一般将建筑物的破坏情况划分为Ⅰ、Ⅱ、Ⅲ、Ⅳ等级。然而，确定开采影响的边界比较复杂[2]。煤矿开采损害鉴定离不开对开采沉陷影响因素、监测手段、沉陷规律、沉陷预计等的研究。蒋旭梓[3]等对济宁井工矿区不同开采条件对沉陷阶段的影响进行了研究；孙国庆[4]等对顶板下沉因素进行了分析；康新亮[5]等对多种地质下开采沉陷规律分析进行了研究；张向阳[6]等对沉陷区内路基沉降变形规律进行了研究；蒋春梅[7]等对地面下沉规律进行了分析；孙庆先[8]等提出了物探技术和开采沉陷学理论相结合的开采损害技术鉴定方法；文献[9-10]等对InSAR测量技术进行了研究；邓喀中[11]等利用D-InSAR技术对开采沉陷进行了监测研究；郭山川[12]等利用DInSAR技术对黄土高原矿区地表形变进行了监测；李柱[13]等利用DS-InSAR技术对乌达煤田火区长时序地表形变进行了监测；张文龙[14]等结合PS-InSAR和DS-InSAR技术对韩城矿区采空区地面塌陷进行了识别；杨隽[15]等将SBAS-InSAR技术对形变特征的监测用于矿区损害鉴定分析；李学良[16]对煤矿老采空区覆岩移动变形监测方法进行了分析；吴群英[17]等对荒漠化矿区的生态环境进行了监测；陈毓[18]对地质灾害风险评估与适应性评价进行了研究；张官进[19]等对矿区的开采沉陷进行了预测；谭志祥[20]等提出在准确掌握采矿资料时，鉴定方法有地表移动变形预计、移动角计算、相似材料模拟以及根据房屋和地表裂缝特征判别等方法；李晓斌[21]等对高强度开采地表损伤程度分类进行了研究。

上述研究为开采影响边界研究提供了有益借鉴，但针对多个工作面较长时间开采导致的房屋损害鉴定问题，传统方法为理论预计结果，缺少实测数据支撑，导致矿村双方对边界位置存在较大争议，而常规InSAR技术由于时空失相干，无法获取长时间的地表沉降信息。

基于此，笔者提出了一种基于DS-InSAR的分时段累加方法，实现对华东地区某村庄2016年7月至2021年10月期间的地表沉降监测，同时结合概率积分法，综合分析划定了该煤矿3个工作面的开采影响边界，得到了矿村双方的认可，解决了矿村矛盾，为开采损害鉴定提供了新的方法。

1 试验原理及流程

1.1 概率积分法

1965年，我国学者刘宝琛、廖国华根据波兰学者J.LITWINISYN的随机介质理论提出了概率积分法，并在我国广泛应用[22]。根据概率积分法，煤矿开采后引起地面某点的下沉值为

式中，L1，L2分别为工作面走向长度和倾向宽度；1l，2l分别为沿走向、倾向计算的开采宽度；m为煤层厚度；S0，Sdown，Sup分别为走向、倾向下山和上山拐点偏移距；Wmax为充分采动时地表最大下沉值；α为煤层倾角；q为下沉系数；θ为开采影响传播角；β为主要影响传播角；r为主要影响半径；H为采深。

1.2 DS-InSAR分时段累加方法

1.2.1 FaSHPS 同质点选取原理

同质像元识别是指识别在地表形变监测区域内SAR影像数据集中具有相同变化特征的像元，包括参数检验和非参数假设检验的两种识别方法。笔者选用参数检验中的快速同质点选取方法(FaSHPS)[23-24]，它是将假设检验问题转化为置信区间估计问题，比较两个像元是否服从相同的分布函数，从而判断它们是否为同质像元。

根据中心极限定理，在SAR影像数目足够的情况下，像元平均振幅A(L) 可以近似服从高斯分布，区间估计可表示为

其中，P{＊}为概率；μ(L) 为像元L的数学期望；z1-α/2为置信度为α时对应的分位点； Var[A(L)]为像元L幅度的真实方差；N为样本数。若像元后向散射特性稳定，变异系数为常量(0.52)，则式(2)变换为

根据式(3)，若平均振幅A(L) 在该置信区间内，即可判定两者为同质像元。

1.2.2 特征值分解相位优化

由于分布式目标像元内地物的后散射特性是相同的，且都不占据像元信号的主导地位，其相位稳定性也较差，因此要对相位进行优化，以达到去除噪声相位的目的。试验选用特征值分解法[25-27]进行相位优化，其通过求解相干矩阵的特征值和特征向量，进而确定不同散射信号对应的特征值，并找到最大特征值对应的特征向量，其所代表的散射体被视为主导散射体，从而实现对相位的优化处理。

定义分布式目标的相干性矩阵T的表达式为

式中，y=[y1,y2,y3,… ,yN]为分布式目标的同质像元在N景影像上的复数观测量经过时间维振幅归一化处理后的复数观测值；NSHPs为同质点数量；Ω为全部分布式目标同质点的集合。

由特征值分解可得：

式中，λi为相干矩阵的特征值，当λi越大时，其对应的散射相位越占优势；μi为对应的特征向量。

因散射体间相互独立，则相干性矩阵可分解为主散射体信号Tsignal和噪声信号Tnoise，即

式中，λ1为最大特征值；μ1为其对应的特征向量，即优化后的分布式目标相位。

在同质像元选取和相位优化完成后，通过计算拟合度，设置时间相干性阈值，从而选取相干性较好的分布式目标点，进行时序InSAR处理，便可获得该时间段内研究区的地表沉降信息。

1.2.3 三次插值法

1959年，DAVIDON首先提出三次插值法，其原理是通过构建多项式来实现插值。

定义一个在给定区间[m,n]上的三次多项式为

求解得到：

假设，＊α为区间[m,n]上的插值点，则

求解得到插值点＊α为

根据三次插值原理，对每个时间段的地表沉降结果进行插值后累加，即可获得整个时间段内研究区的地表沉降结果。

1.3 试验流程

试验的具体处理流程如图1所示。①根据已知地质资料，选取合适的预计参数，利用概率积分法预计鉴定区的沉降情况；②利用快速同质点选取算法和特征值分解相位优化方法获得DS点；③经过相位解缠即可提取地表形变情况，对每个时间段的沉降结果进行插值叠加后可获取鉴定区长时序的地表形变信息，对PIM预计结果进行修正、验证；④确定多个工作面的开采影响边界，完成开采损害鉴定。

图1 试验流程Fig.1 Experimental flow chart

2 研究区概况

研究区域位于华东地区中部，其地表属淮河冲积平原，地形平坦，海拔标高+23.52～+27.58 m，地势西高东低，地表水系密布。需鉴定的村庄北侧分布有3个工作面，分别为121302，121303，121304工作面，工作面位置如图2所示。

图2 工作面与村庄位置Fig.2 Position of the working face and the village

由图2中工作面与村庄位置可知，工作面自北向南推进，采用后退式倾向长壁一次采全高采煤法采煤，全部垮落法管理顶板。

3个工作面的具体参数见表1，鉴定区表土层厚度580 m，工作面上覆岩层以砂岩、泥岩为主，岩性中硬偏软。经过笔者现场勘察和调研，村庄部分房屋受到不同程度的损坏，村民认为是工作面开采导致，但是企业认为是由于村民在村庄北侧大量取土引起房屋部分地面及墙体产生裂痕。因此，需要对该村庄进行开采损害鉴定，以判定村矿之间的责任。

表1 工作面参数Table 1 Parameters of the working face

3 概率积分法预计结果分析

文献[28-30]中给出了该矿首采工作面及附近煤矿10多个地质采矿条件相似工作面的岩层移动参数，结合鉴定区工作面具体情况，通过综合分析，确定采用的概率积分法预计参数为：下沉系数q=0.99，水平移动系数b=0.35，主要影响传播角正切tanβ=1.8，开采影响传播角θ=86°，拐点偏移距s=0.03H。

基于上述参数，首先利用中国矿业大学开采损害及防护研究所研制的矿区开采沉陷预测分析系统(MSAS)，对开采后引起的村庄地表沉陷情况进行了计算，其下沉等值线如图3所示。

图3 研究区预计下沉等值线Fig.3 Anticipated subsidence contour lines in the study area

由图3中的研究区预计下沉等值线结果可知，3个工作面的开采使村庄部分地表产生了沉降，影响最大的位于村庄东北部，开采影响边界位于村庄南部，附近煤矿的开采对村庄影响范围较广。矿村双方认为概率积分法是理论预计模型，不能反映实际地表沉降情况，对开采影响边界具体位置存在较大争议。因此笔者利用DS-InSAR分时段累加方法，获取研究区地表沉降结果，对概率积分法预计结果进行修正和验证。

4 DS-InSAR结果分析

4.1 试验数据

研究使用欧空局发射的Sentinel-1A卫星获取的试验数据，从中选取153景覆盖研究区域的卫星升轨影像，时间跨度为2016年7月至2021年10月。试验采用了90 m分辨率的航天飞机雷达测图计划数字高程模型作为外部DEM数据，以实现地形相位去除的目的。为了降低时空失相干的影响，试验将影像分为2016年7月至2018年6月，2018年6月至2020年6月和2020年6月至2021年10月3个时间段，采用多主影像策略，对数据进行DS-InSAR技术处理。

4.2 DS-InSAR结果

4.2.1 可靠性分析

为验证DS-InSAR监测结果的可靠性，选取2021年3月至10月村庄内23个监测点的水准实测数据，与DS-InSAR监测结果进行对比分析。其水准点与村庄的位置关系如图4所示。DS-InSAR监测结果与水准实测数据对比结果如图5所示。将水准实测数据视为真值，DS-InSAR监测结果作为测量值，计算得到DS-InSAR监测结果的最大误差、平均误差以及均方根误差，结果见表2。

表2 DS-InSAR监测结果精度评价Table 2 Accuracy evaluation of DS-InSAR monitoring results

图4 水准点与村庄位置关系Fig.4 Relationship between the leveling points and the village

图5 InSAR与水准方法监测结果对比Fig.5 Comparison of monitoring results between InSAR and leveling

结合图5和表2分析可知，DS-InSAR获取的地表沉降趋势基本与水准数据一致，误差较小，证明DSInSAR技术在该研究区的监测结果具有较高的可靠性。

4.2.2 研究区时序监测结果分析

对于每个时间段，计算其他SAR影像相对于第一景影像的沉降值，然后将其从卫星的视线方向(Line of Sight，LOS)转换为垂直方向，即可得到研究区域在该时间段内的累积地表沉降值。笔者将累积沉降图与地下工作面的位置进行叠加，以便更加直观展示出沉降区域的形成与工作面开采之间的关系。2016年7月至2018年6月，2018年6月至2020年6月和2020年6月至2021年10月3个时间段的研究区地表累积沉降信息如图6所示。

图6 各时间段研究区地表累计沉降Fig.6 Cumulative surface subsidence in the study area for each time period

由图6(a)可知，2016年7月至2018年6月研究区地表出现了较为明显的沉降区域，其位置基本与121303工作面和121304工作面相对应。由于工作面区域地表沉降较大、中央积水，DS-InSAR无法监测到工作面中心的沉降值，可监测到的最大沉降值为782 mm，位于121303工作面和121304工作面边缘，村庄内监测到的最大沉降值为309 mm，位于村庄北部。由此可见，121303工作面和121304工作面的开采对村庄中部和北部均产生一定影响。由图6(b)可知，2018年6月至2020年6月的DS-InSAR监测结果中也出现了较为明显的沉降区域，其位置基本与121302工作面相对应，监测到的最大沉降值为804 mm。该时间段内村庄内依旧受到工作面开采的影响，地表产生了一定沉降。根据图6(c)可知，2020年6月至2021年10月，121303，121304和121302三个工作面均开采结束，但是工作面附近仍有残余沉降，虽然其沉降值明显减小，但仍对村庄东北部产生了一定影响。

为了更好地展示2016年7月至2021年10月期间，121303，121304和121302工作面的开采对村庄产生的沉降影响，笔者将3个时间段的DS-InSAR监测结果，通过插值处理后，对同一地理位置的沉降值进行叠加，获得了累计沉降情况，如图7所示。根据DS-InSAR叠加沉降结果，3个工作面的开采对地表产生最大沉降值约为1 496 mm，位于工作面边界，对村庄影响较大，村庄地表产生了较为明显的下沉，计算得到的最大沉降值约为576 mm，位于村庄东北部。

图7 2016年7月至2021年10月研究区地表累计沉降Fig.7 Cumulative surface subsidence of the study area from July 2016 to October 2021

4.2.3 监测结果修正与对比分析

由于概率积分法预计结果在采空区边界上方收敛较快[31]，导致其预计得到的下沉边界较实际要小，笔者根据DS-InSAR技术获取的实际地表沉降结果，绘制出研究区地表下沉等值线，据此对概率积分法预计得到的村庄下沉等值线进行向外偏移修正，最终划定出3个工作面对村庄的开采影响边界，位置如图8所示。可见，开采影响边界以内的房屋等建筑均受到3个工作面开采的影响。

图8 研究区DS-InSAR下沉等值线Fig.8 Subsidence contour lines in the study area using DS-InSAR

通常开采沉陷影响边界由岩层移动的边界角确定，但是边界角受采深、采宽、表土层厚度、覆岩岩性、采动程度等诸多因素的影响，通常很难准确选取，用类比法获得的边界角往往存在较大误差，从而导致其确定的开采沉陷边界与实际相差较大；由上述研究结果可知，InSAR技术在确定开采区域边界上方小变形时结果较准确，在条件许可情况下可采用InSAR技术结合PIM方法确定开采沉陷影响边界。

5 结 论

(1) 利用DS-InSAR对时间跨度长达5 a的153景Sentinel-1A数据进行分段处理，并将分段沉降监测结果进行叠加，获取了研究区域长时序大变形的地表沉降信息，有效减小了时空失相干对InSAR监测结果的影响，提升了InSAR技术在长时间地表形变监测方面的能力。

(2) 通过DS-InSAR对概率积分法的预计结果进行验证及修正，划定了开采影响边界，得到了矿村双方的认可，为开采损害技术鉴定提供了新的参考，丰富了开采沉陷学领域的学科内容。