地球物理学进展  2017, Vol. 32 Issue (4): 1584-1590   PDF    
TTI介质弹性波FCT有限差分数值模拟
李立平1,2,3, 何兵寿1,2,3     
1. 中国海洋大学海洋地球科学学院, 青岛 266100
2. 青岛海洋科学与技术国家实验室海洋矿产资源评价与探测技术功能实验室, 青岛 266071
3. 中国海洋大学海底科学与探测技术教育部重点实验室, 青岛 266100
摘要:TTI(Tilted Transversely Isotropic)各向异性是对地下岩石中广泛存在的规则发育的裂缝和层理的一种有效的弹性近似,基于TTI介质的地震波数值模拟技术是分析地震波在复杂各向异性介质中的传播机理的有效工具.同时,高精度的数值模拟算法也能为后续的逆时偏移技术提供重要的技术支撑.由于TTI介质中地震波方程的弹性参数众多且变化复杂,常规有限差分技术在解决TTI介质正演模拟问题时往往会产生严重的数值频散现象,降低了数值模拟精度.通量校正传输(Flux-Corrected Transport,FCT)技术能够有效地压制由空间离散产生的数值频散.本文将FCT技术用于TTI介质中弹性波方程的交错网格高阶精度差分正演,在数值模拟过程中通过对波场进行漫射和反漫射校正实现了空间网格频散的压制.模型模拟结果表明,与常规有限差分算法相比,本文算法能够有效的压制大网格条件下的数值频散,提高模拟精度.
关键词TTI介质    有限差分    交错网格    数值频散    通量校正传输    
Elastic wave FCT finite difference numerical simulation in TTI media
LI Li-ping1,2,3 , HE Bing-shou1,2,3     
1. College of Marine Geo-science Ocean University of China, Qingdao 266100, China
2. Evaluation and Detection Technology Laboratory of marine mineral resources, Qingdao National Laboratory for Marine Science and Technology, Qingdao 266071, China
3. Key Lab of Submarine Geosciences and Prospecting Techniques Ministry of Education, Qingdao 266100, China
Abstract: TTI anisotropy is an effectively elastic approximate to the rule growth fissure and stratification. Seismic wave numerical simulation based on the TTI media is an effective way to analysis propagation mechanisms of seismic wave in complex anisotropic media. High-precision numerical simulation algorithm provides important technical support to reverse time migration. Owing to the elastic parameter of TTI media elastic wave equation are numerous and complex, conventional finite-difference technique have seriously numerical dispersion phenomenon in the forward modeling of TTI media and decrease the numerical simulation accuracy. FCT can suppress numerical dispersion caused by space discrete. The paper applied FCT to the elastic wave staggered grid high-order accuracy difference forward modeling in TTI media, make diffusion and anti-diffusion to wave field in the process of numerical simulation and suppress the numerical dispersion caused by large grid. The result of model simulation shows that compared with conventional finite difference method, algorithm in this paper can suppress numerical dispersion in the large grid and improve simulation accuracy.
Key words: TTI media     finite difference     staggered grid     numerical dispersion     FCT    
0 引言

TTI各向异性介质在岩石中广泛存在,一般认为,它是由沉积岩周期性薄互层和定向排列的裂隙引起的(吴国忱等, 2010).TTI介质的弹性参数随着观测方向变化而发生变化,其弹性特征主要表现为:地震波速度是传播方向的函数、纵横波之间相互耦合以及横波发生分裂(吴国忱, 2006).研究TTI介质中地震波传播规律,对于认识各向异性介质弹性特征、地震波在复杂各向异性介质中的传播机理具有重要意义.

