1 引言

水下光学成像具有重要的应用价值，但也具有挑战性^[1]。受水体及水中粒子等散射介质影响，在水下拍摄的图像存在目标模糊及对比度低等问题，难以对图像进行利用和分析^[2]。

根据文献[3]和文献[4]，人们主要将对比度作为判断图像质量的依据。关于获得高对比度水下图像的技术大体可分为两类：一类是基于数字图像处理技术，主要是水下图像增强^[5-7]与水下图像复原^[8-10]；另一类是基于成像系统，主要是水下距离选通系统^[11]、水下结构光成像系统^[12]等。受海洋动物利用偏振信息来改善视觉^[13]的启发，水下偏振成像技术作为一种解决问题的新思路被提出来。近几十年的相关研究表明，通过偏振技术可以抑制成像中的后向散射光，并通过成像算法的后处理进一步克服后向散射光对图像质量的影响。

Schechner将大气中的偏振去雾算法思想^[14]应用到水下，该方法^[2]采用被动光源进行水下目标成像，利用偏振差分成像得到复原图像，算法过程简单、处理耗时短、易于实现，但模型中的全局参量通过人机交互选取一个常量来代替。Treibitz^[4]重建了在主动光照明条件下的水下偏振成像模型，该算法在不同海域中获取的正交偏振图像的后处理中进行了验证，均取得了较好的实验效果，但算法中将后向散射光的偏振度视为常量，导致复原的水下图像有可能出现明暗不均的情况。Li^[15]在Schechner的研究基础上加入图像增强作为预处理，依然得到了理想的水下去雾图像，但仍没有考虑后向散射光偏振度为全局变量这一因素。此外，还有采用圆偏振照明的水下偏振成像去雾技术^[16]。

为解决现有算法中对参数估计不佳的问题，依据实际课题背景，本文采用主动偏振光进行照明，首先对正交的偏振图像进行图像增强预处理，通过分析水下主动偏振成像模型，推导出后向散射光偏振度特征参量，利用多项式拟合函数估算出全局变量的后向散射光偏振度，再迭代求解出相对于最优图像质量的目标反射光偏振度，最终得到水下去雾图像。通过实际的水下偏振成像实验进行验证，结果表明，本文算法与现有算法相比，可提高图像对比度，丰富图像信息，为水下目标探测和识别提供支持。

1 水下主动偏振成像模型

在水下成像时，探测器接收到的辐射信号有2个来源。第1个是场景目标的辐射，其辐射在水中被吸收和散射减弱，称为目标反射光；第2个是环境照明，环境光经过水中粒子的多次散射照射到探测器上，称为后向散射光^[1]。水下图像的数学表达式为

$I(x,y) = S(x,y) + B(x,y)$ (1)

式中： $ I(x,y)$为水下成像系统获得的图像； $ S(x,y)$为目标反射光，该信号实际由2个分量组成，分别为直接透射光和前向散射光； $ B(x,y)$为后向散射光。

目标反射光包含了水下目标受到水体吸收前的辐射信息（即为直接透射光），又包含了直接透射光向前传播时与介质互相作用而产生的前向散射光，吸收和前向散射均使得图像模糊。前向散射光在水下近距离成像时对图像质量影响较小，因此本文忽略前向散射光的因素。后向散射光是由主动照明光与介质发生多次散射作用而形成的背景光，如图1所示。由于主动光源照度的不均匀性，后向散射光实际是不均匀的。如前所述，对比度影响着图像质量，导致对比度下降的主要原因是后向散射光，因此水下去雾算法就是克服后向散射光引起的“影响”。

图 1. 散射过程示意图

Fig. 1. Schematic diagram of scattering process

下载图片 查看所有图片

利用偏振成像系统，得到光强度最大和最小时的2幅偏振正交图像，其强度分别记为 $ {I_{\max }}(x,y)$和 $ {I_{\min }}(x,y)$，则由式（1）可推得：

${I_{\max }}(x,y) = {S_{\max }}(x,y) + {B_{\max }}(x,y)$ (2)

式中： $ {S_{\max }}(x,y)$和 $ {B_{\max }}(x,y)$为偏振成像系统接收到目标场景“最亮”时的目标反射光和后向散射光的强度。同理：

${I_{\min }}(x,y) = {S_{\min }}(x,y) + {B_{\min }}(x,y)$ (3)

根据米氏散射规律，后向散射光为部分偏振光，因此后向散射光偏振度为

