基于多属性的断溶体特征提取与分带性自动识别

郑笑雪,唐金良

(中国石油化工股份有限公司石油物探技术研究院,江苏南京211103)

鲁新便等[1]在塔河油田2013年的勘探开发中提出了“断溶体”的概念。致密碳酸盐岩或低孔、低渗的碳酸盐岩地层受多期次构造挤压作用,沿断裂带形成了不同规模的破碎带,多期岩溶水沿断裂下渗致使破碎带内断裂、裂缝被溶蚀改造,从而形成板状溶蚀孔、洞储集体,该储集体在上覆泥灰岩、泥岩等盖层封堵以及侧向致密灰岩遮挡下构成的圈闭类型,称为“断溶体圈闭”。这类圈闭后期因油气沿断裂运移、充注后成藏,形成“断溶体油(气)藏”。目前已探明的多个顺北油气富集区均沿大型断裂带分布,断裂破碎带既是油气疏导通道,又是成藏的有利空间。顺北特深断溶体埋深普遍超7000m,构造应力差异导致断裂发育样式多,岩溶作用导致储层非均质性强,因其地震波场特征复杂,故难以通过简单的地球物理方法精细描述和刻画断溶体储层特征,因而制约了断溶体油藏高效勘探开发。因此,有必要对断溶体进行针对性的地球物理方法技术和流程研究。

杨德彬[2]、焦方正[3]对顺北地区特深断溶体的油气富集规律进行了综合分析,所总结的地震反射特征类型及规律对断溶体油藏后续处理解释具有一定的借鉴意义;马乃拜等[4]研究了断溶体的地震反射特征,建立的模型有效指导了断溶体的识别及开发;徐红霞等[5]利用多属性融合技术预测断溶体的连通性;李鹏飞等[6]利用Geoeast解释软件对断溶体进行了描述与刻画。上述技术均取得了良好的应用效果。

顺北奥陶系上部地层的强地震反射同相轴削弱了下部断溶体的地震响应特征,为突出断溶体的地震异常特征,使其便于识别和描述,需采用相应的信号处理技术消除或减弱反映地层结构特征的连续强反射同相轴[7]。1998年,HUANG等[8]提出了经验模态分解(empirical mode decomposition,EMD)算法,将信号分解为多个不同频率特性的分量和一个残余分量;张义平等[9]将其用于地震信号处理,陈东方等[10]利用经验模态分解算法压制信号中的噪声干扰;2003年,二维经验模式分解(bidimensional empirical mode decomposition,BEMD)[11-13]概念被提出并广泛应用于各领域;DAMERVAL等[14]采用三角剖分和分段三次多项式插值改进了EMD算法;宋平舰等[15]采用BEMD算法对遥感图像进行了处理,并取得了一定效果;2005年,刘忠轩等[16]将经验模态分解算法与纹理分析相结合;姜晓宇等[17]利用基于相干技术与广义S变换频谱分解的不连续性检测技术对断溶体油藏进行了检测。采用BEMD信号处理技术和分频滤波等相关处理手段对断溶体地震响应特征进行增强预处理,是实现断溶体精细描述和刻画的前提。

利用地震信息识别和描述断溶体是顺北断溶体油藏勘探和开发的关键环节之一。地震属性技术广泛应用于油气藏表征与储层预测等领域,其中纹理属性和梯度结构张量属性在计算机视觉和图像分类等领域已经取得了较好的应用效果,为将其应用领域拓展至地震数据处理,许多学者进行了研究。GAO[18-20]通过计算纹理属性实现了断层识别;KNUTSSON等[21]提取了地震剖面的结构张量属性;RAO等[22]、FARAKLIOTI等[23]、RANDEN等[24]提取了各向异性介质中结构张量属性并进行了研究;张军华等[25]和崔世凌等[26]分别将纹理属性和结构张量属性提取算法应用于不同工区完成了断层、河道等特殊异常体的识别;龚屹等[27]利用纹理属性的聚类分析算法进行裂缝预测;赵淑琴等[28]将纹理属性提取技术应用于火山岩边界刻画。

综合利用多种属性进行聚类分析,得到构造、岩性和储层等多种信息,可有效降低利用属性进行储层识别的多解性。SELIM等[29]发展了K-means方法并得到了广泛应用;MATOS等[30]利用属性进行聚类分析,并划分了地震相;杨培杰等[31]利用模糊C均值算法进行地震属性聚类分析;蔡涵鹏等[32]利用纹理特征进行了半监督地震相分析;王治国等[33]对地震多属性进行了优化分析;聚类算法在石油工业中展现了广阔的应用前景,具有重要的理论和实用价值[34]。

本文提出的基于多属性的断溶体特征提取与分带性自动识别方法,首先以BEMD信号分解理论为基础,削弱反映地层结构的强反射同相轴,突出断溶体异常响应;然后利用优化梯度结构张量和地震纹理属性算法,提高断溶体储层地震属性识别精度;再利用优选的多种断溶体识别地震属性,采用K-means算法进行降维化聚类;最终实现断溶体空间形态和分带性的自动刻画。

