基于Sentinel-1影像的2019～2020年河北隆尧地裂缝及周边形变监测与分析

冯晓松 杨成生 赵超英 牛玉芬 高 涵

1 长安大学地质工程与测绘学院，西安市雁塔路126号，710054 2 河北工程大学矿业与测绘工程学院，河北省邯郸市太极路19号，056038 3 邯郸市自然资源空间信息重点实验室，河北省邯郸市太极路19号，056038 4 云南省地震局信息中心，昆明市茨坝蓝桉路2号，650201

地裂缝是一种由内部和外部地质动力过程所引起的浅表破裂[1]，通常表现为土壤地表层在各种人为因素(采矿、浇灌、抽水等)[2-4]及其他自然因素(构造断裂[4]、地震[5])影响下，土层的连续性被力的相互作用打断，并以一种重大地裂缝和地表裂缝破坏现象呈现出来[6]。地裂缝直接或间接破坏了人类和野生动物的生存环境[7]，并且地裂缝周边地区的地表往往会出现沉降现象[8]。

隆尧地裂缝位于华北平原隆尧县内，形成于1966年邢台地震[9]，震后裂缝进入非活跃期，2003年裂缝再次暴露[10]，2006年以来裂缝规模及活跃度迅速增加[11]。现今，隆尧地裂缝是华北平原规模最大、灾害程度最为严重的构造地裂缝[12]，分布在近EW向的20多个村庄里[13]。地表裂缝最大宽度为20 cm，可见深度达数m[14]。Yang等[15]研究了2007～2011年间隆尧地裂缝的形变特征，并进行建模处理，但随着《全国地面沉降防治规划(2011～2020年)》、《京津冀平原地面沉降综合防治总体规划》等文件的发布及地下水控采方案和南水北调工程的实施，隆尧地裂缝的特征可能发生了改变，但目前缺乏对其现阶段活动性的研究。

有学者使用全球导航卫星系统GNSS、全站仪等手段对地裂缝进行了大量的形变监测研究，但传统测量技术成本昂贵且数据获取不连续，严重阻碍了对地裂缝的进一步研究[16]。干涉合成孔径雷达InSAR技术是一项基于地球物理学、信号学、计算机等学科发展而来的新型技术，被广泛应用于火山监测、地震监测、滑坡监测、地裂缝监测等地质灾害监测领域。本文采用时序InSAR技术探测隆尧地裂缝及其周边地区的形变，并进一步分析现阶段隆尧地裂缝的活动特征。

1 研究区域及实验数据

1.1 研究区概况

中国大陆位于印度、亚欧和太平洋等板块大量汇聚形成的冲击区和挤压区，华北平原长期遭受拉张应力[17]，因此1 000多条裂缝遍布华北平原[18]。隆尧地裂缝位于华北平原的邢台市，处于华北地块的核心位置，拥有3个地质块：临清地堑、邢-衡隆起、宁晋-衡水断层，形成了一个“三叉戟”的形状。该区域分布有NE向和EW向断层，其中主要裂隙是太行山山前断裂和巨鹿EW向断裂。

隆尧地裂缝是华北平原最长、最活跃的地裂缝，位于隆尧县城与魏家镇之间(图1 红色实线)，全长约36 km，地表裂缝最大宽度为20 cm，可见深度达数m[19]，沿途的道路、河流、农田、房屋等都被其切断，造成了巨大的经济损失。在隆尧地裂缝周边分布有太行山山前断裂、临城断裂、巨鹿西部断裂、巨鹿东部断裂、鸡泽断裂和明华断裂，这些断裂与隆尧地裂缝相互作用，使得隆尧地裂缝周边地区的地表沉降尤为复杂。由于研究区域较为广阔，本文采用InSAR技术对河北隆尧地裂缝及其周边的地表形变进行监测。

图1 隆尧地裂缝及其周围地区地形Fig.1 Topographic of Longyao ground fissure and its surrounding area

1.2 研究数据

