基于限定对比度直方图均衡的超声测井图像增强方法

付青青, 杨立华, 付建栋, 吴爱平

(1.长江大学电信学院, 湖北 荆州 434023; 2.中国石油集团测井有限公司, 陕西 西安 710201)

0 引 言

超声成像测井不仅可以在裸眼井中反映井眼几何形状,识别裂缝、孔洞、层理等地层非均质性,还能在套管井中检查射孔质量、分析套管损坏等。测井过程中,由于超声换能器的形状和工作频率、井孔的大小和形状、井壁的倾角和表面结构、井液泥浆密度等因素会影响超声波的衰减[1-2],因此,一幅超声成像测井图像的灰度分布通常会集中在某一个很窄的动态范围内,使得图像对比度低,导致图像中的细节常常看不清楚,给目标地质体特征的识别和探测带来困难。在成像测井图像处理中,常采用直方图修正的方法增强图像质量[3-6]。胡刚等[3]采用全局直方图均衡对测井图像进行处理[3],虽增强了图像的整体信息,但对局部细节信息增强能力不明显。闫建平等[4]在利用图像直方图局部均衡化实现测井图像动态增强研究的基础上,采用Morphing技术对颜色进行线性插值,消除动态增强产生的台阶问题,效果明显。涂继辉等[5]利用某井段一个局部区域数据统计其直方图,从而实现局部图像的最佳均衡,该方法对局部细节的增强效果优于全局直方图均衡算法,但这种动态直方图增强方法会引入平坦区域少量噪声放大和过度增强等问题。

本文针对超声测井图像对比度低、细节模糊的问题,借鉴一种限定对比度直方图均衡方法对超声测井图像进行增强,达到了突出局部细节,限制平坦区域过度增强的目的,对进一步提高成像测井图像解释精度具有重要意义。

1 直方图均衡化

直方图均衡又称为全局直方图均衡,它是最基本的直方图均衡方法[7],其基本思想是通过对原图像进行某种灰度变换,使变换后图像的直方图能均匀分布,这样就能使原图像中具有相近灰度且占有大量像素点的区域灰度范围展宽,使大区域中的微小变化显现出来,使图像更清晰。对于连续图像,设r代表原图像的灰度级,假定r归一化为0≤r≤1,r=0代表黑(最暗),r=1代表白(最亮),p(r)为原始图像灰度分布的概率分布函数,直方图均衡处理实际上就是寻找一个灰度变换函数T,使变换后的灰度值满足s=T(r),其中,s归一化为0≤s≤1,要求处理后图像灰度分布的概率密度函数p(s)=1,变换函数T(r)必须满足2个条件[8]:

(1) 在0≤r≤1区间内是单值且单调递增函数,确保灰度变换,原始图像的每个灰度级r都对应产生一个灰度级s,且灰度变换前后不倒置(使灰度级保持从黑到白的次序)。

(2) 在0≤s≤1区间内0≤T(r)≤1,确保输出灰度级与输入灰度级有同样的范围。

设原图像的灰度范围为[r,r+dr],包含的像素个数为p(r)dr,经过单调增的一对一变换,变换后灰度范围为[s,s+ds],包含像素个数为p(s)ds,变换前后的像素个数应相等,即

p(r)dr=p(s)ds

(1)

(2)

式(2)右边为直方图p(r)的累积分布函数,也是直方图均衡化对应的灰度变换函数。

对于数字图像离散情况,直方图均衡化的计算过程:

(1) 根据原始图像的统计信息,算出原始直方图,p(i)=ni/n。其中,i=0,1,…,L-1;L为灰度级的个数;ni表示各灰度级的像素数,n表示像素的总数。

(3) 借助累积直方图,映射各像素点的灰度值,运用变换得到公式(3),其中的INT()是取整符号。

j=INT[(L-1)pj+0.5]

(3)

(4) 确定新的映射i→j,据此,将原图像的灰度值f(m,n)=i修正为g(m,n)=j。

(5) 统计变换后各灰度级的像素个数nj,并由此给出处理后图像的直方图p(j)=nj/n。

超声测井图像进行直方图均衡化后,能够扩展图像的动态范围,增加图像的整体对比度,质量得到一定程度的提升,使感兴趣的特征更容易被检测或识别。如果超声测井图像中表示某些特征的像素点较少,当小到一定程度后,在式(3)中运用累积直方图映射各像素点的灰度值时,不同的pj有可能会得到相同的j,这样会造成灰度值合并现象,导致测井图像中细节丢失。另外,对于长达几千m测井数据形成的测井图像,不同的深度有不同的特征,使用全局直方图均衡算法,不可避免地对少数噪声点也进行了放大,如果存在局部区域过亮的部分,会有过度增强的问题[9-10],反而会使这部分图像包含的细节更加不清晰。因此,超声测井图像的增强,需要采用突出局部细节,限制平坦区过度增强的方法。

