角度域共成像点道集(ADCIGs)既包含背景速度的影响, 也包含方位角度反射系数的信息, 因此可用于偏移速度分析和更新速度模型^[1-3]、照明分析^[4]、提取AVA信息^[5], 还可用于介质各向异性分析和提高构造成像质量等。

ADCIGs可以由射线类叠前偏移方法生成(比如克希霍夫偏移方法^[6-8]或者Beam类偏移方法^[9-10]), 也可由波动方程类偏移方法(比如单程波方程偏移^[11-12]或者逆时偏移(RTM)^[13-15])生成。与克希霍夫偏移、Beam偏移和单程波方程偏移相比, 逆时偏移能够精确地处理多波至现象和突破高陡倾角成像限制, 成为了复杂区域的首选成像方法。所以, 逆时偏移的角度道集生成方法研究非常重要^[16-17]。

目前生成逆时偏移角度道集的方法主要有两类:一类是成像后方法, 包括扩展成像条件方法^[18-21]和逆散射成像条件方法^[22-24]; 另一类是成像前方法, 包括波矢量方向估计方法^[25-34]和波场局部平面波分解方法^{[17, 35-39]}。

扩展成像条件方法首先生成地下局部偏移距道集, 然后将其转换为角度域成像道集, 该转换过程需要进行高维傅里叶变换和局部倾斜叠加计算, 存在内存和I/O问题。在进行局部倾斜叠加时, 时窗大小对生成道集的质量影响很大, 选择的局部时窗小时, 角度道集噪声少, 但分辨率低; 选择的时窗大时分辨率高, 但是倾斜叠加产生的假象多, 从而降低了角度域成像道集的质量, 影响AVA特征提取。

逆散射成像条件方法是利用方位角度滤波器对空间求导后的波场进行相关成像, 得到单方位角度的成像结果, 叠加后产生角度道集。尽管该方法可以自动地压制低频噪声, 并且增加的计算量不大, 但是其成像结果的振幅(包括子波相位)保真性及波场方向计算的准确性等方面还存在问题。

波矢量方向估计方法是指在波场传播过程中计算波矢量方向, 利用估计的波矢量方向判断波场的传播方向, 并利用角度域成像条件生成角度道集。由于坡印廷矢量方法计算效率高, 且具有较高的角度分辨率, 因此被广泛应用于逆时偏移角度道集的生成。在角度道集估计过程中, DICKENS等^[27]利用炮点端波场和检波点端波场的传播方向提取角度道集, 并分析了角度道集的振幅特征。但是当检波点端波场比较复杂时, 计算的坡印廷矢量通常会出现不稳定情况, 因此YOON等^[40]和VYAS等^[28]利用坡印廷矢量计算的炮点端波场传播方向和地下反射界面倾角来估计反射角; ZHAO等^[41]利用偏移剖面预测反射层法向向量, 并将其与稳定的炮点端方向波场结合起来计算角度道集。王保利等^[32]采用一阶声波波动方程求解坡印廷矢量, 减少了坡印廷矢量计算时间, 并利用高斯采样函数解决了角度不等间隔问题。虽然坡印廷矢量方法的计算效率比较高, 但是该方法无法解决波前重叠的问题^[42]。VYAS等^[29]和TANG^[43]先将波场进行方向分解, 然后利用分解后的波场来计算坡印廷矢量。在波矢量方向估计方法中, 最重要的是波传播方向的计算, ZHANG等^[30-31]利用P波的位移极化矢量表示波的传播方向; ZHANG^[44]采用光学流计算波场传播方向。

