CN104570085A - 一种纵横波射线参数域联合反演方法 - Google Patents

Abstract

本发明提供了一种纵横波射线参数域联合反演方法，通过将纵波、转换波叠前道集资料转换到射线参数域，利用地震振幅随射线参数变化的信息反演得到地层弹性参数，进而对岩性和流体进行预测。本发明是同点物性联合反演，纵波和转换波来自于同一条入射路径，虽然反射路径彼此分离，但是来自于同一深度同一反射点，各参数之间的联系更加紧密，因此能得到更可靠的地层弹性参数，提高利用纵横波地震资料进行油气预测的精度。

Description

一种纵横波射线参数域联合反演方法

技术领域

本发明属于地震资料解释领域，是一种纵波、转换横波地震资料的联合反演方法。 

背景技术

在油气检测和油藏描述中，介质的弹性参数十分重要，这些参数与岩性和流体成分有关。描述弹性介质的最主要的3个参数是纵波速度、横波速度和密度。目前地球物理勘探的常规技术是利用纵波AVO 

（Amplitude Versus Offset，振幅随炮检距变化）叠前道集反演上述3个参数，反演时考虑了纵波反射振幅随入射角变化的因素。但是大量理论研究和实例应用表明，在有限偏移距的纵波资料中，这3个参数的准确求解存在一定的多解性。随着多波多分量地震勘探技术的不断发展，高质量转换波资料的获取，使得纵波、转换波联合反演技术成为随之发展起来的新技术。利用纵波和转换横波地震数据进行联合反演，可以减少反演结果的多解性，提高反演参数的精度，提供更多油藏描述和监测的潜力。随着近年来多波勘探采集和处理技术的发展，纵波、转换波联合反演的研究得到了较大的发展，出现了多种联合反演方法，通常是基于Zeoppritz方程或其线性近似进行反演，这些方法均利用纵波、转换波地震振幅随入射角的变化信息反演得到地层的弹性参数。纵波、转换波联合反演方法的基本原理是基于Smith和Gidlow（1987）提出的加权叠加方法，通过PP数据的加权叠加来计算纵波和横波波阻抗以及其它的相关参数。Fatti等（1994）做了进一步的研究，取消了对于Gardner公式的依赖并且包含了密度项。Stewart等（2001）、Landro等（2006）将加权叠加的思想用于联合反演，给出了实用的纵横波联合反演方法，反演结果得到的是纵波速度、横波速度和密度的变化率ΔV_p/V_p、ΔV_s/V_s、Δρ/ρ。但这些方法仅使用纵波、转换波地震振幅随入射角的变化信息，没有考虑纵波、转换波地震振幅随射线参数的变化信息，所反演出的结果精度较低，与钻井、测井、试油、生产数据进行对比和统计分析的结果吻合率不高。 

纵波、转换波地震振幅随射线参数的变化信息是多波地震勘探的重要参数，用纵波、转换波地震振幅随射线参数的变化信息进行联合反演方面目前还缺乏类似的研究。由于纵波、转换波地震振幅随射线参数的变化信息中包含着更准确的地震能量传递信息，因此纵横波射线参数域联合反演方法的研究和应用对于更充分地利用纵横波地震信息进行油气检测有着重要意义。 

发明内容

针对现有技术的不足，本发明的目的在于提供一种利用纵波、转换波地震振幅随射线参数的变化信息进行联合反演的方法。 

本发明采用的技术方案是：一种纵横波射线参数域联合反演方法，其特征在于：所述方法通过将纵波、转换波叠前道集资料转换到射线参数域，利用地震振幅随射线参数变化的信息反演得到地层弹性参数，进而对岩性和流体进行预测。 

具体包括以下步骤： 

步骤一：采用射线追踪法将PP波、PS波共成像点道集数据从偏移距域转换到射线参数域，并进行分射线参数叠加，得到分射线参数叠加剖面； 

