2019-07-04美国加利福尼亚州发生M_W6.4地震，34 h后该区域又发生M_W7.1地震，震中位于小湖断裂带及Airport Lake断裂带交汇区域^[1]，该地区构造活动强烈，曾发生多次7级以上地震，包括1952年克恩县M_W7.5地震、1992年兰德斯M_W7.3地震及1999年M_W7.1赫克托矿地震等^[2]。地震发生后，许多学者利用InSAR、GPS数据提取了本次地震的形变场，但由于加州地震的震中区域地表破坏严重，D-InSAR结果在震中区域形变梯度过大，相位解缠失败，无法提取准确的断层破裂轨迹^[3]；此外，基于升降轨InSAR数据的offset-tracking方法虽可避免失相干影响提取方位向、距离向形变场，但在三维解算同震三维形变场时仍受D-InSAR结果控制；GPS数据由于点位稀疏，在地震等大区域形变监测领域受到很大限制。相对而言，光学遥感数据具有覆盖范围全面、分辨率高、获取方便等优势，其局限性在于只适合监测水平向形变，对地表垂向沉降不敏感，同时形变测量精度受限于影像匹配精度及像元分辨率，但在监测以水平位移为主、地表破裂明显的形变领域，高分辨率光学影像同样可以提供高精度地表形变场。目前，该方法广泛应用于冰川、滑坡及同震位移形变场监测领域^[4-7]。

本文利用地震前后的Landsat-8和Sentinel-2光学影像数据提取加州地震同震形变场，通过光学形变后处理算法去除系统误差项，并利用形变误差的方差-协方差函数及融合offset-tracking技术的InSAR同震三维形变场综合验证光学形变测量的精度，最终根据形变信息提取地震的发震断层几何展布及水平位移信息，包括断层长度、最大水平位移量等参数，为该地震的震源参数反演及发震背景研究提供依据。

1 数据处理与方法 1.1 光学影像处理

分别选用地震前后的两景Landsat-8和Sentinel-2光学影像第8波段数据，数据来源于USGS网站，参数如表 1所示。原始影像如图 1所示，选取的两类影像含云量均小于3%，时间基线分别为16 d及10 d，可有效避免地物变化造成的时间失相干及云雾所引起的噪声。Sentinel-2光学数据经大气校正及辐射定标处理，Landsat-8数据属于1级产品，由于两者数据覆盖范围不一致，为减少运算量，将Landsat-8影像裁剪至相同范围。最后，以COSI-Corr为数据处理平台，基于频率域互相关算法，通过逆傅里叶变换计算两景影像的相对位移，经亚像元匹配可达到1/20亚像元精度，即同震位移形变场精度可提高20倍^[8](Landsat-8为0.75 m，Sentinel-2为0.5 m)。假设i₁、i₂分别为震前、震后两景影像，其相对位移量为Δ_x、Δ_y，对其进行傅里叶变换后存在以下关系式：

$ {i_2}(x, y) = {i_1}(x{\rm{ }} - {\Delta _x}, y{\rm{ }} - {\Delta _y}) $

(1)

$ {I_2}({w_x},{w_y}) = {I_1}({w_x},{w_y}){{\rm{e}}^{{\rm{ - j}}({w_x}{\Delta _x} + {w_y}{\Delta _y})}} $

(2)

式中，I₁、I₂为两幅影像的傅里叶变换，w_x、w_y分别为行和列的频率变化。两幅影像的归一化互谱C_{i1, i2}为：

$ \begin{array}{c} {C_{{i_1}, {i_2}}} = \frac{{{I_1}({w_x}, {w_y})I_2^*({w_x}, {w_y})}}{{|{\rm{ }}{I_1}({w_x}, {w_y})I_2^*({w_x}, {w_y})|}} = \\ {{\rm{e}}^{{\rm{j}}({w_x}{\Delta _x} + {w_y}{\Delta _y})}} \end{array} $

(3)

式中，符号*表示复数共轭相乘。基于该函数进行逆傅里叶变化可得到与位移量有关的互相关函数$\mathscr{F}^{ - 1}$，之后取二维脉冲函数的峰值点即可得到两景影像相对位移及信噪比信息：

