融合DInSAR和MAI技术反演九寨沟地震三维形变场

宫熙雯,张正加,陈启浩,刘修国

(中国地质大学(武汉)，湖北 武汉 430070)

2017年8月8日21:19:46，四川省北部阿坝州九寨沟县发生7.0级地震(33.2°N，103.82°E)，震源深度20 km，截至2017年8月18日09:00共记录到地震总数为5 020个(余震总数5 019个)，其中最大余震为：2017年8月9日10:17四川阿坝州九寨沟县4.8级。此次地震，中国地震局徐锡伟等专家判断这次地震与巴颜喀拉块体南东向推挤有关系，发震构造推测为塔藏断裂南侧分支和虎牙断裂北段[1]。

一些学者对该地震进行了相关研究，郑绪君等采用迭代反褶积和叠加方法反演得到破裂模型，用于确定地震的破裂过程[2]，汪建军等采用库仑应力模型研究地震激发的库仑应力变化[3]，文献[4]采用多步定位法对九寨沟地震Ms7.0的主震及其余震进行重新定位，提出震源机制解决方案，并通过CAP法(Cut and Paste Method)反演了九寨沟地震的质心深度。除了对以上地震参数的反演，获取地震同震形变场对揭示地表形变、理解和分析发震构造具有重要意义。由于差分雷达干涉测量(DInSAR)可以获取全天时全天候、大范围、连续空间覆盖、高精度的地震区域同震形变图[5]，在地震监测方面应用广泛，对于本次九寨沟Ms7.0地震，赵强等利用Sentinel-1A和RadarSat-2升轨数据结合两轨法DInSAR反演了地震的同震形变场[6]；熊威等使用Sentinel-1A数据进行同震形变场反演，并确定形变区的相干系数[7]；申文豪等采用DInSAR方法利用升、降轨Sentinel-1数据反演了同震形变场[8]。Zhao基于InSAR和GPS测量，反演地震同震滑动分布和倾角，并获得地震视线向的同震形变场[9]。但由于DInSAR方法的设计本质是一种侧视雷达测量技术，得到的并不是地表真实形变，而是地表真实形变矢量通过一系列坐标系转换和几何转换后在视线向的投影[10]。因此，同震形变的反演仅仅依靠DInSAR技术获取视线向形变是不够的，需要增加观测角度以建立地震三维同震形变场。

1999年，Michael 等人提出了一种基于方位向像素偏移(Azimuth Pixel Offset，AZO)的技术[11]，2002年，Fialko等结合DInSAR和像元偏移量技术(Offset Tracking)成功获取Hector Mine地震的三维位移场[12]。2004年，Wright 等提出通过不同入射角影像进行干涉的多视角InSAR技术并反演Nenana Mountain地震的三维形变场[13]，王家庆等在此基础上结合先验条件的最小二乘迭代逼近法解算改后地震的三维同震形变场[14]，但由于多视角InSAR对数据入射角差异要求较高，获取满足条件的数据源对绝大多数地区较难实现。2006年，Bechor等提出了多孔径干涉技术(Multiple Aperture InSAR，MAI)以获取方位向形变，较AZO技术获取方位向位移矢量精度更高[15]。2014年，Wang等使用该方法并结合DInSAR方法，使用L波段PALSAR数据和C波段ASAR数据反演了玉树地震三维形变场，实验结果表明，使用DInSAR和MAI方法可有效提高震源参数的反演精度[16]。

文中在九寨沟Ms7.0地震同震形变场反演过程中，融合升降轨Sentinel-1A数据，使用DInSAR和MAI组合方法获取雷达视线向和方位向的形变，构建地震三维同震形变场，获取本次地震在竖直方向、南北方向、东西方向的地表位移三分量，较单一雷达视线向更充分地解释了地表真实形变情况。

1 实验区与数据源

1.1 实验区概况

在四川省的北部，阿坝州九寨沟县于2017年8月8日发生了一次Mw7.0级地震，震源中心为(33.20°N，103.82°E)。地震发生在青藏高原与四川盆地两大地貌单元过渡的深切割高山峡谷地带，地处岷山山脉北段的复背斜上，地势西北高，东南低，平均海拔在4 000 m以上，受灾面积约为840 km2(35 km×24 km)。研究区内地形复杂，有多条断裂带纵横交错，震源中心位于塔藏断裂南侧，虎牙断裂延长线北部，与岷江断裂相交[1,17]，如图1所示。

图1 研究区概况

1.2 数据选取

哨兵1号卫星是欧洲航天局哥白尼计划(GMES)中的地球观测卫星，为双星系统，搭载C波段合成孔径雷达，可全天时、全天候不受天气条件影响地提供连续影像。本研究使用数据为Sentinel-1A卫星干涉宽幅模式的斜距单视复数产品(IW SLC)，空间分辨率为5 m×20 m，极化方式为VV。选取了震前震后的升降轨SAR数据共4景用于本次地震同震形变反演，影像基本信息如表1所示。本研究所采用的参考DEM为覆盖震区的SRTM-3 version4 DEM数据。本次实验选取的升降轨震前震后两个影像对的时间基线相近均为12 d，空间基线也相对较小，分别为36.7 m和65.4 m，尽可能减少了时空失相关对最终结果造成的误差。

