CN115856996B - 一种地震数值预测预报方法和系统 - Google Patents

Abstract

本发明公开了一种地震数值预测预报方法和系统，包括获取待预测的监视区的地质、地球物理信息和大地测量信息；根据区域地质、地球物理信息和大地测量信息建立内含断层的有限元数值模型，并获取有限元模型的边界条件、初始应力场；根据有限元数值模型对上述监视区进行反演回溯，模拟上述监视区内历史上发生过的强地震时空序列；在上述最佳模型的边界上进行持续不间断地模拟加载，计算得出监视区未来强震发生的时间、地点及震级。本发明基于力学、数学及地球科学的动力学方法，通过有限元模拟地震的孕育和发生过程，突破目前地震预测只是利用地震前兆的经验方法或基于经验的统计学方法的瓶颈，从而更好地满足实际地震预测预报需求。

Description

一种地震数值预测预报方法和系统

技术领域

本发明属于地震预报技术领域，尤其是涉及一种地震数值预测预报方法和系统。

背景技术

地震是人类面临的重大自然灾害，我国大陆地震活动非常活跃，面临的地震形势十分严峻。目前，国内外地震预测预报都是基于地震前兆的经验方法或基于经验的统计学方法，由于缺乏物理基础，成功率往往很低。

根据弹性回跳理论，构造应力在地震区内缓慢积累，达到一定极限产生地震,地震是岩石的粘滑摩擦失稳过程,大多数构造地震起因于原有断层或板块边界的突然滑动，而很少是突然产生新的剪切破裂。地震从孕育到发生是个非常复杂的力学过程。通常，强震孕育需要数百年甚至几千年的时间；而地震的发生是一刹那，只需数十秒至多上百秒，时间跨度高达10个数量级。可见，数值模拟强震孕育、发生及其循环复发过程非常困难。由于大地震孕育时间长，能量积聚大，而且是在人类目前无法直接探测的地下深处，为此，我们提出一种地震数值预测预报方法来解决上述问题。

发明内容

本发明的目的是针对上述问题，提供一种地震数值预测预报方法。

为达到上述目的，本发明提供如下技术方案：

本发明包括如下步骤：

A获取待预测的监视区的地质、地球物理信息和大地测量信息；

B根据区域地质、地球物理信息和大地测量信息建立内含断层的有限元数值模型，并获取有限元模型的边界条件、初始应力场；

C根据所述有限元数值模型对上述监视区进行反演回溯，模拟上述监视区内历史上发生过的强地震时空序列，获得模拟结果与断裂带上发生的真实历史地震时空序列相一致的最佳模型；

D在上述最佳模型的边界上进行持续不间断地模拟加载后得出所述监视区未来强震发生的时间、地点，再根据模拟给出的断层面上的破裂面积及同震位错分布推出未来地震的震级。

进一步地，所述地质信息包括地震地质构造、断层几何，所述地球物理信息包括地下介质速度结构、密度结构、黏性结构、温度结构、地下流体结构、钻孔测量、应力测量、震源机制解、地震活动性,所述大地测量信息包括GNSS、GPS、InSAR遥感观测信息。

进一步地，所述有限元数值模型的计算方法包括：

a首先利用地质、地球物理及大地测量资料反演有限元模型的初始应力场及边界条件；

b对模型进行加载，计算断裂带上的断层自发破裂过程，模拟给出监视区内历史地震时空分布，获得最佳的数值模型；

c继续对模型加载，计算给出未来断裂带上强震发生的时、空、强三要素。

上述有限元模拟过程，其动力学过程的数学描述如下：

利用有限元方法模拟被断层隔开的两块区域间的破裂动力学过程，其平衡方程如下：

σ_ij,j+f_i-ρu_i,tt-cu_i,t＝0 (1)

上式中σ_ij是应力张量，f_i是体力，ρ是弹性介质的密度，c是阻尼系数，u_i表示位移，u_i,t和u_i,tt分别是u_i对t的一次导数和二次导数，即分别是速度和加速度；-ρu_i,tt和-cu_i,t分别代表惯性力和阻尼力，在线弹性介质中

σ_ij＝λu_k,kδ_ij+μ(u_i,j+u_j,i) (2)

上式中λ和μ为拉梅常数，称为物理方程。

将平衡方程(1)及力的边界条件表示成等效积分形式，可以得到破裂动力学过程的能量表达形式：

