高效的成像域最小二乘逆时偏移

陈生昌， 李代光， 金成玫

浙江大学地球科学学院, 杭州 310027

0 引言

地震数据偏移成像是得到地下构造图像的主要方法技术，特别是地下复杂构造区的构造成像(Gray et al., 2001；Yilmaz, 2001；Etgen et al.，2009；Robein, 2010; Leveille et al.，2011).随着国内“两宽一高”地震数据采集技术的广泛应用和油气勘探开发目标的日趋复杂化和精细化，对于地震数据偏移成像方法技术要求能充分地利用地震数据的波形信息，得到的偏移成像结果不仅要满足复杂区域构造成像的要求，还要具有高分辨和保真性，使之能精细地反映地下反射面岩性变化的特征，以满足地震数据岩性处理解释，特别是复合型复杂油气藏和非常规油气藏地震数据岩性处理解释的需要.逆时偏移是一种基于全波方程的偏移方法(Baysal et al., 1983; MeMECHAN, 1983; Whitmore, 1983)，被认为是当前在理论上最为先进的偏移方法(Gray et al., 2001).相比于通常仅考虑单一射线路径走时的Kirchhoff偏移(Schneider, 1978; Bleistein et al., 2005)，逆时偏移可有效地考虑地震波传播过程中的多路径与多走时.相比于基于地震波方向分解的单程波偏移存在的大角度波(相对于深度方向的大角度波)限制(Ristow and Ruhl, 1994; Zhang et al., 2005)，逆时偏移无地震波传播角度限制.

给定具有地震波运动学准确的光滑偏移速度模型，常规逆时偏移可视为是一种利用一次波信息对有关地下反射率的线性反演问题的粗糙求解(Claerbout, 1985; Schuster, 2017).利用线性反演理论建立的有关地下反射率反演的最小二乘逆时偏移是对常规逆时偏移的改进，它可有效地消除地震子波、地震数据观测系统、成像目标上覆介质空间变化和地震波几何扩散对常规偏移成像的影响，改善地震波弱照明区和深部的偏移成像效果，得到具有高保真性和高分辨率的偏移成像结果(Tarantola, 1984；Yao and Jakubowicz, 2012; 黄建平等, 2013; Yang and Zhang, 2019).由于最小二乘逆时偏移涉及偏移结果的成像域到地震数据的时空域之间线性映射矩阵L的求逆，对于三维炮集地震数据的最小二乘逆时偏移，矩阵L的行数M和列数N，分别有M=Ny×Nx×Nt和列数N=Ny×Nx×Nz(其中，Ny、Nx和Nz分别为成像域的Y、X和Z方向的网格点数，Nt为地震道采样点数)，假设Ny=1000，Nx=1000，Nz=1000，Nt=10000，则L为1010×109的超大规模矩阵，致使利用矩阵L求逆的方法进行最小二乘逆时偏移在现有计算机条件下难以实现.因此，当前对于最小二乘逆时偏移的求解有两类方法，一是数据域方法(Schuster, 1993; Fletcher et al., 2012; Zhang et al., 2015; Zhao et al., 2015; Xue et al., 2016; Fletcher and Cavalca, 2018)；二是成像域方法(Yu et al., 2006; Aoki and Schuster, 2009; Tang, 2009; Fletcher et al., 2016).在数据域最小二乘逆时偏移方法中，以观测的炮集地震记录与数值计算的炮集地震记录间的最佳拟合为准则，采样Landweber迭代方法(陈生昌和周华敏，2018)避免最小二乘逆时偏移中大规模矩阵求逆的计算复杂度.但在数据域最小二乘逆时偏移的迭代中需要进行炮集地震数据的偏移和利用偏移成像结果的反偏移，所以数据域最小二乘逆时偏移方法的计算量巨大，约为常规逆时偏移方法计算量的2×Nite(Nite为迭代次数)倍.在成像域最小二乘逆时偏移方法中，利用常规逆时偏移结果可构建常规逆时偏移成像结果与最小二乘逆时偏移成像结果之间基于Hessian矩阵H(Hessian矩阵H为最小二乘逆时偏移的目标函数关于地下反射率的二阶导数，也是矩阵L的伴随与L的乘积，即H=L*L，H中的列向量可视为成像点的点扩展函数)的空变线性积分关系，使成像域成为了最小二乘逆时偏移的计算空间.对于上述Ny=1000，Nx=1000，Nz=1000的三维最小二乘逆时偏移，矩阵H为109×109的超大规模矩阵，因此在现有计算机条件下也同样难以利用矩阵H求逆的方法实现成像域的三维最小二乘逆时偏移.此外，由矩阵H的表达式可知，形成矩阵H的计算量对于三维最小二乘偏移也是十分巨大的，因此在成像域利用Landweber迭代方法求解也同样存在巨大的计算量问题.鉴于成像域最小二乘逆时偏移的计算量与偏移成像区域大小有关，目前该方法多用于区域较小的目标区最小二乘逆时偏移(Valenciano et al., 2006; Tang and Biondi, 2013; Zhao and Sen, 2018).

