利用信息熵的岩心图像自适应压缩感知重构①

唐国维,刘彦彤,张岩(东北石油大学 计算机与信息技术学院,大庆 163318)



利用信息熵的岩心图像自适应压缩感知重构①

唐国维,刘彦彤,张岩
(东北石油大学 计算机与信息技术学院,大庆 163318)

摘 要:针对BCS-SPL算法对岩心图像进行压缩感知重构的细节模糊的问题,提出一种利用信息熵的岩心图像BCS-SPL压缩感知重构算法.采用小波变换对岩心图像进行稀疏表示,对各子带进行多尺度分块,依据信息熵的大小自适应分配采样率并确定观测矩阵,通过维纳滤波结合Landweber迭代操作实现重构.实验结果表明,在相同采样率下,与原始的BCS-SPL算法相比,该算法的重构质量提高了2-4 dB.

关键词:岩心图像; 压缩感知; 重构; 自适应; 信息熵

岩心是油气田勘探开发中重要的基础地质资料,在推断沉积环境和生储盖组合研究中具有不可替代的作用.将岩心样本通过扫描方式以数字图像形式存储,已成为数字化油田建设的重要组成部分.由于多年的累积和不断新取心,导致岩心数据量及其庞大,因此必须对岩心图像进行压缩处理.通过对大量典型岩心图像分析,发现岩心图像普遍具有纹理信息丰富的特点,并且对比度很弱.因此,对于岩心图像的压缩与普通自然图像压缩有着不同的要求[1].目前,国内外各类岩心图像压缩重构算法均基于Shannon/Nyquist釆样理论[2,3].由于理论框架的原因,基于Shannon/Nyquist采样定理对岩心图像进行压缩重构,其采样数据具有非常大的冗余性,需要耗费大量的处理时间和存储空间,且压缩效果难以保证.

Donoho、Candès及Tao等人建立的压缩感知[4,5](Compressed Sensing,CS)理论指出,只要信号是稀疏的或者在某一变换空间是稀疏的或可压缩的,以远低于奈奎斯特采样率的速率随机采样,仍能够精确地重构原始信号.可见,CS突破了传统信源编码架构已经接近的理论极限,将其用于图像压缩可能获得意想不到的效果.由于直接使用CS方法重整幅图像的计算量相当巨大,Gan将图像分块技术运用到图像中来,即分块压缩感知(Block Compressed Sensing,BCS)[6]方法.BCS方法把整幅图像分成等尺寸的块,独立地对每个图像块进行观测和重构,这样大大降低了存储和计算成本,但是会在低码率下带来块效应.为此,Sungkwang等人提出BCS-SPL(Block Compressed Sensing-Smooth Projected Landweber)[7]算法,该算法通过高斯随机矩阵实现采样,图像的重构策略使用维纳滤波结合Landweber迭代算法实现,虽然改善了块效应,但是在一定程度上降低了重构质量.

根据油田的实际需要,本文将BCS-SPL算法用于解决岩心图像压缩与重构问题.通过大量实验发现,其效果并不理想.分析其原因在于岩心图像普遍包含丰富的目标纹理信息,直接使用通用的BCS-SPL算法必然导致重构图像细节模糊.所以,本文在BCS-SPL框架的基础上,提出一种利用信息熵的岩心图像BCS-SPL压缩感知重构算法.该算法在DWT(Discrete Wavelet Transform)域内,将变换分解后得到的每一级低频和高频子带进行分块,然后根据各个子带的信息熵得到其自适应采样子率进行自适应采样,再通过维纳滤波结合Landweber迭代实现图像的重构,达到进一步提高岩心图像的重构质量和改善视觉效果的目的.

1　BCS-SPL算法

根据压缩感知理论,假设x为从M个采样信号y中获得的长度为N的信号,且M«N,那么,可以从(1)式中恢复信号x .

其中,x∈RN,y∈RM,即x是一个N维向量,y是一个M维向量,Φ是一个具有采样率为S= M N的M´ N维的观测矩阵.由于x的数量远大于观测值y,理论上通过y恢复x是不可能的.然而,如果x足够稀疏,就能够使精确重构成为可能[8].

由于图像数据的多维性,采样过程的维数N会随着图像x的增大而迅速增加,导致存储观测矩阵Φ需要巨大内存空间,并且重构过程会产生相当大的计算量.为此,Gan提出分块CS方法,文献[6]给出了一个二维图像的CS范式.在这个技术中,图像采样是通过应用块到块基的随机矩阵实现的,重构是Landweber迭代结合平滑操作实现的.由于分块CS采样和平滑Landweber迭代重构相结合,所以称之为BCS-SPL技术.

