2.1   张量分解和特征提取

由于现实数据的多维性，使用向量描述显然不能很好反映它们复杂的结构和本质，因此Leopold Kronecker提出了张量的概念。一个$N$阶张量$\mathcal{A}$可以如下表述：

其中：${l_1}, {l_2}, ..., {l_N}$表示每个维度上的元素个数。张量分解作为一种处理高维数据的有效手段，已经在数据分类、图像水印和图像恢复等众多计算机视觉领域上取得了显著效果。张量分解有两种基本形式：CP分解和Tucker分解^[14]。本文将Tucker分解应用到多曝光图像融合领域。Tucker分解是将一个N阶张量$\mathcal{A} \in {\mathfrak{R} ^{{l_1} \times {l_2} \times ... \times {l_N}}}$分解成一个核心张量$\mathcal{B} \in {\mathfrak{R} ^{{l_1} \times {l_2} \times ... \times {l_N}}}$与N个矩阵${{\mathit{\boldsymbol{U}}}^{(n)}} \in {{\mathit{\boldsymbol{R}}}^{{l_n} \times {R_n}}}, (n = 1, 2, ..., N)$的$n$模式乘积形式，即：

其中：${ \times _n}(n = 1, ..., N)$表示张量与矩阵的$n$模式乘积，即张量$\mathcal{B}$沿着第$n$维展开成对应矩阵形式后与${{\mathit{\boldsymbol{U}}}^{(n)}}$的乘积。当${{\mathit{\boldsymbol{U}}}^{(n)}}$为正交矩阵时，N阶张量$\mathcal{A}$与核心张量$\mathcal{B}$存在如下的转换关系：

多曝光图像融合的难点在于融合图像的边缘和纹理等细节特征难以保留，而Tucker分解后的核心张量较大程度地保留了源图像中的所有信息和特征，因此更利于进行特征提取。

首先每幅RGB源图像${\mathit{\boldsymbol{I}}_k}(k = 1, ..., K)$可以表示为三阶张量$\chi \in {\mathfrak{R} ^{{l_1} \times {l_2} \times 3}}$，${l_1}$和${l_2}$分别为图像的宽和高。然后将该三阶张量分别进行模式1、模式2和模式3展开，得到对应的矩阵形式

再对上述矩阵进行SVD分解，得到正交矩阵${{\mathit{\boldsymbol{U}}}^{(1)}}$, ${{\mathit{\boldsymbol{U}}}^{(2)}}$, ${{\mathit{\boldsymbol{U}}}^{(3)}}$，大小分别为${l_1} \times {l_1}$，${l_2} \times {l_2}$和$3 \times 3$。

其中：${\text{SVD}}( \cdot )$表示矩阵的奇异值分解。最后根据式(2)得到对应的核心张量$C \in {\mathfrak{R} ^{{l_1} \times {l_2} \times 3}}$，用于进行特征提取。又由于核心张量$\mathcal{B}$的子张量满足有序性^[15]，即对任意$1 \leqslant n \leqslant N$，有${\mathcal{B}_{{i_n} = 1}} \geqslant {\mathcal{B}_{{i_n} = 2}} \geqslant ... \geqslant {\mathcal{B}_{{i_n} = {l_N}}} \geqslant 0$。所以对于每幅RGB图像的核心张量$C \in {\mathfrak{R} ^{{l_1} \times {l_2} \times 3}}$，满足：

在这里，定义${\mathit{\boldsymbol{C}}_1}$为子带一，${\mathit{\boldsymbol{C}}_2}$为子带二，${\mathit{\boldsymbol{C}}_3}$为子带三。根据核心张量的有序性，${\mathit{\boldsymbol{C}}_1}$携带的信息和能量最多，故在第一子带图上提取特征。又由于人类视觉系统对边缘特征比较敏感，因此采用拉普拉斯算子在第一子带图上提取边缘信息，从而得到${\mathit{\boldsymbol{C}}_1}$的高频分量$\mathit{\boldsymbol{H}}$：

其中：$\mathit{\boldsymbol{L}} = \left[ {0, 1, 0;1, - 4, 1;0, 1, 0} \right]$，该算子相比于一阶微分算子具有更好的细节保留能力，$\mathit{\boldsymbol{C}}_1^k$表示第$k$幅源图像的第一子带图，${\mathit{\boldsymbol{H}}_k}$表示第$k$幅第一子带图上提取的边缘特征。图 2为House序列的多曝光图像组^[9]，图 3为图 2中图像经Tucker分解得到的核心张量图，图 4为在图 3所示的第一子带图上提取的边缘特征。

