0 引言

北京时间2017年8月8日21时19分46秒，四川省阿坝州九寨沟县发生M_S7.0地震，震中位于33.20°N、103.82°E，震源深度20km，最大烈度达Ⅸ度。地震发生后，国内外多家机构发布了此次地震震源机制(表 1)，其中，中国地震台网中心(CENC)测定的面波震级为M_S7.0，断层向NE倾，倾角约55°；美国地质调查局(USGS)和全球矩心矩张量数据中心(GCMT)确定的矩震级为M_W6.5，断层向SW倾，倾角约80°，相对中国地震台网中心的结果较陡。不同机构的震源机制解均表明九寨沟地震是一次左旋走滑型地震事件。

此次地震发生在青藏高原中部巴颜喀拉地块东北角(邓起东等，2003)，向东与四川盆地接壤。巴颜喀拉地块是青藏高原近20年来构造活动最强烈的次级块体(单新建等，2017)。此次地震震中位于塔藏断裂、岷江断裂和虎牙断裂交汇处，这3条断裂均为全新世活动断裂，断层交错，地质构造复杂，地震活动频繁，是南北地震带北段高危区。

近30年来，InSAR技术因其高精度、高时空分辨率、全天时广域监测等特性，广泛应用于地震同震、震间和震后形变监测、城市地面沉降监测、地质灾害隐患早期识别等领域(王超等，2002；张景发等，2008；许强等，2019；王琪等，2020；葛大庆等，2019)。九寨沟M_S7.0地震发生后，许多学者利用InSAR技术提取了同震形变场。单新建等(2017)利用欧空局(ESA)升、降轨Sentinel-1A数据提取了九寨沟地震InSAR同震形变，并以观测结果为近场约束，运用SDM程序反演了破裂断层滑动分布，升、降轨最大雷达视线方向(LOS)形变量分别约22cm和14cm；姚鑫等(2017)利用两景升轨Sentinel数据获取了九寨沟地震的同震形变分布，极震区2000km²范围内LOS向形变在-13cm~28cm之间变化；季灵运等(2017)基于2个升轨和1个降轨Sentinel-1数据，利用InSAR技术获取了九寨沟地震同震形变场，结果显示地震造成的LOS向同震形变最大达20cm，同震形变场具有非对称展布特征。上述研究均基于InSAR雷达视线方向的一维形变监测结果展开，然而任何的形变过程均是三维的，提取地震同震三维形变对认识断层性质、研究震源机理具有更显著的意义。针对单一的InSAR技术只能获取LOS向一维形变，无法提取地表真实的三维形变信息问题，本文基于覆盖九寨沟地震极震区的Sentinel-1升、降轨数据和Radarsat-2升轨数据，利用3个不同入射角的卫星平台进行多视线联合观测，从而解算此次地震的同震三维形变场。

1 InSAR观测的雷达视线向形变与地表三维形变的关系

InSAR是一种侧视成像技术，其不足之处在于InSAR直接观测的形变信息并非地表真实形变，而是地表沿垂直向、SN向、EW向3个形变向量在LOS向投影的矢量叠加和(图 1)，这便是InSAR形变测量固有的视线向模糊问题(洪顺英，2010；高明亮等，2017)。雷达视线方向的形变分量只能定性地解释地表形变的大致状况，不能满足定量形变信息分析的需求。因此，从InSAR雷达视线向一维监测结果中进一步分离出地表形变三分量信息尤为重要。

建立局部UNE直角坐标系，构建地表三维形变场，U、N、E分别取垂直向上、北、东方向为正。InSAR雷达视线向结果可以分解为垂直向形变分量和地距向形变分量，如图 2所示，地距向形变分量又可进一步分解为正东方向形变分量和正北方向形变分量(α₁=α-π)。

InSAR解算的LOS向形变与U、N、E形变之间的转换关系为

$ d_{\mathrm{LOS}}=d_{\mathrm{U}} \cos \theta-d_{\mathrm{gr}} \sin \theta $

(1)

$ d_{\mathrm{gr}}=d_{\mathrm{N}} \sin \alpha_{1}-d_{\mathrm{E}} \cos \alpha_{1}=d_{\mathrm{E}} \cos \alpha-d_{\mathrm{N}} \sin \alpha $

