一种高分辨率的多尺度波动方程反演方法

摘要

一种高分辨率的多尺度波动方程反演方法，本发明提供了一种多尺度波动方程反演方法，其利用地震炮集，基于声波波动方程，对给定的初始体积模量和密度进行反演。反演过程在不同频率、空间尺度中进行，并利用先验模型对目标层进行约束，通过迭代获得体积模量和密度的反演剖面，增强了反演的适定性，提高了反演迭代收敛速度。

说明书

技术领域

本发明涉及地质勘探技术领域，具体的讲是一种高分辨率的多尺度波动方程反演 方法。

背景技术

储层物性特征参数控制着油气的空间分布和油气的储量估计，直接影响油气开采 方案的设计和采收率，一直是油气勘探开发研究的一个重点。随着现阶段油气田勘探 开发难度的增大和计算技术的发展，近年来，储层预测及流体检测地球物理方法与技 术的研究取得了显著进展。地震反演是储层预测及流体检测技术系列之一。

地震反演是利用野外观测的地震数据，通过数学方法及计算技术求取地质目标的 物理参数。双程声波波动方程能够模拟声学介质中的地震波形，与传统的反演（波阻 抗反演、弹性阻抗反演和AVO反演）相比，该方程能处理任意复杂的非均匀介质情 况，提高反演精确性。

声波波动方程反演面临着两大问题：一、计算效率低是实际工业生产面临的瓶颈 性问题。解决这一问题关键是提高计算波动方程模拟的计算效率或者提高反演的收敛 速度。Tarantola等(1988)提出了时间-深度模拟地震波场，Pratt等（1998）提出的频 率空间域的模拟单频波场。这两种方法减少了一次反演迭代所需要时间。另外，高效 的并行计算技术大幅提高了计算速度。二、反演分辨率较低使波动方程反演不能很好 应用到储层反演。地震数据的主频主要为30Hz-50Hz，如果只利用地震数据反演，最 终反演结果分辨率也比较低。

发明内容

为了解决现有技术的问题，提高地震数据反演的效率，提供了一种高分辨率的多 尺度波动方程反演方法，能够快速收敛，提高反演的效率。

本发明实施例提供了一种高分辨率的多尺度波动方程反演方法，包括：

步骤101，采集炮集数据；

步骤102，输入地层层位数据H(x)；

步骤103，建立约束模型m_p(x,z)＝(K_p,ρ_p),其中，输入的模型是所需反演的 储层段的体积模量K_p和密度ρ_p；

步骤104，记每个炮集数据为一个二维数据P_o(t,x)，其中，x代表距离，t代 表采样时间点，设炮集个数为n_s，对炮集数据P_o(t,x)进行多尺度分解，得到在不同 主频尺度下的地震炮集P_o(t,x,ω_i)，ω_i＝i·Δω,多尺度分解公式可表示为：

$G\left(\omega \right)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{{\omega }^2}{{2\sigma }^2}}$

其中ω表示高斯函数的中心频率；σ表示高斯函数的宽度，多尺度下地震炮集可表示 为

P_o(t,x,ω_i)＝FFT^-1[FFT(P(t,x))·G(ω_i)]

对每个尺度的炮集数据P_o(t,x,ω_i)，重复下述步骤105至步骤111；

步骤105，设初始模型的输入m₀(x,z)＝(K₀,ρ₀)，K₀为所需反演储层段的初始体 积模量，ρ₀为所需反演储层段的初始密度；

步骤106，针对频率ω_i，以主频为ω_i的零相位雷克子波，根据声波波动方程及吸 收边界条件计算模拟炮集数据P_s(t,x,ω_i)；

步骤107，计算接收到的地震炮集数据P_o(t,x,ω_i)与模拟炮集数据P_s(t,x,ω_i)的残 差d(t,x)，记为

d(t,x)＝P_o-P_s       （4）

其中，P_o为P_o(t,x,ω_i)，P_s为P_s(t,x,ω_i)；

步骤108，以-d(-t,x)为震源，替换声波波动方程的震源项，模拟残差记录的反 向炮集数据P_s'；

步骤109，根据上述计算得到的炮集数据和反向炮集数据构造反演模型变化的梯 度：

$g=\frac{\partial }{\partial m}E\left(p\right)=<\mathrm{\delta \rho },\frac{\partial }{\partial t}{P}_s·\frac{\partial }{\partial t}{P_s}^{\prime }>+<\mathrm{\delta K},▿P_s·▿{P_s}^{\prime }>---\left(5\right)$

E(p)为炮集数据能量误差，上述式（5）的反演模型参数的共轭修改量可表示为：

$\delta {\rho }_{i_{\omega },j_s}=\frac{1}{{\rho }^2\left(x\right)}{\int }_0^T\frac{\partial }{\partial t}P(x,t;x_s;{\omega }_i)\frac{\partial }{\partial t}{P}^{\prime }(x,t;x_s;{\omega }_i)\mathrm{dt}---\left(6\right)$

$\delta K_{i_{\omega },j_s}=\frac{1}{K^2\left(x\right)}{\int }_0^T▿P_s(x,t;x_s;{\omega }_i)·▿{P_s}^{\prime }(x,t;x_s;{\omega }_i)\mathrm{dt}---\left(7\right)$

其中，x为接收点的位置，t为时间变量，x_s为震源的位置，▽表示梯度算子， T为炮集记录时间长度，K,ρ为反演模型参数，为反演模型参数修 改量，i_ω表示第i个频率，j_s表示第j炮；对该频率下的炮集中所有炮重复步骤106 至步骤109，累加计算整个模型修改量；

步骤110，利用步骤109计算的所有模型修改量，生成新的反演模型，

所述新的反演模型为：

m_n+1＝m_n-α_ng_n         （8）

其中，m＝(K,ρ)为反演模型，n为迭代次数，α_n为第n次迭代步长，g_n＝(δK,δρ) 为第n次迭代的反演模型变化梯度；

步骤111，按照下式计算能量误差，

$E=\frac{1}{2}\Delta d^t{C}_D^{-1}\mathrm{\Delta d}---\left(9\right)$

式中E为能量误差，Δd^t为Δd的转置，Δd为上述步骤108中的残差，C_D为数 据协方差矩阵；

当E＜δ，i_ω＜n_ω时，输出反演结果，即m＝(K,ρ)新的反演模型作为本频段内 下次反演的初始模型，进入重复步骤105；其中i_ω为第i个频率，n_ω为频率采样个数；

当E＞δ，i_ω＜n_ω时，输出反演结果，即m＝(K,ρ)新的反演模型作为下一频段 反演的初始模型，进入重复步骤104；

当E＜δ，i_ω＝n_ω时，输出反演结果，即m＝(K,ρ)新的反演模型作为最终的反 演结果，其中，δ为任意小的实数。

根据本发明实施例所述一种高分辨率的多尺度波动方程反演方法的一个进一步 的方面，所述步骤101中对采集到的炮集数据进行去噪处理，去除环境噪音和系统噪 音。

根据本发明实施例所述一种高分辨率的多尺度波动方程反演方法的再一个进一 步的方面，步骤108中，所述声波波动方程为：

$\left(\begin{array}{c}\rho \frac{\partial v_x}{\partial t}=\frac{\partial P}{\partial x}\\ \rho \frac{\partial v_z}{\partial t}=\frac{\partial P}{\partial z}\\ \frac{\partial P}{\partial t}={\mathrm{\rho \upsilon }}^2\frac{{\partial v}_x}{\partial x}+{\mathrm{\rho \upsilon }}^2\frac{\partial x_z}{\partial z}+s(x_0,z_0,t)\end{array}\right),---\left(1\right)$

其中，s(x₀,z₀,t)为震源项，P为应力向量，v_x和v_z分别为模拟地震剖面 P_s(t,x,ω_i)，及时间-空间波场P(t,x,z)。

根据本发明实施例所述一种高分辨率的多尺度波动方程反演方法的另一个进一 步的方面，所述步骤108中，采用改进的非分裂PML吸收边界条件，即

$\left(\begin{array}{c}\rho \frac{\partial v_x}{\partial t}=\frac{\partial P}{\partial x}-{\Omega }_{\mathrm{xx}}\\ \rho \frac{\partial v_z}{\partial t}=\frac{\partial P}{\partial z}-{\Omega }_{\mathrm{zz}}\\ \frac{\partial P}{\partial t}={\mathrm{\rho \upsilon }}^2\frac{{\partial v}_x}{\partial x}+{\mathrm{\rho \upsilon }}^2\frac{{\partial v}_z}{\partial z}-{\mathrm{\rho \upsilon }}^2{\Psi }_{\mathrm{xx}}-{\mathrm{\rho \upsilon }}^2{\Psi }_{\mathrm{zz}}\end{array}\right)---\left(2\right)$

其中Ω_xx、Ω_zz、Ψ_xx、Ψ_zz、T_xx、T_zz为引入的中间辅助变量，它们对应的 辅助方程为：

$\left(\begin{array}{c}\frac{{\partial \Omega }_{\mathrm{xx}}}{\partial t}+{\sigma }_x{\Omega }_{\mathrm{xx}}={\sigma }_x\frac{\partial T}{\partial x}\mathrm{xx}\\ \frac{\partial {\Omega }_{\mathrm{zz}}}{\partial t}+{\sigma }_z{\Omega }_{\mathrm{zz}}={\sigma }_z\frac{\partial T_{\mathrm{zz}}}{\partial z}\\ \frac{{\partial \Psi }_{\mathrm{xx}}}{\partial t}+{\sigma }_x{\Psi }_{\mathrm{xx}}={\sigma }_x\frac{{\partial V}_x}{\partial x}\\ \frac{{\partial \Psi }_{\mathrm{zz}}}{\partial t}+{\sigma }_z{\Psi }_{\mathrm{zz}}={\sigma }_z\frac{{\partial V}_z}{\partial z}\end{array}\right)---\left(3\right)$

上述的方程组（2）和（3）构成吸收边界条件的方程，

联立声波波动方程和吸收边界条件的方程得到炮集波场P_s。

通过本发明实施例的多尺度反演，使得波动方程反演稳定、收敛快速。理论模型 测试结果表明了该方法的有效性。实例应用为我国四川油气田实际数据反演，经过5 次迭代后，获得结果较好。

附图说明

结合以下附图阅读对实施例的详细描述，本发明的上述特征和优点，以及额外的 特征和优点，将会更加清楚。

图1所示为本发明实施例一种高分辨率的多尺度波动方程反演方法的流程图；

图2所示为本发明实施例野外炮集数据1000炮中的一炮的示意图；

图3a-图3d为本发明实施例地震数据中4个不同频段的示意图；

图4a为本发明实施例速度初始模型的示意图；

图4b为本发明实施例密度初始模型的示意图；

图5a至图5d为本发明实施例根据四个不同频率得到的体积模量的增量示意图；

图6a至图6d为本发明实施例四个不同频率得到的体积模量反演结果示意图；

图7所示为本发明实施例每个频段经过5次迭代得到的最终反演结果示意图。

具体实施方式

下面的描述可以使任何本领域技术人员利用本发明。具体实施例和应用中所提供 的描述信息仅为示例。这里所描述的实施例的各种延伸和组合对于本领域的技术人员 是显而易见的，在不脱离本发明的实质和范围的情况下，本发明定义的一般原则可以 应用到其他实施例和应用中。因此，本发明不只限于所示的实施例，本发明涵盖与本 文所示原理和特征相一致的最大范围。

如图1所示为本发明实施例一种高分辨率的多尺度波动方程反演方法的流程图。

包括步骤101，采集炮集数据，并对该炮集数据进行去噪处理，去除环境噪音和 系统噪音。

步骤102，输入地层层位数据H(x)。在本发明的技术方案中与现有技术中的反 演不同，在本发明中输入的数据为深度域数据。

步骤103，建立约束模型m_p(x,z)＝(K_p,ρ_p),其中，输入的模型是所需反演的 储层段的体积模量K_p和密度ρ_p；该模型为空间中的参数，在现有技术中的全波形 反演（Full Waveform Inversion）中不需要先验模型。

步骤104，记每个炮集数据为一个二维数据P_o(t,x)，其中，x代表距离，t代 表采样时间点，设炮集个数为n_s。

对炮集P_o(t,x)进行多尺度分解，得到在不同主频尺度下的地震炮集P_o(t,x,ω_i)， ω_i＝i·Δω。多尺度分解公式可表示为：

$G\left(\omega \right)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{{\omega }^2}{{2\sigma }^2}}$

其中ω表示高斯函数的中心频率；σ表示高斯函数的宽度，σ越大，高斯函数越宽。 多尺度下地震炮集可表示为

P_o(t,x,ω_i)＝FFT^-1[FFT(P(t,x))·G(ω_i)]

对每个尺度的炮集数据P_o(t,x,ω_i)，重复下述步骤105至步骤111；

步骤105，设初始模型的输入m₀(x,z)＝(K₀,ρ₀)，K₀为所需反演储层段的初始体 积模量，ρ₀为所需反演储层段的初始密度；这两个模型是根据已知的条件建立的粗 糙光滑的深度域模型。反演技术就要是通过不同的迭代次数来更新此模型。

步骤106，针对频率ω_i，以主频为ω_i的零相位雷克子波，根据声波波动方程及吸 收边界条件计算模拟炮集数据P_s(t,x,ω_i)。

步骤107，计算接收到的地震炮集数据P_o(t,x,ω_i)与模拟炮集数据P_s(t,x,ω_i)的残 差d(t,x)，记为

d(t,x)＝P_o-P_s

其中，P_o为P_o(t,x,ω_i)，P_s为P_s(t,x,ω_i)。

步骤108，以-d(-t,x)为震源，替换声波波动方程的震源项，模拟残差记录的反 向炮集数据P_s'。

步骤109，根据上述计算得到的炮集数据和反向炮集数据构造反演模型变化的梯 度：

$g=\frac{\partial }{\partial m}E\left(p\right)=<\mathrm{\delta \rho },\frac{\partial }{\partial t}{P}_s·\frac{\partial }{\partial t}{P_s}^{\prime }>+<\mathrm{\delta K},▿P_s·▿{P_s}^{\prime }>---\left(5\right)$

E(p)为炮集数据能量误差，上述式（5）的反演模型参数的共轭修改量可表示为：

$\delta {\rho }_{i_{\omega },j_s}=\frac{1}{{\rho }^2\left(x\right)}{\int }_0^T\frac{\partial }{\partial t}P(x,t;x_s;{\omega }_i)\frac{\partial }{\partial t}{P}^{\prime }(x,t;x_s;{\omega }_i)\mathrm{dt}---\left(6\right)$

$\delta K_{i_{\omega },j_s}=\frac{1}{K^2\left(x\right)}{\int }_0^T▿P_s(x,t;x_s;{\omega }_i)·▿{P_s}^{\prime }(x,t;x_s;{\omega }_i)\mathrm{dt}---\left(7\right)$

其中，x为接收点的位置，t为时间变量，x_s为震源的位置，▽表示梯度算子， T为炮集记录时间长度，K,ρ为反演模型参数，为反演模型参数修 改量，i_ω表示第i个频率，j_s表示第j炮。

对炮集中所有炮重复步骤106至步骤109，累加计算整个模型修改量；

步骤110，利用步骤109计算的所有模型修改量，生成新的反演模型，

所述新的反演模型为：

m_n+1＝m_n-α_ng_n        （8）

其中，m＝(K,ρ)为反演模型，n为迭代次数，α_n为第n次迭代步长，g_n＝(δK,δρ) 为第n次迭代的反演模型变化梯度；

步骤111，按照下式计算能量误差，

$E=\frac{1}{2}\Delta d^t{C}_D^{-1}\mathrm{\Delta d}---\left(9\right)$

式中E为能量误差，Δd^t为Δd的转置，Δd为上述步骤108中的残差，C_D为数 据协方差矩阵；

当E＜δ，i_ω＜n_ω时，输出反演结果，即m＝(K,ρ)新的反演模型作为本频段内 下次反演的初始模型，进入重复步骤105；其中，i_ω为第i个频率，n_ω为频率采样个数。

当E＞δ，i_ω＜n_ω时，输出反演结果，即m＝(K,ρ)新的反演模型作为下一频段 反演的初始模型，进入重复步骤104；

当E＜δ，i_ω＝n_ω时，输出反演结果，即m＝(K,ρ)新的反演模型作为最终的反 演结果，其中，δ为任意小的实数，例如可以为10^-6或者10^-5的预定数值。

作为本发明的一个实施例，所述步骤108中的声波波动方程可以为：

$\left(\begin{array}{c}\rho \frac{\partial v_x}{\partial t}=\frac{\partial P}{\partial x}\\ \rho \frac{\partial v_z}{\partial t}=\frac{\partial P}{\partial z}\\ \frac{\partial P}{\partial t}={\mathrm{\rho \upsilon }}^2\frac{{\partial v}_x}{\partial x}+{\mathrm{\rho \upsilon }}^2\frac{\partial x_z}{\partial z}+s(x_0,z_0,t)\end{array}\right),---\left(1\right)$

其中，s(x₀,z₀,t)为震源项，P为应力向量，v_x和v_z分别为模拟地震剖面 P_s(t,x,ω_i)，及时间-空间波场P(t,x,z)。

采用改进的非分裂PML吸收边界条件，即

$\left(\begin{array}{c}\rho \frac{\partial v_x}{\partial t}=\frac{\partial P}{\partial x}-{\Omega }_{\mathrm{xx}}\\ \rho \frac{\partial v_z}{\partial t}=\frac{\partial P}{\partial z}-{\Omega }_{\mathrm{zz}}\\ \frac{\partial P}{\partial t}={\mathrm{\rho \upsilon }}^2\frac{{\partial v}_x}{\partial x}+{\mathrm{\rho \upsilon }}^2\frac{{\partial v}_z}{\partial z}-{\mathrm{\rho \upsilon }}^2{\Psi }_{\mathrm{xx}}-{\mathrm{\rho \upsilon }}^2{\Psi }_{\mathrm{zz}}\end{array}\right)---\left(2\right)$

其中Ω_xx、Ω_zz、Ψ_xx、Ψ_zz、T_xx、T_zz为引入的中间辅助变量，可以在初始 阶段全部充零，它们对应的辅助方程为：

$\left(\begin{array}{c}\frac{{\partial \Omega }_{\mathrm{xx}}}{\partial t}+{\sigma }_x{\Omega }_{\mathrm{xx}}={\sigma }_x\frac{\partial T}{\partial x}\mathrm{xx}\\ \frac{\partial {\Omega }_{\mathrm{zz}}}{\partial t}+{\sigma }_z{\Omega }_{\mathrm{zz}}={\sigma }_z\frac{\partial T_{\mathrm{zz}}}{\partial z}\\ \frac{{\partial \Psi }_{\mathrm{xx}}}{\partial t}+{\sigma }_x{\Psi }_{\mathrm{xx}}={\sigma }_x\frac{{\partial V}_x}{\partial x}\\ \frac{{\partial \Psi }_{\mathrm{zz}}}{\partial t}+{\sigma }_z{\Psi }_{\mathrm{zz}}={\sigma }_z\frac{{\partial V}_z}{\partial z}\end{array}\right)---\left(3\right)$

联立以上两式（1）、（2）和（3）可得到PML波动方程问题的求解，得到解P_s， 即为炮集波场P_s。

通过本发明的实施例，通过多尺度反演使得波动方程反演稳定、收敛快速。理论 模型测试结果表明了该方法的有效性。

如图2所示为本发明实施例野外炮集数据1000炮中一炮的示意图，作为反演的 输入；

步骤201，采集实际地震炮集数据（图2），对炮集数据进行静校正、地标一致性 振幅补偿和叠前去噪处理，得到信噪比较高的地震炮集数据，作为反演的输入。

步骤202，输入地层层位H(x)。本例的层位即为深度域的层位，而非常规叠前 反演所需要的深度域层位。

步骤203，建立约束模型m_p(x,z)＝(K_p,ρ_p),本例中，利用测井资料和层位数 据内插出先验模型（体积模量K_p、密度ρ_p），作为约束：先验模型中包含着已知与 地下真实模型一样的信息，用这些信息控制每次反演的数据不会超出先验模型的数值 一定范围，从而使得反演结果获得更准确的结果。

步骤204，对地震数据根据高斯函数进行多尺度分解，具体分解过程可以参考图 1中实施例相应步骤，分解到4个不同的频段（尺度）（图3a-图3d），对每个频段炮 集数据重复步骤205至步骤211。

步骤205，用层析反演获得比较准确的宏观速度初始速度模型，如图4a所示， 其中，层析反演是一种常规的速度建模方法，在这里不做赘述，通过GAEDENAR公 式获得密度初始模型，如图4b所示；通过模型与密度，计算出初始体积模量模型。

步骤206，针对每个频率ω_i，以主频为ω_i的零相位雷克子波，根据声波波动方程 公式（1）及吸收边界条件的方程公式（2）和（3）计算模拟炮集数据P_s(t,x,ω_i)。

步骤207，计算接收到的地震炮集数据P_o(t,x,ω_i)与模拟炮集数据P_s(t,x,ω_i)的残 差d(t,x)。

步骤208，以-d(-t,x)为震源，替换声波波动方程的震源项，模拟残差记录的反 向炮集数据P_s'；

步骤209，根据上述计算得到的炮集数据和反向炮集数据构造反演模型变化的梯 度对炮集中所有炮重复步骤206至步骤209，累加计算整个模型修改 量。图5a至图5d表示四个不同频率得到的体积模量的增量。

步骤210，利用步骤209计算的所有模型修改量，生成新的反演模型（图6a-图 6d）。图6a至图6d表示四个不同频率得到的体积模量反演结果。

步骤211，用四个频率段，每个频率段迭代五次（也可以设定别的整数），如图7 所示。

通过上述本发明实施例的方法，利用地震炮集，基于声波波动方程，对给定的初 始体积模量和密度进行反演。反演过程在不同频率、空间尺度中进行，并利用先验模 型对目标层进行约束，通过迭代获得体积模量和密度的反演剖面，增强了反演的适定 性，提高了反演迭代收敛速度。本发明通过多尺度反演，使得波动方程反演稳定、收 敛快速。理论模型测试结果表明了该方法的有效性。实例应用为我国四川油气田实际 数据反演，经过5次迭代后，获得结果较好。

本发明可以以任何适当的形式实现，包括硬件、软件、固件或它们的任意组合。 本发明可以根据情况有选择的部分实现，比如计算机软件执行于一个或多个数据处理 器以及数字信号处理器。本文的每个实施例的元素和组件可以在物理上、功能上、逻 辑上以任何适当的方式实现。事实上，一个功能可以在独立单元中、在一组单元中、 或作为其他功能单元的一部分来实现。因此，该系统和方法既可以在独立单元中实现， 也可以在物理上和功能上分布于不同的单元和处理器之间。

在相关领域中的技术人员将会认识到，本发明的实施例有许多可能的修改和组合， 虽然形式略有不同，仍采用相同的基本机制和方法。为了解释的目的，前述描述参考 了几个特定的实施例。然而，上述的说明性讨论不旨在穷举或限制本文所发明的精确 形式。前文所示，许多修改和变化是可能的。所选和所描述的实施例，用以解释本发 明的原理及其实际应用，用以使本领域技术人员能够最好地利用本发明和各个实施例 的针对特定应用的修改、变形。