地震形变监测可以提供断层的几何形状、地表形变以及释放的能量等信息,能够更好地了解地震灾害的具体情况。目前,有许多办法可以监测地震形变,使用最广泛的是全球导航卫星系统（global navigation satellite system,GNSS）和干涉合成孔径雷达（interferometric synthetic aperture radar,InSAR）^[1-5]技术。对于有地表破裂的大地震,GNSS技术只能在特定数量的位置进行测量,这会导致GNSS技术测量的地震形变覆盖面很小,空间分辨率较低,不能体现破裂位置附近的具体情况。InSAR技术对于较大地表形变,由于相邻像素之间的形变梯度超过一个相位周期会产生相位失相干^[6],无法获得破裂断层附近的形变。InSAR只能得到卫星视线向上的一维形变,它对沿卫星飞行轨道方向的形变十分不敏感。

Michel等^[7]通过欧洲遥感卫星(Europe remote‑sensing satellite,ERS）的合成孔径雷达（synthetic aperture radar,SAR）影像的亚像素互相关获取了Landers地震的二维形变场。Puymbroeck等^[8]通过SPOT光学影像获取了Landers地震水平形变场,并分析了形变结果的误差。自此以后,基于互相关技术的光学影像相关（optical image correlation,OIC）和SAR像素偏移追踪（pixel offset tracking,POT）被应用到了不同地震的形变监测中^[9-13],文献[9]使用Radarsat-1 SAR数据得到了2002年阿拉斯加Mw 7.9地震距离向形变,文献[10]使用ALOS-2 SAR数据得到了2016年熊本地震的距离向和方位向形变,文献‍[11]使用Landsat 7 ETM+光学影像获取了2001年昆仑山Ms 8.1地震的东西向和南北向形变,文献[12]使用Landsat 8光学影像提取了2013‍年巴基斯坦Mw 7.7地震的水平形变场,并对结果的误差来源和误差处理进行了系统分析。不同于GNSS和InSAR技术,它们对于发生较大地表形变的地震的地表形变监测有较强的适用性,SAR像素偏移追踪可以得到卫星距离向和方位向的形变场,光学影像相关可以得到东西向和南北向的形变场。

无论是SAR像素偏移追踪还是光学影像相关,均只能得到地表二维形变,地震的三维形变场可以观测断层的交叉,拐弯部分的复杂地表形变特征,为了能够探究地表形变最真实的情况,获取三维形变是很有必要的。目前三维形变已经有了较多的方法,其中主要包括：结合不同的入射角和方位角的多方向D-InSAR观测^[14],结合InSAR与多孔径InSAR（multi-aperture InSAR,MAI）^[15],结合InSAR和全球定位系统（global positioning system,GPS）^[16],结合InSAR与像素偏移追踪^[17],结合地应变特征与InSAR^[18],结合升降轨SAR影像像素偏移追踪^[19]等。

2016-11-13T11：02：00（UTC）,新西兰南岛东部凯库拉地区发生Mw 7.8地震,震中位置为42.757°S,173.077°E,震源深度15.1 km。在Mw ‍7.8主震发生的9 h内,发生了多次余震,其中有4次6级以上。新西兰位于太平洋板块和澳大利亚板块交界处,地质构造复杂,板块运动较为活跃,容易发生地震活动。发生本次地震的新西兰南岛东北部地区,太平洋板块以40 mm/a的速度相对于澳大利亚板块向西南偏西移动^[20],如图‍1所示。此次地震产生了大范围的破坏,地震后产生了多个地表破裂,从震中北部的Humps断层开始出现,沿着Hundalee、Jordan Thrust、Papatea和Kekerengu等主要断层向东北传播,破裂总长度超过150 km,带来的浅海破裂在凯库拉附近引发的海啸高度达到1.5 m^[21-22]。图1小图中红线为板块交界线,来自http：//peterbird.name/publications/2003_PB2002/2003_PB2002.htm；白框为哨兵1号与哨兵2号影像范围；黑线为活动断层,蓝线为本次地震的破裂迹线,来自新西兰数据中心（https：//data.gns.cri.nz/af/）；浅红色圆点以及蓝色沙滩球为余震,黑色沙滩球为主震,均来自美国地质调查局（United States Geological Survey,USGS）。

图  1  2016年新西兰凯库拉地震背景图

Figure  1.  Background of 2016 Kaikoura Earthquake, New Zealand

下载: 全尺寸图片 幻灯片

本次地震产生了巨大的地表形变,InSAR相位失相干情况十分严重^[21,23],许多方法都不适于用来获取此次地震的三维形变。对于凯库拉地震,文献[24]利用Planescape、Landsat 8和哨兵‍2‍号通过光学影像相关获取了东西向和南北向形变场。文献[6]将光学影像相关技术与光线追踪方法结合,并应用到Worldview光学影像中获取三维形变。文献[19]结合哨兵1号与ALOS-2像素偏移追踪获取三维形变,文献[23]结合ALOS-2升降轨像素偏移追踪所得到的4个方向形变获取三维形变。这些文献使用了ALOS-2或者高分辨率的光学影像,由于它们对于个人并不那么容易获得,成本较高,而哨兵1号和哨兵2号数据均可以免费获取,且并未有结合光学影像和SAR影像来获取凯库拉地震三维形变的相关研究,本文主要通过结合哨兵2号光学影像相关与哨兵1号像素偏移追踪获取本次地震的同震形变场,它们提供了不同方向的同震形变,这些形变能够组成不同的获取三维形变的组合。为找到一个最好的组合,本文考虑这些形变的不同组合,不加权与加权情况,通过最小二乘方法来解算凯库拉地震破裂地区附近的地表三维形变,预期的结果将有助于了解地震地表破裂区域的形变情况,根据同震三维形变结果对凯库拉地震的形变和地表破裂情况进行分析。

1 同震形变场获取

亚像素互相关技术是确定地震、冰川移动,滑坡和沙丘移动等引起的地表位移有效的技术之一^[25-27]。主要思想是通过将两个影像中的事件前影像作为参考,将事件后影像与其进行比较来评估两幅影像之间的相似性,进而检索同名点之间的位移。通过对事件前后的影像进行亚像素互相关,可以计算形变场^[28]。而光学影像相关与SAR像素偏移追踪正是基于亚像素互相关技术的方法。本文利用哨兵2号光学影像相关和哨兵1号像素偏移追踪来获取凯库拉地震形变场。

哨兵1号与哨兵2号是欧洲空间局（Europe‍an Space Agency, ESA）哥白尼计划下的地球观测任务中的卫星,哨兵1号载有C波段合成孔径雷达,拥有A和B两颗卫星,哨兵2号也是由两颗卫星组成,但其携带多光谱成像仪,它们对于气候变化和土地监测等方面的研究具有重大意义^[29]。本文使用的是来自阿拉斯加卫星设备（Alaska satellite facility, ASF）的哨兵1号数据和来自于USGS的哨兵2号数据如表1所示。

表  1  文中采用的影像信息表

Table  1.  Information of Images Used in This Study

卫星	震前影像采集时间	震后影像采集时间	时间间隔/d	轨道方向
Sentinel-1	2016-09-05	2016-11-16	67	降轨
Sentinel-1	2016-11-03	2016-11-15	12	升轨
Sentinel-2	2016-02-19	2017-02-23	370	-

下载: 导出CSV 

| 显示表格

1.1   光学影像相关

哨兵2号影像拥有4个10 m分辨率的波段,本文选取最适合用来进行光学影像相关的第8波段数据^[22],使用COSI-Corr^[30]来估计东西向和南北向的水平偏移量。并使用文献[22]和文献‍[31]所示方法来处理结果中所存在的主要误差。这种技术一般可以以1/20像素的理论精度获取地表位移^[32]。

如图2所示,形变场中可以清晰地看见几处明显的不连续性,这些不连续性产生的迹线就是地震产生的地表破裂在形变场中的表现,破裂区域与现场测得的地表破裂较为吻合,且与相关文献的哨兵2号光学影像相关结果相似^[22]。在凯库拉地震中,哨兵2号光学影像相关的结果有较为严重的轨道误差,在使用多项式拟合进行去除时,对Kekerengu左侧的形变范围产生了影响,并一定程度上削弱了该部分的形变量。由于形变场远场的形变接近于0,本文假设远场形变为0,取东西向与南北向形变场远场部分区域计算均方根误差,分别为0.361 0 m与0.263 8 m,符合光学影像相关技术的精度水平。

图  2  哨兵2号获得的同震形变

Figure  2.  Coseismic Deformation Obtained from Sentinel-2 Image

下载: 全尺寸图片 幻灯片

1.2   像素偏移追踪

使用像素偏移追踪的主要目的是得到距离向和方位向的形变场,本文使用哨兵1号数据,通过GAMMA软件来实现这一过程^[33]。哨兵1号升轨和降轨影像获得的同震形变如图3和图4所示,与已有研究相似^[34]。

图  3  哨兵1号升轨影像获得的同震形变

Figure  3.  Coseismic Deformation Obtained from Ascending Sentinel-1 Image

下载: 全尺寸图片 幻灯片

图  4  哨兵1号降轨影像获得的同震形变

Figure  4.  Coseismic Deformation Obtained from Descending Sentinel-1 Image

下载: 全尺寸图片 幻灯片

图3和图4中,正值表示靠近卫星的形变和沿卫星飞行方向的形变,负值表示相反方向的形变。距离向形变和方位向形变与地表三维形变有以下转换关系^[35]：

式中,$d_r$和$d_a$分别表示距离向形变和方位向形变；$d_{\mathrm{E}}$、$d_{\mathrm{N}}$、$d_{\mathrm{U}}$分别表示东西方向、南北方向和垂直方向的形变；$\alpha$表示卫星雷达波入射角；$\beta$表示卫星飞行轨道方位角。

与光学影像相关相同,像素偏移追踪的精度很大程度上也取决于SAR影像的分辨率,精度一般可达到1/10~1/30个像素^[36]。如图3和图4所示,距离向形变场中形变明显分布在几个主要破裂附近,噪声水平相对较低,而方位向形变场中形变则集中在东北地区,噪声水平明显高于距离向形变场,距离向形变场显示效果也优于方位向形变场。方位向形变场在西南地区没有明显的不连续性以及形变分布,而距离向形变场则有较大的形变分布,一部分原因是哨兵1号数据方位向分辨率较低。在西南地区距离向形变场中存在较大形变：由式（1）以及哨兵1号方位角可知,距离向形变对垂直向形变较为敏感；由式（2）可知,方位向形变与垂直向形变没有关联,对南北向形变较为敏感,而哨兵1号数据在西南地区检索出的南北向形变较小,原因可能是西南地区主要是垂直形变,而南北向形变较小。取形变场远场部分区域计算均方根误差,哨兵1号升轨数据获得的距离向与方位向形变场分别为0.113 6 m和0.445 9 m,哨兵1号降轨数据获得的距离向和方位向形变场分别为0.249 6 m和0.650 5 m。降轨数据得到的形变场精度和显示效果整体差于升轨,可能是由于降轨数据震前震后时间间隔过长导致配准误差较大,从而影响了像素偏移追踪结果的精度。

2 三维形变场的解算

2.1   方法

SAR像素偏移追踪可以得到距离向和方位向的形变,而升轨和降轨数据所得的距离向和方位向形变不同,相当于4个方向的形变,光学影像相关可以得到东西向和南北向的水平形变。将哨兵1号数据与哨兵2号数据组成能够解算出三维形变的不同的组合,见表2。以哨兵2号数据与哨兵1号升轨数据的组合为例,它表示哨兵2号数据光学影像相关所获取的东西向和南北向形变,以及哨兵1号升轨数据像素偏移追踪所获取的方位向和距离向形变。对于这5种组合,使用最小二乘法来求解,并对各方向形变分别进行不加权与加权操作。

表  2  哨兵1号与哨兵2号的数据组合

Table  2.  Data Combination of Sentinel-1 and Sentinel-2

组合	简称
哨兵2号数据与哨兵1号升轨数据	OIC+POT_As
哨兵2号数据与哨兵1号降轨数据	OIC+POT_Des
哨兵2号数据与哨兵1号升降轨数据	OIC+POT_As_Des
哨兵2号数据与哨兵1号升降轨数据获取的距离向形变	OIC+POT_Range
哨兵1号升降轨数据	POT_As_Des

下载: 导出CSV 

| 显示表格

以OIC+POT_As为例,各方向形变场与地表3个方向形变有以下关系：

式中,$d_{\mathrm{E}\mathrm{W}}$、$d_{\mathrm{N}\mathrm{S}}$、$d_{r\mathrm{A}\mathrm{s}}$和$d_{a\mathrm{A}\mathrm{s}}$分别表示为东西向形变、南北向形变、哨兵1号升轨数据得到距离向形变和方位向形变；$\mathit{W}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mathit{P}}_{\mathrm{E}\mathrm{W}},{\mathit{P}}_{\mathrm{N}\mathrm{S}},{\mathit{P}}_{r\mathrm{A}\mathrm{s}},{\mathit{P}}_{a\mathrm{A}\mathrm{s}})$表示权阵。在不加权情况下,$\mathit{P}$为单位阵；在加权情况下,每个权值代表每个形变场远场部分区域的均方误差（mean square error,MSE）,即P=1/M_MSE,可以写成：