${p_{{\rm{scat}}}} = \frac{{{B_{\max }}(x,y) - {B_{\min }}(x,y)}}{{{B_{\max }}(x,y) + {B_{\min }}(x,y)}}$ (4)

对目标反射光作为非偏振光还是部分偏振光，以往有不同的看法，文献[2]和[3]认为随着传输距离的增加，目标反射光偏振态的影响在图像中的比重降低；文献[17]认为对于一些退偏度低的目标，其目标反射光的偏振态不能忽略。此模型考虑了这两种情况，因此目标反射光偏振度为

${p_{{\rm{obj}}}} = \frac{{{S_{\max }}(x,y) - {S_{\min }}(x,y)}}{{{S_{\max }}(x,y) + {S_{\min }}(x,y)}}$ (5)

下面为表达方便，省略了坐标量（x，y）。由（2）式和（3）式可得到总光强I为

$I = {I_{\max }} + {I_{\min }} = B + S$ (6)

由（2）式~（5）式可得到：

${I_{\max }} - {I_{\min }} = {p_{{\rm{scat}}}}B + {p_{{\rm{obj}}}}S$ (7)

最后，联立（6）式和（7）式，求解得到目标反射光分量S和后向散射光分量B：

$S = \frac{1}{{{p_{{\rm{scat}}}} - {p_{{\rm{obj}}}}}}[{I_{\min }}(1 + {p_{{\rm{scat}}}}) - {I_{\max }}(1 - {p_{{\rm{scat}}}})]$ (8)

$B = \frac{1}{{{p_{{\rm{scat}}}} - {p_{{\rm{obj}}}}}}[{I_{\max }}(1 - {p_{{\rm{obj}}}}) - {I_{\min }}(1 + {p_{{\rm{obj}}}})]$ (9)

经过理论推导可知，在给定后向散射光偏振度p_scat和目标反射光偏振度p_obj的情况下，就可从2张正交偏振图像中分离出后向散射光分量，目标反射光分量S即为复原图像。

2 预处理与全局变量的估计

2.1 偏振图像的预处理

本步骤借鉴文献[15]的算法思想，对正交偏振图像进行预处理，旨在拉伸正交偏振图像的直方图。首先，通过直方图拉伸处理偏振图像 $ {I_{\min }}(x,y)$。直方图拉伸的数学表达式如下

$I_{\min {\rm{pro}}}^{}(x,y) = \frac{{{I_{\min }}(x,y) - \min ({I_{\min }}(x,y))}}{{\max ({I_{\min }}(x,y)) - \min ({I_{\min }}(x,y))}}$ (10)

式中： $I_{\rm {min pro}}^{}(x,y)$是经过处理后的偏振图像； ${I_{\min }}(x,y)$是原图像； $\min ({I_{\min }}(x,y))$是原图像中灰度值最小的像素点； $\max ({I_{\min }}(x,y))$是原图像中灰度值最大的像素点。根据偏振度定义，正交偏振图像的偏振度为

$P(x,y) = \frac{{{I_{\max }}(x,y) - {I_{\min }}(x,y)}}{{{I_{\max }}(x,y) + {I_{\min }}(x,y)}}$ (11)

根据公式（11），正交偏振图像之间的偏振关系可以通过偏振度 $ P(x,y)$表示为

${I_{\max }}(x,y) = \frac{{1 + P(x,y)}}{{1 - P(x,y)}}{I_{\min }}(x,y)$ (12)

因此，根据公式（10）和（12）得到 $ {I_{\max }}(x,y)$的拉伸图像：

$I_{\max {\rm{pro}}}^{}(x,y) = \frac{{1 + P(x,y)}}{{1 - P(x,y)}}I_{\min {\rm{pro}}}^{}(x,y)$ (13)

通过直方图拉伸处理的2张正交偏振图像的偏振度与2张原始正交偏振图像的偏振度相同，如果直接拉伸2张正交偏振图像的直方图，则该偏振特性无法保证。此外，由于混浊水体介质的退偏性大，“最暗”的图像 $ {I_{{\rm{min}}}}\left( {x,{\rm{y}}} \right)$比“最亮”的图像 $ {I_{{\rm{max}}}}\left( {x,y} \right)$受到更少（除光学系统遮蔽外）的后向散射光的影响，目标反射光信息在 $ {I_{{\rm{min}}}}\left( {x,y} \right)$里占比更大，因此，可以认为直接对偏振图像 $ {I_{{\rm{min}}}}\left( {x,y} \right)$进行处理更为有效。