正演模拟是研究TTI介质中地震波传播规律、波场耦合关系以及空间变化特征的有效工具(裴正林和牟永光, 2004; 吴国忱, 2006).有限差分是目前地震波场正演模拟的重要手段,使用有限差分方法进行地震波场正演模拟时,用离散网格逼近连续介质,以差分算子代替微分算子,差分的空间网格间距过大和时间采样不合适、震源主频过高都会引起波形的畸变和重影,使波场发生数值频散.其本质是一种因离散化求解波动方程而产生的伪波动(杨顶辉和滕吉文, 1997).数值频散覆盖了地震波场的有效信息,降低波场数值模拟的精度和分辨率,并影响了后续地震资料处理的效果.

目前,有限差分计算中数值频散的压制思路主要有三种:

(1) 提高有限差分阶数(董良国和李培明, 2004; 吴国忱和王华忠, 2005).一般来说,差分阶数越高数值模拟精度越高.但提高差分阶数不仅会增加计算量,而且达到一定精度以后,差分阶数进一步提高时,数值频散压制效果的提升幅度变得十分缓慢.

(2) 减小网格的剖分步长(Liang et al., 2015).减小网格剖分步长可以明显的压制数值频散,但是网格步长的减少会增加该介质内网格点数,显著的增加计算量.

(3) 采用通量校正传输技术压制数值频散(成景旺等, 2011; 张省等, 2012).

通量校正传输方法最早起源于Boris和Book(1973)等人求解流体力学连续方程,随后被应用于声波方程求解,实现了声波方程正演模拟中数值频散的有效压制.Fei和Larner(1995)杨顶辉和滕吉文(1997)等人将FCT方法应用于弹性波的正演模拟,取得了较好的效果.成景旺等(2011)将这一方法应用于海上黏弹介质盐丘模型,实现了海上黏滞性介质高精度数值模拟.

TTI介质弹性波方程有限差分正演模拟过程中,时间和空间网格离散均使得介质的弹性参数发生变化,并引起数值频散(吴国忱,2006).由于TTI介质弹性参数众多且变化复杂,常规方法无论提高差分阶数还是减小网格剖分步长,都无法完全压制TTI介质有限差分正演模拟产生的数值频散.与以上在频散产生之前进行压制的方法不同,本文使用FCT方法对得到的含频散波场进行压制,对每一个时刻的波场进行漫射和反漫射校正,避免了误差的累积.同时,FCT方法因为适用于大网格剖分步长而减少了计算量,提高了计算效率.

本文研究二维TTI介质有限差分FCT数值频散压制方法,使用漫射通量和反漫射通量对差分计算的波场进行校正,并对反漫射通量进行修正.将FCT技术应用于TTI介质中,使用简单模型和BP模型检验效果,实现TTI介质高精度数值模拟.

1 TTI介质中的弹性波方程

二维TTI介质中,用速度-应力表示的一阶弹性波方程为

(1)

其中,vxvz为速度分量,τxxτzzτxz为应力分量,ρ为介质密度,t为时间,xz代表空间位置,cij(ij=1, 3, 5) 为介质的弹性系数,它由VTI介质中的弹性系数经过Bond变换得来.二维情况下,TTI介质的弹性系数可以表示为

(2)

其中,cv11cv13cv33cv44为VTI介质中的弹性参数,θ为TTI介质的极化角.

2 TTI介质弹性波有限差分算法

图 1所示的交错网格空间内,根据一阶弹性波方程交错网格高阶差分解法(董良国等,2000a),结合完全匹配层(Perfectly Matched Layer, PML)吸收边界条件(Berenger,1994)得到时间2阶,空间2N阶的空间有限差分格式,限于文章篇幅,本文只给出vx分量的差分格式为

图 1 交错网格示意图 Figure 1 Staggered-grid sketch map
(3)

其中,Δt为时间离散步长,Δx和Δz为空间离散步长,ij分别是xz方向的网格序号,k为时间序号,N为差分阶数的一半,CnN为差分系数(董良国等,2000a),dxdz表示x方向和z方向的阻尼因子(韩令贺和何兵寿,2010).

经推导,在TTI介质中,差分格式的稳定性条件(董良国等,2000b)为

(4)
3 数值频散FCT压制方法

