跳至主要內容

SG-PML吸收边界条件（边界反射问题）

原创Xenny约 3322 字大约 14 分钟FWIPMLFWI

SG-PML吸收边界条件（边界反射问题）

  • 本文将对博文 波动方程的有限差分解（地震数据正演） 中提到的边界反射问题进行研究。解决地震正演中的边界反射带来的干扰。

    首先我们要知道波为什么会反射，当从一种介质进入另一种介质时，在介质边界由于阻抗的不连续性，导致一部分能量继续进入新的介质称为透射波，一部分能量返回原介质称为反射波。而我们将需要了解反射部分具体是多少，以及如何消除边界反射问题。

吸收边界条件

  • 自然界中地震波传播到介质边界（中断层）时，界面的波阻抗出现差异导致波发生反射和透射，地球介质可以看作“无限大”的区域，但在模拟计算中受到计算量的限制，模拟的传播区域大小是由界限的，这个人为截取的区域便是人工边界。

    在人工边界处，边界外的波速可以看为0，即边界反射系数绝对值为1，此时表现在数值模拟中便产生了全反射现象，对地震正演中的波场计算将产生影响。此时我们需要引入吸收边界的概念来吸收掉这些传播至边界的波。

  • 吸收边界条件即人为的在模拟计算区域的边界处设置一定厚度的吸收层，让能量在吸收层呈指数曲线衰减，最终在吸收层边界处近似衰减到0而避免边界反射。

    注意这里我们的衰减函数也是有讲究的，因为如果衰减不够连续的话吸收层也将看成断层产生反射，虽然这里使用指数曲线进行衰减，但是实际模拟时还是会产生一定量的反射波。

PML

  • PML边界条件将波场分裂为垂直和平行的PML边界两部分，通过再分裂向量满足的微分方程中添加阻尼项得到PML中的微分方程。同时renger证明了如果不进行离散化PML条件可以完全无反射地吸收所有入射波。当然我们的模拟实现是没办法做到这一点的。

