在图像去噪问题中，由于高斯噪声比较常见，大多数图像去噪算法都假设噪声服从高斯概率模型. 在天文成像、计算机断层扫描、荧光共焦显微镜成像、CCD固态光电检测等光子计数成像系统中，噪声的概率分布服从泊松概率模型^[1]，利用高斯噪声模型不能有效描述其概率分布特征，采用面向高斯噪声的去噪方法难以获得理想的去噪效果.

图像泊松噪声去除的经典方法是通过方差变换将泊松噪声转化为近似高斯分布的噪声^[2]，然后利用成熟的高斯去噪方法，如维纳滤波、全变差模型^[3]、非局部均值^[4]等进行去噪处理，最后经逆方差变换得到去噪后的图像. 这种方法在一定程度上可以去除泊松噪声，但是去噪效果仍然有待提高. 另一类方法是基于贝叶斯最大后验（maximum a posteriori，MAP）^[5]概率的正则化方法，这类方法通常以Richardson-Lucy（RL）算法^[6]为基础. 为了解决RL算法多次迭代时产生的噪声放大问题，研究者提出了很多基于变分正则化的改进方法，常用的有吉洪诺夫（Tikhonov）正则项、全变分（total variation，TV）正则项、高阶正则项等^[7]. Le等^[8]选用TV正则项推导出了经典的泊松全变分去噪模型. 但是，TV模型在去除噪声的同时容易产生阶梯效应，造成部分边缘和细节丢失. 为了解决这一问题，Gilboa等^[9]引入了空间变换；张哲等^[10]提出了四阶去噪模型；白键等^[11]提出了基于积分微分方程的改进方法；胡学刚等^[12]提出了一种基于分数阶导数的泊松去噪算法；张峥嵘等^[13]提出了一种基于非局部全变差正则化的图像去噪模型.

马尔科夫随机场（Markov random fields，MRF）^[14]也是比较有代表性的一种正则项模型. 然而，传统的马尔科夫随机场是基于简单的人为设定的滤波器进行构建的，虽然其实现速度较快，但是对于复杂自然图像而言，模型的准确度欠佳^[15]，会对去噪效果产生影响.

为了构建能够准确描述自然图像先验概率分布特征的MRF模型，Roth等^[16]采用机器学习的方法，建立了一种全新的高阶MRF模型，称之为专家场（fields of experts，FoE）模型，其特点在于构成FoE模型的所有滤波器都是基于自然图像库学习得到的. 因此，FoE模型能够更准确地捕捉复杂自然图像的先验概率分布特点. FoE模型已被应用于一系列视觉图像处理领域，如图像修复、模糊图像复原等，有广阔的研究和应用前景.

本文提出一种基于泊松噪声模型和FoE约束的正则化去噪方法，采用二次惩罚函数对代价函数进行优化求解，并通过仿真实验以及图像评价指标对所提方法的去噪效果进行客观评价.

FoE模型的本质是一种马尔科夫随机场，其邻域系统被定义为包含了K个 $m \times m$的高通滤波器模板 ${{{J}}_k}\;( = 1,2,\cdots,K)$的集合，区域内所有节点都是相邻节点. 假设FoE模型邻域系统中所有的 ${{{J}}_k}$均对应如下式所示的相同势函数：

(1) $V\left( {{{{f}}_i}} \right) = \prod\limits_{k = 1}^K {{\phi _k}\left( {{{J}}_k^{\rm{T}}{{{f}}_i};{{\theta}} } \right)}. $

式中： $K$为邻域系统中滤波器模板总数， ${{{f}}_i}$为 ${{{J}}_k}$遍历图像 ${{f}}$时覆盖的第 $i$个图像块， ${\phi _k}$为专家函数， ${{\theta}} $为专家函数参数集.

根据Hammersley-Clifford理论^[17]，整幅图像在FoE模型下的概率分布可以写成下式：

(2) $P\left( {{f}} \right) = \frac{1}{Z}\prod\limits_{i = 1}^N {V\left( {{{{f}}_i}} \right)} = \frac{1}{Z}\prod\limits_{i = 1}^N {\prod\limits_{k = 1}^K {{\phi _k}\left( {{{J}}_k^{\rm{T}}{{{f}}_i};{{\theta }}} \right)} }. $

式中：N为图像像素总数，Z为归一化函数.

在FoE模型中， ${\phi _k}$采用student-t概率密度函数：

(3) ${\phi _k}\left( {{{J}}_k^{\rm{T}}{{{f}}_i};{{\theta }}} \right) = {\left[ {1 + \frac{1}{2}{{\left( {{{J}}_k^{\rm{T}}{{{f}}_i}} \right)}^2}} \right]^{ - {\alpha _k}}}.$

