利用升降轨Sentinel-1A影像提取2017年九寨沟Ms 7.0地震同震形变场

高二涛，兰艳萍，黄小梅，李 豪，雍 琦

(1.桂林理工大学测绘地理信息学院，桂林 541006；2.西南交通大学地球科学与环境工程学院，成都 611756；3.山西省工业设备安装集团有限公司，太原 030000)

2017年8月8日21点19分46秒，四川省阿坝州九寨沟县发生7.0级地震，震中位于距离九寨沟核心景区西侧约5 km处的比芒村(33°12′N，103°49′E)，震源深度20 km。距离震中200 km的范围以内，近5年期间发生过142次3级以上的地震，此次地震是震级最大、破坏最严重的一次。

随着InSAR(interferometric synthetic aperture radar)技术的深入研究，其在地震形变场提取、火山运动监测、城市地表沉降和滑坡监测等方面表现出了良好的应用前景[1]。洪顺英等人利用3种不同视线(line of sight,LOS)向的ENVISAT ASAR数据进行干涉处理，获取多视线向(multi-LOS)的同震形变场[2]。2018年，曹海坤等[3]利用InSAR观测技术融合GPS(global positioning system)数据，通过添加一个跟点位有关的系统误差函数，对系统误差进行修正，得到更为准确的三维形变场。在地震形变场观测中，D-InSAR(differential synthetic aperture radar interferometric)技术具有不受时相和气象条件限制、可以全天时全天候的探测地球表面细微的变化[4]、覆盖范围广、观测精度高等独特优势。

由于地震发生在九寨沟景区内，植被、森林密布给GPS、水准测量等常规方法带来了较大的困难，InSAR技术无需到实地测量，并且可以获取大范围监测信息，具有较大的优势。基于欧空局Sentienl-1A卫星获取的九寨沟震前和震后的升降轨雷达影像，结合精度较高的AW3D30(ALOS World 3D-30 m)数字高程模型，利用D-InSAR技术分别对升降影像进行数据处理，为了获取九寨沟地震的精确地表同震形变场，进行了升降轨结果的交叉验证分析。结合MATLAB和Surfer绘图软件实现对形变场的精度分析，最后导入Google Earth验证形变场位置的准确性。

1 差分干涉测量基本原理

差分干涉测量基本原理是基于地震前后SAR影像获得地表的形变信息。即利用两幅跨越形变周期的SAR图像对，共轭相乘获取干涉相位信息[5]。通过卫星观测几何参数和轨道信息去除参考椭球相位，使用外部DEM(digital elevation model)去除地形相位，扣除大气和噪声相位，最后得到地表形变相位信息。

合成孔径雷达差分干涉测量工作模式通常有两种，即通过两个传感器同时对同一区域进行观测的单轨模式和单一传感器重复对同一区域进行观测的重复轨道模式[6](图1)。

图1 差分测量原理

S1为卫星第一次成像的位置，S2为卫星第二次成像的位置，R1和R2分别对应两次成像到地面目标点的视线距离，令R2=R1+Δr。B是卫星两次成像的空间基线，α是基线与水平面的夹角，S1位置时卫星的侧视角为θ，S1与参考平面的高度为H，h为目标点P点与参考平面之间的高度。干涉相位可以表示为

Φ=Φflat+Φtopo+Φdef+Φatm+Φnoise

(1)

式(1)中：Φflat为平地相位；Φtopo为地形相位；Φdef为要获取的地表形变相位；Φatm、Φnoise分别为大气和噪声相位。在去除大气和噪声相位后，式(1)可表达为

(2)

(3)

式中：-Δr表示卫星两次对同一目标成像的路程差；λ表示雷达波长。其中，Φflat可以用传感器本身的头文件等相关参数计算后消除，Φtopo可以用外部DEM模拟出地形相位，并从中扣除。这样，地表形变相位可以通过式(3)获得。

2 九寨沟地质背景构造

