CN105572736A - 一种利用地震资料定量预测坡折带的方法 - Google Patents

Abstract

本发明涉及一种利用地震资料定量预测坡折带的方法，A、根据目地质、地震等背景资料，对发育坡折带的层系进行精细三维地震层位解释，得到目的层层位数据；B、求取出目的层的坡度与坡向值；C、根据坡度、坡向平面图上的分布特征初步识别坡折带分布范围；若坡折带边界清晰，则直接求取坡折带几何参数：坡折、坡脚、坡度、高差、水平宽度；若坡折带边界模糊(如断裂影响)，则进行D步骤；D、统计并求取已初步识别的坡折带范围内所有坡向值的平均值θ；E、以θ为中心，δ为步长在θ±ε范围内扫描求取最优视坡度φ；F、根据步骤E求取的最优视坡度，识别坡折带边界，并分别求取折带几何参数。本发明实现坡折带坡度平面刻画；提高坡度计算精度。

Description

一种利用地震资料定量预测坡折带的方法

技术领域

本发明涉及石油勘探技术领域，特别涉及一种利用地震资料定量预测坡折带的方法。

背景技术

坡折带原是地貌学概念，指地形坡度突变的地带，不论在沉积盆地中还是在剥蚀区，都可能发育坡折带。坡折带主要由以下几部分组成：坡折(slopebreak):指地形由缓变陡的转折点连线,为上凸的脊；陡坡(slope):指其上下的地形坡度都迅速变缓的陡坡地带；坡脚(slopetoe):指地形由陡变缓的转折点连线，为下凹的槽。林畅松等通过对半地堑盆地同沉积构造坡折带的研究，认为断陷湖盆中存在的构造坡折带制约着盆地可容空间的变化，控制了层序的发育、沉积体系域及砂体的分布；王英民等通过对准噶尔盆地大型坳陷湖盆坡折带类型的探讨，指出坡折带对层序、沉积和非构造圈闭的发育具有十分重要的控制作用；张善文等提出利用坡折带理论寻找隐蔽油气藏并在济阳坳陷陆相湖盆勘探中取得了显著成果。大量研究表明，坡折带制约着盆地可容纳空间的变化，控制着斜坡区沉积相带的分布和地层层序的发育，对砂体的沉积、横向尖灭和地层岩性圈闭的发育也具有关键控制作用。因此坡折带的识别和几何参数的刻画对斜坡区沉积相划分和地层岩性圈闭发育有利区及圈闭识别都具有重要指导意义。

现有的破折带识别方法概括起来可分为两类，一类是以剖面变化特征来识别，包括地震、地质、沉积相和层序地层等剖面的识别，主要表现为地层厚度的急剧增加以及地层倾角的突然变化；另一类是以坡折带平面变化特征来识别，包括地层厚度、沉积相展布、地层倾角等来识别，其平面变化特征与剖面变化相对应。但是现有方法的不足之处有两点，第一是坡折带的划分多采用人工划分方法，其划分方案影响因素较多，且为对坡折带的定性认识，坡折、坡脚识别不准确；第二是不能刻画坡折带内部及其上下坡度的平面变化。

发明内容

本发明解决上述技术问题所采用的技术方案是：一种利用地震资料定量预测坡折带的方法，包括以下步骤：

A、根据目地质、地震等背景资料，对发育坡折带的层系进行精细三维地震层位解释，得到目的层层位数据；

B、根据精细解释的地震层位数据和坡度计算模型求取出目的层的坡度与坡向值；

C、根据坡度、坡向平面图上的分布特征初步识别坡折带分布范围；若坡折带边界清晰，则直接求取坡折带几何参数：坡折、坡脚、坡度、高差、水平宽度；若坡折带边界模糊(如断裂影响)，则进行D步骤；

D、统计并求取已初步识别的坡折带范围内所有坡向值的平均值θ；

E、以θ为中心，δ为步长在θ±ε范围内扫描求取最优视坡度φ；

F、根据步骤E求取的最优视坡度，识别坡折带边界，并分别求取折带几何参数。

进一步的，所述步骤E中的最优视坡度φ为已初步识别的坡折带范围内∑(β_t-β_r)²最小的视坡度，其中β_t为真坡度，β_r为视坡度，原理是保持坡折带范围内真坡度减小的最少。

