随着互联网的普及、信息高速传输的发展以及智能手机的广泛应用，数字信息处理的需求与日俱增。其中，数字图像以其信息量大、传输速度快、作用距离远等优点成为人们获取信息的重要方式。数字图像在获取、传输以及存储的过程中，往往会因为各种干扰而产生噪声，从而影响图像信息的完整性，因此图像噪声的去除是图像预处理的重要一步。图像噪声主要分为加性噪声与乘性噪声，目前去除图像加性高斯噪声的模型已经比较成熟。在去除图像乘性噪声方面，文献[1-3]提出了几种经典的模型，但是上述模型大都假设乘性噪声服从Gamma分布。一般情况下，图像噪声与图像信息是相关的，例如激光雷达成像、CCD(charge-coupled device)固态光电检测器阵列、卫星遥感成像、医学CT成像等，都是以光量子计数成像的系统，在图像的获取中受到干扰基本都是量子噪声。量子噪声通常服从泊松分布，而不是Gamma分布，并且泊松噪声的强度通常依赖于原始图像信息，从分布的角度看，亮度大的像素受到的干扰更强。相比于一般乘性噪声而言，图像泊松噪声的去除是一个更难解决的问题。

经典的图像泊松噪声去除方法通常是将泊松噪声近似地转化为高斯噪声，然后在空域或者变换域进行。文献[4]与文献[5]分别利用维纳滤波和小波阈值收缩将泊松噪声变换为近似的高斯噪声，此类方法可以在一定程度上去除泊松噪声，但其去噪效果还有待提高。文献[6]中Rudin等人利用BV(bounded variation)空间在图像边缘不连续的特性，引入能量函数，提出了全新的基于全变分(TV)的图像去噪模型，将图像去噪转化成求解极小值的问题，并取得了较好的去噪效果。文献[7]在Bayesian-MAP框架下，使用泊松似然函数的负对数作为去泊松噪声模型的保真项，并结合TV模型中的正则项，提出了去除泊松噪声的全变分模型。该模型将泊松噪声去除问题转化为二阶偏微分方程稳定状态的求解，从而达到去除泊松噪声的目的。该算法虽然能够较好地去除泊松噪声，但它使用的TV正则项属于BV函数空间，在去除噪声的同时容易产生“阶梯效应”。为了解决这一问题，文献[8]将全变分去泊松噪声模型正则项中的微分由一阶变为二阶，该方法能较好地抑制“阶梯效应”，但同时也增加了数值求解的计算量。为了减小数值求解的复杂度，提高模型的去噪效率，文献[9]采用了一种简单高效的交替迭代投影方法来求解全变分去泊松噪声模型。该方法相比于传统的数值解法收敛速度更快，虽然有效地降低了模型的数值计算量，但是去噪效果不是很理想。近年来，新的微分方程去噪方法被相继提出，文献[10]将图像纹理半范数的倒数作为保真项的系数来构造新的微分去噪方程。文献[11]提出了一种带有时间积分的微分方程方法，在去噪效果上，该方法相比于传统微分方法有一定程度的提高，但该方法仍然属于微分方程方法的范围，并没有从根本上解决传统模型存在“阶梯效应”和边缘模糊的问题。

针对上述问题，文中提出了一种基于分数阶的非凸自适应去除泊松噪声新模型。其中非凸自适应正则项相比于原始正则项，能够自适应地根据不同像素点的信息选取不同的参数以更好地保持图像边缘信息。对于传统全变分模型容易产生的“阶梯效应”的问题，文中引入分数阶离散微分向量来替代原去噪模型的一阶离散微分向量。目前分数阶微分已经成功地运用到了图像去加性噪声模型中^[12-14]，克服了原始正则项建立在BV空间的这一缺点，能够有效地抑制“阶梯效应”。实验结果表明，文中提出的自适应分数阶去除泊松噪声模型在去噪的同时能够较好地保护图像边缘和抑制“阶梯效应”。由于新模型正则项的非凸性，使用传统的数值解法存在运算量大、收敛速度慢的缺点，文中结合交替迭代法与加权原始-对偶算法^[15]，提出了一种高效的数值解法。实验结果表明，相比传统的数值解法，文中提出的数值解法不仅收敛速度较快，而且复杂度较低。

