对天文图像进行采集和研究是深空探索的一个重要分支.从获得的天文图像中可以直接得知某一星体的地貌特征, 是否存在未知生命等重要天文信息.由于采集得到的天文图像一般为高分辨率图像, 这些图像的存储会给卫星或者其他探测设备自身携带的存储空间带来较大压力.更为重要的是, 天文图像在采集和传输的过程中经常受到噪声信号的干扰, 从接收到的天文图像中很难分辨出某些重要的特征信息, 大部分接收到的图像都要经过二次处理.然而, 由于高分辨率图像的数据量较大, 目前存在的大部分去噪算法处理高维数据的能力有限, 在图像重构的过程中经常由于很难获得足够的原始信号信息导致重构图像的质量较差^[1-3].

为有效地解决高维信号的重构问题, 国内外学者一直都在努力探索如何从高维信号中获得一种低维的信号结构, 并且同时保证这种低维结构保存了重构原始信号所需要的全部数据, 然后利用这种低维数据精确重构原始高维信号.即, 对高维信号进行压缩, 舍弃信号中非重要信息, 直接使用信号中的重要信息进行信号重构.信号的稀疏表示^[3-4]是目前探索信号内部低维结构的一种比较简单的方式之一.

基于稀疏的思想, Donoho^[5]在2006年提出了著名的信号采样理论:压缩感知理论^[1-6].该理论指出:如果某一信号是稀疏的或者是可压缩的, 则该信号就可以从少量的观测值中精确恢复.换言之, 如果某一信号在某个正交基上是稀疏的, 则可以使用一个低维测量矩阵对原始信号进行采样以获得低维数的信号观测值.最后使用重建算法从低维观测值中精确重构原始信号.

CS理论充分利用了信号的稀疏特性, 并且在信号采样的同时就已经完成了信号的压缩过程.传统的香农采样定理要求信号的采样速率必须大于或者等于2倍的信号带宽, 而CS理论使用低于2倍信号带宽的采样速率仍可精确重构原始信号.

CS理论包含3个重要部分:稀疏基的设计, 测量矩阵的选取以及重构算法的设计.本文主要关注设计高性能去噪重建算法.经过近几年的发展, 学者提出了许多的CS去噪算法, 如迭代收缩阈值类 (iterative shrinkage-thresholding, IST) 算法^[7-8]、匹配追踪类算法^[9]、梯度方法^[10]和Bregman迭代方法^[11]等.通过分析可知, 大部分CS去噪重建算法采用阈值的方法对原始信号进行去噪处理.然而, 由于经常使用的稀疏基如小波变换、轮廓波变换等缺乏平移不变性^[12-13], 采用阈值去噪会导致重构的图像中出现伪吉布斯效应.因此, 本文将小波维纳滤波算子替代阈值算子对图像的小波系数进行筛选, 进而提出了基于小波维纳滤波算子的CS重建算法.使用本文算法对高分辨率天文图像进行去噪重建处理, 实验结果表明, 本文算法能有效地提高天文图像的重构质量, 获得较高的峰值信噪比.

1 压缩感知去噪模型

经典的去噪问题可描述为

$ y = f + e. $

式中: f为任意N×1信号；e为加性噪声, 本文考虑为零均值高斯白噪声；y为含噪观测值.对于高维数信号, N值越大, 则观测值y的维数就越高, 在图像重构时, 从y中寻找有效的原始先验数据就越困难.

本文使用低维M×N(M≪N) 非自适应矩阵 (CS测量矩阵) Φ对信号f进行观测, 进而得到低维含噪观测值y:M×1, 可描述为

$ y = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}f + e. $

(1)

由于M≪N, 因此从y中重构原始信号f是一个病态问题.然而由压缩感知理论可知:如果信号在某个正交基Ψ上是稀疏的, 同时测量矩阵Φ与正交基Ψ不相干, 则原始信号就可以获得高概率重建^[1-4].同时, 式 (1) 可以转化为

$ y = \mathit{\boldsymbol{ \boldsymbol{\varPhi} \boldsymbol{\varPsi} }}s + e = \mathit{\boldsymbol{ \boldsymbol{\varTheta} s}} + e, $

式中s = Ψ^-1f为稀疏系数. Θ = ΦΨ可以作为CS测量矩阵直接对稀疏系数进行观测.本文将正交小波作为信号稀疏基Ψ.

正则化方法^[14-15]经常用来求解以上病态问题, 而这类方法通常转化为满足某些稀疏限制条件的最小l_p范数问题进行求解, 即