2 限定对比度直方图均衡化方法

限定对比度的自适应直方图均衡化算法由Karel Zuiderveld提出[11],基本思想是通过限制局部直方图的高度,从而可以限制噪声放大和削弱局部对比度的过增强,已在许多领域得到了广泛应用[10,12-13]。针对彩色的测井图像,本文将RGB图像变换到HSI(Hue Saturation Intensity)空间,仅对亮度分量进行限定对比度直方图均衡化处理,再将结果作RGB输出,这样保留了图像的颜色信息,处理前后图像的整体色调保持不变。

限定对比度直方图均衡化的基本流程如下[9]:

(1) 将输入原图像划分为若干大小相等的子块,每个子块相互连续并且互不重叠。

(2) 统计各子块内的局部直方图。

(3) 计算子块内所有像素数平均分配到的各灰度级值,并用参数ANum表示,若用GNum表示子块可能的灰度级,则以上过程可表示成

(4)

式中,XP、YP分别表示子块X、Y方向的像素数。

(4) 设定一个剪切系数CV,并假设其具有[0,1]的取值范围。针对不同的图像,通过仿真结果调整到最佳值,则实际剪切极限值NV可用公式(5)表示,其中round是四舍五入函数。

NV=ANum+round[CV(XP×YP-ANum)]

(5)

(5) 利用剪切极限,对局部直方图各个灰度级上的像素进行剪切,多余下来的像素数目重新分配到各直方图的各灰度级中,设已被剪切的像素总数为NClip,于是得到每个灰度级应该分到的剪切像素数NAcp。

NClip=∑{max[H(i)-NV,0]}

(6)

式中,H(i),i=i0,i1,…,in表示子块的局部直方图。

(7)

如果用CH表示重新分配后的直方图,则有

(8)

(6) 设经过以上分配后,像素点剩余数目为NLeft,则分配像素的步长为

(9)

从最小灰度级开始按步长搜索,遇到像素数小于剪切阈值的位置,则分配1个像素,完成从最小灰度级到最大灰度级的循环,直到NLeft为0,直方图分配才算完成,最终得到新的直方图。

(7) 对剪切后的每个子区域的灰度直方图分别进行直方图均衡化。

(8) 把每个子块的中心点作为参考点,获取其灰度值,对图像中的每一个像素进行灰度线性插值,采用双线性插值的方法,每个像素点的映射由其相邻的4个参考点对应区域的映射确定,如图1中黑色小矩形(x,y)代表任意目标点,f(x,y)表示该点待计算的灰度值,与其相邻的4个区域的中心点用黑色圆点表示,分别记为A(x1,y1)、B(x2,y2)、C(x3,y3)、D(x4,y4),其中x1=x3,x2=x4,y1=y2,y3=y4,则目标点f(x,y)的灰度值可由A、B、C、D灰度值的线性组合表示成公式(10),在边界区域的像素,其灰度计算用邻近2个样本点进行线性插值;在图像4个角的点用邻近1个样本点进行计算。

f(x,y)=a[bf(x1,y1)+(1-b)f(x2,y2)]+

(1-a)[bf(x3,y3)+(1-b)f(x4,y4)]

(10)

图1 双线性插值示意图

3 实验结果及分析

3.1 主观评价

图2 长庆宜×井图像增强效果图

图3 长庆宜×井图像的直方图

采用限定对比度的自适应直方图均衡化方法,使用版本为R2008b的Matlab软件对超声测井图像的增强效果进行验证。图2(a)是一幅长庆宜×井超声测井原始彩色幅度图的一段截图;图2(b)是使用全局直方图均衡算法处理的结果;图2(c)是使用限定对比度直方图均衡处理的结果。根据仿真结果,调整参数至最佳,其中分割的子块数目为8×8个,剪切系数CV为0.01。从图2中可以看出,2种方法都一定程度上增强了图像的对比度,但是图2(b)中对应原始图像高亮区域明显过度发亮,淹没了细节信息,造成图像局部区域过度增强,有一定程度失真。图2(c)中对原始图像中高亮局部区域的细节增强效果明显好于图2(b)中的效果,没有出现过度增强的问题。图3分别给出了图2(a)、(b)、(c)对应的灰度直方图信息,横坐标为像素灰度级,纵坐标为像素在某灰度级出现的个数。从图3(a)可以看出原始图像像素集中分布在灰度级50～240的范围,对比度偏低。图3(b)是经过全局直方图均衡算法处理后的直方图,其像素近似分布在0～255之间,同时在低灰度级和高灰度级像素比较集中,而中间部分比较稀疏,整体不均匀。从直方图还可以看出100左右两边附近的有些灰度级的像素个数为0,这是由于较少像素的灰度被合并了,淹没了图像的细节信息,同时在255灰度级附近像素个数为500左右,直方图有尖状化现象,这是导致亮区过度增强的原因。图3(c)经过限定对比度直方图均衡后,其直方图的动态范围被均匀扩展到10～255灰度范围内,没有出现细节退化和过度增强的问题。

3.2 客观评价

(11)

RPSN为最大信号量和噪声强度的比值,它的表达式为

(12)

对于常见的256级灰度图像,fmax=255。

(13)

把信息论中熵值的概率应用到图像信息源,假如图像的灰度级为[1,L],可以通过直方图得到各灰度级的概率ps(sk),k=1,2,…,L,这时,图像的熵为

(14)

图像处理后均方误差(EMS)越小、峰值信噪比(RPSN)越大说明处理效果好。信息熵越大,信息的无序度越高,包含的信息量越大。

表1给出了图2(a)的超声测井图像经过全局直方图均衡和限定对比度直方图均衡算法各参数相应的结果。从表1中可以看出,全局直方图均衡算法处理过的图像,其均方误差大于限定对比度直方图均衡算法处理过的图像,前者的峰值信噪比和信息熵均小于后者,表明采用限定对比度直方图均衡化算法具有较好的图像增强效果。

表1 超声测井图像的客观评价结果

4 结 论

(1) 超声测井图像灰度直方图通常会集中在很窄的动态范围内,导致图像对比度低、细节模糊现象。采用限定对比度直方图均衡方法对超声测井图像进行增强,达到了突出局部细节,限制平坦区域过度增强的目的。

(2) 对图像处理效果及直方图的形态进行分析,通过客观评价参数均方误差、峰值信噪比和信息熵对比2种直方图增强方法,仿真结果显示,限定对比度直方图均衡算法对于超声测井图像质量的改善比全局直方图均衡算法效果更好,对实际的超声测井资料评价具有一定的意义,也为测井图像处理提供了又一种可行的方法。

参考文献:

[1] 李长文, 强毓明, 余春昊, 等. 提高超声成像测井图像质量的途径 [J]. 石油仪器, 2009, 23(5): 51-54.

[2] 付青青, 谢凯. 应用维纳滤波技术改善超声波成像测井图像效果 [J]. 测井技术, 2012, 36(6): 581-584.

[3] 胡刚, 陈凡. 基于直方图均衡化的成像测井彩色图像增强 [J]. 国外测井技术, 2011, 4: 68-70.

[4] 闫建平, 首祥云, 邵在平, 等. 成像测井图像的动态增强及Morphing方法 [J]. 测井技术, 2006, 30(4): 364-366.

[5] 涂继辉, 余厚全, 李长文, 等. 基于超声测井图像的动态直方图均衡算法研究 [J]. 电视技术, 2011, 35(19): 113-114.

[6] 王晓峰, 彭天慈, 雷刚, 等. XRMI动态增强及全井眼成像方法研究与应用 [J]. 测井技术, 2015, 39(4): 432-437.

[7] 曾炜赫, 杨俊炜, 邢宇翔. 动态直方图双向均衡化的图像增强方法 [J]. 中国体视学与图像分析, 2014, 19(2): 129-139.

[8] GONZAKEZ R C, WOODS R E. DigitalImage Processing [M]. 阮秋琦, 阮宇智, 译. 2版. 北京: 电子工业出版社, 2007: 70-74.

[9] 赵静. 彩色图像增强 [D]. 西安: 西安电子科技大学, 2014.

[10] 张璞. 基于视频交通车辆检测系统的图像增强的研究 [D]. 青岛: 青岛大学, 2012.

[11] ZUIDERVELD K. Contrast Limited Adaptive Histogram Equalization [M]. Chapter VIII. 5, Graphics gems IV, San Diego: Academic Press Professional, Inc, 1994: 474-485.

[12] ALAMEEN Z, SULONG G, REHMAN A, et al. An Innovative Technique for Contrast Enhancement of Computed Tomography Images Using Normalized Gamma-corrected Contrast-limited Adaptive Histogram Equalization [J]. EURASIP Journal on Advances in Signal Processing, 2015, 2015(1): 1-12.

[13] 杨有, 李波. CLAHE和细节放大相结合的档案图像增强方法 [J]. 中国图象图形学报, 2011, 16(4): 522-527.

[14] 张德丰. Matlab数字图像处理 [M]. 北京: 机械工业出版社, 2012: 213-218.