0 概述

随着陆地资源逐渐匮乏, 人们开始关注开发深海资源, 因此, 水下图像处理成为相关领域的研究重点。声视觉系统是传统的海洋探测手段, 而目前人们主要通过光学照明来探测深海。海水中的光在成像时, 由于吸收和散射等因素, 会使成像出的图像对比度下降, 而低对比度的图像会对后续的图像特征提取^[1]和目标跟踪识别^[2-3]等产生不良影响。因此, 获取清晰的水下图像具有重要意义。

针对水下图像模糊的问题, 研究者提出较多增强算法, 如文献[4]提出自适应的图像增强算法, 文献[5]提出基于暗原色先验原理的算法, 文献[6]提出基于卷积神经网络的算法, 文献[7]提出基于迭代直方图的算法, 均有效提高了水下图像的清晰度。文献[8]对暗通道模型进行双边滤波操作, 取得了较好的去雾视觉效果, 但算法的实时性较差。文献[9]利用暗通道先验理论处理空气中的带雾图像, 改善了因光折射引起的图像模糊现象, 并且图像的细节保存完整。由于水下图像和雾图像都会因光散射引起图像对比度降低, 因此研究者将暗原色先验直接用于水下图像的预处理。文献[10]指出水下红光衰减得最快, 区别于空气, 因此, 其仅利用蓝色和绿色通道来估计透射率图, 而对红色通道的透射率单独估计, 得到了较好的输出结果。文献[11]指出空气中的雾图像和水下低对比度图像都是因光被吸收和散射所导致的, 因此, 其将空气中的暗原色先验去雾算法直接用于水下图像增强, 但实验证明该算法只能提高图像对比度, 不能避免明显的图像的颜色失真现象。

本文根据水下成像特点, 提出一种结合暗原色先验原理和颜色修正的水下图像增强算法, 以改善颜色失真现象, 同时提高图像清晰度。根据水下光的成像模型, 在求取暗通道图像时考虑红色通道的逆通道并改进全局背景光的估计方法, 使用自定义的模板选取颜色饱和度方差最小的块图像, 分别对块图像RGB通道进行最小滤波, 并选取亮度最大的像素点作为预估背景光。在此基础上, 使用Retinex算法对图像的R通道进行颜色校正, 根据图像各通道的颜色衰减系数比得到颜色校正后的图像。

1 水下图像增强算法 1.1 水下图像成像物理模型

由文献[12-13]提出的水下成像物理模型可知, 到摄像相机的总辐照度是直接光、前向散射光和后向散射光这三部分的线性叠加, 即:

$ {E_{\rm{T}}} = {E_{\rm{D}}} + {E_{\rm{F}}} + {E_{\rm{B}}} $

(1)

直接光是指成像设备接收到物体直接反射的光, 其光的强度表达式为:

$ {E_{\rm{D}}}(x,\theta ) = E(x,\theta ){\rm{exp}}[ - c(\theta )d(x)],\theta \in \{ R,G,B\} $

(2)

其中, E(x, θ)是物体所在位置处的光强, E_D(x, θ)是摄像机接收的直接光光强, c(θ)表示由吸收和散射造成的总的衰减系数, θ表示图像的某一通道, d(x)表示物体到相机的距离。

前向散射光E_F的光强可以表示为:

$ {E_{\rm{F}}} = {E_{\rm{D}}}(x,\theta ) * g(x,\theta ) $

(3)

其中, *表示卷积运算, g(x, θ)为点扩展函数。

背景的散射部分可以表示为:

$ {E_{\rm{B}}}(x,\theta ) = {B_\infty }(\theta )\{ 1 - {\rm{exp}}[ - c(\theta )d(x)] $

(4)

其中, B_∞(θ)为背景光。

将式(2)~式(4)代入式(1)中, 可得:

$ \begin{array}{*{20}{l}} {{E_{\rm{T}}}(x,\theta ) = {E_{\rm{D}}}(x,\theta ) + {E_{\rm{F}}}(x,\theta ) + {E_{\rm{B}}}(x,\theta ) = }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} E(x,\theta ){\rm{exp}}[ - c(\theta )d(x)] + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \{ E(x,\theta ){\rm{exp}}[ - c(\theta )d(x)]\} * g(x,\theta ) + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {B_\infty }(\theta )\{ 1 - {\rm{exp}}[ - c(\theta )d(x)]\} } \end{array} $

