0.   引言

复杂地质构造具有倾角大、断层发育、分布范围小而无规律、埋藏深、速度在横向变化剧烈等特征, 地震勘探的常规处理难以获得明晰的地下复杂地质构造信息.近年来, 在油气工业的实践中, 叠前深度偏移是成像逆冲断裂、深部古潜山内幕等复杂地质构造的关键技术.叠前深度偏移成像的方法, 原理上可分为Kirchhoff积分叠前深度偏移和基于波动方程的叠前深度偏移.Kirchhoff积分叠前深度偏移有计算速度快、适应各种观测系统和易于建立偏移速度场.但存在难以处理多值走时及振幅问题, 以及焦散点和盲区问题.虽然很多研究者对此作了大量努力, 其结果仍不能令人满意.

基于波动方程理论波场延拓的叠前深度偏移方法具有理论完整、成像精度高、保幅和适于处理复杂构造等特点, 一直是研究复杂介质偏移成像的主要内容.波场延拓算子的实现: 一类用有限差分(Ristow and Ruhl, 1994) 来实现, 高角度成像及三维方位角各向异性(Ristow and Ruhl, 1997; Zhang et al., 2000) 的处理需要大量的计算; 另一类用傅里叶变换在频率-波数域实现的显格式方法.其方法尽管在处理横向均匀介质时是准确的(Gazdag, 1978), 但对横向变速介质单程波算子则需进行近似处理.近几年该类方法得到较大的发展, 主要有分裂步相移法(SSF)、广义屏法(GSP)、傅里叶有限差分(FFD)、优化分裂步(OSP)、可分表示法(OSA) 等. Stoffa et al. (1990)提出了分裂步相移法, 将介质慢度(速度) 分解为背景慢度和扰动慢度两部分, 基于小扰动理论, 将相移法推广到横向变速介质中, 该方法适于处理横向弱变速介质.广义屏方法(Le Rousseau and de Hoop, 2001) 在分裂步相移方法基础上对平方根算子作二次泰勒展开, 算子近似精度得到了进一步改善, 但在复平面内支点附近存在不稳定问题.优化分裂步法(Liu and Zhang, 2006) 对泰勒展开系数进行优化近似处理, 而使其数值性能得以较大改善.最优可分表示法(宋健, 2001; Chen and Liu, 2004; 张致付等, 2005) 把波场算子象征用两个单变量函数的乘积求和实现, 来实现算子在波数域和空间域的分离, 波数域内进行相移计算, 在空间域对介质横向变速引起的时移作修正.该方法不用背景速度作相移运算, 对强横向变速介质中各速度值的波场都能很好地归位.二维、三维地震数据的偏移结果表明该方法成像复杂介质中的小断层、透镜体具有较好的能力.

波动方程叠前深度偏移计算量巨大, 在偏移的速度分析中, 一般需要多次偏移计算, 其计算量难以忍受.如何减少叠前的数值计算以提高偏移算法的计算速度也是值得探索的问题.一种方法是面炮技术, 从单炮数据合成叠加接收点的记录来减少偏移的数据量, 该方法在共炮数据上用单平方根方程实现偏移.一种替代的方法是在中点-偏移距域数据上用双平方根方程实现偏移.Popovici (1996)和Jin et al. (2002)分别用分裂步法和广义屏法实现了双平方根方程偏移.对于三维数据而言, Biondi and Palacharla (1996)在共方位角域实现了双平方根波动方程偏移.他们的研究成果显示双平方根方程的共偏移距方法相对于炮域方法可降低波动方程偏移的计算量.

本文将从中点-偏移距坐标双平方根方程出发, 用最优可分表示法来近似双平方根方程波场延拓算子的象征, 采用正反傅立叶变换构造了其波场延拓算子, 并实现成像.用该方法计算非零偏移距的脉冲响应, 以考察可分表示法构造的算子用于强横向变速精度, 并对模型数据作偏移计算.

1.   共偏移距波场延拓

1.1   双平方根算子及象征

横向变速介质中的频率-波数域的双平方根方程(Claerbout, 1985) 如下:

这里, x_s、x_r分别表示炮点和检波点位置, ω是圆频率, 和v (x_r, z) 分别表示深度z上的炮点和检波点位置的速度, k_s、k_r是波数, P (k_s, k_r, z, ω) 是频率-波数域波场, 式(1) 是炮点-检波点坐标下的双平方根方程, 可把它转化为中点-偏移距坐标下的双平方根方程, 假定:

这里, m是中点, h为半偏移距, 把式(2) 代入式(1), 作变换可得中点-偏移距坐标下的双平方根方程为:

当介质均匀时, 该波场延拓算子是一个简单的相移运算, 在这种情况下, 算子由傅立叶变换、相因子的乘积、反傅立叶变换三部分组成, 采用快速傅立叶变换实现.横向非均匀时, 式(3) 的双平方根算子是拟微分算子, 介质中的延拓算子与相移算子有类似的结果, 差别是相因子是依赖空间坐标, 这样需要对不同空间速度值上都要进行傅立叶变换, 这使得大量的计算被要求.从薄板传播子理论(Le Rousseau, 2001), 可以给出该双平方根方程波场延拓算子的象征:

这里E=exp[iΔzk_z (k_m, k_h, z) ]称为双平方根方程波场延拓算子的象征, k_z (k_m, k_h, z) 是与空间坐标有关的垂直波数.当横向速度变化时, 对波场延拓算子的处理即对象征的近似计算, SSF, GSP, FFD等方法对象征的逼近通过泰勒展开, 把象征分解为包含背景速度的相移运算和在空间坐标下对速度变化的时移运算.这种处理各速度值的波场成像精度与背景速度有关.

1.2   象征的可分表示法逼近

根据上述近似方法的思想构造横向变速介质的单程波算子的新方法, 该方法是用最优可分表示逼近双平方根方程延拓算子的象征, 该算子在介质各速度值上都有较高的精度, 而不受背景速度的影响.同时, 在计算中仍然用快速傅立叶变换来实现.

这里, u_r=ω/v (x_r, z), u_s=ω/v (x_s, z), 上式是一个4个变量的方程, 很难求解, 令u= (u_r, u_s), k= (k_m, k_h), 我们可借助于单程波算子象征可分表示的公式, 把波场延拓算子的象征写为(Chen and Liu, 2004; Chen et al., 2007) :

对于β^(l) (k) =β^(l) (k_m, k_h) 有相同的可分近似表达式:

把式(7) 代入式(6), 得

这里λ=λ₁λ₂, 由于炮点和检波点速度对应空间位置上是相同的, 在波场延拓时的空间变量速度项是相同的, 把式(8) 写成:

由单程波算子象征的可分表示逼近的结果可知, 只需取较小的近似阶数s、t, 象征就可获得较高的近似精度.

图 1是3种不同方法(SSF、GSP、OSA) 对双平方根方程的算子象征数值逼近比较, 当偏移距为零, 这样偏移距项函数省略, 象征被简化为E (u, k_m), 取介质速度v在1 500~2 500 m/s变化, k_m在0~π/25变化, ω取40 π时, 用40、100分别对u、k_m (或表示为v、k) 均匀采样.从图 1可以看出3种方法在窄传播角上相位误差都很小, 分裂步相移法与广义屏法随着传播角的增加其相位误差逐渐增加, 广义屏法在90^°角上出现了奇点.5阶最优可分表示法近似在所有传播角上都能很好地逼近算子的象征, 只有在大角度振幅微弱损失.所以用最优可分表示近似算子的象征具有很好的逼近, 所构造的双平方根方程波场延拓算子具有较高的精度.

图  1  不同近似方法数值逼近象征的对比

Fig.  1.  Numerical computation comparison of symbol with different approximate methods

下载: 全尺寸图片 幻灯片

1.3   共偏移距数据的波场延拓

共偏移距双平方根方程偏移是通过波场延拓来实现, 把薄层介质中波场延拓算子的象征用最优可分表示的近似式表示, 同时把该近似式代入式(3), 得:

这里, FFT和IFFT分别表示正反傅立叶变换, 把所有频率的波场延拓后进行求和, 就是在t=0及h=0时的波场, 这就是双平方根方程的偏移波场.

1.4   算子的脉冲响应

为了检查最优可分表示偏移方法在强横向变速介质中的偏移效果, 用横向变速情况下的波场延拓算子来计算脉冲响应.模型采用x方向200个点, 采样间隔10 m, 纵向200个点, 采样间隔5 m, 主频15 Hz的Ricker子波作为点脉冲放置在0.4 s处, 时间采样间隔为4 ms.

