高动态范围图像融合过程中的噪声抑制算法

1 引言

高动态范围(high dynamic range, HDR)图像能准确展现真实场景的动态范围，因此受到了人们的广泛关注。HDR图像是通过合成多幅不同曝光的低动态范围(low dynamic range, LDR)图像得到，它克服了传统LDR图像动态范围小的限制。目前，获取HDR图像的方法主要有两类：第一类方法是通过估计相机响应曲线(camera response function, CRF)将LDR图像映射到辐照域^[1-4]，并在辐照域中进行融合，从而得到一幅动态范围宽的图像。但该方法得到的HDR图像不能在常规图像设备上显示，还需要进行色调映射^[5-7]。第二类方法则是在图像域直接融合多幅不同曝光的LDR图像，从而产生一幅具有HDR效果的LDR图像。由于该方法无需进行CRF估计和色调映射，显得十分简单有效，因此受到了广泛关注。Ma等人^[8]利用结构块分解来融合多曝光图像；Li等人^[9]利用引导滤波，将图像分解为基本层和细节层，在不同层中分别融合图像；Liu等人^[10]利用密集的尺度不变特征变换融合多曝光图像；Mertens等人^[11]利用对比度、饱和度和良好曝光度三个参数构造权重值，融合多曝光图像。

但上述融合方法均没有考虑图像中的噪声。实际中，在低光环境下拍摄多曝光图像序列时，为防止图像模糊，相机会设置为高灵敏度，会产生噪声，且在图像融合后会进一步放大。目前已有的HDR图像噪声抑制方法是在图像融合前，对多曝光图像序列进行预处理。Akyuz等人^[12]根据图像噪声水平与曝光时间成反比的关系，利用高曝光图像去抑制低曝光图像的噪声，但相应高曝光图像的噪声就无法被抑制；孙等人^[13]提出一种基于图像辐照值与曝光时间成线性关系的曲线拟合方法，可以抑制HDR图像中的高斯噪声，但噪声抑制效果受限于CRF的准确标定；Min等人^[14]认为时域上的加权平均不足以降低传感器噪声，提出利用图像的空间相邻和时间相邻信息，对多曝光图像序列使用结构自适应时空平滑滤波器抑制噪声，取得了较好的效果。但他们的方法均需要将LDR图像映射到辐照域进行处理，即需要已知图像序列的曝光信息。在实际中，多数图像的曝光信息很难获得，因此Wei等人^[15-16]利用强度映射函数(intensity mapping function)将Akyuz等人的方法扩展到图像域，其效果与Akyuz等人的方法类似。Lee等人^[17]提出一种HDR图像融合中的噪声抑制方法，首先利用拉普拉斯金字塔分解多曝光图像序列，再利用高斯金字塔融合各分解的子带图像，对融合后的低频子带图像利用双边滤波去噪，对高频子带图像利用小波软阈值滤波去噪，去噪参数设置没有考虑子带图像不同区域的噪声水平，对图像整体设置为常数，进行一致性滤波，导致对噪声水平低的图像区域造成细节损失，对噪声水平高的图像区域去噪不完全。为了有效抑制图像噪声，且较好保留图像细节，Liu等人^[18]提出一种噪声水平估计算法，与现有的噪声抑制算法^[19-21]结合，可以有效提高噪声抑制效果。Chen等人^[22]提出利用神经网络对图像去噪，但特定的网络只针对特定水平的噪声，对于多曝光图像序列来说，噪声水平不一致，很难用单个网络来实现噪声抑制。Zhang等人^[23]提出利用残差学习来克服传统神经网络只针对特定水平的噪声的缺点，但仍假设单幅图像所有区域的噪声水平相同。由于多曝光图像融合，需要保留低曝光图像中高亮区域的细节信息以及高曝光图像中低暗区域的细节信息，且图像噪声与亮度有关，因此直接对整幅图像估计噪声水平，得到的是噪声水平均值，以此指导图像滤波，无法保留多曝光图像序列不同亮度区域的的细节信息，会导致融合图像的视觉质量降低。

