北京时间2016-02-06 03:57:27，台湾省高雄市发生M_S6.7地震。震中位于高雄市美浓区，地理位置为22.94°N、120.43°E，震源深度15 km。地震使地壳积累的能量得到释放，同时引起地表形变。监测地震造成的地表形变对于研究地质灾害触发机制和评估受灾情况十分重要^[1]。过去数10 a，学者们使用D-InSAR(differential synthetic aperture radar interferometry)技术对地震引起的地表形变监测开展研究^[2]。现有的研究成果表明，D-InSAR技术是一种监测地表形变强有力的工具，尤其是在缺乏GPS数据的区域。

新一代SAR卫星Sentnel-1A在地震形变监测以及灾害应急响应中具有巨大优势^[3-4]。众所周知，常规D-InSAR技术只能获取雷达视线方向(line of sight, LOS)的一维形变，为了反演地表三维形变，通常是结合升降轨雷达数据来获取不同视线方向的形变。但是，由于目前所有的SAR卫星都是以近极地轨道方式运行的，雷达视线方向对南北方向的形变不敏感，使获取的南北向形变质量很差^[5]。为了能够改善南北向的形变质量，在升降轨InSAR的基础上，许多学者开始研究获取方位向形变的方法，其中比较常用的方法有多孔径干涉(multi-aperture interferometry, MAI)技术^[6]和偏移量跟踪(offset tracking)方法^[7]。偏移量跟踪技术监测到的方位向形变最多只能达到SAR数据像素大小的1/15~1/30^[8]。而对于多孔径干涉技术而言，采用Sentinel-1A数据的理论形变监测精度仅为27 cm^[9]。因此，对于小量级的地震形变，这两种方法在获取南北向形变(量级较小)方面的能力存在不足。

Sentinel-1A干涉宽幅影像每个子带中相邻的burst在方位向有一定的重叠度。在这个重叠区内，相邻的burst分别从前视和后视两个不同的角度对地面同一目标点进行两次成像。这一成像特点满足沿轨干涉^[10]的条件，可以采用沿轨干涉技术来获取方位向形变。为了准确反演2016年高雄M_W6.7地震的同震三维形变场，本文采用覆盖地震区域的Sentinel-1A升轨(T069)和降轨(T105)干涉宽幅数据(图 1)，提出联合沿轨干涉方法以获取方位向的形变，用D-InSAR技术获取雷达LOS向的形变，采用加权最小二乘方法反演2016年高雄M6.7地震的同震三维形变场, 进而对三维形变场进行分析和讨论。

1 研究区域及数据

图 1显示了本文的研究区域及所采用的Sentinel-1A数据的覆盖情况。采用的SAR数据是两对Sentinel-1A干涉宽幅成像模式影像，包括1对升轨数据和1对降轨数据。每个干涉对由一幅震前和一幅震后影像组成，见表 1。由图可见，采用的SAR数据完全覆盖了地震震中20 km半径内的区域。初步判断，SAR数据可以全面覆盖形变影响范围。

2 方法

为了获取2016年高雄地震的三维地表形变场，需要至少两个不同成像几何结构的影像，即一个升轨和一个降轨，并分别采用D-InSAR技术和沿轨干涉技术对SAR影像进行处理，以获得升轨和降轨在LOS向和方位向的形变量，进而反演地表三维形变场。

2.1 D-InSAR获取雷达LOS向形变

二轨D-InSAR由Massonnet等^[1]提出。先利用研究区域地表变化前后两幅SAR影像生成干涉条纹图φ_int，再获取SRTM DM数据模型地形相位的条纹图φ_{sim_topo}，然后从干涉条纹图中去除地形信息得到地表雷达视线方向形变信息d_los。其一般公式为：

(1)

式中，λ表示雷达波长。

LOS方向的形变标准差可以通过最大似然估计方法获得：

(2)

式中，N表示有效视数，γ表示干涉相位的相干性。

2.2 沿轨干涉技术获取方位向形变

Sentinel-1A是第一颗实现TOPS(terrain observation by progressive scans)为标准获取模式的SAR卫星。该模式从前视和后视两个不同角度对同一地面目标点进行两次成像，可以满足沿轨干涉的条件^[4]。对于Sentinel-1A而言，斜视角的差异在1°左右，这意味着重叠区域是burst长度的10%左右。

