采动区地表异常形变SBAS识别与损害原因揭示

刁鑫鹏，孙全帅，白志辉，吴 侃，杨 静

（1.中国矿业大学 环境与测绘学院，江苏 徐州 221116；2.中国矿业大学 江苏省资源环境信息工程重点实验室，江苏 徐州 221116；3.冀中能源峰峰集团有限公司，河北 邯郸 056000；4.宁夏回族自治区地球物理地球化学调查院，宁夏 银川 750000）

0 引 言

煤炭是我国的基础能源，长期以来为经济社会发展和国家能源的安全稳定供应提供了有力保障[1-3]。煤炭资源大规模开采的同时，对矿区环境亦造成了极大损害，衍生出一系列的地质灾害与社会问题。减轻或防治各类灾害发生的关键在于明确开采沉陷规律，当前对沉陷规律的研究已取得了一系列显著成果[4-8]。然而，因复杂地形及地质构造的存在，使得地表移动和变形分布的正常规律被破坏，地表异常变形现象持续出现，且因此而引起的地表损害较常规情况更为严重[9-13]。

目前，国内外针对地表采动异常损害及其变形规律的研究正逐渐增多，亦取得部分成果[14-16]。如余学义等[17]从矿区地质采矿条件入手，揭示了张坡村出现大量地裂缝和建筑物异常损害的原因为开采诱发的断层滑移；郭迅等[18]分析了建筑物差异性沉降的根本原因是构造和采煤沉陷相互作用下局部区域岩体的“多米诺骨牌”反应；杨隽等[19]、DIAO等[20]多角度分析了昔阳县某村房屋建筑异常损害的原因为采动引起的坡体滑移；DONNELLY等[21]研究发现已稳沉区域的部分断层仍存在“活化”现象，并会引起地表的严重损害，但因缺乏实测数据，未能确定异常形变的持续时间。准确认识和理解地表移动过程的前提是对研究对象进行高精度、高可靠性的监测；然而，地表异常形变及损害相对而言具有特殊性和不可预知性，容易被忽视。同时，由于传统测量手段存在着测点稀疏、测量范围小、监测周期长等不足[22]，使得多维度、长时间序列的地表实测数据相对缺乏，导致目前对地表采动异常形变时、空演变特征的认识尚不充分。随着SAR卫星数据的不断丰富及数据处理手段的发展，InSAR技术获取地表形变的精度已达毫米级；特别是，该技术可利用研究区的存档影像对历史形变进行反演，成为研究采动区地表异常形变规律和揭示损害原因的有效手段。

以河北省邯郸市某城镇的异常损害为研究对象，利用2019-07-14—2020-04-27期间的25景Sentinel-1A SAR影像数据，通过短基线集技术（SBAS-InSAR）反演了异常损害区地表形变的空间分布形态及其演变过程。结合研究区地质资料，揭示了前述损害是开采诱发F15断层“活化”的结果，并得出地表变形存在跨断层传递的结论。研究成果表明，SBAS能够实现对地表异常形变的有效识别和提取，可为异常损害规律研究和损害原因揭示提供数据基础，对沉陷控制理论的发展完善亦具有借鉴意义。

1 研究区及异常损害现象

1.1 研究区概况

研究区位于河北省邯郸市，地理位置处于114°08′～114°11′E与36°36′～36°38′N，图1展示了研究区域的井上下对照情况。其中，地面主要为农田及城镇建筑物；井下工作面开采运用综采工艺，顶板管理方式为全部垮落法；煤层顶、底板岩性以粉砂岩为主。表1列出了研究区周边所开采工作面的具体信息。

表1 研究区周边各工作面开采信息

图1 研究区井上下对照情况Fig.1 Surface-underground contrast plan of study area

1.2 损害特征及沉陷预计分析

1）损害特征。该城镇于2019年10月开始出现损害现象，且随时间推移损害加重。经调查，至2020年5月，建筑墙体裂缝达30 mm，地表张开型裂缝及台阶裂缝宽度达10～20 mm。受损害区域呈条带分布，具体分布位置如图1中红色曲线所示。

