一种计算弹性地球内部同震变形的方法和系统与流程

本申请涉及地球物理学技术领域，特别涉及一种计算弹性地球内部同震变形的方法和系统。

背景技术：

地震的发生一般会伴随着孕震、发震、应力调整等过程，并开始新的地震的孕育，研究地震的孕育机理、破裂过程、震后调整等一直都是科学家们最重视的问题，这可以帮助我们更深刻的认识地震并为减少地震伤害提供科学依据。1958年科学家开始研究地震引起的同震位移和应变变化集中在地球表面的变形计算，科学家对地球模型的认识经历了“均质半无限空间模型—均质球模型—层状半无限空间模型/层状球模型”的过程，对地震引起的地表同震变形计算方法也随着地球模型的演化而改进。

随着地震学的发展以及对孕震特征的研究加深，地震产生的地球内部位移和应变变化也引起了关注，研究地球内部同震变形的特征对地球深部运动及对下一次地震的孕育过程及机理的分析具有重要意义，然而关于计算地球内部同震位移和应变变化的方法却十分少，目前有日本地震学家okada提出的均质半无限空间模型下的解析解计算方法，该模型下的计算非常简便快捷，直接采用解析解的计算方法，就可以快速地计算出任意点源和有限断层引起的地球内部同震变化，但该方法采用的地球模型并不是球形的。日本学者takagi和okada提出的均质球模型下的数值计算方法，该方法是通过应力-应变关系、泊松方程和平衡方程，并引入面球谐函数、缔合勒让德函数，推导出以球谐函数形式表示的地表位移和应变计算公式，它是数值解的计算方法，需要计算到高阶以满足精度要求，尤其对接近于地表的震源需要计算的阶数非常高，从而计算量比较大，比较费时间，而且当最高阶达到一定程度时容易出现数值震荡。两种方法都存在一定的缺陷。

技术实现要素：

(一)申请目的

基于此，为了能真实地计算出任意点源引起的地球内部任意位置处的同震变形，避免纯数值计算方法的高阶震荡问题，保证计算精度的同时节省计算时间，以及更深刻的认识地震并为减少地震伤害提供科学依据，本申请公开了以下技术方案。

(二)技术方案

作为本申请的第一方面，本申请公开了一种计算弹性地球内部同震变形的方法，包括：

将应力场球函数和位移场球函数带入由应力场、位移场组建的平衡方程和泊松方程，化简计算得到球型微分方程组和环型微分方程组，进而算出球型基本解和环型基本解；

将球型基本解和环型基本解带入地表边界条件和独立点源的震源函数中，得到非齐次方程组，进而算出球型解待定系数和环型解待定系数；

根据球型解待定系数和环型解待定系数及计算面深度与震源深度之间的大小关系，算出地球内部不同位置的球型变形值和环型变形值；

将球型变形值和环型变形值带入非齐次方程组，计算不同独立震源的位错love数解析解；

将位错love数解析解带入球函数积分求和公式，得到位移格林函数和应变格林函数；

定义独立点源的震级因子为uds/r²＝1，将实际震级大小带入所述位移格林函数和应变格林函数，得到待计算震源点的弹性地球内部同震变形，包括同震位移和同震应变，其中，u为位移量，s为面积，r为半径；

其中，所述独立点源有四种，包括走滑点源、倾滑点源、水平引张点源和垂直引张点源。

在一种可能的实施方式中，由应力场、位移场组建的泊松方程和平衡方程如下：

其中，f是弹性球表面的单位点力，i是单位张量，上标t表示转置，μ和λ是震源处的弹性介质常数，即拉梅常数。τ为应力张量，u为位移场，ρ为密度。

在一种可能的实施方式中，位移场球函数、应变场球函数、单位点力球函数如下：

