0 引言

自适应光学系统(Adaptive Optics，AO)可以消除因大气湍流导致的光波波前畸变，提高光学系统对环境的适应能力.AO系统通过波前传感器测量波前像差，波前控制器根据测量的波前像差生成控制信号，从而调节波前校正器的面形变化，补偿波前畸变^[1].经自适应光学系统实时像差校正后，大部分低阶像差得到补偿.但一般情况下，受到系统成本、有限带宽及探测噪声等限制，自适应光学系统对大气湍流的补偿不完全，目标的高频信息依然受到抑制和衰减^[2].高频信息的丢失导致目标细节特征不清晰，很难满足对空间精确定位和目标识别的要求.对自适应光学校正后的图像进行处理，以得到高质量成像，称为图像事后处理.图像的成像过程可以看成点扩散函数(Point Spread Function，PSF)和目标的卷积，再加上成像噪声，那么图像的去模糊就是去噪解卷积过程.图像复原根据点扩散函数的先验信息是否已知可以分为解卷积(已知)和盲解卷积(未知)^[3].目前已有许多解卷积方法，如能够缓和振铃效应的RL(Richardson Lucy)算法^[4]；能同时减弱振铃和噪声效应的维纳滤波算法^[5]；以及将问题转化为最小能量函数的优化组合求解问题马尔可夫场理论^[6]等.这些方法都能够有效的对图像进行复原，但在实际应用中，点扩散函数很难直接获取，因此研究盲解卷积复原图像的方法更有意义.

图像的盲解卷积技术在图像事后处理方法中应用非常广泛，盲解卷积的关键部分就是对于点扩散函数的估计，点扩散函数的估计会直接影响到图像复原的效果，因此这也是研究的热点之一，文献[7]提出的椭圆抛物面模型点扩散函数估计方法，可以得到较好的点扩散函数模板.文献[8]对近似点扩展函数估计算法进行改进，改进的方法可以较为精确地获得大气湍流模型的点扩散函数.文献[9]在一种新的改进倒频谱域中采用位平面分解提取和Radon变换相结合的方法，准确估计出PSF的模糊角度.文献[7-9]都可以较好估计点扩散函数，完成图像复原，但算法计算量大并且要求在湍流冻结的时间内，对多幅短曝光图像采集和处理，且仅适用于对实时性没有要求的场合.本文提出一种适应大气湍流动态变化的多通道自适应光学图像恢复方法，以自适应光学校正后不同时刻的天文目标或空间目标成像作为多通道，建立求解系统点扩散函数的线性方程.根据解出的点扩散函数，使用基于超拉普拉斯的交替最小化方案求解待观测目标的估计值，从而恢复出图像.本文首先分析不同时刻的点扩散函数的互质性，然后基于建立的线性方程求解点扩散函数，再利用超拉普拉斯算法解卷积，求出目标图像.此方法利用了大气湍流的动态特性，把不同湍流冻结时间内的目标成像看做不同的通道，建立线性方程组，将图像恢复的过程转化成线性求解，大大减少了计算量.

1 理论基础

1.1 多通道盲识别算法原理介绍

天文或空间目标光学成像可表示为输入光场与光学系统点扩散函数的卷积^[10]，表示为

1 $ \mathit{\boldsymbol{I}}\left( {x,y} \right) = \mathit{\boldsymbol{H}}\left( {x,y} \right) \otimes \mathit{\boldsymbol{O}}\left( {x,y} \right) $

式中，I(x, y)是成像，H(x, y)表示点扩散函数，O(x, y)为目标，$ \otimes $为卷积运算，对目标图像进行恢复，即对式(1)的解卷积过程.

多通道盲识别方法算法原理为：采集不同时刻的两幅成像，第i帧和第j帧成像I_i(x, y)和I_j(x, y)，两者之间满足关系

2 $ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{I}}_i}\left( {x,y} \right) = {\mathit{\boldsymbol{H}}_i}\left( {x,y} \right) \otimes \mathit{\boldsymbol{O}}\left( {x,y} \right)}\\ {{\mathit{\boldsymbol{I}}_j}\left( {x,y} \right) = {\mathit{\boldsymbol{H}}_j}\left( {x,y} \right) \otimes \mathit{\boldsymbol{O}}\left( {x,y} \right)} \end{array} $

利用卷积运算的可交换性，得到

