摘要:实际介质普遍具有黏弹性, 波在传播过程中常伴有能量的耗散、相位畸变和频带变窄等, 对于含有液体和气体的介质, 衰减现象尤为突出. 由于经典的波动理论未考虑介质的黏弹性效应, 基于完全弹性假设的模拟波场和实际传播特征之间的差异明显, 波动理论在工程技术中的应用效果还有提升的空间. 在岩石物理中, 品质因子Q是量化地震衰减强度的参数, 为了研究波在地下介质中的传播规律, 本文从常Q理论出发, 在黏弹性介质频散关系中, 利用多项式拟合和Taylor展开法将频率的分数阶转化为整数阶, 进而推导了时间域复数形式的地下黏弹性介质波动方程. 该近似处理避免了频散关系经域转换后出现分数阶时间微分项, 能有效地降低计算成本. 最后, 采用有限差分法联合伪谱法对均质模型实现了波场的数值模拟, 验证了方程的有效性.
关键词: 黏弹性介质/
波动方程/
数值模拟/
伪谱法

English Abstract

全文HTML

--> --> -->

地下介质普遍具有黏弹性, 波在传播过程中部分能量转变成热, 波场特征会发生明显变化, 如振幅衰减、相位畸变和频带变窄等. 当地层中含有气体或液体时, 黏弹性效应更加突出. 忽略黏弹性的波场模拟易导致波场描述失真, 不利于工程技术中的应用^[1-4], 因此有必要开展黏弹性介质波动理论及波场数值模拟的研究.
在岩石物理中, 品质因子Q用于计算应力与应变之间的相位差, 是度量地震波衰减特征的参数. Liu等^[5]利用广义标准线性体建立了不依赖频率的Q模型, 解释了地震波能量衰减和相位畸变的原因, 理论上Q值与实际介质的匹配度取决于线性体的数目和组合方式^[6,7]. Kjartansson^[8]在衰减系数和频率成线性关系的假设下建立了常Q模型, 由于频散关系的数学形式简单, 在地震资料处理中得到了广泛关注. 从常Q模型的本构关系出发, Carcione等^[9,10]推导了黏弹性介质波动方程, 并采用Grünwald-Letnikov近似实现了地震波场数值模拟, 由于方程中分数阶时间微分项的计算依赖于历史时刻波场, 因此需要保存全时程波场, 计算成本较高. 虽然短时截断法和内部变量法能降低内存的占用^[11,12], 但依然难以应对大尺度模型. Yang和Zhu^[13]推导了分数阶拉普拉斯算子的黏弹性介质波动方程, 并应用于地下构造成像中.
波动方程是双曲线偏微分方程^[14-17], 常用的数值模拟方法有伪谱法、有限元法和有限差分法. 有限差分法无需数据域的转换, 且计算效率高, 是波动方程数值模拟中最常用的算法^[18,19], 为了在相同计算量的情况下提高数值模拟精度, 先后发展了交错网格和旋转交错网格^[20,21]. 伪谱法通过Fourier变换实现, 相当于无限阶差分算子精度, 计算效率不及有限差分法^[22]. 有限元法的优势是网格剖分灵活, 能模拟复杂介质中的波传播, 但计算成本较高^[23].
本文基于Kjartansson提出的常Q模型, 建立了黏弹性介质波动方程, 研究了波动方程的数值模拟方法, 通过数值实验分析了品质因子Q对波场振幅、相位和频带的影响.

在Kjartansson^[8]提出的常Q模型中, 模量$M(\omega )$在频率域的形式为

$M(\omega ) = {{\rm{i}}^{2\gamma }}{\omega ^{2\gamma }}\omega _0^{ - 2\gamma }{M_0}, $

其中${\omega _0}$是参考频率; ${M_0}$是模量; $\omega $是频率; ${\rm{i}}$是虚数单位; $\gamma $与品质因子Q相关, 可通过下式确定^[8]:

$\frac{1}{Q} = \tan (\pi \gamma ).$

式中, $\pi \gamma $的物理意义是介质应力和应变之间相位差.
在地下黏弹性介质中波速${v_{\rm{c}}}$与量$M(\omega )$的关系为^[24]

