CN112230277A

CN112230277A - 基于柱坐标系的隧道地震波传播数值模拟方法及系统

Info

Publication number: CN112230277A
Application number: CN202011065264.7A
Authority: CN
Inventors: 刘斌; 李铎; 任玉晓; 许新骥
Original assignee: Shandong University
Current assignee: Shandong University
Priority date: 2020-09-30
Filing date: 2020-09-30
Publication date: 2021-01-15
Anticipated expiration: 2040-09-30
Also published as: CN112230277B

Abstract

本公开提供了一种基于柱坐标系的隧道地震波传播数值模拟方法及系统，所述方法为：获取待模拟隧道的基础数据；根据获取的基础数据设定初始参数；根据设定的初始参数，基于柱坐标系，在径向和轴向上采取规则网格剖分，在周向上进行变尺度细化网格的剖分，根据变尺度细化网格构建弹性波波动方程的有限差分格式；根据获取的初始参数和得到的弹性波波动方程的有限差分格式进行地震波场延拓，得到模拟波场和地震记录；本公开针对现有柱坐标系数值模拟中存在的牺牲计算量来避免数值频散和中心点数值奇异性问题，分别提出了变尺度加密网格策略和基于波场分离思想的“奇点值”近似算法，实现了隧道真三维地震波传播的较真实模拟。

Description

基于柱坐标系的隧道地震波传播数值模拟方法及系统

技术领域

本公开涉及地震波模拟技术领域，特别涉及一种基于柱坐标系的隧道地震波传播数值模拟方法及系统。

背景技术

本部分的陈述仅仅是提供了与本公开相关的背景技术，并不必然构成现有技术。

目前交通建设发展进入了一个科学、快速的发展阶段，隧道工程作为其中的节点性工程，发展规模不断扩大，但同时，隧道工程建设不断向更加复杂的地质条件区域延伸，如何有效地保障隧道施工过程中的安全问题成为了重要的研究课题。为了提前获知前方不良地质体信息，减少工程地质灾害事故的发生，隧道施工期超前预报技术越来越受到人们的重视。在各类隧道地质超前预报技术当中，地震波法以期探测距离较远、界面识别效果较好的优点，成为了一种在实际工程中应用最广泛的地质预报方法。开展隧道地震波传播数值模拟方法的研究对于实际工程中的应用起到了重要的指导作用。

本公开发明人发现，目前隧道地震波传播数值模拟方法的研究以隧道二维环境为主，且局限于笛卡尔坐标系下，将隧道近似为长方体和矩形，没有考虑隧道“类圆柱体”的实际形状；而在隧道三维地震波传播数值模拟方面，多数研究是将二维笛卡尔坐标系下的地震波模拟方法进行升维处理，存在着隧道角点绕射干扰以及网格精度不足等问题，影响数值模拟精度；而基于柱坐标系下的地震模拟方法目前研究较少，且为了简化处理，均将隧道空间视为轴对称问题，将三维问题简化为二维问题，降低了模拟精度，另外为了避免柱坐标系均匀网格划分中存在的数值频散问题，加密了网格，牺牲了计算量来避免数值频散；且并未针对算法本身特有的“奇点”问题，提出有效的解决办法，造成极轴处绕射影响严重，从而影响成像结果。

发明内容

为了解决现有技术的不足，本公开提供了一种基于柱坐标系的隧道地震波传播数值模拟方法及系统，有效解决了传统笛卡尔坐标系中因采用正方形网格所造成的网格精度不足和掌子面绕射干扰严重的问题，同时针对现有柱坐标系数值模拟中存在的牺牲计算量来避免数值频散和中心点数值奇异性问题，分别提出了变尺度加密网格策略和基于波场分离思想的“奇点值”近似算法，实现了隧道真三维地震波传播的较真实模拟。

为了实现上述目的，本公开采用如下技术方案：

本公开第一方面提供了一种基于柱坐标系的隧道地震波传播数值模拟方法。

一种基于柱坐标系的隧道地震波传播数值模拟方法，包括以下步骤：

获取待模拟隧道的基础数据；

根据获取的基础数据设定初始参数；

根据设定的初始参数，基于柱坐标系，在径向和轴向上采取规则网格剖分，在周向上进行变尺度细化网格的剖分，根据变尺度细化网格构建弹性波波动方程的有限差分格式；

根据获取的初始参数和得到的弹性波波动方程的有限差分格式进行地震波场延拓，得到模拟波场和地震记录。