2.2   结果

取凯库拉地震东北部主要破裂区域各方向形变场,由于光学数据结果与SAR数据结果有不同的分辨率,将它们降采样到相同的分辨率,然后根据文中所述方法计算三维形变,见图5和图‍6,第1行为东西向形变,第2行为南北向形变,第3行为垂直向形变。正值代表东向形变、北向形变和抬升；负值代表西向形变、南向形变和沉降。各组合最小二乘解算的3个方向形变场均方根误差（root mean square error,RMSE）见图7。

图  5  未加权最小二乘解算得到的各组合三维形变

Figure  5.  3D Deformations of Each Combination Obtained by the Least Square Solution Without Weight

下载: 全尺寸图片 幻灯片

图  6  各组合加权最小二乘解算得到的三维形变

Figure  6.  3D Deformation Obtained by the Weighted Least Squares Solution of Each Combination

下载: 全尺寸图片 幻灯片

图  7  各组合最小二乘解算的3个方向形变远场均方根误差

Figure  7.  Far Field Root Mean Square Error of Three Directions Deformation by Least Squares Solution of Each Combination

下载: 全尺寸图片 幻灯片

如图5和图7所示,OIC+POT_As和OIC+POT_Des得到垂直向形变与其他3个组合得到的垂直向形变有所差别,其他3个组合在Papatea断层左侧均存在较大抬升,与已有研究基本一致^[21,34],但OIC+POT_As在Kekerengu断层左侧有较大范围的沉降,而OIC+POT_Des在Kekerengu断层左侧有远超实际情况的抬升量和范围,这是因为仅依靠单个距离向形变场获取的垂直形变会受到组合中东西向形变的较大影响,而哨兵2号所得到的东西向形变在Kekerengu左侧低估了形变范围以及部分形变量,缺乏有效的约束。其他3个组合同时存在升轨数据和降轨数据获取的距离向形变,所得垂直向形变与实际情况较为相符^[21,34]。各组合得到的东西向形变有相对较好的一致性,由于OIC+POT_As和OIC+POT_Des组合所得的东西向形变很大程度上取决于组合中所包含的东西向形变,所以形变范围弱于其他3个同时包含哨兵1号升轨和降轨数据所获取的距离向形变的组合,其中OIC+POT_As_Des所得到的东西向形变精度较高。除OIC+POT_Range外的其他组合得到的南北向形变精度表现较差,这是因为精度较高的距离向形变对南北向形变不敏感,而精度较差的方位向形变对南北向形变影响较大,在有南北向形变约束的组合所得的南北向形变精度比没有南北向形变约束的POT_As_Des的南北向形变好,而仅仅由南北向形变决定的OIC+POT_Range得到的南北向形变精度取决于哨兵2号所获取的南北向形变场精度。从精度上看,包含方位向形变的组合所得到的南北向形变精度较差,由于本次地震哨兵2号所得的南北向形变低估了Kekerengu左侧的形变,故准确度优于不包含南北向形变的组合。