(2)

因此，任意一地面点的InSAR结果与地表三维形变分量可表示为

$ d_{\mathrm{LOS}}=d_{\mathrm{U}} \cos \theta+\left[d_{\mathrm{N}} \sin \alpha-d_{\mathrm{E}} \cos \alpha\right] \sin \theta $

(3)

其中，d_LOS表示地表沿雷达视线方向的形变分量；d_U、d_N、d_E分别代表垂直向、南北向、东西向的形变分量；θ和α分别表示合成孔径雷达视线方向的入射角和卫星运行的轨道方位角；d_gr表示地距向形变分量。

2 地表三维形变解算方法 2.1 InSAR观测与方位向测量联合解算方法

求取地表三维形变分量的一种方法是利用至少2个轨道的雷达视线向和方位向的InSAR观测结果，根据最小二乘法对其进行优化，获取地表三维形变分量。其中雷达视线向的地表形变观测结果可由D-InSAR和时序InSAR(MT-InSAR)方法获得；方位向形变分量可通过偏移量跟踪技术(Offset-tracking)或者多孔径干涉测量(MAI)技术(季灵运等，2009；胡俊等，2010；杨红磊等，2014)获取。

偏移量跟踪技术和多孔径干涉测量技术在方位向的形变监测精度一般在分米级甚至是米级(王晓文等，2014)，大大低于InSAR技术的观测精度，因此一般适用于地震、矿区、火山喷发和冰川运动引起的大地表位移量监测。

2.2 InSAR与GPS数据联合解算方法

InSAR与GPS具有很强的互补性，GPS站网观测具有大区域、连续观测、离散分布、水平精度高等观测优势，联合GPS和InSAR可充分发挥GPS水平精度高以及InSAR视线向观测精度高且对垂直向敏感的优势，较好地将离散的GPS观测点和地表连续的InSAR结合(许才军等，2003；单新建等，2014；伍吉仓等，2020)。

InSAR与GPS联合观测最早由Gudmundsson等(2002)在研究Reykjanes半岛三维形变时提出，应用马尔可夫随机场和模拟退火算法实现了GPS和InSAR的有机融合。Samsonov等(2006)根据贝叶斯统计的随机场理论，对GPS和InSAR联合提取三维形变的算法进行了改进，成功应用于加州地区的三维形变研究。胡俊等(2010)利用GPS与D-InSAR的结果确定加权方案，不依靠先验知识即可提取地表三维形变场。陈威等(2018)联合GNSS连续站和升、降轨Sentinel-1数据，获取了九寨沟地震同震形变场，并基于均匀弹性半无限位错模型，反演了发震断层的滑动分布，计算了同震库伦应力变化。

在InSAR与GPS观测融合过程中，不同观测数据的权重配比对结果影响较大，权重配比问题是当前InSAR与GPS联合解算的热点问题。目前，国内GPS站网的空间密度在数千米至数百千米范围，一般不能满足星载雷达卫星高空间分辨率的监测需求，且GPS站点建设和运维管理费用昂贵，这在一定程度上降低了InSAR与GPS数据联合解算方法的普适性。

2.3 多平台InSAR数据联合观测方法

当前InSAR卫星正在朝多角度、大幅宽、高时空分辨率方向发展，不同源、多轨道的InSAR数据集使得相同覆盖范围内的多源InSAR立体观测成为可能(刘国祥等，2012)。如式(3)所示，假设针对某一地震事件有n个InSAR观测资料，可创建n个相互独立的观测方程(薛腾飞，2013；王家庆，2015；Chang et al，2017)，即

$ d_{\mathrm{LOS}}^{n}=d_{\mathrm{U}} \cos \theta_{n}+\left(d_{\mathrm{N}} \sin \alpha_{n}-d_{\mathrm{E}} \cos \alpha_{n}\right) \sin \theta_{n} \quad n=1, 2, 3, \cdots, n $

(4)

其中，d_LOSⁿ表示第n次观测地表沿LOS方向的形变分量；θ_n和α_n分别表示第n次观测的雷达信号入射角和卫星运行的轨道方位角。

观测方程共有3个未知数，当n=3时，表达式可以改为