射线追踪法是指给定发射点和接收点位置及介质的波速，求从发射点到接收点的射线轨迹。射线追踪方法有很多种，本发明采用弯曲射线追踪法，根据最小走时准则迭代求出接收点处的射线路径，然后将偏移距域共成像点道集中某一确定时间和偏移距的反射同相轴，放置到相应射线参数值的相同时间处，即完成共成像点道集的射线参数域转换。利用得到的射线参数域共成像点道集进行分射线参数叠加，即得到用于反演的分射线参数叠加剖面。 

步骤二：分别对步骤一得到的纵波、转换波分射线参数叠加剖面进行稀疏约束反射系数序列反演，得到PP波、PS波分射线参数的反射系数剖面。 

稀疏约束反射系数序列反演的方法原理是： 

地震反射系数序列反演与地震数据处理中提高数据分辨率的反褶积紧密关联，而反褶积的最佳结果正是地层的反射系数序列。从统计学角度认为反射系数序列呈高斯分布，但是传统的反褶积结果通常较光滑，很难作为反射系数序列来进一步使用。实际地震反射系数序列通常是脉冲状的，它们比高斯分布更稀疏，而与柯西分布更接近。为此，本方法中以迭代反演代替反褶积处理，在反演中引入柯西准则，实现柯西稀疏约束的反射系数反演，提高解的稀疏性，得到的反射系数序列具有较好的高频带特性。 

柯西约束反射系数序列反演的目标函数： 

$J=\frac{1}{2}{(d-C_{w}r)}^T{C}_D^{-1}(d-C_{w}r)+N\ln \left(\mathrm{\π\λ}\right)+\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}\ln (1+\frac{r_i^2}{{\mathrm{\λ}}^2})---\left(1\right)$

式中，r：反射系数，d：地震数据，C_w：子波矩阵，C_D：地震数据协方 

差矩阵，λ：柯西分布参数 

令 $\∂J/\∂r=0,$ 得到的解为 

$r={(C_w^T{C}_D^{-1}{C}_w+\frac{2}{{\mathrm{\λ}}^2}P)}^{-1}{C}_w^T{C}_D^{-1}d---\left(2\right)$

这里P是一个斜对角阵，定义为 

$P=\mathrm{diag}\{{(\frac{r_0^2}{{\mathrm{\λ}}^2}+1)}^{-1},{\left(\frac{r_1^2}{{\mathrm{\λ}}^2}\right)}^{-1},\mathrm{\Λ},{\left(\frac{r_{N-1}^2}{{\mathrm{\λ}}^2}\right)}^{-1}\}$

解（2）式通过迭代方式得到，而柯西分布参数λ则控制了反射系数向量r的稀疏程度。 

步骤三：利用步骤二得到的PP波和PS波的分射线参数反射系数剖面进行联合反演，得到纵波速度、横波速度及密度3个参数。 

联合反演的方法：利用PP波、PS波射线参数域高阶反射系数的线性近似公式，采用广义线性反演算法同步估算纵波速度、横波速度及密度3个弹性参数。 

PP波射线参数域的高阶反射系数： 

$R_{\mathrm{pp}}\left(p\right)\≈R_f-A_1{p}^2+B_1{p}^4---\left(3\right)$

其中， 

$A_1=2\frac{\mathrm{\Δ\μ}}{\mathrm{\ρ}}-(1-R_f)q_{\mathrm{\α}}{q}_{\mathrm{\β}}{\left(\frac{\mathrm{\Δ\μ}}{\mathrm{\ρ}}\right)}^2$

$B_1=2R_f{\left(\frac{\mathrm{\Δ\μ}}{\mathrm{\ρ}}\right)}^2-2q_{\mathrm{\α}}{q}_{\mathrm{\β}}{\left(\frac{\mathrm{\Δ\μ}}{\mathrm{\ρ}}\right)}^3+(1-R_f)q_{\mathrm{\α}}^2{q}_{\mathrm{\β}}^2{\left(\frac{\mathrm{\Δ\μ}}{\mathrm{\ρ}}\right)}^4$

PP波射线参数域的线性近似反射系数： 

$R_{\mathrm{pp}}\left(p\right)\≈R_f-2\frac{\mathrm{\Δ\μ}}{\mathrm{\ρ}}{p}^2---\left(4\right)$

PS波射线参数域的高阶反射系数： 

