利用光学影像相关获取2019年Ridgecrest地震序列同震形变

2021-07-13王乐洋邹阿健许光煜

测绘工程 2021年4期

王乐洋,邹阿健,许光煜

(东华理工大学测绘工程学院，江西南昌 330013)

地震从古至今对人类的危害都十分巨大，而地震同震形变是地震最直观的表现，能用其评估地震灾害等级，研究断层构造和震源机制。获取同震形变的方法有很多，目前主要还是使用GPS技术和InSAR技术[1-4]。GPS技术获取同震形变在远离板块边界和活跃断裂带地区GPS连续观测站分布较为稀疏，空间分辨率较低，不能很好地体现地震产生的地表形变。在获取有地表破裂的地震同震形变时，InSAR技术容易受相位失相干等因素的影响，从而造成靠近断层区域形变数据的缺失。而光学影像在气候良好的情况下，对于监测较大地表形变能够很好地弥补两者的不足。

随着科技的发展，高质量的光学影像逐渐增多，它们可以用来监测由于地质过程、气候变化或者人类活动造成的地表变化。科学家使用光学影像测量地表形变最早是在1991年Bindschadler等用Landsat TM影像观察南极冰川流速[5]。2000年Puymbroeck等首次将光学影像相关应用在地震形变监测中，他们使用SPOT影像获取了1994年Landers地震形变[6]。光学影像的分辨率和质量，成像系统引起的几何畸变，光学影像相关(Optical image Correlation，OIC)技术实行的困难等都是光学影像地表形变监测所需要解决的问题。随着科技的进步，影像处理方法的改进和光学影像的配准和相关(Co-registration of Optically Sensed Images and Correlation,COSI-Corr)工具的开发[7]，这些问题已经得到解决。文献[8]介绍了从一对卫星图像中准确提取地表位移的方法流程，COSI-Corr实现了这些流程。在此之后，使用光学影像监测地表形变变得相对容易，且精度更高[9-10]。

2019年7月6日3时19分(UTC)，美国加利福尼亚Ridgecrest东北处发生Mw7.1地震(震中位于35.770°N，117.599°W)，震源深度8 km，该地震发生前两天(7月4日)发生了一次Mw6.4地震，此次事件被称为Ridgecrest地震序列。尽管加州地区经常发生地震，但已经有20年未发生过如此强烈的地震，根据USGS(United States Geological Survey)地震目录记载，1999年1月1日至2019年7月31日。加利福尼亚东部发生了7次6级以上地震，其中2019年Ridgecrest地震序列是1999年Hector Mine Mw 7.1地震以来震级最大、最强烈的地震。本文首先分析哨兵2号影像相关结果中表现为明显横向条带的姿态角误差的处理，在已有方法上进一步改进；再使用哨兵2号影像相关获取2019年Ridgecrest地震序列同震地表形变，分析地震的地表形变特征；最后对Landsat 8和Landsat 7影像分别进行相关，并与哨兵2号结果进行比较。

1 研究区域和数据处理

1.1 研究区域构造背景

美国加州位于太平洋板块和北美洲板块之间的交汇处，活跃的板块运动导致了加州成为美国地震最活跃的地区。持续的大地测量观测数据表明太平洋板块相对于北美洲板块以约48 mm/a的速度向西北移动[11]。2019年Ridgecrest地震序列的构造背景见图1。

图1 2019年Ridgecrest地震序列构造背景图

图中黑线为板块运动碰撞的边界线，白色方框是哨兵2号影像覆盖区域，黄色五角星表示的是2019年7月4日发生的Mw6.4地震，蓝色五角星表示的是2019年7月5日发生的Mw7.1地震(震源机制球来自USGS)，白色五角星表示的是1999年10月16日发生的Hector Mine Mw7.1地震。红色圆点表示地震前后三周发生的大于4级的地震。黄色正方形表示Ridgecrest所在位置。

1.2 实验数据处理

本文研究使用欧洲空间局(European Space Agency,ESA)的哨兵2号第八波段影像数据，以及美国国家航空航天局(National Aeronautics and Space Administration,NASA)的Landsat 7和Landsat 8卫星第八波段影像数据，影像对的相关信息见表1。

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

