WO2021012592A1 - 一种基于多源sar数据附加对数约束的同震震后时空滑动分布联合反演方法 - Google Patents

一种基于多源sar数据附加对数约束的同震震后时空滑动分布联合反演方法 Download PDF

Info

Publication number
WO2021012592A1
WO2021012592A1 PCT/CN2019/126200 CN2019126200W WO2021012592A1 WO 2021012592 A1 WO2021012592 A1 WO 2021012592A1 CN 2019126200 W CN2019126200 W CN 2019126200W WO 2021012592 A1 WO2021012592 A1 WO 2021012592A1
Authority
WO
WIPO (PCT)
Prior art keywords
slip
fault
earthquake
coseismic
post
Prior art date
Application number
PCT/CN2019/126200
Other languages
English (en)
French (fr)
Inventor
许文斌
刘小鸽
Original Assignee
中南大学
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 中南大学 filed Critical 中南大学
Publication of WO2021012592A1 publication Critical patent/WO2021012592A1/zh

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/885—Radar or analogous systems specially adapted for specific applications for ground probing
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9021—SAR image post-processing techniques
    • G01S13/9023—SAR image post-processing techniques combined with interferometric techniques
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9021—SAR image post-processing techniques
    • G01S13/9027—Pattern recognition for feature extraction
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • G01S7/41—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section

