引言

合成孔径雷达干涉技术（Interferometric Synthetic Aperture Radar，缩写为InSAR）是近30年来发展起来的一种空间对地观测技术，在形变监测方面具有极大的应用潜力。然而，该技术目前仍存在诸多问题未能解决，其中SAR数据干涉对的失相干问题即为其最大障碍之一（刘晓萌等，2007）。失相干会导致失相干区域不能形成干涉图，使得处理出来的同震形变场的部分区域因数据缺失而出现空白（Zebker，Villasenor，1992；王超等，2002）。然而，由于目前的SAR技术尚无法从根源上解决失相干问题，因此考虑从处理已出现失相干的结果入手来恢复失相干，插值法便是一种可行的途径。

2009年意大利拉奎拉M_W6.3地震发生距今已有十年，有关该地震的研究成果相继发表，例如基于DInSAR技术研究拉奎拉地震的同震形变场（杨久东等，2010；Wang et al，2010），基于PSInSAR技术进行该地震的时序分析（Ryder et al，2010；罗三明等，2012），以及多平台融合解算该地震的三维同震形变场（王永哲等，2012）。但是由于植被的影响，获取的拉奎拉地震同震形变场存在严重的失相干现象，对于该地震失相干现象的研究却未见报道。鉴于此，本文使用移动窗口克里金法对拉奎拉地震同震形变场进行插值处理，恢复形变图的失相干区域，并对其结果进行反演验证，以期说明该方法在恢复形变图失相干区域应用中的可行性。本文研究思路如图1所示。

图  1  本文研究思路流程图

Figure  1.  Flow chart of the research ideas in this study

下载: 全尺寸图片 幻灯片

1.   拉奎拉地震同震位移场测量

2009年4月6日，在意大利中部的拉奎拉城附近发生M_W6.3地震，地震震中的地理坐标为（42°25′N，13°23′E），震源深度为8.8 km，该地震造成了巨大的人员伤亡和财产损失（张健，彭军还，2011）。本文以拉奎拉城区为中心选择研究区，图2为研究区域地理位置的示意图，其中粗线框圈定了本次研究所选取的区域，细线框内区域为单景单视复数（single-look complex，缩写为SLC）影像的覆盖范围。

图  2  研究区域（粗框）示意图

Figure  2.  Schematic diagram of the studied area （thick frame）

下载: 全尺寸图片 幻灯片

1.1   数据

拉奎拉地震发生之后，欧洲空间局（European Space Agency，缩写为ESA）陆续共享了覆盖该地震震中区包括ENVISAT、ERS-2以及ALOS等在内的多种卫星观测资料，其中可用于同震处理的ASAR资料有3个轨道近50景的数据。本文采用其中的两景影像（表1），应用ENVI SARscape软件，结合90 m分辨率的SRTM DEM数据进行二轨差分干涉处理，并采用Goldstein滤波（Goldstein，Werner，1998）和最小费用流（minimum cost flow）（王秀萍，2010）解缠算法获取了此次地震的同震形变场。

表  1  ENVI ASAR雷达数据参数列表

Table  1.  Data parameters of ENVI ASAR radar

轨道	升降轨模式	接收日期	时间基线/d	空间垂基线/m
T079	降轨	2008-04-27 2009-04-12	350	40.401

下载: 导出CSV 

| 显示表格

1.2   DInSAR同震位移场

拉奎拉城市位于亚平宁山脉最高段大萨索山西麓，阿泰尔诺河谷地中，地形复杂，且植被较多，这使得所获InSAR同震像对的相干性较低。从所获干涉影像和解缠结果均可看出，干涉条纹有很大部分的缺失，说明InSAR像对存在很严重的失相干现象。

通常情况下，对滤波后的差分干涉图进行相位解缠，得到解缠后的差分干涉图。相位解缠的实质是获取干涉图中各干涉相位整周数之差，当干涉图的质量较好时，使用各种相位解缠方法均能获取令人满意的结果，而当干涉图的质量较差时，各解缠方法皆难以获得高精度解缠结果。由于拉奎拉地震的极震区存在严重的失相干现象，本文就最小费用流算法和3D解缠算法分别尝试对数据进行解缠处理，但解缠结果均不理想。

最终根据解缠后的差分干涉图，将获得的形变信息从雷达图像的斜距投影坐标系中转换至正射的地理坐标系中，进行地理编码，获取了拉奎拉地震降轨数据的同震形变场，如图3所示。可以看出，结果存在严重的失相干现象，这很可能是由于拉奎拉地区植被较多且处理结果受大气影响较大所致。由于失相干现象的存在导致同震形变场大部分区域均出现数据缺失，使得处理出来的同震形变图出现大范围的影像空白，因而不能进行有效的形变分析。

