随着油气勘探开发复杂程度的不断加深, 如何提高储层预测的精度逐渐成为专家学者研究的重点, 而寻找一种适用于复杂介质的反演方法尤其重要。全波形反演(FWI) 利用叠前地震波场中的运动学及动力学信息反演得到地下岩石物理特性, 能够解释复杂地质背景下的构造细节信息^[1-2], 自1984年^[3]问世后得到了飞速发展。PRATT等将全波形反演的思想引入频率域, 推导了频率域内最速下降法全波形反演梯度公式^[4-5]。SHIN等提出Laplace-Fourier域全波形反演方法, 通过引入复频率, 在缺少低频的情况下反演出较好的结果^[6]。此外, SIRGUE等将全波形反演应用于混合域, 得到了较好的反演结果^[7]。在这些方法中, 时间域全波形反演可对数据进行灵活预处理, 适用于多节点并行运算, 已被成功地应用于实际地震数据处理^[8]。

在单参数全波形反演中, 常假设密度随空间位置变化缓慢, 从而将密度设为一已知恒定值, 反演得到速度结果。但在实际地层中, 密度是指示储层流体的重要信息, 在油气解释中起着重要作用。密度的变化会引起反射率等变化, 影响地震波传播及反演结果^[9], 如何从地震数据中得到有效的密度信息越来越成为精细化地震勘探的研究重点。多参数全波形反演是一种强非线性问题, 且速度与密度之间存在串扰现象^[10], 严重影响反演结果的精度, 因此, 选择合适的反演方法和策略是重建高精度参数模型的关键。FORGUES等通过辐射模式分析指出, 很难从全波形反演中得到密度^[11]。JEONG等提出频率域密度反演策略, 先将密度设置为常数, 反演得到弹性参数, 再进行速度、密度联合反演, 提高了反演精度^[12]。PRIEUX等首先采用长偏移距数据反演速度, 然后利用全偏移数据反演速度和密度, 提高了密度反演的稳定性, 但是密度的反演结果仍然不是很理想^[13]。董良国等通过分析目标函数随密度和速度扰动的变化规律, 指出全波形反演对密度的较小尺度摄动具有较好敏感度, 但是无法反演中长波长的密度摄动^[14]。RAKNES等讨论了弹性波全波形反演过程中反演参数的组合方式, 并提出多参数全波形反演策略--逐一反演各个参数^[15]。杨积忠等提出分步变密度声波方程多参数反演策略, 得到了比较好的速度和密度反演结果^[16]。本文将该方法运用于时间域声波全波形反演中，利用照明预处理L-BFGS方法使不同深度的梯度能量分布更为均衡，进一步提高模型参数收敛速度和反演精度。

1 照明预处理L-BFGS

全波形反演是Bayes估计理论在勘探地球物理领域的一个应用范例^[17], 它将波动方程正演模拟得到的地震记录与实际观测地震记录进行匹配, 通过使数据之差最小建立目标函数, 寻找最佳模型参数。时间域全波形反演算法在实现过程中主要包括两个过程:以地震子波为震源的正演过程和以地震记录残差为震源的波场反传过程^[18-19]。其目标函数可以写为:

$E = \frac{1}{2}\delta {\boldsymbol{d}^{\rm{T}}}\delta \boldsymbol{d}$

(1)

式中:E为误差泛函; δd为波场残差。模型参数m可以通过下式进行迭代更新:

${\boldsymbol{m}_{k + 1}} = {\boldsymbol{m}_k}-{\alpha _k}{\boldsymbol{h}_k}$

(2)

式中:α_k是迭代步长; k是迭代次数; h_k是迭代方向。全波形反演的目标函数高度非线性, 在迭代过程中存在强烈的非线性问题, 因此要寻找合适的优化方法提高迭代的稳定性。

牛顿类优化算法中通过对二阶偏导海森算子的计算, 均衡梯度算子修正量, 可以加速迭代收敛效率。但是海森矩阵占用大量的存储空间和计算时间, 对计算机的存储能力和计算能力都提出了很高的要求, 导致大规模全波形反演无法开展。1989年, LIU等提出的L-BFGS方法^[20]是对牛顿类优化方法的改进。该方法对海森矩阵进行近似代替, 不直接保存海森矩阵, 而是在迭代过程中保存前n次迭代的参数修正值及梯度信息, 大大减小了存储量及计算量^{[18, 20-22]}:

$\begin{array}{l} {\mathit{\boldsymbol{H}}_{k + 1}} = \left( {\mathit{\boldsymbol{V}}_k^{\rm{T}} \cdots \mathit{\boldsymbol{V}}_{k-n}^{\rm{T}}} \right)\mathit{\boldsymbol{H}}_k^0\left( {{\mathit{\boldsymbol{V}}_{k-n}} \cdots {\mathit{\boldsymbol{V}}_k}} \right) + \sum\limits_{j = 0}^n {{\mathit{\boldsymbol{\rho }}_{k-n + j}}} \cdot \\ \left( {\prod\limits_{l = 0}^{n - j - 1} {\mathit{\boldsymbol{V}}_{k - l}^{\rm{T}}} } \right) \cdot {\mathit{\boldsymbol{s}}_{k - n + j}}\mathit{\boldsymbol{s}}_{k - n + j}^{\rm{T}}\left( {\prod\limits_{l = 0}^{n - j - 1} {{\mathit{\boldsymbol{V}}_{k - l}}} } \right) \end{array}$

(3)

式中: H_k⁰是简单正定矩阵, s_k=m_k+1-m_k, y_k=∇E(m_k+1), V_k=(I -ρ_ky_ks_k^T), ρ_k=1/(s_k^Ty_k)。

利用近似海森矩阵对梯度的预处理可以校正波场几何扩散对不同深度梯度的影响, 增加参数在模型深部的校正量, 同时降低浅层不合理的低频修正量, 提高建模、成像的分辨率。但是, 在L-BFGS法前n次迭代过程中, 不能计算得到近似海森矩阵, 梯度的能量不平衡, 影响反演结果。全波形反演通过正传波场和反传波场在时间上的零延迟互相关来计算模型参数的更新量, 因此梯度能量的不均衡很大程度上是由地震波传播过程中几何扩散的影响引起的。而照明分析可以研究地震波在地下介质中传播的能量分布情况, 进行能量补偿^[23-24], 同时, 全波形反演需要大量的计算量和存储空间, 为了在得到较好反演效果的同时尽量减少全波形反演的计算负担, 我们采用单向照明强度近似代替双向照明强度的方法, 通过正演模拟得到模型各个点的波场, 将其作为该点单向照明强度。但检波器接收到的地震波是双程旅行的过程^[23], 为了便于计算, 我们将模型中各个点波场值的平方作为双向照明强度, 对前n次L-BFGS法计算的梯度值进行预处理。则N炮照明总强度I_S可用下式表示:

${I_{\rm{S}}}\left( {x, z} \right) = \sum\limits_t {\sum\limits_{i = 1}^N {{u_i}\left( {x, z} \right) \cdot {u_i}\left( {x, z} \right)} } $

(4)

式中:u_i(x, z) 为t时刻在模型坐标(x, z) 处的波场值。因此, 参数m的前n次迭代可以用下式来表示:

${\boldsymbol{m}_{k + 1}} = {\boldsymbol{m}_k}-{\alpha _k}\nabla E\left( \boldsymbol{m} \right)I_{\rm{S}}^k$

(5)

2 多参数全波形反演

时间域二阶变密度声波方程如下:

$\frac{1}{{\kappa \left( {x, z} \right)}}\frac{{{\partial ^2}p}}{{\partial {t^2}}}- \nabla \cdot \left[{\frac{1}{{\rho \left( {x, z} \right)}}\nabla p} \right] = \delta \left( {x -{x_{\rm{s}}}} \right)s$

(6)

式中:κ表示体积模量; ρ表示密度; p表示压力场; x_s为震源坐标; s为震源函数。则速度和密度的梯度可分别用下式表示:

$\frac{{\partial E}}{{\partial v}} = \frac{2}{{\rho {v^3}}}{p_{\rm{f}}}{p_{\rm{b}}}$

(7)

$\frac{{\partial E}}{{\partial \rho }} = \frac{2}{{{\rho ^2}{v^2}}}{p_{\rm{f}}}{p_{\rm{b}}} + \frac{1}{{{\rho ^2}}}\nabla {p_{\rm{f}}} \cdot \nabla {p_{\rm{b}}}$

(8)

式中:p_f是正传波场; p_b是反传波场; v表示速度。

多参数全波形反演因速度和密度之间的串扰影响非常不稳定, 难以获得精确的反演结果^[10]。速度的变化会引起地震波传播的旅行时和振幅发生扰动, 而密度的变化会影响地震波的振幅大小, 且随着角度的变大, 波场变化逐渐减小^[9]。同时, 地震反射波对密度的较小尺度摄动非常敏感, 而中长波长的密度摄动对反射波基本没有影响^[16]。因此, 多参数全波形反演虽然可以重建较好的速度模型, 但只能重建密度的高频成分^{[9, 14, 16]}。与此同时, 全波形反演, 尤其是多参数全波形反演, 对初始模型的依赖性较强, 初始速度模型越准确, 密度反演结果越精确。