上式中t为时间，ρ为弹性介质的密度，v为速率矢量，U单位质量的内能,f为体力，T^l为边界分布荷载，T^qb为吸收边界上吸收辐射能量的阻尼力，T^f为摩擦力，

为外力做功功率，/>

为吸收边界上阻尼耗散功率，/>

为摩擦耗散吸收功率。

离散后的位移表达式代入(3)式，得到求解域内有限元节点位移满足的运动方程：

上式中

和u(t)分别是系统的节点加速度向量、节点速度向量和节点的位移向量。M，C，K和Q(t)分别是系统的质量矩阵、阻尼矩阵、刚度矩阵、和节点荷载向量。式中：

式中F_t(t)和

分别表示构造力和断层面上的接触力，D表示断层上的相对滑移距离，/>

表示断层上的相对滑移速率。实际计算中，运用有限元中的接触技术来模拟断层的运动。

在断层两侧法向位移假设满足连续条件：

u_--u₊＝0 (6)

此外假设断层面上应力满足：

上式中σ_n为法向接触应力，τ为切向应力，μ为摩擦系数。

利用方程(4)～(7)就可以得到地震位移场。

进一步地，根据介质的构造模型，三维有限单元的剖分采用了的8节点6面体单元。

一种地震数值预测预报系统，包括：

数据获取模块，用于获取待预测的监视区的地质、地球物理信息和大地测量信息，

模型建立模块，利用所述区域地质、地球物理信息和大地测量信息，建立有限元数值模型，通过有限元计算，反演获取有限元模型的边界条件及初始应力场；

模拟模块，用于根据所述有限元数值模型对上述监视区进行反演回溯，模拟上述监视区内历史上发生过的强地震时空序列；

预测模块，用于将最佳模型的边界上进行持续不间断地模拟加载后得出所述监视区断层的自发破裂产生的时间和位置，获得未来地震的地点与时间，通过破裂面积及破裂位错推演地震震级。

一种电子设备，包括：一个或多个处理器；存储装置，用于存储一个或多个程序，其中，当所述一个或多个程序被所述一个或多个处理器执行时，使得所述一个或多个处理器执行根据权利要求1～6中任一项所述的基于地震数值预测预报方法。

一种计算机可读存储介质，其上存储有可执行指令，该指令被处理器执行时使处理器执行根据权利要求1～4中任一项所述的基于地震数值预测预报方法。

与现有技术相比，本发明提供了一种地震数值预测预报方法，具备以下有益效果：

本发明完全基于地学、力学和数学方法，通过有限元数值模拟给出未来强震发生的时、空、强三要素，从而从本质上改变了目前地震预测预报利用地震前兆的经验方法，使得地震预报由经验走向科学，预测的准确度将大大提高，能够更好的满足实际预测需求。

附图说明

图1为本发明地震数值预测预报方法的流程示意图；

图2为本发明有限元模拟给出的龙门山断裂带上缓倾角断层及陡倾角断层上的位错随着时间变化分布图。

图3为本发明地震数值预测预报系统的流程示意图；

具体实施方式

以下实施例仅处于说明性目的，而不是想要限制本发明的范围。

如图1所示，本发明包括如下步骤：

A获取待预测的监视区的地质、地球物理信息和大地测量信息；

为了建立包含监视区的有限元数值模型，首先需要获取监视区的多种观测数据，这些数据包括：活动构造、地球物理、大地测量、应力测量、断层几何、波速结构、地形高程，等。此外还要收集前人研究所获得的成果中所给出的结果数据，如：介质的密度结构，黏性结构，电性结构，温度结构，流体分布，等。

B根据区域地质、地球物理信息和大地测量信息建立内含断层的有限元数值模型，并获取有限元模型的边界条件、初始应力场；

根据监视区的空间几何范围，建立包含监视区的三维有限元模型(立方体)，通过约束条件获取三维数值核型的初始应力场及边界条件。为了能够比较准确地模拟应力、应变的时空演化过程以及分辨破裂的凝聚区(cohesive zone),主要断裂带附近的网格尺寸为50米，距离断层远处的网格大小为150米。有限元模型中各个单元中介质的物理性质(如：杨氏模量，泊松比，介质密度，介质孔隙率，含水性，温度，黏度，电导率，屈服强度，热膨胀系数，等等)。这些都是根据对监视区收集的数据或前人的研究结果中获得。

