地球物理学报  2020, Vol. 63 Issue (6): 2400-2414   PDF    
基于Remez迭代算法求解差分系数的双相介质自适应有限差分正演模拟
王文化1,2, 文晓涛1,2, 张波3, 张懿疆1,2, 王威4     
1. 成都理工大学油气藏地质及开发工程国家重点实验室, 成都 610059;
2. 成都理工大学地球物理学院, 成都 610059;
3. 美国阿拉巴马大学地质科学系, 塔斯卡卢萨 35401;
4. 中国石油化工股份有限公司勘探分公司, 成都 610041
摘要:利用传统有限差分方法对基于Biot理论的双相介质波动方程进行数值求解时,由于慢纵波的存在,数值频散效应较为明显,影响模拟精度.相对于声学近似方程及普通弹性波方程,Biot双相介质波动方程在同等数值求解算法和精度要求条件下,其地震波场正演模拟需要更多的计算时间.本文针对Biot一阶速度-应力方程组发展了一种变阶数优化有限差分数值模拟方法,旨在同时提高其正演模拟的精度和效率.首先结合交错网格差分格式推导Biot方程的数值频散关系式.然后基于Remez迭代算法求取一阶空间偏导数的优化差分系数,并用于Biot方程的交错网格有限差分数值模拟.在此基础上把三类波的平均频散误差参数限制在给定的频散误差阈值和频率范围内,此时优化有限差分算子的长度就能自适应非均匀双相介质模型中的不同速度区间.数值频散曲线分析表明:基于Remez迭代算法的优化有限差分方法相较传统泰勒级数展开方法在大波数范围对频散误差的压制效果更明显;可变阶数的优化有限差分方法能取得与固定阶数优化有限差分方法相近的模拟精度.在均匀介质和河道模型的数值模拟实验中将本文变阶数优化有限差分算法与传统泰勒展开算法、最小二乘优化算法进行比较,进一步证明其在复杂地下介质中的有效性和适用性.
关键词: 正演模拟      Remez迭代算法      最小二乘      双相介质      自适应有限差分     
Adaptive finite difference forward modeling of two-phase media based on the Remez iterative algorithm for solving differential coefficients
WANG WenHua1,2, WEN XiaoTao1,2, ZHANG Bo3, ZHANG YiJiang1,2, WANG Wei4     
1. State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Chengdu University of Technology, Chengdu 610059, China;
2. College of Geophysics, Chengdu University of Technology, Chengdu 610059, China;
3. Department of Geological Sciences, The University of Alabama, Tuscaloosa 35401, USA;
4. Exploration Company of SINOPEC, Chengdu 610041, China
Abstract: When the finite difference method (FDM) is used to numerically solve the two-phase medium wave equation based on Biot theory, the numerical dispersion effect is obvious due to the existence of slow compressional waves (P-wave), which affects the accuracy of numerical simulation. Compared with the acoustic approximation equation and the ordinary elastic wave equation, the Biot two-phase medium wave equation requires more calculation time for the seismic wave-field forward modeling under the same numerical solution algorithm and accuracy requirements. To solve this problem, this work developed a variable-order optimal FDM (VOFDM) for Biot first-order velocity-stress equations to improve the accuracy and efficiency of the forward modeling. Firstly, the numerical dispersion relation of Biot equation is derived by combining the staggered-grid FD (SFD) scheme. Then, based on the Remez iterative (RI) algorithm, the optimal differential coefficients of the first-order spatial partial derivatives are obtained and used for the SFD of the Biot equation. On this basis, the average dispersion error parameters of three kinds of waves are limited to a given dispersion error threshold and frequency range, so that the length of the FD operator can be adaptive to different velocity intervals in the inhomogeneous two-phase media model. The numerical dispersion curve analysis shows that the optimized FDM based on the RI algorithm is more effective in suppressing the dispersion error in the large wave number range than the traditional Taylor-series expansion (TE) method. The VOFDM can achieve similar simulation accuracy with the fixed-order optimized FDM (FOFDM). The numerical simulation experiments on homogeneous medium and river models prove VOFDM's effectiveness and applicability in complex underground media compared with the TE algorithm and least squares (LS) optimization algorithm.
Keywords: Forward modeling    Remez iterative algorithm    Least squares    Two-phase medium    Adaptive finite difference    
0 引言

实际地下介质, 特别是油气储层, 往往是由固体骨架和孔隙流体组成的双相介质或多相孔隙介质.因此, 相比常规的声学或单相弹性介质近似, 研究双相介质中地震波的传播规律更符合当今地球物理勘探需求.Biot理论(Biot, 1956a, b, 1962)被认为是进行孔隙介质中地震波传播研究的理论基础, 从理论角度证明了双相介质中的慢纵波的存在.Plona(1980)Berryman(1980)在岩石物理实验中观测到了慢纵波的存在, 从而证明了Biot理论的正确性.此后, 国内外学者针对双相介质的波场数值模拟问题进行了大量的研究工作, 发展了多种数值求解算法.比如有限差分法(Zhu and McMechan, 1991; Dai et al., 1995; Carcione and Helle, 1999; Sheen et al., 2006; 裴正林, 2006; O′Brien, 2010; 刘财等, 2014)、伪谱法(Sidler et al., 2010)、有限元法(张金波等, 2018)、有限体积法(Lemoine et al., 2013)等.已经用于双相介质波动方程正演模拟的有限差分方法有:交错网格方法(Özdenvar and McMechan, 1997; Masson et al., 2006; O′Brien, 2010; Zhang et al., 2018a); 变网格差分法(Liu et al., 2014); 时空域差分法(Zhang and Gao, 2014b); 隐式差分法(Itzá et al., 2016)等.相比其他方法, 有限差分法具有易于实施、灵活多变和计算高效的优点, 因此广泛应用于各类波动方程的正演模拟(Alterman and Karal, 1968; Dablain, 1986; Robertsson et al., 1994; 董良国等, 2000; Saenger and Bohlen, 2004; Zhang and Gao, 2014a; 刘财等, 2018), 但是因网格离散带来的数值频散是不可避免的, 因此波场模拟的精度会受到影响.此外, 对于双相介质波动方程而言, 由于慢纵波比其他两类波(横波、快纵波)的传播速度小, 使得有限差分方法在此类波动方程中的空间网格数值频散问题更为突出.

