图像恢复问题是一类典型的不适定问题, 它可以不满足一个或多个Hadamard条件, 原因在于成像系统本身的特性或观测数据的误差, 所以图像恢复的过程无论是理论分析还是数值模拟都存在着一定的困难.目前众多学者对图像恢复问题进行了深入的研究, 主要的恢复算法分为直接法和迭代法两大类.由于直接法存储需求和计算量都比较大, 且难以控制, 所以迭代法就成为首选.其中与正则化技术相结合的Krylov子空间方法备受广大学者青睐, 当这种Krylov子空间方法的天然正则性并不足以获取足够的信息时, 通过添加额外的正则化来获得足够好的近似解, 这就构成了混合正则化方法的基本思想. 1981年, O'Leary和Simmons^[1]最先提出了混合方法, 先利用Lanczos双对角化过程将问题投影成低维问题, 再对小问题进行TSVD正则化. 2001年, Hanke^[2]也讨论了基于Lanczos方法的混合正则化问题. 2011年, 赵苗苗^[3]等提出了图像恢复的正则化GMRES方法. 2014年, Novati等^[4]提出用于图像恢复的自适应Arnoldi-Tikhonov正则化方法.

考虑到Tikhonov正则化方法^[5]在求解不适定问题中具有明显的优势, 并充分利用LSQR算法在数值计算上收敛快、稳定的优点, 本文将其与正则化技术相结合应用于图像恢复问题中.

图像在形成、传输和记录过程中, 由于很多原因使图像的质量变坏, 称为图像退化^[6].第一类Fredholm积分方程常用于描述图像的线性退化模型, 表达式为

$ \begin{equation} \int_a^b\int_c^{d}K(u, v, x, y)f(x, y)dxdy+\eta (u, v) = g(u, v), \end{equation} $

(2.1)

其中$ g(u, v) $是观察到的图像, $ f(x, y) $为原始图像, $ \eta(u, v) $为加性随机噪声, 积分方程的核函数$ K(u, v, x, y) $称为点扩散函数(PSF).点扩散函数有多种类型, 本文中所使用的核函数为高斯核, 形如

$ \begin{equation*} K(u, v, x, y) = \frac{1}{2\pi\sigma^{2}}\mathrm{e}^{-\frac{(u-x)^2+(v-y)^2}{2\sigma^{2}}} = \eta_{\sigma}(u-x)\eta_{\sigma}(v-y), \end{equation*} $

其中$ \eta_{\sigma}(t) = \frac{1}{\sqrt{2\pi}\sigma}\mathrm{e}^{-\frac{t^2}{2\sigma^2}} $, $ \sigma = 0.02 $.由于很难求出(2.1)式的解析解, 本文考虑求其数值解, 为此将(2.1)积分方程离散化.所谓积分方程离散化指的就是将积分方程转化为可以数值求解的线性方程组, 从而得到积分方程的数值解.文献[7]中较为完整的叙述了用于积分方程离散的多种求积法则(如矩形公式、梯形公式等等), 在这里, 选用中点公式对(2.1)进行离散,

$ \begin{equation} \int_0^{1}\varphi (t)dt\approx \sum\limits_{j = 1}^{n}w_{j}\varphi (t_j), \end{equation} $

(2.2)

其中$ t_1, t_2, \cdots, t_n $表示对应求积法则的横坐标, $ w_1, w_2, \cdots, w_n $为对应的权重.对于中点公式有

$ \begin{aligned} t_j = \frac{j-0.5}{n}, w_j = \frac{1}{n}, j = 1, 2, \cdots, n. \end{aligned} $

遵循该求积法则, 可以将(2.1)二维第一类Fredholm积分方程离散化.令$ x, y\in[a, b]\times[c, d]$, $ u, v\in[a_1, b_1]\times[c_1, d_1]$, 将积分区间$[a, b]\times[c, d]$分为$ m1 \times n1 $等份, $[a_1, b_1]\times[c_1, d_1]$分为$ m2 \times n2 $等份, 可得到