通量校正传输是在求解一维流体中的冲击波前阵面极大的浓度梯度过程中提出的,该技术满足物理上的守恒定律.这就要求在构造交错网格差分格式时,网格交界面处的通量可以相互抵消,即网格交界面处从一个方向流入的通量等于反方向流出的通量.使用这种方式,抵消掉了构造差分格式过程中产生的截断误差,避免了误差的累积.FCT假设所有的极值点都是由于数值频散引起的,对所有的网格点进行漫射通量校正,然后使用反漫射通量校正消除漫射通量校正产生的误差,同时在计算过程中对反漫射项进行了修正,从而保证解的正确性.具体算法如下所示:

(1) 首先建立一阶速度-应力方程交错网格有限差分格式,由k时刻的速度分量vxk(ij)、vzk(ij)和应力分量τxxk(ij)、τzzk(ij)、τxzk(ij)分别求取k+1时刻的速度分量和应力分量.为了节省计算时间,提高运算效率,通量校正传输只对速度分量vxvz进行处理.根据交错网格差分格式可知,速度分量校正完成后,应力分量也会得到相应的校正.

(2) 漫射运算

第一步:以速度分量为例,计算k时刻漫射通量,公式为

(5)

其中,η1为漫射因子,其值满足0 < η1 < 1,具体大小取决于差分过程数值频散情况.它可以是一个常数,也可以是一个算子.

第二步:使用漫射通量PQ对速度分量vx进行修正,公式为

(6)

(3) 反漫射运算

第一步:计算k+1时刻的漫射通量,公式为

(7)

其中,η2为反漫射因子,其选取准则与η1一致,但在实际应用中,其值应略大于η1,这是因为反漫射运算不仅要补偿漫射运算造成的损失,而且还要补偿有限差分运算过程带来的振幅损失.

第二步:求解k+1时刻的反漫射通量,公式为

(8)

第三步:利用反漫射通量修正,得到压制数值频散后的正确波场值,公式为

(9)

其中:

(10)
4 数据模型试算 4.1 均匀TTI介质模型

模型的大小水平方向为400 m,垂直方向为400 m,使用的网格大小为Δxz=10.0 m,采样间隔Δt=0.5 ms,定义在VTI中的弹性矩阵为:C11=2.5×109C13=9.0×109C33=18.4×109C44=5.6×109,对其进行二维Bond变换,得到在TTI介质中的弹性矩阵C11C13C33C15C35C55,极化角为30°,震源采用主频为35 Hz的雷克子波,在模型中心用集中力源激发,采用时间二阶,空间二阶差分对上述模型做正演,得到波前快照如图 2所示.

图 2 FCT校正前波前快照 (a)为水平分量;(b)为垂直分量. Figure 2 Wavefront before FCT correction (a)Is horizontal component; (b)Is vertical component.

使用FCT方法与有限差分相结合,取漫射参数η1=0.001,反漫射参数η2=0.002,得到同一时刻的波前快照如图 3所示.

图 3 FCT校正后波前快照 (a)为水平分量;(b)为垂直分量. Figure 3 Wavefront after FCT correction (a)Horizontal component; (b)Vertical component.

通过比较FCT校正前后的波前快照,可以明显的看出:传统有限差分方法在空间步长较大的情况下,模拟波传播的过程出现了明显的伪波动,即数值频散现象;FCT方法压制了有限差分正演模拟过程中产生的数值频散,使波场更加平滑.可以看出,数值模拟的结果得到明显的改善,常规有限差分中产生数值频散得到较好的压制效果.

4.2 BP模型试算

本文使用局部的BP模型进行正演模拟来检验算法的适用性.选用的模型大小横向为800 m,垂直深度为1200 m,模型表层为海水,底部盐丘向上渗入到岩层当中.图 4为使用Thomsen参数表示的BP模型.

图 4 BP模型示意图 (a)δ分量; (b)ε分量; (c)γ分量; (d)vp0分量. Figure 4 Schematic diagram of BP model (a)δ component; (b)ε component; (c)γ component; (d)vp0 component.

