1.3　图像复原非线性迭代算法

基于逆滤波和Tikhonov正则化的图像复原方法是线性的，存在封闭解，然而，这两种方法均存在缺陷，逆滤波的解不稳定，而Tikhonov正则化方法的解又过于平滑。基于全变差或是高阶变分、小波框架理论、稀疏理论和随机场的非线性正则化复原方法更易得到良好的复原结果，然而这些方法往往并不存在封闭解，其求解需要借助数值迭代算法。事实上，迭代求解方式更利于将关于解的先验知识融入求解过程，也更有利于实现对复原过程的“监控”［114］。

尽管所采用的正则化方法不同，但当前图像复原的基本实现途径都是建立包含正则项和保真项的最小化函数，然后求取正则化函数的极小点，并以此作为图像复原的结果。目标函数的凸性对于求解过程的快速实现和解的稳定性都极为重要。若目标函数为非凸，则结果很难保证为目标函数的全局最优解。在以稀疏性为基础的正则化方法中，通常会将NP难的l0优化问题凸松弛为l1的凸优化问题，能够证明，在十分宽松的条件下，l1凸优化问题的解趋向于相应的l0非凸优化问题的解［95］。因此，基于l1范数的正则化在当前的图像反问题中应用最为广泛。

1.3.1　传统方法

TV模型是最具代表性的l1正则化子（regularizer），因此，以TV模型的求解为例说明图像反问题中非线性函数求解算法的发展历程。事实上，早期的一些方法通常是针对具体正则化模型专门设计的。

当噪声为Gauss白噪声时，基于TV的图像复原有以下目标函数：

　　 （1-16） 　　

其中为一阶差分算子，是l1型的凸函数，因而上式为典型l1-l2最小化问题。由于TV范数的不可微性和非线性式（1-16）的求解并非易事。尽管TV模型被引入图像处理已有超过二十年的历史，迄今，最小化函数式（1-16）仍然是检验众多新算法的试金石。

最早用于求解TV去噪（K=I）的方法是Rudin等人提出的时间推进（time-marching）算法［6］，该方法将时间变量引入函数式（1-16）的Euler-Lagrange方程，不仅速度缓慢，计算精度也不尽人意。随后，出现了求解TV去噪模型的滞后扩散不动点法［115］，一定程度上克服了时间推进方法的不足。该方法在处理TV的不可微点时需在方程分母上引入一个小的常数，这一策略在当前的一些文献之中依然可见［28］。2004年Chambolle［116］基于TV模型的对偶模型提出了经典的对偶梯度下降法，并严格证明了算法的收敛性以及所需的收敛性条件。该方法大体思路是，先将模型式（1-16）转化为如下原始-对偶（primal-dual）形式：

　　 （1-17） 　　

其中K=I，div为散度算子，其Hilbert伴随算子为。设置“min”目标函数梯度为零可得u=f-λ-1divv，将其代入式（1-17）得到模型式（1-16）的对偶模型：

　　 （1-18） 　　

其最优化必要条件为：

　　 （1-19） 　　

随后Chambolle采用了如下隐式人工推进方法迭代求解对偶变量：

　　 （1-20） 　　

最终根据上述的原始变量和对偶变量的关系解出原始变量，即为最终的去噪图像。该方法通过求解TV去噪模型式（1-16）的Fenchel对偶［5］模型巧妙规避了TV模型的不可微问题，成为迄今效率最高的图像去噪方法之一。此外，它也常作为嵌套算法出现在一些图像反卷积算法之中［15］。但是这种对偶思想很难直接推广到其他图像反问题中去，原因在于当退化矩阵K不是单位阵时，式（1-16）的对偶模型中会出现K-1。不幸的是K可能是奇异的，即其逆矩阵可能并不存在。

此外，用于非线性正则化模型求解的方法还有二阶锥法［117］、正交投影法［118］、内点法［119］和预优法［120］等方法。这些传统方法虽能够针对某一特定问题给出合理的解，但其也存在近似求解、无法充分挖掘问题本身的结构、不利于大规模并行计算等故有缺陷，这些都限制了它们在当前图像大数据背景下的应用。