通常用来压制数值频散效应的方法有以下几种:其一是采用较小的空间间隔、时间步长(Liu et al., 2014)或者高阶差分算子(Liu and Wei, 2008), 但是这些方法需要较大的计算量(刘洋, 2014; 辛维等, 2015);其二是结合FCT(Flux-Corrected Transport)技术缓解数值频散现象(杨宽德等, 2002; 2011; 周竹生和唐磊, 2012; 李立平和何兵寿, 2017); 其三是利用数值优化算法求取新的差分系数进行波场正演模拟(Holberg, 1987; Zhang and Yao, 2013a; Liu, 2014; Yang et al., 2015; 李振春等, 2016).目前用于双相介质波动方程的有限差分方法主要采用前面两种方法压制数值频散, 其都采用传统泰勒级数展开(TE)求解差分系数,TE方法可以在较小波数区域取得很高的精度, 但是其频散误差随波数的增加而显著增大.相同阶数下, 常系数优化算法能在较大的波数区间取得更高的模拟精度.尽管常系数优化有限差分方法已经广泛用于提高各类波动方程的正演模拟精度, 但是针对双相介质波动方程的研究较少涉及.张杰和吴国忱(2015)采用最小二乘法实现了基于优化差分系数的双相介质交错网格正演模拟.Itzá等(2016)将基于最小二乘优化的隐式交错网格有限差分方法应用到双相介质的波场数值模拟中, 提高了数值模拟的精度, 但是计算过程较为耗时并且降低了正演模拟的稳定性.

常系数优化有限差分方法主要用频散关系近似优化偏导数的有限差分系数(Liu, 2013).典型的两种方法是采用最小二乘(LS)近似和极大极小近似(MA)方法来建立目标函数求解有限差分系数, 它们分别属于二范数和最大范数范畴(Yang et al., 2017).其中基于LS建立的目标函数形式较为固定, 这在一定程度上限制了它的应用(Zhang and Yao, 2013b).基于MA建立目标函数较为灵活, 这导致它的误差相对较小(Holberg, 1987; Zhang and Yao, 2013b).但是此类方法的难点在于目标函数的求解, 并且计算量较大.另外一种常见方法是利用频散关系式推导线性化方程组求解优化差分系数.基于此类线性化数值优化算法的有限差分方法可以在大波数区间取得较高精度的同时基本不会带来额外的计算量, 因此近年来备受关注.辛维等(2015)对比了模拟退火、最小二乘和采样近似(SA)这三种优化有限差分算法的频散压制效果和计算效率,认为SA方法有着与其他两种方法相近的压制频散效果, 并且计算量远远小于其他两种方法.Yang等(2016)基于采样近似方法实现了优化隐式交错网格弹性波数值模拟.Yan等(2016)提出一种新的常系数优化算法(TS), 其结合泰勒级数展开和采样近似方法求解差分系数, 数值频散分析表明这种方法能提高频散带宽, 但是其在小波数区域的数值频散误差并未变小.Yang等(2017)采用最大范数建立目标函数, 用Remez迭代算法(Remez, 1934)求解一阶空间偏导数的优化交错网格差分系数, 数值频散分析表明此方法相比采样近似法在较大波数时的频散误差更小.此方法在每次迭代求取差分系数的过程中采用采样近似方法, 只需迭代三到四次就能获得满意的计算结果.He等(2019)认为Yang等(2017)在目标函数中引入的不恰当零波数约束条件会导致数值频散曲线收敛不到等波纹解,因此引入了新的约束条件去抑制零波数附近的数值频散误差.徐世刚和刘洋(2018)Yan等(2016)研究的基础上利用Remez算法求解目标函数获得二阶空间偏导数的优化有限差分系数, 数值频散曲线分析表明此方法精度略高于LS方法, 同时将其成功用于三维VTI介质声波和弹性波方程的正演模拟.本文将此算法推广到一阶偏导数的有限差分系数优化中, 并且用于双相介质波动方程的相关数值模拟工作中.

