﻿ 基于光学卫星影像的同震形变提取技术及形变梯度的地学意义
 文章快速检索 高级检索
 大地测量与地球动力学  2024, Vol. 44 Issue (5): 485-490  DOI: 10.14075/j.jgg.2023.08.156

### 引用本文

HUANG Yu, LIU Xiaoli, MO Xinyu, et al. Coseismic Deformation Extraction Method and Geodynamic Significances of Displacement Gradients Based on Optical Imagery[J]. Journal of Geodesy and Geodynamics, 2024, 44(5): 485-490.

### Foundation support

Scientific Research Fund of Institute of Seismology, CEA and National Institute of Natural Hazards, MEM, No. IS20226325; The Spark Program of Earthquake Technology of CEA, No. XH22003C; Open Fund of Wuhan Gravitation and Solid Earth Tides, National Observation and Research Station, No. WHYW202201.

### Corresponding author

LIU Xiaoli, associate researcher, majors in disaster reduction by remote sensing and earthquake geodesy, E-mail: liuxl.j@163.com.

### 第一作者简介

HUANG Yu, postgraduate, majors in earthquake geodesy, E-mail: 2315839404@qq.com.

### 文章历史

1. 中国地震局地震研究所，武汉市洪山侧路40号，430071;
2. 湖北省地震局，武汉市洪山侧路48号，430071

1 光学同震形变处理技术

1.1 亚像素相关性匹配技术

 $f_{\text {post }}(x, y)=f_{\text {pre }}(x-\Delta x, y-\Delta y)$ (1)

 $\begin{gather*} F_{\text {post }}(u, v)= \\ F_{\text {pre }}(u, v) * \exp (-\mathrm{j} 2 {\mathsf{π}}(u \Delta x+v \Delta y)) \end{gather*}$ (2)

 $\begin{gather*} \boldsymbol{H}(u, v)=\frac{\boldsymbol{F}_{\text {post }}(u, v) \boldsymbol{F}_{\text {pre }}(u, v)^{*}}{\left|\boldsymbol{F}_{\text {post }}(u, v) \boldsymbol{F}_{\text {pre }}(u, v)^{*}\right|}= \\ \exp (-\mathrm{j} 2 {\mathsf{π}}(u \Delta x+v \Delta y)) \end{gather*}$ (3)

 $\begin{gather*} \boldsymbol{H}(u, v)=\exp (-\mathrm{j} 2 {\mathsf{π}} u \Delta x) \exp (-\mathrm{j} 2 {\mathsf{π}} u \Delta y)= \\ \boldsymbol{q}_{\Delta x}(u) \boldsymbol{q}_{\Delta y}^{\mathrm{G}}(v) \end{gather*}$ (4)

 $\boldsymbol{H} \approx \boldsymbol{u}_{\max } * \sigma_{\max } * \boldsymbol{v}_{\max }^{\mathrm{G}}$ (5)

 $\Delta x=\frac{-k_{x}}{2 {\mathsf{π}}}, \Delta y=\frac{-k_{y}}{2 {\mathsf{π}}}$ (6)

1.2 光学同震形变提取的优化流程

1) 失相关区域掩膜。对同震形变场中信噪比(signal-noise ratio，SNR)低于设定经验阈值的像素点进行掩膜，抑制失相关噪声点对后续误差处理的影响。

2) 条带误差处理。为消除CCD线阵列姿态变化引起的条带误差，将每列条带区域内像元形变值叠加求平均，再由每个像元形变值减去平均值，去除条带误差。

3) 非局部均值滤波(NL-means)。相对于传统的中值滤波，实验选择非局部均值滤波(NL-means)平滑噪声、减少突跳点。NL-means滤波是基于图像自相似性的邻域滤波方法，充分利用图像中的冗余信息，在去噪的同时最大程度地保持图像的细节特征。在条带误差处理后，进行NL-means滤波，可进一步平滑噪声、突出同震形变场细节特征。其核心是计算图像中所有像素与当前像素间的相似性，一般会设定两个固定大小的窗口——一个大的搜索窗口和一个小的邻域窗口，邻域窗口在搜索窗口中进行滑动，根据邻域间的相似性修改权值来确定对应中心像素对当前像素的影响度。

 $\tilde{m}(a)=\sum\limits_{b \in I} w(a, b) * m(b)$ (7)

 \begin{aligned} w(a, b)=\frac{1}{\sum\limits_{b} \exp \left(-\frac{\|M(a)-M(b)\|^{2}}{h^{2}}\right)} \\ \;\;\;\;\;\;\;\; \exp \left(-\frac{\|M(a)-M(b)\|^{2}}{h^{2}}\right) \end{aligned} (8)
 $\begin{array}{*{20}{c}} \|M(a)-M(b)\|^{2}= \\ \frac{1}{d^{2}} \sum\limits_{c_{a} \in M(a), c_{b} \in M(b)}\left\|m\left(c_{a}\right)-m\left(c_{b}\right)\right\|^{2} \end{array}$ (9)

