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

本发明属于基于雷达遥感技术的大地测量和地球物理领域，尤其是一种基于多源合成孔径雷达(sar)数据附加对数约束的同震震后时空滑动分布联合反演方法。

背景技术：

星载合成孔径雷达干涉测量(interferometricsyntheticapertureradar，insar)自上世纪90年代发展至今，已经被成功地应用在震间、同震和震后整个地震周期的形变监测中。特别是当地面数据(如gps,水准等)无法获取时，insar技术为地震研究提供了独一无二的形变数据源。雷达卫星得到的sar(syntheticapertureradar)数据可以通过insar技术获取地震周期形变，进而反演得到有限断层滑动模型。这不仅有助于提高我们对断层几何、断层滑动以及断层摩擦性质的认识，而且有利于评估区域地震风险。然而，由于sar卫星有限的重访周期以及数据获取政策，insar技术得到的同震滑动通常包含几天或几个月的震后形变贡献，相反得到的震后余滑缺失了这部分震后形变的影响，这必然增加了地震周期断层滑动估计的不确定性。如何有效分离基于insar技术反演得到的断层同震滑动和震后余滑对推动地震学的研究具有重要的科学指导意义。

然而实际研究应用中，由于sar卫星重访周期限制造成的基于insar技术反演得到的同震滑动分布常常受到震后余滑的严重污染，这严重阻碍了人类关于地震孕震机制认识和理解，并且严重影响了地震风险评估。为了更好地认识地球内部断层构造特性，需要有效地分离同震断层滑动和震后余滑。然而根据上述分析，目前相关同震震后滑动分离的研究具有一定的局限性，基于sar数据的断层滑动反演技术难以有效地分离同震断层滑动和震后余滑。

技术实现要素：

本发明的目的在于，克服现有基于sar数据的断层滑动反演技术中的不足和局限性，提供一种基于多源sar数据附加对数约束的同震震后时空滑动分布联合反演方法。

一种基于多源sar数据附加对数约束的同震震后时空滑动分布联合反演方法，包括以下步骤：

步骤一：利用insar技术处理多源sar数据得到地震区域的同震和震后雷达视线向(即斜距向)观测值，并通过地理编码将insar测量值从雷达坐标系下转换到通用的横轴墨卡托(universaltransversemercator,utm)投影坐标系；

步骤二：四叉树降采样insar观测值，利用多峰值粒子群算法求解发震断层倾角、发震断层走向、均一滑动量等主要断层几何参数，然后固定断层倾角和走向并将其细分为若干有限元断层；

步骤三：基于弹性位错理论构建地表形变和断层时空滑动分布之间的格林函数模型：

其中x为断层时空滑动向量，d为insar技术得到的地表形变，g为通过okada弹性位错模型得到的格林函数矩阵，b为设计系数矩阵，ε为模型残差；t表示获取的总计震后sar数据数量，和分别是断层有限元i在震后时刻j时的同震滑动向量和时空余滑累积向量，其中p为断层有限元个数，和分别表示不同时间节点的同震形变场和震后形变场。

震后余滑遵循对数式衰减的特征，由下式表征：

xpost＝a·log(1+(t-t0)/τ)

其中，其中，a为描述震后余滑大小的幅度系数，t为震后sar数据获取的时间节点，t0为地震发生时的时间节点，τ为描述震后余滑衰减速度的时间常数。

步骤四：联合上述两公式得到附加对数约束的同震震后时空滑动分布联合模型。利用附加约束的非线性最小二乘算法求解联合反演模型，得到所有断层有限元上的同震滑动分布以及每个断层有限元上的余滑幅度系数和衰减时间常数，进而得到研究时间段内整个发震断层上的时空滑动分布。

与现有技术相比，本发明的有益效果如下：

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测量值从雷达坐标系下转换到通用的横轴墨卡托(universaltransversemercator,utm)投影坐标系；

然后，四叉树降采样insar观测值，利用多峰值粒子群算法求解发震断层倾角、发震断层走向、均一滑动量等主要断层几何参数，然后固定断层倾角和走向并将其细分为若干有限元断层；

根据弹性位错理论，断层滑动x就会造成地表形变d，并且二者的关系满足：

d＝g·b·x+ε

其中g为通过okada弹性位错模型得到的格林函数矩阵，b为设计系数矩阵，ε为模型残差；t表示获取的总计震后sar数据数量，和分别是断层有限元i在震后时刻j时的同震滑动向量和时空余滑累积向量，其中p为断层有限元个数，和分别表示不同时间节点的同震形变场和震后形变场。

而对于震后余滑而言，根据其遵循对数函数式衰减的本质，可以将任意断层有限元上任意时刻的余滑写成：

xpost＝a·log(1+(t-t0)/τ)(1)

联合公式(1)和(2)可得如下矩阵形式联合模型：

其中，b1＝[1…1]1×p和b0＝[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数据和对数函数约束可以很好地分离同震滑动和震后余滑，因此说明本发明是可行可靠的。

对于本领域技术人员而言，显然本发明不限于上述示范性实施例的细节，而且在不背离本发明的精神或基本特征的情况下，能够以其他的具体形式实现本发明。因此，无论从哪一点来看，均应将实施例看作是示范性的，而且是非限制性的，本发明的范围由所附权利要求而不是上述说明限定，因此旨在将落在权利要求的等同要件的含义和范围内的所有变化囊括在本发明内。不应将权利要求中的任何附图标记视为限制所涉及的权利要求。