多尺度水下偏振成像方法∗

韩平丽刘飞张广陶禹邵晓鹏

1)(西安电子科技大学物理与光电工程学院,西安 710071)

2)(中国科学院长春光学精密机械与物理研究所,应用光学国家重点实验室,长春 130033)

3)(陆军炮兵防空兵装备技术研究所,北京 100012)

4)(北京理工大学,北京 100081)

(2017年9月11日收到;2017年11月14日收到修改稿)

1 引 言

水下光学成像技术因结构简单、图像分辨率高、易解析等优点,在海洋资源勘探、海底环境监测、海洋搜救及军事侦察领域具有重要应用价值[1−3].由于水中悬浮的微小粒子、可溶性有机物等混沌介质及气泡等非均匀性因素对光波产生严重的散射作用[4],背景散射光叠加在目标弹道光上形成噪声,使得图像产生一层“雾”的感觉,导致图像对比度大幅降低,细节信息严重丢失[5].因此,有效去除背景散射光,提供清晰的高对比度图像在水下光学成像领域具有重要研究价值.

近年来,国内外许多学者均对水下成像技术开展了研究[6−9],其中水下偏振成像技术由于结构简单,性价比高等特点备受关注[10−12].以Huang等[13],Schechner等[14]以及Gilbert和Pernicka[15]的研究为代表,利用散射光的偏振特性分离场景中目标信息光和背景散射光,可获取清晰的目标图像.水下偏振成像技术的研究重点普遍集中于背景散射光和目标信息光分离算法的改进,但重建图像受限问题长期困扰水下偏振成像技术的广泛应用,究其原因,与上述两部分光分离过程中产生的噪声放大现象不无关系.本研究针对以上问题,利用图像分层思想,结合小波变换的多尺度特性,全面考虑散射光的空间特性和频谱特性,提出多尺度水下偏振成像方法,用以解决传统水下偏振成像技术中存在的噪声放大问题及其导致的成像质量受限问题.

多尺度水下偏振成像方法的核心在于对于目标和背景差异较明显的高对比度基础层利用联合双边滤波进行处理;而对目标本身对比度较低但细节信息丰富的细节层进行多尺度小波收缩处理,最后通过小波逆变换进行图像重建.大量水下偏振成像实验结果证明该方法不仅有效抑制了偏振成像过程中的噪声放大,而且能够显著提高重建图像对比度,丰富图像细节,提供更多有效的目标信息.

2 水下偏振成像方法

2.1 水下偏振成像模型

水体的吸收作用主要导致光波能量的损失,散射才是影响成像质量的主要原因[16].水下偏振成像技术通过提取光波偏振信息并利用其分离背景散射光和目标信息光,获得清晰的目标图像[17].由图1所示水下偏振成像原理示意图可知,探测器所接收到的总光强ETotal由两部分构成,分别为经水中粒子散射并被探测器接收到的背景散射光EBS和经目标反射最终到达探测器的目标信息光ETI,二者之间存在如下所示的物理关系:

其中目标信息光ETI的图像形式为目标场景清晰图像,可知获取水下清晰图像的关键在于将背景散射光EBS从总图像中分离,得到目标信息光ETI.结合光波的偏振特性,ETotal,EBS和ETI间的关系可以细化为如(2)式所示的形式[18]:

结合(1)—(3)式,可以将背景散射光EBS和目标信息光ETI的分布情况分别表示为偏振度的函数:

图1 水下偏振成像原理示意图Fig.1.Schematic of underwater polarization imaging.

因此,目标信息光与背景散射光能够有效分离取决于二者偏振度的准确计算.针对这一问题,本文通过改变偏振片透光轴方向采集两幅偏振子图像后,选取场景中背景散射光均匀且不含目标的背景区域估算背景散射光偏振度pBS[11].由(4)式可知,任何对目标信息光偏振度pTI的错误估计,都会使背景散射光EBS和目标信息光ETI的分离效果变差,使得分离出的背景散射光存在部分目标信息光,增加了背景散射光和目标信息光的相关性,进而影响重建图像质量.为保证目标信息光偏振度的准确计算,采用互信息IM(mutual information)来衡量背景散射光和目标信息光之间的相关性[20]:

其中,m和n分别为背景散射光图像EBS和目标信息光图像ETI的灰度级;prob(m,n)表示联合概率分布函数;prob(m)和prob(n)为边缘分布函数.根据互信息的定义,当互信息取最小值时,背景散射光和目标信息光之间相关性最低,此时分离出的背景散射光EBS中将不再包含目标信息光ETI,达到二者的最佳分离效果.在此条件下目标信息光的偏振度pTI能够取得最佳估计值.图2(a)为利用互信息估计目标信息光偏振度的示例,图2(b)为实验采集的原始水下强度图像,通过上述背景散射光偏振度估算方法计算得pBS为0.51,此时根据图2(a)所示互信息计算可得目标信息光偏振度pTI为0.24.在此基础上,利用水下偏振成像技术分离背景散射光和目标信息光,重建图像如图2(c)所示.图2(d)和(e)分别为原始强度图像和重建图像的局部放大结果,相比原始图像,重建图像对比度显著提高,背景散射光得到有效移除,包括细节纹理信息及清晰度在内的图像整体成像质量得到有效提升.

图2 目标信息光偏振度计算及重建结果示例 (a)互信息随目标信息光偏振度的变化;(b)水下原始强度图像;(c)水下偏振成像算法重建图像;(d)原始强度图像局部放大结果;(e)偏振重建图像局部放大结果Fig.2.pTIcalculation and an example of reconstructed images:(a)MI of EBSand ETIversus pTI;(b)raw image of an underwater scene;(c)reconstructed result of(b)by polarization imaging method;(d)zoomed-in view of the region of interest in(b)marked with a red rectangle;(e)zoomed-in view of the region of interest in(c)marked with a red rectangle.

2.2 水下偏振成像模型噪声分析

水下偏振成像方法中图像重建效果依赖于两幅相互正交的偏振子图像如(4)式所示,而正交子偏振图像作为两个独立的测量值,与重建图像呈线性关系.因此,重建图像中的噪声方差与两偏振子图像的噪声方差存在如下关系[5]:

式中ETI为重建图像,为重建图像的噪声方差,分别为两偏振子图像的噪声方差.

图3 噪声放大倍率随偏振度pBS和pTI变化示意图Fig.3.Noise amplification ratio as a function of pBSand pTI.

当背景散射光偏振度和目标信息光偏振度相等时,重建图像ETI的噪声方差趋于无穷大,此时水下退化图像无法得到有效恢复.假设噪声只与成像系统有关,即σa=σi=σ0成立,此时重建图像的噪声方差为

重建图像噪声放大倍率σT/σ0和σB/σ0随偏振度pBS和pTI的变化情况如图3所示.总体来看,噪声放大倍率均随pBS和pTI的增大而减小,但当背景散射光偏振度pBS=1且目标信息光偏振度pTI=0时,重建图像的噪声放大倍率达到最小值随着二者偏振度值越来越接近,噪声放大倍率逐渐接近,表现为重建图像中的噪声增大,目标信号湮没于噪声中.因此,水下偏振成像过程中,重建所得的清晰目标图像中的噪声方差σT是恒被放大的,且在pBS=pTI时,重建图像噪声达到无穷大,此时清晰场景图像将无法被复原.以图4所示的原始强度图像和传统偏振处理结果的频谱分布为例,处理后图像的整体频谱分布波动明显增强,这表明在目标得到有效突出的同时,重建图像中的噪声被严重放大.

图4 (a)原始图像(图2(b))的频谱分布图;(b)传统水下偏振成像算法重建图像(图2(c))的频谱分布图Fig.4.(a)Spectral distribution of a raw underwater image(2(b));(b)spectral distribution of the reconstructed image(2(c)).

3 多尺度水下偏振成像方法分析

根据2.2节的分析,水下偏振成像中噪声放大现象是一个不可避免的问题.分析发现,水下偏振图像主要分为包含目标但与背景差异较明显的高对比度区域,以及目标本身对比度较低但含有丰富细节信息的低对比度区域两部分.基于两部分信息的差异性,本文采用图像分层处理的思想建立多尺度水下偏振成像方法,重建高对比度、高信噪比的水下图像.在高对比度的基础层,利用联合双边滤波法来对噪声进行相应抑制[21];而在低对比度细节层,结合小波变换的多分辨率特性[22],对其进行多尺度小波收缩处理,最后结合小波逆变换对水下退化图像进行高质量重建.

