[发明专利]一种“M(2，4)FDTD+FDTD”的低色散低频地波传播时延预测方法

摘要

本发明公开了一种“M(2，4)FDTD+FDTD”的低色散低频地波传播时延预测方法，该方法通过对整个计算区域进行分层处理，将粗网格M(2，4)FDTD算法(用于上层区域)与细网格传统FDTD算法(用于下层区域)相结合，利用修正因子以及亚网格技术，进行低频地波传播时延的高精度快速预测，本发明解决了传统全区域FDTD算法在减小数值色散误差与降低计算机消耗之间的矛盾，提高了低频地波传播时延预测精度的同时，减少了计算机内存资源占用，提高了计算速度。

主权项

1.一种M(2，4)FDTD+FDTD的低色散低频地波传播时延预测方法，其特征在于，通过对低色散低频地波计算区域进行分层处理，将粗网格M(2，4)FDTD与细网格传统FDTD方法相结合，进行低频地波传播时延的预测，从而达到在保证预测精度的同时，提高计算速度、减少计算机内存占用，具体步骤如下：步骤1、设置并生成模型文件、输入模型文件；步骤2、参数设置及初始化；步骤3、添加场源；步骤4、更新整个计算区域的ρ向电场分量Eρ；步骤5、更新整个计算区域的z向电场分量Ez；步骤6、更新整个计算区域的φ向磁场分量Hφ；步骤7、更新激励源；步骤8、判断结束条件，循环；步骤9、观测点观测量计算并输出；步骤1具体为：步骤(1.1)、设模型文件的上层区域大小Nρ1×Nz1，下层区域大小Nρ2×Nz2，CFS‑PML层数为NPML，其中Nρ为ρ方向网格数，Nz为z方向网格数，标号1和2分别表示上层和下层；步骤(1.2)、设置空间、时间步长：上层粗网格空间步长为δ1，其中，Δρ1＝Δz1＝δ1，下层细网格空间步长为δ2，其中，Δρ2＝Δz2＝δ2，其中Δρ为ρ方向网格大小，Δz为z方向网格大小，标号1和2分别表示上层和下层，时间步长为Δt，粗网格M(2，4)FDTD的时间步长设置为与细网格FDTD同样的时间步长；步骤(1.3)、设置迭代次数为NTCalc；步骤(1.4)、设置传播路径电参数：ρ方向起始网格位置为ρStart、ρ方向结束网格位置为ρEen、z方向起始网格位置为zStart、z方向结束网格位置为zEen、大地电导率为σ、相对介电常数为εr；步骤(1.5)、设置吸收边界：CFS‑PML层数为NPML，相关参数为κηmax，αηmax，σηmax，其中η＝ρ,z；步骤(1.6)、设置场源：源的个数NSource、位置[ρStart，ρEen]和[zStart，zEen]，源的种类：有两种激励源供选择：单频正弦源、Loran_C源，激励场形式：有三种激励形式：Ez、Eρ、Hφ，源的类型：软源或硬源，幅度A，频率f，单频正弦源中的常数t0，高斯脉宽Tp，时延/包周差τ；步骤(1.7)、设置观测点：观测点个数NVPoint，位置[ρStart，ρEen]和[zStart，zEen]，输出场量类型Ez、Eρ或Hφ；步骤4具体为：步骤(4.1)、首先根据所述步骤1中模型文件中定义的上层区域采用粗网格M(2，4)FDTD方法，对该区域内的电磁场分量Eρ1进行更新，其中，上层区域不包含PML层，Eρ1中的下标1表示上层区域，具体更新公式为：式中，n表示时间步，标号i1和k1分别表示上层区域中ρ向和z向的空间位置，ε₀为真空中介电常数，k₁为环路系数；步骤(4.2)、根据所述步骤1中模型文件中定义的下层区域采用传统FDTD方法，对该区域内的电磁场分量Eρ2进行更新，其中，下层区域不包含PML层，Eρ2中的下标2表示下层区域，具体更新公式为：式中，n表示时间步，标号i2和k2分别表示下层区域中ρ向和z向的空间位置；步骤(4.3)、对边界区域的场量传递进行更新：边界区域场量的传递分为如下几种情况：a、当上层粗网格中Eρ1的计算需用到下层细网格中的Hφ2时，由于上下粗细网格比为奇数，场量重叠，故直接取相应细网格的场量即可；b、当下层细网格中边界网格处Eρ2的计算需用到边界线上的Hφ时，直接采用边界线上细网格的Hφ2即可；步骤(4.4)、对CFS‑PML中的电场分量进行更新：上方吸收边界中的网格空间大小与粗网格相同，右侧吸收边界中网格的空间大小与相邻左侧计算区域的网格空间大小相同，也分为上下两层，上层与粗网格一致，下层与细网格一致，吸收边界中的场分量的差分公式为：其中，其中，n表示时间步，i和k分别表示计算区域中ρ向和z向的空间位置，kz、σz和αz为吸收边界参数；所述步骤5具体为：步骤(5.1)、先根据所述步骤1中模型文件中定义的上层区域采用粗网格M(2，4)FDTD方法，对该区域内的电磁场分量Ez1进行更新，其中，上层区域不包含PML层，Ez1中的下标1表示上层区域，具体更新公式为：式中，n表示时间步，标号i1和k1分别表示上层区域中ρ向和z向的空间位置，ε₀为真空中介电常数，k₁为环路系数；步骤(5.2)、根据所述步骤1中模型文件中定义的下层区域采用传统FDTD方法，对该区域内的电磁场分量Ez2进行更新，其中，下层区域不包含PML层，Ez2中下标2表示下层区域，具体更新公式为：式中，n表示时间步，标号i2和k2分别表示下层区域中ρ向和z向的空间位置；步骤(5.3)、交界线上Ez均为细网格的Ez2，上层粗网格中Ez1的更新不包括交界线上的点；步骤(5.4)、对CFS‑PML中的电场分量进行更新：吸收边界中的场分量的差分公式为：式中，标号m＝(i,k+1/2)，且其中，n表示时间步，i和k分别表示计算区域中ρ向和z向的空间位置，αρ、kρ、σρ、kρmax和σρmax为吸收边界参数，ρ＝i×δ表示ρ向长度，ρ0是PML层与非PML层的分界位置；所述步骤6具体为：步骤(6.1)、首先根据所述步骤1中模型文件中定义的上层区域采用粗网格M(2，4)FDTD方法，对该区域内的电磁场分量Hφ1进行更新，其中，上层区域不包含PML层，具体更新公式为：其中，n表示时间步，i1和k1分别表示上层区域中ρ向和z向的空间位置，k1和k2为环路系数，c为光速，μr为磁导率；步骤(6.2)、根据所述步骤1中模型文件中定义的下层区域采用传统FDTD方法，对该区域内的φ向磁场分量Hφ2进行更新，不包括PML层，具体更新公式为：其中，n表示时间步，i2和k2分别表示下层区域中ρ向和z向的空间位置，c为光速，μr为磁导率；步骤(6.3)、边界区域的场量传递：边界区域场量的传递分为如下2种情况：a、当上层粗网格中Hφ1的计算需用到下层细网格中的Eρ2时，由于上下粗细网格比为奇数，场量重叠，故直接取相应细网格的场量即可，同时粗网格中Hφ1的更新不包括交界线上的点；b、当交界线上的Hφ均取细网格的Hφ2，且交界线上Hφ2的计算需用到上层细网格中的Eρ1时，采用线性插值的方法获取，具体过程如下：的差分公式为：式中，根据粗网格中的电场分量和采用线性插值的方法得到即步骤(6.4)、对CFS‑PML中的电场分量Hφ进行更新：吸收边界中的场分量Hφ的差分公式为：式中，其中，n表示时间步，i和k分别表示计算区域中ρ向和z向的空间位置，αη、kη和ση为吸收边界参数。