1.3.2　算子分裂方法

近些年来，为更好地应对高维、海量、高品质要求的图像大数据处理问题，一类强大的、通用的、灵活的、并行的算子分裂方法被引入到了图像反问题领域，其基本思想是“化繁为简，分而治之”。它们共同的数学基础是由Fenchel（1905—1988）、Moreau（1923—2014）和Rockafellar（1935—）等先驱者所奠定的现代凸优化分析。次微分（subdifferential）、临近映射（proximal mapping）、卷积下确界（infimal convolution）等概念被频繁地运用。这类方法可以更好地应对目标函数的非线性和非光滑性以及各种繁杂的约束条件。通常算子分裂方法只需要挖掘目标函数的一阶信息，因而其计算实现又是足够简单的。

很多信号处理问题，如图像复原问题，通常可建模为以下最小化模型：

　　 （1-21） 　　

其中f1，…，fm为从映射到（-∞，+∞］的凸函数。求解这一模型通常会碰到的一个难题是某些函数项是不可微的，这就使得一些传统的光滑优化技术无用武之地。而算子分裂方法则通过“分离”并分别求解这些函数项，得出可行的求解算法。通常的一个假定是非光滑函数fi是“可临近的”（proximal），即其临近算子（proximity operator）存在封闭解或可以被方便地求得，临近分裂（proximal splitting）是算子分裂方法的基础。实际应用表明，这一假设是足够宽松的。尽管临近算法被引入图像处理领域的时间并不长，但它的推广散布却异常迅速。

在介绍常用的几个算子分裂方法之前，首先给出要用到的几个凸分析概念。

记为N维的Euclidean空间，记<·，·>为内积符号，记凸函数的域为；记正常凸函数集合映射到（-∞，+∞］的域是非空的，且下半连续［5］的凸函数。f的fenchel共轭（fenchel conjugate）定义为：

　　 （1-22） 　　

f的次微分为点集（set-valued）映射：

　　 （1-23） 　　

通常所讲的次微分是集合，而次梯度（subgradient）指的是其中的某一个元素，显然，次梯度概念是对光滑函数梯度概念的推广。应用次梯度可以得到非光滑凸函数最小化的Fermat法则［5］（Fermat’s rule）。

Fermat法则：若，则有：

　　 （1-24） 　　

次梯度的一个重要性质是极大单调性［5］（关于极大单调算子和极大单调性质的介绍见2.4.3节），即它满足：

　　 （1-25） 　　

且算子I+∂f的值域为。该性质在有关分裂算法的收敛性证明中常常用到。

f关于点x与x'的Bregman距离定义为：

　　 （1-26） 　　

Bregman距离是广义距离，显然它不满足对称性，其非负性可由次微分的定义得到。

设Ω为中的非空集合，其示性函数（indicator fuction）定义为：

　　 （1-27） 　　

示性函数的fenchel共轭为支撑函数（support fuction），其定义为：

　　 （1-28） 　　

容易验证如果Ω为非空凸集，则有。通过引入示性函数，可以将一个约束优化问题转化为等价的无约束优化问题，从而更方便于算子分裂类方法的应用。

记的临近算子为：

　　 （1-29） 　　

临近算子是次微分算子的预解算子（定义见第2章），即，prox是单值映射的，且是固定非扩张的（firmly nonexpansive）（定义见第2章），即：

　　 （1-30） 　　

此外，次微分的反射算子（reflection operator，或是反射预解算子，reflected resolvent）2proxf-I也是非扩张的（定义见第2章）［5］。运用临近算子proxf的固定非扩张性通常可以将一个复杂的函数极小化问题转化为不动点问题，而其单值性则对相应算法的稳定性有着至关重要的意义。容易验证凸集Ω示性函数ιΩ的临近算子为投影到Ω的投影算子，因此，临近算子被认为是投影算子的推广［121］。事实上若proxf存在，则f（x）的极小点可通过如下临近不动点迭代求得：

　　 （1-31） 　　

上式又可写为或。这就是所谓的临近点算法（proximal point algorithm，PPA）。接下来，简要介绍当前流行的几种以临近点方法为基础的算子分裂方法。

