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 研究区周边各工作面开采信息
Table 1 Mining information of working faces around study area
序号工作面开采起止日期开采厚度m/m开采深度h/m煤层倾角α/（°）16102015-01—2017-024.0650～7001921582019-08—2020-074.0750～9002436192020-01—2020-104.0600～75017
图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时间段内地表
图2 概率积分法预计的研究区地表下沉等值线
Fig.2 Contour of surface subsidence in study area predicted by probability integral method
视线向形变相位；δφε为目标区地形相位；δφα、δφ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为最终解译的地表垂直位移的分布情况，显示了不同时间段内研究区地表变形的空间分布形态及其随时间的演变过程。
图4 SBAS解译研究区不同时间段内地表垂直位移分布情况
Fig.4 Distribution of surface vertical displacements in different time periods interpreted with SBAS method
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形变解译结果可对上述特征做进一步分析。
图5 监测点位InSAR与水准测量变形与对比曲线
Fig.5 Deformation and comparison curve of InSAR and leveling results at each observation point
图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断层“活化”的结果，断层赋存环境下采动影响范围会远超以移动角划定的影响边界，且地表变形存在跨断层传递现象。研究成果对沉陷控制理论的发展完善具有借鉴意义。
[1] 刘 峰,曹文君,张建明,等. 我国煤炭工业科技创新进展及“十四五”发展方向[J]. 煤炭学报,2021,46（1）:1-15.
LIU Feng, CAO Wenjun, ZHANG Jianming, et al. Current technological innovation and development direction of the 14th Five -Year Plan period in China coal industry[J]. Journal of China Coal Society,2021,46（1）:1-15.
[2] 易四海. 特厚急倾斜煤层水平分层开采地表沉陷规律研究[D]. 北京：中国矿业大学（北京）, 2009.
YI Sihai. The research on surface subsidence due to mining horizontal slice of extra-thick and steep coal seam[D]. Beijing: China University of Mining and Technology-Beijing,2009.
[3] 李春意. 覆岩与地表移动变形演化规律的预测理论及实验研究[D]. 北京：中国矿业大学（北京）, 2010.
LI Chunyi. Prediction theory and experiment research of overburden and surface movement deformation evolution regulation[D]. Beijing:China University of Mining and Technology-Beijing,2010.
[4] 胡青峰. 特厚煤层高效开采覆岩与地表移动规律及预测方法研究[D]. 北京：中国矿业大学（北京）,2011.
HU Qingfeng. Overburden and surface movement due to high efficient mining extra-thick coal seam and prediction method[D]. Beijing:China University of Mining and Technology-Beijing,2011.
[5] 胡海峰. 不同土岩比复合介质地表沉陷规律及预测研究[D]. 太原：太原理工大学, 2012.
HU Haifeng. Research on surface subsidence regularity and prediction under composite mediums with different thickness rations of loose and bedrock[D]. Taiyuan:Taiyuan University of Technology,2012.
[6] 邓喀中,谭志祥,姜 岩,等. 变形监测及沉陷工程学[M]. 徐州:中国矿业大学出版社,2014.
[7] 胡海峰，廉旭刚，蔡音飞，等. 山西黄土丘陵采煤沉陷区生态环境破坏与修复研究[J]. 煤炭科学技术, 2020, 48（4）:70-79.
HU Haifeng, LIAN Xugang, CAI Yinfei,et al. Study on ecological environment damage and restoration for coal mining-subsided[J]. Coal Science and Technology, 2020, 48（4）:70-79.
[8] 胡炳南,颜丙双. 充填采煤开采沉陷主控因素及其影响规律研究[J]. 煤矿开采,2016,21（2）:57-59.
HU Bingnan,YAN Bingshuang. Main control factors and influence rule of mining subsidence by backfill mining[J]. Coal Mining Technology,2016, 21（2）: 57-59.
[9] 张华兴. 倾斜煤层开采地表移动变形等价计算方法研究[J]. 煤炭科学技术, 2020, 48（10）:119-123.
ZHANG Huaxing. Study on equivalent calculation method of surface subsidence in inclined coal seam mining [J ].Coal Science and Technology，2020，48（10）: 119- 123.
[10] GUO W, ZHAO G, BAI E, et al. Effect of overburden bending deformation and alluvium mechanical parameters on surface subsidence due to longwall mining[J]. Bulletin of Engineering Geology and the Environment, 2021, 80（3）: 2751-2764.
[11] 崔希民,张 兵,彭 超. 建筑物采动损害评价研究现状与进展[J].煤炭学报,2015,40（8）:1718-1728.
CUI Ximin, ZHANG Bing, PENG Chao. Research status and progress in damage assessment of buildings due to mining [J]. Journal of China Coal Society, 2015,40（8）:1718-1728.
[12] 张华兴，仲惟林. 受断层影响的地表移动计算[J].煤炭学报, 1995, 20（2）: 163-166.
ZHANG Huaxing, ZHONG Weilin. Calculation of surface movement incurred by fault[J]. Journal of China Coal Society, 1995,20（2）:163-166.
[13] DIAO X P, BAI Z H, WU K, et al. Assessment of mining-induced damage to structures using InSAR time series analysis: a case study of Jiulong mine, China[J]. Environmental Earth Sciences, 2018, 77（5）:166.
[14] 郭文兵. 断层影响下地表裂缝发育范围及特征分析[J]. 矿业安全与环保,2000, 27（2）:25-27.
GUO Wenbing. Analysis on the development range and characteristics of surface cracks under the influence of faults[J]. Mining Safety & Environmental Protection,2000, 27（2）:25-27.
[15] 孔宪森. 峰峰矿区采动断层滑移数值模拟研究[D]. 北京：中国矿业大学, 2017.
KONG Xiansen. Numerical simulation study of the fault slip induced by coal mining in the Feng-Feng Mining Area[D]. Beijing: China University of Mining and Technology,2017.
[16] 徐乃忠,高 超. 正断层存在的地表沉陷特殊性规律研究[J]. 采矿与岩层控制工程学报,2020,2（1）:101-106.
XU Naizhong,GAO Chao. Study on the special law of surface subsidence in normal fault[J]. Journal of Mining And Strata Control Engineering,2020,2（1）:101-106.
[17] 余学义,黄国耀,赵兵朝,等. 断层对采动变形影响及村庄保护煤柱留设方法研究[C]//煤炭开采新理论与新技术——中国煤炭学会开采专业委员会2010年学术年会论文集, 2010.
YU Xueyi,HUANG Guoyao,ZHAO Bingchao,et al.Study on the influence of fault on mining deformation and the retaining method of village protective coal pillar[C]//New theory and technology of coal mining—Proceedings of 2010 academic annual meeting of Mining Professional Committee of China Coal Society,2010.
[18] 郭 讯,戴君武. 采煤沉陷与断层相互作用引起地表建筑破坏特点分析[J]. 辽宁工程技术大学学报,2006,25（6）:851-854.
GUO Xun,DAI Junwu. Analysis of building damage caused by interaction between faults and coal mining subsidence[J]. Journal of Liaoning Technical University（Natural Science）, 2006, 25（6）: 851-854.
[19] 杨 隽,吴 侃,刁鑫鹏,等. SBAS-InSAR形变特征监测在矿区损害鉴定分析中的应用[J]. 金属矿山,2020,51（6）:174-178.
YANG Jun, WU Kan, DIAO Xinpeng, et al. Application of SBAS-InSAR deformation characteristic monitoring in cause analysis of building damage[J]. Safety in Coal Mines,2020,51（6）:174-178.
[20] DIAO X, WU K, CHEN R, et al. Identifying the cause of abnormal building damage in mining subsidence areas using InSAR Technology[J]. IEEE Access, 2019, 7: 172296-172304.
[21] DONNELLY L J, CULSHAW M G, BELL F G. Longwall mining-induced fault reactivation and delayed subsidence ground movement in British coalfields[J]. Quarterly Journal of Engineering Geology and Hydrogeology, 2008, 41（3）: 301-314.
[22] DIAO X P, WU Kan Wu, ZHOU D W, et al. Integrating the probability integral method for subsidence prediction and differential synthetic aperture radar interferometry for monitoring mining subsidence in Fengfeng, China[J]. Journal of Applied Remote Sensing, 2016, 10（1）: 016028.
[23] BERARDINO P, FORNARO G, LANARI R, et al. A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms[J]. IEEE Transactions on Geoscience and Remote Sensing, 2002,40（11）:2375-2383.
[24] LANARI R, MORA O, MANUNTA M, et al. A small-baseline approach for investigating deformations on full-resolution differential SAR interferograms[J]. IEEE Transactions on Geoscience and Remote Sensing, 2004,42（7）:1377-1386.
[25] HOOPER A. A multi-temporal InSAR method incorporating both persistent scatterer and small baseline approaches[J]. Geophysical Research Letters, 2008,35: L16302.
[26] 武丽梅, 金丽华, 季碗栎. 结合SBAS技术的鲁西南地区地面沉降监测[J]. 测绘通报, 2019（1）:108-113.
WU Limei, JIN Lihua, JI Wanyue. Monitoring ground subsidence in southwestern Shandong based on SBAS technology [J]. Bulletin of Surveying and Mapping, 2019（1）:108-113.
[27] 何秀凤, 何 敏. InSAR对地观测数据处理方法与综合测量[M]. 北京: 科学技术出版社, 2012.
[28] 杨成生, 张 勤, 赵超英, 等. 短基线集InSAR技术用于大同盆地地面沉降、地裂缝及断裂活动监测[J]. 武汉大学学报（信息科学版）, 2014,39（8）:945-950.
YANG Chengsheng, ZHANG Qin, ZHAO Chaoying, et al. Small baseline subset InSAR technology used in Datong basin ground subsidence, fissure and fault zone monitoring[J]. Geomatics and Information Science of Wuhan University, 2014,39（8）:945-950.
