河流相储层是一类重要的含油气储集层，其精细描述的核心是有效识别出河道。但是此类储层纵向上常表现为砂体多期叠置，可识别厚度小于地震分辨率；横向上物性变化快，非均质性特征强，河道边缘特征不清晰。为了提高此类储层的识别精度，众多专家学者进行了深入研究，相干体^[1-2]、地震波形结构属性^[3]、甜点属性^[4]、谱分解和分频属性融合^[5]、概率神经网络^[6]、相层双控智能识别描述^[7]等方法技术得到广泛应用。

上述研究主要集中在河道构型、地震相整体特征、含油气分布预测等相对宏观的信息，对河道边缘特征的精细刻画还较少，在河道信号微弱、强屏蔽干扰或者横向不连续等情况下难以精准刻画储层。而图像边缘检测可通过捕捉图像局部特征的细微突变实现边缘的定位^[8]，其本质是通过突出图像的有效边缘信息、弱化非边缘信息来提取边缘细节。在地震解释领域，已有学者将此项技术应用于断层与裂缝预测中^[9-11]以及砂坝、砂体厚度的描述中^[12-14]。然而现有相关文献对河道边缘特征的数学意义以及算子处理的本质含义分析不够深入，还有待精细解剖。

本文从河道边缘几何特征与数学含义出发，从图像处理角度分析相干体边缘特征识别优势及存在的不足，以Sobel算子为例探究算子处理边缘识别的本质。对于图像灰度不均导致小河道识别困难的问题，提出基于直方图均衡化和模糊集理论的图像增强方法。上述技术方法均用三维河道理论模型和实际资料加以检验，应用效果显著，相关认识与结论对开展河流相储层精细描述有一定的借鉴作用。

1.   河道边缘特征的几何表征与三维建模

1.1   河道边缘几何特征与数学含义

河流相储层中的河道，其垂向厚度一般小于调谐厚度，其正演模型的左右两边缘振幅为零，中心位置振幅最大（图1（a））。若提取振幅切片，河道内有幅值，呈“实心”特征。对图1（a）截面计算一阶导数，导数从左到河道中心为正，边缘特征为“右负左正”或者“左正右负”。如果对一阶导数求模，左右两边界均为正值，边缘呈“左右双正”特征（图1（c））。二阶导数是一阶导数的导数，所以它的值与原始振幅相比，符号相反，但河道空间形态会更加锐化（图1（d））。

图  1  河道边缘几何特征与数学含义

Figure  1.  Geometric characteristics and mathematical implication of the channel edge

下载: 全尺寸图片 幻灯片

1.2   河流相储层三维建模

为便于对本文方法进行验证，构建200×300×100（线×道×时间）三维模型。模型有5个水平层，设计了4条河道，速度结构如图2（a）所示，河道1到河道4的砂岩速度分别为3600、3200、3150和3300 m/s，围岩速度为3100 m/s。从速度结构上看，各支河道反射系数由强到弱分别为河道1、河道4、河道2、河道3。河道1和河道2的埋深设计在47～57 ms之间，河道3和河道4在55～60 ms之间，河道1和河道4宽度为10道，河道3宽度为15道，河道2宽度为5道。正演模型采样间隔为1 ms，子波采用主频为30 Hz的雷克子波，加入最大振幅5％的随机噪声（图2（b））。

图  2  用于方法测试的三维河道模型

Figure  2.  Three-dimensional channel model for method testing

下载: 全尺寸图片 幻灯片

2.   基于相干属性的经典边缘识别技术

2.1   方法原理概述

对于地质异常体的边缘特征识别，常用的地震属性有均方根振幅、相干、曲率、甜点等，尤以相干属性为最常用^[15]。对于以特征值为基础的第3代相干，其表征公式常有以下几种（表1）^[2]，物理意义最明确、应用最广的是最大特征值相干。通过相干属性计算，河道特征比原始切片会更加清晰。

表  1  特征值相干不同表征公式

Table  1.  Different characterization formulas of eigenvalue coherence