$ {\mathscr{F}^{ - 1}}\{ {{\rm{e}}^{{\rm{j}}({\omega _x}{\Delta _x} + {\omega _y}{\Delta _y})}}\} = \delta (x{\rm{ }} + {\Delta _x}, y{\rm{ }} + {\Delta _y}) $

(4)

图 2为Sentinel-2、Landsat-8融合EW向、NS向形变场及信噪比的三波段融合影像，其中area1、area2、area3为失相干区、轨道区及条带误差区，area4为本文选定的精度评定区域。

1.2 误差后处理

1) 光学影像特征追踪法通过计算同名点像元值的偏移量提取形变信息，在植被、建筑、雨雪覆盖区域影像DN值随时间变化剧烈，容易丢失相关性。如图 2(a)和2(b)所示，area1位于Ridgecrest小镇湖泊及农田附近，受时间失相干影响，在远离震中区域仍存在大量异常形变值，在形变图上表现为信噪比均低于0.9。本文通过掩膜处理去除整景形变场中存在的低相干区域，并对该区域进行插值平滑处理。

2) 卫星成像过程中由于拍摄位置的变化，形变场全局存在一个固有偏移量，即轨道误差。如图 2(a)所示，受轨道误差影响，形变场存在一个明显的线性趋势，常采用多项式曲面拟合模型方法去除轨道误差^[9]。在非形变区域建立参考点，根据最小二乘原理解算待求参数，用原始形变场减去模拟的整幅图像轨道误差项，即可得到地表真实形变场：

$ {\varphi _{{\rm{orbit}}}} = {\rm{ }}{a_0} + {a_1}x{\rm{ }} + {a_2}y{\rm{ }} + {a_3}xy $

(5)

式中，φ_orbit为轨道趋势项误差，x为水平向坐标，y为垂直向坐标，a₀、a₁、a₂、a₃为待求参数。

3) 推扫式卫星成像过程中传感器保持固定，CCD(电荷耦合器件)线性阵列探测器所导致的非线性失真小于1 mm，在形变监测的精度范围之内，因此误差主要来源于飞行平台抖动导致的影像几何变形。如图 2(b)所示，由于Landsat-8、Sentinel-2均为处理后的一级产品，无法获取传感器轨道参数及姿态参数进行校正，沿形变场交叉轨道方向存在大量平行排列的带状条纹。常规的均值相减法容易引入地形、植被造成的异常形变值，因此本文采用改进均值相减法^[10]，将旋转后的形变图均分成12等份，使条带误差呈竖直分布，再以每一等份的均值代表该行像素所受到的姿态误差影响水平，最终减去该误差即可去除形变场中带状条纹。

处理后的影像如图 3所示，整体形变场信噪比均高于0.96，Sentinle-2的EW向形变场中异常低值及Landsat-8影像中残存的条带误差均被有效去除。由于地震发生于巴拿明山谷及死亡谷附近，受地形误差影响，在Sentinel-2形变场左上区域仍存在部分异常低值；Landsat-8影像由于只经过粗略的辐射校正，且数据质量、分辨率较低，同时卫星姿态误差并未完全去除，远场区域存在大量异常形变值。

2 同震形变场分析与评定 2.1 地表形变场分析

通过数据处理流程，可获取整个研究区EW向及NS向的同震形变场。如图 3所示，加州地震由AB、CD两段地表形变场中的断层破裂迹线组成，其中Landsat-8因受分辨率及数据质量制约，远场区域存在异常形变值，但所揭露的断层轨迹与Sentinel-2一致。在地表形变场中，AB段发震断层迹线整体呈NNW向展布，断层南端距加洛克断层仅3 km，在南部三分之一处走向转变为NW，由南向北逐渐延伸。该段发震断层最大水平滑移量为2.82 m，长度达55 km，且上盘形变明显大于下盘，在远离断层方向形变值逐渐减小，断层呈右旋走滑运动特征，推测M_W7.1主震是由NNW向的断层运动所引起的。距NNW破裂迹线南部18 km处存在一个与之相交的NE向破裂迹线，该断层迹线整体走向NE，最大滑移量为1.05 m，地表破裂轨迹15 km，运动学特征以左旋走滑为主，推测M_W6.4前震是由NE向断层运动所致。由于本次地震震源机制复杂，地震序列中包含多达20个2 km以上的断层破裂，受到M_W7.1主震及余震影响，NS向形变场中CD段断层上、下两盘均表现为正N向运动，但断层上盘形变值远低于下盘；EW向形变场中断层上下两盘仍呈相对EW向运动趋势，说明实际断层呈现走滑特征，滑动方向为左旋。