图 2.  House序列的多曝光图像组

Figure 2.  Multi-exposure image stack of House sequence

下载: 全尺寸图片 幻灯片

图 3.  House序列的核心张量。

(a)第一子带图；(b)第二子带图；(c)第三子带图

Figure 3.  The core tensor of House sequence.

(a) The first sub-band; (b) The second sub-band; (c) The third sub-band

下载: 全尺寸图片 幻灯片

图 4.  第一子带的边缘特征

Figure 4.  The edge feature of the first sub-band

下载: 全尺寸图片 幻灯片

2.2   卷积稀疏表示和字典学习

当给定一个原始信号${\mathit{\boldsymbol{S}}} \in {{\mathit{\boldsymbol{R}}}^n}$和一个过完备字典矩阵${\mathit{\boldsymbol{D}}} \in {{\mathit{\boldsymbol{R}}}^{n \times m}}(n < m)$，SR的目标是为了估计一个只有几个非零项的稀疏向量${\mathit{\boldsymbol{x}}} \in {{\mathit{\boldsymbol{R}}}^m}$，例如${\mathit{\boldsymbol{S}}} \approx {\mathit{\boldsymbol{Dx}}}$。因此大多数稀疏编码问题可以表述为以下优化问题^[6]：

其中：$|| \cdot |{|_1}$表示L1范数，用来计算矩阵中非零项的数目。${\mathit{\boldsymbol{D}}}$表示训练好的字典矩阵，${\mathit{\boldsymbol{x}}}$表示稀疏系数，${\mathit{\boldsymbol{S}}}$表示原始图像，$\lambda$表示正则化参数。

传统的SR方法是通过对一组重叠块进行稀疏表示，其结果往往是多值的且没有对整幅图像进行优化。所以Zeiler^[16]提出用带卷积操作的稀疏编码来代替传统的稀疏表示，即利用字典滤波器与相应特征响应的卷积形式来替代字典矩阵与稀疏编码系数的乘积形式，该表示是单值的且对整幅图像进行了联合优化。类似于式(7)中的基追踪去噪^[17](basis pursuit denoising, BPDN)问题，相应的卷积基追踪去噪(convolutional basis pursuit denoising, CBPDN)问题定义如下：

其中：$\{ {{\mathit{\boldsymbol{d}}}_m}\}$是由$M$个滤波器组成的卷积字典，$\otimes$表示卷积操作，$\{ {{\mathit{\boldsymbol{x}}}_m}\}$是$M$维稀疏系数图。

目前，已经提出了很多算法来解决式(8)中的卷积稀疏分解问题。本文是用改进的交替方向乘法器算法(alternating direction method of multipliers, ADMM) ^[18]来解决该优化问题，该算法能较大地提高计算效率。除此之外，相应的卷积字典学习问题同样具有挑战性。大多数卷积字典学习是在稀疏编码步骤和字典更新步骤之间交替进行工作，当训练数据很大时，通常会使计算成本急剧增大。另外，传统的基于批量处理的字典学习算法是直接对所有训练数据进行处理，这还会造成大量的内存消耗，因此本文采用在线卷积字典学习^[19]算法来解决这两个限制。

CSR作为一种有效的图像表示模型，在多曝光图像融合领域主要有两个优点：1) CSR是对整幅图像进行稀疏表示，而不是基于块的多值表示，图像细节保留问题可以得到有效的解决；2) CSR是一个移位不变的图像表示模型，当融合误配准的区域时，融合质量也可以得到显著提高。

基于上述分析，本文采用在线字典学习方法，用50张尺寸为256×256的高清自然图像训练字典滤波器。为了在算法复杂度和融合质量之间取得平衡，本文设置每个字典滤波器的空间大小为$8 \times 8$，字典滤波器的数量为32。图 5是训练得到的卷积字典${\mathit{\boldsymbol{D}}}$。

图 5.  训练的卷积字典

Figure 5.  A trained convolution dictionary

下载: 全尺寸图片 幻灯片

2.3   亮度权重图生成