进一步的，所述步骤B中坡度与坡向值的求取公式如下：

$\begin{array}{c}{f}_x=\frac{1}{13968w}\{496(z_4+z_{24}-z_2-z_{22})+455(z_1+z_{21}-z_5-z_{25})+1960(z_{19}+z_9-z_7-z_{17})\\ +50(z_{15}-z_{11}+88(z_{14}-z_{12}\left)\right)+152(z_6+z_{16}-z_{10}-z_{20})\}\end{array}$

$\begin{array}{c}{f}_y=\frac{1}{13968w}\{496(z_6+z_{10}-z_{16}-z_{20})+455(z_{21}+z_{25}-z_1-z_5)+1960(z_7+z_9-z_{17}-z_{19})\\ +50(z_3-z_{23}+88(z_8-z_{18}\left)\right)+152(z_{22}+z_{24}-z_2-z_{24})\}\end{array}$

$\mathrm{\β}=\arctan \sqrt{f_x^2+f_y^2}$

其中β为坡度，α为坡向，fx为x方向的偏导数，fy为y方向的偏导数，z_i为5×5局部移动窗口内对应的属性值(如高程，地层深度等)。

本发明的有益效果：(1)实现坡折带从定性识别到定量刻画，对坡折带坡折、坡脚、坡度、水平宽度、高差等几何参数进行定量刻画；(2)实现坡折带坡度平面刻画，对坡折带斜坡内坡度、坡折之上、坡脚之下的坡度进行平面展布刻画；(3)提高坡度计算精度，减小断裂、局部构造变形等非坡折因素的影响，提高坡折带识别和参数刻画精度。

附图说明

图1是坡度、坡向示意图；

图2是5×5局部移动窗口及权值分布图；

图3是实施例中目的层坡折带地震剖面；

图4是实施例中目的层坡度图；

图5是实施例中目的层坡向图；

图6是实施例中最优视坡度图；

图7是实施例中坡折带与断裂系统叠合图；

图8是坡度带参数定量刻画流程图；

具体实施方式

本发明的一种利用地震资料定量预测坡折带的方法，包括以下步骤：A、根据目地质、地震等背景资料，对发育坡折带的层系进行精细三维地震层位解释，得到目的层层位数据；B、根据精细解释的地震层位数据和坡度计算模型求取出目的层的坡度与坡向值；C、根据坡度、坡向平面图上的分布特征初步识别坡折带分布范围；若坡折带边界清晰，则直接求取坡折带几何参数：坡折、坡脚、坡度、高差、水平宽度；若坡折带边界模糊(如断裂影响)，则进行D步骤；D、统计并求取已初步识别的坡折带范围内所有坡向值的平均值θ；E、以θ为中心，δ为步长在θ±ε范围内扫描求取最优视坡度φ；F、根据步骤E求取的最优视坡度，识别坡折带边界，并分别求取折带几何参数。

优选的实施方案，所述步骤E中的最优视坡度φ为已初步识别的坡折带范围内∑(β_t-β_r)²最小的视坡度，其中β_t为真坡度，β_r为视坡度，原理是保持坡折带范围内真坡度减小的最少。

优选的实施方案，所述步骤B中坡度与坡向值的求取公式如下：

$\begin{array}{c}{f}_x=\frac{1}{13968w}\{496(z_4+z_{24}-z_2-z_{22})+455(z_1+z_{21}-z_5-z_{25})+1960(z_{19}+z_9-z_7-z_{17})\\ +50(z_{15}-z_{11}+88(z_{14}-z_{12}\left)\right)+152(z_6+z_{16}-z_{10}-z_{20})\}\end{array}$

$\begin{array}{c}{f}_y=\frac{1}{13968w}\{496(z_6+z_{10}-z_{16}-z_{20})+455(z_{21}+z_{25}-z_1-z_5)+1960(z_7+z_9-z_{17}-z_{19})\\ +50(z_3-z_{23}+88(z_8-z_{18}\left)\right)+152(z_{22}+z_{24}-z_2-z_{24})\}\end{array}$

$\mathrm{\β}=\arctan \sqrt{f_x^2+f_y^2}$

其中β为坡度，α为坡向，fx为x方向的偏导数，fy为y方向的偏导数，z_i为5×5局部移动窗口内对应的属性值(如高程，地层深度等)。