作为可能的一些实现方式，待模拟隧道的基础数据至少包括隧道的形状、围岩和隧道前方区域的弹性参数以及采用的观测装置形式、震源主频、采样频率和采样时长。

作为可能的一些实现方式，在周向上进行变尺度细化网格的剖分，具体为：

根据剖分网格在径向的长度和波场中的最小波长确定最里层初始网格数N1；

判断当前周向的空间差分间隔是否符合频散条件，如果符合则保持当前网格剖分数；

如果不符合，则网格加密系数n加1，第i层的网格数为N_i＝2ⁿN₁。

作为可能的一些实现方式，构建柱坐标系下的初始弹性波方程，利用位移对时间的二阶偏导转换为质点振动速度对时间的一阶偏导数，将初始弹性波方程降阶为一阶速度应力方程组，随后在空间和时间上均采用二阶的交错网格离散方法，得到差分格式。

作为可能的一些实现方式，对网格加密分界线两侧区域的参数进行插值处理，转换为常规网格差分。

作为进一步的限定，对第一参数进行差值处理，具体为：

切取k+1/2平面内的部分网格，先利用相邻两个网格中的第二参数分量进行线性插值，获得虚位的第二参数分量，得到差值后的第二参数分量，进而得到细化交界线处的第一参数差值结果。

作为可能的一些实现方式，将柱坐标系下纵波波场中的奇点值和相近网格上的纵波波场在周向的速度分量建立联系，利用环向速度矢量来插值计算出奇点值。

作为可能的一些实现方式，利用柱坐标系电磁波有限差分中的奇点处理方法计算横波波场的奇点值。

本公开第二方面提供了一种基于柱坐标系的隧道地震波传播数值模拟系统。

一种基于柱坐标系的隧道地震波传播数值模拟系统，包括：

数据获取模块，被配置为：获取待模拟隧道的基础数据；

初始参数设定模块，被配置为：根据获取的基础数据设定初始参数；

弹性波波动方程构建模块，被配置为：根据设定的初始参数，基于柱坐标系，在径向和轴向上采取规则网格剖分，在周向上进行变尺度细化网格的剖分，根据变尺度细化网格构建弹性波波动方程的有限差分格式；

传播数值模拟模块，被配置为：根据获取的初始参数和得到的弹性波波动方程的有限差分格式进行地震波场延拓，得到模拟波场和地震记录。

本公开第三方面提供了一种介质，其上存储有程序，该程序被处理器执行时实现如本公开第一方面所述的基于柱坐标系的隧道地震波传播数值模拟方法中的步骤。

本公开第四方面提供了一种电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的程序，所述处理器执行所述程序时实现如本公开第一方面所述的基于柱坐标系的隧道地震波传播数值模拟方法中的步骤。

与现有技术相比，本公开的有益效果是：

1、本公开所述的基于柱坐标系的隧道地震波传播数值模拟方法及系统，能够解决笛卡尔坐标系下的三维数值模拟中存在的掌子面绕射波干扰严重的问题，同时针对目前柱坐标系下隧道地震波传播数值模拟方法中常采用简化的2.5D模拟方法、牺牲计算换取数值计算稳定性和“奇点”绕射干扰问题，从三维柱坐标系下弹性波波动方程出发，未作任何简化处理，同时引入了变尺度细化网格方法和“奇点”值计算方法，避免了柱坐标轴线的绕射干扰，能够以尽可能小的计算量，实现对隧道环境下地震波的传播更加真实的模拟。

2、本公开所述的基于柱坐标系的隧道地震波传播数值模拟方法及系统，以柱坐标系下的真三维弹性波波动方程的差分格式为基础，能够解决笛卡尔坐标系下的三维数值模拟中存在的掌子面绕射波干扰严重的问题，同时进行真三维的地震波数值模拟，相较于常用的简化后的2.5D数值模拟方法，模拟精度更高。

3、本公开所述的基于柱坐标系的隧道地震波传播数值模拟方法及系统，采用了变尺度细化网格剖分方式，能够有效缓解牺牲计算量换取计算稳定性的问题，能够以尽可能小的计算量实现基于柱坐标系的真三维地震波场的准确模拟。

4、本公开所述的基于柱坐标系的隧道地震波传播数值模拟方法及系统，借鉴了电磁和流体领域对奇点的处理方式，避免了极轴上的奇点产生的绕射干扰假象，进一步提高了地震波波场模拟精度，同时该方法可以直接用于柱坐标系弹性波有限差分计算中。