假定被噪声污染的图像为$f$，原始无噪图像为$u$，图像的紧支撑域为$\Omega$，并且设噪声服从泊松分布，即

$ {P_\mu }\left( \eta \right) = \frac{{{\rm e^{-\delta }}{\mu ^\eta }}}{{\eta !{\rm{ }}}},\eta \ge 0 $

(1)

式中，$\delta$为泊松分布的均值，根据泊松分布的特点，其方差与均值相等。

根据对无噪图像$u$的最大后验估计分析，由贝叶斯定理可从观察到的图像$f$中估计出原始无噪图像$u$。

$ P\left( {u|f} \right) = {\rm{ }}\frac{{P\left( {f|{\rm{ }}u} \right)P\left( u \right){\rm{ }}}}{{P\left( f \right)}} $

(2)

式中，由于$f$是已观察到的图像，因此$P(f)$通常为一固定值，所以最大化$P(u|f)$的问题等效于最大化$P(f|u)P(u)$。在泊松噪声环境下，对于每一个像素$x_i∈\Omega$，都有

$ P\left( {f\left( {{x_i}} \right)|u} \right) = {P_{u({x_i})}}\left( {f\left( x \right)} \right) = \frac{{{{\rm{e}}^{ - u({x_i})}}u{{\left( x \right)}^{f({x_i})}}}}{{f({x_i})!}} $

(3)

假设在图像紧支撑域$\Omega$中像素之间都是彼此相互独立的，由式(3) 可得

$ P\left( {f|u} \right) = \prod\limits_i {\frac{{{{\rm {e}}^{ - u({x_i})}}u{{({x_i})}^{f({x_i})}}}}{{f({x_i})}}} {\rm{ }} $

(4)

通常，原始图像$u$的先验分布服从吉布斯分布，即

$ P\left( u \right) = \prod\limits_i {\lambda {\rm{exp}}( - \lambda \left| {\nabla {u_i}} \right|)} $

(5)

式中，$\lambda$为参数。

因此，极大化$P(f|u)P(u)$问题可以转换为极小化问题$ - \ln (P(f|u)P(u))$，将式(4)(5) 代入其中，可得

$ \mathop {\min }\limits_{u \in \Omega } = \sum\limits_i {\{ - \ln \lambda - \ln ({{\rm {e}}^{ - u({x_i})}}u{{({x_i})}^{f({x_i})}}({{\rm {e}}^{ - \lambda \left| {\nabla u} \right|}}))\} } $

(6)

忽略常数项$\ln \lambda$，式(6) 可进一步转化为

$ \mathop {\min }\limits_{u \in \Omega } = \sum\limits_i {\left( {u\left( {{x_i}} \right) - f\left( {{x_i}} \right)\ln u\left( {{x_i}} \right)} \right) + \lambda \sum\limits_i {\left| {\nabla u} \right|} } $

(7)

由式(7) 可得文献[7]中离散化后的传统去泊松噪声全变分模型，即

$ E\left( u \right) = \mathop {\min }\limits_{u \in \Omega } \left\| {u - f\ln u} \right\|_2^2 + \lambda \left\| {\nabla u} \right\|{_1} $

(8)

式中，$\left\| {u - f\ln u} \right\|_2^2$为保真项，主要起保留原图像特性和降低图像失真度的作用。$\lambda {\left\| {\nabla u} \right\|_1}$为正则项，起到平衡去噪与平滑的作用，$\lambda$为正则化参数。

模型式(8) 中的正则项是建立在BV空间上的函数，而BV空间上的函数特性会造成分段平滑，所以模型的稳态解中始终无法解决“阶梯效应”的问题。另外，由吉布斯分布推导得到的正则项参数$\lambda$是一个固定的值，不能根据像素点的变动而改变，因此在保护图像的边缘信息上存在一定的缺陷。