其中，为位移场，为应力张量，是应力在径向的分量，球谐函数和是缔合勒让德函数的函数，θ是余纬，是经度，r是地球半径，代表单位点力f的位置，“*”表示复共轭，v是单位矢量，t代表环型变形。y1至y4是球型变形因子，y1和y3是位移的径向与水平向分量，y2和y4是应力的径向和水平向分量，和是环型变形因子中的水平位移和应力分量。

在一种可能的实施方式中，得到球型微分方程组和环型微分方程组，进而算出球型基本解和环型基本解，其中，得到的球型微分方程组如下：

其中，y1至y4是球型变形因子，y1和y3是位移的径向与水平向分量，y2和y4是应力的径向和水平向分量；r为地球半径，μ和λ是震源处的弹性介质常数，β＝λ+2μ，δ代表变化量，n代表阶数，r0是单位点源的位置，

得到的环型微分方程组如下：

其中，和是环型变形因子中的水平位移和应力分量，上标t表示环型解部分，是缔合勒让德函数的函数，代表单位点力f的位置，*表示复共轭，v是单位矢量，

球型基本解如下：

环型基本解如下：

在一种可能的实施方式中，地表边界条件和独立点源的震源函数如下：

y2(r)|r＝r＝y4(r)|r＝r＝0

y(r)|r＝0＜+∞

其中，s代表4种独立点源的震源函数，包括：走滑点源、倾滑点源、水平引张点源和垂直引张点源。

在一种可能的实施方式中，非齐次方程组如下：

其中，r是地球半径，rs＝(r-r0)/r代表正则化的震源深度，βi(i＝1,2,…,6)是球型解系数；

其中，βi^t(i＝1,2,3)是环型解系数。

在一种可能的实施方式中，球型变形值和环型变形值如下：

其中，当计算面深度小于震源深度时，球型变形值和环型变形值如下：

当计算面深度大于震源深度时，球型变形值和环型变形值如下：

在一种可能的实施方式中，利用球型变形值和环型变形值计算位错love数后再利用球函数积分求和得到位移格林函数和应变格林函数位移格林函数和应变格林函数；

求得的独立震源的位错love数解析解如下：

其中，ij＝12,32,220,33，代表4种独立震源；

得到的位移格林函数如下：

得到的应变格林函数如下：

作为本申请的第二方面，本申请还公开了一种计算弹性地球内部同震变形的系统，其特征在于，包括：

基本解计算模块用于将应力场球函数和位移场球函数带入由应力场、位移场组建的平衡方程和泊松方程，化简计算得到球型微分方程组和环型微分方程组，进而算出球型基本解和环型基本解；

待定系数计算模块用于将球型基本解和环型基本解带入地表边界条件和独立点源的震源函数中，得到非齐次方程组，进而算出球型解待定系数和环型解待定系数；

变形值计算模块用于根据球型解待定系数和环型解待定系数及计算面深度与震源深度之间的大小关系，算出地球内部不同位置的球型变形值和环型变形值；

位错love数计算模块用于将球型变形值和环型变形值带入非齐次方程组，计算不同独立震源的位错love数解析解；

格林函数计算模块用于将位错love数解析解带入球函数积分求和公式，得到位移格林函数和应变格林函数，

同震变形计算模块用于定义独立点源的震级因子为uds/r²＝1，将实际震级大小带入所述位移格林函数和应变格林函数，得到待计算震源点的弹性地球内部同震变形，包括同震位移和同震应变，其中，u为位移量，s为面积，r为半径；

其中，所述独立点源有四种，包括走滑点源、倾滑点源、水平引张点源和垂直引张点源。

在一种可能的实施方式中，由应力场、位移场组建的泊松方程和平衡方程如下：

其中，f是弹性球表面的单位点力，i是单位张量，上标t表示转置，μ和λ是震源处的弹性介质常数，即拉梅常数。τ为应力张量，u为位移场，ρ为密度。

在一种可能的实施方式中，位移场球函数、应变场球函数、单位点力球函数如下：

