大而复杂形变区域中的精确三维形变反演——基于偏移的相位解缠和改进的多孔径SAR干涉测量集成技术在2016年熊本地震中的应用

Won-Kyung Baek, Hyung-Sup Jung*

Department of Geoinformatics, University of Seoul, Seoul 02504, Republic of Korea

1. 引言

常规合成孔径雷达（synthetic aperture radar, SAR）干涉测量（InSAR）功能强大，可以在大于1000 km2的大面积上测量精确的表面形变。该技术已成功应用于地震形变[1-6]、火山形变[7-11]、冰川运动[12-15]、地面沉降[16-19]、临时和季节性湿地水位变化[20,21]等。但是，由于InSAR方法只能观察到雷达视线（line-ofsight, LOS）方向上的一维（one-dimensional, 1D）形变，因此使用此方法几乎无法测量沿轨迹（along-track, AT）方向的形变。

然而，Bechor和Zebker [22]提出了多孔径SAR干涉技术（multiple-aperture SAR interferometry, MAI），Jung等[23-26]也进一步改进了该方法。MAI法可以精确测量沿AT方向的一维表面形变。该方法克服了InSAR方法的局限性。MAI处理包括以下步骤：①使用前视多普勒频谱生成前视差分干涉图；②使用后视多普勒频谱产生后视差分干涉图；③通过估计两个干涉图的相位差来创建MAI干涉图。众所周知，MAI性能要比方位角偏移跟踪性能高出[22-26]好几倍。

通过集成InSAR和MAI技术，可以用一组干涉对精确测量二维（two-dimensional, 2D）表面形变。而且，该方法可以使用升轨和降轨干涉对来观测三维（three-dimensional, 3D）形变。SuperSAR多方位SAR系统已通过精确3D形变测量可行性测试[27]。3D反演已广泛用于测量火山形变[8,9,28-32]、同震形变[4,5,33-38]、冰川运动[39,40]等。使用地中海盆地观测小卫星-SkyMed（Constellation of Small Satellites for the Mediterranean Basin Observation-SkyMed, COSMO-SkyMed）的X波段星座InSAR和MAI集成方法在东、北和垂直方向的3D测量精度分别约为0.86 cm、1.04 cm和0.55 cm [29]。但是，由于MAI相关损失和InSAR相位展开误差，3D形变反演在大而复杂形变区域中的应用仍然具有挑战性。

最近，Jung等[26]引入了前视和后视残差干涉图以减少相位噪声水平。结果显示，该方法极大地增强了MAI相关性，特别是对于由于大而复杂的形变而失去相关性的区域。Baek等[41]引入了多核偏移跟踪方法[13,15,42]，以最大限度地减小高陡形变区域的相位展开误差。他们指出，使用多内核偏移跟踪方法估算的ALOS-2 PALSAR-2偏移图与用于最小化相位展开误差的偏移图一样精确。在大而复杂形变区域中，InSAR测量不可用，但是偏移图是有效的，因此可以使用多核偏移跟踪方法来测量LOS形变。即使在由于大而复杂的表面形变而失去相关性的区域中，改进的InSAR和MAI测量的集成也可以精确绘制3D形变。但是，这种集成方法从未应用于生成3D形变图。

在这项研究中，我们发现即使在由于大而复杂的表面形变而失去相关性的区域中，基于偏移相位展开的InSAR [41]和改进的MAI方法的集成[26]也可以观测到3D形变场形变。通过集成改进的InSAR和MAI方法我们观测到了2016年的熊本地震3D形变场。2016年4月14日，日本熊本地震发生，东北—西南走向的走滑断层沿地表形变。在断层线附近，与地震有关的地表形变很大且很复杂。因此，尚未有研究尝试使用InSAR和MAI集成来进行精确的3D表面形变反演。本研究获取并处理了两个升轨对和一个降轨ALOS-2 PALSAR-2干涉对，以反演2016年熊本地震3D形变场。我们使用了11个全球定位系统（global positioning system, GPS）站的观测来评估通过改进的InSAR和MAI集成来进行3D形变测量的效果。

