一阶声波方程建模的混合显式隐式交错网格有限差分格式

www.nature.com/scientificreports一种混合显式隐式交错网格有限差分方案对于一阶声波方程建模打开梁文全1、王艳飞2、曹静杰3 & Ursula Iturrarán‑Viveros4*隐式交错网格有限差分 (SGFD) 方法广泛用于一阶声学波动方程模型。相同的隐式 SGFD 运算符通常用于所有第一个一阶声波方程中的阶空间导数。在本文中，我们提出了一种混合显式隐式 SGFD (HEI‑SGFD) 方案，可以同时保留波动方程仿真精度，提高波动方程仿真速度。我们使用二阶显式一阶声波方程中一半一阶空间导数的 SGFD 算子。同时，我们使用隐式 SGFD 算子，在对角线方向添加点一阶声波方程中的其他一阶空间导数。拟议的 HEI-与隐式 SGFD 相比，SGFD 方案的波动方程模拟速度几乎翻倍计划。本质上，所提出的 HEI-SGFD 方案等价于二阶 FD 方案使用普通的网格格式。然后，我们通过以下方式确定时空域中的 HEI-SGFD 系数使用最小二乘法最小化相速度误差。最后，该方案的有效性通过色散分析和数值模拟证明了所提出的方法。地震波方程在地震勘探、震源定位等方面有着广泛的应用。其他1-4。有限差分时域（FDTD）方法是最流行的求解方法之一地震波方程 5-7。交错网格有限差分 (SGFD) 方法是 FDTD 方法的一部分，因其精度高、内存要求低而常用于地震波外推并且易于实施8-11。为了提高方法的准确性和效率，在对角线上添加点的模板方向被采纳12-14。与之前的高阶 FD 钢网相比，使用这些新的 FD 钢网可以提高波动方程的模拟效率，同时通过使用更大的时间来保持高精度步。显式和隐式 SGFD 方案都广泛用于波动方程建模。隐含的与显式 SGFD 方案相比，SGFD 方案可以用更短的 FD 算子实现更高的精度15,16。提高 FDTD 方法的准确性和效率的另一种方法是利用优化的 FD 算子系数 17-21。地球物理成像需要大量的波动方程建模22。通常需要很长时间才能进行地球物理成像。这些计算必须部署在超级计算机上或依赖 GPU（图形处理单元）以便在可容忍的等待时间 23,24 内获得结果。使用有效波方程建模算法可以减少运行时间，从而节省能源并潜在地减少碳排放。相同的 SGFD 算子是用于所有一阶空间导数的最常用算子一阶声波方程。我们建议对不同的空间导数使用不同的 SGFD 算子在声波方程 25,26 中的 tives 并在本文中提出了 HEI-SGFD 方案。我们使用二阶一阶声波方程中一些一阶空间导数的显式SGFD格式。同时我们使用高阶隐式SGFD方案，其余的在对角线方向上添加点1龙岩大学资源工程学院,龙岩 364000 2钥匙中国科学院地质与地球物理研究所石油资源研究室，中华人民共和国北京 100029。 3智能检测与装备重点实验室京津冀城市群地下空间，自然资源部，河北省地质勘探局石家庄大学, 050031, 中华人民共和国4国立大学数学系墨西哥自治市，墨西哥城，墨西哥。邮箱：yfwang@mail.iggcas.ac.cn*科学报告|(2022) 12:10967|https://doi.org/10.1038/s41598-022-15112-x1 www.nature.com/scientificreports/的一阶空间导数。提出的 FD 方案称为：HEI-SGFD 方案。 TIS HEI-SGFD与隐式 SGFD 方案相比，该方案可以节省约 45% 的计算时间，同时它仍然保持高精度。理论Te 一阶速度-应力声波方程可以描述如下：∂吨 = −ρv2∂五∂五 zX(1)(2)(3)+,∂z∂ X∂五X∂吨1 ∂p,ρ ∂ X= −= −∂五z∂吨1 ∂p,ρ ∂z在哪里p是压力，ρ是密度和v是波的传播速度，五和五是粒子速度。Xz本文中的方程假设密度恒定。Te 声波方程 (1) 到 (3) 可以离散化为ķ+1/2z(一世+1/2,j)ķ+1/2,ķķ−1一世,j2− v τ � 五(4)(5)(6)p = p+ �X五X(一世,j+1/2)一世,jzķ+1/2X(一世,j+1/2)ķ−1/2X(一世,j+1/2)− τ pķ ,五= 五X 一世,jķ+1/2z(一世+1/2,j)ķ−1/2z(一世+1/2,j)ķ 五= 五− τ � p,z 一世,j时间步长在哪里，ςn = ς(z + 小时, X + jh,吨 + nτ), ς = [五 ,五 , p]. Te 符号 Δ 表示隐含空间导数的离散形式，例如我们在以下等式中使用它：米,jzXδ2 � pķ米X 一世,j 1 2+ /11+ Cķp一世,j+米− pķ + pķ− pķ ,一世+1,j 一世−1,j+1 一世−1,j1+bh2=C米− pķpķ一世,j−米+1米+1 一世+1,j+1δX2HH米=1(7)(8)在哪里2δ qδX2q(X + H) + q(X − H) − 2q(X).=H2为了使用更大的时间和空间网格间隔，在方程。 (8),米是 SGFD 算子的长度，C米和b是 SGFD 系数和H是网格大小，米 = 3 和 C= 0本文采用。什么时候b和米+1C等于零，方程式。 (4) 到 (8) 是传统的显式高阶 SGFD 方案。米+1十我们可以隐式地得到空间导数的近似值q ≈ ∂p∂ X 作为 米1ķ+(1−2b)q我知道,j+1/2+bqķķķ) + C米+1(pķķ+ pķ− pķ ) .bq=C (p− p− p一世,j+3/2一世,j−1/2米一世,j+米一世,j−米+1一世+1,j+1一世+1,j一世−1,j+1一世−1,jH米=1(9)此后，我们表示方程式。 (4) 到 (6) 作为隐式 SGFD 方案。可以使用更大的时间步长系数C.同时，可以使用较短的算子长度和系数b27-30。我们建议米+1对方程中的空间导数使用最简单的显式二阶 SGFD 算子。 (2) 和 (3)。提议的HEI-SGFD 方案由下式给出：ķ+1/2z(一世+1/2,j)ķ+1/2X(一世,j+1/2),ķ一世,jķ−1一世,j2− v τ � 五(10)(11)(12)p = p+ �X五 ķzτķ+1/2X(一世,j+1/2)ķ−1/2X(一世,j+1/2)− p我知道,j,五= 五−p一世,j+1Hτķ+1/2z(一世+1/2,j)ķ−1/2z(一世+1/2,j) ķ− p我知道,j.五= 五−p一世+1,jH方程式中的 Te SGFD 方案。 （10）是隐含的，而方程式中的 SGFD 方案。 (11) 和 (12) 是明确的并且二阶。因此，我们将提出的 SGFD 方案称为 HEI-SGFD 方案。 Te HEI-SGFD 方案是一个隐式 SGFD 方案，尽管一些空间导数是用显式 SGFD 方案近似的。似乎方程式。 (11) 和 (12) 不准确，但事实并非如此。方程中的 SGFD 系数。 (10)通过考虑方程式确定。 (11) 和 (12)。很容易观察到，所提出的 HEI-SGFD 方案在科学报告|(2022) 12:10967 |https://doi.org/10.1038/s41598-022-15112-x2 www.nature.c om/scientificreports/方程。与方程式中给出的隐式 SGFD 方案相比，（10）到（12）需要更少的浮点计算。(4)-(6)。方程式所需的计算资源。 (11) 和 (12) 小于等式。 (5) 和 (6)。替换方程。 (11) 和 (12) 进入等式。 (10) 我们得到：2v τ22 2v τķ+1一世,jķķ−1一世,jķ′′(一世,j)zķ′′(一世,j)X(13)p=2p − p+p+p,一世,jH2H2在哪里δ2p′′(一世,j)z米−11 + bh2= −2C p +ķ(C米 − C米+1)(pķ+ p我知道−米,j) + C米(pķ+ pķ一世+米,j 一世−米,j))1 一世,j一世+米,jδz2(14)米=1+C米+1(pķ+ pķ+ pķ + p我知道−1,j−1),一世+1,j−1一世+1,j+1一世−1,j+1δ2p′′(一世,j)X米−11 + bh2= −2C p +ķ(C米 − C米+1)(pķ+ p我知道,j−米) + C米(pķ + pķ一世,j+米 一世,j−米1 一世,j一世,j+米δX2(15)米=1+C米+1(pķ+ pķ+ pķ+ p我知道−1,j−1).一世−1,j+1一世+1,j+1一世+1,j−1理论上，方程式中的一阶 HEI-SGFD 方案。 (10) 到 (12) 等价于二阶隐式具有普通网格格式的 FD 方案，如方程式所示。 (13)。所提出的 HEI-SGFD 方案的色散关系由下式给出米−2sin(0.5kXH）2c 罪（（米 − 0.5)k H） + 4c米余弦(k h)(cos(k H） − 1z XX米1米=1米[1 + 2b(cos(k H） − 1 + −2sin(0.5k H）z2c 罪（（米 − 0.5)k H） + 4c米余弦(k h)(cos(k H） − 1X zzz米1米=11[1 + 2b(cos(kXH） − 1 = [1 + 2b(cos(k H） − 1+ 2b(cos(k H） − 1余弦(kvτ)− 2].Xzr 2(16)方程中的 SGFD 系数。 (16) 应仔细确定优化方法。相似的对于 Wang 等人 31 提出的方法，我们从方程得到目标函数。 (16) 通过最小化误差真实速度和相速度之间：2r2G2[1+2b(cos(k H）−1+2b(cos(k H）−1ķπ/4ķπ/ 阿科斯 1 +v2 fdzX（三） =− 1=− 1,v千伏ķ=3e- 12θ=0ķ= 3e- 12θ=0(17)在哪里米G = −2sin(0.5kXH）2c 罪（（米 − 0.5)k H） + 4c米余弦(k h)(cos(k H） − 1 [1 + 2b(cos(k H） − 1zX米+1Xz米=1米+ −2sin(0.5k H）z2c 罪（（米 − 0.5)k H） + 4c米余弦(k h)(cos(k H） − 1 [1 + 2b(cos(k H） − 1Xz米+1zX米=1(18)方程中只有未知数。 (17) 是C (米 = 1,· · · ,米 + 1） 和b.波数ķ在等式。 (17) 应该开始米从零开始。然而，ķ出现在分母中，导致不稳定。因此，我们让ķ从一个非常小数，例如，3.0× 10−.我们假设波的传播速度、时间等参数step 和空间网格间隔已经给出（然后r =H会知道的）。其他两个参数是ķ和θ.传播角θ是从 0 到π/4π/16π/16π/4)。从方程式。 (17)，我们发现在频波数域 SGFD 系数与 Courant 比有关r.当一个人考虑一个具有固定时间和空间步长网格间隔的波动方程模拟，SGFD 系数不同对于不同的速度，模板形成一个大的 3D 矩阵，用于 2D 复杂的速度模型。为了得到我们应用 MATLAB 函数的混合显式-隐式 SGFD 系数lsqnonlin解决非线性方程中的最小二乘问题。 (17)。提出的 HEI-SGFD 方案可以很容易地扩展到 3D：ķ+1/2z(l+1/2,米,n)ķ+1/2X(l,米+1/2,n)ķ+1/2是的(l,米,n+1/2),ķpl,米,nķ−1l,米,n2− v τ � 五(19)(20)= p+ �X五+ �是的五zτķ+1/2X(l,米+1/2,n)ķ−1/2X(l,米+1/2,n) ķ− pķ ,l,米,n五= 五−pl,米+1,nH科学报告|(2022) 12:10967 |https://doi.org/10.1038/s41598-022-15112-x3 www.nature.c om/scientificre