随着人类日益增多的海洋活动，水下成像技术已经成为人类了解海洋的一种重要手段。水下成像技术面临的主要挑战是如何克服水体对于光线的吸收、散射等作用产生的独特的噪声对于图像质量的影响。常用的水下成像技术包括普通光源成像，水下激光成像。由于水中介质对自然光的衰减作用较强，普通光源成像的距离较短，形成的水下图像的质量较差。采用具有较好的通透性的激光光源能够在一定程度上抑制介质对于光线的衰减以及后向散射的作用，水下图像的传输距离、成像精度、噪声抑制能力均优于普通光源成像。

虽然激光水下成像系统能够在一定程度上降低了水下介质对于散射作用，但是图像仍然受到由于后向散射所引起的散斑噪声的影响^[1]。又由于像增强器的非线性特征在一定程度上改变了散斑噪声的统计特性，因此激光水下图像的噪声效果受到了像增强器非线性响应的限制。本文研究了像增强器的非线性响应函数(camera response function,CRF)以及逆响应函数(inverse camera response function,ICRF)^[2]对于去噪算法的影响；在此基础上提出了一种对于光晕具有鲁棒性的辐射标定算法用于提高水下激光图像的噪声抑制效果。

1 研究背景 1.1 像增强器的非线性特性

像增强器是光路接收器的核心器件，主要由光阴极、微通道板(microchannel plate,MCP)、荧光屏、光纤、CCD组成。理想的像增强器是线性器件，然而在实际环境中，MCP、光纤都会对输入信号产生一定程度的非线性影响，从而导致像增强器的ICRF具有非线性特性。实际情况和理想状态下，像增强器的响应如图 1所示。

1.2 非线性响应对于抑制散斑噪声的影响

激光水下图像主要受到与光路同向的后向散射的影响，后向散射在数字图像中表现为与入射激光强度有关的散斑噪声。散斑噪声的统计特性易受到非线性响应的影响，从而影响噪声的抑制效果。令原始的激光图像是x_i，u_i的均值为u，方差是${\tilde u}$的白噪声，g是非线性响应函数，z_i是经过g作用后的图像。散斑噪声的表达式^[3]为

$ g\left( {{z_i}} \right) = {x_i}u $

(1)

对式(1)两边同时取期望，由于散斑噪声与有用信号统计独立，则x_i的局部均值为

$ {{\bar x}_i} = \frac{{\bar g\left( {{z_i}} \right)}}{{\bar u}} $

(2)

对式(1)两边同时乘方再取期望，根据x_i的局部均值表达式得到x_i的局部方差：

$ {{\tilde x}_i} = \frac{{{{\tilde g}^2}\left( {{z_i}} \right) + {{\bar g}^2}\left( {{z_i}} \right)}}{{\tilde u + {{\bar u}^2}}} - \frac{{{{\bar g}^2}\left( {{z_i}} \right)}}{{{{\bar u}^2}}} $

(3)

得到x_i和$ {{\tilde x}_i} $后，根据Kalman滤波，则x_i的最优估计可以表示为以下迭代形式：

$ {{\hat x}_i} = {x_i} + \frac{{\bar u{{\tilde x}_i}\left[ {g\left( {{z_i}} \right) - \bar u{{\bar x}_i}} \right]}}{{{{\bar x}_i}^2\tilde u + {{\tilde x}_i}{{\bar u}^2}}} $

(4)

通过式(4)可知，如果直接对输出的灰度图像进行Kalman滤波，则图像复原结果与理论值存在一定的理论误差。

1.3 辐射度标定

传统的辐射标定方法通过查表以及不断改变入射光源辐射度，标定出入射光子以及系统响应之间的对应关系^[4-5]。由于制造工艺以及水下物质的衰减作用导致真实的辐照度与查表得到的理论值存在误差。另一种常用的辐射度标定方法是基于多曝光图像的辐射标定方法。这类方法的前提条件是相机位置需要保持绝对静止，如果相机出现微小的移位则会出现光晕现象。

2 标定算法原理 2.1 灰阶映射

设函数f是g的反函数，辐射标定的目标是求解g(0),g(1),...,g(255)。令l_i是灰阶i的辐射度，t_j是积分时间，y_i.j是灰阶i在j次积分的灰阶值：

$ g\left( {{y_{i,j}}} \right) = {l_i}{t_j}\; $

(5)

这里将l_i在积分时间t_j、t_j+1对应的灰阶y_i,j、y_i,j+1定义为灰阶映射。根据式(5)可知y_i,j、y_i.j+1满足以下约束关系：