（1）Bregman迭代与线性Bregman迭代方法

Osher等于2005年将Bregman迭代方法引入了图像处理领域，并将其应用到了基于TV的图像去噪与去模糊［122］。相比于此前的方法，该方法有着较好的通用性和较高的计算效率。Bregman迭代法可以用来求解以下类型的图像反问题：

　　 （1-32） 　　

式中，算子A可以是非线性的，Bregman迭代方法的迭代规则为：

　　 （1-33） 　　

其中pk为f在xk处的次梯度，β∈（0，+∞）为惩罚参数，为Bregman距离。如果ϕ=A是线性的（这种情况在图像反问题中更为普遍），迭代规则式（1-33）可以转化为下列更为紧凑的形式［123］：

　　 （1-34） 　　

针对同样类型的问题，Yin等则将Bregman迭代法加以改进，得到了线性Bregman迭代法［123，124］，并将其应用到了压缩感知的基追踪问题上（）。令式（1-33）中ϕ=A为线性，其基本思想是将二次项在xk附近做Taylor展开：

　　 （1-35） 　　

线性Bregman迭代法针对问题式（1-32）（ϕ=A）的迭代规则为：

　　 （1-36） 　　

需要指出的是，只有当δ满足上述给定条件时，线性Bregman迭代法才是收敛的，该方法通常具有比基本Bregman迭代更高的执行效率，因为它往往可以避免反问题中常见的矩阵求逆环节。

问题式（1-32）显然是针对“干净”数据的，当观测数据含有噪声时，两种Bregman方法则采用作为停机准则，这需要根据噪声水平预先估计出c的值。文献［122］和［124］分别提供了Bregman迭代法和线性Bregman法的收敛性证明。Bregman方法在处理更复杂的问题时显得捉襟见肘，但它们为后续分裂Bregman方法的提出建立坚实的理论基础。

（2）分裂Bregman法

2009年，Goldstein和Osher［125］在变量分裂（variable splitting，VS）［126］和Bregman迭代法的基础上提出了分裂Bregman算法（splitting bregman algorithm，SBA），并将其应用到了l1正则化的反问题中。该方法可以用来求解以下类型的问题：

　　 （1-37） 　　

首先通过引入辅助变量d，可将上式转化为下列等价线性约束优化问题：

　　 （1-38） 　　

令凸函数E（x，d）=f1（x）+f2（d），参考Bregman迭代的思想可以得到以下迭代形式：

　　 （1-39） 　　

这即是所谓的分裂Bregman方法。第一步中的xk+1和dk+1的求解可以交替进行。看上去第一步需要引入一个嵌套迭代，但实际上它并不需要精确求解，且精确求解的精度会被后续变量的更新浪费掉。实验表明，xk+1和dk+1只需交替求解一次即可，且该情况下的收敛性依然可以严格证明［127，128］。与Bregman迭代不同的是，这里等式约束的两端都是变量。当ϕ=A为线性算子时，迭代规则式（1-39）同样可以简化为：

　　 （1-40） 　　

上式中，xk+1和dk+1的求解仅交替进行了一次。辅助变量的引入是十分重要的，它使得两个凸函数项在求解时可以分离开来，而关于辅助变量d的求解是一个临近点问题。

运用线性Bregman的思想也可以方便地将分裂Bregman方法转换为线性分裂Bregman方法［129］。

（3）交替方向乘子法与线性交替方向乘子法

与Bregman分裂法类似，交替方向乘子法（alternating direction method of multipliers，ADMM）［15，130］，又称为交替方向法（alternating direction method，ADM），同样融入了变量分裂思想。它通过求解式（1-38）（ϕ=A）的增广Lagrange函数的鞍点来求解式（1-37）。式（1-38）的增广Lagrange函数为：

　　 （1-41） 　　

其中b为Lagrange对偶变量，又称为Lagrange乘子，β∈（0，+∞）为惩罚参数。增广Lagrange方法（augmented lagrangian method，ALM）［15］可以通过以下迭代规则求解式（1-41）的鞍点：

　　 （1-42） 　　