一维声波方程

  • 这里我们先以一维声波为例说明PML边界条件思想，一维声波在均匀介质中的波动方程为[1]

    1v2(x)∂p(x,t)∂t2=∂2p(x,t)∂x2(1) \frac{1}{v^2(x)}\frac{\partial p(x,t)}{\partial t^2} = \frac{\partial^2 p(x,t)}{\partial x^2}\tag{1}

    其中vv代表速度，pp代表波信息（即位移），我们添加一个中间变量A(x,t)A(x,t)将上述方程改写为两个一阶偏微分方程

    ∂p(x,t)∂t=v2(x)∂A(x,t)∂x∂A(x,t)∂x=∂p(x,t)∂x(2) \begin{aligned} \frac{\partial p(x,t)}{\partial t} = v^2(x)\frac{\partial A(x, t)}{\partial x}\\ \frac{\partial A(x,t)}{\partial x} = \frac{\partial p(x,t)}{\partial x} \end{aligned}\tag{2}

    显然式2再对tt求偏导即得式1，此时我们可以得到一个PML控制方程组

    ∂p∗(x,t)∂t+d(x)p∗(x,t)=v2(x)∂A∗(x,t)∂x∂A∗(x,t)∂t+d(x)A∗(x,t)=∂p∗(x,t)∂x(3) \begin{aligned} \frac{\partial p^{*}(x,t)}{\partial t} + d(x)p^{*}(x,t) = v^2(x)\frac{\partial A^{*}(x,t)}{\partial x}\\ \frac{\partial A^{*}(x,t)}{\partial t} + d(x)A^{*}(x,t) = \frac{\partial p^{*}(x,t)}{\partial x} \end{aligned}\tag{3}

    其中p∗p^{*}代表式3方程的解，d(x)d(x)代表衰减系数，我们后续步骤的目的便是研究这个PML控制方程的解p∗p^{*}和原解pp的关系。对式3做关于时间的傅里叶变换得

    iωp^∗(x,ω)+d(x)p^∗(x,ω)=v2(x)∂A^∗(x,ω)∂xiωA^∗(x,ω)+d(x)A^∗(x,ω)=∂p^∗(x,ω)∂x(4) \begin{aligned} \mathrm{i}\omega\hat{p}^{*}(x, \omega) + d(x)\hat{p}^{*}(x,\omega) = v^2(x)\frac{\partial \hat{A}^{*}(x,\omega)}{\partial x}\\ \mathrm{i}\omega\hat{A}^{*}(x,\omega) + d(x)\hat{A}^{*}(x,\omega) = \frac{\partial \hat{p}^{*}(x,\omega)}{\partial x} \end{aligned}\tag{4}

    其中p^∗(x,ω), A^∗(x,ω)\hat{p}^{*}(x,\omega),\ \hat{A}^{*}(x,\omega)分别代表p∗(x,ω), A∗(x,ω)p^{*}(x,\omega),\ A^{*}(x,\omega)的傅里叶变换。再将实数空间域转换到复数空间域后的坐标得

    x~(x)=x−iω∫0xd(s)ds(5) \widetilde{x}(x) = x - \frac{\mathrm{i}}{\omega}\int_{0}^x{d(s)}\mathrm{d}s\tag{5}

    则有

    x~(x)∂x=1−id(x)ω∂A^∗(x)∂x=∂A^∗(x)∂x~x~(x)∂x=(1−id(x)ω)∂A^∗(x~)∂x~∂p^∗(x)∂x=∂p^∗(x)∂x~x~(x)∂x=(1−id(x)ω)∂p^∗(x~)∂x~(6) \begin{aligned} \frac{\widetilde{x}(x)}{\partial x} &= 1 - \frac{\mathrm{i}d(x)}{\omega}\\ \frac{\partial \hat{A}^{*}(x)}{\partial x} = \frac{\partial \hat{A}^{*}(x)}{\partial \widetilde{x}}\frac{\widetilde{x}(x)}{\partial x} &= (1 - \frac{\mathrm{i}d(x)}{\omega})\frac{\partial \hat{A}^{*}(\widetilde{x})}{\partial \widetilde{x}}\\ \frac{\partial \hat{p}^{*}(x)}{\partial x} = \frac{\partial \hat{p}^{*}(x)}{\partial \widetilde{x}}\frac{\widetilde{x}(x)}{\partial x} &= (1 - \frac{\mathrm{i}d(x)}{\omega})\frac{\partial \hat{p}^{*}(\widetilde{x})}{\partial \widetilde{x}} \end{aligned}\tag{6}

    将式6带入式4得

    iωp^∗(x~,ω)=v2(x~)∂A^∗(x~)∂x~iωA^∗(x~,ω)=∂p^∗(x~)∂x~(7) \begin{aligned} \mathrm{i}\omega \hat{p}^{*}(\widetilde{x}, \omega) = v^2(\widetilde{x})\frac{\partial \hat{A}^{*}(\widetilde{x})}{\partial \widetilde{x}}\\ \mathrm{i}\omega \hat{A}^{*}(\widetilde{x}, \omega) = \frac{\partial \hat{p}^{*}(\widetilde{x})}{\partial \widetilde{x}} \end{aligned}\tag{7}

    再对此式子做傅里叶逆变换到时间域上得

    ∂p∗(x~,t)∂t=v2(x~)∂A∗(x~,t)∂x~∂A∗(x~,t)∂t=∂p∗(x~,t)∂x~(8) \begin{aligned} \frac{\partial p^{*}(\widetilde{x}, t)}{\partial t} = v^2(\widetilde{x})\frac{\partial A^{*}(\widetilde{x}, t)}{\partial \widetilde{x}}\\ \frac{\partial A^{*}(\widetilde{x}, t)}{\partial t} = \frac{\partial p^{*}(\widetilde{x}, t)}{\partial \widetilde{x}} \end{aligned}\tag{8}

    可以发现式8与式2形式完全相同，设式2的解为p(x,t)p(x,t)，式8的解为p(x~,t)p(\widetilde{x},t)，则式3的解p∗p^{*}为

    p(x−iω∫pxd(s)ds,t) p(x - \frac{\mathrm{i}}{\omega}\int_p^x{d(s)}\mathrm{d}s, t)

  • 这代表的含义为我们可以先求式2的解（有限差分或者其他方法），再对解做坐标变换便可以得到式3的解。同时对于式2在均匀介质中有如下特解

    u(x,t)=p0exp⁡[−i(kxx−ωt)](9) u(x,t) = p_0 \exp[-\mathrm{i}(k_xx-\omega t)]\tag{9}

    则对应的式3中的解为

    u∗(x,t)=p0exp⁡[−i(kxx~−ωt)](10) u^{*}(x,t) = p_0 \exp[-\mathrm{i}(k_x\widetilde{x}-\omega t)]\tag{10}

    将式5带入得

    u∗(x,t)=p0exp⁡[−i(kxx−ωt)−kxω∫0xd(s)ds](11) u^{*}(x,t) = p_0\exp[-\mathrm{i}(k_xx - \omega t)-\frac{k_x}{\omega}\int_0^x{d(s)}\mathrm{d}s]\tag{11}

    可得上述解得振幅比（阻尼系数）为

    ∥p∗(x,t)∥∥p(x,t)∥=exp⁡[−kxω∫0xd(s)ds](12) \frac{\lVert p^{*}(x, t)\rVert}{\lVert p(x,t)\rVert} = \exp\left[-\frac{k_x}{\omega}\int_0^x{d(s)}\mathrm{d}s \right]\tag{12}

    这代表引入d(x)d(x)后波信号将按传播距离指数衰减。

