数字图像采集设备受制于传感器有限的动态范围，无法捕捉真实环境的光线变化，故高动态范围(HDR)图像技术应运而生。它能捕捉整个场景动态范围的光照信息，获得更多的细节与视觉体验，具有重要的理论研究价值与商业价值。

目前主要采用硬件与软件两种方式采集HDR图像。其中，硬件方式是采用HDR传感器直接采集HDR图像，但代价过高难以实现。故当前普遍采用融合多幅不同曝光时间的低动态范围(LDR)图像序列，即多曝光图像融合的方法，以获取一幅具有完整场景的动态范围图像。传统的HDR融合一般是通过标定相机响应函数(CRF)来重建场景HDR光辐射图。

相机响应函数是一条真实场景的照度信息经过相机成像后映射到数字图像像素值的关系非线性函数，通过求解CRF的逆函数，再由多幅已知精确曝光时间的LDR图像像素值加权平均而得到HDR场景辐照图，即可还原真实场景的照度信息。Debevec等人^[1]根据多曝光图像序列采样得到数据而构建相机响应函数的标定方程，求解方程得到相机响应函数，再从多幅已知曝光时间的多曝光图像序列逆映射得到HDR场景辐射图。Mitsunaga等人^[2]以多项式的形式表示相机响应函数，只需获得相邻曝光图像的曝光时间之比即可标定相机响应函数。方华猛等人^[3]提出了快速标定相机响应函数的方法，通过降低标定方程的维数并采用稳健QR分解算法求解，引入高斯加权函数，减少融合明暗区域时引入的噪声。上述基于标定相机响应函数的多曝光图像融合方法虽能获取HDR图像，但是抗干扰能力较弱，并且求解相机响应函数的方法复杂度过高，实时性较低。

因此，多曝光融合即从图像融合的角度出发的融合方法成为近年的研究热点。多曝光图像的像素级融合算法省略了标定相机响应函数的过程，将一组同一场景但曝光时间不同的图像序列以一定的权值限定，再融合而成的一幅能够直接用于显示设备的高质量LDR图像，而不需要经过HDR图像生成并色调映射为LDR图像的过程。

Goshtasby^[4]对多曝光图像分块并基于信息熵融合图像最优的曝光图像块，利用高斯函数将最优块进行融合，以消除块与块之间存在的边界不连续现象。Mertens等人^[5]对输入图像序列进行拉普拉斯金字塔分解并根据对比度、饱和度与曝光良好率求得融合图像的权值，将图像拉普拉斯金字塔与权重高斯金字塔进行融合，由融合后的金字塔重构HDR图像。Shen等人^[6]提出一种基于曝光质量测量以及最小可察觉失真获取多曝光图像融合权重的方法，同时利用快速的拉普拉斯金字塔加速分解图像纹理层与细节层，分别对各层给予对应权重进行融合，融合后的图像具有保持细节与色彩的效果。Bruce^[7]提出了一种基于局部熵获取像素融合权重的方法，并把图像转换到对数域进行处理，通过设置权值的大小突出局部细节的对比度。由于直接融合多曝光图像像素的方法除了能压缩动态范围，还能作为HDR图像重现方法。付争方等人^[8]将多曝光图像直接融合结合Sigmoid函数拟合，直接从获得的LDR图像中提取每个像素位置的亮度序列，利用最小二乘法来拟合适应视觉的S曲线，从而建立起亮度信息的数学模型，给出最佳成像亮度的判决方式，快速有效地合成HDR图像。

由于多曝光图像序列不可避免会存在物体移动与噪声污染的问题，特别是弱光源场景，而上述方法对这些问题的处理能力不足。

低秩矩阵恢复方法可从稀疏噪声污染的数据矩阵中提取出低秩的矩阵部分，故Peng等人^[9]从批量线性相关的图像中分解出对齐后线性相关的低秩图像与稀疏噪声。这为多曝光图像序列的对齐与融合提供了很好的思路，可在解决多曝光图像序列对齐(低秩)问题的同时，减少因噪声污染与运动物体而产生的伪影(稀疏的噪声)。

因此，本文从低秩矩阵恢复的角度出发，展开多曝光HDR图像融合方法的研究，并对低秩矩阵恢复算法进行改进，提出一种新的HDR图像融合方法。

随着维度的增加，数据的冗余度与相关性也会随之增加，从而引发维度灾难，增加了处理的困难与成本。在实际应用中，常假设高维数据近似存在于一个低维线性子空间中，即可利用低维线性子空间的性质对数据进行特征提取、噪声移除等处理。传统的线性子空间方法对含有高斯噪声数据具有较好的鲁棒性，但对稀疏噪声数据却非常敏感^[10]。

根据压缩感知理论^[11-12]，稀疏的数据对稀疏噪声具有更好的鲁棒性^[13]。而低秩矩阵恢复(low-rank matrix recovery)又称为鲁棒主成分分析(RPCA)^[14]，指的是将数据矩阵分解为低秩矩阵与稀疏矩阵之和。

它以矩阵的秩作为目标函数，其奇异值构成的向量具有稀疏性。于是，数据的稀疏表示便可拓展到矩阵的低秩上。从多曝光的LDR图像序列所形成的多维度线性空间便可还原场景真实光照的低维度空间，同时抑制场景中运动物体所产生的伪影，即稀疏噪声。因此，本文将从RPCA角度解决HDR图像融合问题。