Fielding等^[3]利用升降轨ALOS-2和Sentinel-1影像，通过联合offset-tracking及升降轨D- InSAR结果解算加州地震三维形变场，最后基于GNSS结果验证同震形变场精度，提取出的形变场EW向形变范围为－1.5~1.5 m，NS向形变范围为－2~2 m。该结果与本文利用光学影像所提取的同震形变场具有很好的一致性，同时由于D-InSAR提取的LOS向形变在震中区域出现相位缠绕，沿发震断层附近区域存在大量空值区，本文所提取的断层轨迹可为后续震源参数分析及滑动模型反演提供约束条件。

2.2 精度评定

为估算形变场中存在的误差，本文选取远离形变区分布的area4开展精度评定，该区域远离震中分布，去除误差项后实际形变值应为0。根据地统计学中的变异函数理论，分别采用方差、协方差函数对该区域形变场取样，以推算整体形变场的误差水平^[11]。其中，变差函数及协方差函数可表示为：

$ \gamma ({h_i}) = \frac{1}{{2N}}\sum\limits_{i = 1}^N {{{[f({x_i}) - f({x_i} + {h_i})]}^2}} $

(6)

$ C({h_i}) = \frac{1}{{2N}}\sum\limits_{i = 1}^N {f({x_i})\cdot f({x_i} + {h_i})} $

(7)

式中，x_i为影像中随机采样点的位置，f(x_i)为在影像中随机采样点x_i上的形变值，h_i为距离间隔，N为给定距离h时采样点对的数目。

从图 4可以看到，EW向形变协方差在15 km处逐渐趋于稳定，NS向协方差在20 km才开始稳定，表明形变场在大于20 km处的两点空间上不存在相关性。由于Sentinel-2光学影像未借助卫星辅助姿态文件进行校正，同震形变场始终存在因地形因素造成的误差形变值，在协方差曲线上反映为影像存在微弱的相关性。根据变差函数显示，EW向及NS向的方差分别为35.17 cm和74.96 cm，说明NS向形变场受地形影响更为明显。

3 结语

本文以2019年加州M_W7.1和M_W6.4地震为例，详细论述Landsat-8、Sentinel-2光学影像同震形变数据处理流程及误差后处理方法，开展了地表同震形变场分析，并结合InSAR同震三维形变及形变误差的方差-协方差函数综合验证光学形变精度。结果显示，两次地震的发震断层分别为NNW向的左旋走滑断裂及NE向的右旋走滑断裂，形成了共轭关系，其中M_W6.4前震以左旋走滑为主，最大滑移量约为1.05 m，断层长度为15 km；M_W7.1主震以右旋走滑为主，地表最大滑移量达到2.82 m，断层长度约55 km。由于D-InSAR结果在缺少断层近场数据的情况下常导致错估滑动分布^[12]，利用GPS融合InSAR的方法常导致断层反演的最大滑动量小于实际值^[13]。基于高分辨率光学影像获取的同震形变场可有效避免D-InSAR形变梯度过大时出现的失相干现象，准确揭示该地震的地表形变、发震断层几何展布和水平位移信息，从而提高加利福尼亚地震序列的断层滑动分布和震源参数反演的准确性，并为本次地震的发震断层及运动学背景分析提供依据，同时为国内开展基于光学遥感影像的地表同震形变场研究奠定基础。

致谢: 本文所采用的Landsat-8、Sentinel-2光学影像来源于美国USGS，在此表示感谢。