2. 北京桔灯地球物理勘探股份有限公司, 北京 102200;
3. 吉林油田地球物理勘探研究院, 松原 138006
2. Beijing Orangelamp Geophysical Exploration Company Limited, Beijing 102200, China;
3. Geophysical Exploration and Research Institute of Jilin Oilfield, Songyuan 138006, China
基于弹性波传播方程高频近似的射线理论类偏移成像技术(Hill, 1990, 2001; 高成等, 2015; Gao et al., 2017;韩建光等,2017),其在实现偏移成像的过程当中,采用射线理论方法求解运动学和动力学射线追踪系统获得波场传播所需的走时振幅和相位信息,然后进行地震波场的延拓外推从而达到最终成像的目的.其实现过程和实现方式相对于波动方程类偏移成像具有很大的时效性和经济性,因此一直是工业界应用的主流成像方法和技术(岳玉波, 2011;杨珊珊等,2015;袁茂林等,2016).基于体波、面波等地震弹性波和电磁波的有效联合,解译了华南深部结构(Di et al., 2021).那么在应用射线理论实现波场延拓成像的过程当中,难免会遇到依赖速度模型进行的射线追踪系统的求解.在实际应用过程当中,一般给定的速度模型都比较复杂,而且是以离散的形式只给定网格节点上的数值,非网格节点上的值只能通过构建的插值算法来进行插值计算.在应用射线理论进行运动学和动力学追踪计算中,波场随着时间的不断外推增加,每一步计算所应用到的速度场信息不一定都落在给定的原离散网格上,如果计算所需要的速度场不在给定的离散节点之上,就需要对速度模型进行插值计算,从而解决在波前外推当中速度在x、z方向上非网格节点处的插值问题.众所周知,插值算法以及应用插值算法所得的插值计算结果的准确性直接影响着速度模型上计算的走时、相位和振幅,最终影响到偏移成像的质量(韩复兴, 2009).在偏移成像的过程当中,给定的速度模型往往都具有一定的速度梯度,如何根据速度梯度的空间变化选择一种高效、稳定的插值算法对于实现基于高频近似理论的射线类偏移成像有着十分重要的作用(Gao et al., 2017).
截止到目前为止,在数值计算当中常用的插值算法主要分为两大类型:第一大类型主要为通过给定的网格节点上的数值直接进行插值计算,例如双线性插值、分片线性插值、临近域插值、距离加权插值等.直接计算这种类型的插值算法经过几十年的发展和完善,算法相对成熟和稳定,而且简单实用、计算速度快.在实际插值计算处理当中,该类插值算法只需要4个已知网格节点上的属性值进行插值计算,但缺点在于只能达到二阶近似的计算精度,确保了计算的效率而忽略了计算的精度问题;第二大类型主要为基于构建插值核函数的插值计算方法,例如卷积类插值(韩复兴等, 2008a, b )、样条插值、双立方厄米特插值(王德人和杨忠华, 1990; Keys, 1981)、双三次多项式插值等等,该类插值算法依据所构建的核函数不同,所选取插值点周围的网格节点数也不尽相同,二维三次卷积、双三次多项式插值等需要插值点周围的16个网格节点上的属性值进行插值,二维四次卷积依据其核函数需要插值点周围的36个网格节点上的属性值进行插值.由于采用的网格节点数比较多,该类核函数插值算法计算精确度相对比较高,并且根据实际需求可以满足插值结果在网格点处的一阶导数、二阶导数连续,在实现的过程当中考虑算法的优化能够同时保证计算具有较高的计算效率.但该类算法存在的问题在于由于插值过程当中应用的插值点周围的网格节点数较多,虽然保证了插值结果的计算精度和导数的连续性,但对于速度梯度在空间上变化剧烈的模型来说,忽略了原始速度模型横向、纵向空间结构的变化,从而造成插值结果后的速度模型不满足原始速度模型空间梯度的连续性.也就是说基于插值核函数的插值算法的缺点在于其插值函数一旦固定下来对于整个模型空间来说就是不变的,这种不变性对于速度梯度变化不明显的模型来说是可以适应的,但对于速度梯度在横向和纵向变化较大的模型来说,这个固定核函数的插值算法就很难适应.如果应用基于固定核函数的插值算法在速度梯度变化剧烈的模型上进行插值计算就会造成梯度变化大的地方边缘模糊和锯齿现象,最终造成射线路径、走时和振幅相位的误差,影响偏移成像的结果.
为了解决上述射线类偏移成像求解射线追踪系统当中遇到的速度在波场延拓过程当中的不连续问题,本研究的主要内容就是考虑原始速度模型的空间结构,基于变分原理引入基于速度梯度构建的偏微分方程插值算法,使模型空间被插值点处速度梯度的空间变化尽可能与原始模型保持一致(Perona and Malik, 1990; Catté et al., 1992; Alvarez et al., 1992; Jiang and Moloney, 2002a, b ),从而实现基于原始速度模型空间梯度变化的插值计算.由于偏微分方程方法本身所固有的特性(局部特征不变性、解的唯一性和线性叠加性),自偏微分方程算法提出以来,随着时间的推移,目前已经应用到各个领域(Perona and Malik, 1990; Catté et al., 1992; Alvarez et al., 1992; 仵冀颖, 2009; 肖义男, 2005).因此,在射线类偏移成像当中应用基于速度梯度构建的偏微分方程插值算法可以最终实现不破坏原始速度模型空间结构前提下的插值计算,从而获得较为准确的走时、振幅和相位信息,最终提高偏移成像的质量.
1 偏微分方程(Partial Differential Equation)的插值算法基于速度梯度构建的偏微分方程插值算法在进行插值计算的过程当中引入总变分最小原则,使被插值点的速度值v(x, z)满足(Jiang and Moloney, 2002a, b ):
(1) |
也就是要求需要插值地方的速度的梯度尽可能与原始速度模型的空间结构保持一致,从而达到对原始速度模型空间结构的保护性.其约束条件为采用基于速度梯度构建的偏微分方程在网格节点上所得的插值计算结果必须等于原始速度模型网格节点上的速度值,即:
(2) |
其中v(i, j)是原始速度模型空间网格节点上的速度值,X和Z是原始速度模两个方向上的空间尺度,Δh、Δk是X方向和Z方向上的空间采样间隔,φ(·)≥0是一个递减函数.
对于公式(1)来说,初始条件的选择直接决定着该问题的求解难易程度.为了使公式(1)的求解简单化,引入对被插值点速度方向角的条件限制(Jiang and Moloney, 2002a, b ):
(3) |
式中α(x, z)=arg[∇v(x, z)],是速度v(x, z)的梯度角,ϕ() 是一个单调递减函数.该函数必须满足两个条件:(1) ϕ(B)=1;(2)
假设已知原始速度模型网格节点处的梯度角,且在应用基于速度梯度构建的偏微分方程插值计算过程当中原已知的梯度角不随插值条件发生改变,即限定条件为:
(4) |
根据公式(3)可知ϕ (‖∇α(x, z)‖) 是一个单调递减函数,由于cos(α(x, z)) 的取值范围就在0和1之间,所以选取ϕ (‖∇α(x, z)‖)=1-cos(‖∇α(x, z)‖),因此可以采用如下方向模型估计实际公式(3)中的
(5) |
在实际插值计算当中梯度角α(x, z) 可以采用数值计算的方法获得.
基于以上假设,在偏微分方程插值计算的过程当中,首先需要求解梯度角α(x, z),即(3)式,接下来求解基于约束条件(2)和约束条件(4)的变分问题,即可解决基于变分原理的插值计算问题.
关于插值点处梯度角的计算问题如下,对应(5)式的Euler方程为:
(6) |
由于在实际计算过程当中只使用插值点4个斜线方向上的相邻速度值,所以(6)式可以简化为:
(7) |
其中k∈{(i, j), (i, j+1), (i+1, j), (i+1, j+1)},αk表示对应于第k个速度点上的梯度角,因此得到非线性方程(7)的解为:
(8) |
其中
(9) |
将公式(9)展开并代入方向角限定条件化简得到:
(10) |
公式(10)是关于插值点处的梯度角以及插值点处速度偏导数的计算公式,将该公式采用数值计算的方法进行离散处理,然后应用被插值点周围的4个已知点上的速度值就可以进行插值计算求得被插值点处的速度值.其具体的计算步骤如图 1所示.
图 1中G(x, z)为需要插值的点,h为速度模型X方向上的网格间距,k为速度模型Z方向上的网格间距,v(i, j)、v(i+1, j)、v(i, j+1)、v(i+1, j+1)为已知网格节点上的速度值,A1、A2、B1、B2为计算过程当中的待求点,偏微分方程插值计算步骤如下:
(1) 首先应用四个已知网格节点上的速度值v(i, j)、v(i+1, j)、v(i, j+1)、v(i+1, j+1)求解待定点A1、A2、B1、B2,即:
(11) |
(2) 然后应用不等距差分和泰勒级数展公式求解被插值点处的速度导数值vxx(x, z)、vxz(x, z)、vzz(x, z),即:
(12) |
(3) 把通过第二步解出的值vxx(x, z)、vxz(x, z)、vzz(x, z)、vzx(x, z) 带入公式(10),就可以获得待求点处的速度值G(x, z),其公式为:
(13) |
其中:
(14) |
下面以地震数据处理当中常用的Marmousi速度模型和Sigsbee速度模型为例,分析说明在原始速度模型上进行插值计算后速度模型空间结构特别是模型边缘位置的变化情况.Marmousi速度模型横向384个网格节点,纵向122个网格节点,原模型横向和纵向的网格间距Δx=Δz=24 m,现在采用dΔx=dΔz=4 m的间距进行插值计算,插值后速度模型横向网格点数为2298,纵向网格点数为726,在相同的运算环境下,采用卷积类插值算法插值计算耗时18.899 s,采用偏微分方程插值计算耗时16.204 s,原Marmousi模型和经过卷积插值、PDE插值计算后的结果如图 2所示.为了进一步分析说明不同插值算法对速度模型空间结构和速度梯度变化区域边缘的影响,将原模型和插值计算后模型红色局部区域进行放大,如图 3所示.通过图 3可以看出,采用固定核函数的卷积类插值算法在原始速度模型空间结构梯度变化剧烈的区域,由于采用的插值点数较多(16个),插值造成的结果就是破坏了原模型的空间梯度结构,因此造成插值后速度模型红色部分同原始速度相比,在边缘区域变得模糊,而基于速度梯度构建的偏微分方程插值算法,由于去插值函数的单调性和可变性,插值结果基本保持和原始速度模型空间结构的一致性.
Sigsbee速度模型横向2400个网格节点,纵向800个网格节点,横向网格间距Δx=37.5 m,纵向网格间距为Δz=25 m.采用x方向间距为dΔx=12.5 m,z方向间距为dΔz=5 m进行插值计算,插值后速度模型横向网格点数为7197,纵向网格点数为3995,在相同的运算环境下,采用卷积类插值算法完成计算耗时408.972 s,采用基于速度梯度的偏微分方程插值完成计算耗时374.696 s.原模型以及插值计算后所得到的速度模型如图 4所示.通过对比插值前后的速度模型可以看出,对于速度梯度变化比较剧烈的区域,同样由于卷积类插值算法采用的插值点数较多,忽略了原模型梯度的空间变化,插值计算后整个模型的边缘区域变得模糊,而采用基于速度梯度的偏微分方程插值后的速度模型考虑和原模型速度梯度的空间变化,边缘区域都得到的较好的保持.
另外,通过插值前后的模型对比也可以看出,在速度模型底部红色圆点区域,基于速度梯度的偏微分方程插值能够较好的保持原始速度的局部特征,红色圆点清晰可见,而卷积类插值后红色圆点区域变得模糊.其主要是因为Sigsbee盐丘模型内外速度差异比较大,速度梯度在局部变化特别明显,卷积类插值核函数的空间不变性造成了其难以适应速度梯度变化剧烈区域的插值计算,从而造成对原模型空间结构的改变.
3 射线路径计算为了进一步说明偏微分方程插值的优越性,下面给出一些实际的速度模型,分别用卷积类插值和PDE插值计算射线路径并分析其差异.
速度模型1:起伏海底速度模型,模型横向网格点数2000,纵向网格点数1600,横向和纵向网格间距都为2.5m,层速度为:海水层
速度模型2:倾斜层速度模型,模型横向网格点数1600,纵向网格点数1500,横向和纵向网格间距都为4.0 m,层速度为:海水层
速度模型3:Marmousi速度模型:模型横向网格点数为384,纵向网格节点数为122,横向纵向网格间距都为24 m,射线路径计算时间步长为0.004 s,应用二维三次卷积插值以及PDE插值计算的射线路径如图 7所示.
从图 5-图 7的射线路径计算结果可以看出:对于起伏海底以及倾斜层速度模型,由于每一层层速度都为常数,插值计算只发生在层与层之间,由于不同插值算法应用的插值点数以及对边界的处理不同,应用二维三次卷积插值以及PDE插值得到的射线路径会在在第一个分界面后存在差异;而对于速度横向纵向剧烈变化的Marmousi速度模型来说,由于其速度模型横向纵向梯度变化十分明显,因此插值点速度由于采用的插值算法不同而造成图 7b中射线路径的差异.
4 偏移成像结果对比下面以南海复杂的崎岖海底大陡坡速度模型为例,验证在射线类偏移成像中插值算法的选取对成像结果的影响.复杂的大陡坡模型横向网格点数为3001,纵向网格点数为1501,横向和纵向网格间距都为2.5 m,如图 8a所示.数值模拟采用炮间距500 m,共16炮,每炮601道,道间距12.5 m,时间采样间隔2 ms.分别采用卷积类插值算法和偏微分方程插值算计算射线路径和走时并应用高斯波束实现最终的偏移成像,偏移结果剖面如图 8b、c所示.
从图 8的偏移成像剖面结果可以看出,对于复杂的南海某区域速度梯度变化较大的大陡坡速度模型,在射线类偏移成像当中采用基于速度梯度构建的偏微分方程插值算法对原始速度模型进行插值处理后获得成像结果其海水底部层位清晰可见,所获得成像效果要明显优于采用二维三次卷积插值所获得的偏移成像效果.
5 结论本文针对基于弹性波传播方程高频近似的射线理论类偏移成像过程当中求解运动学和动力学追踪系统当中所遇到的模型非网格节点处速度以及导数的插值问题,基于原始速度模型空间梯度变化,引入总变分最小原则,使插值点处的速度沿x方向和z方向的梯度与原始模型梯度保持一致的要求构建偏微分方程插值算子实现非网格节点处速度以及速度导数的插值计算.通过在两个速度模型上的插值对比分析可以得出,由于偏微分方具有很好的局部特征保持特性,应用基于速度梯度构建的偏微分方程插值算法对速度模型进行插值可以较好的保持原速度模型的空间结构和边缘特征,同时通过插值的耗时可以看出,偏微分方程插值速度要快于卷积类插值;通过不同模型的射线路径计算可以得出,由于不同的插值算法采用的插值点不同以及对速度模型的边缘处理不同,偏微分方程插值能更好的保持原始速度模型的空间结构,从而造成相同模型射线路径存在较大的差异.通过图 8南海复杂的大陡坡模型的成像结果也能够很好的说明应用基于速度梯度构建的偏微分方程插值算法对原始速度模型进行插值计算能够在保持原始速度模型空间结构的同时可以获得较好的成像结果.
Alvarez L, Lions P L, Morel J M. 1992. Image selective smoothing and edge detection by nonlinear diffusion. Ⅱ. SIAM Journal on Numerical Analysis, 29(3): 845-866. DOI:10.1137/0729052 |
Catté F, Lions P L, Morel J M, et al. 1992. Image selective smoothing and edge detection by nonlinear diffusion. SIAM Journal on Numerical Analysis, 29(1): 182-193. DOI:10.1137/0729012 |
Di Q Y, Tian F, Suo Y, et al. 2021. Linkage of deep lithospheric structures to intraplate earthquakes: A perspective from multi-source and multi-scale geophysical data in the South China Block. Earth-Science Reviews, 214(1): 103504. DOI:10.1016/j.earscirev.2021.103504 |
Gao C, Sun J G, Qi P, et al. 2015. 2-D Gaussian-beam migration of common-shot records in time domain. Chinese Journal of Geophysics (in Chinese), 58(4): 1333-1340. DOI:10.6038/cjg20150420 |
Gao Z H, Sun J G, Sun X, et al. 2017. A fast algorithm for depth migration by the Gaussian beam summation method. Journal of Geophysics and Engineering, 14(1): 173-183. DOI:10.1088/1742-2140/aa52dd |
Han F X. 2009. On some computational problems in wavefront construction method[Ph. D. thesis] (in Chinese). Changchun: Jilin University.
|
Han F X, Sun J G, Yang H. 2008a. Ray-Tracing of wavefront construction by bicubic convolution interpolation. Journal of Jilin University (Earth Science Edition) (in Chinese), 38(2): 336-340, 346. |
Han F X, Sun J G, Yang H. 2008b. Interpolation algorithms in wavefront construction. Chinese Journal of Computational Physics (in Chinese), 25(2): 197-202. |
Han J G, Wang Y, Chen C X, et al. 2017. Multiwave Gaussian beam prestack depth migration for 2D anisotropic media. Progress in Geophysics (in Chinese), 32(3): 1130-1139. DOI:10.6038/pg20170324 |
Hill N R. 1990. Gaussian beam migration. Geophysics, 55(11): 1416-1428. DOI:10.1190/1.1442788 |
Hill N R. 2001. Prestack Gaussian-beam depth migration. Geophysics, 66(4): 1240-1250. DOI:10.1190/1.1487071 |
Jiang H, Moloney C. 2002a. A new direction adaptive scheme for image interpolation.//Proceedings. International Conference on Image Processing. Rochester, NY, USA, USA: IEEE.
|
Jiang H, Moloney C. 2002b. Error concealment using a diffusion based method.//Proceedings. International Conference on Image Processing. Rochester, NY, USA, USA: IEEE.
|
Keys R. 1981. Cubic convolution interpolation for digital image processing. IEEE Transactions on Acoustics, Speech, and Signal Processing, 29(6): 1153-1160. DOI:10.1109/TASSP.1981.1163711 |
Munk W H. 1974. Sound channel in an exponentially stratified ocean, with application to SOFAR. The Journal of the Acoustical Society of America, 55(2): 220-226. DOI:10.1121/1.1914492 |
Perona P, Malik J. 1990. Scale-space and edge detection using anisotropic diffusion. IEEE Transactions on Pattern Analysis and Machine Intelligence, 12(7): 629-639. DOI:10.1109/34.56205 |
Wang D R, Yang Z H. 1990. Numerical Approximation Introduction (in Chinese). Beijing: Higher Education Press: 314-338.
|
Wu J Y. 2009. Research on partial differential equation-based methods for image processing[Ph. D. thesis] (in Chinese). Beijing: Beijing Jiaotong University.
|
Xiao Y N. 2005 Adaptive digital image interpolation and its application[Master's thesis] (in Chinese). Chongqing: Chongqing University.
|
Yang S S, Huang J P, Li Z C, et al. 2015. The review of imaging methods based on gaussian beam. Progress in Geophysics (in Chinese), 30(3): 1235-1242. DOI:10.6038/pg20150332 |
Yuan M L, Huang J P, Deng C, et al. 2016. Gaussian-beam based diagonal Hessian for amplitude-preserved imaging. Progress in Geophysics (in Chinese), 31(3): 1268-1273. DOI:10.6038/pg20160346 |
Yue Y B. 2011. Study on Gaussian beam migration methods in complex medium[Ph. D. thesis] (in Chinese). Qingdao: China University of Petroleum (East China).
|
高成, 孙建国, 齐鹏, 等. 2015. 2D共炮时间域高斯波束偏移. 地球物理学报, 58(4): 1333-1340. DOI:10.6038/cjg20150420 |
韩复兴. 2009. 论波前构建法中的几个计算问题[博士论文]. 长春: 吉林大学.
|
韩复兴, 孙建国, 杨昊. 2008a. 基于二维三次卷积插值算法的波前构建射线追踪. 吉林大学学报(地球科学版), 38(2): 336-340, 346. |
韩复兴, 孙建国, 杨昊. 2008b. 不同插值算法在波前构建射线追踪中的应用与对比. 计算物理, 25(2): 197-202. |
韩建光, 王赟, 陈传绪, 等. 2017. 二维各向异性多波高斯束叠前深度偏移. 地球物理学进展, 32(3): 1130-1139. DOI:10.6038/pg20170324 |
王德人, 杨忠华. 1990. 数值逼近引论. 北京: 高等教育出版社: 314-338.
|
仵冀颖. 2009. 偏微分方程在图像处理中应用的研究[博士论文]. 北京: 北京交通大学.
|
肖义男. 2005. 数字图像自适应插值技术及其应用研究[硕士论文]. 重庆: 重庆大学.
|
杨珊珊, 黄建平, 李振春, 等. 2015. 高斯束成像方法研究进展综述. 地球物理学进展, 30(3): 1235-1242. DOI:10.6038/pg20150332 |
袁茂林, 黄建平, 邓聪, 等. 2016. 基于高斯束的对角Hessian在保幅成像中的应用. 地球物理学进展, 31(3): 1268-1273. DOI:10.6038/pg20160346 |
岳玉波. 2011. 复杂介质高斯束偏移成像方法研究[博士论文]. 青岛: 中国石油大学(华东).
|