清晰的目标场景图像ETI可以表示为

其中BTI为图像基础层,DTI为图像细节层.基础层BTI利用

所示的联合双边滤波来进行提取[23],其中ETI(ξ)表示目标场景图像ETI在ξ处的像素值;w(x,ξ)为双边滤波的内核函数,

σs和σr分别表示空间域和变换域的核尺寸,通过控制核尺寸的变化来获取图像基础层;g为导向图[24].由于双边滤波综合考虑了图像空间域核和变换域核的信息,故当图像中像素点x与ξ的欧氏距离或像素值距离增大时,双边滤波内核函数都会随之改变,进而影响图像基础层的提取精度.因此,为保证在基础层的良好保边去噪效果,同时有效保留原始图像的纹理信息,通过导向图g实现对原始图像的迭代处理,以精确提取图像基础层.

通过(8)—(10)式获取图像的基础层和细节层后,结合小波变换的能量集中性和多分辨时频特性,对细节层图像进行小波变换处理[25],

其中a为尺度因子,b为反应因子.通过改变尺度因子和反应因子可以获得任意中心频率、任意带宽和任意时间所对应的信号.此外,由于小波变换后的小波系数中幅值较大的部分包含了图像的重要信息;而幅值较小的部分则被认为是噪声变换后的系数,故通过选择双边滤波器的核函数w(x,f)作为窗口函数对导向图和细节层图像进行变换,保留或收缩小波变换系数,去除图像噪声,得到

ψ为小波基函数,GT(x,f)与是定义在时频窗口Fx中频率f、尺度因子a及反应位移b的函数;窗口Fx的大小和Np的大小相同.此外,由于直接使用变换域内核会保留均值为零的噪声信号,丢弃频率幅值较大的有用信号.而与无用的噪声信号相比,有用的目标信号在某些频谱处的幅值比较高.因此,为实现丢弃噪声信号、保留目标信号的目的,对频谱系数距离取倒数,定义收缩因子w(x,f),

式中γf主要用来调控变换域像素值核形状,而在(14)式所示的收缩因子中,导向谱G(x,f)主要起引导作用.对于幅值较高的频谱(对应原图像有用的目标信息),导向谱会引导收缩因子保留该部分变换系数;而对于幅值较低的频率(对应的噪声信号),导向谱就会引导收缩因子丢弃该信号.通过在变换域窗口Fx内平均所有收缩变换后的小波系数,并通过逆变换获得低对比度细节层信号结合高对比度基础层信号BTI,获取下次迭代的导引图像g,

通过对导引图像的迭代更新,修正小波变换系数,在有效抑制目标场景图像噪声的同时对受损细节信息进行不断完善,使最终图像逐渐逼近原始清晰目标场景图像.

4 实验结果与分析

为验证和说明多尺度水下偏振成像方法的有效性,根据图1所示的水下偏振成像原理搭建实验,采集偏振子图像.通过脱脂牛奶和自来水溶液模拟实际海洋环境中悬浮粒子对光波的散射和吸收情况[14],溶液配比为30 mL脱脂牛奶(分子大小约为0.040.3µm)和60 L自来水,二者均匀混合于体积为55 cm×55 cm×30 cm的水槽中.

为直观地表示实验介质的光学特性,计算其光学厚度.根据光学厚度的定义,光学厚度为介质透过率的自然对数,

式中τ为介质光学厚度,T为介质透过率.光学厚度的计算可以通过介质散射系数和成像距离实现[26],

式中τ为介质光学厚度,µs为介质散射系数(其大小与介质粒子类型和浓度有关[27]),d为物体与探测器之间的物理距离.结合(16)及(17)式可知,光学厚度可以反映介质透过率与成像距离的关系.经计算,实验中所用介质光学厚度为1.155.

图5为实验结果,其中图5(a)为实验采集的原始水下图像,图5(b)为传统偏振成像结果,相比原始图像,图5(b)背景散射光得到有效移除,对比度明显提高.但图5(h)中的图像频谱图频谱波动较大,说明处理结果中图像细节信息增加的同时,噪声信息也被放大.图5(c)所示的多尺度水下偏振成像结果不仅视觉效果明显改善,对比度显著提升,细节信息增加,图5(i)所示频谱图表明图像噪声得到明显抑制,图5(d),(e)和(f)所示的局部放大结果进一步对图像细节进行了直观展示.