如图6和图7所示,对几个方向的形变场进行简单加权后（在包含方位向和南北向的组合中,根据升降轨所得方位向形变的精度增加一定的权重来计算的南北向形变）,OIC+POT_As_Des和OIC+POT_Range所得东西向形变在左上角有小部分明显的范围扩大,而其他组合并无明显变化。因为升轨哨兵1号数据所得方位向形变精度优于降轨数据所得,相对于等权时相对增加了这部分权重,而降低了降轨数据所得权重,所以提高了同时包含南北向形变和方位向形变的组合所得的南北向形变精度。垂直向形变在OIC+POT_As_Des以及OIC+POT_Range组合结果中左侧有明显的小部分沉降,这与文献‍[21]和文献[34]该部分相似,其他组合同样并未有明显变化。综上所述,在凯库拉地震中,OIC+POT_As_Des拥有最好的约束。由于哨兵‍2号光学影像相关的结果中低估了Kekerengu左侧的南北向形变范围和形变量,导致OIC+POT_Range组合数据约束性虽并未降低太多,而且有着相对较少的数据种类,但所得南北向形变的准确度不如OIC+POT_As_Des。而相对于POT_As_Des,它所得到的南北向形变精度更高,在Kekerengu断层右侧以及Papatea断层左侧的形变表现也更好。所以在哨兵2号光学影像相关和哨兵1号升降轨数据像素偏移追踪所得到的各类形变场组合中,OIC+POT_As_Des比较适合用来检测三维形变。