设计如下观测系统:集中力源激发,震源使用主频为30 Hz的雷克子波,网格大小Δxz=5 m, 炮点位置定义在(400 m, 0 m), 即在海水表面处激发.接收点位于海水表面,道间距为5 m.FCT漫射因子的值为0.002,反漫射因子的值为0.003.采用时间二阶,空间二阶精度的交错网格有限差分格式对上述模型进行正演,得到FCT校正前后的单炮记录如图 5所示.

图 5 FCT校正前后单炮地震记录 (a)FCT校正前水平分量;(b)FCT校正前垂直分量;(c)FCT校正后水平分量;(d)FCT校正后垂直分量. Figure 5 Horizontal and vertical components of seismic record before and after FCT (a)Horizontal component before FCT correction; (b) Vertical component before FCT correction; (c)Horizontal component after FCT correction; (d) Vertical component after FCT correction.

为了更好的对比观察FCT的效果,以x分量为例,对正演得到的单炮地震记录重新抽道如图 6所示.通过比较可以看到FCT同样适用于复杂介质的正演模拟,由于空间网格过大引起的数值频散得到了明显压制,地震记录中同相轴清晰的显示出来,有效提高了地震记录的分辨率.

图 6 抽道后的地震记录 (a)FCT校正前的水平分量;(b)FCT校正后水平分量. Figure 6 Seismic record after line shot (a)Horizontal component before FCT correction; (b)Horizontal component after FCT correction.

研究发现,FCT算法当中漫射因子η1和反漫射因子η2的取值对数值模拟的结果影响很大,二者取值过大或者过小都会在一定程度上影响频散的压制效果.取值过大,平滑过于强烈,会使地震记录中的同相轴变粗,同时丢失许多有效信息.取值过小,则不能有效的压制数值频散.只有漫射因子η1和反漫射因子η2都在适当范围内取值时,FCT算法才可以既取得较好的频散压制的效果,又能保证较为真实的波场.

5 结论

FCT方法能够有效地压制各向异性介质有限差分正演模拟过程中由于网格剖分等原因造成的数值频散,可以在大空间网格情况下实现高精度的数值模拟.本文分别在均匀的TTI介质和BP模型进行了模型试算,验证了FCT的可行性.计算过程中,FCT既保证了较高的计算精度,同时该方法因适应较大的采样间隔而节省了计算量,从而具有较高的计算速度,实现了TTI介质中的高精度弹性波正演模拟.