$ \left[\begin{array}{c} d_{\mathrm{L} . \mathrm{S}}^{1} \\ d_{\mathrm{LOS}}^{2} \\ d_{\mathrm{LOS}}^{3} \end{array}\right]=\left[\begin{array}{ccc} \cos \theta_{1} & \sin \theta_{1} \sin \alpha_{1} & -\sin \theta_{1} \cos \alpha_{1} \\ \cos \theta_{2} & \sin \theta_{2} \sin \alpha_{2} & -\sin \theta_{2} \cos \alpha_{2} \\ \cos \theta_{3} & \sin \theta_{3} \sin \alpha_{3} & -\sin \theta_{3} \cos \alpha_{3} \end{array}\right]\left[\begin{array}{c} d_{\mathrm{U}} \\ d_{\mathrm{N}} \\ d_{\mathrm{E}} \end{array}\right] $

(5)

令矩阵$ \boldsymbol{A}=\left[\begin{array}{ccc}\cos \theta_{1} & \sin \theta_{1} \sin \alpha_{1} & -\sin \theta_{1} \cos \alpha_{1} \\ \cos \theta_{2} & \sin \theta_{2} \sin \alpha_{2} & -\sin \theta_{2} \cos \alpha_{2} \\ \cos \theta_{3} & \sin \theta_{3} \sin \alpha_{3} & -\sin \theta_{3} \cos \alpha_{3}\end{array}\right]$，其逆矩阵为$\boldsymbol{A}^{-1}, \boldsymbol{D}=\left[\begin{array}{lll}d_{\mathrm{LOS}}^{1} & d_{\mathrm{LOS}}^{2} & d_{\mathrm{LOS}}^{3}\end{array}\right]^{T}$，$\boldsymbol{X}=\left[\begin{array}{lll}d_{\mathrm{U}} & d_{\mathrm{N}} & d_{\mathrm{E}}\end{array}\right]^{T}$。通过矩阵运算可求解出X

$ \boldsymbol{X}=\boldsymbol{A}^{-1} \boldsymbol{D} $

(6)

当有3个以上观测方程，考虑InSAR干涉测量观测误差E，误差方程可表示为

$ \mathit{\boldsymbol{E}}{{ = }}\mathit{\boldsymbol{A X}}{{ - }}\mathit{\boldsymbol{D}} $

(7)

根据最小二乘原理，未知向量X的最优无偏估计为

$ X=\left(\boldsymbol{A}^{T} P \boldsymbol{A}\right)^{-1} \boldsymbol{A}^{T} P \boldsymbol{D} $

(8)

其中，A为系数矩阵，由InSAR卫星观测几何决定；P为观测值权重；D为雷达视线向观测矩阵，E为对应的观测误差矩阵。

九寨沟M_S7.0地震发生后，我们收集了覆盖九寨沟地震极震区的升、降轨Sentinel-1数据和Radarsat-2升轨数据(图 3)，理论上可实现九寨沟地震的三维立体观测。鉴于目前尚未有相关学者利用3个轨道的InSAR数据对九寨沟地震进行同震三维形变提取研究，故本文采用3种观测模式的卫星平台进行InSAR多视线观测，联合解算九寨沟地震同震三维形变场。

3 InSAR同震形变场提取 3.1 InSAR影像

为提取九寨沟M_S7.0地震的同震雷达视线向形变结果，选取欧空局的Sentinel-1数据和加拿大MDA公司的Radarsat-2(RS-2)InSAR数据(表 2)。其中，Sentinel-1数据为C波段、图像分辨率5m×20m、VV极化的SLC数据，幅宽为250km；Radarsat-2数据为超宽精细模式、分辨率5m×5m、HH极化的SLC数据，幅宽为125km。

3.2 差分干涉处理