通过训练得到的卷积字典${\mathit{\boldsymbol{D}}}$，对核心张量中第一子带的边缘特征图${\mathit{\boldsymbol{H}}_k}$进行卷积稀疏分解，它的稀疏系数图${\mathit{\boldsymbol{X}}_{k, m}}, m \in \{ 1, ..., M\}$通过求解式(8)中的CSR模型得到。

其中：${\mathit{\boldsymbol{X}}_{k, 1:M}}(x, y)$表示第$k$幅源图像在像素$(x, y)$位置处的稀疏系数，$M$是稀疏系数图的数量，${\mathit{\boldsymbol{D}}_m}$是通过在线学习方法得到的卷积字典，$\lambda$为0.01。由于源图像中各像素点的活跃水平是衡量权重的依据，它反映了源图像的特征信息，因此活跃水平${\mathit{\boldsymbol{A}}_k}(x, y)$可用稀疏系数${\mathit{\boldsymbol{X}}_{k, 1:M}}(x, y)$的L1范数来获取。

同时，为了减少噪声的影响，对${\mathit{\boldsymbol{A}}_k}(x, y)$执行基于窗口的平均策略来获得最终的活跃水平图${\mathit{\boldsymbol{A}}_k}(x, y)$：

其中r是窗口大小。若r很大时，这种方法会对噪声更加鲁棒，但是可能会丢失一些细节。在多曝光图像融合中，源图像的边缘和纹理这些特征信息会对最终的融合结果有很大影响，故在本文中，窗口大小r设定为3。最后采用"赢者取全"策略得到初始权重图$\mathit{\boldsymbol{P}}_k^i$，图 6是House序列的初始权重图。

图 6.  House序列的初始权重图

Figure 6.  Initial weight map of the House sequence

下载: 全尺寸图片 幻灯片

然而，如图 7所示，初始权重图中存在大量噪声和边缘匹配不准的情况，这可能会在融合图像中产生边缘伪影。例如图 7(a)中的红色框部分的场景十分暗，该处的像素点处于欠曝光区域，理论上不应分配权重，但图 7(b)的红色框部分却表明分配了权重，说明该初始权重图需要进一步精细化。而使用空间一致性是解决这个问题的有效途径，这里的空间一致性是指如果两个相邻像素点具有相似的亮度，它们将具有相似的权重。而引导滤波器作为一种边缘保留滤波器，可以促使滤波后的图像保留引导图像的边缘信息，因此Kang等^[10]将初始权重图和多曝光序列分别作为滤波输入图像和引导图像，使得滤波输出图像与引导图像的边缘进行对准，从而解决了空间一致性的问题。本文在此基础上，利用引导滤波器对初始权重图进行精细化，即将初始权重图${\mathit{\boldsymbol{P}}_k}$作为滤波输入图像，相应的多曝光图像的亮度分量${\mathit{\boldsymbol{Y}}_k}$作为引导图像。同时，为了更好地保留${\mathit{\boldsymbol{Y}}_k}$中的细节信息，采用基本层和细节层分别融合的方式，而基本层图像表示亮度分量的低频信息，其看起来在空间上平滑，那么对应的权重图也应平滑，否则会出现边缘伪影现象。因此，对于基本层来说，优先选择较大的滤波半径和正则化参数，而对于细节层而言，优先选择较小的滤波半径和正则化参数。

图 7.  (a) 放大后的House序列的低曝光图像；(b)图 7(a)所对应的初始权重图

Figure 7.  (a) The enlarged low exposure image of House sequence; (b) The initial weight map of Fig. 7(a)

下载: 全尺寸图片 幻灯片

其中：${{\text{G}}_{r, \varepsilon }}( \cdot )$代表引导滤波器，${r_1}$，${r_2}$和${\varepsilon _1}$，${\varepsilon _2}$分别是引导滤波器的窗口半径和正则化参数，文中设置为45，7，0.3，10^-6。$\mathit{\boldsymbol{W}}_k^{\text{B}}$和$\mathit{\boldsymbol{W}}_k^{\text{D}}$分别是滤波后的基本层和细节层的权重图。最后，将$K$个精细化的权重图归一化，使得它们在每个像素点上的权重和为1，图 8是精细化后的权重图。

图 8.  精细化后的权重图。(a)基本层权重图；(b)细节层权重图

Figure 8.  Refined weight map. (a) The weight map of basic layer; (b) The weight map of detailed layer

下载: 全尺寸图片 幻灯片

2.4   亮度融合

