波场分解算法与逆时偏移角道集

2018-08-01王美霞张晓慧

石油物探 2018年4期

王美霞,张晓慧,张钋,唐冰,徐昇

(1.Statoil Gulf Services,Houston 77042;2.前斯塔托尔北京技术服务有限公司,北京100022)

复杂构造区的地震勘探,如大角度地层或盐丘等需要高质量的速度模型和高精度的成像方法。逆时偏移[1-4]利用双程波传播方程,能够准确模拟反射波、折射波、多次波等各种地震波,是目前最为精确的成像技术之一。

逆时偏移算法在20世纪80年代就已经提出,但由于其巨大的计算量和超量的计算资源需求,直到21世纪初才有实际应用。逆时偏移应用面临着理论和计算两方面的挑战。理论方面主要有两点:①逆时偏移低频噪声问题[4]。逆时偏移成像是通过对震源波场和检测器波场应用传统的零延迟互相关成像条件[5]得到的,如果存在强反射,这个过程不仅在实际反射界面处产生图像,也会沿着波路径产生严重的低频噪声。这些噪声可能掩盖真实的构造,给解释工作造成困难。将这个传统的成像条件和真振幅偏移方法[6-11]做比较,我们发现主要区别在于依赖于震源和检波器波场夹角的权重函数。实际上,真振幅偏移中关于反射角的函数可以压制很多反向散射噪声,关键在于如何计算逆时偏移中的反射角。②快速生成偏移共成像点道集问题。共成像点道集是质量控制、速度建模、属性分析等的重要输入数据,被广泛应用于复杂地区数据处理。常用的共成像点道集是共偏移距道集,也就是使用地面偏移距作为索引的道集。然而,在复杂情况下,共偏移距道集具有很多由于多路径传播而引入的偏移假像,而且很难直接由逆时偏移得到高质量的共偏移距道集,需要进行特殊处理[12]。相比之下,由逆时偏移产生角道集[13-15]就比较自然。此时角道集由地下反射角作为索引,偏移噪声也相对较弱。如果使用真振幅逆时偏移,角道集还可以为AVA分析提供输入数据。角道集可以使用拓展成像条件的方法得到[16-19],也可以直接由波场产生[11-20]。这些方法在二维情况下计算效率都比较高,但是在三维情况下计算量都很大,所以在实现逆时偏移算法时首先要考虑算法效率以及算法的并行计算和优化等。计算效率方面,引入Poynting矢量方法提高了生成角道集的计算效率,但是由于检波器波场比较复杂,一般Poynting矢量仅被用于震源波场,使用逆时偏移成像或反射界面倾角而不是检波器端波传播方向信息[21-22]。在三维道集逆时偏移应用中,面对存储量大、计算量高等挑战,Poynting矢量方法更具吸引力。YOON等[23]利用Poynting矢量方法估计波传播方向,该方法的缺点是在一个时空点处,Poynting矢量只能给出一个方向,因此,当波场复杂时,不同传播方向的波混叠在一起,Poynting矢量就不能准确计算方向。若用于角道集计算,就会使得有些能量被成像在错误的角度上,进而降低角道集质量或导致角道集存在部分能量缺失,这严重影响Poynting矢量方法在检波器端波场的应用。

本文进一步分析了Poynting矢量方法。我们认为只需要在使用Poynting矢量计算波传播方向之前对波场进行分解,便可以将其同时用于震源和检波器波场。波场分解越精细,Poynting矢量就越精确,但计算量也会随之增加。对于一般的地面观测系统,初至波主要是震源的下行波,反射波主要是检波器的下行波。如果模型中存在强速度差界面,波场中初至波和反射波就会叠加,因此进行简单的上、下行波场分解[2-4]就可以显著提高Poynting矢量的精度。

传统的波场分解算法[2-4]一般是通过傅里叶变换[25]在频率-波数域进行的。分解计算效率取决于快速傅里叶变换的实现效率。然而,傅里叶变换是全局运算,计算耗时。另一种方法是使用希尔伯特变换(HT)[24,26-28]进行波场分解。下行波和上行波可以通过原始波场和在时间-深度方向进行希尔伯特变换后的波场组合得到[2-7]。这和利用傅里叶变换得到的结果一致。优点是希尔伯特变换的实现比较灵活;可以通过快速傅里叶变换实现;也可以在时间-空间域使用有限长度卷积算子实现,以便进一步进行算法优化以提高计算效率。本文使用SSE和多线程优化,提高利用希尔伯特变换进行波场分解的计算效率。

