波阻抗反演是LINDSETH^[1]于1979年首次提出的方法, 通过测井数据合成的地震记录进行地震资料解释, 并且把波阻抗改写成反射系数的积分形式。经过数十年的发展, 波阻抗反演已成为岩石属性、储层表征和沉积相分析等领域中最重要的工具^[2]。反演方法从确定性反演发展为随机性反演, 从线性反演发展为非线性反演, 从叠后反演发展为叠前反演^[3-4]。其中, 叠前反演是根据纵横波速度和密度计算弹性参数, 达到划分岩性和检测流体的目的。传统的叠前反演是在确定性反演的基础上, 根据已有的测井信息建立初始模型, 利用Zoeppritz方程的线性近似公式, 结合叠前地震记录估算纵横波阻抗或弹性阻抗^[5-6]。早期的确定性反演存在抗噪性较差和过度依赖初始模型等问题^[7], 而基于贝叶斯框架下的叠前三参数同步反演可以在地震数据和测井数据等先验信息约束下对反演的结果进行不确定性评估^[8-9]。在反演过程中, 允许先验信息加入贝叶斯函数, 然后得到最大似然分布以及后验分布曲线。将得到的最大似然分布和后验分布曲线应用于概率分析和储层描述^[10-11]。通常的贝叶斯反演在优化过程中, 采用共轭梯度或者拟牛顿法在空间中线性地寻找全局最优位置。这种方法很容易陷入局部最小值, 使得在优化过程中不能跳出局部最小值寻找全局最小, 因而得不到正确结果。基于随机搜索的叠前同步三参数反演可以克服共轭梯度法的缺陷, 充分利用先验信息, 输出若干符合地震地质信息的结果^[12]。反演目标函数中的正演算子, 早期多采用线性近似方程, 而目前多采用Zoeppritz方程, 进而保证叠前反演的精度。在随机反演中, 关键问题之一就是优化目标函数的效率。尽管人们提出了许多非线性随机算法解决优化问题, 但这些方法效率较低, 不能满足实际生产需求。

自适应粒子群优化算法(adaptive particle swarm optimization, APSO)是一种改进的非线性优化方法, 已经被证明比标准粒子群优化更具成本效益^[13]。在保证反演精度的同时, 针对叠前反演中非线性优化的效率问题提出了利用自适应粒子群优化算法。该方法可以提供与地震反演最佳解一致的有限模型空间中全局最小值解, 并且具有良好的抗噪性。克服了确定性反演的问题, 即过度拟合和良好初始模型的需求。此外, 该方法允许贝叶斯统计和马尔可夫链法在数据域每一次迭代中更新粒子群。由于在宽角度区域, 传统的线性近似公式与Zoeppritz公式存在误差, 在反演中不能够满足精度条件。因此, 在正演过程中选择了Zoeppritz方程而非其近似公式。最后, 利用NS工区的二维海洋地震资料和测井资料, 估算纵横波速度和密度, 并且与传统的随机反演方法对比, 验证了该方法的有效性和可行性。

1 方法技术 1.1 叠前三参数同步反演

叠前三参数同步反演, 即利用叠前地震道集中包含的振幅和角度信息同时反演纵波速度、横波速度和密度3个参数, 进而用于指示储层的物性和含油气性^[14]。基于完整形式的Zoeppritz方程叠前三参数反演, 可以精确地描述振幅随入射角变化, 并且在入射角接近临界角时依然可以保证反演的精度^[15-16]。由于波阻抗反演问题总是不适定的, 并且地震数据的带宽有限以及缺少低频分量, 因此几乎不可能从地震道中直接估算绝对阻抗^[17-18]。因此, 在贝叶斯框架下的叠前三参数同步反演, 可以根据生成的概率分布, 减少反演结果的不确定性^[19-20]。在贝叶斯框架的基础上, 结合随机反演的特点, 建立了叠前三参数同步反演流程, 如图 1所示。首先, 对叠前地震数据进行处理得到角道集, 对测井数据进行处理得到低频背景模型。此外, 利用测井信息建立纵、横波速度以及密度的扰动先验分布。利用先验分布结合Metropolis-Hasting方法实现随机抽样, 建立初始粒子群。初始粒子群经过不断优化, 最终得到3个参数的全局最优值和后验分布。

