CN105374054A - 基于空谱特性的高光谱图像压缩方法 - Google Patents

Abstract

本发明涉及高光谱遥感图像压缩技术领域，提出了一种基于空谱特性的高光谱图像压缩方法。首先将高光谱图像分成若干块，然后利用高光谱图像的强谱间相关性对每一块高光谱图像的所有波段进行自适应聚类分组，再对每一个高光谱图像块中的每一个波段组分别进行PCA运算,得到若干经过降维后的主成分，最后对所有主成分进行编解码，得到重构图像。本发明提供的高光谱遥感图像压缩方案，针对现有压缩算法对高光谱图像的特性利用不够充分以致压缩性能受限，算法复杂度偏高等问题，充分利用高光谱图像的空间及谱间特性，再结合PCA算法，可以降低高光谱数据维度，减少了算法的计算量，并且前面所诉对波段进行聚类分组也可以降低PCA的运算量。

Description

基于空谱特性的高光谱图像压缩方法

技术领域

本发明属于图像压缩技术领域，具体涉及一种基于空谱特性并结合JPEG2000的图像数据压缩方法。

背景技术

高光谱遥感(Hyperspectralremotesensing)是近几十年才发展起来的一门新兴学科。成像光谱仪在航空、航天载体上对地物进行观测，记录下数百个波段的高光谱遥感图像数据。和传统的遥感图像比，高光谱遥感图像具有图谱合一的特点，同时具有高空间分辨力和光谱分辨力，在农业、地质、军事侦察等方面具有广泛的应用前景。

随着成像光谱技术的飞跃发展，高光谱数据的光谱分辨力和空间分辨力越来越高，数据量化的深度也有所增加，同时卫星不间断周期性扫描次数也在不断增加。这些因素都导致了高光谱遥感图像数据量的激增。典型AVIRIS图像大小为614×512，成像波段数为224，如果每像素灰度值用16比特来表示，一组图像的数据量则约为140M字节。高光谱遥感图像巨大的数据量给高光谱遥感图像的传输和存储带来巨大压力，严重制约了高光谱遥感技术的发展。对高光谱遥感图像进行压缩可以：(1)较快地传输各种信源，降低信道占用费用。(2)在现有的通信干线上开通更多的并行业务。(3)降低发射机的功率。(4)紧缩数据存储容量，降低存储费用。因此，高光谱图像压缩技术成为高光谱遥感领域的研究热点。

由于高光谱图像要用于诸如卫星勘探、制导等领域，另一方面受到卫星上所用计算机以及数据传输设备等限制，故高光谱遥感图像的压缩技术必须满足：(1)高保真度，要求解压后得到的恢复图像与原始图像相比，不影响分析和识别的应用。(2)高压缩比，能够适应分辨力与采样率不断提高的情形下高光谱遥感图像的实时传输与存储的要求。(3)算法简便、耗内存低，高光谱遥感图像数据量巨大，压缩算法的设计力争做到软件高效、硬件低成本，以便设计出实时压缩处理系统。

与传统遥感图像和二维图像相比，高光谱遥感图像存在着更多的冗余信息，这样就使高光谱图像数据压缩成为了可能，为高光谱数据的压缩提供了发挥空间。针对高光谱遥感图像而言主要存在两种冗余——谱间冗余和空间冗余。因此对高光谱图像数据的压缩主要就是有效去除光谱数据间存在的空间冗余信息和谱间冗余信息。

目前，高光谱图像压缩算法基本上分为三类：基于预测的编码方法、基于矢量量化技术的编码方法、基于变换的编码方法。其中基于预测的压缩方法主要用于无损压缩，压缩比受到较大限制，对传输带宽依然有很高要求，不利于数据的实时传输；基于矢量量化压缩方法的算法复杂度太高，计算量随着矢量维数的增加呈指数增长。基于变换的压缩方法对高光谱图像的特性利用还不够充分。

综上所述，现有技术存在以下问题：

1)对高光谱图像的空间及谱间特性利用不够充分，压缩性能无法进一步提升。2)算法复杂度偏高，不利于硬件实现。

发明内容

本发明针对现有高光谱图像有损压缩算法对图像特性利用不够充分，图像恢复质量低以及算法复杂度高的缺点，提出了一种充分利用高光谱图像空间特性和谱间特性，并结合JPEG2000(基于小波变换的图像压缩标准，由JointPhotographicExpertsGroup组织创建和维护)的高光谱图像压缩方法。