假设一幅大小为N´ N的图像x被分成大小为B´ B的图像块,第i个图像块的向量表示记为xi,使用观测矩阵ΦB进行采样,得到观测值.

其中,B的大小根据图像重构的速率和重构的质量要求综合决定: 当B较小的时候,内存占用少且计算速度快; 当B较大的时候,图像的重构效果比较好.i =1Kn ,n =N2B2,ΦB是大小为MB×B2的正交观测矩阵,M=(M×B2)N2,M为对整幅图像的观测采样数.

在文献[6]中,维纳滤波被纳入基本Landweber迭代框架中,目的是为了去除块效应.本质上,这个操作对于Landweber迭代来说,除了固有稀疏性还能够起到平滑作用.具体地,维纳滤波步骤被插入到公式(3)、(4)中的Landweber迭代,可以看出,第i+1次的迭代图像x(i+1)近似值是通过x(i)得到的.

BCS-SPL算法将图像分块进行观测采样,从根本上减少了观测矩阵的存储量,使重构图像的效率明显提高,但对图像分块进行观测采样会割裂图像的整体信息,而其投影迭代过程中会产生块效应,去除块效应会带来额外的资源消耗和图像信息丢失,该算法采用维纳滤波来去除块效应,导致丢失图像的边缘和细节信息.在BCS-SPL算法的基础上提出的MS-BCS -SPL算法[9]在DWT域内,对变换分解后的每一级的子带进行分块采样,再通过平滑迭代重构图像.该算法兼顾了CS的计算开销和图像的重构质量,但是直接将其应用到岩心图像压缩感知重构的应用中,重构的岩心图像的细节信息仍然模糊.

2　岩心图像压缩感知采样与重构

2.1岩心图像频谱分析与信息熵计算

图像可以看作是一个离散的二维函数,其频谱│F(U,V)│可以由该图像矩阵作二维离散傅里叶变换得到.而二维离散傅里叶变换能够描述图像纹理近似周期模式的分布规律[10],所以本文使用基于傅里叶变换的频谱分析对岩心图像的频谱进行分析.图1(a)和(b)分别给出了岩心图像及其傅里叶频谱.

从图1(b)可以看出,在岩心图像的傅里叶频谱中,相对于远离坐标原点的边缘位置,其靠近坐标原点的中间部分并不是特别亮,也就是说,岩心图像在低频部分并没有显著的优势,而是在高频部分同样占有较重要的比例.这说明岩心图像在经过小波变换后,除了低频子带,每一级高频子带必然会包含更多的、不同方向的重要信息.

图1　岩心图像的傅里叶频谱

由于岩心图像的高频成分占有相对较多的比例,那么图像分块后,不同子块纹理不同,即包含信息量不同.由Mallat塔式小波分解理论可知,图像经分解层数为L的二维离散小波变换以后,分为3L+ 1个子频带,即1个低频子带和3个高频子带.理想的采样方法是信息量少的块少采样,信息量多的块多采样,在总采样率不变的情况下,将有限的资源有效地分配给纹理相对复杂的图像块.因此,本文改进了所有块都使用相同的采样率的BCS采样方法,根据子块间纹理结构不同引起信息量的差异的特点,采用信息熵作为纹理信息的度量.通过计算经DWT变换后各级高频子带的信息熵,得到自适应采样率.

信息熵反映了岩心图像中平均信息量的多少,即表示了岩心图像中灰度分布的聚集特征所包含的信息量,将岩心图像灰度值进行数学统计,便可得到每个灰度值出现的次数及概率.一般情况下,信息熵值越大表明图像信息保留的程度越好,其携带信息量的能力越强[11].定义岩心图像信息熵的计算公式为:

其中,pi表示岩心图像的概率密度函数,可利用直方图近似计算.利用信息熵的计算公式便可计算出岩心图像的信息熵.

2.2基于信息熵的自适应采样

在本文的改进算法中,把观测矩阵Φ分成两个部分: 一个是DWT多尺度变换矩阵Ω,而另一个是多尺度分块自适应观测矩阵Φ″,即Φ=Φ″Ω.假设Ω为L 级DWT分解,那么,Φ″是由3L+1个不同的观测矩阵组成.这时,被分成大小为Bl×Bl的图像块的图像x在l级的低频和高频子带分别通过自适应观测矩阵Φ″进行采样.自适应采样的实现步骤如下:

1)计算分解层数为L的DWT变换的l级的采样子率Sl:

在l级,根据其分块大小Bl使用矩阵Φl进行采样,会产生采样子率Sl[12].其中,设DWT基带子率S0为全采样率,即S0= 1.若l级的采样子率Sl,可以得到公式(6):

这里,Wl为l级的采样子率Sl的加权系数.加权系数Wl可由下式得到.

那么,整个图像的采样率为:

由此可知,当已知图像的目标采样率S和加权系数Wl后,由(8)式很容易求出S' ,再通过(6)式得到l级的采样子率Sl.表1给出了在不同的目标采样率S下,分解层数L=3的DWT变换实现的l级的采样子率Sl统计;

表1　L3=级DWT变换实现的采样子率统计

2)计算l级h ,v ,d子带的信息熵Hlh,Hlv和Hld;

3)计算l级h ,v ,d子带的自适应采样子率Sla,见公式(9).

根据自适应采样的实现步骤可以计算出自适应采样子率.表2给出了岩心图像1在不同的目标采样率S 下,分解层数L=3时的DWT变换实现的l级h、v和d子带的自适应采样子率Sla统计.

表2　L3=级DWT变换实现的高频子带自适应采样子率统计

根据表2的l级的h、v和d子带的自适应采样子率Sla,可以计算出相应方向的高频子带的观测采样数M .设MBl表示块大小为Bl时观测采样数,由计算.表3给出了岩心图像1在不同的目标采样率S下,L=3级DWT变换实现的每一级的h、v 和d子带的自适应观测采样数MBl统计.

表3　L3=级DWT变换实现的高频子带的自适应观测采样数统计

根据表3的每一级的h、v和d子带的自适应观测采样数MBl,可以看出观测采样数MBl会根据分解级数不同和高频子带方向不同而自适应改变,充分体现了其根据图像块所包含信息量不同观测采样数不等的自适应性,同样的工作量却能够保留更多的边缘和细节信息,从而提高岩心图像的重构质量.

2.3多尺度BCS自适应重构算法

在图像 DWT 稀疏变换域内,结合图像边缘的3´ 3维纳滤波和稀疏提升阈值处理[13]实现重构.维纳滤波在空间域实现,而平滑和阈值操作在变换域进行.该算法能够实现图像的快速重构,在DWT分解的每一级中的每个块都使用自适应观测矩阵Φ″和Landweber迭代操作.重构算法具体步骤如下:

步骤1: 利用MMSE(Minimum Mean Square Error)估计得到xi的近似解xi,从而得到图像的初始解xi;

1)用3´ 3邻域的自适应维纳滤波器去除图像分块所带来的块效应;

2)将滤波后的图像投影在凸集上,可由下式得到:

当Φ是正交矩阵的时候,即ΦΦT=1,上式可简化为:

3)用小波域双变量阈值[14]对投影结果进行滤波;

4)将滤波后的图像再次投影到凸集上;

5)判断并终止迭代,直到得到最优解.

3　实验仿真与结果分析

实验中用到的测试图像是512×512的两幅岩心图像,对其进行利用信息熵的分块自适应采样与重构,并与BCS-SPL和MS-BCS-SPL(Multiscale Block Compressed Sensing with Smoothed Projected Landweber)算法进行比较.本文算法、BCS-SPL和MS-BCS-SPL算法均使用双树复小波变换(Dual-Tree Complex Wavelet Transform,DTCWT)[15]作为稀疏基,采样时使用97双正交3级DWT作为多尺度变换矩阵Ω.Ω进行l级分解时,使用大小为Bl×Bl的图像块采样.该采样过程使用随机DCT(Discrete Cosine Transform)SRM(Site Recovery Manager)观测矩阵[16]实现.所有实验都在MATLAB R2013b环境下完成.本文算法中,当l=1,2,3时,图像块的大小分别为Bl=16,32,64,l级的每个方向的高频子带的采样子率Sl都使用表2的计算结果,然后根据表3的l级的每个方向的高频子带的自适应观测采样数MBl得到自适应观测矩阵Φl″.MS-BCS-SPL算法中,当l=1,2,3时,块的大小分别为Bl=16,32,64,每一级的采样子率都使用表1的计算结果.BCS-SPL算法中,B=32.图2、3给出了三种算法重构的岩心图像的部分实验结果.由图可见,当采样子率S=0.2时,本文提出的利用信息熵的分块自适应采样和多尺度重构的图像质量优于BCS-SPL算法约4 dB,也优于MS-BCS-SPL算法约2 dB.表3给出三种算法对岩心图像1和岩心图像2重构结果的峰值信噪比.

图2　岩心图像1的3种重构算法效果对比(S= 0.2)

