1 引 言

高光谱影像反映地物空间分布的同时提供高维的光谱特性，为地物的准确分类提供依据^[1]。然而，高光谱影像在采集、转换、传输、压缩和储存等过程中，常因成像设备与外部环境等因素的影响而引入噪声^[2]，导致影像质量的降低，限制了影像的应用^[3]。近年来国内外学者提出了诸多有效的常规图像降噪方法：文献^[4]提出了小波变换软阈值降噪方法；文献^[5]考虑到像元间的空间关系，利用邻域信息实现降噪。高光谱遥感影像包含了比常规图像更多的波段与信息，受到更多因素的影响。文献^[6]探究了利用正交子空间估计高光谱影像中条带噪声的方法；文献^[7]研究了影像各波段间的相关性，求取全局阈值进行小波阈值降噪。在国内，文献^[8]提出用于测深信号处理的最佳小波阈值滤波法；文献^[9]提出了基于PDE的混合遥感影像去噪模型；基于移不变全方向角提升小波^[10]、多指标融合小波^[11]、调制传递函数^[12]等算法也被学者用于去噪处理。

上述影像降噪方法的处理对象通常是一维或二维信号，处理高维数据时则将其视为一系列离散的低维信号，再利用相应的降噪算法进行处理，但人为分割并独立处理各波段导致原影像的相关性信息丢失，无法保证滤波前后影像的空间、光谱一致性。张量^[13]作为描述多维数组的数据模型，在模式分析和机器智能领域有一定应用^{[13, 14]}。相关研究表明：张量方法能够更好地保持遥感影像中空间与光谱的信息，可用于影像降维提高分类精度和目标识别精度等工作^[15]。本文主要介绍一种基于张量的多维滤波算法，针对该算法处理高光谱影像噪声问题进行讨论，通过同两种常规降噪算法进行对比，验证本文算法去除高光谱影像中噪声的同时能够较好地保持影像的空间、光谱信息。

高光谱影像通常被视为多维数组，可以表示成张量X∈R^{I₁×I₂×…×I_N}，其中，N为张量的维度数，也称为阶数；I_i表示第i维的大小。张量中每一个元素表示为X_i1,i2,…,iN(1≤i_n≤I_n且1≤n≤N)，其中，i_n表示该元素在第n阶上的位置。张量是向量、矩阵在维度上的扩展，而后者也是前者在特殊情况下的表达：一阶张量X∈R^I₁表示维度大小为I₁的向量；二阶张量X∈R^I₁×I₂表示行列大小分别为I₁、I₂的矩阵；本文采用三阶张量X∈R^{I₁×I₂×I₃}表示高光谱影像，如图 1所示。I₁、I₂表征影像空间维(第一、二维)的大小，I₃则反映光谱维(第三维)的大小。

图 1 张量模型与高光谱影像的张量表达 Fig. 1 Tensor modeling and hyperspectral image representation

图选项

张量的d维度展开(d-mode unfolding)是将张量X∈R^{I₁×I₂×…×I_N}在第d(d=1,2,…,N)维上展开成矩阵R_d∈R^I_d×M_d的过程。其中，M_d=I₁×I₂×…×I_d-1×I_d+1×…×I_N。

张量的外积(outer product)是一种对张量进行维度扩展的运算，张量X∈R^{I₁×I₂×…×I_N}和Y∈R^{I′₁×I′₂×…×I′_M}的外积为N+M阶的张量，定义为

点击浏览公式

d维度与矩阵相乘(d-mode product)要求张量X∈R^{I₁×I₂×…×I_d×…×I_N}与矩阵U∈R^J×I_d在d阶上维度大小相等，其结果具有参加运算的两个元素的性质，具体数值表达式为：Δ=X_d×U，张量中的元素为

点击浏览公式

本文算法假设无噪声高光谱影像S被高斯加性白噪声N所影响，且噪声与信号相互独立，根据文献^[2]得到含噪影像R的数学模型

点击浏览公式

本文算法对含噪影像R的3个维度分别施加滤波器，获得滤波后影像，表达式如下

点击浏览公式