$ \mathop {\min }\limits_f \left\| {y - \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}f} \right\|_2^2 + \lambda {\left\| {{\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}^{ - 1}}f} \right\|_{{l_p}}}. $

式中:λ为正则化参数; 第1项为惩罚项, 用来估计观测值与计算值之间的偏差; 第2项为正则化项, 它代表原始信号的先验信息；Ψ^-1可以看作是正则化运算算子, 根据针对不同先验信息的Ψ^-1, 学者已经提出了许多的正则化方法, 如全变差正则化.本文将CS重构问题转化为最小l₁范数问题来解决, 可表示为

$ \mathop {\min }\limits_f \left\| {y - \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}f} \right\|_2^2 + \lambda {\left\| {{\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}^{ - 1}}f} \right\|_{{l_1}}}. $

2 本文重构算法

在压缩感知中, IST算法的设计过程比较简单且易于实现, 其迭代过程可描述为

$ {f_{k + 1}} = {H_{\rm{W}}}\left( {{f_k} + {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\rm{T}}}\left( {y - \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{f_k}} \right)} \right), $

式中H_W(·) 为阈值算子, 通常考虑为硬阈值算子, 可表示为

$ {H_{\rm{W}}}\left( {{f_k}} \right) = \left\{ \begin{array}{l} {f_k},\left| {{f_k}} \right| \ge T;\\ 0,\left| {{f_k}} \right| < T. \end{array} \right. $

(2)

式中：$T = \sigma \sqrt {2\log \left( {N \times N} \right)} $为硬阈值；σ为噪声标准差.对于高分辨率图像, 阈值T的值会由于N×N的值较大而变的过大, 在迭代重构的过程中丢失较多的图像信息.

此外, 从式 (2) 可以看出硬阈值函数在±T处存在断点, 在图像重构的过程中会产生较大的震荡效应, 最终导致重构的图像中出现伪吉布斯现象.为消除阈值函数的缺陷, 本文将小波维纳滤波算子代替硬阈值算子对天文图像信息进行筛选.

2.1 小波维纳滤波算子

对原始图像进行随机投影得到的观测值为独立同分布的高斯信号, 可以认为是在原始清晰图像中加入了高斯白噪声^[16-17].小波系数通常被认为是条件独立的高斯随机变量, 而噪声信号则是统计独立的零均值高斯变量^[18].因此原始含噪图像的小波系数a_n可表示为

$ {a_n} = {a_p} + \zeta . $

式中：a_p为纯净图像的小波系数, ζ为噪声的小波系数.由于维纳滤波器是一种最小均方误差滤波器, 因此滤波器系数c可以表示为

$ c = \frac{{{\theta ^2}\left( k \right)}}{{{\theta ^2}\left( k \right) + \hat \sigma _n^2}}, $

式中θ²(k) 为小波系数方差, 由于小波系数与噪声信号两者相互独立, 可以采用5×5 LAWMAP方法^[18]进行计算, 即

$ \begin{array}{l} {\theta ^2}\left( k \right) = \max \left( {\frac{Q}{{4\lambda }}\left[ { - 1 + } \right.} \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\left. {\sqrt {1 + \frac{{8\lambda }}{{{Q^2}}}\sum\limits_{i \in N\left( k \right)} {a_n^2\left( i \right)} } } \right] - \hat \sigma _n^2,0} \right). \end{array} $

式中:Q为N(k) 中的系数数量；N(k) 为以a_n(k) 为中心的5×5邻域.参数λ可以采用文献[18]中的方法进行计算.

${{\hat \sigma }_n}$为噪声小波系数的标准差, 在实际当中, 噪声是未知的, 本文采用鲁棒性中值估计方法进行估算, 可表示为

$ {{\hat \sigma }_n} = \frac{{{\rm{Median}}\left( {\left| {{a_n}} \right|} \right)}}{{0.6745}}. $

本文设计的小波维纳滤波算子H_S(·) 为

$ {H_{\rm{S}}}\left( {{a_n}} \right) = \frac{{{\theta ^2}\left( k \right)}}{{{\theta ^2}\left( k \right) + \hat \sigma _n^2}}{a_n}. $

2.2 去噪算法

在图像去噪中, 全变差 (total variation, TV) 方法^[19-20]能够有效地抑制图像中的伪吉布斯现象, 提高重构图像质量.因此, 这里将全变差方法应用到本文的去噪算法当中, 对重构的图像进行调节.

本文算法的具体步骤总结如下：

1) 初始化迭代次数m=1和重构图像f_m=0.

2) 对原始图像进行随机投影观测, 得到含噪观测值y= Φf+e.