式中： ${\alpha _k}$为滤波器参数.

在实际问题中，通常采用吉布斯分布形式对式（2）所示FoE模型进行表示：

(4) $P\left( {{f}} \right) = \frac{1}{Z}\exp\; \left[ {\sum\limits_{i = 1}^N {\sum\limits_{k = 1}^K {\ln\; {\phi _k}\left( {{{J}}_k^{\rm{T}}{{{f}}_i};{{\theta }}} \right)} } } \right].$

邻域系统 ${{{J}}_k}\;(k = 1,2,\cdots,K)$由Roth等^[16]采用对比散度（contrast divergence，CD）算法从自然图像库中学习得到. 由于CD算法的随机性，每次学习的结果不尽相同，但是实验发现，CD学习的结果对于随机性参数的确切值不是很敏感，其典型结果^[16]如图1所示. 图中，每个 $3 \times 3$大小的方格代表一个滤波器，每个小格的灰度代表该点的数值大小， $\alpha $为该滤波器student-t函数的指数.

噪声图像的退化模型可以表示为

(5) ${{g}} = {{f}} + {{n}}.$

式中： ${{g}}$为噪声图像， ${{f}}$为清晰图像， ${{n}}$为噪声.

基于噪声的随机特性，根据贝叶斯定律可得

(6) $P({{f}}|{{g}}) = \frac{{P({{f}})P({{g}}|{{f}})}}{{P({{g}})}}.$

图像去噪问题可以表述为噪声图像已知条件下的最大后验概率估计问题：

(7) ${{f}} = \mathop {\max }\limits_{{f}} \;P({{f}}|{{g}}) = \mathop {\max }\limits_{{f}}\; P({{f}})P({{g}}|{{f}}).$

对式（7）取负自然对数，上述求最大值问题可转为求最小值问题：

(8) ${{f}} = \mathop {\min }\limits_{{f}}\; \left[ { - \ln\; P({{f}}) - \ln\; P({{g}}|{{f}})} \right].$

由式（8）可知，求解清晰图像 ${{f}}$的主要步骤如下：1）找到能反应图像固有概率分布特征的先验概率模型 $P\left( {{f}} \right)$；2）找到描述观测噪声概率分布特征的概率模型 $P\left( {{{g}}|{{f}}} \right)$. 构建图像的先验模型和噪声模型后将其代入式（8），得到相应的代价函数. 通过求解代价函数即可从噪声图像中恢复出清晰图像.

若 ${{{g}}_i}$代表噪声图像中的某个像素， ${{{f}}_i}$为清晰图像中对应的像素，并且噪声服从泊松概率模型，假设所有像素独立同分布，则对于噪声图像可得：

(9) $P\left( {{{g}}|{{f}}} \right) = \prod\limits_i {\frac{{{{f}}_i^{{{{g}}_i}}}}{{{{{g}}_i}!}}} \exp \;\left( { - {{{f}}_i}} \right).$

对式（9）进行负对数操作后有

(10) $ - \ln\; P\left( {{{g}}|{{f}}} \right) = \sum\limits_{i = 1}^N {\left[ {{{{f}}_i} - {{{g}}_i}\ln\; {{{f}}_i}} \right]} .$

对 $P\left( {{f}} \right)$的建模则采用FoE模型. 由式（4）可得

(11) $ - \ln\; P\left( {{f}} \right) = \sum\limits_{i = 1}^N {\sum\limits_{k = 1}^K {{\alpha _k}\ln \;\left[ {1 + \frac{1}{2}{{\left( {{{J}}_k^{\rm{T}}{{{f}}_i}} \right)}^2}} \right]} } .$

将式（10）和（11）代入式（8），可得代价函数：

(12) $\begin{split}{{f}} = &\mathop {\min }\limits_{{f}} \;\frac{\lambda }{2}\sum\limits_{i = 1}^N {\left( {{{{f}}_i} - {{{g}}_i}\ln\; {{{f}}_i}} \right)} + \\ &\sum\limits_{i = 1}^N {\sum\limits_{k = 1}^K {{\alpha _k}\ln\; \left[ {1 + \frac{1}{2}{{\left( {{{J}}_k^{\rm{T}}{{{f}}_i}} \right)}^2}} \right]} } + {l_ + }\left( {{f}} \right)\end{split}.$

式中： $\lambda $为正则系数， ${l_ + }\left( {{f}} \right)$为引入的非负性约束，

