全波形反演(FWI) 根据观测地震数据构建地下速度模型, 利用了地震波波形中的所有信息, 包括振幅和相位信息。常规的全波形反演^[1-3]和反射波波形反演^[4-5]都是根据最小化观测数据和正演数据之间的波形差异, 沿着波路径反传数据误差, 得到模型更新量。FWI最大的问题是相减型目标函数与速度模型呈强非线性关系。只要预测数据和观测数据中存在相位差异大于半个周期的情况, FWI就会遇到严重的周波跳跃问题。在这种情况下, 反演容易陷入局部极小值。从小偏移距到大偏移距进行FWI是解决周波跳跃问题的一种实用化方法; 另外, 从低频到高频进行FWI也能够缓解这一问题。为了避免观测波形和正演波形之间存在太大的旅行时差异, 仍然需要一个尽可能准确的初始速度模型。

相对于波形, 旅行时和长波长速度变化之间更符合线性关系, 因此工业界常采用旅行时层析反演地下结构信息。地震层析只利用了旅行时信息, 但是能够提供满足构造成像要求的速度模型, 也即运动学意义上正确的速度模型。射线层析是工业界最常用的旅行时反演方法, 该方法通过旅行时信息反演得到平滑的偏移速度模型。在射线层析中, 需要在偏移后的共成像点道集拾取剩余时差, 然后将拾取的剩余时差沿着射线路径反传得到速度更新量。射线层析方法计算效率高，当地震波传播过程可用射线理论近似时, 利用射线层析反演能够得到合理的速度模型^[6-8]。但是, 在复杂地质构造条件下, 由于射线多路径的影响及射线理论高频近似的局限, 射线层析方法很难反演出正确的速度模型。

通过引入有限频地震波传播算子, 波动方程旅行时层析可以克服射线层析的多路径问题。在处理运动学频散和采集孔径有限的数据时, 波动方程层析优于射线层析^[7]。波动方程层析中一个关键的步骤是提取正演数据和观测数据之间的旅行时差异。在数据域中, 最直接的旅行时时差估算方法是在波形数据中拾取两个相关的波形, 然后选择观测波形和计算波形互相关最大时所对应的时移量^[9]。在成像域中, 旅行时时差可以在共成像点道集(CIGs) 中拾取同相轴的剩余时差或剩余深度差^[10-11]。上述旅行时反演方法的目标函数是最小化拾取的相对旅行时时差。尽管工业界存在多种自动拾取技术, 但是在多数情况下, 为了降低自动拾取带来的误差, 质量控制和人工干预仍然必不可少。

为了避免拾取带来的误差, VAN LEEUWEN等^[12]采用互相关加权范数作为目标函数, 该目标函数隐式地度量正演数据和观测数据之间的旅行时差, 不需要做任何拾取, 就可以进行全自动反演。但是, 由于观测数据中不同反射波到达时之间容易发生串扰, 因此这个目标函数对反射波数据容易出现多走时问题。一些成像域的速度建模方法, 如波动方程偏移速度分析^[13], 可以通过偏移成像克服这个问题。这类方法通过最小化成像扰动, 使得成像空间中的共成像点道集拉平或者地下偏移距道集聚焦在零偏移距。在波动方程偏移速度分析类方法中, 微分相似度最优化(DSO) 方法将震源波场及接收波场互相关之后的结果乘以地下偏移距作为误差函数^[14]。DSO通过隐式度量地下偏移距道集的聚焦性或者角道集的拉平程度自动更新速度模型。但是, 常规的DSO梯度往往出现明显的近垂直的条带状噪声, 这类噪声会降低反演的收敛速度^[15-18]。通过改进成像剩余量的提取方法能够避免这种噪声, 比如水平收缩算子提取无限小的成像扰动^[18], 但是这类方法的速度更新量已经不是原目标函数的梯度。2009年, ZHANG等^[19]提出了一种基于旅行时的波动方程反演方法, 该方法反传延迟炮记录, 并通过最大化反传波场在震源位置处对应激发时间的能量, 更新速度模型, 它只依赖于地震数据中的走时信息, 而对地震数据中存在的振幅误差不敏感。2011年, LUO等^[20]提出反褶积型波动方程反演, 该方法最小化正演数据和观测数据之间加权反褶积的范数, 全自动, 对振幅误差不敏感, 而且对带限(或非脉冲) 震源函数也不敏感。2014年, WARNER等^[21]提出自适应波形反演(AWI), 该方法首先利用维纳滤波器匹配预测数据和观测数据, 然后通过最小化维纳滤波器非零延迟系数的能量, 更新速度模型。如果维纳滤波器和输入数据长度相同, AWI和反褶积型波动方程反演类似。实际上, 设计一个正演数据和观测数据每个地震道都完全匹配的维纳滤波器十分困难。