2.2 全局变量p_scat的估计

现行算法对后向散射光偏振度的估计是基于手动所选的区域值，近似地将后向散射光偏振度作为常量来处理。为提高去雾算法的稳定性，避免复原图像出现明暗不均的现象，本文重新分析了后向散射光偏振度，将其视为全局变量，并利用二元多项式函数进行数据拟合，拟合过程是基于最小二乘法。

不同物体具有不同的偏振特性，如石块的退偏性较强，金属的退偏性较弱^[18-19]。当目标反射光的偏振特性不可忽略时，无法直接利用偏振图像估算目标处后向散射光偏振度，但可以依托于目标周围区域的后向散射光偏振度的数据，通过二元多项式函数来估算目标处的后向散射光偏振度的值。实验表明，后向散射光偏振度p_scat受主动光源的影响，在整张图像上确实是有变化的，如图2所示。因此本文将后向散射光偏振度作为全局变量来处理是合理的。

图 2. 偏振度图

Fig. 2. Schematic diagram of polarization degree

下载图片 查看所有图片

将偏振图像的后向散射光偏振度值设为关于（x，y）的二元多项式函数，以图像左下角为原点建立三维坐标系，像素点分别为x轴和y轴，偏振度值为z轴，则二元多项式拟合函数数学表达式为

$\begin{split} {\textit{z}}(x,y) = & \sum\limits_{i = 0,j = 0}^{i = {n_1},j = {n_2}} {{a_{ij}}{x^i}{y^j}} \\ & (i + j \leqslant \max ({n_1},{n_2})) \end{split} $ (14)

式中：n₁为x的最高阶数；n₂为y的最高阶数。将（14）式展开推导，即：

$\begin{split} & {a_{00}} + {a_{10}}{x_1} + \cdot\cdot\cdot + {a_{{n_1}0}}x_1^{{n_1}} + {a_{01}}{y_1} +\cdot\cdot\cdot + \\ &\quad{a_{0n}}_2y_1^{{n_2}} +\cdot\cdot\cdot+ {a_{ij}}x_1^iy_1^j = {{\textit{z}}_1} \\ &{a_{00}} + {a_{10}}{x_2} +\cdot\cdot\cdot + {a_{{n_1}0}}x_2^{{n_1}} + {a_{01}}{y_2} +\cdot\cdot\cdot + \\ & \quad{a_{0n}}_2y_2^{{n_2}} + \cdot\cdot\cdot + {a_{ij}}x_2^iy_2^j = {{\textit{z}}_2} \\ &... \\ & {a_{00}} + {a_{10}}{x_n} +\cdot\cdot\cdot+ {a_{{n_1}0}}x_n^{{n_1}} + {a_{01}}{y_n} + \cdot\cdot\cdot + \\ & \quad{a_{0n}}_2y_n^{{n_2}} + \cdot\cdot\cdot + {a_{ij}}x_n^iy_n^j = {{\textit{z}}_n} \end{split} $ (15)

以上线性方程组写成矩阵形式：

$Aa{\bf{ = }}Z$ (16)

因为方程组无解，要找到一组所设的二元多项式系数 ${a_{00}},{a_{10}},\cdots,{a_{{n_1}0}},{a_{01}},\cdots,{a_{0{n_2}}},{a_{ij}}$，使平方偏差

$ \begin{split} \delta = & \sum {({a_{00}} + {a_{10}}{x_n} +\cdots + {a_{{n_1}0}}x_n^{{n_1}} + } \\ & {a_{01}}{y_n} + \cdots + {a_{0n}}_2y_n^{{n_2}} + {a_{ij}}x_n^iy_n^j - {{\textit{z}}_n}{)^2} \end{split} $ (17)

最小。经推导，二元多项式系数所满足的代数方程为

$ {{{A}}^{\bf{T}}}{{Aa}} = {{{A}}^{\bf{T}}}{{Z}} $  (18)

因此，系数最优解为

${{a}} = {({{{A}}^{\bf{T}}}{{A}})^{ - 1}}{{{A}}^{\bf{T}}}{{Z}}$ (19)

由此得到二元多项式拟合函数，可估算出全局变量的后向散射光偏振度。

关于x与y阶数的确定，以决定系数R²为依据进行迭代求得。决定系数R²可用来判断拟合函数的解释力，其正常取值范围为0～1，决定系数越接近1，表明拟合函数的解释能力越强，即对数据拟合也越好。决定系数R²定义为回归平方和SS_reg与总平方和SS_tot的比值，数学表达式为

