CN104516021B - 一种同时提高解析式稳定性和精度的射线弹性参数反演方法 - Google Patents

Abstract

本发明预测和油气检测中同时提高解析式稳定性和精度的射线弹性参数反演方法，通过提取子波，用纵波速度、地层密度制作合成地震记录，确定地震剖面上地震反射层位；在地震反射层位的约束下，由测井数据分别建立三个射线弹性阻抗模型，滤波得到低频模型；由三个入射角的叠加剖面递推反演得到三个相对弹性阻抗剖面；在已知井的射线弹性阻抗约束下，将射线弹性阻抗低频模型和对应入射角的相对弹性阻抗剖面相加获得三个角度的射线弹性阻抗剖面；并由此反演得到纵横波速度比、纵波阻抗、横波阻抗，然后求出其他储层弹性参数。本发明略去射线弹性阻抗公式的高次项，提高了解析式稳定性，提高了入射角大于25°的弹性反演的精度。

Description

一种同时提高解析式稳定性和精度的射线弹性参数反演方法

技术领域

本发明属于油气田勘探、开发过程中的储层预测和油气检测技术，是一种同时提高解析式稳定性和精度的射线弹性参数反演方法。

背景技术

目前叠前弹性反演较多采用的是Connolly(1999)提出的基于Zoeppritz方程近似的、与角度有关的弹性阻抗反演方法(Elastic Inversion，简称EI反演)，该方法假设整个地震剖面不同深度上纵、横波速度比为常数，不符合实际；假设地震波在地层中按照常入射角度传播；随着纵波入射角度的增大，EI反演的误差逐渐增大；在实际应用中弹性阻抗(EI)的幅值随入射角度变化很大，用的归一化方法难以控制(Whitcombe，2002；马劲风，2003)，幅值的变化易导致岩性和流体识别的错误。王仰华(2003)和马劲风(2004)分别公开了一种射线路径弹性阻抗反演方法，将弹性阻抗表示为纵波阻抗和纵横波速度比的函数关系式，射线弹性阻抗不需要假设纵横波速度比为常数，也无需进行入射角归一化处理，但在求解弹性参数时，解析式不稳定，抗噪能力差，因此限制了射线弹性阻抗的实际应用。刘力辉等(2011)对射线弹性阻抗方法进行二项式展开，将射线弹性阻抗表达成纵波阻抗和横波阻抗的函数，通过误差推导和试验分析，在小入射角度(小于25°)的情况下，与射线弹性阻抗方法相比，刘力辉的近似射线弹性阻抗方法在稳定性和抗噪能力方面有所改善，但是在入射角大于25°时，与射线弹性阻抗误差大于5％，而且随着入射角的增大误差越来越大。而对于入射角大于25°的弹性反演，这种方法不能有效区分储层物性特征。上述射线弹性阻抗反演具有诸多优势，但都存在解析式不稳定的问题。

发明内容

本发明目的是提供一种提高解析式稳定性和抗噪能力，同时提高入射角大于25°的弹性反演的精度，满足有些地区需要入射角大于25°弹性反演的同时提高解析式稳定性和精度的射线弹性参数反演方法。

本发明通过以下步骤实现：

1)在探区采集地震资料，处理得到地震剖面；

2)在探区井中激发地震纵波和横波，得到纵波速度、横波速度数据；

3)在探区测井得到地层密度；

4)选择标准地震子波或者从井旁地震剖面提取子波，用纵波速度、地层密度合成地震记录；

所述的标准地震子波是雷克子波或克劳德子波。

5)用合成地震记录对步骤1)地震剖面进行层位标定，确定目的层段，并拾取地震反射层位；

6)采用下式计算随入射角变化的射线弹性阻抗合成记录模型，确定射线弹性阻抗随入射角变化的类型；

公式中：θ为入射角；α为纵波速度；β为横波速度；ρ为密度；m为调节系数；REI(θ)为入射角为θ的射线弹性阻抗；