使用基于ENVI平台的COSI-Corr软件包，通过光学影像相关技术获取同震地表形变场。由于欧空局的哨兵2号数据是L2A级数据，USGS的Landsat 7和Landsat 8数据是L1C级数据，都是经过粗略的几何校正和正射校正的产品，预处理只需保证影像包含整个震区以及影像对大小一致，相关参数设置如下: 配准运算窗口为32像素×32像素，偏移量步长为6像素×6像素(最终的分辨率结果为60 m)；掩膜阈值设为0.9；迭代次数为2次。哨兵2号影像对(2019-06-13—2019-07-03)相关的结果见图2。其中，图2(a)中正值代表东向形变，负值代表西向形变；图2(b)中正值代表北向形变，负值代表南向形变。

图2 哨兵2号影像对(2019-6-13—2019-7-03)得到的未处理误差的东西向和南北向形变场

如图2所示，处理数据得到的水平形变场中存在各种误差，其中包括失相关噪声、轨道误差、条带误差、卫星姿态角误差、地形阴影误差。由于地形阴影误差较为复杂，本文未进行考虑。为减弱或消除以上误差带来的影响，一般情况下采用以下方法：对于失相关噪声，在软件中掩膜掉形变场中信噪比(Signal Noise Ratio，SNR)低于经验阈值0.9的区域，将过大和过小的值替代为空值(Nan值)。由于下载的影像没有进行严格的校正，得到的形变场中依旧会存在轨道误差，对于轨道误差，用一次多项式曲面拟合来消除。在卫星成像过程中，传感器的电荷耦合器件(Charge Coupled Device，CCD)阵列的错位引起卫星条带伪影即条带误差[12]，对于条带误差，旋转影像使条带呈竖直分布，对列采用均值相减法[10]。在结果中由于卫星成像过程中姿态角的变化，导致形变结果中出现卫星姿态角误差，对于卫星姿态角误差，一般表现为微弱的横向带状信号，与条带误差相互垂直，旋转相同角度使得条带呈水平分布，再对行采用均值相减法。最后再用3像素×3像素的中值滤波窗口进行进一步降噪。

对于哨兵2号光学影像相关的结果，姿态角误差表现为明显的横向条带，使用传统的均值相减法并不能去除，文献[10]提出了改进的均值相减法，将研究区域分成12等份再进行均值相减，很好地去除了凯库拉地震区域哨兵2号影像相关结果中的横向条带，但使用此方法对本次加州地震区域哨兵2号影像相关结果中横向条带的去除效果不是很理想(见图3(a))，本文对其进行进一步的改进。由于地形起伏，植被覆盖，存在云和河流湖泊等因素，形变场结果中各区域相干性有差别，导致传统的均值相减并不能完全去除横向条带，所以把结果分成若干区域，增强均值相减区域的相干性，再对各区域进行均值相减来去除条带。从图3(a)中可以看出，虽然将横向条带分成12份使得各自区域内相干性变强，但处理后还存在一些横向的条带，仔细观察发现这些未去除的条带存在于各竖向条带的两旁，可能是由于各竖向条带两旁区域相干性较差，但两竖向条带之间区域的横向条带相干性较好。基于这点，将横向条带按照竖向条带所处位置进行划分，然后对划分后的各区域进行均值相减，并按照去除情况对竖向条带内的区域进行一定数量的均分。

图3 哨兵2号影像对(2019-6-13—2019-7-03)得到的处理完误差的东西向和南北向形变场

由于结果中各处相干性有差别，各处的值是不均匀的，很难依靠旋转后的列均值比较等方法来分辨图中各竖向条带所处的位置。本文根据旋转前两竖向条带顶端经差来确定两竖向条带之间区域所占比例，确定旋转前形变场竖向条带顶端各点在形变场所代表矩阵中的位置，然后得到旋转后形变场中竖向条带的大致位置。相关过程可表示为：

(1)

式中:L1，L2表示旋转前第一个和第二个竖向条带顶点的经度；l表示旋转前第一个和第二个竖向条带顶点的经度差；k表示旋转前形变场所代表矩阵的列数；L表示整个结果图中最右端与最左端的经度差；k1表示旋转前两个竖向条带顶点处的列数差；round(·)表示按照四舍五入原则将结果转化成整数；按旋转前形变场设立坐标系，令左下角坐标为(0,0)，(x0,y0)表示形变场的中心位置坐标；(x1,y1)表示旋转前第一个竖向条带顶点的坐标，(x2,y2)表示旋转后该点的坐标；b表示逆时针旋转使得横向条带呈水平分布的角度。按求得k1相同的方法可得到形变图中竖向条带所划分的各区域的列数差，从而得到旋转前各竖向条带顶点在所设坐标系中的坐标，根据求得(x2,y2)相同的方法可以得到旋转后各竖向条带顶点坐标。