$v_c^{} = \sqrt {\frac{{M(\omega )}}{\rho }} ,$

式中, $\rho $为介质密度. 相速度${v_{\rm{p}}}$可利用下式计算,

${v_{\rm{p}}} = \frac{\omega }{{{k_{\rm{R}}}}}, $

式中, ${k_{\rm{R}}}$是波数$k$的实部, 利用等价关系${\rm{i}} \!=\! \cos ({\text{π}}/{2}) $$ + {\rm{i}}\sin {({\text{π}}} / 2)/2$, $k$的实部${k_{\rm{R}}}$和虚部${k_{\rm{I}}}$根据(5)式和(6)式确定,

${k_{\rm{R}}} = \operatorname{Re} \left( {\frac{\omega }{{{v_{\rm{c}}}}}} \right) = {\omega ^{1 - \gamma }}\omega _0^\gamma v_0^{ - 1}, $

${k_{\rm{I}}} = \operatorname{Im} \left( {\frac{\omega }{{{v_{\rm{c}}}}}} \right) = - \tan \left( {\frac{{\gamma {\rm{\pi }}}}{2}} \right){\omega ^{1 - \gamma }}\omega _0^\gamma v_0^{ - 1},$

式中, Re和Im分别为取实部和虚部算子; ${v_0} = $$ \sqrt {{{{M_0}} / \rho }} {\cos ^{ - 1}}({{\gamma {\text{π}}} / 2})$, 是在参考频率处的相速度. 将(5)式代入(4)式可得

${v_{\rm{p}}} = {v_0}{\omega ^\gamma }\omega _0^{ - \gamma }.$

类似于弹性介质, 在黏弹性介质中波的频散关系可表示为

$\frac{{{\omega ^2}}}{{v_{\rm{c}}^2}} = {k^2}.$

将(1)式和(3)式代入(8)式有

${\omega ^{2 - 2\gamma }}\omega _0^{2\gamma }\cos (\gamma \pi ) - {\rm{i}}{\omega ^{2 - 2\gamma }}\omega _0^{2\gamma }\sin (\gamma {\rm{\pi )}} = v_{}^2{k^2}.$

式中, ${v^2} = {{{M_0}} / \rho }$, $v$是Q趋于无穷大时的相速度, 即介质趋于弹性介质时波的相速度. 若将(9)式从频率-波数域变换至时间-空间域, 可得含有分数阶时间微分的波动方程, 虽然采用Grünwald-Letnikov近似能实现分数阶微分项的计算, 但该方法需保存历史时间波场. 为此, 将(9)式等号左侧第一项进行多项式拟合有

${\omega ^{2 - 2\gamma }}\omega _0^{2\gamma }\cos (\gamma {\rm{\pi )}} = {a_1}{\omega ^2} + {a_2}\omega + {a_3}.$

式中, ${a_1}$, ${a_2}$和${a_3}$为拟合系数, 表1提供了参考频率为50 Hz时不同品质因子Q对应的拟合系数. 当Q = 30时, (10)式拟合曲线如图1所示, 图中实心点是理论数值, 实线为拟合曲线, 以上拟合均是将频率转换至圆频率所得.

Q	5	10	20	30	60	100	5000
a1	0.8144	0.9081	0.9545	0.9698	0.9850	0.991	0.9998
a2	58.452	30.728	15.662	10.499	5.276	3.1716	0.0636
a3	–2078.2	–1133.1	–587.58	–396.08	–200.14	–120.57	–2.4255

表1不同Q值拟合系数
Table1.Fitting coefficients of different Q values.

图 1 拟合曲线 (Q = 30)
Figure1. Fitting curve (Q = 30).

(9)式中等号左侧第二项中${\rm{i}}\omega $经Fourier变换是时间的一阶微分算子, 为了使剩余的${\omega ^{1 - 2\gamma }}$易于计算, 先对${k^{2\gamma - 1}}$进行处理,

