卢欣奇，李学峰，张勤斌，黄海斌

(1.北京矿冶科技集团有限公司，北京 100160； 2.广西大学资源环境与材料学院，广西 南宁 530004)

传统的矿区开采沉陷监测方法存在监测成本高、效率低、受天气影响以及无法实现实时监测等问题，并且各类开采沉陷监测方法与开采沉陷预计算法难以进行有机集成[1-6]。合成孔径雷达干涉测量(interferometric synthetic aperture radar,InSAR)是通过对相同研究区域内的两景或多景SAR影像数据进行相干性处理，获得研究区域内高精度的三维地形信息[7]。对其进行时间序列分析可以得到研究区域内地形微小变形信息的一种区别于传统测量的新技术。该技术可以同时达到对研究区域的大范围、低成本、高精度、高效率的变形监测，其理论精度可达到毫米级，这些优势是传统监测手段无法达到的。卫星雷达干涉测量技术在地面沉降监测方面已显示了独特的优势和高精度，无疑将成为未来监测地表形变的主要技术之一[8]。

2001年FERRETTI等[9]提出PS-InSAR(permanent scatters InSAR)技术，通过永久散射体具有较高相干性，受空间、时间影响比较小的特性，为降低信号噪声的影响提供了新的解决方式。2017年白泽朝等[10]利用Sentinel-1A雷达影像采用PS-InSAR技术对天津地区地面沉降进行了监测，综合水文地质情况对天津地区的地表沉陷情况及原因进行了详细的分析说明，其监测效果较为突出，为非商业SAR影像进行地表沉降监测提供了新的思路。李曼等[11]等通过PS-InSAR技术对唐山市地面沉降灾害较严重的地区之一曹妃甸新区的地面沉降发育特征及其影响因素进行了详细的分析，得出了该地区地面累积沉降量分布及演化状况。地下水超采、大规模工程扰动是诱发和加剧地面沉降的外在动力的结论，为PS-InSAR研究地表沉降与实际扰动情况相结合提供了新思路。黄佳璇[12]对乌东德区域蠕动型滑坡进行检测中使用了PS-InSAR技术与ArcGIS软件结合的方式，对该地区的沉降变化规律进行总结，使PS-InSAR技术测量的结果直观化、整体化，使沉降规律的总结数字化，对沉降的发展情况更具预测性。邹昊等[13]应用PS-InSAR技术的对老采空区地表沉降监测进行了研究，将PS-InSAR技术应用到了煤矿的采空区监测当中，综合工作面的布置情况对采空区地表变形规律进行了总结和预测。PS-InSAR技术监测大范围内采空区地表沉降具有较大优势，为采空区地表沉降的预防、规律总结及治理提供重要的数据基础。

本文针对广西平南锡基坑铅锌矿矿区范围进行研究，利用2015年6月26日至2018年1月4日期间覆盖概况区部分的21景Sentinel-1影像开展基于小数据集的PS-InSAR技术矿区监测地面沉降研究，分析其变形规律，为预测和评价老采空区残余变形提供基础。

1 PS-InSAR技术

PS-InSAR技术是20世纪末由意大利学者首先提出的，以解决常规干涉中大气影响、失相干、DEM误差等问题，极大地拓展了InSAR技术的应用前景，为精确研究地壳形变提供了强有力工具[14]。一般情況来说，PS(永久散射体)点位数量城市区域每平方千米可在数十个点以上，而郊区部分也可达每平方千米内有几个点，这样的资料密度，已经远超过多数地区GPS的站位密度。

其原理为利用永久散体在长时间内保持稳定反射的特性[15]，其反射信号大大高于信号噪声，将这些永久散体提供的信号在时间序列上的相位进行排列，在一定特征尺度范围内大气效应一致的假设前提下可以将InSAR结果中的大气延迟误差和DEM误差从信号中消除，从而达到利用永久散射体对地表变形的精密观测。

2 实验区概况与数据处理

试验区位于贵港市某矿区西北部，矿区范围内有罗马村、农昌村等村庄。全矿区内均为第四系所覆盖，早在1998年进行开采活动，主要为出露地表的矿体及至地表以下30～40 m之间的铅锌矿体，矿区范围内+20 m以上主要矿体已基本采完，目前地下存在大面积采空区。矿石以铅锌矿为主呈层状、似层状产于中泥盆统东岗岭组下段白云岩中，顶板、底板均为白云岩。

