高噪声遥感图像稀疏去噪重建

删除或更新信息，请邮件至freekaoyan#163.com(#换成@)

本站小编哈尔滨工业大学/2019-10-24

高噪声遥感图像稀疏去噪重建

张杰¹,史小平²,张焕龙¹,耿盛涛¹

(1.郑州轻工业大学电气信息工程学院,郑州 450002; 2.哈尔滨工业大学航天学院,150080 哈尔滨)

摘要:

高噪声遥感图像去噪一直是遥感领域研究的一个重要难题, 为进一步提高高噪声遥感图像的重建质量,在经典的压缩感知迭代小波阈值算法的基础上,提出了一种改进迭代小波阈值算法. 首先,提出一种自适应小波滤波算子在图像稀疏变换过程中对获取的遥感图像小波系数进行筛选,去除图像中的部分噪声信息；其次,使用提出的下降BayesShrink阈值在每次迭代过程中对获取的小波系数进行二次筛选过程；最后,使用改进的块稀疏全变差方法对获得的重建图像进行调整以进一步提高重建遥感图像的质量. 试验结果表明,该算法的去噪重建性能优于经典的压缩感知迭代小波阈值算法,可以从高噪声图像中重建一幅高质量的遥感图像,验证了该算法的有效性. 此外,该算法能够有效地保护遥感图像的边缘和纹理等重要特征信息. 在低压缩采样比情况下,该算法也能够获得相对较高的峰值信噪比和视觉质量. 在卫星地面接收站,该算法可直接使用获取的少量含噪遥感图像数据重建一幅清晰的遥感图像.

关键词: 高噪声遥感图像去噪压缩感知小波阈值改进的块稀疏全变差

DOI：10.11918/j.issn.0367-6234.201806161

分类号:TN911.73

文献标识码:A

基金项目:国家自然科学基金(9,7,61873246)

High noise remote sensing image sparse denoising reconstruction

ZHANG Jie¹,SHI Xiaoping²,ZHANG Huanlong¹,GENG Shengtao¹

(1.School of Electrical and Information Engineering, Zhengzhou University of Light Industry, Zhengzhou 450002, China; 2.School of Astronautics, Harbin Institute of Technology, Harbin 150080, China)

Abstract:

High noise remote sensing image denoising was a difficult problem in the field of remote sensing research. In this paper, to improve the reconstruction quality of remote sensing images, an improved iterative wavelet thresholding (IWT) algorithm is proposed on the foundation of the classical compressed sensing iterative wavelet thresholding (IWT-CS) algorithm. First, an adaptive wavelet filtering operator was proposed to remove parts of image noise, which selects wavelet coefficients of the remote sensing image in the process of image sparsity transform. Second, a descending BayesShrink threshold was put forward to select the wavelet coefficients obtained again in each iteration. Finally, an improved group sparse total variation (IGSTV) method was utilized to adjust the obtained reconstructed image to further improve the reconstruction quality of the remote sensing image. The experimental result demonstrates that the proposed algorithm was superior to the IWT-CS algorithm in terms of denoising reconstruction, which could recover a high quality remote sensing image from high noise image, and the effectiveness of the proposed algorithm was validated. In addition, the proposed algorithm could effectively protect the edges, textures, and other important feature information in remote sensing image. Under low compression sampling ratio, the proposed algorithm could still obtain relatively high peak signal to noise ratio (PSNR) and visual quality. In the satellite receive station, the proposed algorithm can be directly used to reconstruct a clear remote sensing image using a small amount of received-noise remote sensing image data.

Key words: high noise remote sensing image denoising compressed sensing wavelet threshold IGSTV

张杰, 史小平, 张焕龙, 耿盛涛. 高噪声遥感图像稀疏去噪重建[J]. 哈尔滨工业大学学报, 2019, 51(10): 47-54. DOI: 10.11918/j.issn.0367-6234.201806161. 复制到剪切板

ZHANG Jie, SHI Xiaoping, ZHANG Huanlong, GENG Shengtao. High noise remote sensing image sparse denoising reconstruction[J]. Journal of Harbin Institute of Technology, 2019, 51(10): 47-54. DOI: 10.11918/j.issn.0367-6234.201806161. 复制到剪切板

基金项目国家自然科学基金(61427809, 61074127, 61873246) 作者简介张杰(1986—), 男, 讲师;
史小平(1965—), 男, 教授, 博士生导师通信作者史小平, sxp@hit.edu.cn 文章历史收稿日期: 2018-06-26

Contents            -->Abstract            Full text            Figures/Tables            PDF

高噪声遥感图像稀疏去噪重建
张杰¹, 史小平², 张焕龙¹, 耿盛涛¹
1. 郑州轻工业大学电气信息工程学院, 郑州 450002;
2. 哈尔滨工业大学航天学院，150080 哈尔滨

收稿日期: 2018-06-26
基金项目: 国家自然科学基金(61427809, 61074127, 61873246)
作者简介: 张杰(1986—), 男, 讲师;
史小平(1965—), 男, 教授, 博士生导师
通信作者: 史小平, sxp@hit.edu.cn