所述的调节系数取2～6。

以上公式是用标准砂岩模型正演测试改进的射线弹性阻抗公式，式中略去了射线弹性阻抗公式的高次项，提高了解析式稳定性；与刘力辉等人提出的方法比较，提高了入射角大于25°的弹性反演的精度。

7)由步骤5)层位标定得到的时深关系，用均匀采样或者非均匀采样把测井数据由深度域转换为时间域；

8)在时间域采用下式由测井数据计算三个入射角对应的射线弹性阻抗；

公式中：θ为入射角；α为纵波速度；β为横波速度；ρ为密度；m为调节系数；REI(θ)为入射角为θ的射线弹性阻抗。

所述的三个入射角可以选择为5°、15°、30°。

所述的调节系数取2～6。

9)在地震反射层位的约束下，由三个入射角的射线弹性阻抗值分别建立三个射线弹性阻抗模型；

所述的建立三个射线弹性阻抗模型采用反距离内插建模算法、克里金地质建模算法。

10)用数字低通滤波器对三个射线弹性阻抗模型进行滤波，得到三个射线弹性阻抗低频模型，根据地震剖面的频谱中低频缺失的范围确定数字低通滤波器的低通频率；

所述的确定数字低通滤波器的低通频率选择5～8hz。

11)从地震叠前道集数据提取模型对应的三个入射角的角度道集叠加剖面，根据叠前道集数据信噪比的大小选择每个角度叠加的范围；

所述的角度叠加的范围提取入射角为5°的角度道集叠加时，可选择3°～7°的道集进行叠加。

所述的三个入射角可以选择为5°、15°、30°。

12)对三个入射角的角度道集叠加剖面进行时差、振幅的一致性校正，通过时差调整消除三个剖面上主要目的层段反射同相轴的时差，通过振幅匹配方法使三个剖面上主要目的层段的总能量一致，进一步消除球面扩散等因素的影响；

13)对三个入射角的角度道集叠加剖面采用脉冲反褶积算法压缩地震子波，获得反射系数剖面，在此基础上采用递推反演算法得到三个相对弹性阻抗剖面；

14)在已知井的射线弹性阻抗约束下，将射线弹性阻抗低频模型和对应入射角的相对弹性阻抗剖面相加；获得三个角度的射线弹性阻抗剖面；

15)在已知三个角度的射线弹性阻抗的情况下，由下式计算β/α横纵波速度比；由1/β/α计算纵横波速度比α/β：

公式中：θ1、θ2、θ3为入射角；α为纵波速度；β为横波速度；ρ为密度；m为调节系数；REI(θ1)、REI(θ2)、REI(θ3)分别为入射角θ1、θ2、θ3对应的射线弹性阻抗,σ1min为入射角θ1、θ2时的目标函数，σ2min为入射角θ2、θ3时的目标函数，使σ1min和σ2min同时达到最小值时的横纵波速度比β/α即为该方程组的解；

16)将横纵波速度比β/α带入下式中的任意一个公式，求解纵波阻抗αρ；(创新技术)

公式中：θ1、θ2、θ3为入射角；α为纵波速度；β为横波速度；ρ为密度；m为调节系数；REI(θ1)、REI(θ2)、REI(θ3)分别为入射角θ1、θ2、θ3对应的射线弹性阻抗。

17)将横纵波速度比β/α和纵波阻抗αρ带入下式求出横波阻抗βρ；

βρ＝αρ*βρ/αρ＝αρ*β/α

公式中：α为纵波速度；β为横波速度；ρ为密度；

18)由纵横波速度比、纵波阻抗、横纵波阻抗求出其他储层弹性参数，完成射线弹性参数反演。

其他储层弹性参数为泊松比、杨氏模量、剪切模量。