$k_{}^{2\gamma - 1} = {(k_{\rm{R}}^{} + {\rm{i}}k_{\rm{I}}^{})^{2\gamma - 1}} = k_{\rm{R}}^{2\gamma - 1}{\left( {1 + {\rm{i}}\frac{{{k_{\rm{I}}}}}{{{k_{\rm{R}}}}}} \right)^{2\gamma - 1}}.$

将(5)式和(6)式代入(11)式得

$k_{}^{2\gamma - 1} = k_{\rm{R}}^{2\gamma - 1}{\left( {1 - {\rm{i}}\tan \frac{{\gamma \pi }}{2}} \right)^{2\gamma - 1}}.$

对于弱衰减介质有$\gamma \ll 1$, 利用Taylor展开可将${[1 - {\rm{i}}\tan ({{\gamma \pi } / 2})]^{2\gamma - 1}}$表示为

${\left( {1 - {\rm{i}}\tan \frac{{\gamma \pi }}{2}} \right)^{2\gamma - 1}} = 1 + {\rm{i}}\frac{{\gamma \pi }}{2} + o{(r)^2}.$

忽略(13)式的高阶无穷小量后代入(12)式有

$k_{\rm{R}}^{1 - 2\gamma } = k_{}^{1 - 2\gamma }\left( {1 + {\rm{i}}\frac{{\gamma \pi }}{2}} \right).$

(4)式的等价形式为

${\omega ^{1 - 2\gamma }} = v_{\rm{p}}^{1 - 2\gamma }k_{\rm{R}}^{1 - 2\gamma }.$

根据(7)式, 当参考频率接近于波的频率时, 有${v_{\rm{p}}} \approx $$ {v_0}$, 将(14)式代入(15)式为

${\omega ^{1 - 2\gamma }} = v_0^{1 - 2\gamma }k_{}^{1 - 2\gamma }\left( {1 + {\rm{i}}\frac{{\gamma \pi }}{2}} \right).$

将(10)式和(16)式代入(9)式有

$\begin{split}& {a_1}{\omega ^2} + {a_2}\omega + {a_3} - {\rm{i}}\tau \omega k_{}^{1 - 2\gamma } \\& + {{\gamma \pi \tau \omega k_{}^{1 - 2\gamma }}}/{2} = v_{}^2{k^2}, \end{split} $

其中$\tau = \omega _0^{2\gamma }v_0^{1 - 2\gamma }\sin (\gamma \pi )$.

$\begin{split}&F\begin{array}{*{20}{c}} {[{\omega ^2}}&\omega &{{\rm{i}}\omega }&{{k^\xi }} \end{array}] \\=\;& \left[ {\begin{array}{*{20}{c}} { - \dfrac{{{\partial ^2}}}{{\partial {t^2}}}}&{ - {\rm{i}}\dfrac{\partial }{{\partial t}}}&{\dfrac{\partial }{{\partial t}}}&{{{\left( { - {\nabla ^2}} \right)}^{\tfrac{\xi }{2}}}} \end{array}} \right].\end{split}$

根据Fourier变换的微分性质, (18)式给出了算子在频率-波数域和时间-空间域的对应关系, F为Fourier变换算子, $\nabla $是拉普拉斯算子. 将(17)式乘频率-波数域波场后, 经域转换可得时间-空间域地下黏弹性介质波动方程:

$\begin{split}&{a_1}\frac{{{\partial ^2}u}}{{\partial {t^2}}} + {\rm{i}}{a_2}\frac{{\partial u}}{{\partial t}} - {a_3}u + \tau {\left( { - {\nabla ^2}} \right)^{0.5 - \gamma }}\frac{{\partial u}}{{\partial t}}\\&+ {\rm{i}}\frac{{\gamma \pi \tau {{\left( { - {\nabla ^2}} \right)}^{0.5 - \gamma }}}}{2}\frac{{\partial u}}{{\partial t}} = v_{}^2{\nabla ^2}u\end{split},$

(19)式中, $u$是复数波场, 取其实部即物理意义上的波场. 为了在数值模拟中研究波动方程中各项对波场的主控因素, 可将(19)式改写为