4) 轨道误差处理。为剔除图像中的多项式趋势，在非形变区域均匀选取若干个像素点，通过最小二乘平差求得4个参数和多项式轨道误差，即将一个多项式曲面拟合到图像上，并将其从图像中移除[5]

5) 赋空值-中值滤波。当形变场的空值较多时，会在一定程度上影响形变表达和解读，传统方法中对此不加考虑。实验通过中值滤波法将失相关掩膜空值区域赋值，只对失相关区域进行中值滤波，可以在不失真的情况下尽可能减小噪声影响。

 图 1 光学同震形变提取优化流程 Fig. 1 Optimization process for optical coseismic deformation extraction
1.3 主流光学同震形变软件对比分析

COSI-Corr是一款依赖ENVI平台的非开源免费软件。COSI-Corr软件采用频域相关器完成像素偏移检测。首先对固定尺寸的初始窗口内逐像素偏移进行粗略估计，并以此设定截止窗口进行亚像素相关性检测。当影像噪声较大或预期形变较大时，初始窗口尺寸设置较大为宜，而截止窗口尺寸应设计成能容纳一定噪声的最小尺寸。COSI-Corr软件根据对数-交叉频谱振幅大小对噪声频率进行屏蔽，筛选噪声值，修改每次重新计算的频率屏蔽次数，减少测量中的噪声值。

MicMac是完全开源的数字摄影测量和三维重建软件，采用正则化相关器，通过检测窗口尺寸和移动步长来调整检测精度和速度，从而完成亚像素相关性匹配[10]；采用非线性相关性系数(式(10))在像素对偏移检测的同时限制噪声对形变检测结果的影响：

 $\text { Cost }=C_{1} *\left(1-C_{3}\right)$ (10)
 $C_{1}=1-C^{\min }$ (11)
 $C_{2}=\frac{\operatorname{Max}\left(\text { Cor }, C^{\min }\right)-C^{\min }}{C_{1}}$ (12)
 $C_{3}=C_{2}^{\gamma}$ (13)

2 基于光学同震形变场的定量化形变算子

 $\begin{array}{*{20}{c}} {\partial u(x, y) / \partial x=}\\ {\frac{2 u(x+1, y)+u(x+1, y-1)+u(x+1, y+1)-(2 u(x-1, y)+u(x-1, y-1)+u(x-1, y+1))}{8\;\mathrm{px}}} \end{array}$ (14)
 $\begin{gather*} \partial u(x, y) / \partial y= \\ \frac{2 u(x, y+1)+u(x-1, y+1)+u(x+1, y+1)-(2 u(x, y-1)+u(x-1, y-1)+u(x+1, y-1))}{8\;\mathrm{px}} \end{gather*}$ (15)

 $\omega=(\nabla \times u)_{z}=\frac{\partial u_{\mathrm{EW}}}{\partial y}-\frac{\partial u_{\mathrm{NS}}}{\partial x}$ (16)
 $\delta=\nabla_{h} \cdot u=\frac{\partial u_{\mathrm{EW}}}{\partial x}+\frac{\partial u_{\mathrm{NS}}}{\partial y}$ (17)
 $\gamma=(\Delta \times u)_{z}=\frac{\partial u_{\mathrm{EW}}}{\partial y}+\frac{\partial u_{\mathrm{NS}}}{\partial x}$ (18)

3 2021年玛多MW7.4地震实例分析

2021-05-22青海玛多发生MW7.4地震，是继2008年汶川地震之后中国大陆发生的震级最高的一次地震，造成沿昆仑山口-江错断层158 km长的地表破裂[12]，并在东、西尾端发生分叉破裂[13]。玛多地震发生后，我们选用完全覆盖震区的10 m分辨率的Sentinel-2影像提取同震形变场。考虑到气候、时间、拍摄角度等影响，选择地震前、后同一时间段(震前：2020-10-12，震后：2021-10-17)、太阳方位角相近(159°~ 163°)、含云量低的各3景数据。考虑到短波红外波段(Band8)形变值标准差最小，选择Band8作为输入进行亚像素相关性匹配处理。震区水系分布较广，为最大程度抑制水体噪声，对其进行掩膜和空值赋值处理。其中将COSI-Corr中的初始窗口和截止窗口分别设置32×32 px和16×16 px、迭代次数为2、掩膜阈值为0.9；MicMac中的滑动窗口设置为3×3 px、相关性最小阈值为0.5、相关影响程度为2、正则化系数为0.3。

 图 2 EW方向同震形变场 Fig. 2 EW coseismic displacement field

 图 3 玛多地震精细化后处理 Fig. 3 Refined post-processing of the Maduo earthquake

 图 4 旋度、膨胀量、剪切应变 Fig. 4 Curl, dilatation, and shear strain