局部平面波分解方法是先将炮点端波场和检波点端波场分解为各自的局部平面波分量, 然后应用角度域成像条件来提取地下局部角度域成像道集, 通常在时间-空间域或频率-空间域中实现。XIE等^[35]利用局部倾斜叠加方法提取偏移波场中的角度域信息。YAN等^[37]提出在时间-空间域利用局部慢度方向提取炮点端和检波点端的波场方向来构建角度道集的方法, 这种方法可以解决波前重叠的问题, 但是存在角度分辨率低的问题, 而且需要储存足够的时间片来计算波传播方向。XU等^[17]给出了一种在频率-波数域直接构建角度道集的方法, 并采用抗泄露傅里叶变换^[45-46]解决了角度分辨率问题, 但是该方法需要在逆时偏移中将整个波场全部储存下来进行关于时间维的傅里叶变换, 计算代价很大。HU等^[47]提出基于解析时间波场外推与局部波场分解的逆时偏移角度道集生成方法, 解析时间波场只包含正频率的波场, 不需要进行时间维的傅里叶变换即可确定波场的传播方向, 解决了储存波场的问题, 但是该方法需要对震源和波场记录实施Hilbert变换, 分别进行正反传计算, 导致一次逆时偏移过程至少需要两次正传和两次反传, 计算代价还是较大。

本文结合波矢量方向估计方法和局部平面波分解方法的优点, 提出一种混合构建角度道集的方法。首先在常规逆时偏移基础上利用局部平面波分解类方法提取局部平面波, 然后利用波矢量方向估计方法计算的波场方向来确定最终的波前传播方向, 最后采用角度域成像条件生成角度道集。同时, 利用结构张量特征值快速判断波前是否交叉, 只在波前交叉的情况下使用本文提出的混合方法, 而在单一的波前面采用高效的波矢量估计方法, 进一步提高了计算效率。

1 波矢量方向估计方法

波矢量方向估计方法是指利用波场的振幅梯度和相位梯度等信息来计算地下任意时刻、任意点的波传播方向, 其最主要的难点是波前方向矢量的计算。计算波前方向矢量的方法通常有:坡印廷矢量方法、极化矢量方法、瞬时波数矢量方法以及光学流方法等。

坡印廷矢量原本是指电磁场中的能流密度矢量, 用于计算单位时间内穿过垂直于矢量方向的单位面积的电磁场能量。2006年, YOON等^[26]首次将坡印廷矢量引入到逆时偏移中, 用于确定波的传播方向。坡印廷矢量可由下式计算:

$ {\mathit{\boldsymbol{p}}}\left( {{\mathit{\boldsymbol{x}}}, t} \right) =-\nabla {\mathit{\boldsymbol{u}}}\left( {{\mathit{\boldsymbol{x}}}, t} \right)\mathit{\boldsymbol{\dot u}}\left( {{\mathit{\boldsymbol{x}}}, t} \right) $

(1)

其中, p(x, t)为地下x=(x, y, z)位置t时刻的坡印廷矢量, ▽u表示波场的空间导数, $\mathit{\boldsymbol{\dot u}}$表示波场的时间导数。经过正演和逆时外推后, 可分别得到炮点端波场u_s(x, t, p_s)和检波点端波场u_r(x, t, p_r), 利用相关成像条件获得角度道集:

$ \mathit{\boldsymbol{I}}\left( {\mathit{\boldsymbol{x}}, \theta } \right) = \smallint {\mathit{\boldsymbol{u}}_{\rm{s}}}(\mathit{\boldsymbol{x}}, t, {\mathit{\boldsymbol{p}}_{\rm{s}}}){\mathit{\boldsymbol{u}}_{\rm{r}}}(\mathit{\boldsymbol{x}}, t, {\mathit{\boldsymbol{p}}_{\rm{r}}}){\rm{d}}t $

(2)

其中, p_s和p_r分别代表炮点端波场和检波点端波场在地下x点处t时刻的波场传播方向; θ为地下反射角度, θ=arccos(p_s－p_r)。反射角的计算方法有多种, 如DICKENS等^[27]利用炮点端波场和检波点端波场传播方向计算反射角; YOON等^[40]和VYAS等^[28]利用炮点端波场传播方向和反射截面倾角计算反射角。

坡印廷矢量方法只需要计算波场的空间导数和时间导数, 计算效率高, 易实现, 具有较高的角度分辨率, 但易受噪声影响。特别是通过地表数据逆时延拓得到的检波点端波场, 由于只有地表的上边界数据作为波动方程的边值条件, 在进行逆时外推时, 波前无法干涉闭合, 加上存在噪声, 因而坡印廷矢量方法计算的方向矢量连续性差, 稳定性和准确性难以保证。另外, 坡印廷矢量方法假设地下任一点任一时刻只有一个传播方向, 当波场出现交叉或存在界面时, 波前方向的计算不准确, 不能得到正确的角度道集。