2. 研究领域和数据

2016年4月14日，日本熊本县发生了两次震级分别为6.2级和6.0级的地震。在主震之后，共观察到140多场余震。随后，大约28 h后发生7.0级的主震。据报道，2016年熊本地震导致了东北—西南走向的走滑断裂发生大而复杂的形变[38,43,44]。在断层带附近观测到约2.1 m的沉降，测得的隆起约为0.3 m，还观测到约2.0 m的东向形变[40]。因此与2016年熊本地震有关的地表形变既大又复杂，尤其是在断裂带附近。

本研究使用了两组升轨对和一组降轨ALOS-2 PALSAR-2同震干涉对。一组升轨对是在2015年11月19日和2016年6月16日（20151119_20160616）获得的，另一组是在2016年2月11日和2016年6月2日（20160211_20160602）获得的。降轨对是在2016年3月7日和2016年4月18日获得的（20160307_20160418）。这三组的垂直基线分别短至30 m、-74 m和-121 m，时间基线分别为210 d、140 d和42 d。表1总结了用于本研究的同震干涉测量对的干涉测量参数。图1是日本熊本市的地貌晕渲图，其中包括两组升轨对和一组降轨干涉对的边界。如图1所示，用实心三角形标记的GPS站测量值是用于评估所测量的3D形变的精度。

3. 方法

由于InSAR处理中的相位解缠误差和MAI处理中的相关性损耗，在大而复杂的形变区域中，3D形变反演很难进行。基于偏移的相位缠绕方法已经成功地应用于ALOS-2 PALSAR-2干涉图[41]。因此，有研究表明通过参考文献[41]中提出的改进的基于偏移的相位解缠的InSAR方法，即使在大而复杂的形变区域中，人们也可以观测到InSAR测量的LOS形变。与此同时，自Bechor和Zebker [22]首次提出MAI处理方法以来，该方法就一直被用于降低MAI相关性损失。最近，人们通过残差的前视和后视干涉图计算MAI干涉图，对MAI方法进行了改进[26]。如参考文献[26]中所述，改进的MAI方法可以有效降低MAI相关性损失。因此，即使在大而复杂的形变区域中，也可以通过MAI精确测量AT形变[26]。这些技术改进使得InSAR和MAI集成可以应用于大而复杂3D形变的精确测量。图2展示了改进的InSAR和MAI集成用于3D形变反演的详细工作流程。Baek等[41]和Jung等[26]改进的InSAR和MAI方法可以用于生成精确的2D形变图。如图2所示，给定升轨和降轨干涉对就可以从升轨和降轨2D形变图精确地创建3D形变图。

表1 从升轨和降轨上的ALOS-2 PALSAR-2数据中获取的同震干涉对

3.1. 改进的合成孔径雷达干涉测量方法简介

改进的InSAR方法包括5个主要步骤：①缠绕差分干涉图的创建；②距离偏移图像的计算与滤波；③通过从解缠的差分干涉图中减去滤波后的距离偏移图像来生成残差干涉图；④解缠的残差干涉图生成；⑤通过将解缠后的残差相位与滤波后的距离偏移量相加来创建解缠后的差分干涉图。有关该过程的更多详细信息请参见参考文献[41]。在此InSAR方法中，有效的残差干涉图的生成是关键。通常情况下，有效的残差干涉图很难生成，因为偏移方法的精度远低于InSAR方法的精度。但是，使用ALOS-2 PALSAR-2干涉对以及对该干涉对[41,42]使用改进的多核偏移量追踪法可以克服这一局限性。由于更高的空间分辨率和经过改进的偏移跟踪测量[13,41,42]，L波段ALOS-2 PALSAR-2干涉对的偏移跟踪性能很高。即使在大而复杂的形变区域中，这种InSAR方法也可以精确测量LOS形变。