2）沉陷预计分析。为确定常规情况下工作面开采对该镇的影响，采用概率积分法进行地表沉陷预计。预计参数基于该矿地表移动观测站实测数据并参照相关规范指南确定，具体为：下沉系数q=0.78，主要影响角正切tanβ=1.70，开采影响传播角θ=90°-0.4α，水平移动系数b=0.30，拐点偏移距S=0。图2显示了基于概率积分法预计的地表下沉等值线。

经预计计算得到，该城镇范围内地表下沉量、倾斜和水平变形的最大值分别为：200 mm、3.0 mm/m和2.5 mm/m；井下工作面开采会造成该区域房屋建筑一定程度上的损害。但基于《建筑物、水体、铁路及主要井巷煤柱留设与压煤开采规范》中的变形控制指标，地下开采所致建筑物的损害程度不会超过II级，这与实际损害情况不符。此外，井下工作面开采的直接影响边界（图2中绿色曲线），并未波及至现场调研的严重损害区，两者间水平投影距离超130 m。表明该城镇范围内房屋建筑的损害属异常损害。

2 SBAS原理方法与数据准备

为进一步分析地表变形的演变过程，准确揭示异常损害的原因，利用InSAR技术可通过存档影像回溯地表历史形变的特性，采用SBAS时序分析方法对研究区建筑受损害前、后的地表变形情况进行反演。

2.1 原理与方法

SBAS技术的提出是为削弱时空失相干影响，提高InSAR地表形变解译精度。BERARDINO、LANARI等[23-24]先后探讨了该技术在地表形变监测中的应用，而后HOOPER等[25]又对处理过程中高相干点的选取方法进行了优化。

与常规差分干涉测量相比，SBAS技术的基本思想是：通过设置时、空基线条件，将同一研究区的多期SAR影像分组成若干个短基线集合；利用最小二乘准则，获取每个短基线集合的地表形变时间序列；通过奇异值分解方法再将分组后的短基线集进行联合求解；进而得到研究区整个观测时间段内时序的地表形变信息[26]。具体[27-28]可表述为：假定同一研究区成像时间为t0，t1，t2，…，tn-1的N幅SAR影像，在短基线距的条件下可形成M幅差分干涉图，且满足N/2≤M≤N（N-1）/2。对于干涉图i，在去除平地相位后，干涉图中任意位置（x，y）的干涉相位δφ（x,y）可表示为

δφi（x,y）=δφi（tB,x,y）-δφi（tA,x,y）≈

δφdef,i（x,y）+δφε,i（x,y）+δφα,i（x,y）+δφn,i（x,y）

（1）

其中，i为干涉图序号；tA、tB为干涉图i所对应SAR影像的获取时间；δφdef为tA～tB时间段内地表

视线向形变相位；δφε为目标区地形相位；δφα、δφn分别为大气相位和噪声相位。其中前3项可以表示为：

（2）

式中，di（tB,x,y）和di（tA,x,y）分别为tB、tA时刻地表沿雷达视线向的变形量；λ为雷达波长；R为斜距；B⊥为垂直基线；Δz为DEM高程差；θ为入射角。若不同干涉图间的形变速率为vk,k+1，则tA～tB间的累积形变量可表达为：

（3）

对M幅干涉图进行三维时空相位解缠和地理编码后，便可求出不同SAR影像获取时间的形变速率。

2.2 数据准备

用于研究区形变信息提取的SAR数据为2019-07-14—2020-04-27期间的25景Sentinel-1A影像，C波段，地面分辨率为5 m×20 m，相关参数如下：

成像时间 2019-07-14、2019-07-26、2019-08-07、2019-08-19

2019-08-31、2019-09-12、2019-09-24、2019-10-06

2019-10-18、2019-10-30、2019-11-11、2019-11-23

2019-12-05、2019-12-17、2019-12-29、2020-01-10