(5)

可以把物体与成像设备的距离看成一个相对很小的数, 这样就可以忽略前向散射造成的图像模糊, 则式(5)表示为:

$ \begin{array}{*{20}{l}} {{E_{\rm{T}}}(x,\theta ) = E(x,\theta ){\rm{exp}}[ - c(\theta )d(x)] + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {B_\infty }(\theta )\{ 1 - {\rm{exp}}[ - c(\theta )d(x)]\} } \end{array} $

(6)

假如水下透射率t(x, θ)=exp[-c(θ)d(x)], 则式(6)可表示为:

$ \begin{array}{*{20}{l}} {{E_{\rm{T}}}(x,\theta ) = E(x,\theta )t(x,\theta ) + {B_\infty }(\theta )\{ 1 - t(x,\theta )\} \} }\\ {\theta \in \{ R,G,B\} } \end{array} $

(7)

1.2 基于R逆通道暗原色先验的水下增强

通过实验可知, 可以直接用图像中的暗像素来估计雾和光的透射率, 暗通道公式表示为:

$ {J_{{\rm{ dark }}}}(x) = \mathop {{\rm{min}}}\limits_{y \in \varOmega (x)} [\mathop {{\rm{min}}}\limits_{c \in \{ R,G,B\} } {J_c}(y)] $

(8)

其中, J_c(y)为J(x)的某个通道, 表示以像素为中心的一个窗口。由式(8)可知, 暗通道图是分别对图像3个颜色通道取最小值运算和局部区域取最小值(它们的运算不分先后)得到的。根据暗原色先验原理, 可得:

$ {J_{{\rm{ dark }}}}(x) = \mathop {{\rm{min}}}\limits_{y \in \varOmega (x)} [\mathop {{\rm{min}}}\limits_{c \in \{ R,G,B\} } {J_c}(y)] \to 0 $

(9)

水下图像和雾图像都是被浑浊的传输介质所退化的图像, 虽然它们具有相似性, 但是因为水下光的衰减率和大气光衰减率不同, 所以不能用空气中的暗原色先验直接恢复E(x, θ)。因此, 在求取水下图像的暗原色通道时考虑红色通道的逆通道^[14]和蓝绿颜色通道, 用J^wdark表示水下暗原色通道, 表达式为:

$ {J^{{\rm{wdark}}}}(x) = \mathop {{\rm{min}}}\limits_{y \in \varOmega (x)} \{ {\rm{min}}[E(y,{R_{{\rm{inv}}}}),E(y,G),E(y,B)]\} $

(10)

其中, R_inv表示红色通道的逆通道, 其值为1和归一化红色通道像素点的亮度差值。

图 1分别是采取空气中的暗原色通道方法和本文方法求取的暗原色通道对比。根据水下光的衰减率可知, 红光在水下先衰减完, 图像中的红色通道像素值低, 水下图像主要呈现出偏蓝色或者绿色, 如图 1(a)所示。采用文献[9]方法直接求取水下图像暗原色通道, 将导致整体偏暗, 如图 1(b)所示, 其是考虑红色通道逆通道后求取的结果, 可以看出离相机远即意味着光的衰减程度大, 因此, 其暗通道区域的亮度也就越高。

1.3 改进的背景光估算方法

由于水下图像十分模糊, 图像的暗通道整体亮度偏高, 因此文献[9]方法对背景光估计不准确, 对图像的恢复会产生不利的影响。在水下采集图像时, 会有人工照明和白色鱼群的存在, 若用这种方式选取, 会使增强后的图像偏白, 而且亮度不均匀。因此, 本文提出一种改进的全局背景光估计方法, 以提高背景光估计的精确度, 具体如下:

1) 设计一个计算颜色饱和度方差的模板, 图像的大小为M×N, 经试验得出最佳的模板大小是M/2×N/2, 设定步长2, 从模糊图像的左上方开始移动, 每次移动一步就计算该模板区域的颜色饱和度方差, 选择全局最小的区域颜色饱和度方差, 因为方差值越大则此区域受光源的影响就越大, 该区域不适合做背景光估计, 所以待完全移动完成后选取颜色饱和度方差最小的区域作为背景光预估计区域。