利用希尔伯特变换进行波场分解时,需要考虑时间和空间方向的希尔伯特变换。有两种方法可以实现波场在时间方向的希尔伯特变换。一种是利用希尔伯特变换后的震源函数进行波传播得到一个新的波场[27],该方法会加倍波传播计算量和存储量,因此读写时间也会增加。另一种办法是只使用原始波场,在计算道集的过程中边计算边对原始波场进行希尔伯特变换。若效率和算法不是瓶颈,后者需要的数据读写较少。因此,本文采用后者。

本文首先简单回顾真振幅逆时偏移成像和角道集理论以及Poynting矢量方法。其次,我们证明通过傅里叶变换或希尔伯特变换都能够进行波场分解,且理论上等价,进一步使用单指令多数据流(SIMD)优化希尔伯特变换算法,并使用数值例子验证了优化后算法的精度和效率。然后,将优化的希尔伯特算法用于波场分解,进一步使用多线程优化并行算法。利用数值计算例子说明这种方法在得到和传统傅里叶方法一致结果同时,能够显著提高数值计算效率。最后,将这种加速的波场分解算法应用于逆时偏移产生角道集,利用二维Sigsbee2A和三维SEAM TTI计算实例说明本文方法的有效性和稳定性。

1 方法理论

利用双程波方程进行波传播的数值计算[29-35],可以得到震源波场us(x,t)和检波器波场ur(x,t)。这里t表示时间,向量x表示空间位置。在二维介质中,x=(z,x)。在三维介质中,x=(z,x,y)。下面我们首先回顾成像条件和角道集理论,之后引出一种快速的波场分解算法。

1.1 成像条件

偏移成像的传统互相关成像条件[5]可以表示为:

(1)

式中:xs和ys表示震源位置;ω为角频率。此成像条件在反射点处产生像,但同时也在满足成像条件的路径上产生假像。假像和真像的区别在于在假像处震源波场和检波器波场的传播方向相异。YOON等[23]提出只有当震源和检波器波传播方向所成角度在一定范围内时才使用互相关条件成像,且波传播方向由Poynting矢量计算得到,该方法的一个缺点是很难选取用于成像的最大角度。COSTA等[36]提出了依赖于角度的光滑函数作为偏移权重,在一定程度上克服了最大成像角度选取的困难。与真振幅偏移[6-7,10-11]相比,其和传统成像条件的区别在于依赖于震源和检波器波场夹角的权重函数。因此,一个更加适当的成像条件可以表达为[6,11,36-38]:

(2)

式中:θ表示反射角。由于偏移噪声是当震源和检波

器波传播方向成较大角度时产生的,所以权重函数cos2θ可以有效压制偏移中的低频噪声。

1.2 使用Poynting矢量产生真振幅逆时偏移角道集

共炮逆时偏移的真振幅三维角度域共成像点道集[11,39,40]可以按下式计算:

(3)

式中:R(x,θ0,φ0)表示在空间位置x处方位角φ0和反射角θ0方向上的反射系数;δ(θ-θ0)和δ(φ-φ0)指狄拉克函数,它们将积分限制在要求的角度θ0和φ0,v(x)为速度。

反射角和方位角计算如下:

其中,p=1/2(ps+pr),x,z为坐标轴的单位向量;ps和pr分别表示震源和检波器的波传播方向。震源波传播方向可以按下式利用Poynting矢量进行计算:

(6)

类似地,检波器波传播方向pr也可以使用Poynting矢量计算得到。

Poynting矢量计算是在时间-空间域进行的,当波场比较复杂时,Poynting矢量并不精确。图1a展示了从Sigsbee2A模型得到的一个波场快照,可以看到波场很复杂,几种不同类型的波交织在一起,这使得Poynting矢量很难给出一个正确的方向。然而,如果我们仅关注下行波(图1b),波场就变得简单很多,Poynting矢量可以较准确地计算波传播方向。

图1 由Sigsbee2A得到的波场快照(a)及其下行波场分量(b)