本文介绍一种全自动的波动方程旅行时反演方法, 该方法能够利用旅行时信息反演得到运动学意义上准确的速度模型。首先介绍透射波和反射波情形下全走时反演(FTI) 的基本原理, 并分析FTI的梯度和Born近似下的梯度存在的差异, 然后通过理论模型和实际资料测试, 验证FTI方法的正确性和有效性。

1 基本原理 1.1 透射波FTI

首先介绍数据域透射波FTI的基本原理。为方便起见, 本文假定:①地下模型为常密度声介质模型；②假定模型扰动Δs只产生地震波场的时间变化Δτ。基于第②个假设, 地震波同相轴可以完全由旅行时来表征。令P_obs(x_r, t; x_s) 和P_cal(x_r, t; x_s) 分别表示炮点位于x_s, 接收点位于x_r, 时刻为t的观测地震道和计算地震道。定义两道之间的互相关函数:

$c\left( {{x_{\rm{r}}}, \tau ;{x_{\rm{s}}}} \right) = \int {{P_{{\rm{obs}}}}\left( {{x_{\rm{r}}}, t + \tau ;{x_{\rm{s}}}} \right){P_{{\rm{cal}}}}\left( {{x_{\rm{r}}}, t;{x_{\rm{s}}}} \right){\rm{d}}} t$

(1)

式中:τ为观测地震道和计算地震道之间的时移量。目标函数定义为^[12]:

$\varepsilon = \frac{1}{2}\sum\limits_{{x_{{\rm{r }}}}, {x_{\rm{s}}}} {\int {{{\left[{\tau c\left( {{x_{\rm{r}}}, \tau ;{x_{\rm{s}}}} \right)} \right]}^2}{\rm{d}}} \tau } $

(2)

式中:τ∈[-T, T]是线性加权函数, 用于压制互相关c(x_r, τ; x_s) 在非零时移处的能量, 其中, T是估算观测地震道和计算地震道之间旅行时时差的最大值。最小化目标函数((2) 式) 意味着选择速度模型使得互相关能量聚焦在零时移处。与LUO等^[9]提出的互相关波动方程旅行时反演方法相比, 该方法不需要显式拾取互相关时移量, 而是自动计算出来, 且对振幅谱的差异也更加不敏感^[12]。

经过一系列的公式推导^[22], 目标函数((2) 式) 关于模型的Born梯度为:

$\begin{array}{l} \frac{{\partial \varepsilon }}{{\partial s\left( {x'} \right)}} =-2{s_0}\left( {x'} \right)\sum\limits_{{x_{\rm{r}}}, {x_{\rm{s}}}} {\int {{\tau ^2}c\left( {{x_{\rm{r}}}, \tau ;{x_{\rm{s}}}} \right)} } \int {{P_{{\rm{obs}}}}\left( {{x_{\rm{r}}}} \right., } \\ \left. {t + \tau ;{x_{\rm{s}}}} \right)g\left( {x', t;{x_{\rm{r}}}} \right)*{{\ddot P}_{{\rm{cal}}}}\left( {x', t;{x_{\rm{s}}}} \right){\rm{d}}t{\rm{d}}\tau \end{array}$

(3)

这里, g(x′, t; x_r) 为格林函数，表示单位源在地表x_r处激发后在地下x′处所观测的波场。Born近似表明速度扰动和波场扰动之间呈线性关系。波场扰动包含了振幅和相位信息。但是, 互相关目标函数((2) 式) 用于度量相对时移量, 和基于Born近似的梯度公式((3) 式) 的意义不一致, 这会导致常规的Born近似最优化算法无法获得目标函数所对应的正确梯度, 使得反演过程收敛很慢或者不收敛。

目标函数隐式地最小化拾取观测和计算的透射波到达时差。当梯度和目标函数意义一致时, 目标函数的这种误差度量方式才能产生有意义的模型更新量。

假设模型扰动Δs只产生波传播的旅行时变化Δτ。, 对于单个炮点-检波点对, 当慢度模型从s₀变化到s₀+Δs, 计算地震道P_cal(x_r, t; x_s)_s0仅存在时间变化Δτ, 即:

${P_{{\rm{cal}}}}{\left( {{x_{\rm{r}}}, \tau ;{x_{\rm{s}}}} \right)_{{s_0} + \Delta s}} = {P_{{\rm{cal}}}}{\left( {{x_{\rm{r}}}, t + \Delta \tau ;{x_{\rm{s}}}} \right)_{{s_0}}}$

(4)

利用公式(4) 并经过一系列公式推导^[22], 得到目标函数((2) 式) 关于模型的FTI梯度:

$\begin{array}{l} \frac{{\partial \varepsilon }}{{\partial s\left( {x'} \right)}} = 2{s_0}\left( {x'} \right)\sum\limits_{{x_{\rm{r}}}, {x_{\rm{s}}}} {\int {\tau c{{\left( {{x_{\rm{r}}}, \tau ;{x_{\rm{s}}}} \right)}^2}{\rm{d}}\tau } } \\ \frac{{\int {{{\dot P}_{{\rm{cal}}}}\left( {{x_{\rm{r}}}, t;{x_{\rm{s}}}} \right)g\left( {x', t;{x_{\rm{r}}}} \right)*{{\ddot P}_{{\rm{cal}}}}\left( {x', t;{x_{\rm{s}}}} \right){\rm{d}}t} }}{{\int {{{\dot P}_{{\rm{cal}}}}\left( {{x_{\rm{r}}}, t;{x_{\rm{s}}}} \right){{\dot P}_{{\rm{cal}}}}\left( {{x_{\rm{r}}}, t;{x_{\rm{s}}}} \right){\rm{d}}t} }} \end{array}$

(5)

与Born近似下的梯度公式((3) 式) 相比, FTI梯度仅包含旅行时信息, 这和目标函数的误差度量方式一致。注意, 对单个炮点-检波点对, 伴随源(或FTI反传剩余量) 为:

$\delta S = {\dot P_{{\rm{cal}}}}\left( {{x_{\rm{r}}}, t;{x_{\rm{s}}}} \right)\frac{{\int {tc{{\left( {{x_{\rm{r}}}, \tau ;{x_{\rm{s}}}} \right)}^2}\tau } }}{{\int {{{\dot P}_{{\rm{cal}}}}\left( {{x_{\rm{r}}}, t;{x_{\rm{s}}}} \right){{\dot P}_{{\rm{cal}}}}\left( {{x_{\rm{r}}}, t;{x_{\rm{s}}}} \right){\rm{d}}} t}}$

(6)

从(6) 式可以看出, 这个剩余量为合成地震道的一阶导数Ṗ_cal(x_r, t; x_s) 乘以加权平均时移量∫τc(x_r, τ; x_s)²dτ, 再除以照明因子∫Ṗ_cal(x_r, t; x_s)Ṗ_cal(x_r, t; x_s) dt。与射线层析类似, FTI波路径只取决于速度模型, 而与观测数据无关。因此, 该方法对噪声不敏感, 且较为稳定。需要注意的是, FTI梯度包含一个惩罚因子τ, 而Born梯度包含的惩罚因子是τ², 这是FTI能够在迭代过程中计算出正确梯度的关键所在。

基于速度扰动只产生旅行时变化的假设, 可以将公式中的P_cal(x_r, t; x_s) 替换为P_obs(x_r, t+Δτ; x_s), 其中，Δτ为合成地震道匹配观测地震道的时移量。伴随源可以改写为:

$\begin{array}{l} \delta S = \int {\tau c{{\left( {{x_{\rm{r}}}, \tau ;{x_{\rm{s}}}} \right)}^2}{\rm{d}}\tau \cdot } \\ \frac{{{{\dot P}_{{\rm{obs}}}}\left( {{x_{\rm{r}}}, t + \Delta \tau ;{x_{\rm{s}}}} \right)}}{{\int {{{\dot P}_{{\rm{obs}}}}\left( {{x_{\rm{r}}}, t + \Delta \tau ;{x_{\rm{s}}}} \right){{\dot P}_{{\rm{obs}}}}\left( {{x_{\rm{r}}}, t + \Delta \tau ;{x_{\rm{s}}}} \right){\rm{d}}t} }}\\ \approx \int {\tau c{{\left( {{x_{\rm{r}}}, \tau ;{x_{\rm{s}}}} \right)}^2}{\rm{d}}\tau \cdot } \\ \frac{{{{\dot P}_{{\rm{obs}}}}\left( {{x_{\rm{r}}}, t + \tau ;{x_{\rm{s}}}} \right)}}{{\int {{{\dot P}_{{\rm{obs}}}}\left( {{x_{\rm{r}}}, t + \tau ;{x_{\rm{s}}}} \right){{\dot P}_{{\rm{obs}}}}\left( {{x_{\rm{r}}}, t + \tau ;{x_{\rm{s}}}} \right){\rm{d}}t} }} \end{array}$

(7)

实际上, 满足方程(7) 的条件是地震子波是Delta函数, 这样对任意的τ≠Δτ都满足c(x_r, τ; x_s)=0。在此情况下, 积分项中只有τ=Δτ时不为0。和直接拾取时移量Δτ不同, 伴随源是通过最大化互相关函数c(x_r, τ; x_s) 计算出“最佳”的时移量。

1.2 反射波FTI

反射波FTI在成像空间可以由透射波的波路径来表征, 其目标是将扩展成像结果聚焦在反射层深度上。地震波场由平滑的背景速度模型和反射系数正演得到。这种方式预测的地震数据包含所有的一次反射波, 而不含任何多次波。

震源波场P_s(x, t; x_s) 和反向传播得到的接收波场P_r(x, t; x_s) 在成像点x位置处的时移互相关^[23]为:

$I\left( {x, \tau } \right) = \int {\int {{P_{\rm{s}}}\left( {x, t + \tau ;{x_{\rm{s}}}} \right){P_{\rm{r}}}\left( {x, t;{x_{\rm{s}}}} \right){\rm{d}}t{\rm{d}}{x_{\rm{s}}}} } $

(8)

其中, τ(x, x_s)∈[-T, T]是时移量。沿着时移坐标轴方向对偏移成像进行扩展, 可以线性模拟速度扰动产生的时移量。时移量弥补了速度误差并使得成像结果在非零时移处形成反射能量。正确模型下的成像结果聚焦在零时移处。为了自动度量反射能量在不同深度位置的时移量, 定义目标函数为:

$\varepsilon = \frac{1}{2}\sum\limits_x {\int {{{\left[{\tau I\left( {x, \tau } \right)} \right]}^2}} } {\rm{d}}\tau $

(9)

时移量τ用于压制非零时移处的能量。最小化上述目标函数意味着选择速度模型使得扩展成像结果聚焦在零延迟时刻。

对于DSO方法, 速度更新量(梯度) 表示为:

$\frac{{\partial \varepsilon }}{{\partial s\left( {x'} \right)}} = \sum\limits_x {\int {{\tau ^2}I\left( {x, \tau } \right)\frac{{\partial I\left( {x, \tau } \right)}}{{\partial s\left( {x'} \right)}}{\rm{d}}\tau } } $

(10)

公式(10) 的详细推导过程参考文献[14]。在原来的DSO方法中, 扩展成像道集定义在地下偏移距域。为简单起见, 这里仅分析时移成像道集。地下空移量和时移量可以相互转换^[24]。与公式(5) 类似, DSO梯度包含了振幅和相位信息, 在Born近似下无法解耦。但是目标函数度量的是相对时移量, 因此, DSO梯度和目标函数的意义不一致。FEI等^[15]指出反射层截断或畸变会对DSO梯度引入噪声, 降低收敛速度。SHEN等^[18]提出了一种地下偏移距域水平收缩方法, 可以产生无噪声的速度更新量。该方法通过水平校正算子无限小地改进当前成像结果。成像扰动则是校正后的和原来的成像结果之差。尽管由此生成的速度更新量不会产生噪声, 但是它已经不是原来DSO目标函数的梯度。

和透射波情形类似, FTI假定模型扰动只产生时移成像道集聚焦能量的时间变化, 即:

$I{\left( {x, \tau } \right)_{{s_0} + \Delta s}} = I{\left( {x, \tau + \Delta \tau } \right)_{{s_0}}}$

(11)

和数据域推导类似, 目标函数(9) 式关于慢度的梯度公式^[22]为:

$\begin{array}{l} \frac{{\partial \varepsilon }}{{\partial s\left( {x'} \right)}} = \frac{{2{s_0}\left( {x'} \right)}}{E}\sum\limits_x {\int {\tau I{{\left( {x, \tau } \right)}^2}} } \left\{ {\int {\left[{{{\ddot P}_{\rm{s}}}\left( {x', t;{x_{\rm{s}}}} \right)*} \right.} } \right.\\ \left. {g\left( {x, t + \Delta \tau ;x'} \right)} \right]{{\dot P}_{\rm{r}}}\left( {x, t;{x_{\rm{s}}}} \right) + {{\dot P}_{\rm{s}}}\left( {x, t + \Delta \tau ;{x_{\rm{s}}}} \right) \cdot \\ \left[{{{\ddot P}_{\rm{r}}}\left( {x', t;{x_{\rm{s}}}} \right)*g\left( {x, -t;x'} \right){\rm{d}}t} \right]{\rm{d}}\tau \end{array}$

(12)

其中, 分母E为:

$E = \iint{\ddot P_{\rm{s}}}\left( {x, t + \Delta \tau ;{x_{\rm{s}}}} \right){P_{\rm{r}}}\left( {x, t;{x_{\rm{s}}}} \right){\rm{d}}t{\rm{d}}{x_{\rm{s}}}$

(13)

旅行时剩余量Δτ度量的是聚焦性, 可以在时移成像道集中拾取得到。由于成像道集I(x, τ) 只聚焦在Δτ位置, 因此公式(12) 在任何τ≠Δτ处生成的波路径可以忽略, 即:

$\begin{array}{l} \frac{{\partial \varepsilon }}{{\partial s\left( {x'} \right)}} = \frac{{2{s_0}\left( {x'} \right)}}{E}\sum\limits_x {\int {\tau I{{\left( {x, \tau } \right)}^2}} } \left\{ {\int {\left[{{{\ddot P}_{\rm{s}}}\left( {x', t;{x_{\rm{s}}}} \right)*} \right.} } \right.\\ g\left( {x, t + \tau ;x'} \right){{\dot P}_{\rm{r}}}\left( {x, t;{x_{\rm{s}}}} \right) + {{\dot P}_{\rm{s}}}\left( {x, t + \tau ;{x_{\rm{s}}}} \right) \cdot \\ \left. {\left[{{{\ddot P}_{\rm{r}}}\left( {x', t;{x_{\rm{s}}}} \right)*g\left( {x, -t;x'} \right)} \right]{\rm{d}}t} \right\}{\rm{d}}\tau \end{array}$

(14)

(14) 式不需要拾取Δτ, 而是采用积分自动提取Δτ对应的主要能量。

2 数值试验 2.1 Wedge模型

首先采用二维理论模型对本文方法进行测试。图 1a是真实速度模型。采用线性Born正演方法合成反射地震数据。震源是主频为10 Hz的雷克子波, 共正演126炮, 炮间距为80 m, 每炮以固定接收排列方式接收数据, 共251道, 道间距为40 m。初始速度模型为常速(2 000 m/s) 模型。图 1b为FTI反演迭代10次后的速度模型。对比图 1a和图 1b可以看出, 反演结果较好地恢复了真实速度模型的变化趋势。图 2各个子图分别显示初始速度模型、反演速度模型和真实速度模型下的叠前深度偏移成像结果和地下偏移距道集。对比不同速度模型下的成像结果和道集可以看出, 采用反演后的速度模型, 反射层基本偏移到真实的空间位置, 地下偏移距道集聚焦在零偏移距处。理论模型测试结果表明:FTI从常速度模型开始反演, 得到的背景速度模型, 可以满足构造成像的精度要求。

2.2 实际资料

采用某二维实际地震数据对本文方法进行测试。该实际资料共801炮, 炮间距25 m, 每炮以固定接收排列方式接收数据, 共1 026道, 道间距25 m。采用常速模型(4 500 m/s) 作为初始速度模型进行反演。图 3是FTI反演迭代10次后的速度模型。图 4a和图 4b分别为初始速度模型下的叠前深度偏移成像剖面和地下偏移距道集。图 5a和图 5b分别为FTI反演速度模型下的叠前深度偏移成像剖面和地下偏移距道集。对比图 4和图 5可以看出, 采用FTI反演得到的速度模型, 可以改善成像质量, 提高偏移成像道集的聚焦程度。

3 结论

FTI是一种利用旅行时信息构建运动学意义上正确速度模型的反演方法，其核心思想是减少振幅对旅行时反演的影响。FTI可以在数据域也可以在成像域实现, 分别用于接收端的透射波和反射层上的反射波旅行时反演。数据域全走时反演的目标函数是观测和计算的透射波数据之间互相关的线性加权范数。成像域FTI的目标函数是时间或者空间延迟扩展成像道集的线性加权范数。两个目标函数都可以隐式地度量旅行时相对误差。基于速度扰动只产生波场的旅行时变化这一假设, 本文分别介绍了透射波和反射波情形下的FTI梯度(波路径)。与Born近似意义下的梯度相比, FTI梯度只包含旅行时信息, 和目标函数极小化旅行时时差的意义一致。FTI不需要准确的初始速度模型和地震数据中的低频信息。理论模型和实际资料测试结果表明, FTI具有很好的收敛性, 采用十分简单的初始速度模型, 在很少的迭代次数内也能收敛到满足构造成像精度的要求。