另外, 由于双相介质波动方程相比普通声波和弹性波方程参数更多、形式更为复杂, 需要更大的计算量, 涉及到三维问题的时候更是如此(Zhang et al., 2018b), 所以采用高效的有限差分方法解决此类波动方程的正反演问题显得尤为重要.Liu和Sen(2011)提出一种用于声波方程的自适应变阶数有限差分方法.其核心思想是地下介质不同的速度区域对应不同阶数的有限差分算子, 在较低的速度区域采用较高阶数的差分算子, 在较高的速度区域采用较低阶数的差分算子(Yao et al., 2016).数值算例证明这种方法可以在减小计算量的同时确保模拟精度.Gao和Zhang(2013)将此方法用于二维双相介质波动方程的数值模拟中, 并且后来成功推广至三维情况下的数值模拟(Zhang et al., 2018b).本文在此基础上结合新的常系数优化算法进行双相介质波动方程变阶数有限差分正演模拟.即在给定的频散误差阈值和频率范围内, 优化有限差分算子的不同差分阶数可以自适应非均匀双相介质模型中的不同速度区段.此方法相比变网格或者变时间步长算法更容易实现, 并且不需要特殊的边界条件处理.数值频散对比分析和正演模拟算例证明变阶数优化有限差分方法不但能减小计算量, 而且能提高数值模拟精度.

1 理论与方法 1.1 双相介质波动方程及其数值频散关系式

根据Dai等(1995)的论文, 基于Biot理论的二维一阶速度-应力弹性波方程为

(1)

此方程组的系数矩阵见附录A.对于方程组(1)的频率波数域频散关系式可以表示为(Zhang and Gao, 2014a):

(2)

(3)

(4)

其中vS, vsPvfP分别表示横波、慢纵波和快纵波的速度.这三种波的速度分别等于:

(5)

(6)

(7)

其中,, Y=P2R222-4PQR12R22-2PRR11R22+4PRR122+4Q2R11R22+R2R112-4RQR11R12.P, Q, R为Biot弹性常数.

同等模拟参数条件下,基于交错网格差分格式的有限差分方法相对于常规网格方法具有更高的精度, 本文采用此方法求解方程组(1).在此方法中, 二阶精度的时间偏导数和2M阶精度的空间偏导数分别如式(8)、式(9)所示:

(8)

(9)

其中τ是时间步长, h是网格步长, u可以是质点速度、固体压力或者流体压力, cm(m=1, 2, …, M)是交错网格有限差分系数.

根据平面波理论, 假设:

(10)

其中u0是一个常量, , k表示波数.将公式(10)分别代入公式(8)和(9), 可得:

(11)

(12)

同理可得:

(13)

将以上3式代入(3)式可得快纵波的数值频散关系式:

(14)

类似可得慢纵波和横波的数值频散关系式.令kx=kcosθ, kz=ksinθ, 其中θ是和波传播方向有关的角度.则二维双相介质波动方程的数值频散关系式可以表示为

(15)

其中ri=viτ/h, i=1, 2, 3.i=1表示慢纵波, i=2表示横波, i=3表示快纵波.根据公式(15)可定义频散误差为

(16)

其中,

(17)

1.2 基于Remez迭代算法的变阶数优化有限差分方法

徐世刚和刘洋(2018)年提出的二阶空间偏导数优化差分方法(MN)推广到一阶空间导数交错网格有限差分系数求解中.假设:

(18)

则方程(12)或(13)可以表示为

(19)

其中.基于公式(19)利用采样近似方法可以得到M个线性方程组(辛维等, 2015; Yang et al., 2015):

(20)

其中βm(m=1, 2, …, M)是均匀分布在区间[0, b]上的一系列参考点.

另外利用方程(19)采用最大范数定义如下的目标函数(Zhang and Yao, 2013b):

(21)

其中,

(22)

根据Remez算法, 令:

(23)

其中E是一个可变的常数.

结合传统泰勒展开方法和以上公式可以得到如下的M阶线性方程组:

(24)

在给定波数范围[0, b]内,基于公式(24)可以利用Remez算法使得频散误差参数ε的绝对值的最大值极小化, 同时解得空间一阶偏导数的优化交错网格有限差分系数cm.具体算法流程如下:

(1) 在[0, b]上选取一系列均匀分布的采样点βm(m=1, 2, …, M), 再代入方程(24)利用高斯消元法求得cmE的值.

(2) 利用现有cm值(通过步骤(1)或者步骤(5)得到)和方程(21)计算A在特定波数区间[0, b]中得最大值.

(3) 对比当下A和|E|的值.如果|E|≥A, 此时获得最优有限差分系数cm, 结束循环迭代; 否则在步骤(4)中将现有采样点与计算得到极值点对应的位置进行交换.

(4) 由方程(23)可知M个参考点上ε(β)的值正负交替出现, 这样必然存在(M-1)个根使得ε(β)=0.通过计算每个极值点对应的β值得到新的βm(m=1, 2, …, M)值.需要注意的是βM始终等于b.

(5) 利用新的交换点βm求解方程(24)获得新的cmE.

(6) 重复步骤(2-5)直到找到最优差分系数cm.

按照以上步骤就可获得优化后的差分系数, 将其代入公式(16)对比三种波的频散误差曲线(如图 13所示,其中τ=1 ms,h=10 m,, vsP=1180 m·s-1, vs=1790 m·s-1, vfP=3210 m·s-1)可以看出:利用传统TE方法和MN方法得到的频散误差都随着差分阶数的增加而减小; 相同阶数下, 采用MN方法得到频散曲线显示在大波数条件下的空间数值频散误差小于TE方法.Liu和Sen(2011)认为有限差分算子的长度可以随速度而改变, 低速区域可以采用较长的空间差分算子(高阶), 高速区域可以采用较短的差分算子(低阶).为了提高双相介质波动方程的计算效率, 可以在优化有限差分系数的基础上进行变阶数有限差分数值模拟, 即变阶数优化有限差分方法.