(13) ${l_ + }\left( {{f}} \right) = \left\{ \begin{aligned} &0,\quad{{f}} \in {\bf{R}}_{\rm{N}}^ + ;\\ &+ \infty ,\quad {\text{其他}}. \end{aligned} \right.$

采用二次惩罚函数对如式（12）所示的代价函数进行求解. 引入2组辅助变量： ${{u}}$ 和 ${w_k}\;(k = 1,2,\cdots,K)$. 其中， ${{u}}$对应于清晰图像 ${{f}}$， ${w_k}$对应于滤波器 ${{{J}}_k}$和图像 ${{f}}$的卷积，例如： ${w_1}$对应于 ${{J}}_1^{\rm{T}}{{f}}$，依此类推. 引入2个惩罚系数 ${\beta _1}$和 ${\beta _2}$，将式（12）中的问题转化为下式：

(14) $ \begin{split} \left( {{{f}},{w_k},{{u}}} \right) =& \mathop {\min }\limits_{\left( {{{f}},{w_k},{{u}}} \right)} \frac{\lambda }{2}\sum\limits_{i = 1}^N {\left( {{{{u}}_i} - {{{g}}_i}\ln\; {{{u}}_i}} \right)} + \frac{{{\beta _1}}}{2}{\left( {{{u}} - {{f}}} \right)^2} +\\ &\sum\limits_{i = 1}^N {\sum\limits_{k = 1}^K {{\alpha _k}\ln\; \left( {1 + \frac{1}{2}w_{ki}^2} \right) + } }\\ &\frac{{{\beta _2}}}{2}\sum\limits_{k = 1}^K {{\alpha _k}{{\left( {{w_k} - {{J}}_k^{\rm{T}}{{f}}} \right)}^2}} + {l_ + }\left( {{f}} \right).\\[-17pt] \end{split} $

根据二次惩罚函数原理，当 ${\beta _1}$、 ${\beta _2}$的值趋近于无穷大时，式（14）的解收敛至式（12）的解，如图2所示，图中， $n$为迭代次数， $e$为误差. 当迭代次数增大到一定程度时，误差收敛到一个稳定值，式（14）的解趋近于式（12）的解. 又由于式（14）可以利用优化迭代的方法对各变量进行分离求解，在计算复杂度方面远低于式（12），因此下面根据式（14）来求解清晰图像. 式（14）的求解步骤如下所示，其中上标t为迭代次数.

1) 初始化： $\,\beta _1^{\left( 0 \right)} = 1,\;\beta _2^{\left( 0 \right)} = 1,\;{{{f}}^{\left( 0 \right)}} = {{g}}$；

2) 循环；

$\begin{aligned}3) \; w_k^{\left( {t + 1} \right)} = \mathop {\min }\limits_{{w_k}} \;\sum\limits_{i = 1}^N {\sum\limits_{k = 1}^K {{\alpha _k}\ln\; \left( {1 + \frac{1}{2}{{\left( {w_{ki}^{\left( t \right)}} \right)}^2}} \right)} } + \\ \frac{{{\beta _2}}}{2}\sum\limits_{k = 1}^K {{\alpha _k}{{\left( {w_k^{\left( t \right)} - {{J}}_k^{\rm{T}}{{{f}}^{\left( t \right)}}} \right)}^2}} {\text{；}} \quad\quad\end{aligned}$

$\begin{aligned}4) \;{{{u}}^{\left( {t + 1} \right)}} = \mathop {\min }\limits_{{u}}\; \frac{\lambda }{2}\sum\limits_{i = 1}^N {\left( {{{{u}}_i}^{\left( t \right)} - {{{g}}_i}\ln\; {{{u}}_i}^{\left( t \right)}} \right)} +\\ \frac{{{\beta _1}}}{2}{\left( {{{{u}}^{\left( t \right)}} - {{{f}}^{\left( t \right)}}} \right)^2}{\text{；}}\quad\;\;\quad\quad\quad\end{aligned}$

$\begin{aligned}5)\; {{{f}}^{\left( {t + 1} \right)}} = \mathop {\min }\limits_{{f}}\; \frac{{{\beta _1}}}{2}{\left( {{{{u}}^{\left( {t + 1} \right)}} - {{{f}}^{\left( t \right)}}} \right)^2}+\quad\quad\quad\quad \end{aligned}$ 　　　　　　 $\begin{aligned}\quad\frac{{{\beta _2}}}{2}\sum\limits_{k = 1}^K {{\alpha _k}{{\left( {w_k^{\left( {t + 1} \right)} - {{J}}_k^{\rm{T}}{{{f}}^{\left( t \right)}}} \right)}^2}} + {l_ + }\left( {{f}} \right){\text{；}}\end{aligned} $