其中，为位移场，为应力张量，是应力在径向的分量，球谐函数和是缔合勒让德函数的函数，θ是余纬，是经度，r是地球半径，代表单位点力f的位置，“*”表示复共轭，v是单位矢量，t代表环型变形。y1至y4是球型变形因子，y1和y3是位移的径向与水平向分量，y2和y4是应力的径向和水平向分量，和是环型变形因子中的水平位移和应力分量。

在一种可能的实施方式中，得到球型微分方程组和环型微分方程组，进而算出球型基本解和环型基本解，其中，得到的球型微分方程组如下：

其中，y1至y4是球型变形因子，y1和y3是位移的径向与水平向分量，y2和y4是应力的径向和水平向分量；r为地球半径，μ和λ是震源处的弹性介质常数，β＝λ+2μ，δ代表变化量，n代表阶数，r0是单位点源的位置，

得到的环型微分方程组如下：

其中，和是环型变形因子中的水平位移和应力分量，上标t代表环型解部分，是缔合勒让德函数的函数，代表单位点力f的位置，“*”表示复共轭，v是单位矢量，

球型基本解如下：

环型基本解如下：

在一种可能的实施方式中，地表边界条件和独立点源的震源函数如下：

y2(r)|r＝r＝y4(r)|r＝r＝0

y(r)|r＝0＜+∞

其中，s代表4种独立点源的震源函数，包括：走滑点源、倾滑点源、水平引张点源和垂直引张点源。

在一种可能的实施方式中，非齐次方程组如下：

其中，r是地球半径，rs＝(r-r0)/r代表正则化的震源深度，βi(i＝1,2,…,6)是球型解系数；

其中，βi^t(i＝1,2,3)是环型解系数。

在一种可能的实施方式中，球型变形值和环型变形值如下：

其中，当计算面深度小于震源深度时，球型变形值和环型变形值如下：

当计算面深度大于震源深度时，球型变形值和环型变形值如下：

在一种可能的实施方式中，利用球型变形值和环型变形值计算位错love数后再利用球函数积分求和得到位移格林函数和应变格林函数位移格林函数和应变格林函数；

求得的独立震源的位错love数解析解如下：

其中，ij＝12,32,220,33，代表4种独立震源；

得到的位移格林函数如下：

得到的应变格林函数如下：

(三)有益效果

本申请公开的计算弹性地球内部同震位移和应变的方法及系统，通过一种半解析半数值的计算方法，能真实地计算出任意点源引起的地球内部任意位置处的同震变形，避免了纯数值计算方法的高阶震荡问题，保证计算精度的同时节省了计算时间，同时更深刻的认识了地震并为减少地震伤害提供了科学依据。

附图说明

以下参考附图描述的实施例是示例性的，旨在用于解释和说明本申请，而不能理解为对本申请的保护范围的限制。

图1是本申请公开的计算弹性地球内部同震变形方法实施例的流程图。

图2是4种独立震源的模型图。

图3是本申请公开的计算弹性地球内部同震变形方法的具体流程图。

图4是走滑点源(strike-slip)在地球内部不同深处的变形值。

图5是倾滑点源(dip-slip)在地球内部不同深处的变形值。

图6是水平引张点源(horizontaltensile)在地球内部不同深处的变形值。

图7是垂直引张点源(verticaltensile)在地球内部不同深处的变形值。

图8是4种独立点震源在地球内部0km–40km处产生的位移和应变变化。

图9是本申请公开的计算弹性地球内部同震变形系统实施例的结构框图。

具体实施方式

为使本申请实施的目的、技术方案和优点更加清楚，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行更加详细的描述。

下面参考图1详细描述本申请公开的一种计算弹性地球内部同震变形的方法实施例。如图1所示，本实施例公开的方法主要包括有以下步骤100至步骤600。

步骤100，将应力场球函数和位移场球函数带入由应力场、位移场组建的平衡方程和泊松方程，化简计算得到球型微分方程组和环型微分方程组，进而算出球型基本解和环型基本解。

