﻿ 起伏地表下基于改进BISQ模型双相介质中曲线网格有限差分法波场模拟
 地球物理学报  2018, Vol. 61 Issue (8): 3356-3373 PDF

1. 长安大学地质工程与测绘学院地球物理系, 西安 710054;
2. 长安大学计算地球物理研究所, 西安 710054;
3. 阿联酋哈利法科技大学石油学院, 阿布扎比 2533

Wavefield modeling in two-phase media including undulated topography based on reformulated BISQ model by Curvilinear Grid FD method
YANG ShangBei1, BAI ChaoYing1,2, ZHOU Bing3
1. Dept. of Geophysics, School of Geology Engineering and Geomatics, Chang'an University, Xi'an 710054, China;
2. Institute of Computational Geophysics, Chang'an University, Xi'an 710054, China;
3. The Petroleum Institute, Khalifa University of Science and Technology, Po BOX 2533, Abu Dhabi, UAE
Abstract: In this article, the corresponding equations in curvilinear coordinate system based on the first-order velocity stress equation of two-phase isotropic media based on the improved BISQ model are derived, then the equations are numerically solved by an optimized high-order non-staggered finite difference scheme, that is, DRP/opt MacCormack scheme. To implement the undulating free-surface topography, we derive the analytical relationship between derivatives of velocity components and use the compact finite-difference scheme and traction-image method. The curvilinear grid finite difference method uses body-conforming grid to describe the free surface, thus avoids the numerical approximation caused by scattering. In the undulating free surface and the undulating interface of two-phase medium, the complex reflection wave and transmission wave can be clearly demonstrated by the numerical simulation results. The simulation results show that the curvilinear grid finite-difference method can accurately solve the propagation problem of seismic waves in a two-phase isotropic medium with undulating surface.
Key words: Undulating surface    Reformulated BISQ mechanism    Porous isotropic media    Curvilinear grid    Finite difference method
0 引言

1 基于改进BISQ模型双相介质中曲线网格有限差分算法 1.1 曲线坐标系与直角坐标系的转换关系

 图 1 曲线网格和计算网格的网格变换关系 Fig. 1 Transform a boundary-conforming grid into a uniform grid

 (1)

 (2)

 (3)

 (4)

1.2 曲线坐标系下基于改进BISQ模型各向同性双相介质的一阶速度-应力方程

 (5)

λμ为固体骨架的拉梅常数，φ表示孔隙度，ρsρf分别是固体和流体的密度，ρa表示固流耦合附加质量密度，∂为有效应力之孔隙弹性系数，KsKf分别表示固体骨架和固体颗粒的体积模量，F表示Biot流动系数密度，η表示孔隙流体的黏滞系数，k表示孔隙介质的渗透率.

 (6a)

 (6b)

 (6c)

 (6d)

 (6e)

 (6f)

 (6g)

 (6h)

 (7)

 (8)

 (9)

 (10)

 (11)

 (12)

1.3 空间和时间差分格式

Hixon和Turkel(2000)采用Tam和Webb(1993)提出的DRP(Dispersion Relation Preserving)的方法对MacCormack格式的系数进行了优化，得到了低频散、低耗散的DRP/opt MacCormack格式.DRP/opt MacCormack是一种高阶优化的MacCormack类有限差分格式，相对于传统的2~4阶MacCormack格式，它能够更好地抑制数值频散, 祝贺君等(2009)将此方法应用于二维各向异性介质的地震波场数值模拟中，验证了DRP/opt MacCormack差分格式既保留了同位网格的优势又避免了它的缺点, 是一种有效地研究复杂各向异性介质中地震波传播规律的工具.本文采用DRP/opt MacCormack有限差分算子来离散方程.

 (13)

 (14)

 (15)

 (16)

 (17)

 (18)

 (19)

 (20)

1.4 自由地表边界条件

 (21)

 (22)

 (23)

 (24)

 (25)

 (26)

 (27)

