基于InSAR技术对陕西矿区塌陷地震早期预警的研究

吴 明，严 珊，王金硕，杨 敏，张晨蕾，王玥明

（1.陕西省地震局,西安710038；2.中国地震局第二监测中心,西安710054）

0 引言

由于经济发展带来的资源需求，陕西地区矿产资源开采持续增长，同时也带来地面塌陷等伴生灾害，根据陕西地震台网记录到的塌陷地震来看，其发生的区域基本与矿区位置吻合。随着陕北能源基地的建设和煤炭资源的大量开采，采空区的塌陷灾害越来越严重，伴生的地面形变、地下水流失、土地沙化、房屋开裂等问题对当地居民的生产生活造成了困扰，自2004年以来陕西数字地震台网记录到的塌陷地震事件越来越多[1]。榆林地区位于鄂尔多斯地块，构造稳定，地震活动较弱，自公元前780年至1960年记载的地震共9次，均小于6级，1990年至今未发生过除塌陷类型以外的地震，未见有地震形成的灾害记录[2]。而从中国地震台网记录的2004—2019年地震监测数据来看，陕西省有记录的塌陷地震多达324起，震级最大的是2004年10月14日发生在陕西省榆林市神木县的ML4.2塌陷地震。因此采空区的塌陷地震早期预警研究具有重要意义。

塌陷监测预警分为地下监测和地表监测，地下监测如应力应变监测[3]、土壤含水率监测、爆破震动监测等方法[4]，通常需要在地下埋设大量传感器与设备，监测预警成本高且应用范围小；地表监测包括全球定位系统（GPS）、热红外（TI）、合成孔径雷达（InSAR）、大地水准测量等。传统的地面形变监测手段如大地水准测量、GPS等虽然方法成熟且单点精度较高，但是其仍存在成本高和空间分辨率低等问题。而InSAR技术因具有覆盖范围广、全天时、全天候等特点被广泛应用于区域形变监测研究中[5-6]。1996年Carnec等人[7]首次利用D-InSAR技术监测法国GARDANNE附近煤矿开采引起的地表形变，验证了InSAR技术监测矿区形变的可行性，2001年Berardino等[8]提出了SBAS-InSAR方法。InSAR用于矿区地面沉降监测已有20多年的发展历史，在滑坡监测、同震形变监测等领域获得了很多成果[9-13]。随着InSAR技术的发展，Ps-InSAR、SBAS-InSAR等基于时序的InSAR技术开始应用，在克服形变监测中的大气效应、时空失相干等方面效果更加显著。为了更科学地研究塌陷地震预警，我们综合应用InSAR以及沉降规律模型开展研究。

1 数据源及研究区域

地面沉降是一个时空演化的过程，一些发展比较快速的沉降可能会演化为塌陷，塌陷大多发生于沉降的发展期中，塌陷的发展也是一个时空渐变的演化过程，在这个过程中地表形变速率有其发展的演化规律，那么，在前期速率由慢到快的发展阶段中，如果能够使用InSAR技术手段监测到塌陷发生前的速率特征，那么或许可以为塌陷地震的前期预警提供帮助。塌陷地震引发的地面垂直形变较大，塌陷区相位梯度过大导致失相干较严重，研究最终获得的信息较少，目前基于InSAR的塌陷地震研究相对较少，对于塌陷地震早期预警的研究也存在困难。通过比较发现C波段SAR数据是用于观测塌陷过程活跃阶段沉降边缘变形带的较理想数据，用X波段数据观测会因为相位梯度太小导致可观测变形带过窄，而使用L波段数据则因为波长过长，边缘形变量不够，导致无法在干涉图上形成明显的环状干涉条带。而如果使用时间基线过长的SAR数据做干涉，则会难以判断沉降所处的阶段，所以建议使用时间基线尽量短的数据来观测，时间基线越短，可做出的判断越准确。基于观测效果以及时间成本综合考虑，我们选择C波段Sentinel-1A卫星数据进行处理（数据源Sentinel-1A来自欧空局网站：https://scihub.esa.int/）。