本发明解决上述技术问题的技术方案为：一种基于空谱特性的高光谱图像压缩方法，包括步骤：利用高光谱图像的空间相关性以及地物在空间上分布的连续性对高光谱图像分块；对每个高光谱图像块进行自适应波段聚类分组；对每个高光谱图像块的每个波段组分别进行主成分分析PCA(PrincipalComponentAnalysis)运算，得到经过若干降维后的高光谱图像主成分，对所有高光谱图像主成分进行编解码，得到重构图像。

本发明其中一个实施例包括，所述对高光谱图像分块进一步包括：读取高光谱图像，根据空间维上像元的行数M、列数N、波段数得到一个高光谱数据立方体D；用像元间的欧式距离衡量像元间的相似度，遍历空间维上第一行像元，如果相邻像元间的欧式距离小于等于阈值μ，则两像元属于同一块，如果相邻像元间的欧式距离大于μ，则两像元属于不同的块；遍历空间维上非第一行像元，比较相邻像元的相似度，选择相似度最大的像元以及对应像元之间的欧式距离，若该欧式距离小于μ，则两像元属于同一块。

本发明其中一个实施例包括，所述对高光谱图像块进行自适应波段聚类分组进一步包括，根据所有波段之间的相似度组成相似度矩阵，根据公式确定聚类的能量函数，其中，C_l为波段l的聚类中心，S(l,C_l)为波段l与其聚类中心波段C_l的相似度；根据任意两波段l、k的相似度调用公式 $p=\underset{l\≠k,l,k=1,...,N}{min}S(l,k)$ 确定加权偏好参考度p。

本发明其中一个实施例包括，读取每个高光谱图像块中的高光谱数据，将高光谱数据重新排列为矩阵X，矩阵X中的每一行表示一个波段的数据，矩阵X中的每一列为一个L维矢量；根据公式计算每个波段的样本均值，根据公式C_l,i＝E[(X_i-m_X)(X_i-m_X)^T]计算任意光谱的协方差矩阵，根据公式获得所有光谱矢量的协方差矩阵；计算协方差矩阵的特征值和特征向量，对矩阵X进行降维，选取大于阈值的特征值所对应的特征向量组成特征矩阵，用矩阵X乘以特征矩阵得到经过PCA运算后的高光谱图像主成分。

存在正交矩阵U＝[u₁,u₂,…,u_L]，使得所有光谱矢量的协方差矩阵C_l可对角化，即对角矩阵的对角元素{λ₁,λ₂,…,λ_L}为矩阵C_l的特征值，u_i为矩阵C_l的第i个特征值对应的特征向量。

本发明其中一个实施例包括，遍历整个空间范围的高光谱图像，建立一个M行、N列的分块矩阵A，将矩阵A中元素数值相同的位置所对应的高光谱数据立方体D中的像元块归为同一块，将高光谱数据立方体D分成若干块，其中，矩阵A的元素A_ij表示在空间坐标(i,j)处光谱D_ij所属的块数。