令不同曝光参数拍摄$N$幅($m$×$n$)静止的LDR图像构成的序列为{${\mathit{\boldsymbol{I}}_1}$, …, ${\mathit{\boldsymbol{I}}_N}$}，向量化每一幅图像作为数据矩阵的列向量，则$N$幅$m$×$n$的数据矩阵$\mathit{\boldsymbol{L}}$转换为一个低秩矩阵，即

$ \mathit{\boldsymbol{L}} = \left[ {vec\left( {{\mathit{\boldsymbol{I}}_1}} \right), \cdots ,vec\left( {{\mathit{\boldsymbol{I}}_N}} \right)} \right] $

(1)

式中，$vec$(.)表示列向量函数，${\mathit{\boldsymbol{I}}_i}$为第$i$幅图像。

在拍摄过程中，难免存在相机抖动，或者场景突然出现运动物体等情况。图像序列未对齐，便对其进行融合处理，必然会导致伪影与模糊现象。令伪影与模糊部分为噪声${\mathit{e}_i}$，根据相机成像原理，则第$i$幅图像${\mathit{\boldsymbol{I}}_i}$为

$ {\mathit{\boldsymbol{I}}_i} = f\left( {\left( {{R_i} + {e_i}} \right)\Delta {t_i}} \right) $

(2)

式中，${R_i}$表示场景的辐射照度，$f$(·)为相机响应函数，Δ${t_i}$为曝光时间，则式(2)的逆函数为

$ {f^{ - 1}}\left( {{\mathit{\boldsymbol{I}}_i}} \right) = {R_i}\Delta {t_i} + {e_i}\Delta {t_i} = {\mathit{\boldsymbol{L}}_i} + {\mathit{\boldsymbol{E}}_i} $

(3)

式中，${\mathit{\boldsymbol{L}}_i}$为第$i$幅图像的场景辐射照度矩阵，${\mathit{\boldsymbol{E}}_i}$为第$i$幅图像的噪声矩阵。经过相机响应函数逆映射后，亮度满足线性关系，场景的辐射照度$\mathit{\boldsymbol{R}}$将具有低秩特点。