在贝叶斯框架下的叠前三参数同步反演中, 后验分布由先验分布表示, 根据给定的叠前地震数据D_obs, 那么三参数的后验分布p(m| D_obs)可以表示为:

$ p\left(m \mid \boldsymbol{D}_{\mathrm{obs}}\right)=\frac{p\left(\boldsymbol{D}_{\mathrm{obs}} \mid m\right) p(m)}{p\left(\boldsymbol{D}_{\mathrm{obs}}\right)} $

(1)

式中: p(m)是扰动模型m的先验概率分布。在三参数叠前反演的问题中, m为纵横波速度和密度。假定纵横波速度和密度服从高斯分布, 那么模型的先验分布表达式为:

$ p(m)=\frac{1}{\sqrt{2 \pi \sigma_{m}^{2}}} \exp \left[-\frac{\left(m-\mu_{m}\right)^{2}}{2 \sigma_{m}^{2}}\right] $

(2)

式中: σ_m是模型的协方差; μ_m是模型的期望值。在贝叶斯框架中, 构造似然函数是检验模型选择和评估模型不确定性考虑的另一种方法。在此函数中, 如果假设噪声均值为零, 且具有噪声协方差的高斯噪声分布, 那么噪声的分布表达式为:

$ \varepsilon \sim N\left(0, \sigma_{N}^{2}\right) $

(3)

因此, 根据叠前角道集D_obs和噪声协方差σ_N, 那么似然函数p(D_obs|m)的表达式为:

$ \begin{array}{c} p\left(\boldsymbol{D}_{\mathrm{obs}} \mid m\right)=p_{\varepsilon}\left(\boldsymbol{D}_{\mathrm{obs}}-R(m) S_{w}\right)=\frac{1}{\sqrt{2 \pi \sigma_{N}^{2}}} \exp \cdot \\ {\left[-\frac{\left(\boldsymbol{D}_{\mathrm{obs}}-\boldsymbol{R}(m) \boldsymbol{S}_{w}\right)^{\mathrm{T}}\left(\boldsymbol{D}_{\mathrm{obs}}-\boldsymbol{R}(m) \boldsymbol{S}_{w}\right)}{2 \sigma_{N}^{2}}\right]} \end{array} $

(4)

式中: R(m)是Zoeppritz纵波反射系数方程; S_w是利用角道集提取的各个入射角子波。

1.2 自适应粒子群

在随机反演的优化过程中, 粒子群优化算法(article swarm optimization, PSO)是被广泛应用于目标函数的优化和机器学习的一种方法^[21]。该算法通过研究群体中个体之前的信息共享和协作在限定的空间中寻找最优解。虽然粒子群算法已被广泛应用于解决各种复杂的工程问题, 但在实际问题中仍需研究克服粒子群算法计算量大和精度低等缺点^[22]。尤其在地球物理反演中, 针对高维度参数空间的反演, 需要更为高效的算法。其中, 自适应粒子群(APSO)算法加入了自适应参数和精英学习策略(elite learning strategy, ELS)。与传统粒子群优化算法相比, 自适应粒子群算法在求解精度、收敛速度和算法可靠性等方面更具优势。由于收敛速度比传统的粒子群算法更快, 因而自适应粒子群算法在全局搜索上具有更高的计算效率。在叠前三参数同步反演中, 自适应粒子群算法在优化目标函数和寻找地震反演问题的全局最优解方面比粒子群优化算法有更好的计算能力^[23]。

在粒子群优化的过程中, 将粒子的当前最优值和全局最优值分别表示为p_Best[i]和p_Best[g]。假设粒子群中粒子的个数为N, 那么将第i次迭代中某个粒子的当前空间位置表示为x_i。图 2显示了空间Ω中第i次迭代粒子的当前位置和更新后的位置。由图 2可见, 更新速度取决于粒子的当前最优值和全局最优值以及上次的更新速度。在传统粒子群优化算法的迭代中, 粒子更新速度和更新后的新位置可以写成:

$ \begin{array}{c} \boldsymbol{v}_{i}=\omega \boldsymbol{v}_{i-1}+c_{1} * \operatorname{rand}_{1} *\left(p_{\text {Best }[i]}-x_{i}\right)+ \\ c_{2} * \operatorname{rand}_{2} *\left(p_{\text {Best }[g]}-x_{i}\right) \end{array} $

(5)

$ x_{i+1}=\boldsymbol{v}_{i}+x_{i} $

(6)

式中: x_i为当前粒子的位置; x_i+1为迭代更新速度后粒子的位置; v为更新速度; c₁和c₂为学习因子; rand₁和rand₂为[0, 1]之间的随机数; ω为惯性因子。

而自适应粒子群算法在粒子群优化算法的基础上, 修正了惯性因子, 并且增加了精英学习策略。图 3描述了在反演过程中, 自适应粒子群算法的具体工作流程, 在迭代过程中每个粒子会不断移动, 直到满足迭代的终止条件。在叠前三参数同步反演中, 实现自适应粒子群算法包括以下几个步骤。

首先是粒子的初始化阶段。对于每一个粒子, 里面包含有纵波速度、横波速度以及密度3个分量。根据纵波速度、横波速度以及密度的先验分布, 构建粒子的初始状态。针对粒子中每一道中的采样点, 粒子的初始状态可以表示为:

$ v_{\mathrm{P}}^{0}=L_{v_{\mathrm{P}}}+\Delta_{v_{\mathrm{P}}} $

(7)

$ v_{\mathrm{S}}^{0}=L_{v_{\mathrm{S}}}+\Delta_{v_{\mathrm{S}}} $

(8)

$ \rho^{0}=L_{\rho}+\Delta_{\rho} $

(9)

式中: {v_P⁰, v_S⁰, ρ⁰}为单个初始的粒子状态; {L_vP, L_vS, L_P}为低频背景模型, {Δ_vP, Δ_vS, Δ_ρ}为模型的扰动变量。那么在第i次迭代过程中, 第k个粒子可以表示为:

$ x_{i}^{k}=\left\{v_{\mathrm{P}}^{k}, v_{\mathrm{S}}^{k}, \rho^{k}\right\} $

(10)

在粒子群演化的过程中, 用参数f_d表示粒子群的演化状态。在迭代中, 计算粒子k到其它粒子j的距离可以表示为:

$ d_{\mathrm{i}}=\frac{1}{N-1} \sum\limits_{j=1, j \neq k}^{N} \sqrt{\sum\limits_{s=1}^{D}\left[x_{i}^{k}(s)-x_{i}^{j}(s)\right]^{2}} $

(11)

式中: N代表粒子群中总的粒子数; D是单个地震道参与反演的采样点数。之后比较所有粒子的距离, 对距离d_i从大到小进行排序, 计算演化因子f_d, 定义如下:

$ f_{\mathrm{d}}=\frac{d_{\mathrm{g}}-d_{\min }}{d_{\max }-d_{\min }} \in[0, 1] $

(12)

式中: d_g为全局最佳位置所对应的距离。f_d在初始阶段的默认值为0, 可以通过演化因子判断当前粒子群的状态。图 4为不同的演化因子所对应的演化状态, 共4个状态, 即跳出阶段, 探索阶段, 开发阶段以及收敛阶段。不同的演化因子对应着不同的探索阶段。

然后, 计算自适应惯性参数, 该参数可以平衡搜索局部最佳和全局最佳位置的能力。

$ \begin{array}{c} \omega\left(f_{\mathrm{d}}\right)=\frac{1}{1+1.5 e^{-2.6 f_{\mathrm{d}}}} \in[0.4, 0.9] \\ \forall f_{\mathrm{d}} \in[0, 1] \end{array} $

(13)

式中: ω(f_d)在每次迭代中从0.9减小到0.4。因此, 由于演化因子的存在, ω(f_d)使得粒子群优化算法变化更为有效。由于f_d初始值为0, 因此ω(f_d)的初始值为0.4。

在迭代过程中, 需要定义惯性系数c₁和c₂, 如:

$ c_{i}(i+1)=c_{i}(i) \pm \delta_{c}, i=1, 2 $