$R_{\mathrm{ps}}\≈-\frac{{\mathrm{\α}}_1}{{\mathrm{\β}}_1}\frac{\mathrm{\Δ\μ}}{\mathrm{\ρ}}{q}_{\mathrm{\α}}p(1-A_2{p}^2+B_2{p}^4)---\left(5\right)$

$A_2=\frac{1}{q_{\mathrm{\α}}{q}_{\mathrm{\β}}}-2\frac{\mathrm{\Δ\μ}}{\mathrm{\ρ}}+q_{\mathrm{\α}}{q}_{\mathrm{\β}}{\left(\frac{\mathrm{\Δ\μ}}{\mathrm{\ρ}}\right)}^2$

$B_2=\frac{2}{q_{\mathrm{\α}}{q}_{\mathrm{\β}}}\frac{\mathrm{\Δ\μ}}{\mathrm{\ρ}}-{\left(\frac{\mathrm{\Δ\μ}}{\mathrm{\ρ}}\right)}^2-2q_{\mathrm{\α}}{q}_{\mathrm{\β}}{\left(\frac{\mathrm{\Δ\μ}}{\mathrm{\ρ}}\right)}^3+q_{\mathrm{\α}}^2{q}_{\mathrm{\β}}^2{\left(\frac{\mathrm{\Δ\μ}}{\mathrm{\β}}\right)}^4$

PS波射线参数域的线性近似反射系数： 

                          (6) 

$P^{R_{\mathrm{ps}}}\≈-\frac{{\mathrm{\α}}_1}{{\mathrm{\β}}_1}\frac{\mathrm{\Δ\μ}}{\mathrm{\ρ}}{q}_{\mathrm{\α}}p(1-\frac{p^2}{q_{\mathrm{\α}}{q}_{\mathrm{\β}}})$ 线参数，θ_i和分别代表了入射角和反射角；α_i,β_i和ρ_i分别代表了纵波速度、横波速度和密度，α、β和ρ代表了上述三个参数在相邻介质的平均值。纵波相关的反射系数R_f可表示为 $R_f=\frac{{\mathrm{\ρ}}_{i+1}{q}_{\mathrm{\αi}}-{\mathrm{\ρ}}_i{q}_{\mathrm{\α}(i+1)}}{{\mathrm{\ρ}}_{i+1}{q}_{\mathrm{\αi}}+{\mathrm{\ρ}}_i{q}_{\mathrm{\α}(i+1)}},$ 其中 $q_{\mathrm{\αi}}=\sqrt{1/{{\mathrm{\α}}_i}^2-p^2}$ 和 

是纵波垂直慢度，q_α和q_β分别代表了纵波和横波的平均垂直慢度。Δμ=ρ_i+1β_i+1 ²-ρ_iβ_i ²是剪切模量差。 

广义线性反演方法：通过用泰勒展开把非线性反问题化成线性反问题进行反演求解。将PP、PS波线性近似反射系数公式对弹性参数模型求一阶泰勒展开式，可得 

$R_{\mathrm{PP}/\mathrm{PS}}-R_{\mathrm{PP}/\mathrm{PS}}\left(m\right)\≈\frac{{\∂R}_{\mathrm{PP}/\mathrm{PS}}\left(m\right)}{\∂m}(m-m_0)---\left(7\right)$

其中R_PP/ps是由步骤二求得的实际数据的射线参数域反射系数值，R_PP/ps(m)是由PP、PS波线性近似反射系数公式求得的理论反射系数值，m₀=[α₀,β₀,ρ₀]是初始模型向量，m=[α,β,ρ]是待解的模型向量。 

通过迭代计算，不断求解(5)式并修改模型参数，通过一定的准则（本方法采用模型修改量的模达到极小）得到模型的逼近值，即最终的反演结果。 

步骤四：利用步骤三得到的纵波速度、横波速度及密度三个参数，计算纵横波速度比、泊松比等对岩性和流体敏感的参数，对储层及流体分布进行预测。 

射线追踪的理论基础是,在高频近似条件下,地震波场的主能量沿射线轨迹传播。由于地震波的振幅主要和地震波场的主能量有关，因此本发明利用地震振幅随射线参数变化的信息反演，相比现有技术中利用纵波、转换波地震振幅随入射角的变化信息反演，得到的地层弹性参数更加准确，能有效的提高反演精度。 

由于采用的是PP波、PS波共成像点道集数据，因此纵波和转换波来自于同一条入射路径，虽然反射路径彼此分离，但是来自于同一深度同一反射点。因此，射线参数域反演是同点物性联合反演。采用这种同点物性联合反演方法，各参数之间的联系更加紧密，因此能更准确地对纵波、转换波资料进行联合反演，得到更可靠的地层弹性参数，提高纵横波地震资料油气预测的精度。 

本发明实现了纵波、转换波射线参数域的联合反演，反演精度较高，为从多波地震信息中提取岩性和流体参数提供了一种新的解决途径。 

附图说明

图1罗家地区PP波偏移距域共成像点道集图。 

图2罗家地区PP波射线参数域共成像点道集图。 

图3罗家地区PS波偏移距域共成像点道集图。 

图4罗家地区PS波射线参数域共成像点道集图。 

图5罗家地区射线参数值为60的PP波地震剖面图。 

图6罗家地区射线参数值为60的PS波地震剖面图。 

图7罗家地区PP波射线参数值为60的反射系数剖面图。 

图8罗家地区PS波射线参数值为60的反射系数剖面图。 

图9罗家地区射线参数域联合反演得到的纵波速度剖面图。 

图10罗家地区射线参数域联合反演得到的横波速度剖面图。 

图11罗家地区射线参数域联合反演得到的密度剖面图。 

图12罗家地区射线参数域联合反演纵横波速度比平面图。 

具体实施方式

现在将在下文中详细描述本发明的示例性实施例；然而，本发明的示例性实施例可以以许多不同的形式来实施，且不应该限于这里阐述的示例。相反，提供这些示例使得本公开将是彻底并完整的，并将向本领域技术人员充分地传达本发明的范围。 

根据本发明的示例性实施例，可以预先对采集到的地震资料进行处理，然后将处理后的地震数据输入到例如计算机中。例如，可以以本领域公知的方法进行二维多分量或三维多分量野外勘探采集，以得到原始采集的多分量数据。然后，经过二维或三维纵波、转换波地震资料处理，以得到包括能够进行储层纵横波速度比反演的纵波、转换波叠前时间偏移道集数据的地震数据。然后，可以将得到的地震数据输入到计算机中，从而在后面的步骤中对这样的地震数据进行处理。可选择地，也可以直接将采集到的地震资料输入到计算机中并在计算机中进行相应地处理。如此，输入的地震数据可以包括与二维纵波和转换波相关的数据和与三维纵波和转换波相关的数据中的至少一种。 

本发明的方法在胜利油田罗家多波工区进行了实验，取得了较好的应用效果。本实验利用研究区三维纵波、转换波地震资料进行射线参数域叠前联合反演，并对目的层进行含油气性预测。 

首先是射线追踪及分射线参数剖面的求取。如图1到图4所示，采用射线追踪方法将纵波、转换波共成像点道集从偏移距域转换到射线参数域，并进行分射线参数叠加，其射线参数从0到140ms/km，间隔为20ms/km。图5、图6为得到的射线参数值为60ms/km的PP、PS波地震剖面。 

然后分别对纵波、转换波分射线参数叠加剖面进行稀疏约束反射系数序列反演，反射系数反演的主要目的是去除子波的影响，图7、图8为得到的PP波、PS波分射线参数的反射系数剖面，可以看到，稀疏约束反演结果反映了更为精细的地层特征。 

利用得到的PP波和PS波的分射线参数反射系数剖面，运用三参数同步反演方法完成联合反演。图9到图11为得到的纵波速度、横波速度及密度剖面图。 