令数据矩阵为$\mathit{\boldsymbol{D = }}\left[{vec\left( {{{I'}_1}} \right), \cdots, \mathit{vec}\left( {{{\mathit{\boldsymbol{I'}}}_N}} \right)} \right] \in {{\bf{R}}^{M \times N}}$，其中，${{\mathit{\boldsymbol{I'}}}_i} = {f^{-1}}\left( {{\mathit{\boldsymbol{I}}_i}} \right)$，$M$与$N$分别表示一幅图像像素个数与多曝光图像序列的数量。令低秩矩阵为$\mathit{\boldsymbol{L = }}\left[{vec\left( {{\mathit{\boldsymbol{L}}_1}} \right), \cdots, \mathit{vec}\left( {{\mathit{\boldsymbol{L}}_N}} \right)} \right]$，稀疏噪声矩阵为$\mathit{\boldsymbol{E = }}\left[{vec\left( {{\mathit{\boldsymbol{E}}_1}} \right), \cdots, \mathit{vec}\left( {{\mathit{\boldsymbol{E}}_N}} \right)} \right]$，则式(3)变换为

$ \mathit{\boldsymbol{D}} = \mathit{\boldsymbol{L}} + \mathit{\boldsymbol{E}} $

(4)

若场景中不存在伪影，即矩阵$\mathit{\boldsymbol{E}}$为零矩阵，对齐后的数据矩阵$\mathit{\boldsymbol{D}}$将与场景辐射照度矩阵$\mathit{\boldsymbol{L}}$相等，则秩为1。但噪声使得数据矩阵$\mathit{\boldsymbol{D}}$的秩不可能为1。假设噪声矩阵$\mathit{\boldsymbol{E}}$是稀疏的，则可分解数据矩阵$\mathit{\boldsymbol{D}}$为秩接近1的低秩矩阵$\mathit{\boldsymbol{L}}$与稀疏的噪声矩阵$\mathit{\boldsymbol{E}}$，即

$ \begin{array}{*{20}{c}} {\mathop {\min }\limits_{\mathit{\boldsymbol{L}},\mathit{\boldsymbol{E}}} rank\left( \mathit{\boldsymbol{L}} \right) = \lambda {{\left\| \mathit{\boldsymbol{E}} \right\|}_0}}\\ {{\rm{s}}.\;{\rm{t}}.\;\;\mathit{\boldsymbol{D}} = \mathit{\boldsymbol{L}} + \mathit{\boldsymbol{E}}} \end{array} $

(5)

式中，$\mathit{\boldsymbol{D}}$, $\mathit{\boldsymbol{L}}$, $\mathit{\boldsymbol{E}}$∈${{\bf{R}}^{{n_1} \times {n_2}}}$，$\lambda $>0是权重参数，$rank$($\mathit{\boldsymbol{X}}$)是集合$\left\{ {\mathit{\boldsymbol{X}} \in {{\bf{R}}^{{n_1} \times {n_2}}}\left| {\left\| \mathit{\boldsymbol{X}} \right\| \le 1} \right.} \right\}$上的凸包(即$\mathit{\boldsymbol{X}}$的核范数)。但因求秩函数与l₀范数都是非凸的，求解该优化问题是NP难题。但在满足一定的条件下，松弛目标函数可以较大概率得到近似低秩的最优解，以及相对稀疏的噪声矩阵。利用凸松弛将求秩函数转换为核函数，l₀范数转换为l₁范数，可求得最优解，则目标函数可重写为

$ \begin{array}{*{20}{c}} {\mathop {\min }\limits_{\mathit{\boldsymbol{L}},\mathit{\boldsymbol{E}}} {{\left\| \mathit{\boldsymbol{L}} \right\|}_ * } + \lambda {{\left\| \mathit{\boldsymbol{E}} \right\|}_1}}\\ {{\rm{s}}.\;{\rm{t}}.\;\;\mathit{\boldsymbol{D}} = \mathit{\boldsymbol{L}} + \mathit{\boldsymbol{E}}} \end{array} $

(6)

该凸优化问题又称为主成分追踪(principal component pursuit)^[15]。只要低秩矩阵$\mathit{\boldsymbol{L}}$的奇异值分布合理并且稀疏矩阵$\mathit{\boldsymbol{E}}$的非零元素均匀分布，RPCA即可以接近1的概率恢复出低秩矩阵。

实际上，多曝光图像序列通常较少，通常只有3幅或5幅图像，无法令RPCA求解得到低秩的最优解，且核函数作为目标函数虽使得恢复的矩阵具有低秩的性质，却对数据的重要程度不敏感。根据Oh等人^[16-17]提出的以部分奇异值和(PSSV)作为优化目标函数，尽可能减少对最大奇异值的影响，不断减少其余奇异值所代表的数据比重，从而可保证即使在只有少量数据的条件下，也能够得到更低秩的矩阵。于是，将式(5)转换为

$ \begin{array}{*{20}{c}} {\mathop {\min }\limits_{\mathit{\boldsymbol{L}},\mathit{\boldsymbol{E}}} \left| {rank\left( \mathit{\boldsymbol{L}} \right) - N} \right| + \lambda {{\left\| \mathit{\boldsymbol{E}} \right\|}_0}}\\ {{\rm{s}}.\;{\rm{t}}.\;\;\mathit{\boldsymbol{D}} = \mathit{\boldsymbol{L}} + \mathit{\boldsymbol{E}}} \end{array} $

(7)

该优化问题求解得到的低秩矩阵$\mathit{\boldsymbol{L}}$秩接近目标大小$N$。定义投影算子${P_r}\left( \mathit{\boldsymbol{X}} \right) = \mathit{\boldsymbol{U}}_{1:\mathit{r}}^{\rm{T}}\mathit{\boldsymbol{X}}{\mathit{\boldsymbol{V}}_{1:\mathit{r}}}$，其中${\mathit{\boldsymbol{U}}_{1:\mathit{r}}}$与${\mathit{\boldsymbol{V}}_{1:\mathit{r}}}$分别为矩阵$\mathit{\boldsymbol{X}}$的左奇异值矩阵与右奇异值矩阵中前$r$个最大奇异值所对应的向量构成的矩阵。当$rank$($\mathit{\boldsymbol{X}}$)≥$r$时，${P_r}$($\mathit{\boldsymbol{X}}$)的秩等于$r$。利用投影算子，可将式(7)的第1项进行松弛，得

$ \begin{array}{*{20}{c}} {\left| {{{\left\| \mathit{\boldsymbol{L}} \right\|}_ * } - {{\left\| {{P_N}\left( \mathit{\boldsymbol{L}} \right)} \right\|}_ * }} \right| = }\\ {\left| {\sum\limits_{i = 1}^{rank\left( L \right)} {{\sigma _i}\left( \mathit{\boldsymbol{L}} \right)} - \sum\limits_{i = 1}^N {{\sigma _i}\left( \mathit{\boldsymbol{L}} \right)} } \right| = }\\ {\sum\limits_{i = N + 1}^{rank\left( L \right)} {{\sigma _i}\left( \mathit{\boldsymbol{L}} \right)} = {{\left\| \mathit{\boldsymbol{L}} \right\|}_{p = N}}} \end{array} $

(8)

式中，$\sigma $(·)为矩阵的奇异值函数；${\left\| \cdot \right\|_{p = N}}$就是令目标矩阵秩为$N$的PSSV范数，即截断的核范数(truncated nuclear norm)^[18]。则基于PSSV的多曝光HDR图像融合的低秩矩阵模型可表示为

$ \begin{array}{*{20}{c}} {\mathop {\min }\limits_{\mathit{\boldsymbol{L}},\mathit{\boldsymbol{E}}} {{\left\| \mathit{\boldsymbol{L}} \right\|}_{p = 1}} + \lambda {{\left\| \mathit{\boldsymbol{E}} \right\|}_1}}\\ {{\rm{s}}.\;{\rm{t}}.\;\;\mathit{\boldsymbol{D}} = \mathit{\boldsymbol{L}} + \mathit{\boldsymbol{E}}} \end{array} $

(9)

对式(9)进行低秩矩阵恢复求解，即可得到一幅融合后的HDR图像。

目前，求解低秩矩阵恢复的主要方法有迭代阈值法(iterative thresholding)^[19]，加速近端梯度算法(accelerated proximal gradient)^[20]，对偶方法(dual method)^[20]以及增广拉格朗日乘子法(ALM)^[21]等。

由于增广拉格朗日乘子法效率高、精度较高与占用存储空间较少等优点，故本文采用精确增广拉格朗日乘子法(EALM)求解以PSSV为目标函数的低秩矩阵。

对于一般的带约束优化问题

$ \min f\left( \mathit{\boldsymbol{X}} \right){\rm{s}}.\;{\rm{t}}.\;\;h\left( \mathit{\boldsymbol{X}} \right) = 0 $

(10)

式中，$f:{{\bf{R}}^n} \to {\bf{R}}, \mathit{h} :{{\bf{R}}^n} \to {{\bf{R}}^m}$。构造增广拉格朗日函数

$ \begin{array}{*{20}{c}} {L\left( {\mathit{\boldsymbol{X}},\mathit{\boldsymbol{Y}},\mu } \right) = f\left( \mathit{\boldsymbol{X}} \right) + < \mathit{\boldsymbol{Y}},h\left( \mathit{\boldsymbol{X}} \right) > + }\\ {\frac{\mu }{2}\left\| {h\left( \mathit{\boldsymbol{X}} \right)} \right\|_{\rm{F}}^2} \end{array} $

(11)

式中，$\mathit{\boldsymbol{Y}}$是拉格朗日乘子向量，$\mu $>0是惩罚因子。<$A$, $B$>=tr(${\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{B}}$)为矩阵内积。${\left\| \mathit{\boldsymbol{X}} \right\|_{\rm{F}}} = \sqrt {\sum\limits_i {\sum\limits_j {x_{ij}^2} } } $为Frobenius范数。带约束的优化问题式(10)便可转换为无约束的优化问题min $L$($\mathit{\boldsymbol{X}}$, $\mathit{\boldsymbol{Y}}$, $\mu $)。令$k$为迭代步数，即

$ \begin{array}{l} {\mathit{\boldsymbol{X}}_{k + 1}} = \arg \mathop {\min }\limits_\mathit{\boldsymbol{X}} L\left( {\mathit{\boldsymbol{X}},{\mathit{\boldsymbol{Y}}_k},{\mu _k}} \right)\\ {\mathit{\boldsymbol{Y}}_{k + 1}} = {\mathit{\boldsymbol{Y}}_k} + {\mu _k}h\left( {{\mathit{\boldsymbol{X}}_{k + 1}}} \right)\\ {\mu _{k + 1}} = \rho {\mu _k} \end{array} $

(12)

式中，$ρ$>0为加快收敛的权重参数。迭代执行使得$\mathit{\boldsymbol{X}}$最终收敛到最优解。采用ALM求解以PSSV为目标函数的低秩矩阵恢复问题，需做改变

$ \begin{array}{l} \mathit{\boldsymbol{X}} = \left( {\mathit{\boldsymbol{L}},\mathit{\boldsymbol{E}}} \right)\\ f\left( \mathit{\boldsymbol{X}} \right) = {\left\| \mathit{\boldsymbol{L}} \right\|_{p = N}} + \lambda {\left\| \mathit{\boldsymbol{E}} \right\|_1}\\ h\left( \mathit{\boldsymbol{X}} \right) = \mathit{\boldsymbol{D}} - \mathit{\boldsymbol{L}} - \mathit{\boldsymbol{E}} \end{array} $

(13)

则求解该问题的拉格朗日方程

$ \begin{array}{*{20}{c}} {L\left( {\mathit{\boldsymbol{L}},\mathit{\boldsymbol{E}},\mathit{\boldsymbol{Y}},\mu } \right) = {{\left\| \mathit{\boldsymbol{L}} \right\|}_{p = N}} + \lambda {{\left\| \mathit{\boldsymbol{E}} \right\|}_1} + }\\ { < \mathit{\boldsymbol{Y}},\mathit{\boldsymbol{D}} - \mathit{\boldsymbol{L}} - \mathit{\boldsymbol{E}} > + }\\ {\frac{\mu }{2}\left\| {\mathit{\boldsymbol{D}} - \mathit{\boldsymbol{L}} - \mathit{\boldsymbol{E}}} \right\|_{\rm{F}}^2} \end{array} $

(14)

根据交替方向乘子法(ADMM)的思想，优化问题又可分解为2个求解矩阵$\mathit{\boldsymbol{L}}$与$\mathit{\boldsymbol{E}}$的子优化问题

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{L}}_{k + 1}} = \arg \mathop {\min }\limits_\mathit{\boldsymbol{L}} L\left( {\mathit{\boldsymbol{L}},{\mathit{\boldsymbol{E}}_k},{\mathit{\boldsymbol{Y}}_k},{\mu _k}} \right) = }\\ {\arg \mathop {\min }\limits_\mathit{\boldsymbol{L}} \mu _k^{ - 1}{{\left\| \mathit{\boldsymbol{L}} \right\|}_{p = N}} + }\\ {\frac{1}{2}\left\| {\mathit{\boldsymbol{L}} - \left( {\mathit{\boldsymbol{D}} - {\mathit{\boldsymbol{E}}_k} + \mu _k^{ - 1}{\mathit{\boldsymbol{Y}}_k}} \right)} \right\|_{\rm{F}}^2}\\ {{\mathit{\boldsymbol{E}}_{k + 1}} = \arg \mathop {\min }\limits_\mathit{\boldsymbol{E}} L\left( {{\mathit{\boldsymbol{L}}_{k + 1}},\mathit{\boldsymbol{E}},{\mathit{\boldsymbol{Y}}_k},{\mu _k}} \right) = }\\ {\arg \mathop {\min }\limits_\mathit{\boldsymbol{E}} \lambda \mu _k^{ - 1}{{\left\| \mathit{\boldsymbol{E}} \right\|}_1} + }\\ {\frac{1}{2}\left\| {\mathit{\boldsymbol{E}} - \left( {\mathit{\boldsymbol{D}} - {\mathit{\boldsymbol{L}}_{k + 1}} + \mu _k^{ - 1}{\mathit{\boldsymbol{Y}}_k}} \right)} \right\|_{\rm{F}}^2} \end{array} $

(15)

定义部分奇异值阈值算子(PSVT)^[17]为

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{F}}_{N,\tau }}\left[ \mathit{\boldsymbol{Y}} \right] = \mathit{\boldsymbol{U}}\left( {{\mathit{\boldsymbol{S}}_{{Y_1}}} + {\mathit{\boldsymbol{H}}_\tau }\left[ {{\mathit{\boldsymbol{S}}_{{Y_2}}}} \right]} \right){\mathit{\boldsymbol{V}}^{\rm{T}}} = }\\ {\arg \mathop {\min }\limits_\mathit{\boldsymbol{X}} \frac{1}{2}\left\| {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{Y}}} \right\|_{\rm{F}}^2 + \tau {{\left\| \mathit{\boldsymbol{X}} \right\|}_{p = N}}} \end{array} $