实验选取该矿1～6号采空区进行观测，观测期间6个工作面均已停采，该实验是对老采空区地表沉降的变形监测。其中6号采空区面积最大，3号采空区距离地表最近，PS-InSAR技术数据处理流程见图1。实验数据为21景高分辨率Sentinel-1A卫星IW模式影像数据，影像选择及时间基线见表1。采用PS-InSAR技术对数据进行分析得出超级主影像，对21幅SAR单视复数影像，经配准、辐射定标、PS探测和干涉处理，并借助分辨率为90 m、高程精度约为10 m的SRTM-3的DEM文件进行差分干涉处理，得到20幅干涉和差分干涉图、研究区域内的5 200个PS点以及各PS点的差分干涉相位集。在考虑地表变形、高程误差、大气误差、DEM误差等其他失相干因素影响的前提下，得到每个PS点在每幅差分干涉图上的差分干涉相位组成，其中，对形变速率增量和高程误差增量积分，可以得到每个PS点相对于主参考点的形变速率和高程误差。对影像进行灵活的时序分析叠加，将相邻影像进行干涉处理生成影像干涉对以保证影像的相干性。

图1 PS-InSAR技术数据处理流程Fig.1 PS-InSAR technology data processing flow

表1 影像的选择及时间基线Table 1 SAR selection and SAR time baseline

序号时间时间基线序号时间时间基线12015-06-25-122017-02-143622015-08-1248132017-04-034832015-09-2948142017-05-214842015-11-1648152017-06-021252016-01-0348162017-07-083662016-05-26144172017-08-254872016-08-0672182017-09-182482016-09-2348192017-10-122492016-11-1048202017-11-2948102016-12-2848212018-01-0436112017-01-0912

最后将雷达坐标系下的PS矢量点信息导入ArcGIS进行坐标变换、重新投影坐标系转化、克里金面生成、与采空区分布平面图进行叠加，最终将PS点位信息可视化，得到研究区域内采空区的沉降量分布云图、平均沉降速率云图。

3 试验结果与讨论

3.1 试验结果

通过对覆盖该矿区的Sentinel-1数据进行处理取得了较好结果，在实验区域内有5 200多个相位稳定PS点被选取出来。将时序处理得到的沉降信息结果文档导入到ArcGIS中并绘制成整个区域的沉降信息图，调节沉降数值的显示分布，叠加强度影像图作为底图，得到各采空区沉降量分布云图，如图2所示，并根据累计沉降量得出研究期间平均沉降速率图，如图3所示。

由图2可知，在采空区范围内，最大沉降量主要集中在6号采空区北部，最大沉降量为205 mm；其次为4号采空区和5号采空区，最大沉降量为195 mm。与其他采空区的沉降量相比，1号采空区沉降量较小，沉降量为127～80 mm；3号采空区和6号采空区交接也出现了部分沉降，沉降量为146 mm。其他部分区域表现出地表抬升现象，抬升区域主要是3号采空区中部(该矿区主井口附近)，为该矿区的废石堆积处。

由图3可知，1号采空区沉降速率小，低于4.27 mm/a；3号采空区中部抬升速率高于10 mm/a属于人为堆积造成；1号采空区、2号采空区、4号采空区、5号采空区相比6号采空区的其他区域沉降速率较小，均低于11.73 mm/a；6号采空区的北部区域主要表现出较为严重的沉降情况，沉降速率较大，可达到13 mm/a左右，最大处可达到15.46 mm/a。整体来看，6号采空区沉降速率较大，沉降区域范围较大且连成一片，形成一个典型的沉降漏斗。

图2 2015～2018年矿区沉降量分布云图Fig.2 The settlement distribution nephogram of the mining area from 2015 to 2018

图3 2015～2018 年矿区平均沉降速率云图Fig.3 The average settlement rate nephogram of the mining area from 2015 to 2018

3.2 剖面线情况讨论

为了详细研究沉降漏斗区域的分布特征及其时序规律，在图2中选取矿区内采空区沉降漏斗中心处分别绘制平行于勘探线的剖面线，分析沉降漏斗中心线两侧区域的沉降变化情况，6号采空区北部沉降量最大，1号采空区沉降量最小，3号采空区与6号采空区交接位置也有沉降产生，故对1号采空区、6号采空区及3号采空区与6号采空区交接位置进行详细分析，其他剖面情况在以下讨论中不予列出。

由图4可知，1号采空区剖面线西侧呈现三“漏斗状”并列分布，最大沉降点东侧“漏斗面”较为平滑。由采空区沉降量分布云图可知，该沉降区域的最大沉降点A为1号采空区的东南方向，应属于1号采空区顶板的薄弱部位，在最大沉降点西侧有两个明显的小型漏斗区域存在，沉降量为51 mm、43 mm，为最大沉降量点A的47%、39%。