3 $ \begin{array}{l} {\mathit{\boldsymbol{H}}_j}\left( {x,y} \right) \otimes {\mathit{\boldsymbol{I}}_i}\left( {x,y} \right) = {\mathit{\boldsymbol{H}}_j}\left( {x,y} \right) \otimes \left[ {{\mathit{\boldsymbol{H}}_i}\left( {x,y} \right) \otimes \mathit{\boldsymbol{O}}\left( {x,y} \right)} \right] = {\mathit{\boldsymbol{H}}_i}\left( {x,y} \right) \otimes \\ \;\;\;\;\left[ {{\mathit{\boldsymbol{H}}_j}\left( {x,y} \right) \otimes \mathit{\boldsymbol{O}}\left( {x,y} \right)} \right] = {\mathit{\boldsymbol{H}}_i}\left( {x,y} \right) \otimes {\mathit{\boldsymbol{I}}_j}\left( {x,y} \right) \end{array} $

显然式(3)中并不含有目标信息，两个通道输出的信息仅与两个点扩散函数有关，因此利用这一特点可以得到关于点扩散函数的线性方程.由数学知识可知，卷积运算可以表示为托普列兹(toeplize)矩阵相乘，将式(3)进行变换，可得一个齐次线性方程为

4 $ \left[ {{\mathit{\boldsymbol{B}}_i} - {\mathit{\boldsymbol{B}}_j}} \right]\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{H}}_j^{\rm{T}}}\\ {\mathit{\boldsymbol{H}}_i^{\rm{T}}} \end{array}} \right] = 0 $

式中，B_i和B_j是由模糊图像I_i(x, y)和I_j(x, y)的toeplize矩阵构成，T表示矩阵的转置，H_j^T和H_i^T为待求解的未知量.在众多求解齐次方程的方法中，选用奇异值分解求解这个方程，那么最小奇异值对应的右奇异向量即为误差最小的解，求出的点扩散函数误差最小.假设待求矩阵R=[B_i－B_j]大小为m×n，矩阵R乘以转置R^T，得到方阵RR^T(m×m)，再对其进行奇异值分解即可求得奇异值和右奇异向量，约需要进行2×m²次乘法和2×m次加法运算.

式(4)有唯一解的前提是点扩散函数互质，因此在将此理论用于自适应光学图像事后处理之前要先判断不同时刻的点扩散函数是否具有互质性，下面利用点扩散函数的相关Sylvester矩阵判断它们之间的互质性.

1.2 不同时刻点扩散函数的互质性分析

文献[11]中给出Sylvester矩阵判断互质性原理：设d(s)和n(s)为两个多项式，且d(s)≠0，不存在多项式p(s)和q(s)使得

5 $ - q\left( s \right)d\left( s \right) + p\left( s \right)n\left( s \right) = - q\left( s \right)p\left( s \right)\left[ {\begin{array}{*{20}{c}} {d\left( s \right)}\\ {n\left( s \right)} \end{array}} \right] = 0 $

且多项式次数degp(s) < degd(s)，即可得d(s)和n(s)互质.

本文利用式(5)判定点扩散函数是否互质，证明如下

设