因此,波场分解,特别是上、下行波场分解(对于地表观测系统),能显著提高Poynting矢量计算方向的精度,这也可以进一步提高角道集的质量。

1.3 波场分解

1.3.1 利用傅里叶变换进行波场分解

传统方法利用傅里叶变换在频率-波数域进行波场分解。约定时间-深度方向的二维傅里叶变换如下:

(7)

于是下行波在频率-波数域可表示为:

(8)

上行波在频率-波数域可表示为:

UU(kz,ω)=U(kz,ω)-UD(kz,ω)

(9)

在公式(8)和公式(9)中,对于满足ω·kz=0的U(kz,ω),我们使用了系数1/2以便使得UD(kz,ω)+UU(kz,ω)=U(kz,ω)。这是合理的,因为满足ω·kz=0的波既非上行波也非下行波。

1.3.2 利用希尔伯特变换进行波场分解

SHEN等[27]指出上、下行波可以使用希尔伯特变换进行波场分离。他们使用一对波场实现分解:震源函数传播得到的原始波场和进行希尔伯特变换后的震源函数传播得到的波场。这种方法很有效,但计算效率不高,它需要两倍的波传播计算量。本文提出另一种实现方法,仅使用原始波场,在计算道集过程中进行希尔伯特变换。重点以上、下行波分解来描述算法,并将该方法拓展到其它方向的波场分解。

通过傅里叶变换和希尔伯特变换之间的关系,最终下行波可以表示为:

(10)

上行波可以表示为:

(11)

式中:Ht[u(x,t)]表示波场u(x,t)在时间方向的希尔伯特变换;Hz[u(x,t)]表示深度方向的希尔伯特变换。公式(10)和公式(11)与SHEN等[27]提出的公式一致,具体推导过程见附录A。

公式(10)和公式(11)提供了另一种波场分解的办法。希尔伯特变换可以简单通过傅里叶变换实现,但是这样计算效率与传统频率-波数域的波场分解一样。为了提高计算效率,我们提出直接在时间-空间域通过卷积进行希尔伯特变换,并进一步使用SSE优化该算法。

对于一维输入数组aj,j=0,1,…,n-1,使用(2m+1)长度卷积算子进行的离散希尔伯特变换如下:

(12)

我们注意到进行希尔伯特变换的卷积算子是反对称的,即:hk=-h2m-k,k=0,1,…,m。而且大约一半的系数hk都是0。基于这两条性质,公式(12) 可以进一步改写为如下形式以减少计算量:

(13)

式中:k+=2表示k的步长为2。

以下的实现和数值计算实例都是基于公式(13)。

1.3.3 希尔伯特变换波场分解算法的优化

若将波场分解用于逆时偏移中进行上、下行波分解,我们需要时间和深度方向的希尔伯特变换(如公式(10)和公式(11))。约定波场在内存中的存储规律为最快维为深度,最慢维为时间。因此,提高波场分解算法速度的关键是分别优化快维和慢维方向的希尔伯特变换。下面先说明如何优化二维数组的快维和慢维的希尔伯特变换,然后,将其用于波场分解。

给定二维数组ai,j,i=0,1,…,n2-1,且j=0,1,…,n1-1(下角标i表示慢维索引,j表示快维索引),记数组{ai,j}的希尔伯特变换为数组{bi,j}。传统方法在快维方向的希尔伯特变换可直接表示为:

(14)

若要计算慢维方向的希尔伯特变换,传统方法是先通过转置将慢维变为快维,再由公式(14)计算,最后再做一次数组转置。对数据做两次转置增加了计算量。我们提出一种不用转置的算法,即按照公式(15)直接对慢维进行数值运算:

(15)

由于我们使用相对较短的卷积算子,k相对于慢维长度n2较小,就可以使用SIMD沿快维j方向加速算法,即将每一列ai-m+k,*-ai+m-k,*看作一个元素,同一个系数h2m-k与该元素相乘,这样就可以使用SSE同时处理多个数据。我们将这种优化的算法记为慢维方向的向量化希尔伯特变换(vectorized Hilbert transform,VHT)。

通过使用SSE优化,慢维方向的希尔伯特变换不再需要对数组进行转置。理论上,数值计算效率可以提高到4倍(SSE)或8倍(AVX,Advanced Vetor Extensions)。但是考虑到从缓存中获取数据也需要时间,实际提高水平可能略低于理论值。

