CN109270590A - 非均匀椭球地球地震和地表载荷库伦应力计算方法 - Google Patents

Abstract

本发明公开了一种非均匀椭球地球地震和地表载荷库伦应力计算方法，利用将地震剪切位错和地表载荷分别等效为有限元计算中的体力项和面力项，再通过局部网格自适应加密技术保证计算精度，此方法可计算地震包括震间、同震和震后，以及地表载荷包括地形剥蚀、冰川消融和水库抽蓄水的加卸载引起的库伦应力变化；具体步骤是：地震剪切位错与地表荷载的等效替代、网格自适应加密、粘弹性maxwell体积分型本构递推实现和横向非均匀椭球形地球的实现，最后通过程序实现自动计算。本发明的有益效果是：其计算更为全面和准确的体现了地球浅表圈层地震和火山库伦应力变化的全过程。

Description

非均匀椭球地球地震和地表载荷库伦应力计算方法

技术领域

本发明涉及地球物理学领域，具体地说，涉及一种用于真实的考虑地形、考虑横向非均匀介质的椭球形地球的地震、地表载荷等过程的弹性-粘弹性库伦应力变化计算方法。主要应用于大地测量学、地震危险性领域，其功能是根据地震静态滑动分布或者地表冰川消融、水库蓄水等地球物理资料，考虑真实地球的地表地形起伏、介质横向非均匀性，计算这些加载/卸载过程引起的周边断层库伦应力变化。

背景技术

地球系统永久变形可以改变周围弹性介质的应力分布，最典型的就是大地震、冰川消融、水库蓄水等过程引起的应力分布往往会对周围断层产生应力加载，库伦应力就是描述断层上应力加载的一个物理量。近年来静态库仑应力通常被用来解释余震的分布、预测地震危险性，讨论地震的触发机制，但目前的库伦应力变化计算方法仍然存在改进空间。

目前国际上广泛使用的地震引起的库伦应力变化计算方法主要基于地震位错理论，将计算得到应力降投影到接收面上得到库伦应力变化。Steketee(1958) 最早把位错理论引入地震学，导出了弹性半无限空间中走滑断层的位移格林函数。Okada(1985,1992)给出了经典的半无限空间介质剪切与张拉断层引起的地表和内部的位移、应变表达式，简称为半空间平面模型。Wang等(2003)发布程序 EDGRN/EDCMP用于计算层状地球模型中地震位错引起的位移和应力，简称为平面分层模型。位错理论的进一步发展是能够考虑地球曲率和分层的球形位错理论，例如，Sun(1992)、Sun和Okubo(1993)发展了考虑重力位的球状分层地球模型位错理论，简称为球面分层模型。但基于这些理论的库伦应力变化计算因为未能考虑实际地球的地表地形起伏、介质横向非均匀性等特性，因此计算均距离实际地球有较大差距，还有很大改进空间。

同时，国际上广泛使用的地表载荷过程引起的库伦应力变化计算主要采用谱方法(Haskell,1935,1936；van Bemmelen and Berlage,1935；McConnell,1968； Peltier,1974；Cathles,1975；Wu and Peltier,1982)。谱方法的基本假设是地球模型径向分层、横向均匀，但当地球粘弹性具有横向非均匀性时会发生模的耦合现象。耦合现象的出现使得谱方法不能处理具有不同类型和幅度横向非均匀的地球模型(汪汉胜和许厚泽，2009)。然而，地震学结果显示壳幔的密度、剪切模量、粘滞度以及厚度等性质均具有横向非均匀性(Vasco and Johnson,1998；Mé gnin and Romanowicz,2000；Kozlovskaya et al.,2008)，这与谱方法的基本假设(横向均匀)相矛盾。

综上所述，位错理论和谱方法的局限性带来的库伦应力变化计算局限性，迫切要求有一种新的计算方法能够计算真实地球地质条件下的库伦应力变化。除此之外，目前的库伦应力变化计算方法，例如Coulomb3.3或者其他有限元等数值计算方法往往只能考虑同震库伦应力变化、或者同震震后库伦应力变化，迫切需要发展一种同时考虑地震循环的同震、震后、震间库伦应力变化和地表载荷过程引起的粘弹性库伦应力变化的计算方法。因此，我们通过等效体力和等效面力的方式集成考虑地震和地表载荷过程，同时通过粘弹性积分型本构关系计算地震循环的同震、震后、震间库伦应力变化和地表载荷过程引起的粘弹性引起的库伦应力变化，并且利用数值方法的优点综合考虑了地球椭率、地形起伏和介质的横向非均匀性。