本发明其中一个实施例包括，确定波段间的相关系数为波段间的相似度，以任一图像块E_n所有波段之间的相似度组成相似度矩阵S，根据公式 $R\[l,k\]=\frac{\underset{i=1}{\overset{M}{\Σ}}\underset{j=1}{\overset{N}{\Σ}}\left(X_k\right(i,j)-{\mathrm{\μ}}_k)\left(X_l\right(i,j)-{\mathrm{\μ}}_l)}{\sqrt{\underset{i=1}{\overset{M}{\Σ}}\underset{j=1}{\overset{N}{\Σ}}{\left(X_k\right(i,j)-{\mathrm{\μ}}_k)}^2\underset{i=1}{\overset{M}{\Σ}}\underset{j=1}{\overset{N}{\Σ}}{\left(X_l\right(i,j)-{\mathrm{\μ}}_l)}^2}}$ 获得波段l和波段k之间的谱间相关系数，其中，X_k(i,j)和X_l(i,j)分别为图像块中第k波段和第l波段上空间坐标(i,j)的像素，μ_k和μ_l分别为波段k和波段l上所有像素的均值，M和N分别表示图像块空间维上的行数和列数；根据公式：确定加权偏好参考度λ，根据公式： $r(l,k)\←s(l,k)-\underset{k^{\′}\≠k}{max}\{a(l,k^{\′})+s(l,k^{\′})\},$ r(l,k)←λ*r(l,k)_old+(1-λ)*r(l,k)_new计算波段间的吸引度并更新；根据公式： $a(l,k)\←\min \{0,r(k,k)+\underset{l^{\′}\∉\{l,k\}}{\Σ}\max \{0,r(l^{\′},k)\left\}\right\},$ a(l,k)←λ*a(l,k)_old+(1-λ)*a(l,k)_new计算波段间的归属度并更新，其中，s(l,k)表示波段l与波段k间的相似度，s(l,k')表示波段l与波段k'间的相似度，a(l,k')表示波段l对波段k'的归属度，r(l,k)表示波段k对波段l的吸引度，r(l,k)_old表示波段k对波段l前一次的吸引度，r(l,k)_new代表波段l与波段k当前的吸引度，a(l,k)_old表示波段k对波段l前一次的归属度，a(l,k)_new代表波段k对波段l当前的归属度，不断进行迭代更新，根据当前得到的吸引度与归属度，当r(l,k)+a(l,k)最大时，获得波段l的中心波段k，根据中心波段将所有波段分成若干组。

加权偏好参考度用来反映各波段最终被选作类代表波段的统一偏好程度，聚类中心波段与该类中其他任何波段的相似度必须大于加权偏好参考度，其数值越大，得到的类越多，反之越少。算法初始时将高光谱图像块的所有波段都视为潜在的聚类中心，通过一种迭代循环不断地更新，产生若干个聚类中心波段及其对应的聚类，同时聚类的能量函数也得到了最小化。将各波段分配到对应中心波段所包含的类，即将所有波段分成了若干组。

本发明首先利用地物在空间上分布连续性这一信息对高光谱图像分块，相比传统压缩方案，充分利用了高光谱图像的空间相关性，使得每一块高光谱图像的地物空间信息非常相似，有利于后续的压缩处理。其次，针对现有压缩算法对高光谱图像谱间相关性利用不够充分的情况，对高光谱图像分块后，再对每块图像的所有波段进行聚类分组，充分利用了高光谱图像的强谱间相关性，使得每一块中的每一组高光谱数据在空间和谱间信息上都具有很高的相似性。利用了PCA对高光谱图像进行降维处理，而PCA不仅可以用作特征的降维，还可以用来消除数据中的噪声，并且能与被处理数据的统计特性很好地匹配，是统计特征基础上的最优线性正交变换，使得经过PCA后的数据能量更为集中。总之利用PCA对高光谱图像进行处理，不仅可以在一定程度上去除高光谱图像的谱间冗余及噪声，还能降低整个压缩算法的复杂度。

本发明对高光谱图像波段聚类分组后的每个波段组进行PCA，不仅进一步利用了高光谱图像的谱间相关性，还显著降低了PCA的运算量。充分利用高光谱图像的空间相关性以及谱间相关性，对高光谱图像进行分块、波段聚类分组以及PCA等一系列处理，最后再结合JPEG2000算法对图像进行编码压缩，相比传统算法，不仅在失真性能上有所改善，而且在重构图像的异常检测以及光谱角性能方面均有所改善，使得图像中的某些细节信息得以保留，重构图像更加接近原始图像。

附图说明

图1为本发明的整个方案流程图。

具体实施方式

为了使本发明的目的、技术方案以及优点更加清楚明白，下面结合附图及实施例进一步阐明本发明的内容，但本发明的内容不仅仅局限于下面的实施例。本领域技术人员可以对本发明作各种改动或修改，这些等价形式同样在本申请所列权利要求书限定范围之内。

本发明提出的一种基于空谱特性的高光谱图像压缩方案。具体包括，对高光谱图像分块；对每一块高光谱数据的所有波段进行聚类分组；对每一块中的每一组高光谱数据做PCA运算；利用JPEG2000算法对经过PCA运算后的所有主成分进行编码压缩，得到重构图像。本发明流程图如图1所示。

本方案以典型的AVIRIS高光谱遥感图像为实例对本发明的实施作具体说明,大小为512×614×224，为处理方便，在每个波段中截取大小为256×256的图像块，最终形成大小为256×256×224的高光谱图像。

一、高光谱图像的分块

传统的高光谱图像分块只是简单地利用地物在空间上分布的连续性，通常认为如果某点处为C类地物，那么在此点的领域上也是C类地物，因此将高光谱图像立方体均分或者近似均分成若干块，但这样划分误差很大。

本方案的高光谱图像分块主要是利用了高光谱图像的空间相关性以及地物在空间上分布的连续性，用D代表高光谱数据立方体，则D_ij表示在空间坐标(i,j)处的光谱(i＝1,2,…,Μ；j＝1,2,…,N)，其中M、N分别表示高光谱图像在空间维上的行数和列数。引入一个M行、N列的分块矩阵A以及一个分块阈值μ，其中的元素A_ij表示在高光谱图像块空间坐标(i,j)处D_ij所属的块数，D_ij表示高光谱图像在空间坐标(i,j)处的像元即光谱曲线，每个像元都是该像素点对应的所有波段，分块矩阵A用来标记高光谱图像每一个像元分别属于哪一块。用F(x,y)表示相邻两条光谱曲线x和y的相似程度(欧式距离)，x和y分别表示高光谱图像在某空间坐标处的像元。若F(x,y)越小，则两条光谱曲线的相似程度越大。即当F(x,y)≤μ时，x和y相似程度比较大，则把x和y归为同一块；反之，当F(x,y)>μ时，x和y相似程度比较小，则x和y属于不同的块。

具体的处理步骤如下：

读取高光谱图像，根据空间维上像元的行数M，列数N，波段数得到一个高光谱数据立方体D，D_ij表示D在空间坐标(i,j)处的像元即光谱曲线，并对数据进行归一化处理。

对于空间维上第一行像元，从第一列至第N列开始遍历，如果相邻像元的距离小于等于阈值μ，则两像元属于同一块，如果相邻像元的距离大于μ，则两像元属于不同的块。其中μ的选取依赖具体的高光谱图像，当高光谱图像中地物类别有较明显的界限时，μ可以选得比较大；否则，μ应比较小。一般可取所有相邻像元距离的平均值。

对于空间维上非第一行像元，从第一列至第N列开始遍历，将与其相邻的像元的相似度进行比较，选出相似度最大的像元以及对应像元之间的距离，若该距离小于μ，则两像元属于同一块，否则不属于同一块。

遍历整个空间范围，建立一个M行、N列的分块矩阵A，其中的元素A_ij表示在空间坐标(i,j)处光谱D_ij所属的块数，

将矩阵A中元素数值相同的位置所对应的高光谱数据立方体D中的像元块归为同一块，即可将高光谱数据立方体D分成若干块E_n(n＝1,2,3,…)。

二、每一块高光谱图像的所有波段聚类分组

传统的波段聚类分组一般采用k-means聚类，其基本思想是依据波段相关性，从大量数据对象中随机选取k个对象作为初始聚类中心,根据最小距离准则对数据对象进行分类。对于不同的初始聚类中心,该算法容易产生不同的结果,聚类结果带有很大的随机性；此外,该算法的聚类结果不能根据高光谱图像本身的变化自动改变k值,用户需指定聚类数目、初始质心和阈值大小,其本质上不是根据聚类元素的分布进行自适应聚类,且聚类精度和效率依赖于初始质心的选取等预设参数,这不仅增加了选取k值的难度,也影响了聚类结果的准确性。

本方案将高光谱图像分成若干块后，以每一块高光谱图像E_i的所有波段作为聚类对象，以波段间的相关性作为输入数据的相似度。将所有波段之间的相关性组成的相似度矩阵S，根据公式确定聚类的能量函数，其中，C_l为波段l的聚类中心，S(l,C_l)为波段l与其聚类中心波段C_l的相似度，聚类的能量函数为各波段与其聚类中心波段的相似度之和；根据公式确定加权偏好参考度，其中，l、k为任意两波段，加权偏好参考度用来反映各波段最终被选作类代表波段的统一偏好程度，聚类中心波段与该类中其他任何波段的相似度必须大于加权偏好参考度，其数值越大，得到的类越多，反之越少。算法初始时将高光谱图像块的所有波段都视为潜在的聚类中心，通过一种迭代循环不断地更新，产生若干个聚类中心波段及其对应的聚类，同时聚类的能量函数也得到了最小化。将各波段分配到对应中心波段所包含的类，即将所有波段分成了若干组。

具体处理步骤如下：

选取第n块高光谱图像块E_n，设l和k为其中任意两波段；

根据公式s(l,k)＝R[l,k]计算波段l和波段k间的相似度s(l,k)(相关系数)，获得图像块E_n所有波段的相似度，以所有波段之间的相关系数组成相似度矩阵S，其中R[·]为谱间相关系数。

$R\[l,k\]=\frac{\underset{i=1}{\overset{M}{\Σ}}\underset{j=1}{\overset{N}{\Σ}}\left(X_k\right(i,j)-{\mathrm{\μ}}_k)\left(X_l\right(i,j)-{\mathrm{\μ}}_l)}{\sqrt{\underset{i=1}{\overset{M}{\Σ}}\underset{j=1}{\overset{N}{\Σ}}{\left(X_k\right(i,j)-{\mathrm{\μ}}_k)}^2\underset{i=1}{\overset{M}{\Σ}}\underset{j=1}{\overset{N}{\Σ}}{\left(X_l\right(i,j)-{\mathrm{\μ}}_l)}^2}}$

其中，X_k(i,j)和X_l(i,j)分别为图像块中第k波段和第l波段上空间坐标(i,j)的像素；μ_k和μ_l分别为波段k和波段l上所有像素的均值，M和N分别表示图像块空间维上的行数和列数。即， ${\mathrm{\μ}}_k=\frac{1}{MN}\underset{i=1}{\overset{M}{\Σ}}\underset{j=1}{\overset{N}{\Σ}}{X}_k(i,j)$

令r(l,k)为波段k对波段l的吸引度，用来描述波段k适合作为波段l的类代表程度；a(l,k)为波段l对波段k的归属度，用来描述波段l选择波段k作为其类代表的适合程度。初始化，对所有的l，k，r(l,k)＝0，a(l,k)＝0。

根据公式：确定加权偏好参考度λ。

根据公式： $r(l,k)\←s(l,k)-\underset{k^{\′}\≠k}{max}\{a(l,k^{\′})+s(l,k^{\′})\},$ r(l,k)←λ*r(l,k)_old+(1-λ)*r(l,k)_new计算吸引度信息与更新。其中，s(l,k)表示波段l与波段k间的相似度，s(l,k')表示波段l与波段k'间的相似度，a(l,k')表示波段l对波段k'的归属度，初始时设为0，r(l,k)表示波段k对波段l的吸引度，初始时设为0，r(l,k)_old表示前一次计算出的吸引度，r(l,k)_new代表本次计算出的波段l与波段k的吸引度。

根据公式： $a(l,k)\←\min \{0,r(k,k)+\underset{l^{\′}\∉\{l,k\}}{\Σ}\max \{0,r(l^{\′},k)\left\}\right\},$ a(l,k)←λ*a(l,k)_old+(1-λ)*a(l,k)_new计算归属度信息与更新。其中，a(l,k)表示波段l对波段k的归属度，a(l,k)_old表示前一次计算出的归属度，a(l,k)_new代表本次计算出的波段l与波段k的归属度。

不断进行迭代更新，执行样本决策：根据得到的吸引度与归属度信息，当r(l,k)+a(l,k)最大时，可得对于波段l的中心波段k，由此反推可得到若干中心波段及其对应的类，即可将所有波段分成若干组。

三、高光谱图像中的高光谱数据的PCA运算

将每块高光谱图像的所有波段分成若干组后，将每个图像块的每个组的高光谱数据重新排列为二维矩阵X，使其每一行表示一个波段上所有像素点，每一列表示一个像素点对应的所有波段。计算矩阵X的协方差矩阵C；计算协方差矩阵C的特征向量e₁,e₂,…,e_N及其对应的特征值λ₁,λ₂,…,λ_n，把特征值从大到小排序，取前r位特征值对应的特征向量组成投影矩阵W；最后投影数据到W组成的空间之中，即可得到矩阵X经过PCA后得到的主成分。其中，r表示最终所保留的主成分个数。r越小，则最终所保留数据的维数越低，便于分析，同时也降低了噪声，但信息的损失也越大；相反r越大，则最终所保留数据的维数越高，但所保存的信息也越完整。具体包括：

读取高光谱图像每一块中的每一组高光谱数据，将每个波段的数据作为矩阵的行元素，矩阵中每一列为一个L维矢量，将高光谱数据重新排列为矩阵X，

$X=\left[\begin{array}{cccc}{X}_1(1,1)& X_1(1,2)& ...& X_1(M,N)\\ X_2(1,1)& X_2(1,2)& ...& X_2(M,N)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ X_L(1,1)& X_L(1,2)& ...& X_L(M,N)\end{array}\right]$

X中的每一行表示一个波段的数据，而X中的每一列为一个L维矢量，即：X_i＝[x_1,i,x_2,i,…,x_L,i]^T，i＝1,2,…,M×N，X_i表示矩阵中的某一行，x_L,i表示矩阵中具体某一行某一列的元素。显然，X＝[X₁,X₂,…,X_MN]^T，M和N分别表示某一块中的某一组高光谱数据空间维上的行数和列数。

根据公式 $m_r=\frac{1}{MN}\underset{j=1}{\overset{MN}{\Σ}}{X}_{r,j},r=1,2...,L$ 计算每个波段的样本均值，假设均值向量为m_X＝[m₁,m₂,…m_L]。m_X是由每个波段的样本均值m_r组成的均值向量。

计算高光谱数据矩阵X的协方差矩阵。根据公式C_l,i＝E[(X_i-m_r)(X_i-m_r)^T]，i＝1,2,…,M×N计算任意光谱的协方差矩阵，其中X_i表示高光谱数据矩阵中一个波段的数据，m_r表示对应波段的样本均值，E表示求方括号中向量的期望。所有光谱矢量的协方差矩阵为

计算协方差矩阵的特征值和特征向量，其中C_l为L×L的对称矩阵，存在正交矩阵U＝[u₁,u₂,…,u_L]，使得C_l可对角化。

其中，对角矩阵的对角元素{λ₁,λ₂,…,λ_L}为C_l的特征值，且λ₁≥λ₂≥…≥λ_L≥0，正交矩阵中u_i为对称矩阵C_l的第i个特征值对应的特征向量(i＝1,2,…,L)。

对高光谱数据矩阵X进行降维。根据经验设定阈值，选取大于阈值的特征值所对应的特征向量组成特征矩阵，用X乘以特征矩阵即可得到经过PCA运算后的主成分。

四、利用JPEG2000进行最终的编码压缩

传统的编码压缩部分可以采用熵编码、JPEG等方式，但JPEG2000是最新一代的静态图像压缩标准，其编码压缩效果非常好。

得到每一块中的每一组高光谱数据经过PCA运算后的主成分后，再利用JPEG2000压缩标准对所有主成分进行编码压缩得到重构图像。

本发明的优点由以下仿真实例和所得实验结果进行具体说明。

实验中使用美国喷气推进实验室提供的免费高光谱数据“LowAltitude”和“LunarLake”作为实验数据，大小均为512×614×224，为处理方便，在每个波段中截取大小为256×256的图像块，最终形成大小为256×256×224的高光谱图像，均在MATLAB2013a软件平台下进行仿真，这些源数据被绝大多数研究者使用，实验结果具有代表性。

高光谱遥感图像的压缩性能主要用图像的信噪比(SNR)、图像的压缩比(CR)、重构图像的异常检测、重构图像的光谱角以及算法的计算复杂度来进行评价。在同一压缩比下，图像的信噪比越高，重构图像的异常检测性能越好，以及重构图像的光谱角越小，则重构图像的质量越好，重构图像越接近于原始图像，故该算法的压缩性能越好。但图像的信噪比与重构图像的异常检测以及光谱角性能并不是成正比，即图像的信噪比越高，重构图像的异常检测以及光谱角性能不一定更好。因此，针对某一方案压缩性能的评价要从多方面综合考虑。

1、SNR计算公式

$SNR=10\×\lg \frac{mean\left(I^2\right)}{MSE+\frac{1}{12}}$

I表示高光谱原始图像，I表示重构图像，I(x,y,λ)表示一个像素的像素值，其中x，y，λ分别表示行、列和波段序号。n_x，n_y，n_λ分别表示行列和波段数。mean(I²)表示高光谱图像所有像素点的均方值，MSE表示原始高光谱图像与重构图像的均方误差。

$MSE=\frac{1}{n_x{n}_y{n}_{\mathrm{\λ}}}\underset{x,y,\mathrm{\λ}}{\Σ}{\left(e_I\right(x,y,\mathrm{\λ}\left)\right)}^2$

其中，原始数据与重构数据的误差表示为

2、CR计算公式

$CR=\frac{B}{B^{*}}$

其中，B表示原始高光谱数据的大小，B^*表示重构高光谱数据的大小。

3、衡量重构图像异常检测性能计算公式

$\mathrm{\ρ}(l,k)=\frac{\underset{x=1}{\overset{M}{\Σ}}\underset{y=1}{\overset{N}{\Σ}}\[f(x,y)-{\mathrm{\μ}}_f\]\[f(x+k,y+l)-{\mathrm{\μ}}_f\]}{\underset{x=1}{\overset{M}{\Σ}}\underset{y=1}{\overset{N}{\Σ}}{\[f(x,y)-{\mathrm{\μ}}_f\]}^2}$

其中，用空间相关系数ρ来反映某一波段图像异常检测性能，其值越接近于1.0，表明异常检测的性能越好，即图像中的一些小的目标或物质在压缩过程中保留了下来，也即图像在压缩中的损失越小。(l,k)为像素间距，f(x,y)为某一波段图像上空间坐标为(x,y)灰度函数，μ_f为某一波段图像灰度平均值，M×N为图像空间尺寸，最终求得重构图像所有波段的空间相关系数ρ的均值，即可衡量重构图像异常检测性能。

4、平均光谱角计算公式

${\mathrm{\μ}}_{SAM}=arccos\frac{T\×R}{\left|T\right|\×\left|R\right|}$

其中，T和R分别为原始高光谱图像与重构图像中同一像元所对应的光谱矢量。根据所有光谱矢量求得的平均光谱角越小，表明重构图像与原始高光谱图像光谱之间的相似程度越高，即重构图像与原始图像越接近，也即重构图像的质量越好。

通过上述压缩性能的评价方法对DWT+JPEG2000算法、PCA+JPEG2000算法、波段上分组+PCA+JPEG2000算法与本发明方案进行比较，表1列出了四种方法在相同压缩比CR的条件下，信噪比SNR、异常检测性能ρ以及平均光谱角μ_SAM的比较。

表1.五种算法在信噪比、异常检测以及光谱角性能方面的比较

从实验结果可以看出：对于高光谱图像LowAltitude，本发明所提方案在相同压缩比的情况下，与DWT-JPEG2000算法相比，信噪比平均提高了7.3dB，异常检测性能ρ平均提高了0.198，平均光谱角减少了0.654°；与PCA-JPEG2000算法相比，信噪比平均提高了1.4dB，异常检测性能ρ平均提高了0.162，平均光谱角减少了0.441°；与Group-PCA-JPEG2000算法相比，信噪比基本没有提高，但异常检测性能ρ平均提高了0.036，平均光谱角减少了0.190°。而对于高光谱图像LunarLake，本发明所提方案在相同压缩比的情况下，与DWT-JPEG2000算法相比，信噪比平均提高了6.4dB，异常检测性能ρ平均提高了0.209，平均光谱角减少了0.375°；与PCA-JPEG2000算法相比，信噪比平均提高了0.6dB，异常检测性能ρ平均提高了0.096，平均光谱角减少了0.172°；与Group-PCA-JPEG2000算法相比，信噪比基本没有提高，但异常检测性能ρ平均提高了0.028，平均光谱角减少了0.116°。

综上所述，本方案以及Group-PCA-JPEG2000算法与其他算法相比，无论是信噪比还是异常检测和光谱角性能均有较大改善，而本方案与Group-PCA-JPEG2000算法相比，虽然在信噪比方面没有明显提升，但在异常检测性能以及光谱角方面均有一定改善，说明本方案更加充分地利用了高光谱图像的空谱特性，有效地保留了图像中一些细节信息，改善了压缩性能，使得重构图像的质量进一步提升。随着高光谱技术的发展和广泛应用，高光谱遥感图像的波段数不断增加，成像光谱仪获取的数据量越来越大，这给数据的存储和传输带来了巨大的压力，对高光谱数据进行有效压缩显得非常迫切，因此本发明具有很高的应用价值。

Claims (7)

1.一种基于空谱特性的高光谱图像压缩方法，其特征在于，包括步骤：利用高光谱图像的空间相关性以及地物在空间上分布的连续性对高光谱图像分块；对每个高光谱图像块进行自适应波段聚类分组；对每个高光谱图像块的每个波段组分别进行PCA运算，得到经过若干降维后的高光谱图像主成分，对所有高光谱图像主成分进行编解码，得到重构图像。

2.根据权利要求1所述的方法，其特征在于，所述对高光谱图像分块进一步包括：读取高光谱图像，根据空间维上像元的行数M、列数N、波段数得到一个高光谱数据立方体D；用像元间的欧式距离衡量像元间的相似度，遍历空间维上第一行像元，如果相邻像元间的欧式距离小于等于阈值μ，则两像元属于同一块，如果相邻像元间的欧式距离大于μ，则两像元属于不同的块；遍历空间维上非第一行像元，比较相邻像元的相似度，选择相似度最大的像元以及对应像元之间的欧式距离，若该欧式距离小于μ，则两像元属于同一块。

3.根据权利要求1所述的方法，其特征在于，所述对高光谱图像块进行自适应波段聚类分组进一步包括，根据所有波段之间的相似度组成相似度矩阵，根据公式确定聚类的能量函数，其中，C_l为波段l的聚类中心，S(l,C_l)为波段l与其聚类中心波段C_l的相似度；根据任意两波段l、k的相似度调用公式确定加权偏好参考度p。

4.根据权利要求1所述的方法，其特征在于，读取每个高光谱图像块中的高光谱数据，将高光谱数据重新排列为矩阵X，矩阵X中的每一行表示一个波段的数据，矩阵X中的每一列为一个L维矢量；根据公式计算每个波段的样本均值，根据公式C_l,i＝E[(X_i-m_X)(X_i-m_X)^T]计算任意光谱的协方差矩阵，根据公式获得所有光谱矢量的协方差矩阵；计算协方差矩阵的特征值和特征向量，对矩阵X进行降维，选取大于阈值的特征值所对应的特征向量组成特征矩阵，用矩阵X乘以特征矩阵得到经过PCA运算后的高光谱图像主成分。

5.根据权利要求2所述的方法，其特征在于，遍历整个空间范围的高光谱图像，建立一个M行、N列的分块矩阵A，将矩阵A中元素数值相同的位置所对应的高光谱数据立方体D中的像元块归为同一块，将高光谱数据立方体D分成若干块，其中，矩阵A的元素A_ij表示在空间坐标(i,j)处光谱D_ij所属的块数。

6.根据权利要求3所述的方法，其特征在于，确定波段间的相关系数为波段间的相似度，以任一图像块E_n所有波段之间的相似度组成相似度矩阵S，根据公式 $R\[l,k\]=\frac{\underset{i=1}{\overset{M}{\Σ}}\underset{j=1}{\overset{N}{\Σ}}\left(X_k\right(i,j)-{\mathrm{\μ}}_k)\left(X_l\right(i,j)-{\mathrm{\μ}}_l)}{\sqrt{\underset{i=1}{\overset{M}{\Σ}}\underset{j=1}{\overset{N}{\Σ}}{\left(X_k\right(i,j)-{\mathrm{\μ}}_k)}^2\underset{i=1}{\overset{M}{\Σ}}\underset{j=1}{\overset{N}{\Σ}}{\left(X_l\right(i,j)-{\mathrm{\μ}}_l)}^2}}$ 获得波段l和波段k之间的谱间相关系数，其中，X_k(i,j)和X_l(i,j)分别为图像块中第k波段和第l波段上空间坐标(i,j)的像素，μ_k和μ_l分别为波段k和波段l上所有像素的均值，M和N分别表示图像块空间维上的行数和列数；根据公式：确定加权偏好参考度λ，根据公式： $r\left(l,k\right)\←s\left(l,k\right)-\underset{k^{\′}\≠k}{max}\left\{a\left(l,k^{\′}\right)+s\left(l,k^{\′}\right)\right\},$ r(l,k)←λ*r(l,k)_old+(1-λ)*r(l,k)_new计算波段间的吸引度并更新；根据公式：a(l,k)←λ*a(l,k)_old+(1-λ)*a(l,k)_new计算波段间的归属度并更新，其中，s(l,k)表示波段l与波段k间的相似度，s(l,k')表示波段l与波段k'间的相似度，a(l,k')表示波段l对波段k'的归属度，r(l,k)表示波段k对波段l的吸引度，r(l,k)_old表示波段k对波段l前一次的吸引度，r(l,k)_new代表波段l与波段k当前的吸引度，a(l,k)_old表示波段k对波段l前一次的归属度，a(l,k)_new代表波段k对波段l当前的归属度，不断进行迭代更新，根据当前得到的吸引度与归属度，当r(l,k)+a(l,k)最大时，获得波段l的中心波段k，根据中心波段将所有波段分成若干组。

7.根据权利要求4所述的方法，其特征在于，存在正交矩阵U＝[u₁,u₂,…,u_L]，使得所有光谱矢量的协方差矩阵C_l可对角化，即对角矩阵的对角元素{λ₁,λ₂,…,λ_L}为矩阵C_l的特征值，u_i为矩阵C_l的第i个特征值对应的特征向量。