图3　岩心图像2的3种重构算法效果对比(S= 0.2)

表4　三种算法重构结果的峰值信噪比PSNR(dB)

4　结论

借鉴DWT域的多尺度分块压缩感知技术,本文提出一种利用信息熵的岩心图像BCS-SPL压缩感知重构算法.该算法针对岩心图像特性,利用DWT的多分辨率和多尺度特性以及信息熵计算使每级分解层上每个方向的高频子带上的采样具有自适应性,所以其观测结果能够充分表示岩心图像的结构特点.因此,使用本文算法对岩心图像进行压缩重构,在重构质量和视觉效果方面都有所提升.

参考文献

1Zhan X,Zhang R,Yin D,Huo C.SAR image compression using multiscale dictionary learning and sparse representation.IEEE Geoscience and Remote Sensing Letters,2013,10(5): 1090–1094.

2Xiao C.Reconstruction of bandlimited signal with lost samples at its Nyquist rate—the solution to a nonuniform sampling problem.IEEE Trans.on Signal Processing,1995,43(4): 1008–1009.

3Yang F,Hu J,Li SQ.A total least squares reconstruction algorithm of UWB signals based on sub-nyquist sampling.Journal of Electronics and Information Technology,2010,32(6): 1418–1422.

4Donoho DL.Compressed sensing.IEEE Trans.on Information Theory,2006,52(4): 1289–1306.

5Candès E,Tao T.Near-optimal signal recovery from random projections: Universal encoding strategies.IEEE Trans.on Information Theory,2006,52(12): 5406–5425.

6Lu G.Block compressed sensing of natural images.2007 15th International Conference on Digital Signal Processing.Cardiff.IEEE.2007.403–406.

7Mun S,Fowler JE.Block compressed sensing of images using directional transforms.2009 16th International Conference on Image Processing.Cairo.IEEE.2009.3021– 3024.

8Candès E,Romberg J,Tao T.Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information.IEEE Trans.on Information Theory,2006,52(2): 489–509.

9Fowler JE,Mun S,Tramel EW.Multiscale block compressed sensing with smoothed projected landweber reconstruction.2011 19th European Signal Processing Conference.Barcelona.IEEE.2011.564–568.

10Chen J,Deng M,Xiao PF,Yang MH,Mei XM,Liu HM.Optimal spatial scale choosing for high resolution imagery based on texture frequency analysis.Journal of Remote Sensing,2011,15(3): 492–511.

11徐洁,张国海.基于熵的图书馆网络信息安全的风险评估.科技情报开发与经济,2008,18(14):3–5.

12Candès EJ,Wakin MB.An introduction to compressive sampling.IEEE Signal Processing,2008,25(2): 21–30.

13蒋业文,于昕梅.基于DWT的多尺度分块变采样率压缩感知图像重构算法.中山大学学报:自然科学版,2013,52(3): 30–33.

14Jia J,Jiao LC,Xiang HL.Using bivariate threshold function for image denoising in NSCT domain.Journal of Electronics and Information Technology,2009,31(3): 532–536.

15Selesnick IW,Baraniuk RG,Kingsbury NC.The dual-tree complex wavelet transform.IEEE Signal Processing,2005,22(6): 123–151.

16Do TT,Tran TD,Lu G.Fast compressive sampling with structurally random matrices.IEEE International Conference on Acoustics,Speech and Signal Processing.Las Vegas,NV.IEEE.2008.3369–3372.

Adaptive Compressed Sensing Reconstruction of Core Images Using Information Entropy

TANG Guo-Wei,LIU Yan-Tong,ZHANG Yan
(School of Computer and Information Technology,Northeast Petroleum University,Daqing 163318,China)

Abstract:Aimed at the details vague problem of Block Compressed Sensing-Smooth Projected Landweber compressed sensing reconstruction of core images,a Block Compressed Sensing-Smooth Projected Landweber compressed sensing reconstruction of core images using information entropy is proposed.The method introduces discrete wavelet transform into the sparse representation and conducts multiscale block for each subband,and then adaptively allocates the sampling rates and determines the measurement matrix.The reconstruction can be achieved by Wiener filter combined with Landweber iterative.The experimental results show that the reconstruction quality is improved by 2-4dB compared with that of Block Compressed Sensing-Smooth Projected Landweber under the same sampling rates.

Key words:core images; compressed sensing; reconstruction; adaptive; information entropy

基金项目:①东北石油大学研究生创新科研项目(YJSCX2015-034NEPU);黑龙江省教育厅科学技术研究项目(12521050)

收稿时间:2015-09-24;收到修改稿时间:2015-11-11