弹性地球内部产生位错从而引起地表及地球内部不同深度面上的震荡变形，弹性球表面的单位点力产生同震位移场和应力场，将位移场应力张量写成球函数，将其以球函数形式带入泊松方程和平衡方程，

泊松方程为：

其中，δ代表拉普拉斯算子，而和f可以是在流形上的实数或复数值的方程。

平衡方程是力系平衡条件的数学形式。空间任意力系的平衡条件是，力系的主矢和主矩都等于零，即r＝0，mo＝0

平面任意力系的平衡方程为：

∑fx＝0，∑fy＝0，∑mo(f)＝0

平面汇交力系的平衡方程为：

∑fx＝0，∑fy＝0

在至少一种实施方式中，由应力场、位移场组建的泊松方程和平衡方程如下：

其中，f是弹性球表面的单位点力，i是单位张量，上标t表示转置，μ和λ是震源处的弹性介质常数，即拉梅常数。τ为应力张量，u为位移场，ρ为密度。

在至少一种实施方式中，位移场球函数、应变场球函数、单位点力球函数如下：

其中，为位移场，为应力张量，是应力在径向的分量，球谐函数和是缔合勒让德函数的函数，θ是余纬，是经度，r是地球半径，代表单位点力f的位置，“*”表示复共轭，v是单位矢量，t代表环型变形。y1至y4是球型变形因子，y1和y3是位移的径向与水平向分量，y2和y4是应力的径向和水平向分量，和是环型变形因子中的水平位移和应力分量。

在至少一种实施方式中，得到球型微分方程组和环型微分方程组，进而算出球型基本解和环型基本解，其中，得到的球型微分方程组如下：

其中，y1至y4是球型变形因子，y1和y3是位移的径向与水平向分量，y2和y4是应力的径向和水平向分量；r为地球半径，μ和λ是震源处的弹性介质常数，β＝λ+2μ，δ代表变化量，n代表阶数，r0是单位点源的位置，

得到的环型微分方程组如下：

其中，和是环型变形因子中的水平位移和应力分量，上标t表示环型解部分，是缔合勒让德函数的函数，代表单位点力f的位置，*表示复共轭，v是单位矢量，

球型基本解如下：

环型基本解如下：

步骤200，将球型基本解和环型基本解带入地表边界条件和独立点源的震源函数中，得到非齐次方程组，进而算出球型解待定系数和环型解待定系数。

边界条件是指在求解区域边界上所求解的变量或其导数随时间和地点的变化规律。边界条件是控制方程有确定解的前提，对于任何问题，都需要给定边界条件。边界条件的处理，直接影响了计算结果的精度。而解微分方程要有定解，就一定要引入条件，这些附加条件称为定解条件。

在至少一种实施方式中，地表边界条件和独立点源的震源函数如下：

y2(r)|r＝r＝y4(r)|r＝r＝0

y(r)|r＝0＜+∞

其中，s代表4种独立点源的震源函数，包括：走滑点源、倾滑点源、水平引张点源和垂直引张点源。

4种独立震源的模型图如图2所示，从左到右依次为走滑点源、倾滑点源、水平引张点源、垂直引张点源。

在至少一种实施方式中，非齐次方程组如下：

其中，r是地球半径，rs＝(r-r0)/r代表正则化的震源深度，βi(i＝1,2,…,6)是球型解系数；

其中，βi^t(i＝1,2,3)是环型解系数。

步骤300，根据球型解待定系数和环型解待定系数及计算面深度与震源深度之间的大小关系，算出地球内部不同位置的球型变形值和环型变形值。

震源深度是指从震源到地面(震中)的垂直距离。

计算面深度与震源深度之间的大小关系分为两种，第一种为计算面深度小于震源深度，即计算面在震源上方，处于地表与震源之间；第二种为计算面深度大于震源深度，及计算面在震源下方，处于震源与核幔边界之间。不同位置处的内部同震变形计算表达式是不同的，需要分情况计算地球内部的球型和环型变形。

在至少一种实施方式中，球型变形值和环型变形值如下：