利用联合反演得到的纵波速度、横波速度等参数，计算速度比、泊松比等对岩性和流体敏感的参数，对储层及流体分布进行预测。图12为罗家地区目的层的射线参数域联合反演纵横波速度比平面图，图中用不同的井符号代表实际井的含油气情况（其中黑色圆点代表油井，白色圆点代表水井，圆圈代表干井及无生物灰岩井），对比纵横波速度比参数可以看出，大部分钻遇油层的井位于纵横波速度比低值区（红、黄、绿色区），大部分钻遇水层和未钻遇该套储层的井位于纵横波速度比高值区（蓝色区）；纵横波速度比的分布规律性较好，低值区主要分布于罗37至罗9-1井一带，与钻井揭示的油气发育区对应。利用区内56口井的钻井、测井、试油、生产数据进行对比和统计分析，吻合率达到85.7%，因此，反演取得了良好的效果。 

Claims (4)

1.一种纵横波射线参数域联合反演方法，其特征在于：所述方法通过将纵波、转换波叠前道集资料转换到射线参数域，利用地震振幅随射线参数变化的信息反演得到地层弹性参数，进而对岩性和流体进行预测。 

2.如权利要求1所述的一种纵横波射线参数域联合反演方法，其特征在于所述方法包括以下步骤： 

步骤一：采用射线追踪法将PP波、PS波共成像点道集数据从偏移距域转换到射线参数域，并进行分射线参数叠加，得到分射线参数叠加剖面； 

步骤二：分别对步骤一得到的纵波、转换波分射线参数叠加剖面进行稀疏约束反射系数序列反演，得到PP波、PS波分射线参数的反射系数剖面； 

步骤三：利用步骤二得到的PP波和PS波的分射线参数反射系数剖面进行联合反演，得到纵波速度、横波速度及密度3个参数； 

步骤四：利用步骤三得到的纵波速度、横波速度及密度三个参数，进一步计算得到纵横波速度比、泊松比等对岩性和流体敏感的参数，对储层及流体分布进行预测。 

3.如权利要求2所述的一种纵横波射线参数域联合反演方法，其特征在于：步骤二中所述稀疏约束反射系数序列反演的方法是：以迭代反演代替反褶积处理，在反演中引入柯西准则，实现柯西稀疏约束的反射系数反演，提高解的稀疏性，得到的反射系数序列具有较好的高频带特性； 

柯西约束反射系数序列反演的目标函数： 

式中，r：反射系数，d：地震数据，C_w：子波矩阵，C_D：地震数据协方差矩阵，λ：柯西分布参数 

令得到的解为 

这里P是一个斜对角阵，定义为 

解（2）式通过迭代方式得到，而柯西分布参数λ则控制了反射系数向量r的稀疏程度。 

4.如权利要求2所述的一种纵横波射线参数域联合反演方法，其特征在于：步骤三中所述联合联合利用PP波、PS波射线参数域高阶反射系数的线性近似公式，采用广义线性反演算法同步估算纵波速度、横波速度及密度3个弹性参数； 

所述PP波射线参数域的高阶反射系数为： 

其中， 

所述PP波射线参数域的线性近似反射系数为： 

所述PS波射线参数域的高阶反射系数为： 

所述PS波射线参数域的线性近似反射系数为： 

P是地下第i层介质的射线参数，θ_i和分别代表了入射角和反射角；α_i,β_i和ρ_i分别代表了纵波速度、横波速度和密度，α、β和ρ代表了上述三个参数在相邻介质的平均值。纵波相关的反射系数R_f可表示为其中和 是纵波垂直慢度，q_α和q_β分别代表了纵波和横波的平均垂直慢度；Δμ=ρ_i+1β_i+1 ²-ρ_iβ_i ²是剪切模量差； 

所述广义线性反演方法为：通过用泰勒展开把非线性反问题化成线性反问题进行反演求解；将PP、PS波线性近似反射系数公式对弹性参数模型求一阶泰勒展开式，可得 

其中R_PP/ps是由步骤二求得的实际数据的射线参数域反射系数值，R_PP/ps(m)是由PP、PS波线性近似反射系数公式求得的理论反射系数值，m₀=[α₀,β₀,ρ₀]是初始模型向量，m=[α,β,ρ]是待解的模型向量； 

通过迭代计算，不断求解(5)式并修改模型参数，得到模型的逼近值，即最终的反演结果。