上述步骤的流程为图8所示，且原理是利用了坡度计算模型；其中坡度模型中空间曲面的坡度和坡向是互相联系的两个点位函数，坡度反映曲面的倾斜程度；坡向是斜坡面对的方向，一般以北方向为起始方向，并按顺时针方向度量(如图1所示)，地表某点的坡度β、坡向α是地形曲面z＝f(x，y)分别沿地图坐标x、y方向求取偏导f_x、f_y的函数，关系式为：

$\mathrm{\β}=\arctan \sqrt{f_x^2+f_y^2}---\left(1\right)$

上述公式(1)和(2)中，对偏导数f_x、f_y的计算是准确求取坡度的关键。Evan(1980)首次提出局部曲面拟合法，即在规则网格的局部对称窗口(如3×3，5×5窗口)利用二次曲面函数拟合地表高程，通过函数系数与函数偏导数对应关系求取f_x、f_y；Zevenbergen和Thornev(1987)等用类似的方法在3×3局部窗口中提出利用不完全四次曲面拟合；Florinsky(2009)在5×5局部窗口提出利用三次曲面拟合，拟合函数：

$z=\frac{1}{6}{\mathrm{ax}}^3+\frac{1}{6}{\mathrm{dy}}^3+\frac{1}{2}{\mathrm{bx}}^{2}y+\frac{1}{2}{\mathrm{cxy}}^2+\frac{1}{2}{\mathrm{rx}}^2+\frac{1}{2}{\mathrm{ty}}^2+sxy+px+qy+u---\left(3\right)$

其中，函数的各项系数与函数偏导数存在以下关系：

$a=\frac{{\∂}^{3}z}{\∂x^3}d=\frac{{\∂}^{3}z}{\∂y^3}b=\frac{{\∂}^{3}z}{\∂x^2\∂y}c=\frac{{\∂}^{3}z}{\∂x\∂y^2}r=\frac{{\∂}^{2}z}{\∂x^2}t=\frac{{\∂}^{2}z}{\∂y^2}s=\frac{{\∂}^{2}z}{\∂x\∂y}p=\frac{\∂z}{\∂x}q=\frac{\∂z}{\∂y}---\left(4\right)$

在5×5局部窗口(图2)中，设定各点坐标分别为：(-2w,2w,z₁),(-w,2w,z₂),(0,2w,z₃),(w,2w,z₄),(2w,2w,z₅),(-2w,w,z₆)，(-w,w,z₇),(0,w,z₈),(w,w,z₉),(2w,w,z₁₀),(-2w,0,z₁₁),(-w,0,z₁₂),(0,0,z₁₃),(w,0,z₁₄),(2w,0,z₁₅),(-2w,-w,z₁₆),(-w,-w,z₁₇),(0,-w,z₁₈),(w,-w,z₁₉),(2w,-w,z₂₀),(-2w,-2w,z₂₁),(-w,-2w,z₂₂),(0,-2w,z₂₃),(w,-2w,z₂₄)和(2w,-2w,z₂₅)。将各个点坐标带入方程可求解各系数，由于方程个数大于未知数个数，因此采用最小二乘法求解各系数。根据最小二乘法原理，可以得到25个线性方程：

α＝Fβ(5)

其中，α为25个已知高程z的列矩阵：

β为包含10个未知数的列向量：

$\mathrm{\β}=\left|\right|\begin{array}{c}a\\ d\\ b\\ c\\ r\\ t\\ s\\ p\\ q\\ u\end{array}\left|\right|---\left(7\right)$

F为已知坐标带入方程后与未知数β、已知数α对应的25×10的常数矩阵：

考虑不同距离的格网点对格网中心点参数提取的不同影响，即权的概念，常采用如下形式的权函数：m_i为对应格网点的权值，d_i为该格网点离中心点的距离，n一般取1或者2。以n＝2为例，权矩阵M为25×25的方阵：

$M=\left|\right|\begin{array}{c}\frac{1}{8}000000000000000000000000\\ 0\frac{1}{5}00000000000000000000000\\ 00\frac{1}{4}0000000000000000000000\\ 000\frac{1}{5}000000000000000000000\\ 0000\frac{1}{8}00000000000000000000\\ .\\ .\\ .\\ 0000000000000000000000000\\ .\\ .\\ .\\ 00000000000000000000\frac{1}{8}0000\\ 000000000000000000000\frac{1}{5}000\\ 0000000000000000000000\frac{1}{4}00\\ 00000000000000000000000\frac{1}{5}0\\ 000000000000000000000000\frac{1}{8}\end{array}\left|\right|---\left(9\right)$