图1. 日本熊本市的RGB图（前哨2），其中包括两个升轨对和一个降轨干涉对的边界。灰线表示断裂，黄色五角星表示2016年熊本地震的震中。白色和黑色三角形表示GPS站的位置。其中，白色三角形用于验证3D测量。标签“Asc1”和“Asc2”分别表示20151119_20160616和20160211_20160602的干涉对，“Dsc”是指20160307_20160418干涉对。

图2. 3D形变反演的详细工作流程。SLC：单视复数。

3.2. 改进的合成孔径雷达干涉测量方法的简介

改进的MAI方法主要包括五个步骤：①分别使用主图像和副图像的全孔径带宽及前视和后视孔径带宽，通过带通滤波创建三个单视复数（single-look complex,SLC）图像；②从三个SLC主图像和副图像生成三个差分干涉图；③通过从前视和后视差分相位中减去全孔径差分相位来创建两个残差干涉图；④通过计算两个剩余干涉图之间的相位差来创建一个MAI干涉图；⑤校正残差相位，包括地形变化和平地效应相位校正，以及MAI干涉图的自适应滤波。在第一步中，三个SLC图像必须具有相同的图像位置，这可以通过应用相同的距离像元徒动校正来实现。在第二步中，通过使用全孔径干涉对估算共配准参数，然后将其应用于三组干涉对。在第三步中，应使用诸如Goldstein滤波器之类的自适应滤波器对全孔径差分干涉图进行硬滤波。可以通过迭代应用具有大窗口内核的滤波器来执行硬滤波。最好在第四步进行相位差计算之前对两个残差干涉图进行轻微滤波。有关处理的更多详细信息请参见参考文献[26]。在MAI处理中，从两个残留干涉图生成MAI干涉图是关键。重要的是，应小心地生成硬滤波全孔径差分干涉图，因为它是无噪声的干涉图。为了生成更好的硬滤波干涉图，可以将8×8或16×16像素的小内核大小迭代地应用于自适应滤波器。即使在大而复杂的形变区域中，这种MAI方法也可以精确测量AT形变。

3.3. 3D形变反演

通过组合升轨和降轨2D 形变测量值来反演3D形变。使用参考文献[3,8,43]提供的InSAR升轨和降轨形变以及MAI升轨和降轨形变来定义形变矢量（r），如式（1）所示：

式中，d是3D表面形变矢量。U可以定义如下：

式中，uInSAR和uMAI分别是InSAR和MAI测量的单位向量。InSAR向量单位定义如下：

式中，[ ]T表示向量的转置，θ和φ分别是雷达的入射角和方位角，单位MAI向量可以定义如下：

最后，可以通过最小二乘解如下生成3D表面形变，其中，∑是InSAR和MAI测量[8]的协方差矩阵：

InSAR测量值既包含水平形变分量，也包含垂直形变分量，但包含较少的北向形变分量。但是，东向和垂直方向的MAI测量值几乎没有关系。另外，从上升轨道和下降轨道获得的MAI几何形状是相反但相似的。因此，如果两个MAI测量值中有一个不好，则不使用其反演北向形变。因此，为了最好地区分东、北和垂直形变，总测量值应大于三个，其中包括至少两个来自上升轨道和下降轨道的InSAR测量值和一个MAI测量值，如等式（2）所示。

4. 结果

使用InSAR和MAI集成技术的三维形变反演被应用于2016年熊本地震，目的是为了测试在大而复杂的形变区域通过改进的InSAR和MAI集成技术精确测量三维地表形变的可行性。众所周知，2016年熊本地震是沿一条右旋走滑断层发生的。据报道，此次地震相关形变在垂直方向上约为2.1 m，在水平方向上约为2.0 m [40,44]。

我们获得了两组升轨对和一组降轨同震干涉对来分析2016年熊本地震（表1）。通过如下7个处理步骤可以生成三个缠绕差分干涉图：①方位共带滤波；②偏差参数估计；③副SLC图重采样；④干涉图的生成；⑤使用航天飞机雷达地形探测任务数字高程模型生成合成干涉图；⑥缠绕差分干涉图的生成以及在方位向和距离向上使用15×12像素外观的多视；⑦内核大小为32的缠绕差分干涉图滤波。图3展示了升轨对20151119_20160616和20160211_20160602、降轨对20160307_20160418的缠绕差分干涉图。如图3所示，由于存在较大的同震形变，InSAR相位梯度较大。梯度最大的形变出现在断裂带附近（图3）。