k=k0-1和k=k0-2网格层不满足DRP/opt MacCormack有限差分算子的使用条件，为使边界处的差分格式精度与内部高阶格式相匹配，k=k0-1和k=k0-2网格层的速度关于η轴的导数由单边4/4阶紧致MacCormack差分算子计算得到：

 (28)

1.5 震源和人工边界条件

 (29)

 (30)

2 数值模拟算例及分析 2.1 地表水平半空间模型和水平两层模型

 图 2 模型网格剖分示意图 Fig. 2 Mesh generation of models

 图 3 地表水平半空间模型0.35 s时刻的波场快照 (a)流相Vx; (b)流相Vz; (c)固相Vx; (d)固相Vz分量. Fig. 3 Wavefield snapshots at time of 0.35s for half-space model with horizontal surface (a) Vx in fluid phase; (b) Vz in fluid phase; (c) Vx in solid phase; (d) Vz in solid phase.

 图 4 水平两层模型0.30s时刻的波场快照 (a)流相Vx; (b)流相Vz; (c)固相Vx; (d)固相Vz分量.图中黑色实线为界面位置, 黑星为震源位置. Fig. 4 Wavefield snapshots at time of 0.30s for horizontal two-layer model (a) Vx in fluid phase; (b) Vz in fluid phase; (c) Vx in solid phase; (d) Vz component in solid phase. In figure the black lines and star represent the subsurface interface and the source location, respectively.
2.2 曲线网格与矩形网格方法模拟结果对比

 图 5 含起伏地表半空间模型0.30 s时刻的波场快照 (a)曲线网格流相Vz; (b)矩形网格流相Vz; (c)曲线网格固相Vz; (d)矩形网格固相Vz分量. Fig. 5 Wavefield snapshots at time of 0.30s for half-space model with undulating surface (a) Curvilinear grid Vz in fluid phase; (b) Rectangular grid Vz in fluid phase; (c) Curvilinear grid Vz in solid phase; (d) Rectangular grid Vz component in solid phase.
 图 6 含起伏地表半空间模型的单炮地震记录 (a)曲线网格流相Vz; (b)矩形网格流相Vz; (c)曲线网格固相Vz; (d)矩形网格固相Vz分量. Fig. 6 Common source point gather seismic sections for half-space model with undulating surface (a) Curvilinear grid Vz in fluid phase; (b) Rectangular grid Vz in fluid phase; (c) Curvilinear grid Vz in solid phase; (d) Rectangular grid Vz component in solid phase.
 图 7 含起伏地表半空间模型的单道合成理论地震图(检波器位置：1278 m，-250 m) (a)固相Vz分量; (b)液相Vz分量.图中实线、虚线分别为曲线、矩形网格方法模拟结果. Fig. 7 Synthetic seismogram for half-space model with undulating surface (receiver location:1278 m, -250 m) (a) Vz in solid phase; (b) Vz in fluid phase. In figure the solid and dotted lines indicate the results of curvilinear grid and rectangular grid, respectively.