为了获得方位向的地表形变，可以对burst重叠区域采用沿轨干涉，该方法类似于多孔径干涉技术。通过以下步骤可以获取重叠区域差分相位的二次差分φ_d：①对主、从SAR影像进行配准，并提取前视和后视的burst重叠区域SAR影像；②利用前视主、从SAR影像生成前视几何结构的差分干涉图φ_forward；③利用后视主、从SAR影像生成后视几何结构的差分干涉图φ_back；④计算前、后视差分干涉图的二次相位差φ_d=φ_forward－φ_back。

burst重叠区域差分干涉图的二次差φ_d与方位向形变d_az的关系可以表示为^[11]：

(3)

式中，表示平行于卫星轨道的水平单位矢量，Δψ_d表示一个子带内相邻burst的斜视角的差。

burst重叠区域的沿轨干涉形变监测的标准差为^[11]：

(4)

式中，N表示有效视数，γ表示干涉相位的相干性，Δf表示重叠区的频谱间隔，Δt_s为方位向采样间隔。

2.3 三维形变场求解

理论上，反演地面三维形变需要3个以上不同方向的观测量。本文采用的是Sentinel-1A升轨和降轨2个干涉对数据，因此，可以通过D-InSAR技术获取升轨和降轨2个雷达LOS形变d_los^A、d_los^D，以及采用升降轨沿轨干涉技术获取2个方位向形变d_az^A、d_az^D，进而反演三维形变场d_n、d_e、d_u。具体理论推导可以参考文献[5]。三维形变的误差方程如下：

(5)

式中，误差向量，观测向量，参数向量 ，且方程系数矩阵为：

(6)

通过加权最小二乘方法可以求得参数X的解：

(7)

式中，Σ为观测向量L的先验方差矩阵。

3 实验结果及讨论 3.1 数据处理及结果

数据处理包括以下几个主要步骤。

1）D-InSAR处理。为了使SLC影像配准的精度达到1/1 000像元，我们采用文献[12]的配准方法。SLC影像配准之后，后续的D-InSAR数据处理与传统的Stripmap模式数据一致。采用精密轨道数据去除干涉图的平地相位，以及30 m分辨率的SRTM DEM模拟地形贡献相位并去除，然后进行相位解缠。最后通过差分相位到形变转化以及地理编码, 得到升降轨雷达视线方向的形变（图 2）。

2) 沿轨干涉处理。对burst重叠区域采用沿轨干涉方法获取方位向形变场。为了保证获取的方位向形变与D-InSAR结果大小一致，在burst重叠区的前视和后视干涉中，对距离向和方位向进行了相同视数(10 :2)的多视处理。相干性低于0.65的地方进行掩膜，丢弃不可靠的相位值。使用低通滤波方法(将1 km范围内的最大似然估计值作为每个像素的值)对干涉图进行滤波。经过去平地相位、去地形相位及相位解缠等步骤得到前视和后视差分干涉图，然后对前视和后视差分干涉图作差，得到干涉图的二次差分，根据式(3)可以求得burst重叠区的方位向形变。对burst非重叠区域，利用拉普拉斯算子进行内插，最后得到台湾高雄M_W6.7地震升降轨方位向的形变场(图 3)。

3) 三维形变反演。运用式(5)和式(6)，采用加权最小二乘法对其求解，最终得到高雄地震东西、南北和垂直方向的地表三维形变场(图 4)。

3.2 分析与讨论

图 2(a)和图 2(b)显示了升轨和降轨数据雷达LOS向的形变。可以看出，此次地震引起的地表形变主要集中在A、B两瓣区域。两升轨LOS向最大的正形变量达到+14 cm(图 2(a)中区域A)，降轨的为+10 cm(图 2(b)中区域A)。尽管形变数值存在差异，但升降轨雷达视线方向的最大形变位置比较一致，都是在震中西北方向约20 km，台南市东南方向约17 km处。而对于区域B，在升轨几何结构中反映出来的形变值接近于0，在降轨的几何结构中形变量达到-8 cm。结合SAR成像几何结构以及地表形变在升降轨LOS向上的数值差异推断，此次地震导致A瓣抬升、B瓣下沉，同时A、B瓣均往西移动。根据这一推断，可以进一步解释升降轨LOS向的形变差异：由于A瓣抬升，其在升轨几何LOS向表现为正的形变(即靠近卫星)，在降轨几何也表现为正的形变。同时A瓣往西移动，这一信号也会在升轨几何LOS向表现为正的形变，而其在降轨几何中表现为负的形变，因而在升轨几何中，A瓣的抬升和往西移动同时造成正形变(靠近卫星)，正的形变信号得到放大；而在降轨几何中，往西移动引起负的形变(远离卫星)，与抬升引起的正形变信号抵消。因此，A瓣在升轨几何中表现出更大的形变量，而在降轨几何中表现出来的形变量较小。同理，对于B瓣进行分析可以得到类似的结果。图 3显示了升降轨burst重叠区域方位向形变及插值之后的方位向形变结果，可以看出，方位向的形变量较小，在-2~2 mm之内。