对于快维方向的希尔伯特变换,我们发现如下使用SSE优化的算法比直接利用公式(14)计算更加高效:

(16)

在以上过程中,二维转置可以使用优化的SSE转置函数(_MM_TRANSPOSE4_PS)实现。尽管使用了转置,SSE向量化仍然可以大幅提高计算效率。我们称此为快维方向的VHT。

由于重点是上、下行波分解,所以进行波场分解算法的核心是二维波场的分解。对于三维情况,y方向是最慢维,只需对每一个y切面波场进行二维波场分解。在使用SSE向量化的同时,还可以使用OpenMP进行多线程并行计算以提高效率。对于地震成像应用,目前的缓存一般足够容纳二维子波场。因此,基于公式(10)和公式(11),我们提出以下流程进行高效波场分解:

1) 对每一个ix,将二维子波场u(z,ix,t)从内存加载到缓存中,并使用慢维方向的向量化希尔伯特变换计算Ht[u(z,ix,t)];

2) 对每一个it,将二维子波场Ht[u(z,x,it)]从内存加载到缓存中,并使用快维方向的VHT计算HzHt[u(z,x,it)];

3) 使用公式(10)(或公式(11))得到下行(或上行)波场。

对于三维波场的分解,我们只需要对每一个iy切片进行以上二维波场分解。

1.3.4 数值算例

先以一个二维波场快照u(z,x)为例说明VHT的计算效率及其在波场分解中的应用。输入波场是由常速度v=1790m/s模型通过伪谱法[41]模拟产生,震源函数是主频为14Hz的Ricker子波,我们对子波进行了3Hz的高通滤波。波传播网格为dx=dz=15m,nz=nx=701,且深度z为快维。震源在区域中心。我们使用一个相对较长的卷积算子,m=40,即81个点,进行希尔伯特变换以便得到较好的精度。

首先,我们验证与传统方法相比,慢维方向的VHT能够提高计算效率。图2是t=2s时刻的波场快照。图3a是使用向量化希尔伯特变换计算得到的二维波场在慢维x方向的希尔伯特变换波场。图3b比较了慢维方向的VHT和传统希尔伯特变换的计算效率。为了得到可靠的计算时间统计,我们重复了5000次计算。从图3b可以看到,对于慢维方向的希尔伯特变换,向量化算法比传统算法效率高三倍多。为检验VHT的精度,我们取图3a 中一条深度方向的记录与傅里叶方法得到的结果进行比较,见图3c。由于傅里叶方法对于所有波数都是精确的,可以认为其得到的结果为理想结果。图3c中VHT结果与傅里叶方法的结果十分吻合,说明本文的VHT方法准确有效。在进行VHT时,卷积算子越长,精度越高。在实际应用中,如果不需要特别高的精度,可以使用一个较短的卷积算子,这样可以进一步提高计算效率。

图2 t=2s时刻的二维波场快照

其次,我们再验证快维方向的VHT。图4a是图3a 所示波场快照在深度方向的希尔伯特变换。图4b 比较了VHT方法和传统HT方法用于计算快维方向希尔伯特变换的计算效率。我们同样进行了5000次重复计算。可以看到向量化方法可以减少大概一半的计算时间。

图3 使用VHT得到的波场快照(a)、VHT和传统希尔伯特变换(HT)计算效率对比(b)以及深度z=5250m处VHT结果和傅里叶变换结果(FT)比较(c)

图4 快维(z方向)的希尔伯特变换a 使用VHT得到的波场快照在z方向的希尔伯特变换; b z 方向的VHT和传统希尔伯特变换(HT)的效率比较(计算重复了5000次)

从以上两个实例(图3b和图4b),我们可以看到:传统HT方法在计算慢维方向希尔伯特变换时比计算快维方向希尔伯特变换时需要时间更多;相反地,向量化方法在计算快维方向希尔伯特变换时比计算慢维方向希尔伯特变换时需要更多时间。这均是由于是否需要数据转置所致。

