0 引言

基于AVO理论的叠前地震反演技术利用角度道集数据，以地质、测井资料为约束，挖掘地下介质的纵波、横波速度及密度等多种弹性参数信息，进而为流体识别^[1]、岩性分类^[2]、储层预测^[3]等提供数据基础，是地震数据定量分析最重要的手段之一。

早期的叠前反演技术采用Zoeppritz方程近似式(如Aki-Richards公式等)构建线性正演算子，其反问题可通过线性反演方法解析计算^[4-5]，然而近似式存在近入射角及弱阻抗差等假设条件，势必影响反演结果的准确性。为获取更可靠的弹性参数结果，后期的叠前反演技术通常采用精确Zoeppritz方程^[6-8]，但基于精确方程的地震正演算子非线性程度高，其反问题难以解析表达。此外，叠前反演是欠定问题，需依赖先验约束以稳定反演过程。传统Tikhonov正则化易造成反演结果过度平滑^[9]，地震反演适合于采用边界保护类正则化方法^[10-13]，但此类非二次约束项会造成目标函数偏导矩阵求解困难，致使常规梯度类优化算法难以适用。因此，全局优化算法是求解此类叠前反演问题的重要途径^[14-16]。

模拟退火算法(Simulated Anneding, SA)是一种蒙特卡洛类启发式算法，作为最主要的全局优化算法之一，它具有不依赖初始模型且易跳出局部极值的优势^[17]。Rothman^[18]最先将模拟退火算法引入地球物理反演领域，以解决地震资料的剩余静校正问题。此后，模拟退火算法凭借其可靠性迅速拓展至地震反演领域，如：李景叶等^[19]使用快速模拟退火算法进行叠后时移地震数据的反演；Misra等^[20]使用非常快速模拟退火算法，并结合保边算子获取了具有块状特征的弹性参数模型；周琦等^[21]采用模拟退火算法进行横波速度预测及地震岩石物理建模。然而，模拟退火算法涉及多种优化参数(如初始温度、模型扰动范围等)，这些参数设置的准确性对反演结果具有重要影响。目前，模拟退火参数优化主要采用试算方法，通过合成数据或井旁地震道进行参数调试^[22-23]，但此类方法经验性强，易造成算法不稳定。在优化参数的定量估算方面，Basu等^[24]选用较高的初始温度及较小的衰减系数以保证算法收敛成熟，但这增加了迭代次数，降低了计算效率；姚姚^[25]发展了基于局部势能确定终止温度的策略；Aarts等^[26]提出依据首次接收概率确定初始温度的思路。然而，上述方法的结果在具体实现中仍依赖初始模型、扰动范围等基本参数。

鉴于叠前地震反演技术面向的研究区地质条件日益复杂，降低人为因素干扰、提高算法的稳定性及适用性尤为重要。为此，本文提出了一种基于自适应优化参数设置模拟退火的叠前地震联合反演方法。首先，联合贝叶斯线性反演与模拟退火非线性反演方法，通过线性反演结果驱动后续优化参数的估算及先验模型的构建；其次，针对不同地震道的数据差异，逐道计算适用的初始温度及扰动范围，有效提升了算法的适用性及稳定性；最后，通过合成数据测试及实测资料应用证明该方法提高了弹性参数反演结果的准确度。

1 方法原理 1.1 叠前地震反演技术

作为叠前地震反演技术的理论基础，Zoeppritz方程描述了地层分界面处反射及透射系数随入射角的变化规律，其表达式为^[27]

$ \left[\begin{array}{cccc} \sin \theta_1 & \cos \varphi_1 & -\sin \theta_2 & \cos \varphi_2 \\ -\cos \theta_1 & \sin \varphi_1 & -\cos \theta_2 & -\sin \varphi_2 \\ \sin 2 \theta_1 & \frac{V_{\mathrm{P} 1}}{V_{\mathrm{S} 1}} \cos 2 \varphi_1 & \frac{\rho_2 V_{\mathrm{S} 2}^2 V_{\mathrm{P} 1}}{\rho_1 V_{\mathrm{S} 1}^2 V_{\mathrm{P} 2}} \sin 2 \theta_2 & -\frac{\rho_2 V_{\mathrm{S} 2} V_{\mathrm{P} 1}}{\rho_1 V_{\mathrm{S} 1}^2} \cos 2 \varphi_2 \\ \cos 2 \varphi_1 & -\frac{V_{\mathrm{S} 1}}{V_{\mathrm{P} 1}} \sin 2 \varphi_1 & -\frac{\rho_2 V_{\mathrm{P} 2}}{\rho_1 V_{\mathrm{P} 1}} \cos 2 \varphi_2 & -\frac{\rho_2 V_{\mathrm{S} 2}}{\rho_1 V_{\mathrm{P} 1}} \sin 2 \varphi_2 \end{array}\right] \left[\begin{array}{c} R_{\mathrm{PP}} \\ R_{\mathrm{PS}} \\ T_{\mathrm{PP}} \\ T_{\mathrm{PS}} \end{array}\right]=\left[\begin{array}{c} -\sin \theta_1 \\ -\cos \theta_1 \\ \sin 2 \theta_1 \\ -\cos 2 \varphi_1 \end{array}\right] $