(14)

每一次迭代中, c₁和c₂的值都在变化。为了收敛, c₁越大, c₂的值越小; 反之, c₁的值越小, c₂的值越大。c₁和c₂在[1.5, 2.5]范围内变化, 它们的总和应在[3.0, 4.0]范围内, 以便粒子在局部最小值和全局最小值之间进行调整, δ_c是在[0.05, 0.10]范围内均匀生成的随机值。而在传统的粒子群算法中, c₁和c₂等于2, 不具备调整的能力。在迭代过程中, 采取精英学习策略(ELS)寻找全局最佳位置。标准的粒子群优化算法如果陷入局部最优, 就不能帮助全局最优从当前最优中跳出来。因此, 提出利用精英学习策略来克服这一不足。当从每一次迭代中找出全局最佳时, 它将增加轻微的扰动量, 通过高斯扰动利用学习策略在空间中移动, 如:

$ p^{k}=p^{k}+\left(X_{\max }^{k}-X_{\min }^{k}\right) * \operatorname{Gauss}\left(\mu, \sigma^{2}\right) $

(15)

式中: Gauss(μ, σ²)表示期望μ和标准差σ的高斯分布。X_max^k和X_min^k可以从模型空间中计算。在叠前三参数反演中, 我们通常假设变量服从高斯分布, 即:

$ \sigma=\sigma_{\max }-\left(\sigma_{\max }-\sigma_{\min }\right) *(i / T) $

(16)

式中: σ_max=1.0和σ_min=0.1是从实证研究测试得出的。i为当前的迭代次数, T为总的迭代次数。如果新的p^k比当前位置更适合, 那么它将被接受。当满足预定最小值或达到最大迭代次数时的收敛准则时, 输出全局最优解和若干次优解。

1.3 基于粒子群优化-马尔科夫链蒙特卡洛的随机模拟反演实现过程

在叠前三参数同步反演中, 已经假设模型的扰动分布服从高斯分布。根据三参数扰动先验分布抽样获取初始粒子群^[24]。首先, 利用测井曲线中的三参数信息, 去掉低频成分, 生成三参数扰动曲线。通过对扰动进行统计, 建立三参数样本点x和样本集X。

典型的蒙特卡洛方法以相对较小的步长围绕着分布移动, 而这些步长没有具备朝着相同方向进行的趋势。尽管这些方法易于实现和分析, 但可能需要很长时间才能探索空间中的最佳位置^[25]。因此为了加速马尔科夫链蒙特卡罗算法(MCMC), 提出了与自适应粒子群优化算法相结合的马尔科夫链蒙特卡罗算法, 即自适应粒子群优化-马尔科夫链蒙特卡罗算法(APSO-MCMC)。这样在保证随机性的同时, 进一步提高全局搜索的计算效率。在随机模拟的过程中, 马尔科夫链蒙特卡罗算法需要一个实用的解决方案来解决生成跃迁矩阵的问题。跃迁矩阵与自适应粒子群优化算法相关联, 在粒子群中根据每个粒子的目标函数值优选后, 将优选得到的粒子按照概率p₁继续在空间中搜索, 而剩下的粒子则以一定概率p₂加入随机分量, 如果在加入随机分量后目标函数值小于之前, 则保留当前结果; 否则停留在先前结果。按照这个方法构建跃迁矩阵, 在加速搜索过程的同时, 保持了其自身的收敛性。在本方法中, 跃迁矩阵的表达式为:

$ \boldsymbol{P}=\left[\begin{array}{cc} p_{1} & 1-p_{1} \\ 1-p_{2} & p_{2} \end{array}\right] $

(17)

式中: p₁和p₂为自身转移的概率, 其大小可以根据实际情况定义。当以p₁的概率搜索时, 表达式为:

$ x_{i}=x_{i}+p_{1} * \operatorname{rand}(0, 1) * V_{i} $

(18)

式中的V_i由自适应粒子群优化算法在每次迭代过程中获得。图 5描述了自适应粒子群优化-马尔科夫链蒙特卡罗算法的简要技术流程。基于目标函数和过渡状态, 可以计算新概率, 然后确定过程是否继续。主要步骤为:

1) 从三参数的扰动分布中随机采样, 生成初始状态X₀, 并对初始粒子群进行优选;

2) 在迭代时, 对每个粒子在转化过程中应用跃迁矩阵P, 然后计算目标函数的值。按目标函数值降序对粒子排序, 然后输出最佳粒子;

3) 粒子群应用自适应粒子群优化算法方法, 获得局部最优p_Best[i]和全局最优p_Best[g]。计算自适应惯性系数以确定演化状态ω(f_d), 然后更新粒子迭代后的位置X_M+1和速度V_M+1;

4) 根据目标函数的值确定每个粒子移动到下一个位置或停止移动;

5) 确定收敛准则。如果全局最佳p_Best[g]小于预设定值ε或达到最大迭代次数后, 则终止迭代并转到步骤7)。如果未达到收敛准则, 进行步骤6);

6) 计算更新粒子的分布, 增加迭代次数, 然后返回步骤2);

7) 输出最佳粒子, 然后获得最大似然结果以及后验分布。

2 应用实例 2.1 NS工区概述

利用本方法对NS工区二维海洋地震数据进行反演。由于海洋地震数据包含大量多次波, 因而在反演前需要对数据进行预处理。在预处理后, 通过叠前时间偏移得到共反射点道集和含有角度信息的叠前数据体。测井解释的结果表明, 在2 000 m左右含有丰富的含气显示。反演之前, 根据实际地下情况建立含气饱和地层模型, 然后利用Zoeppritz方程分析工区含气储层的AVO效应。根据测井信息将模型简化成一个双层模型, 上层为泥岩, 下层为含气砂岩, 并且假定双层模型的各层是均匀介质。根据模型得到的正演模拟结果如图 6所示, 其曲线特征符合第三类AVO响应, 即地震波的负振幅随着反射角的增加而显著增加。将此特征关系与实际的叠前道集中的含气储层区段的振幅比较后, 振幅特征与AVO曲线特征一致, 说明工区内发育有丰富的含气储层。而后利用井旁道的低频信息, 结合叠前数据, 对三参数最大后验估计进行计算, 并和实际的测井曲线对比, 估算反演的精度和误差。完成井旁道反演后, 再对整个剖面进行反演。此外, 利用同样的数据, 将本方法与传统的粒子群算法(PSO)和遗传算法(GA)对比, 观察同等迭代次数的情况下不同方法的计算效率。

2.2 反演结果分析

在反演过程中, 首先对井旁道进行反演, 以验证反演结果的有效性。通过提取井旁道的道集, 结合测井曲线中的低频信息, 建立扰动分布, 生成初始粒子群后进行叠前三参数反演。图 7为纵波速度的扰动统计结果和先验模型, 利用统计结果建立针对纵波速度扰动的高斯分布。用同样的方法建立横波速度和密度的扰动分布。利用该高斯分布建立先验模型, 然后对该分布进行随机采样, 生成若干初始扰动模型, 形成初始粒子群。粒子群由40个粒子构成。通过对整个粒子群的反演, 不断估算演化因子, 判定收敛状态。图 8显示了演化因子随迭代次数变化的情况。由图 8可见, 随着迭代次数的增加, 演化因子在初始阶段逐步增加, 而后逐渐变小。说明在反演初期阶段, 粒子群中的粒子在空间中不断搜索全局最优位置。缩小搜索的范围后, 状态逐渐稳定, 演化因子变小, 进入收敛状态。