图 1 不同差分阶数下慢纵波的频散曲线对比 (a)传统方法(TE); (b)本文方法(MN). Fig. 1 Comparison of dispersion curves of the slow P-wave with different finite difference orders (a) Conventional TE method; (b) Proposed MN method of this work.
图 2 不同差分阶数下横波的频散曲线对比 (a)传统方法(TE); (b)本文方法(MN). Fig. 2 Comparison of dispersion curves of the S-wave with different finite difference orders (a) Conventional TE method; (b) Proposed MN method of this work.
图 3 不同差分阶数下快纵波的频散曲线对比 (a)传统方法(TE); (b)本文方法(MN). Fig. 3 Comparison of dispersion curves of the fast P-wave with different finite difference orders (a) Conventional TE method; (b) Proposed MN method of this work.

为了求取自适应差分算子的阶数, 先定义数值频散误差为

(25)

其中ε=0时, 无数值频散, ε≠0时, 会出现数值频散现象(ε<0时对应由于相速度变大而出现的时间数值频散;ε>0时对应由于相速度变小而出现的空间数值频散).

接下来定量分析双相介质中的数值频散, 这里定义一个有关慢纵波、横波和快纵波的频散参数如下:

(26)

图 4表示相同空间差分阶数和不同地震波传播速度时的频散误差(log10|ξ|)曲线.从图中可以看出相同频率时, 地震波传播的速度越大, 频散误差越小, 精度越高.本文采用Liu和Sen(2011)给出的方法决定不同速度情况下的差分阶数.在给定频散误差阈值η和最大频率值fmax的情况下令:

(27)

图 4 不同速度下频散误差log10|ξ|随频率的变化曲线(相同差分阶数,M=8时) Fig. 4 Dispersion parameters log10|ξ| varied with frequencies for different velocities with the same FD order (M=8)

公式(27)提供了确定差分阶数的准则,利用它来确定给定条件下不同速度对应的M最小值,就能使有限差分阶数(2M)自适应地下非均匀介质的速度区间变化.图 5表示不同阶数和不同速度时的频散误差(log10|ξ|)变化曲线, 可以看出在特定频散误差阈值(η=10-5)及频率范围(0~25 Hz)内, 速度为7000 m·s-1时采用八阶差分算子的模拟精度与速度1000 m·s-1时采用十六阶差分算子的精度相当.对比图 4图 5, 可以看出当地下介质速度较大时(比如4000 m·s-1),采用较小的差分阶数(2M=10)可以取得与较大的差分阶数(2M=16)时相近的模拟精度.这证明了本文变阶数优化有限差分数值模拟方法的可行性.图 6是不同速度值、M值和不同频散误差阈值间的关系变化.可以看到有限差分阶数(2M)取决于频散误差阈值和速度值, 并且它随频散误差阈值和速度值的减小而增大.

图 5 不同速度下频散误差log10|ξ|随频率的变化曲线(变阶数时) Fig. 5 Dispersion parameters log10|ξ| varied with frequencies for different velocities with the different FD orders
图 6 不同速度区间和频散误差阈值下的M Fig. 6 The SFD orders for different velocities and error thresholds
2 数值算例 2.1 均匀介质模型数值模拟

这里为了验证本文提出的常系数优化算法的有效性, 设计一个均匀双相介质模型(如图 7所示), 其参数见表 1.网格数300×300, 网格大小10 m×10 m, 时间采样间隔为Δt=1 ms, 选用20 Hz的雷克子波作为震源, 震源激发点位于模型中央, 接收点位于(x=1.25 km, z=1.35 km)处.为了简单起见, 本文只给出流相z分量的数值模拟结果.本文运用PML(完全匹配层)边界条件处理来自模型边界的地震波反射.

图 7 均匀介质模型 Fig. 7 Homogeneous medium model
表 1 均匀介质模型相关参数 Table 1 Parameters of the homogeneous medium model

图 8给出了不同方法得到的流相z分量t=0.8 s时的波场快照.可以看出采用十阶TE方法得到的波场快照(图 8b)中存在明显的数值频散现象(如白色箭头所指); 对比LS优化算法得到的波场快照(图 8c), 基于MN算法的常系数优化有限差分方法(图 8d)可以取得同样的压制数值频散效果.图 9是位于(x=1.25 km, z=1.35 km)的接收点得到的地震记录波形对比图.可以看出相同阶数时(2M=10), 优化有限差分算法(LS和MN)的数值模拟精度明显优于传统TE方法, 它们可以取得与高阶TE方法(2M=60)相近的精度, 红色虚线方框内的慢纵波(图 9b所示)尤为明显.在实际数值模拟中, 在较低差分阶数情况下, 两种优化有限差分方法计算量大致相当, 但随着阶数的增大, LS方法的计算量会逐渐大于MN方法.

图 8 不同方法得到的流相z分量t=0.8 s时的波场快照 (a) TE方法,2M=60; (b) TE方法,2M=10; (c) LS方法,2M=10; (d) MN方法,2M=10. Fig. 8 Snapshots of the z component of the fluid phase obtained by different methods (t=0.8 s) (a) TE method, 2M=60; (b) TE method, 2M=10; (c) LS method, 2M=10; (d) MN method, 2M=10.
图 9 三种方法(TE、LS和MN)的波形对比图 (a)原始波形图; (b)慢纵波波形放大显示图. Fig. 9 Waveforms comparison of the three methods (TE, LS and MN) (a) Original waveform; (b) Slow P-wave amplification diagram.
2.2 河道模型数值模拟