针对上述分析，本文提出一种基于亮度分区、噪声水平估计的高动态范围图像融合过程中的噪声抑制算法，也适用于多曝光图像融合。新方法主要从两方面提高了噪声抑制算法的性能：一方面通过亮度分区，保留多曝光图像不同亮度区域的细节信息，避免对图像整体处理导致图像模糊；另一方面通过估计噪声水平，指导噪声抑制算法，相比于传统的盲去噪具有更好的去噪效果。

2 高动态范围图像融合过程中的噪声抑制算法

现有的HDR图像融合过程中的噪声抑制算法没有考虑图像在不同亮度区域的噪声水平，容易出现细节损失、模糊等不好的视觉现象。因此，本文提出了一种基于亮度分区、噪声水平估计的高动态范围图像融合过程中的噪声抑制算法，可以在有效抑制图像低暗区域噪声的同时，保留图像的细节信息。

首先对输入的LDR图像进行自适应亮度分区，根据图像的亮度中值，将原图像分解为三幅不同亮度范围的图像，这可以有效保留所有曝光图像的中间亮度细节信息，以及低曝光图像的高亮度区域细节信息和高曝光图像的低亮度区域细节信息，再将三幅图像通过块滤波器分解为重叠图像块；其次对每个重叠图像块利用结构张量提取弱纹理区域，对该区域利用图像的协方差矩阵估计噪声水平，以此准确指导图像去噪，不会对图像造成过度平滑；由于每个图像块估计的噪声水平可靠性不一致，且噪声水平与亮度有关，故对低亮度区域图像块采用“选择最大”策略，对中间亮度和高亮度区域图像块采用“平均”策略，得到最终的噪声水平；然后根据得到的噪声水平，指导图像的稀疏去噪；对去噪后的三幅不同亮度图像进行重构；最后将处理后的多幅LDR图像融合成HDR图像，如图 1所示。

[Image omitted: See PDF]

2.1 图像亮度分区

HDR图像具有比LDR图像更宽的亮度范围，其特点是图像对比度大，保留了多曝光图像序列的信息，即低曝光图像的高亮区域信息，高曝光图像的低暗区域信息以及所有曝光图像的中间亮度区域信息，因此需要对多曝光图像序列进行亮度分区，避免整体处理造成细节损失。将彩色图像转换颜色空间，从RGB转到XYZ空间，其中Y是亮度分量，计算亮度中值，利用线性归一化函数，将亮度中值归一化至0~1之间，作为图像中间亮度值T_h。根据T_h值对图像进行亮度分区，如下：

\( {L_{{\rm{dark}}}} = {L_{{\rm{max}}}}-\;(0.9 + 0.1\;{T_{\rm{h}}})\;({L_{{\rm{max}}}}-{L_{{\rm{min}}}}), \)(1) \( {L_{{\rm{bright}}}} = {L_{{\rm{min}}}}\;\; + \;[\;0.6 + 0.4\;(1-{T_{\rm{h}}})]\;({L_{{\rm{max}}}} -{L_{{\rm{min}}}}\;), \)(2)

其中： \({L_{{\rm{dark}}}}\) 为图像低亮度区域的阈值，称为暗阈值； \({L_{{\rm{bright}}}}\) 为图像高亮度区域的阈值，称为亮阈值； \({L_{{\rm{max}}}}\) 和 \({L_{{\rm{min}}}}\) 分别是图像灰度化并归一化到0~1之间的最大值和最小值。若图像亮度小于暗阈值，则认为是低亮度区域；若图像亮度大于亮阈值，则认为是高亮度区域；介于两者之间，则认为是中间亮度区域。

作为一个实例，图 2给出了对多曝光图像序列进行亮度分区的结果，图 2(a)是7个不同曝光得到的原始LDR图像，图 2(b)是分割得到的对应的低亮度图像，图 2(c)是对应的中间亮度图像，图 2(d)是对应的高亮度图像；图中，白色表示属于该亮度区域，黑色表示不属于该亮度区域。由实验结果分析，可以发现一个现象：亮度分区在图像局部存在相似性，即对于单个亮度区间，若某像素属于该区间，则存在一定尺寸的邻域也属于该区间。利用该现象可以对图像不同亮度区域估计噪声水平。

[Image omitted: See PDF]