成像域最小二乘逆时偏移中成像点的点扩散函数的局部分布形态主要由地震数据观测系统和用于偏移的速度模型决定，而地震数据观测的多次覆盖特点和偏移速度模型的空间光滑性，所以成像域的点扩散函数具有一定的相似性.针对上述存在的问题和认识，本文把常规逆时偏移成像结果与最小二乘逆时偏移成像结果之间在成像域的空变积分关系近似为局部区域的空不变褶积关系，提出基于局域空不变反褶积的成像域最小二乘逆时偏移方法.通过利用求解空不变反褶积的波数域伪广义逆方法(陈生昌和王芳,1993；陈生昌和肖鹏飞，2007)，可得到一种稳定高效的成像域最小二乘逆时偏移方法.把本文提出的成像域高效最小二乘逆时偏移方法应用于模拟和实际地震数据均取得了理想的偏移成像结果.

1 最小二乘逆时偏移简介

给定具有地震波运动学准确的光滑偏移速度模型，作为一种线性反演的地震数据偏移所对应的地震数据线性正演方程可写为

d(xr,xs,ω)=L(xr,x,xs,ω)m(x),

(1)

式中，d(xr,xs,ω)为炮道集的频率域一次波地震数据；L(xr,x,xs,ω)为偏移成像结果到地震数据的线性映射算子；m(x)为偏移成像结果；xr为接收点坐标；xs为震源点坐标；ω为频率；x为地下成像点坐标.关于方程(1)中线性算子L(xr,x,xs,ω)的具体表达式，当前有三种不同的版本：一是平面波无穷大平反射面的经典反射版本(Claerbout, 1971; Bleistein et al., 2001)；二是基于Born近似的散射版本(Schuster, 2017；Yang and Zhang, 2019)；三是局部平面波局部平反射面的局部反射版本(陈生昌和周华敏，2018).对于经典反射版本的线性算子L(xr,x,xs,ω),其具体表达式有

L(xr,x,xs,ω)=Gr(xr,x,ω)Gi(x,xs,ω)s(ω)，

(2)

对于散射版本的线性算子L(xr,x,xs,ω),其具体表达式有

(3)

对于局部反射版本的线性算子L(xr,x,xs,ω),其具体表达式有

(4)

上述式中，vm(x)表示用于偏移的光滑速度模型；Gr(xr,x,ω)表示散射(或反射)传播的Green函数；Gi(x,xs,ω)表示入射传播的Green函数；s(ω)表示震源子波时间函数的Fourier变换.

把方程(1)简写为矩阵方程形式，有

d=Lm，

(5)

式中，d为地震数据d(xr,xs,ω)所对应的矩阵；L为线性算子L(xr,x,xs,ω)所对应的矩阵；m为偏移成像结果m(x)所对应的矩阵.利用线性最小二乘反演方法，由方程(5)可得到有关m的最小二乘解估计，即

mls=(L*L)-1L*d,

(6)

式中L*为L矩阵的伴随.对于方程(6)的解表达，由于矩阵L的超大规模，致使难以通过直接求矩阵逆的方式得到解mls.针对这样的大规模矩阵方程，在计算数学中常采用近似迭代的方式进行求解，如Landweber迭代(Kirsch, 2011).对于方程(6)，有如下的Landweber迭代格式，

(7)

式中，k为迭代次数；τ为解的修正步长.由迭代格式(7)可知，其反演运算主要是在数据域进行，因此应用迭代格式(7)的最小二乘偏移被称为数据域最小二乘偏移方法.

对于方程(6)可写为

mls=H-1mc,