(16)

$ {\mathit{\boldsymbol{S}}_{{Y_1}}} = {\rm{diag}}\left( {{\sigma _1}, \cdots ,{\sigma _N},0, \cdots ,0} \right) $

$ {\mathit{\boldsymbol{S}}_{{Y_2}}} = {\rm{diag}}\left( {0, \cdots ,0,{\sigma _{N + 1}}, \cdots ,{\sigma _l}} \right) $

式中，$\tau $>0，$l$=$rank$($\mathit{\boldsymbol{Y}}$)，${\sigma _i}$为奇异值，$\mathit{\boldsymbol{U}}$、$\mathit{\boldsymbol{V}}$与$\mathit{\boldsymbol{S}} = {\mathit{\boldsymbol{S}}_{{Y_1}}}+{\mathit{\boldsymbol{S}}_{{Y_2}}}$分别是矩阵$\mathit{\boldsymbol{Y}}$的奇异值分解结果，即

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{Y}} = \mathit{\boldsymbol{US}}{\mathit{\boldsymbol{V}}^{\rm{T}}} = \mathit{\boldsymbol{U}}\left( {{\mathit{\boldsymbol{S}}_{{Y_1}}} + {\mathit{\boldsymbol{S}}_{{Y_2}}}} \right){\mathit{\boldsymbol{V}}^{\rm{T}}},}\\ {{\mathit{\boldsymbol{H}}_\tau }\left[ x \right] = {\mathop{\rm sgn}} \left( x \right)\max \left\{ {\left| x \right| - \tau ,0} \right\}} \end{array} $