有限元模拟中，除了上述介质的物性参数外，模型的边界条件及初始条件将影响着计算结果。但是，由于震源区通常位于地下10多千米的地下深度，目前人类对震源区及其附近的应力场无法直接观测；此外，地表介质的移动速度可以根据大地测量的结果，但地下深处，模型边界的运动速度无法探测。所以，根据震源机制解、应力测量、地质调查以及大地测量结果(如：GNSS、InSAR、GPS,等)作为约束条件来反演有限元模型的初始应力场及模型的边界运动速度。

C根据所述有限元数值模型对上述监视区进行反演回溯，模拟上述监视区内历史上发生过的强地震时空序列，获得模拟结果与断裂带上发生的真实历史地震时空序列相一致的最佳模型；

基于上述的有限元模型，以边界上反演的速度对模型实施加载，这样模型中的应力、应变会不断变化。因为在不同区域，介质的物理性质不同，应力、应变在空间是不均匀分布的，特别是在断层附近应力容易集中、积累。当断层上的应力达到极限时，产生自发破裂，形成大地震。地震之后，积累的应力得到释放，断层愈合、重新闭锁，但随着边界上的不断加载，势必会再次孕育新的地震。模拟中就是通过这种方式，模拟1000年以来，监视区破裂带上的强震事件(M≥7.0)的时空序列。

通常，上述模拟给出的强震序列与真实的历史上发生的强地震时空分布是不一致的，为此通过调整断层面上的摩擦本构关系及其参数，让模型模拟给出的强震序列与历史地震序列相一致(震级之差小于0.5级；时间之差小于3个月；空间距离相差不超过40千米)。并将此时的有限元模型视为最佳模型。

D在上述最佳模型的边界上进行持续不间断地模拟加载后得出所述监视区断层面上的破裂开始的时间为发震时间，其位置为震源；断层上破裂的面积与位错的积分可以推演地震的震级。

基于上述对历史地震拟合得最好的最佳模型，让模型边界持续不断地加载，这样模型中的应力场在断层附近集中、积累，断层最终会破裂而形成大地震(M≥7.0)。若模拟的地震其震级小于7.0级，则让模型继续计算，直至产生一个M≥7.0的大地震。显然，这个大地震是完全通过有限元模拟计算给出的，其发生的地点、时间、震级都是通过科学计算得到的，与目前运用地震前兆的经验方法或基于前兆的统计学方法完全不同。是一种基于力学和数学的物理方法。有望对未来的大地震进行准确的预测预报，减轻地震造成的灾害和人员伤亡。

由于地震孕育过程非常缓慢、变形速率很小，该过程可视为静态或准静态过程，可用描述静态力学的方程来求解(也可以用动力学方法，视加速度为0)；但地震发生时，断层产生错动，发射地震波，地面产生强烈震动，只能用动力学过程来描述；另外地震在孕育过程中，由于变形小，可以认为介质是连续的，但是在地震发生过程中断层发生相互错动，介质的连续性遭到破坏，此时只能用不连续介质力学的理论来表达。在地震孕育、发生及震后调整过程的连续数值模拟中，不仅涉及力学状态的改变，同时，如何来选择计算中的时间步长，也面临着很大的困难：若按照漫长的地震孕育过程来选择时间步长计算步长为半年或一个月，则无法模拟短暂的地震发生过程；若以短暂的地震发生过程来选择时间步长，则很难模拟漫长的地震孕育过程数千年。因此，模拟过程中，还必须改变计算的时间步长。计算过程能够在静态与动态之间自动转换，计算时间步长可以根据需要自动控制，因此既能计算长期缓慢的地震孕育过程通常百年尺度或更长，又能对地震发生的同震过程即瞬时的断层破裂过程，计算时间步长可短至0.001秒或更小的时间进行模拟。

模拟中，只要有断层面上摩擦本构关系的数学表达式，就可以方便地嵌入到有限元计算的主程序中，在模拟地震破裂过程中，无需进行任何人工干预，如给定成核区域、减小摩擦系数、改变应力状态等进行连续计算。这样利用改进后的有限元方法，将研究地震破裂过程与地震孕育过程综合考虑，减少了人工干预，按照实际地质过程，由缓慢的构造加载获得断层破裂之前的应力状态，克服了人为给定初始应力场的限制；