3) 计算小波维纳滤波算子H_S(·).

4) 使用得到含噪观测值y迭代重建图像, 可描述为

$ {f_w} = {f_{m - 1}} + {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\rm{T}}}\left( {y - \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{f_{m - 1}}} \right). $

5) 将得到的图像f_w进行小波变换, 并使用小波维纳滤波算子H_S(·) 对图像的小波信息进行筛选, 可描述为

$ {f_s} = {H_S}\left( {{\mathit{\Psi }^{\rm{T}}}{f_w}} \right). $

6) 通过小波逆变换f_m= Ψf_s, 重建去噪后的图像.

7) 如果满足迭代停止条件‖f_m-f_m-1‖≦, 则进入步骤8);否则令m=m+1, 重复步骤4)~6), 进入下一次迭代过程.

8) 利用式 (3) 对去噪重建的天文图像f_m进行TV调整.

$ {f_m} = {f_m} - \frac{{\partial {\rm{TV}}\left( {{f_m}} \right)}}{{\partial {f_m}}}.\frac{q}{m}, $

(3)

其中q表示TV调整步长.同时,

$ {\rm{TV}}\left( {{f_m}} \right) = \\ \sum\limits_{i,j} {\sqrt {{{\left[ {{f_m}\left( {i + 1,j} \right) - {f_m}\left( {i,j} \right)} \right]}^2} + {{\left[ {{f_m}\left( {i,j + 1} \right) - {f_m}\left( {i,j} \right)} \right]}^2}} } . $

9) 输出调整后的重构图像.

3 实验及分析

本文首先设定压缩比为3, 高斯白噪声标准差σ_n=10, 并引入峰值信噪比 (PSNR) 对重构图像的质量进行评价.由于重构一幅高分辨图像需要较长的时间, 为节省重构时间 (reconstruction time, RT), 本文从大小为5 120×3 840的高分辨率图像中取出部分1 024×1 024的图像作为实验图像.为测试本文算法的去噪能力, 将它与迭代硬阈值算法 (IHT)、改进的正交匹配追踪算法 (IOMP)^[21]和快速重建去噪算法^[22]对比.首先测试拍摄的月球图像, 图 1显示了不同算法的重构结果.

从图 1可以看出, 与其他3种算法相比, 本文算法能够恢复一幅清晰的天文图像, 同时获得的PSNR值较高, 但是花费的重构时间较长.

设定σ_n=10, 图 2显示了不同算法随压缩比增加得到的PSNR和RT变化曲线.

从图 2可以看出, 随着压缩比的增大, 不同算法的重构性能逐渐降低.但与其他算法相比, 本文算法能获得相对较高的PSNR值, 同样也需要相对较长的重构时间.从图 2中也可以注意到, 当压缩比较高时, 本文算法仍可获得相对较高的PSNR值.

测试1 024×1 024土星图像, 设定压缩比为10和σ_n=10, 图 3显示了不同算法重构得到的图像.从图 3可以看出, 当压缩比较高时, 本文算法仍具有相对较好的去噪重建性能.

测试更多的天文图像, 表 1显示了当σ_n=10时, 不同算法随压缩比降低时得到的PSNR.表 2显示了当压缩比为3时, 不同重构算法随σ_n增加时重构得到的PSNR值.

从表 1也可以看出, 针对不同的天文图像, 本文算法仍然具有较好的去噪重建能力.从表 2可以得知, 随着σ_n的增加, 不同算法的去噪能力逐渐降低, 但与其他算法相比, 本文算法获得的PSNR值相对较高.

测试本文算法对天文图像局部特征的重构能力, 从高分辨率月球图像中截取大小为512×512的局部特征图像作为实验图像.设定压缩比为3和σ_n=10, 图 4显示了文献[22]算法和本文算法的重构结果.从重构图像可以看出, 本文算法能够恢复较多的细节和纹理特征.

4 结论

1) 本文将小波维纳滤波算子替代压缩感知迭代阈值算法中的小波阈值对高分辨率天文图像进行筛选, 并使用全变差方法对获得的重构图像进行调整, 进而提出了基于小波维纳滤波的迭代阈值改进算法.

2) 本文算法与IHT、IOMP和文献[22]算法的对比结果表明, 本文算法能够获得较优的重建性能, 同时可以有效地保护高分辨率天文图像的细节特征.当压缩比较高时, 该算法也能获得相对较高的峰值信噪比.

3) 尽管本文算法具有较优的去噪性能, 但是消耗的重构时间相对较长, 因此如何降低算法的重构时间是以后研究的重点.