时间	作者	表征公式	物理含义
1999	Gersztenkorn 和Marfurt	${C_{31}} = \displaystyle \dfrac{{{\lambda _1}}}{{\sum\limits_{j = 1}^J {{\lambda _j}} }}$	用最大特征值在所有特征值中的占比来表示相干
2000	Randen等	${C_{32} } = \displaystyle \dfrac{ {2{\lambda _2} } }{ { {\lambda _1} + {\lambda _3} } } - 1 = \dfrac{ { {\lambda _2} - {\lambda _3} - ( { {\lambda _1} - {\lambda _2} } )} }{ { {\lambda _1} + {\lambda _3} } }$	被称为“chaos”的相干属性
2002	Bakker	${C_{33} } = \displaystyle \dfrac{ {2{\lambda _2}( { {\lambda _2} - {\lambda _3} })} }{ {( { {\lambda _1} + {\lambda _2} } )( { {\lambda _2} + {\lambda _3} } )} }$	重点考虑第2特征值和第3特征值的差异
2007	Donias等	${C_{34} } = 1 - \displaystyle \dfrac{3}{2}\dfrac{ { {\lambda _2} + {\lambda _3} } }{ { {\lambda _1} + {\lambda _2} + {\lambda _3} } } = \dfrac{ { {\lambda _1} - {\lambda _2} + {\lambda _1} - {\lambda _3} } }{ {2( { {\lambda _1} + {\lambda _2} + {\lambda _3} } )} }$	“disorder”相干属性
2017	Wu	${C_{35}} = \displaystyle \dfrac{{{\lambda _1}{{ - }}{\lambda _2}}}{{{\lambda _1}}}$	利用第1特征值和第2特征值的差异

下载: 导出CSV 

| 显示表格

2.2   模型测试与效果分析

对图2（b）所示模型，计算C₃₁相干，提取时间为55 ms处的相干切片与原始切片进行比较（图3）。可以看到：①原始切片（图3（a））展示的是振幅信息，所以河道内是包含振幅信息的，河道显示是“实心”的；②相干（图3（b））突出的是异常信息，所以在相干体上展示的是河道边缘信息，河道是“空心”的；③无论是原始切片还是相干切片，显示的幅值范围都需要人工设置比较合理的范围，才能做到基本清晰（如图3（c）调整了色标幅值范围，河道边缘对比度有了一定的提高，但受背景噪声影响，即使精细调整色标范围仍难以精准刻画能量较弱的河道边缘）。

图  3  模型相干切片与原始振幅切片比较

Figure  3.  Comparison between model coherence and original amplitude slices

下载: 全尺寸图片 幻灯片

2.3   实际资料应用与存在问题分析

选取胜利油田罗家高密度资料作为实例数据，数据范围400线和600道，时间950～1050 ms。研究区曲流河、辫状河等河流相储层比较发育，河道相互叠置，横向变化大，砂体厚度小，河道砂体的精细描述是勘探开发的重点。

通过观察原始振幅切片（图4（a），T0=984 ms）可以看到：①研究区中间偏右侧南北向发育有一条宽度较大的河道；②北部1条河道虽然不是很宽，但振幅较强。研究区大量发育的是河道较小的曲流河，它们的相干值差异不大，大量值分布在白色的小值上（图4（b），显示范围0.45～1.0），此时需要人为设置显示范围才能较好地展示全区河道分布（图4（c），显示范围0.8～1.0）。且即使色标做了调整，对河道的的刻画效果仍然欠佳，个别振幅较小的河道的边缘特征不明显，进行基于图像处理的目标处理具有客观必要性。

图  4  实际资料相干切片与原始振幅切片比较

Figure  4.  Comparison between coherence and original amplitude slices of the actual data

下载: 全尺寸图片 幻灯片

3.   基于算子处理的边缘增强技术

3.1   方法基本原理

计算一个点的导数，容易出现异常值。利用多个点特殊的组合，进行数值或导数计算，可以得到具有某些特殊特征的图像，这就是算子处理的功效。在图像处理领域，常用算子包括Canny、Prewitt、Sobel、Roberts、Laplace、LoG算子等，其中Sobel算子在梯度运算中应用了局部平均的思想，使其具备较强的抗噪能力，检测结果不易受噪声干扰，在地震边缘检测中已有应用^[12-14]。