哨兵一号Sentinel-1采用的是C波段，其中心频率为5.405 GHz，C波段在X波段和S波段之间，兼具两者的特性。本文采用31景VV极化的Sentinel-1数据，所选时段为2019-01-03～2020-12-23，如表1所示。数据覆盖范围如图1所示，图中蓝色实线代表裁剪后的雷达影像范围，即本文的研究区域。本文使用的外部DEM数据是由美国地质勘探局USGS提供的SRTM DEM，分辨率为30 m。

2 研究方法

2.1 SBAS-InSAR

由于研究区域植被覆盖度较高，失相干现象较为严重，为了通过不同时间基线干涉图提供的多余观测来降低失相干的影响，并去除干涉图中的各种误差(如大气误差、DEM残差等)，本文首先采用SBAS-InSAR技术获得隆尧地裂缝及周边地区的时间形变序列。SBAS方法是在利用小基线干涉图集进行形变测量的前提下，合成符合条件的小基线干涉图对，然后基于形变速率的最小范数准则，通过SVD方法获得系数矩阵的伪逆，最终获得相干目标的形变速率及时间序列[20]。计算公式如下：

(1)

表1 实验数据

2.2 Stacking InSAR

为获取研究区域含有较小误差的形变速率，使用Stacking InSAR技术对解缠后的差分干涉图进行处理。Stacking InSAR技术假设大气在时间上的扰动是随机的，且地表形变为线性形变，因此可以将所有符合运算条件、经过解缠的干涉图进行线性叠加，从而减少信号穿过对流层、电离层所产生的大气误差对监测结果的影响，提高监测精度[21]。公式如下：

(2)

式中，Vi为第i幅干涉图的相位变化速率，φi为第i幅干涉图的解缠相位，ΔTi为第i幅干涉图的成像时间间隔。

3 数据处理及结果分析

3.1 数据处理

数据处理流程如下：1)首先对原始数据进行预处理，转化为SLC格式，然后进行基线估计，空间基线设置为200 m，时间基线设置为100 d，时空基线如图2所示；2)将得到的干涉对进行共轭相乘，得到差分干涉图，并对干涉图进行滤波以提高影像的信噪比；3)通过最小费用流法对滤波后的干涉图进行相位解缠，得到地表像元的真实相位；4)分别通过SBAS-InSAR技术和Stacking InSAR技术对研究区域的数据进行时序处理。经过上述数据处理得到研究区域的形变时间序列和形变速率，其中形变速率如图3所示，2019-01-03～2020-12-23的累积形变量如图4所示。

图2 时空基线Fig.2 Spatial-temporal baseline

图3 2019～2020年隆尧地裂缝周边地区形变速率Fig.3 Deformation rate of Longyao ground fissure from 2019 to 2020

图4 2019～2020年隆尧地裂缝周边地区累积沉降量Fig.4 Cumulative subsidence around Longyao groundfissure from 2019 to 2020

3.2 隆尧地裂缝及其周边地区形变特征分析

3.2.1 隆尧地裂缝

由图3和4可见，隆尧地裂缝的南北两侧有一定的形变速率差异，为了更加直观形象地观察隆尧地裂缝南北两侧的速率差，对隆尧地裂缝南北两侧进行剖线分析。为与Yang等[15]2007～2011年的ALOS PALSAR结果进行比较，本文选择相近位置的剖线(图3黑色虚线EF)。隆尧地区形变速率剖线所得结果如图5所示，图中浅灰色实线为隆尧地裂缝，实线左侧为地裂缝北侧、右侧为地裂缝南侧。由图可见，相对于隆尧地裂缝南侧，隆尧地裂缝北侧呈抬升趋势，二者的形变速率差达4 cm/a。

图5 沿剖线EF隆尧地裂缝南北两侧形变速率差异Fig.5 Deformation rate of both sides of Longyao ground fissure along profile EF