则(5)式可化为：

α＝F^TMβ(10)

可解得系数β：

β＝(F^TMF)^-1F^TMα(11)

以及不含权时，系数β：

β＝(F^TF)^-1F^Tα(12)

由于计算坡度只需要偏导数因此这里只给出f_x、f_y表达式，采用反距离平方权：

$\begin{array}{c}{f}_x=\frac{1}{13968w}\{496(z_4+z_{24}-z_2-z_{22})+455(z_1+z_{21}-z_5-z_{25})+1960(z_{19}+z_9-z_7-z_{17})\\ +50(z_{15}-z_{11}+88(z_{14}-z_{12}\left)\right)+152(z_6+z_{16}-z_{10}-z_{20})\}\end{array}---\left(13\right)$

$\begin{array}{c}{f}_y=\frac{1}{13968w}\{496(z_6+z_{10}-z_{16}-z_{20})+455(z_{21}+z_{25}-z_1-z_5)+1960(z_7+z_9-z_{17}-z_{19})\\ +50(z_3-z_{23}+88(z_8-z_{18}\left)\right)+152(z_{22}+z_{24}-z_2-z_{24})\}\end{array}---\left(14\right)$

上述优选的方案是，所述步骤D中的最优视坡度φ为最能反映坡折带边界的视坡度。在实际地质条件中，地层构造往往包含一定的方向性信息，即构造特征在一定范围内是沿着某个方向分布或延伸的，为了突出某一特定方向上的构造信息，在计算地层坡度时引入任意方向的视坡度β_r：

β_r＝βcos(α-θ)(15)

式中，θ为选定的特定角度，是与正北方向顺时针的夹角。定义最能反映地质工作者(或地震解释人员)感兴趣的构造或地质体的视坡度为该构造或体质体的最优视坡度φ，一般情况下，最优视坡度与该构造或地质体的坡向有关，对于单一倾向的构造或地质体，可通过下式确定：

$\mathrm{\Ω}=\underset{i=1}{\overset{n}{\Σ}}{{\mathrm{\β}}_i}^2{(1-\cos ({\mathrm{\α}}_i-\mathrm{\θ}\left)\right)}^2---\left(16\right)$

式中，β_i、α_i分别为该构造或地质体范围内每一点的真坡度和坡向，n为该范围内所有格网点；当Ω取得最小值时，此时的θ为最优视坡度φ。

上述坡度模型的精度分析：由于用实际资料进行精度分析没有明确的可对比对象或真值，因此采用一种与实际地层构造相似的数学曲面进行分析，该数学曲面的坡度真值可通过表达式直接求取，与利用模型求取的坡度对比可进行精度、误差分布等分析。该曲面函数表达式：

$\begin{array}{c}z=30\[1-{\left(\frac{x}{500}\right)}^2\]e^{-{\left(\frac{x}{500}\right)}^2-{\left(\frac{y}{500}+1\right)}^2}-100\[0.2\left(\frac{x}{500}\right){\left(\frac{x}{500}\right)}^3-{\left(\frac{y}{500}\right)}^5\]e^{-{\left(\frac{x}{500}\right)}^2-{\left(\frac{y}{500}\right)}^2}-\\ \frac{10}{3}{e}^{-{\left(\right(\frac{x}{500})+1)}^2-{\left(\frac{y}{500}\right)}^2}\end{array}---\left(16\right)$

将曲面按不同间隔(表1)离散并利用该算法求取坡度，精度指标采用中误差：

$\mathrm{\σ}=\sqrt{\frac{{\Σ}_{i=1}^n{({\mathrm{\β}}_{iF}-{\mathrm{\β}}_{iT})}^2}{n}}---\left(17\right)$

式中，β_iF为利用(13)、(14)式求取的坡度，β_iτ为曲面真坡度，n为样点总数。

表1不同格网间隔中误差

下面结合附图和实施例对发明做进一步详细的说明。