最后,我们进一步检验VHT用于波场分解的综合计算效率,并将其和传统的利用傅里叶变换的波场分解进行比较。为了合理地进行比较,在进行傅里叶变换时使用了FFTW(the Faster Fourier Transform in the West)[42]。波场时间方向采样为:nt=660,dt=4ms。图5a和图5b分别是由本文算法和傅里叶方法得到的下行波场在t=2s时的波场快照。图5c 比较了本文方法和傅里叶变换方法的计算效率。两种方法都是用16核计算机进行计算。为得到可靠的计算时间统计,运算重复了100次。可以看到本文提出的使用VHT的分解方法能够很大地提高计算效率,比传统使用傅里叶变换进行分解的方法快了6倍多。为进行更详细的精度比较,图6比较了从图5a和图5b中取出的x=4950m处一条深度方向的记录,可以看到两种方法得到的结果非常一致。

图5 上、下行波分解示例(采用t=2s 时刻的下行波分量)a VHT结果; b FFTW计算结果; c 两者效率比较(计算重复100次)

图6 下行波场中x= 4950m 处深度方向记录的比较

1.3.5 波场分解拓展

以上主要阐述了如何进行上、下行波分解,我们也可以将该方法拓展到其他方向的波场分解。为了理论的完整性,附录B列出了其他方向波分量的公式,包括左行波(ul)、右行波(ur)、左上行波(uul)、右上行波(uur)、左下行波(udl)、右下行波(udr)、左上后行波(uulb)、左上前行波(uulf)、右上后行波(uurb)、右上前行波(uurf)、左下后行波(udlb)、左下前行波(udlf)、右下后行波(udrb)、右下前行波(udrf)。

同理,用于上、下行波分解的优化方法也可以用于其它方向波场的分解。这里给出一个左、右行波场分解的例子来说明其有效性。我们同样使用图2所使用的模型。图7a和图7b分别是t=2s时刻的左行波和右行波。图7c比较了两种方法的计算效率。可以看到,本文提出的VHT方法能够显著提高左、右行波分解的计算效率。值得注意的是本文方法用于左、右行波分解时不需要任何数据转置,这一点也优于传统方法。

图7 向量化希尔伯特变换得到的t=2s 时刻的左、右行波分解结果a 左行波分量; b 右行波分量; c 两者效率比较(计算重复100次)

2 应用实例

将波场分解算法用于逆时偏移成像来检验其在逆时偏移和角道集生成中的作用。偏移像可以通过叠加角道集得到,这使得成像更加灵活,例如可以只叠加某一角度范围内的道集或使用合适的关于角度的权重函数。下面介绍了二维Sigsbee2A模型和三维SEAM TTI 模型例子。

2.1 二维Sigsbee2A例子

首先使用Sigsbee2A模型和数据。偏移中共使用500炮的合成地震数据。最大频率为36Hz。为了准确模拟波场,我们使用伪谱法[41]进行波传播,并对时间频散进行了校正[34]。图8a是不使用波场分解时得到的角道集,可以观察到盐丘下边界之上道集并不连续也不清晰。图8b为使用VHT进行波场分解得到的结果,道集连续性变好,噪声也得到压制。这说明波场分解能够提高Poynting矢量方向计算的精度。图9是使用互相关成像条件(方程(1))和未进行分解的波场得到的叠加图像,从图上可以看到严重的低频噪声,以致掩盖了真实的反射界面。图10a是使用包含角度权重的成像条件(方程(2))和未进行分解的波场得到的叠加图像,图10b是使用包含角度权重的成像条件(方程(2))和下行波得到的叠加图像。比较图9和图10a,可以看到依赖于角度的权重函数可以压制很多噪声,盐丘以下的反射层也变得清楚。但是图10a中还存在一些噪声,使盐丘下边界之上的图像比较模糊,如图中红色箭头所示区域。若进一步使用波场分解,叠加图像就变得比较清晰,质量得到提高,如图10b 所示。这也说明波场分解之后使用包含角度权重的成像条件(方程(2))可以有效降低偏移噪声。

图8 Sigsbee2A模型的角道集(道集抽取位置6004.56～25816.56m，间隔1524m)a 不使用波场分解得到的角道集; b 使用VHT进行波场分解得到的角道集

图9 利用互相关成像条件得到的二维Sigsbee2A模型叠加图像(没有使用波场分解)

图10 利用0～60°角道集叠加得到二维Sigsbee2A模型叠加图像(成像条件为依赖于角度的权重函数)a 未使用波场分解; b 采用VHT对下行波场进行波场分解