对于横向不均匀介质, 假定介质横向速度从1 500~2 500 m/s随机变化, 对u以间隔40均匀采样, 用5阶最优可分近似不同速度单程波算子的象征, 并构造双平方根方程的波场延拓算子.图 2从上到下分别是2 500 m/s介质中理论解、SSF法、OSA法计算的脉冲响应比较, 可以看出, 在小传播角上, SSF法和OSA法都很好地逼近理论解, 对于横向变速介质, 这两种方法在小角度都能对波场较好地成像, 但是在大传播角, SSF法计算的脉冲响应偏离了理论值, 这因为SSF算子的象征数值近似在大角度相位误差较大; 而OSA法的脉冲响应在大传播角几乎接近于理论值, 只有能量略有减弱.说明该方法对大角度波场成像优于SSF法.图中也显示共中点-偏移距的脉冲响应是一个椭圆, 且随着偏移距的增加, 椭圆率在增大.

图  2  脉冲响应的比较

a, b分别是offset=200 m和400 m的理论解; c, d分别是offset=200 m和400 m的SSF法脉冲响应; e, f分别是offset=200 m和400 m的OSA法脉冲响应

Fig.  2.  Comparison of impulse responses

下载: 全尺寸图片 幻灯片

2.   模型与实测数据计算

为了进一步检验本文成像方法, 用SEG/Marmousi模型数据对OSA方法进行共偏移距双平方根方程叠前深度偏移试验.模型速度场是IFP基于Cuanza盆地的地质构造剖面给出的, 单炮记录用2D声波有限差分法模拟.如图 3a所示, 横向737个采样点, 纵向750个采样点, 速度场水平采样间隔12.5 m, 最大深度是3 000 m, 深度采样间隔为4 m.数据采样点数是750, 采样率是4 ms, 总采样长度3 000 ms, 总共240炮, 每炮96道.图 3b是第100、101炮的模拟地震记录, 由于模型构造复杂, 多断层, 地层也基本非水平, 单炮记录上绕射波明显, 反射波同相轴不是双曲线状, 都是倾斜的, 且同相轴不连续.

图  3  Marmousi模型与单炮记录

Fig.  3.  Marmousi model and shot records

下载: 全尺寸图片 幻灯片

图 4是采用了5阶共偏移距及炮域最优可分表示叠前深度偏移结果, 从剖面上看, 成像较明晰, 剖面显示3条断层面清晰, 断点收敛基本到位, 深部含气低速度体的成像准确, 对于倾角小的同相轴, 其成像能力较好.炮域法的成像效果要略好于共偏移法.因为对于倾斜地层, 抽取的共偏移距数据对应的CMP与对应的不完全一致.但该方法可选取任意偏移距的记录, 可根据需要来选取目标数据偏移, 这样能够减少计算量, 能用于叠前深度偏移速度分析.

图  4  Marmousi模型偏移结果

Fig.  4.  Marmousi data test with OSA migration method

下载: 全尺寸图片 幻灯片

3.   结论

最优可分表示方法是区域逼近的表达方法, 故对区域内各速度离散点上的双平方根方程波场延拓算子的象征都能很好地逼近, 用可分表示法构造的双平方根方程波场延拓算子具有较高的精度, 且波数域变量与空间(速度) 域变量分离.在波数域内进行相移计算, 在空间域对因介质横向变速引起的时移作修正.该方法在中点-偏移距坐标中同时向下延拓炮点和检波点波场, 地震数据偏移不需要逐炮计算, 具有较高的效率.对于三维数据, 地震数据要抽成共方位角数据道集, 其构造算子的方法类似, 这样可以大大减少计算量(Biondi and Palacharla, 1996).

该方法计算了变速情况下算子的脉冲响应, 其脉冲计算值与理论值基本一致, 特别在大传播角仍然能够很好地逼近理论波前面, 从而说明共偏移距最优可分表示叠前深度偏移能够适应强变速介质内的构造成像.实现了二维Marmousi模型偏移, 其结果显示该方法成像复杂构造的能力.由于该方法的偏移数据是中点-偏移距道集, 可针对目标选取数据成像, 这样可以用于波动方程的叠前偏移速度估计.