$ \begin{equation} \begin{split} g(u_p, v_q) & = \int_a^{b} \int_c^{d}K(u_p, v_q, x, y)f(x, y)dxdy \approx \sum\limits_{j = 1}^{n1} \sum\limits_{i = 1}^{m1} K(u_p, v_q, x_i, y_j)f(x_i, y_j)w_{i}w_{j}^{'}, \end{split} \end{equation} $

(2.3)

其中$ i = 1, 2, \cdots, m1 $; $ j = 1, 2, \cdots, n1 $, $ p = 1, 2, \cdots, m2 $; $ q = 1, 2, \cdots, n2 $.当用矩阵表示这些方程时, 得到了线性系统

$ \begin{equation} Af = g, \end{equation} $

(2.4)

其中$ A $是$ (m2 \times n2) \times (m1 \times n1) $矩阵, $ f $和$ g $分别表示为$ m1 \times n1 $维, $ m2 \times n2 $维的列向量.这意味着$ A $的行将取决于参数$ u $和$ v $的离散, $ A $的列将取决于$ x $和$ y $的离散. $ A $, $ f $和$ g $的元素分别为$ A = K(u_p, v_q, x_i, y_j)w_iw_j^{'} $, $ f = f(x_i, y_j) $以及$ g = f(u_p, v_q), $其中$ i = 1, 2, \cdots, m1 $; $ j = 1, 2, \cdots, n1 $, $ p = 1, 2, \cdots, m2 $; $ q = 1, 2, \cdots, n2 $.

由于$ A $是大规模奇异矩阵, 具有严重的病态性, 使得计算大规模矩阵的奇异值分解(SVD)或逆需要很大的运算量和存储量, 利用常见的截断奇异值分解方法(TSVD)、求逆等直接法求解大规模的问题是不实际的.目前求解$ Af = g $的大规模稀疏线性方程组的首选方法是Krylov子空间方法^[8], 其基本思想是在一个具有更小维数的子空间$ K \subset R^n $中寻找满足精度要求的近似解.因此, 这类方法也可以看作是投影方法, 即寻找解析解在子空间中的投影.

求解线性系统(2.4)等价于求解最小二乘问题^[9]

$ \begin{equation} \min\limits_{f\in R^n}\lVert Af-g \rVert, \end{equation} $

(3.1)

核心思想是求出的解向量$ f $能够使矩阵方程两边的误差平方和最小化.事实上, 最小二乘问题(3.1)也等价于下面的正规方程

$ \begin{equation} A^TAf = A^Tg. \end{equation} $

(3.2)

将原方程$ Af = g $转化为正规方程的一个最大缺点是:条件数是原问题的平方.如果原问题的条件数本身就比较大, 则正规方程的条件数就更大, 迭代求解时收敛会更慢, 而且舍入误差的影响也更明显.所以, 如果系数矩阵条件数比较大, 则不建议直接求解正规方程.在这种情况下, 可以采用Paige和Saunders^[10]提出的LSQR方法.

注意到, 原线性方程组等价于下面的鞍点问题

$ \begin{equation} \left[ \begin{array}{cccc} I & A\\ A^T & 0\\ \end{array} \right]\left[ \begin{array}{cccc} r\\ f\\ \end{array} \right]= \left[ \begin{array}{cccc} g\\ 0\\ \end{array} \right]\Leftrightarrow A\widetilde{f} = \widetilde{g}, \end{equation} $

(3.3)

其中$ r = g-Af $.选取初值$[0, 0]^{T} $, 然后采用Lanczos方法构造Krylov子空间$ K(A, g_0) = $ span$ \{g_0, Ag_0, \cdots, A^{m-1}g_0\} $的一组标准正交基, 其中$ g_0 =[g, 0]^{T} $.通过直接计算可知

$ \begin{eqnarray*} v_1 = \frac{1}{\lVert g \rVert_2}\left[ \begin{array}{cccc} g\\ 0\\ \end{array} \right], v_2 = \frac{1}{\lVert A^Tg \rVert_2}\left[ \begin{array}{cccc} 0\\ A^Tg\\ \end{array} \right]. \end{eqnarray*} $

依次类推, 可以发现这样一个规律:基向量$ v_{2k-1} $和$ v_{2k} $分别具有下面的形式

$ \begin{eqnarray*} v_{2k-1} = \left[ \begin{array}{cccc} \ast\\ 0\\ \end{array} \right], v_{2k} = \left[ \begin{array}{cccc} 0\\ \ast\\ \end{array} \right]. \end{eqnarray*} $

因此, 引入另外两组单位向量$ \{u_k\} $和$ \{w_k\} $来表示这组基, 即

$ \begin{eqnarray*} v_{2k-1} = \left[ \begin{array}{cccc} u_k\\ 0\\ \end{array} \right], v_{2k} = \left[ \begin{array}{cccc} 0\\ w_k\\ \end{array} \right], k = 1, 2, \cdots. \end{eqnarray*} $

记$ \beta = \lVert g \rVert_2 $, 则对称鞍点问题(3.3)的Lanczos过程可描述为

$ \begin{eqnarray*} \beta u_1 = g, \alpha_1w_1 = A^Tu_1 , \end{eqnarray*} $

$ \begin{equation} \left\{ \begin{aligned} \beta_ku_{k+1} & = Aw_k-\alpha_ku_k \\ \alpha_{k+1}w_{k+1} & = A^Tu_{k+1}-\beta_kw_k \end{aligned} \right., \quad k = 1, 2, \cdots, \end{equation} $

(3.4)

其中$ \beta_k>0 $和$ \alpha_{k+1}>0 $的选取是确保$ \lVert u_{k+1} \rVert_2 = \lVert w_{k+1} \rVert_2 = 1 $, 即

$ \begin{eqnarray*} \beta_k = \lVert Aw_k-\alpha_ku_k \rVert_2, \alpha_{k+1} = \lVert A^Tu_{k+1}-\beta_kw_k \rVert_2 . \end{eqnarray*} $

上面的这个过程实际上就是一个Lanczos双对角化过程, 是由Golub和Kahan^[11]提出来的.记$ U_m =[u_1, u_2, \cdots, u_m]$, $ U_m =[w_1, w_2, \cdots, w_m]$, 则有

$ \begin{equation} AW_m = U_m\widetilde{H}_m, AU = W_m\widetilde{H}_m^{T}+\alpha_mw_me_{m+1}^{T}, \end{equation} $

(3.5)

其中$ e_{m+1} =[0, 0, \cdots, 1]^T\in R^{m+1} $,

$ \begin{eqnarray*} \widetilde{H}_m = \left[ \begin{array}{cccc} \alpha_1 & \quad & \quad & \quad \\ \beta_1 & \alpha_2 & \quad & \quad\\ \quad & \beta_2 & \ddots & \quad\\ \quad & \quad & \ddots & \alpha_m\\ \quad & \quad & \quad & \beta_m \end{array} \right]\in R^{(m+1)\times m}, \end{eqnarray*} $

且$ U_{m+1}^{T}U_{m+1} = I_{m+1} $, $ W_{m}^{T}W_{m} = I_{m} $.所以矩阵

$ \begin{eqnarray*} \left[ \begin{array}{ccccccc} u_1 & 0 & u_2 & 0 & \cdots & 0 & u_m \\ 0 & w_1 & 0 & w_2 & \cdots & w_{m-1} & 0 \end{array} \right], \end{eqnarray*} $

或者

$ \begin{eqnarray*} \left[ \begin{array}{cccc} U_{m+1} & 0 \\ 0 & W_m \end{array} \right]\end{eqnarray*} $

的列向量组构成了$ K_{2m+1}(A, g_0) $的一组基.

下面, 在子空间$ K_{2m+1}(A, g_0) $中寻找方程组(3.3)的近似解, 记为$[\widetilde{f}, \widetilde{r}]^T $, 使得原方程组的残量范数(即$ \lVert g-A\widetilde{f} \rVert_2 $)达到最小.

首先, 将$[\widetilde{f}, \widetilde{r}]^T $写成

$ \begin{eqnarray*} \left[ \begin{array}{cccc} \widetilde{r}\\ \widetilde{f} \end{array} \right]= \left[ \begin{array}{cccc} U_{m+1} & 0 \\ 0 & W_m \end{array} \right]\left[ \begin{array}{cccc} t_{m+1} \\ y_m \end{array} \right], \end{eqnarray*} $

或者

$ \begin{eqnarray*} \left\{ \begin{aligned} \widetilde{r}& = U_{m+1}t_{m+1}, \\ \widetilde{f}& = W_my_m. \end{aligned}\right. \end{eqnarray*} $

由(3.5)式可知

$ \begin{equation} \begin{aligned} g-A\widetilde{f} & = \beta u_1-AW_my_m \\ & = U_{m+1}(\beta e_1)-U_{m+1}\widetilde{H}_my_m\\ & = U_{m+1}(\beta e_1-\widetilde{H}_my_m). \end{aligned} \end{equation} $

(3.6)

设$ T_m $是Lanczos方法^[12]作用到正规方程$ A^TAf = A^Tg $上得到的对称三对角矩阵, 则可以证明, $ \widetilde{H}_m $就是$ T_m $的双对角化因子, 即$ T_m = \widetilde{H}_m^T\widetilde{H}_m $.由于$ U_{m+1}^{T}U_{m+1} = I_{m+1} $, 因此对(3.6)式取范数有

$ \begin{equation} \lVert g-A\widetilde{f} \rVert_2 = \lVert U_{m+1}(\beta e_1-\widetilde{H}_my_m) \rVert_2 = \lVert \beta e_1-\widetilde{H}_my_m \rVert_2, \end{equation} $

(3.7)

所以极小化残量范数就等价于求解下面的最小二乘问题

$ \begin{equation} \min\limits_{y\in R^m}\lVert \beta e_1-\widetilde{H}_my_m \rVert_2. \end{equation} $

(3.8)

可以采用基于Givens变换的QR分解^[13]来求解该问题.因为$ \widetilde{H}_m $具有特殊形态, QR分解格外简单, 仅耗费$ O(m) $个运算.这就形成了完整的LSQR方法.求解$ Af = g $时观测数据误差的放大因子是奇异值的倒数, 求解$ A^TAf = A^Tg $时观测数据误差的放大因子是奇异值平方的倒数.所以, 对于病态问题LSQR算法的效果要好, 但当数据误差较大时LSQR算法仍会发散.为此, 根据Tikhonov正则化的思想, 通过对非满秩的矩阵$ A $的协方差矩阵$ A^TA $的每一个对角元素加入一个很小的扰动$ \lambda $, 使得奇异的协方差矩阵$ A^TA $求逆变为非奇异矩阵的$ A^TA+\lambda I $求逆, 从而大大改善求解非满秩矩阵$ Af = g $的数值稳定性, 同时降低了条件数的大小.于是得到最小二乘问题

$ \begin{equation} \min \lVert Af-g \rVert_2^2+\lambda \lVert f \rVert_2^2, \end{equation} $

(3.9)

式中$ \lambda>0 $为正则化参数, 用于控制求解后图像的精度和平滑程度之间的平衡, 对图像恢复质量的好坏有直接的影响.在本文中, 由于原始数据中$ g $的相对误差参数$ \varepsilon $可以近似计算出, 所以用Morozov偏差原理^{[14, 15]}来确定正则化参数$ \lambda $.知道(3.9)等价于增广法方程

$ \begin{equation} (A^TA+\lambda I)f = A^Tg, \end{equation} $

(3.10)

简记为

$ \begin{equation} A_1f = g_1, \end{equation} $

(3.11)

其中$ A_1 = A^TA+\lambda I $, $ g_1 = A^Tg $.然后利用上述LSQR算法求解方程(3.11), 如此便将正则化技术和LSQR方法结合在一起, 各取所长, 构成了求解大规模离散不适定问题的新算法—混合正则化LSQR算法.

根据以上叙述给出混合正则化LSQR算法的具体实现过程.

(1) 赋初值, 输入$ A $和$ g $与$ M $, 计算矩阵$ A $的条件数cond$ (A) $, 如果cond$ (A)\geq M $ (可取$ M = 10^8 $), 则转向(2), 否则$ A_1 = A $, $ g_1 = g $转向(3);

(2) 完成正则化过程:选取正则化参数$ \lambda $ (本文利用Morozov偏差原理确定$ \lambda $的值), 正则化$ A_1 = A^TA+\lambda I $, $ g_1 = A^Tg $;

(3) 进行初始化

$ \begin{aligned} x_0& = 0, \beta_1 = \lVert g_1 \rVert_2, u_1 = g_1/\beta_1, \alpha_1 = \lVert A_1^Tu_1 \rVert_2, \\ v_1& = A_1^Tu_1/\alpha_1, w_1 = v_1, \bar{\phi}_1 = \beta_1, \bar{\rho}_1 = \alpha_1\; ; \end{aligned} $

(4) 对$ i = 1, 2, \cdots $, 重复(5)–(8);

(5) Lanczos双对角化

$ \begin{aligned} w& = A_1v_i-\alpha_iu_i, \beta_{i+1} = \lVert w \rVert_2, u_{i+1} = w/\beta_{i+1}, \\ z& = A_1^Tu_{i+1}-\beta_{i+1}v_i, \alpha_{i+1} = \lVert z \rVert_2, v_{i+1} = z/\alpha_{i+1}\; ; \end{aligned} $

(6) 进行QR分解

$ \begin{aligned} \rho_i& = (\bar{\rho}_i^2+\beta_{i+1}^2)^{1/2}, c_i = \bar{\rho}_i/\rho_i, s_i = \beta_{i+1}/\rho_i, \\ \theta_{i+1}& = s_i\alpha_{i+1}, \bar{\rho}_{i+1} = -c_i\alpha_{i+1}, \phi_i = c_i\bar{\phi}_i, \bar{\phi}_{i+1} = s_i\bar{\phi}_i\; ; \end{aligned} $

(7) 更新迭代求解

$ x_i = x_{i-1}+(\phi_i/\rho_i)w_i\; , \; w_{i+1} = v_{i+1}-(\theta_{i+1}/\rho_i)w_i\; ; $

(8) 进行收敛判别.

为了验证所提方法的可行性, 下面给出数值算例, 模拟时不妨假设原始图像$ f $已知, $ g $可通过离散线性系统$ Af = g $获得, 并对$ g $加入高斯白噪声, 即$ g_n = g+\varepsilon \cdot $ randn$ (n, 1) $, 其中randn$ (n, 1) $是服从方差为1均值为0的高斯分布, 本文取$ \varepsilon = 0.01 $.另外, 用解的相对误差$ \delta $和改进信噪比^[16]描述恢复图像的逼真度, 计算公式分别为

$ \begin{aligned} \delta = \frac{\lVert f-\widetilde{f}\rVert_2}{\lVert f \rVert_2}, ISNR = 10\cdot \text{lg}\left( \frac{\lVert g-f\rVert^2}{\lVert\widetilde{f}-f\rVert^2} \right), \end{aligned} $

其中$ g $为退化图像, $ f $为原始图像, $ \widetilde{f} $为恢复图像. $ \delta $越小, $ ISNR $越大, 图像恢复的逼真度越高, 恢复效果越好.基于上述评价标准, 给出两个测试算法性能的算例.

算例1 选取原始图像大小$ 41\times41 $, 本例中的矩阵$ A $是利用中点公式离散积分公式所得到的, 由Morozov偏差原理确定正则化参数$ \lambda = 3.725e-09 $.图 1给出了原始图像、模糊噪声图像以及用混合正则化LSQR方法去模糊去噪图, 图 2给出了LSQR算法、文献[3]中的正则化Gmres算法($ tol = 0.01 $)和本文新算法的$ ISNR $变化曲线图, 所对应的数值实验结果见表 1.

图 1显示的结果说明混合正则化LSQR方法有效地去除了噪声和模糊, 从视觉效果来看更为清晰.从图 2可以得知, 提出的新算法迭代6次后的$ ISNR $值高于LSQR算法且趋于平稳, 这也就达到了对已有LSQR算法改进的目的, 而$ ISNR $值的提高意味着降低了噪声对信号的影响, 从而得到最优的恢复图像.

通过表 1, 显而易见的是, 本文所提出的新算法$ ISNR $值最高, 相对误差最小, 耗时较短.所以新算法的数值表现总体上均优于LSQR算法和文献[3]所提出的的正则化Gmres算法.

算例2 选取大小为$ 71\times71 $的卫星图作为原始图像, 退化图像的模糊矩阵$ A $依然利用中点公式离散积分公式所得到的, 不同的正则化参数$ \lambda $在图像恢复时会得到不同的改进信噪比$ ISNR $, 具体见图 3.并用LSQR算法、文献[3]中的正则化Gmres算法(tol$ = 0.01 $)和本文新算法的恢复效果如图 4所示, 所对应的数值实验结果见表 2.

由图 3可以确定参数$ \lambda $为$ 10^{-2} $时信噪比最高.在图 4中, 绘制了各算法的$ ISNR $随迭代次数的变化情况.显然, LSQR算法表现出半收敛现象, 而本文算法迭代35次后的$ ISNR $值趋于平稳, 不仅有效地改善了LSQR算法的半收敛性, 还提高了$ ISNR $值, 表现出数值稳定性的优点.图 5可以直观地看出LSQR算法恢复的图像背景仍然被噪声污染严重, 导致其$ ISNR $值很低, 正则化Gmres算法恢复的卫星部分较为模糊, 相比之下本文提出的混合正则化LSQR方法更接近原始图像, 得到结果的图像质量较好, 表现较为稳定.

从表 2的实验结果来看, 虽然本文提出的算法耗时略微慢点, 但是在提高$ ISNR $较和降低相对误差方面的表现突出, 进一步说明了本文提出的混合正则化LSQR算法是相对有效的.

本文在基于投影原理的Krylov子空间方法的基础上, 结合Tikhonov正则化技术提出了一种新算法——混合正则化LSQR算法, 该算法不仅发挥了各自优势, 还有效克服了各自的弊端, 图像的恢复效果也令人满意.实验结果表明, 新算法得到的图像逼真度均优于LSQR算法和正则化Gmres算法, 尤其是在改进信噪比$ ISNR $和相对误差方面有明显的优势, 为求解图像去模糊去噪问题提供了一个行之有效的方法.