3 三维形变分析

对于本文所求的形变场范围,哨兵1号降轨影像的InSAR结果基本处于失相干的情况,而哨兵1号升轨影像的InSAR结果有部分未失相干,根据OIC+POT_As_Des的组合,并将哨兵1号升轨影像InSAR结果中未失相干部分替代哨兵1‍号升轨影像像素偏移追踪结果,计算整个凯库拉地震三维形变如图8所示,其中,图8（c）中黑色箭头表示水平形变,按照图8北部区域两个剖面线绘制对应的剖面图如图9所示。

图  8  凯库拉地震的三维形变

Figure  8.  3D Deformation of Kaikoura Earthquake

下载: 全尺寸图片 幻灯片

图  9  沿图8(b)剖面的形变

Figure  9.  Deformation Along Profiles of Fig.8(b)

下载: 全尺寸图片 幻灯片

选取所求形变场内的22个GNSS点计算均方根误差作为精度评定的指标,除去各形变场中几个噪声较大区域的GNSS点,其中东西向形变场有18个点,精度为0.56 m,主要在Kekerengu右侧以及Humps南侧的形变较大程度低于区域内3‍个GNSS点所测形变；南北向形变场有22个点,精度为0.36 ‍m；垂直向形变场有22个点,精度为0.45 m,对于本文所获取的近断层形变范围,以及哨兵2号光学影像相关和哨兵1号像素偏移追踪的精度水平较为合理。组合所计算的形变与18‍个GNSS点所测形变见图10。黑色箭头代表组合计算形变,棕色箭头代表GNSS形变,圆点颜色和形变场颜色表示垂直形变。