其中，当计算面深度小于震源深度时，球型变形值和环型变形值如下：

当计算面深度大于震源深度时，球型变形值和环型变形值如下：

步骤400，将球型变形值和环型变形值带入非齐次方程组，计算不同独立震源的位错love数解析解。

位错love数：位错love数是描述点震源作用下的地球弹性变形而引入的无量纲参数，包括h(与径向位移有关)、l(与水平向位移有关)、k(与引力位有关)，本方法中采用的地球模型不考虑重力的影响，故未有位错love数k的相关计算。

根据计算面深度和震源深度之间不同关系对应的球形变形值和环型变性值，计算不同独立震源的位错love数解析解。

在至少一种实施方式中，求得的独立震源的位错love数解析解如下：

其中，ij＝12,32,220,33，代表4种独立震源。

步骤500，将位错love数解析解带入球函数积分求和公式，得到位移格林函数和应变格林函数。

根据步骤400中得到的位错love数解析解，再利用球函数积分求和得到地球内部4种独立点源在任意层面上产生的位移格林函数和应变格林函数。

在至少一种实施方式中，得到的位移格林函数如下：

得到的应变格林函数如下：

步骤600，定义独立点源的震级因子为uds/r²＝1，将实际震级大小带入所述位移格林函数和应变格林函数，得到待计算震源点的弹性地球内部同震变形，包括同震位移和同震应变，其中，u为位移量，s为面积，r为半径。

该套格林函数是由独立点源直接引起的，我们定义该独立点源的震级因子为uds/r²＝1(即，位移量×面积/半径²＝1)。根据步骤500中得到的位移格林函数和应变格林函数，带入实际震级大小，可得到待计算震源点的弹性地球内部同震位移变形和应变变形。

具体的，参照图3，该方法完整计算步骤为：

根据弹性球表面的单位点力产生的同震位移场、应力场建立应力-应变关系、泊松方程和平衡方程，引入位移场球函数和应力场球函数组建微分方程组，随后求解上述方程组，得到球型基本解(解析解)和环型基本解(解析解)。将球型基本解和环型基本解引入地表边界条件和独立点源的震源函数中，组建非齐次方程组，求解球型解和环型解未知系数,。将解得的球型解和环型解带回非齐次方程组中，计算位移分量和应力分量，通过判断计算面深度(h)与震源深度(ds)的关系，对于地球内部震源上方(h<ds)和震源下方(h>ds)不同位置处的内部同震变形值利用不同的表达式分别计算，再根据计算结果得到4种独立震源的位错love数解析解。最后根据得到的位错love数解析解利用球函数积分求和得到地球内部4种独立点源在任意层面上产生的位移格林函数和应变格林函数。利用得到的位移格林函数和应变格林函数计算待计算震源点的弹性地球内部同震位移变形和应变变形。

以地球内部20km处的走滑点震源为例进行计算：

为展示4种独立点震源在地球内部不同深度产生的地球内部变形，我们画出位移和应变在不同震中角距(θ)范围的计算结果，分别计算它们在地表(h0)，地下2km(h2)、12km(h12)、28km(h28)和40km(h40)处的内部变形值，如图4-7所示；

为展示地球内部变形随深度变化的结果，我们分别给出4种独立震源的位移(displacement)和应变(strain)分量计算结果，地球内部变形值随震中距的增加而衰减，故我们只给出了地表(0km)到40km深处的计算结果，如图8所示，

下面参考图9详细描述本申请公开的一种计算弹性地球内部同震变形系统实施例。本实施例用于实施前述一种计算弹性地球内部同震变形方法实施例。如图9所示，本实施例公开的系统包括：

基本解计算模块用于将应力场球函数和位移场球函数带入由应力场、位移场组建的平衡方程和泊松方程，化简计算得到球型微分方程组和环型微分方程组，进而算出球型基本解和环型基本解；

待定系数计算模块用于将球型基本解和环型基本解带入地表边界条件和独立点源的震源函数中，得到非齐次方程组，进而算出球型解待定系数和环型解待定系数；