如上所述，式(5) 中原始图像$u$的先验分布服从吉布斯分布，参数$\lambda$是一个固定的值，即整幅图像的不同像素点对应的都是一个相同的参数$\lambda$。由于一幅图像中不同元素之间存在着一定的关联性，因此，将式(5) 中的参数$\lambda$固定为一个常数显然是不尽合理的。考虑到图像像素之间的关联性，文中假设图像中不同像素点对应不同的参数$\lambda_i$。因此，式(5) 可以改写为

$ P\left( u \right) = \prod\limits_i {{\lambda _i}{\rm{exp}}( - {\lambda _i}\left| {\nabla {u_i}} \right|)} $

(9)

式中，$i=1，2，…，n$，$n$为图像像素的个数。将式(9) 与式(4) 代入极小化问题$-\ln (P(f|u)P(u))$中，可得

$ \mathop {\min }\limits_{u \in \Omega } = \sum\limits_i {\left\{ { - \ln{\lambda _i} - \ln \left( {{{\rm {e}}^{ - u\left( {{x_i}} \right)}}u{{\left( {{x_i}} \right)}^{f({x_i})}}\left( {{{\rm {e}}^{ - {\lambda _i}\Delta u}}} \right)} \right)} \right\}} $

(10)

式(10) 可进一步简化为

$ \mathop {\min }\limits_{u \in \Omega } = \sum\limits_i {\left( {u\left( {{x_i}} \right) - f\left( {{x_i}} \right)\ln u\left( {{x_i}} \right)} \right)} + \lambda \sum\limits_i {\left| {\nabla u} \right|} - \ln {\lambda _i} $

(11)

由式(11) 可得离散化后的自适应去泊松噪声全变分模型

$ E(u,{\lambda _i}) = \mathop {\min }\limits_{u \in \Omega } \left\| {u - f\ln u} \right\|_2^2 + {\lambda _i}\left\| {\nabla u} \right\|{_1} - \ln {\lambda _i} $

(12)

对$\lambda_i$求导，可得

$ \left| {\nabla u} \right| - \frac{1}{{{\lambda _i}}} = 0 \Rightarrow {\lambda _i} = \frac{1}{{\left| {\nabla u} \right|}} $

(13)

为了避免分母为零的情况出现，$\lambda_i$可取为

$ {\lambda _i}{\rm{ = }}\frac{1}{{1 + k\left| {\nabla u} \right|}} $

(14)

式中，参数$k>0$。

使用式(14) 中的$\lambda_i$替换文献[7]中的正则项参数，可以得到如下新的自适应去除泊松噪声的全变分模型

$ \mathop {\min }\limits_{u \in \Omega } \left\| {u - f\ln u} \right\|_2^2 + g\left( u \right)\left\| {\nabla u} \right\|{_1} $

(15)

式中，$g\left( u \right){\rm{ = 1/}}\left( {1 + k\left| {\nabla u} \right|} \right)$。

新模型(15) 中的自适应系数$g(u)$能够在图像的不同区域自适应的选取不同系数，从而得到更好的保边效果。然而在式(15) 的新模型中，正则项的系数虽然是自适应的，但是正则项中的$\nabla u$是一阶离散微分向量，其正则项仍然是建立在BV空间上的函数，所以新模型的稳态解中始终会存在“阶梯效应”。并且一阶离散微分向量$\nabla u$只能结合图像中两个像素点的信息，具有明显的局限性。由于图像像素之间存在着一定的关联性，因此在设计新模型的时候应当更充分地结合图像的局部信息。式(15) 中的一阶离散微分向量是局部算子，显然，它不能较好地满足这一要求。