附图说明

构成本公开的一部分的说明书附图用来提供对本公开的进一步理解，本公开的示意性实施例及其说明用于解释本公开，并不构成对本公开的不当限定。

图1为本公开实施例1提供的径向速度分量在网格细化交界线处的插值方法。

图2为本公开实施例1提供的中心点附近纵波波场切片图。

图3为本公开实施例1提供的不同时刻纵横波波场二维切片波场快照示意图。

图4为本公开实施例1提供的不同坐标系下两条测线地震记录对比结果示意图。

具体实施方式

下面结合附图与实施例对本公开作进一步说明。

应该指出，以下详细说明都是示例性的的，旨在对本公开提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本公开所属技术领域的普通技术人员通常理解的相同含义。

需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本公开的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。

在不冲突的情况下，本公开中的实施例及实施例中的特征可以相互组合。

实施例1：

如图1所示，本公开实施例1提供了一种基于柱坐标系的隧道地震波传播数值模拟方法，主要包括以下步骤：

A：初始参数设定。

B：基于变尺度细化网格的弹性波波动方程有限差分格式的构建：

a)变尺度细化网格的剖分；

b)建立基于变尺度细化网格的弹性波波动方程有限差分格式；

c)求取柱坐标系中轴线上的“奇点值”。

C：进行波场延拓，获得模拟波场和地震记录。

一般在进行地震波正演模拟时，首先需要设定一些用于正演模拟的固定参数，比如说：模拟的地层参数信息，选择的虚拟震源函数及震源主频，以及剖分的网格数量和吸收边界条件等。这是正演模拟的一个必要条件；另外正演模拟需要采用一定的模拟方法，这里采用的是有限差分方法，需要构建有限差分格式，才能够进行迭代计算。对于步骤C，是一个具体的实施过程，即根据A的参数和B的计算方法进行计算，模拟波场在隧道条件下的传播。

进一步的，对于步骤A，所述初始参数设定包括：正演模拟所需的地震参数模型、震源主频、正演所采用的时间与空间步长及地震记录时长，CPML吸收边界条件的参数。

进一步的，对于步骤B，基于变尺度细化网格的弹性波波动方程有限差分格式的构建过程如下：

对于a)变尺度细化网格的剖分，在径向和轴向上采取常规的规则网格剖分，在周向上，针对牺牲计算量换取数值计算稳定性的问题，综合考量计算量和计算稳定性，采用如下网格剖分方法：

1)最里层初始网格数N1设置：

其中，Δr为剖分网格在r方向的长度，ceil函数表示向大取整(N1为偶数)； λ_min为波场中的最小波长；

2)判断当前θ方向的空间差分间隔是否符合频散条件：

其中，i为层数，如果上式成立，则保持当前网格剖分数。如果不满足条件，则n＝n+1并进入下一步；

3)第i层的网格数N_i：

N_i＝2ⁿN₁ (3)

其中，n为网格加密系数，n的初始值为0；

进一步的对于步骤b)，柱坐标系下的弹性波方程如下：

其中，r、θ、z为柱坐标系三个分量方向，分别对应径向、周向和轴向，λ、μ 为拉梅常数，ρ为介质密度，u_i(i＝r,θ,z)为振幅。

利用位移对时间的二阶偏导转换为质点振动速度对时间的一阶偏导数，将上式降阶为一阶速度-应力方程组，如式(5)，随后在空间和时间上均采用2阶的交错网格离散方法，得到其差分格式。

上式(5)中，r、θ、z为柱坐标系三个分量方向，λ、μ为拉梅常数，ρ为介质密度，V_i(i＝r,θ,z)为速度分量，τ_ij、σ_ij(i,j＝r,θ,z)分别为剪应力和正应力分量。

由于加密网格会造成分界线两侧的网格参数点错位半个网格间距，阻碍了差分算法的正常进行。因此，必须对分界线两侧区域的参数进行插值处理，使其转换为常规网格差分。网格加密分界线两侧的各个参数的计算思想基本一致，都是利用相邻节点上的已知量计算插值，下面以V_r的插值计算具体说明，如附图1。

切取k+1/2平面内的部分网格，根据差分公式可知，细化交界线处的Vr在径向上缺少区域A中对应的值。因此需要先利用相邻两个网格中的σ_rr值进行线性插值，获得一个虚位的

值，

此外：

从而可得细化交界线处V_r值：

公式(6)-(8)中，r、θ、z为柱坐标系三个方向分量，i,j,k为节点位置坐标，t为时间，ρ为介质密度，