对于复杂边界条件及非平面断层问题都能够进行灵活处理，可以模拟复杂的介质物性结构如：黏弹、黏塑、孔隙介质等。此外，在计算过程中网格剖分可以实现自适应，即在计算过程中，根据精度要求，实施网格的自动加密。模拟结果显示改进后的有限元方法在模拟地震的孕育和地震破裂方面有独到的优点。

如图2所示，有限元模拟给出的龙门山断裂带上缓倾角断层及陡倾角断层上的位错随着时间变化分布图。图中的每个台阶表示一次地震事件产生突然位错(一次地震)，红色箭头对应着一次大地震，中间的数字是地震的复发间隔，平均约为3298年，该值与古地震探测及地质调查结果有很好的一致性。此外，模拟给出在大地震之间，龙门山断裂带上缓倾角断层上的小地震事件的发生服从滑移可预测模型(slip-predictable)，而缓倾角断层上的较大地震服从时间可预测模型(time-predictable)，但是陡倾角上发生的强地震很复杂，不服从任何预测模型。

本公开的实施例还包括一种计算机程序产品，其包括计算机程序，该计算机程序包含用于执行流程图所示的方法的程序代码。当计算机程序产品在计算机装置中运行时，该程序代码用于使计算机装置实现本公开实施例所提供的预测方法。

在一种实施例中，如图3所示，一种地震数值预测预报系统，包括：

数据获取模块，用于获取待预测的监视区的地质、地球物理信息和大地测量信息，

模型建立模块，用于根据所述区域地质、地球物理信息和大地测量信息建立内含断层的有限元数值模型，并获取有限元模型的边界条件、初始应力场；

模拟模块，用于根据所述有限元数值模型对上述监视区进行反演回溯，模拟上述监视区内历史上发生过的强地震时空序列；

预测模块，用于对最佳模型的边界上进行持续不间断地模拟加载后得出所述监视区地震时断层面上的断层自发破裂，给出发震时间、地点及地震的震级。

在该计算机程序被处理器执行时执行本公开实施例的装置/装置中限定的上述功能。根据本公开的实施例，上文描述的装置、模块、单元等可以通过计算机程序模块来实现。

在一种实施例中，该计算机程序可以依托于光存储器件、磁存储器件等有形存储介质。在另一种实施例中，该计算机程序也可以在网络介质上以信号的形式进行传输、分发，并通过通信部分被下载和安装，和/或从可拆卸介质被安装。该计算机程序包含的程序代码可以用任何适当的网络介质传输，包括但不限于：无线、有线等等，或者上述的任意合适的组合。

在这样的实施例中，该计算机程序可以通过通信部分从网络上被下载和安装，和/或从可拆卸介质被安装。在该计算机程序被处理器执行时，执行本公开实施例的装置中限定的上述功能。根据本公开的实施例，上文描述的设备、装置、模块、单元等可以通过计算机程序模块来实现。

根据本公开的实施例，可以以一种或多种程序设计语言的任意组合来编写用于执行本公开实施例提供的计算机程序的程序代码，具体地，可以利用高级过程和/或面向对象的编程语言、和/或汇编/机器语言来实施这些计算程序。程序设计语言包括但不限于诸如Fortran,Java，C++，Python，“C”语言或类似的程序设计语言。程序代码可以完全地在用户计算设备上执行、部分地在用户设备上执行而部分在远程计算设备上执行、或者完全在远程计算设备或服务器上执行。在涉及远程计算设备的情形中，远程计算设备可以通过任意种类的网络，包括局域网(LAN)或广域网(WAN)，连接到用户计算设备，或者，可以连接到外部计算设备(例如，利用因特网服务提供商来通过因特网连接)。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种地震数值预测预报方法，其特征在于：包括如下步骤：

A获取待预测的监视区的地质、地球物理信息和大地测量信息；

B根据区域地质、地球物理信息和大地测量信息建立内含断层的有限元数值模型，并获取有限元数值模型的边界条件、初始应力场,所述有限元数值模型的计算方法包括：

a首先利用地质、地球物理及大地测量资料反演有限元数值模型的初始应力场及边界条件；

b对有限元数值模型进行加载，计算断裂带上的断层自发破裂过程，模拟给出监视区内历史地震时空分布，获得最佳的数值模型；