近年来，分数阶微分理论在图像处理中的应用得到了人们的关注。分数阶离散微分向量是全局算子，它特有的“弱导数”性质可以结合更多的局部信息，能够较好地抑制“阶梯效应”。因此，文中将新模型中的一阶离散微分向量用分数阶离散微分向量替换，从而更好地抑制“阶梯效应”和提升保护图像边缘的能力。文中选择图像处理领域应用较多的Grünwald -letnikov(G-L)^[16]分数阶离散微分向量来构造新的正则项。对于大小为$M×N$的图像$u$，新模型中离散形式的分数阶离散微分向量为

$ {\nabla ^\alpha }u = \left( {{{\left( {\Delta _x^\alpha u} \right)}_{i,j}},{{\left( {\Delta _y^\alpha u} \right)}_{i,j}}} \right)\;i = 1,2, \ldots ,M;j = 1,2, \ldots ,N $

(16)

式中，$α$是分数阶的阶数，$1<α<2$，且

$ \left\{ \begin{array}{l} {(\Delta _x^\alpha u)_{i,j}} = {\left( { - 1} \right)^k}\sum\limits_{k = 0}^{K - 1} {C_k^\alpha {u_{i - k,j}}} \\ {(\Delta _y^\alpha u)_{i,j}} = {\left( { - 1} \right)^k}\sum\limits_{k = 0}^{K - 1} {C_k^\alpha {u_{i,j - k}}} \end{array} \right. $

(17)

式中，$K≥3$为整数，$C_k^\alpha = {\rm{ }}\frac{{\mathit{\Gamma }\left( {\alpha + 1} \right)}}{{\mathit{\Gamma }\left( {k + 1} \right)\mathit{\Gamma }\left( {\alpha - k + 1} \right)}}{\rm{ }}$，$\mathit{\Gamma }$(·)表示Gamma函数。

基于上述考虑，文中提出了一种基于分数阶变分的自适应去泊松噪声新模型，即

$ \mathop {\min }\limits_{u \in \Omega } \left\| {u - f\ln u} \right\|_2^2 + g\left( u \right)\left\| {{\nabla ^\alpha }u} \right\|{_1} $

(18)

式中

$ \left\{ {\begin{array}{*{20}{l}} {{{\left\| {{\nabla ^\alpha }u} \right\|}_1} = \sum\limits_{\begin{array}{*{20}{c}} {1 \le i \le M}\\ {1 \le j \le N{\rm{l}}} \end{array}} {\left| {{{\left( {{\nabla ^\alpha }u} \right)}_{i,j}}} \right|} }\\ {\left| {{{({\nabla ^\alpha }u)}_{i,j}}} \right| = \sqrt {{{({{(\Delta _x^\alpha )}_{i,j}})}^2} + {{({{(\Delta _y^\alpha )}_{i,j}})}^2}} } \end{array}} \right. $

(19)

另外，对于自适应系数$g\left( u \right){\rm{ = 1/}}\left( {1 + k\left| {{\nabla ^\alpha }u} \right|} \right){\rm{ }}$，在图像的光滑区域${{\nabla ^\alpha }u}$值较小，$g(u)$值较大，其惩罚力度加大，噪声将得到更好的抑制。在图像的边缘区域${{\nabla ^\alpha }u}$值较大，$g(u)$值较小，其惩罚力度减小，图像边缘将得到更好的保持。

由上可知，新模型的自适应系数$g(u)$会根据不同的区域自适应地选取较优的惩罚力度，在去噪与边缘保护之间得到一个合适的阈值。另外，新模型中采用分数阶离散微分向量替代TV模型中的一阶离散微分向量，能有效地抑制“阶梯效应”，从而取得更好的去噪效果。

由于新模型的正则项为非凸函数，因此很难使用经典的欧拉-拉格朗日方程方法进行求解。另外，传统的偏微分方程方法求解变分去噪模型时存在着收敛速度慢、运行时间长等缺点。为了能够快速有效地求解新模型，文中结合交替迭代法和加权原始-对偶法，提出了一种高效的数值解法。

参考文献[17]的方法，文中引入新的辅助变量${\rm{ }}v \in {\boldsymbol{{\rm{R}}}^n}$将式(18) 的无约束问题转换成等价的有约束问题，等价模型如下

$ \left\{ \begin{array}{l} \mathop {\min }\limits_{u,v} \left\| {u - f\ln u} \right\|_2^2 + g\left( v \right)\left\| {{\nabla ^\alpha }v} \right\|{_1}\\ s.t.{\nabla ^\alpha }u = {\nabla ^\alpha }v \end{array} \right. $

(20)

对于式(20)，进一步引入容许变量${\rm{ }}z \in {\boldsymbol{{\rm{R}}}^n}$，将式(20) 的有约束项问题转换成求解如下最小值的问题

$ \mathop {\min }\limits_{u,v} \left\| {u - f\ln u} \right\|_2^2 + g\left( v \right){\left\| {\nabla v} \right\|_1} + L{\left\| {u - v + z} \right\|^2} $

(21)

$ {z^k} = {z^{k - 1}} + {u^k} - {v^k} $

(22)

对于式(21)，文中采用交替迭代方法进行求解，即先固定$v$求解$u$，然后再固定$u$求解$v$。为此将式(21) 分解成下面两个模型

$ \hat u = \mathop {\arg \min }\limits_u \left\| {u - f\ln u} \right\|_2^2 + L{\left\| {u - v + z} \right\|^2} $

(23)

$ \hat v = \mathop {\arg \min }\limits_v g\left( v \right)\left\| {{\nabla ^\alpha }v} \right\|{_1} + L{\left\| {u - v + z} \right\|^2} $

(24)

对式(23) 中的$u$进行求导，可得方程

$ \left(1 - \frac{f}{u} \right)+ 2L\left( {u - v + z} \right) = 0 $

(25)

对于式(25) 可以采用牛顿迭代法进行求解，即

$ {u^{k + 1}} = {u^k} - {\rm{ }}\frac{{\left( {1 - \frac{f}{{{u^k} + \varepsilon }}} \right) + 2L({u^k} - {v^k} + {z^k})}}{{f/{{({u^k} + \varepsilon )}^2} + 2L}} $

(26)

在式(26) 中，为避免分母为零，添加极小值$\varepsilon = 10{{\rm{e}}^{ - 5}}$。

对于式(24)，其与鞍点结构优化模型的原始模型在形式上较为相似，故可采用原始-对偶方法进行求解。原始-对偶方法已经广泛应用于图像处理的各个领域，它具有求解容易、收敛快等特点。具有鞍点结构的优化模型为

$ \mathop {\min }\limits_{x \in \boldsymbol{X}} \mathop {\max }\limits_{p \in \boldsymbol{Y}} \left\langle {Ax,p} \right\rangle + G\left( x \right) - {F^*}\left( p \right) $

(27)

式中，$\boldsymbol{X},\boldsymbol{Y}$表示有限维度实向量空间，$\left\langle { \cdot , \cdot } \right\rangle $表示内积，$A$表示任意的线性算子，$G$和$F$表示函数，$F^*$代表$F$的拓扑对偶。在式(27) 中，可以把$p$看做对偶变量，$x$看做原始变量，于是具有鞍点结构的优化模型可以描述为原始模型与对偶模型的原始-对偶模型。其中，原始模型为

$ \mathop {\min }\limits_{x \in X} F\left( {Ax} \right) + G\left( x \right) $

(28)

对偶模型为

$ \mathop {\max }\limits_{p \in Y} - \left( {{G^*}\left( { - {A^*}p} \right)} \right) + {F^*}\left( p \right) $

(29)

由上面的分析可以得到模型式(24) 与原始模型式(28) 的映射关系，即

$ \left\{ \begin{array}{l} Ax = g\left( v \right){\left| {{\nabla ^\alpha }v} \right|_1}\\ G\left( u \right) = L{\left\| {u - v + z} \right\|^2} \end{array} \right. $

(30)

对于模型式(27) 的求解，文献[15]中给出了一种基于预解式的原始—对偶求解方法。其中基于预解式的求解是先令$x$固定，对$p$进行求导，然后固定$p$对$x$进行求导。可以得到预解式

$ \left\{ \begin{array}{l} p = {\left( {I + \nabla {F^*}} \right)^{ - 1}}\left( {p + Ax} \right)\\ x = {\left( {I + \nabla G} \right)^{ - 1}}\left( {x - Ap} \right) \end{array} \right. $

(31)

根据原始—对偶算法，当$F^*$和$G$中至少有一个是凸函数的时候，基于预解式的原始-对偶算法求解式(27)，即

$ \left\{ \begin{array}{l} {p^{k + 1}} = {\left( {I + \delta {F^*}} \right)^{ - 1}}\left( {{p^k} + \delta A{{\overline x }^k}} \right)\\ {x^{k + 1}} = {\left( {I + \tau G} \right)^{ - 1}}\left( {{x^k} - \tau {A^*}{p^{k + 1}}} \right)\\ {\overline x ^{k + 1}} = {\overline x ^{k + 1}} - {x^k} \end{array} \right. $

(32)

式中，参数$\delta ,\tau > 0$。

式(24) 模型中的$L{\left\| {u - v + z} \right\|^2}$为凸函数，满足原始-对偶算法的前提条件，因此使用原始-对偶算法求解式(24) 是完全合理的。

利用原始-对偶算法求解式(24)，需先求出预解算子${{{\left( {I + \nabla {F^*}} \right)}^{ - 1}}\left( {p + Ax} \right)}$、${{{\left( {I + \nabla G} \right)}^{ - 1}}\left( {x - Ap} \right)}$与线性算子$A$。由于${F^*}\left( p \right) = F_p^*\left( p \right)$、$G\left( u \right) = L{\left\| {u - v + z} \right\|^2}$，于是可以得到

$ \left\{ {\begin{array}{*{20}{l}} {p = {{\left( {I + \delta \nabla {F^*}} \right)}^{ - 1}}\left( {\tilde p} \right) \Leftrightarrow }\\ {p = \tilde p/\max \left( {1,\frac{{\left| {\tilde p} \right|}}{{g\left( v \right)}}} \right)} \end{array}} \right. $

(33)

$ \left\{ {\begin{array}{*{20}{l}} {v = {{\left( {I + \tau \nabla G} \right)}^{ - 1}}\left( {\tilde v} \right) \Leftrightarrow }\\ {v = \left( {\tilde v + {\tau _n}\lambda \left( {v + r} \right)} \right)/\left( {1 + {\tau _n}\lambda } \right)} \end{array}} \right. $

(34)

式中，$\tilde p = p + \delta A\bar v$，$\tilde v = v - \tau {A^*}p$，$A = \overline {{{\left( { - 1} \right)}^\alpha }} {\rm{di}}{{\rm{v}}^\alpha }$。

综上所述，文中给出利用原始—对偶方法求解式(24) 的具体算法(算法1) 如下：

1) 初始化$g=u+z^0，v^0=f，p^0=0，\delta，α，τ，\varepsilon，L>0。$

2) 迭代

$ \left\{ \begin{array}{l} {p^{k + 1}} = pro{j_p}\left( {{p^k} + \delta {\nabla ^\alpha }{{\overline v }^k}} \right)\\ {v^{k + 1}} = \frac{{{v^k} + \tau {\rm{di}}{{\rm{v}}^\alpha }\left( {{p^{k + 1}}} \right) + \tau Lg}}{{1 + \tau L}}\\ {\overline v ^{k + 1}} = 2{v^{k + 1}} - {v^k} \end{array} \right. $

(35)

式中，$pro{j_p}\left( {\tilde p} \right) = \frac{{{{\tilde p}^k}}}{{\max \left( {1,\frac{{\left| {{{\tilde p}^k}} \right|}}{{g\left( v \right)}}} \right)}}$，$\rm{div}^α$表示${{\nabla ^\alpha }}$的负共轭。

3) 当$\left\| {{u^{k + 1}} - {u^k}} \right\|_2^2/\left\| {{u^k}} \right\|_2^2 \le \varepsilon $时，结束迭代。

结合式(26) 与式(24) 的求解过程，文中给出非凸分数阶的泊松去噪模型的具体算法(算法2) 如下：

1) 初始化$z^0=0，u^0=v^0=f，p^0=0，\delta，τ，\varepsilon，L>0；$

2) 重复迭代；

3) 应用牛顿迭代求解式(25)；

4) 应用算法1求解式(24)；

5) 当$\left\| {{u^{k + 1}} - {u^k}} \right\|_2^2/\left\| {{u^k}} \right\|_2^2 \le \varepsilon $时，结束迭代。

为了验证新模型及新数值算法的有效性，文中从泊松噪声去除模型和数值算法两个方面进行了相关测试。所有程序均在MATLAB R2012b上编程实现，CPU为i5-3210M 2.5 GHz。如图 1所示，文中选取Lena、SynImage、Peppers与Phantom作为测试图像，图像的大小均为256×256像素。

模型分别与TV模型、自适应TV模型进行比较以验证文中新模型的有效性。实验中，选取图 1所示的4幅图像作为测试图像，并加入泊松噪声。为了客观评价实验效果，文中用峰值信噪比(PSNR)、均方误差(MSE)、信噪比(SNR)及图像相似度(SSIM)等指标来衡量去噪后的结果图像。

$ {\rm{SNR}} = 101g\left( {\frac{{\left\| {{u_0} - {\mu _u}} \right\|_2^2}}{{\left\| {{u_0} - u} \right\|_2^2}}} \right) $

(36)

$ {\rm{MSE}} = \int_\Omega {\left( {u - {u_0}} \right){\rm{d}}x} $

(37)

$ {\rm{PSNR = }}101g\left( {\frac{{{n^2}}}{{\left\| {u - {u_0}} \right\|_2^2}}} \right) $

(38)

$ {\rm{SSIM}} = \frac{{\left( {2{\sigma _{uu_0}} + c1} \right)\left( {2{\mu _u}{\mu _{{u_0}}} + c2} \right)}}{{\left( {\sigma _u^2 + \sigma _{{u_0}}^2 + c1} \right)\left( {\mu _u^2 + \mu _{{u_0}}^2 + c2} \right)}} $

(39)

式中，$u$、$u_0$分别是原始图像和处理后的图像，$σ_u$、$σ_{u0}$、$μ_u$、$μ_{u0}$分别为$u$、$u_0$的标准差和均值，$σ_{uu_0}$为$u$与$u_0$的协方差。参考文献[18]，文中将SSIM中的维稳系数$c_1$，$c_2$分别设定为0.002与0.001。在实验中，文中新模型的参数设定为：$τ=0.03$，$\delta=0.3$，$k=2.0$，$L=7.5$，$α=1.3$。

图 2所示是含噪图像及不同模型的去噪结果，从上至下依次为Lena、SynImage、Peppers、Phantom。图 2 (a)是被泊松噪声污染的图像，图 2 (b)是用经典的TV模型去噪后的结果图像，图 2(c)是自适应TV模型去噪后的结果图像，其中自适应系数使用的是式(15) 中使用的系数$g(u)$，图 2 (d)是使用文中新模型去噪后的结果图像。从图 2(b)的TV模型去噪实验结果可列以看出，其去噪结果图像存在着图像边缘模糊现象。另外，在图像的平坦区域也出现了“阶梯效应”，尤其以Lena图像最为明显。从图 2(c)的自适应TV模型去噪实验结果可列以看出，其去噪结果图像的边缘部分得到了一定的保持，但是在部分图像边缘处的噪声也被当作边缘被保护起来，从而形成了锯齿状的边缘。从图 2 (d)的本文所提模型去噪实验结果图像可以看出，图像细节较为清晰，边缘部分得到了较好的保持，且在图像的平坦区域较好地抑制了“阶梯效应”。

表 1所示是用TV模型、自适应TV模型及文中模型对含泊松噪声的不同图像进行去噪后的实验结果。从表 1可以看出，对4幅典型图像的去除泊松噪声结果中，文中提出的新模型无论是在峰值信噪比、均方误差、信噪比还是图形相似度上都明显优于传统的TV模型和自适应TV模型。

选取Lena图像作为测试图像，图 3所示是不同模型的迭代次数与SNR、MSE、PSNR、SSIM的关系曲线比较结果。从图 3所示的实验结果可以看出文中提出的新的去除泊松噪声的模型在SNR与MSE、PSNR、SSIM曲线中，明显优于传统的TV与自适应TV模型。其中从图 3(d)中可以明显看出，随着迭代次数的增加，传统的TV模型会造成图像的模糊，从而导致图像结构相似度的大幅度下滑。

图 4所示是不同模型去噪结果的局部3维维图。图 4(a)为未被噪声污染的原始图像，为了更形象地说明不同模型的性能，文中以红线所围区域的局部3维图来描述TV模型、自适应TV模型和文中新模型的去噪效果。图 4(b)是未被噪声污染的原始区域3维图，图 4(c)是被噪声污染的原始区域3维图，图 4(d)是经过TV模型处理之后的3维图，图 4(e)是经过自适应TV模型处理之后的3维图，图 4(f)是经过本文模型处理之后的3维图，从图 4(d)中可以看出TV模型的处理结果出现了明显的过度平滑，并且在3维图像的顶部，出现了“阶梯效应”。从图 4(e)中可以看出自适应模型处理后的图在顶部与底部之间的交界处的边缘位置，出现明显的参差不齐，误将一些噪声当成图像的边缘保留了下来，使得噪声不能被完全去除。而从图 4 (f)中的本文模型3维局部效果图中可以看出，在顶部与底部平滑区域无明显的“阶梯效应”，且在底部与顶部之间的边缘区域颜色层次分明，并没有出现自适应TV模型中的参差不齐的现象，从图 4的局部3维图中可以看出，本文模型在边缘保护和抑制“阶梯效应”上均很好的效果。

为了验证本文数值解法的快速优势，将与其他一些经典的算法进行比较，分别是原始偏微分方程算法^[7]和chambolle投影算法^[19]。实验中，选取Lena图像作为测试图像，并加入泊松噪声，在表 2中，给出了当$k=0，α=1.0$，解的均方根误差$\varepsilon≤10^{-3}，10^{-4}$时几种算法的迭代次数与CPU运行时间的比较。从表 2中可以看出，文中采用的交替迭代算法和加权原始-对偶算法在运行速度上明显优于文献[7]与文献[19]的算法。

图 5所示是不同算法的迭代次数与均方根误差的关系曲线比较结果。从图 5中可以看出，在迭代次数为8时，本文算法的均方误差曲线趋于平稳，文献[7]的均方误差曲线趋于稳定时的迭代次数为45，而文献[19]的均方误差曲线趋于稳定时的迭代次数为21。实验结果表明，本文中算法能有效快速地收敛于鞍点。

在去除泊松噪声的同时，为了更好的保留图像的边缘信息以及有效抑制“阶梯效应”，提出了一种分数阶自适应去泊松噪声新模型与数值解法。针对新的模型的求解，文中采用交替迭代与加权原始-对偶相结合的全新数值解法，新模型相比与传统模型，能更好的保留边缘信息与抑制“阶梯效应”的效果，同时新的数值解法想比于传统的数值解法，具有收敛速度快、复杂度低等优点。此外，通过对多幅图像的实验数据结果表明，文中的新模型与新的数值解法具有一定的普适性。

虽然新的模型与新的数值求解算法可以有效的抑制“阶梯效应”与降低模型求解的复杂度，但是模型中分数阶的阶次是固定的值，众所周知，图像的不同区域使用不同的分数阶阶次会有不同的效果，对整幅图像来说使用固定的分数阶阶次处理的效果还有进一步提高的空间，下一步的研究中将考虑如何做到分数阶阶次的自适应，以达到最好的处理效果。