摘要: 高噪声遥感图像去噪一直是遥感领域研究的一个重要难题, 为进一步提高高噪声遥感图像的重建质量，在经典的压缩感知迭代小波阈值算法的基础上，提出了一种改进迭代小波阈值算法.首先，提出一种自适应小波滤波算子在图像稀疏变换过程中对获取的遥感图像小波系数进行筛选，去除图像中的部分噪声信息；其次，使用提出的下降BayesShrink阈值在每次迭代过程中对获取的小波系数进行二次筛选过程；最后，使用改进的块稀疏全变差方法对获得的重建图像进行调整以进一步提高重建遥感图像的质量.试验结果表明，该算法的去噪重建性能优于经典的压缩感知迭代小波阈值算法，可以从高噪声图像中重建一幅高质量的遥感图像，验证了该算法的有效性.此外，该算法能够有效地保护遥感图像的边缘和纹理等重要特征信息.在低压缩采样比情况下，该算法也能够获得相对较高的峰值信噪比和视觉质量.在卫星地面接收站，该算法可直接使用获取的少量含噪遥感图像数据重建一幅清晰的遥感图像.
关键词: 高噪声遥感图像    去噪    压缩感知    小波阈值    改进的块稀疏全变差
High noise remote sensing image sparse denoising reconstruction
ZHANG Jie¹, SHI Xiaoping², ZHANG Huanlong¹, GENG Shengtao¹
1. School of Electrical and Information Engineering, Zhengzhou University of Light Industry, Zhengzhou 450002, China;
2. School of Astronautics, Harbin Institute of Technology, Harbin 150080, China

Abstract: High noise remote sensing image denoising was a difficult problem in the field of remote sensing research. In this paper, to improve the reconstruction quality of remote sensing images, an improved iterative wavelet thresholding (IWT) algorithm is proposed on the foundation of the classical compressed sensing iterative wavelet thresholding (IWT-CS) algorithm. First, an adaptive wavelet filtering operator was proposed to remove parts of image noise, which selects wavelet coefficients of the remote sensing image in the process of image sparsity transform. Second, a descending BayesShrink threshold was put forward to select the wavelet coefficients obtained again in each iteration. Finally, an improved group sparse total variation (IGSTV) method was utilized to adjust the obtained reconstructed image to further improve the reconstruction quality of the remote sensing image. The experimental result demonstrates that the proposed algorithm was superior to the IWT-CS algorithm in terms of denoising reconstruction, which could recover a high quality remote sensing image from high noise image, and the effectiveness of the proposed algorithm was validated. In addition, the proposed algorithm could effectively protect the edges, textures, and other important feature information in remote sensing image. Under low compression sampling ratio, the proposed algorithm could still obtain relatively high peak signal to noise ratio (PSNR) and visual quality. In the satellite receive station, the proposed algorithm can be directly used to reconstruct a clear remote sensing image using a small amount of received-noise remote sensing image data.
Keywords: high noise remote sensing image    denoising    compressed sensing    wavelet threshold    IGSTV
由于高分辨率遥感卫星拍摄的遥感图像具有广泛的应用价值，如环境监测、城市规划等.获取的高分辨率遥感图像通常含有丰富的数据信息，但是由于拍摄环境以及长距离传输的影响，接收到的图像数据通常含有大量的噪声，给图像的分析带来了影响.遥感图像的去噪问题一直是遥感应用领域研究的热点，尤其是高噪声遥感图像去噪.图像去噪的目的在于去除噪声的同时最大程度地保护图像的某些重要特征信息.目前，学者提出的遥感图像去噪方法主要分为:1)图像域去噪.该类方法直接对整幅图像数据进行处理，典型代表算法为PDE去噪算法^[1-2]、Nonlocal去噪算法^[3-4]等.2)变换域去噪.该类算法先将图像数据投影到某个变换域，将图像数据与噪声信息分离开；随后使用设计的高性能阈值或者滤波器去除变换域的噪声系数进而达到去噪的目的.典型代表算法为：轮廓波变换去噪^[5]、曲波去噪^[6]和小波去噪^[7]等.3)压缩感知(compressed sensing, CS)去噪^[8-9].这类方法将信号经过稀疏变换，然后使用与稀疏变换不相干的测量矩阵从获得稀疏系数选取少量的重要信息进行原始信号的重建.典型的CS重建算法为：迭代收缩阈值(IST)算法^[10-11]、Bayes去噪算法^[12]、梯度下降算法^[13]等.

由于目前采集到的遥感图像通常为高分辨率图像，通常含有海量的数据信息.通过对比可知，前两种去噪方法在重建过程中需要对所有的数据信息进行处理，不仅需要消耗大量的时间，而且在某种情况下很难重建原始图像.CS去噪重建方法仅使用少量的观测数据就可以实现信号的高质量重建，可以有效地解决高维数信号的重建问题.因此，为解决高噪声高分辨率遥感图像的去噪重建问题，本文基于CS理论对高性能的遥感图像去噪重建算法进行设计.