4 结语

 [1] Ferrario M F, Livio F. Conditional Probability of Distributed Surface Rupturing during Normal-Faulting Earthquakes[J]. Solid Earth, 2021, 12(5): 1 197-1 209 DOI:10.5194/se-12-1197-2021 (0) [2] Xu X H, Tong X P, Sandwell D T, et al. Refining the Shallow Slip Deficit[J]. Geophysical Journal International, 2016, 204(3): 843-862 (0) [3] Massonnet D, Feigl K L. Radar Interferometry and Its Application to Changes in the Earth's Surface[J]. Reviews of Geophysics, 1998, 36(4): 441-500 DOI:10.1029/97RG03139 (0) [4] Van Puymbroeck N, Michel R, Binet R, et al. Measuring Earthquakes from Optical Satellite Images[J]. Applied Optics, 2000, 39(20): 3 486 DOI:10.1364/AO.39.003486 (0) [5] 贺礼家, 冯光财, 冯志雄, 等. 哨兵-2号光学影像地表形变监测: 以2016年MW7.8新西兰凯库拉地震为例[J]. 测绘学报, 2019, 48(3): 339-351 (He Lijia, Feng Guangcai, Feng Zhixiong, et al. Coseismic Displacements of 2016 MW7.8 Kaikoura, New Zealand Earthquake, Using Sentinel-2 Optical Images[J]. Acta Geodaetica et Cartographica Sinica, 2019, 48(3): 339-351) (0) [6] Leprince S, Ayoub F, Klinger Y, et al. Co-Registration of Optically Sensed Images and Correlation(COSI-Corr): An Operational Methodology for Ground Deformation Measurements[C]. 2007 IEEE International Geoscience and Remote Sensing Symposium, Barcelona, 2008 (0) [7] Foroosh H, Zerubia J B, Berthod M. Extension of Phase Correlation to Subpixel Registration[J]. IEEE Transactions on Image Processing, 2002, 11(3): 188-200 DOI:10.1109/83.988953 (0) [8] Hermas E, Leprince S, Abou El-Magd I. Retrieving Sand Dune Movements Using Sub-Pixel Correlation of Multi-Temporal Optical Remote Sensing Imagery, Northwest Sinai Peninsula, Egypt[J]. Remote Sensing of Environment, 2012, 121: 51-60 DOI:10.1016/j.rse.2012.01.002 (0) [9] Leprince S, Barbot S, Ayoub F, et al. Automatic and Precise Orthorectification, Coregistration, and Subpixel Correlation of Satellite Images, Application to Ground Deformation Measurements[J]. IEEE Transactions on Geoscience and Remote Sensing, 2007, 45(6): 1 529-1 558 DOI:10.1109/TGRS.2006.888937 (0) [10] Rosu A M, Pierrot-Deseilligny M, Delorme A, et al. Measurement of Ground Displacement from Optical Satellite Image Correlation Using the Free Open-Source Software MicMac[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2015, 100: 48-59 DOI:10.1016/j.isprsjprs.2014.03.002 (0) [11] Milliner C W D, Dolan J F, Hollingsworth J, et al. Quantifying Near-Field and Off-Fault Deformation Patterns of the 1992 MW7.3 Landers Earthquake[J]. Geochemistry, Geophysics, Geosystems, 2015, 16(5): 1 577-1 598 DOI:10.1002/2014GC005693 (0) [12] 姚文倩, 王子君, 刘静, 等. 2021年青海玛多MW7.4地震同震地表破裂长度的讨论[J]. 地震地质, 2022, 44(2): 541-559 (Yao Wenqian, Wang Zijun, Liu Jing, et al. Discussion on Coseismic Surface Rupture Length of the 2021 MW7.4 Madoi Earthquake, Qinghai, China[J]. Seismology and Geology, 2022, 44(2): 541-559) (0) [13] 刘小利, 夏涛, 刘静, 等. 2021年青海玛多MW7.4地震分布式同震地表裂缝特征[J]. 地震地质, 2022, 44(2): 461-483 (Liu Xiaoli, Xia Tao, Liu Jing, et al. Distributed Characteristics of the Surface Deformations Associated with the 2021 MW7.4 Madoi Earthquake, Qinghai, China[J]. Seismology and Geology, 2022, 44(2): 461-483) (0)
Coseismic Deformation Extraction Method and Geodynamic Significances of Displacement Gradients Based on Optical Imagery
HUANG Yu1     LIU Xiaoli1,2     MO Xinyu1     DENG Debei'er1     RUAN Qiaozhe1     JIA Zhige1,2
1. Institute of Seismology, CEA, 40 Hongshance Road, Wuhan 430071, China;
2. Hubei Earthquake Agency, 48 Hongshance Road, Wuhan 430071, China
Abstract: This paper introduces the key technology, refined process, and differences in mainstream software packages for optical coseismic deformation extraction. Geoscientific operators such as displacement gradient, curl, dilatation, and shear strain are derived based on the optical coseismic deformation field. Taking the MW7.4 earthquake in Maduo, Qinghai province in 2021 as an example, we test the refined process of coseismic deformation extraction based on Sentinel-2 optical images. The results indicate post-processing of subpixel correlation matching results can effectively eliminate white noise and additive noise caused by orbit error, strip error, scenes and sensor factors, highlight deformation details, and more objectively describe the coseismic deformation spatial distribution.
Key words: coseismic deformation; optical imagery; displacement gradient; Sobel operator; 2021 Maduo earthquake