将本文采用Sentinel-1卫星影像获取的2019～2020年隆尧地裂缝南北侧形变速率剖面，与采用日本地球观测卫星合成孔径雷达ALOS PALSAR-1获取的2007～2011年形变速率剖面[14]进行对比分析发现，前者的隆尧地裂缝南侧沉降速率比后者的相对沉降速率约低1 cm/a。搜集资料后发现，国务院在2012年发布了《全国地面沉降防治规划(2011～2020年)》，随后自然资源部又组织编制了《京津冀平原地面沉降综合防治总体规划》。随着地下水控采措施和南水北调工程的实施，京津冀地区的沉降状况在2016年以后得到了大幅度的缓解，沉降范围和沉降中心形变量级的降低可能影响了隆尧地裂缝的活动性。此外，由剖线图可以明显看出，隆尧断层两侧形变特征存在明显差异，这一现象与Yang等[15]的研究中L波段获取的形变特征类似。

3.2.2 地裂缝周边地区

由图3和4可见，除巨鹿县县城所在位置有轻微的抬升外，境内其他地区均处于沉降状态。其中，固城店镇、魏家镇、官庄镇是明显的沉降中心，3个区域年沉降速率均达6～8 cm/a。对固城店镇、魏家庄镇、官庄镇的形变时间序列进行统计分析(图6)可知，固城店镇、魏家庄镇及官庄镇处于周期沉降状态，即当年2～9月处于快速沉降状态，当年9月至次年2月又存在一定程度的抬升，但其总体形变趋势仍是沉降。

图6 固城店镇、魏家庄镇及官庄镇单点形变时序Fig.6 Time series subsidence of Guchengdian, Weijiazhuang and Guanzhuang

调查发现，固城店镇、魏家庄镇、官庄镇都种植了大面积的农作物，包括冬小麦、玉米、棉花、蔬菜等，每年需要抽取大量的地下水进行灌溉。由于河北独特的气候条件，当地的主要农作物(小麦、玉米)为一年两收模式，冬小麦通常在当年10月播种并进行初次灌溉，次年开春时进行第2次灌溉。由此推测，当地对地下水的大量抽取及周边断层的影响造成了2～10月的沉降状态，之后随着灌溉抽水用量的急剧减少及冬季降雪和周围地下水回流，区域表现为缓慢的地表抬升[22]。此外，隆尧地裂缝东南端有一个较大的形变区域(图3和4中区域D)，查阅资料发现，当地于2003-11成立了隆尧工业园，主要生产方便面和饮料，其工业过程消耗大量的地下水，可能是造成当地地面沉降的主要原因。

4 隆尧地裂缝模型构建

研究表明，弹性均匀空间中的蠕变位错[6,23-24]可用来解释活动断层的大地观测。一个活动断层的地表变形可以被建模为反正切，为断层滑动速率和锁定深度提供了估计方法。垂直走滑断层运动的一个简单模型是均质弹性半空间的，可以描述为:

(3)

式中，V(x)为沿断层垂直剖面估计的各点平行断层速度，x为离断层的距离，S为断层滑动速率，D为位错深度，α为静态偏移量。

为更好地了解隆尧地裂缝的地面变形，使用上述模型对InSAR形变剖面进行建模。通过计算得知，地裂缝深度为5 m，几乎破裂至地表，滑移速率为27 mm/a，相比于2007～2011年的48 mm/a，滑移速率减缓约21 mm/a，说明现阶段隆尧地裂缝虽仍处于活跃状态，但随着对京津冀地区地下用水的宏观调控，地裂缝活跃度降低。

5 结 语

本文采用SBAS-InSAR技术，基于31景Sentinel-1卫星影像对邢台市隆尧地裂缝及其周边区域进行时间序列分析和物理模型构建。结果表明，随着对京津冀地区地下水使用的宏观调控，隆尧地裂缝及其周边地区形变趋势有所减缓，活跃度降低。具体结论如下：

1)现阶段隆尧地裂缝两侧相对形变有所降低，相较于2007～2011年降低约1 cm/a；

2)隆尧地裂缝周边地区如固城店镇、魏家庄镇、官庄镇是明显的沉降中心，3个区域年沉降速率均达6～8 cm/a，且呈现周期性的形变特征；

3)物理模型构建结果表明，隆尧地裂缝几乎破裂至地表，滑移速率约为27 mm/a，相较于2007～2011年减缓约21 mm/a。