(8)

其中，H=L*L,根据线性最小二乘反演理论可知，H矩阵为线性反演问题(5)的Hessian矩阵；对于mc=L*d，L*d表示对地震数据d进行常规逆时偏移运算，mc表示常规逆时偏移结果.

由公式(8)可知，把Hessian矩阵H的逆作用于常规逆时偏移结果就可得到最小二乘逆时偏移结果.由于公式(8)的运算仅在成像域进行，通常把下述三步运算组成的最小二乘逆时偏移称为成像域最小二乘逆时偏移方法，即

(9)

如果偏移成像的范围比较大，直接按公式(9)进行最小二乘逆时偏移存在形成Hessian矩阵H和求取逆矩阵H-1的两大计算困难，它们涉及的计算量巨大.因此，由公式(9)组成的成像域最小二乘逆时偏移比较适合偏移成像范围比较小的区域，如面向目标区的最小二乘逆时偏移(Tang and Biondi, 2013; Zhao and Sen, 2018).为了进行大区域三维地震数据的成像域最小二乘逆时偏移，我们需要研究可实现公式(9)所表示的成像域最小二乘逆时偏移的近似算法.

2 高效的成像域最小二乘逆时偏移方法

公式(8)所表示的反演求解形式有与其对应的正演表示式，即

mc=Hmls.

(10)

在公式(10)中，mc和mls分别表示由常规逆时偏移结果和最小二乘逆时偏移结果组成的列矩阵，即mc=[mc1,mc2,…,mcN]T，mls=[mls1,mls2,…,mlsN]T,它们分别是正演表达式中的数据和模型，N代表成像域的网格点数；H为N×N的方阵.

对于矩阵方程(9)，如果矩阵mls中仅第i个元素值等于1，其他的元素值都等于0，则矩阵mc的值等于矩阵H中第i个列向量的值，即

(11)

因此，H中的列向量值也称为成像域中成像点的点扩散函数(PSF)值.根据矩阵H的表达式，H=L*L，利用公式(2)或(3)或(4)，可得到下述多炮多道的PSF计算式(利用公式(3))，

(12)

式中，x为PSF在成像域中的计算点位置；x0代表成像域中的成像点位置；PSF(x,x0)为成像点x0的点扩散函数；Re表示取实部运算.

式(12)所表示的计算实质是在偏移速度模型vm(x)下对成像点x0处单位脉冲式散射体的地震波场多炮模拟和多炮逆时偏移成像，即从炮点xs激发的震源波场传播到成像点x0，与该点处的单位脉冲式散射体发生作用形成散射虚源激发散射波，散射波场传播到接收点xr，得到记录波场，然后多道记录波场经伴随传播到PSF计算点x，同时从炮点xs激发的震源波场也经伴随传播到计算点x，最后再利用成像条件得到偏移成像结果.因此，式(12)所得到的PSF(x,x0)即为成像点x0处单位脉冲式散射体的逆时偏移成像结果.由地震波的传播理论和偏移成像理论可知，PSF(x,x0)具有以成像点x0为中心的局部分布特征.该局部分布特征为多成像点的PSF同时计算提供了可能.

根据地震波传播与成像理论、矩阵H的表达式和PSF的计算式(12)，公式(10)对应一个线性空变积分，因此可把式(10)写为下述的空变积分形式，即

(13)

式中k(x,y)为与式(10)中矩阵H对应的空变积分核.与式(10)对应的反演转化为式(13)的空变线性积分方程求解.由于积分核k(x,y)的空变，对于积分方程(12)的求解没有快速算法.

k(x,y)的空变特征与偏移速度模型vm(x)的空间变化和地震数据观测系统范围的有限性、覆盖的不均匀性密切相关.如果vm(x)为均匀模型、观测系统具有无限大的观测范围且均匀覆盖，则k(x,y)退化为空不变，相应的式(13)退化为线性空不变褶积型积分方程，即

(14)

对于式(14)的线性空不变褶积型积分方程可利用波数域伪广义逆方法(陈生昌和王芳，1993；陈生昌和肖鹏飞，2007)进行稳定高效地求解，即

(15)

由于偏移速度模型vm(x)为空间宏观变化的光滑速度模型和地震数据采集的规则多次覆盖特性，我们认为光滑模型vm(x)和多次覆盖下k(x,y)的空间变化是很缓慢的，即k(x,y)在一定的空间范围内可近似为空不变的.对应的在一定空间范围内的各成像点的PSF(x,x0)的变化形态具有很大的相似性.