2 局部平面波分解方法

局部波场分解方法主要是利用局部傅里叶变换或局部倾斜叠加等方法将局部波场变换到其它域(空间-时间域、频率-波数域或时间-波数域)中实现角度分解。在各种局部波场分解方法中, 比较高效的是HU等^[47]提出的基于解析时间波场外推与局部波场分解的角度道集生成方法。由于解析时间波场中只包含正频率的波场信息, 该方法只需要在时间域逆时外推中进行空间维的傅里叶变换即可实现波场的角度分解。

假设波场方向分解借助傅里叶变换来实现, 将波场实施时间维和空间维的傅里叶变换后, 根据波场的频率和波数即可定义波场的上下行波^[48], 从而实现波场的方向分解。在时间域中实现这种波场分离时, 需要存储地下整个波场, 进行时间维的傅里叶变换, 其I/O量和计算量都很大。如果能够在外推过程中固定波场的频率, 那么波场的方向分解只需要利用波场的空间波数信息即可, 而空间波数却很容易获得。在信号处理中有一类时间信号叫做解析时间信号, 其实部为信号本身, 而虚部为信号的Hilbert变换, 它只包含正频率的信号。HU等^[47]将其引入到波场计算中, 提出解析时间波场的概念。解析时间波场是复波场, 其实部可以通过波动方程的传播过程得到, 其虚部可通过对地下波场进行Hilbert变换得到。但是由于Hilbert变换是关于时间方向的卷积, 因此需要存储整个波场并进行卷积, 不利于计算。Hilbert变换为一个90°相移算子, 借助于波动方程的震源项与波场具有线性关系^[49], 因此可以将波场的Hilbert变换转化为震源项的Hilbert变换, 利用新的‘震源’进行传播, 从而可以获得解析时间波场的虚部。

以炮点端为例, 解析时间波场的外推过程可表示为:

$ \left\{ \begin{array}{l} \frac{1}{{{v^2}\left( \mathit{\boldsymbol{x}} \right)}}\frac{{{\partial ^2}{\mathit{\boldsymbol{u}}_{\rm{s}}}(\mathit{\boldsymbol{x}},t;{\mathit{\boldsymbol{s}}_\mathit{\boldsymbol{x}}})}}{{\partial {t^2}}} - {\nabla ^2}{\mathit{\boldsymbol{u}}_{\rm{s}}}\left( {\mathit{\boldsymbol{x}},t;{\mathit{\boldsymbol{s}}_\mathit{\boldsymbol{x}}}} \right) = {\mathit{\boldsymbol{f}}_{\rm{s}}}\left( {t\mathit{\boldsymbol{,}}{\mathit{\boldsymbol{s}}_\mathit{\boldsymbol{x}}}} \right)\\ \frac{1}{{{v^2}\left( \mathit{\boldsymbol{x}} \right)}}\frac{{{\partial ^2}\mathit{\boldsymbol{u}}_{\rm{s}}^{\rm{H}}(\mathit{\boldsymbol{x}},t;{\mathit{\boldsymbol{s}}_\mathit{\boldsymbol{x}}})}}{{\partial {t^2}}} - {\nabla ^2}\mathit{\boldsymbol{u}}_{\rm{s}}^{\rm{H}}(\mathit{\boldsymbol{x}},t;{\mathit{\boldsymbol{s}}_\mathit{\boldsymbol{x}}}) = \mathit{\boldsymbol{f}}_{\rm{s}}^{\rm{H}}\left( {t,{\mathit{\boldsymbol{s}}_\mathit{\boldsymbol{x}}}} \right)\\ {{\mathit{\boldsymbol{\tilde u}}}_{\rm{s}}} = {\mathit{\boldsymbol{u}}_{\rm{s}}} + {\rm{i}}\mathit{\boldsymbol{u}}_{\rm{s}}^{\rm{H}} \end{array} \right. $

(3)

式中:f_s为地震子波震源项; t为时间; s_x=(s_x, s_y)为震源位置; x=(x, y, z)为空间位置; v为地下介质速度; u_s是炮点端波场, 即解析时间波场的实部; f_s^H为f_s的Hilbert变换; u_s^H为u_s的Hilbert变换; 即解析时间波场的虚部; ${\mathit{\boldsymbol{\tilde u}}_{\rm{s}}}$为解析时间波场。获得解析时间波场后, 可以在逆时外推的每个时间片上对解析时间波场施加空间傅里叶变换实现波场的角度方向分解。

由于在局部的空间窗内仅包含几个局部的平面波, 因此可以采用局部波场的稀疏分解方法来高效提取波场的局部角度分量^[47]。在实现局部平面波稀疏分解时, 定义角度分解的误差泛函为:

$J = {\left\| {{\mathit{\boldsymbol{u}}^{{\rm{win}}}}\left( {t, \mathit{\boldsymbol{x}}} \right)-\sum\limits_\theta {{\mathit{\boldsymbol{u}}^{{\rm{win}}}}} \left( {t, \mathit{\boldsymbol{x}}, \theta } \right)} \right\|^2} $

(4)

其中, u^win(t, x)代表局部解析空间波场, u^win(t, x, θ)代表提取的局部平面波, θ为局部平面波的传播角度。实现波场稀疏角度分解的过程为:首先对局部解析时间波场实施空间维的傅里叶变换, 并投影到波数域, 然后在波数域中将同一传播角度的波数进行反傅里叶变换, 实现局部平面波的分解。通过使目标泛函(4)最小化, 可将局部窗内的全部平面波提取出来。

对稀疏角度分解后提取的炮点端波场和检波点端波场的角度分量应用角度域成像条件, 提取成像结果, 得到角度道集。采用解析时间波场外推与局部波场稀疏分解的方法不需要进行关于时间维的傅里叶变换, 避免了整个外推波场的存储, 因此在逆时偏移中可以很方便地实现。同时该方法可以处理复杂的波前交叉情况, 得到的角度道集物理意义更加明确, 能够比较好地反映地下界面的AVA特征。但是该方法需要对解析波场虚部进行波场外推, 额外增加的波场外推计算量非常大, 特别是在各向异性介质中, 实现波场外推需要付出巨大的计算代价。

3 混合方法

波矢量方向估计方法可以高效地实现, 但是在检波点端存在方向估计不稳定的问题, 并且无法处理波场交叉的情况。利用局部平面波分解方法可以准确地提取波传播方向, 即使在波场交叉、检波点波场比较复杂的情况下, 仍然能够很好地提取波传播方向。但是该方法需要进行高维傅里叶变换, 需要大量内存和I/O, 效率低下, 即使采用解析时间波场, 也需要对震源和波场记录实施Hilbert变换, 分别进行波场外推, 一次逆时偏移至少需要做4次正演, 计算量还是比较大。本文结合局部平面波场分解的准确性和波矢量方向估计方法的高效性, 提出了一种混合生成RTM-ADCIGs的方法, 其实现过程如下:

1) 在正演和逆时波场外推的每个时间片上, 选择局部波场进行空间维傅里叶变换, 将局部波场从空间域变换到波数域。

2) 在波数域采用HU等^[47]提出的稀疏波场分解方法, 按照能量谱大小提取局部平面波(由于在波数域中能量谱是对称的, 因此只需提取正波数方向), 此时可以得到该平面波的两个波前方向φ或-φ(φ与-φ相差180°)。由于缺乏时间维的信息, 因此无法确定该波前传播方向是φ还是-φ, 但是提取的局部波前矢量具有稳健性, 其方向是正确的。

3) 利用坡印廷矢量方法计算局部波场矢量方向, 此时计算的波场方向在波场交叉和复杂情况下是不准确的。

