在地震数据信噪比低的情况下, 建立准确的背景速度很困难, 而透射波信噪比高, 受道间时差的影响小, 只要观测方式合适, 利用透射波旅行时层析可以建立比较准确的表层(0~200m)、浅层(200~500m)甚至中层(500~1500m)的速度模型。对于面向高分辨和高保真的精确成像, 即使浅层速度只存在较小的相对误差, 也会导致中深层成像达不到同相位叠加, 从而严重影响高分辨率和高保真成像。

如今, 地震数据采集技术进步明显, “两宽一高”(宽频带、宽方位、高密度)地震数据采集逐渐成为主流, 研究与采集技术相匹配, 能够充分挖掘地下介质弹性参数信息的地震数据分析方法越来越重要。特别是, 在“两宽一高”地震数据采集中, 偏移距的增大使透射波信息更加丰富, 除了来自浅层的透射波, 中深层的透射波信息也能被采集到, 从而使透射波具备了反演中深层弹性参数的能力。宽方位采集方式还能优化透射波的照明范围, 降低了反问题的病态性。而高密度的地震采集则提供了更多的数据, 配合适当的数据预条件方法, 可以很好地降低层析核函数的稀疏性, 使层析反演结果更加稳定。同时, 相比反射波层析成像算法, 透射波层析更容易实现自动化, 对工业生产更具现实意义。因此非常有必要发展一套高精度的透射波层析反演方法与技术流程。

利用透射波信息进行层析成像, 属于数据域反演方法, 根据利用的数据属性, 透射波层析成像技术分为走时层析、相位层析和波形层析。与波形层析相比, 走时扰动对速度扰动的非线性程度较弱, 对初始模型的要求较低, 走时层析在实际生产中得到广泛应用。

按正问题的不同, 透射波层析则可分为射线层析、波束层析和波动方程层析。射线层析基于高频射线理论, 对初始模型精度的要求较低, 但反演精度不高^[1]。由于正问题基于射线理论, 受制于焦散、多路径和阴影区等射线理论的固有缺陷, 且层析矩阵十分稀疏, 病态性强, 求解时收敛较慢, 反演结果受射线照明的影响严重^[2]。波动方程层析反演方法中, 由于采用波动方程作为正演算子, 因而可以同时预测地震数据的运动学和动力学信息。经典的数据域波动方程层析方法为全波形反演方法^[3-4]。但是, 正是因为采用波动方程计算子波的正传播波场与残差的反传播波场, 即两次正演计算, 所以计算量较大, 加之海量的采集数据、复杂地区实际资料的信噪比低导致波形匹配困难等原因使得波动方程层析基本处于理论研究阶段, 很少能够得到实际应用。

可见, 射线层析与波动层析是层析反演的两个极端, 前者效率高但其效果受射线相关假设的制约, 后者理论上完善但效率及初始模型限制了方法的实际应用。射线层析明显与波传播规律不符, 可以从直观上改进射线层析, 引入波束的概念^[5-8], 采用Beam束模拟波传播路径, 将走时残差投影到Beam束路径上, 这样做在不增加太多走时层析计算量的基础上可以大幅度降低射线“盲区”数量及范围。更重要的是, Beam束层析摆脱了高频近似的射线理论, 考虑了波的频率信息, 使得走时层析更加稳定, 分辨率与精度也得到了显著提高。因此, 本文主要论述如何利用Beam束层析成像方法从透射波旅行时中挖掘与弹性参数相关的信息。

本文从层析成像的理论框架出发, 系统论述了波动方程层析、波束层析以及射线层析的异同, 剖析了参数迭代公式中各项的含义, 总结出利用透射波旅行时层析成像方法挖掘“两宽一高”地震数据信息时应注意的若干核心问题(层析正问题、数据的预条件与模型的正则化), 并重点阐述了加权走时的含义和自动化的获取加权地震走时信息的方法。最后, 论述了Beam层析算法理论以及对观测系统的适应性, 结合MPI+OpenMP策略在某探区实际资料上进行了应用, 其处理结果说明了透射波旅行时Beam层析处理流程在“两宽一高”海量地震数据处理中的可行性。

1 透射波层析成像的理论框架

基于地表观测的地震波数据估计地下介质的弹性参数进而进行精确的油藏描述是勘探地震的核心问题。这是一个病态的、强非线性的问题, 求解这样的问题需要有充分体现其本质的高度抽象的理论框架。地震数据是随机的, 确定性的反(散)射波场(同相轴)分布在满足不同统计特征的随机噪声中。基于不完美的地表观测数据和不完美的地震波预测模型进行Bayes框架下的地震波反演成像就成了必然的理论选择。Bayes估计理论框架下透射波反演成像方法的理论基础描述如下。

1.1 透射波层析反问题数学表达与求解

基于Bayes反演理论^[9-10], 在L2范数下定义误差泛函为:

$\begin{array}{l} J(\mathit{\boldsymbol{m}}) = \frac{1}{2}\left\| {\mathit{\boldsymbol{K}}(\mathit{\boldsymbol{m}}) - {\mathit{\boldsymbol{d}}_{{\rm{obs}}}}} \right\|_D^2 + \frac{1}{2}\left\| {\mathit{\boldsymbol{m}} - {\mathit{\boldsymbol{m}}_{{\rm{prior}}}}} \right\|_M^2\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\; \frac{1}{2}{\left[ {\mathit{\boldsymbol{K}}(\mathit{\boldsymbol{m}}) - {\mathit{\boldsymbol{d}}_{{\rm{obs}}}}} \right]^{\rm{T}}}\mathit{\boldsymbol{C}}_D^{ - 1}\left[ {\mathit{\boldsymbol{K}}(\mathit{\boldsymbol{m}}) - {\mathit{\boldsymbol{d}}_{{\rm{obs}}}}} \right] + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{1}{2}{\left( {\mathit{\boldsymbol{m}} - {\mathit{\boldsymbol{m}}_{{\rm{prior}}}}} \right)^{\rm{T}}}\mathit{\boldsymbol{C}}_M^{ - 1}\left( {\mathit{\boldsymbol{m}} - {\mathit{\boldsymbol{m}}_{{\rm{prior}}}}} \right) \end{array} $

(1)

(1)式为经典数据域FWI的目标泛函。式中:K(m)是依赖模型参数m的非线性函数(正问题), 是模型空间M向数据空间D映射的有界非线性算子; C_D表示与数据相关的协方差矩阵; C_M则是与模型相关的协方差矩阵; m_prior为先验模型信息; d_obs为观测数据。

通常而言, 上述泛函采用局部最优化方法求解。若给定初始模型m_k-1, 牛顿迭代算法写为:

$\begin{array}{l} {\mathit{\boldsymbol{m}}_k} = {\mathit{\boldsymbol{m}}_{k - 1}} - {\left( {\mathit{\boldsymbol{C}}_D^{ - 1}{\mathit{\boldsymbol{H}}_{{m_{k - 1}}}} + \mathit{\boldsymbol{C}}_M^{ - 1}} \right)^{ - 1}}\left\{ {\mathit{\boldsymbol{F}}_{{m_{k - 1}}}^{\rm{T}}} \right. \cdot \\ \mathit{\boldsymbol{C}}_D^{ - 1}\left[ {\mathit{\boldsymbol{K}}\left( {{\mathit{\boldsymbol{m}}_{k - 1}}} \right) - {\mathit{\boldsymbol{d}}_{{\rm{obs}}}}} \right] + \mathit{\boldsymbol{C}}_M^{ - 1}\left( {{\mathit{\boldsymbol{m}}_{k - 1}} - {\mathit{\boldsymbol{m}}_{{\rm{ prior }}}}} \right)\} \end{array} $

(2)

式中:${\mathit{\boldsymbol{F}}_{{m_{k - 1}}}} = {\left. {\frac{{\partial \mathit{\boldsymbol{K}}(\mathit{\boldsymbol{m}})}}{{\partial \mathit{\boldsymbol{m}}}}} \right|_{{m_{k - 1}}}}$为波场对模型的Fréchet导数; ${\mathit{\boldsymbol{H}}_{{m_{k - 1}}}} = {\left. {\frac{{{\partial ^2}\mathit{\boldsymbol{K}}(\mathit{\boldsymbol{m}})}}{{\partial {\mathit{\boldsymbol{m}}^2}}}} \right|_{{m_{k - 1}}}}$称为Hessian算子, 一般由线性和非线性两项构成。若假定目标函数在初始模型m_k-1处局部为二次型函数, 则可以对Fréchet微商引入线性化近似${\mathit{\boldsymbol{F}}_{{m_{k - 1}}}} = {\left. {\frac{{\partial \mathit{\boldsymbol{K}}(\mathit{\boldsymbol{m}})}}{{\partial \mathit{\boldsymbol{m}}}}} \right|_{{m_{k - 1}}}} \approx {\mathit{\boldsymbol{K}}_{{m_{k - 1}}}}$, 同时, 仅考虑Hessian的线性部分, 此时, 迭代算法为高斯-牛顿反演算法:

${\mathit{\boldsymbol{m}}_k} = {\mathit{\boldsymbol{m}}_{k - 1}} - {\left( {\mathit{\boldsymbol{K}}_{{m_{k - 1}}}^{\rm{T}}\mathit{\boldsymbol{C}}_D^{ - 1}\mathit{\boldsymbol{K}}_{k - 1}^m + \mathit{\boldsymbol{C}}_M^{ - 1}} \right)^{ - 1}}\left\{ {\mathit{\boldsymbol{K}}_{{m_{k - 1}}}^{\rm{T}}\mathit{\boldsymbol{C}}_D^{ - 1} \cdot }\\ \;\;\;\;\;\;\;\;\;\right.\left[ {\mathit{\boldsymbol{K}}\left( {{\mathit{\boldsymbol{m}}_{k - 1}}} \right) - {\mathit{\boldsymbol{d}}_{{\rm{obs}}}}} \right] + \mathit{\boldsymbol{C}}_M^{ - 1}\left( {{\mathit{\boldsymbol{m}}_{k - 1}} - {\mathit{\boldsymbol{m}}_{{\rm{prior}}}}} \right) $

(3)

大多数情况下, Hessian的逆(K_mk-1^TC_D^-1K_mk-1+C_M^-1)^-1很难求解, 一般我们并不用(3)式来更新参数扰动, 而是忽略Hessian项, 采用最速下降法求解, 即:

${\mathit{\boldsymbol{m}}_k} = {\mathit{\boldsymbol{m}}_{k - 1}} - {\mu _{k - 1}}\left[ {{\mathit{\boldsymbol{C}}_M}\mathit{\boldsymbol{K}}_{{m_{k - 1}}}^{\rm{T}}\mathit{\boldsymbol{C}}_D^{ - 1}\delta \mathit{\boldsymbol{d}} + \left( {{\mathit{\boldsymbol{m}}_{k - 1}} - {\mathit{\boldsymbol{m}}_{{\rm{prior}}}}} \right)} \right] $

(4)

式中:μ_k-1为最速下降迭代步长。

当δd表示地震采集数据中的透射波数据残差时, 公式(4)为常用的透射波反问题的迭代求解公式, 其物理意义为对透射波数据残差的加权处理结果(数据预条件), 通过由模型正则化项C_M作用后的K_mk－1^T算子反投影到模型空间, 用其修正慢度。需要指出的是, 在实际应用中, 我们一般简单地选择初始模型m₀=m_prior, 先验信息一般来自于重磁电信息、构造解释层位信息、测井信息以及近地表调查信息等。

1.2 透射波数据残差与预条件(C_D^－1δd)

如公式(4)所述, 我们根据不同的需求, 灵活选择不同的透射波属性(δd)进行反传, 属性可以是走时、相位、振幅或波形等。理论上, 为了充分挖掘透射波携带的信息, 应该从透射波走时反演出发, 逐渐过渡到波形反演。结合前文分析, 本文重点论述如何充分挖掘透射波走时(δt)中携带的地下弹性参数信息。

C_D^－1为数据的预条件项, 可看作是数据协方差的逆, 本质上是对数据进行变换, 也可以选一个基函数特征表达数据, 实质上是解开矩阵向量各分量之间的相关性, 得到稀疏的特征表达。根据C_D^－1不同的数学表达形式, 其物理含义十分丰富。HU等^[2]认为其为照明均衡的算子, 最终使反传的梯度分布更加均匀。刘玉柱等^[11]将其定义为偏移距加权算子, 数据的权重系数是偏移距的光滑函数, 可以区别利用所有偏移距数据, 实现由浅到深的逐步反演。ZHANG等^[12]将C_D^－1定义为基于相似系数的数据筛选矩阵, 在每一次迭代过程中, 通过度量观测数据与模拟数据的相似性, 保留相关性较强的数据残差进行反传, 逐次迭代克服FWI中的周期跳跃现象。可见, 在选定某种透射波属性后, 引入C_D^－1对数据进行预处理, 能显著地提升透射波层析反问题的精度与收敛效率。

针对透射波旅行时层析反演, 数据的不确定性应该综合考虑初至波同相轴的空间相关性、信噪比、子波一致性等因素确定到达时的可靠性程度, 可靠性大的到达时, 由C_D^－1引入的加权系数大; 可靠性小的到达时, 由C_D^－1引入的加权系数小。这是数据预条件最基本的处理方法。

1.3 透射波旅行时层析正问题(K(m))

为了求解公式(4), 必须研究透射波走时层析的正问题t=K(m)。最经典的透射波层析方法是波形层析方法, 譬如初至波或早至波波形层析。低频长偏移数据下的FWI本质上实现了波形层析, 叠前数据Laplace变换下的FWI是一种典型的做法^[13-14], 通过调节Laplace变换中的系数控制参与层析反演的早至波成分的多少。如此进行透射波波形层析的优势是基本不需要人工干预, 但初至波及后续一段时间内的波现象非常复杂, 成因不明, 加上初始模型不准, 数据中缺乏低频成分, 透射波波形层析很难收敛。因此, 针对陆上数据, 最好是进行波动理论的初至波旅行时层析, 如果早至波到达时检测问题解决得好, 还可以进行波动理论的早至波旅行时层析。

建立透射波旅行时扰动与速度扰动之间的关系是透射波速度层析反演的最关键环节。

通过波场的Rytov近似可以导出旅行时扰动与速度扰动之间的线性关系, 得到层析反演的线性矩阵方程, 然后在最小二乘意义下进行层析反演。冯波基于波动方程的Rytov近似并采用模拟数据的归一化能谱作为单频走时扰动的加权函数, 导出带限地震波的走时敏感度核函数^[15-16]:

${\mathit{\boldsymbol{K}}_{{\rm{ wave }}}} = \frac{2}{{{v^3}(\mathit{\boldsymbol{x}})}}\int {{u_0}} \left( {\mathit{\boldsymbol{x}}, t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)p\left( {\mathit{\boldsymbol{x}}, t;{\mathit{\boldsymbol{x}}_{\rm{r}}}} \right){\rm{d}}t $

(5)

从(5)式可以看出, 对于任意一对炮检(x_r, x_s), 均需计算一次波场正演(计算u₀(x, t; x_s))以及一次波场逆时传播(计算p(x, t; x_r), 因此计算量远高于射线追踪或波束层析。

相应地, 地震波走时层析正问题矩阵形式可表示为:

${\mathit{\boldsymbol{K}}_{{\rm{ wave }}}}\delta \mathit{\boldsymbol{v}} = \delta \mathit{\boldsymbol{t}} $

(6)

式中:δt为透射波的走时残差; δv为速度参数扰动。

从程函方程的角度也可以建立速度扰动和旅行时扰动之间的关系进行射线理论的透射波旅行时层析反演。为此, 首先给出波动方程的平面波解:

$p\left( {{x_i}, t} \right) = P\left( {{x_i}} \right)\exp \left\{ { - {\rm{i}}\omega \left[ {t - T\left( {{x_i}} \right)} \right]} \right\} $

(7)

式中:p(x_i, t)为空间某一点x_i在t时刻的压力场; P(x_i)为振幅项; T(x_i)为旅行时; ω=2πf为圆频。(7)式也是波动方程WKBJ近似的零阶解。将(7)式回代到声波波动方程中, 得到:

$\begin{array}{*{20}{c}} { - {\omega ^2}P\left[ {{{(\nabla T)}^2} - \frac{1}{{{v^2}}}} \right] + {\rm{i}}\omega (\nabla P \cdot \nabla T + }\\ {P{\nabla ^2}T) + {\nabla ^2}P = 0} \end{array} $

(8)

当圆频ω→∞(高频近似)时, (8)式中第三项可以忽略, 同时为了使公式(8)成立, 需要满足如下两个方程:

${(\nabla T)^2} = {\left( {\frac{{\partial T}}{{\partial x}}} \right)^2} = \frac{1}{{{v^2}}} $

(9)

$2\nabla P \cdot \nabla T + P{\nabla ^2}T = 0 $

(10)

公式(9)和公式(10)便是波动方程在高频近似下推导出来的程函方程和输运方程。公式(9)是关于旅行时T(x_i)的一阶非线性偏微分方程, 是射线旅行时理论的基础; 公式(10)则为振幅项P(x_i)的一阶线性偏微分方程。

求解程函方程可以进行射线追踪, 建立起如下的沿射线路径的速度和旅行时之间的积分关系:

${t^i} = {g^i}(\mathit{\boldsymbol{s}}) = \int_{Ri(s)} {\rm{d}} T = \int_{Ri(s)} {\rm{d}} {l^i}\left( {{\mathit{\boldsymbol{x}}^j}, \mathit{\boldsymbol{s}}} \right)s\left( {{\mathit{\boldsymbol{x}}^j}} \right) $

(11)

式中:s表示慢度场模型; s(x^j)表示在点x^j处的慢度值, 与速度场的关系为s(x^j)=v^－1(x^j); dlⁱ(x^j, s)表示在慢度模型s中第i根射线在网格x^j处的射线长度; Rⁱ(s)代表在模型s中第i根射线的射线路径; tⁱ是第i根射线旅行时。(11)式属于Fredholm第一类积分方程(Radon变换也属于此类方程)。

在此基础上, 引入介质慢度扰动量δs, 并对旅行时进行了Taylor展开并只保留一阶项, 假设介质的小扰动不会引起射线路径太大的变化, 即假设射线在光滑介质中穿过, 可以进一步得到慢度扰动量和旅行时扰动量之间的近似线性关系:

$\delta {t^i} = {g^i}(\mathit{\boldsymbol{s}} + \delta \mathit{\boldsymbol{s}}) - {g^i}(\mathit{\boldsymbol{s}}) = \int_{{R^i}(s)} {\rm{d}} {l^i}\left( {{\mathit{\boldsymbol{x}}^j}, \mathit{\boldsymbol{s}}} \right)\delta s\left( {{\mathit{\boldsymbol{x}}^j}} \right) $

(12)

公式(12)引入射线层析核函数矩阵K_Ray, 则射线层析正问题的矩阵方程形式可写为:

${\mathit{\boldsymbol{K}}_{{\rm{Ray}}}}\delta \mathit{\boldsymbol{s}} = \delta \mathit{\boldsymbol{t}} $

(13)

基于公式(13), 当引入波束的概念时, 传统的射线类走时层析核函数矩阵K_Ray和正问题的矩阵形式可重写为:

${\mathit{\boldsymbol{K}}_{{\rm{ Beam }}}}(l, q) = \frac{1}{{\sigma (l)\sqrt {2\pi } }}{{\rm{e}}^{ - \frac{{{q^2}}}{{2\sigma {{(l)}^2}}}}} $

(14)

${\mathit{\boldsymbol{K}}_{{\rm{ Beam }}}}\delta \mathit{\boldsymbol{s}} = \delta \mathit{\boldsymbol{t}} $

(15)

式中:K_Beam为引入波束概念后的核函数, 以射线交点为中心呈高斯衰减, 表示Beam-ray计算的波场振幅可以作为加权因子, 来提高透射波层析的反演精度; (l, q)是射线中心坐标系, l是沿射线的弧长, q代表垂直于射线的距离; σ是标准差, 与频率相关。这里需要指出的是, 根据对K_Beam的不同定义, 在保证沿垂直射线的椭圆(三维)或线段(二维)的积分为1的前提下, 可以形成不同形式的波束(Beam)。

结合公式(6)、公式(13)、公式(15)以及公式(4)可以看出, 从波动方程出发, 通过不同的假设近似, 透射波旅行时层析反演从波动理论逐渐退化到射线理论, 复杂程度逐步下降, 挖掘透射波走时能力也在递减, 其计算耗时也呈指数形式下降。因此, 折中选取研究波束走时层析, 在精度和效率两者间做个平衡, 更具有现实意义。

1.4 透射波层析模型正则化(C_M)

地震波反演问题是一个典型的病态问题, 但求解这一病态问题具有现实意义。求解病态反问题的基本思想是施加各种正则化方法, 使得反问题的解能解决实际问题。

从目前的实际情况来看, 促进层析成像技术走向实用化的关键之一就是反演过程中模型正则化的施加。本质上, 模型正则化是利用模型的先验信息约束模型估计本身及其分布范围, 降低模型估计的多解性。通过模型参数协方差(逆)算子对反演解进行正则化约束是非常根本的方法。MANDELLI等^[17]系统对比了不同模型正则化(倾角滤波、TV以及各向异性扩散)的差异, 认为把地质构造约束定量化地引入到层析成像中, 将增加反演的稳定性。KAZEI等^[18]则从另一个角度出发, 将模型正则化转化为惩罚函数选择问题, 从模型的L2范数问题过渡到Lp范数, 从而实现由低波数向高波数的过渡。ROMAHN等^[19]利用已知井信息构建模型正则化算子C_M, 相当于构建了一个匹配滤波器, 作用于反演速度扰动上, 达到用井信息修正梯度的目的。

在“两宽一高”地震数据采集基础下, 实用化的透射波旅行时层析成像技术核心问题有:旅行时(差)的拾取或计算; 旅行时差质量评价以及关于旅行时残差的数据预条件(C_D^－1δd); 观测系统、模型正则化与矩阵(核函数)形成(C_MK^T)以及大型稀疏矩阵的存储与求解。总之, 我们应该将精力聚焦于以上几点, 在正则化/预条件思想的约束下, 发展一套实用化的模型参数与数据之间的映射关系, 用以实现高精度的参数建模。基于此, 我们将重点从透射波到达时的(自动)测量和线性层析成像方程中的Kernel定义(Beam束)以及对观测系统的适应性三方面阐明透射波旅行时Beam层析成像算法流程。

2 透射波层析成像中到达时测量问题

地震波走时的拾取是透射波层析成像中重要的一步。一般意义下, 地震波走时拾取表达一个物理实体的到达, 数学上抽象为一个质点的到达。虽然数学上一个质点的到达概念很明确, 但是不具可测性, 因此, 实体的到达时必然引起不确定性。同理, 子波的到达时测量也存在不确定性。理论意义下, 在时间域中, 测理论脉冲的到达时是准确的, 测带限波形的到达时存在不确定性; 在频率域中, 用单频波的相位同相轴测时差是可行的, 而测单频波的到达时存在困难。传统的旅行时拾取算法大体可分为两类:滑动时窗方法和相干类方法。在滑动时窗方法^[20-21]中, 地震信号序列的属性在连续或重叠的移动窗口中计算。相干类方法^[22]依赖于使用一些相似度测量技术, 比如利用互相关来比较单个或多个波形, 从而得到走时信息。

本文在相速度和群速度的基础上, 定义了基于相位延迟和群延迟^[23]的两种加权带限子波到达时。并在文献[24]和文献[25]的基础上, 讨论了基于群延迟的自动走时拾取的迭代优化策略。通过引入一个滤波算子, 迭代拾取的策略能够自动拾取走时, 用于透射波层析反演算法。

2.1 相位延迟

在真空中, 所有波长的电磁波均以相同的相速度——真空中的光速传播。在色散介质中, 由于介质的折射率与光的波长(频率)有关, 不同波长(频率)的电磁波具有不同的相速度。程乾生^[23]给出了相应“相位延迟”的概念, 一般的滤波器频谱S(ω)为:

$S(\omega ) = A(\omega ){{\rm{e}}^{{\rm{i}}\Phi (\omega )}} $

(16)

式中:A(ω)=|S(ω)|为振幅谱; Φ(ω)为相位谱。滤波器的相位延迟T_p(ω)为:

${T_{\rm{p}}}(\omega ) = - \frac{1}{{2\pi }}\frac{{\mathit{\Phi} (\omega )}}{\omega } $

(17)

相位延迟表示频率为ω的单一正弦波经过滤波器后的延迟时间。如果将检波器接收到的信号看成若干频率的谐波经过地球滤波后被检波器接收, 则带限子波走时可以看作相位延迟的一个加权效应。XIE等^[26]指出带限地震信号的走时扰动可以用单频相位延迟扰动加权叠加表示:

${T_{\rm{p}}} = \int_{{\omega _{\min }}}^{{\omega _{\max }}} P (\omega ){T_{\rm{p}}}(\omega ){\rm{d}}\omega $

(18)

式中:P(ω)为加权算子。冯波^[15]详细论述了公式(18)的物理含义, 在此不再赘述。

2.2 群延迟及瞬时走时

实际存在的波不是单频的, 传播介质对这个(或这些)波必然是色散的, 那么, 传播中的波由于各不同频率的成分运动快慢不一致, 会出现扩散, 但假若(注意这个假设)这个波由一系列频率差别不大的简谐波组成, 这时在相当长的传播过程中整个波组仍将维持为一个整体, 以固定的速度传播, 这个特殊的波群称为“波包”, 这个速度称为群速度, 也称包络延迟。与相速度不同, 群速度比波包的中心相速度小, 并且二者的差值同中心相速度随波长而变化的平均率成正比。群速度是波包的能量传播速度, 也是波包所表达信号的传播速度。滤波器的群延迟T_g(ω)为:

${T_{\rm{g}}}(\omega ) = - \frac{1}{{2\pi }}\frac{{{\rm{d}}\mathit{\Phi} (\omega )}}{{{\rm{d}}\omega }} $

(19)

群延迟(包络延迟)反映某一频率的包络经过滤波器后的延迟时间, 体现了在频率ω的邻域内局部延迟性质。如果我们将检波器接收到的信号看成若干频率的地震波经过地球滤波后被检波器接收, 则带限子波走时可以看作群延迟在一定频率内的平均效应(平均时间), 可写为:

${T_{\rm{g}}} = \frac{{\int {{T_{\rm{g}}}} (\omega ){\rm{d}}\omega }}{{{\omega _{{\rm{max}}}} - {\omega _{{\rm{min}}}}}} $

(20)

LUO等^[27-28]阐明了瞬时旅行时的具体意义, 其等价于群延迟。在时-频域内瞬时走时τ(t)的表达形式为:

$\tau (t) = {M_\omega }\{ \tau (t, \omega )\} $

(21)

$\tau (t, \omega ) = {\mathop{\rm Im}\nolimits} \left\{ {\frac{{{\rm{d}}S(t, \omega )}}{{{\rm{d}}\omega }}/S(t, \omega )} \right\} $

(22)

${M_\omega } = \frac{{\Delta \omega \sum\limits_{\omega = {\omega _{\min }}}^{{\omega _{\max }}} \tau (t, \omega )}}{{{\omega _{\max }} - {\omega _{\min }}}} $

(23)

式中:S(t, ω)为信号s(t)的时频变换; Im{·}代表取虚部。(22)式中, 计算导数dS(t, ω)/dω可能不稳定, 因此采用FOMEL^[29]提出的平滑除法以减小除法带来的不稳定。

根据公式(21), 计算带限信号的瞬时走时可以分为互不干涉的两步:首先, 使用公式(22)将单道地震数据映射到时间-频率域; 然后, 通过算子M_ω(公式(23))从时频域重新映射回时间域。我们通过图 1和图 2进一步说明如何计算瞬时走时。图 1a展示的是由两个地震事件(0.32s, 1.28s)合成的某一地震信号。图 1b和图 1c分别为合成信号的时频谱S(t, ω)和采用公式(22)计算得到的时频图τ(t, ω)。最终计算的结果如图 2所示。其中图 2a表示将τ(t, ω)通过映射算子作用, 也就是在有效带宽内对频率平均后的结果, 即τ(t), 当τ(t)=t时, 其对应的走时t则为瞬时旅行时。而旅行时的自动拾取结果t位于0.32s和1.28s处, 如图 2b中脉冲位置所示(红线处)。

2.3 基于瞬时走时拾取的迭代优化策略

这里, 我们引入一个阻尼滤波算子:

$w(i, m) = - a\left\{ {t - \left[ {\tau (i, m) + \frac{1}{a}} \right]} \right\}{{\rm{e}}^{a\left\{ {t - \left[ {\tau (i, m) + \frac{1}{a}} \right]} \right\}}} $

(24)

$\tau (i, m) = \eta (i) + \frac{{{T_0} - \eta (i)}}{M}m $

(25)

$\begin{array}{*{20}{c}} {(m = 0, 1, \cdots , M - 1)}\\ {\eta (i) = \{ t|t \in {\rm{Method}}({\rm{STA}}/{\rm{LTA}}, {\rm{etc}}.)\} } \end{array} $

(26)

式中:w(i, m)为定义的迭代滤波算子, 该算子作用于每一道; a为阻尼因子, 决定了滤波器的宽度; T₀和M分别代表地震信号的采样时间和最大迭代次数; η(i)相当于一个粗略的初至走时。

这种方法可以分为以下几个步骤:①通过公式(26)计算每一道的η(i); ②为每一次迭代设计滤波器算子w(i, m)并作用于每一道上; ③利用公式(21)和公式(23)计算每一道的瞬时旅行时; ④重复步骤②和步骤③, 直到迭代结束。

2.4 数值试验

我们使用YILMAZ^[30]提供的Data6作为实际测试数据。该地震数据为炸药震源激发, 48道接收, 道间距为100m, 采样间隔为4ms。该实际数据存在多个接近初至波信号的来自更深层的折射波信号, 两次反射波信号约1.5s和2.7s。

我们分别采用Saragiotis的方法和基于瞬时走时拾取的迭代优化策略对该实际数据进行处理, 结果如图 3所示。对比图 3a与图 3b可以看出, 改进后的方法能够更加准确地拾取初至以及反射走时信息, 在此基础上, 配合一定的质量控制(QC)准则, 可以逐层选取合适的透射旅行时信息用于层析成像。

3 透射波旅行时Beam层析成像技术 3.1 Beam波路径模拟

当引入波束的概念时, 可以用公式(11)和(12)模拟Beam波路径。图 4显示了三维梯度速度场的Beam波路径模拟情况。该模型X方向3750m, 共301个网格采样点, 网格大小为12.50m;Y方向2500m, 共201个网格采样点, 网格大小为12.50m;深度h=1875m, 共301个网格采样点, 网格大小为6.25m。梯度速度模型可表示为:v(h)=600+(h/6.25)·3, 单位m/s。

将炮点放置在(100m, 800m, 25m)处, 检波点放置在(1625m, 800m, 25m)处。图 5为不同频率的三维Beam波路径, 可以看出, 沿着半径方向能量呈高斯衰减, 随着频率的增加Beam波半径变小, 当频率趋近于无穷大时, Beam束近似为射线。

3.2 三维透射波旅行时Beam层析与射线层析对比分析

为了具体说明Beam束相对于射线层析的优势, 采用如图 6a所示的速度模型进行数值试验。该模型在近地表附近存在两个高速异常体。该模型X方向3750m, 共301个网格采样点, 网格大小为12.50m;Y方向2500m, 共201个网格采样点, 网格大小为12.50m;深度为1875m, 共301个网格采样点, 网格大小为6.25m。背景梯度速度场可表示为:v(h)=600+(h/6.25)·3, 单位m/s, 高速异常体速度为1600m/s。每条炮线上有11炮, 总共16条炮线, 共计176炮, 炮线距为250m, 炮间距为250m, 道间距125m。在此基础上, 只选择满覆盖左侧异常体的4条线数据进行反演。

以图 6b所示的背景梯度场为初始速度场, 经过5次反演迭代的结果如图 7所示。图 7a为传统射线类三维透射波层析反演结果(未加入任何约束)。图 7b为三维Beam束层析反演结果。对比图 7a和图 7b可见, 三维Beam束层析结果分辨率更高。为了进一步说明效果, 抽取分别穿过两个异常体的纵向速度曲线进行对比, 结果如图 8所示, 可以看出, 经过5次迭代的Beam层析(绿线)结果明显优于经过5次迭代的Ray层析(红线)结果。

相对射线走时层析而言, Beam束走时层析在不增加太多计算量的基础上可以大幅度降低射线“盲区”数量; 同时, 该方法考虑了波的频率信息, 使得走时层析更加稳定, 分辨率与精度也得到了显著提高。在同等观测系统条件(不完备)下, Beam束参数反演纵横分辨率以及反演深度都比射线类方法要好。一般, Beam束层析耗时高于射线层析, 但远低于传统波动方程层析, 因此可以在生产上, 特别是处理“两宽一高”地震数据时发挥一定作用。

3.3 观测系统与透射波层析成像精度的关系分析

随着“两宽一高”采集的出现, 地震波中包含的信息越来越丰富。ZHOU等^[31]通过对实际资料的处理评价了透射波层析成像技术的能力, 如图 9所示, 采用不同偏移距数据进行透射波层析反演得到的速度模型精度存在一定差异。可见, 在充分利用中、深层折射波以及回转波信息(长/超长偏移距生成的初至/早至波)时, 透射波层析反演的精度能够得到一定程度的提升。

首先, 截取BP模型的部分数据, 采用不同道间距数据进行反演。图 10显示了真实速度模型以及初始速度模型。图 11为不同道间距的反演结果。图 12为横向速度的抽道对比结果。通过分析可知, 对道间距的基本要求是:反演目标越接近地表、异常体尺度越小, 道间距越小。同时, 考虑到成本以及效率的原因, 在保证较好反演分辨率和消除道间距时差的基础上, 可以适当选取较大的道间距进行反演。

接着, 设计如图 13a所示的宽方位观测系统(蓝线为炮线, 红线为检波线), 5条炮线, 炮线距200m, 炮线长度为16km, 每条炮线上有639炮, 炮间距25m, 因此存在长偏移距信息。37条检波线, 检波线间距为200m, 检波线长度为4800m, 检波点间距为50m, 相当于1炮激发, 3552道接收。并以图 13b所示真实速度场正演数据, 使用梯度速度场作为初始速度。通过对比不同形式地震数据的反传梯度的区别, 简要总结了宽方位观测系统对Beam透射波层析成像方法技术的影响。

图 14为二炮线与五炮线数据速度反传梯度的对比结果, 可以看出, 随着炮线的增加, 有效方位增加, 从而射线覆盖更好, 某些在二炮线时未反演出来的区域在五炮线情况下出现。图 15为不同偏移距情况下的反传梯度。可见, 随着偏移距的逐渐增加, 深部照明能量也在不断增加。

结合前面的数值分析结果可知, 透射波层析反演的精度与观测系统密切相关。在各向同性介质假设下, 宽方位的采集模式能够优化照明范围, 提升反演精度。同时, 需要采集长偏移距/超长偏移距数据, 这样可以加大反演深度, 得到更深层的反演结果。在此基础上, 道间距并非越小越好, 而是应该根据处理数据的规模、反演目标的尺度/深度以及生产成本等信息, 权衡分辨率与成本, 在高密度数据中抽取合适的道间距。总之, 合适的高密度宽方位长偏移距观测系统, 能够提高透射波层析成像的精度以及深度。

4 实际资料测试 4.1 瞬时旅行时拾取结果

海上某“两宽一高”工区的实际地震记录如图 16所示, 共639道, 道间距25m, 最大偏移距为13625m, 纵向采样4100个, 采样间隔为1ms。该数据在红框处存在严重的吸收衰减现象, 信噪比很低; 在绿框处由于偏移距较大导致深层折射波与初至波相邻很近, 严重影响了拾取精度。分别采用传统方法和本文迭代优化策略进行瞬时旅行时拾取, 其结果分别如图 16中的绿线和红线所示, 图 17a和图 17b分别为图 16中红框处和绿框处的放大显示。对比图 17a和17b可以发现, 采用本文方法能够得到更加精确的初至旅行时。

按第3节中介绍的方法, 半自动拾取了该工区共101309153道的瞬时旅行时。通过QC处理保留大约85961691道, 有效道数占比为84.85%。需要指出的, QC剔除的道, 主要集中在吸收衰减区域以及超长偏移距处(图 16), 这是因为采用半自动化(加入一些人为QC策略)算法时同样的参数可能不适合所有的数据, 所以拾取质量有差异。

4.2 反演结果

获取合适的透射波旅行时信息后, 我们从图 18a所示初始速度场出发, 采用透射波旅行时Beam层析算法进行了两次迭代, 其迭代更新量如图 18b所示, 图 18c为更新后的速度场, 共处理大约1.6TB数据。分析迭代更新量可以发现, 在速度场的东北方向存在一个高速异常界面, 与地震数据上出现远偏移距吸收衰减严重以及生成折射波现象相吻合。图 19为走时差对比结果(红色箭头指示为拾取走时与更新后速度场正演走时作差的结果, 而绿色箭头则为拾取走时与初始速度场正演走时作差的结果)。可见, 采用更新后的速度正演模拟走时的结果更趋近于拾取结果。分别对更新前、后的速度场进行偏移成像, 图 20为采用双平方根叠前深度偏移得到的X方向的偏移剖面切片, 对比图 20a和图 20b的绿色椭圆处可以看出, 采用更新后速度得到的偏移结果同相轴更加连续聚焦。结合迭代更新量的速度异常分布、走时差对比以及偏移结果, 说明透射波旅行时Beam层析速度反演有能力处理大规模的“两宽一高”地震数据。

4.3 大型稀疏矩阵的求解性能分析

该工区数据生成层析敏感度矩阵, 即使是稀疏存储, 仍然会占用计算机约378GB内存。显而易见, 直接采用传统LSQR算法求解, 需要单节点至少拥有400GB内存以上, 否则难以实现。因此, 本文采用MPI+OpenMP策略, 即以MPI模式将灵敏度矩阵存储分摊到多个节点上, 单节点则采用OpenMP模式加速计算, 以此降低对单节点内存的要求。我们使用了7个节点, 每个节点上分摊大约54GB内存。采用共轭梯度法求解, 共耗时约17.5h, 主要是数据IO以及多节点归约等待。

可见, 面对“两宽一高”采集数据, 特别是TB级规模的数据, 层析核函数的稀疏存储与求解问题仍然是工业生产需要面对的难题。一方面, 需要研究无需显式计算和存储核函数的求解层析反问题的算法, 比如隐式矩阵向量乘方法、伴随状态法^[32]等。这类算法将层析反问题转化为单道并行处理, 降低对计算机内存的要求。另一方面, “两宽一高”的地震数据为我们提供了更加充足完备的样本, 我们可以借鉴机器学习中的优化算法, 引入随机梯度的概念, 将地震数据分成若干样本子集, 每次只使用部分数据进行反演。总之, 如何在现有计算机条件下, 完成“两宽一高”地震数据的大规模层析成像处理, 是值得我们关注且具有现实意义的。

5 结论

复杂介质情况下, “两宽一高”采集的地震数据, 需要有针对性的参数估计方法。高密度采集加上炮集中透射波具有较高的信噪比, 利用透射波(初至波)进行层析成像对于提高浅中层速度(包括Q值)建模的精度具有重要意义。

透射特征波场的有效提取是利用透射波信息进行层析反演的关键环节。初至带限子波的到达时测量(计算瞬时走时), 应该分为四步:首先, 提取数据某种特征; 然后, 定义走时与该特征的关系; 再次, 度量观测数据和模拟数据该种特征的距离; 最后, 对拾取结果进行质量监控, 剔除精度较差的结果。同时, 随着地震数据量的日益增大, 孤立地讨论单道的地震事件提取而忽视数据的横向空间连续性, 这样没有充分利用数据的特征信息, 从而降低特征提取的精度。因此, 在机器学习框架下, 将地震数据拓展到高维空间, 基于多域、多属性特征的自动地震旅行时拾取技术, 或许是一个很好的思路。

相对射线走时层析而言, Beam束走时层析在不增加太多计算量的基础上可以大幅度降低射线“盲区”数量, 并且考虑了波的频率信息, 使得走时层析更加稳定, 分辨率与精度也得到了显著提高。数值试验结果表明, 相同情况下Beam束层析反演纵横分辨率以及反演深度都比射线类方法好。经过海上实际资料测试, 透射波(初至和早至波)Beam层析反演的深度并不局限于浅层, 随着采集偏移距的增加, 其反演深度将随之增加(达到浅层甚至中层)。采用宽方位数据进行透射波层析反演, 增加了波场的照明范围, 使反演解更加稳定。高密度采集则为我们提供了充足的数据信息, 允许我们灵活设计合适的求解策略。因此, 和透射波层析反演算法相匹配的观测系统应该是宽方位高密度长偏移距观测系统。同时, Beam束层析耗时远低于传统波动方程层析, 但仍高于射线层析。为了更好地在“两宽一高”采集数据中利用透射波信息(特别是TB级数据), 透射波旅行时Beam层析反演的精度与效率问题仍需继续探究。

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