式中,U、V、W分别为高光谱影像第1、2、3维上对应的多维滤波器，d维度与矩阵相乘则是作用滤波器的具体过程，为实现影像降噪首先需要构建各维度上的滤波器。

降噪算法的目的在于使滤波后影像最大限度地接近无噪声影像，本文利用最小均方误差准则建立降噪模型，使估计影像与无噪声影像间的均方误差最小

点击浏览公式

将式(4)代入展开，表达式可变为

点击浏览公式

为使均方误差最小，可求取式(6)的梯度并令该值为0，得

点击浏览公式

3个滤波器相互关联，其求解可以通过迭代计算实现，不妨假定矩阵V、W已知，来求取滤波器U，结合式(6)和式(7)，于是在第一维上存在关系

点击浏览公式

式中，S₁、R₁分别为张量S、R的第一维展开矩阵。依据矩阵迹的特性，简化式(6)可获得

点击浏览公式

视(V⊗W)^T为权矩阵Ω₁，可以发现基于张量的多维滤波是二维维纳滤波在维度上的扩展。但无噪声影像S通常不可预知，因此有必要研究影像中信号与噪声的能量分布。

模型中噪声与信号相互独立，因此含噪影像在第n维上可以被唯一地分解为噪声和信号两部分^[16]，即含噪影像的能量子空间由相互正交的无噪影像和噪声的能量子空间直和获得

点击浏览公式

无噪影像作为含噪影像的一部分，可以由含噪影像的子空间进行重构获得^[17]

点击浏览公式

式中，S_n是无噪声影像S在维度n上的展开，大小为I_n×M_n；V_s是一组由I_n个大小为1×K_n正交向量构成的特征矩阵；Δ_s是大小为K_n×M_n的权矩阵。K_n是矩阵S_n的秩，决定信号在含噪影像中的比重。本文根据Akaike信息论(AIC)准则^[18]结合噪声水平决定K_n的取值。依据文献^[16]提到的特征值分解重构模型，相应维度上的滤波器Hⁿ在特征空间中可表示成

点击浏览公式

式中，Hⁿ(1≤n≤3)分别对应滤波器U、V、W，Λ_s是大小为K_n×K_n的权矩阵。而影像对应的协方差矩阵的特征值反映了影像中各部分能量的分布^[16]，因此，笔者规定高于某特定阈值的特征值对应影像中有效信息部分，低于该阈值的特征值划归为噪声。在降噪模型中进行特征值分解，可以分别获得E(S_n,R_nΩ_n)和E(R_nΩ_n,R_nΩ_n)的特征值λ_SR、λ_RR。以降序方式重排λ_SR，规定后I_n-K_n个特征值对应于影像中的噪声，取其均值作为噪声强度的估值

点击浏览公式

于是可得含噪影像中的无噪信号对应的特征值序列为；结合式(9)，根据两组协方差包含能量的不同构造权矩阵

点击浏览公式

再取协方差矩阵E(R_nΩ_n,R_nΩ_n)按强度重排后的前K_n个特征向量为V_s，代入式(9)即可求解相应维度上的滤波器。至此完成由两个已知的滤波器求解未知滤波器，再利用新求得的滤波器去更新另外两个滤波器，完成一次迭代。算法的目的是将滤波结果尽可能地接近无噪影像，因此，用上一次迭代获得的滤波影像来替代S，迭代计算滤波器直至收敛。

由上可知，利用张量多维滤波算法进行降噪处理，算法将影像视为整体，分别在3个维度上进行张量展开获得展开矩阵，再根据信号强度在子空间的分布进行特征提取，达到滤波的目的。展开矩阵分别保留了对应维度的信息，迭代过程综合地保持了影像的各维度的信息，因此算法可以更好地保持影像的空间、光谱一致性。

基于张量的多维滤波算法用于高光谱影像降噪的基本流程如下。

(1) 初始化参数：滤波器初始化为单位对角阵，即有U⁰=I₁,V⁰=I₂,W⁰=I₃；初始无噪影像估计为，即有Σ⁰=P。计算第一次滤波影像Σ¹=P₁×U₂⁰×V₃⁰×W⁰。

(2) 若相邻两次迭代的滤波结果的均方差值大于规定阈值ε，则计算滤波器：

1. 利用V^iter-1、W^iter-1构成权矩阵Ω₁^iter。在第一维上展开S^iter为S₁^iter，R为R₁。

2. 计算协方差矩阵SR=E(S₁^iter,R₁Ω₁^iter)，进行特征值分解后从大到小重排得到特征值序列λ_SR，取后I₁-K₁项的均值作为噪声能量。

3. 计算协方差矩阵RR=E(R₁Ω₁^iter,R₁Ω₁^iter)，进行特征值分解后从大到小重排得到特征值序列λ_RR，取前K₁项λ_RR及对应的特征向量构成V_s。

4. 构建权矩阵Λ₁^iter，求得第iter次迭代的第一维滤波器U^iter。

5. 重复以上4步，根据U^iter和W^iter-1计算V^iter，W^iter由U^iter和V^iter计算获得。

6. 依据本次迭代获取的滤波器获得滤波结果，Σ^iter=P₁×U₂^iter×V₃^iter×W^iter。

(3) 若均方误差小于规定阈值ε则判定算法收敛，输出滤波后影像S^iter。

图 2 基于张量的多维滤波算法流程图 Fig. 2 The flow chart of tensor based multiway filtering algorithm

图选项

为验证本文算法的有效性，本文进行了两组试验：试验1选用的高光谱影像数据Indiana Pine是AVIRIS于1992年6月获取，位于Indiana西北100 km²区域，影像大小为128像素×128像素，波段数为220，波段范围为0.4~2.5 μm，光谱分辨率为10 nm，影像含有比较明显的噪声。试验2选用的高光谱影像研究区域位于常州市夏桥，采用的是国产PHI(推扫式光谱成像仪)遥感影像(342像素×581像素)，波段数为80，成像波段范围为0.417~0.854 μm。试验中分别添加了图像信噪比为30 dB、25 dB、20 dB、15 dB、10 dB的高斯白噪声。

对比试验选用逐波段二维维纳滤波算法^[19]和逐波段小波阈值降噪算法^[3]，采用的分析工具为Matlab 7.11.0 ，ENVI 4.7；对比项目主要包括降噪后影像的目视效果、PSNR值、SSIM值、影像分类精度以及影像降噪前后的光谱角与光谱差值曲线等。

在试验中，逐波段维纳滤波采用文献^[19]中的方法；逐波段小波阈值降噪采用atrou小波自适应软阈值降噪；基于张量的多维滤波算法设置维度参数(K₁,K₂,K₃)为(85,85,135)。影像前4个波段含明显噪声，目视对比图如图 3所示。

图 3 Indiana影像降噪结果目视对比 Fig. 3 The visual comparison of denoised images of Indiana

图选项

对比目视效果：基于张量的多维滤波算法可以有效去除影像中大部分噪声，并较好地保留了边缘与细节信息。两种对比算法的降噪结果滤除了噪声的同时模糊了影像的细节。

为检验滤波结果对遥感影像分类的影响，本文选用支持向量机(SVM)为分类器，对降噪后影像在同一组样本下进行分类处理，并对分类结果进行精度评价。在分类试验中，为了避免训练样本选择对分类精度的影响，采取100组独立的分类试验，每次随机选择一定数量的训练样本进行分类，统计平均分类精度。类别、训练样本和测试样本的数量如表 1所示。

各降噪方法的分类结果如图 4所示，分类后精度评价如图 5所示。

图 4 Indiana影像监督分类影像 Fig. 4 The supervised classification result of Indiana

图选项

图 5 降噪后影像分类的总精度和Kappa系数比较 Fig. 5 The comparison of overall accuracy and kappa value of denoised images of Indiana

图选项

由图 5可知，原始影像由于受到噪声影响导致分类精度偏低，经逐波段维纳滤波和小波降噪处理后分类精度有所提高，而经基于张量的多维滤波算法处理后的影像取得了最高的分类精度，达到了88.01％，比原始影像提高了11.04％，Kappa系数也从原始影像的0.728 7提高到了0.857 9。证明本文算法能够有效去除高光谱影像中的噪声，帮助影像提高分类精度。