2.2 三维SEAM TTI例子

实际上,我们更加重视三维应用的计算效率。为此,我们使用三维SEAM (SEG Advanced Modeling) TTI 模型来说明本文算法在三维逆时偏移应用中的计算效率。SEAM 模型区域为40km×35km×15km。包含一个复杂的盐丘[43],这在墨西哥湾地区非常典型。不仅模型大小符合实际情况,合成地震数据也比较符合真实情形。由于计算资源有限,我们选取了215炮数据进行测试。图11展示了炮的位置。在算法中,我们使用了XU 等[33]提出的拟P波方程,这样不仅效率高而且没有S波产生的噪声。测试中,偏移最高频率为18Hz。以下我们选取位置y=32250m

处观测偏移结果。图12a和图12b分别是不使用和使用波场分解得到的结果。可以看到:图12a中,特别是红色椭圆标记区域,丢失了很多成像细节;图12b相对清晰很多,也有更多细节信息,特别是在盐丘之上。这个比较说明即便对于非常复杂的三维模型,波场分解也能够提高Poynting矢量的精度,从而提高角道集的质量。图13是不使用波场分解得到的叠加偏移图像(使用成像条件(1)),可以看到严重的噪声,真实的构造被掩盖。图14是使用成像条件(2)得到的图像,图14a和图14b分别使用波场分解前后得到的结果。比较图14a 和图13,可见关于角度的权重函数能够压制一部分噪声,但是图像质量还不是很高。进一步使用波场分解,可以得到更加清晰、细节更丰富的图像(图14b)。因此依赖于角度的权重函数可以压制一些噪声,波场分解可以进一步提高成像质量,去除偏移噪声。

图11 三维 SEAM TTI模型实例(选取215 炮的位置)

图12 三维 SEAM TTI模型角道集(角道集位置11～25km，间隔为1km)a 未进行波场分解得到的角道集; b 利用VHT波场分解后得到的角道集

图13 三维 SEAM TTI模型叠加偏移图像(采用互相关成像条件(方程(1)),未进行波场分解)

图14 利用0～60°的角道集叠加得到三维 SEAM TTI模型叠加偏移图像(采用依赖于角度的权重函数 (公式(2)成像)a 未使用波场分解; b 采用VHT进行波场分解

3 结论

角度域共成像点道集对复杂构造的地震成像十分重要,可以使用Poynting矢量高效计算波传播方向,但是当波场复杂时,传统的Poynting方法精度较低。为提高Poynting矢量的精度,可以在计算Poynting矢量之前将波场进行分解,再进行Poynting矢量计算,得到高精度的传播方向,从而得到高质量的角道集和偏移图像。我们提出了VHT和多核并行的波场分解算法。经过此优化,波场分解算法加快大约5倍。二维Sigsbee2A和三维SEAM TTI的数值计算实例均说明本文的波场分解算法能够有效去除逆时偏移图像和角道集中的噪声,提高精度和计算效率。该方法并不局限于上、下行波场分解,可以很容易拓展到其它方向的波场分解。本文算法也可以进一步应用于最小二乘偏移和全波形反演以消除成像噪声和提高成像质量。

致谢:感谢Hongbo Zhou 和Jun Mu 的帮助和讨论。感谢MIT的FFTW,SMAART JV的Sigsbee2A模型和数据,以及SEG的SEAM模型和数据。

附录A 公式(10)和公式(11)的推导

我们首先定义如下两个辅助函数:

通过使用辅助函数(A1)式和 (A2)式,可以将频率-波数域下行波的定义公式(8)改写成如下形式:

(A3)

同理,使用辅助函数(A1)式和 (A2)式,可以将频率-波数域上行波的定义公式(9)改写成如下形式:

(A4)

对于任意信号,其希尔伯特变换和傅里叶变换之间存在如下关系:

(A5)

式中:ξ表示频率或波数,F(ξ)表示信号的傅里叶变换,Fh(ξ) 表示信号的希尔伯特变换。根据我们关于希尔伯特变换和傅里叶变换的定义,σ(ξ) 是如下函数:

(A6)

根据函数定义(A1)式、(A2)式和 (A6)式,我们可得到以下g1(ξ),g2(ξ)关于σ(ξ)的表示: