基于Zoeppritz方程近似式的油气预测方法及系统

摘要

本发明提供一种基于Zoeppritz方程近似式的油气预测方法及系统，方法包括：采集当前目标区的测井数据以及叠前角道集数据；根据测井数据以及叠前角道集数据提取地震子波；根据叠前角道集数据、地震子波以及Zoeppritz方程近似式得到反演目标函数；根据第一关系、第二关系以及Zoeppritz方程近似式确定目标Zoeppritz方程近似式；根据反演目标函数以及目标Zoeppritz方程近似式确定纵波模量以及横波模量；根据当前目标区地层的纵波模量以及横波模量确定流体因子；对流体因子进行综合解释，得到当前目标区地层的油气藏预测结果。可稳定求取纵波、横波模量信息，提升了流体因子的精度，进而可降低油气藏预测的风险。

说明书

技术领域

本发明关于石油地球物理勘探技术领域，特别是关于石油、天然气能源 的预测技术，具体的讲是一种基于Zoeppritz（佐普利兹）方程近似式的油气 预测方法及系统。

背景技术

在石油地球物理勘探技术领域中，振幅随角度变化（Amplitude versus  angle，AVA）反演是利用反射地震波振幅随入射角的变化规律反演地下岩性 参数，进而估算储层岩石及流体的性质并进行油气预测的技术。Zoeppritz （1919）利用反射界面两侧位移和应力的连续性作为边界条件，推导出了反 射系数和透射系数，它们是入射角θ和介质弹性参数（如：纵波速度V_p、横波 速度V_s、密度ρ、纵波模量M、横波模量μ、拉梅参数λ和泊松比σ等）的函 数。

Zoeppritz方程的线性近似是AVA反演的理论基础。Aki和Richard（阿基 和理查德）（1980）给出了以纵横波速度和密度相对变化表示的反射系数近似 方程，可用于反演地下岩层纵横波速度和密度。Shuey（舒伊）（1985）研究 了泊松比对反射系数的影响，提出了反射系数的AVA截距和梯度概念，可用 于反演泊松比参数。Gidlow（盖德露）（1992）给出了以纵横波阻抗反射系数 表示的近似方程，用于反演纵横波阻抗。Zong（宗兆云）和Yin（印兴耀）（2012） 给出了以纵横波模量反射系数表示的近似方程，以反演纵横波模量。将上述 各种Zoeppritz方程近似式反演，得到各种地层弹性参数，进而可求取多种流 体因子并进行油气藏预测。

油气流体因子的求取过程通常可分为反射系数形式和阻抗形式。常用的 阻抗形式流体因子有Goodway（古德威）（1997）给出的和其中I_p和I_s分别为纵波和横波阻抗。λρ可很好地识别岩性，而μρ可预测流体 变化。Russell（鲁赛尔）（2003，2011）给出了一种新的流体因子其中c是干燥岩石的(V_p/V_s)²系数。Quakenbush（曲肯布什）等（2006）给出 了一种泊松比阻抗流体因子PI=I_p-γ·I_s。Liu（刘苗）和Gong（龚同举）（2012） 给出了一种流体检测因子该流体检测因子其实是 Russell和Quakenbush流体因子的一种组合。上述提及的几种流体因子，均不 能通用于所有实际情况油气藏预测。

因此，对于实际地震勘探工区，如何反演多种地下岩性参数，计算多种 流体因子，并联合起来进行油气预测，降低油气预测的风险是本领域亟待解 决的技术难题。

发明内容

为了克服现有技术存在的上述问题，本发明提供了一种基于Zoeppritz方 程近似式的油气预测方法及系统，在测井和岩石物理分析前提下，充分考虑 了多种敏感流体因子，得到了新的与纵横波模量反射系数有关的Zoeppritz方 程近似式，进一步确定出高精度的流体因子，达到降低油气藏预测风险的目 的。

本发明的目的之一是，提供一种基于佐普利兹Zoeppritz方程近似式的油 气预测方法，包括：采集当前目标区的测井数据以及叠前角道集数据；根据 所述的测井数据以及叠前角道集数据提取地震子波；根据所述的叠前角道集 数据、所述的地震子波以及Zoeppritz方程近似式得到反演目标函数；根据纵 波速度、纵波模量以及密度的弹性参数关系推导纵波模量反射系数、纵波速 度反射系数与密度反射系数的关系，称为第一关系；根据横波速度、横波模 量以及密度的弹性参数关系推导横波模量反射系数、横波速度反射系数与密 度反射系数的关系，称为第二关系；根据所述的第一关系、第二关系以及 Zoeppritz方程近似式确定目标Zoeppritz方程近似式；根据所述的反演目标函 数以及目标Zoeppritz方程近似式确定纵波模量以及横波模量；根据当前目标 区地层的纵波模量以及横波模量确定流体因子；对所述的流体因子进行综合 解释，得到当前目标区地层的油气藏预测结果。

本发明的目的之一是，提供了一种基于佐普利兹Zoeppritz方程近似式的 油气预测系统，包括：数据采集装置，用于采集当前目标区的测井数据以及 叠前角道集数据；地震子波确定装置，用于根据所述的测井数据以及叠前角 道集数据提取地震子波；目标函数确定装置，用于根据所述的叠前角道集数 据、所述的地震子波以及Zoeppritz方程近似式得到反演目标函数；第一关系 推导装置，用于根据纵波速度、纵波模量以及密度的弹性参数关系推导纵波 模量反射系数、纵波速度反射系数与密度反射系数的关系，称为第一关系； 第二关系推导装置，用于根据横波速度、横波模量以及密度的弹性参数关系 推导横波模量反射系数、横波速度反射系数与密度反射系数的关系，称为第 二关系；目标方程近似式确定装置，用于根据所述的第一关系、第二关系以 及Zoeppritz方程近似式确定目标Zoeppritz方程近似式；流体因子确定装置， 用于根据所述的反演目标函数以及目标Zoeppritz方程近似式确定纵波模量以 及横波模量，根据当前目标区地层的纵波模量以及横波模量确定流体因子； 综合解释装置，用于对所述的流体因子进行综合解释，得到当前目标区地层 的油气藏预测结果。

本发明的有益效果在于，在测井和岩石物理分析前提下，充分考虑了多 种敏感流体因子，通过建立新的与纵横波模量反射系数有关的Zoeppritz方程 近似式以及反演目标函数，可稳定求取纵横波模量信息，提升了流体因子的 计算精度，进而可降低油气藏预测的风险。

为让本发明的上述和其他目的、特征和优点能更明显易懂，下文特举较 佳实施例，并配合所附图式，作详细说明如下。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实 施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面 描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲， 在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例提供的一种基于Zoeppritz方程近似式的油气预测方 法的流程图；

图2为图1中的步骤S103的具体流程图；

图3为本发明实施例提供的一种基于Zoeppritz方程近似式的油气预测系 统的结构框图；

图4为本发明实施例提供的一种基于Zoeppritz方程近似式的油气预测系 统中目标函数确定装置的结构框图；

图5(a)为实施方式一中的纵波速度楔形模型图；

图5(b)为实施方式一中的纵波速度反射系数楔形模型图；

图5(c)为实施方式一中的横波速度楔形模型图；

图5(d)为实施方式一中的横波速度反射系数楔形模型图；

图5(e)为实施方式一中的无噪声情况下的叠前角道集示意图；

图5(f)为实施方式一中的含20%随机噪声的叠前角道集示意图；

图5(g)为实施方式一中的时间厚度20ms的无噪叠前角道集示意图；

图5(h)为实施方式一中的时间厚度20ms的含噪叠前角道集示意图；

图6(a)为实施方式一中的Zhang和Sen的叠前AVA反演结果中的纵波反 射系数反演结果示意图；

图6(b)为实施方式一中的Zhang和Sen的叠前AVA反演结果中纵波反射 系数反演残差剖面图；

图6(c)为实施方式一中的Zhang和Sen的叠前AVA反演结果中横波反射 系数反演结果示意图；

图6(d)为实施方式一中的Zhang和Sen的叠前AVA反演结果中横波反射 系数反演残差剖面图；

图6(e)为实施方式一中的本发明的高分辨率叠前AVA反演结果中纵波反 射系数反演结果示意图；

图6(f)为实施方式一中的本发明的高分辨率叠前AVA反演结果中纵波反 射系数反演残差剖面图；

图6(g)为实施方式一中的本发明的高分辨率叠前AVA反演结果中横波反 射系数反演结果示意图；

图6(h)为实施方式一中的本发明的高分辨率叠前AVA反演结果中横波反 射系数反演残差剖面图；

图7(a)为实施方式一中的Zhang和Sen的叠前AVA反演结果中纵波反射 系数反演结果示意图；

图7(b)为实施方式一中的Zhang和Sen的叠前AVA反演结果中纵波反射 系数反演残差剖面图；

图7(c)为实施方式一中的本发明的高分辨率叠前AVA反演结果中纵波反 射系数反演结果示意图；

图7(d)为实施方式一中的本发明的高分辨率叠前AVA反演结果中纵波反 射系数反演残差剖面图；

图8(a)为实施方式二中的纵波速度楔形模型图；

图8(b)为实施方式二中的纵波模量反射系数示意图；

图8(c)为实施方式二中的横波速度楔形模型图；

图8(d)为实施方式二中的横波模量反射系数示意图；

图8(e)为实施方式二中的密度楔形模型图；

图8(f)为实施方式二中的密度反射系数示意图；

图8(g)为实施方式二中的无噪声情况下的叠前角道集示意图；

图8(h)为实施方式二中的含10%随机噪声的叠前角道集示意图；

图8(i)为实施方式二中的时间厚度10ms的无噪叠前角道集示意图；

图8(j)为实施方式二中的时间厚度10ms的含噪叠前角道集示意图；

图9(a)为实施方式二中采用Zong和Yin近似式反演的纵波模量反射系数 示意图；

图9(b)为实施方式二中采用本发明近似式反演的纵波模量反射系数示意 图；

图9(c)为实施方式二中采用Zong和Yin近似式反演的横波模量反射系数 示意图；

图9(d)为实施方式二中采用本发明近似式反演的横波模量反射系数示意 图；

图9(e)为实施方式二中采用Zong和Yin近似式反演的密度反射系数示意 图；

图9(f)为实施方式二中采用本发明近似式反演的密度反射系数示意图；

图9(g)为实施方式二中采用Zong和Yin近似式反演的纵波模量示意图；

图9(h)为实施方式二中采用本发明近似式反演的纵波模量示意图；

图9(i)为实施方式二中采用Zong和Yin近似式反演的横波模量示意图；

图9(j)为实施方式二中采用本发明近似式反演的横波模量示意图；

图9(k)为实施方式二中采用Zong和Yin近似式反演的密度示意图；

图9(l)为实施方式二中采用本发明近似式反演的密度示意图；

图10(a)为反演输入的叠前AVA角道集数据示意图；

图10(b)为反演的纵波速度图；

图10(c)为反演的横波速度图；

图10(d)为反演的密度图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行 清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而 不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做 出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明提出了一种新的与纵横波模量反射系数有关的Zoeppritz方程近似 式，主要用于提升计算Russell流体因子的精度。在测井和岩石物理分析前提 下，充分考虑多种敏感流体因子，包括基于本发明的Zoeppritz方程近似式计 算的流体因子，达到降低油气藏预测风险的目的。

图1为本发明实施例提供的一种基于Zoeppritz方程近似式的油气预测方 法的流程图，由图1可知，该方法具体包括：

S101：采集当前目标区的测井数据以及叠前角道集数据。

在油气勘探中，测井数据以及叠前角道集数据中保留了丰富的信息，如 何利用这些信息来获取目标区地层参数为油气检测和识别提供数据支撑是地 震勘探最为重要的工作之一。目前勘探中主要利用反射地震波振幅随入射角 的变化规律反演地下岩性参数，估算储层岩石及流体的性质并进行油气预测， 即AVA反演技术。

S102：根据所述的测井数据以及叠前角道集数据提取地震子波。

其中，叠前角道集数据用s(t,θ)来表示，提取的地震子波用w(t,θ)来表示， 其中t为地震数据采样时间序列，θ为叠前角道集所对应的角度。在具体的实 施方式中，该步骤首先由叠前角道集数据确定地震子波的振幅谱，然后由叠 前角道集数据中的井旁道数据以及测井数据确定地震子波的相位谱，最后， 根据振幅谱和相位谱确定地震子波。

S103：根据所述的叠前角道集数据、所述的地震子波以及Zoeppritz方程 近似式得到反演目标函数。

基于Zoeppritz方程的线性近似式众多，本发明采用的Zoeppritz方程近似 式是较为经典的一种Aki-Richards近似式，其具体形式如下

$\frac{1}{2}(1+{\tan }^2\bar{\theta })\frac{\Delta V_p}{{\bar{V}}_p}-4{\sin }^2\bar{\theta }\frac{{\bar{V}}_s^2}{{\bar{V}}_p^2}\frac{\Delta V_s}{{\bar{V}}_s}+\frac{1}{2}(1-4{\sin }^2\bar{\theta }\frac{{\bar{V}}_s^2}{{\stackrel{-}{V}}_p^2})\frac{\mathrm{\Delta \rho }}{\bar{\rho }}=r_{\mathrm{pp}}\left(\bar{\theta }\right)---\left(1\right)$

其中，为地层分界面处上下两弹性介质的平均纵波速度； 为地层分界面处上下两弹性介质的平均横波速度；为地层分界面处上下两弹性介质的平均密度；ΔV_p=V_p2-V_p1为上下两弹性介质 纵波速度的差值；ΔV_s=V_s2-V_s1为上下两弹性介质横波速度的差值；Δρ=ρ₂-ρ₁为上下两弹性介质密度的差值；为入射角θ_i和透射角θ_t的平均值， 由于实际入射角θ_i和透射角θ_t相差很小，通常采用入射角θ_i来代替表示地层的横-纵波速度比值；为随角度变化的反射纵波反射系数。

定义纵波反射系数为r_p、横波反射系数为r_s和密度反射系数为r_ρ，则有

$r_p\left(t\right)=\frac{\Delta V_p\left(t\right)}{{\bar{V}}_p\left(t\right)},r_s\left(t\right)=\frac{\Delta V_s\left(t\right)}{{\bar{V}}_s\left(t\right)},r_{\rho }\left(t\right)=\frac{\mathrm{\Delta \rho }\left(t\right)}{\bar{\rho }\left(t\right)}$

${}^{·}·^{·}{r}_p=\frac{{\mathrm{\Delta V}}_p\left(t\right)}{\overline{V_p}\left(t\right)}=\frac{V_{p2}-V_{p1}}{\frac{1}{2}(V_{p2}+V_{p1})}=\frac{{2V}_{p2}-{2V}_{p1}}{V_{p2}+V_{p1}}$

${}_{·}·_{·}{r}_p(V_{p2}+V_{p1})={2V}_{p2}-{2V}_{p1}{}_{·}·_{·}{V}_{p2}\left({2-r}_p\right)=V_{p1}\left({2+r}_p\right)$

${}_{·}·_{·}{V}_{p2}=V_{p1}\frac{{2+r}_p}{{2-r}_p}.$同理可得，$V_{s2}=V_{s1}\frac{2+r_s}{2-r_s},$${\rho }_2={\rho }_1\frac{2+r_{\rho }}{2-r_{\rho }}.$

以上推导表明，由反演得到的反射系数r_p、r_s和r_ρ可递推得到速度V_p、V_s和密度ρ。为了便于介绍本发明的原理，采用如下的简记符号：

$C_p\left(\theta \right)=\frac{1}{2}(1+{\tan }^2\theta ),$$C_s(t,\theta )=-4{\sin }^2\theta \frac{{\bar{V}}_s^2}{{\bar{V}}_p^2\left(t\right)},$$C_{\rho }(t,\theta )=\frac{1}{2}-2{\sin }^2\theta \frac{{\bar{V}}_s^2\left(t\right)}{{\bar{V}}_p^2\left(t\right)}$

此时，Zoeppritz方程Aki-Richards近似式，即公式(1)，可简写为：

C_p(θ)r_p(t)+C_s(t,θ)r_s(t)+C_ρ(t,θ)r_ρ(t)=r_pp(t,θ)  (2)

鉴于大部分Zoeppritz方程近似式具有相同的形式，因此，(2)式可写为

C_a(θ)r_a(t)+C_b(θ)r_b(t)+C_c(θ)r_c(t)=r_pp(t,θ)  (3)

针对某一入射角θ_i，其反射纵波反射系数r_pp(t,θ_i)可表示为：

[C_a(θ_i)  C_b(θ_i)  C_c(θ_i)][r_a(t)  r_b(t)  r_c(t)]^T=r_pp(t,θ_i)  (4)

将反射纵波反射系数r_pp(t,θ_i)同地震子波w(t)褶积，进而可得到这一入射角 θ_i所对应的角道集合成记录s_pp(t,θ_i)：

[C_a(θ_i)  C_b(θ_i)  C_c(θ_i)][r_a(t)  r_b(t)  r_c(t)]^T*w(t)=r_pp(t,θ_i)*w(t)=s_pp(t,θ_i)  (5)

对所有入射角θ_i，θ_i∈[θ₁ θ_N]，N为入射角个数。此时，公式(5)可详写成 矩阵向量方程组的形式，有

$w\left(t\right)*\left(\begin{array}{ccc}{C}_a\left({\theta }_1\right)& C_b\left({\theta }_1\right)& C_c\left({\theta }_1\right)\\ C_a\left({\theta }_2\right)& C_b\left({\theta }_2\right)& C_c\left({\theta }_2\right)\\ ·& ·& ·\\ ·& ·& ·\\ ·& ·& ·\\ C_a\left({\theta }_N\right)& C_b\left({\theta }_N\right)& C_c\left({\theta }_N\right)\end{array}\right)\left(\begin{array}{c}{r}_a\left(t\right)\\ r_b\left(t\right)\\ r_c\left(t\right)\end{array}\right)=\left(\begin{array}{c}{s}_{\mathrm{pp}(t,{\theta }_1)}\\ s_{\mathrm{pp}(t,{\theta }_2)}\\ ·\\ ·\\ ·\\ s_{\mathrm{pp}(t,{\theta }_N)}\end{array}\right)---\left(6\right)$

上述(6)式是时间域AVA反演的基本公式。Zong和Yin（2012）、Alemie 和Sacchi（2011）、Buland和Omre（2003）、Zhang和Sen（2013）等广泛 采用了公式(6)进行AVA反演。在时间域，由于有效信号和干扰噪声信号耦合 在一起，在提高分辨率的同时，不可避免地放大了噪声的影响；另一方面， 其利用了全频段信息，计算成本较高。

所以，本发明给出了下面的基于叠前地震资料部分谱的AVA反演方法， 其利用信噪比较高的优势频段，反演结果具有全息性，计算成本也大大降低。 在原理上也摈弃了Zhang和Sen（2013）等所采用的反射系数奇偶分解理论， 即不是先求取反射系数的奇偶分量再重构反射系数，而是直接反演最终目标 反射系数，原理上更加简捷。

图2为图1中的步骤S103的具体流程图，由图2可知，该步骤具体包括：

S201：对所述的叠前角道集数据以及所述的地震子波进行傅里叶变换。 对公式(5)做快速傅里叶Fourier变换，有

$\mathrm{FFT}\left\{\left(\begin{array}{ccc}{C}_a\left({\theta }_i\right)& C_b\left({\theta }_i\right)& C_c\left({\theta }_i\right)\end{array}\right)\left(\begin{array}{c}{r}_a\left(t\right)\\ r_b\left(t\right)\\ r_c\left(t\right)\end{array}\right)\right\}\times W\left(f\right)=R_{\mathrm{pp}}(f,{\theta }_i)\times W\left(f\right)=S_{\mathrm{pp}}(f,{\theta }_i)---\left(7\right)$

其中，S_pp(f,θ_i)、W(f)和R_pp(f,θ_i)分别是某一入射角θ_i的角道集地震记录 s_pp(t,θ_i)、地震子波w(t)和反射纵波反射系数r_pp(t,θ_i)对应的复数谱。

S202：对傅里叶变换后的叠前角道集数据以及地震子波进行频谱分析， 并选定优势频带范围。

对于某一特定入射角θ_i，按时间t_j,j∈[1,K]，对所选优势频带范围内的频 率成分ω_k=2πf_k，ω_k∈[ω_Low,ω_High]，ω_Low即为所选取的低截角频率，ω_High为所选 的高截角频率，展开有

其中，C_a,b,c(t_k,θ_i)=[C_a(t_k,θ_i) C_b(t_k,θ_i) C_c(t_k,θ_i)]，r_a,b,c(t_k)=[r_a(t_k) r_b(t_k) r_c(t_k)]^T；

为子波频谱矩阵；

为Zoeppritz方程近似式中某一 入射角θ_i所对应的三参数系数矩阵。

S203：根据所述的优势频带范围以及Zoeppritz方程近似式确定频率域反 演模型。

将公式(8)简记为：

${\mathrm{WC}}_{{\theta }_i}m=d_{{\theta }_i}$

联立所有入射角θ_i∈[θ₁ θ_N]，则有：

WCm=d   (9)

其中，m=[r_a,b,c(t₁),r_a,b,c(t₂),…,r_a,b,c(t_K)]^T为要反演的模型参数，比如纵波速度、 横波速度和密度反射系数；$d_{{\theta }_i}={[S({\omega }_1,{\theta }_i),S({\omega }_2,{\theta }_i),···,S({\omega }_M,{\theta }_i)]}^T$为已知的某一入射 角θ_i角道集地震记录的复数谱信息。此时，频率域AVA反演就简化为求上面 的矩阵方程组(9)式。

S204：对所述的频率域反演模型加入1范数正则化，得到反演目标函数。 由于公式(9)是一个高度病态的反问题，因此本发明考虑对模型加入1范数正 则化，即得到的反演目标函数为：

min‖d-WCm‖₂+λ‖m‖₁  (10)

其中，‖d-WCm‖₂用来控制反演过程中模型合成数据WCm同实际观测数据 d之间的匹配，‖m‖₁用来控制反演反射系数的稀疏度，λ用于调节信噪比和分 辨率之间的折中。本发明中，假设地下反射系数序列是由主反射系数r_main（对 应于幅值明显大于0的情况）和背景弱反射系数r_background（对应于幅值近乎为0 的情况，实际地下反射系数一般不为0值）组成，此处的强反射系数个数反 映了稀疏度。

本发明提供的基于Zoeppritz方程近似式的油气预测方法还包括：

S104：根据纵波速度、纵波模量以及密度的弹性参数关系推导纵波模量 反射系数、纵波速度反射系数与密度反射系数的关系，称为第一关系。

由地震波动力学理论，纵波速度V_p同纵波模量M、密度ρ满足弹性参数 关系：

$V_p=\sqrt{\frac{M}{\rho }}---\left(a\right)$

由高等代数全微分相关定理：如果函数z=f(x,y)在点(x,y)处可微分，则 函数在点(x,y)处的偏导数必定存在，且函数在点(x,y)的全微分为：

$\mathrm{dz}=\frac{\partial z}{\partial x}\mathrm{dx}+\frac{\partial z}{\partial y}\mathrm{dy}---\left(b\right)$

所以，利用(b)式，对(a)式两边取微分，得到：

$d{V}_p={\rho }^{-\frac{1}{2}}d{M}^{\frac{1}{2}}+M^{\frac{1}{2}}d{\rho }^{-\frac{1}{2}}---\left(c\right)$

根据反射系数的定义形式，对(c)式两边同除纵波速度V_p，得到：

$\frac{d{V}_p}{V_p}=\frac{{\rho }^{-\frac{1}{2}}d{M}^{\frac{1}{2}}}{V_p}+\frac{M^{\frac{1}{2}}d{\rho }^{-\frac{1}{2}}}{V_p}=\frac{{\rho }^{-\frac{1}{2}}d{M}^{\frac{1}{2}}}{M^{\frac{1}{2}}{\rho }^{-\frac{1}{2}}}+\frac{M^{\frac{1}{2}}d{\rho }^{-\frac{1}{2}}}{M^{\frac{1}{2}}{\rho }^{-\frac{1}{2}}}=\frac{d{M}^{\frac{1}{2}}}{M^{\frac{1}{2}}}+\frac{d{\rho }^{-\frac{1}{2}}}{{\rho }^{-\frac{1}{2}}}---\left(d\right)$

对(d)式进行展开细化和整理，可得：

$\frac{d{V}_p}{V_p}=\frac{\frac{1}{2}{M}^{-\frac{1}{2}}\mathrm{dM}}{M^{\frac{1}{2}}}+\frac{-\frac{1}{2}{\rho }^{-\frac{3}{2}}\mathrm{d\rho }}{{\rho }^{-\frac{1}{2}}}=\frac{1}{2}(\frac{\mathrm{dM}}{M^{\frac{1}{2}}{M}^{\frac{1}{2}}}-\frac{\mathrm{d\rho }}{{\rho }^{-\frac{1}{2}}{\rho }^{\frac{3}{2}}})=\frac{1}{2}(\frac{\mathrm{dM}}{M}-\frac{\mathrm{d\rho }}{\rho })---\left(e\right)$

所以，纵波模量反射系数同纵波速度反射系数和密度反射系数 之间的关系即第一关系为：

$\frac{\Delta V_p}{\overline{V_p}}=\frac{1}{2}(\frac{\mathrm{\Delta M}}{\bar{M}}-\frac{\mathrm{\Delta \rho }}{\bar{\rho }})---\left(11\right)$

S105：根据横波速度、横波模量以及密度的弹性参数关系推导横波模量 反射系数、横波速度反射系数与密度反射系数的关系，称为第二关系。

由地震波动力学理论，横波速度V_s同横波模量μ、密度ρ满足弹性参数关 系：

$V_s=\sqrt{\frac{\mu }{\rho }}---\left(f\right)$

所以，利用(b)式，对(f)式两边取微分，得到：

$d{V}_s={\rho }^{-\frac{1}{2}}d{\mu }^{\frac{1}{2}}+{\mu }^{\frac{1}{2}}d{\rho }^{-\frac{1}{2}}---\left(g\right)$

根据反射系数的定义形式，对(g)式两边同除横波速度V_s，得到：

$\frac{d{V}_s}{V_s}=\frac{{\rho }^{-\frac{1}{2}}d{\mu }^{\frac{1}{2}}}{V_s}+\frac{{\mu }^{\frac{1}{2}}d{\rho }^{-\frac{1}{2}}}{V_s}=\frac{{\rho }^{-\frac{1}{2}}d{\mu }^{\frac{1}{2}}}{{\mu }^{\frac{1}{2}}{\rho }^{-\frac{1}{2}}}+\frac{{\mu }^{\frac{1}{2}}d{\rho }^{-\frac{1}{2}}}{{\mu }^{\frac{1}{2}}{\rho }^{-\frac{1}{2}}}=\frac{d{\mu }^{\frac{1}{2}}}{{\mu }^{\frac{1}{2}}}+\frac{d{\rho }^{-\frac{1}{2}}}{{\rho }^{-\frac{1}{2}}}---\left(h\right)$

对上式进行展开细化和整理，可得：

$\frac{d{V}_s}{V_s}=\frac{\frac{1}{2}{\mu }^{-\frac{1}{2}}\mathrm{d\mu }}{{\mu }^{\frac{1}{2}}}+\frac{-\frac{1}{2}{\rho }^{-\frac{3}{2}}\mathrm{d\rho }}{{\rho }^{-\frac{1}{2}}}=\frac{1}{2}(\frac{\mathrm{d\mu }}{{\mu }^{\frac{1}{2}}{\mu }^{\frac{1}{2}}}-\frac{\mathrm{d\rho }}{{\rho }^{-\frac{1}{2}}{\rho }^{\frac{3}{2}}})=\frac{1}{2}(\frac{\mathrm{d\mu }}{\mu }-\frac{\mathrm{d\rho }}{\rho })---\left(i\right)$

所以，横波模量反射系数同横波速度反射系数和密度反射系数之间的关系即第二关系为：

$\frac{\Delta V_s}{\overline{V_s}}=\frac{1}{2}(\frac{\mathrm{\Delta \mu }}{\bar{\mu }}-\frac{\mathrm{\Delta \rho }}{\bar{\rho }})---\left(12\right)$

S106：根据所述的第一关系、第二关系以及Zoeppritz方程近似式确定目 标Zoeppritz方程近似式。

将公式(11)和(12)代入公式(1)，经整理即可得本发明的目标Zoeppritz方程 近似式：

$\frac{1}{4}(\frac{\mathrm{\Delta M}}{\bar{M}}+\frac{\mathrm{\Delta \rho }}{\bar{\rho }})+[\frac{1}{4}(\frac{\mathrm{\Delta M}}{\bar{M}}-\frac{\mathrm{\Delta \rho }}{\bar{\rho }})-2{\left(\frac{\overline{V_s}}{V_p}\right)}^2\frac{\mathrm{\Delta \mu }}{\bar{\mu }}]{\sin }^2\theta +\frac{1}{4}(\frac{\mathrm{\Delta M}}{\bar{M}}-\frac{\mathrm{\Delta \rho }}{\bar{\rho }}){\tan }^2\theta {\sin }^2\theta =r_{\mathrm{pp}}\left(\theta \right)---\left(13\right)$

公式(13)直接建立了纵波反射系数与纵波模量反射系数、横波模量反射系 数以及密度反射系数一种近似的线性关系。

S107：根据所述的反演目标函数以及目标Zoeppritz方程近似式确定纵波 模量以及横波模量。在求解公式(10)时有很多种算法，下面以基追踪（Basis  pursuit，BP）算法为例，进行说明。

针对本发明的近似式(13)，采用下面符号简记：

$A=\frac{1}{4}(\frac{\mathrm{\Delta M}}{\bar{M}}+\frac{\mathrm{\Delta \rho }}{\bar{\rho }})$$B=\frac{1}{4}(\frac{\mathrm{\Delta M}}{\bar{M}}-\frac{\mathrm{\Delta \rho }}{\bar{\rho }})-2{\left(\frac{\overline{V_s}}{\overline{V_p}}\right)}^2\frac{\mathrm{\Delta \mu }}{\bar{\mu }}$$C=\frac{1}{4}(\frac{\mathrm{\Delta M}}{\bar{M}}-\frac{\mathrm{\Delta \rho }}{\bar{\rho }})$

将上式A、B和C代入公式(13)，有

A+Bsin²θ+Ctan²θsin²θ=R_pp(θ)  (14)

此时，公式(9)中的C=[1 sin²θ tan²θsin²θ]，m=[A B C]^T。求解公式 (10)，可获得系数A、B和C。然后通过下列的转换：

$\frac{\mathrm{\Delta M}}{\bar{M}}=2(A+C)$$\frac{\mathrm{\Delta \rho }}{\stackrel{-}{\rho }}=2(A-C)$$\frac{\mathrm{\Delta \mu }}{\bar{\mu }}=\frac{1}{2}{\left(\frac{V_p}{V_s}\right)}^2(C-B)$

即可得到纵波模量反射系数和横波模量反射系数最后，联同低 频初始模型，即可获得纵波模量M和横波模量μ。此处提及的低频初始模型 可由步骤S101采集的测井数据以及叠前角道集数据得到，此过程为现有技术， 故在本发明中不再赘述。

S108：根据当前目标区地层的纵波模量以及横波模量确定流体因子。表1 展示了几种不同的弹性参数对油气检测的敏感度。其中，AI(p)=ρV_p为纵波阻 抗，AI(s)=ρV_s为横波阻抗，为饱和岩石纵波模量，与介质抗压缩性 和硬度直接相关，体现储层骨架和流体信息；为饱和岩石横波模量， 与介质抗剪切性和刚度直接相关，体现储层骨架信息。(V_p/V_s)_dry为干燥岩石纵 横波速度比。

从表1可明显看出，湿砂岩和气砂岩的ρf信息相比于其它参数来说差异 较大。因此，本发明根据当前目标区地层的纵横波模量信息计算油气高敏感 度流体因子ρf，进行油气藏预测。其中，流体因子

表1

若利用Zong和Yin（2012）近似式

$\frac{1}{4}(1+{\tan }^2\theta )\frac{\mathrm{\Delta M}}{\bar{M}}-2{\sin }^2\theta {\left(\frac{\overline{V_s}}{\overline{V_p}}\right)}^2\frac{\mathrm{\Delta \mu }}{\bar{\mu }}+\frac{1}{4}(1-{\tan }^2\theta )\frac{\mathrm{\Delta \rho }}{\bar{\rho }}=r_{\mathrm{pp}}\left(\theta \right)---\left(15\right)$

此时，公式(9)中的$C=\left(\begin{array}{ccc}\frac{1}{4}(1+{\tan }^2\theta )& {-2\sin }^2\theta {(\overline{V_s}/\overline{V_p})}^2& \frac{1}{4}(1-{\tan }^2\theta )\end{array}\right),$$m={\left(\begin{array}{ccc}\frac{\mathrm{\Delta M}}{\bar{M}}& \frac{\mathrm{\Delta \mu }}{\bar{\mu }}& \frac{\mathrm{\Delta \rho }}{\bar{\rho }}\end{array}\right)}^T.$求解公式(10)，可获得纵波模量、横波模量和密度反射 系数。然后结合低频初始模型，也可反演得到纵波模量M和横波模量μ。

利用本发明的公式(13)可直接求取纵波模量M和横波模量μ，而不用先采 用Aki-Richards近似式求取纵波速度V_p、横波速度V_s和密度ρ，再进一步得 到纵波模量和横波模量这在一定程度上减少了累积误差。而 在不给定情况下，利用公式(15)则无法进行反演，但利用本发明近似 式(14)可以先反演得到A、B和C，这使得A和C中的和以及B的反 演结果不受估计误差的影响。而利用公式(15)，若有误差，其 对和均有影响。这也是本发明的优异之处。

S109：对所述的流体因子进行综合解释，得到当前目标区地层的油气藏 预测结果。利用前面反演所得的地层弹性参数：纵波速度V_p、横波速度V_s、 密度ρ、纵波模量M和横波模量μ等计算Goodway（1997）给出的流体因子 和λρ可很好地识别岩性，而μρ可预测流体变化；计算 Russell（2003，2011）给出的流体因子计算Quakenbush（2006） 给出的泊松比阻抗流体因子PI=I_p-γ·I_s；计算Liu（刘苗）和Gong（龚同举） （2012）给出的流体因子$F=(I_p-C_1·I_s)(I_p^2-C_2·I_s^2)$等。对不同的流体因子属性 剖面进行解释分析，联合不同流体因子属性剖面标定具有共性的油气异常区 域，确定油气储层位置，降低利用单一流体因子进行油气藏预测的风险。

图3为本发明实施例提供的一种基于Zoeppritz方程近似式的油气预测系 统的结构框图，由图3可知，所示系统具体包括：

数据采集装置100，用于采集当前目标区的测井数据以及叠前角道集数 据。在油气勘探中，测井数据以及叠前角道集数据中保留了丰富的信息，如 何利用这些信息来获取目标区地层参数为油气检测和识别提供数据支撑是地 震勘探最为重要的工作之一。目前勘探中主要利用反射地震波振幅随入射角 的变化规律反演地下岩性弹性参数，估算储层岩石及流体的性质并进行油气 预测，即AVA反演技术。

地震子波确定装置200，用于根据所述的测井数据以及叠前角道集数据提 取地震子波。其中，叠前角道集数据用s_pp(t,θ)来表示，提取的地震子波用w(t,θ) 来表示，其中t为地震数据采样时间序列，θ为叠前角道集所对应的角度。在 具体的实施方式中，首先由叠前角道集数据确定地震子波的振幅谱，然后由 叠前角道集数据中的井旁道数据以及测井数据确定地震子波的相位谱，最后， 根据振幅谱和相位谱确定地震子波。

目标函数确定装置300，用于根据所述的叠前角道集数据、所述的地震子 波以及Zoeppritz方程近似式得到反演目标函数；

第一关系推导装置400，用于根据纵波速度、纵波模量以及密度的弹性参 数关系推导纵波模量反射系数、纵波速度反射系数与密度反射系数的关系， 称为第一关系，第一关系如公式(11)所示。

第二关系推导装置500，用于根据横波速度、横波模量以及密度的弹性参 数关系推导横波模量反射系数、横波速度反射系数与密度反射系数的关系， 称为第二关系，第二关系如公式(12)所示。

目标方程近似式确定装置600，用于根据所述的第一关系、第二关系以及 Zoeppritz方程近似式确定目标Zoeppritz方程近似式，目标Zoeppritz方程近 似式如公式(13)所示。

流体因子确定装置700，用于根据所述的反演目标函数以及目标Zoeppritz 方程近似式确定纵波模量以及横波模量，根据当前目标区地层的纵波模量以 及横波模量确定流体因子。在求解公式(10)时有很多种算法，诸如基追踪（Basis  pursuit，BP），即可获得纵波模量M和横波模量μ。计算流体因子 $\mathrm{\rho f}=\rho [M_{\mathrm{sat}}-{(V_p/V_s)}_{\mathrm{dry}}^2\times {\mu }_{\mathrm{sat}}],$其中，$M_{\mathrm{sat}}=\rho V_p^2,$${\mu }_{\mathrm{sat}}=\rho V_s^2.$

综合解释装置800，用于对所述的流体因子进行综合解释，得到当前目标 区地层的油气藏预测结果。利用反演所得的地层弹性参数：纵波速度V_p、横 波速度V_s、密度ρ、纵波模量M和横波模量μ等；计算Goodway（1997）给 出的流体因子和λρ可很好地识别岩性，而μρ可预测流体 变化；计算Russell（2011）给出的流体因子计算Quakenbush （2006）给出的泊松比阻抗流体因子PI=I_p-γ·I_s；计算Liu（刘苗）和Gong （龚同举）（2012）给出的流体因子等。对不同的流 体因子属性剖面进行解释分析，联合不同流体因子属性剖面标定具有共性的 油气异常区域，确定油气储层位置，降低利用单一流体因子进行油气藏预测 的风险。

图4为本发明实施例提供的一种基于Zoeppritz方程近似式的油气预测系 统中目标函数确定装置300的结构框图，由图4可知，目标函数确定装置具 体包括：

傅里叶变换单元301，用于对所述的叠前角道集数据以及所述的地震子波 进行傅里叶变换。

分析提取单元302，用于对傅里叶变换后的叠前角道集数据以及地震子波 进行频谱分析，并选定优势频带范围；

反演模型确定单元303，用于根据所述的优势频带范围以及Zoeppritz方 程近似式确定频率域反演模型，频率域反演模型如公式(9)所示。

目标函数确定单元304，用于对所述的频率域反演模型加入1范数正则化， 得到反演目标函数。反演目标函数如公式(10)所示。

下面结合具体的实施例，详细介绍本发明提供的基于Zoeppritz方程近似 式的油气预测方法以及系统的可行性。

图5为本实施例提供的原始模型与合成数据，模型和数据主要用于比较 本发明的叠前频率域AVA反演同Zhang和Sen（2013）提出的方法反演纵、 横波反射系数和密度反射系数的效果。其中，图5(a)为纵波速度楔形模型图， 图5(b)为纵波速度反射系数楔形模型图，图5(c)为横波速度楔形模型图，图 5(d)为横波速度反射系数楔形模型图，图5(e)为无噪声情况下的叠前角道集示 意图，图5(f)为含20%随机噪声的叠前角道集示意图，图5(g)为时间厚度20ms 的无噪叠前角道集示意图，图5(h)为时间厚度20ms的含噪叠前角道集示意图。 为公平比较，模型数据源同于Zhang和Sen（2013）Geophysics论文。

图6为无噪情况下基于本发明的叠前频率域AVA反演结果同Zhang和Sen （2013）反演结果的对比，其中，图6(a)为Zhang和Sen的叠前AVA反演结果 中的纵波反射系数反演结果示意图，图6(b)为Zhang和Sen的叠前AVA反演 结果中纵波反射系数反演残差剖面图，图6(c)为Zhang和Sen的叠前AVA反 演结果中横波反射系数反演结果示意图，图6(d)为实施方式一中的Zhang和 Sen的叠前AVA反演结果中横波反射系数反演残差剖面图，图6(e)为本发明 的高分辨率叠前AVA反演结果中纵波反射系数反演结果示意图，图6(f)为本 发明的高分辨率叠前AVA反演结果中纵波反射系数反演残差剖面图，图6(g) 为本发明的高分辨率叠前AVA反演结果中横波反射系数反演结果示意图，图 6(h)为本发明的高分辨率叠前AVA反演结果中横波反射系数反演残差剖面图。 Zhang和Sen反演结果图源同于其Geophysics论文。通过对比分析，无噪情 况下，本发明的反演结果分辨薄层能力较好，反演残差剖面较为干净，适用 于目前石油工业高分辨率勘探需求。

图7为加入20%随机噪声情况下，本发明的叠前AVA反演结果同Zhang 和Sen（2013）反演结果的对比。图7(a)为实施方式一中的Zhang和Sen的叠 前AVA反演结果中纵波反射系数反演结果示意图；图7(b)为实施方式一中的 Zhang和Sen的叠前AVA反演结果中纵波反射系数反演残差剖面图；图7(c) 为实施方式一中的本发明的高分辨率叠前AVA反演结果中纵波反射系数反演 结果示意图；图7(d)为实施方式一中的本发明的高分辨率叠前AVA反演结果 中纵波反射系数反演残差剖面图。Zhang和Sen反演结果图源同于其 Geophysics论文。通过对比分析，含噪情况下，本发明的反演残差剖面较为 干净。

图8为原始模型和合成数据。该模型用于比较基于本发明的Zoeppritz方 程近似式同基于Zong和Yin（2012）近似式反演的纵、横波模量和密度的效 果。图8(a)为为纵波速度楔形模型图，图8(b)为纵波模量反射系数示意图，图 8(c)为横波速度楔形模型图，图8(d)为横波模量反射系数示意图，图8(e)为密 度楔形模型图，图8(f)为密度反射系数示意图，图8(g)为无噪声情况下的叠前 角道集示意图；图8(h)为含10%随机噪声的叠前角道集示意图；图8(i)为时间 厚度10ms的无噪叠前角道集示意图，图8(j)为时间厚度10ms的含噪叠前角 道集示意图。模型数据来源于第一类AVO纵横波速度和密度数据，第一类 AVO异常也是目前石油工业勘探过程中遇到较多的油气异常。

图9为加入10%随机噪声情况下基于本发明的Zoeppritz方程近似式同基 于Zong和Yin（2012）近似式反演结果的对比。图9(a)为采用Zong和Yin近 似式反演的纵波模量反射系数示意图；图9(b)为采用本发明近似式反演的纵 波模量反射系数示意图；图9(c)为采用Zong和Yin近似式反演的横波模量反 射系数示意图；图9(d)为采用本发明近似式反演的横波模量反射系数示意图； 图9(e)为采用Zong和Yin近似式反演的密度反射系数示意图；图9(f)为采用 本发明近似式反演的密度反射系数示意图；图9(g)为采用Zong和Yin近似式 反演的纵波模量示意图；图9(h)为采用本发明近似式反演的纵波模量示意图； 图9(i)为采用Zong和Yin近似式反演的横波模量示意图；图9(j)为采用本发 明近似式反演的横波模量示意图；图9(k)为采用Zong和Yin近似式反演的密 度示意图；图9(l)为采用本发明近似式反演的密度示意图。本模型数据中，纵 波模量反射系数为-0.0840，幅值最强；密度反射系数为-0.0486，幅值次之； 横波模量反射系数为0.0101，幅值最弱。加10%随机噪声情况下，横波模量 反射系数最易受噪声的影响，对比分析采用本发明的Zoeppritz方程近似式结 果和采用Zong和Yin（2012）近似式反演结果，可知本发明近似较为稳定， 特别是横波模量和密度信息反演结果。

图10为本发明采用叠前地震数据部分谱信息进行AVA反演并油气预测 的二维剖面效果图。图10(a)为反演输入的叠前AVA角道集数据示意图，图 10(b)为反演的纵波速度图，图10(c)为反演的横波速度图，图10(d)为反演的 密度图。由图10(a)叠前角道集数据，由于地震子波的干涉效应，分辨率较低， 顶底界面难以识别，油气藏边界范围难以划定；而本发明AVA反演的纵横波 速度和密度剖面，去掉了地震子波的干涉效应，提高了分辨率，综合利用各 种流体因子进行油气检测，油气储层刻画较为清楚。

综上所述，本发明首先给出了基于叠前地震资料部分谱的高分辨率AVA 反演框架。主要包括：详细推导了采用部分谱信息进行AVA反演的原理公式， 给出了AVA反演的目标函数矩阵方程组，介绍了地下岩性参数反演的1范数 正则化方法。

其次，本发明给出了纵横波模量反射系数同纵横波反射系数关系的详细推 导过程，发明了一种新的由纵横波模量反射系数和密度反射系数表示的 Zoeppritz方程近似式。该发明可稳定求取纵横波模量信息，提升了计算Russell 流体因子的精度。

最后，利用本发明给出的AVA反演框架，选用本发明的Zoeppritz方程 近似式和前面所介绍的Zoeppritz方程近似式，反演计算不同流体因子所需要 的地下岩石弹性参数，包括：纵波速度V_p、横波速度V_s、纵波阻抗I_p、横波 阻抗I_s、密度ρ、纵波模量M、横波模量μ和泊松比σ，等。综合利用多种流 体因子对薄层、薄互层油气藏进行预测，以达到降低油气藏预测风险的目的。

本发明的有益成果是：

1）本发明给出了纵横波模量反射系数同纵横波反射系数关系的详细推导 过程，发明了一种新的由纵横波模量反射系数和密度反射系数表示的 Zoeppritz方程近似式。该发明近似式可稳定求取纵横波模量信息，提升了计 算Russell流体因子的精度。

2）选用本发明的Zoeppritz方程近似式可以稳定求取纵横波模量信息，提 升了计算Russell流体因子的精度。

3）在测井和岩石物理分析前提下，充分考虑多种敏感流体因子，包括基 于本发明的Zoeppritz方程近似式计算的流体因子，可降低油气藏预测的风险。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流 程，可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于一 般计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施 例的流程。其中，所述的存储介质可为磁碟、光盘、只读存储记忆体（Read-Only  Memory，ROM）或随机存储记忆体（Random Access Memory，RAM）等。

本领域技术人员还可以了解到本发明实施例列出的各种功能是通过硬件 还是软件来实现取决于特定的应用和整个系统的设计要求。本领域技术人员 可以对于每种特定的应用，可以使用各种方法实现所述的功能，但这种实现 不应被理解为超出本发明实施例保护的范围。

本发明中应用了具体实施例对本发明的原理及实施方式进行了阐述，以 上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于 本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上 均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。