2020-01-22、2020-02-03、2020-02-15、2020-02-27

2020-03-10、2020-03-22、2020-04-03、2020-04-15

2020-04-27

数据类型 IW

像元大小（m×m） 2.3×14.1

Sentinel-1影像的轨道参数需根据影像的成像时间下载，欧空局提供了哨兵影像的Precise Orbit Ephemerides（POD精密定轨星历数据）。此外，为消除地形相位，外部DEM采用的是美国宇航局SRTM3数据。

3 形变解译与特征分析

3.1 InSAR形变解译

SBAS时序分析主要包括：差分干涉对生成、点目标选取、差分干涉图解缠，以及时间形变序列获取等几个步骤，详细的处理流程如图3所示。

时序分析过程中相关参数设置如下：①连接图生成时，为充分利用影像数据，时间基线设置为50 d，空间基线设置为临界基线的10%；共连接成90个像对，其中最长时间基线为48 d，最长空间基线为175.56 m（平均59.78 m）。②干涉处理过程中距离向与方位向的视数比设置为4∶1，影像配准利用精密轨道数据采用强度互相关算法，滤波方法选用自适应滤波法，解缠方法采用最小费用流法（相干性阈值设置为0.3）。③轨道精炼和重去平处理过程中，地面控制点选择相干性高、相位好的点，并且所选控制点位置的形变需为0。④第一次解译反演选用较为稳定的线性模型。其余参数均选用处理Sentinel-1数据的推荐参数。

图3 SBAS地表形变解译流程Fig.3 Flow of surface deformation with SBAS method

图4为最终解译的地表垂直位移的分布情况，显示了不同时间段内研究区地表变形的空间分布形态及其随时间的演变过程。

3.2 InSAR解译可靠性评定

为验证InSAR解译结果的可靠程度，收集了48个地面观测点在2020-01-12—2020-05-03期间的水准测量数据，水准点位的具体分布位置如图2所示。因水准测量时间与InSAR解译时间段并不完全一致，为降低因时间差异而造成的对比误差，用于参与对比的InSAR数据选用时间段最为接近的2020-01-10—2020-04-27期间的解译结果。图5为2类数据各监测点位的变形曲线及对比情况。

基于图5中地表水准实测数据，可分析得出：①对于Q测线，Q13点以东地表变形不明显，Q13与Q14点间地表变形存在突变现象，Q14点以西的地面沉降先减小后又逐渐增加。②对于P测线，P5点以东地表变形不明显，P5与P7点间地表变形存在突变现象，P7至P18测点地面沉降有减小趋势。③地表变形突变点位连线，与城镇严重受损区位置具有一致性。

经统计得到，对于上述48个监测点位，水准测量与InSAR解译的平均差值为+5.32 mm，均方根差为11.49 mm。两者监测结果存在差异的主要原因有：①观测时间段不同（水准测量较InSAR解译的时间间隔长），导致上述对比数据中水准测量变形值较InSAR解译结果偏大；②InSAR解译过程中，因控制点选取及大气效应去除等因素导致解译结果存在误差。

虽然2类数据的监测结果间存在一定偏差，但平均差值占最大下沉值的比例仅为5%；同时，2类监测数据的变形趋势具有一致性。因此，我们可以基于InSAR解译结果分析地表变形的时、空演变特征。

3.3 InSAR解译形变特征分析

图5中的水准实测曲线一定程度上展示了采动地表的变形情况，但因时间和空间分辨率低的缺陷，对于地表突变出现的时机、变形的动态演变过程和空间分布特征难以准确表达，而基于SBAS形变解译结果可对上述特征做进一步分析。

图4a显示，2019-07-14—2019-09-12期间，周边区域地下工作面的开采并未引起城镇相关位置出现明显的变形现象。图4b显示，至2019-09-24，受损区域地表的零星位置开始出现变形，量值约为10 mm。此时，618工作面开采已接近结束，158工作面已进入村庄煤柱进行开采，但建筑物尚未出现明显损害。表明周期性的InSAR形变监测可实现对地表异常损害的预警。