图  3  拉奎拉地震同震形变场（a）及其部分放大区域（b）

Figure  3.  Coseismic deformation field of the L′Aquila earthquake （a） and part of its magnified area （b）

下载: 全尺寸图片 幻灯片

2.   移动窗口克里金插值法实现相干恢复

前已述及，获得的同震形变场由于失相干现象的存在导致震区的形变图出现部分数据缺失，难以对该地震进行有效的形变分析。针对已缺失数据的恢复，本文使用移动窗口克里金插值法（moving window Kriging method）对失相干区域进行插值，实现相干恢复（Harris et al，2010）。该方法的处理步骤如图4所示。

图  4  插值法处理形变场流程图

Figure  4.  Flow chart of interpolation process for deformation field

下载: 全尺寸图片 幻灯片

2.1   数据预处理

首先对所获得的同震形变图（图3）进行掩膜处理以去除形变异常值，然后对图像进行裁剪、降采样等处理以减少计算量，其中，行向、列向的重采样因子均设为0.5，采样后的像元数为216万9 222，满足软件的操作要求。降采样使用像素聚合（pixel aggregate）法，通过对所有输出像元值有贡献的像元进行加权平均来采样。经过上述处理后获得结果形变图（图5），之后对形变图进行彩色渲染显示，显示为白色的区域均为空值区，即失相干区，可以看出拉奎拉地震存在严重的失相干现象。

图  5  经过裁剪、降采样处理所获得的结果形变图

Figure  5.  Resulting deformation image obtained after the cropping and downsampling processing

下载: 全尺寸图片 幻灯片

2.2   移动窗口克里金法实现同震形变场失相干恢复

克里金插值法（Kriging relax）是一种对预测点最优无偏估计的空间插值方法，根据研究区域已知控制点的空间分布构建半方差函数并求取预测点的权系数，并基于此进行高程数据的内插。该方法是对预测点的线性无偏估计，即满足最小二乘原理，主要分析已知控制点的大小、形状及空间位置上的相互关系，实测点与预测点之间的空间位置关系，以及变异函数提供的结构信息。相对而言，地形变化多样，克里金插值法能够更好地反映这种变化（王艳妮等，2008）。

克里金技术是一种广泛使用于点样的光栅数据插值技术，如数字高程模型和地球化学图的生成（Harris et al，2010）。使用该技术所获结果的质量取决于采样值的空间分布和半方差模型的性质，这与样本数据集在未知位置预测值的经验全局函数相吻合。然而，这种半方差图模型可能不适用于空间分布上局部情况复杂多变的数据集，如在大地震中观测到的差分干涉合成孔径雷达（DInSAR）数据。为此，需引进移动窗口克里金法来解决那些空间分布上局部趋势复杂的数据集。

移动窗口克里金法的基本原理是根据较小的邻域重新计算变程、块金和偏基台半变异函数参数。为了避免使用全局半方差模型，首先引入移动窗口机制，在每个像素位置的窗口区域内导出局部半方差图，将该窗口内所有参考点的距离加权函数作为中心像素位置的克里金插值的一个最小二乘解；然后基于窗口的克里金插值过程在每个像素位置重复，直到整个形变图均被插值。本文应用Surfer软件，使用移动窗口克里金法对上述经过裁剪、降采样预处理的形变图（图5）进行插值处理（王子龙，2011），插值结果的优劣主要取决于插值参数的设置，参数主要包括滤波参数、重采样后像元大小和移动搜索窗口大小等。进行参数设置时，移动搜索的扫描窗口要足够大，以确保有足够多的高质量数据点（参考点）进行插值。通过不断调整插值参数进行试验，最终得到较为合理的插值形变图，将结果展示在ArcGIS中（叶宝莹等，2009），如图6所示。

图  6  使用移动窗口克里金法插值后的形变图

Figure  6.  Interpolated deformation map by moving window Kriging

下载: 全尺寸图片 幻灯片

将插值前（图5）、后（图6）的形变图进行对比，很明显可见形变场中的失相干区域得到了很好的恢复，表明移动窗口克里金法能够恢复InSAR像对失相干区域的数据缺失问题，达到了对有数据缺失影像的相干恢复的目的。从插值后形变图（图6）可以看出，极震区的上下盘形变图像已较完整，上、下盘最大视线向上的位移分别约为−27 cm和8 cm，呈椭圆状沿NW−SE方向展布。

3.   拉奎拉地震同震形变场插值结果验证

