[发明专利]一种基于弥散黏滞性波动方程的FCT‑FDM正演模拟方法

摘要

本发明公开了一种基于弥散黏滞性波动方程的FCT‑FDM正演模拟方法，从数值模拟的角度研究了弥散黏滞性波的传播特性。首先阐述了FCT‑FDM的计算步骤，设计一种流体饱和质模型利用该方法模拟弥散黏滞性波的传播，并将数值结果与声波模拟结果做对比分析。结果表明弥散黏滞性波相对于声波在振幅上有明显的衰减，相位也会发生显著的变化。弥散黏滞性波的衰减程度主要取决于弥散衰减系数和黏滞性衰减系数的取值。本发明为地震波的数值模拟提供了一种有效的工具，尤其对于含有流体的介质，它能用于刻画地震波在流体介质中的衰减和相移特性。此外，该技术方案易于实现，可操作性强。

主权项

一种基于弥散黏滞性波动方程的FCT‑FDM正演模拟方法，其特征在于，包括以下步骤：1)根据实际地质背景条件、岩石物理测试及其测井资料建立介质模型的参数；2)对所设计的介质模型进行网格离散；3)震源函数在空间上采用高斯函数，时间上采用Ricker子波，形式为：s(x,z,t)＝g(x,z)·f(t)式中：f(t)＝(1‑2(πf0t)2)exp(‑(πf0t)2)g(x,z)=exp{-[(x-x0)2+(z-z0)2]}β2]]>式中：f0表示Ricker子波的中心频率，模型计算中f0＝30Hz；β为常数；(x0,z0)表示震源的空间位置；4)将二维弥散黏滞性波动方程中的微分用差分近似替代，得到其相应的有限差分格式；其中空间采样步长和时间采样步长必须满足以下数值格式的稳定性条件：0≤Δt≤min(64hυ,h2γh2+8η)]]>其中，h为空间网格步长；5)对二维弥散黏滞性波动方程进行FCT‑FDM计算，得到每一时刻的波场值，其具体方法如下：首先，给出弥散黏滞性波动方程的数学描述，二维弥散黏滞性波动方程其形式为：∂2u∂t2+γ∂u∂t-η(∂3u∂x2∂t+∂3u∂z2∂t)-υ2(∂2u∂x2+∂2u∂z2)=0---(1)]]>式中：u为波场函数；γ，η分别为弥散和黏滞性衰减系数，它们与岩石的孔隙度、渗透率以及流体的密度、粘度等有关；υ为非频散介质中波的传播速度；x，t分别为空间和时间变量；式(1)中左端第一项表示惯性项，第二项为弥散耗损力即扩散项，第三项表示黏滞性阻尼，第四项为波动方程的弹性部分；弥散黏滞性波动方程FCT‑FDM方法的具体的计算步骤如下：5‑1)定义式(1)的数值解为给定初值为震源函数；其中，表示第n时间步在网格点(xj,zk)处的波场值；n为时间采样点，j为空间x方向的采样点，k为空间z方向的采样点；5‑2)利用弥散黏滞性波动方程有限差分格式计算uj,kn+1=[(2-4a-γ(Δt)-4b]uj,kn-[1-γ(Δt)-4a]uj,kn-1-a(uj+1,kn-1+uj-1,kn-1+uj,k+1n-1+uj,k-1n-1)+(a+b)[(uj+1,kn+uj-1,kn+uj,k+1n+uj,k-1n)]---(2)]]>式中：h为空间网格步长；Δt为时间步长；5‑3)计算漫射通量：a)计算第n时间步的漫射通量Pj,kn=λ1(uj+1,kn-uj,kn-uj+1,kn-1+uj,kn-1),0≤λ1≤1Qj,kn=λ1(uj,k+1n-uj,kn-uj,k+1n-1+uj,kn-1),0≤λ1≤1---(3)]]>其中，λ1为参数；b)利用漫射通量平滑差分方程的解u~j,kn+1=uj,kn+1+(Pj+12,kn-Pj-12,kn)+(Qj,k+12n-Qj,k-12n)---(4)]]>其中，为平滑以后的解；5‑4)计算反漫射通量：a)计算第n+1时间步的漫射通量P~j,kn+1=λ2(uj+1,kn+1-uj,kn+1-uj+1,kn+uj,kn),0≤λ1≤1Q~j,kn+1=λ2(uj,k+1n+1-uj,kn+1-uj,k+1n+uj,kn),0≤λ2≤1---(5)]]>式中：λ2与λ1类似，为参数；b)利用和计算反漫射通量Xj+1,kn+1=u~j+1,kn+1-uj,kn+1-(u~j+1,kn-uj,kn)Yj,k+1n+1=u~j,k+1n+1-uj,kn+1-(u~j,k+1n-uj,kn)---(6)]]>c)利用反漫射通量修正平滑以后的解Uj,kn+1=u~j,kn+1-(X~j+12,kn-X~j-12,kn)-(Y~j,k+12n-Y~j,k-12n)---(7)]]>式中：为修正以后的解，且X~j+12,kn=(sx·max{0,min[sx·Xj-12,kn,|P~j+12,kn+1|,sx·Xj+32,kn]})Y~j,k+12n=(sy·max{0,min[sy·Yj,k-12n,|Q~j,k+12n+1|,sy·Yj,k+32n]})---(8)]]>Xj-12,kn=12(Xj-1,kn+Xj,kn),Xj+32,kn=12(Xj+1,kn+Xj+2,kn),Yj,k-12n=12(Yj,k-1n+Yj,kn),Yj,k+32n=12(Yj,k+1n+Yj,k+2n)---(9)]]>P~j+12,kn+1=12(P~j+1,kn+1+P~j,kn+1),Q~j,k+12n+1=12(Q~j,k+1n+1+Q~j,kn+1)---(10)]]>sx=sign(P~j+12,kn+1),sy=sign(Q~j,k+12n+1)---(11)]]>(11)式中“sign”为符号函数；(8)式中max(·)，min(·)分别表示物理量的最大值和最小值；5‑5)重复步骤5‑2)～5‑4)，直到最大的计算时间步，得到整个计算区域内不同时间的波场值；6)在计算区域的边界处使用吸收边界条件并对其进行处理，模拟实际地下介质中波的传播。