下面以Sobel算子为例，阐述其物理含义与计算公式。如图5，对于中心点（x, y），将第3列的和（局部平均）减去第1列的和，构建的是该点x方向的梯度G_x。求和运算可以均衡幅值，起到平滑的作用，计算时距离中心点最近的两个点（x+1, y）和（x − 1, y）给予2倍的权重；相减（差分）操作可以突出异常的边缘信息，同时起到去除相似背景的作用。类似地，可以得到y方向的梯度G_y。上述操作既可以用计算公式（1）实现，也可以核函数卷积实现，核函数设置见图5。

图  5  Sobel算子及核函数解析图

Figure  5.  Analytical diagram of the Sobel operator and kernel function

下载: 全尺寸图片 幻灯片

实际应用时既可以用x、y方向的梯度单独表示，使其具有方向性，也可以用1范数、2范数或$ \infty $范数（极大模）表示。

3.2   模型测试与效果分析

为检验本文方法有效性，对图6（a）所示切片，用Sobel算子进行处理，结果如图6（b）～图6（d）所示。梯度处理后河道两侧出现了不同极性的边缘，以能量较强的河道1为例，边缘右负左正、下负上正，这与前文中图1的理论分析一致。经G_x方向处理后，线方向的河道边缘特征得到较好展示；G_y方向处理后，道方向的河道边缘特征有所加强，但南北向的河道空间特征变得模糊。图6（d）是计算G_x和G_y模极大值结果，可以看到由于已计算了梯度的绝对值，河道出现类似于相干检测的边缘特征。总的来看，原始切片在采用Sobel算子做相关处理后相比于传统相干属性其边缘检测效果获得一定的增强，但对于能量较弱的河道边缘检测效果仍有待提高。

图  6  理论模型Sobel算子处理及效果比较

Figure  6.  Sobel operator processing and effect comparison of the theoretical model

下载: 全尺寸图片 幻灯片

3.3   实际资料应用与尚存在的问题

对图7（a）实际资料进行处理，结果如图7（b）～图7（d）所示。对比分析可知，不同方向Sobel算子对南北向和东西向河道识别能力不一样，这与模型讨论认识一致；相对方向梯度，使用模值表征更利于刻画细小河道。

图  7  实际资料Sobel算子处理及效果比较

Figure  7.  Sobel operator processing and effect comparison of the actual data

下载: 全尺寸图片 幻灯片

需要指出的是，以上只是对切片进行Sobel卷积操作，它没有构造导向和挖掘数据子体内部特征的功能，效果还比较受限，分辨率还有待提高。

4.   基于图像灰度处理的边缘增强技术

经相干体处理后的切片其数值分布范围很不均匀，且该经典方法与本文后续提出的基于算子处理的边缘增强技术对能量较弱地质异常体边缘显示效果仍然不佳。针对这样的问题，我们采用基于图像灰度处理的边缘增强技术来加以改善。

图像边缘的抽象定义为图像灰度发生空间突变的像素集合，以灰度值表征的图像在均匀区域之间的幅值突变被称为边缘。若将地震数据视为图像，地震振幅视为灰度，地震边缘对应的幅值突变区域往往是波阻抗界面产生的位置，它是地下构造或岩性信息的一种反映。

4.1   方法基本原理

4.1.1   灰度均衡图像处理技术

为解决相干处理后切片数值分布范围不均匀的问题，文章引入图像处理中的灰度直方图均衡化技术。首先统计灰度图像中元素具有的相同灰度值i出现的次数n_i，计算灰度值i出现的概率：

式中，n为图像像素总数，L₀为图像灰度矩阵中非重复灰度级的个数。再计算灰度值i的累计概率：

按以下映射关系对原始像素进行处理，其中$ {C_{\text{p}}}(i) $为当前元素灰度值的均衡化结果，Le为图像处理预期的总灰度级（取值256）：

round为四舍五入取整函数，用${C_{\text{p}}}(i)$替代原始灰度值为i的像素，得到最终均衡结果。

4.1.2   基于模糊集理论的边缘增强技术