2.3 含起伏地表半空间模型和起伏两层模型

 图 8 含自由地表水平半空间模型0.30 s时刻的波场快照(图中黑色星为震源位置) (a)流相Vx; (b)流相Vz; (c)固相Vx; (d)固相Vz. Fig. 8 Wavefield snapshots at time of 0.30 s for half-space model with horizontal free surface (In figure black star indicates source location) (a) Vx in fluid phase; (b) Vz in fluid phase; (c) Vx in solid phase; (d) Vz component in solid phase.
 图 9 含自由地表水平半空间模型的单炮地震记录 (a)流相Vx; (b)流相Vz; (c)固相Vx; (d)固相Vz分量. Fig. 9 Common source point gather seismic sections for half-space model with horizontal surface (a) Vx in fluid phase; (b) Vz in fluid phase; (c) Vx in solid phase; (d) Vz component in solid phase.

 图 10 含起伏地表半空间模型0.30 s时刻的波场快照(图中黑色星为震源位置) (a)流相Vx; (b)流相Vz; (c)固相Vx; (d)固相Vz分量. Fig. 10 Wavefield snapshots at time of 0.30 s for half-space model with undulating surface (In figure the black star indicates the source location) (a) Vx in fluid phase; (b) Vz in fluid phase; (c) Vx in solid phase; (d) Vz component in solid phase.
 图 11 含起伏地表半空间模型的单炮地震记录 (a)流相Vx; (b)流相Vz; (c)固相Vx; (d)固相Vz分量. Fig. 11 Common source point gather seismic sections for half-space model with undulating surface (a) Vx in fluid phase; (b) Vz in fluid phase; (c) Vx in solid phase; (d) Vz component in solid phase.

 图 12 起伏两层模型流相Vz分量不同时刻的波场快照(图中黑色星为震源位置) (a) 0.15 s; (b) 0.20 s; (c) 0.25 s; (d) 0.30 s. Fig. 12 Wavefield snapshots of Vz component in fluid phase at different time steps for undulated two-layered model (In figure black star indicates the source location)
 图 13 起伏两层模型中的单炮地震记录 (a)流相Vx; (b)流相Vz; (c)固相Vx; (d)固相Vz分量. Fig. 13 Common source point gather seismic sections for undulated two-layered model (a) Vx in fluid phase; (b) Vz in fluid phase; (c) Vx in solid phase; (d) Vz component in solid phase.
3 结果与讨论

(1) 通过与前人结果的对比分析，验证了基于改进的BISQ模型双相各向同性介质中采用曲线网格有限差分法进行波场数值模拟的正确性；

(2) 当流体为非黏滞性时，固相和流相波场中均可以观察到快纵波、横波和慢纵波；在模型介质分界面和自由地表处，每种波都会发生反射(或反射转换)、或透射(或透射转换)，从而产生多种波，使得波场信息非常丰富；通过分析各种波的传播特征，可以帮助我们加深对双相介质地震响应的认识；

(3) 当存在起伏地表或起伏介质分界面时，曲线网格有限差分法采用贴体网格来剖分模型，从而避免矩形网格阶梯状近似带来的数值频散；对于研究地震波在包含起伏地表或起伏界面的双相各向同性介质中的传播问题，曲线网格有限差法是一种有效且精确的方法.但是，当介质包含断层、地层尖灭等复杂地质构造时，曲线网格有限差法无法用贴体网格来精确的剖分模型，此方法不再适用.

(4) 曲线网格有限差法可以直接推广到三维各向同性双相介质，但是对于各向异性双相介质，由于波动方程的复杂性，本文施加自由地表的方法不再适用.