采用GAMMA软件进行数据处理，基于二轨法差分干涉生成差分干涉图。对于Sentinel-1数据，利用欧空局提供的精轨数据去除轨道相位，利用30m分辨率的SRTM DEM数据去除地形相位，为使Sentinel干涉像对保持较高的空间分辨率，并同时达到去除斑点噪声的目的，选取8 ︰ 2的距离向、方位向视数比，利用最小费用流法(Eineder et al，1998；祝杰等，2020)进行相位解缠；对于Radarsat-2数据，由于两期数据成像时间间隔较长(72天)，降低了SAR干涉对的相干性，故采用较大的20 ︰ 10距离向、方位向视数比进行多视处理，采用加权功率谱法对干涉图进行3次滤波，滤波窗口依次设置为128、64、32，以提升干涉图的相干性，相位解缠同样采用最小费用流法。基于表 2中3个干涉组合各自的数字高程模型建立大气延迟模型，从原始差分干涉相位中去除大气延迟相位，获取去除各种相位误差的清晰干涉条纹。最后，将解缠后的干涉相位转换为地表沿雷达视线方向形变。

3.3 同震形变场结果分析

利用Sentinel-1升、降轨和Radarsat-2升轨数据提取的九寨沟M_S7.0地震雷达视线向同震形变结果，如图 4所示。由图可见，3个形变场的空间分布大体一致，同震形变主要发生在塔藏断裂、虎牙断裂、雪山梁子断裂、岷江断裂的交汇处，总体呈“果仁状”(姚鑫等，2017；Liu et al，2019)，同震形变场长轴方向为NW向，形变区范围约55km×45km。

图 4中的形变正值表示靠近卫星飞行方向，即LOS向缩短，形变负值表示远离卫星飞行方向，即LOS向拉伸。Sentinel-1升轨、Radarsat-2升轨和Sentinel-1降轨LOS向缩短的最大形变量分别为8.5cm、8.6cm、12.9cm；LOS向拉伸的最大形变量分别为19.5cm、14.7cm、10.1cm。升轨和降轨数据提取的形变结果运动方向相反、量级相当，表明断层的运动性质具有典型的走滑变形特征。

3个InSAR雷达视线向形变监测结果稍有差异，主要原因是雷达信号入射角、卫星轨道方位角差异导致地表真实形变在雷达视线向具有不同的敏感度，由此证明单一的InSAR监测不能提取地表真实的形变信息。本文在空间后方交会思想(刘国祥等，2012)的基础上，联合雷达视线向同震形变结果(图 4)，提取九寨沟M_S7.0地震同震三维形变场。

4 多平台三维形变场提取与结果分析 4.1 三维形变场提取与结果分析

使用Sentinel-1升、降轨和Radarsat-2升轨三种观测模式的数据提取了九寨沟地震同震LOS向形变，3种观测数据构成的空间交会关系如图 5所示。Sentinel-1和Radarsat-2升轨数据的雷达信号入射角分别为43.86°和34.99°，Sentinel-1降轨数据的雷达信号入射角为39.25°，较大的入射角差异为联合解算空间三维形变场提供了有益的观测几何条件(刘国祥等，2012；杨梦诗等，2017)。

根据式(6)，联合3个轨道InSAR雷达视线方向观测结果，再结合卫星飞行的轨道参数，可精确计算出地表垂直向、SN向和EW向形变。图 4给出了九寨沟地震Sentinel-1升、降轨和Radarsat-2升轨的LOS向形变，由表 2给出的3个干涉对卫星飞行入射角和轨道方位角的准确值，可计算系数矩阵$\boldsymbol{A}=\left[\begin{array}{ccc}0.7210 & -0.1545 & -0.6755 \\ 0.7744 & -0.1408 & 0.6168 \\ 0.8193 & -0.1109 & -0.5626\end{array}\right]$，其逆矩阵$\boldsymbol{A}^{-1}=\left[\begin{array}{ccc}-2.5097 & 0.2039 & 3.2367 \\ -15.9974 & -2.5112 & 16.4532 \\ -0.5014 & 0.7919 & -0.3072\end{array}\right]$，因此得到九寨沟M_S7.0地震同震三维形变结果，即

