2. 青岛海洋科学与技术国家实验室海洋矿产资源评价与探测技术功能实验室, 青岛 266071
3. 中国海洋大学海底科学与探测技术教育部重点实验室, 青岛 266100
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
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) |
其中,vx和vz为速度分量,τxx、τzz、τxz为应力分量,ρ为介质密度,t为时间,x和z代表空间位置,cij(i,j=1, 3, 5) 为介质的弹性系数,它由VTI介质中的弹性系数经过Bond变换得来.二维情况下,TTI介质的弹性系数可以表示为
(2) |
其中,cv11、cv13、cv33、cv44为VTI介质中的弹性参数,θ为TTI介质的极化角.
2 TTI介质弹性波有限差分算法在图 1所示的交错网格空间内,根据一阶弹性波方程交错网格高阶差分解法(董良国等,2000a),结合完全匹配层(Perfectly Matched Layer, PML)吸收边界条件(Berenger,1994)得到时间2阶,空间2N阶的空间有限差分格式,限于文章篇幅,本文只给出vx分量的差分格式为
(3) |
其中,Δt为时间离散步长,Δx和Δz为空间离散步长,i和j分别是x和z方向的网格序号,k为时间序号,N为差分阶数的一半,CnN为差分系数(董良国等,2000a),dx和dz表示x方向和z方向的阻尼因子(韩令贺和何兵寿,2010).
经推导,在TTI介质中,差分格式的稳定性条件(董良国等,2000b)为
(4) |
通量校正传输是在求解一维流体中的冲击波前阵面极大的浓度梯度过程中提出的,该技术满足物理上的守恒定律.这就要求在构造交错网格差分格式时,网格交界面处的通量可以相互抵消,即网格交界面处从一个方向流入的通量等于反方向流出的通量.使用这种方式,抵消掉了构造差分格式过程中产生的截断误差,避免了误差的累积.FCT假设所有的极值点都是由于数值频散引起的,对所有的网格点进行漫射通量校正,然后使用反漫射通量校正消除漫射通量校正产生的误差,同时在计算过程中对反漫射项进行了修正,从而保证解的正确性.具体算法如下所示:
(1) 首先建立一阶速度-应力方程交错网格有限差分格式,由k时刻的速度分量vxk(i,j)、vzk(i,j)和应力分量τxxk(i,j)、τzzk(i,j)、τxzk(i,j)分别求取k+1时刻的速度分量和应力分量.为了节省计算时间,提高运算效率,通量校正传输只对速度分量vx、vz进行处理.根据交错网格差分格式可知,速度分量校正完成后,应力分量也会得到相应的校正.
(2) 漫射运算
第一步:以速度分量为例,计算k时刻漫射通量,公式为
(5) |
其中,η1为漫射因子,其值满足0 < η1 < 1,具体大小取决于差分过程数值频散情况.它可以是一个常数,也可以是一个算子.
第二步:使用漫射通量P、Q对速度分量vx进行修正,公式为
(6) |
(3) 反漫射运算
第一步:计算k+1时刻的漫射通量,公式为
(7) |
其中,η2为反漫射因子,其选取准则与η1一致,但在实际应用中,其值应略大于η1,这是因为反漫射运算不仅要补偿漫射运算造成的损失,而且还要补偿有限差分运算过程带来的振幅损失.
第二步:求解k+1时刻的反漫射通量,公式为
(8) |
第三步:利用反漫射通量修正,得到压制数值频散后的正确波场值,公式为
(9) |
其中:
(10) |
模型的大小水平方向为400 m,垂直方向为400 m,使用的网格大小为Δx=Δz=10.0 m,采样间隔Δt=0.5 ms,定义在VTI中的弹性矩阵为:C11=2.5×109,C13=9.0×109,C33=18.4×109,C44=5.6×109,对其进行二维Bond变换,得到在TTI介质中的弹性矩阵C11、C13、C33、C15、C35、C55,极化角为30°,震源采用主频为35 Hz的雷克子波,在模型中心用集中力源激发,采用时间二阶,空间二阶差分对上述模型做正演,得到波前快照如图 2所示.
使用FCT方法与有限差分相结合,取漫射参数η1=0.001,反漫射参数η2=0.002,得到同一时刻的波前快照如图 3所示.
通过比较FCT校正前后的波前快照,可以明显的看出:传统有限差分方法在空间步长较大的情况下,模拟波传播的过程出现了明显的伪波动,即数值频散现象;FCT方法压制了有限差分正演模拟过程中产生的数值频散,使波场更加平滑.可以看出,数值模拟的结果得到明显的改善,常规有限差分中产生数值频散得到较好的压制效果.
4.2 BP模型试算本文使用局部的BP模型进行正演模拟来检验算法的适用性.选用的模型大小横向为800 m,垂直深度为1200 m,模型表层为海水,底部盐丘向上渗入到岩层当中.图 4为使用Thomsen参数表示的BP模型.
设计如下观测系统:集中力源激发,震源使用主频为30 Hz的雷克子波,网格大小Δx=Δz=5 m, 炮点位置定义在(400 m, 0 m), 即在海水表面处激发.接收点位于海水表面,道间距为5 m.FCT漫射因子的值为0.002,反漫射因子的值为0.003.采用时间二阶,空间二阶精度的交错网格有限差分格式对上述模型进行正演,得到FCT校正前后的单炮记录如图 5所示.
为了更好的对比观察FCT的效果,以x分量为例,对正演得到的单炮地震记录重新抽道如图 6所示.通过比较可以看到FCT同样适用于复杂介质的正演模拟,由于空间网格过大引起的数值频散得到了明显压制,地震记录中同相轴清晰的显示出来,有效提高了地震记录的分辨率.
研究发现,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. |