致谢 感谢审稿专家提出的修改意见和编辑部的大力支持!
参考文献
[] Berenger J P. 1994. A perfectly matched layer for the absorption of electromagnetic waves[J]. Journal of computational Physics, 114(2): 185–200. DOI:10.1006/jcph.1994.1159
[] Boris J P, Book D L. 1973. Flux-corrected transport. I. SHASTA, a fluid transport algorithm that works[J]. Journal of Computational Physics, 11(1): 38–69. DOI:10.1016/0021-9991(73)90147-2
[] Cheng J W, Gu H M, Liu L, et al. 2011. Flux correction transport (FCT) finite difference forward modeling for marine viscoelastic media[J]. Journal of Oil and Gas Technology, 33(12): 83–87.
[] Dong L G, Ma Z T, Cao J Z, et al. 2000a. A staggered-grid high-order difference method of one-order elastic wave equation[J]. Chinese Journal of Geophysics, 43(3): 411–419.
[] Dong L G, Ma Z T, Cao J Z. 2000b. A study on stability of the staggered-grid high-order difference method of first-order elastic wave equation[J]. Chinese Journal of Geophysics, 43(6): 856–864. DOI:10.3321/j.issn:0001-5733.2000.06.015
[] Dong L G, Li P M. 2004. Dispersive problem in seismic wave propagation numerical modeling[J]. Natural Gas Industry, 24(6): 53–56.
[] Fei T, Larner K. 1995. Elimination of numerical dispersion in finite-difference modeling and migration by flux-corrected transport[J]. Geophysics, 60(6): 1830–1842. DOI:10.1190/1.1443915
[] Han L H, He B S. 2010. First-order Quasi-P wave equation forward modeling and boundary conditions in VTI medium[J]. Oil Geophysical Prospecting, 45(6): 819–825.
[] Liang W Q, Wang Y F, Yang C C. 2015. Comparison of numerical dispersion in acoustic finite-difference algorithms[J]. Exploration Geophysics, 46(2): 206–212. DOI:10.1071/EG13072
[] Pei Z L, Mu Y G. 2004. Numerical simulation of seismic wave propagation[J]. Progress in Geophysics, 19(4): 933–941. DOI:10.3969/j.issn.1004-2903.2004.04.038
[] Wu G C. 2006. Seismic Wave Propagation and Imaging in Anisotropy Media[M]. Dongying: Press of China University of Petroleum.
[] Wu G C, Liang K, Yin X Y. 2010. The analysis of phase velocity and polarization feature for elastic wave in TTI media[J]. Chinese Journal of Geophysics, 53(8): 1914–1923. DOI:10.3969/j.issn.0001-5733.2010.08.017
[] Wu G C, Wang H Z. 2005. Analysis of numerical dispersion in wave-field simulation[J]. Progress in Geophysics, 20(1): 58–65. DOI:10.3969/j.issn.1004-2903.2005.01.012
[] Yang D H, Teng J W. 1997. FCT finite difference modeling of three-component seismic records in anisotropic medium[J]. Oil Geophysical Prospecting, 32(2): 181–190.
[] Zhang S, He B S, Wang Y F. 2012. Staggered-grid FCT finite difference numerical simulation in VTI media[J]. Chinese Journal of Engineering Geophysics, 9(5): 565–571.
[] 成景旺, 顾汉明, 刘琳, 等. 2011. 海上粘弹介质FCT有限差分正演模拟[J]. 石油天然气学报, 33(12): 83–87. DOI:10.3969/j.issn.1000-9752.2011.12.017
[] 董良国, 马在田, 曹景忠, 等. 2000a. 一阶弹性波方程交错网格高阶差分解法[J]. 地球物理学报, 43(3): 411–419.
[] 董良国, 马在田, 曹景忠. 2000b. 一阶弹性波方程交错网格高阶差分解法稳定性研究[J]. 地球物理学报, 43(6): 856–864. DOI:10.3321/j.issn:0001-5733.2000.06.015
[] 董良国, 李培明. 2004. 地震波传播数值模拟中的频散问题[J]. 天然气工业, 24(6): 53–56.
[] 韩令贺, 何兵寿. 2010. VTI介质一阶准P波方程正演模拟及边界条件[J]. 石油地球物理勘探, 45(6): 819–825.
[] 裴正林, 牟永光. 2004. 地震波传播数值模拟[J]. 地球物理学进展, 19(4): 933–941. DOI:10.3969/j.issn.1004-2903.2004.04.038
[] 吴国忱. 2006. 各向异性介质地震波传播与成像[M]. 东营: 中国石油大学出版社.
[] 吴国忱, 梁锴, 印兴耀. 2010. TTI介质弹性波相速度与偏振特征分析[J]. 地球物理学报, 53(8): 1914–1923. DOI:10.3969/j.issn.0001-5733.2010.08.017
[] 吴国忱, 王华忠. 2005. 波场模拟中的数值频散分析与校正策略[J]. 地球物理学进展, 20(1): 58–65. DOI:10.3969/j.issn.1004-2903.2005.01.012
[] 杨顶辉, 滕吉文. 1997. 各向异性介质中三分量地震记录的FCT有限差分模拟[J]. 石油地球物理勘探, 32(2): 181–190.
[] 张省, 何兵寿, 王玉凤. 2012. VTI介质交错网格FCT有限差分数值模拟[J]. 工程地球物理学报, 9(5): 565–571.