由图5可知，3号采空区与6号采空区交接部位剖面线东侧呈现两“漏斗状”梯度分布，西侧呈“V”型阶梯状分布。由采空区沉降量分布云图可知，该抬升区域的最大抬升点C为3号采空区与6号采空区交接部位，应属于3号采空区与6号采空区连接部位，在最大沉降点西侧有明显的小型漏斗区域存在，沉降量为46 mm，为最大沉降量点C的32%。

图4 Ⅰ-Ⅰ等值线剖面沉降情况Fig.4 Ⅰ-Ⅰ contour section settlement

图5 Ⅲ-Ⅲ等值线剖面沉降情况Fig.5 Ⅲ-Ⅲ contour section settlement

由图6可知，6号采空区剖面线西侧较为平缓东侧较陡且漏斗面较为平滑，整体呈现“V”形沉降梯度分布，由采空区沉降量分布云图可知，该沉降区域的最大沉降点F为6号采空区的东南方向，应属于6号采空区顶板的薄弱部位。在最大沉降点西侧平缓部位，并未出现明显漏斗区域，平均沉降量为87 mm，为最大沉降量点F的41%。

3.3 最大沉降点沉降速率情况讨论

为了详细研究沉最大沉降点沉降速率情况，对最大沉降点进行沉降速率分析，其各点速率变化曲线如图7所示。

图6 Ⅵ-Ⅵ等值线剖面沉降情况Fig.6 Ⅵ-Ⅵ contour profile settlement

图7 最大沉降点沉降速率曲线Fig.7 Settlement rate curve of the maximum settlement point

由点A的沉降速率变化曲线可知，A点的沉降速率出现了三次峰值。该点沉降速率一直缓慢增加，达到极值后逐渐下降，持续一定时间又出现了加速状态，沉降速率均匀，在达到第二次峰值后又出现了沉降速率减小的情况。大体上呈现“平衡状态-非平衡状态-新平衡状态”这一开采沉陷规律，当采空区顶板在沉陷过程中达到平衡状态后，又会出现新的失稳状态，当顶板及矿柱经历过失稳过程后又一次达到平衡时即进入到新的平衡状态。

由点C的沉降速率变化曲线可知，C点的沉降速率出现了三次峰值。该点沉降速率最初缓慢减少，第一次达到极值后沉降速率急速上升，达到第二次极值后沉降速率又减小，持续一定时间到达第三次极值后沉降速率又开始增加。造成这种情况的原因为达到第一次峰值后受到了下部开采扰动的影响，使原本即将到达稳定状态的地表产生了变化，产生了一段时间沉降速率较大的变形，达到第二次峰值。但这种变化并未持续较长时间，就产生了逐渐趋于稳定的状态，直到第三次峰值的出现。在第三次峰值之后因受到其他扰动影响沉降速率先快速增加，后缓慢增加，且有即将达到峰值的趋势。

由点F的沉降速率变化曲线可知，该点沉降速率最初为较小，后沉降速率逐渐增大，达到极值后逐渐下降，持续一定时间又出现了加速状态，且加速度逐渐增加。

C点与F点均受到扰动影响，其沉降速率不同的原因是下部矿体开采工作面的推进导致距离这两点的距离逐渐改变导致的。当下部开采扰动距离上部采空区较近时则该区域先出加速状态，在扰动过后，采空区整体的稳定性下降，新的平衡状态将会后延，当开采扰动影响过后采空区依旧会达到新的平衡。

4 数值模拟验证

4.1 岩体力学参数及矿柱平均强度

岩体力学参数是依据岩石力学参数特性测试结果，并考虑了岩体的结构效应、地下水、节理裂隙等因素，对岩石力学参数按照POPOV的RMR岩体质量分类法进行适当的修正[16]。根据长沙矿山研究院在“金属矿山安全技术国家重点实验室”对矿体及围岩采用MTS实验系统进行实验得出的岩石力学参数，选用饱和岩样实验结果，使用加拿大Rocscience公司开发的Roclab1.0软件计算矿体及顶板、底板围岩的岩体力学参数，从而确定基于数值模拟的矿山采空区稳定性分析所需的岩体力学参数，岩体及后期采用的充填体力学参数见表2。

考虑到矿柱尺寸效应的影响，确定了上述的岩体力学参数后，采用式(1)进行矿柱强度估算。

(1)

式中：QP为矿柱强度，MPa；Qr为岩体强度，MPa；B/H为宽高比。

各采空区的平均矿柱强度见表3。

表2 矿岩介质的力学参数Table 2 Mechanical parameters of rock medium

表3 平均矿柱强度统计表Table 3 StatisticalTable of average strength of pillar

4.2 计算分析

4.2.1 最大主应力分析

由于篇幅限制故只对以上沉降较为明显的位置进行分析，其他采空区模拟结果不予列出。

图8所示为各采空区顶板及矿柱纵剖面、横剖面的最大主应力云图。