为了进一步分离出相干等切片上的有效边缘信息，对均衡化后的图像做进一步模糊增强处理。方法核心思想是利用隶属度函数将空间域图像变换为模糊域内的模糊特征平面，在模糊特征平面上对模糊特征进行非线性变换，最后将其变换回空间域，得到增强处理的图像。

在模糊集理论中，将一幅灰度级为L的M×N图像u视为一个模糊集，通过隶属度函数将其由空间域转换到模糊特征域，并由模糊矩阵P进行表征：

式中，${\mu _i}_j\left( {0 \leq {\mu _i}_j \leq 1} \right)$表示图像中第（i，j）个像素点的灰阶，$ {u_i}_j $为特定灰度级的隶属度，一般这个特定灰度级取图像最大灰度$ {u_{\max }} $，$ {\mu _i}_j $组成模糊特征平面$ \left\{ {{\mu _i}_j} \right\} $。采用Pal-King算法中的隶属函数^[16]：

式中，F_d为倒数模糊因子，F_e为指数模糊因子，这里取F_d=128，F_e=2。对$ {\mu _i}_j $进行模糊增强，非线性增强变换函数为：

式中，$ {E_r} $表示对函数E的r次迭代运算，$ {\mu _c} $为渡越点，一般取0.5。利用$ {F^{ - 1}} $对$ {\mu _i}_j' $进行逆变换即可得到增强后的图像$ {l_{ij}} $：

4.2   模型测试与效果分析

4.2.1   灰度均衡直方图对比及分析

从前面分析已知，相干切片上能量较弱的河道2与河道3的边缘特征依然不够清晰（图8（a）），视觉分辨率不高。分析其原因主要是灰度值过于集中（图8（c），灰度值大多聚焦在白色的大值上），所以与背景值接近的河道不易识别。图8（b）为灰度均衡化后结果，此时直方图主灰度级已分散到6个以上（图8（d）），由此河道弱边缘特征得到较大增强，不过背景噪声也被同步放大，需要后续进一步处理。

图  8  理论模型灰度均衡化处理前后切片及直方图对比

Figure  8.  Comparison of slices and histograms before and after grayscale equalization of the theoretical model

下载: 全尺寸图片 幻灯片

4.2.2   模糊边缘增强及迭代次数评价

针对均衡后噪声同步放大的问题（图9（a）），用模糊增强方法做进一步处理，图9（b）～图9（e）为1～4次模糊增强结果。可以看到：①1次模糊增强处理后，背景噪声得到较大程度的压制，弱边缘得到进一步增强；②模糊次数增加到 4次，噪声确实已得到较好去除，但边缘信息损失也很大。总的来看，1～2次模糊增强效果较好。

图  9  理论模型图像模糊增强处理及效果对比

Figure  9.  Image fuzzy enhancement processing and effect comparison of the theoretical model

下载: 全尺寸图片 幻灯片

4.3   实际资料应用与效果评价

利用以上方法对实际资料相干切片，做直方图均衡化和模糊增强处理，结果如图10所示。可以看到：①原始相干切片的灰度范围较为集中，对比度不高，直方图均衡化后河道边缘特征变得比较清楚，由此也说明这样的图像目标处理十分必要；②与模型测试类似，模糊增强处理1～2次，能消除部分较强能量的背景信息，使河道边缘得到更好展示，有利于储层精细描述与油气预测。

图  10  实际资料图像处理前后相干切片对比

Figure  10.  Comparison of coherent slices before and after image processing of the actual data

下载: 全尺寸图片 幻灯片

5.   结语

对于频繁分叉改道、多期砂体叠置的复杂河流相储层，直接用振幅切片，识别能力有限。多道相干计算，利用河道的空间信息，河道的边缘特征得到了有效的提取。但由于背景噪声的影响，相干体上能量较弱的细小河道还是得不到很好成像。

利用图像处理中的灰度直方图均衡化技术，可以有效提高相干属性切片的对比度，凸显微小边缘特征。在此基础上进一步结合模糊增强，分离有效边缘和无效背景，可以进一步增强相干属性的边缘表征能力。

图像处理中，Canny、Prewitt、Sobel等算子处理，物理意义明确，操作简单、快捷，可用于河流相储层的边缘检测。