2.2 对不同亮度区域估计噪声水平

噪声水平是指反映信号中噪声强度的参数，它能够被应用于指导噪声抑制算法的设计，准确估计的噪声水平可以有效提高噪声抑制算法的性能，在抑制图像噪声的同时，保留图像细节信息。

图像的噪声和纹理细节都属于图像的高频信息，直接对原图像估计噪声水平会导致估计结果偏高，从而使处理后的图像纹理信息丢失，因此本文利用结构张量估计图像纹理强度，提取图像的弱纹理区域估计噪声水平。假设噪声为加性高斯白噪声，其噪声模型如下：

\( {y_i} = {x_i} + n\;\;, \)(3)

其中： \({y_i}\) 是图像局部区域的观测像素值， \({x_i}\) 是图像局部区域的真实像素值，n是均值为0、方差为 \(\sigma _n^2\) 的高斯噪声。将图像数据投影到最小方差方向上，由于噪声与信号无关，所以投影后噪声方差不变。最小方差方向根据PCA计算可得，由协方差矩阵的最小特征值所对应的特征向量确定，而方差等于最小特征值，即：

\( V = {\lambda _{\min }}\;(\frac{1}{M}\sum\limits_{i = 1}^M {{\mathit{\boldsymbol{y}}_i}\mathit{\boldsymbol{y}}_i^{\rm{T}}} ) = {\lambda _{\min }}\;(\frac{1}{M}\sum\limits_{i = 1}^M {{\mathit{\boldsymbol{x}}_i}\mathit{\boldsymbol{x}}_i^{\rm{T}}} ) + \sigma _n^2, \)(4)

其中：M是图像局部区域的总数。 \({\lambda _{\min }}\) 表示矩阵的最小特征值。由于高斯噪声在各个方向上具有相同的强度，即所有特征值都相同，所以估计的噪声水平为

\( \sigma _n^2 = {\lambda _{\min }}\;(\frac{1}{M}\sum\limits_{i = 1}^M {{\mathit{\boldsymbol{y}}_i}\mathit{\boldsymbol{y}}_i^{\rm{T}}} )-{\lambda _{\min }}\;(\frac{1}{M}\sum\limits_{i = 1}^M {{\mathit{\boldsymbol{x}}_i}\mathit{\boldsymbol{x}}_i^{\rm{T}}} ), \)(5)

结构张量定义为图像梯度的协方差矩阵，如下：

\( \begin{array}{l} \mathit{\boldsymbol{E}} = {\left[{\begin{array}{*{20}{c}} {{\partial _{\rm{h}}}}&{{\partial _{\rm{v}}}} \end{array}} \right]^{\rm{T}}}\left[{\begin{array}{*{20}{c}} {{\partial _{\rm{h}}}}&{{\partial _{\rm{v}}}} \end{array}} \right]\\ \;\; = \left[{\begin{array}{*{20}{c}} {\partial _{\rm{h}}^2}&{{\partial _{\rm{h}}}{\partial _{\rm{v}}}}\\ {{\partial _{\rm{v}}}{\partial _{\rm{h}}}}&{\partial _{\rm{v}}^2} \end{array}} \right], \end{array} \)(6)

其中： \({\partial _{\rm{h}}}\) 和 \({\partial _{\rm{v}}}\) 分别表示图像的水平梯度和垂直梯度。该矩阵的迹可以衡量图像的纹理强度，迹越大，表示纹理越强，对应于图像局部区域 \({y_i}\) 的纹理强度 \({H_i}\) 为

\( {H_i} = {\rm{tr}}(\mathit{\boldsymbol{E}}) \)(7)

其中：tr表示矩阵的迹。定义阈值τ，用于提取图像的弱纹理区域：