${R^2} = \frac{{S{S_{{\rm{reg}}}}}}{{S{S_{{\rm{tot}}}}}}$ (20)

z假设后向散射光偏振度中包括 $ {{\textit{z}}_1},\cdots,{{\textit{z}}_n}$，共n个观察值，相对应的拟合值分别为 $ {f_1},\cdots,{f_n}$，则平均观察值为 $\overline {\textit{z}} = \dfrac{1}{n}\displaystyle\sum\limits_{i = 1}^n {{{\textit{z}}_i}} $。（20）式中回归平方和SS_reg为 $\displaystyle\sum\limits_i {{{({f_i} - \overline {\textit{z}} )}^2}} $，总平方和SS_tot为 $\displaystyle\sum\limits_i {{{({{\textit{z}}_i} - \overline {\textit{z}} )}^2}} $。

2.3 目标反射光偏振度p_obj的确定

本文将目标反射光也视为部分偏振光。通过水下主动偏振成像模型可知，目标反射光的偏振度为

${p_{{\rm{obj}}}} = \frac{{{S_{\max }}(x,y) - {S_{\min }}(x,y)}}{{{S_{\max }}(x,y) + {S_{\min }}(x,y)}}$ (21)

由（8）式可知，目标反射光偏振度p_obj对目标反射光分量S的影响比重较小，将目标反射光偏振度作为全局变量处理意义不大，而且会增加算法的复杂性。为适应实际应用的要求，通过迭代法选取一个相对于水下图像质量的最优p_obj值。为此，采用图像增强测量值^{[15, 17, 20]}（the value of measure of enhancement，EME）来量化图像质量。EME的数学表达式为

${\rm{EME}} = \left| {\frac{1}{{{k_1}{k_2}}}\sum\limits_{l = 1}^{{k_2}} {\sum\limits_{k = 1}^{{k_1}} {20\log \frac{{i_{\max ;k,l}^\omega (x,y)}}{{i_{\min ;k,l}^\omega (x,y)}}} } } \right|$ (22)

式中，数字图像按像素被划分成k₁×k₂个单位块，其二维序号为（k，l）， $i_{\max ;k,l}^\omega (x,y)$和 $i_{\min ;k,l}^\omega (x,y)$为序号为（k，l）的ω块中的最大灰度值和最小灰度值。通过迭代搜索来求取EME最大化的目标反射光偏振度p_obj，继而得到图像质量最优的去雾图像。

3 实验结果与分析

在实际的水下场景中进行成像实验，验证算法的可行性，实验装置如图3所示。光源为连续光谱LED光源，光束发散角为35°～37°，波长范围为435 nm～656 nm，照明装置内置线偏振片。成像系统是Teledyne DALSA公司的Genie Nano M2450偏振相机，相机的传感器是Sony公司的IMX250MZR CMOS单色传感器，该相机可输出一张同时包含4个偏振方向（0°；45°；90°；135°）的灰度图像，输出图像的分辨率为2 464×2 056像素。实验中，在PMMA水箱中将清水和牛奶混合，使水浑浊以达到改变水体介质的目的。实验目标为目标靶（材料如表1所示）和钢索，如图4所示。

图 3. 实验装置示意图（俯视视角）

Fig. 3. Schematic diagram of experimental device (top view)

下载图片 查看所有图片

图 4. 实验目标

Fig. 4. Experimental target

下载图片 查看所有图片

将PMMA水槽（1 180 mm×580 mm×310 mm）注满清水，在清水中加入1 000 μL染料墨水和4 000 μL全脂牛奶改变水体介质。实验过程中，主动光为唯一光源，避免其他杂散光干扰。调整光源和相机，使得相机0°偏振传感器方向获得的结果“最亮”。拍摄实验目标，成像系统获得一张含有正交偏振关系的原始合成图像，如图5所示。子图像中，偏振方向0°的图像为光强度最大的图像I_max，偏振方向90°的图像为光强度最小的图像I_min。在该水体条件下，原始图像的对比度低，图像中目标靶及钢索细节处无法有效识别。

图 5. 原始图像

Fig. 5. Original image

下载图片 查看所有图片

