全变分2范数有限差分小波域成像测井图像缺陷修复算法

林光明,张 健,张 顺

(1.武汉理工大学 计算机科学与技术学院, 湖北 武汉 430070; 2.武汉理工大学 信息工程学院,湖北 武汉 430070; 3.华北油田测井公司,河北 任丘 062550)



全变分2范数有限差分小波域成像测井图像缺陷修复算法

林光明1,张 健2,张 顺3

(1.武汉理工大学 计算机科学与技术学院, 湖北 武汉 430070; 2.武汉理工大学 信息工程学院,湖北 武汉 430070; 3.华北油田测井公司,河北 任丘 062550)

针对超声成像测井中因各种不确定因素而导致测井图像中出现带状异常区域的现象,通过分析修复异常区域所采用的常规方法尤其是小波域修复方法的不足,文章提出了一种能够有效解决该类问题的全变分2范数有限差分小波域超声测井图像修复算法。首先,设计基于图像相邻行相关值与绝对误差值相结合的“链条”异常区域自动识别算法,实现对异常区域的准确定位;然后，针对全变分小波域图像修复过程中出现的“分片常数”效应,设计一种全变分2范数有限差分小波域修复算法来有效消除该影响。通过实际测井资料处理结果与标准测井图像效果的比对，证明该算法进一步提升了针对“链条”异常区域的修复效果,为后期图像资料中该类状况的修复提供了有效的解决方案。

2范数;有限差分;改进;全变分小波域;缺陷修复

0 引 言

超声成像测井技术在工程上应用较为广泛。当今世界上最先进的测井技术是美国研发的三维特征测井方法,该方法利用井下阵列传感器沿井眼纵向、径向扫描从而采集大量的层数据信息,再根据采集到的数据使用图像处理方法获得井眼附近的三维图像或井壁的二维图像[1],因此,所获得的目标图像比曲线方式更加直观、精确和方便。

但由于超声成像测井设备在移动测绘成像过程中,受制于井下阻力不均匀的干扰,会出现骤停现象,而由于测井电缆弹性形变的存在,导致在骤停深度上测井设备会反复进行扫描成像,相应的深度测量仪器也不停地在匀速运动,造成测井数据与深度位置的不匹配,将会在测井图像中出现带状分布区域,称为“链条效应”[2]。

为解决该问题,文献[3]采用偏微分方程方式对超声成像测井获得的图像进行修复;文献[4]利用总体变分模型(TV模型)对测井图像进行修复;文献[5]利用成像纹理设计了测井成像的修复算法。后续学者也提出了许多改进方法并取得了一定效果,但总体上这些修复方案在针对曲线或面积较小图像区域修复时效果较好,而当图像异常或区域较大时则效果不佳。

因此,本文提出了一种基于图像相邻行绝对误差的带状区域自动识别算法,在实现对该类区域准确定位的同时,借助于小波域在图像修复方面的优势,设计了一种全变分2范数有限差分小波域成像测井修正算法,以期能够较好地解决该类问题。

1 “链条效应”分布区域自动识别方法

在进行超声成像测井图像修复过程中,“链条”区域的确定是该修复过程中的重要环节,若无法准确定位“链条”的异常区域,则图像修复算法将无法准确应用。由于“链条”区域是测井成像仪器在同一深度进行重复扫描所形成的异常区域[2],因此在该区域中相邻两行数据的相关性非常好;根据这个特性,对于链条区域的自动识别可采用基于相邻行相关值与绝对误差值相结合的自动识别方法。

“链条”异常区域如图1所示。

图1 “链条”异常区域

假设成像测井图像资料的起始深度为D1st,而结束深度为D1end,则从超声成像测井图像资料的起始深度D1st开始计算。“链条”异常区域的自动识别算法步骤如下：

(1) 读取当前图像深度行D1cur的超声成像测井图像数据A(D1cur,k),然后计算与其相邻m行图像数据A(D1cur+m,k)的相关值,即:

(1)

其中,N为成像测井仪器对井壁扫描1圈的采样点数;A(D1cur+p,k)为与深度D1cur相邻p行的图像数据。

(2) 判断步骤(1)中计算得到的m个相关值是否都大于相关值阈值Rmax,即

(2)

若(2)式中R(D1cur)≤Rmax,则说明深度行D1cur不在链条异常区域,跳转步骤(1),令D1cur=D1cur+1,继续对下一行图像进行判别;若(2)式成立,则记录满足条件的m′≤m值,跳转步骤(3)。

