一种基于单轨InSAR观测的三维地壳形变转换方法与流程

本发明属于大地测量与地壳形变监测技术领域，并具体涉及一种稳健的基于单轨insar观测的三维地壳形变转换方法。

背景技术：

合成孔径干涉雷达(interferometricsyntheticapertureradar,insar)技术以其空间分辨率高，覆盖范围广，对形变敏感度高等优势成为了高密度面状形变主要观测手段，已广泛应用于监测大尺度连续的地表形变，如地震形变、火山运动、地面沉降、矿山形变、火山活动、冰川漂移、山体滑坡等领域。

然而，insar探测得到的只是地表三维形变北(n)、东(e)、上(u)方向在雷达视线向的投影，属一维形变，存在方向模糊问题，视线方向上表现的正负并不意味着地表的沉陷和隆升。因此，不能根据视线方向观测值识别和分析地表水平和垂直具体形变特征，

近年来，更多的学者开始致力于insar三维形变重建方面的研究，如：基于互相关技术的偏移轨道(offset-tracking)、多孔径(multi-apertureinsar，mai)技术和升降轨技术等。这些方法需要不同几何结构的多轨观测数据来支撑，且没有稳定的专门适用于三维地壳形变的转换方法；尽管通过位错模型的正演可以得到连续的三维形变场，然而其结果主要取决与模型的精细程度和断层先验信息，通常情况是不同学者用不同手段得到的断层滑动参数大相径庭，各执己见，导致生成完全不同的地表三维形变场。因此，融合大地测量与地形模型，研究一个稳健的三维形变场重建方法势在必行。

技术实现要素：

本发明的目的是提供一种基于单轨insar观测、使用于地壳形变的三维转换方法，利用简化的模型亦能够恢复出可靠三维形变场。

本发明提出了一种基于单轨insar数据的三维地壳形变转换方法，包括以下步骤：

1)、基于弹性位错模型，利用单轨insar视线向一维数据和初步或简化的断层几何参数反演断层滑动分布；然后正演得到地表空间连续的三维形变方向；

2)、以步骤1)所得的三维形变方向为约束，将insar视线向一维形变数据转换为n、e、u三维数据，进而得到地壳形变三维形变场。

所述的步骤1)具体包括以下步骤：

11)、确定断层初步的几何参数：根据震源机制解得到；

12)、反演得到断层面滑动分布；

13)、地表三维形变方向参数估计。

所述的步骤12)具体如下所述：

根据弹性均匀、各向同性、半无限空间位错理论和步骤11)得到的断层几何参数，将位错模型反演转化为视线向观测一维数据与断裂滑动分布的线性反演问题，即：

式(1)中，yi(i＝1,...,n)为n个insar采样点视线向向形变观测数据；为m个离散的子断层沿走向和倾向上的滑移分量；g为格林函数；εi为insar视线向观测误差；

首先，根据已知的最大滑移量smax和断层滑动性质先验信息，采用不等式约束来控制断层滑移的大小和方向，使之始终保持在一定范围和一定方向内变动：

式(2)中，为最大、最小滑动角；

其次：将整个断层面离散为许多子断层，并且认为每一个子块内的位错是均匀的，但不同的子块滑动大小和方向不同，地表上任一点的位移是各个子块单元位错造成的位移的叠加；采用拉普拉斯算子约束断层面上的滑动在空间分布上的平滑程度，见式(3)：

联合式(1)、(2)、(3)，建立断层表面上的滑动量与同震地表形变之间的反演模型为：

式(4)中，τ为平滑因子，h为平滑算子；

根据约束最小二乘便可反演得到断层面滑动参数s。

所述的步骤13)具体如下所述：根据得到的断层滑动参数，由位错模型按式(6)正演地表三维形变值dn,de,du，进而根据式(7)计算断层上方地表空间连续的三维形变方向参数cosα，cosβ，cosγ：