与分裂Bregman迭代法类似，第一步并不需要精确求解，若xk+1和dk+1交替求解的迭代次数为1，则可得到如下交替方向乘子法的迭代规则：

　　 （1-43） 　　

式（1-40）与式（1-43）仅相差一个常数，而这一常数通过变量替换可以消除，这说明分裂Bregman与交替方向法在线性约束条件下是完全等价的。

如果对式（1-43）第一步中的二次项在xk附近做Taylor展开则可以导出适用范围更广的线性交替方向乘子法（linearized ADMM，LADMM）［92，131-135］：

　　 （1-44） 　　

有关ADMM类算法的发展历程，可以参阅综述文献［136］～［138］。

（4）前向-后向分裂法

前向-后向分裂（forward-backward splitting，FBS）法［139］用以解决以下类型问题：

　　 （1-45） 　　

其中f1是“可临近的”，即其临近算子存在闭合形式或可以方便地求解，f2则是可微的，且其梯度为Lipschitz连续，即：

　　 （1-46） 　　

取ε∈（0，min{1，1/γ}），则式（1-45）可通过如下迭代规则求解：

　　 （1-47） 　　

式中，涉及梯度下降的第一步被称为前向步骤，而涉及临近算子的第二步被称为后向步骤。式（1-47）的导出运用了临近算子的非扩张性［5］，具体的是：

　　 （1-48） 　　

其中Fix表示不动点。事实上，FBS算法可以看作PPA的一种推广，FBS的（次）梯度形式为：

　　 （1-49） 　　

即为分离f1与f2将（PPA形式）换成了。

当λk≡1且f1=ιΩ时，记PΩ为投影到凸集Ω的投影算子，式（1-49）转化为：

　　 （1-50） 　　

即经典的梯度投影（gradient projection）算法。

运用前向-后向分裂法的一个经典例子是Beck和Teboulle的快速迭代收缩/阈值算法（fast iterative shrinkage/thresholding algorithm，FISTA）［140］，其迭代规则为：

　　 （1-51） 　　

在文献［140］中，Beck和Teboulle将FISTA应用到了TV去噪问题式（1-16）的对偶问题上，但是在将FISTA应用到TV去模糊时，则需要嵌套使用该算法，原因是TV去模糊的对偶问题涉及模糊矩阵求逆，而这一过程通常是无法实现的。

（5）Douglas-Rachford分裂法与Peaceman-Rachford分裂法

上述的FBS算法需要式（1-45）两函数中的一项为可微，这一点对于很多实际应用是较为苛刻的，Douglas-Rachford分裂（Douglas-Rachford splitting，DRS）法［141］则只要求式（1-45）中两函数的临近算子存在，它通过以下迭代实现两函数项的解耦：

　　 （1-52） 　　

记􀳱，则DRS的次梯度形式为：

　　 （1-53） 　　

即：

　　 （1-54） 　　

根据临近算子prox的非扩张性和反射预解算子的非扩张性（见2.4.3节），TPRS是非扩张的，故（I+TPRS）/2为1/2平均算子，是固定非扩张的（见2.4.2节）。将式（1-54）中的1/2平均进行松弛，则可以得到下列Peaceman-Rachford分裂（Peaceman-Rachford splitting，PRS）算法［142，143］：

　　 （1-55） 　　

其次梯度形式为：

　　 （1-56） 　　

（6）原始-对偶分裂法

原始-对偶分裂（primal-dual splitting，PDS）法［144］可以通过求解Lagrange问题的鞍点，同时求解原始问题和对偶问题。相比于前述的分裂方法，原始-对偶分裂灵活性更强，所能解决的问题也更宽泛，因而逐步成为新的研究热点。

Chambolle和Pock两人于2011年提出了一种用于求解

　　 （1-57） 　　

的一般性原始-对偶算法［145］，掀起了原始-对偶理论在图像反问题中的应用热潮。

利用f2（Ax）的Fenchel共轭可以得到问题式（1-57）的Lagrange函数：

　　 （1-58） 　　

再次运用f1（x）的Fenchel共轭可以得到问题式（1-57）的Fenchel-Rockafellar对偶问题：

　　 （1-59） 　　