按照上述方法得到两个竖向条带旋转后的横坐标之差可以得到两竖向条带的间距，根据此间距可以对旋转后的形变场进行划分，再对每个区域进行均值相减。本次事件覆盖区域哨兵2号影像相关结果旋转后最左端至第一个竖向条带间距与两竖向条带间距相同，将其与其他4个两竖向条带之间区域形变场提取出来，按照下述过程去除横向条带。

(2)

式中:φ表示横向条带误差水平；δ表示原形变场所代表矩阵；β表示减去横向条带误差的形变场所代表矩阵；m,n表示提取出的形变场所代表矩阵的行数和列数；a,b表示分成的每个列等份的初始像素点和最后的像素点在矩阵中的位置；Dk表示该像素点的值；ceil(·)表示大于某一整数的最小值;c表示将形变场分成5个列等份，可以按照去除情况对各区域内部进一步均分，即可以将其分成5的倍数等份，但不宜过多。具体流程见图4，结果见图3(c)、图3(d)。

图4 本文方法去除横向条带数据处理流程

在图3(c)、图3(d)中，没有看到明显的横向条带，从视觉上看去除效果较为良好。选取竖向条带旁区域作剖面线(见图3(a)黑色线条)，其各点误差值的情况见图5。可以看出东西向形变场中直接均分的均值相减法得到的值明显较为分散，而按本文方法则基本集中在-0.2～0.2 m之间，两者精度分别为0.215 m和0.133 m，本文方法改善较大。而南北向形变场中可能由于其竖向条带两旁相干性差距不大，直接均分的均值相减法得到的结果中也没有看到明显的横向条带，两种方法精度分别为0.148 m和0.138 m，从精度上看，本文方法略好于直接均分的均值相减法。

图5 直接均分的均值相减法与本文方法得到的形变场剖面

2 同震地表形变场特征分析

在2019年Ridgecrest地震序列研究中，本文采用基于ENVI的COSI-Corr软件包，通过光学影像相关技术，使用震前震后的哨兵2号影像对(2019-06-28—2019-07-18)，使用前节所述处理方法获取该次地震的同震地表形变场如图6(c)、图6(d)所示，具体破裂区域如图7所示。从图中可以看出Ridgecrest Mw7.1地震的破裂迹线从东南到西北呈近直线形状，Mw6.4地震的地表形变主要表现在Mw7.1地震破裂的东南方向，破裂迹从东北到西南呈近直线形状，两次破裂相互交叉。在东西形变场中，Mw7.1地震的破裂迹以北形变值为正，向东位移，破裂迹以南为负，向西位移，两侧形变较为平均，呈现拉伸趋势，最大形变值在1.7 m左右；Mw6.4地震的破裂迹东西两侧均为负。南北形变场中，Mw7.1地震破裂线北侧形变为负，向南位移，而南侧形变为正，向北位移，南侧形变主要集中在中部，北侧形变比较平均，北侧形变明显高于南侧，呈挤压趋势，最大形变值在2.3 m左右；Mw6.4地震的破裂迹东西两侧均为正。Mw6.4破裂迹两侧形变情况可能受在其之后发生的Mw7.1地震的影响，导致其破裂迹两侧位移方向相同，但还能看到明显破裂迹线，说明虽然位移方向相同，但两侧大小还是有一定程度的落差。在结果图中还可以看到在主破裂西方向的北部和南部均有一个小破裂。形变结果不仅能提供形变数据，还可以提供地震的断层几何参数。从本次地震的形变结果来看，可以得到该地震序列的主破裂的走向角大约为320°，破裂长度约为45 km。次破裂的走向角约为230°，破裂长度约为15 km。相关文献中使用PlanetScope影像[13]，worldview影像[14]得到的Ridgecrest地震序列同震形变场结果，在形变大小的量级，以及断层迹线细节上均与本文相似，几者间存在着良好的一致性，可见本文结果较为合理。

图6 哨兵2号影像获取的Ridgecrest地震序列的同震形变场