图5 (a)原始强度图像;(b)传统偏振成像方法处理结果;(c)多尺度水下偏振图像处理结果;(d),(e)和(f)分别为图(a),(b)和(c)的局部放大图;(g),(h)和(i)分别为图(a),(b)和(c)的频谱强度分布图Fig.5.(a)A raw underwater image;(b)and(c)reconstructed image by traditional polarization imaging method and multi-scale polarization imaging method;(d),(e)and(f)zoomed-in view of the region of interest in(a),(b)and(c);(g),(h)and(i)spectral distributions of(a),(b)and(c).

图6 (a),(b)和(c)分别为图5(a),(b)和(c)的R,G,B三通道像素强度值统计;(d)为图5(a),(b)和(c)第210行像素强度值统计曲线Fig.6.(a),(b)and(c)are the pixel intensity distribution of channel R,G,and B of Fig.5(a),(b)and(c);(d)horizontal line plot at the vertical position pixel 210 for Fig.5(a),(b)and(c).

图6对应原始水下图像(即图5(a)),传统水下偏振重建图像(即图5(b))和多尺度水下偏振方法处理图像(即图5(c))的R,G,B三个色彩通道的像素强度统计值.对比发现,图6(c)中像素强度统计分布较其他两组结果均有明显提高,对应图像动态范围增大,即对各色彩通道均起到了一定的拉伸作用,增强了图像的层次感.此外,图6(d)所示为分别选取图5(a),(b)和(c)中穿过背景和目标的第210行(由上至下)像素的强度分布曲线,该曲线直观地表征了成像场景中背景散射光和目标信息光之间的差异情况及图像对比度变化情况.经多尺度水下偏振成像方法重建结果(蓝色曲线)与传统水下偏振成像结果(红色曲线)在背景区域强度的像素强度值(接近于0)基本相当,均小于原始强度图像(绿色曲线)的像素强度值(约为130).与此同时,三幅图像在目标位置的像素强度值大小相近(约为180),说明传统偏振成像结果对比度较原始强度图像得到有效提升.此外,如图6(d)右上角局部放大图所示,在目标位置处,经多尺度水下偏振成像方法重建后图像的像素强度值曲线波动明显增大,说明图像的细节轮廓及纹理信息得到了凸显.

在此基础上,为验证所述方法的普遍适用性和有效性,实验中选取了多个不同类型的目标,结果如图7所示.以图7(a)为例,目标物体上“我们d”几个字受光波强散射的影响,细节不易分辨;而重建图像中同一目标上“e一辈”几个字却能直接辨识.图7(e)所示的书签左侧强度图像中三列汉字几乎完全无法被识别,但图7(e)右侧经多尺度水下偏振成像方法处理之后的图像中右边三列汉字清晰可见.图7中其余图像存在类似的结论,也充分证明本文所述方法的普遍适用性.

图7 不同目标的多尺度水下偏振成像方法重建图像Fig.7.Reconstructed images of different underwater scenes utilizing multi-scale underwater polarization imaging method.

表1所列为四种常用的图像质量客观评估参数计算数据,其中平均梯度可以反映图像的边缘和细节等高频信息,其值越大,表明图像边缘越清晰,细节越明显;标准差反映图像灰度等级的离散程度,其值越大,表明图像反差越大,图像越清晰;图像对比度则表示图像整体亮暗区域之间的比率,其值越大,亮暗渐变层次越多,图像信息越丰富;而峰值信噪比则直观地反映出了图像信噪比的变化情况.总体看来,多尺度水下偏振成像图像的质量显著提升.如平均梯度较原始强度图像普遍提升了3.5倍左右,图像标准差提升了1.2倍左右;与此同时,图像对比度相比于原始强度图像提升了4倍左右,图像信噪比也有较为显著的提升.这都表明经多尺度水下偏振成像方法重建后,图像质量尤其在图像对比度和图像细节的改善以及噪声抑制方面有显著提升,且与图像主观评价和分析结果一致.

表1 水下成像结果的客观评价参数Table 1.Objective evaluations of underwater images.