为了测试本文提出的变阶数优化有限差分方法在非均匀双相介质中的适应能力, 设计一个河道地质模型(如图 10所示), 其包括六个水平层和一个低速的河道模型, 相关模型参数如表 2所示(参考Dai, 1993).模型网格点个数400×300, 网格大小10 m×10 m, 时间采样间隔为t=1 ms.采用雷克子波作为震源子波(在图 10中用星号表示), 其位于模型(x=2 km, z=0.5 km)处.400个接收排列组合位于z=20 m处.

图 10 地质模型(这里只给出快纵波速度模型) Fig. 10 Geological model (only fast P-wave velocity model is given here)
表 2 河道模型相关参数 Table 2 Parameters of the river channel model

运用公式(27)可以获得河道模型中不同速度体对应的M值(如表 3所示).图 11a11b11c11d11e分别为高阶TE方法、低阶TE方法、LS方法、MN方法和变阶数MN(V-MN)方法得到的t=1.8 s时流相z分量的波场快照.图 12a12b12c12b12e是分别对应的单炮地震记录.对比图 11a12a图 11b12b, 可以看出低阶TE方法(M=8)存在明显的数值频散现象; 优化有限差分方法(图 11c12c图 11d12d)和变阶数优化有限差分方法(图 11e图 12e)能取得相同精度的数值频散压制效果.为了更清楚的对比不同方法数值模拟后的波形变化特征, 从图 12a12b12c12b12e中分别提取第200道的部分波场记录作对比显示(如图 13所示), 结果表明本文提出的变阶数优化有限差分方法(2M≤16)能取得与高阶TE方法(2M=60)相近的数值模拟精度, 且能取得与固定阶数(2M=16)优化有限差分方法(LS和MN)类似的压制频散效果. 图 14是对图 12不同方法的计算时间测试结果, 对于双相介质波动方程而言, 比较其他方法, 采用V-MN有限差分方法数值模拟更为高效.

表 3 地质模型不同速度体对应的M Table 3 M values for different velocities of the geological model
图 11 流相z分量t=1.8 s时的波场快照 (a) TE方法,2M=60; (b) TE方法,2M=16; (c) LS方法,2M=16; (d) MN方法,2M=16; (e) V-MN方法. Fig. 11 Snapshots of the z component of the fluid phase (t=1.8 s) (a) TE method, 2M=60; (b) TE method, 2M=16; (c) LS method, 2M=16; (d) MN method, 2M=16; (e) V-MN method.
图 12 流相z分量的合成地震记录 (a) TE方法,2M=60; (b) TE方法,2M=16; (c) LS方法,2M=16; (d) MN方法,2M=16; (e) V-MN方法. Fig. 12 Seismograms of the z components of the fluid phase (a) TE method, 2M=60; (b) TE method, 2M=16; (c) LS method, 2M=16; (d) MN method, 2M=16; (e) V-MN method.
图 13 对比图 12不同方法得到的单道合成地震记录 第1、2、3、4和5道分别截取于图 12(a)(b)(c)(d)(e)中的第200道. Fig. 13 Comparison of traces obtained by different methods from Fig. 12 Traces 1, 2, 3, 4 and 5 with x coordinate of 2000 m are some traces from Fig. 12 (a), (b), (c), (d) and (e), respectively.
图 14 测试时间对比 Fig. 14 Computation time test
3 讨论

针对双相介质的波动方程稳定性条件可以表示为

(28)

一般而言, 稳定性条件会随有限差分算子的差分阶数增加而变的更加严格(Liu and Sen, 2011).所以, 当常规固定阶数的有限差分算子满足稳定性条件时, 本文变阶数有限差分算子也必然满足稳定性条件, 因为变阶数有限差分方法中的差分阶数都不会大于常规固定阶数的有限差分方法.值得注意的是, 对于几乎所有的常系数优化有限差分方法而言, 尽管它们能一定程度上提高数值模拟的精度, 但是相比传统TE方法稳定性要求更为严格, 本文所使用的MN算法也不例外.

本文所使用的数值优化算法(MN)主要针对空间域的数值频散问题, 根据公式(12)和(13)可以用参数δ定义空间数值频散:

(29)

将不同优化算法(SA, TS, LS和MN)所得的有限差分系数代入上式就能进行空间数值频散分析.图 15是四种常系数优化算法与传统TE方法的空间频散曲线.在不加任何约束条件的情况下,可以看出四种优化方法在大波数区段都能取得相比传统TE方法更高的精度, 其中MN优化算法的大波数覆盖范围略高于其他三种方法.这些优化算法都存在一定的频散误差问题, 低波数范围内MN方法的频散误差较大, 但其在大波数区段的精度较其他方法高.在真正实际应用时, 考虑到各种优化方法的数值频散误差积累性, 还需要进一步对各种优化算法进行大规模模型和更长时程条件下的数值模拟精度测试, 以便进一步对比其优劣.

图 15 不同优化方法空间频散曲线对比 (a) 2M=10; (b) 2M=16. Fig. 15 Comparison of spatial dispersion curves of different optimization methods (a) 2M=10; (b) 2M=16.
4 结论