在计算机上进行图像处理验证本文算法。从原始图像中得到图像I_min和图像I_max（分辨率1 232×1 028像素），并对两正交偏振图像进行预处理，原偏振图像及预处理后的偏振图像如图6所示。根据本文算法，对偏振图像进行处理，得到水下去雾图像，如图7所示。与强度图像相比，处理后的图像中目标靶的轮廓信息更加明显，淹没的信息得以重现，图像对比度得到改善，表明该算法可以去除后向散射光对水下成像的影响，改善水下图像质量。

图 6. 原偏振图像和预处理后的偏振图像

Fig. 6. Original polarized image and preprocessed polarized image

下载图片 查看所有图片

图 7. 强度图像和处理后的图像

Fig. 7. Intensity image and processed image

下载图片 查看所有图片

通过二元多项式拟合出的后向散射光偏振度如图8(b)所示。x的阶数为2，y的阶数为3，决定系数R²为0.909 7。对比图8(a)和图8(b)可知，拟合函数依托于目标周围区域的后向散射光偏振度，估算出目标处的后向散射光偏振度，去除掉目标反射光的偏振特性影响，实现将后向散射光偏振度作为全局变量。目标反射光偏振度p_obj取值0.1，EME达到最大值，取值示例如图9所示。

图 8. 偏振度图和后向散射光偏振度拟合图

Fig. 8. Polarization degree diagram and backscattering light polarization degree fitting diagram

下载图片 查看所有图片

图 9. 目标反射光偏振度取值示例

Fig. 9. Example of target reflected light polarization degree value

下载图片 查看所有图片

利用图像增强测量值EME和直方图对处理后的图像进行评估，将本文算法与其他算法处理的水下图像进行对比分析，如图10所示。其中，图10（a）为强度图像，图像受后向散射光的影响，目标与背景之间对比度下降，图像质量受到很大影响；图10（b）为Schechner方法处理得到的图像，该方法在非均匀主动光照明条件下适用性降低，由于将后向散射光特征参量作为常量处理，使得处理后的图像中出现失真现象；图10（c）为Treibitz方法处理得到的图像，该方法建立偏振主动成像模型，手动取一块背景区域作后向散射光特征参量，处理后的图像出现了明暗不均的现象；图10（d）为Li方法处理得到的图像，该方法为经典Schechner方法的延展，预处理虽然提升了图像对比度，但也使后向散射光特征参量增大，处理后的图像中失真现象更加明显；图10（e）为本文算法处理得到的图像，将后向散射光偏振度作为全局变量处理，可以有效地解决复原图像的明暗不均及失真现象。从EME值来看，图10（e）的EME值较以往算法约提高70%；从直方图来看，图10（e）的直方图灰度值分布范围更宽广且相对均匀，意味着图像的对比度更大；从主观视觉来看，图10（e）克服了明暗不均和失真的现象，图像质量最优。通过对比分析，可看出本文算法和以往算法相比有一定的提高和改善。

图 10. 各算法处理的图像及其直方图对比

Fig. 10. Comparison of images processed by each algorithm and histograms

下载图片 查看所有图片

本文还采用两种常用的客观评价参数对实验图像进行比较，如表2所示。对比结果表明，本文算法处理的图像依然具有优势，与主观评价和EME参数评价结论相同。

为证明本文算法在不同浓度介质中的普遍适用性，在相同实验条件下加入1 000 μL染料墨水和6 000 μL全脂牛奶。如图11所示，与上述实验相比，由于介质浓度增大，强度图像中后向散射光的影响也随之增加，经过本文算法处理的图像表现出更高的对比度，图像质量优于强度图像，证明本文算法适用于不同浓度的介质，图像处理效果稳定。

图 11. 强度图像和处理后的图像

Fig. 11. Intensity image and processed image

下载图片 查看所有图片

4 结论

本文在发散型主动光照明的实验条件下，通过建立水下主动偏振成像模型，提出一种基于二元多项式拟合函数的后向散射光参数估计算法，并有效结合图像增强技术作为预处理，设计实验进行了算法的验证。实验结果表明，本文算法可以改善后向散射光对水下图像质量的影响，提高图像对比度，图像复原效果较以往算法相比有一定的提高，几种客观评价参数均有优势，直方图灰度值分布更广，且对于不同浓度的介质，算法有良好的复原稳定性。与此同时，本文算法受限于光源视场范围和图像噪声，图像噪声与后向散射光偏振度和目标反射光偏振度紧密相关，复原图像的直方图中像素点曲线波动变大，说明在处理的过程中噪声也被放大，这些问题还有待进一步讨论。