图1为用于成像点PSF计算的Marmousi模型的光滑偏移速度模型.图2为9个不同网格点((340,100)、(370,100)、(400,100)，(340,130)、(370,130)、(400,130)，(340,160)、(370,160)、(400,160))位置成像点的PSF分布图.计算PSF的观测系统为地面均匀全覆盖.由图2所示的PSF分布图可看出，成像点的PSF分布图不仅具有很好的以成像点为中心的局部分布特征，而且x和z方向距离30个网格点(网格点间距dx=dz=10 m)的成像点的PSF还具有很强的相似性，因此可认为在30×30网格点的区域内各个成像点的PSF形态变化是很小的，也即在该范围内的k(x,y)可近似为空不变的.图3为图1 Marmousi模型上间距30×30网格的所有成像点的PSF分布图(图3的显示加了perc=99的增益).由图3同样可以看出图中相邻成像点PSF分布图的相似性.

图1 用于成像点PSF计算的Marmousi模型的光滑速度模型Fig.1 The smooth velocity model of Marmousi model for the PSF calculation of imaging points

图2 9个不同网格点位置成像点的PSF分布图Fig.2 The PSF distribution of 9 imaging points at different grid-point positions

图3 图1速度模型上网格点间距30×30的所有成像点的PSF分布图Fig.3 The PSF distribution of all imaging points with 30×30 grid-points interval in velocity model of Fig.1

根据上述有关k(x,y)的局部分布特性和光滑速度模型、规则多次覆盖下k(x,y)在一定空间范围内可近似为空不变的认识，式(13)的空变积分方程可近似为局域空不变褶积的叠加，即

(16)

式中，Lb表示可把k(x,y)近似为空不变的局部区块；kLb(x)表示以局部区块的中心点xb为成像点的PSF近似的空不变核函数.把式(15)表示的线性空不变褶积型积分方程求解方法，应用于式(16)，有

(17)

(18)

(19)

式中，变量β、γ和φ的意义见示意图4.

图4 区块边界拼接示意图Fig.4 The sketch of block boundary tapering

把上述的常规逆时偏移、多炮多道的PSF计算和式(17)的分块局域空不变褶积型积分方程求解方法结合起来，就可得到一种稳定高效的成像域最小二乘逆时偏移方法，我们也把这种最小二乘逆时偏移方法称之为成像域的局域空不变反褶积方法.常规数据域最小二乘逆时偏移的计算量约为常规逆时偏移的2×Nite(一般情况下，Nite>10)倍，而本文的最小二乘逆时偏移的计算量约为常规逆时偏移的2倍，所以本文方法相对于常规数据域最小二乘逆时偏移方法约提高计算效率Nite倍.

在本文提出的成像域最小二乘逆时偏移方法中，局域区块Lb尺寸的大小不仅关系到方法的计算量，也关系到方法的最终偏移成像效果.局域区块Lb尺寸大，虽然方法的计算量小，但对k(x,y)的空变性考虑不足，会影响方法的最终偏移成像效果.如果局域区块Lb取为整个偏移成像区域，虽然方法的计算量最小，但相当于把k(x,y)近似为全成像域空不变(如方程(13)).局域区块Lb尺寸小，虽然方法的计算量大，但能较好地考虑k(x,y)的空变性，有利于最终的偏移成像结果.如果局域区块Lb取为一个点，则方法就等同于第1节公式(8)的成像域最小二乘逆时偏移成像方法.因此，区块Lb尺寸大小的选择原则为，在保证区块中心点xb的PSF(x,xb)完整的条件下，使区块Lb的尺寸尽量小，即在兼顾计算量的同时尽可能考虑k(x,y)的空变性以保证偏移成像效果.PSF(x,xb)的完整也就是褶积算子的完整，这是得到高分辨率反褶积结果的必要条件之一.

3 数值试验