References
 Biot M A. 1956a. Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low-frequency range.. J. Acoust. Soc. Am., 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. J. Acoust. Soc. Am., 28(2): 179-191. Diallo M S, Appel E. 2000. Acoustic wave propagation in saturated porous media:reformulation of the Biot/Squirt flow theory. J. Appl. Geophys., 44(4): 313-325. DOI:10.1016/S0926-9851(00)00009-4 Diallo MS, Prasad M, Appel E. 2003. Comparison between experimental results and theoretical predictions for P-wave velocity and attenuation at ultrasonic frequency. Wave Motion, 37(1): 1-16. DOI:10.1016/S0165-2125(02)00018-5 Dvorkin J, Nur A. 1993. Dynamic poroelasticity:a unified modelwith the squirt and the Biot mechanisms. Geophysics, 58(4): 524-533. DOI:10.1190/1.1443435 Gist G A. 1994. Interpreting laboratory velocity measurements in partially gas-saturated rocks. Geophysics, 59(7): 1100-1109. DOI:10.1190/1.1443666 Gurevich B, Lopatnikov S L. 1995. Velocity and attenuation of elastic waves in finely layered porous rocks. Geophys. J. Int., 121(3): 933-947. DOI:10.1111/gji.1995.121.issue-3 Han C, Chen X F, Jiang D, et al. 2007. Application of the improved BISQ model to seismic numerical simulation of two-phase media. Natural Gas Industry (in Chinese), 27(10): 49-52. Hayashi K, Bums D R, Toksöz M N. 2001. Discontinuous-grid finite-difference seismic modeling including surface topography. Bull. Seismol. Soc. Am., 91(6): 1750-1764. DOI:10.1785/0120000024 Hestholm S O, Ruud B O, Husebye E S. 1999. 3-D versus 2-D finite-difference seismic synthetics including real surface topography. Phys. Earth Planet. Int., 113(1-4): 339-354. DOI:10.1016/S0031-9201(99)00010-2 Hestholm S O, Ruud B O. 2000. 2D finite-difference viscoelastic wave modelling including surface topography. Geophys. Prosp., 48(2): 341-373. DOI:10.1046/j.1365-2478.2000.00185.x Hestholm S O, Ruud B O. 1994. 2D finite-difference elastic wave modelling including surface topography. Geophys. Prosp., 42(5): 371-390. DOI:10.1111/gpr.1994.42.issue-5 Hixon R. 1997. Evaluation of a high-accuracy MacCormack-type scheme using benchmark problems. J. Comput. Acous., 6(3): 291-305. Hixon R, Turkel E. 2000. Compact implicit MacCormack-type schemes with high accuracy. J. Comput. Phys., 158(1): 51-70. Jiang L L, Sun J G. 2008. Source terms of elliptic system in grid generation. Global Geology (in Chinese), 27(3): 298-305. Jiang L L. 2010. Body-fitted grid generation for the geological conditions[Ph. D. thesis]. Changchun: Jilin University. Jih R, McLaughlin K L, Der Zoltan A. 1988. Free-boundary conditions of arbitrary polygonal topography in a two-dimensional explicit elastic finite-difference scheme. Geophysics, 53(8): 1045-1055. DOI:10.1190/1.1442541 Jones T D. 1986. Pore fluids and frequency-dependent wave propagation in rocks. Geophysics, 56(10): 1939-1953. Lan H T, Liu C, Guo Z Q. 2014. Modelling of seismic wave propagation in two-phase medium based on reformulated BISQ model using time-splitting staggered pseudospectral method. Global Geology (in Chinese), 33(1): 190-199. Li N. 2012. Finite difference seismoelectric wave modeling in 2D porous media with surface topography[Ph. D. thesis] (in Chinese). Hefei: University of Science and Technology of China. Liu C, Lan H T, Guo Z Q, et al. 2003. Pseudo-spectral modeling and feature analysis of wave propagation in two-phase HTI medium based on reformulated BISQ mechanism. Chinese J. Geophys. (in Chinese), 56(10): 3461-3473. DOI:10.6038/cjg20131021 Mavko G M, Jizba D. 1991. Estimating grain-scale fluid effects on velocity dispersion in rocks. Geophysics, 56(12): 1940-1949. DOI:10.1190/1.1443005 Opršal I, Zahradník J. 1999. Elastic finite-difference method for irregular grids. Geophysics, 64(1): 240-250. DOI:10.1190/1.1444520 Peng C Z, Li C M, Wang M C. 2007. Seismic wave-field modeling based on the modefied BISQ model. Computing Techniques for Geophysical & Geochemical Exploration (in Chinese), 29(1): 15-18. Sun Y C, Zhang W, Chen X F. 2016. Seismic-wave modeling in the presence of surface topography in 2D general anisotropic media by a curvilinear grid finite-difference method. Bull. Seism. Soc. Am., 106(3): 1036-1054. DOI:10.1785/0120150285 Tam C K W, Webb J C. 1993. Dispersion-relation-preserving finite difference schemes for computational acoustics. J. Comput. Phys., 107(2): 262-281. DOI:10.1006/jcph.1993.1142 Tarrass I, Giraud L, Thore P. 2011. New curvilinear scheme for elastic wave propagation in presence of curved topography. Geophys. Prosp., 59(5): 889-906. DOI:10.1111/gpr.2011.59.issue-5 Tessmer E, Kosloff D, Behle A. 1992. Elastic wave propagation simulation in the presence of surface topography. Geophys. J. Int., 108(2): 621-632. DOI:10.1111/gji.1992.108.issue-2 Thompson J F, Warsi Z U A, Mastin C W. 1985. Numerical Grid Generation: Foundations and Applications. New York: Elsevier North-Holland, Inc. Winkler K W. 1985. Dispersion analysis of velocity and attenuation in Berea sandstone. Journal of Geophysical Research, 90(B8): 6793-6800. DOI:10.1029/JB090iB08p06793 Yang D H, Zhang Z J, Teng J W, et al. 2000. The study of two-phase anisotropy, questions and applied prospects. Progress in Geophysics (in Chinese), 15(2): 7-21. Yang D H. 2002. Finite Element method of the elastic wave equation and wavefield simulation in two-phase anisotropic media. Chinese J. Geophys. (in Chinese), 45(4): 575-583. Zhang W, Chen X F. 2006. Traction image method for irregular free surface boundaries in finite difference seismic wave simulation. Geophys. J. Int., 167(1): 337-353. DOI:10.1111/gji.2006.167.issue-1 Zhang W, Zhang Z G, Chen X F. 2012. Three-dimensional elastic wave numerical modelling in the presence of surface topography by a collocated-grid finite-difference method on curvilinear grids. Geophys. J. Int., 190(1): 358-378. DOI:10.1111/gji.2012.190.issue-1 Zhang Y, Xu Y X, Xia J H, et al. 2015. Characteristics and application of surface wave propagation in fluid-filled porous media. Chinese J. Geophys. (in Chinese), 58(8): 2759-2778. DOI:10.6038/cjg20150812 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 H J, Zhang W, Chen X F. 2009. Two-dimensional seismic wave simulation in anisotropic media by non-staggered finite difference method. Chinese Journal of Geophysics (in Chinese), 52(6): 1536-1546. DOI:10.3969/j.issn.0001-5733.2009.06.015 韩翀, 陈雪菲, 蒋东, 等. 2007. 改进BISQ模型在双相介质地震波数值模拟中的应用. 天然气工业, 27(10): 49-52. DOI:10.3321/j.issn:1000-0976.2007.10.013 蒋丽丽, 孙建国. 2008. 基于Poisson方程的曲网格生成技术. 世界地质, 27(3): 298-305. 蒋丽丽. 2010. 面向地质条件的贴体网格生成技术[博士论文]. 长春: 吉林大学. 兰慧田, 刘财, 郭智奇. 2014. 利用时间分裂的错格伪谱法模拟地震波在基于改进BISQ模型的双相介质中的传播. 世界地质, 33(1): 190-199. 李娜. 2012. 含起伏地表的二维双相介质中地震波及电磁波传播有限差分算法[博士论文]. 合肥: 中国科学技术大学. 刘财, 兰慧田, 郭智奇, 等. 2013. 基于改进BISQ机制的双相HTI介质波传播伪谱法模拟与特征分析. 地球物理学报, 56(10): 3461-3473. DOI:10.6038/cjg20131021 彭传正, 李才明, 王明春. 2007. 基于改进BISQ模型的地震波场数值模拟. 物探化探计算技术, 29(1): 15-18. 杨顶辉, 张中杰, 滕吉文, 等. 2000. 双相各向异性研究、问题与应用前景. 地球物理学进展, 15(2): 7-21. 杨顶辉. 2002. 双相各向异性介质中弹性波方程的有限元解法及波场模拟. 地球物理学报, 45(4): 575-583. 张煜, 徐义贤, 夏江海, 等. 2015. 含流体孔隙介质中面波的传播特性及应用. 地球物理学报, 58(8): 2759-2778. DOI:10.6038/cjg20150812 周竹生, 唐磊. 2012. 改进BISQ模型的双相介质地震波场数值模拟及频散校正. 中南大学学报:自然科学版, 43(4): 1411-1418. 祝贺君, 张伟, 陈晓非. 2009. 二维各向异性介质中地震波场的高阶同位网格有限差分模拟. 地球物理学报, 52(6): 1536-1546. DOI:10.3969/j.issn.0001-5733.2009.06.015