这些区域的干涉相位不能被很好地观测到，因为在陡峭和复杂的形变中相位无法解缠。为了克服这一缺陷，我们采用多核偏移量追踪法对三对图像进行处理。使用从32×32到256×256像素的16种不同内核尺寸，总共得到16个偏移量测量值。使用内核尺寸为3×3×16像素的三维中值滤波估计每个像素点的最终偏移量。将最终距离偏移图由像素转换为弧度，然后用内核尺寸为11×11像素的非局部（non-local, NL）均值滤波器使之平滑，以生成缠绕残差干涉图。因此，利用常规最小费用流（minimum cost flow, MCF）算法，可以方便、准确地对缠绕残差干涉图进行解缠[45]。

图4为升轨对20151119_20160616、20160211_20160602和降轨对20160307_20160418的解缠差分干涉图。值得注意的是，图4中的一条条纹与图3中的不同。图4中的这个条纹对应4π。从图4可以看出，在断层附近的高条纹率区域，干涉相位可以正确地解缠。即使在断层线附近没有任何条纹图的区域，也发现了正确解缠的相位。在这些区域，干涉相位没有条纹图，而距离偏移量具有有效的测量。这是由于该地区形变大而复杂，而去相关系数低。因此，由于较低的去相关性，尽管它不会产生任何条纹图，测量的干涉相位仍是有效的。结果表明，改进的InSAR方法能较好地测量断裂带附近的大而复杂的形变。测量结果与GPS衍生的LOS形变匹配良好。对于Asc1、Asc2和Dsc对，最终干涉图的精度约为2.88 cm、1.96 cm和1.90 cm [41,44]。因此，我们认为额外的电离层校正并不是必要的。

为了生成MAI干涉图，对ALOS-2 PALSAR-2 SLC图像进行方位向傅里叶变换，利用汉明窗将多普勒频谱分割为前视和后视带宽。在此步骤中，由于通过使用窗口函数应用了多普勒频谱，我们首先删除应用的窗口，然后对多普勒频谱应用分带滤波。然后，将全孔径差分干涉图应用于内核尺寸分别为256×256、128×128、64×64像素的硬滤波器，生成残差干涉图。假设硬滤波差分干涉图为无噪声干涉图，我们可以通过在硬滤波干涉图与前视和后视干涉图之间进行相位相减，得到前视和后视残差干涉图。用Goldstein滤波器对前视和后视残差干涉图进行轻微滤波后，再由两幅残差干涉图的复共轭得到MAI干涉图。由于MAI干涉图具有地形相和地平相，我们采用参考文献[23]中提出的方法对相位进行校正。一些MAI干涉图存在电离层效应，这些误差应被消去；因此，我们将定向滤波方法应用到MAI干涉图中。更多关于定向滤波的细节可以在参考文献[13,15,46,47]中找到。

图3. 由20151119_20160616 (a)、20160211_20160602 (b)、20160307_20160418 (c)干涉对生成的缠绕InSAR干涉图。1 rad=180°/π。

图4. 使用改进的InSAR处理器从20151119_20160616 (a)、20160211_20160602 (b)和20160307_20160418 (c)干涉对中创建的展开InSAR干涉图。白框显示了3D测量的覆盖范围。