(1)

式中：V_P、V_S和ρ分别表示纵波、横波速度和密度，下标1、2对应界面上、下地层；θ₁和θ₂分别表示纵波的入射角和透射角；φ₁和φ₂分别表示横波的反射角和透射角。基于褶积模型，PP波地震道集可表示为

$ \boldsymbol{d}=G(\boldsymbol{m})=\boldsymbol{R}_{\mathrm{PP}}(\boldsymbol{m}, \theta) * \boldsymbol{w}+\boldsymbol{e} $

(2)

式中：G表示地震正演算子；m表示地层弹性参数；w表示角度子波；e为随机误差；符号*表示卷积运算。

为有效引入先验约束并描述反演结果的不确定性，将地震反演问题置于贝叶斯框架下，则弹性参数的后验概率P(m|d)可表示为^[4]

$ P(\boldsymbol{m} \mid \boldsymbol{d}) \propto P(\boldsymbol{d} \mid \boldsymbol{m}) \times P(\boldsymbol{m}) $

(3)

式中：P(d|m)为似然函数，用以描述正演模型的不确定性；P(m)为先验概率分布，用以引入先验约束。地震反演是欠定问题，需借助先验约束以稳定反演结果，借鉴边界保护正则化思想^{[13, 22]}，添加两种先验约束项，即P(m)=P₁(m)×P₂(m)，其分别表示为

$ P_1(\boldsymbol{m})=\frac{1}{Z} \exp \left\{-\sum\limits_{\mathit{\Omega}} \mathit{\Phi}\left[D_c(\boldsymbol{m})\right]\right\} $

(4)

$ \begin{aligned} P_2(\boldsymbol{m})= & \frac{1}{\left(2 \rm{ \mathsf{ π} }\left|\mathit{\boldsymbol{ \boldsymbol{\varSigma}}}_{\boldsymbol{m}}\right|\right)^{\frac{L}{2}}} \times \\ & \exp \left[-\left(\boldsymbol{m}-\boldsymbol{\mu}_{\boldsymbol{m}}\right)^{\mathrm{T}} \mathit{\boldsymbol{ \boldsymbol{\varSigma}}}_\boldsymbol{m}^{-1}\left(\boldsymbol{m}-\boldsymbol{\mu}_{\boldsymbol{m}}\right)\right] \end{aligned} $

(5)

式中：Φ为势函数；D_c为c阶差分算子；Ω为邻域系统；Z为归一化常数；μ_m和Σ_m分别为弹性参数的期望与协方差；L为弹性参数的维度。本文采用一阶邻域(即四近邻)及一阶差分算子(c=1)，并采用具有边界保护特性的势函数^[22]

$ \mathit{\Phi }(\boldsymbol{m})=\frac{\boldsymbol{m}^2}{1+\boldsymbol{m}^2} $

(6)

在反演过程中，P₁(m)通过描述弹性参数的邻域梯度以保护模型的边界特征，P₂(m)通过期望均值及多参数间统计关系以施加低频约束。结合式(3)~式(5)，将反演问题等价于最大后验概率估算，则目标函数可表示为

$ \begin{aligned} O(\boldsymbol{m})= & {[\boldsymbol{d}-G(\boldsymbol{m})]^{\mathrm{T}} \mathit{\boldsymbol{ \boldsymbol{\varSigma}}}_\boldsymbol{e}^{-1}[\boldsymbol{d}-G(\boldsymbol{m})]+} \\ & \eta_1 \sum\limits_{\mathit{\Omega}} \mathit{\Phi }\left[D_c(\boldsymbol{m})\right]+ \\ & \eta_2\left(\boldsymbol{m}-\boldsymbol{\mu}_\boldsymbol{m}\right)^{\mathrm{T}} \mathit{\boldsymbol{ \boldsymbol{\varSigma}}}_\boldsymbol{m}^{-1}\left[\boldsymbol{m}-\boldsymbol{\mu}_{\boldsymbol{m}}\right) \end{aligned} $

(7)