本发明略去了射线弹性阻抗公式的高次项，提高了解析式稳定性；同时，与刘力辉等人提出的方法比较，提高了入射角大于25°的弹性反演的精度。本方法对标准模型I(图9)，在50°入射角以下的误差小于5％；对标准模型II和(图10)和III(图11)，在任意入射角，与射线弹性阻抗方法的误差小于3％，同时与Zoeppritz方程的正演结果误差也比较小。而刘力辉等人提出的方法在25°入射角以下的误差小于5％，大于25°入射角时误差大于5％，而且随着入射角的增大越来越大。

附图说明

图1 S087017测线地震剖面；

图2 S116井测井曲线显示；

图3 S087017测线地震剖面层位标定及层位解释；

图4 S116井射线弹性阻抗合成记录模型；

图5a入射角5°射线弹性阻抗模型；

图5b入射角15°射线弹性阻抗模型；

图5c入射角30°射线弹性阻抗模型；

图6a入射角5°角度道集叠加剖面；

图6b入射角15°角度道集叠加剖面；

图6c入射角30°角度道集叠加剖面；

图7a入射角5°射线弹性阻抗剖面；

图7b入射角15°射线弹性阻抗剖面；

图7c入射角30°射线弹性阻抗剖面；

图8本发明射线弹性参数反演剖面；

图9标准模型I的正演结果比较；

图10标准模型II的正演结果比较；

图11标准模型III的正演结果比较。

具体实施方式

以下结合附图和实验例详细说明本发明。

1)在探区采集地震资料，处理得到地震剖面(见图1)；

2)在探区井中激发地震纵波和横波，得到纵波速度、横波速度数据(见图2)；

3)在探区测井得到地层密度(见图2)；

4)选择25hz雷克子波，用纵波速度、地层密度合成地震记录(见图3)；

5)用合成地震记录对步骤1)地震剖面进行层位标定，确定目的层段，并拾取地震反射层位Tp7、Tp8、Tc2(见图3)；

6)采用下式计算随入射角变化的射线弹性阻抗合成记录模型，确定目的层段(见图4中红虚线范围)射线弹性阻抗随入射角θ的增大而增大属于第三类砂岩类型；

公式中：θ为入射角；α为纵波速度；β为横波速度；ρ为密度；m为调节系数；REI(θ)为入射角为θ的射线弹性阻抗；

所述的调节系数取2～6。

7)由步骤5)层位标定得到的时深关系，用均匀采样或者非均匀采样把测井数据由深度域转换为时间域；表1中左侧框为深度域测井曲线数据，右侧框转换得到的时间域测井曲线数据。

表1测井曲线深时转换及重新采样数据表

8)在时间域采用下式由测井数据计算三个入射角对应的射线弹性阻抗；见表2，其中REI(5)、REI(5)、REI(5)为入射角为5°、15°、30°对应的射线弹性阻抗；

公式中：θ为入射角；α为纵波速度；β为横波速度；ρ为密度；m为调节系数；REI(θ)为入射角为θ的射线弹性阻抗。

所述的三个入射角可以选择为5°、15°、30°。

所述的调节系数取2～6。

表2弹性阻抗数据表

9)在地震反射层位的约束下，由三个入射角的射线弹性阻抗值分别建立三个射线弹性阻抗模型，图5a为入射角5°射线弹性阻抗模型；图5b入射角15°射线弹性阻抗模型；图5c入射角30°射线弹性阻抗模型。

所述的建立三个射线弹性阻抗模型采用反距离内插建模算法、克里金地质建模算法。

10)用数字低通滤波器对三个射线弹性阻抗模型进行滤波，得到三个射线弹性阻抗低频模型，根据地震剖面的频谱中低频缺失的范围确定数字低通滤波器的低通频率；所述的确定数字低通滤波器的低通频率选择5～8hz。

11)从地震叠前道集数据提取模型对应的三个入射角的角度道集叠加剖面，图6a为入射角5°射线弹性阻抗模型；图6b入射角15°射线弹性阻抗模型；图6c入射角30°射线弹性阻抗模型，角度叠加的范围选择：提取入射角为5°的角度道集叠加时，选择3°～7°的道集进行叠加；提取入射角为15°的角度道集叠加时，选择13°～17°的道集进行叠加；提取入射角为30°的角度道集叠加时，选择28°～32°的道集进行叠加。

12)对三个入射角的角度道集叠加剖面进行时差、振幅的一致性校正，通过时差调整消除三个剖面上主要目的层段反射同相轴的时差，通过振幅匹配方法使三个剖面上主要目的层段的总能量一致，进一步消除球面扩散等因素的影响；

13)对三个入射角的角度道集叠加剖面采用脉冲反褶积算法压缩地震子波，获得反射系数剖面，在此基础上采用递推反演算法得到三个相对弹性阻抗剖面；

14)在已知井的射线弹性阻抗约束下，将射线弹性阻抗低频模型和对应入射角的相对弹性阻抗剖面相加；获得三个角度的射线弹性阻抗剖面；

图7a为入射角5°射线弹性阻抗剖面；图7b入射角15°射线弹性阻抗剖面；图7c入射角30°射线弹性阻抗剖面，可以看出入射角为30°射线弹性阻抗剖面在S116井、T26井的弹性阻抗有明显差异；而试气结果S116井4.13万方/日、T26井无产量，所以在较大入射角弹性阻抗也可以识别储层含气性。

15)在已知三个角度的射线弹性阻抗的情况下，由下式计算β/α横纵波速度比；由1/β/α计算纵横波速度比α/β：

公式中：θ₁、θ₂、θ₃为入射角；α为纵波速度；β为横波速度；ρ为密度；REI(θ₁)、REI(θ₂)、REI(θ₃)分别为入射角θ₁、θ₂、θ₃对应的射线弹性阻抗，σ¹ _min为入射角θ₁、θ₂时的目标函数，σ² _min为入射角θ₂、θ₃时的目标函数，使σ¹ _min和σ² _min同时达到最小值时的横纵波速度比β/α即为该方程组的解；

16)将横纵波速度比β/α带入下式中的任意一个公式，求解纵波阻抗αρ；

公式中：θ₁、θ₂、θ₃为入射角；α为纵波速度；β为横波速度；ρ为密度；REI(θ₁)、REI(θ₂)、REI(θ₃)分别为入射角θ₁、θ₂、θ₃对应的射线弹性阻抗。

17)将横纵波速度比β/α和纵波阻抗αρ带入下式求出横波阻抗βρ；

βρ＝αρ*βρ/αρ＝αρ*β/α

公式中：α为纵波速度；β为横波速度；ρ为密度；

18)由纵横波速度比、纵波阻抗、横纵波阻抗求出其他储层弹性参数，完成射线弹性参数反演。

其他储层弹性参数为泊松比、杨氏模量、剪切模量。

图8为射线弹性阻抗参数反演得到的纵波阻抗、横波阻抗、纵横波速度比剖面，可以看出盒8目的层116井附近纵波阻抗低、T26井附近纵波阻抗高；横波阻抗在两口井附近变化不大；段纵横波速度比剖面在S116井附近为低值，T26井附近为高值，具有含气后纵波速度下降造成的纵横波速度比低特点，与试气结果S116井4.13万方/日、T26井无产量相吻合。

模型正演测试如下：

表3是标准砂岩模型物性参数表。

标准砂岩模型物性参数表表3

图9-11是正演结果比较。图中①为Zoeppritz方程的正演结果；②为马劲风的射线弹性阻抗方法的正演结果；③为本方法的正演结果，④为刘力辉等人的近似射线弹性阻抗方法的正演结果。

Claims (5)

1.一种同时提高解析式稳定性和精度的射线弹性参数反演方法，特点是通过以下步骤实现：

1)在探区采集地震资料，处理得到地震剖面；

2)在探区井中激发地震纵波和横波，得到纵波速度、横波速度数据；

3)在探区测井得到地层密度；

4)选择标准地震子波或者从井旁地震剖面提取子波，用纵波速度、地层密度合成地震记录；

5)用合成地震记录对步骤1)地震剖面进行层位标定，确定目的层段，并拾取地震反射层位；

6)采用下式计算随入射角变化的射线弹性阻抗合成记录模型，确定射线弹性阻抗随入射角变化的类型；

$REI\left(\mathrm{\θ}\right)=\frac{\mathrm{\α}\mathrm{\ρ}}{\cos \mathrm{\θ}}(1-4\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}{\sin }^2\mathrm{\θ}+m\frac{{\mathrm{\β}}^4}{{\mathrm{\α}}^4}{\sin }^4\mathrm{\θ});$

公式中：θ为入射角；α为纵波速度；β为横波速度；ρ为密度；m为调节系数；REI(θ)为入射角为θ的射线弹性阻抗；调节系数m取2～6；

7)由步骤5)层位标定得到的时深关系，用均匀采样或者非均匀采样把测井数据由深度域转换为时间域；

8)在时间域采用下式由测井数据计算三个入射角对应的射线弹性阻抗；

$REI\left(\mathrm{\θ}\right)=\frac{\mathrm{\α}\mathrm{\ρ}}{\cos \mathrm{\θ}}(1-4\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}{\sin }^2\mathrm{\θ}+m\frac{{\mathrm{\β}}^4}{{\mathrm{\α}}^4}{\sin }^4\mathrm{\θ});$

公式中：θ为入射角；α为纵波速度；β为横波速度；ρ为密度；m为调节系数；REI(θ)为入射角为θ的射线弹性阻抗；

9)在地震反射层位的约束下，由三个入射角的射线弹性阻抗值分别建立三个射线弹性阻抗模型；

10)用数字低通滤波器对三个射线弹性阻抗模型进行滤波，得到三个射线弹性阻抗低频模型，根据地震剖面的频谱中低频缺失的范围确定数字低通滤波器的低通频率；

11)从地震叠前道集数据提取模型对应的三个入射角的角度道集叠加剖面，根据叠前道集数据信噪比的大小选择每个角度叠加的范围；

12)对三个入射角的角度道集叠加剖面进行时差、振幅的一致性校正，通过时差调整消除三个剖面上主要目的层段反射同相轴的时差，通过振幅匹配方法使三个剖面上主要目的层段的总能量一致，进一步消除球面扩散因素的影响；

13)对三个入射角的角度道集叠加剖面采用脉冲反褶积算法压缩地震子波，获得反射系数剖面，在此基础上采用递推反演算法得到三个相对弹性阻抗剖面；

14)在已知井的射线弹性阻抗约束下，将射线弹性阻抗低频模型和对应入射角的相对弹性阻抗剖面相加；获得三个角度的射线弹性阻抗剖面；

15)在已知三个角度的射线弹性阻抗的情况下，由下式计算β/α横纵波速度比；由1/β/α计算纵横波速度比α/β：