表1 Sentinel-1A影像参数信息

2 研究方法

文中使用DInSAR和MAI结合的方法获取地表位移分量。由于MAI所获取方位向位移的精度在方位向形变获取中精度较高，且和DInSAR方法在一个数量级[18]，因此可结合上述数据处理得到的形变图数据，对九寨沟地区的地表形变三分量进行解算，方法流程图如图2所示。

图2 方法流程图

2.1 DInSAR

DInSAR技术是通过对同一区域不同时相两景影像产生的干涉图差分处理，借助外部DEM消除地形影响而获得形变信息一种技术手段[19]。DInSAR获取视线向形变的主要步骤为：选取主、辅影像生成干涉图，并借助外部DEM进行自动配准并进行去平地效应；采用Goldstein方法进行前置滤波以提高干涉条纹的清晰度；采用Delaunay MCF方法进行解缠以消除低相干像元的影响，只对高相干像元进行解缠，最小化相位突变的影响，减少“孤岛”情况的存在；对于山区可能存在的大气延迟，主要表现为与地形相关的相位斜坡，用线性模型进行校正[20]，对于轨道误差引起的相位斜坡，用双二次模型进行校正[21]，后经轨道精炼、重去平及地理编码等得到视线向形变。

2.2 MAI

MAI技术是基于分光束InSAR进行数据处理，生成前向和后向干涉图，根据两幅干涉图上的相位差获取轨道方位向的位移分量的一种干涉测量方法[22]。MAI技术获取方位向形变的主要步骤为：选取主、辅影像生成干涉图，主、辅影像分别通过带通滤波，得到各自的前视和后视的单视复数影像，主辅影像的前视和后视数据对应进行干涉处理，分别得到前视干涉图和后视干涉图，并将二者进行干涉处理得到MAI干涉图，通过Goldstein滤波、区域增长法解缠，对干涉图中存在的轨道误差和大气延迟进行校正后，经轨道精炼、重去平及地理编码后得到方位向形变。

2.3 三维形变场

文中使用的升轨数据的轨道运行方向为东南—西北，降轨数据的轨道运行方向为东北—西南，由于Sentinel-1卫星为右视雷达，升降轨数据与真实地表形变三维分量的几何关系分别如图3、图4所示。

图3 升轨与地表三维向几何关系

图4 降轨与地表三维向几何关系

图中dU，dN，dE分别为竖直方向、南北方向、东西方向上的形变，规定向上、正北、正东为正；dlos，daz分别为视线向、方位向上的形变，规定沿卫星轨道方位向为正，统一正向可得到地震三维形变的观测方程为：

(1)

式中:θ为雷达波入射角，α为轨道方位角，下标A，D分别表示升轨和降轨。观测矢量D:

D=(dlos,A,dlos,D,daz,A,daz,D)T.

(2)

得到误差方程：

V=CX-D.

(3)

其中，向量V为每次观测的误差，即：

V=(Vlos,A,Vlos,D,Vaz,A,Vaz,D)T.

(4)

由加权最小二乘原理[23]可计算三维位移场：

X=(ATPA)-1ATPD.

(5)

式中:C为观测方程的系数矩阵，P为观测权阵，由于DInSAR测量误差和MAI的方位向形变误差都在cm级，近似认为4次观测的精度是相等的，即权阵P为单位矩阵。

3 结果分析

3.1 InSAR处理结果

融合升降轨SAR数据使用DInSAR和MAI技术反演的地震升降轨视线向、方位向同震形变如图5所示，同震位移剖面如图6所示。

综合分析图5、图6可知，本次九寨沟Ms7.0地震，形变主要在西北区域和东南区域，存在非对称性分布特征。升降轨数据的形变特征显示出不一致性：升轨数据形变反演结果显示，震中西北侧区域发生了沿视线向最大值为0.19 m、沿方位向最大值为0.25 m的沉降形变，震中东南侧区域发生了沿视线向最大值为0.10 m、沿方位向最大值为0.11 m的隆升形变；而在降轨数据形变反演结果中，震中西北侧区域发生了沿视线向最大值为0.14 m、沿方位向最大值为0.14 m的隆升形变，震中东南侧区域发生了沿视线向最大值为0.10 m、沿方位向最大值为0.14 m的沉降形变。这种不一致性主要与升降轨不同的成像几何特征以及地表形变方向有关，而且升降轨呈现出相反的形变态势，表明九寨沟Ms7.0地震的地表形变以水平形变为主，与地震学结论保持一致[24]。

由方位向升降轨同震形变图可知，九寨沟Ms7.0地震不仅引发了地表竖直方向上的位移，同时在地表平面上也存在位移且形变量和视线向的形变量为同一个数量级，在以水平形变为主的九寨沟地震中，构建三维同震形变场对于其发震机理及构造分析有着重要的意义。

图5 九寨沟Ms7.0地震同震形变场

图6 九寨沟Ms7.0地震同震形变剖面

3.2 三维形变场