2) 对预估计区域每个通道进行最小值滤波, 得到在暗原色通道中的亮度值排在前0.1%的点, 对这些点的亮度值进行排序, 选取最大值作为大气光的值B_{θ, ∞}。

本文通过实验对比不同模板尺寸下所得背景光恢复图片的信息熵和PNSR, 实验结果如表 1所示。

表 1反映了不同尺寸大小模板估计出的背景光对图片质量的影响, 可以看出, 以尺寸为M/2×N/2模板估计出的背景光恢复图像质量最好。

1.4 水下透射率估算

通过R通道的逆通道来求水下图像的暗通道, 在求取水下透射率时, 图像的通道要分开考虑。对式(7)等号两边取最小值运算得出:

$ \begin{array}{*{20}{l}} {\mathop {{\rm{min}}}\limits_{y \in \varOmega (x)} [1 - {E_{\rm{T}}}(y,\theta )] = }\\ {t(x,\theta )\mathop {{\rm{min}}}\limits_{y \in \varOmega (x)} [1 - E(y,\theta )] + [1 - {B_\infty }(\theta )][1 - t(x,\theta )]}\\ {\theta \in \{ R\} } \end{array} $

(11)

$ \begin{array}{*{20}{l}} {\mathop {{\rm{min}}}\limits_{y \in \varOmega (x)} [{E_{\rm{T}}}(x,\theta )] = }\\ {\mathop {{\rm{min}}}\limits_{y \in \varOmega (x)} [E(x,\theta )]t(x,\theta ) + {B_\infty }(\theta )[1 - t(x,\theta )]}\\ {\theta \in \{ G,B\} } \end{array} $

(12)

分别对红通道和蓝绿通道取最小值运算, 对式(11)和式(12)分别除以1-B_∞(θ)和B_∞(θ), 可以得到:

$ \begin{array}{*{20}{l}} {\mathop {{\rm{min}}}\limits_{\theta \in \{ R\} } \left\{ {\mathop {{\rm{min}}}\limits_{y \in \varOmega (x)} \left[ {\frac{{1 - {E_{\rm{T}}}(y,\theta )}}{{1 - {B_\infty }(\lambda )}}} \right]} \right\} = }\\ {t(x,\theta )\mathop {{\rm{min}}}\limits_{\theta \in \{ R\} } \left\{ {\mathop {{\rm{min}}}\limits_{y \in \varOmega (x)} \frac{{[1 - E(y,\theta )]}}{{1 - {B_\infty }(\lambda )}}} \right\} + \mathop {{\rm{min}}}\limits_{\theta \in \{ R\} } [1 - t(x,\theta )]} \end{array} $

(13)

$ \begin{array}{*{20}{l}} {\mathop {{\rm{min}}}\limits_{\theta \in \{ G,B\} } \left\{ {\mathop {{\rm{min}}}\limits_{y \in \varOmega (x)} \left[ {\frac{{{E_{\rm{T}}}(x,\theta )}}{{{B_\infty }(\theta )}}} \right]} \right\} = }\\ {\mathop {{\rm{min}}}\limits_{\theta \in \{ G,B\} } \left\{ {\mathop {{\rm{min}}}\limits_{y \in \varOmega (x)} \left[ {\frac{{E(x,\theta )}}{{{B_\infty }(\theta )}}} \right]} \right\}t(x,\theta ) + \mathop {{\rm{min}}}\limits_{\theta \in \{ G,B\} } [1 - t(x,\theta )]} \end{array} $

(14)

水下暗原色通道J^wdark→0, 即:

$ \begin{array}{*{20}{l}} {{J^{{\rm{ wdark }}}}(x) = \mathop {{\rm{min}}}\limits_{y \in \varOmega (x)} \{ {\rm{min}}[E(y,{R_{{\rm{ inv }}}}),E(y,G)]}\\ {E(y,B)\} = 0} \end{array} $

(15)

由式(13)~式(15)可得:

$ \begin{array}{*{20}{l}} {\mathop {{\rm{min}}}\limits_\theta [t(x,\theta )] = 1 - \mathop {{\rm{min}}}\limits_\theta \left[ {\mathop {{\rm{min}}}\limits_{y \in \varOmega (x)} \left( {\frac{{1 - E(y,R)}}{{1 - {B_\infty }(R)}}} \right)} \right],}\\ {\left. {\mathop {{\rm{min}}}\limits_{y \in \varOmega (x)} \left( {\frac{{E(y,G)}}{{{B_\infty }(G)}}} \right),\mathop {{\rm{min}}}\limits_{y \in \varOmega (x)} \left( {\frac{{E(y,B)}}{{{B_\infty }(B)}}} \right)} \right]}\\ {\theta \in \{ R,G,B\} } \end{array} $

(16)

由于红光在水下衰减最快, 因此红光的透射率最小, 即:

$ \mathop {{\rm{min}}}\limits_\theta [t(x,\theta )] = t(x,R) $

(17)

文献[15]实验得出图像的全局背景光的值与衰减系数和散射系数的关系式为:

$ {B_{\theta ,\infty }} \propto {b_\theta }/{c_\theta } $

(18)

在2类不同的海水中进行9组不同波长对应散射系数的数据的实验得出散射系数与波长的关系式:

$ b(\lambda ) = ( - 0.001{\kern 1pt} {\kern 1pt} {\kern 1pt} 13\lambda + 162{\kern 1pt} {\kern 1pt} {\kern 1pt} 517)b({\lambda _r}) $

(19)

其中, λ为波长, 本文选择红光波长620 nm、绿光波长540 nm、蓝色波长450 nm。为参照波长的散射系数, 此处以R通道为参考, 计算其他2个通道衰减系数比为:

$ {\frac{{{c_G}}}{{{c_R}}} = \frac{{{b_G}{B_{R,\infty }}}}{{{b_R}{B_{G,\infty }}}}} $

(20)

$ {\frac{{{c_B}}}{{{c_R}}} = \frac{{{b_B}{B_{R,\infty }}}}{{{b_R}{B_{B,\infty }}}}} $

(21)

则水下各通道的透射率为:

$ \begin{array}{*{20}{l}} {t(x,R) = 1 - \left\{ {\mathop {{\rm{min}}}\limits_R \left[ {\mathop {{\rm{min}}}\limits_{y \in \varOmega (x)} \left( {\frac{{1 - E(y,R)}}{{1 - {B_\infty }(R)}}} \right)} \right],} \right.}\\ {\left. {\mathop {{\rm{min}}}\limits_{y \in \varOmega (x)} \left( {\frac{{E(y,G)}}{{{B_\infty }(G)}}} \right),\mathop {{\rm{min}}}\limits_{y \in \varOmega (x)} \left( {\frac{{E(y,B)}}{{{B_\infty }(B)}}} \right)} \right\}} \end{array} $

(22)

$ t(x,G) = [{t_R}(x)]\frac{{{c_G}}}{{{c_R}}} $

(23)

$ t(x,B) = [{t_R}(x)]\frac{{{c_B}}}{{{c_R}}} $

(24)

本文采用文献[16]中的导向滤波算法对各通道透射率进行细化, 该算法避免了文献[8]利用区域为单位估计出的透射率, 会产生块效应, 如图 2所示。

最后通过透射率图和估计出的全局背景光的值可恢复出清晰的图像, 其光强表达式为:

$ \begin{array}{*{20}{l}} {E(x,\theta ) = \frac{{{E_{\rm{T}}}(x,\theta ) - {B_\infty }(\theta )}}{{{\rm{min}}[t(x,\theta ),{t_0}]}} + {B_\infty }(\theta )}\\ {\theta \in \{ R,G,B\} } \end{array} $

(25)

其中, t₀可以防止求取出的透射率低导致图像的颜色过度失真。图 3为上述方法恢复后的图像。

1.5 颜色较正