图5为20151119_20160616、20160211_20160602升轨对和20160307_20160418降轨对的MAI干涉图。当把图5中的MAI干涉图与图4中的InSAR干涉图进行对比时，MAI干涉图的空间分辨率要低得多，因为它们是通过子孔径处理生成的。值得注意的是，图5中的一个条纹对应的是0.2π。从图5可以看出，在断裂带附近的高条纹率区域以及在断裂带附近没有条纹图的区域，MAI相位的测量都是正确的。这是由于前视和后视残差干涉图因低去相关性而具有有效的测量值。结果表明，用改进的MAI方法可以较好地测量断裂带附近的大型复杂形变。形变时两组升轨测量精度约为8.13 cm和9.87 cm[图 5（a）、（b）]，降轨测量精度约为3.36 cm [图5（c）]。前两次测量的较大误差是由于升轨的MAI干涉图包含了严重的电离层失真，如图6（a）、（b）所示。为了削弱如图6（a）所示的电离层失真，我们在电离层条纹横向和纵向迭代应用旋转角度为45°、窗口内核为151×63像素的定向中值滤波器。此外，对图6（b）所示的MAI干涉图迭代应用旋转角度为50°、窗口内核为751×63像素的定向中值滤波器。有关电离层失真削弱的更多细节见参考文献[13,46]。图5（a）、（b）显示了通过定向中值滤波削弱电离层失真的效果。通过这种削弱方法，我们将20151119_20160616和20160211_20160602升轨对的精度从约52.29 cm和47.55 cm分别提高到约8.13 cm和9.87 cm。但形变时升轨的精度远低于降轨的精度。因此，在三维形变场的反演中不采用升轨AT形变法。

干涉去相关是评估InSAR和MAI干涉图测量精度的主要因素。特别是MAI的测量精度比InSAR对去相关性更敏感。因此，在使用MAI干涉图之前，需要考虑去相关因素[26]。空间的、时间的、热量的和体积的去相关性是众所周知的去相关组分[48]。相干性可以很好地描述相位去相关，通常用于评估理论误差水平[25,48]。利用相干性计算相位信号的空间稳定性时，其在高梯度形变区可能会被低估[41]。也就是说，相干性不能用来判断在一个大而复杂的形变区域内是否存在有效的形变。在大而复杂的形变区域，应更谨慎地进行分析。

图7显示了由降轨前视残差干涉图和相干值小于0.5的降轨MAI干涉图所估计的相干图。使用5×5像素的移动窗口计算相干值。大部分区域的相干值接近1.0；断层线附近的相干值为0.5 ～ 0.7。如图7（b）所示，低相干值主要分布在极陡形变区。在这些地区，算出的相干值不能用来评估MAI精度；因此，利用相干性测量这些区域的形变是否有效还不能确定。因此，偏移量追踪信息被另外用来确定一个有效的测量。

图5. 使用改进的MAI处理器创建的来自20151119_20160616 (a)、20160211_20160602 (b)和20160307_20160418 (c)干涉对的电离层校正的MAI干涉图。白框显示了3D测量的覆盖范围。

图6. 电离层校正前的MAI干涉图。（a）20151119_20160616；（b）20160211_20160602。

图7. （a）降轨前视残差干涉图的相干图；（b）在降轨MAI干涉图上以0.5为标准的相干阈值图。黑色像素表示相干值低于0.5的区域；ρ即相干值。

图8. 采用InSAR和MAI方法集成的2016年熊本地震地表三维形变场。（a）向东；（b）向北；（c）垂直形变。（a）～（c）上的彩色菱形表示来自GPS站的地表形变。

图8为采用改进的InSAR和MAI方法集成的2016年熊本地震地表三维形变场。我们使用了由两组升轨和一组降轨测量的LOS形变图和由降轨测量的AT形变图。如图8所示，即使是在大而复杂的形变中，三维形变场的反演效果也很好。正负最大形变向东约为1.78 m和-1.81 m [图8（a）]，向北约为1.57 m和-1.04 m [图8（b）]，垂直形变约为2.49 m和-0.56 m [图8（c）]。

图9为结合东向和北向形变计算得到的水平形变矢量场。如图9所示，矢量场的底图来源于垂直的形变。矢量场表明，2016年熊本地震是沿右旋走滑断层产生的。此外，断层西段观测到约2.49 m的隆起，东段测量到约0.56 m的下沉。这意味着2016年熊本地震具有右旋走滑和正断层运动的双重特征。