图  10  组合计算形变与GNSS形变

Figure  10.  Deformation of Combination and GNSS

下载: 全尺寸图片 幻灯片

由图10可知,在凯库拉地震震中东北方向,两个相距较近的地区检测到了非常复杂且巨大

的地表位移,它们中间有一小段未发生较大位移的区域,其中一个在东北地区,Kekerengu、Jordan Thrust及Papatea断层附近发生巨大的地表形变与十分明显的不连续性。在Kekerengu断层与Papatea断层周围有着最大超过9 m的水平形变,是本次地震发生最大地表形变的区域。Kekerengu断层滑动主要表现为右旋走滑兼少许逆冲,左侧有少量向上抬升。Jordan Thrust断层形变主要发生在右侧,表现为右旋走滑,右侧向上抬升的运动。这种异常的最大超过10 m的垂直方向的位移是由于Papatea右侧区域向西南运动造成的,Papatea断层表现出强烈的逆冲趋势以及少量的左旋走滑。在西南地区,产生的地表位移要小于东北地区,形变主要发生在Hope断层以南,Humps断层以及Hundalee断层区域,地表产生多个破裂,形变模式复杂,在Hope断层的南部以及Hundalee断层的北部有较大的形变分布。西南地区垂直位移与水平位移相比较大,表现为垂直位移,最大抬升超过5 m,整体显示出较强的逆冲运动并带有少量右旋走滑。从整体三维形变上看,凯库拉地震垂直形变表现为抬升,运动在东北地区以右旋走滑为主,在西南地区以逆冲为主。

4 结 语

本文利用哨兵2号光学数据和哨兵1号SAR数据,基于光学影像相关和像素偏移追踪技术,获取了2016年凯库拉Mw 7.8地震的同震形变场。根据获得的东西向、南北向、距离向和方位向形变场,将它们组成能够解算出三维形变的不同的组合,使用最小二乘方法计算东北部主要破裂区域三维形变。结果表明,OIC+POT_Des和OIC+POT_As获得的垂直形变与实际情况符合较差；对各类数据进行简单加权能够增强组合中的南北向形变精度及改善部分组合的东西向和垂直向形变；POT_As_Des所得的南北向形变效果和精度较差；OIC+POT_As_Des和OIC+POT_Range结合能够得到较好的三维形变。考虑数据约束性和结果的准确性,对于哨兵‍1‍号和哨兵2号数据,OIC+POT_As_Des的组合最适合获取凯库拉地震三维形变。使用OIC+POT_As_Des获取凯库拉地震整个区域同震三维形变,结果揭示了凯库拉地震巨大且复杂的破裂和地表形变,形变区域主要集中在东北地区的Kekerengu、Jordan Thrust和Papatea断层和西南地区的Hope、Humps和Hundalee断层附近,垂直形变主要表现为抬升,东北地区形变整体大于西南地区,该地震是以右旋走滑为主带有逆冲的地震。