根据中国地震台网数据统计，2004—2019年陕西地区发生的324起塌陷地震，榆林地区达293起，占据了陕西省塌陷地震的绝大多数，同时考虑到本研究采用的哨兵卫星数据从2014年开始升空监测，我们从2014—2019年塌陷地震中选择了震级最大的2次地震作为研究对象，分别为2017年9月5日陕西彬州ML3.2塌陷地震（以下简称彬州ML3.2塌陷地震）和2018年8月16日陕西府谷ML3.2塌陷地震（以下简称府谷ML3.2塌陷地震）（图1）。ML3.2塌陷地震震中位于陕西省北部矿区，陕北地区地层松软，煤矿埋藏较浅，上覆岩层普遍薄。陕北位于鄂尔多斯地块内部，这个块体主要以缓慢的升降震荡运动为主，具有整体性和稳定性的特征，除了块体周缘发育不同程度的褶皱与断裂，块体内部完整且稳定，构造活动较少，没有发生过6级以上地震[1]，在该区域记录到的地震主要是采矿引发的塌陷地震。彬州ML3.2塌陷地震震中位于陕西省彬州市与咸阳市的交界处。

图1 研究区域示意图

2 沉降模型

对于采空区塌陷的理论模型也有很多学者做了相关研究，20世纪90年代钱鸣高院士首次提出了关键层理论用于揭示采空区上覆岩层力学机理，该理论认为采空区上覆岩层中较为坚硬的厚层状岩层控制整个上覆岩层和地表的移动与变形，这种较厚层状的坚硬岩层称为关键层[14]。20世纪50年代波兰学者Knothe对地下采矿引起的的地表塌陷进行了研究，基于地表下沉速度和该时刻下沉量与最终下沉量的差值成正比的假设前提，精确建立了下沉量时间序列模型。由于Knothe提出的模型与很多实际的监测资料拟合的结果不符，李远耀[15]在Knothe模型的基础上，提出了修正的时间序列模型，在Knothe的时间函数上加了一个幂指数k，k为地面下沉系数，增加参数后，时间函数为：

李远耀的模型模拟的沉降过程大致分加速度正负2个阶段。郭增长[16]按照下沉速度对地表的影响程度不同，将地表下沉分为3个阶段：①开始阶段：下沉量达到10 mm的时刻为地表下沉开始阶段，从下沉开始到沉降漏斗下沉速度达到1.67 mm/d为开始阶段。②活跃阶段：沉降漏斗下沉速度大于1.67 mm/d的阶段，该阶段的下沉量占总下沉量的85%～90%，地表的损坏主要发生在这个阶段。③衰退阶段：沉降漏斗下沉速度降到1.67 mm/d以下，并且半年内各点的累计沉降量不超过30 mm。

通过研究分析了许多学者们用各种手段监测得到的塌陷沉降时间序列曲线后，我们发现了一些有趣的规律，总结成如图2的沉降规律模型，根据各时期的沉降量，将沉降分为4个阶段，即开始阶段（stage A）、前活跃阶段（stage B）、后活跃阶段（stage C）、稳定阶段（stage D）。图中2条黑色曲线属于开始阶段（stage A），这个阶段整个沉降漏斗沉降速率都低于1.67 mm/d；紫色曲线属于前活跃阶段（stage B），这一阶段的特点是，整体下沉速度急剧加大，伴随着沉降漏斗的边缘水平扩展，这一时期漏斗边缘基本成型，stage B的发展速度很快，持续时间通常也比较短；红色曲线属于后活跃阶段（stage C），该阶段的特点是，漏斗中央沉降速度依然很大，但是漏斗边缘扩展的速度急剧下降，漏斗的中央逐渐发生超临界沉降，随着开采的推进，发展成为盆地，盆地中心的最大沉降值保持不变，沉降中心底部逐渐变平，stage C的持续时间通常受开采深度、开采方式、上覆岩层厚度及坚硬程度等影响；蓝色曲线属于稳定阶段（stage D），这个阶段整个沉降漏斗速度都低于1.67 mm/d，由于在前一阶段发生了塌陷，应力快速释放并形成新的平衡，漏斗内的活动速率迅速降低。塌陷绝大多数都发生在stage C，后活跃阶段采空区上覆关键层全部断裂，压力平衡拱达到极限状态，此时意外的外力因素极易导致应力失稳形成塌陷。陕西的煤矿埋藏深度普遍较浅，上覆黄土和岩层，岩层一般较薄。