c继续对有限元数值模型加载，计算给出未来断裂带上强震发生的时空强三要素；

上述有限元模拟过程，其动力学过程的数学描述如下：

利用有限元方法模拟被断层隔开的两块区域间的破裂动力学过程，其平衡方程如下：

σ_ij,j+f_i-ρu_i,tt-cu_i,t＝0 (1)

上式中σ_ij，j是应力张量，f_i是体力，ρ是弹性介质的密度，c是阻尼系数，u_i表示位移，u_i,t和u_i,tt分别是u_i对t的一次导数和二次导数，即分别是速度和加速度；-ρu_i,tt和-cu_i,t分别代表惯性力和阻尼力，在线弹性介质中

σ_ij＝λu_k,kδ_ij+μ(u_i,j+u_j,i) (2)

上式中λ和μ为拉煤常数；

将平衡方程(1)及力的边界条件表示成等效积分形式，可以得到破裂动力学过程的能量表达形式：

上式中t为时间，ρ为弹性介质的密度，v为速率矢量，U为单位质量的内能,f为体力，T^l为边界分布荷载，T^qb为吸收边界上吸收辐射能量的阻尼力，T^f为摩擦力，

为外力做功功率，

为吸收边界上阻尼耗散功率，/>

为摩擦耗散吸收功率；

离散后的位移表达式代入(3)式，得到求解域内有限元节点位移满足的运动方程：

上式中

和u(t)分别是系统的节点加速度向量、节点速度向量和节点的位移向量，M，C，K和Q(t)分别是系统的质量矩阵、阻尼矩阵、刚度矩阵、和节点荷载向量，式中：

式中F_t(t)和F_f(D(t),

分别表示构造力和断层面上的接触力，D表示断层上的相对滑移距离，/>

表示断层上的相对滑移速率，实际计算中，运用有限元中的接触技术来模拟断层的运动，

在断层两侧法向位移假设满足连续条件：

u_--u₊＝0 (6)

此外假设断层面上应力满足：

上式中σ_n为法向接触应力，τ为切向应力，μ_f为摩擦系数；

利用方程(4)～(7)就可以得到地震位移场；

C根据所述有限元数值模型对上述监视区进行反演回溯，模拟上述监视区内历史上发生过的强地震时、空序列，获得模拟结果与断裂带上发生的历史地震时空序列相一致的最佳模型；

D在上述最佳模型的边界上进行持续不间断地模拟加载后得出所述监视区未来强震发生的时间、地点，再根据模拟给出的断层面上的同震位错分布推出未来地震的震级。

2.根据权利要求1所述的一种地震数值预测预报方法，其特征在于：地质信息包括地震地质构造、断层几何，所述地球物理信息包括地下介质速度结构、密度结构、黏性结构、温度结构、地下流体结构、钻孔测量、应力测量、震源机制解、地震活动性,所述大地测量信息包括GNSS、GPS、InSAR遥感观测信息。

3.根据权利要求1所述的一种地震数值预测预报方法，其特征在于：根据介质的构造与断层的几何模型，三维有限单元的剖分采用了8节点6面体单元。

4.根据权利要求1所述的一种地震数值预测预报方法，其特征在于，包括：所述地震数值预测预报方法的预报系统包括

数据获取模块，用于获取待预测的监视区的地质、地球物理信息和大地测量信息，

模型建立模块，用于根据区域地质、地球物理信息和大地测量信息建立内含断层的有限元数值模型，并获取有限元数值模型的边界条件、初始应力场；

模拟模块，用于根据所述有限元数值模型对上述监视区进行反演回溯，模拟上述监视区内历史上发生过的强地震时、空序列；

预测模块，用于将最佳模型的边界上进行持续不间断地模拟加载后得出所述监视区内，发震断层上发生断层自发破裂，其破裂的时间即为地震的发震时刻，破裂位置为发震地点，由破裂面积和破裂位错可以推算出地震震级。

5.一种电子设备，其特征在于，包括：一个或多个处理器；存储装置，用于存储一个或多个程序，其中，当所述一个或多个程序被所述一个或多个处理器执行时，使得所述一个或多个处理器执行根据权利要求1～4中任一项所述的地震数值预测预报方法。

6.一种计算机可读存储介质，其上存储有可执行指令，该指令被处理器执行时使处理器执行根据权利要求1～4中任一项所述的地震数值预测预报方法。

Images