\( {p_i} = \left\{ {\begin{array}{*{20}{c}} {0, \;\quad {\rm{if}}\;{H_i} \le \tau }\\ {1, \;\quad {\rm{otherwise}}} \end{array}} \right.。\)(8)

当p_i=0，则表示该图像局部区域属于弱纹理区域，p_i=1表示该图像局部区域属于强纹理区域。对于图像的弱纹理区域，由于纹理信息很少，其原始无噪图像的协方差矩阵的特征值可以假定为零，因此式(5)改为

\( \sigma _n^2 = {\lambda _{\min }}\;(\frac{1}{N}\sum\limits_{i = 1}^N {{\mathit{\boldsymbol{y}}_i}\mathit{\boldsymbol{y}}_i^{\rm{T}}} ), \)(9)

其中： \(\sigma _n^2\) 是估计的噪声水平，N是图像弱纹理区域的总数， \(N \le M\) 。由于不同图像纹理强度有差异，直接固定阈值τ会导致提取方法不能自适应图像内容，所以采用迭代方式。具体过程如下：初始阈值直接取原始图像的迹，即初始的弱纹理区域是原图，再进行噪声水平估计，用新得到的噪声水平修改阈值。在本文中，经过大量实验，迭代次数取3，可以准确地估计出噪声水平。

由于需要分别对图像不同亮度区域估计噪声水平，根据2.1节中的观察：亮度分区在图像局部存在相似性，即对于单个亮度区间，若某像素属于该区间，则存在一定尺寸的邻域也属于该区间，如图 2所示，因此只要将噪声图像分解为合适尺寸的重叠块，使该重叠块的区域均属于低亮度区域，则对该重叠块估计的噪声水平可作为低亮度区域的噪声水平。同理，中间亮度和高亮度区域的噪声水平可以同样测得。定义重叠块尺寸为w×w，非重叠区域尺寸为d，给定一幅图像，尺寸为m×n，则分解的重叠块数目为

\( {N_{{\rm{over-b}}}} = floor\;(\;\frac{{m-w}}{d}{\rm{ + 1}}\;) \times floor\;(\;\frac{{n-w}}{d}{\rm{ + 1}}\;), \)(10)

其中floor()表示向下取整。由于要保证每个重叠块尽可能属于单个亮度区间，所以非重叠区域不能过大，但过小的非重叠区域会产生冗余的重叠块，导致算法效率降低。根据大量的实验经验，在这里，w=100，d=20可以准确地分解图像的重叠块。对分解的所有重叠块检测亮度区间并估计噪声水平。其中，每个重叠块存在三个亮度区间，计算三个亮度区间在重叠块中的比例：

\( R(k, i) = \frac{{{S_{k, i}}}}{{{S_{k, {\rm{sum}}}}}}, \)(11)

其中：k表示重叠块索引，即第k个重叠块；i=1, 2, 3，表示三个亮度区间，定义i=1表示低亮度区间；i=2表示中间亮度区间；i=3表示高亮度区间。S_{k, i}表示第k个重叠块中属于i亮度区间的像素点个数，S_{k, sum}表示该重叠块的像素点总数。由于分解规则是针对图像空间的，即每个重叠块均有可能包含多个亮度区间的像素，但像素比例不同，因此为区分重叠块属于某个亮度区间，定义阈值 \({T_{\rm{s}}}\) ，当 \(R(k, i) \ge {T_{\rm{s}}}\) 时，表示该重叠块中大部分像素都属于i亮度区间，即认为该重叠块属于i亮度区间。为保证图像的每个亮度区间均有对应的多个重叠块，阈值 \({T_{\rm{s}}}\) 不能取1，经过大量的实验，在这里取值0.9，可以准确地区分重叠块。

由于噪声主要存在于图像的低亮度区域，即低亮度区域的噪声水平高，因此对该区域采用“选择最大”策略估计噪声水平。首先找出对应R(k, 1)大于 \({T_{\rm{s}}}\) 的所有重叠块；然后根据噪声水平对选出的重叠块排序；最后选出重叠块的噪声水平最大值作为最终估计的噪声水平。对图像中间亮度区域和高亮度区域采用“平均”策略，可以解决由低亮度区域过多而造成的噪声水平估计过高的问题，如下：

\( \left\{ \begin{array}{l} {{\hat \delta }_{n1}} = \max \;({\delta _k}){\rm{, }}\;\;{\rm{if}}\;R(k, 1) \ge {T_{\rm{s}}}\\ {{\hat \delta }_{n2}} = \frac{1}{P}\sum\limits_{} {{\delta _k}} {\rm{,...