实施例中的研究区是一个形成于晚海西期并长期继承发育的古凸起，整体隆升不均衡；图3该研究区L880测线地震剖面，目的层为北西高、南东低以斜坡形式过渡大型不整合面，其上大部分地层逐层沿该不整合面上超尖灭，其下部分地层直接超覆于石炭系之上并伴随地层削蚀现象；受石炭系古凸起和深部大型逆掩断裂的控制，目的层发育挠曲坡折带。

(1)根据精细的地震解释层位求取目的层坡度和坡向，并根据坡度和坡向平面分布特征初步识别坡折带分布范围；

图4为目的层坡度，图5为目的层坡向；据图4所示，看出坡折带呈“S”型分布贯穿研究区，坡脚线十分清晰，但坡折带中部由于断裂发育，坡折线难以划分；坡折带的坡向分布与坡度分布有一定的相关性，易于区分其他构造形态，同时，坡度和坡向能不同程度的反应断裂的分布。

(2)统计坡折带坡向值平均值θ，θ＝89°。

(3)以θ为中心，δ为步长在θ±ε范围内扫描求取最优视坡度，其中δ＝3°，ε＝15°，得φ＝95°。

图6为最优视坡度，是根据(2)、(3)步求取最能反应坡折带坡折和坡脚的视坡度；图7是已识别的坡折带和断裂系统叠加图，对比图3可以看出，最优视坡度减小了坡折带附近断裂的坡度值，而其本身的坡度基本不变，这种坡度的减小与构造的坡向十分相关，与坡折带得坡向差异越大，坡度值减小越多。断裂坡度的减小使得坡折线变得十分清晰，提高了坡折带边界识别精度。同时，坡折带的分布与实际勘探井位十分吻合。

(4)根据以后上步骤的成果，求取坡折带几何参数：坡折、坡脚、坡度、高差、水平宽度。根据坡折带形态，将坡折带分为三段，分别为坡折带北段、坡折带中段和坡折带南段，各参数如下表2：

表2坡折带参数表

Claims (3)

1.一种利用地震资料定量预测坡折带的方法，其特征在于，包括以下步骤：

A、根据目地质、地震等背景资料，对发育坡折带的层系进行精细三维地震层位解释，得到目的层层位数据；

B、根据精细解释的地震层位数据和坡度计算模型求取出目的层的坡度与坡向值；

C、根据坡度、坡向平面图上的分布特征初步识别坡折带分布范围；若坡折带边界清晰，则直接求取坡折带几何参数：坡折、坡脚、坡度、高差、水平宽度；若坡折带边界模糊(如断裂影响)，则进行D步骤；

D、统计并求取已初步识别的坡折带范围内所有坡向值的平均值θ；

E、以θ为中心，δ为步长在θ±ε范围内扫描求取最优视坡度φ；

F、根据步骤E求取的最优视坡度，识别坡折带边界，并分别求取折带几何参数。

2.根据权利要求1所述的一种利用地震资料定量预测坡折带的方法，其特征在于，所述步骤E中的最优视坡度φ为已初步识别的坡折带范围内∑(β_t-β_r)²最小的视坡度，其中β_t为真坡度，β_r为视坡度，原理是保持坡折带范围内真坡度减小的最少。

3.根据权利要求1或2所述的一种利用地震资料定量预测坡折带的方法，其特征在于，所述步骤B中坡度与坡向值的求取公式如下：

$\begin{array}{c}{f}_x=\frac{1}{13968w}\{496(z_4+z_{24}-z_2-z_{22})+455(z_1+z_{21}-z_5-z_{25})+1960(z_{19}+z_9-z_7-z_{17})\\ +50\left(z_{15}-z_{11}+88\left(z_{14}-z_{12}\right)\right)+152\left(z_6+z_{16}-z_{10}-z_{20}\right)\}\end{array}$

$\begin{array}{c}{f}_y=\frac{1}{13968w}\{496(z_6+z_{10}-z_{16}-z_{20})+455(z_{21}+z_{25}-z_1-z_5)+1960(z_7+z_9-z_{17}-z_{19})\\ +50\left(z_3-z_{23}+88\left(z_8-z_{18}\right)\right)+152\left(z_{22}+z_{24}-z_2-z_{24}\right)\}\end{array}$

$\mathrm{\β}=\arctan \sqrt{f_x^2+f_y^2}$

其中β为坡度，α为坡向，fx为x方向的偏导数，fy为y方向的偏导数，z_i为5×5局部移动窗口内对应的属性值(如高程，地层深度等)。