CS理论主要包括：稀疏变换、测量矩阵和重建算法.针对高噪声遥感图像，测量矩阵对经过稀疏变换后获得的稀疏系数进行采样时，获取的观测数据中仍含有较多的噪声信息.因此为了提高去噪效果，本文在图像进行稀疏变换的过程中，加入提出的自适应去噪算子对图像信息进行筛选.测量矩阵对筛选的滤波系数进行观测时，获取的观测数据中噪声信息会相应减少，进而提高了去噪重建效果.

通过对上述CS去噪重建算法进行分析可知，IST算法的实现方式较简单，而且大部分稀疏基都能较容易地融入到IST框架中.因此，本文在IST算法的基础上对提出的去噪重建算法进行设计.为了提高IST算法的高噪声遥感图像的重建效果，经过稀疏变换后，本文提出一种下降BayeShink阈值在迭代过程中对稀疏变换的遥感图像系数进行“二次筛选”过程.阈值去噪虽然能够获得较优去噪效果，但在迭代过程中会导致重建图像出现边缘效应，进而降低重建图像的质量^[14].GSTV方法^[15]可以有效地消除边缘效应，提高图像的重建效果.为提高遥感图像的重建质量，本文在GSTV方法的基础上提出了一种IGSTV方法在迭代过程中对重建的遥感图像进行调节.基于以上技术，本文提出了一种高性能CS高噪声遥感图像去噪重建算法.实验结果表明，本文算法具有较优的高噪声图像去噪重建性能，使用少量的观测值就可重建一幅高质量的遥感图像，验证了本文算法的有效性.

1 CS去噪重建模型对于任一K-稀疏的N×1原始信号f，本文考虑不完全观测的CS去噪重建问题为

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

式中：Φ为M×N(M < < N)测量矩阵; e为噪声信号.由于测量矩阵Φ的行数M远远小于列数N, 导致从观测值y中重建原始信号x是一个欠定线性系统.但是CS理论指出：如果信号f在某一正交基Ψ上是可压缩的或者稀疏的，且测量矩阵Φ与正交基Ψ不相关，则原始信号f可以获得高质量重建.

由于常用的随机矩阵与大部分稀疏基Ψ均不相关，本文选择高斯随机矩阵作为测量矩阵Φ, 同时选取小波变换作为稀疏基Ψ.将小波变换基Ψ应用到式(1)中可得：

$y = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}f + e = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}^{ - 1}}\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}f + e = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}^{ - 1}}s + e,$ (2)

式中，s=Ψf.

l₁范数最小化方法^[16-17]常用来求解式(1)或者式(2)的病态问题, 可表示为

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

或者转化到稀疏变换域Ψ下进行求解:

$\mathop {\min }\limits_s \left\{ {\frac{1}{2}\left\| {y - \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}^{ - 1}}s} \right\|_2^2 + \lambda {{\left\| s \right\|}_{{l_1}}}} \right\},$ (4)

式中:λ为正则化参数.式(3)、(4)中的第1项用来计算观测值与估计值之间的偏差，第2项为原始信号的先验信息.

2 算法研究与设计 2.1 自适应小波滤波算子在CS理论中，如果信号在某一稀疏基Ψ是稀疏的，则可以使用与该稀疏基不相关的测量矩阵Φ对获取的稀疏信号进行采样.针对高噪声遥感图像，本文在图像的稀疏变换过程中加入了自适应滤波算子在不同尺度下对图像信号进行筛选.自适应滤波算子定义η(f)为