Definitions

  • the invention belongs to the field of geodesy and geophysics based on radar remote sensing technology, in particular to a joint inversion method of coseismic post-coseismic sliding distribution based on multi-source synthetic aperture radar (SAR) data with logarithmic constraints.
  • SAR synthetic aperture radar
  • Spaceborne Synthetic Aperture Radar Interferometry Interferometric Synthetic Aperture Radar, InSAR
  • InSAR Interferometric Synthetic Aperture Radar
  • ground data such as GPS, leveling, etc.
  • SAR Synthetic Aperture Radar
  • SAR data obtained by radar satellites can be used to obtain seismic periodic deformation through InSAR technology, and then inversion to obtain a finite fault slip model. This not only helps to improve our understanding of fault geometry, fault slip, and fault friction properties, but also helps assess regional seismic risk.
  • the coseismic slip obtained by InSAR technology usually contains a few days or months of post-earthquake deformation contribution.
  • the post-earthquake afterslip obtained lacks this part of the post-earthquake deformation. This will inevitably increase the uncertainty of fault slip estimation during the earthquake period.
  • How to effectively separate the fault co-seismic slip and post-earthquake afterslip obtained by the InSAR technology inversion has important scientific guiding significance for the advancement of seismological research.
  • the purpose of the present invention is to overcome the shortcomings and limitations of the existing fault slip inversion technology based on SAR data, and provide a joint inversion method of coseismic post-coseismic slip distribution based on multi-source SAR data with additional logarithmic constraints.
  • a joint inversion method for post-coseismic spatio-temporal slip distribution based on multi-source SAR data with additional logarithmic constraints including the following steps:
  • Step 1 Use InSAR technology to process multi-source SAR data to obtain co-seismic and post-seismic radar line-of-sight observations (ie oblique range) in the seismic area, and use geocoding to convert the InSAR measurements from the radar coordinate system to the universal horizontal Axis Mercator (Universal Transverse Mercator, UTM) projection coordinate system;
  • UTM Universal Transverse Mercator
  • Step 2 Quadtree downsampling the InSAR observations, using the multi-peak particle swarm algorithm to solve the main fault geometric parameters such as the inclination of the seismogenic fault, the strike of the seismogenic fault, and the uniform slip, and then fix the inclination and strike of the fault and subdivide them into Several finite element faults;
  • Step 3 Based on the elastic dislocation theory, construct a Green's function model between the surface deformation and the temporal and spatial slip distribution of the fault:
  • X is the time-space sliding vector of the fault
  • D is the surface deformation obtained by InSAR technology
  • G is the Green's function matrix obtained by Okada elastic dislocation model
  • B is the design coefficient matrix
  • ⁇ is the model residual
  • T is the total post-earthquake obtained SAR data quantity, with They are the co-seismic slip vector and time-space residual slip accumulation vector of the fault finite element i at the post-earthquake time j, where P is the number of fault finite elements, with Respectively represent the coseismic deformation field and post-earthquake deformation field at different time nodes.
  • the after-slip after the earthquake follows the characteristics of logarithmic attenuation, which is characterized by the following formula:
  • A is the amplitude coefficient describing the magnitude of the afterslip after the earthquake
  • t is the time node of the SAR data acquisition after the earthquake
  • t 0 is the time node when the earthquake occurs
  • ⁇ is the time constant describing the attenuation rate of the afterslip after the earthquake.
  • Step 4 Combine the above two formulas to obtain a combined model of post-coseismic slip distribution with additional logarithmic constraints.
  • the joint inversion model is solved by the nonlinear least squares algorithm with additional constraints, and the coseismic slip distribution on all fault finite elements and the after-slip amplitude coefficient and attenuation time constant on each fault finite element are obtained, and then the research time period is obtained.
  • the entire invention process is clear, simple to implement, not limited by region, and does not rely on any other geodetic data (such as GPS) and seismic wave data, which provides a low level of effective separation of coseismic sliding distribution and post-earthquake afterslip. Cost, high time resolution and effective and feasible method.
  • geodetic data such as GPS
  • seismic wave data which provides a low level of effective separation of coseismic sliding distribution and post-earthquake afterslip. Cost, high time resolution and effective and feasible method.
  • the present invention breaks through the bottleneck that it is difficult to strictly separate co-seismic slip and post-earthquake afterslip when inverting fault slip based on InSAR observations, and expands the understanding and understanding of the mechanism of the earthquake rupture process and the distribution of after-earthquake afterslip.
  • the internal structure research has important scientific value and guiding significance; and the present invention has the compatibility of multi-source SAR satellite data, and actively promotes the practical and market-oriented development of SAR data and InSAR technology.
  • the present invention mainly uses the after-earthquake afterslip to follow the nature of logarithmic function attenuation, and effectively inverts and separates co-seismic sliding and post-earthquake afterslip based on InSAR technology, which greatly reduces the amount of model waiting. At the same time, it makes full use of multi-source SAR data to make up for the limitation of the revisit period of single-orbit SAR satellites to improve the time resolution of time-space slip distribution.
  • Figure 1 is a flow chart of a joint inversion method for post-coseismic spatio-temporal slip distribution based on multi-source SAR data with logarithmic constraints.
  • Figure 2 is a schematic diagram of simulated SAR data for different satellite orbits.
  • Figure 3 shows the distribution of time-space slip after coseismic (the first row) and inversion (the second row) and its relative residuals (the third row).
  • the first column represents the coseismic sliding results; the other columns represent the post-seismic sliding results.
  • the black contours with an interval of 2.5m in the figure indicate the coseismic sliding distribution of simulation (first row) and inversion (second row).
  • Figure 4 shows the simulated ascending orbit co-seismic (first row), after ascending orbit earthquake (second) and after descending orbital earthquake (third row) deformation field.
  • ac represents the original simulated deformation field
  • df represents the deformation field after corresponding downsampling
  • gi represents the deformation field after downsampling and adding noise
  • jl represents the deformation field of the optimal model forward
  • mo represents the original simulated deformation field and the optimal model Forward the residuals between the deformation fields
  • pr represents the corresponding residual histogram.
  • Figure 5 is a correlation diagram between simulated slip distribution and estimated fault slip.
  • the quadtree downsamples the InSAR observations, uses the multi-peak particle swarm algorithm to solve the main fault geometric parameters such as the seismogenic fault dip, the seismogenic fault trend, and the uniform slip, and then fixes the fault dip and trend and subdivides them into several Finite element fault;
  • fault sliding X will cause surface deformation D, and the relationship between the two satisfies:
  • G is the Green's function matrix obtained by Okada elastic dislocation model
  • B is the design coefficient matrix
  • ⁇ is the model residual
  • T is the total number of post-earthquake SAR data obtained, with They are the co-seismic slip vector and time-space residual slip accumulation vector of the fault finite element i at the post-earthquake time j, where P is the number of fault finite elements, with Respectively represent the coseismic deformation field and post-earthquake deformation field at different time nodes.
  • the afterslip at any time on any fault finite element can be written as:
  • ⁇ represents the L2 norm criterion, combined with the optimal logarithmic function parameters to restore and reconstruct the time-space residual slip, and finally get the time-space slip distribution model in the study period.
  • the ascending orbit includes the data before the earthquake of 1 scene and the data after the earthquake of 5 scenes, and the orbiting data only contains the data after the earthquake of 5 scenes.
  • the lifting rails can form 15 and 10 InSAR interference pairs respectively. Assuming that the strike angle of the fault is 350°, the dip angle is 15°, the coseismic slip angle and post-seismic slip angle are 150° and 120°, respectively, and the after-slip decay time constant is 10 days. Then divide the fault into 120 rectangular finite elements with sides of 5 kilometers. Then, according to the simulated SAR data, the time-space slip distribution is simulated on the fault plane at all times (the first row in Figure 3).
  • the deformation field of the InSAR observation value is forwarded ( Figure 4a-c).
  • Figure 4d-f the quadtree downsampling method is used to simulate the observation value for downsampling processing.
  • Gaussian white noise with a mean value of zero and a standard deviation of 10% of the maximum simulated deformation is added to the simulated observations ( Figure 4g-i). It should be noted that Fig. 4 only shows a representative deformation map after an ascending orbit coseism, after an ascending orbit earthquake, and after a descending orbital earthquake.
  • the above-mentioned simulated noisy InSAR observation value inversion can be used to separate co-seismic slip and post-seismic spatio-temporal residual slip.
  • the second row of Figure 3 shows the spatiotemporal fault slip distribution obtained by the inversion of the present invention. It can be seen that the model inversion results and the simulated fault slip have good consistency in the temporal and spatial distribution.
  • the residual root mean square error (RMS) of the coseismic sliding distribution is calculated based on the residual of the sliding distribution (the third row of Figure 3) is only 8.5% of the coseismic sliding, and the residual RMS of the residual The maximum is 7.1cm, which is only 3.5% of the simulated after slip. Therefore, the model can well recover the slip distribution of the separated coseismic faults and the time-space residual slip, thus indicating that the present invention is feasible and reliable.
  • Fig. 4j-1 shows the deformation field obtained by the inversion sliding forward modeling. It can be seen that its magnitude and spatial distribution are very consistent with the original simulated deformation field.
  • Figure 4m-o shows the residuals between the forward deformation results and the simulated values. It can be seen that the relative deformation values of the residuals of the three interference pairs are very small, and the main residuals appear on the northeast of the fault. This is mainly due to the design of the fault. It is caused by tilting to the northeast and adding noise to the down-sampling deformation field, and the residual distribution shown in Figure 4p-r is similar to the Gaussian distribution and the RMS is less than the added noise level, so the residuals are mainly derived from the added noise.
  • Figure 5 is a graph of the correlation coefficient between the simulated slip distribution and the inverted slip distribution.
  • the correlation coefficient can reach 0.92, and the co-seismic slip angle and residual
  • the sliding angle and the decay time constant are 150.1° ⁇ 0.4°, 119.4° ⁇ 0.8°, and 10.6 ⁇ 0.9 days, respectively, which are slightly different from the model input values, indicating that the joint model has good performance.
  • the present invention can well separate co-seismic slip and post-earthquake after-slip based on multi-source SAR data and logarithmic function constraints, thus indicating that the present invention is feasible and reliable.

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