从图 3中可以看出, 估算水下透射率后恢复的图像仍然会偏蓝绿色, 而且呈现出光照不均的现象, 因此, 还需要对增强后的图像进行颜色校正。本文首先采用Retinex算法^[17]处理增强后的图像, 得到亮度和反射分量, 然后用高斯核对原图做卷积操作得到低通滤波后的图像D(x, y), 最后与原图像做差分得到高频图像G(x, y), 即得到R(x, y)。将其转化到对数域分别处理, 则有:

$ \begin{array}{*{20}{l}} {{R_i}(x,y) = {\rm{ln}}[{r_i}(x,y)] = {\rm{ln}}[{s_i}(x,y)] - {\rm{ln}}[{l_i}(x,y)]}\\ {i \in \{ R,G,B\} } \end{array} $

(26)

本文用r_R(x, y)根据各颜色通道的衰减系数比求出r_G(x, y)和r_B(x, y), 再将r_i(x, y)从对数域转换为实数域R_i(x, y)。

$ {r_B}(x,y) = [{r_R}(x,y)]\frac{{{c_B}}}{{{c_R}}} $

(27)

$ {r_G}(x,y) = [{r_R}(x,y)]\frac{{{c_G}}}{{{c_R}}} $

(28)

图 4是颜色校正后恢复的图像, 可以看出, 图像的色偏以及光照不均得到了较好的改善, 并且突出了图像的细节, 清晰度得到提高。

1.6 算法流程

本文基于改进全局背景光估算法, 结合暗原色先验算法恢复出增强后的图像, 并采用Retinex理论进行颜色校正, 得到清晰的图像, 本文方法流程如图 5所示。

2 实验与结果分析

实验采用3幅水下原始图像, 分别用HE算法^[18]、同态滤波算法^[19]、基于暗原色先验的算法^[14]、MSRCR算法^[20]和本文算法进行处理, 然后对上述算法进行实验结果对比, 如图 6~图 8所示。

观察图 6~图 8可知:HE算法处理的结果使图像失真, 经过同态滤波的图像亮度被降低, 有明显偏色现象, 基于暗原色先验算法不能很好地区分前景和背景, 也不能避免图像的偏色现象, 虽然MSRCR算法去除了基于暗原色先验算法的弊端, 但图像颜色淡化且失真, 经过本文算法处理的图像前景和背景的区分度高, 色偏校正显著改善, 明显地优化了图像的视觉效果。

上述图像主观对比是建立在本文做了颜色校正的基础上进行的, 为论证本文算法的有效性和公平性, 对上述经典算法处理的图像和本文算法得出未经过颜色校正的图像进行PSNR计算, 如表 2所示。可以看出, 本文算法第1组和第2组的PNSR值低于MSRCR算法, 这是因为缺少清晰的图像来对比, 水下的图像降质十分严重, 而算法处理效果好, 与原始图像的方差较大, 从而使PSNR值降低。

为进一步验证本文算法的有效性, 对上述图像的清晰度(平均梯度)以及信息熵进行分析与统计, 对处理结果进行客观评价, 如表 3和表 4所示。

可以看出, 本文算法处理得到的图像清晰度高于其他算法。信息熵是度量图像信息量多少的指标, 信息熵越高说明图像所所含的信息越多, 本文算法信息熵值也高于其他算法。从主观和客观性上分析都可以得出本文算法更具优越性, 图像恢复的清晰度较高。

3 结束语

本文提出的水下图像增强算法基于水下成像模型, 根据暗原色先验原理分别考虑红色通道的逆通道和蓝绿色通道求出暗通道图像, 在此基础上利用模板分块对区域颜色饱和度方差进行计算, 得出方差最低的块图像进行背景光估计, 避免了光源、白色物体等不利因素对图像全局背景光的影响, 最后结合透射率图得到增强后的图像。该算法利用图像颜色通道的衰减系数比和Retinex理论对增强后的图像进行颜色校正, 改善了水下图像的偏色现象, 并且提高了图像整体的清晰度。实验结果表明, 本文算法对于水下严重降质图像效果明显, 处理后的图像不易过度失真且细节得到突出。由于每次做图像增强时需要人工提取特征, 此步骤较为繁琐, 因此下一步将采用深度学习的网络框架提取特征映射图, 设计一个端到端的图像增强算法, 从而提高鲁棒性和实时性。