4) 在选取的局部平面波空间窗内, 沿着垂直于该平面波的波前矢量方向选择带状的空间窗。在这个带状的空间窗内计算的波矢量方向和平面波的波前方向矢量大体一致。因此, 在选择的带状空间窗内统计波矢量方向落在计算的两个波前方向φ或-φ的个数, 统计个数最多的方向就代表了波前的方向(可选择一定的角度统计范围), 利用统计得到的波前方向来约束第2步中提取的波矢量方向的符号, 从而得到该平面波准确的波前扩散方向。空间窗的宽度以能覆盖有效波矢量方向为准则, 一般为局部平面波宽度的1/3~1/2除以局部平面波个数。角度统计范围根据具体的平面波个数限定, 本文选用的范围为30°除以局部平面波个数。

5) 利用提取的炮点端局部平面波和检波点端局部平面波, 应用角度域成像条件提取该时刻的ADCIGs。

6) 重复上述2)~5)步可提取下一个局部平面波, 直到空间窗内没有平面波为止。

下面通过一个简单的模型来论述角度道集的生成过程。如图 1a所示, 该模型从上往下速度值依次为3000, 3200, 3400m/s。对应1s时刻的波场快照如图 1b所示。选取该时刻某一局部空间波场如图 2a所示。采用HU等^[47]提出的解析时间波场方法计算的波场虚部如图 2b所示。将图 2a所示的局部空间波场进行傅里叶变换, 投影到波数域, 其能量谱如图 3a所示。将图 2a和图 2b组合形成的解析时间波场变换到波数域, 其能量谱如图 3b所示。由于解析时间波场中只包含正频率的信息, 两个局部平面波只含有两个波数域能量团, 依次提取对应的最大能量团即可获得该空间窗内的局部平面波方向。而采用常规波场得到的波数域能量谱中两个局部平面波包含4个能量团, 无法直接根据波数域能量团确定局部平面波的传播方向。分析发现, 图 3a中的波数域能量团是关于中心点对称的, 同一局部平面波的两个波数能量团的方向相反, 其夹角为180°, 如图 3a中红色箭头和蓝色箭头所示。利用对称性, 只在正波数区域内采用稀疏局部平面波分解方法提取最大波数能量团对应的局部平面波, 如图 4a所示。由于缺乏时间维信息, 一个局部平面波在波数域中对应两个能量团, 无法确定该局部平面波的传播方向, 只能得到两个可能的传播方向, 如图 4a中红色实箭头和黑色虚箭头所示。重复上述过程提取该区域内的第二个局部平面波如图 4b所示, 两个可能的传播方向分别为图 4b中蓝色实箭头和黑色虚箭头所示。图 5为采用坡印廷矢量方法提取的局部波场方向, 可以看到, 在波场交叉的区域, 波矢量方向的估计是错误的, 如图 5中橙色圆圈所示。根据图 4a中提取的局部平面波确定带状的空间窗如图 5中红色线框所示, 根据图 4b中提取的局部平面波确定带状的空间窗如图 5中蓝色线框所示。在红色线框内统计波矢量方向, 落在图 4a红色实箭头和黑色虚箭头所指方向的个数多的方向即为该局部平面波的传播方向(图 4a中红色实箭头所指的方向)。同样, 统计蓝色方框内波矢量方向落在图 4b中蓝色实箭头和黑色虚箭头所指方向的个数, 根据统计结果判断图 4b中蓝色实箭头方向为其局部平面波的传播方向。

接下来采用图 6所示的层状模型测试本文混合方法提取波场方向的有效性。正演采用的观测系统参数如下:震源位置为(3000m, 0), 震源为30Hz的雷克子波, 检波点间距10m, 偏移距范围为(-3000m, 3000m), 时间采样点数n_t=5000, 时间采样间隔Δt=0.5ms。检波点端波场在0.505s时刻采用不同方法获得的波场方向如图 7所示, 其中图 7a为采用坡印廷矢量方法计算的波场方向(为了与局部平面波显示一致, 只选取了局部空间窗中心点的矢量方向代表该空间窗内波场方向); 图 7b为对常规波场进行空间维傅里叶变换, 根据波数域能量谱大小获得的局部平面波波场方向; 图 7c为在图 7b基础上施加坡印廷矢量约束后计算的波场方向(即本文提出的混合方法); 图 7d为采用解析时间波场方法获得的波场方向。对比图 7a至图 7d可见, 只采用坡印廷矢量方法会在波场分解面处出现方向计算不准确的情况(图 7a); 而对常规波场进行空间维傅里叶变换方法由于缺乏时间维的信息, 也无法得到准确的波场方向(图 7b); 施加坡印廷矢量方向约束后, 可以获得准确的波场方向(图 7c), 其精度也与解析时间波场方法提取的波场方向相当(图 7d)。