1 断溶体边界刻画技术

1.1 基于BEMD信号处理的断溶体特征增强分析

采用BEMD[11-13]算法对断溶体剖面数据进行预处理,假定地震数据满足相应条件[11-13],具体分析过程如下。

1) 定义初始剖面参数,将叠后地震剖面数据进行灰度化处理后,作为输入数据:

r1,1(x,y)=f(x,y)

(1)

式中:f(x,y)为地震剖面;r1,1(x,y)为初始输入剖面。

2) 计算输入剖面rp,k(x,y)的极大值和极小值,其中p代表从地震剖面计算得到的第p个分解分量,k为筛分次数(P=1,2,,p;K=1,2,,k)。

(2)

4) 将断溶体剖面rp,k(x,y)减去计算的曲面均值mp,k(x,y),得到:

hp,k(x,y)=rp,k(x,y)-mp,k(x,y)

(3)

5) 重复上述步骤,当hp,k(x,y)不能再被分解时,即满足约束条件SD:

(4)

分解得到的第p个IMF分量为Cp(x,y),Cp(x,y)=hp,k(x,y)。

6) 计算残余量:

rp+1,k(x,y)=rp,1(x,y)-Cp(x,y)

(5)

7) 经过p次分解计算得到残余量rp(x,y)=rp+1,k(x,y),目标信号为:

(6)

对如图1a所示的地质模型进行BEMD分解测试。该模型通过在断裂带内添加不同大小和发育密度的缝洞体来模拟断溶体的分带性,其中断溶体核部为大尺度缝洞体且缝洞体发育密集,翼部缝洞体尺度相对较小且发育密度低,断溶体上覆地层为泥岩层。从地震剖面可以看出,断溶体发育的碳酸盐岩地层上部存在强能量地震同相轴,影响了断溶体的特征识别,此外,断溶体弱破碎带(相当于断溶体翼部)地震响应呈现弱反射特征,当破碎程度较强(相当于断溶体核部)时,地震响应表现为杂乱异常反射特征(图1b)。BEMD分解结果(图1c)在一定程度上使异常更聚焦,同时对背景数据进行了压制。对模拟地震数据进行处理,通过BEMD分解得到单频或多频段内的不同组合剖面,分解得到的剖面分别代表IMF1,IMF2,IMF3分量。分析其有效频段数据组合方式获取有效的BEMD提取结果(图2),IMF1分量主要反映断溶体空间异常特征,IMF2分量主要反映地层结构信息,IMF3分量代表算法分解的残差部分,IMF1分量增强了断溶体地震响应特征信息,为实际工区效果分析奠定了基础。

图1 模型的BEMD分解a 地质模型; b 地震剖面; c BEMD分解结果

图2 模型BEMD提取结果a 地震剖面; b IMF1分量; c IMF2分量; d IMF3分量

1.2 基于地震纹理分析的断溶体识别研究

对剖面进行BEMD处理后,通过计算纹理属性和结构张量属性,探讨地震纹理分析在断溶体识别和预测中的应用。

1.2.1 基于地震纹理属性的断溶体识别

HARALICK等[34]利用灰度共生矩阵计算纹理属性,并将其应用于地球物理领域。首先对地震剖面振幅数据进行灰度化计算,然后统计一定距离和方向上地震数据的相似性,得到地震振幅数据在不同方向、不同角度等方面的变化特征。进行灰度化计算时,沿着某一方向计算距离为δ的任意两个剖面振幅满足一定条件的概率为:

p(i,j)=∑{g(x,y)=i,g(x+Δx,y+Δy)=j}i,j=0,1,,L-1,L

(7)

1) 能量。

能量统计了纹理灰度的变化程度,反映地震剖面振幅的分布均匀程度。

(8)

2) 熵。

熵可以度量剖面中纹理的随机性。当共生矩阵的值大小都相等时,熵为最大值;反之,其值较小。

(9)

3) 对比度。

对比度为灰度共生矩阵主对角线附近的惯性矩,可以反映剖面矩阵的清晰度和纹理的沟纹深浅。

(10)

选取如图1a所示模型进行纹理属性提取分析,为了更好地识别断溶体发育特征,提取不同方向、不同灰度级、不同窗口的纹理属性,结果如图3所示。其中对比度属性、熵属性、能量属性均可表征剖面上的断溶体,优化后的纹理属性连续性得到了增强,其中熵属性更为连贯,说明断溶体灰度值变化较大。

图3 纹理属性提取结果a 某商业软件的对比度属性; b 某商业软件的能量属性; c 某商业软件的熵属性; d 采用本文方法得到的对比度属性; e 采用本文方法得到的能量属性;f采用本文方法得到的熵属性

1.2.2 基于地震梯度结构张量属性的断溶体识别研究

梯度结构张量(GST)主要用于计算地震构造属性和预测裂缝发育区,本文在现有算法的基础上研发优化了GST属性的提取方法,具体包括以下步骤。

1) 在计算梯度结构张量之前利用BEMD进行断溶体特征增强处理,然后采用特定滤波器算子对剖面进行处理,增强断溶体边界。