式中，α,β,γ是空间形变向量分别与n、e、u三轴的夹角，用三个夹角余弦cosα，cosβ，cosγ表示三维形变方向参数。

所述的步骤2)具体如下所述：

21)、构建insar视线向一维观测与n、e、u三维形变之间的几何模型；

根据insar的成像特点，其形变观测结果是n、e、u三个方向的形变在雷达视线方向上的投影叠加，α为雷达飞行坐标方位角为，θ为雷达侧视角，这些参数由sar影像头文件获取；

根据各向空间几何位置关系，列出insar视线向形变与n、e、u三维形变之间的转换关系式，即：

式(8)中，dn、de、du为地面上任一点n、e、u向三维形变分量大小。dlos，εlos分别代表insar视线向一维形变量和观测误差；

22)、得到转换后的三维形变分量dn、de、du：

根据insar成像几何关系，三维形变分量dn、de、du由地表三维空间总向量s表示为：

式中，将式(8)和式(9)联立得：

由式(10)求得s为：

因此，根据insar一维观测值dlos和获取的地表三维形变方向参数cosα，cosβ，cosγ，结合式(9)和(11)即得到转换后的三维形变分量dn、de、du。

本发明融合insar观测与地学模型于一体，充分利用了断层破裂引起的地表形变具有高空间相关性的特点，尤其是地表形变方向在一定范围内具有较强的空间一致性。即使断层几何参数有所差异，但是一旦断层的运动特性确定，其地表形变整体趋势特征将确定，地表形变方向差异不大。如，走滑断裂引起的地表形变方向平行于断裂，逆冲断裂引起的地表形变方向垂直与断裂。因此，基于断层破裂的地学模型，以稳健的三维形变方向为约束，将insar视线向一维观测转换至三维分量，从而实现了单轨insar数据的三维地壳形变转换。

该三维地壳形变转换方法，稳健可靠。相比传统方法，在初步简化的断层模型情况下，也能够得到可靠的三维形变场。该方法不仅体现了地壳运动的内在的物理意义，而且还反映了断层形变方向的局部细节特征，为地球内部构造变形的反演提供了很重要的边界与约束条件，也为震间连续地表三维形变场的获取与地壳形变特征的进一步分析提供了一个新方法。

附图说明

图1是本发明的方法流程图；

图2是insar视线向观测与n、e、u三维方向的空间几何关系图；

图3是本方法获取的地表三维形变场示意图；

图4是本方法获取的三维形变结果与gnss观测方法的对比图。

具体实施方式

如图1所示，本发明公开了一种基于单轨insar数据的三维地壳形变转换方法，具体包括以下步骤：

一、基于弹性位错模型，利用单轨insar视线向一维数据和初步或简化的断层几何参数反演断层滑动分布；然后正演得到地表空间连续的三维形变方向；

二、以步骤一所得的三维形变方向为约束，将insar视线向一维形变数据转换为n、e、u三维数据，进而得到地壳形变三维形变场。

其中步骤一具体实施时，包括以下步骤：

11)、确定断层初步的几何参数；

断层几何参数例如断层长度、宽度、深度、位置、走向、倾角等，本领域技术人员均可以根据震源机制解得到，在此不再赘述。

12)、反演得到断层面滑动分布；

由于okada弹性位错模型(okada,1985)是描述震源位移与地表位移之间的物理方程，广泛用于断层参数反演。根据弹性均匀、各向同性、半无限空间位错理论和已知的断层几何参数，将位错模型反演转化为视线向观测一维数据与断裂滑动分布的线性反演问题，即

式(1)中，yi(i＝1,...,n)为n个insar采样点视线向向形变观测数据；为m个离散的子断层沿走向和倾向上的滑移分量；g为格林函数；εi为insar视线向观测误差。

首先：为确保估计结果的可靠性和稳定性，附加可靠的先验信息和约束条件至关重要。本实施例可以根据已知的最大滑移量smax和断层滑动性质如走滑、倾滑，左旋、右旋等先验信息，采用不等式约束来控制断层滑移的大小和方向，使之始终保持在一定范围和一定方向内变动：