反演是基于同震表面观测地震数据，将已知的地质定律和数据作为约束条件来模拟地下断层的空间结构和物理性质，其目的是找到模拟数据与地球表面观测数据之间所对应的差值，当所对应的差值最小时即对应于反演模型的参数解。基本步骤为：先假设一个均匀滑动的分布源，再设定所研究地震的断层几何参数，然后根据地表观测值与断层滑动量呈线性相关，通过得到的线性关系反演求解断层位错量。实际上，观测数据量和观测数据的误差是有限的，所以得到的反演结果并不唯一，故求得震源参数后，需正演求得模拟形变值，再与真实形变值相比较，以评定反演结果的精度及可靠性，同时将反演得到的同震形变场与上节中插值后的形变图作比较，来判定插值方法是否可靠。常用方法包括基于梯度下降法（supervised descent method，缩写为SDM）反演和PSCMP程序正演。

3.1   形变场重采样

InSAR数据具有高分辨率的特点，观测数据量巨大，而用于反演的观测数据量充足即可，没必要将所有的InSAR数据作为约束进行反演（冯万鹏，2010；王永哲等，2012；王乐洋等，2017），因此反演时需要对形变场数据进行合理的降采样处理。由于此次地震引起的地表形变梯度较小，影响范围不大，但是东南角区域的形变值异常，再者近场数据对断层参数的影响较远场数据更大，故使用均匀选点法（刘唐志，郭小宏，2003），同时结合局部特征进行重采样选点，重采样后最终用于反演的地面观测点总共2 919个（图7）。

图  7  选点图

Figure  7.  Point-choosing map

下载: 全尺寸图片 幻灯片

3.2   梯度下降法反演处理流程

采用梯度下降法对同震位移场进行反演，获得目标区的地球物理模型（刘智荣等，2012）。本文选取均匀介质的弹性半空间位错模型，介质的泊松比设为0.25，在确定断层的几何参数之后，基于梯度下降法（Wang et al，2013）反演拉奎拉地震的断层参数，模型拟合度可达95%，反演参数结果列于表2，反演得到的模拟形变结果如图8a所示。由表2和图8a可见：断层走向结果非常稳定，约为141.8°，断层倾向西南；该断层未出露地表，上边界埋深约为2.9 km。断层倾角约为55°。所有反演模型的模拟残差基本控制在1.2 cm左右，保持了较低的残差水平（图8b），说明了反演结果的可靠性。

表  2  拉奎拉地震断层模型几何参数列表

Table  2.  List of geometric parameters of L′Aquila seismic fault model

采样点个数	断层参数								MW	误差/cm
	东经/（°）	北纬/（°）	走向/°	倾角/°	断层长度/km	顶深/km	底深/km	滑动角/°
2 919	13.470	42.367	141.8°	55°	13.2	2.9	12.6	−118.7	6.3	1.2

下载: 导出CSV 

| 显示表格

图  8  模拟形变（a）及其残差（b）

Figure  8.  Simulated deformation （a） and its residual （b）

下载: 全尺寸图片 幻灯片

反演本身不会对失相干的部分进行处理，但是在通过使用基于梯度下降法（SDM）反演获得地球物理模型后，假定地表形变与地下破裂所产生的位移存在线性关系，便可根据反演所获的震源模型，通过正演获取地表的形变量（Wang et al，2003，2006）。由这个过程获得的形变场，就是完整的形变场。图8a即是利用基于梯度下降法反演得到表2中的模型参数，而后通过PSSMP正演模拟所获的同震位移场。将图8a与图6的结果作残差计算，得到残差图，如图9所示，可见插值结果与反演结果在极震区吻合得很好，说明通过插值法得到的同震形变结果具有很高的可信度。

图  9  形变残差图

Figure  9.  Deformation residual map

下载: 全尺寸图片 幻灯片

4.   讨论与结论

本文利用二轨法对2009年意大利拉奎拉M_W6.3地震前后的ENVISAT ASAR数据进行干涉处理，成功获得了此次地震的差分干涉图，进而估计了同震形变场，但是由于植被、大气等失相干源的存在导致形变图出现了失相干现象。利用两景ENVISAT ASAR数据，使用移动窗口克里金插值法实现了研究拉奎拉地震极震区的失相干重构，并通过正反演模拟的形变场结果对插值结果进行了验证。所得结论如下：

1） 拉奎拉地震极震区的上下盘形变图像已较完整，上、下盘最大视线向上的位移分别约为−27 cm和8 cm，呈椭圆状沿NW−SE方向展布。

2） 拉奎拉地震的断层走向结果非常稳定，约为141.8°，断层倾向西南。该断层未出露地表，上边界埋深约为2.9 km，断层倾角约为55°。

3） 反演结果与插值结果吻合得很好，验证了移动窗口克里金插值法恢复失相干的可靠性，说明插值法可以成为实现地震同震形变场的失相干恢复的一种途径。