$ t_j^{ - 1}g\left( {{y_{i,j}}} \right) = t_{j + 1}^{ - 1}g\left( {{y_{i,j + 1}}} \right) $

(6)

设函数H(·)表示概率分布函数，C(·)表示累积概率分布函数，则存在以下引理和定理。

引理1 对于任意严格单调函数h，等式Hh(l_i)=H(l_i)恒成立。

证明 因为h具有严格的单调性，则l_i与f(l_i)存在一一映射的关系。因此l_i与对应的f(l_i)出现的概率相等。证毕。

定理1 同一光照辐射度在相邻两个曝光图像中对应灰阶的累积概率分布相等，即C(y_i,j)=C(y_i,j+1)。

证明 在灰度序列中对y_i,j进行累积分布运算得到以下表达式：

$ C\left( {{y_{i,j}}} \right) = C\left( {f\left( {{l_i} \cdot {t_j}} \right)} \right)\; $

(7)

根据累积分布函数的定义，式(7)可以变换为以下形式：

$ C\left( {{y}_{i,j}} \right)=\sum{_{f\left( {{t}_{j}}\cdot {{l}_{k}} \right)\le {{y}_{i,j}}}H\left( f\left( {{l}_{i}}\cdot {{t}_{j}} \right) \right)} $

(8)

由于f是严格单调增加的函数，根据引理1，式(8)可以变换为以下形式：

$ C\left( {{y}_{i,j}} \right)=\sum{_{f\left( {{t}_{j}}\cdot {{l}_{k}} \right)\le {{y}_{i,j}}}H\left( {{t}_{j}}\cdot {{l}_{k}} \right)} $

(9)

设h(k)=t·k，显然h是一个严格单调增加的函数，根据引理1，式(9)可以转化为

$ C\left( {{y}_{i,j}} \right)=\sum{_{f\left( {{t}_{j}}\cdot {{l}_{k}} \right)\le {{y}_{i,j}}}H\left( {{l}_{k}} \right)} $

(10)