6) $\,{\beta _1} = 2{\beta _1}$，直到 $\,{\beta _1} = {\beta _{1{\rm{max}}}}$；

7) $\,{\beta _2} = 2{\beta _2}$，直到 $\,{\beta _2} = {\beta _{2{\rm{max}}}}$；

8) 输出： $\widehat {{f}} = {{{f}}^{\left( {t + 1} \right)}}$， $\widehat { f}$是清晰图像结果.

以下依次对求解步骤中涉及的几点问题进行解决，为了便于描述，略去上标 $t$.

由于步骤3)中的问题不存在解析解，采用Newton-Raphson方法^[18]进行求解. 以 $k = 1$为例，假设：

(15) $ \begin{split} {w_1} = \mathop {\min }\limits_{{w_1}}\; H\left( {{w_1}} \right) = & \mathop {\min }\limits_{{w_1}} \;{\alpha _1}\ln \;\left( {1 + \frac{1}{2}{w_1}^2} \right) +\\& \frac{{{\beta _2}}}{2}{\alpha _1}{\left( {{w_1} - {{J}}_1^{\rm{T}}{{f}}} \right)^2}. \end{split} $

根据Newton-Raphson方法可得

(16) ${w_1}^{\left( {n + 1} \right)} = {w_1}^{\left( n \right)} - \frac{{{H'}\left( {{w_1}} \right)}}{{{H{''}}\left( {{w_1}} \right)}}.$

对于每一个 $k$，根据以上方法即可求解得到 ${w_k}$.

针对步骤4)中的问题，设定损失函数 $L\left( {{x}} \right)$：

(17) $L\left( {{x}} \right) = \frac{\lambda }{2}{{x}} - \frac{\lambda }{2}{{g}}\ln \;{{x}} + \frac{{{\beta _1}}}{2}{\left( {{{x}} - {{f}}} \right)^2}.$

则 $L\left( {{x}} \right)$一阶导数为

(18) $\frac{{{\rm{d}}\left( {L\left( {{x}} \right)} \right)}}{{{\rm{d}}{ x}}} = \frac{\lambda }{2} - \frac{{\lambda {{g}}}}{{2{{x}}}} + {\beta _1}\left( {{{x}} - {{f}}} \right).$

当 ${{{\rm{d}}\left( {L\left( {{x}} \right)} \right)}/ {{\rm{d}}{{x}}}} = 0$时，有

(19) $\frac{\lambda }{2}{{x}} - \frac{{\lambda {{g}}}}{2} - {\beta _1}{{fx}} + {\beta _1}{{{x}}^2} = 0.$

由于 $ - {{\lambda {{g}}} / {\rm{2}}} \leqslant 0$，且 ${{{{\rm{d}}^2}\left( {L\left( {{x}} \right)} \right)} / {{\rm{d}}{{{x}}^2}}} > 0$，式（19）的非负根即为 ${{x}} = \mathop {\min }\limits_{{x}}\; L\left( {{x}} \right)$的最优解，即

(20) $ {{u}} = \frac{{\rm{1}}}{{2{\beta _1}}}\left\{ {{\beta _1}{{f}} - \frac{\lambda }{2} + {{\left[ {{{\left( {{\beta _1}{{f}} - \frac{\lambda }{2}} \right)}^2} + 2\lambda {\beta _1}{{g}}} \right]}^{{\rm{1/2}}}}} \right\}. $

针对步骤5)中的问题，可在频域中进行求解. 首先，在不考虑非负约束 ${l_ + }\left( {{f}} \right)$时，其解为

(21) ${{f}} = {F^{ - 1}}\left( {\frac{{{\beta _1}{{U}} + {\beta _2}\displaystyle\sum\limits_{k = 1}^K {{\alpha _k}{{{\varGamma}} _k}^*{W_k}} }}{{{\beta _1} + {\beta _2}\displaystyle\sum\limits_{k = 1}^K {{\alpha _k}{{{\varGamma}} _k}^*{{{\varGamma}} _k}} }}} \right).$

式中： ${{U}}{\text{、}}{{{\varGamma}} _k}{\text{、}}{W_k}$分别为 ${{u}}{\text{、}}{{{J}}_k}{\text{、}}{w_k}$的离散傅里叶变换，*为共轭运算， ${F^{ - 1}}$为离散傅里叶逆变换.

由于 ${\beta _1} > 0$， ${\beta _2} > 0$， ${\alpha _k} > 0$， ${{{\varGamma}} _k}^*{{{\varGamma}} _k} \geqslant 0$，式（21）括号中的分母不存在零点. 当不存在非负约束时，步骤5)中问题的代价函数是凹函数，并且由式（21）所得的结果 ${{f}}$是其唯一解， ${{f}}$的任意变化都将导致代价函数值的增加. 因此，当步骤5)中的问题带有非负约束时，只需要对式（21）的结果增加逐点非负操作即可得到最优解，如下式所示：