(3) 步骤(2)计算得到m′值,则图像区域[D1cur,D1cur+m′]为链条异常区域的预估值,跳转步骤(4),利用相邻行绝对误差值对图像链条的异常区域进行准确定位。

(4) 读取当前图像深度行D1cur的成像测井图像数据A(D1cur,k),然后计算与其相邻p行图像数据A(D1cur+p,k)的列误差绝对值和,即

(3)

其中,m′的值为上一步计算得到的异常区域的预估值。

(5) 判断步骤(4)中计算得到的m′个相关值是否都小于阈值Smin，即

(4)

若(4)式中S(D1cur)>Smin,则说明深度行D1cur不在异常区域,跳转步骤(4),令D1cur=D1cur+1,继续对下一行图像进行判别;若(4)式成立,则记录满足条件的m″≤m′值。

(6) 输出D1cur和m″,即为“链条”区域的限定值。

上述步骤通过采用基于相邻行相关值的方式对异常区域进行预估,然后利用绝对误差值作为补充对异常区域进一步细化,使得定位区域更加准确可靠,为修复算法奠定良好基础。

2 小波域图像修复算法

“链条效应”产生的主要原因是由于井下阻力不均匀而造成测井数据与深度位置不匹配,工程上常用的方法是根据“链条”区域上下边缘的测量数据进行插值得到[6]。实践证明,这种插值方法在处理较窄“链条”时能够得到很好的效果,但当“链条”区域过大时则修复效果不佳。文献[7]提出了一种全变分小波域图像修复算法,但是在实际应用该修复算法时发现,在小波稳态解中会出现“分片常数”效应,影响图像修复的效果。本文结合2范数采用有限差分方法对该小波域图像修复算法进行改进。

2.1 全变分小波域图像修复算法模型

假设图像可表述为:

(5)

其中,u0(x)为成像测井原始图像;n(x)为图像噪声,一般实验环境下取为高斯白噪声。假设图像大小为n×m,则z(x)的标准正交小波变换形式为:

(6)

在小波变换图像表示中,β={βj,k}为小波系数。需要注意的是,该模型中并未区分系数的高、低频,这样处理的原因为:① 简化计算;② 该模型在高、低频范围内允许出现一定系数的丢失。在实际应用中,噪声与数据丢包情况不可避免,造成的后果是小波系数的丢失;而剩余的系数又要受到噪声的干扰。为准确表述小波域图像,文献[8]提出含噪声的全变分小波域图像修复算法模型为:

(7)

其中,当不存在噪声时，λ=0;若j、k∈D,其中D属于图像待修复区域,则λj,k=0;若j、k∉D,则λj,k为正常数。引入时间变量t，采用梯度下降法可求得小波系数为:

(8)

图2 全变分特性

图2表征的是3条单调递增的曲(直)线函数。若u∈[a,b],并且函数在[a,b]范围内单调递增,满足u(a)=α,u(b)=β,那么上述3种曲线函数的全变分结果都相同,即

(9)

从图像形态学的角度上分析这3条线段是完全不同的,特别是利用(7)式进行修复时,小波域系数与曲率的计算不在同一个域中。因而,采用全变分方式修复图像时会导致小波系数稳态值中出现“分片常数”效应,进而影响图像修复效果。

2.2 改进型全变分小波域图像修复算法模型

图像修复在像素域上的主要工作是修正异常区域的像素值[8-9],而在小波域上的主要工作是补充经过小波变换之后缺失的系数。如前所述,全变分小波修复方式在稳态时存在“分片常数”效应,这会直接影响修复效果。而全变分小波域修复方式因其具备较强的边缘保持能力,适应于异常“链条”区域较窄情况下的修复;但当异常“链条”区域较宽时,该算法平滑性的不足容易产生边缘阶梯效应。而2范数∫|u(β,x)|2dx同全变分∫|u(β,x)|dx相比正好相反,图2中3条曲线的2范数关系为:

(10)

其各向同性扩散系数与沿梯度的扩散系数一致,不会产生阶梯效应，但是会导致边缘模糊。综合这2种方式的优点,对小波域模型进行改进,得到的模型如下:

(11)

其中,j、k∈D,D属于图像待修复区域;τ∈[0,1]为权重系数,用来调整全变分与2范数的比重。

2.3 改进算法模型的有限差分求解

假设βD={βj,k|j,k∈D}为“链条”异常区域待修复小波系数,采用有限差分方式求解如下:

(12)

若母波函数φ满足(11)式,那么(12)式由分部积分可得:

(13)

为简化公式表示,令

(14)