因为g具有单调增加性，则式(10)的取值范围转化为l_k≤t_j^-1·g(y_i,j。将l_k≤t_j^-1·g(y_i,j代入式(10)中得到以下表达式：

$ C\left( {{y}_{i,j+1}} \right)=\sum{_{{{l}_{k}}\le t_{j}^{-1}\cdot g\left( {{y}_{i,j}} \right)}H\left( {{l}_{k}} \right)} $

(11)

同理可得y_i,j+1的累积分布如式(12)所示：

$ C\left( {{y_{i,j + 1}}} \right) = \sum {_{{l_k} \le t_{j + 1}^{ - 1} \cdot g\left( {{y_{i,j + 1}}} \right)}H\left( {{l_k}} \right)} $

(12)

根据约束条件(6)可知式(11)与(12)的累加范围相等，因此C(y_i,j)=C(y_i,j+1)。证毕。

根据定理1可知，只要在j+1幅图像中找到对应的y_i,j+1，使y_i,j+1的累积分布与y_i,j相等，则y_i,j+1是y_i,j的灰阶映射。

2.2 递推表达式

假设积分时间的初始条件是y_i.1=i，根据灰阶映射满足等式(6)的约束关系，则初次积分与二次积分的灰阶映射如式(13)所示：

$ g\left( i \right) = g\left( {{y_{i,1}}} \right) = {t_1} \cdot t_2^{ - 1} \cdot g\left( {{y_{i,2}}} \right) $

(13)

同理g(y_i,2)如式(14)所示：

$ g\left( {{y_{i,2}}} \right) = {t_2} \cdot t_3^{ - 1} \cdot g\left( {{y_{i,3}}} \right) $

(14)

将f_CRF函数作用于式(14)的等号两侧

$ {y_{i,2}} = f\left( {{t_2} \cdot t_3^{ - 1} \cdot {f^{ - 1}}\left( {{y_{i,3}}} \right)} \right) $

(15)

将式(15)代入式(13)得到经过三次曝光之后灰阶i对应的辐射度如式(16)所示：

$ g\left( i \right) = {t_1} \cdot t_3^{ - 1} \cdot g\left( {{y_{i,3}}} \right) $

(16)

以此类推可以得到N次曝光后灰阶i对应的灰阶映射的通项公式：

$ g\left( i \right) = g\left( {{y_{i,1}}} \right) = {t_1} \cdot t_N^{ - 1} \cdot g\left( {{y_{i,N}}} \right)\; $

(17)

根据式(17)可知只要g(y_i,N)的表达式已知，则能够得到g(i)的标定结果。

2.3 积分时间增量

将多曝光图像按照积分时间从大到小的顺序进行排列，当积分时间足够大时存在以下关系y_i,1>y_i,2>…>y_i,j≈0。根据图 1可知当灰阶y_i,j足够小时，可以将ICRF近似为线性函数，所以这里将ICRF分解为非线性曲线和线性曲线两部分，利用递推的思想将响应强度较高的非线性部分转化为响应强度较低的线性部分。

假设$ {\tilde{y}} $_i,2是距离浮点型的灰阶y_i,2最近的整数型灰阶，t₂+△t_i,2是$ {\tilde{y}} $_i,2对应的积分时间；同理，$ {\tilde{y}} $_i,3是距离y_i,3最近的整型灰阶，t₃+△t_i,3是$ {\tilde{y}} $_i,3对应的积分时间。以此类推，可以得到一个时间序列△t_i,j以及整型灰阶序列$ {\tilde{y}} $_i,j。这里将△t_i,j定义为灰阶i在第j次积分时间对应的积分时间增量，则$ {\tilde{y}} $_i,j，t_j，△t_i,j满足如下函数关系：

$ {{{\tilde{y}}}_{i,j}}=f\left( \left( {{t}_{j}}+\vartriangle {{t}_{i,j}} \right)\cdot {{l}_{j}} \right) $

(18)

将等式(18)在(t_j+△t_i,j)·l_i=t_j·l_i进行一阶泰勒级数展开：

$ {{{\tilde{y}}}_{i,j}}=f\left( {{t}_{j}}\cdot {{l}_{i}} \right)+{f}'\left( {{t}_{j}}\cdot {{l}_{j}} \right)\cdot {{l}_{i}}\cdot \vartriangle {{t}_{i,j}} $

(19)

同理j+1时刻的表达式：

$ {{y}_{i,j+1}}=f\left( {{t}_{j}}\cdot {{l}_{i}} \right)+{f}'\left( {{t}_{j}}\cdot {{l}_{j}} \right)\left( {{t}_{j+1}}-{{t}_{j}} \right){{l}_{i}} $

(20)

联立式(19)、(20)得到△t_i,j的表达式：

$ \vartriangle {{t}_{i,j}}=\frac{\left( {{{\tilde{y}}}_{i,j}}-{{y}_{i,j}} \right)\cdot \left( {{t}_{j}}-{{t}_{j+1}} \right)}{{{y}_{i,j+1}}-{{y}_{i,j}}} $

(21)

通过式(21)可以得到灰阶i对应的积分时间增量序列，即△t_i,1,△t_i,2,...,△t_i,N。根据式(17)得到引入时间增量后ICRF表达式：

$ g\left( i \right)=g\left( {{y}_{i,1}} \right)=\frac{{{t}_{1}}\cdot g\left( {{{\tilde{y}}}_{i,N}} \right)}{{{t}_{N}}+\vartriangle {{t}_{i,N}}} $

(22)

2.4 线性部分标定

因为ICRF曲线端点的灰阶值分别是0和1，所以采用0和255作为参考点，通过线性插值对ICRF的线性部分进行估计。设灰阶i与0距离最近的灰度映射是y_i,mi，m_i是y_i,mi的曝光次数，则ICRF为

$ \begin{matrix} g\left( {{{\tilde{y}}}_{i,{{m}_{i}}}} \right)=\frac{g\left( {{{\tilde{y}}}_{255,{{m}_{255}}}} \right)-g\left( {{{\tilde{y}}}_{0,{{m}_{0}}}} \right)}{{{{\tilde{y}}}_{255,{{m}_{255}}}}-{{{\tilde{y}}}_{0,{{m}_{0}}}}}\times \\ \left( {{{\tilde{y}}}_{i,{{m}_{i}}}}-{{{\tilde{y}}}_{0,{{m}_{0}}}} \right)+{{{\tilde{y}}}_{0,{{g}_{0}}}} \\ \end{matrix} $

(23)

根据先验知识可知$ {\tilde{y}} $_0,m0=0，g($ {\tilde{y}} $_0,m0)=0，因此只有g($ {\tilde{y}} $_255,m255)的值未知。由于灰阶255对应值为1，通过式(22)可以得到灰阶255对应的ICRF的值可以表示为以下形式：

$ g\left( 255 \right)=1=\frac{{{t}_{1}}\cdot g\left( {{{\tilde{y}}}_{255,{{m}_{255}}}} \right)}{{{t}_{{{m}_{255}}}}+\vartriangle {{t}_{255,{{m}_{255}}}}}\ $

(24)

整理式(24)得到g($ {\tilde{y}} $_255,m255)的值并将g($ {\tilde{y}} $_255,m255)以及ICRF的已知条件代入式(23)中得到灰阶i线性部分的ICRF表达式：

$ g\left( {{{\tilde{y}}}_{i,{{m}_{i}}}} \right)=\frac{{{{\tilde{y}}}_{i,{{m}_{i}}}}\cdot \left( {{t}_{{{m}_{255}}}}+\vartriangle {{t}_{255,{{m}_{255}}}} \right)}{{{t}_{1}}\cdot {{{\tilde{y}}}_{255,{{m}_{255}}}}} $

(25)

将式(25)代入式(23)中，将ICRF的递推关系重新写为以下形式：

$ g\left( i \right)=\frac{{{{\tilde{y}}}_{i,{{m}_{i}}}}\cdot \left( {{t}_{{{m}_{255}}}}+\vartriangle {{t}_{255,{{m}_{255}}}} \right)}{{{{\tilde{y}}}_{255,{{m}_{255}}}}\left( {{t}_{{{m}_{i}}}}+\vartriangle {{t}_{i,{{m}_{i}}}} \right)},\left( 0,255 \right) $

(26)

由于式(26)中的$ {\tilde{y}} $_i,mi和$ {\tilde{y}} $_255,m255可以通过式(22)计算出来，根据式(21)可以计算曝光增量△t_255,m255和△t_i,mi，因此通过式(26)可以得到任意灰阶值对应的光照辐射度。

3 实验结果与数据分析

实验部分通过脉宽为50 ns的激光器和最小门宽为15 ns的增强型CCD构成距离选通激光水下成像系统。

3.1 辐射标定噪声抑制效果的影响

实验部分首先对图像Lena进行散斑噪声进行污染；通过如图 1所示的ICRF进行非线性变换；通过对比辐射标定前、后，恢复图像的主观效果以及峰值信噪比(peak signal noise ratio,PSNR)、结构相似度(structural similarity index measurement,SSIM)^[6]，验证了本文算法的有效性。图 2~4分别是通过Sarode算法^[3]、Kaur算法^[7]、Hoque算法^[8]获得的图像恢复结果。

从图 2~4可以看出同种噪声抑制算法对于辐射标定后的散斑噪声抑制效果好于未加入辐射标定的对应的恢复图像。相比于未加入辐射标定的图像复原结果，辐射标定后对应的恢复图像局部的颗粒效果以及灰度的起伏得到明显的抑制，图像整体的噪声抑制能力得到显著的提高。表 1是图 2~4对应的SSIM以及PSNR，从表 1可知经过辐射标定后，SSIM以及PSNR均高于标定前，说明辐射标定有助于改善散斑噪声的抑制效果，加入辐射标定后图像的复原结果与原始图像的相似度较高。

3.2 辐射标定对于噪声抑制算法的影响

由于多曝光图像对应同一场景，因此线性图像亮度分布规律需要尽可能的保持一致。本文采用线性多曝光图像的亮度方差来测量辐射标定方法对于光晕的抑制程度，线性多曝光图像的方差为

$ \gamma =\frac{1}{N}\sum\limits_{j=1}^{N}{\left[ \frac{1}{M\cdot {{t}_{j}}}\sum\limits_{i=1}^{M}{{{g}^{2}}\left( {{z}_{i,j}} \right)-{{\varphi }^{2}}} \right]}\ $

(27)

式中：φ是线性多曝光图像的均值；γ描述了图像的亮度分布与亮度均值的偏移程度，如果γ的值越小，说明辐射标定后图像的亮度分布距离亮度均值较近，标定函数抑制光晕的效果越好，反之亦然。图 5是通过调整像增强器的积分时间得到的多曝光图像。

图 6、7分别是通过Sarode^[3]算法得到的立方体、椭圆形的恢复图像。图 8是通过线性方法、Debevec算法^[9]、Granados算法^[10]、本文算法标定的ICRF。

表 2是4种ICRF标定方法对应的γ值，可以看出经过辐射标定后，本文算法对应的亮度的方差小于其他辐射标定法，说明辐射标定后图像序列能够最大程度地保证同一场景的亮度分布规律，因此本文算法对于光晕具有较高的鲁棒性。

4 结束语

本文研究了像增强器的非线性响应对于散斑噪声抑制算法的影响，提出一种对于光晕具有鲁棒性的辐射标定算法。通过实验，验证了在噪声抑制算法之前先对增强器进行辐射标定有助于提高激光水下图像的散斑噪声抑制效果。