为了说明拉普拉斯算子内插技术能够恢复非重叠区的形变场，设计一个实验进行验证。由于非重叠区缺少可靠的形变数据，这里选取burst重叠区作为研究对象。使用重叠区两端各10%的区域，借助拉普拉斯算子内插方法得到重叠区其余80%范围的形变场，然后将内插得到的形变场与沿轨干涉获取的形变场进行对比分析，最终得到的形变内插精度为2 mm左右，这说明拉普拉斯算子内插方法可以恢复形变场。

由于前视和后视的二次差分可以很大程度上抵消地形和对流层的影响，所以沿轨干涉相位质量优于D-InSAR，从而使其在获取小量级形变南北方向分量方面具有优势。而对于传统的条带模式，沿轨干涉无法使用。所以，对Sentinel-1A影像采用沿轨干涉方法，是目前获取小量级南北向形变的一种新手段。从图 4(d)中三维形变拟合的均方根误差可以清楚看到，在形变较大的区域其拟合误差较大。对比图 3获取的方位向形变及图 4的三维形变结果可以发现，形变大的区域没有落在burst重叠区，这些区域的小尺度方位向形变插值结果存在较大误差，且这些误差会通过三维形变拟合均方根误差反映出来。图 4(d)中的均方根误差代表了模型(7)的拟合误差。

联合升降轨LOS向和方位向形变观测反演的三维形变及拟合均方根误差列于图 4。由东西方向、南北方向以及垂直方向的形变可以看到，此次地震造成的地表形变以东西方向和垂直方向为主，A、B两瓣均发生往西的位移，同时A瓣抬升而B瓣下沉，南北方向有少量形变。在震中区域(B瓣)产生了大约4 cm的下沉，而在A瓣区域发生显著的抬升，最大抬升量达到12 cm；同时也发现，A瓣往西移动的位移最大值也比B瓣大，差异值达到4 cm。通过垂直方向形变可以初步推断，此次地震的发震断层位于A瓣和B瓣之间(图 4中的虚线位置)，且未出露地表，断层线左侧地表发生向上的抬升运动，而右侧地表发生向下的沉降运动，这一特征与逆冲断层的性质相符。由此可知，此次地震的发震断层具有逆冲性质，并可以判断断层线左侧为上盘、右侧为下盘，断层走向南南东。同震过程中，断层下盘向上盘方向推挤，使得下盘向西运动，而上盘向上逆冲，也有向西的运动。这一特征与本文所得的东西向形变结果非常吻合。图 3中沿轨干涉获取的方位向形变场显示有小的形变发生，表明断层具有一定的走滑特征。由上可知，此次地震的发震断层以逆冲为主兼具走滑特征。通过三维形变推断的断层性质与CMT利用地震波形数据反演所得的震源机制解^[13]符合得较好。

4 结语

本文采用Sentinel-1A数据监测2016年高雄M_S6.7地震同震形变，提出使用沿轨干涉技术获取方位向的形变，并利用D-InSAR技术获取升降轨LOS向形变。最后联合LOS向形变和方位向形变，通过加权最小二乘方法获取了此次地震同震形变的三维形变场，并对地震的发震断层性质进行分析。结果显示，高雄地震的发震断层为南南东走向，倾向西南，以逆冲为主兼具走滑性质。同震形变以垂直和东西方向为主，且形变同时分布在断层上盘和下盘。断层上盘抬升，最大12 cm；下盘下沉，最大8 cm；同时上下两盘均伴随着向西的位移，最大6 cm。通过联合LOS向形变和方位向形变解算2016年高雄地震的同震三维形变场结果表明，本文提出的沿轨干涉技术具有获取方位向形变的能力，Sentinel-A数据能够监测地震中小量级同震形变。