式(2)中，为最大、最小滑动角。

其次：为得到断层面精细滑动，要考虑到断层的不均匀变形特征，于是将整个断层面离散为许多子断层，并且认为每一个子块内的位错是均匀的，但不同的子块滑动大小和方向不同，地表上任一点的位移是各个子块单元位错造成的位移的叠加。

为避免滑动分布解的振荡，避免相邻子断层滑动量在大小和方向上存在显著差异，保持断层滑动分布在空间上的平稳性，可以采用拉普拉斯算子约束断层面上的滑动在空间分布上的平滑程度，见式(3)：

联合式(1)、(2)、(3)，建立断层表面上的滑动量与同震地表形变之间的反演模型为：

式(4)中，τ为平滑因子，h为平滑算子。

根据约束最小二乘便可反演得到断层面滑动参数s。

13)、地表三维形变方向参数估计；

根据上述得到的断层滑动参数，由位错模型按式(6)正演地表三维形变值dn,de,du，进而根据式(7)计算断层上方地表空间连续的三维形变方向参数cosδ，cosβ，cosγ。

式中，δ,β,γ是空间形变向量分别与n、e、u三轴的夹角，用三个夹角余弦cosδ，cosβ，cosγ表示三维形变方向参数。

本发明所述方法的步骤二具体实施时，包括以下步骤：21)、构建insar视线向一维观测与n、e、u三维形变之间的几何模型；

根据insar的成像特点，其形变观测结果是n、e、u三个方向的形变在雷达视线方向上的投影叠加。如图2所示，α为雷达飞行坐标方位角为，θ为雷达侧视角，这些参数可由sar影像头文件获取。

根据各向空间几何位置关系，可以写列出insar视线向形变与n、e、u三维形变之间的转换关系式，即

式(8)中，dn、de、du为地面上任一点n、e、u向三维形变分量大小。dlos，εlos分别代表insar视线向一维形变量和观测误差。

显然，由式(8)可以容易的将地面三维形变分量dn、de、du转换为insar视线向一维形变dlos。然而，没有外部数据的支持，难以由insar视线向一维形变量dlos转换到地表水平和垂直方向上的三维形变分量dn、de、du。

22)、得到转换后的三维形变分量dn、de、du

根据图2所示的insar成像几何关系，三维形变分量dn、de、du可由地表三维空间总向量s表示为：

式中，

将式(8)和式(9)联立得：

由式(10)可求得s为：

因此，根据insar一维观测值dlos和步骤一所获取的地表三维形变方向参数cosδ，cosβ，cosγ，结合式(9)和(11)即可得到转换后的三维形变分量dn、de、du。

为进一步证实本发明所提方法的有效性，我们用alos/palsar数据，以2008年汶川ms8.0地震同震三维地表形变场的转换为例，从简化-到精细构建了四种截然不同的断层几何模型(m1-m4)，其详细参数如表1所示。

表1

利用本发明提出的方法获取的四种模型的同震三维形变场如图3所示。显然，不同模型之间呈现了较强的一致性，无论模型的精细与否。重建的三维形变场与全球卫星导航定位系统(globalnavigationsatellitesystem，gnss)观测结果比较完全吻合，见图4。试验结果表明，利用本方法，即使采用简化的断层模型也能获取可靠的三维形变场，正确反映出断层破裂所呈现的逆冲兼走滑的运动特性。从重建的三维形变场结果可以清楚地判断断层表面局部细节运行的方向和性质及各种分量(逆冲和走滑)的大小和分布，分析特定方向上更为精细的形变特征。证实了本方法是一种稳健的三维形变估计方法。本方法为揭示大区域近、远场真实的形变特征，理解断层运动方式和地震触发机制提供了有价值的参考，也为震间断层三维蠕滑形变场的建立提供了一个好方法。