变形值计算模块用于根据球型解待定系数和环型解待定系数及计算面深度与震源深度之间的大小关系，算出地球内部不同位置的球型变形值和环型变形值；

位错love数计算模块用于将球型变形值和环型变形值带入非齐次方程组，计算不同独立震源的位错love数解析解；

格林函数计算模块用于将位错love数解析解带入球函数积分求和公式，得到位移格林函数和应变格林函数；

同震变形计算模块定义独立点源的震级因子为uds/r²＝1，将实际震级大小带入所述位移格林函数和应变格林函数，得到待计算震源点的弹性地球内部同震变形，包括同震位移和同震应变，其中，u为位移量，s为面积，r为半径；

其中，所述独立点源有四种，包括走滑点源、倾滑点源、水平引张点源和垂直引张点源。

在至少一种实施方式中，由应力场、位移场组建的泊松方程和平衡方程如下：

其中，f是弹性球表面的单位点力，i是单位张量，上标t表示转置，μ和λ是震源处的弹性介质常数，即拉梅常数。τ为应力张量，u为位移场，ρ为密度。

在至少一种实施方式中，位移场球函数、应变场球函数、单位点力球函数如下：

其中，为位移场，为应力张量，是应力在径向的分量，球谐函数和是缔合勒让德函数的函数，θ是余纬，是经度，r是地球半径，代表单位点力f的位置，*表示复共轭，v是单位矢量，t代表环型变形。y1至y4是球型变形因子，y1和y3是位移的径向与水平向分量，y2和y4是应力的径向和水平向分量，和是环型变形因子中的水平位移和应力分量。

在至少一种实施方式中，得到球型微分方程组和环型微分方程组，进而算出球型基本解和环型基本解，其中，得到的球型微分方程组如下：

其中，y1至y4是球型变形因子，y1和y3是位移的径向与水平向分量，y2和y4是应力的径向和水平向分量；r为地球半径，μ和λ是震源处的弹性介质常数，β＝λ+2μ，δ代表变化量，n代表阶数，r0是单位点源的位置，

得到的环型微分方程组如下：

其中，和是环型变形因子中的水平位移和应力分量，上标t表示环型解部分，是缔合勒让德函数的函数，代表单位点力f的位置，*表示复共轭，v是单位矢量，

球型基本解如下：

环型基本解如下：

在至少一种实施方式中，地表边界条件和独立点源的震源函数如下：

y2(r)|r＝r＝y4(r)|r＝r＝0

y(r)|r＝0＜+∞

其中，s代表4种独立点源的震源函数，包括：走滑点源、倾滑点源、水平引张点源和垂直引张点源。

在至少一种实施方式中，非齐次方程组如下：

其中，r是地球半径，rs＝(r-r0)/r代表正则化的震源深度，βi(i＝1,2,…,6)是球型解系数；

其中，βi^t(i＝1,2,3)是环型解系数。

在至少一种实施方式中，球型变形值和环型变形值如下：

其中，当计算面深度小于震源深度时，球型变形值和环型变形值如下：

当计算面深度大于震源深度时，球型变形值和环型变形值如下：

在至少一种实施方式中，利用球型变形值和环型变形值计算位错love数后再利用球函数积分求和得到位移格林函数和应变格林函数位移格林函数和应变格林函数，

求得的独立震源的位错love数解析解如下：

其中，ij＝12,32,220,33，代表4种独立震源。

得到的位移格林函数如下：

得到的应变格林函数如下：

本文中的模块的划分仅仅是一种逻辑功能的划分，在实际实现时可以有其他的划分方式，例如多个模块可以结合或集成于另一个系统中。作为分离部件说明的模块在物理上可以是分开的，也可以是不分开的。

以上，仅为本申请的具体实施方式，但本申请的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本申请揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本申请的保护范围之内。因此，本申请的保护范围应以权利要求的保护范围为准。