基于以上分析, 我们将杨积忠等提出的分步反演策略^[16]用于多参数时间域全波形反演, 得到速度和密度反演结果。第一步是利用多参数全波形反演得到较为精确的速度模型, 此时密度反演结果较差; 第二步是将第一步反演的速度结果作为初始速度模型, 结合初始密度模型进行速度、密度联合反演, 此时由于初始速度模型精度的提高, 密度反演的精度也有很大提高。为了进一步提高反演精度, 可将第二步得到的速度反演结果作为初始模型, 利用初始密度模型再进行一轮速度、密度联合反演。

3 模型测试 3.1 照明预处理L-BFGS方法测试

为了测试照明预处理L-BFGS方法的有效性, 我们将其应用于Marmousi模型，进行速度参数的全波形反演, 并与最速下降法和L-BFGS法进行对比。实际速度模型如图 1a所示, 网格点数为260×130, 网格间距10 m。地表放置38个震源, 震源间距70 m, 并放置260个检波器, 检波器间距10 m。所用震源为主频30 Hz的雷克子波。图 1b为初始速度模型。

图 2是最速下降法、L-BFGS法、照明预处理L-BFGS法速度反演结果, 可以看到:最速下降法反演效果不理想, 虽然浅层速度反演得很好, 但是深层速度值偏差较大; L-BFGS方法在反演出较好的浅层速度的基础上, 由于自身可以校正波场几何扩散的影响, 深部速度的校正量增加, 速度反演更为准确; 而经过照明预处理后的L-BFGS方法进一步弥补了波场传播时深层能量的缺失, 得到了更好的深部速度反演结果。图 3是三种方法在不同水平位置处的反演速度随深度变化的曲线, 可以看到照明预处理L-BFGS法得到的速度结果与真实值拟合最好, 反演精度最高。图 4是三种反演方法的目标函数收敛曲线, 可以看到照明预处理L-BFGS法收敛最快, 失配函数最小, 说明照明预处理L-BFGS方法可以提高收敛速度及精度, 减小运算量及运算时间。

3.2 分步多参数全波形反演测试

比较最速下降法、传统L-BFGS法及照明预处理L-BFGS法的效果, 我们选择反演精度最高的照明预处理L-BFGS方法, 采用与单参数反演模型相同的参数, 进行速度、密度的联合反演测试。图 5a是实际速度模型, 图 5b是用Gardner公式和实际速度模型计算得到的实际密度模型; 图 6a和图 6b是初始速度模型和初始密度模型, 图 7a和图 7b是速度和密度第一步反演结果。从图 7可以看出, 速度和密度的反演精度都不高, 浅层速度反演较好, 构造也得到很好的重建, 但深部速度和实际值仍有很大偏差, 密度的反演结果精度不高, 不能满足反演精度的要求。

为了得到更好的反演结果, 我们舍弃密度反演结果, 将第一步的速度反演结果(图 7a) 作为初始速度模型, 联合初始密度模型再进行一次多参数同时反演(第二步反演), 结果如图 8所示。可以看到, 速度反演结果精度有了明显提高, 中深层速度值与真实值吻合较好, 密度反演精度较之前也有了提高, 证明分步多参数全波形反演在提高反演精度方面是有效的。但是, 受到初始速度模型精度的制约, 密度精度的提高有限(图 8b)。要进一步重建更为精确的密度模型, 需要更为准确的初始速度模型, 因此我们将第二步的速度反演结果作为初始速度模型, 结合初始密度模型, 又进行了一次多参数联合反演(最终反演), 结果如图 9所示。由图 9可以看到，速度和密度反演的精度都有了明显改善。比较分步多参数全波形反演过程中不同水平位置处速度、密度值纵向剖面图(图 10), 可以看到速度、密度联合反演的精度每一次都在提高, 最终的模型反演结果与实际模型值拟合较好。

4 结论

在多参数全波形反演中, 由于速度和密度的耦合关系影响了反演的精度和效率, 密度反演结果很不理想, 选择合适的反演算法和反演策略对全波形反演至关重要。将分步多参数全波形反演策略应用于时间域多参数声波全波形反演中, 通过建立更为准确的速度初始模型, 减小速度和密度间的串扰影响, 得到了更为准确的密度反演结果。L-BFGS是一种平衡梯度能量十分有效的方法, 但是得到的反演结果分辨率仍有上升空间, 通过对L-BFGS法前n次的计算结果进行照明预处理, 可以在不增加过多计算量的基础上更有效地平衡模型修正量的能量, 提高反演精度和收敛效率, 减少反演计算量和计算时间。