5 结 论

水下偏振成像技术能够有效利用光的偏振特性解决水下图像质量退化问题,提高图像对比度.但在传统偏振处理方法中,目标信息光和背景信息光偏振度估算对重建图像质量存在显著影响,尤其在两者偏振度相等时,由于传统偏振处理产生的噪声放大,几乎无法有效复原图像.本文利用图像分层思想,结合小波变换的多尺度特性,提出了多尺度水下偏振成像方法,用以抑制噪声影响,重建高对比度的清晰图像.实验结果表明,该方法能够有效提升图像的对比度,抑制噪声放大,获得高对比度、高信噪比的清晰的水下场景图像,为水下偏振成像技术的应用提供了理论支持.与此同时,本文方法也存在不足,该方法实现过程中仍需对同一场景先后采集两幅偏振方向正交的图像,对运动目标的可行性较低.

[1]Panetta K,Gao C,Agaian S 2016IEEE J.Ocean.Eng.41 541

[2]George M J R,Kattawar W 1999Appl.Opt.38 6431

[3]Harvey E S,Shortis M R 1998Mar.Technol.Soc.J.32 3

[4]Zhao X W,Jin T,Chi H,Qü S 2015Acta Phys.Sin.64 104201(in Chinese)[赵欣慰,金韬,池灏,曲嵩 2015物理学报64 104201]

[5]Schechner Y Y,Averbuch Y 2007IEEE Trans.Pattern Anal.Mach.Intell.29 1655

[6]Guan J G,Zhu J P,Heng T 2015Chin.Phy.Lett.32 074201

[7]Lewis G D,Jordan D L,Roberts P J 1999Appl.Opt.38 3937

[8]Schechner Y Y,Karpel N 2005IEEE J.Ocean.Eng.30 570

[9]Miller D A,Dereniak E L 2012Appl.Opt.51 4092

[10]Kattawar W,Gray D J 2003Appl.Opt.42 7225

[11]Treibitz T,Schechner Y Y 2009IEEE Trans.Pattern Anal.Mach.Intell.31 385

[12]Han J,Yang K,Xia M,Sun L,Cheng Z,Liu H,Ye J 2015Appl.Opt.54 3294

[13]Huang B J,Liu T G,Han H F,Han J H,Yu M X 2016Opt.Express24 9826

[14]Dubreuil M,Delrot P,Leonard I,Alfalou A,Brosseau C,Dogariu A 2013Appl.Opt.52 997

[15]Gilbert G D,Pernicka J C 1967Appl.Opt.6 741

[16]Schettini R,Corchs S 2010EURASIP J.Adv.Signal Process2010 1

[17]Liu F,Shao X P,Gao Y,Xiang L B,Han P L,Li G 2016J.Opt.Soc.Am.A33 237

[18]Han P L,Liu F,Yang K,Ma J Y,Li J J,Shao X P 2017Appl.Opt.56 6631

[19]Liu F,Shao X P,Xiang L B,Gao Y,Han P L,Wang L 2015Chin.Phys.Lett.32 114203

[20]Zhao L Y,Lü B Y,Li X R,Chen S H 2015Acta Phys.Sin.64 124204(in Chinese)[赵辽英,吕步云,厉小润,陈淑涵2015物理学报64 124204]

[21]Knaus C,Zwicker M 2013Proceedings of the 20th IEEE International Conference on Image ProcessingNew Jersey,USA,September 15–18,2013 p440

[22]Liu F,Cao L,Shao X P,Han P L,Xiang L B 2015Appl.Opt.54 8116

[23]Li J C,Huang S X,Peng Y X,Zhang W M 2012Acta Phys.Sin.61 119501(in Chinese)[李金才,黄思训,彭宇行,张卫民2012物理学报61 119501]

[24]Papari G,Idowu N,Varslot T 2017IEEE Trans.Image Process26 251

[25]Myint S W,Zhu T,Zheng B J 2015IEEE Geosci.Remote Sens.12 1232

[26]Dubreuil M,Delrot P,Leonard P,Alfalou A,Brosseau C,Dogariu A 2013Appl.Opt.52 997

[27]Piederrière Y,Boulvert F,Cariou J,Jeune B L,Guern Y,Brun G L 2005Opt.Express13 5030