图2 采空区沉降模型图

如果我们假设的沉降模型合理，那么使用重访周期12 d的Sentinel-1A卫星C波段SAR数据对各个沉降阶段的地面进行干涉处理，滤波和去平之后，可能得到的干涉的结果如下：开始阶段由于整体下沉速率低，在12 d累计最大沉降不超过20 mm，干涉图上应无明显现象。前活跃阶段，漏斗中央速度迅速增大，超越C波的相位梯度，而漏斗边缘也在水平和垂直向有所运动，但是相位梯度相对中央平缓许多，在漏斗外边缘和内边缘之间可能形成一个不超越C波段形变梯度的沉降环，所以这一阶段得到的干涉图可能会在沉降漏斗外边缘附近形成若干干涉条纹，而条纹内由于沉降速度超越形变梯度阈值，导致干涉条纹发生混叠。后活跃阶段，漏斗中央沉降速度依然很快，但漏斗边缘速度却急剧下降，所以中央由于梯度过大无法形成干涉，而边缘由于速度太低也无法形成干涉条纹。稳定阶段整体沉降速度都低于1.67 mm/d，也无法形成干涉条纹。也就是说，如果此沉降模型符合实际，那么使用重访周期12 d的C波段SAR影像来不断对整个沉降过程的地面做干涉的话，应该会在第2阶段前活跃阶段观测到沉降漏斗边缘的活跃，表现在干涉图像上应该是在沉降中心的边缘地带形成若干环状的干涉条纹。陕西的矿区大部分位于黄土覆盖区，原涛[17]通过相似材料模拟实验和FLAC数值模拟揭示了黄土区开采沉陷的基本规律，认为厚黄土层矿区地表移动具有发展快、稳定快、活跃期短的特征，且黄土层厚度是影响地表动态移动规律的重要因素。因此在黄土覆盖的薄岩层矿区第3阶段的持续时间较短，观测到边缘的环状干涉条带距离塌陷发生时间也相应较短。

3 数据处理与验证

在研究中我们选择了最近几年2次震级相对较大的塌陷地震事件，2次地震分别发生在府谷县三道沟附近和彬州市新民镇东务村附近，府谷和彬州都拥有丰富的煤矿资源，也长期受到采空区地面沉降的困扰。使用短周期时序InSAR对2次地震事件进行分析，验证提出的模型是否符合预期。

3.1 2018年8月16日府谷ML3.2塌陷地震

据榆林地震台网测定2018年8月16日13:09分在陕西榆林市府谷县三道沟镇（39.16°N，110.73°E）发生ML3.2地震（地震类型塌陷）（图3）。我们尝试用重访周期12 d的C波段升轨哨兵卫星数据，对发生塌陷的区域做干涉处理，影像从2018年1月1日到2018年8月29日，包含了共21景SAR数据及精密轨道数据，获取以12 d为间隔的20对干涉图像，干涉后的影像进行滤波和去平处理。

图3 府谷塌陷地震位置示意图（底图来源于Google earth）

处理结果：

从处理后获得的干涉序列（图4）上来看，2018年1月1日—7月24日间的17幅干涉图属于沉降的stageA阶段，2018年7月24日—8月5日的干涉图出现了非常明显的环状干涉条纹，表示沉降进入stageB，沉降速度明显加大且沉降范围迅速扩大。塌陷发生于2018年8月16日，说明在该区域stageC时间跨度较短，长度大约为10～20 d，在塌陷发生后，沉降迅速进入stageD，沉降速率快速下降，在干涉图上无明显干涉现象。出现在stageB的干涉环最大半径达到2 km。

图4 府谷塌陷地震雷达干涉序列图

SAR影像基本参数见表1，参考DEM为SRTM3 DEM数据，地面分辨率90 m。

表1 哨兵2A影像基本参数列表

3.2 2017年9月5日彬州市ML3.2级塌陷地震