每幅源图像的亮度分量${\mathit{\boldsymbol{Y}}_k}$首先被分解成基本层$\mathit{\boldsymbol{Y}}_k^{\text{b}}$和细节层$\mathit{\boldsymbol{Y}}_k^{\text{d}}$，而获取基本层和细节层的方法有很多，例如双边滤波、引导滤波等。本文通过解决以下优化问题^[20]来获得基本层：

该问题等同于利用$\eta$作为正则化参数，离散梯度作为算子的Tikhonov正则化问题，可以通过快速傅里叶变换解决。其中${{\mathit{\boldsymbol{g}}}_x} = [ - 1\;\;1]$和${{\mathit{\boldsymbol{g}}}_y} = {[ - 1\;\;1]^{\text{T}}}$分别是水平和垂直梯度算子，$\eta$是正则化参数，在本文中设置为5。通过式(15)得到基本层$\mathit{\boldsymbol{Y}}_k^{\text{b}}$后，细节层$\mathit{\boldsymbol{Y}}_k^{\text{d}}$可以通过下式得到：

最后根据式(14)和式(16)重构得到融合亮度${\mathit{\boldsymbol{Y}}_{\text{f}}}$。

2.5   色度融合

与亮度融合过程不同的是，色度融合过程直接决定了最终融合图像的颜色特征，所以本文采取不同的融合策略来处理色度分量。首先将每幅RGB源图像${\mathit{\boldsymbol{I}}_k}$转换颜色空间，得到两个色度分量${\mathit{\boldsymbol{Cb}}_k}$和${\mathit{\boldsymbol{Cr}}_k}$，过程如下：

而式(18)的逆过程可以如下表示：

根据式(19)，当$\mathit{\boldsymbol{Cb}}$、$\mathit{\boldsymbol{Cr}}$值越接近128时，图像的视觉效果越接近灰度图像，因此携带了较少的色度信息，而当$Cb$、$\mathit{\boldsymbol{Cr}}$值越远离128时，图像的颜色会越丰富。基于上述特性，Paul等^[21]利用$Cb$、$\mathit{\boldsymbol{Cr}}$值与128之间的L1范数来构造线性色度融合权重，即当色度分量越接近128时，则分配较少的权重，从而保证融合图像具有丰富的色彩信息。本文在此基础上，对该权重进行了简单的改进，即采用非线性高斯曲线进行加权。首先令${\mathit{\boldsymbol{Cb}}_k}$和${\mathit{\boldsymbol{Cr}}_k}(k = 1, 2, ..., K)$表示归一化后的源图像${\mathit{\boldsymbol{I}}_k}$的色度分量，则融合后的色度分量${\mathit{\boldsymbol{Cb}}_{\text{f}}}$和${\mathit{\boldsymbol{Cr}}_{\text{f}}}$可根据下式得到。

其中：

$\mu$是衡量色度信息的标准值，在本文中设置为0.5，即当源图像中某个像素的色度值等于0.5时，对应的权重值为0。$\sigma$控制高斯曲线的平坦程度，在本文中设置为0.5。图 9为线性加权与高斯加权方式得到的色度分量间的直方图对比结果，可以看出，用高斯加权方式得到的色度分量值在0.5附近的像素个数明显少于用线性加权方式得到的结果，而图 10(a)和10(b)分别为Paul等^[21]和本文所提出色度融合方法得到的最终融合图像。可以看出，图 10(b)在窗外的树木、室内的窗帘和墙壁上具有更加丰富的颜色，而图 10(a)在红色框内的墙壁区域有明显的颜色退化现象，因此本文所提出的高斯色度加权方法能更好地处理多曝光源图像序列的色度信息。

图 9.  高斯与线性方式加权得到的色度分量间的直方图对比

Figure 9.  Histogram comparison for chrominance components between Gaussian weighting and linear weighting algorithm

下载: 全尺寸图片 幻灯片

图 10.  两种色度融合方法得到的融合图像对比。

(a) Paul等^[21]；(b)本文提出的方法

Figure 10.  Comparison of fused images obtained by two chrominance fusion methods.

(a) Paul, etc.^[21]; (b) The proposed algorithm

下载: 全尺寸图片 幻灯片

2.6   融合图像生成

最后，根据上述求取的融合亮度和融合色度分量，通过颜色空间转换得到最终的融合图像${I_{\text{f}}}$。