混合方法的角度道集估计过程还可以扩展到不同类型的方法组合, 比如局部平面波的提取可以选择傅里叶变换、局部倾斜叠加和τ-p变换等方法, 波矢量方向的计算方法也可扩展到极化矢量方法、瞬时波数矢量方法以及光学流方法等。同时, 该方法很容易扩展到各向异性介质中, 与解析时间波场方法相比, 增加的计算量更少。

4 改进的混合方法

分析地下介质不同时刻的波场快照可知, 不同时刻的波场复杂程度不同, 在波场简单情况下, 大部分区域由单一的波前面组成, 只在局部区域存在波前交叉情况; 随着波场复杂程度的增加, 地下波场的交叉情况也变得越来越复杂, 即使如此, 在局部空间范围内也存在单一的波前面。因此, 可以根据波前的复杂情况, 选择合适的角度道集提取方法, 在波前简单的情况下采用快速的波矢量估计方法, 而在波前出现重叠、交叉的情况下采用局部波场分解方法。

考虑到波场在不同时空范围内的特征, 本文借助于结构张量特征值, 快速判断波场交叉的情况, 从而只在波场能量交叉的情况下使用混合方法, 估计一个或多个局部波场的方向, 而在其它区域采用快速的波矢量估计方法提取波场方向, 以进一步提高角度道集的计算效率。

结构张量, 又称二阶矩矩阵, 在图像处理中非常流行^[50-51], 也常用于地震数据分析。BAKKER^[52]对结构张量在地震数据滤波方面的应用做了详细的描述, HALE^[53]利用结构张量实现了地震剖面沿构造方向的平滑滤波处理, 增强了地下结构特征, 保持了断层结构特征。

以二维为例, 给定一波场u, 其结构张量矩阵可表示为:

$ \mathit{\boldsymbol{G}} = \nabla \mathit{\boldsymbol{u}} \cdot \nabla {\mathit{\boldsymbol{u}}^{\rm{T}}} = \left[{\begin{array}{*{20}{c}} {u_x^2}&{{u_x}{u_z}}\\ {{u_x}{u_z}}&{u_z^2} \end{array}} \right] $

(5)

其中, ▽u表示波场的梯度; T代表转置; u_x, u_z分别为x方向和z方向的方向导数。对结构张量矩阵进行特征值分解得:

$ \mathit{\boldsymbol{G}} = {\lambda _e}\mathit{\boldsymbol{e}}{\mathit{\boldsymbol{e}}^{\rm{T}}} + {\lambda _v}\mathit{\boldsymbol{v}}{\mathit{\boldsymbol{v}}^{\rm{T}}} $

(6)

其中, λ_e和λ_v分别是结构张量矩阵G的特征矢量e和v的特征值。结构张量矩阵特征值可以展开为:

$ {\lambda _e} = \frac{1}{2}(u_x^2 + u_z^2) + \frac{1}{2}\sqrt {{{(u_x^2-u_z^2)}^2} + 4{{({u_x}{u_z})}^2}} $

(7a)

$ {\lambda _v} = \frac{1}{2}(u_x^2 + u_z^2)-\frac{1}{2}\sqrt {{{(u_x^2-u_z^2)}^2} + 4{{({u_x}{u_z})}^2}} $

(7b)

结构张量矩阵特征矢量可以展开为:

$ \mathit{\boldsymbol{e}} = \left[{\begin{array}{*{20}{c}} {2{u_x}{u_z}}\\ {u_z^2-u_x^2 + \sqrt {{{(u_x^2-u_z^2)}^2} + 4{{({u_x}{u_z})}^2}} } \end{array}} \right] $

(8a)

$ \mathit{\boldsymbol{v}} = \left[{\begin{array}{*{20}{c}} {-u_z^2-u_x^2-\sqrt {{{(u_x^2 - u_z^2)}^2} + 4{{({u_x}{u_z})}^2}} }\\ {2{u_x}{u_z}} \end{array}} \right] $

(8b)

式中, λ_e为最大的特征值, 代表了第一个特征矢量e方向的张量能量, 特征矢量e的方向与梯度的方向一致; λ_v为最小的特征值, 代表了最小特征矢量v方向的张量能量, 特征矢量v代表了平行于梯度的方向。当波场中只存在单一波前面时, 利用结构张量计算的特征值只有λ_e有值, λ_v为零; 当波场中存在交叉的波前面时, 利用结构张量计算的特征值λ_e和λ_v均有值。因此, 可以利用λ_e和λ_v值的大小来判断波场是否存在交叉现象。为了方便起见, 对特征值λ_e和λ_v进行了归一化处理。图 8a为图 6所示层状模型1s时刻的波场快照, 图 8b为归一化后垂直于局部波前的特征矢量e对应的特征值λ_e, 图 8c为归一化后平行于局部波前的特征矢量v对应的特征值λ_v。可以看出, 在图 8a波前面单一的区域, 特征值λ_v为零, 在波前面存在交叉的区域, 特征值λ_v非零, 说明可以根据归一化后特征值λ_v的大小来判断是否存在波场交叉情况。

另外, 分析公式(1)和公式(7)可以发现, 利用公式(7)计算特征值时需要的u_x和u_z已在波矢量方向估计中进行了计算, 额外增加的计算量非常少, 因此可以在波矢量估计方法的框架中高效地实现波场交叉判断。

下面以局部Sigsbee模型(图 9)为例, 说明方法对不同复杂程度波场的交叉判断能力。图 10和图 11分别为Sigsbee模型1.6s时刻和3.2s时刻的波场快照、归一化后垂直于局部波前的特征向量对应的特征值λ_e和平行于局部波前的特征向量对应的特征值λ_v。对比图 10和图 11可以发现, 当波场比较简单时, 特征值λ_v全是零, 此时可用波矢量估计方法准确地计算波场的传播方向(图 10); 当波场比较复杂时, 特征值λ_v在波前单一区域为零, 在波场交叉区域非零, 于是可在零值区域采用波矢量估计方法计算波场方向, 在非零区域利用本文提出的混合方法计算波场方向, 从而提高了计算效率(图 11)。另外, 改进后的混合方法是一种自适应生成角度道集的方法, 在极端情况下, 即地下波场特别简单时, 可退化成波矢量估计方法; 当地下波场特别复杂时, 则转化为第3节介绍的混合方法。

5 数值实验

首先采用局部Sigsbee 2A模型(地表范围0~6000m)测试了本文混合方法的有效性(图 12)。检波点间隔10m, 偏移距范围为(-3000m, 3000m)。图 13a为正演计算的原始单炮炮集, 其加入信噪比(S/N)为5的高斯随机噪声后如图 13b所示。采用图 13b所示含噪声炮集计算角度道集, 提取其中的3道如图 14所示, 提取位置如图 12中蓝色竖线所示。图 14a为采用坡印廷矢量估计方法计算的角度道集, 图 14b为采用本文局部平面波分解+坡印廷矢量混合的方法得到的角度道集, 图 14c是采用结构张量特征值约束改进后的混合方法得到的角度道集, 图 14d是采用解析时间波场得到的角度道集。上述4种方法对应的0~60°角度道集叠加结果如图 15所示。由图 14a可以看到, 由于随机噪声的加入, 坡印廷矢量方法估计的方向不准确, 角度道集中含有较多噪声, 模糊情况严重, 对应的叠加成像结果(图 15a)也非常差, 整个剖面淹没在噪声中, 有效同相轴不易识别。而采用本文混合方法计算的角度道集(图 14b和图 14c)质量与采用解析波场计算的角度道集(图 14d)质量相当, 改进后的混合方法产生的角度道集与混合方法产生的角度道集质量也基本一致, 只在局部区域质量略低于混合方法。0~60°角度道集叠加结果也说明本文混合方法及改进方法具有很好的抗噪性。