图4c—图4f显示，自2019-9-24日起（158工作面推进距离约130 m），城镇相关位置已经开始出现明显的地表变形现象，且地表变形发育的停止位置与镇内受损严重房屋的分布位置具有很好的一致性。图4g—图4i显示，在158、618、619工作面作面开采引起常规的地表沉陷以外，地表变形有向东扩展的趋势，且地表变形空间发育的终止位置近似成直线分布，而该终止位置恰为现场调研确定的城镇异常损害的位置。

自图4b显示2019-9-24受损区域地表的零星位置开始出现变形后，图4g、图4h及图4i进一步显示上述位置地表的变形量逐渐增大、变形范围逐渐增加，但与工作面开采引起的主要变形区域在空间上相对独立。图4k、图4l亦表明，地下工作面开采在引起地表的常规变形以外，在沉陷盆地以东存在较为明显的地表变形，且该异常形变与常规沉陷盆地在空间分布上具有较为明显的独立性。图4m及图4n则展示出，随地下开采范围的不断扩展，上述两独立的沉陷区域逐渐发育为一体，发育的终止位置城镇异常损害的分布具有一致性。

上述地表变形的空间分布形态及特征，一定程度上符合断层赋存环境下开采沉陷的规律。

4 损害原因揭示

综合地表实测数据与InSAR动态解译结果的分析情况，研究区地表异常形变与损害大概率与某断层的“活化”有关。此外，现场调研发现损害区地裂缝两侧存在较为明显落差（图1b），符合断层露头处地表损坏特征。

井上下对照情况（图2）显示，研究区周边发育有3条较为明显的断层，分别为：F49（H=140 m，∠65°）、F13（H=100 m，∠70°）和F15（H>600 m，∠70°），图6为8号勘探线的地质剖面图。

图6 勘探线地质剖面Fig.6 Geological section of exploration line

由图6可得，地下工作面的开采容易扰动到F13断层，引起断层的“活化”；以45°作为松散层移动角所划定的F13断层露头位置并未出现在城镇屋建筑的受损区域，距离损害区约306 m。若F13断层露头出现在建筑物受损区域，则松散层移动角约为26°，这与常规情况下松散层移动角的取值范围不符；即工作面开采引起F13断层“活化”，而直接导致建筑物受损的可能性较低。

当以56°作为松散层移动角时，F15断层的露头位于城镇的异常受损区，且移动角的选择满足常规取值要求；此外，F15断层落差大于600 m，属特大型断层，相对而言更易受到扰动。因此，结合前节地表实测数据及InSAR形变反演的分析结果，认为F15断层受到扰动而发生“活化”是研究区城镇建筑发生异常损害的主要原因。基于上述结论，可进一步得到断层赋存环境下，采动影响范围会远超常规沉陷理论划定的影响边界，地表变形存在跨断层传递现象。

5 结 论

1）现场调研及水准实测显示，研究区受损害位置及程度超出常规沉陷理论划定的采动影响边界和级别，地表变形存在突变现象，建（构）筑物损害属异常损害。

2）通过InSAR解释结果与水准数据的对比得到，两者监测地表变形的趋势具有一致性，监测变形的平均差值为5.32 mm，证明SBAS形变解译结果具有可靠性，时序分析技术能够有效识别并提取采动区地表的异常变形。

3）InSAR解译地表变形空间分布的最终形态特征，一定程度上符合断层存在时开采沉陷的规律；但从地表变形的动态发育过程看，异常形变初始发育位置与常规沉陷盆地在空间上相对独立，而后才逐渐发育为一体；进一步揭示了断层受开采扰动的过程。

4）地表异常损害为开采诱发F15断层“活化”的结果，断层赋存环境下采动影响范围会远超以移动角划定的影响边界，且地表变形存在跨断层传递现象。研究成果对沉陷控制理论的发展完善具有借鉴意义。