利用本文算法进行试验，需要设置降噪维度K₁、K₂、K₃，在AIC准则计算获得降噪维度(K₁=K₂=85，K₃=135)的基础上，本文探索了降噪维度和降噪后影像分类精度间的关系，由图 6可知：固定两个维度，分类精度会同第三维的增加而快速上升，达到最大值后基本保持稳定，后期会有缓慢下降。这个结果也反映了当降噪维度取值较小时，影像大部分信息被视为噪声信号，原始信息被大量滤除导致重建质量变差，分类精度降低；当降噪维度取值较大时，一部分噪声信号被当成有效信号，导致滤波后影像包含较多噪声，分类精度下降，只要降噪维度在一定区间内取值就可以保证取得满意的降噪结果。而试验结果也显示使用AIC准则能将每一个维度的值K_n估计到比较准确的区间。

图 6 降噪维度K1、K2、K3与降噪后影像总体分类精度关系曲线 Fig. 6 The relationship between the denoising rank K1,K2,K3  and the overall accuracy of the denoised images

图选项

模拟试验研究了基于张量的多维滤波算法在不同强度噪声影响下的降噪结果，图 7显示了信噪比强度为10 dB下的夏桥影像的目视降噪结果。

图 7 夏桥影像(SNR=10 dB)降噪结果第4波段目视对比< Fig. 7 The visual comparison on the 4th band of denoised images of imnoised Xiaqiao(SNR=10 dB)

图选项

对降噪前后的各组影像进行分类试验，分类后精度评价见表 2。

为评价各算法对影像的空间、光谱信息的保持能力，本文选用4种国际常用的评估参量。峰值信噪比(PSNR)和结构相似性指数(SSIM)^[20]是衡量降噪后影像与无噪声影像空间相似性的经典参数，试验统计了夏桥影像各算法降噪后的PSNR值和SSIM值，由表 3可知：较对比算法，本文算法结果更接近理想影像，有助于保持影像空间信息。

本文采用光谱差值曲线和光谱角比较各降噪结果与无噪声影像的光谱一致性^[3]：在影像各类地物中各取一像元，获得该位置在无噪声影像和各降噪后影像中的光谱，再以横轴为波段，纵轴为差值计算降噪结果和无噪声影像间的光谱差值，图 8展示了在信噪比为10 dB强度噪声影响下的光谱差值曲线，该曲线可以定性衡量降噪算法对影像光谱特征的保持能力；计算各降噪结果与无噪声影像的光谱角来定量衡量光谱一致性，结果见表 4。

图 8 夏桥影像(SNR=10 dB)降噪后影像与无噪影像间的光谱差值曲线 Fig. 8 The difference between spectrum of denoisd images of Xiaqiao(SNR=10 dB) and the noise-free image

图选项

分析图 8和表 4可以得出如下结论：对比常规逐波段降噪算法，基于张量的多维滤波算法取得的降噪结果更加接近理想影像，在保持影像光谱一致性上具有较大优势。

算法也被试验证明具有较好的收敛特性，图 9(a)和(b)分别描述了Indiana影像和添加信噪比强度为10 dB的夏桥影像在各自较优的3个降噪维度参数下的算法收敛情况，图中纵轴为相邻两次迭代获取的滤波影像的差值，横轴为迭代次数。从图 9(a)和(b)可知，本文采用的迭代算法能在5次迭代以内收敛。进一步的试验证明，本文算法能够适应不同的高光谱数据和不同的参数设置情况，迭代计算迅速收敛达到最优解。

图 9 降噪试验迭代关系图 Fig. 9 The convergence of denoising experiments during iterations

图选项

基于张量的多维滤波算法可以有效地用于高光谱影像降噪处理，相对于经典的逐波段降噪算法，该方法在滤波过程中可以有效保持原始影像的空间、光谱一致性，从而取得相对较好的目视降噪效果、帮助提升影像PSNR值以及降噪后影像的分类精度。