其中，A*为A的伴随算子（adjoint operator，又称共轭算子，Hilbert空间中有<Ax，y>=<x，A*y>，Euclidean空间中，A*即为AT）。

Chambolle和Pock的方法可以求解式（1-58）的鞍点，在不存在对偶间隙的情况下，式（1-57）和式（1-59）的目标函数值相同。该方法的迭代规则为：

　　 （1-60） 　　

改进格式为：

　　 （1-61） 　　

此外，Chambolle和Pock还针对不同情况对式（1-61）作了改进。

事实上，早在2008年，Zhu和Chan已将式（1-61）在θ=0时的特殊情况［146］（笔者将其命名为原始-对偶联合梯度法，primal-dual hybrid gradient，PDHG）应用于基于TV的图像复原模型式（1-16），并以此来解决TV原始模型不可微和对偶问题（TV去噪）解不唯一的情况，但当时笔者并未将该算法一般化，也未严格地分析其收敛性。

当前，原始-对偶分裂方法多基于极大单调算子理论和非扩张算子理论来构建，而并非简单地应用次梯度的一些性质，这样可以导出一些性能更为强大、适用范围更广的算法［147-151］。

事实上，以上所讨论的算子分裂方法在某些特定的情况下可以建立等价关系。如分裂Bregman算法（SBA）、交替方向乘子法（ADMM）和Douglas-Rachford分裂算法（DRSA）在应用于原始问题式（1-57）（P）和对偶问题式（1-59）（D）时具有如图1-2所示的等价性关系［152］。

以上算子分裂方法也存在一些共同缺陷，如多数基本算法仅针对目标函数包含两个函数项的情况；很多分裂方法在设计过程中会引入辅助变量，这使得这些方法在求解上相对烦琐复杂［153］。发展并行的算子分裂方法［154-158］，并将其与分布式计算相结合来应对图像大数据问题，是今后很长一段时期内学术界研究的热点问题。此外，原始-对偶算法与其他方法之间的联系也是值得深入研究的学术问题［153］。

1.3.3　分裂算法的收敛性分析

目前，有关算子分裂方法的一个重要研究热点和难点是其渐进收敛行为分析。算子分裂方法之所以发展迅速，应用广泛，一个很重要的原因是其理论基础坚实，这既体现在算法的设计上，也体现在收敛性的分析上。通常，算法的收敛性分析包含两个方面，一个是收敛性证明，即算法是否能够准确找到目标函数的解；另一个则是收敛速率分析，即算法能以多快速率逼近问题的最优解。收敛速率可以通过不动点残差（fixd-point resudual，FPR，即连续两步迭代结果的Euclidean距离）、目标函数值偏差（objecitve error）和对偶间隙（duality gap）来描述［143，159］。

目前算子分裂的收敛性分析有两种主要途径，变分不等式（variation inequality，VI）和非扩张算子［5］。较先推广应用的Bregman类算法和ADMM算法大都基于变分不等式进行算法收缩性和收敛性分析，而前向-后向分裂法、Douglas-Rachford分裂法和最新的原始-对偶分裂法因其导出与极大单调算子和非扩张算子直接相关，大都会通过非扩张算子不动点迭代的收缩性进行收敛性分析。相比于基于非扩张算子收缩性的方法，基于变分不等式的方法更为繁杂，但对理论基础要求较低。我国南京大学知名学者何炳生教授在其课件《凸优化和单调变分不等式的收缩算法》（个人网站）中详细介绍了有关变分不等式的基础。

Yin Wotao等［160-163］（美国UCLA）以及何炳生团队［127，128］对Bregman类算法和ADMM类算法进行了系统的收敛性分析。Yin等在文献［143］和［159］中指出，对于现有的诸多算子分裂方法而言，在不加附加条件的情况下通常可以证得收敛速率O（1/k）；在文献［161］中则指出，如果式（1-57）中的某项具备强凸性和Lipschitz连续的梯度，则ADMM算法可以获得收敛速率O（1/ck）（c为某个大于0的常数）。何炳生等在文献［127］和［128］中分别证明了DRS/ADMM遍历和非遍历的O（1/k）收敛速率。此外，Goldstein等在文献［164］中通过采用Nesterov加速法［165］（同样被FISTA［140］所采用）使ADMM获得了O（1/k2）的收敛速率；Chambolle和Pock则证明了其原始-对偶分裂算法［145］的O（1/k）收敛性，在采用Nesterov加速法时证明了其O（1/k2）收敛性，并在假设两函数项均为强凸的情况下证明了算法的O（1/ck）收敛性。