根据九寨沟地震震中位置(图2)，结合已有的地质资料可知，此次地震位于巴颜喀拉块体东缘岷江断裂、塔藏断裂和虎牙断裂周围，是巴颜喀拉块体发生激烈碰撞的结果。东昆仑断裂向东部延伸形成塔藏断裂，在巴颜喀拉体块东北方向上，大体上呈现出西北走向的趋势，它涵盖了东北村、塔藏、九寨沟、沙尕里等广泛地区[7]。对塔藏断裂西段的研究，主要以左旋剪切走滑为主兼挤压活动为主，滑动速率为2.43～2.89 mm/a[8-9]。

图2 地震地质背景

塔藏块体东北部在青藏高原的东北方向上，是柴达木地块、扬子地块、巴颜喀拉地块的交汇地区，是我国南北地震带地震中段。此处地形极其复杂，板块之间的运动剧烈，发生过许多强震。岷江断裂北起九寨沟县，南至松潘县，断层走向主要以左旋为主，倾向为西北方向，倾角在60°～70°之间，显示出该断裂强烈活动的现象。

虎牙断裂是这次地震的发震断裂，呈现出由西北向南北方向的延伸，为上冲兼左旋走滑方式[10]。断裂的走向方向表现出上部和下部不同的现象，在上部，发生从北西方向向南北方向过渡，向东倾斜80°左右；在下部由南北方向向南东方向过渡，倾向向西，倾角在30°～70°范围变化。在这样的情况和地势下，使得巴颜喀拉块体向东部运动激烈，导致九寨沟地震的发生。

3 九寨沟地震形变场提取

3.1 实验数据

利用Sentinel-1A卫星获取LOS方向形变场。其中使用的4幅影像参数如表1所示。

表1 Sentinel-1A升降轨影像参数

DEM数据使用的是AW3D30。它是日本宇航局对地观测中心于2016年5月发布的全球数字高程模型，它是由ALOS(advanced land observation satellite)卫星上搭载的全色立体遥感测绘仪PRISM(panchromatic remote-sensing instrument for stereo mapping)生成，覆盖了北纬82°到南纬82°之间的全球陆地面积，空间分辨率为30 m，高程精度为5 m。AW3D30是当今用户获取的覆盖范围最广，精度最高的开源数字高程模型之一[11]。通过MATLAB软件对DEM数据进行拼接、裁剪得到研究区DEM，运用其去地形相位处理。

3.2 数据处理过程

主要使用ENVI软件中的SARscape模块，结合AW3D30 DEM数据分别对Sentinel-1A升降轨影像对进行差分干涉处理。技术路线如图3所示。

图3 D-InSAR技术路线

3.2.1 影像配准

首先找出两景影像的同名点，计算它们的坐标函数关系，按照函数关系把从影像经过简单的平移，使两幅影像有统一的分辨单元。配准的步骤分为粗配准和精配准[6]，得到配准信息如表2所示。

表2 主从影像配准信息

从表2可以看出，图像配准精度较好，距离向、方位向、信噪比和相关系数都比较稳定。

3.2.2 生成干涉图

对配准后的主影像和从影像进行共轭乘法运算来生成干涉图，设置距离向和方位向的多视比例为19∶5。评估所生成的干涉图的质量，通常用的方法是使用相干系数对干涉图进行评价[12]。相关系数的范围为[0,1]，一般相干系数大于0.3，说明干涉图质量较好。

对于平地相位使用传感器本身的头文件携带的轨道参数计算并减去；利用AW3D30 DEM模拟地形相位信息，并从干涉结果中去除。经过以上两步处理后，利用MATLAB对干涉图进行裁剪，去平、去地形相位后的干涉图如图4所示。

图4 扣除地形相位后的干涉图

从图4可以看出地震的干涉图条纹不是很清晰，因为存在噪声、近场失相关等的影响，需对其进行滤波处理，如图5所示为自适应滤波并编码的结果，通过滤波处理之后，干涉条纹信息明显清晰了很多。

图5 D-InSAR差分干涉图滤波

3.2.3 相位解缠

相位解缠的本质是恢复整周相位周期的过程。九寨沟景区，植被密度较高，使用最小费用流(minimum cost flow, MCF)的方法进行解缠可以提高解缠效率。这种方法对像元进行了充分的考虑，对于一些相干性系数小的像元，进行了掩膜处理，提高了解缠的精度。如图6所示，利用MATLAB裁剪出的相位解缠结果。