针对Biot一阶速度-应力方程组, 本文发展了一种变阶数优化有限差分数值模拟方案.将基于Remez迭代算法的二阶偏导数有限差分系数计算方法推广到一阶偏导数的差分系数求取中, 并结合交错网格有限差分格式和变阶数有限差分算法求解双相介质波动方程.采用Remez算法求取有限差分系数只需要几次循环迭代就能获得满意的结果, 其每次迭代过程中用于求解差分系数的线性方程组并不需要过多的计算量, 并且相比其他优化算法能覆盖更广的频带范围.对比可变阶数和固定阶数的频散误差变化曲线, 发现在一定频散误差阈值和频率范围内, 变阶数方法能取得和固定阶数有限差分方法相同的数值模拟精度.在数值实验方面, 相比传统泰勒展开方法, 变阶数优化有限差分数值模拟方案压制数值频散效果好, 计算效率更高, 适用于双相介质波动方程的正演模拟; 对比最小二乘优化算法, 本文提出的变阶数优化有限差分算法能在保证数值模拟精度的前提下进一步节省计算时间成本.

附录A

未知向量定义为

(A1)

其中(vx, vz)代表固相速度分量, (Vx, Vz)代表流相速度分量, (σxx, σzz, σxz)固相应力分量, s和流体压力p有关, 孔隙度ϕ定义如下

(A2)

公式(1)中的矩阵系数可以表示如下

(A3)

(A4)

(A5)

系数P, Q, RN表示Biot弹性常数, 公式(A3)—(A5)中的其他参数定义如下Rij=ρij/(ρ11ρ22-ρ122), B1=b(R11+R12), B2=b(R12+R22), , , , , ρ11ρ22分别表示固体和流体的质量系数, ρ12表示固体和流体之间的耦合质量系数, b与达西渗透系数κ有关:

(A6)

其中μ是孔隙流体的黏滞系数.