6 $ \left\{ {\begin{array}{*{20}{l}} {d\left( s \right) = {d_0} + {d_1}s + {d_2}{s^2} + \cdots {d_n}{s^n}}\\ {n\left( s \right) = {n_0} + {n_1}s + {n_2}{s^2} + \cdots {n_m}{s^m}} \end{array}} \right. $

和

7 $ \left\{ {\begin{array}{*{20}{l}} {p\left( s \right) = {p_0} + {p_1}s + {p_2}{s^2} + \cdots {p_{n - 1}}{s^{n - 1}}}\\ {q\left( s \right) = {q_0} + {q_1}s + {q_2}{s^2} + \cdots {q_{m - 1}}{s^{m - 1}}} \end{array}} \right. $

定义多项式d(s)和n(s)的(m+n)阶相关Sylvester矩阵

8

并令s相同次数的系数等于零，可得一个线性齐次方程组为

9 $ \left[ { - {q_0} - {q_1} - {q_2} \cdots - {q_{m - 1}} \vdots {p_0}{p_1}{p_2} \cdots {p_{n - 1}}} \right]\mathit{\boldsymbol{S}} = [ - \mathit{\boldsymbol{QP}}]\mathit{\boldsymbol{S}} = 0 $

此方程有唯一平凡解的前提是Sylvester矩阵非奇异，相应的p(s)=q(s)=0，则满足上述条件，所以d(s)和n(s)互质.

由此可以得出d(s)和n(s)互质的充要条件是它们的相关Sylvester矩阵非奇异.将这一性质运用到式(4)中，那么当H_i和H_j的相关Sylvester矩阵非奇异时就可以认为这两个点扩散函数互质，式(4)有唯一解，即可准确的恢复出目标图像.

1.3 超拉普拉斯解卷积算法

在估测点扩散函数之后，本文采用超拉普拉斯算法解卷积复原图像，Krishnan^[12]提出的超拉普拉斯算法采用交替最小化方案，又称为半二次分裂，一部分用来处理问题的非凸部分，另一部分在频域中解决二次方程问题.算法模型为

10 $ \mathop {\min }\limits_{\mathit{\boldsymbol{x}},\mathit{\boldsymbol{w}}} \sum\limits_i {\left( {\frac{\lambda }{2}\left( {\mathit{\boldsymbol{x}} \otimes \mathit{\boldsymbol{k}} - \mathit{\boldsymbol{y}}} \right)_i^2} \right)} + \frac{\beta }{2}\left( {\left\| {\mathit{\boldsymbol{F}}_i^1\mathit{\boldsymbol{x}} - \mathit{\boldsymbol{w}}_i^1} \right\|_2^2 + \left\| {\mathit{\boldsymbol{F}}_i^2\mathit{\boldsymbol{x}} - \mathit{\boldsymbol{w}}_i^2} \right\|_2^2} \right) + {\left| {\mathit{\boldsymbol{w}}_i^1} \right|^\alpha } + {\left| {\mathit{\boldsymbol{w}}_i^2} \right|^\alpha } $

式中，定义F_i^jx=(x$ \otimes $f_j), j=1, …, J.x为原始图像，y为模糊图像，k是点扩散函数，$ \otimes $为卷积，f为一阶导数滤波器，λ是控制正则化的权值，i是像素，w为每个像素处引入的辅助变量，β是一个在正则化过程变化的权值，α是补偿函数的指数，一般取0~2之间.算法分成两部分进行，给定x通过查LUT表的方法快速求解w，给定w解方程得x.在β增大之前，对x和w交替给定，直到β取得最大值.此算法采用交替迭代最小化方案来优化能量方程，通过查LUT表，对图像快速解卷积，在保证速度的同时也可以达到高质量图像.

1.4 图像质量评价

图像复原后，采用可靠的图像质量评价方法对复原结果进行评价非常有必要^[13].评价方法一般有主观评价和客观评价.主观评价即直接人眼判断图像恢复的质量好坏，这种方法评价的结果因人而异，所以客观评价图像质量的评定方法应用更广泛.常用的有均方误差(Mean Squared Error，MSE)、归一化均方误差(Normalized Mean Square Error, NMSE)、信噪比(Signal Noise Ratio，SNR)和峰值信噪比(Peak Signal to Noise Ratio，PSNR)等^[14].本文采用均方误差对图像质量进行评价.

设两个m×n单色图像I和K，如果一个为另外一个的近似，那么它们的均方误差定义为

11 $ {\rm{MSE}} = \frac{1}{{mn}}\sum\limits_{i = 0}^{m - 1} {\sum\limits_{j = 0}^{n - 1} {{{\left\| {\mathit{\boldsymbol{I}}\left( {i,j} \right) - \mathit{\boldsymbol{K}}\left( {i,j} \right)} \right\|}^2}} } $

均方误差是对有参照物图像复原质量评价的一个最常用的方法，其值越小代表图像误差越小，恢复质量就越好.

2 结果与分析

采用RODDIER N提出的方法^[15]生成不同湍流条件下的畸变波前，分别对大气湍流强度为D/r₀=5，D/r₀=10，D/r₀=15和D/r₀=20时的四种波前像差生成的点扩散函数进行互质性分析(其中D为望远镜孔径，r₀为大气相干长度).基于不同时刻下点扩散函数的互质性，选取两幅原始图像作为目标，对在不同的大气湍流强度下的成像采用超拉普拉斯算法进行恢复，对恢复的结果使用均方误差(MSE)判断图像的恢复质量.

2.1 不同时刻点扩散函数互质性分析

在大气湍流强度为D/r₀=5的情况下生成100帧像差，随机选取两帧波前像差，分别生成大小为32×32的点扩散函数h₁和h₂，如图 1，其中(a)、(b)为波前像差，(c)、(d)为对应的点扩散函数.由1.2节分析可知点扩散函数互质的充要条件是对应的Sylvester矩阵非奇异.将h₁和h₂的每一行生成大小为64×32的toeplize矩阵，再对toeplize矩阵转置形成32×64的矩阵，h₁的每一行与h₂的所有行依次构成64×64的Sylvester矩阵，并验证是否奇异，即构成的Sylvester矩阵行列式是否为0，结果显示构成的所有Sylvester矩阵行列式均不等于0，因此得到这两帧随机两帧像差对应的点扩散函数互质.采用上述互质性分析方法，得到图 2中的两帧随机像差对应的点扩散函数互质.图 2给出了D/r₀=15时随机选取的两帧波前像差和对应的点扩散函数，其中(a)、(b)为波前像差，(c)、(d)为对应的点扩散函数.

图 1. D/r₀=5时随机选取的两帧波前像差和对应的PSFRandomly selected two frames of wavefront aberrations and corresponding PSF under D/r₀=5

Fig. 1. D/r₀=5时随机选取的两帧波前像差和对应的PSFRandomly selected two frames of wavefront aberrations and corresponding PSF under D/r₀=5

下载图片 查看所有图片

图 2. D/r₀=15时随机选取的两帧波前像差和对应的PSFRandomly selected two frames of wavefront aberrations and corresponding PSF under D/r₀=15

Fig. 2. D/r₀=15时随机选取的两帧波前像差和对应的PSFRandomly selected two frames of wavefront aberrations and corresponding PSF under D/r₀=15

下载图片 查看所有图片

将上述测试方法推广到D/r₀=5，D/r₀=10，D/r₀=15和D/r₀=20四种湍流情况下的各100帧随机像差，进行两两测试，运行结果没有出现行列式为0的情况，因此可以认为不同大气湍流强度下的点扩散函数也存在互质关系.

2.2 图像恢复

本文实验基于PC机，在Matlab编程环境中实现.主机CPU型号为Intel(R)Core(TM) i5-4590，主频3.30 GHz，6 GB内存，64位操作系统.采集自适应光学校正后不同时刻的成像，建立齐次线性方程组求解点扩散函数，即式(4)，将图像复原问题转化为线性求解问题，然后采用超拉普拉斯变换方法恢复目标.实验中，超拉普拉斯算法的正则化权值为λ=2×10¹¹，β₀=1，β_max=256，比列增长因子取2$\sqrt 2 $.原始目标如图 3所示，其中(a)为Lena图像，(b)为satellite图像，像素大小为64×64.

图 3. 原始图像Original image

Fig. 3. 原始图像Original image

下载图片 查看所有图片

考虑2个通道的退化图像，在大气湍流强度为D/r₀=5情况下，随机选取两帧波前像差生成点扩散函数，对目标成像，成像结果如图 4所示，其中(a)、(e)是随机选取的两帧像差，(b)、(f)是像差对应的点扩散函数，(c)、(d)是(a)像差下的成像，(g)、(h)是(e)像差下的成像.

图 4. D/r₀=5时随机选取的两帧像差、对应的PSF和目标成像Randomly selected two frames of aberration, corresponding PSF and target imaging under D/r₀=5

Fig. 4. D/r₀=5时随机选取的两帧像差、对应的PSF和目标成像Randomly selected two frames of aberration, corresponding PSF and target imaging under D/r₀=5

下载图片 查看所有图片

复原结果如图 5所示.其中PSF(a)、(d)对应图 4原PSF(b)、(f)；恢复的目标(b)、(c)对应图 4成像(c)、(d)；恢复的目标(e)、(f)对应图 4中成像(g)、(h).从图中可以明显看出，求解出的点扩散函数和退化图像细节已被有效恢复，目视效果接近原始图像.其中求解点扩散函数用时14 s，图像复原2 s.

图 5. 对应图 4D/r₀=5时恢复的PSF和目标Corresponding PSF and target recovered in Fig. 4 under D/r₀=5

Fig. 5. 对应图 4D/r₀=5时恢复的PSF和目标Corresponding PSF and target recovered in Fig. 4 under D/r₀=5

下载图片 查看所有图片

调整大气湍流强度为D/r₀=10，D/r₀=15和D/r₀=20，使用同样的方法对不同大气湍流强度下的成像复原.图 6给出大气湍流强度D/r₀=15时的两帧随机像差，对应的点扩散函数以及成像.

图 6. D/r₀=15时随机选取的两帧像差、对应的PSF和目标成像Randomly selected two frames of aberration, corresponding PSF and target imaging under D/r₀=15

Fig. 6. D/r₀=15时随机选取的两帧像差、对应的PSF和目标成像Randomly selected two frames of aberration, corresponding PSF and target imaging under D/r₀=15

下载图片 查看所有图片

图 7为图 6中PSF和成像复原的结果，其中PSF(a)、(d)对应图 6原PSF(b)、(f)；恢复的目标(b)、(c)对应图 6成像(c)、(d)；恢复的目标(e)、(f)对应图 6中成像(g)、(h).实验中求解点扩散函数用时14 s，图像复原2 s.从复原的效果来看，目视接近原始图像，同样在大气湍流强度为D/r₀=10和D/r₀=20时，也得到了极好的复原结果，初步验证了本文提出的多通道盲识别图像复原方法的可行性.

图 7. 对应图 6D/r₀=15时恢复的PSF和目标Corresponding PSF and target recovered in Fig. 6 under D/r₀=15

Fig. 7. 对应图 6D/r₀=15时恢复的PSF和目标Corresponding PSF and target recovered in Fig. 6 under D/r₀=15

下载图片 查看所有图片

为了评价本文提出方法的有效性，对不同大气湍流条件下求解出的点扩散函数和复原图像采用均值误差(MSE)指标评价其恢复效果.分别对两个目标在大气湍流强度为D/r₀=5，D/r₀=10，D/r₀=15和D/r₀=20时复原，在每一种大气湍流强度下实验20次取MSE平均值，结果如表 1和表 2所示.

从表 1中数据可以看出：点扩散函数的均方误差在10^－30~10^－27量级，证明在不同大气湍流强度下本文建立的求解点扩散函数方程组均可以准确地恢复出点扩散函数.由表 2可得：图像复原的MSE值均在8.2×10^－5~1.1×10^－4之间，随着大气湍流强度的增加，图像复原的误差也在增大，但是波动的误差都非常小，可以认为复原的图像质量极好，能够准确地保持图像基本信息和视觉特性，同时也验证了本文所提出的多通道盲识别图像复原算法在不同大气湍流强度下均适用.实验中点扩散函数计算时长均值约为14 s，解卷积平均用时2 s.以上结果验证了本文提出的多通道盲识别图像复原方法的合理性和可行性，为后续在噪声情况下图像复原提供了理论基础.

3 结论

盲解卷积在自适应光学图像事后处理中凭借不需要知道模糊图像的先验信息而被广泛应用.如何有效地估测点扩散函数，对之后的解卷积具有重要意义.本文利用不同时刻的点扩散函数的互质性，提出多通道自适应光学图像恢复方法，采集自适应光学校正后不同时刻的天文目标成像或空间目标成像，建立求解点扩散函数的线性方程，根据解出的点扩散函数求解待观测目标的估计值，采用超拉普拉斯的交替最小化方案解卷积，从而恢复出图像.分别选择大气湍流强度为D/r₀=5，D/r₀=10，D/r₀=15和D/r₀=20的四种湍流条件下的波前像差生成的点扩散函数分析是否互质，利用Sylvester矩阵分析互质性，超拉普拉斯算法解卷积，并对复原的图像进行图像质量评价.实验结果表明：不同大气湍流强度下的像差所产生的点扩散函数均互质，满足理论要求.基于点扩散函数互质关系所采用的多通道盲识别图像复原方法在不同大气湍流强度下均可以准确地恢复出目标，恢复出的目标与原始目标之间的均方误差在10^－5~10^－4量级，图像质量评价结果验证了本文所提出的多通道盲识别图像复原方法的合理性和可行性.本方法将图像恢复的过程转化成线性求解，大大减少了计算量，并且以上结论为后续在噪声情况下采用多通道盲识别方法进行图像复原提供了理论基础.