发明内容

本发明正是为了解决上述传统库伦应力变化计算功能单一、半空间或球分层的基本假设不符合真实地球的局限性，而设计的一种非均匀椭球地球地震和地表载荷库伦应力计算方法。

本发明解决其技术问题所采用的技术方案是：

分别通过等效体力和等效面力的方式考虑地震和地表载荷过程，进一步把体力项装入单元右端项，然后组装总体刚度矩阵和总体右端项，求解大型线性方程组，然后根据位移解的梯度利用网格自适应加密技术保证计算精度。采用粘弹性积分型本构关系粘弹性考虑地震的同震和震后过程，采用backslip方式考虑地震的震间过程。并且通过移动节点的方式实现地球曲率和地表地形起伏，高斯积分点赋值的方法实现介质横向非均匀性，最后利用开源的有限元库和并行求解库进行有限元并行代码编写，实现非均匀椭球地球地震和地表载荷库伦应力计算方法。

非均匀椭球地球地震和地表载荷库伦应力计算方法，利用将地震剪切位错和地表载荷分别等效为有限元计算中的体力项和面力项，再通过局部网格自适应加密技术保证计算精度，此方法可计算地震包括震间、同震和震后，以及地表载荷包括地形剥蚀、冰川消融和水库抽蓄水的加卸载引起的库伦应力变化，具体步骤如下：

步骤一：地震剪切位错与地表荷载的等效替代

包含两类情况的等效替代，对于地震剪切位错采用体力项进行等效，对于地表载荷采用面力项进行等效；

a、地震剪切位错的等效体力替代如下：

对于与地震错动等效的剪切位错，与其等效的体力项如下：

其中，为等效体力，Φ_m表示单元K中第m个形函数，v_j是位错面法向， u(r)_i为位错，C_ijpq为弹性系数，dΣ是单元内部位错面；

b、地表载荷的等效面力替代如下：

地表载荷根据获取的岩石、冰川或者水体的质量数据，和总面积，按照质量载荷除以底面积的方式，得到地表非均匀的面载荷，其中质量载荷是质量乘以重力加速度所得；

步骤二：地震剪切位错与地表荷载处网格自适应加密

网格自适应加密后，引入悬挂节点:

为了满足解的连续性，悬挂节点8和9的解需要满足式[2]：

上述u₈和u₉需要从有限元线性系统里消去，待求解完线性系统后通过回带获得它们的值；悬挂点的约束方程写成式[3]：

u_i ^constrained＝C_iju_j [3]

其中，

此时，有限元线性系统，Au＝b [5]

需要转换为：

u_i ^constrained＝C_iju_j [7]

待求解式[6]后，代入式[7]得到悬挂节点的解；

步骤三：粘弹性maxwell体积分型本构递推实现

各向同性材料的三维粘弹性积分型本构关系写为：

式中，E_ijkl为松弛函数，ε_kl为应变，σ_ij为应力；在小应变假设条件下，在足够小时间增量Δt_n内应变率视为常数，经推导可得应力递推公式:

假设体积形变和剪切形变解耦，则对由1个弹簧和N个Maxwell体链并联的广义Maxwell体，体积模量和剪切模量松弛函数分别用Prony级数表示：

式中，K^∞ _ijkl和G^∞ _ijkl分别为时间无穷时刻的体积模量和剪切模量，N为 Maxwell链数目，λ为松弛时间；认为地球介质体积模量不随时间衰减，即体积变形为弹性，剪切变形为粘弹性；对于Maxwell体，式[10]-[11]表示的剪切模量和体积模量写为：

K_ijkl(t)＝K_ijkl [12]

G_ijkl(t)＝G_ijklexp(-t/λ) [13]

体积形变和剪切形变解耦，公式[5]写为应力球张量和偏张量的形式：

