云南漾濞5.21地震诱发地表形变的D-InSAR快速探测分析

朱智富，甘淑，袁希平，张晓伦，黄志豪

(1.昆明理工大学国土资源工程学院，云南 昆明 650093；2.云南省高校高原山地空间信息测绘技术应用工程研究中心，云南 昆明 650093；3.滇西应用技术大学云南省高校山地实景点云数据处理及应用重点实验室，云南 大理 671006)

1 引 言

根据中国地震台网(http：//www.ceic.ac.cn/)报道，2021年5月21日，中国云南省大理州漾濞县(北纬25.67°，东经99.87°)发生MS6.4级地震，震源深度 8 km。本次地震灾害发生在夜间，全州十二县市均有明显震感，目前漾濞县受灾最严重，大量房屋遭到破坏，多处道路中断，造成了严重的经济损失。据历史地震目录记载，漾濞地区处于红河断裂带[1]，属于地震多发地区，在该断裂带的中段，即下关到漾濞这一段，历史上发生的地震比较多，因此，针对漾濞震后的地表变化进行研究具有现实指导意义。

传统的地表形变监测方式，由于地形条件和仪器自身局限性的影响，只能获取到研究区域内有限的地表信息，不能在短时间内，得到大范围、监测点连续的形变信息，并且这些监测点都是以离散点的形式存在，不能很好地反映出研究区域的地表变化。与之相比，差分合成孔径雷达干涉测量技术(Differential Synthetic Aperture Radar Interferometry)具有全天时、全天候、大面积和高精度的优势，能够监测厘米级甚至毫米级的地表形变，可以为地震监测提供空间连续地表形变信息[2]。目前，D-InSAR技术已经成为监测地震地表形变的重要手段。文献[3]利用合成孔径干涉测量对九寨沟震后的滑坡隐患点进行识别与形变监测，共探测出13处滑坡隐患，验证了基于InSAR技术对滑坡隐患点进行早期识别的可靠性。文献[4]利用D-InSAR技术对2017年伊朗地震进行反演，获取了发震断层的精确几何参数和最优断层滑动分布。文献[5]利用D-InSAR结合GIS技术对矿区开采的沉陷形变进行监测，很好地反映出煤炭在开采过程中的沉降区域位置和分布范围。文献[6]利用二轨法D-InSAR技术对北京城区地面沉降进行监测，发现北京市地面沉降较为明显的两个沉降中心，表明了D-InSAR技术在城市地区的地面不均匀沉降方面具有明显优势。文献[7]利用升降轨D-InSAR技术对新疆伽师地震的同震形变场进行提取，成功获取了此次地震的LOS向形变。基于D-InSAR的一系列优势，本文将利用二轨法D-InSAR技术对漾濞震后的地表形变进行快速探测分析。

2 研究区及数据源

2.1 研究区概况

研究区位于漾濞彝族自治县如图1所示，隶属于云南省大理州，位于大理州中部，东与大理市、巍山彝族回族自治县毗邻，西与永平、云龙二县接壤，南交保山市昌宁县，北连洱源县。境跨北纬25°12′～25°54′、东经99°36′～100°07′，属亚热带和温带高原季风气候区，地势北高南低，最高海拔 4 122 m，最低海拔 1 174 m。2021年5月21日，漾濞县(北纬25.67°，东经99.87°)发生MS6.4级地震，震源深度 8 km。

图1 研究区及其裁剪后的DEM数据

2.2 实验数据

本次实验数据来自开源的欧空局提供的Sentinel-1A C波段SAR数据(https：//scihub.copernicus.eu/)，模式为干涉测量宽幅模式(IW)，其极化方式为VV极化，距离向和方位向分辨率为 5 m×20 m。实验数据的具体参数如表1所示。为了提高轨道精度，采用阿拉斯加(https：//s1qc.asf.alaska.edu/)提供的精密AUX_POEORB轨道数据文件对原始的SAR数据进行轨道精炼。采用美国太空总署(NASA)和国防部国家测绘局(NIMA)联合测量的空间分辨率为 30 m的STRM1-DEM数据，作为参考高程数据，用于地形相位的去除。

Sentinel-1A干涉参数 表1

3 研究方法与技术

3.1 D-InSAR探测技术

InSAR技术是利用雷达系统获取同一地区两幅SAR影像所提供的相位信息进行干涉处理，来获取地表的形变[8]。D-InSAR是在传统InSAR基础上发展起来，利用同一地区不同时刻获取SAR数据进行差分干涉处理，通过差分处理去除两次观测相位中的平地效应、地形相位和大气延迟等，得到形变相位，进而获得地表形变信息[9]。目前的D-InSAR方法包括了二轨法，三轨法和四轨法，这主要是由所选影像的数量来决定的。综合对比不同方法性能差异，二轨法干涉方法简单而最为可靠，因此在本文实证研究中采用了二轨法。

研究采用的二轨法D-InSAR技术的基本原理及流程如图2所示。在形变监测中，二轨法D-InSAR的主要技术核心是使用跨越形变期的两幅SAR影像和该区域的外部DEM来完成差分干涉[10]。

图2 二轨法D-InSAR原理示意图

将两幅SAR影像配准重采样并共轭相乘形成一幅干涉图，干涉图中每个像元的涉相位包含该点相对于参考位置的高程和形变信息。考虑到大气环境和地表的影响，干涉图中的相位φint可以表示为：

φint=φref+φtopo+φdef+φatm+φnoise

(1)

在(1)式中，φref表示参考椭球相位，φtopo表示地形相位，φdef表示形变相位，φdef表示差分干涉相位，φatm表示为大气相位，φnoise表示噪声相位[11]。

特别需要说明的是，由于后面两项的大气和噪声的影响是无法消除的，所以这里可以不再考虑大气和噪声的影响，从而相位φint可简化为：

φint=φref+φtopo+φdef

(2)

通过轨道数据和外部DEM去除参考椭球相位和地形相位后，便可得最终的形变相位φdef，表示为：

(3)

其中λ为波长，△r为传感器斜距方向的形变量。

3.2 地震前后D-InSAR数据处理

本次实验采用了二轨法D-InSAR差分技术，利用SARscape软件对震前和震后的SAR数据进行处理，其基本的处理技术流程如图3所示，包括：①原始影像的预处理；②基线估算；③干涉相位图生成；④滤波与相干性计算；⑤相位解缠；⑥轨道精炼与重去平；⑦相位转形变及地理编码。

图3 D-InSAR处理流程

试验研究中的具体操作步骤如下：

(1)原始影像的预处

每景Sentinel-1影像数据覆盖范围为 400 km×400 km，而研究区域远远小于一景影像的覆盖范围，为了节省时间和提高处理数据的工作效率，首先对原始影像进行裁剪，得到震前和震后的研究区影像。

(2)基线估算

这一步主要是估算主从数据的基线情况，包括时间基线、空间基线、多普勒偏移、一个相位周期变化代表的高程变化等信息[12]，其中空间基线的垂直分量不能超过临界基线的阈值，否则就会丢失相位信息，从而无法进行干涉。因此，垂直基线是判断所获取的SAR数据能否做干涉的重要依据。

对裁剪后的研究区进行基线估算，得到空间基线为 49.061 m，远远小于临界基线 6 262.828 m，满足D-InSAR处理要求。

(3)干涉相位图生成

干涉相位图生成主要是基于轨道信息将两景影像进行配准，然后用SLC数据所对应的像元进行复数共轭相乘，所得到的干涉图像(图4)包含了平地效应的影响，接下来引入外部DEM数据去除平地效应[13]，得到去平后的干涉图像(图5)。

图4 干涉图 图5 去平后的干涉图

(4)滤波与相干性计算

经过干涉处理后的图像中还包含形变相位、大气延迟相位和噪声相位所带来的影响，这些影响在差分干涉图中表现为周期不明显，条纹不清晰和斑点噪声，为后续的相位解缠带来困难[14]。通过相位滤波，削弱噪声信号，提高相位解缠的准确性。常用的滤波方法为Goldstein自适应滤波，该方法可以在保证条纹清晰度的情况下有效的滤出噪声，经过SARscape软件处理后得到滤波后的干涉图像(图6)和相干性图(图7)。相干性系数的范围为0-1，相干性系数越大，相干性越高，解缠效果也越好，反之亦然。

图6 滤波后的干涉图 图7 相干性图

(5)相位解缠

差分干涉图中的形变相位被缠绕在[-π，π)区间内，不能直接反映出真实的地表形变信息。经过相位解缠处理可以恢复每一像素的相位整周数，得到连续的地形变化[15]。本实验采用了SARscape软件提供的Minimum Cost Flow(最小费用流法)进行相位解缠处理，得到解缠后的干涉图像(图8)。该方法可以考虑图像上所有的像元，对相干性小于阈值的像元做掩膜处理。

图8 解缠后的干涉图

(6)轨道精炼与重去平

以滤波后和解缠后的干涉图为参考，选取GCP地面控制点。利用选取的GCP点进行轨道精炼和相位偏移计算，消除可能的相位坡道和残余的轨道相位，此时就只剩下形变的相位。在GCP控制点的选取过程中要远离干涉条纹密集和解缠错误的区域，选择相干性较高和平地区域。本实验选取了11个GCP点，用Automatic Refinement(自动优化法)进行轨道精炼与重去平。

(7)相位转形变及地理编码

把经过绝对校准和解缠的相位，结合成合成相位，将仅包含形变信息的差分干涉相位转为雷达视线方向的形变量[16]，然后通过地理编码，把SAR坐标系下的形变图转化到WGS-84坐标下，最终得到LOS向的地理坐标形变图(图9)。

图9 LOS方向上的地表形变图

4 D-InSAR形变监测结果分析

4.1 形变监测结果的沉降凸显斑块识别分析

根据图9的形变空间分异所示，首先可以看出此次地震造成4处形变比较明显的区域，且4个区域呈现出以EE′为分界线的对称分布，其中分界线EE′的西侧区域的形变量为负值(远离卫星方向)，表现为沉降，而东侧区域的形变量为正值(靠近卫星方向)，表现为抬升。实验中为了进行具体对比分析，对4个形变量较大集中区进行多边形圈定如图10所示，并分别命名为区域A、B、C、D，即得出区域A和区域B的最大沉降值为 8.0 cm，区域C和区域D的最大抬升值为 8.8 cm。区域A沉降区域面积较大，但沉降数值相对较小，而区域B虽然沉降区域面积较小，但它的沉降值相对较大，二者在沉降面积和沉降数值上形成鲜明对比。区域C和区域D的抬升面积和抬升数值基本一致，其中抬升面积较沉降区域面积而言，相对较小，说明此次地震造成的地表形变主要还是以沉降为主。最后结合4个形变区与最大震中位置的空间相关关系进行初步解读分析，表明此次地震造成地表形变最大的区域并不是出现在最大震级发生的震源中心位置，而是发生在偏离最大震级震中的东南侧，即形变区A、B、C、D4个区域。

图10 典型形变区的多边形圈定示意图

4.2 形变监测结果与震源空间关联分析

为了进一步对利用D-InSAR地震监测结果的形变空间分异特性开展深化探讨，研究中收集整理了对中国地震台网发布的关于漾濞县在2021年5月21日～5月22日的4级以上地震记录资料如表2。初步统计分析可知，在此期间，漾濞县一共发生了14次4级以上的地震，而将这14次地震点的位置信息与D-InSAR监测的沉降形变图进行空间叠加可得到图11。分析图示空间关联可以发现：记录的14次4级以上的余震中心点位置，有10次分布于沉降A区和B区，有2次分布于抬升C区，亦即有12/14次的频度的4级以上的地震中心分布与利用D-InSAR技术探测出的地表形变较大区域吻合，该结果一定程度上较好地实证了利用D-InSAR技术进行定量、定位而连续探测地震诱发地表形变空间分异的可行性及可靠性。

4级以上地震记录 表2

图11 4级余震分布图

4.3 典型形变监测剖面分析

最后，为了更加深入地揭示D-InSAR形变监测结果空间分异特性，进行不同方向的形变大小量化认识分析，针对D-InSAR形变监测结果空间图层，研究中特别设置了与分界线EE′近似平行的，分别穿过沉降B区和抬升C区的两条剖面线BB′和CC′，以及一条与分界线EE′近似垂直且穿过抬升D区中心的剖面线B′C′如图12(a)所示，并绘制出对应剖面线的形变剖面图(图12(b)、(c)、(d))。根据三条典型形变剖面曲线图的解析可以看出：

图12 典型形变区剖面线图

(1)剖面线BB′横跨沉降区域A和B，主要是以沉降为主，剖面线整体呈现凹形，局部呈起伏变化状态。主要集中在 15 km～20 km、22 km～25 km和26 km～28 km之间发生了3次剧烈沉降，形成3个较大的“漏斗型”沉降区，其中沉降的最大值达到 5.5 cm左右。

(2)剖面线CC′主要位于抬升区C，剖面线整体呈现凸形，局部呈起伏变化状态，在剖面线 0 km～6 km之间变化迅速，抬升量达到了 7 cm左右，之后又相对变得平缓，在 11 km～15 km之间虽然曲线表现为下降趋势，但形变值仍然没有出现负值，说明该区域的主体形变仍然还是以抬升为主。

(3)剖面线B′C′横跨抬升和沉降区，呈正负分布，在对称线EE′的西侧表现为沉降，东侧表现为抬升，且抬升趋势远远大于沉降趋势，大约在 13 km位置刚好为沉降与抬升的分界点，可以大概推测它为分界线EE′的位置，剖面线B′C′经过分界线EE′从沉降区B进入抬升区C。剖面线在分界点处没有出现断裂现象，说明此次地震的破裂没有到达地表。在 13 km～18 km之间上升趋势剧烈，达到了抬升的峰值，之后的变化相对比较连续平稳，基本没有“跳点”出现。

5 结 论

本文利用二轨法D-InSAR技术对5.21漾濞地震诱发的山地形变空间分异特征进行快速探测，识别出4个地表形变较大的区域，并且发现地表形变变化最大的位置并不是在最大震级的地震中心，而是分布在最大震级震中东南侧的次级地震频发的空间区域。结合中国地震台网发布的大于4级余震发生的震中位置记录资料整理，空间叠加D-InSAR形变探测结果进行初步关联分析，发现这些点位与D-InSAR探测出的形变集中区域具有较大的空间关联性，一定程度上实证了D-InSAR技术用于震后地表形变快速探测的实用性及可行性。最后，基于D-InSAR监测出的地震后形变空间分异图及所识别出的4个地表形变集中地块，通过三个典型剖面图绘制及曲线特征发现，初步分析认为此次地震造成的主要地表形变变化以沉降为主。