式中：η₁和η₂均为正则化参数，用以平衡两种约束项的关系；Σ_e表示地震数据随机误差的协方差。

1.2 模拟退火算法

模拟退火算法是根据固体物质退火过程与组合、优化问题间的相似性，基于Metropolis概率接收准则，实现全局寻优^[17]。使用模拟退火算法求解地震反演问题具有以下优势：①不依赖初始模型且易跳出局部极值；②避开目标函数偏导矩阵计算，适用于高度非线性反问题。鉴于式(7)非线性强，偏导矩阵求解困难，常规梯度类反演算法难以适用，本文采用模拟退火算法对目标函数全局寻优。

模拟退火算法采用随机扰动以更新模型参数，采用非常快速模拟退火算法^[28-29]的似Cauchy分布扰动公式

$ \begin{aligned} \boldsymbol{m}^{(k)}= & \boldsymbol{m}^{(k-1)}+t^{(k-1)} \operatorname{sign}(\xi-0.5) \times \\ & \left\{\left[1+\frac{1}{t^{(k-1)}}\right]^{|2 \xi-1|}-1\right\} \Delta \boldsymbol{x} \end{aligned} $

(8)

式中：t为温度；sign表示符号函数；ξ属于[0, 1]的随机数；Δx=[ΔV_P, ΔV_S, Δρ]为模型扰动范围；k表示迭代次数。根据Metropolis准则判别是否接收更新后的模型参数

$ P\left[\boldsymbol{m}^{(k)} \rightarrow \boldsymbol{m}^{(k-1)}\right]= \begin{cases}1 & \Delta O \leqslant 0 \\ \exp \left[-\frac{\Delta O}{t^{(k-1)}}\right] & \Delta O>0\end{cases} $

(9)

式中ΔO为更新前、后目标函数的变化量。可见更新后若目标函数值降低，该模型参数无条件被接收；倘若目标函数值增加，该模型参数(恶化解)仍有一定概率被接收，从而赋予算法跳出局部极值的特性。温度的降低按照冷却进度表执行

$ t^{(k)}=t^{(0)} \exp \left(-\beta k^{\frac{1}{L}}\right) $

(10)

式中：t⁽⁰⁾为初始温度；β为衰减因子。

虽然模拟退火算法对求解非线性反问题具有优势，但涉及多种优化参数的预先设置，其中初始温度t⁽⁰⁾及扰动范围Δx作为重要的优化参数，对收敛效果及反演结果具有重要影响。

1.3 自适应优化参数

关于地震反演中模拟退火优化参数的设置，常规试算经验性方法易产生误差，而现有定量设置方法仍依赖初始参数；此外，二维地震数据反演通常采用固定的优化参数，然而地震道之间存在数据差异，因此采用不同的优化参数更符合实际情况。

对此，本文提出一种基于自适应优化参数模拟退火的叠前地震联合反演方法，流程(图 1)主要包括：①采用贝叶斯线性反演获取弹性参数模型(线性反演结果)；②利用线性反演结果逐道估算扰动范围及初始温度；③设置线性反演结果为先验约束项的期望均值；④逐道更新弹性参数并由终止条件确定总迭代次数；⑤输出最终弹性参数反演结果。具体实施过程中，线性反演采用Aki-Richard近似式作为地震正演算子，地震正演模型可简化为

$ \boldsymbol{d}=\boldsymbol{F m}^*+\boldsymbol{e}=\boldsymbol{W S D m}^*+\boldsymbol{e} $

(11)

式中：F为正演算子矩阵；S为含Aki-Richard公式系数的分块矩阵；W为子波矩阵；D为一阶差分矩阵；m^*表示弹性参数的自然对数形式。采用贝叶斯线性反演解析计算弹性参数反演结果^[4]

$ \begin{aligned} \boldsymbol{\mu}_{\boldsymbol{m} \mid \boldsymbol{d}}= & \boldsymbol{\mu}_{\boldsymbol{m}}+\mathit{\boldsymbol{ \boldsymbol{\varSigma}}}_{\boldsymbol{m}} \boldsymbol{F}^{\mathrm{T}}\left(\boldsymbol{F} \mathit{\boldsymbol{ \boldsymbol{\varSigma}}}_{\boldsymbol{m}} \boldsymbol{F}^{\mathrm{T}}+\mathit{\boldsymbol{ \boldsymbol{\varSigma}}}_{\boldsymbol{e}}\right)^{-1} \times \\ & \left(\boldsymbol{d}-\boldsymbol{F} \boldsymbol{\mu}_{\boldsymbol{m}}\right) \end{aligned} $

(12)