$\left\{\begin{array}{c}{\mathrm{\σ}}_{\min }^1={(\frac{REI\left({\mathrm{\θ}}_1\right)}{REI\left({\mathrm{\θ}}_2\right)}-\frac{{\mathrm{cos\θ}}_2}{{\mathrm{cos\θ}}_1}\frac{\[1-4\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}{\sin }^2{\mathrm{\θ}}_1+m\frac{{\mathrm{\β}}^4}{{\mathrm{\α}}^4}{\sin }^4{\mathrm{\θ}}_1\mathrm{\]}}{\[1-4\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}{\sin }^2{\mathrm{\θ}}_2+m\frac{{\mathrm{\β}}^4}{{\mathrm{\α}}^4}{\sin }^4{\mathrm{\θ}}_2\mathrm{\]}})}^2\\ {\mathrm{\σ}}_{\min }^2={(\frac{REI\left({\mathrm{\θ}}_3\right)}{REI\left({\mathrm{\θ}}_2\right)}-\frac{{\mathrm{cos\θ}}_2}{{\mathrm{cos\θ}}_3}\frac{\[1-4\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}{\sin }^2{\mathrm{\θ}}_3+m\frac{{\mathrm{\β}}^4}{{\mathrm{\α}}^4}{\sin }^4{\mathrm{\θ}}_3\mathrm{\]}}{\[1-4\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}{\sin }^2{\mathrm{\θ}}_2+m\frac{{\mathrm{\β}}^4}{{\mathrm{\α}}^4}{\sin }^4{\mathrm{\θ}}_2\mathrm{\]}})}^2\end{array}\right.$

公式中：θ₁、θ₂、θ₃为入射角；α为纵波速度；β为横波速度；ρ为密度；m为调节系数；REI(θ₁)、REI(θ₂)、REI(θ₃)分别为入射角θ₁、θ₂、θ₃对应的射线弹性阻抗，σ¹ _min为入射角θ₁、θ₂时的目标函数，σ² _min为入射角θ₂、θ₃时的目标函数，使σ¹ _min和σ² _min同时达到最小值时的横纵波速度比β/α即为该方程组的解；

16)将横纵波速度比β/α带入下式中的任意一个公式，求解纵波阻抗αρ；

$\left\{\begin{array}{c}REI\left({\mathrm{\θ}}_1\right)=\frac{\mathrm{\ρ}\mathrm{\α}}{{\mathrm{cos\θ}}_1}\[1-4\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}si{n}^2{\mathrm{\θ}}_1+m\frac{{\mathrm{\β}}^4}{{\mathrm{\α}}^4}si{n}^4{\mathrm{\θ}}_1\]\\ PEI\left({\mathrm{\θ}}_2\right)=\frac{\mathrm{\ρ}\mathrm{\α}}{{\mathrm{cos\θ}}_2}\[1-4\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}{\sin }^2{\mathrm{\θ}}_2+m\frac{{\mathrm{\β}}^4}{{\mathrm{\α}}^4}{\sin }^4{\mathrm{\θ}}_2\]\\ REI\left({\mathrm{\θ}}_3\right)=\frac{\mathrm{\ρ}\mathrm{\α}}{{\mathrm{cos\θ}}_3}\[1-4\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}{\sin }^2{\mathrm{\θ}}_3+m\frac{{\mathrm{\β}}^4}{{\mathrm{\α}}^4}si{n}^4{\mathrm{\θ}}_3\]\end{array}\right.$

公式中：θ₁、θ₂、θ₃为入射角；α为纵波速度；β为横波速度；ρ为密度；REI(θ₁)、REI(θ₂)、REI(θ₃)分别为入射角θ₁、θ₂、θ₃对应的射线弹性阻抗；

17)将横纵波速度比β/α和纵波阻抗αρ带入下式求出横波阻抗βρ；

βρ＝αρ*βρ/αρ＝αρ*β/α

公式中：α为纵波速度；β为横波速度；ρ为密度；

18)由纵横波速度比、纵波阻抗、横纵波阻抗求出泊松比、杨氏模量和剪切模量，完成射线弹性参数反演。

2.根据权利要求1的方法，特点是步骤8)所述的三个入射角为5°和15°和30°；调节系数取2～6。

3.根据权利要求1的方法，特点是步骤10)所述的确定数字低通滤波器的低通频率选择5～8hz。

4.根据权利要求1的方法，特点是步骤11)所述的角度叠加的范围提取入射角为5°的角度道集叠加时，可选择3°～7°的道集进行叠加。

5.根据权利要求1的方法，特点是步骤11)所述的三个入射角为5°、15°、30°。