由于高分辨率遥感卫星拍摄的遥感图像具有广泛的应用价值，如环境监测、城市规划等.获取的高分辨率遥感图像通常含有丰富的数据信息，但是由于拍摄环境以及长距离传输的影响，接收到的图像数据通常含有大量的噪声，给图像的分析带来了影响.遥感图像的去噪问题一直是遥感应用领域研究的热点，尤其是高噪声遥感图像去噪.图像去噪的目的在于去除噪声的同时最大程度地保护图像的某些重要特征信息.目前，学者提出的遥感图像去噪方法主要分为: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所示.

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

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

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

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

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

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

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

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

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

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

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

4 结论

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

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

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