为确保采空区经治理后达到长期稳定，矿柱安全系数采用保守方式计算，计算公式见式(2)。

F=QP/σPmax

(2)

式中：F为矿柱安全系数；QP具体数值详见表2；σPmax为矿柱的最大垂直应力，MPa，由数值模拟给出。

矿柱的最佳强度安全系数由式(3)确定[5]。

n=1+t×σ

(3)

式中：n为矿柱最佳安全系数；t为相当于规定的可靠性系数，根据该矿山实际情况，确定t值为3.48；σ为强度安全系数的均方差，根据该矿山地质条件，计算得出σ=0.25。将数值带入(3)式得n值约为1.9，即当F

采空区顶板、矿柱的应力状态、不稳定矿柱个数等信息详见表4。

表4 采空区应力状态统计表Table 4 StatisticalTable of stress state in goaf

图8 最大主应力云图Fig.8 Maximum principal stress nephogram

由图8(c)和图8(f)可知，在矿柱上的最大主应力呈现压应力状态，且最大主应力出现于矿柱的中上部位，因各采空区矿柱众多，矿柱的最大主应力云图将按上述规律在矿柱的中上部位做横剖面图，各采空区横剖面的最大主应力云图见图8(b)、图8(d)和图8(e)。

由最大主应力可知，在1号采空区的中部、3号采空区整体及6号采空区北部及南部均有大量不稳定的矿柱集中，因此这些部位为各区域的薄弱位置。

4.2.2 最小主应力分析

由图9可知，最小主应力呈现出压应力及拉应力状态。1号采空区、2号采空区、4号采空区、5号采空区及6号采空区仅顶板局部小范围内出现拉应力，拉应力范围是0～0.20 MPa；3号采空区顶板绝大部分为拉应力状态，拉应力范围是0～0.27 MPa，而顶板围岩抗拉强度仅0.19 MPa，说明3号采空区顶板有较大几率出现拉伸破坏。

4.3 结果对比

PS-InSAR测量得出的沉降量较大区域与数值模拟得出的薄弱区域对照信息详见表5。

由表5可知，该矿区1号采空区和6号采空区通过PS-InSAR技术测量得到的沉降量较大区域与数值模拟所得的薄弱区域基本一致；3号采空区数值模拟得出的薄弱区域为整个采空区范围，而通过PS-InSAR技术测量得出的较大沉降量的部位为3号采空区的北部。该情况产生的原因为3号采空区南部的废渣堆放导致了该部位地表变形出现上升。在废渣堆放处边缘有明显下沉边界，该边界外部为3号采空区与6号采空区交接部位(即3号采空区与6号采空区交界位置)，在该区域出现了沉降量较大的区域，综合PS-InSAR测量得出的沉降量较大区域与数值模拟得出的薄弱区域进行对比可知，通过PS-InSAR技术对采空区沉降监测的数据具有较高的真实性，可以为采空区地表沉陷防治及沉降规律总结、采空区治理提供数据支持。

图9 最小主应力云图Fig.9 Minimum principal stress nephogram

表5 数值模拟及实测沉降区域对照表Table 5 Numerical simulation and measured settlement area check list

采空区编号PS-InSAR测量得出的沉降量较大区域数值模拟得出的薄弱区域1号中部中部3号北部整体6号北部及南部北部及南部

5 结 论

1) 利用PS-InSAR技术对地表沉降进行检测，绘制出地表沉降量分布云图、等值线剖面、观测点沉降速率变化曲线。与传统方法相比大大提高了地表沉降研究的直观性、整体性及预测性。

2) 通过数值模拟对PS-InSAR沉降监测结果进行验证，对比结果基本一致，因此PS-InSAR技术获得的采空区沉降数据可以为地表沉陷防治及沉降规律总结、采空区治理提供支持。

3) 1号采空区没有受到下部采动影响，其沉降趋势较为稳定，符合“平衡状态-非平衡状态-新平衡状态”这一开采沉陷规律。

4) 3号采空区与6号采空区交接位置受到开采扰动影响出现沉降速率较快增加的现象，达到峰值后整体沉降趋势较为稳定，但在数据截至时沉降速率依旧在加快。

5) 基于分析结果，1号采空区正处于新的平衡状态后的减速下沉状态。6号采空区在受到下部开采扰动的影响后出现了加速沉降状态，沉降的极限情况有待研究。

6) 1号采空区整体沉降变化较为稳定，暂时不需要进行处理。但3号采空区与6号采空区交接位置及6号采空区沉降部位在数据截止时沉降速率持续增加，需对这两个部位进行追踪研究，为避免发生大规模采空区塌陷等灾害，应对这两个沉降速率持续增加的区域进行加固处理。