$\begin{split}&\frac{{{\partial ^2}u}}{{\partial {t^2}}} + \lambda \left[ {({a_1} - 1)\frac{{{\partial ^2}u}}{{\partial {t^2}}} + {\rm{i}}{a_2}\frac{{\partial u}}{{\partial t}} - {a_3}u} \right]\\&+ \mu \left[ {\tau {{\left( { - {\nabla ^2}} \right)}^{0.5 - \gamma }}\frac{{\partial u}}{{\partial t}} + {\rm{i}}\frac{{\gamma \pi \tau {{\left( { - {\nabla ^2}} \right)}^{0.5 - \gamma }}}}{2}\frac{{\partial u}}{{\partial t}}} \right] \\=\;& v_{}^2{\nabla ^2}u.\\[-10pt]\end{split}$

(20)式中参数$\lambda $和$\mu $的取值为1或0, 两参数组合可形成四个方程, 当$\lambda = 1,\; {\rm{ }}\mu = 1$时, (20)式和(19)式等价; 当$\lambda = 0,\; {\rm{ }}\mu = 0$时, (20)式退化为声波方程.

本文采用二阶精度有限差分法求解(19)式中的时间微分项, 其离散格式为

$\frac{{{\partial ^2}}}{{\partial {t^2}}}{u^n} = \frac{{{u^{n + 1}} - {\rm{2}}{u^n} + {u^{n - 1}}}}{{\Delta {t^2}}}, $

$\frac{\partial }{{\partial t}}{u^n} = \frac{{{u^n} - {u^{n - 1}}}}{{\Delta t}}.$

式中, $n$为时间节点; $\Delta t$是时间采样间隔.
(19)式中的空间微分项采用伪谱法计算, 先通过Fourier变换将时间-空间域波场转换至波数域, 与波数$k$和拉普拉斯算子的阶数作用后, 再经反Fourier变换即可实现求解^[25,26]. 以上过程可用(23)式表示:

${\left( { - {\nabla ^2}} \right)^\xi }u = {F^{ - 1}}\left[ {{k^{2\xi }}F\left( u \right)} \right].$

结合(21)式、(22)式和(23)式, 波动方程(19)式的波场递推公式为

$\begin{split} {u^{n + 1}} =\;& \frac{{\Delta {t^2}}}{{{a_{\rm{1}}}}}\bigg\{ v^2{F^{-1}} [-K^2F(u)] - {\rm{i}}\Big({a_2}\frac{{{u^n} - {u^{n - 1}}}}{{\Delta t}} \\& + \frac{{\gamma \pi \tau }}{2}{F^{-1}}\left( {{K^{1 - 2\gamma }}F\left( {\frac{{{u^n} - {u^{n - 1}}}}{{\Delta t}}} \right)} \right) \Big) \\& +{{a_3}{u^n} - \tau {F^{-1}}\left( {{K^{1 - 2\gamma }}F\left( {\frac{{{u^n} - {u^{n - 1}}}}{{\Delta t}}} \right)} \right) }\bigg\} \\& + 2{u^n} - {u^{n - 1}}. \\[-12pt] \end{split} $

为研究地下黏弹性介质中波的传播规律, 设计水平和垂直方向的距离为3.8 km的均质模型, 空间网格间距为10 m, 波在参考频率50 Hz处的相速度为2216 m/s. 利用主频是20 Hz的Ricker子波作为激发震源, 其振幅峰值为1 Pa, 时间采样间隔1 ms, 子波的函数形式为

${\rm{Ricker}} (t) \!=\! [1 - 800{\pi ^2}{(t - 0.1)^2}]{{\rm{e}}^{ - 400{\pi ^2}{{(t - 0.1)}^2}}}.$

震源位于模型中央位置, 距离震源0.6 km处布置检波器, 数值实验模型如图2所示. 设置模型品质因子Q为20, 编程验证(19)式的有效性和数值模拟方案的可行性. 图3为600 ms的波场, 由于介质中不含波阻抗界面, 因此波场特征是外扩的圆. 在模型垂直方向1.9 km处, 水平方向0—3.8 km范围分别提取500和600 ms波场, 如图4所示, 随着时间推移, 波场的能量由1.9 km处向外传递, 受到介质的吸收衰减, 波场的振幅明显变小.
图 2 数值模拟模型
Figure2. Numerical simulation model.