分析本文混合方法及改进方法计算的角度道集(图 14b和图 14c)可知, 位于炮点下方的角度道集主要集中在小角度范围, 说明该炮下方的照明能量主要集中在小角度位置; 随着地表偏移距的增大, 该炮对地下层位照明的入射角越来越大。说明本文混合方法生成的角度道集可以准确地反映地下层位的角度照明情况, 为后续的角度照明分析等奠定基础。

分析不同角度道集生成方法的单炮计算时间, 如图 16所示。可以看出, 坡印廷矢量方法计算效率最高, 解析时间波场方法效率最差, 本文提出的混合方法与解析时间波场方法相比, 计算效率有所提升, 改进的混合方法在此基础上有进一步提升。综合效率和质量看, 本文改进的混合方法能在满足一定精度的基础上, 高效地实现逆时偏移角道集的生成。

接下来采用不同方法计算了图 17所示Sigsbee 2A模型右侧(地表范围6000~10000m)的角度道集。该模型的观测系统如下:炮点间距10m, 检波点间距10m, 炮点范围0~4000m;每一炮的检波器接收排列方式为震源两侧各2000m。图 18a~图 18d分别是坡印廷矢量方法、本文混合方法、改进的混合方法和解析时间波场方法每隔10个CDP点提取的角度道集, 其中角度范围为0~60°。由于坡印廷矢量方向在复杂区域计算不准确, 无法获得正确的角度值, 导致图 18a中角度道集出现较多噪声, 且分布间断、不连续。而本文提出的混合方法及改进方法(图 18b和图 18c)计算的波场方向准确, 角度道集连续性较好, 剖面噪声也少, 成像质量基本与解析时间波场方法一致(图 18d)。

对图 6所示水平层状模型不同方法生成的角度道集进行了校正^[27], 得到的AVA关系(浅层)如图 19所示。其中蓝色实线为利用佐普利兹方程近似得到的角度道集随角度变化的振幅值, 红色十字线为采用坡印廷矢量方法计算的角度道集随角度变化的振幅值; 蓝色圆圈线为采用解析时间波场计算的角度道集随角度变化的振幅值; 绿色五星线为采用本文混合方法计算的角度道集随角度变化的振幅值。可以看到, 随着角度的增大, 振幅值也不断增加, 采用混合方法得到的AVA关系基本上符合实际情况。这说明简单介质和简单波场情形下计算的角度道集振幅保真性是有保障的。在复杂介质情况下, 需结合保真RTM成像才能产生保真的角度成像道集。

6 结论

本文结合波矢量方向估计方法和局部平面波分解方法的特点, 提出一种混合构建角度道集的方法。首先利用局部波场分解方法获得波前方向, 然后利用坡印廷矢量方法确定波的传播方向, 最后采用角度域成像条件提取成像道集, 从而高效地获取高质量角度道集。该方法对波场方向的计算比较准确, 具有一定的抗噪性, 可解决波场交叉等问题, 并且很容易扩展到各向异性介质。在此基础上, 利用结构张量的特征值快速判断波场交叉情况, 只在波场能量交叉情况下使用混合方法, 而在其它区域采用快速的波矢量估计方法提取波场方向, 进一步提高了角度道集的计算效率。然而, 由于地下介质的复杂性、叠前地震数据采集时观测系统的不完备以及地下照明角度的非均匀性等因素影响, 真正产生保真的角度道集还非常困难, 需要进行深入的研究。

致谢: 感谢中国石油天然气股份有限公司勘探开发研究院及西北分院、中海油研究总院和湛江分公司、中国石油化工股份有限公司石油物探技术研究院和胜利油田分公司对波现象与反演成像研究组(WPI)研究工作的资助与支持。