图6 DInSAR相位解缠

3.2.4 轨道精炼和重去平

轨道精炼是通过摄取控制点来校正细小的轨道误差，从而提高精度。在解缠和干涉图平坦的地方选点，震中区域不可以选点，由于震中失相关，会给重去平造成很大的误差。

3.2.5 相位转形变和地理编码

把轨道精炼后的结果转换成形变量，它的本质是把雷达的形变相位转换成一个形变量(以米为单位)，即乘以λ/(4π)以获得LOS方向形变场，并与外部DEM结合通过地理编码投影到WGS84(World Geodetic System 1984)坐标系。

利用MATLAB软件裁剪形变场，使用freadbk命令读取结果并裁剪出所需要的形变场，并且通过格式转换，导入Surfer软件，如图7所示，实验获取的升降轨形变场。

图7 D-InSAR获取的地表形变场

为了验证形变场位置的准确性，在Surfer软件中，把形变场转化为Google Earth能识别的KML(keyhole markup language)文件，将实验得到的升降轨形变场叠加显示在Google Earth地图上。从图8中可以看出形变出现的区域与实际地震发生的区域吻合度较高。

4 形变结果及分析

基于Sentinel-1A影像，利用InSAR技术，提取九寨沟地震升降轨形变场，对形变场的范围、时空分布以及形变大小进行分析，为地震断层反演打下良好基础。

首先，从图7可以看出，形变场的分布不对称，左边(西边)的形变程度较严重，右边(东边)形变相对较弱。其主要原因是该地发震断层处于塔藏断层、虎牙断层和岷江断层这三个断层交汇的地震带上，断层结构相对复杂所致。断层活动表现出不确定特性，结合历史资料，该地震带上发生过很多次大地震，显示出地震群体的特征。断层活动的变化发生的方式不同，变形场的分布特征也遵循不同的模式。2003年巴姆地震是一次右旋走滑形式，它在LOS形变上呈现出正四面体分布[13]；2008年汶川地震也是一次右旋走滑形式，它在LOS形变分布不明显[6]。九寨沟地形复杂，在提取地震形变场时，需要考虑断裂山脉的滑坡影响[14]。九寨沟地震形变场与巴姆地震形变场有类似之处，都处于四象限分布，不同之处在于九寨沟地震形变场是不对称分布，而巴姆地震形变场是对称分布。

其次，从图5可以看出，震中区域的相干性较低。其主要原因有两个，一是因为九寨沟景区植被森林密集和地形复杂所致。植被森林密集在雷达反射过程中，由于绿色植被中水分含量较多，吸收率变大，因此反射率就大大下降了。雷达波接收不到反射光信号，在干涉处理过程中就出现了干涉失相关的现象。此外，震中造成的破坏力度强，导致干涉失相关。西北和东南方向的干涉条纹显示出不同形状。在西北方向干涉条纹表现出彩虹渐变色变化趋势；而在东南方向上变化较弱。从另一个层面上说明地形、地物复杂，变形对该地的破坏程度不同。

最后，从整体上看，它的总体形状是一个“核仁”状，长度大约有24 km，宽度大约有29 km，并且呈现出西南向东北走向。从MATLAB中精确的读取到九寨沟地震LOS向形变的相关参数，形变场的面积大约为696 km2，上升的最大形变值达到22.4 cm，约有8个干涉条纹；下降极值约为11.2 cm，约有4个干涉条纹。说明该地区变形较大，地震破坏程度较大。

5 结论

利用升降轨Sentinel-1A数据，基于D-InSAR技术进行SAR数据处理，通过图像配准等一系列数据处理，结合已有的资料和谷歌地球等数据，利用MATLAB、Surfer软件实现对形变场的裁剪和显示，得到了2017年九寨沟7.0级地震精确同震形变场。最后本文对形变场的位置精度、形变量级大小等进行了系统的分析，为下一步地震断层位置反演打好了良好的前期基础。

致谢:感谢欧空太空局在科学数据中心网站公布的Sentinel-1A卫星SAR影像；感谢日本宇宙航空航天局地球观测研究中心提供的AW3D30 DEM数据。