为在t+1/2时刻节点在r方向上的速度分量，

分别节点在t时刻对应的剪应力和正应力分量。

对于V_θ、V_z、σ_rr、σ_θθ、σ_zz、τ_rθ、τ_zθ、τ_rz可采用相同的思想进行求解。

进一步地，对于步骤c)，奇点问题的处理方法如下：

通过参考柱坐标系下流场的有限差分模拟，提出了针对纵波波场

值的插值算法。首先，柱坐标系中单位向量的径向分量可以通过在笛卡尔坐标系下的正交单位向量表示：

e_r＝e_x cosθ+e_y sinθ. (9)

其中，e_i(i＝r,x,y)为单位向量，x,y为笛卡尔坐标系下的方向分量，r,θ为柱坐标系下的径向和周向分量。

如附图2，将柱坐标系下纵波波场中的“奇点值”

和相近网格上的v_θ建立联系，从而利用环向速度矢量来插值计算出“奇点值”

插值公式的具体推导过程如下：

通过式(10)，可以建立起1/2网格上v_x、v_y同

之间的关系，然后在笛卡尔坐标系下，分别求取x＝y＝1/2上v_x、v_y值的平均值，即可获取中心点处的

和

最后，联立式(9)和式(10)，可得

的插值近似计算公式：

其中，θ_j+1/2＝(j+1/2)Δθ。

式(9)～(11)中，上标代表时间步，下标代表空间网格坐标(分别对应x， y，z轴)。在柱坐标中，定义逆时针方向为正，x,y为笛卡尔坐标系下的方向分量，θ为柱坐标系下的周向分量，t为时间，

为纵波波场速度在对应方向上的速度分量。

对于横波波场，横波波场同电磁场一样都是一个无散场，利用柱坐标系电磁波有限差分中的处理方法来计算横波波场的“奇点值”(

和

)。得到横波波场“奇点值”的计算公式：

式(12)-(14)中，上标代表时间步，下标代表空间网格坐标(分别对应x， y，z轴)。在柱坐标中，定义逆时针方向为正，x,y为笛卡尔坐标系下的方向分量，θ为柱坐标系下的周向分量，t为时间，

为纵波波场速度在对应方向上的速度分量和剪应力分量。V_S是横波波速，N是当前层的网格数量， N须为4的倍数。

对于步骤C，利用步骤A中设定的初始参数和步骤B中的于变尺度细化网格的弹性波波动方程有限差分方法进行波场延拓，获得模拟波场和地震记录包括，采用雷克子波和CPML边界条件对波长进行延拓，获得模拟波场和地震记录

本实施例在MATLAB2016b平台上实现，以具体实例进行分析，包括以下步骤：

获取初始参数设定包括：正演模拟所需的地震参数模型、震源主频、正演所采用的时间与空间步长及地震记录时长、CPML吸收边界条件的参数；

本实例采用模型大小40m×40m×140m，网格大小为0.5m。隧道长30m，宽度(直径)为12m。

在观测方式方面，设置一个震源点，震源点位于隧道顶端中心点，距离掌子面5m，震源主频120Hz。检波器测线有两条，测线1设置在掌子面上，围绕隧道一圈；测线2设置在掌子面前方20m至50m，沿着z轴布置，间隔1m，共31 个。

在围岩介质设置方面，除隧道区域设置为空气介质外，围岩区域设置为均匀介质，纵波波速4000m/s，横波波速2300m/s。

采用变尺度细化网格方法，在径向和轴向上采取常规的规则网格剖分，在周向上，针对牺牲计算量换取计算稳定性的问题，考虑频散条件和计算量，采用如下网格剖分方法：

(1)最里层初始网格数N1设置：

其中ceil函数表示向大取整(N1为偶数)；λ_min为波场中的最小波长；

(2)判断当前θ方向的空间差分间隔是否符合频散条件：

i 为层数，如果符合，则保持当前网格剖分数。如果不符合，则n＝n+1并进入下一步；

(3)建立基于变尺度细化网格的弹性波波动方程有限差分格式；

(4)借鉴流体和电磁领域对于“奇点”值的处理方法，分别对横波和纵波奇点值进行求取；

(5)借鉴流体和电磁领域对于“奇点”值的处理方法，分别对横波和纵波奇点值进行求取；

(6)采用雷克子波和CPML边界条件对波长进行延拓，获得模拟波场和地震记录，如图3和4。