2) 计算三维地震数据体每一点的梯度矢量:

(11)

(12)

式中:u(x,y,z)为地震数据。

3) 构建梯度结构张量T[33]:

(13)

4) 采用不同的高斯核参数对梯度结构张量的每一个成分进行平滑测试,选取最优结果:

(14)

式中:G(x,σT)是高斯核函数;σT是尺度参数。

(15)

5) 计算梯度结构张量的特征值和特征向量:

|Tυ-λυ|=0

(16)

式中:λ为高斯核函数的相关参数。利用不同高斯核函数的相关参数对模型(图1a)提取梯度结构张量属性,结果如图4所示,从左至右所选用的参数λ分别

图4 模型梯度结构张量属性提取结果a λ1分量; b λ2分量; c λ3分量

为0.1,0.5,2.0,2.5。不同高斯核参数下得到的结构张量属性λ1,λ2差异明显,λ3差异较小;λ1提取结果相对较差,整体信息虽然较丰富,但是背景噪声较严重;λ2相对较好,但是仍存在部分模糊现象,整体提取效果较理想且不受参数影响。对比剖面可知,当参数λ为0.1和2.5时,属性计算效果存在模糊性(图4a);参数λ为0.5时,背景存在较多噪声;参数λ为2时的效果较为理想。

1.3 基于多属性聚类的断溶体分带性自动识别研究

属性解释的多解性造成储层的准确识别存在困难。属性提取结果是同一地质体的不同方面的响应,地震属性聚类分析是对相互独立参数进行分析,因此可以在无井约束的情况下完成储层分类预测。本文利用K-means原理将多个地震数据进行聚类,在实现属性降维的同时,实现断溶体分带性自动识别。本文将属性集分为c个类别,算法描述如下。

1) 适当选择c个类的中心作为初始剖面聚类中心。

2) 在第k次迭代中,对任意一个属性数据,求其到c个中心剖面的距离,将该样本归到距离最短的中心所在的类:

(17)

式中:xi代表属性剖面;μj代表聚类中心剖面。

3) 利用均值等方法更新该类的中心值:

(18)

4) 对于所有的c个聚类中心,如果利用步骤2)和3)迭代更新后,中心值保持不变,则迭代结束,否则继续迭代。

采用上述方法计算相关属性,模型优选的多种属性如图5a所示,从左到右分别是为能量性、熵属性、以断溶体地震响应特征增强预处理方法为基础,开展以纹理属性和GST属性为核心的多属性断溶体识别分析,结合空间描述和分带性识别技术,提出了针对断溶体的一体化处理识别方法,具体流程见图6。

图5 模型属性聚类结果a 模型优选多种属性; b 属性分带性聚类结果; c 聚类结果与原剖面叠合显示

图6 基于多属性的断溶体特征提取与分带性自动识别算法流程

2 实际资料处理

对顺北油田B工区三维地震资料进行处理,结果如图7所示,顺北油田发育一系列沿断裂分布的碳酸盐岩断溶体油藏,其中典型断溶体地震响应特征如图7中红框所示。对实际地震资料进行BEMD预处理,结果如图8所示。对预处理后的实际地震资料提取纹理属性(图9a)和梯度结构张量属性(图9b),可见利用优化算法的纹理属性所识别断溶体形态更加连续,梯度结构张量属性在一定程度上可以指示断溶体边界轮廓。将不同过井曲线的纹理属性、结构张量属性等作为输入数据进行断溶体分带性聚类分析,结果如图10所示,断溶体分带性与钻井所揭示的储层信息特征相吻合。图11为多属性聚类断溶体分带性自动识别的时间切片和三维可视化结果,断溶体分带性刻画精度高,与井资料吻合程度高,证明了方法的有效性。

图7 实际地震资料剖面(含断溶体油藏)

图8 实际地震资料BEMD预处理结果a IMF1分量; b IMF2分量; c IMF3分量

图9 对预处理后的实际地震资料提取的纹理属性(a)和梯度结构张量属性(b)

图10 断溶体分带性聚类分析结果a 过井1-5曲线; b 过井1-1曲线

图11 多属性聚类断溶体分带性聚类结果时间切片(a)和三维可视化结果(b)

3 结论

1) 断溶体的检测很大程度上依赖于地震剖面处理结果的质量。基于BEMD的地震信号分解技术有效削弱了与地层结构相关的地震强反射同相轴,增强了断溶体的地震响应特征。

2) 通过优化纹理属性和GST属性算法,有效强化了纹理属性识别断溶体的能力,提高了属性识别断溶体的精度,获得了突出的断溶体识别应用效果。综合利用多种地震属性,充分挖掘单一属性内涵信息,利用K值聚类算法可以实现属性的降维,降低多属性预测储层的多解性,进一步提高储层预测精度。

3) 地震属性处理结果,关系到断溶体边界及内幕缝洞储层刻画精度,随着越来越多智能化优化算法的出现以及混合算法的改进,智能化自动识别断溶体描述和刻画算法具有广阔的应用前景,将其应用于断溶体边界识别并结合内幕储层特征研究是下一步研究方向。