σ(t_n)＝(α(Δt)π_d+π_V)Δε(t_n)+β(Δt)σ_d(t_n-1)+σ_V(t_n-1) [14]

其中，

β(t_n)＝exp(-Δt_n/λ)，

其中，G和K分别为剪切模量和体积模量；

将[10]改写成σ(t_n)＝E'ε'+σ₀' [16]

其中，

E'＝α(Δt)π_d+π_V

ε'＝Δε(t_n)

σ₀'＝β(Δt)σ_d(t_n-1)+σ_V(t_n-1) [17]

式[12]在形式上与带初应力的线弹性本构关系完全相同，不同的是模量和初应力随时间变化；这样，粘弹性问题就转化为带初应力的线弹性本构关系；tn 时刻，用应变增量表示的虚功方程可表示为：

通过式[18]的第3项可以按照等效体力考虑地震位错，第4项可以按照等效面力考虑地表剥蚀、冰川消融等地球表面载荷过程；

步骤四：横向非均匀椭球形地球的实现

在计算模型中添加地形起伏以及介质的横向不均匀性，二者采用不同的方法实现；

对于地形起伏，首先生成没有地形的球形地球，然后选取合适的地形数据，计算出地表每个节点需要挪动的距离，将这些节点挪动到目标位置处；为了避免挪动节点造成的网格畸变，同时对地球内部的节点施加适当的挪动量；具体方法是将地表节点的挪动量看做位移，将这个位移作为地球模型的边界条件，以地球模型为求解区域，求解一个Poisson方程，见式[19]，得到地球内部各节点需要的位移量；Poisson方程解的光滑性保证了使用这种方式添加地形不会造成太大的网格畸变；

对横向不均匀性以及Moho面起伏的处理，不采用移动网格的方式，而是在组建刚度矩阵时，根据所采用的介质参数数据，针对不同的Gauss积分点位置选取对应的参数进行计算。

步骤五：程序实现