从附图3可以观察到，在笛卡尔坐标系中，掌子面出存在突变点，波场会存在绕射波干扰，在柱坐标系中，可以很好的避免这个问题。但是因为沿隧道两侧传播横波在隧道底部相遇会产生干涉，因此会产生轻微的频散现象。更进一步对比测线1和测线2的地震记录，如附图4所示，同样可以看出柱坐标在压制突变点绕射波干扰上要明显优于笛卡尔坐标系。

实施例2：

本公开实施例2提供了一种基于柱坐标系的隧道地震波传播数值模拟系统，包括：

数据获取模块，被配置为：获取待模拟隧道的基础数据；

所述系统的工作方法与实施例1提供的基于柱坐标系的隧道地震波传播数值模拟方法相同，这里不再赘述。

实施例3：

本公开实施例3提供了一种介质，其上存储有程序，该程序被处理器执行时实现如本公开实施例1所述的基于柱坐标系的隧道地震波传播数值模拟方法中的步骤，所述步骤为：

获取待模拟隧道的基础数据；

根据获取的基础数据设定初始参数；

详细步骤与实施例1提供的基于柱坐标系的隧道地震波传播数值模拟方法相同，这里不再赘述。

实施例4：

本公开实施例4提供了一种电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的程序，所述处理器执行所述程序时实现如本公开实施例1所述的基于柱坐标系的隧道地震波传播数值模拟方法中的步骤，所述步骤为：

获取待模拟隧道的基础数据；

根据获取的基础数据设定初始参数；

本领域内的技术人员应明白，本公开的实施例可提供为方法、系统、或计算机程序产品。因此，本公开可采用硬件实施例、软件实施例、或结合软件和硬件方面的实施例的形式。而且，本公开可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器和光学存储器等)上实施的计算机程序产品的形式。

本公开是参照根据本公开实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于一计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。其中，所述的存储介质可为磁碟、光盘、只读存储记忆体(Read-Only Memory， ROM)或随机存储记忆体(RandomAccessMemory，RAM)等。

以上所述仅为本公开的优选实施例而已，并不用于限制本公开，对于本领域的技术人员来说，本公开可以有各种更改和变化。凡在本公开的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本公开的保护范围之内。

Claims

1.一种基于柱坐标系的隧道地震波传播数值模拟方法，其特征在于，包括以下步骤：

获取待模拟隧道的基础数据；

根据获取的基础数据设定初始参数；

2.如权利要求1所述的基于柱坐标系的隧道地震波传播数值模拟方法，其特征在于，待模拟隧道的基础数据至少包括隧道的形状、围岩和隧道前方区域的弹性参数以及采用的观测装置形式、震源主频、采样频率和采样时长。

3.如权利要求1所述的基于柱坐标系的隧道地震波传播数值模拟方法，其特征在于，在周向上进行变尺度细化网格的剖分，具体为：

4.如权利要求1所述的基于柱坐标系的隧道地震波传播数值模拟方法，其特征在于，构建柱坐标系下的初始弹性波方程，利用位移对时间的二阶偏导转换为质点振动速度对时间的一阶偏导数，将初始弹性波方程降阶为一阶速度应力方程组，随后在空间和时间上均采用二阶的交错网格离散方法，得到差分格式。

5.如权利要求1所述的基于柱坐标系的隧道地震波传播数值模拟方法，其特征在于，对网格加密分界线两侧区域的参数进行插值处理，转换为常规网格差分。

6.如权利要求5所述的基于柱坐标系的隧道地震波传播数值模拟方法，其特征在于，对第一参数进行差值处理，具体为：

7.如权利要求1所述的基于柱坐标系的隧道地震波传播数值模拟方法，其特征在于，将柱坐标系下纵波波场中的奇点值和相近网格上的纵波波场在周向的速度分量建立联系，利用环向速度矢量来插值计算出奇点值；

或者，

利用柱坐标系电磁波有限差分中的奇点处理方法计算横波波场的奇点值。

8.一种基于柱坐标系的隧道地震波传播数值模拟系统，其特征在于，包括：

数据获取模块，被配置为：获取待模拟隧道的基础数据；

9.一种介质，其上存储有程序，其特征在于，该程序被处理器执行时实现如权利要求1-7任一项所述的基于柱坐标系的隧道地震波传播数值模拟方法中的步骤。

10.一种电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的程序，其特征在于，所述处理器执行所述程序时实现如权利要求1-7任一项所述的基于柱坐标系的隧道地震波传播数值模拟方法中的步骤。