为软阈值算子。采用PSVT算子更新矩阵${\mathit{\boldsymbol{L}}_{k + 1}}$以及矩阵${\mathit{\boldsymbol{E}}_{k + 1}}$, 即

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{L}}_{k + 1}} = {\mathit{\boldsymbol{P}}_{N,\mu _k^{ - 1}}}\left[ {\mathit{\boldsymbol{D}} - {\mathit{\boldsymbol{E}}_k} + \mu _k^{ - 1}{\mathit{\boldsymbol{Y}}_k}} \right]}\\ {{\mathit{\boldsymbol{E}}_{k + 1}} = {\mathit{\boldsymbol{H}}_{\lambda \mu _k^{ - 1}}}\left[ {\mathit{\boldsymbol{D}} - {\mathit{\boldsymbol{L}}_{k + 1}} + \mu _k^{ - 1}{\mathit{\boldsymbol{Y}}_k}} \right]} \end{array} $

(17)

更新矩阵${\mathit{\boldsymbol{L}}_{k + 1}}$与矩阵${\mathit{\boldsymbol{E}}_{k + 1}}$后，更新拉格朗日乘子矩阵${\mathit{\boldsymbol{E}}_{k + 1}}$以及惩罚因子${\mu _k}$

$ \begin{array}{l} {\mathit{\boldsymbol{Y}}_{k + 1}} = {\mathit{\boldsymbol{Y}}_k} + {\mu _k}\left( {\mathit{\boldsymbol{D}} - {\mathit{\boldsymbol{L}}_{k + 1}} - {\mathit{\boldsymbol{E}}_{k + 1}}} \right)\\ {\mu _{k + 1}} = \rho {\mu _k} \end{array} $

(18)

EALM方法与Bhardwaj等人^[22]所提出的融合方法不同之处在于，文献[22]中的目标函数为凸松弛后的RPCA问题，且求解方法为非精确拉格朗日乘子法(IALM)。文献[21]证明了在EALM算法中(${\mathit{\boldsymbol{L}}_k}$, ${\mathit{\boldsymbol{E}}_k}$)任意一个聚点($\mathit{\boldsymbol{L}}$, $\mathit{\boldsymbol{E}}$)都是RPCA问题的最优解，并且收敛速度至少是${\rm{O}}\left( {\mu _k^{-1}} \right)$。如果不采用反复迭代的方法求解矩阵$\mathit{\boldsymbol{L}}$与$\mathit{\boldsymbol{E}}$，即为IALM算法，在${\mu _k}$增长不过快时，可使得$\sum\limits_{k = 1}^{ + \infty } {\mu _k^{-2}{\mu _{k + 1}} < \infty } $且$\mathop {{\rm{lim}}}\limits_{k \to + \infty } {\mu _k}\left( {{\mathit{\boldsymbol{E}}_{k + 1}}-{\mathit{\boldsymbol{E}}_k}} \right) = 0$，它也可保证RPCA的收敛到最优解$\left( {{\mathit{\boldsymbol{L}}_k}, {\mathit{\boldsymbol{E}}_k}} \right)$。

此外，求解低秩矩阵恢复问题时，要求低秩矩阵最大的奇异值能够涵盖绝大部分数据，而代表噪声的奇异值则要尽可能小。奇异值代表噪声的方差，数据矩阵中的最大奇异值表示数据在对应的奇异值向量上具有最大方差，即数据在该向量上具有最大投影。

在数据较为充分的情况下RPCA得到的最优解足够低秩，但数据较少时结果并不理想。而求解PSSV问题使得最优解主要数据部分，即秩等于$N$的数据不受到影响，提高了数据的稳定性，但它只保证数据在其余奇异值向量上方差不大，并不能保证低秩部分的数据在奇异值向量上有最大方差。

对不同的奇异值设置自适应的惩罚因子$\tau $，将使得最优解尽量集中在最大奇异值的空间。可定义Sigmoid函数为

$ {S_{a,b}} = {\left( {1 + {{\rm{e}}^{ - a\left( {x - b} \right)}}} \right)^{ - 1}} $

(19)

该函数设置不同的参数$a$与$b$时，函数形状如图 1所示。

根据Sigmoid函数定义一种新的自适应软阈值算子

$ {{S'}_\tau } = \left\{ {\begin{array}{*{20}{c}} \begin{array}{l} \sigma - \tau {S_{a,b}}\left( \sigma \right)\\ 0 \end{array}&\begin{array}{l} \sigma > \tau + \tau {S_{a,b}}\left( \sigma \right)\\ \sigma \le \tau + \tau {S_{a,b}}\left( \sigma \right) \end{array} \end{array}} \right. $

(20)

根据输入的奇异值$\sigma $来控制优化方程中惩罚因子$\tau $的大小，那么新的PSVT算子${{F'}_{N, \tau }}\left[\cdot \right]$为

$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{F'}}}_{N,\mu _k^{ - 1}}}\left[ \mathit{\boldsymbol{Y}} \right] = \mathit{\boldsymbol{U}}\left( {{\mathit{\boldsymbol{S}}_{{Y_1}}} + {{\mathit{\boldsymbol{H'}}}_\tau }\left[ {{\mathit{\boldsymbol{S}}_{{Y_2}}}} \right]} \right){\mathit{\boldsymbol{V}}^{\rm{T}}} = }\\ {\arg \mathop {\min }\limits_\mathit{\boldsymbol{X}} \frac{1}{2}\left\| {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{Y}}} \right\|_{\rm{F}}^2 + \tau {{\left\| \mathit{\boldsymbol{X}} \right\|}_{p = N}}} \end{array} $

(21)

式中，${\mathit{\boldsymbol{S}}_{{Y_1}}} = {\rm{diag}}\left( {{\sigma _1}, \cdots, {\sigma _N}, 0, \cdots, 0} \right)$

${\mathit{\boldsymbol{S}}_{{Y_2}}} = {\rm{diag}}\left( {{\sigma _1}, \cdots, {\sigma _{N + 1}}, 0, \cdots, {\sigma _l}} \right)$