图 3 600 ms波场快照
Figure3. Wavefield snapshot at 600 ms.

图 4 不同时刻的1维波场
Figure4. One dimensional wavefield at different times.

介质的黏弹性不仅衰减波场的振幅, 而且影响相位和频带. 现采用不同Q值实施波场模拟, 检波器从震源激发时开始采集信息, 直至1000 ms结束, 在模拟区域外侧添加吸收衰减层以避免边界产生虚假反射. 图5是在相同传播路径下, 不同Q值所对应的记录, 由于检波器接收的信息集中在380 ms, 排除无效信息仅对200至600 ms记录显示, Q值由30变为10后, 波场的振幅明显减小, 相位发生畸变. 图5中两记录的归一化频谱如图6所示, Q值影响了记录频带的特征, 主频由初始的20 Hz分别降至15和10 Hz. 在Kjartansson提出的常Q模型中, Q的减小对应着应力和应变之间的滞后变大, 进而介质对波场的吸收衰减变强, 此外该模型的假设是衰减系数与波场的频率成正比关系, 即低频成分的衰减比高频成分弱, 通过数值实验表明, 建立的波动方程能够准确地描述波场的衰减规律, 且数值模拟方案可行.
图 5 波场记录
Figure5. Wavefield record.

图 6 记录频谱
Figure6. Record spectrum.

为了对比分析(20)式中各项对波场的主控因素, 采用控制参数$\lambda $和$\mu $的方式对方程实施数值模拟, 图7(a)—(d)分别对应$\lambda = 0,\; {\rm{ }}\mu = 0$, $\lambda = 0, $$ \; {\rm{ }}\mu = 1$, $\lambda = 1,\; {\rm{ }}\mu = 0$和$\lambda = 1,\; {\rm{ }}\mu = 1$的截取波场, 图7(a)与图7(b)和图7(c)的差异分别在于振幅和相位, 图7(a)和图7(d)在振幅和相位上均有差异. 基于以上研究, (20)式中等号左侧第二项主控相位畸变, 第三项主控振幅衰减. 图8为4个波场的局部(图2中A点至B点)单道显示, 该图展示了弹性介质波场经振幅衰减和相位畸变后演变为黏弹性介质波场的过程.
图 7 组合波场
Figure7. Combined wavefield.

图 8 部分波场
Figure8. Part wavefield.

为了测试数值方法对于复杂模型的有效性, 进行了波场数值模拟. 采用模型在水平和垂直方向的网格点数分别为401和301, 网格间距为10 m, 即水平距离为4 km, 垂直深度为3 km, 模型速度分布在2500至3500 m/s区间, 以层状构造为主, 如图9(a)所示, 品质因子Q取25. 震源布置在模型中央, 三个不同时刻的波场如图9(b),(c)所示. 随着时间递增, 波场能量逐步外扩, 由于介质的黏滞性效应, 波场振幅出现了衰减特征. 通过测试表明, 构建的黏弹性介质波动方程能刻画波传播至波阻抗界面所产生的透射波和反射波.
图 9 复杂模型波场模拟　(a) 速度模型; (b) 0.3 s波场; (c) 0.4 s波场; (d) 0.5 s波场
Figure9. Wavefield simulation for the complex model: (a) Velocity; (b) wavefield at 0.3 s; (c) wavefield at 0.4 s; (a) wavefield at 0.5 s.

本文在常Q理论基础上, 建立了黏弹性介质波动方程, 该方程能解耦振幅衰减和相位畸变效应, 有助于分析介质中的复杂波场现象, 可为地震资料处理提供理论基础. 需要注意的是方程中拉普拉斯算子的阶数是Q的函数, 伪谱法仅能处理Q场为常数的情况, 对于Q场较平滑的介质, 可将Q场的平均应用至全局. 为了应对Q场具有强非均质性的波场模拟, 如何将Q从拉普拉斯算子的阶数中分离值得深入研究.