1.3.4　正则化参数的自适应估计

尽管图像复原技术已走过50多年的发展历程，但到目前为止，实时图像复原技术仍鲜见报道，其原因在于，一是图像复原计算量大，对硬件技术要求高；二是算法的自动化实现问题仍未得到很好的解决。图像复原算法自动化实现的一个关键问题是正则化参数的自适应选取。

正则化参数起到平衡保真项与正则项的作用。在正则化函数式（1-7）中，若正则化参数选取过小，则容易致使复原结果过分偏离观测数据，导致过平滑的结果；若正则化参数选取过大，则又会导致复原结果中含有噪声，不够平滑。

选取正则化参数最简单的方法是在求解目标函数前人为选定，这也是当前大多数文献中的做法［84，125，126，166，167］。然而，这种人为选取的方式不仅耗时过长，且不利于图像复原的全自动化实现。此外，诸多因素如图像噪声水平、模糊函数类型和尺寸、图像类型等，都会对正则化参数的选取产生影响［168］，这对于执行者的经验是一个不小的挑战。

当前，带有正则化参数自动更新的自适应图像复原算法越来越受到学术界的重视，成为图像处理的一个热点问题。从现有的研究成果来看，多种手段可以用于实现正则化参数的自适应选取，这些方法包括Morozov偏差原理［15，30，168-172］（需要关于噪声水平的先验知识）、广义交叉确认法［173，174］（generalized cross-validation，GCV）（无需输入噪声水平）、L曲线法［175，176］、无偏预先风险估计法［177］（unbiased predictive risk estimator，UPRE）、变分Bayes方法［106］和参数下降法［178］等。

GCV方法在正则项具有二次形式时可以方便地应用，然而，这在实际中并不容易满足（如TV正则化的图像复原），GCV公式极小点的求得也并不容易，此外，该方法易导致欠平滑的结果［30］。L曲线法通过找到关于正则项和保真项对数曲线的角点来确认正则化参数，但若曲线本身较为光滑则角点会难以确定，且该方法计算量庞大。变分Bayes方法的正则化参数可以表示为图像和噪声模型超参数的函数，在参数估计的过程中，可以自然地得到正则化参数，但变分Bayes方法自身的求解就是一个难题。参数下降法首先选择一个较大的正则化参数对目标函数进行求解，再采用某种策略减小正则化参数，直至满足一定的停止准则，则最终的参数即为所要求取的正则化参数。然而，参数下降准则和算法停止准则的选取是该方法所要面对的一个难题，这也使得算法的收敛性难以保证。

当观测图像的噪声水平可估计时，Morozov偏差原理是实现正则化参数自适应选取的可行方案。该原理通过匹配残差到某一上限来确定正则化参数。根据该原理，复原结果的可行域为：

　　 （1-62） 　　

其中c为噪声相关的常数。若观测噪声为Gauss噪声，则c=τmnσ2，其中τ为预先确定的噪声依赖的常数（当噪声不为Gauss白噪声时，不等式具有不同形式）。该方法的本质是直接求解约束优化问题式（1-6），并进行正则化参数估计。

当前，基于偏差原理的自适应图像复原所存在的主要问题是，在实现正则化参数自适应选择的同时，需要在基本迭代算法中引入内部迭代［15，30，170-172］。此外，自适应图像复原算法的收敛性并不能直接由非自适应算法的收敛性保证，需要有严格的理论证明，而这在大多数的文献中并未涉及。

在过去的十年间，算子分裂研究无论在理论方法，还是在实际应用上都取得了显著进步，特别是近几年，国内外学者在基于算子分裂的图像反问题求解上取得了丰硕成果。当前，面对高维、海量的图像大数据处理问题，发展一类通用灵活、并行快速的算子分裂方法成为必然。