计算过程中需要在小波域与像素域之间进行转换。首先在像素域中进行计算,对于“链条”异常区域的每个点(i,j)有:

(15)

然后将像素域转换至小波基投影,即

(16)

改进算法的有限差分求解步骤如下：

(2) 若i

(3) 令βold=βnew,并根据(14)～(16)式计算“链条”异常区域像素域到小波基投影。

(4) 对于每个像素点(j,k),更新其小波系数，即

(17)

其中

(18)

(5) 更新算法偏差E=‖βnew-βold‖2,并令i=i+1,转步骤(2)。

(6) 终止条件满足,终止算法并输出修复后的小波参数βD={βj,k|j,k∈D}。

3 仿真结果与分析

仿真数据资料来自阿尔7井的超声成像测井图像。阿尔7井位于内蒙古自治区二连盆地,隶属于华北油田二连公司。原始图像、“链状区域”自动定位识别仿真图以及采用3种不同类型小波算法修复效果的对比如图3所示。仿真参数采用Daubechines 9/7双正交小波;图像噪声参数选取λ=0.05[10-11];仿真硬件平台配置为AMD A10-5 800 K的CPU 、8 G内存、Windows7旗舰版操作系统(由于图内空间有限,故图中“全变分2范数有限差分”简称为“全变分2范数”)。

从图3可以看出,原始测井图像存在较明显的“链条”异常区域,这在实际成像测井图像中是普遍存在的;通过“链条”区域自动识别算法准确辨别出的异常区域,黑线位置即为“链条”区域。如前所述,2范数小波趋向于图像修复的平滑处理而全变分小波侧重于边缘保持,2范数小波给出的修复结果非常平滑,纹理间的边界相对模糊;全变分小波则注重于边缘的保持,因此修复后的图像纹理特征更加清晰;全变分2范数小波的修复结果介于两者之间,较好地平衡了图像修复的平滑性与边缘保持的要求。综上所述,全变分2范数小波在原有小波域特性的基础上较好地提升了图像修复的效果。

图3 异常区域自动识别及3种小波算法效果比对

下面从客观角度对比全变分2范数小波、全变分小波以及2范数小波修复方法[12]在成像测井图像修复过程中,峰值信噪比(PeakSingaltoNoiseratio,PSNR)与小波系数丢失率[13]的特征关系,如图4所示。

图4 PSNR与小波系数丢失率特征关系

由图4可知,全变分2范数小波修复算法的峰值信噪比值明显优于其他2种修复算法的,而受损原始图像的PSNR值最低。这说明应用该3种算法修复受损图像都达到了一定的效果,而全变分2范数小波算法的修复效果明显优于另外2种算法的。同时,随着小波系数丢失率的上升,3种算法所得PSNR值都有所下降,说明图像受损越严重,其PSNR值也越低。

在对比算法运行效率方面,得到在同等条件下算法运行时间随小波系数丢失率增加而变化的关联曲线,如图5所示。为稳定仿真结果,算法运行时间的取值为每种算法各自运行40次的平均值[14-15]。

从图5可以看出,随着小波系数丢失率的上升,3种算法的运行时间都在增加，说明图像受损越严重,算法的运行时间越长,这符合客观现实。2范数小波修复算法在小波系数丢失率低于20%时,其算法的运行时间皆短于其他2种算法;当小波系数丢失率处于20%～30%范围之内,其运行时间高于全变分2范数小波修复时间而低于全变分小波修复时间;当小波系数丢失率达30%以上时,其运行时间比另外2种算法大幅增加;而后两者的变化趋势相似,但全变分2范数小波修复算法的运行时间总体低于全变分小波修复算法,原因是改进的算法提高了修复效率。

图5 算法运行时间对比结果

4 结 论

本文针对成像测井图像中存在带状异常区域现象且利用常规修复算法效果不理想的情况,引入并改进了原有全变分小波域图像修复算法,提出了一种基于全变分2范数有限差分的小波域图像修复算法并进行了仿真与对比实验。该算法考虑了小波系数的方向选择性、传递性和聚集性;较好地解决了拥有较大带状异常区域的成像测井图像缺陷修复问题;同时在充分考虑图像结构性的基础上,有效去除了“分片常数”效应;实验结果验证了该算法的有效性与鲁棒性。实践证明，本文算法在修复效果与运行效率等方面均优于常规算法,能更好地提高真实超声图像的修复性能，为超声成像测井缺陷图像修复提供了一种新的有效思路,是一种有效可行的修复方法。