由于从底层开始编写有限元程序，工作量巨大，并且调试十分繁复，因此该发明方法的有限元程序的编写主要是在C++开源有限元库deal II平台 (https://www.dealii.org/)上进行。对硬件系统的要求主要是一些常用科学计算库：C++编译器、CMake编译配置工具、线性代数库BLAS和LAPACK、MPI 库以支持程序并行计算、大型线性方程组求解器Trilinos和PETSc、并行网格分区库Metis、C++有限元库deal II、后处理库Paraview，这些库均是开源科学计算库，用户可以很方便的下载并安装。在编写有限元代码时，可以分为以下几个专用类，根据有限元的流程，逐一调用这些类，就可以完成有限元代码的编写。

所述非均匀椭球地球地震和地表载荷库伦应力计算方法，从时间格式和时间增量步上保证数值计算的稳定性；时间格式上，把介质粘弹性问题通过应力递推公式转化为形式上带初应力的线弹性本构关系，其中初应力和弹性模量随时间变化并更新；合理时间增量步选取上，通过遍历所有单元的剪切模量和粘滞系数，得到单元最小松弛时间，并将时间增量设为单元最小松弛时间的0.1倍。

所述非均匀椭球地球地震和地表载荷库伦应力计算方法，横向非均匀椭球地球的实现从网格和介质参数上分两步实现；第一步网格上，实现地表的椭率和地形起伏，先生成球形网格，根据WGS-84椭球地球数据和高程数据计算表面网格节点需要挪动的距离，将这些节点挪动到目标位置处；为了避免挪动节点造成网格畸变，求解Poisson方程进一步得到内部网格节点的新位置；二、实现介质计算参数的非均匀性，根据地球内部波速结构、拉美系数和粘滞系数，针对Gauss积分点对应的经纬度和深度坐标，利用三线性插值公式，直接对单元赋予弹性模量、剪切模量、粘滞系数和松弛时间参数，实现地球介质的三维非均匀性。

所述非均匀椭球地球地震和地表载荷库伦应力计算方法，完整地震循环库伦应力变化的计算通过同震、震后、震间三步实现，第一步，利用地震剪切位错的等效体力方法计算地震剪切位错发生瞬间的同震库伦应力变化；第二步，利用粘弹性应力递推关系，实现粘弹性介质的松弛效应，计算震后库伦应力变化；第三步，利用backslip方法，根据区域断层闭锁反演资料，对断层施加一个负位错速率，负位错速率等于地震发生的位错除以地震的复发周期，其方向与地震发生时刻的同震剪切位错相反，以此得到每年的震间库伦应力变化。

本发明的有益效果是：通过等效体力和等效面力的方式集成考虑地震和地表载荷过程，同时通过粘弹性积分型本构关系计算地震循环的同震、震后、震间库伦应力变化和地表载荷过程引起的粘弹性引起的库伦应力变化，并且利用数值方法的优点综合考虑了地球椭率、地形起伏和介质的横向非均匀性。其计算更为全面和准确的体现了地球浅表圈层地震和火山库伦应力变化的全过程。

附图说明

图1为本发明工作流程框图。

图2为网格自适应加密之后悬挂节点示意图。

图3为计算震间库伦应力变化的backslip方法示意图。

海洋板片(1)、大陆地壳(2)、闭锁区(3)，震间地形变形，地壳出出现缩短和隆升。

图4为Neumann边界条件下粘弹性本构的验证示意图之边界条件描述。

图5为Neumann边界条件下粘弹性本构的验证示意图。

图6为汶川地震网格示意图。

图7为汶川地震后在芦山地震处库仑应力、剪应力和正应力随时间变化示意图。

图8为汶川地震后在九寨沟地震处库仑应力、剪应力和正应力随时间变化示意图。

具体实施方式

下面结合附图和实施例对本发明进一步说明。

如图1-5所示，本发明非均匀椭球地球地震和地表载荷库伦应力计算方法，其特征在于：利用将地震剪切位错和地表载荷分别等效为有限元计算中的体力项和面力项，再通过局部网格自适应加密技术保证计算精度，此方法可计算地震包括震间、同震和震后，以及地表载荷包括地形剥蚀、冰川消融和水库抽蓄水的加卸载引起的库伦应力变化，具体步骤如下：

步骤一：地震剪切位错与地表荷载的等效替代

包含两类情况的等效替代，对于地震剪切位错采用体力项进行等效，对于地表载荷采用面力项进行等效；

a、地震剪切位错的等效体力替代如下：

对于与地震错动等效的剪切位错，与其等效的体力项如下：

其中，为等效体力，Φ_m表示单元K中第m个形函数，v_j是位错面法向， u(r)_i为位错，C_ijpq为弹性系数，dΣ是单元内部位错面；

b、地表载荷的等效面力替代如下：

地表载荷根据获取的岩石、冰川或者水体的质量数据，和总面积，按照质量载荷除以底面积的方式，得到地表非均匀的面载荷，其中质量载荷是质量乘以重力加速度所得；

步骤二：地震剪切位错与地表荷载处网格自适应加密

网格自适应加密后，引入悬挂节点:

为了满足解的连续性，悬挂节点8和9的解需要满足式[2]：

上述u₈和u₉需要从有限元线性系统里消去，待求解完线性系统后通过回带获得它们的值；悬挂点的约束方程写成式[3]：

u_i ^constrained＝C_iju_j [3]

其中，

此时，有限元线性系统，Au＝b [5]

需要转换为：

u_i ^constrained＝C_iju_j [7]

待求解式[6]后，代入式[7]得到悬挂节点的解；

步骤三：粘弹性maxwell体积分型本构递推实现

各向同性材料的三维粘弹性积分型本构关系写为：

式中，E_ijkl为松弛函数，ε_kl为应变，σ_ij为应力；在小应变假设条件下，在足够小时间增量Δt_n内应变率视为常数，经推导可得应力递推公式:

假设体积形变和剪切形变解耦，则对由1个弹簧和N个Maxwell体链并联的广义Maxwell体，体积模量和剪切模量松弛函数分别用Prony级数表示：

式中，K^∞ _ijkl和G^∞ _ijkl分别为时间无穷时刻的体积模量和剪切模量，N为Maxwell链数目，λ为松弛时间；认为地球介质体积模量不随时间衰减，即体积变形为弹性，剪切变形为粘弹性；对于Maxwell体，式[10]-[11]表示的剪切模量和体积模量写为：

K_ijkl(t)＝K_ijkl [12]

G_ijkl(t)＝G_ijklexp(-t/λ) [13]

体积形变和剪切形变解耦，公式[5]写为应力球张量和偏张量的形式：

σ(t_n)＝(α(Δt)π_d+π_V)Δε(t_n)+β(Δt)σ_d(t_n-1)+σ_V(t_n-1) [14]

其中，

β(t_n)＝exp(-Δt_n/λ)，

其中，G和K分别为剪切模量和体积模量；

将[10]改写成σ(t_n)＝E'ε'+σ₀' [16]

其中，

E'＝α(Δt)π_d+π_V

ε'＝Δε(t_n)

σ₀'＝β(Δt)σ_d(t_n-1)+σ_V(t_n-1) [17]

式[12]在形式上与带初应力的线弹性本构关系完全相同，不同的是模量和初应力随时间变化；这样，粘弹性问题就转化为带初应力的线弹性本构关系；tn 时刻，用应变增量表示的虚功方程可表示为：

通过式[18]的第3项可以按照等效体力考虑地震位错，第4项可以按照等效面力考虑地表剥蚀、冰川消融等地球表面载荷过程；

步骤四：横向非均匀椭球形地球的实现

在计算模型中添加地形起伏以及介质的横向不均匀性，二者采用不同的方法实现；

对于地形起伏，首先生成没有地形的球形地球，然后选取合适的地形数据，计算出地表每个节点需要挪动的距离，将这些节点挪动到目标位置处；为了避免挪动节点造成的网格畸变，同时对地球内部的节点施加适当的挪动量；具体方法是将地表节点的挪动量看做位移，将这个位移作为地球模型的边界条件，以地球模型为求解区域，求解一个Poisson方程，见式[19]，得到地球内部各节点需要的位移量；Poisson方程解的光滑性保证了使用这种方式添加地形不会造成太大的网格畸变；

对横向不均匀性以及Moho面起伏的处理，不采用移动网格的方式，而是在组建刚度矩阵时，根据所采用的介质参数数据，针对不同的Gauss积分点位置选取对应的参数进行计算。

所述非均匀椭球地球地震和地表载荷库伦应力计算方法，从时间格式和时间增量步上保证数值计算的稳定性；时间格式上，把介质粘弹性问题通过应力递推公式转化为形式上带初应力的线弹性本构关系，其中初应力和弹性模量随时间变化并更新；合理时间增量步选取上，通过遍历所有单元的剪切模量和粘滞系数，得到单元最小松弛时间，并将时间增量设为单元最小松弛时间的0.1倍。

所述非均匀椭球地球地震和地表载荷库伦应力计算方法，横向非均匀椭球地球的实现从网格和介质参数上分两步实现；第一步网格上，实现地表的椭率和地形起伏，先生成球形网格，根据WGS-84椭球地球数据和高程数据计算表面网格节点需要挪动的距离，将这些节点挪动到目标位置处；为了避免挪动节点造成网格畸变，求解Poisson方程进一步得到内部网格节点的新位置；二、实现介质计算参数的非均匀性，根据地球内部波速结构、拉美系数和粘滞系数，针对Gauss积分点对应的经纬度和深度坐标，利用三线性插值公式，直接对单元赋予弹性模量、剪切模量、粘滞系数和松弛时间参数，实现地球介质的三维非均匀性。

所述非均匀椭球地球地震和地表载荷库伦应力计算方法，完整地震循环库伦应力变化的计算通过同震、震后、震间三步实现，第一步，利用地震剪切位错的等效体力方法计算地震剪切位错发生瞬间的同震库伦应力变化；第二步，利用粘弹性应力递推关系，实现粘弹性介质的松弛效应，计算震后库伦应力变化；第三步，利用backslip方法，根据区域断层闭锁反演资料，对断层施加一个负位错速率，负位错速率等于地震发生的位错除以地震的复发周期，其方向与地震发生时刻的同震剪切位错相反，以此得到每年的震间库伦应力变化。

实施例

2008年汶川矩震级7.9级地震位于巴颜喀拉块体东边界的龙门山断裂带具有明显地震活动分段性，其中南段不仅是川西高原与四川盆地剧烈地形高差过渡带，也是历史和现今地震活动带，汶川地震的同震-震后库伦应力变化计算是验证本发明方法有效性的一个最佳且重要实施例子。

图6为汶川地震的计算网格，很好地显示了模型计算对地表地形和横向非均匀性的考虑。计算时根据中国地震局和USGS(美国地质调查局)等机构发布的断层静态分布，使直接计算地震引起的同震和震后的应力降，然后带入公式得到同震-震后库伦应力变化，采用本方法得到的汶川地震引起同震-震后库伦应力变化计算值与前人方法具有可比性，验证了方法的可靠：与单斌等(2009)以及 Wan和Shen(2010)相比，在鲜水河断裂、东库仑断裂、西秦岭断裂上库仑应力一致。Nalbant和McCloskey(2010)的计算是同震库仑应力变化约0.014MPa，5 年的粘弹性库仑应力变化值约为0.003MPa，和本方法视摩擦系数取为0.8的结果一致。

由于采用粘弹性介质，因此断层面上的库伦应力变化会随着时间的推移发生变化，图7和图8分别给出了芦山地震和九寨沟地震震源附近库仑应力变化随时间的演化关系。

本发明不局限于上述最佳实施方式，任何人在本发明的启示下得出的其他任何与本发明相同或相近似的产品，均落在本发明的保护范围之内。

Claims (4)

1.非均匀椭球地球地震和地表载荷库伦应力计算方法，其特征在于：利用将地震剪切位错和地表载荷分别等效为有限元计算中的体力项和面力项，再通过局部网格自适应加密技术保证计算精度，此方法可计算地震包括震间、同震和震后，以及地表载荷包括地形剥蚀、冰川消融和水库抽蓄水的加卸载引起的库伦应力变化，具体步骤如下：

步骤一：地震剪切位错与地表荷载的等效替代

包含两类情况的等效替代，对于地震剪切位错采用体力项进行等效，对于地表载荷采用面力项进行等效；

a、地震剪切位错的等效体力替代如下：

对于与地震错动等效的剪切位错，与其等效的体力项如下：

其中，为等效体力，Φ_m表示单元K中第m个形函数，v_j是位错面法向，u(r)_i为位错，C_ijpq为弹性系数，dΣ是单元内部位错面；

b、地表载荷的等效面力替代如下：

地表载荷根据获取的岩石、冰川或者水体的质量数据，和总面积，按照质量载荷除以底面积的方式，得到地表非均匀的面载荷，其中质量载荷是质量乘以重力加速度所得；

步骤二：地震剪切位错与地表荷载处网格自适应加密

网格自适应加密后，引入悬挂节点:

为了满足解的连续性，悬挂节点8和9的解需要满足式[2]：

上述u₈和u₉需要从有限元线性系统里消去，待求解完线性系统后通过回带获得它们的值；悬挂点的约束方程写成式[3]：

u_i ^constrained＝C_iju_j [3]

其中，

此时，有限元线性系统，Au＝b [5]

需要转换为：

u_i ^constrained＝C_iju_j [7]

待求解式[6]后，代入式[7]得到悬挂节点的解；

步骤三：粘弹性maxwell体积分型本构递推实现

各向同性材料的三维粘弹性积分型本构关系写为：

式中，E_ijkl为松弛函数，ε_kl为应变，σ_ij为应力；在小应变假设条件下，在足够小时间增量Δt_n内应变率视为常数，经推导可得应力递推公式:

假设体积形变和剪切形变解耦，则对由1个弹簧和N个Maxwell体链并联的广义Maxwell体，体积模量和剪切模量松弛函数分别用Prony级数表示：

式中，和分别为时间无穷时刻的体积模量和剪切模量，N为Maxwell链数目，λ为松弛时间；认为地球介质体积模量不随时间衰减，即体积变形为弹性，剪切变形为粘弹性；对于Maxwell体，式[10]-[11]表示的剪切模量和体积模量写为：

K_ijkl(t)＝K_ijkl [12]

G_ijkl(t)＝G_ijklexp(-t/λ) [13]

体积形变和剪切形变解耦，公式[5]写为应力球张量和偏张量的形式：

σ(t_n)＝(α(Δt)π_d+π_V)Δε(t_n)+β(Δt)σ_d(t_n-1)+σ_V(t_n-1) [14]

其中，

β(t_n)＝exp(-Δt_n/λ)，

其中，G和K分别为剪切模量和体积模量；

将[10]改写成σ(t_n)＝E'ε'+σ₀' [16]

其中，

E'＝α(Δt)π_d+π_V

ε'＝Δε(t_n)

σ₀′＝β(Δt)σ_d(t_n-1)+σ_V(t_n-1) [17]

式[12]在形式上与带初应力的线弹性本构关系完全相同，不同的是模量和初应力随时间变化；这样，粘弹性问题就转化为带初应力的线弹性本构关系；tn时刻，用应变增量表示的虚功方程可表示为：

∫_V(δε)^TE'dV+∫_V(δε)^Tσ₀'dV-∫_V(δu)^TfdV-∫_S(δu)^TbdS＝0 [18]

通过式[18]的第3项可以按照等效体力考虑地震位错，第4项可以按照等效面力考虑地表剥蚀、冰川消融等地球表面载荷过程；

步骤四：横向非均匀椭球形地球的实现

在计算模型中添加地形起伏以及介质的横向不均匀性，二者采用不同的方法实现；

对于地形起伏，首先生成没有地形的球形地球，然后选取合适的地形数据，计算出地表每个节点需要挪动的距离，将这些节点挪动到目标位置处；为了避免挪动节点造成的网格畸变，同时对地球内部的节点施加适当的挪动量；具体方法是将地表节点的挪动量看做位移，将这个位移作为地球模型的边界条件，以地球模型为求解区域，求解一个Poisson方程，见式[19]，得到地球内部各节点需要的位移量；Poisson方程解的光滑性保证了使用这种方式添加地形不会造成太大的网格畸变；

对横向不均匀性以及Moho面起伏的处理，不采用移动网格的方式，而是在组建刚度矩阵时，根据所采用的介质参数数据，针对不同的Gauss积分点位置选取对应的参数进行计算。

2.根据权利要求1所述的非均匀椭球地球地震和地表载荷库伦应力计算方法，其特征在于：从时间格式和时间增量步上保证数值计算的稳定性；时间格式上，把介质粘弹性问题通过应力递推公式转化为形式上带初应力的线弹性本构关系，其中初应力和弹性模量随时间变化并更新；合理时间增量步选取上，通过遍历所有单元的剪切模量和粘滞系数，得到单元最小松弛时间，并将时间增量设为单元最小松弛时间的0.1倍。

3.根据权利要求1所述的非均匀椭球地球地震和地表载荷库伦应力计算方法，其特征在于：横向非均匀椭球地球的实现从网格和介质参数上分两步实现；第一步网格上，实现地表的椭率和地形起伏，先生成球形网格，根据WGS-84椭球地球数据和高程数据计算表面网格节点需要挪动的距离，将这些节点挪动到目标位置处；为了避免挪动节点造成网格畸变，求解Poisson方程进一步得到内部网格节点的新位置；二、实现介质计算参数的非均匀性，根据地球内部波速结构、拉美系数和粘滞系数，针对Gauss积分点对应的经纬度和深度坐标，利用三线性插值公式，直接对单元赋予弹性模量、剪切模量、粘滞系数和松弛时间参数，实现地球介质的三维非均匀性。

4.根据权利要求1所述的非均匀椭球地球地震和地表载荷库伦应力计算方法，其特征在于：完整地震循环库伦应力变化的计算通过同震、震后、震间三步实现，第一步，利用地震剪切位错的等效体力方法计算地震剪切位错发生瞬间的同震库伦应力变化；第二步，利用粘弹性应力递推关系，实现粘弹性介质的松弛效应，计算震后库伦应力变化；第三步，利用backslip方法，根据区域断层闭锁反演资料，对断层施加一个负位错速率，负位错速率等于地震发生的位错除以地震的复发周期，其方向与地震发生时刻的同震剪切位错相反，以此得到每年的震间库伦应力变化。