式中μ_m|d为弹性参数线性反演结果(后验均值)，此处结果为对数形式，需转换后得到弹性参数反演结果。

将各地震道的采样点范围作为局部势能估算区域，以线性反演结果作为扰动前模型参数，根据扰动前、后目标函数值变化量均值，逐道计算模拟退火初始温度

$ t_i^{(0)}=-\frac{\frac{1}{N} \sum\limits_j^N\left[O_{i j}\left(\boldsymbol{m}_i\right)-O_i\left(\boldsymbol{\mu}_{\boldsymbol{m} \mid \boldsymbol{d}}\right)\right]}{\ln \left(P_{\text {init }}\right)} $

(13)

式中：P_init为初始接收概率，设置P_init=0.9；O_ij(m_i)表示第i道模型参数在第j次扰动时的目标函数值(i=1, 2, ..., I; j=1, 2, ..., N)，其中I和N分别表示总地震道数和最大扰动次数；O_i(μ_m|d)表示以线性反演结果作为模型参数时第i道的目标函数值。由于温度t与扰动范围Δx均未知，因此无法直接采用式(8)进行模型扰动以估算初始温度，对此采用如下简化扰动公式

$ \boldsymbol{m}_i=\boldsymbol{m}_i^{(0)}+\operatorname{sign}(\xi-0.5) \Delta \boldsymbol{x}_i $

(14)

式中：m⁽⁰⁾表示弹性参数初始模型；Δx_i表示第i道的扰动范围。模型参数扰动范围与搜索区间相关，设置扰动范围处于模型参数最大差值范围内，以保证扰动过程的稳定性，利用线性反演结果进行估算

$ \Delta \boldsymbol{x}_i=\frac{1}{2}\left[\max \left(\boldsymbol{\mu}_{\boldsymbol{m} \mid \boldsymbol{d}}\right)-\min \left(\boldsymbol{\mu}_{\boldsymbol{m} \mid \boldsymbol{d}}\right)\right] $

(15)

式中max(μ_m|d)和min(μ_m|d)分别表示第i道线性反演结果的最大值和最小值。

此外，为提高算法的计算效率，设置连续20次更新的解都不被接收或达到最大迭代次数为终止条件，则各地震道具有不同的迭代次数(终止温度)。

2 合成数据测试

采用部分Marmousi 2模型(图 2)进行合成数据测试。该模型共有126道，每道有121采样点(采样间隔为2 ms)。根据褶积模型，采用Ricker子波(主频为50 Hz)合成地震道集数据，角度范围为5°~40°(图 3a)。由图 3b可见，初始模型仅能反映纵波速度的低频背景趋势。

为说明优化参数对反演结果的影响，测试了不同初始温度及扰动范围。图 4展示了第80道地震数据的误差(目标函数值)随迭代次数的收敛情况。由图 4a可见，初始温度对收敛效果的影响不可忽略。若给定相同的迭代次数，较高的初始温度(灰色曲线)难以收敛至较小的误差值，较低的初始温度(蓝色曲线)收敛不充分，而适中的初始温度(红色曲线)的收敛效果最佳。为便于说明扰动范围对结果的影响，图 4b仅测试了不同纵波速度扰动范围下的收敛情况，每道各采样点使用相同的扰动范围。由图 4b可见，扰动范围对收敛效果有显著影响。若选择的扰动范围偏小(蓝色曲线)，则易造成收敛不充分。

图 5展示了在上述优化参数设置下弹性参数的反演结果。由图可见，采用适中初始温度的反演结果(红色曲线)与真实模型吻合度最高(图 5a)；采用较大扰动范围的反演结果(灰色曲线)与真实模型吻合度较高(图 5b)，但在部分区域异常值偏大，需要可靠的先验约束。因此，在后续反演测试及应用中，除利用线性反演结果估算初始温度及扰动范围外，还将其作为目标函数约束项的期望均值以稳定反演过程。

对Marmousi 2二维剖面进行反演，利用贝叶斯线性反演方法获取弹性参数线性反演结果(图 6)。由于线性反演采用Aki-Richards近似式，对大入射角模拟的精度较低并受限于弱阻抗差假设，且在反演过程中未施加边界保护约束，其反演结果(图 6a)对构造细节和地层边界的识别能力较差，但相较于初始模型(图 3b)，该结果能较好地反映模型的中、低频趋势，因而可用于定量估算初始温度及扰动范围，并可作为目标函数约束项的期望均值以稳定后续非线性反演过程。