图7 哨兵2号数据获取的Ridgecrest地震序列的同震形变场具体区域

本文通过选取远场部分区域(见图6(c)黑色方框)，计算其标准差作为精度评定的指标。由于结果中的大部分误差都已经去除，理论上选定区域的形变值为0。但误差的处理会受到软件亚像素相关过程的精度的限制，光学影像用COSI-Corr进行相关性匹配监测形变精度至少为1/20像素[8]，对于哨兵2号结果来说形变监测精度为50 cm。按图6(c)方框分别选取范围内的东西向和南北向形变，得到它们东西向和南北向的均方根误差分别为0.287 m和0.3 m，精度水平比较合理[9-10]。

加州地区有数以千计的GNSS站，但是尽管有这么多精度较高的GNSS站，但其空间覆盖率仍较低。目前Ridgecrest 地震序列周围地区GNSS站点之间的平均间距为20～30 km[15]，与加州西部板块边界区域相比，本次地震周边的GPS站点更少，站点之间的间距更远，与这次地震的断裂长度相当，不足以用来研究本次地震形变的细节。如果相邻像素间的应变梯度超过一个相位周期，InSAR测量通常会在破裂迹线附近产生失相关，相关文献中哨兵1号与ALOS-2 PLAYSAR-2数据得到的Ridgecrest地震序列地表形变，在地表破裂迹线附近区域失相干情况较为严重，且不能清晰的看到前震破裂和主震破裂之外的破裂迹线[16]。相较于GPS和InSAR结果，哨兵2号影像相关的结果虽然损失了一定精度，但它既具有较高的空间分辨率，又能够弥补InSAR结果中失相干的情况，更能体现该次地震的同震形变情况。

3 Ridgecrest地震序列的不同光学影像相关结果

为了探究使用其他开源光学影像获取的本次地震的形变场情况，本文选取Landsat 7和Landsat 8两种光学影像分别进行相关，其结果见图8和图9。

图9 Ridgecrest地震序列的Landsat 7影像相关结果

从图6和图8中可以看出，Landsat 8与哨兵2号相关结果中均能明显看出地表破裂迹线，均能得出大致断层走向和长度。但Landsat 8相关结果中破裂迹线没有那么清晰，整个形变场中异常值更多，存在的噪声误差也更加明显。选取图中黑色方框区域，东西向和南北向均方根误差为0.345 m和0.334 m，显示效果和精度均稍差于哨兵2号结果。

图8 Landsat 8获取的Ridgecrest地震序列处理完误差的同震形变场

从图8和图9可以看出，Landsat 7结果中两图均未发现Landsat 8结果中出现的明显破裂迹线。在Landsat 8结果中，东西向形变场和南北向形变场中异常值和空值主要集中在瑟尔斯湖以及一些植被茂密、地形起伏的区域；在Landsat 7结果中，只有东西向形变场中充满空值的瑟尔斯湖区域与Landsat 8结果同样明显，但其他区域充满了异常值，在南北向形变场中则是充满了空值。由此可见在本次地震中，Landsat 8结果与Landsat 7结果存在着较大的差距，且用Landsat 7影像相关无法获取本次地震的同震形变。从结果可以看出，光学影像的分辨率对结果的显示效果和精度有一定的影响，影像分辨率越高，相关结果的显示效果和精度越好。而即使是同种分辨率的光学影像，结果也会有差别。Landsat 7和Landsat 8同样为15 m空间分辨率，但两者的传感器不同，Landsat 7携带增强型专题制图仪(Enhanced Thematic Mapper, ETM+)传感器，而Landsat 8则携带陆地成像仪(Operational Land Imager, OLI)和热红外传感器(Thermal Infrared Sensor, TIRS)，两者发射时间相差14年，各项技术Landsat 8均较为成熟。与Landsat 7第八波段相比，Landsat 8的第八波段范围较窄，可以更好地区分植被和无植被特征，数据较好。而且2003年之后Landsat 7卫星受到损坏，导致之后的数据产生缺失，虽然已有方法对其数据进行补充，但质量依旧与损坏前有一定的差距。由于Landsat 7影像数据质量相对较差，横向的条带误差大小间隔过大，地形起伏、河流湖泊、植被覆盖等因素造成的失相关情况更加严重，导致误差不能很好地削弱。2021年3月将发射的Landsat 9卫星将为光学影像地表监测提供更加有效且丰富的数据。