从全部40个粒子中选择5个最优解。图 9为5个粒子的纵波反演结果。尽管粒子之间的反演结果不一致, 但这些最优解的粒子在搜索空间中彼此接近, 并且反演结果都符合岩石物理模型。从输出模型的后验分布中选择最大后验估计(MAP), 并将其与测井曲线进行比较。图 10为反演得到的最大后验估计与实际真值和低频模型的对比。纵波速度和密度显示出与实际测井信息对比有良好的一致性, 总体误差小于5%, 并且反演结果不过分依赖于低频模型。虽然反演得到的横波速度的误差相比纵波速度和密度大, 但是整体趋势仍然与实际真值保持一致。将测井解释结果与泊松比曲线对比发现, 三参数反演计算得到的泊松比, 在气层处显示有低泊松比的特征。因此利用这一特征, 可以帮助识别含气储层。此外, 提取了粒子群中的某个粒子的目标函数值和全局最优值, 其对比情况如图 11所示。粒子在不断收敛, 最后粒子群中的所有粒子都会在搜索空间中趋于一个方向。在反演过程中追踪粒子群中的第30个粒子可以发现, 在迭代初期阶段目标函数值为34.12, 随着迭代次数的增加, 目标函数的值受全局最小值的影响不断减少, 并且逐步收敛。而目标函数的全局最小值在反演初期从最初的23.91经过500次迭代后逐渐降低至7.36, 进一步证实了该算法的有效性。图 12为反演后的合成地震记录(蓝色)与实际地震记录(红色)。对比发现, 相关系数达到了0.97, 反演结果符合测井信息和地震信息。在验证了井旁道反演的有效性后, 而后针对工区的整条剖面进行反演。图 13、图 14和图 15分别为反演得到的二维剖面上的纵、横波速度和密度。从二维剖面上可以看出, 反演结果在保留低频信息的同时, 增加了很多细节, 并且与地震波形的变化特征保持一致。图 16为根据三参数反演结果计算的泊松比剖面, 结合井上含气储层的泊松比特征, 可以进一步利用泊松比指示目的层含气储层的存在。

此外, 通过与传统的粒子群算法和遗传算法对比后发现, 该算法具有更快的收敛速率, 表 1对比了3种方法在同一目标函数下不同迭代次数的耗时情况, G_best代表每次迭代过程中目标函数的值; t代表迭代所需要的时间。在迭代初期, 3种方法差距不大, 并且粒子群优化算法的单次迭代时间要少于遗传算法的单次迭代时间。随着迭代次数的增加, 可以发现, 自适应粒子群优化-马尔科夫链蒙特卡罗算法的优化速率明显快于其它两种算法。虽然自适应粒子群优化-马尔科夫链蒙特卡罗算法的计算时间比粒子群优化算法所需要的时间略长, 但该方法能够在固定的迭代次数内更加快速收敛, 并且在固定的迭代次数下, 自适应粒子群优化-马尔科夫链蒙特卡罗算法花费的时间比遗传算法花费的时间少。与常规粒子群优化算法相比, 虽然每一步中增加了演化因子的计算, 潜在地增加了时间, 但是, 与粒子群优化算法相比, 自适应粒子群优化-马尔科夫链蒙特卡罗算法的收敛更快, 能够在空间中更快地搜索全局最小值。针对高维度的反演问题, 需要更加快速的收敛, 因此自适应粒子群优化-马尔科夫链蒙特卡罗算法相比粒子群优化算法更高效。

3 结论

本文讨论了传统的叠前波阻抗地震反演方法存在的不足, 在此基础上分析了同步随机反演的优势。在利用Zoeppritz方程保证随机反演精度的同时, 实现了基于自适应粒子群优化-马尔科夫链蒙特卡罗反演算法, 该方法用于对NS地区实际地震资料的测试结果表明:

1) 在构建正演算子过程中利用Zoeppritz公式而非近似公式, 保证了反演的精度, 反演过程中不断随机寻找全局最优值而非确定性方法, 确保了反演的可靠性;

2) 基于自适应粒子群优化-马尔科夫链蒙特卡罗算法的反演算法能够在有限的空间中寻找全局最优位置, 方法中包含的演化因子可以根据其大、小判定搜索状态, 方法中的精英策略能够帮助粒子跳出局部最优位置;

3) NS地区二维叠前数据的应用结果表明该地区的含气储层符合第三类AVO特征响应, 经过本方法的测试验证了应用的有效性和可行性。将本方法与传统的粒子群算法和遗传算法进行对比测试后发现, 新的方法具有更高的计算效率。

未来将岩性、深度和地层反射系数等约束信息加入到方法中, 可以使反演的结果与沉积相、地震相和岩相获得的认知更为一致, 进一步提高反演的可靠性和实用性。