一种基于对数约束和多源SAR数据的同震震后时空滑动分布联合反演方法：首先获取地震区域多源SAR数据并利用InSAR技术获取该区域的同震和震后时空形变场；随后基于震后余滑遵循对数函数式衰减的本质和弹性位错理论构建InSAR多源观测值与断层同震震后时空滑动之间的非线性函数模型；然后利用非线性求解方法反演出有限元断层上的同震滑动分布，震后余滑衰减常数以及相应的滑动幅度系数；进而利用非线性模型参数求解震后滑动时空分布，同时联合多源SAR数据提高断层时空滑动的时间分辨率；最后根据反演的断层时空滑动分布计算地震释放能量、应力状态以及构造摩擦性质等重要地球物理参数。

Description

一种基于多源SAR数据附加对数约束的同震震后时空滑动分布联合反演方法 技术领域
本发明属于基于雷达遥感技术的大地测量和地球物理领域，尤其是一种基于多源合成孔径雷达(SAR)数据附加对数约束的同震震后时空滑动分布联合反演方法。
背景技术
星载合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar，InSAR)自上世纪90年代发展至今，已经被成功地应用在震间、同震和震后整个地震周期的形变监测中。特别是当地面数据(如GPS,水准等)无法获取时，InSAR技术为地震研究提供了独一无二的形变数据源。雷达卫星得到的SAR(Synthetic Aperture Radar)数据可以通过InSAR技术获取地震周期形变，进而反演得到有限断层滑动模型。这不仅有助于提高我们对断层几何、断层滑动以及断层摩擦性质的认识，而且有利于评估区域地震风险。然而，由于SAR卫星有限的重访周期以及数据获取政策，InSAR技术得到的同震滑动通常包含几天或几个月的震后形变贡献，相反得到的震后余滑缺失了这部分震后形变的影响，这必然增加了地震周期断层滑动估计的不确定性。如何有效分离基于InSAR技术反演得到的断层同震滑动和震后余滑对推动地震学的研究具有重要的科学指导意义。
然而实际研究应用中，由于SAR卫星重访周期限制造成的基于 InSAR技术反演得到的同震滑动分布常常受到震后余滑的严重污染，这严重阻碍了人类关于地震孕震机制认识和理解，并且严重影响了地震风险评估。为了更好地认识地球内部断层构造特性，需要有效地分离同震断层滑动和震后余滑。然而根据上述分析，目前相关同震震后滑动分离的研究具有一定的局限性，基于SAR数据的断层滑动反演技术难以有效地分离同震断层滑动和震后余滑。
发明内容
本发明的目的在于，克服现有基于SAR数据的断层滑动反演技术中的不足和局限性，提供一种基于多源SAR数据附加对数约束的同震震后时空滑动分布联合反演方法。
一种基于多源SAR数据附加对数约束的同震震后时空滑动分布联合反演方法，包括以下步骤：
步骤一：利用InSAR技术处理多源SAR数据得到地震区域的同震和震后雷达视线向(即斜距向)观测值，并通过地理编码将InSAR测量值从雷达坐标系下转换到通用的横轴墨卡托(Universal Transverse Mercator,UTM)投影坐标系；
步骤二：四叉树降采样InSAR观测值，利用多峰值粒子群算法求解发震断层倾角、发震断层走向、均一滑动量等主要断层几何参数，然后固定断层倾角和走向并将其细分为若干有限元断层；
步骤三：基于弹性位错理论构建地表形变和断层时空滑动分布之间的格林函数模型：
D＝G·B·X+ε
Figure PCTCN2019126200-appb-000001
Figure PCTCN2019126200-appb-000002
其中X为断层时空滑动向量，D为InSAR技术得到的地表形变，G为通过Okada弹性位错模型得到的格林函数矩阵，B为设计系数矩阵，ε为模型残差；T表示获取的总计震后SAR数据数量，
Figure PCTCN2019126200-appb-000003
和
Figure PCTCN2019126200-appb-000004
分别是断层有限元i在震后时刻j时的同震滑动向量和时空余滑累积向量，其中P为断层有限元个数，
Figure PCTCN2019126200-appb-000005
和
Figure PCTCN2019126200-appb-000006
分别表示不同时间节点的同震形变场和震后形变场。
震后余滑遵循对数式衰减的特征，由下式表征：
x post＝A·log(1+(t-t 0)/τ)
其中，其中，A为描述震后余滑大小的幅度系数，t为震后SAR数据获取的时间节点，t 0为地震发生时的时间节点，τ为描述震后余滑衰减速度的时间常数。
步骤四：联合上述两公式得到附加对数约束的同震震后时空滑动分布联合模型。利用附加约束的非线性最小二乘算法求解联合反演模型，得到所有断层有限元上的同震滑动分布以及每个断层有限元上的余滑幅度系数和衰减时间常数，进而得到研究时间段内整个发震断层上的时空滑动分布。
与现有技术相比，本发明的有益效果如下：
1、整个发明流程清晰，实现简单，不受区域的限制，且不依赖于任何其它大地测量数据(如GPS)和地震波数据，对有效分离同震 滑动分布和震后余滑提供了一种低成本、高时间分辨率且有效可行的方法。
2、本发明突破了基于InSAR观测值反演断层滑动时难以严格分离同震滑动和震后余滑的瓶颈，拓展了人类对于地震破裂过程和震后余滑分布的机理认识和理解，对地球内部构造研究具有重要的科学价值和指导意义；而且本发明具有多源SAR卫星数据的兼容性，积极推动了SAR数据及InSAR技术的实用化和市场化发展。
3、本发明主要利用震后余滑遵循对数函数式衰减的本质，基于InSAR技术有效的反演分离同震滑动和震后余滑，大大减少了模型待求量。与此同时，充分利用多源SAR数据弥补单轨SAR卫星重访周期的限制来提高时空滑动分布的时间分辨率。
附图说明
图1是一种基于多源SAR数据附加对数约束的同震震后时空滑动分布联合反演方法的流程图。
图2是不同卫星轨道的模拟SAR数据示意图。
图3是模拟(第一行)和反演(第二行)同震震后时空滑动分布图以及其相对残差(第三行)，第一列表示同震滑动结果；其他列表示震后滑动结果，图中间隔2.5m的黑色等值线表示模拟(第一行)和反演(第二行)的同震滑动分布。
图4是模拟升轨同震(第一行)、升轨震后(第二)和降轨震后(第三行)形变场。a-c表示原始模拟形变场；d-f表示相应降采样后 的形变场；g-i表示降采样并加噪后的形变场；j-l表示最优模型正演的形变场；m-o表示原始模拟形变场与最优模型正演形变场之间的残差；p-r表示相应的残差直方图。
图5是模拟滑动分布和估计断层滑动之间的相关图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
在为了便于理解本发明，首先提供本发明的理论基础：
首先，利用InSAR技术处理多源SAR数据得到地震区域的同震和震后雷达视线向(即斜距向)观测值，并通过地理编码将InSAR测量值从雷达坐标系下转换到通用的横轴墨卡托(Universal Transverse Mercator,UTM)投影坐标系；
然后，四叉树降采样InSAR观测值，利用多峰值粒子群算法求解发震断层倾角、发震断层走向、均一滑动量等主要断层几何参数，然后固定断层倾角和走向并将其细分为若干有限元断层；
根据弹性位错理论，断层滑动X就会造成地表形变D，并且二者的关系满足：
D＝G·B·X+ε
Figure PCTCN2019126200-appb-000007
Figure PCTCN2019126200-appb-000008
Figure PCTCN2019126200-appb-000009
其中G为通过Okada弹性位错模型得到的格林函数矩阵，B为设计系数矩阵，ε为模型残差；T表示获取的总计震后SAR数据数量，
Figure PCTCN2019126200-appb-000010
和
Figure PCTCN2019126200-appb-000011
分别是断层有限元i在震后时刻j时的同震滑动向量和时空余滑累积向量，其中P为断层有限元个数，
Figure PCTCN2019126200-appb-000012
和
Figure PCTCN2019126200-appb-000013
分别表示不同时间节点的同震形变场和震后形变场。
而对于震后余滑而言，根据其遵循对数函数式衰减的本质，可以将任意断层有限元上任意时刻的余滑写成：
x post＝A·log(1+(t-t 0)/τ)        (1)
联合公式(1)和(2)可得如下矩阵形式联合模型：
Figure PCTCN2019126200-appb-000014
其中，B 1＝[1…1] 1×P和B 0＝[0…0] 1×P。
从公式(3)可以看出，模型参数个数仅与断层有限元个数有关，因此大大减少了待求参数量。公式(3)所示联合模型参数仅有2P+3，其中包括P个同震有限元断层滑动和P个震后滑动幅度系数，一个均 一断层滑动衰减时间常数τ以及嵌套在格林矩阵中的两个断层滑动角，该模型待求参数远远少于公式(1)中的模型待求参数。此外由于加入了对数函数约束，该联合模型将不存在矩阵质亏的问题，然而由于该系统是一个非线性系统，因此需要利用非线性方程求解方法来求解，本项目拟采用附加约束的非线性最小二乘方法来求解该联合反演模型得到模型最优解：
||D-G·B·X||＝min
其中‖·‖代表L2范准则，结合最优对数函数参数恢复重建时空余滑，最终得到研究时间段内时空滑动分布模型。
本发明效果可以通过以下模拟仿真实验作为实施例进一步说明。
首先模拟两个不同SAR卫星轨道的升降轨数据，如图2所示，升轨包含了1景地震前的数据和5景震后数据，降轨仅包含5景震后数据。升降轨分别可以形成15和10个InSAR干涉对。假设断层走向角350°、倾角15°，同震滑动角和震后滑动角分别为150°和120°，余滑衰减时间常数为10天。然后将断层分割为120个边长为5千米的矩形有限元。然后根据模拟SAR数据时刻在断层面上模拟时空滑动分布(图3第一行)。进而根据滑动分布正演InSAR观测值形变场(图4a-c)。考虑到InSAR得到的单个形变场就有数百万观测量，在降低数据量的情况下最大程度地保留形变细节，采用四叉树降采样方法来模拟观测值进行降采样处理(图4d-f)。为了让模拟观测值具有真实性，在模拟观测值中加入均值为零、标准偏差为模拟形变最大值10％的高斯白噪声(图4g-i)。需要说明的是图4仅给出了升轨同震、升轨震后和 降轨震后各一个代表性形变图。
通过本发明所提出的联合反演方法，就可以利用上述模拟的含噪InSAR观测值反演分离出同震滑动和震后时空余滑。图3第二行显示的就是本发明反演得到的时空断层滑动分布。可以看出，模型反演结果和模拟断层滑动在时空分布上具有良好的一致性。为了定量验证本发明的效果，根据滑动分布残差(图3第三行)计算得到同震滑动分布残差均方根误差(RMS)仅为同震滑动的8.5％，且余滑残差RMS最大为7.1cm，仅相当于模拟余滑量的3.5％。因此模型可以很好地恢复分离同震断层滑动分布和时空余滑，从而说明本发明是可行的和可靠的。为了进一步说明本发明的效果，图4j-l给出了反演滑动正演得到的形变场，可以看出其量级和空间分布均与原始模拟形变场非常一致。图4m-o显示的正演形变结果与模拟值之间的残差，可以看出，三个干涉对残差相对形变值非常小，主要残差出现在断层东北边，这主要是由于设计断层向东北方向倾斜和对降采样形变场加噪所致，且图4p-r显示的残差分布近似高斯分布且RMS小于所加噪声水平，因此残差主要来源于所加噪声。图5是模拟滑动分布和反演滑动分布之间的相关系数图，在噪声水平为最大形变值10％的情况下，相关系数可以达到0.92，且联合模型反演得到的同震滑动角、余滑滑动角和衰减时间常数分别为150.1°±0.4°、119.4°±0.8°和10.6±0.9天，这与模型输入值相差微乎其微，说明联合模型具有良好的性能。总体而言，本发明基于多源SAR数据和对数函数约束可以很好地分离同震滑动和震后余滑，因此说明本发明是可行可靠的。
对于本领域技术人员而言，显然本发明不限于上述示范性实施例的细节，而且在不背离本发明的精神或基本特征的情况下，能够以其他的具体形式实现本发明。因此，无论从哪一点来看，均应将实施例看作是示范性的，而且是非限制性的，本发明的范围由所附权利要求而不是上述说明限定，因此旨在将落在权利要求的等同要件的含义和范围内的所有变化囊括在本发明内。不应将权利要求中的任何附图标记视为限制所涉及的权利要求。