则更新后的低秩矩阵$\mathit{\boldsymbol{L}}$可表示为

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{L}}_{k + 1}} = \arg \mathop {\min }\limits_\mathit{\boldsymbol{L}} L\left( {\mathit{\boldsymbol{L}},{\mathit{\boldsymbol{E}}_k},{\mathit{\boldsymbol{Y}}_k},{\mu _k}} \right)}\\ { = \arg \mathop {\min }\limits_\mathit{\boldsymbol{L}} \mu _k^{ - 1}{{\left\| \mathit{\boldsymbol{L}} \right\|}_{p = N}} + }\\ {\frac{1}{2}\left\| {\mathit{\boldsymbol{L}} - \left( {\mathit{\boldsymbol{D}} - {\mathit{\boldsymbol{E}}_k} + \mu _k^{ - 1}{\mathit{\boldsymbol{Y}}_k}} \right)} \right\|_{\rm{F}}^2 = }\\ {{{\mathit{\boldsymbol{F'}}}_{N,\mu _k^{ - 1}}}\left[ {\mathit{\boldsymbol{D}} - {\mathit{\boldsymbol{E}}_k} + \mu _k^{ - 1}{\mathit{\boldsymbol{Y}}_k}} \right]} \end{array} $

(22)

当矩阵$\mathit{\boldsymbol{X}}$奇异值${\sigma _i} \to 0$, $i \in \left\{ {1, \cdots, rank\left( \mathit{\boldsymbol{X}} \right)} \right\}$，控制参数$a$与$b$使得${{F'}_{N, \tau }}\left[\mathit{\boldsymbol{X}} \right] \to \mathit{\boldsymbol{X}}$；当${\sigma _i} \to + \infty $, ${{\mathit{\boldsymbol{F'}}}_{N, \tau }}\left[\mathit{\boldsymbol{X}} \right] \to {\mathit{\boldsymbol{F}}_{N, \tau }}\left[\mathit{\boldsymbol{X}} \right]$，仍以PSVT求解该优化问题；当处于Sigmoid函数变化的部分时，自适应惩罚因子使得较小的奇异值有较大的惩罚因子。本文取$a$=4，$b$=0。

于是，求解低秩矩阵的工作流程如图 2所示。先对输入的数据矩阵构造增广拉格朗日函数，然后计算奇异矩阵的奇异值。对于奇异值趋于0的情况，采用新的PSVT算子以更新低秩矩阵与噪声矩阵；否则采用原PSVT算子更新低秩矩阵与噪声矩阵。同时更新拉格朗日乘子向量与惩罚因子。若此时算法已收敛，则输出低秩矩阵与噪声矩阵；否则重新计算奇异值，进入下一次迭代。

已知基于ALM求解低秩矩阵恢复问题的收敛性与{${\mu _k}$}有关，如果{${\mu _k}$}有界该求解方法线性收敛，如果{${\mu _k}$}无界该求解方法超线性收敛。

从实验的角度观察收敛性，随机生成大小为10 000×50像素、秩为$r$的数据矩阵，并在数据矩阵中添加比例为$p$的稀疏噪声。第1组实验令生成矩阵的秩$r$=1，测试多组不同的噪声比例$p$，结果如图 3(a)所示。第2组实验令稀疏噪声比例为0.2，测试多组不同秩$r$，结果如图 3(b)所示。由恢复低秩矩阵每一步的迭代结果与原始矩阵的误差可知，本文提出方法有着较好的收敛性。

本节将对低秩矩阵恢复以及本文所提出的多曝光HDR图像融合方法进行仿真，并与现有的算法比较分析。对于图像质量本文将从主观与客观指标两方面来进行评价。

由于HDR图像存储的是与绝对亮度有关的信息，而人眼是在感受相对亮度来察觉不同亮度等级下的细节信息，因此对于以相对亮度作为参考的峰值信噪比(PSNR)与结构相似度指标(SSIM)图像质量标准并不适用于存储绝对亮度的HDR图像的质量评价。

TO等人^[23]通过研究对比度感知函数，提出了一种感知一致性(PU)映射，利用PU映射函数可以将10^-5~10⁸ cd/m²的亮度值转换到一致的感知尺度中，令PSNR以及SSIM评价的结果更加可靠。故本文采用基于PU映射的PSNR与SSIM对融合后的HDR图像进行客观评价。

实验生成数据的方法是构造1个左矩阵${\mathit{\boldsymbol{A}}_{\rm{L}}} \in {{\bf{R}}^{M \times r}}$，1个右矩阵${\mathit{\boldsymbol{A}}_{\rm{R}}} \in {{\bf{R}}^{r \times N}}$，${\mathit{\boldsymbol{A}}_{\rm{L}}}$与${\mathit{\boldsymbol{A}}_{\rm{R}}}$满足独立同分布，其值满足均值为0方差为1的正态分布，令秩为$r$的数据矩阵$H \in {{\bf{R}}^{M \times N}}$为$H = {\mathit{\boldsymbol{A}}_{\rm{L}}} \cdot {\mathit{\boldsymbol{A}}_{\rm{R}}}$。设置稀疏噪声比例为$p$∈[0, 1]，在数据矩阵$\mathit{\boldsymbol{A}}$中以均匀分布的概率抽选出$p$×$M$×$N$个位置并置换为稀疏噪声，噪声大小同样满足均值为0方差为1的正态分布，此时得到观测数据矩阵$\mathit{\boldsymbol{D}}$。通过求解RPCA、PSSV以及本文的低秩矩阵恢复算法得到的低秩矩阵$\mathit{\boldsymbol{L}}$，由此观察矩阵$\mathit{\boldsymbol{L}}$与原始数据矩阵$\mathit{\boldsymbol{A}}$之间的客观指标，包括归一化均方误差以及低秩矩阵$\mathit{\boldsymbol{L}}$的奇异值。

令生成矩阵的秩$r$=1，尺寸为10 000×50像素，通过改变噪声比例$p$，比较RPCA、PSSV以及本文算法的归一化均方误差，如图 4所示。本文算法比PSSV具有稍好的性能，且均优于RPCA。

令数据矩阵尺寸为10 000×50像素，设置多组噪声比例$p$与生成矩阵的秩$r$，由图 5可以看出数据矩阵随着噪声比例的增加，其低秩矩阵恢复误差逐渐增大；随着数据矩阵的秩$r$变大，数据矩阵非线性相关的数据增多，使得低秩矩阵恢复的误差逐渐增大。

令数据矩阵的秩$r$=1，设置多组噪声比例$p$与矩阵尺寸，由表 1可以看出, 本文算法恢复的低秩矩阵更加低秩，其占主要信息的奇异值大小与原始矩阵的奇异值比较，不会出现非常大的衰减，还能保证非重要数据的奇异值足够小。但由表 1最后一组数据再次证明，随着矩阵尺寸与噪声比例的增大，秩也跟着增大，非线性相关的数据增多，衰减也会随着增大。

PSSV与本文算法在较少数据的情况下恢复的低秩结果都要比RPCA的结果好，非常适用于多曝光图像融合数据集并不多的情况，而本文提出的算法对比于PSSV能够在较低重构误差的前提下得到更加低秩的结果。

通过对多曝光图像序列组成的低秩矩阵进行求解，得到了场景的低秩矩阵$\mathit{\boldsymbol{L}}$，它包含了多曝光图像序列完整的场景光照信息。根据式(3)的关系，可由低秩矩阵的平均值，得到场景的辐射照度

$ R = \frac{1}{N}\sum\limits_{i = 1}^N {\frac{{{\mathit{\boldsymbol{L}}_i}}}{{\Delta {t_i}}}} $

(23)

从而，求解低秩矩阵融合多曝光图像序列以获取HDR图像的工作流程如图 6所示。首先将获得的多曝光图像序列通过标定的相机响应函数逆映射到线性的亮度空间，接着将向量化得到的数据矩阵$\mathit{\boldsymbol{D}}$进行RGB处理，即分为红、绿、蓝3个通道依次进行低秩矩阵的恢复，采用EALM即精确拉格朗日乘子法来求解低秩矩阵恢复，最终将3个通道得到的低秩矩阵$\mathit{\boldsymbol{L}}$求平均即可获得目标HDR图像。

由多曝光HDR图像融合的低秩模型，以Memorial Church^[1]与Arch^[24]多曝光图像序列为实验数据集，进行多曝光HDR图像融合，并将仿真结果分别与文献[1]、RPCA与PSSV方法进行主观与客观的比较。此外，由于无法直接通过展示HDR图像进行效果比对，故本文一致采用${\mathit{\boldsymbol{L}}_d} = {\left( {\mathit{\boldsymbol{L}}/{\mathit{\boldsymbol{L}}_{{\rm{avg}}}}} \right)^\gamma }$压缩HDR图像的动态范围。其中，参数$\gamma $表示压缩因子，本文取值为$\gamma $=1/2.2，${{\mathit{\boldsymbol{L}}_{{\rm{avg}}}}}$是平均亮度，${{\mathit{\boldsymbol{L}}_d}}$是压缩了动态范围的LDR图像。

对图 7中的Memorial Church图像序列作HDR图像融合，首先采用传统的HDR融合方法，即文献[1]计算出多曝光图像序列的相机响应函数，并将该图像序列进行融合得到一幅没有任何噪声的HDR参考图像。如图 8为融合后的各个算法的结果图以及相对参考图像图 8(a)所计算出的PSNR与SSIM(01之间变化，值越大越相似)指标。

从主观角度看，上述方法的融合结果基本能展示出真实场景的光线变化。但文献[1]的图像整体亮度偏暗，对比度相对其余算法都要差。RPCA方法在恢复屋顶的圆形部分时出现了错误(色彩出现异常，细节完全丢失)，如图 9为局部细节放大图。

从图 9可以进一步观察到上述算法对左侧屋顶细节部分都有着较好的融合效果，扩展了图像的动态范围。但文献[1]结果偏暗。从圆形屋顶的局部放大图可明显看到RPCA的融合效果相对较差，文献[1]与PSSV以及本文算法能正确融合多曝光图像序列，文献[1]隐约可见圆形屋顶的线条，PSSV与本文算法，比文献[1]丢失了不少细节。

接着，在Memorial Church图像序列中添加均匀分布的噪声污染。添加稀疏噪声后的图像序列如图 10所示，对该序列进行HDR融合来检验各算法的性能。添加噪声后的图像序列如图 10所示，对该序列进行HDR融合来检验各算法的性能。

融合后的结果以及PSNR与SSIM指标如图 11所示，文献[1]的方法并不能对噪声进行处理，导致融合的结果出现了噪声，如图 11 (a)所示。低秩矩阵恢复不受稀疏噪声的影响，消除了稀疏噪声对融合结果的影响，基于低秩矩阵恢复的多曝光图像融合方法具有较好的鲁棒性。

从客观指标上来分析，本文算法融合结果的PSNR与SSIM指标均高于其余算法，效果较好。基于低秩矩阵恢复的多曝光图像融合方法在处理原图像序列与加噪图像序列的评价指标基本相同，本文方法在低秩对齐的同时保证了融合图像的质量。

图 12(a)为Arch图像序列，属于一般的自然场景，在拍摄过程中存在运动的物体，直接融合会导致伪影的出现。目前还未有成熟的去伪影效果客观评价指标，去伪影效果评价方法仍有待检验，根据文献[25]提出去伪影的评价方法，测试的数据集必须是复杂的真实场景，包含快速变化的运动、剧烈的变形运动、手持摄像设备、各类遮蔽物体、剧烈与微小的位移运动、静态场景等。因此，本文主要从融合后的图像效果，如去伪影效果、细节效果等来对仿真结果进行评价。

图 12(b)—(e)分别为各个算法的融合结果图以及局部细节放大图，从中可观察到文献[1]的方法融合图像存在明显的伪影现象。在低秩矩阵恢复的融合结果中，RPCA并不能很好地抑制伪影，而PSSV以及本文方法都能够抑制伪影的出现且有较好的主观视觉效果。从地砖上的细节来看，RPCA损失的细节最多，其余算法都能较好地保留真实场景中的细节。基于低秩矩阵恢复的多曝光图像融合不仅能消除动态场景的伪影，并且在一定程度上也提高了融合图像的质量。

本文将多曝光HDR图像融合问题与低秩最优化理论相结合，提出了一种改进的多曝光图像HDR融合方法，该方法能消除动态场景的伪影，提高了融合图像的质量，有较好的主观视觉效果，比传统融合方法具有更好的鲁棒性。

仿真结果表明，该方法在较少的数据情况下恢复的低秩结果比RPCA的要好。RPCA需要足够多的曝光图像才能保证矩阵的低秩，但大多数情况下无法通过RPCA求解得到低秩的最优解，而本文提出的算法适用于多曝光HDR图像融合数据集较少的情况；PSSV以部分奇异值作为优化目标函数，虽然能在较少数据的情况下得到较好的低秩矩阵恢复，但它只保证数据在其余奇异值向量上方差不大，并不能保证低秩部分的数据在奇异值向量上有最大方差。本文算法对不同奇异值设置自适应的惩罚因子，能在较低重构误差的前提下得到更加低秩的恢复矩阵。

本文方法依赖于相机响应函数，由于图像采集设备获取的数字图像都是非线性映射的结果，所以标定相机响应函数的准确性会影响融合的质量。因此，在研究基于低秩矩阵恢复的多曝光HDR图像融合方法中应该弱化相机响应函数的条件，同时使得求解的结果更加精确。此外，由于多曝光图像序列不一定是静止的，可能存在因采集设备抖动导致图像位移，或者场景中存在运动的物体。如果不消除这类运动，会导致融合后的HDR图像出现伪影与模糊。

因此，下一步拟对动态场景融合过程中产生的鬼影(ghost)问题进行去伪影的融合处理，以进一步得到更多亮、暗区域的细节信息并且避免失真的融合图像。