References
Alterman Z, Karal F C Jr. 1968. Propagation of elastic waves in layered media by finite difference methods. Bulletin of the Seismological Society of America, 58(1): 367-398.
Berryman J G. 1980. Confirmation of Biot's theory. Applied Physics Letters, 37(4): 382-384. DOI:10.1063/1.91951
Biot M A. 1956a. Theory of propagation of elastic waves in a fluid-saturated porous solid. Ⅰ. Low-frequency range. The Journal of the Acoustical Society of America, 28(2): 168-178. DOI:10.1121/1.1908239
Biot M A. 1956b. Theory of propagation of elastic waves in a fluid-saturated porous solid. Ⅱ. Higher frequency range. The Journal of the Acoustical Society of America, 28(2): 179-191. DOI:10.1121/1.1908241
Biot M A. 1962. Mechanics of deformation and acoustic propagation in porous media. Journal of Applied Physics, 33(4): 1482-1498. DOI:10.1063/1.1728759
Carcione J M, Helle H B. 1999. Numerical solution of the poroviscoelastic wave equation on a staggered mesh. Journal of Computational Physics, 154(2): 520-527. DOI:10.1006/jcph.1999.6321
Dablain M A. 1986. The application of high-order differencing to the scalar wave equation. Geophysics, 51(1): 54-66.
Dai N. 1993. Finite-difference simulation and imaging of seismic wave in complex media[Ph.D. thesis], Alberta: University of Alberta.
Dai N, Vafidis A, Kanasewich E R. 1995. Wave propagation in heterogeneous, porous media; a velocity-stress, finite-difference method. Geophysics, 60(2): 327-340.
Dong L G, Ma Z T, Cao J Z, et al. 2000. A staggered-grid high-order difference method of one-order elastic wave equation. Chinese Journal of Geophysics (in Chinese), 43(3): 411-419.
Gao J H, Zhang Y J. 2013. Staggered-grid finite difference method with variable-order accuracy for porous media. Mathematical Problems in Engineering, 2013: 157071. DOI:10.1155/2013/157071
He Z, Zhang J H, Yao Z X. 2019. Determining the optimal coefficients of the explicit finite-difference scheme using the Remez exchange algorithm. Geophysics, 84(3): S137-S147. DOI:10.1190/geo2018-0446.1
Holberg O. 1987. Computational aspects of the choice of operator and sampling interval for numerical differentiation in large-scale simulation of wave phenomena. Geophysical Prospecting, 35(6): 629-655. DOI:10.1111/j.1365-2478.1987.tb00841.x
Itzá R, Iturrarán-Viveros U, Parra J O. 2016. Optimal implicit 2-D finite differences to model wave propagation in poroelastic media. Geophysical Journal International, 206(2): 1111-1125. DOI:10.1093/gji/ggw180
Lemoine G I, Ou M Y, LeVeque R J. 2013. High-resolution finite volume modeling of wave propagation in orthotropic poroelastic media. SIAM Journal on Scientific Computing, 35(1): B176-B206. DOI:10.1137/120878720
Li L P, He B S. 2017. Elastic wave FCT finite difference numerical simulation in TTI media. Progress in Geophysics (in Chinese), 32(4): 1584-1590.
Li Z C, Yang F S, Wang X D. 2016. High-precision numerical simulation of first-order qP-waves in the transversely isotropic (TI) medium optimized by the LS-RSGFD method. Chinese Journal of Geophysics (in Chinese), 59(4): 1477-1490.
Liu C, Yang Q J, Lu Q, et al. 2014. Simulation of wave propagation in two-phase porous media using a frequency-space domain method. Chinese Journal of Geophysics (in Chinese), 57(9): 2885-5899. DOI:10.6038/cjg20140914
Liu C, Hu N, Guo Z Q, et al. 2018. Numerical simulation of the wavefield in a viscous fluid-saturated two-phase VTI medium based on the constant-Q viscoelastic constitutive relation with a fractional temporal derivative. Chinese Journal of Geophysics (in Chinese), 61(6): 2446-2458.
Liu X X, Yin X Y, Li H S. 2014. Optimal variable-grid finite-difference modeling for porous media. Journal of Geophysics and Engineering, 11(6): 065011. DOI:10.1088/1742-2132/11/6/065011
Liu Y, Wei X C. 2008. Finite-difference numerical modeling with even-order accuracy in two-phase anisotropic media. Applied Geophysics, 5(2): 107-114. DOI:10.1007/s11770-008-0014-6
Liu Y, Sen M K. 2011. Finite-difference modeling with adaptive variable-length spatial operators. Geophysics, 76(4): T79-T89. DOI:10.1190/1.3587223
Liu Y. 2013. Globally optimal finite-difference schemes based on least squares. Geophysics, 78(4): T113-T132. DOI:10.1190/geo2012-0480.1
Liu Y. 2014. Research progress on time-space domain finite-difference numerical solution and absorption boundary conditions of wave equation. Oil Geophysical Prospecting (in Chinese), 49(1): 35-46.
Liu Y. 2014. Optimal staggered-grid finite-difference schemes based on least-squares for wave equation modelling. Geophysical Journal International, 197(2): 1033-1047. DOI:10.1093/gji/ggu032
Masson Y J, Pride S R, Nihei K T. 2006. Finite difference modeling of Biot's poroelastic equations at seismic frequencies. Journal of Geophysical Research:Solid Earth, 111(B10): B10305. DOI:10.1029/2006JB004366
O'Brien G S. 2010. 3D rotated and standard staggered finite-difference solutions to Biot's poroelastic wave equations:stability condition and dispersion analysis. Geophysics, 75(4): T111-T119. DOI:10.1190/1.3432759
Ozdenvar T, McMechan G A. 1997. Algorithms for staggered-grid computations for poroelastic, elastic, acoustic, and scalar wave equations. Geophysical Prospecting, 45(3): 403-420. DOI:10.1046/j.1365-2478.1997.390275.x
Pei Z L. 2006. Staggering grid high-order finite-difference modeling of elastic wave transmission in biphase anisotropic medium. Oil Geophysical Prospecting (in Chinese), 41(2): 137-143.
Plona T J. 1980. Observation of a second bulk compressional wave in a porous medium at ultrasonic frequencies. Applied Physics Letters, 36(4): 259-261. DOI:10.1063/1.91445
Remez E Y. 1934. Sur un procédé convergent d'approximation successives pour de éterminer les polynoômes d'approximation. Comptes Rendus de l'Acadé mie des Sciences, 198: 2063-2065.
Robertsson J O A, Blanch J O, Symes W W. 1994. Viscoelastic finite-difference modeling. Geophysics, 59(9): 1444-1456. DOI:10.1190/1.1443701
Saenger E H, Bohlen T. 2004. Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid. Geophysics, 69(2): 583-591. DOI:10.1190/1.1707078
Sheen D H, Tuncay K, Baag C E, et al. 2006. Parallel implementation of a velocity-stress staggered-grid finite-difference method for 2-D poroelastic wave propagation. Computers & Geosciences, 32(8): 1182-1191.
Sidler R, Carcione J M, Holliger K. 2010. Simulation of surface waves in porous media. Geophysical Journal International, 183(2): 820-832. DOI:10.1111/j.1365-246X.2010.04725.x
Xin W, Yan Z C, Liang W Q, et al. 2015. Methods to determine the finite difference coefficients for elastic wave equation modeling. Chinese Journal of Geophysics (in Chinese), 58(7): 2486-2495. DOI:10.6038/cjg20150724
Xu S G, Liu Y. 2018. 3D acoustic and elastic VTI modeling with optimal finite-difference schemes and hybrid absorbing boundary conditions. Chinese Journal of Geophysics (in Chinese), 61(7): 2950-2968.
Yan H Y, Yang L, Li X Y. 2016. Optimal staggered-grid finite-difference schemes by combining Taylor-series expansion and sampling approximation for wave equation modeling. Journal of Computational Physics, 326: 913-930. DOI:10.1016/j.jcp.2016.09.019
Yang K D, Song G J, Li J S. 2011. FCT compact difference simulation of wave propagation based on the Biot and the squirt-flow coupling interaction (in Chinese), 54(5): 1348-1357.doi: 10.3969/j.issn.0001-5733.2011.05.024.
Yang K D, Yang D H, Wang S Q. 2002. Wave-field simulation based on the Biot-Squirt equation. Chinese Journal of Geophysics (in Chinese), 45(6): 853-861.
Yang L, Yan H Y, Liu H. 2015. Optimal rotated staggered-grid finite-difference schemes for elastic wave modeling in TTI media. Journal of Applied Geophysics, 122: 40-52. DOI:10.1016/j.jappgeo.2015.08.007
Yang L, Yan H Y, Liu H. 2016. Optimal implicit staggered-grid finite-difference schemes based on the sampling approximation method for seismic modeling. Geophysical Prospecting, 64(3): 595-610. DOI:10.1111/1365-2478.12325
Yang L, Yan H Y, Liu H. 2017. Optimal staggered-grid finite-difference schemes based on the minimax approximation method with the Remez algorithm. Geophysics, 82(1): T27-T42.
Yao G, Wu D, Debens H A. 2016. Adaptive finite difference for seismic wavefield modelling in acoustic media. Scientific Reports, 6: 30302. DOI:10.1038/srep30302
Zhang H X, Yin X Y, Zong Z Y. 2018a. Rock moduli estimation of inhomogeneous two-phase media with finite difference modeling algorithm. Journal of Geophysics and Engineering, 15(4): 1517-1527. DOI:10.1088/1742-2140/aaa045
Zhang J, Wu G C. 2015. Seismic modeling method of two phase medium utilizing staggered grid with optimal difference coefficients. Progress in Geophysics (in Chinese), 30(5): 2301-2311.
Zhang J B, Yang D H, He X J, et al. 2018. Discontinuous Galerkin method for solving wave equations in two-phase and viscoelastic media. Chinese Journal of Geophysics (in Chinese), 61(3): 926-937.
Zhang J H, Yao Z X. 2013a. Optimized finite-difference operator for broadband seismic wave modeling. Geophysics, 78(1): A13-A18.
Zhang J H, Yao Z X. 2013b. Optimized explicit finite-difference schemes for spatial derivatives using maximum norm. Journal of Computational Physics, 250: 511-526. DOI:10.1016/j.jcp.2013.04.029
Zhang Y J, Gao J H. 2014a. A 3D staggered-grid finite difference scheme for poroelastic wave equation. Journal of Applied Geophysics, 109: 281-291. DOI:10.1016/j.jappgeo.2014.08.007
Zhang Y J, Gao J H. 2014b. Time-space domain high-order staggered-grid finite difference method for porous media. Journal of Porous Media, 17(9): 785-795. DOI:10.1615/JPorMedia.v17.i9.30
Zhang Y J, Gao J H, Peng J G. 2018b. Variable-order finite difference scheme for numerical simulation in 3-D poroelastic media. IEEE Transactions on Geoscience and Remote Sensing, 56(5): 2991-3001. DOI:10.1109/TGRS.2017.2789159
Zhou Z S, Tang L. 2012. Seismic modeling and numerical dispersion correcting in double-phase medium based on reformulated BISQ model. Journal of Central South University (Science and Technology) (in Chinese), 43(4): 1411-1418.
Zhu X, McMechan G A. 1991. Numerical simulation of seismic responses of poroelastic reservoirs using Biot theory. Geophysics, 56(3): 328-339. DOI:10.1190/1.1443047
董良国, 马在田, 曹景忠, 等. 2000. 一阶弹性波方程交错网格高阶差分解法. 地球物理学报, 43(3): 411-419. DOI:10.3321/j.issn:0001-5733.2000.03.015
李立平, 何兵寿. 2017. TTI介质弹性波FCT有限差分数值模拟. 地球物理学进展, 32(4): 1584-1590. DOI:10.6038/pg20170423
李振春, 杨富森, 王小丹. 2016. 基于LS-RSGFD方法优化的横向各向同性(TI)介质一阶qP波高精度数值模拟. 地球物理学报, 59(4): 1477-1490. DOI:10.6038/cjg20160428
刘财, 杨庆节, 鹿琪, 等. 2014. 双相介质中地震波的频率-空间域数值模拟. 地球物理学报, 57(9): 2885-2899. DOI:10.6038/cjg20140914
刘财, 胡宁, 郭智奇, 等. 2018. 基于分数阶时间导数常Q黏弹本构关系的含黏滞流体双相VTI介质中波场数值模拟. 地球物理学报, 61(6): 2446-2458. DOI:10.6038/cjg2018M0063
刘洋. 2014. 波动方程时空域有限差分数值解及吸收边界条件研究进展. 石油地球物理勘探, 49(1): 35-46.
裴正林. 2006. 双相各向异性介质弹性波传播交错网格高阶有限差分法模拟. 石油地球物理勘探, 41(2): 137-143. DOI:10.3321/j.issn:1000-7210.2006.02.005
辛唯, 闫子超, 梁文全, 等. 2015. 用于弹性波方程数值模拟的有限差分系数确定方法. 地球物理学报, 58(7): 2486-2495. DOI:10.6038/cjg20150724
徐世刚, 刘洋. 2018. 基于优化有限差分和混合吸收边界条件的三维VTI介质声波和弹性波数值模拟. 地球物理学报, 61(7): 2950-2968. DOI:10.6038/cjg2018L0250
杨宽德, 宋国杰, 李静爽. 2011. Biot流动和喷射流动耦合作用下波传播德FCT紧致差分模拟. 地球物理学报, 54(5): 1348-1357. DOI:10.3969/j.issn.0001-5733.2011.05.024
杨宽德, 杨顶辉, 王书强. 2002. 基于Biot-Squirt方程的波场模拟. 地球物理学报, 45(6): 853-861. DOI:10.3321/j.issn:0001-5733.2002.06.013
张杰, 吴国忱. 2015. 基于优化差分系数的双相介质交错网格正演模拟. 地球物理学进展, 30(5): 2301-2311. DOI:10.6038/pg20150543
张金波, 杨顶辉, 贺茜君, 等. 2018. 求解双相和黏弹性介质波传播方程的间断有限元方法及其波场模拟. 地球物理学报, 61(3): 926-937. DOI:10.6038/cjg2018L0095
周竹生, 唐磊. 2012. 改进BISQ模型的双相介质地震波场数值模拟及频散校正. 中南大学学报(自然科学版), 43(4): 1411-1418.