(22) $\widehat {{f}} = \max\; ({{f}},{ 0}).$

式中： $\max\;({{f}},{\bf 0})$表示将 ${{f}}$逐点与0比较大小，并取较大值.

由此，求解过程中的所有问题得到解决，按照求解步骤即可解得清晰图像 ${{f}}$.

为了验证提出的基于FoE模型的泊松噪声去除算法的有效性，针对泊松噪声出现环境，选择大小为 $512 \times 512$的星云图像，使用8种方法作为对比算法. 这8种算法可以分为2类：1）面向高斯噪声的去噪方法，包括高斯滤波（Gaussian filter，GF）、各项同性总变分^[19]（Gaussian total variation L2，GTVL2）、自然图像梯度稀疏约束^[20]（Gaussian gradient sparse，GGS）；2）面向泊松噪声的去噪方法，包括最大似然估计（Poisson Richardson-Lucy，PRL）、各项同性总变分（Poisson total variation L2，PTVL2）、自然图像梯度稀疏约束（Poisson gradient sparse，PGS）、加权核范数最小化^[21]（weighted nuclear norm minimization，WNNM）以及文献[22]中提到的（fields of experts by using Anscombe transform，ATFoE）方法. 由于图像不含泊松噪声，在实验中通过仿真方法添加泊松噪声. 实验中采用的FoE模型的邻域系统包含8个滤波器模板，窗口大小为 $3 \times 3$，如图1所示.

以星云测试图像为例，采用上述不同方法对图像进行去噪的结果分别如图3所示，图中左上角为每幅图像的峰值信噪比（peak signal-noise-ratio，PSNR）. 从图3中可以看出，采用GF方法处理后图像整体模糊；采用GTVL2与GGS方法处理后的图像在亮度较高的区域（如：星云中心处）仍含有部分噪声，噪声去除不彻底；PRL方法对噪声的去除效果较差，结果图像中仍含有大量噪声；采用PTVL2、WNNM以及ATFoE方法处理后图像边缘模糊，细节信息丢失；采用PGS方法处理后在高亮度区域（如：星云中心点处）产生不规则花纹，会对接下来的处理产生影响. 采用提出的方法处理后的图像边缘清晰，纹理信息完整，图像中噪声基本得到去除，如图3（k）所示.

为进一步对所提方法进行客观的评价，选用峰值信噪比以及结构相似度（structural similarity，SSIM）作为评价标准. 使用8幅标准测试图像进行对比实验，其中大小为 $1\;024 \times 1\;024$的图像4幅（“baboon”、“boat”、“cameraman”、“man”），大小为 $512 \times 512$的图像4幅（“couple”、“peppers”、“lake”、“lena”）. 8幅测试图像在9种去噪算法下的PSNR和SSIM分别如图4和图5所示.

由图4和图5可以看出，除少数个别图像以外，使用所提出的方法得到的去噪后图像的PSNR和SSIM值显著优于其他方法. 相较于PGS方法和WNNM方法，所提出的算法在PSNR上平均提升了0.48和0.18 dB. 与文献[22]中同样使用FoE模型的ATFoE算法相比，在PSNR结果上，本文方法优势不明显，但是在SSIM结果上，本文方法具有明显的优越性. 与传统的GF、PRL及TVL2方法相比，本文算法能够更加有效地消除噪声，获得的图像失真度小，边缘纹理结构保留较好，与原始图像相似度高. 以上对比结果充分证明了本文算法的有效性.

在峰值信噪比方面，本文提出的方法与传统方法相比均有提升，至少提升了0.18 dB. 在结构相似度方面，由于现有的去噪算法大多是直接或间接地利用整数阶积分函数改造的滤波器进行滤波，这样的滤波器会破坏图像中的高频信息致使图像边缘或纹理丢失. 提出的算法所采用的FoE模型主要是利用图像中相邻像素之间的关系，因此边缘纹理得到了保留. 本文在贝叶斯最大后验概率估计的框架下，构建了图像泊松噪声去噪模型. 该模型以泊松概率分布为保真项，以马尔科夫专家场为正则项构造代价函数，利用二次惩罚函数优化求解，能够有效地还原清晰图像，有着不错的泊松噪声去除效果. 但由于求解过程需要多次迭代计算，运算复杂度较高，耗时较长，未来需要考虑精简计算以提高去噪速度.