2017年9月5日01:14:22，据中国地震台网测定，陕西省彬州市（35.09°N，108.18°E）发生 ML3.2 地震（塌陷），震源深度0 km（图5）。选择从2017年4月2日—10月31日的15幅哨兵升轨数据，干涉处理后获得14张间隔周期为12 d的干涉影像，滤波和去平处理后如图5所示。塌陷发生于2017年9月5日，说明在彬州市的这次塌陷stageC时间跨度也较短，为大约5～15 d。

图5 彬州塌陷地震位置示意图（底图来源于Google earth）

处理结果：

从干涉的结果（图6）来看，2017年4月22日—9月25日之间的8幅干涉影像都属于stage A，2017年8月20日—9月1日的干涉影像上出现明显环状干涉，表示沉降进入stage B，塌陷发生之后，迅速进入stage D，干涉图上表现的沉降活动微弱。从stage B的干涉环测量得到的沉降区域最大半径约为0.6 km，规模相对较小。

图6 彬州塌陷地震干涉序列图

SAR影像基本参数见表2，参考DEM为SRTM3 DEM数据，地面分辨率90 m。

表2 哨兵2A影像基本参数列表

4 讨论与分析

从2次塌陷地震的干涉时间序列来看，发生于黄土覆盖区的矿区塌陷地震，短周期的C波段雷达数据干涉均可以在塌陷地震发生前明显地识别出前活跃阶段（stage B）（图7）。2次塌陷地震从stage B到地震发生之间的时间间隔较短，且地震发生后沉降活动迅速减弱，印证了原涛[17]关于厚黄土层矿区地表移动形变具有发展快、稳定快、活跃期短的观点。影响stage C时间长短的因素有很多，主要可以分为自然地质因素和采矿技术因素[16]，自然地质因素如煤层厚度、煤层倾角、埋藏深度、松散层厚度、断层密度、地层不整合、岩石坚硬系数等，采矿技术因素如采空区位置、开采方法、工作面推进速度、顶板控制方法等。如果要准确地估算可能的塌陷地震与stage B的时间间隔，需要综合前面的自然地质因素与采矿技术因素来综合考量。通过短周期的时序InSAR干涉测量判断潜在塌陷地震区位置、所处阶段，结合现场考察或资料搜集等获得塌陷区的自然地质因素以及采矿技术因素，可以为塌陷地震的早期预警提供有力的支撑。

图7 采空区在阶段B的三维剖面示意图

2017年9月5日彬州ML3.2级塌陷地震在SAR的干涉影像上不如2018年8月16日府谷ML3.2级塌陷地震表现得明显，主要是塌陷范围小导致，从InSAR观测结果上分析，彬州塌陷的最大直径约600 m，而府谷塌陷的最大直径为2 000 m。如果使用波长较短的X波段SAR数据来研究彬州塌陷地震效果会更加明显。

需要注意的是，即使在短周期的干涉图像中发现明显的环状干涉区域，该区域在未来一段时间是否会发生塌陷地震以及什么时候发生地震，这仍然需要结合区域实际的地质地貌条件、采空区深度等多种因素来综合判断。考虑到哨兵雷达影像的分辨率不高，如果塌陷的范围过小，会影响到整体的监测效果。小范围的塌陷使用波长更短分辨率更高的雷达数据或许可以获得更佳的监测效果。

5 结论

结合前人对地表下沉速率与周期的研究以及大量现场监测得到的沉降数据，总结沉降发育的时空规律经验模型，配合C波段SAR数据的特性获得的短周期干涉序列，可以大范围、低成本地精确筛选出处于stage B的潜在塌陷区域，结合野外实际调查以及资料收集获得的采空区深度、上覆岩层厚度及坚硬系数、开采方法等信息，可以综合判断塌陷是否会发生以及可能发生在什么时间，从而为塌陷地震的早期预警提供支撑。

通过在黄土覆盖矿区开展塌陷地震预警研究工作，可为未来在更广泛的区域监测塌陷地震、做塌陷风险评估、提前预警等开展更多研究提供经验。

致谢感谢袁志祥研究员在论文写作过程中的悉心指导，给本文写作提供了建设性的意见，在此表示感谢。