$ \left[ {\begin{array}{*{20}{l}} {{d_{\rm{V}}}}\\ {{d_{\rm{N}}}}\\ {{d_{\rm{E}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} { - 2.5097}&{0.2039}&{3.2367}\\ { - 15.9974}&{ - 2.5112}&{16.4532}\\ { - 0.5014}&{0.7919}&{ - 0.3072} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {d_{{\rm{LOS}}}^1}\\ {d_{{\rm{LOS}}}^2}\\ {d_{{\rm{LOS}}}^3} \end{array}} \right] $

(9)

最终得到的形变结果如图 6所示，对比图 6中垂直向、SN向和EW向地表形变结果发现，SN向形变量的极值约为垂直向和EW向的6倍。对于相同的雷达视线向形变监测误差，三维形变解算中SN向形变解算的精度最差，垂直向和EW向的解算精度相当(洪顺英，2010)，其主要原因是InSAR技术对地表SN向形变测量敏感度最差，导致SN向形变求解模型中系数值最大，放大了测量误差。

此次九寨沟地震造成的极震区地表位移在垂直向上方向达23.4cm，垂直向下方向达17.6cm；北向位移达141.8cm，南向位移达100.2cm；东向位移达22.0cm，西向位移达10.7cm，三维形变场的解算结果与前人研究结果(张庆云，2019；王燕燕等，2019)基本一致。地表形变以水平向形变为主，与震源机制解显示的左旋走滑型地震的运动特征相吻合。

从三维形变的空间分布特征看，垂直向、SN向形变的空间分布相对较杂乱。由EW向形变结果可见，树正断裂两端地块呈水平相对运动(易桂喜等，2017)，断层破裂主要发生在树正断裂上盘，且沿断层面向东形变剧烈，与前人给出的断层滑动分布反演结果(单新建等，2017；季灵运等，2017)一致。

4.2 解算过程精度评价

对比中国大陆构造环境监测网络和北斗地基增强系统GNSS连续站分布(王阅兵等，2018；陈威等，2018)，在3个InSAR影像重叠区域未找到上述2个GNSS观测网络的连续站点，因此无法以GNSS站点的形变量作为基准进行精度验证。本文将解算后的地表三维形变投影至卫星观测的雷达视线方向，并将投影后的雷达视线向形变量与InSAR技术直接提取的雷达视线向形变量对比(洪顺英，2010)，以此评价本文解算过程的精度。

图 6给出的九寨沟地震三维形变结果是利用3个雷达视线向观测值，运用代数运算直接解算出来的，各视线向的形变量相互独立，且均有厘米至毫米级的监测精度，相对于同震米级的大尺度形变而言，InSAR干涉测量的误差很小。通过比较InSAR直接提取的3个LOS向形变结果与三维形变场在Sentinel-1升、降轨和Radarsat-2观测下雷达视线向投影值(图 7)，发现3种观测的残差均小于0.01cm，即解算过程引起的误差小于0.01cm，证明通过3个InSAR平台联合观测方法直接解算九寨沟地震地表三维形变的过程具有较高的精度。

5 结论

利用覆盖九寨沟M_S7.0地震的Sentinel-1升、降轨和Radarsat-2升轨3个独立观测的InSAR数据，分别提取了九寨沟地震地表沿雷达视线方向形变。在空间立体交会思想下，运用多平台联合观测方法解算了九寨沟地震沿地表真实的垂直向、SN向和EW向形变，得到初步结论如下：

(1) 当前InSAR技术受限于其固有的视线向模糊，多平台InSAR数据联合方法在不借助其他观测手段的前提下，直接利用多个独立的InSAR观测即可解算出地表真实的三维形变信息，体现出该方法的技术优势。

(2) 3个轨道InSAR测量均清晰监测到九寨沟M_S7.0地震同震雷达视线向形变，主要发生在塔藏断裂、虎牙断裂、雪山梁子断裂、岷江断裂交汇处，呈“果仁状”，长轴方向为NW向，形变区范围约55km×45km。LOS向缩短最大形变量为12.9cm(降轨)，LOS向拉伸最大形变量为19.5cm(升轨)。

(3) 九寨沟M_S7.0地震三维形变结果显示，极震区地表位移在垂直向上方向达23.4cm，垂直向下方向达17.6cm；北向位移达141.8cm，南向位移达100.2cm；东向位移达22.0cm，西向位移达10.7cm。水平向形变表明树正断裂两端地块呈水平相对运动，上盘沿断层面向东向形变较剧烈，符合九寨沟地震左旋走滑的运动特征。

致谢: 审稿专家对本文提出宝贵的修改意见，欧洲航天局提供了Sentinel-1数据，环球星云遥感技术有限公司对Radarsat-2数据处理提供了技术支持，在此一并表示感谢。