Claims (4)

  1. 一种基于多源SAR数据附加对数约束的同震震后时空滑动分布联合反演方法，其特征在于，包括以下步骤：
    步骤一：利用InSAR技术处理多源SAR数据得到地震区域的同震和震后雷达视线向(即斜距向)观测值，并通过地理编码将InSAR测量值从雷达坐标系下转换到通用的横轴墨卡托投影坐标系；
    步骤二：通过四叉树降采样InSAR观测值，利用多峰值粒子群算法求解包括发震断层倾角、发震断层走向、均一滑动量的断层几何参数，然后固定断层倾角和走向并将其细分为若干有限元断层；
    步骤三：基于弹性位错理论构建地表形变和断层时空滑动分布之间的格林函数模型，其表示如下
    Figure PCTCN2019126200-appb-100001
    其中X为断层时空滑动向量，D为InSAR技术得到的地表形变，G为通过弹性位错模型得到的格林函数矩阵，B为设计系数矩阵，ε为模型残差；T表示获取的总计震后SAR数据数量，
    Figure PCTCN2019126200-appb-100002
    和
    Figure PCTCN2019126200-appb-100003
    分别是断层有限元i在震后时刻j时的同震滑动向量和时空余滑累积向量，其中P为断层有限元个数，
    Figure PCTCN2019126200-appb-100004
    和
    Figure PCTCN2019126200-appb-100005
    分别表示不同时间节点的同震形变场和震后形变场；
    而震后余滑遵循对数式衰减的特征，由下式表征：
    x post＝A·log(1+(t-t 0)/τ)    (2)
    其中，A为描述震后余滑大小的幅度系数，t为震后SAR数据获取的时间节点，t 0为地震发生时的时间节点，τ为描述震后余滑衰减 速度的时间常数；
    联合公式(1)和(2)将该公式(1)模型中待求的震后时空余滑转化为对数函数模型，可得如下矩阵形式的附加对数约束的同震震后时空滑动分布联合反演模型：
    Figure PCTCN2019126200-appb-100006
    其中，B 1＝[1…1] 1×P和B 0＝[0…0] 1×P；
    步骤四：求解联合反演模型，得到所有断层有限元上的同震滑动分布以及每个断层有限元上的余滑幅度系数和衰减时间常数，进而得到研究时间段内整个发震断层上的时空滑动分布。
  2. 根据权利要求1所述的基于多源SAR数据附加对数约束的同震震后时空滑动分布联合反演方法，其特征在于，使用附加约束的非线性最小二乘求解所述联合反演模型。
  3. 根据权利要求1所述的基于多源SAR数据附加对数约束的同震震后时空滑动分布联合反演方法，其特征在于，公式(3)所示联 合反演模型参数仅有2P+3，其中包括P个同震有限元断层滑动和P个震后滑动幅度系数，一个均一断层滑动衰减时间常数τ以及嵌套在格林矩阵中的两个断层滑动角。
  4. 根据权利要求1-3所述的一种基于多源SAR数据附加对数约束的同震震后时空滑动分布联合反演方法，其特征在于，所述地震形变为任意区域任意地震类型造成的形变。