为了验证本文提出的把空变的PSF近似为局域空不变PSF的可行性、正确性和局域空不变反褶积成像域最小二乘逆时偏移方法对实际数据的有效性，我们进行了基于二维Marmousi模型的合成地震数据、二维和三维实际地震数据的数值试验.在试验中，要求已得到地震数据的常规逆时偏移结果和用于逆时偏移的光滑偏移速度模型.根据PSF的局部分布特性，在给定的光滑偏移速度模型上，按照确定的间距规则地放置单位脉冲式点散射体，点散射体间的间距所确定区域即认为是可把k(x,y)近似为空不变的局部区块，点散射体(即成像点)位于局部区块中心；然后根据地震数据采集的观测系统利用公式(12)计算点散射体(成像点)在对应的局部区块的PSF分布；最后把常规逆时偏移结果和各个局部区块的PSF分布(对应的Fourier变换谱)应用于本文提出的稳定高效的成像域高效最小二乘逆时偏移计算式(17)，得到最小二乘逆时偏移成像结果.由于各个局部区块的PSF仅分布于其对应的区块内，因此利用公式(12)的一次计算就可以计算出成像域内各个局部区块的PSF，计算代价约为常规逆时偏移的2倍.

3.1 Marmousi模型合成地震数据

Marmousi模型为二维模型，试验中，用于成像点PSF计算的局部区块尺寸为30×30网格点，即成像点在x方向和z方向的间距均为30个网格点.图3为利用图1所示的光滑偏移速度模型和公式(12)计算得到的全成像域间距30×30的所有成像点的PSF分布图.由于k(x,y)的本质空变特征，所以各个局部区块的成像点PSF还是存在差异的，尤其是深部和边部的成像点PSF.图5为Marmousi模型合成数据利用图1所示的光滑偏移速度模型得到的常规逆时偏移成像结果.

图5 Marmousi模型数据的常规逆时偏移成像结果Fig.5 The conventional RTM result of Marmousi model data

图7的试验结果表明，即使k(x,y)在理论上是空变的，由于用于偏移和PSF计算的速度模型的光滑性和多次覆盖的数据采集，在实践中把k(x,y)近似为局域空不变，得到的成像域最小二乘逆时偏移效果是可行的和正确的.图6的试验结果表明，用某个局部区块的PSF分布作为全成像域的PSF分布，虽然也能得到不错的成像效果，计算量也很小，但由于不能考虑k(x,y)的空变特性，我们不建议使用全域空不变反褶积进行成像域最小二乘逆时偏移.

图6 全域空不变反褶积得到的 Marmousi模型数据的成像域最小二乘逆时偏移成像结果Fig.6 The IDLSRTM result of Marmousi model data by global spatial-invariant de-convolution

图7 局域空不变反褶积得到的 Marmousi模型数据的成像域最小二乘逆时偏移成像结果Fig.7 The IDLSRTM result of Marmousi model data by local spatial-invariant de-convolution

3.2 二维实际地震数据

为了进一步验证本文提出的局域空不变反褶积最小二乘逆时偏移方法在实际地震数据上的效果，我们首先用某地的一条二维线进行试验.图8为该线的常规逆时偏移成像结果.试验中计算各个局部区块中心点PSF分布的速度模型是用于该线逆时偏移成像的偏移速度模型，观测系统采用实际地震数据采集的观测系统.用于成像点PSF计算的局域区块尺寸为20×40网格点(网格点间距dx=20 m，dz=10 m)，即成像点在x方向和z方向的间距分别为20和40个网格点.图9为计算得到的成像域网格点间距20×40的所有成像点的PSF分布图.

图8 二维实际数据的常规逆时偏移成像结果Fig.8 The conventional RTM result of 2D field data

图9 成像域网格点间距20×40的所有成像点的PSF分布图Fig.9 The PSF distribution of all imaging points with 20×40 grid-points interval in imaging domain

图10 二维实际数据的成像域最小二乘逆时偏移成像结果Fig.10 The IDLSRTM result of 2D field data

3.3 三维实际地震数据

为了验证本文方法对三维实际地震数据最小二乘逆时偏移成像的有效性，我们对某工区的三维实际地震数据进行了试验.与上面的二维实际数据试验一样，我们首先有该实际数据的三维逆时偏移成像结果，并利用实际地震数据采集的观测系统和逆时偏移的偏移速度模型计算三维全成像域各个局部区块中心点的PSF分布.局域区块尺寸为20×20×50网格点(网格点间距dx=30 m，dy=30 m，dz=10 m)，即成像点在x方向和y方向的间距均为20个网格点，在z方向的间距为50个网格点.