交错网格

  • 在求解二维声波方程的PML控制方程前，我们先来说明一下交错网格[2]的概念，在求解高维时，我们会将波场分裂为各个维度的分量，例如二维中的xx方向和zz方向，将二阶微分方程变成两个一阶速度-应力方程进行求解。

    对于二维声波方程，设px,pzp^x,p^z分别为x,zx,z方向的波场位移分量，则对xx和zz方向的中心差分方程为

    {pi,jx=1Δx∑n=1NCn(N)(pi+2n−12,jx−pi−2n−12,jx)pi,jz=1Δz∑n=1NCn(N)(pi,j+2n−12z−pi,j−2n−12z)(13) \begin{cases} p_{i,j}^x = \frac{1}{\Delta x}\sum_{n=1}^N C_n^{(N)}(p^x_{i+\frac{2n-1}{2}, j} - p^x_{i-\frac{2n-1}{2},j})\\ p_{i,j}^z = \frac{1}{\Delta z}\sum_{n=1}^N C_n^{(N)}(p^z_{i, j+\frac{2n-1}{2}} - p^z_{i,j-\frac{2n-1}{2}})\tag{13} \end{cases}

    其中Cn(N)C_n^{(N)}为差分系数（通过泰勒展开推导），2N2N为差分空间阶数。

  • 由式3有

    (∂∂t+dx)A=∂p∂x(14) (\frac{\partial}{\partial t}+d_x)A = \frac{\partial p}{\partial x}\tag{14}

    设右式经过空间差分后近似结果为■k\blacksquare_k，其中kk代表第kΔtk\Delta t时刻。

    则在交错网格下的PML控制方程近似方案如下

    (∂∂t+dx)Ak=Ak+1/2−Ak−1/2Δt+dxAk−1/2(15) (\frac{\partial}{\partial t}+d_x)A_k = \frac{A_{k+1/2}-A_{k-1/2}}{\Delta t} + d_xA_{k-1/2}\tag{15}

    带入式14得

    Ak+1/2=(1−Δt⋅dx)Ak−1/2+Δt⋅■∣k(16) A_{k+1/2} = (1-\Delta t\cdot d_x)A_{k-1/2} + \Delta t\cdot \blacksquare|_k\tag{16}

二维声波方程

  • 那对于二维也是同理，对于二维声波方程

    1v2(x,z)∂2p(x,z,t)∂t2=∂2p(x,z,t)∂x2+∂2p(x,z,t)∂z2(17) \frac{1}{v^2(x, z)} \frac{\partial^2 p(x, z, t)}{\partial t^2} = \frac{\partial^2 p(x,z,t)}{\partial x^2} + \frac{\partial^2 p(x,z,t)}{\partial z^2}\tag{17}

    将pp分解为两个方向得分量p=px+pzp = p_x + p_z，转为一阶速度——应力方程为

    ∂px∂t=−ρv2(x,z)∂Ax∂x∂pz∂t=−ρv2(x,z)∂Az∂z∂Ax∂t=−1ρ∂p∂x∂Az∂t=−1ρ∂p∂z(18) \begin{aligned} \frac{\partial p_x}{\partial t} &= -\rho v^2(x,z)\frac{\partial A_x}{\partial x}\\ \frac{\partial p_z}{\partial t} &= -\rho v^2(x,z)\frac{\partial A_z}{\partial z}\\ \frac{\partial A_x}{\partial t} &= -\frac{1}{\rho}\frac{\partial p}{\partial x}\\ \frac{\partial A_z}{\partial t} &= -\frac{1}{\rho}\frac{\partial p}{\partial z} \end{aligned}\tag{18}

    对应的PML控制方程为

    ∂px∂t+dxpx=−ρv2∂Ax∂x∂pz∂t+dzpz=−ρv2∂Az∂z∂Ax∂t+dxAx=−1ρ∂p∂x∂Az∂t+dzAz=−1ρ∂p∂z(19) \begin{aligned} \frac{\partial p_x}{\partial t} + d_xp_x &= -\rho v^2\frac{\partial A_x}{\partial x}\\ \frac{\partial p_z}{\partial t} + d_zp_z &= -\rho v^2\frac{\partial A_z}{\partial z}\\ \frac{\partial A_x}{\partial t} + d_xA_x &= -\frac{1}{\rho}\frac{\partial p}{\partial x}\\ \frac{\partial A_z}{\partial t} + d_zA_z &= -\frac{1}{\rho}\frac{\partial p}{\partial z} \end{aligned}\tag{19}

  • 空间上按照式13，时间上按照式16进行差分近似，得到交错网格有限差分声波方程时间递推方程为

    {Ax∣i+1/2,jk=(1−Δt⋅dx)Ax∣i+1/2,jk−1−ΔtρΔx∑n=1NCn(N)[p∣i+1/2+2n−12,jk−1/2−p∣i+1/2−2n−12,jk−1/2]Az∣i+1/2,jk=(1−Δt⋅dz)Az∣i,j+1/2k−1−ΔtρΔz∑n=1NCn(N)[p∣i,j+1/2+2n−12k−1/2−p∣i,j+1/2−2n−12k−1/2]px∣i,jk+1/2=(1−Δt⋅dx)px∣i,jk−1/2−ρv2ΔtΔz∑n=1NCn(N)[Az∣i,j+2n−12k−Az∣i,j−2n−12k]pz∣i,jk+1/2=(1−Δt⋅dz)pz∣i,jk−1/2−ρv2ΔtΔz∑n=1NCn(N)[Az∣i,j+2n−12k−Az∣i,j−2n−12k](20) \begin{cases} A_x\vert_{i+1/2, j}^k = (1-\Delta t\cdot d_x)A_x\vert_{i+1/2,j}^{k-1} - \frac{\Delta t}{\rho\Delta x}\sum_{n=1}^N{C_n^{(N)}[p\vert_{i+1/2+\frac{2n-1}{2},j}^{k-1/2}-p\vert_{i+1/2-\frac{2n-1}{2},j}^{k-1/2}]}\\ A_z\vert_{i+1/2, j}^k = (1-\Delta t\cdot d_z)A_z\vert_{i,j+1/2}^{k-1} - \frac{\Delta t}{\rho\Delta z}\sum_{n=1}^N{C_n^{(N)}[p\vert_{i,j+1/2+\frac{2n-1}{2}}^{k-1/2}-p\vert_{i,j+1/2-\frac{2n-1}{2}}^{k-1/2}]}\\ p_x\vert_{i, j}^{k+1/2} = (1-\Delta t\cdot d_x)p_x\vert_{i,j}^{k-1/2} - \frac{\rho v^2\Delta t}{\Delta z}\sum_{n=1}^N{C_n^{(N)}[A_z\vert_{i,j+\frac{2n-1}{2}}^{k}-A_z\vert_{i,j-\frac{2n-1}{2}}^{k}]}\\ p_z\vert_{i, j}^{k+1/2} = (1-\Delta t\cdot d_z)p_z\vert_{i,j}^{k-1/2} - \frac{\rho v^2\Delta t}{\Delta z}\sum_{n=1}^N{C_n^{(N)}[A_z\vert_{i,j+\frac{2n-1}{2}}^{k}-A_z\vert_{i,j-\frac{2n-1}{2}}^{k}]}\tag{20} \end{cases}

衰减系数

  • 在PML中衰减系数dd满足d>0d>0，且在模拟时我们可以将PML区域分为三种不同情况

    图1. PML吸收边界
    图1. PML吸收边界

    对于图1中区域2满足dx(x)=0,dz(z)≠0d_x(x)=0,d_z(z)\not=0，即垂直衰减水平方向不衰减，对于区域3则满足dx(x)≠0,dz(z)=0d_x(x)\not=0,d_z(z)=0。而对于区域1则满足dx(x)≠0,dz(z)≠0d_x(x)\not =0,d_z(z)\not= 0。

  • 一般选择[3]

    d(x)=d0(xδ)n(22) d(x) = d_0\left(\frac{x}{\delta}\right)^n\tag{22}

    其中σ\sigma为PML的厚度，d0d_0是理论反射系数RR中入射角θ=0\theta = 0的情况，有

    R(θ)=e−(2/(n+1))(d0δ/ϵ0Vp)cos⁡θ(23) R(\theta) = e^{-(2/(n+1))(d_0\delta/\epsilon_0 V_p)\cos{\theta}}\tag{23}

    其中VpV_p为波速，ϵ0\epsilon_0是啥我也不太清楚。当nn取n=2n=2时则有

    d0=−log⁡(R)3Vp2δ(24) d_0 = -\log(R)\frac{3V_p}{2\delta}\tag{24}

正演实验

  • 先通过波场快照的动画，看看吸收效果。这里依然使用的是FCNVMB中的速度模型

    结果
    图2. 波场快照动画
    图2. 波场快照动画

    通过图2可以看出吸收效果还是非常好的，然后我们在地面设置检波器进行正演。

    图3. 单炮数据正演
    图3. 单炮数据正演

    可以看到效果相比[4]中没有吸收边界的情况好太多了。

参考

  1. 王守东. 声波方程完全匹配层吸收边界[J].石油地球物理勘探,2003(01):31-34+115-112. ↩︎

  2. SG-PMLopen in new window ↩︎

  3. Berenger J P. A perfectly matched layer for the absorption of electromagnetic waves[J]. Journal of computational physics, 1994, 114(2): 185-200. ↩︎

  4. 波动方程的有限差分解（地震数据正演） ↩︎

Xenny
https://xenny.wiki/posts/fwi/experiment/pml.html
本博客所有文章除特别声明外，均采用 CC BY-NC-SA 4.0 许可协议。转载请注明来自 Xenny 的博客！
2024 © Xenny All rights reserved. 著作权归 Xenny 所有