为了测试InSAR和MAI集成得到的地表三维形变场的精度，我们将SAR观测的形变与来自11个站点的GPS原位形变数据进行了对比（图10）。利用三次插值法提取11个站位的来自SAR的形变。SAR观测的形变和原位GPS形变的均方根误差（root-mean-square error,RMSE）在东、北、垂直三个方向分别约为2.96 cm、3.75 cm和2.86 cm。InSAR和MAI集成结果与原位GPS测量对北形变场的测量结果的一致性较差，而对东、垂直形变场的测量结果的一致性较好，因为用InSAR的LOS形变比用MAI的AT形变更精确。

图9. 结合东向和北向形变计算水平形变矢量场。矢量场的底图是垂直形变。

结果表明，在大而复杂的形变区域，在去相关系数较低的条件下，采用改进的InSAR和MAI集成方法可以获得精确的三维形变场。为了评估在大而复杂的形变区域内去相关性是否较低，可以使用偏移量追踪方法。如果能从该区域获得有效的偏移量，那么由于高相关性，我们可以获得有效的三维形变测量。否则，在该区域测量的三维形变是无效的。精确测量的3D形变可以帮助人们更好地理解地震和火山爆发等地质事件。

图10. 东（a）、北（b）、垂直（c）三个方向的SAR观测的形变与GPS原位形变的比较。

5. 结论

我们测试了通过集成改进的InSAR和MAI方法在大而复杂形变区域内获得精确三维形变测量的可行性。为此，我们采用集成方法对2016年熊本地震地表三维形变场进行了观测。2016年熊本地震在断裂带附近的形变大而复杂，因此目前还没有采用InSAR和MAI集成方法对其进行精确的地表三维形变反演。本研究使用两组升轨的和一组降轨的ALOS-2 PALSAR-2干涉对进行了三维形变反演。

对SLC干涉对采用常规InSAR处理生成三幅缠绕差分干涉图，并且采用多核偏移量追踪方法生成三幅距离偏移图。距离偏移量追踪法较好地测量了断层线附近大而复杂的形变场。这意味着断层线附近的去相关系数很低。因此，我们尝试使用基于偏移的解缠方法解缠三幅差分干涉图；解缠方法在三幅干涉图中得到了很好的应用。这表明利用基于偏移测量的解缠方法成功地观测到了断层附近的大而复杂的形变。利用前视和后视残差干涉图生成了三幅MAI干涉图，可以很好地测量断层线附近大而复杂的形变场。由于断层线的去相关性很低，两幅升轨MAI干涉图出现严重的电离层失真，但已经被校正。然而，经校正的升轨MAI干涉图的RMSE为降轨MAI干涉图的1/3。因此，没有使用两幅升轨MAI干涉图来反演三维形变。

利用三幅解缠的差分干涉图和降轨MAI干涉图反演2016年熊本地震地表三维形变场。三维形变图表明，即使在大而复杂的形变区域，三维形变场的反演效果也很好。从三维形变图上可以看出，2016年熊本地震具有右旋走滑和正断层运动的双重特征。通过将从SAR提取的结果与GPS原位形变对比，我们对三维形变场进行了精度评估。东、北、垂直三个方向的RMSE分别约为2.96 cm、3.75 cm和2.86 cm。这些结果表明，在大而复杂的形变区域中，在去相关性较低的条件下，采用改进的InSAR和MAI集成方法可以获得一个精确的三维形变场。精确测量的3D形变可以更好地理解地震和火山爆发等地质事件。

Acknowledgements

This study was funded by the Korea Meteorological Administration Research and Development Program(KMI2017-9060) and the National Research Foundation of Korea funded by the Korea government (NRF-2018M1A3A3A02066008). In addition, the ALOS-2 PALSAR-2 data used in this study are owned by the Japan Aerospace Exploration Agency (JAXA) and were provided through the JAXA’s ALOS-2 research program (RA4, PI No. 1412). The GPS data were provided by the Geospatial Information Authority of Japan.

Compliance with ethics guidelines

Won-Kyung Baek and Hyung-Sup Jung declare that they have no con fl ict of interest or fi nancial con fl icts to disclose.