$\eta \left( f \right) = \left\{ {\begin{array}{*{20}{l}}{f,f \ge T;}\\{0,f < T.}\end{array}} \right.$ (5)

式中：$T=σ \sqrt {2\log(N×N)} $为统一阈值；σ为噪声标准差.对于高分辨率图像去噪，由于N值过大导致阈值T的值变大，进而影响了图像的去噪效果.遥感图像经过小波多尺度变换后，能量主要集中在第1层，最后一层的能量可以忽略不计.为了能够获得较优的去噪效果，需要尽量减少中间层的能量损失.结合遥感图像经过小波变换后的能量规律，本文提出了尺度调整参数ρ对统一阈值T进行调整：

$\rho = {2^{1 - j}}\sqrt {{{\log }_2}\left( {\frac{{{L_k}}}{J}} \right)} .$

式中：L_k为k层子带长度；J为分解总层数且j∈(1, 2, …, J).因此，本文提出的自适应阈值$T=ρσ \sqrt {2\log(N×N)} $.

综上所述，本文遥感图像的小波稀疏变换过程可表示为：s=η(〈f, Ψ〉).

2.2 下降BayesShrink阈值BayesShrink阈值是Chang于2000年提出的一种收缩阈值去噪方法，常用于含有加性噪声图像的去噪当中.设定σ_f²为原始图像方差，σ²为噪声信号方差，则含噪声图像的方差σ_n²可表示为

$\sigma _n^2 = \sigma _f^2 + {\sigma ^2}.$

对于N×N的原始含噪图像，σ_n²可以通过下式的计算方法得到：

$\sigma _n^2 = \frac{1}{{N \times N}}\sum\limits_{i,j} {s_{i,j}^2} ,$ (6)

式中s_{i, j}为原始含噪图像的小波系数.σ²通常采用经典的鲁棒中值估计方法^[18]进行估计：

${\sigma ^2} = {\left[ {\frac{{{\rm Median}\left| {{s_{i,j}}} \right|}}{{0.674 \;5}}} \right]^2}.$

综上所述，BayesShrink阈值可表示为

$\hat T = \frac{{{\sigma ^2}}}{{{\sigma _f}}} = \frac{{{\sigma ^2}}}{{\sqrt {\left| {\sigma _n^2 - {\sigma ^2}} \right|} }}.$

由于本文处理对象为高分辨率遥感图像，随着图像维数的变大，式(6)计算所得到的噪声方差会变得过小，进而导致BayesShrink阈值变得过大，进而影响了重建算法的去噪能力.本文提出一种调节因子ρ=1/ $\sqrt {n} $在迭代过程中对BayesShrink阈值进行调节.因此，本文提出的下降BayesShrink阈值可表示为

$\hat T = \rho \frac{{{\sigma ^2}}}{{{\sigma _f}}} = \frac{1}{{\sqrt n }}\frac{{{\sigma ^2}}}{{{\sigma _f}}},$

式中n为迭代次数.

2.3 IGSTV方法TV模型在重建过程中，除了能够去除噪声信息之外，还能最大限度地保留图像的边缘和纹理特征信息.TV模型可表示为

${\rm{TV}}\left( f \right) = \sum\limits_{i,j} {\sqrt {{{\left( {{f_{i,j}} - {f_{i - 1,j}}} \right)}^2} + \left( {{{\left( {{f_{i,j}} - {f_{i,j - 1}}} \right)}^2}} \right)} } ,$ (7)

式中:i、j分别为图像像素的位置.从式(7)可知，该TV模型并没有结合图像自身所具有的结构信息，进而导致对图像的细节特征信息重建效果较差.由于大部分图像都具有“块稀疏”特性，因此结合遥感图像的特点，将TV模型与遥感图像的稀疏性质相结合，进而提出了一种IGSTV模型.与传统的TV模型不同，IGSTV模型分别对图像的离散梯度进行水平和垂直方向分组，重建性能得到进一步提高.具体IGSTV模型可表示为

$\begin{aligned} \operatorname{IGSTV}(f)=& \sum_{i, j}\left[\sum_{l=0}^{L}\left|\boldsymbol{D}_{h f}(i+l, j)\right|^{2}+\right.\\ &\left.\left|\boldsymbol{D}_{t f}(i+l, j)\right|^{2}\right]^{1 / 2} .\end{aligned}$

式中:L为图像稀疏块个数, 当L=1时，可以看出IGSTV模型与TV模型等价; D_hf、D_vf分别为垂直、水平方向的稀疏算子，对于二维图像，其微分形式定义为：

${\mathit{\boldsymbol{D}}_{hf}}\left( {i,j} \right) = \left\{ {\begin{array}{*{20}{l}}{{f_{i + 1,j}} - {f_{i,j}},i < N - 1;}\\{0,i = N - 1,}\end{array}} \right.$

${\mathit{\boldsymbol{D}}_{if}}(i,j) = \left\{ {\begin{array}{*{20}{l}}{{f_{i,j + 1}} - {f_{i,j}},j < N - 1;}\\{0,j = N - 1.}\end{array}} \right.$

原始图像的一阶微分表示为：

${\mathit{\boldsymbol{D}}_h} = \left[ {\begin{array}{*{20}{c}}{ - 1}&1&{}&{}&{}\\{}&{ - 1}&1&{}&{}\\{}&{}& \ddots & \ddots &{}\\{}&{}&{}&{ - 1}&1\end{array}} \right],$

${\mathit{\boldsymbol{D}}_v} = \left[ {\begin{array}{*{20}{c}}{ - 1}&{}&{}&{}\\1&{ - 1}&{}&{}\\{}&1& \ddots &{}\\{}&{}&{}&{ - 1}\\{}&{}&{}&1\end{array}} \right].$

令(Df)_{i, j}=(D_hf, D_vf)表示为二维向量.在TV方法当中，当前迭代所使用的TV步长μ通常是由前一次迭代获得的D$\hat f$所确定，其计算表达式为

$\mu = \chi \left( {\mathit{\boldsymbol{D}}\hat f} \right) = M \times \frac{{\sum\limits_{1 \le i,j \le n} {{{\left\| {{{\left( {\mathit{\boldsymbol{D}}f} \right)}_{i,j}}} \right\|}^2}} }}{{{{\left( {\sum\limits_{1 \le i,j \le n} {\left\| {{{\left( {\mathit{\boldsymbol{D}}\hat f} \right)}_{i,j}}} \right\|} } \right)}^2}}},$ (8)

其中(D$\hat f$)_{i, j}的离散梯度矩阵表示如下：

${\left( {\mathit{\boldsymbol{D}}\hat f} \right)_{i,j}} = \left[ {\begin{array}{*{20}{c}}{\left\| {{{\left( {\mathit{\boldsymbol{D}}\hat f} \right)}_{1,1}}} \right\|}& \cdots &{\left\| {{{\left( {\mathit{\boldsymbol{D}}\hat f} \right)}_{1,j}}} \right\|}& \cdots &{\left\| {{{\left( {\mathit{\boldsymbol{D}}\hat f} \right)}_{1,n}}} \right\|}\\ \vdots & \ddots & \vdots & \ddots & \vdots \\{\left\| {{{\left( {\mathit{\boldsymbol{D}}\hat f} \right)}_{i,1}}} \right\|}& \cdots &{\left\| {{{\left( {\mathit{\boldsymbol{D}}\hat f} \right)}_{i,j}}} \right\|}& \cdots &{\left\| {{{\left( {\mathit{\boldsymbol{D}}\hat f} \right)}_{i,n}}} \right\|}\\ \vdots & \ddots & \vdots & \ddots & \vdots \\{\left\| {{{\left( {\mathit{\boldsymbol{D}}\hat f} \right)}_{n,1}}} \right\|}& \cdots &{\left\| {{{\left( {\mathit{\boldsymbol{D}}\hat f} \right)}_{n,1}}} \right\|}& \cdots &{\left\| {{{\left( {\mathit{\boldsymbol{D}}\hat f} \right)}_{n,n}}} \right\|}\end{array}} \right].$

式(8)中的M表示测量矩阵Φ的行向量维数，当M取值很大时，χ(D$\hat f$)则用来收缩权值.当D$\hat f$中存在较多非零元素时，χ(D$\hat f$)则可以产生较大权值.因此，为了避免由于M值过小导致TV步长μ的取值过小，本文给出μ的选取规则如下：

$\mu = \max \left\{ {1,\chi \left( {\mathit{\boldsymbol{D}}\hat f} \right)} \right\}.$

因此，本文提出的IGSTV方法可表示为

${{\hat f}_i} = {f_i} - \frac{{\partial {\rm IGSTV}\left( {{f_i}} \right)}}{{\partial {f_i}}} \cdot \frac{i}{\mu },$

式中i为迭代次数.

2.4 本文算法本文定义遥感图像的去噪重建过程如下：

$\hat f = S\left( f \right) = \theta \left[ {\sum \eta (\langle f,\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}\rangle )\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}} \right].$

式中：$\hat f$为重建图像；〈f, Ψ〉为图像f的稀疏系数；θ(f)为滤波算子，如硬阈值算子(式(5))或者软阈值算子(下式).

$\theta (f) = \left\{ {\begin{array}{*{20}{l}}{f - T,f \ge T;}\\{0,\left| f \right| < T;}\\{f + T,f \le - T.}\end{array}} \right.$

如何选择合适的阈值算子θ(f)以及相对应的阈值T主要取决于所解决的问题，同时如何选择最佳的阈值T是一个亟待解决的问题.经过多次实验，本文选择硬阈值算子和提出的下降BayesShrink阈值$\hat T$在迭代过程中对获得的遥感图像小波系数进行筛选.即

$\theta (f) = \left\{ {\begin{array}{*{20}{l}}{f,f \ge \hat T;}\\{0,f < \hat T.}\end{array}} \right.$

令f_k+1=f_k+Φ^T(y－Φf_k)，则遥感图像的迭代重建过程可表示为

${{\hat f}_{k + 1}} = S\left( {{f_{k + 1}}} \right) = \theta \left[ {\sum \eta \left( {\left\langle {{f_{k + 1}},\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}} \right\rangle } \right)\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}} \right],$

式中k为迭代次数.

综上所述，本文算法实现的具体步骤概括如下.

Step1??初始化过程.迭代索引k=0, 含噪观测值y=Φf+e以及重构遥感图像f_k=0.

Step2??使用获得的含噪观测值y进行迭代计算过程为

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

Step3??将获得的重建图像f_k+1进行稀疏变换过程，同时使用设计的自适应滤波算子定义η(f)在稀疏变换过程中对获得的重建稀疏系数进行筛选：

$s' = \eta \left( {\left\langle {{f_{k + 1}},\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}} \right\rangle } \right),$

式中s′为筛选后的重建稀疏系数.

Step4??经过小波变换过程Ψ重建图像为

${{f'}_{k + 1}} = \sum \eta \left( {\left\langle {{f_{k + 1}},\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}} \right\rangle } \right)\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}.$

Step5??计算下降BayesShrink阈值$\hat T$.

Step6??使用设计的阈值算子θ(f)对重建图像进行筛选为

${f_{k + 1}} = \theta \left( {{{f'}_{k + 1}}} \right).$

Step7??使用IGSTV方法对重建图像f_k+1进行调整为

${f_{k + 1}} = {f_{k + 1}} - \frac{{\partial {\rm IGSTV}\left( {{f_{k + 1}}} \right)}}{{\partial {f_{k + 1}}}} \cdot \frac{k}{\mu }.$

Step8??如果‖f_k+1-f_k‖≤ε, 则进入Step9;否则令k=k+1, 返回Step2，重复Step2~8过程.

Step9??输出重建结果.

3 实验及分析为了验证本文算法的去噪重建性能，将它与指数小波迭代收缩阈值(IST-EW)算法^[19]、GMCA算法^[20]和NLB算法^[21]进行实验对比.由于采集到的遥感图像通常为高分辨率图像，为节省算法的运行时间(running time, RT), 本文从获取的高分辨率遥感图像中取出部分2048×2048的图像作为实验图像.同时，选取峰值信噪比(PSNR)、TEI^[22]、SSIM^[15]和MI^[23]等评价指标对重建图像的质量进行评价.TEI主要用来描述算法对图像边缘信息的重建能力，其值越大表明图像重建质量最优；SSIM用来对去噪前后两幅图像的相似度进行评价，它的值越大表明两幅图像越相似；MI主要用来衡量原图像与重建图像之间的相互性，其值越大表示重建图像保留的原始图像信息特征越多.

设定压缩采样比(compression sampling ratio, CSR)为0.3，高斯白噪声标准差σ_s=25, 对快鸟卫星拍摄的青海玉树震后图像进行重建实验，获得的重建结果如图 1所示.

Fig. 1
点击查看原图

图 1 不同算法的重构结果 Fig. 1 Comparison of the results reconstructed by different algorithms

从重建结果可以看出，针对高噪声遥感图像，本文算法和NLB算法均能获得较优的视觉质量，但本文算法获得的PSNR值相对较高.与IST-EW算法和NLB算法花费的运行时间相比，本文算法的运行时间相对较长，但低于GMCA算法的运行时间.

图 2(a)显示了各个算法随CSR增加获得的PSNR值，可以看出随着CSR的增加，各个算法获得的PSNR值不断增大.与其他算法相比，本文算法能够获得相对较高的PSNR值.图 2(b)显示了不同算法随CSR变化获得的SSIM的值.

Fig. 2
点击查看原图

图 2 不同算法得到的曲线 Fig. 2 Reconstructed curves from different algorithms

可以看出，随着CSR的增加，重建图像与原始图像的相似度越高，也即获得的重建质量也越高.与其他3种对比算法相比，本文算法获得的SSIM值相对较高.

图 3(a)、(b)分别显示了不同算法随CSR增加获得的TEI和MI的值.可以看出，与其他算法相比，本文算法均能获得相对较高的TEI和MI值，也即本文算法能够保留较多的图像的边缘和特征细节信息.

Fig. 3
点击查看原图

图 3 各个算法得到的重建结果 Fig. 3 Reconstructed results of different algorithms

图 4显示了各个算法重建青海玉树震后图像所花费的运行时间.从图 4可以看出，本文算法的运行时间高于IST-EW算法和NLB算法，但低于GMCA算法的运行时间.也即，本文算法的时间复杂度高于IST-EW和NLB算法，低于GMCA算法.但在低压缩采样比情况下，可以看出本文算法的时间复杂度与GMCA和NLB算法的时间复杂度差距较小.

Fig. 4
点击查看原图

图 4 各个算法花费的运行时间 Fig. 4 Running time cost by different algorithms

设定CSR为0.1，σ_s=25对高分1号卫星拍摄的2048×2048西安城图像进行重建，获得的重建结果如图 5所示.

Fig. 5
点击查看原图

图 5 不同算法重构得到的图像 Fig. 5 Images reconstructed by different algorithms

可以看出，在低CSR情况下，本文算法仍可获得较优的重建性能，获得的PSNR相对较高, 此时与其他算法所花费的运行时间的差距较小.因此，CSR在一定范围内越低，本文算法的优势越明显.

图 6显示了当CSR为0.3时，对西安城图像的重建结果.可以看出，随着σ_s的不断增大，图像中的噪声信息不断增加，各个算法获得的PSNR值逐渐降低.与其他算法相比，本文算法随σ_s的增加获得的PSNR变化差距相对较小，证明了本文算法针对高噪声图像具有较优的重建性能.

Fig. 6
点击查看原图

图 6 各个算法得到的重建结果 Fig. 6 Reconstructed results of different algorithms

对更多的2 048×2 048卫星图像进行测试，表 1显示了当噪声标准差为25时，各个算法随CSR变化获得的PSNR值; 表 2显示了当采样率为0.1时，各个算法随噪声标准差变化获得的PSNR值.从表 1、2获得的重建结果可知，针对不同的遥感图像，本文算法仍可获得相对较优的去噪重建性能，进一步验证了本文算法针对高噪声遥感图像去噪的有效性.

表 1

表 1 不同CSR下各个算法获得的PSNR Tab. 1 PSNR obtained from different algorithms under different CSRs dB

图像算法CSR

0.10.20.30.40.50.60.7

高分2号卫星图像GMCA23.7226.5728.9430.2233.1135.4736.38

IST-EW24.0326.8529.3231.4934.2836.6137.25

NLB25.2227.2630.5532.1335.2737.6938.83

本文算法25.9427.8531.0333.4636.2338.7839.06

高分4号卫星图像GMCA23.8126.4329.2231.2933.5734.7636.83

IST-EW23.9726.6229.9531.6333.8935.1437.08

NLB24.3627.4330.3432.5534.2635.6737.68

本文算法24.9727.8630.7533.4635.3337.1238.26

表 1 不同CSR下各个算法获得的PSNR Tab. 1 PSNR obtained from different algorithms under different CSRs

表 2

表 2 表 2不同σ_s下各个算法获得的PSNR Tab. 2 PSNR obtained from different algorithms under different σ_s dB

图像算法σ_s

102025303540

高分1号卫星图像GMCA25.5724.8524.2622.1820.3517.71

IST-EW26.3325.2624.6323.0521.2217.82

NLB26.3725.6825.1223.7421.4718.33

本文算法27.6526.9126.3425.1924.3622.87

高分2号卫星图像GMCA24.6724.1623.7221.0818.3516.66

IST-EW25.6224.8324.0323.2720.0817.14

NLB26.5125.7625.2224.2621.4318.29

本文算法27.6426.8725.9425.4324.0522.07

高分4号卫星图像GMCA25.1224.5523.8122.1620.7118.33

IST-EW25.5924.7623.9722.4820.8818.85

NLB26.1225.4324.3622.6721.0919.36

本文算法26.8725.7524.9724.1622.4820.38

表 2 表 2不同σ_s下各个算法获得的PSNR Tab. 2 PSNR obtained from different algorithms under different σ_s

测试本文算法对遥感图像局部特征图像的重建能力，从拍摄的西安城图像中取出大小为1 024×1 024的西安体育场图像作为测试图像.设定CSR为0.1，σ_s=25，获得的重建结果如图 7所示.从重建结果可以看出，本文算法能够获得相对较多的遥感图像特征，同时花费的运行时间与NLB算法的运行时间相差较小.

Fig. 7
点击查看原图

图 7 两种算法重构结果 Fig. 7 Reconstructed results of two algorithms

同时结合图 6、7可以看出，针对低分辨率遥感图像，本文算法能够花费较短的运行时间重建一幅高质量图像.分辨率越低，本文算法的优势越明显.

4 结论1) 本文提出了一种自适应滤波算子在遥感图像的稀疏表示阶段对图像系数进行筛选，同时使用设计的IGSTV方法对重建遥感图像进行调整.结合以上技术，在迭代小波阈值算法的基础上，提出了一种改进算法.

2) 本文算法与GMCA、IST-EW和NLB算法的对比实验结果表明, 本文算法能够获得相对较优的高噪声遥感图像去噪重建性能，能够获得较高峰值信噪比和保留较多的遥感图像特征信息.

3) 尽管本文算法具有较优的去噪性能, 但是花费的运行时间相对较长, 也即算法的时间复杂度相对较高.因此如何提高算法的重建速度是下一步亟待解决的问题.

参考文献
[1] KIM S. PDE-based image restoration: A hybrid model and color image denoising[J]. IEEE Transactions on Image Processing, 2006, 15(5): 1163. DOI:10.1109/TIP.2005.864184

[2] ZHANG Changjiang, CHEN Yuan, DUANMU Chunjiang, et al. Image denoising by using PDE and GCV in tetrolet transform domain[J]. Engineering Applications of Artificial Intelligence, 2016, 48: 204. DOI:10.1016/j.engappai.2015.10.008

[3] LAUS F, NIKOLOVA M, PERSCH J. A nonlocal denoising algorithm for manifold-valued images using second order statistics[J]. SIAM Journal on Imaging Sciences, 2017, 10(1): 416. DOI:10.1137/16M1087114

[4] HU Xianjun, ZHANG Weiming, LI Ke, et al. Secure nonlocal denoising in outsourced images[J]. ACM Transactions on Multimedia Computing, Communications, and Applications, 2016, 12(3): 1. DOI:10.1145/2886777

[5] LIU Jian, LI Tong, XU Ke, et al. An improved image denoising method based on contourlet transform and neighshrink algorithm[J]. International Journal of Computer Applications in Technology, 2018, 57(2): 206. DOI:10.1504/IJCAT.2018.091646

[6] ANSARI R A, BUDHHIRAJU K M. A comparative evaluation of denoising of remotely sensed images using wavelet, curvelet and contourlet transforms[J]. Journal of the Indian Society of Remote Sensing, 2016, 44(6): 843. DOI:10.1007/S12524-016-0552-y

[7] DING Yin, SELESNICK I W. Artifact-free wavelet denoising: Non-convex sparse regularization, convex optimization[J]. IEEE Signal Processing Letter, 2015, 22(9): 1364. DOI:10.1109/LSP.2015.2406314

[8] ZHANG Jie, SHI Xiaoping. An improved curvelet thresholding denoising algorithm for astronomical image[J]. International Journal of Innovative Computing, Information and Control, 2017, 13(2): 509. DOI:10.24507/ijicic.13.02.509

[9] 李洋.压缩感知在天文图像中的应用研究[D].南京: 南京信息工程大学, 2014: 69
LI Yang. Research of compressed sensing in astronomyimages[D]. Nanjing: Nanjing University of Information & Technology, 2014: 69 http://cdmd.cnki.com.cn/Article/CDMD-10300-1015580062.htm

[10] 张杰, 朱奕, 史小平. 压缩感知的天文图像去噪算法[J]. 哈尔滨工业大学学报, 2017, 49(10): 78.
ZHANG Jie, ZHU Yi, SHI Xiaoping. Compressed sensing denoising algorithm for astronomical image[J]. Journal of Harbin Institute of Technology, 2017, 49(10): 78. DOI:10.11918/j.issn.0367-6234.201609002

[11] CHAMBOLLE A, DOSSAL C. On the convergence of the iterates of the "fast iterative shrinkage/thresholding algorithm"[J]. Journal of Optimization Theory and Applications, 2015, 166(3): 968. DOI:10.1007/s10957-015-0746-4

[12] BASELICE F, FERRAIOLI G, PASCAZIO V. A 3D denoising algorithm based on Bayesian theory[J]. Biomedical Engineering Online, 2017, 16(1): 25. DOI:10.1186/s12938-017-0319-x

[13] SNYDER J C, RUPP M, MVLLER K R. Nonlinear gradient denoising: Finding accurate extrema from inaccurate functional derivatives[J]. International Journal of Quantum Chemistry, 2015, 115(16): 1102. DOI:10.1002/qua.24937

[14] 王蓓, 张根耀, 李智, 等. 基于新阈值函数的小波阈值去噪算法[J]. 计算机应用, 2014, 34(5): 1499.
WANG Pei, ZHANG Genyao, LI Zhi, et al. Wavelet threshold denoising algorithm based on new threshold function[J]. Journal of Computer Applications, 2014, 34(5): 1499. DOI:10.11772/j.issn.1001-9081.2004.05.1499

[15] 赵扬, 汤敏. 基于不同全变差的医学图像压缩感知重构[J]. 计算机工程与设计, 2017, 38(9): 2443.
ZHAO Yang, TANG Min. Compressed sensing reconstruction of medical images based on different total variation[J]. Computer Engineering and Design, 2017, 38(9): 2443. DOI:10.16208/j.issn1000-7024.2017.09.028

[16] TIKHONOV A N, GONCHARSK A, STEPANOV V V, et al. Numerical methods for the solution of ill-posed problems[M]. New York: Springer, 2013: 32. DOI:10.1007/978-94-015-8480-7

[17] CANDES E J, TAO T. Decoding by linear programming[J]. IEEE Transactions on Information Theory, 2005, 51(12): 4203. DOI:10.1109/TIT.2005.858979

[18] KANG Xianggui, STAMM M C, PENG Anjie, et al. Robust median filtering forensics using an autoregressive model[J]. IEEE Transactions on Information Forensics and Security, 2013, 8(9): 1456. DOI:10.1109/TIFS.2013.2273394

[19] ZHANG Yudong, DONG Zhengchao, PHILLIPS P, et al. Exponential wavelet iterative shrinkage thresholding algorithm for compressed sensing magnetic resonance imaging[J]. Information Sciences, 2015, 322(20): 115. DOI:10.1016/j.ins.2015.06.017

[20] YU Chong, CHEN Xiong. Remote sensing image denoising application by generalized morphological component analysis[J]. International Journal of Applied Earth Observation and Geoformation, 2014, 33: 83. DOI:10.1016/j.jag.2014.04.004

[21] MASSE A, LEFEVRE S, BINET R. Fast and accurate denoising method applied to very high resolution optical remote sensing images[C]//Proceeding of 23th Image and Signal Processing for Remote Sensing. Warsaw, Poland: SPIE Remote Sensing, 2017: 1. DOI: 10.1117/12.2277705

[22] QU Xiaobo, ZHANG Weiru, GUO Di, et al. Iterative thresholding comprssed sensing MRI based on contourlet transform[J]. Inverse Problems in Science and Engineering, 2010, 18(6): 737. DOI:10.1080/17415977.2010.492509

[23] QU Guihong, ZHANG Dali, YAN Pingfan. Information measurement for performance of image fusion[J]. Electronics Letters, 2002, 38(7): 313. DOI:10.1049/el:20020212