将升轨雷达波入射角θA、降轨雷达波入射角θD、升轨轨道方位角αA、降轨轨道方位角αD带入观测方程的系数矩阵C：

(6)

进而求得三维形变场X的系数矩阵Z：

(7)

从而得到地表的三维同震形变分量，即：

dU= 0.645 8IlosA+0.646 3IlosD+

0.0934IAzA-0.093 5IAzD.

(8)

dN=-0.000 1IlosA+0.000 1IlosD+

0.512 9IAzA-0.512 9IAzD.

(9)

dE=-0.716 9IlosA+ 0.716 5IlosD-

0.258 6IAzA-0.258 4IAzD.

(10)

式中：IlosA，IlosD，IAzA，IAzD分别表示升轨视线向形变结果、降轨视线向形变结果、升轨方位向形变结果和降轨方位向形变结果。根据本次地震的影响范围和升降轨数据的重叠情况对形变结果进行了裁剪，选取103.49°E—104.11°E、32.93°N—33.55°N的震中区域，并借助Matlab编程计算求解九寨沟地震三维同震形变场，得到地表位移三分量，地震三维同震形变图如图7所示。实验结果表明，本次地震的断裂带为西北—东南(NW-SE)走向，位置大致与地震实际破裂带相符[25]。在竖直方向的同震形变量较小且震源中心附近上下形变贡献较为平均，存在震源中心东南侧最大向上0.09 m、西北侧最大向下0.11 m的形变；在南北方向存在较为明显的位移分界，存在震源中心东南侧最大向北0.09 m、西北侧最大向南0.18 m的形变；在东西方向的同震形变量较大，且存在明显的位移分界，存在震源中心西北侧最大向东0.25 m、东南侧最大向西0.12 m的形变。综合地表位移三分量可以判断，地表在震源中心的西北侧表现为沉降运动，在震源中心的东南侧表现为隆升运动，但地表运动纵向起伏较小，地表平面地表位移运动较为明显，与李渝生等在震后实地考察测量结果相符[17]。

图7 九寨沟Ms7.0地震三维同震形变场

九寨沟地震的震源中心被塔藏断裂、岷江断裂和虎牙断裂的延长线3条较大的断裂带包围，图7中黑色虚线位置为本次地震中地表出露的位置[25]，即断层所在地表部分的破裂区域，该断层两侧在图中的位移在3个方向均较为明显。地震断裂面西北侧(下降盘)错动位移较大，主要向东南偏东方向滑动；东南侧(上升盘)较西北侧错动位移小且3个方向位移的贡献较为平均，大致向正西北方向滑动，总体表现为震中东南方向、西北方向向内聚敛，震中西南方向、东北方向向外拉张，与中国地震局地壳应力研究所申文豪等的强震记录解算同震形变结果[8]及中国地震局第二监测中心刘芳等得到的九寨沟Ms7.0地震的重力场结果保持一致[26]。综合断裂带邻近区域的同震形变场及断层位移形式，可以认为引发本次地震的断层性质为左旋走滑断层，与单新建[25]、李勇[27]、任俊杰[28]等对本次地震的研究结果吻合。

3.3 精度分析

由于震中位于山区地带，以地面观测站数据为基础的常规观测技术在本次地震中难以实施，震中附近观测站点分布稀疏，文中研究区内只有GNSS的J416、JB33流动观测站和GPS的51JZZ、51JZB测站，点位数目较少[8,29]。很多学者在九寨沟地震研究过程中选择了InSAR方法获取视线向的形变数据[8,24-25,30-31]，其反演结果如表2所示。由于选用卫星、运行方向、时空基线、数据处理方式不同等原因，形变结果存在一定的差异，整理了地震视线向升降轨结果得到：升轨数据得到的隆升形变区间为0.07～0.11 m、沉降形变区间为 0.18～0.22 m；降轨数据得到的隆升形变区间为0.10～0.16 m、沉降形变区间为0.05～0.10 m，本次实验结果均在区间内，可以认为由此得到的三维形变结果是可信的。

表2 不同研究者对九寨沟地震视线向形变的反演结果

4 结 论

以四川省九寨沟7.0级地震为研究对象，融合升降轨Sentinel-1A数据，使用DInSAR和MAI结合的方法对该地震的雷达视线向和轨道运行方位向形变进行了提取，构建了九寨沟Ms7.0地震三维同震形变场，结果表明：九寨沟7.0级地震在震源中心的西北侧表现为沉降运动，在震源中心的东南侧表现为隆升运动，地表运动以水平形变为主，纵向起伏较小，其中，东西方向对本次地震的同震形变贡献最大，最大形变值为0.25 m。综合断裂带邻近区域的同震形变场及断层位移形式，可判断引发本次九寨沟Ms7.0地震的断层性质为左旋走滑断层。另外，得到的地表位移三分量，更加准确地反演了地震真实地表的同震形变情况，也可为地震发震机理和地震烈度判断、震后重建等后续工作提供三维数据支持，具有一定的实用价值。文中将观测精度近似认为相等并选用了单位矩阵进行权重评定，这种近似存在一定的误差，三维形变场的提取模型的优化有待进一步研究。