PCT/CN2019/126200 2019-07-19 2019-12-18 一种基于多源sar数据附加对数约束的同震震后时空滑动分布联合反演方法 WO2021012592A1 (zh)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN201910655738.4A CN110333508B (zh) 2019-07-19 2019-07-19 基于多源sar数据的同震震后时空滑动分布联合反演方法
CN201910655738.4 2019-07-19

Publications (1)

Publication Number Publication Date
WO2021012592A1 true WO2021012592A1 (zh) 2021-01-28

Family

ID=68145982

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2019/126200 WO2021012592A1 (zh) 2019-07-19 2019-12-18 一种基于多源sar数据附加对数约束的同震震后时空滑动分布联合反演方法

Country Status (2)

Country Link
CN (1) CN110333508B (zh)
WO (1) WO2021012592A1 (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117111057A (zh) * 2023-08-23 2023-11-24 首都师范大学 一种煤矿采空区形变敏感性评价方法

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110333508B (zh) * 2019-07-19 2021-02-19 中南大学 基于多源sar数据的同震震后时空滑动分布联合反演方法
CN112147610B (zh) * 2020-09-28 2022-08-05 武汉大学 一种融合拓扑分析技术的InSAR干涉图降采样方法
CN112395789A (zh) * 2020-10-23 2021-02-23 马培峰 一种耦合InSAR和数值模拟分析城区滑坡形变的方法
CN113253270A (zh) * 2021-06-11 2021-08-13 中国测绘科学研究院 基于InSAR和Okada模型反演地下采矿参数的方法和系统

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101770027A (zh) * 2010-02-05 2010-07-07 河海大学 基于InSAR与GPS数据融合的地表三维形变监测方法
US20130144574A1 (en) * 2011-12-06 2013-06-06 Korea Institute Of Geoscience And Mineral Resources Method of Analyzing 3D Geological Structure Using Structure Index
CN104123464A (zh) * 2014-07-23 2014-10-29 中国国土资源航空物探遥感中心 一种高分辨率InSAR时序分析反演地物高程与地面沉降量的方法
CN105158760A (zh) * 2015-08-10 2015-12-16 中南大学 一种利用InSAR反演地下流体体积变化和三维地表形变的方法
CN108446516A (zh) * 2018-04-04 2018-08-24 长安大学 一种基于力源反演的构造断裂与地面沉降形变分解法
CN109738892A (zh) * 2019-01-24 2019-05-10 中南大学 一种矿区地表高时空分辨率三维形变估计方法
CN110333508A (zh) * 2019-07-19 2019-10-15 中南大学 一种基于多源sar数据附加对数约束的同震震后时空滑动分布联合反演方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9207318B2 (en) * 2011-06-20 2015-12-08 California Institute Of Technology Damage proxy map from interferometric synthetic aperture radar coherence
US10145972B2 (en) * 2014-08-15 2018-12-04 California Institute Of Technology Systems and methods for advanced rapid imaging and analysis for earthquakes
CN109444879A (zh) * 2018-10-19 2019-03-08 西南交通大学 一种DInSAR近断层同震形变场提取方法

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101770027A (zh) * 2010-02-05 2010-07-07 河海大学 基于InSAR与GPS数据融合的地表三维形变监测方法
US20130144574A1 (en) * 2011-12-06 2013-06-06 Korea Institute Of Geoscience And Mineral Resources Method of Analyzing 3D Geological Structure Using Structure Index
CN104123464A (zh) * 2014-07-23 2014-10-29 中国国土资源航空物探遥感中心 一种高分辨率InSAR时序分析反演地物高程与地面沉降量的方法
CN105158760A (zh) * 2015-08-10 2015-12-16 中南大学 一种利用InSAR反演地下流体体积变化和三维地表形变的方法
CN108446516A (zh) * 2018-04-04 2018-08-24 长安大学 一种基于力源反演的构造断裂与地面沉降形变分解法
CN109738892A (zh) * 2019-01-24 2019-05-10 中南大学 一种矿区地表高时空分辨率三维形变估计方法
CN110333508A (zh) * 2019-07-19 2019-10-15 中南大学 一种基于多源sar数据附加对数约束的同震震后时空滑动分布联合反演方法

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117111057A (zh) * 2023-08-23 2023-11-24 首都师范大学 一种煤矿采空区形变敏感性评价方法

Also Published As

Publication number Publication date
CN110333508B (zh) 2021-02-19
CN110333508A (zh) 2019-10-15

Similar Documents

Publication Publication Date Title
WO2021012592A1 (zh) 一种基于多源sar数据附加对数约束的同震震后时空滑动分布联合反演方法
Feng et al. The 2011 MW 6.8 Burma earthquake: Fault constraints provided by multiple SAR techniques
Funning et al. Surface displacements and source parameters of the 2003 Bam (Iran) earthquake from Envisat advanced synthetic aperture radar imagery
Delouis et al. Joint inversion of InSAR, GPS, teleseismic, and strong-motion data for the spatial and temporal distribution of earthquake slip: Application to the 1999 Izmit mainshock
Wright et al. Constraining the slip distribution and fault geometry of the M w 7.9, 3 November 2002, Denali fault earthquake with interferometric synthetic aperture radar and global positioning system data
Simons et al. Coseismic deformation from the 1999 M w 7.1 Hector Mine, California, earthquake as inferred from InSAR and GPS observations
Feng et al. A slip gap of the 2016 M w 6.6 Muji, Xinjiang, China, earthquake inferred from Sentinel‐1 TOPS interferometry
Zhu et al. Interseismic slip rate and locking along the Maqin–Maqu Segment of the East Kunlun Fault, Northern Tibetan Plateau, based on Sentinel-1 images
Fielding et al. Surface deformation of north‐central Oklahoma related to the 2016 M w 5.8 Pawnee earthquake from SAR interferometry time series
Shen et al. A spatially varying scaling method for InSAR tropospheric corrections using a high‐resolution weather model
CN112556660B (zh) 基于卫星测高数据的海域重力异常反演方法及系统
Cai et al. A new algorithm for landslide dynamic monitoring with high temporal resolution by Kalman filter integration of multiplatform time-series InSAR processing
Kimata et al. Understanding the 2007–2008 eruption of Anak Krakatau Volcano by combining remote sensing technique and seismic data
Lyu et al. Overall subshear but locally supershear rupture of the 2021 Mw 7.4 Maduo earthquake from high-rate GNSS waveforms and three-dimensional InSAR deformation
Lauer et al. Fault geometry and slip distribution of the 2013 Mw 7.7 Balochistan earthquake from inversions of SAR and optical data
Eshagh On integral approach to regional gravity field modelling from satellite gradiometric data
Takada et al. Interseismic crustal deformation in and around the Atotsugawa fault system, central Japan, detected by InSAR and GNSS
De Natale et al. Abruzzo, Italy, earthquakes of April 2009: heterogeneous fault-slip models and stress transfer from accurate inversion of ENVISAT-InSAR data
Zhang et al. Research on deformation characteristics of the 2021 Qinghai Maduo MS7. 4 earthquake through coseismic dislocation inversion
Morishita et al. Complex crustal deformation of the 2016 Kaikoura, New Zealand, earthquake revealed by ALOS‐2
He et al. The 2022 M w 6.6 Menyuan Earthquake in the Northwest Margin of Tibet: Geodetic and Seismic Evidence of the Fault Structure and Slip Behavior of the Qilian–Haiyuan Strike‐Slip Fault
Li et al. The 1998 Mw 5.7 Zhangbei‐Shangyi (China) earthquake revisited: A buried thrust fault revealed with interferometric synthetic aperture radar
Bagherbandi et al. Improving gravimetric–isostatic models of crustal depth by correcting for non-isostatic effects and using CRUST2. 0
Zhou et al. A new GRACE filtering approach based on iterative image convolution
Du et al. InSAR-based rapid damage assessment of urban building portfolios following the 2023 Turkey earthquake

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 19938582

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 19938582

Country of ref document: EP

Kind code of ref document: A1

32PN Ep: public notification in the ep bulletin as address of the adressee cannot be established

Free format text: NOTING OF LOSS OF RIGHTS PURSUANT TO RULE 112(1) EPC (EPO FORM 1205A DATED 23/06/2022)

122 Ep: pct application non-entry in european phase

Ref document number: 19938582

Country of ref document: EP

Kind code of ref document: A1