图11为其中一个局部区块中心点的三维PSF分布的xz平面切片图.由于该PSF的成像点位于成像域的中浅部，PSF的完整性比较好.

图11 三维PSF分布的xz平面切片图Fig.11 The xz plane slice of 3D PSF distribution

图12、13分别为纵测线131的常规逆时偏移成像结果和本文提出的成像域最小二乘逆时偏移成像结果，图14、15分别为纵测线231的常规逆时偏移成像结果和本文提出的成像域最小二乘逆时偏移成像结果，图16、17分别为横测线100的常规逆时偏移成像结果和本文提出的成像域最小二乘逆时偏移成像结果.通过比较图12与图13、图14与图15和图16与图17可看出，成像域最小二乘逆时偏移成像结果相对于常规逆时偏移成像结果在空间分辨率方面有很大的提高.

图12 纵测线131的常规逆时偏移成像结果Fig.12 The conventional RTM result of In-Line 131

图13 纵测线131的成像域最小二乘逆时偏移成像结果Fig.13 The IDLSRTM result of In-Line 131

图14 纵测线231的常规逆时偏移成像结果Fig.14 The conventional RTM result of In-Line 231

图15 纵测线231的成像域最小二乘逆时偏移成像结果Fig.15 The IDLSRTM result of In-Line 231

图16 横测线100的常规逆时偏移成像结果Fig.16 The conventional RTM result of Cross-Line 100

图17 横测线100的成像域最小二乘逆时偏移成像结果Fig.17 The IDLSRTM result of Cross-Line 100

Marmousi模型合成数据的成像域最小二乘逆时偏移结果，验证了把空变PSF近似为局域空不变PSF的可行性与正确性；而在后续的两个实际数据的成像域最小二乘逆时偏移中, 对成像域进行了分块，使用各个局部区块中心点的PSF(即多点的PSF)，验证了本文提出的局域空不变反褶积的最小二乘逆时偏移的有效性.局域空不变反褶积对于不同位置的局域有不同的PSF，计算效率不如全域空不变反褶积，但它可以有效地考虑PSF在成像域的空变特性，有助于提高成像域最小二乘逆时偏移的成像效果,这是一种既兼顾PSF的(缓慢)空变特性，又具有计算高效性的方法.

4 结论

给定具有地震波运动学准确的光滑偏移速度模型，地震数据的一次地震波偏移成像可视为一个有关地下反射率分布的线性反演问题.利用最小二乘反演方法，在成像域常规逆时偏移结果是最小二乘反演的Hessian矩阵与最小二乘逆时偏移结果的空变积分.Hessian矩阵中的列向量也称为成像点的点扩散函数，是成像点单位脉冲体所产生的地震响应的逆时偏移结果，因此点扩散函数具有空间局部性.由于偏移速度模型的光滑性和地震数据采集多次覆盖，使得一定空间范围内成像点的点扩散函数具有很大的相似性.利用三维常规逆时偏移结果和以所有成像点的点扩散函数为列向量的Hessian矩阵求取最小二乘逆时偏移结果，不仅涉及形成Hessian矩阵的巨大计算量，而且还是一个超大规模的线性反演问题.利用一定空间范围内成像点的点扩散函数所具有的相似性，可把该空间范围内的常规逆时偏移结果近似为Hessian矩阵与最小二乘逆时偏移结果的空不变褶积，再利用求解空不变反褶积的波数域伪广义逆方法可稳定高效地得到最小二乘逆时偏移结果.利用点扩散函数的空间局部性，可同时计算成像域内以一定空间间隔均匀分布的多个成像点的点扩散函数，这也是成像域高效最小二乘逆时偏移方法的组成部分.在本文所提出的最小二乘逆时偏移方法中，局域区块Lb尺寸的大小是计算效率与偏移成像效果之间的折中.由于深部成像点的点扩散函数的分布范围大于浅部的点扩散函数的分布范围，因此局域区块Lb尺寸大小的选择应以深部点扩散函数的完整性为原则.Marmousi模型数据和实际数据的数值试验结果验证了本文所提出的成像域最小二乘逆时偏移方法的有效性.

致谢本文的研究工作得到了中石化科技部的资助，也得到了中石化石油物探技术研究院段心标、白英哲、陶永慧、孙敏傲等专家的大力帮助，在此表示感谢.