Numerical Simulation of Seismic Wave Propagation Using Convolutional Differentiator
-
摘要: 为了解决地震波场数值模拟中的速度与精度的匹配及局部信息与全局信息的耦合等问题, 该文借助新推出的基于Forsyte广义正交多项式的褶积微分算子, 将计算数学中的Forsyte多项式褶积微分算子应用到地震波传播的数值模拟中.复杂非均匀介质模型数值模拟结果说明了该方法的可行性和优越性.该方法同时具有广义正交多项式褶积微分算子的高精度和有限差分短算子算法的高速度.通过对算子长度的调节及算子系数的优化, 可同时兼顾波场解的全局信息与局部信息.Abstract: To improve the accuracy and the efficiency of seismic wave simulation and to couple the local and global information better, this paper develops a novel modeling approach referring to the convolutional differentiator based on generalized Forsyte orthogonal polynomial, which applies optimal convolutional operators for spatial differentiation in wave equation. The numerical experiment of complex heterogeneous model demonstrates that the algorithm can bring reliable results with high precision and can be extended to seismic wave simulation in anisotropic media. This method is highly precise in generalized orthogonal polynomial convolutional differentiator and also highly efficient in finite difference short operator method. The local and global information can be considered at the same time by optimizing the coefficients of the operator and adjusting the operator length.
-
地震波场数值模拟算法历来都是地震学研究的一个热门领域.迄今为止, 地震波场正演模拟的数值算法已有很多.常用的地震波场数值模拟方法有射线追踪法、积分方程法、基于波动方程直接解法的有限差分法(Alterman and Karal, 1968; Magnier et al., 1994; Graves, 1996)、伪谱法(Kosloff and Baysal, 1982; Kosloff et al., 1984; Lin and Zhao, 2007)、有限元法(周辉等, 1997)和谱元法(Komatitsch and Tromp, 1999)等, 它们都是比较成熟、比较常用的传统方法.这些传统的方法中, 精度较好地球科学———中国地质大学学报第33卷的, 其计算速度和效率难以忍受, 而速度和效率尚可的, 其精度又满足不了要求.在数值计算中, 前人提出的褶积微分算子都利用傅立叶正反变换技术来构造, 而傅立叶方法是一种全局算法, 即它实际上使用了长度是空间采样点数的长褶积滤波器, 虽然它精度高, 占用内存小, 但难以很好地耦合局部信息(张中杰等, 1996; 滕吉文等, 1995).常规有限差分法由于使用了短算子而具有算法简单快速的特点, 但难以克服频散效应; 而要解决频散问题, 须加密数值计算的网格, 这势必导致计算量增加, 效率下降.
为了兼顾波场模拟的精度与效率, 本文以计算数学中的佛尔塞(Forsyte)多项式插值函数为基础(谢靖, 1981), 构造一个新的褶积微分算子, 并将之应用于弹性波传播的数值模拟中.该方法同时具有广义正交多项式迭积微分算子的高精度和有限差分短算子算法的高速度.通过对算子长度的调节及算子系数的优化, 可同时兼顾波场解的全局信息与局部信息.数值实验说明了该方法的可行性和优越性.
1. Forsyte多项式褶积微分算子
首先我们给出Forsyte多项式褶积微分算子. Forsyte多项式是一个广义正交多项式, 其插值函数可写为(谢靖, 1981):
(1) 其中:
f(xi)为被插值函数f(x)在点xi处的值.式(1)中P0(x), …, Pj+1(x)定义为Forsyte多项式系统.
对(1)式中的x求导, 可得:
(2) 其中:
Forsyte多项式微分算子可写为:
(3) 其中:
将上述微分算子(3)离散化可得:
(4) 这里, i为采样指标; Δx为沿着x轴的采样间隔.就实际应用而言, 为了兼顾插值精度和计算量, 须将微分算子截成短算子.这势必引起Gibbs现象.另外, 多项式的引入还将引起Runge现象.为了消除这些现象, 必须采用窗函数以截断长微分算子.本文采用的是下列Gaussian窗函数:
(5) 其中: mx为单边截断长度的采样数; c为常数; a (0.1≤ a≤ 0.75)为衰减因子.将微分算子(4)用式(5)截断并锯齿化后, 可得一阶褶积微分算子:
(6) 通过对算子长度的调节及算子系数的优化, 可同时兼顾波场解的全局信息与局部信息, 本文运用算子长度为9点的一阶褶积微分算子求解弹性波场的一阶速度-应力方程.通过试验选取的9点褶积算子的最优权系数为: - 0.002 34, 0.008 74, - 0. 027 4, 0. 085, 0. 0, - 0. 085, 0. 027 4, - 0.008 74, 0.002 34, 可以看出, 算子的系数是反对称的.图 1是9点褶积微分算子振幅随采样点的变化曲线, 可以看出, 9点算子振幅已经可以很快地衰减.
图 1 9点褶积微分算子振幅随采样点的变化曲线Fig. 1. Amplitude versus sampling point for nine point convolutional differentiator以f(x)=e-a2x2, 且a=1/Δx, Δx=10作为检验函数, 检验本文所构造的Forsyte广义正交多项式褶积微分算子的计算精度.图 2是基于Forsyte褶积微分算子法的解与基于错格伪谱微分算子法的解之间的精度比较.从图 2可以看出, 9点Forsyte褶积微分算子法的解比错格伪谱微分算子法的数值解要精确.因此, 9点Forsyte褶积微分算子法可以较精确地模拟地震波场.
2. 波动方程的褶积微分算子离散形式
对模型区间离散后, 设n, m, k分别是沿着空间x轴、z轴、时间t轴的采样点数, Δx、Δz、Δt分别是沿x轴、z轴、t轴的采样间隔, mx、mz是沿x轴、z轴采样数的半算子长度, 则二维非均匀各向同性介质弹性波一阶速度-应力方程, 其离散化的时间错格差分褶积微分算子法格式为(体力为零):
(7) 
(8) 
(9) 
(10) 
(11) 地震波数值模拟中常用的子波函数有Ricker子波和Gauss子波, 本文运用的是Ricker子波, 其震源时间函数形式为:
震源空间函数形式(俞寿朋, 1996)为:
其中, f0为Ricker子波的主频; (x0, z0)为震源中心位置; 常数α的取值与脉冲宽度有关, 它决定了震源脉冲在空间的衰减速度.
本文数值模拟试验中采用的边界条件是吸收边界(Cerjan et al., 1985).该方法需要在计算网格的边界上设置波场过渡带, 在过渡带内波场逐渐衰减至零, 这样在边界上将没有波场被反射.这需要在过渡带内设一个衰减因子:
(12) 上式中: β是边界振幅衰减系数或衰减率(attenuation rate); N是过渡带的网格数; i是过渡带的网格编号(1 < i < N); G是衰减因子且是上述几个量的函数.
3. 数值实验
3.1 Marmousi经典模型检验
Marmousi模型是检验地震波成像效果及估计速度精度的标准模型.该模型在构造上十分复杂, 模型的空间网格为384×122, 网格间距Δx=Δz= 10 m.利用基于Forsyte广义正交多项式褶积微分算子法对模型进行数值模拟, 采样间隔Δt=1ms, 震源坐标为(192, 20), Ricker子波主频为25 Hz.在原始Marmousi模型的纵波速度基础上获得其他介质参数, 即: 用
图 3 Marmousi速度模型及500 ms时的波场快照a.速度模型; b.x分量; c.z分量Fig. 3. Snapshots for the Marmousi velocity model at the time of 500 ms从Marmousi模型的波场快照来看, 震源发出的地震波在介质界面和突变点处产生的反射和绕射等现象非常明显, 基于Forsyte广义正交多项式褶积微分算子的数值模拟方法能很好地模拟地震波在速度变化比较剧烈的介质中的传播.
3.2 与伪谱法的模拟效果比较
我们利用一个倾斜界面模型对本文发展的方法和经典伪谱法作比较, 检验两种方法在刻画波场精细细节和抑制频散方面的效果.模型的几何结构如图 4所示.介质密度是2 650 kg/m3, 界面上方纵、横波速度分别为2 000 m/s与1 155 m/s, 界面下方纵、横波速度分别为3 000 m/s与1 730 m/s.网格大小为256×256, 网格间距为dx=dz=10 m, 时间步长为1 ms, 震源位置坐标为(128, 120), 子波主频为25 Hz.
图 5是利用本文发展的方法得到的400 ms时波场快照, 图 6是利用伪谱法得到的400 ms时波场快照.对比图 5与图 6可以看到, 本文发展的方法对数值频散的抑制能力跟伪谱法基本相当, 利用褶积算子法与伪谱法的计算时间比为M/(4log2N), 计算得到二者的计算时间比为9/32, 即, 其计算效率明显优于伪谱法.
图 5 Forsyte褶积微分算子法波场快照a.x分量; b.z分量Fig. 5. Snapshots obtained by Forsyte convolutional differentiator method4. 与有限差分法及伪谱法的比较
4.1 与有限差分法和伪谱法的精度比较
在检验波动方程基础上模拟方法准确性的最精确实验就是均匀介质中脉冲响应的计算.因为在数值计算中, 单脉冲能够产生最宽的频段.Zhou and Greenhalgh(1992)、Zhou et al.(1993)计算了不同的计算方案所对应的零炮检距脉冲响应.他这些计算结果表明, 伪谱法可以产生很好的效果, 子波延续时间很短, 并且没有二次扰动.四阶有限差分的准确性要相对差一些.子波开始时是正确的, 由于网格频散, 子波延续时间太长, 并有二次扰动产生.9点褶积微分算子计算结果接近于伪谱法, 比四阶有限差分法要好得多(戴志阳等, 2005).
4.2 与伪谱法的时间比较
假设二维空间计算区间的大小为N×N, 褶积微分算子的长度为M, 对声波方程进行数值模拟时, 若使用伪谱法, 需要N2log2N次复数乘法; 而若使用褶积微分算子法, 只需要M×N2次实数乘法. 显然, 褶积算子法与伪谱法的计算时间比为M/(4log2N), 只要M < (4log2N), 褶积微分算子法将比伪谱法的计算速度快(Zhou et al., 1993).例如, 在本文256×256的计算区域内, 9点褶积微分算子与伪谱法的计算时间比为9∶ 32, 即伪谱法的计算速度此时是本文方法的3.5倍.
5. 结论与讨论
本文所构造的方法是将褶积微分算子应用于地震波场正演模拟的又一个尝试.该微分算子具有快速衰减和反对称性质.通过引入高斯窗, 极大程度地压制了算子截断所引起的吉谱斯效应.理论模型计算结果表明, 该方法计算速度快, 效率高.这种方法是继戴志阳等(2005)发展的褶积微分算子法之后的又一种弹性波数值模拟的颇具潜力的工具.基于该理论所开发的时间错格有限差分-褶积微分算子法地震波场数值模拟程序, 能有效地对均匀和非均匀介质中的地震波传播进行数值模拟.将本文方法与伪谱法进行比较, 发现本文发展的方法对数值频散的抑制能力跟伪谱法相当, 但计算时间却明显占优势, 可以更为高效地模拟地震波在复杂介质中的传播过程.
-
图 1 9点褶积微分算子振幅随采样点的变化曲线
Fig. 1. Amplitude versus sampling point for nine point convolutional differentiator
图 3 Marmousi速度模型及500 ms时的波场快照
a.速度模型; b.x分量; c.z分量
Fig. 3. Snapshots for the Marmousi velocity model at the time of 500 ms
图 5 Forsyte褶积微分算子法波场快照
a.x分量; b.z分量
Fig. 5. Snapshots obtained by Forsyte convolutional differentiator method
-
[1] Alterman, Z., Karal, C., 1968. Propagation of seismic waves in layered media by finite difference methods. Bull. Seism. Soc. Am. , 58(1): 367-398. [2] Cerjan, C., Kosloff, D., Kosloff, R., et al., 1985. A nonreflecting boundary condition for discrete acoustic and elastic wave equations. Geophysics, 50(4): 705-708. doi: 10.1190/1.1441945 [3] Dai, Z.Y., Sun, J.G., Zha, X.J., 2005. Seismic wave modeling with hybrid order convolution differentiator. Computing Techniques for Geophysical and Geochemical Exploration, 27(2): 111-114(in Chinese with English abstract). [4] Graves, R.W., 1996. Simulating seismic wave propagation in 3D elastic media using staggered grid finite difference. Bull. Seism. Soc. Am. , 86(4): 1091-1106. [5] Komatitsch, D., Tromp, J., 1999. Introduction to the spectral-element method for 3-D seismic wave propagation. Geophys. J. Int. , 139: 806-822. doi: 10.1046/j.1365-246x.1999.00967.x [6] Kosloff, D., Baysal, E., 1982. Forward modeling by a Fourier method. Geophysics, 47(10): 1402-1412. doi: 10.1190/1.1441288 [7] Kosloff, D., Reshef, M., Loewenthal, D., 1984. Elastic wave calculation by the Fourier method. Bull. Seism. Soc. Am. , 74(3): 875-891. doi: 10.1785/BSSA0740030875 [8] Lin, S.H., Zhao, L.Y., 2007. High precision time domain forward modeling for crosshole electromagnetic tomography. Journal of China University of Geosciences, 18 (4): 320-325. doi: 10.1016/S1002-0705(08)60012-6 [9] Magnier, S., Mora, P., Tarantola, A., 1994. Finite difference on minimal grids. Geophysics, 59(9): 1435-1443. doi: 10.1190/1.1443700 [10] Teng, J.W., Zhang, Z.J., Yang, D.H., et al., 1995. Digital simulation of tricomponent seismic data in anisotropic media through convolution differential operator. Geophysical Prospecting for Petroleum, 34(3): 14-22(in Chinese with English abstract). [11] Xie, J., 1981. Methematical methods for the data processing of geophysical prospecting. Geological Publishing Press, Beijing, 97-104(in Chinese). [12] Yu, S.P., 1996. Wide-band Ricker wavelet. Oil Geophysical Prospecting, 31(5): 605-615(in Chinese with English abstract). [13] Zhang, Z.J., Teng, J.W., Yang, D.H., 1996. Simulations of the acoustic and elastic wavefield by convolutional differentiator. Acta Seismologica Sinica, 18(1): 63-69 (in Chinese). [14] Zhou, B., Greenhalgh, S.A., 1992. Seismic scalar wave equation modeling by a convolutional differentiator. Bull. Seism. Soc. Am. , 82(1): 289-303. [15] Zhou, B., Greenhalgh, S.A., Zhe, J., 1993. Numerical seismogram computations for inhomogeneous media using a short, variable length convolutional differentiator. Geophysical Prospecting, 41: 751-766. doi: 10.1111/j.1365-2478.1993.tb00882.x [16] Zhou, H., Xu, S.Z., Liu, B., et al., 1997. Modeling of wave equations in anisotropic media using finite element method and its stability condition. Chinese Journal of Geophysics, 40(6): 833-841(in Chinese with English abstract). [17] 戴志阳, 孙建国, 查显杰, 2005. 地震波混合阶褶积算法模拟. 物探化探计算技术, 27(2): 111-114. https://www.cnki.com.cn/Article/CJFDTOTAL-WTHT200502004.htm [18] 滕吉文, 张中杰, 杨顶辉, 等, 1995. 各向异性介质中褶积微分算子法三维地震资料的数字仿真. 石油物探, 34(3): 14 -22. https://www.cnki.com.cn/Article/CJFDTOTAL-SYWT199503001.htm [19] 谢靖, 1981. 物探数据处理的数学方法. 北京: 地质出版社, 97 -104. [20] 俞寿朋, 1996. 宽带Ricker子波. 石油地球物理勘探, 31(5): 605-615. https://www.cnki.com.cn/Article/CJFDTOTAL-SYDQ199605000.htm [21] 张中杰, 滕吉文, 杨顶辉, 1996. 声波与弹性波场数值模拟中的褶积微分算子法. 地震学报, 18(1): 63-69. https://www.cnki.com.cn/Article/CJFDTOTAL-DZXB601.008.htm [22] 周辉, 徐世浙, 刘斌, 等, 1997. 各向异性介质中波动方程有限元法模拟及其稳定性. 地球物理学报, 40(6): 833-841. doi: 10.3321/j.issn:0001-5733.1997.06.012 -
下载:
下载:
点击查看大图
计量
- 文章访问数: 3768
- HTML全文浏览量: 113
- PDF下载量: 41
- 被引次数: 0