[1] 倪路桥,余厚全，李长文，等.基于纹理的超声成像测井图像“城墙效应”修复研究[J].测井技术,2010,34(5):428-432.

[2] 张健.一种连续带状分布区域修复方法在成像测井资料处理中的应用[J].科学技术与工程,2014,14(10):145-148.

[3] BERTALMIO M,SAPIRO G,CASELLES V,et al.Image inpainting[C]//Proceedings of the 27th annual Conference on Computer Graphics and Interactive Techniques.New York:ACM Press/Addison-Wesley Publishing Co,2002:417-424.

[4] CHAN T F,SHEN J H.Mathematical models for local nontexture inpaintings[J].SIAM Journal on Applied Mathematics,2002,62(3):1019-1043.

[5] CHAN T F,SHEN J H.Non-texture inpainting by curvature driven diffusions[J].Journal of Visual Communication and Image Representation,2001,12(4):436-449.

[6] 贾茜,易本顺,肖进胜.基于结构成分双向扩散的图像插值算法[J].电子与信息学报,2014,36(11):2541-2548.

[7] KIM G D,YOON C H,KYE S B.A single FPGA-based portable ultrasound imaging system for point-of-care application[J].IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control,2012,59(7):1386-1394.

[8] 张聚,王陈,程芸.小波与双边滤波的医学超声图像去噪[J].中国图象图形学报,2014,19(1):126-132.

[9] 谭小容,陈朝峰,查代奉.数字图像多小波逆变换及后置滤波算法[J].计算机应用,2014,34(5):1486-1490.

[10] TASINKEVYCH Y,TROTS I,NOWICKI A,et al.Modified synthetic transmit aperture algorithm for ultrasound imaging[J].Ultrasonics,2012,52(2):333-342.

[11] 胡文瑾,刘仲民,李战明.一种改进的小波域图像修复算法[J].计算机科学,2013,41(5):299-303.

[12] TANG Z J,HUANG L Y,DAI Y M,et al.Robust image hashing based on multiple histograms [J].International Journal of Digital Content Technology and its Applications,2012,23(6):39-47.

[13] 郑驰超,彭虎.基于编码发射与自适应波束形成的超声成像[J].电子与信息学报,2010,32(4):959-962.

[14] ZADEH P B,AKBARI A S.Image resolution enhancement using multi-wavelet and cycle-spinning[C]//Proceedings of 2012 UKACC International Conference on Control.Piscataway:IEEE,2012:789-792.

[15] WANG X.Image de-noising based-on multi-wavelet[C]//Proceedings of the 2009 International Forum on Information Technology and Applications.Piscataway:IEEE,2009:523-525.

(责任编辑 胡亚敏)

A total variational wavelet domain logging image defect inpainting algorithm based on two norm and finite difference method

LIN Guangming1,ZHANG Jian2,ZHANG Shun3

(1.College of Computer Science and Technology, Wuhan University of Technology, Wuhan 430070, China; 2.College of Information Engineering, Wuhan University of Technology, Wuhan 430070, China; 3.Huabei Oilfield Well Logging Company, Renqiu 062550, China)

Aiming at the phenomenon that in the process of ultrasonic image logging, the zonal abnormal region will appear in the images due to several uncertain factors, the inefficiency of normal ways including wavelet domain restoration way especially to clear away the zonal abnormal region is analyzed, and an effective total variational two norm finite difference wavelet domain ultrasonic logging images inpainting algorithm is put forward to solve this problem. Firstly, a kind of automatic recognition algorithm about chain abnormal region based on adjacent rows and absolute error value is proposed, which realizes the accurate location of imaging zonal region effectively. Then, aiming at the presence of “piecewise constant” effect in the process of the total variational wavelet domain image restoration, the total variational two norm finite difference wavelet domain ultrasonic logging images inpainting algorithm is proposed to eliminate the impact effectively. By comparing the experimental results of the actual logging image data processing with the standard logging image data, it is shown that this algorithm improves the inpainting effect and provides an effective solution for the inpainting of similar problems in image materials in the later stage.

two norm; finite difference; improvement; total variational wavelet domain; defect inpainting

2015-04-27；

2015-09-16

国家自然科学基金资助项目(41374148;41372155);湖北省自然科学基金青年基金资助项目(2014CFB248)和湖北省教育厅科技计划资助项目(Q20111306)

林光明(1985-),男,山东青岛人,武汉理工大学硕士生;

张 健(1981-),男,湖北荆州人,博士,武汉理工大学讲师.

10.3969/j.issn.1003-5060.2016.10.010

TP391

A

1003-5060(2016)10-1347-06