图 6b和图 6c是分别采用常规模拟退火算法和本文方法的弹性参数反演结果，其中常规算法采用固定优化参数(t⁽⁰⁾=0.5，ΔV_P=0.05、0.03、0.02 km/s)，衰减因子β与初始温度无关，均设为0.95。相比线性反演结果(图 6a)，采用自适应优化参数设置的反演结果(图 6c)边界保护效果更好，对地层及构造边界有更好的恢复能力(见图 6c左箭头指示处)；相比常规模拟退火算法结果(图 6b)，本文方法反演结果层间异常值明显减少(见图 6c中、图 6c右箭头指示处)，反演结果更为稳定。分别选取第90道及200 ms深度反演结果进行对比(图 7)，并列出第90道反演结果与真实模型的相关系数(表 1)，结果表明本文方法的弹性三参数结果与真实模型在纵向和横向上均具有较高的吻合度。

在实际反演中，观测误差或噪声也可引起道间数据差异，为使分析结果更符合实际情况，选用加噪数据进行测试。对合成数据分别添加信噪比SNR=10、5、2 dB的高斯白噪声。对比常规模拟退火与本文方法的反演结果，加噪数据的纵波速度反演结果(图 8)表明，噪声幅度增大会明显导致反演效果变差，但在相同噪声情况下本文方法相比常规方法均有所改善(见箭头指示处)。由反演结果的相关系数(表 2)可知，数据SNR由10 dB降至2 dB，常规方法结果准确度(相关系数)降低7.54%，而本文方法仅降低5.69%，这表明本文方法在一定程度上能克服加噪数据产生的道间差异，提高了反演结果的抗噪性。然而，需要指出的是，引起道间数据差异的原因除观测噪声外，还包括波场干扰、处理参数等因素，因此差异性分析存在一定局限。

3 实际资料应用

为进一步验证本文方法的有效性，选取了二维实际地震数据进行应用测试。实际资料来自于中国南海M探区，该二维剖面共317个地震道集，每个地震道集的角度范围为3°~36°(角度间隔为3°)，采样间隔为2 ms。图 9展示了该二维地震剖面及井旁地震道集。井口位于第98道(图 9a三角箭头指示处)，已探明目的层深度位于1640~1725 ms，为碳酸盐岩(以石灰岩为主)含气层。

采用本文方法对实际地震数据进行反演。对比图 10可见，本文方法的纵波速度反演结果(图 11a)能够恢复部分层间薄层(箭头指示处)；密度反演结果(图 11b)的异常值(箭头指示处)明显减小，横向连续性有明显的改善。本文方法的纵波速度和密度井旁反演结果与测井曲线相关系数分别为0.723和0.668，相比常规方法井旁结果相关系数(0.706和0.641)有所提高。

采用常规模拟退火算法对实际地震数据进行反演。各地震道采用固定优化参数(t⁽⁰⁾=0.9，ΔV_P=0.06、0.03、0.02 km/s)，最大迭代次数为2×10⁴次。由图 10可见，纵波速度(图 10a)和密度(图 10b)反演结果与测井曲线基本吻合，但道间连续性较差，特别是密度反演结果在第230道与第290道附近有明显的异常值，这是固定优化参数在相应地震道处设置不合理所致。

在实际数据反演中，本文方法(含线性反演)和常规方法分别耗时1.85、2.25 h。由于本文方法的迭代终止条件因地震道而异，相比常规方法，总迭代次数有所减少，因而计算效率有一定程度的提升。然而，鉴于地震资料与测井资料的采样率存在差异(两者分别为2、1 ms)，另外还存在地震数据带限因素的影响，本文方法反演结果难以还原测井曲线中的高频信息。此外，由于全局优化算法的随机性，反演结果在目的层底部仍有明显的异常值。

4 结束语

本文联合贝叶斯线性反演与模拟退火非线性反演，提出一种基于自适应优化参数模拟退火的叠前地震联合反演方法。通过线性反演结果，根据局部势能估算初始温度和扰动范围，实现模拟退火优化参数逐道自适应设置，提高了算法的适用性及稳定性。模型数据测试与实测资料应用表明，该方法有效提升了模拟退火算法的收敛效果，比常规方法的弹性参数反演结果更准确。

自适应优化参数设置难以解决全局寻优过程的随机性，本文方法反演结果仍存在一定程度的不稳定，可考虑采用增加总迭代次数或有效调节正则化参数等方式加以缓解。虽然本文测试分析表明扰动范围对收敛效果的影响大于初始温度，但两者对优化过程呈现耦合作用，其影响规律还有待进一步深入研究。此外，作为一种联合反演方法，本文选用常规贝叶斯线性反演驱动后续非线性反演，但线性反演方法对最终结果的影响并未讨论，例如线性结果需满足何种要求才可有效提升最终结果精度，仍有待进一步研究。