快速压缩高光谱信号的矢量量化方案的制作方法

本发明属于高光谱图像信号处理领域，具体涉及一种快速压缩高光谱信号的矢量量化方案。

背景技术：

高光谱遥感技术是20世纪80年代兴起的新型对地观测技术。这一技术将确定物质或地物性质的光谱特性与把握其空间和几何关系的特性革命性地结合在一起。成像光谱仪在空间成像的同时，以同一谱间分辨率记录下几十到几百组光谱数据，将它们叠合在一起，就可以得到一个高光谱图像立方体。从高光谱图像立方体的每个像元均可提取一条连续的光谱曲线。光谱特性为其获取的遥感数据提供了丰富的地物细节，在海洋遥感、植被研究、精细农业、地质勘察、大气和环境遥感以及军事侦察等方面获得了广泛的应用，并显现出巨大潜力。

随着成像光谱技术的不断发展，人们对遥感对象的了解不断深入。与此同时，新的问题产生，给高光谱数据应用带来巨大挑战，即海量数据的存储与传输问题。高光谱数据的光谱分辨力和空间分辨力越来越高，数据量化深度也不断增加，同时卫星不间断周期性扫描次数也在不断增加。这些因素都导致了高光谱遥感图像数据量的急剧增长。如果每个像素灰度值用2字节存储，那么一幅大小为614×512×224的高光谱图像，大约需要140M字节。对于当今越来越紧张的信道资源来说，为了进行有效的传输，研究出快速有效的高光谱图像压缩方法是十分重要的。

高光谱图像压缩的可行性在于图像数据的高相关性。数字图像压缩主要就是去除数据间存在的相关性或冗余度，二维静止图像的压缩主要是去除数据的空间冗余，而高光谱图像立方体相当于一个三维图像，在具有空间相关性的同时，还具有谱间相关性，高光谱图像的压缩主要在于去除这两类相关性。

矢量量化技术压缩比高、编解码简单，是一种有效的快速数据压缩手段。矢量量化的原理是将若干个标量数据组构成一个矢量，然后在矢量空间给以整体量化，从而压缩了数据而不损失多少信息，能在高压缩率和平均最小失真间获得最佳折衷。

矢量量化由两大关键技术组成：码书设计和码字搜索。码书设计指寻找最优码书，使恢复图像与原始图像之间的失真达到最小，以保证重构图像较好的质量；码字搜索指快速地找到与输入矢量误差最小的码字。矢量量化的主要问题是其较高的运算复杂度，尤其是码书训练的过程，其计算量随着矢量维度的增加呈指数级增长。因此，需要寻找快速的码书搜索算法，在保证码书质量几乎不受影响的前提下，尽可能的降低运算量

矢量量化中码书设计的经典算法是1980年由Linde、Buzo、Gray提出的LBG算法。它不需要事先知道输入矢量的概率分布，通过训练矢量集和一定的迭代算法来逼近最优码书。LBG算法的思想是：①随意选取N个训练矢量作为初始码书；②由这N个码矢量对所有的训练矢量进行划分，即分成N个集合，使每个集合中的矢量，都是与码书中各码字矢量欧氏距离中最小的；③由这N个集合的质心，得到新的码书；④如果这些码矢量与原来的码矢量变化小于某个阈值，就完成码书的训练，否则重新进行②、③步。LBG算法的思想比较简单，但存在两大缺点：一是对初始码书非常敏感，初始码书的选取，不但影响递推的收敛速度，而且影响码书的性能；二是训练时间长，不能自适应的跟踪图像信源的统计特性。

技术实现要素：

本发明针对现有的高光谱图像压缩方案，在获得较高的图像恢复质量的同时，计算量较大的问题，提出了一种快速压缩高光谱图像的矢量量化方案。该方案能在保证图像恢复质量较降低运算复杂度之前的恢复质量相比几乎不变的情况下，大大降低运算复杂度，大大减少压缩所需时间。具体为：

种快速压缩高光谱图像的矢量量化方法，读取高光谱图像数据构建初始码书Y，将初始码书做PCA运算，获得PCA特征矩阵U，从U的分量Ur中提取初始码书的主成分得到初始码书主成分YU；在训练矢量中排除异常像素，从排除异常像素的训练矢量中提取训练矢量主成分BU；记录训练矢量主成分对应的关系索引IB；根据初始码书主成分YU对训练矢量主成分BU进行聚类，所有矢量聚类后根据索引IB找到对应训练矢量，利用每个胞腔中找到的训练矢量的质心更新码书，获得最终训练矢量码书，将最终训练矢量码书和训练矢量码字索引进行存储并传输。构建初始码书Y具体为：读取三维高光谱图像数据，转化为二维矢量数据，将二维矢量数据的所有维度相加得到和值，根据和值将二维矢量数据排序，根据码书大小将排序后的矢量平均分成多份，取每份的质心作为初始码书。

本发明进一步具体包括，从U的分量Ur中提取初始码书的初始码书的主成分具体包括：将同一波段的数据作为矩阵的相同列元素，矩阵中每一行为一个k维码子；根据码书Y的元素Yjr,j＝1,2,…,N,r＝1,2,…,k调用公式计算每个波段的样本均值，其中N为码书大小；根据各波段的样本均值和高光谱数据矩阵中各个波段的数据计算任意两个波段的协方差组成协方差矩阵C，根据协方差矩阵C的特征向量大小对C的特征向量排序得PCA特征矩阵，取特征矩阵U中预定分量与初始码书相乘得到初始码书主成分YU。

本发明进一步具体包括，利用改进的RXD算法对训练矢量进行异常检测，记录异常像素的矢量值BRXD1,BRXD2,…,BRXDM和位置索引IRXD，从训练矢量中排除异常像素，对异常像素的矢量不压缩，直接传输其矢量值。

所述将初始码书Y做PCA运算具体包括：将Y的每一行作为一个码子，Y的每一列为一个波段的数据，得到一个由N个k维行矢量组成的向量组，即根据公式r＝1,2…,k计算每个波段的样本均值，根据公式Cij＝E[(Yi-mi)(Yj-mj)^T]，i,j＝1,2,…,k计算任意两个波段的协方差Cij，根据协方差构建协方差矩阵：寻找使得C可对角化的正交矩阵U作为PCA特征矩阵，其中，Yi，Yj表示高光谱数据矩阵中i,j波段的数据，mi(mj)表示对应i(j)波段的样本均值，E[·]表示求方括号中向量的数学期望。

本发明进一步具体包括，将特征矩阵中预定分量Ur与训练矢量B1相乘得到训练矢量矩阵主成分BU，获取训练矢量和训练矢量主成分对应关系索引IB，依次将矩阵BU中的矢量利用高效排除不等式算法在初始码书主成分YU中快速搜索最佳匹配码字，依次将BU中的训练矢量划分到最佳匹配码字所对应的胞腔中，记录对应的码字索引，根据索引IB找到对应的原训练矢量。

本发明进一步具体包括，根据第i个训练矢量的和值αi计算和值的累积分布函数Pr(a)；给定可靠系数γ，找到满足Pr(a)≤γ的最大和值α0；找到最大和值α0对应的训练矢量r0，根据公式：δ＝r0^TU^-1r0计算阈值，如训练矢量大于阈值，则对应像素为异常像素。

本发明实现过程可包括以下几个阶段：

初始化阶段：读取高光谱图像三维数据，将同一位置的所有波段的像素按序排列成一个矢量，即把高光谱图像的每一个像元对应的光谱矢量作为一个矢量数据。根据矢量和值(一范数)将所有矢量按序排列，根据初始码书大小，将排序后的矢量平均分成若干类，取每一类的质心组成初始码书。

数据预处理阶段：对初始码书做主成分PCA分析(Principle Component Analysis)，得到PCA特征矩阵，通过初始码书和PCA特征矩阵获得码书主成分。具体可采用方法包括：取PCA特征矩阵中的部分分量(如可取前5％)，与初始码书相乘，得到码书主成分。

异常检测阶段：对PCA特征矩阵进行异常像素检测RXD(Reed-Xiaoli Detection)，对矩阵中的每个矢量进行异常检测。记录异常像素值和对应的位置索引，并从训练矢量中排除。对异常像素不做矢量量化，按照其原有像素值进行存储和传输。

训练阶段：根据训练矢量和特征矩阵获取训练矢量主成分。将训练矢量与PCA特征矩阵中的预定部分分量(最优可取前5％的分量)相乘得到训练矢量主成分，记录训练矢量主成分与原训练矢量对应关系的索引。利用高效排除不等式算法快速搜索训练矢量主成分对应的最佳匹配主成分码字，并将训练矢量主成分划分到对应的胞腔中，记录对应主成分码字的索引，根据关系索引找到训练矢量主成分对应的原训练矢量，计算原训练矢量胞腔的质心，以各胞腔的质心代替包腔对应码字更新码书。完成所有码书的更新。

具体可采用如下步骤：

(1)读取高光谱图像数据，根据图像尺寸和光谱波段数构建三维数据矩阵，以同一空间位置对应的各波段像素值作为一个矢量数据，即把高光谱图像每一个像元对应的光谱矢量对应一个矢量数据，把三维数据转化为二维矩阵数据B，B的每一行代表一个矢量数据，B的行数为矢量个数，B的列数为光谱波段数；

(2)计算B的每一行矢量的和值(一范数)，将每一行矢量根据和值排序，记录排序索引I。根据I将排序后的B根据码书大小N平均分为N组，取每组的质心组成初始码书Y。

(3)将初始码书Y做PCA运算得到PCA特征矩阵U，取特征矩阵的部分分量Ur＝{u1,u2,…,ur}对初始码书进行主成分提取,其中ui,i＝1,…,r为特征矩阵的第i个特征向量，得到初始码书主成分YU。

(4)预设一异常像素判决阈值δ，对PCA特征矩阵U进行RXD异常检测，具体可为：对二维数据矩阵B中每一行矢量做RXD异常检测，若检测结果大于δ，则该像素就被判为异常像素。记录异常像素的矢量值BRXD1,BRXD2,…,BRXDM和其位置关系索引IRXD，并在B中排除上述矢量值获得排除异常像素后的训练矢量矩阵B1。异常像素不进行压缩，按照本身的矢量值进行存储和传输。

(5)取PCA特征矩阵U的部分分量Ur与排除异常像素后的训练矢量矩阵B1中所有非异常像素做PCA运算，得到训练矢量主成分矩阵BU，记录训练矢量和主成分其对应位置关系索引IB。当达到设置的迭代次数ite或失真阈值e时结束。

(6)根据B1、BU、IB、YU训练获得新的码书，具体可采用如下方法：依次将训练矢量主成分矩阵BU中的矢量在初始码书主成分YU中快速搜索最佳匹配码字，依次将BU中的训练矢量划分到其最佳匹配码字所对应的胞腔中，根据索引IB在训练完成的胞腔中找到对应的原训练矢量，根据原训练矢量的质心获得新的码字，构成新的码书。

(7)判断迭代次数是否已经达到预先设定的值ite，或者判断本次迭代产生的平均失真与前一次迭代产生的平均失真之间的相对误差是否小于预设的失真阈值e。如果满足上述迭代结束条件，则终止训练，并保存最后一次迭代产生的码书和训练矢量码字索引作为最终压缩结果，进行存储和传输。

附图说明

图1为传统的LBG算法流程图；

图2为本发明矢量量化方法流程图。

具体实施方式

以下结合附图和具体实例对本发明作进一步说明。本发明提出的快速压缩高光谱图像的矢量量化方案，主要包括以下几个阶段：

在量化处理之前可对矢量数据初始化获得初始码书。

初始化阶段：读取高光谱图像三维数据，以同一空间位置各波段像素值作为一个矢量数据，即把高光谱图像每一个像元对应的光谱矢量对应一个矢量数据，将三维的原始数据转化为二维的矢量矩阵数据。求矢量矩阵中所有矢量的和值(一范数)，将矢量根据和值排序，将排序后的矢量平均分为N组，取每组的质心构成初始码书，其中，N为码书大小。

数据预处理阶段：对初始码书做主成分分析PCA(Principle Component Analysis)，得到PCA特征矩阵，将初始码书的每一维看作一个随机变量，获取初始码书的协方差矩阵，再对协方差矩阵做特征分解，将分解得到的特征向量按照特征值从大到小排序。取特征矩阵中部分特征向量(如可取前5％的分量)与初始码书相乘得到码书主成分。

异常检测阶段：将预处理中得到的PCA特征矩阵进行异常检测(可采用异常检测RXD(Reed-Xiaoli Detection)算法对矩阵中每个像素进行检查获得异常检测值)，预设一异常像素判决阈值，如果异常检测值大于阈值，则该像素就被判为异常像素。异常像素不进行压缩，按照实际的矢量值完全传输。记录异常像素的矢量值和位置索引，并从矢量量化训练矢量中排除异常像素矢量。

训练阶段：将排除了异常像素的训练矢量与PCA特征矩阵中的预定分类相乘(可取特征矩阵前5％分量)，得到训练矢量主成分，记录训练矢量主成分与原训练矢量对应关系索引。快速搜索训练矢量主成分对应的最佳匹配主成分码字(可利用高效排除不等式算法或本领域技术人员了解的方法)，并将训练矢量划分到对应的胞腔中。直到所有训练矢量都已经迭代完成时，根据索引找到对应的原训练矢量。根据原训练矢量的质心代替原来胞腔对应码字，得到新的码书。

中断检查阶段：如果迭代次数达到预先设定次数ite，或者本次迭代产生的平均失真与前一次迭代产生的平均失真之间的相对误差小于设定的失真阈值e，则达到迭代结束的条件，若满足迭代结束条件，则终止训练，并保存最后一次迭代产生的码书和训练矢量码字索引作为最终压缩结果；否则，将码书重新进行主成分提取，求得新的主成分码书，进行新一轮训练矢量的训练。

以下具体描述本发明的详细实现方式：

(1)构造矢量数据

读取高光谱图像，截取每个波段中的同一空间位置的一块图像块作为输入信号，根据截取的图像块尺寸和波段数构建三维高光谱图像矩阵X，X的大小为D1×D2×D3，其中D1×D2为单个波段的二维图像，D3为波段数。以同一空间位置各波段像素值作为一个矢量数据，即把截取的图像块每一个像素对应的光谱矢量作为一个矢量数据，把三维数据转化为二维D×D3矩阵数据B，D＝D1×D2，记录矢量在原图像中的位置索引IX。

(2)构造初始码书

将每个训练矢量的每一维相加，得到训练矢量的和值(一范数)，根据和值的大小将训练矢量排序。根据初始码书的大小N将排序后的训练矢量B平均分为N份B1,B2,…,BN，取每一份的质心作为初始码书Y。Y＝{Y1^T,Y2^T,…,YN^T}，Xi∈Bi,i＝1,2,…,N。size(·)表示每份的训练矢量个数。

(3)求初始码书的主成分

将初始码书Y做PCA运算求主成分，具体可采用如下步骤：

将Y的每一行构成一个码子，得到一个由N个k维行矢量组成的向量组。Y的每一列为一个波段的数据。

根据公式r＝1,2…,k计算每个波段的样本均值，设mX＝[m1,m2,…mk]，mX是由每个波段的样本均值mr组成的均值向量。

计算初始码书矩阵Y的协方差矩阵。根据公式Cij＝E[(Yi-mi)(Yj-mj)^T]，i,j＝1,2,…,k计算任意两个波段的协方差Cij，其中，Yi表示高光谱数据矩阵中一个波段的数据，mi(mj)表示对应波段的样本均值，E[·]表示求方括号中向量的数学期望。得到协方差矩阵

计算协方差矩阵C的特征值和特征向量，其中C为k×k的对称矩阵。存在正交矩阵U＝[u1,u2,…,uk]，使得C可对角化。

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

取正交矩阵U＝[u1,u2,…,uk]中的预定部分分量Ur＝[u1,u2,…,ur]，将初始码书Y乘以Ur即可得到初始码书主成分YU。

(4)异常像素的提取

异常像素是指与周围像素特性相差较大且无法计算先验概率(无法事先判断像素出现的可能性)的像素，这一类像素出现概率小，信号能量低。与整幅图像光谱矢量特性偏差较大。这类像素会严重增加矢量量化的量化噪声。因此，在矢量量化前进行异常像素提取是必要的。提取的异常像素不做量化，按其本身的矢量值传输。

异常像素提取通常可采用RXD算法，通过计算矢量的马氏距离是否大于某一阈值来判断是否为异常像素，令r为一个高光谱图像的L维矢量，则RXD公式计算矢量r的阈值为：

δRXD(r)＝(r-μ)^TKL×L^-1(r-μ)

其中，μ为所有像素的均值，KL×L为图像的协方差矩阵。

本发明可采用如下改进的RXD算法，考虑高光谱图像中一阶统计量对异常检测的影响，同时考虑异常像素在整个图像中所占的能量非常少，因此，根据公式计算矢量r的阈值，其中，r^T为r的转置，对每一个训练矢量进行异常像素检测，如训练矢量大于阈值，则对应像素为异常像素。

阈值δ可采用以下方法确定：一般寻找的是在像素相关度较大的区域内灰度值突然增大的像素，灰度值越高的像素，越有可能是异常像素。将所有像素的灰度值做统计，得到灰度值的累积分布函数。设一可靠系数γ，求得累积分布函数等于γ的灰度值对应的像素。

具体步骤如下：根据每个训练矢量的和值αi计算和值的累积分布函数Pr(a)；给定可靠系数γ，找到满足Pr(a)≤γ的最大值α0；找到和值α0对应的训练矢量r0，根据公式：δ＝r0^TU^-1r0计算阈值。检测完成后，记录异常像素的矢量值BRXD1,BRXD2,…,BRXDM，在矢量在原图像中的位置索引IX中找到异常像素的位置组成位置索引IRXD并从训练矢量中排除，剩余的非异常像素记作B1。

(5)非异常像素的主成分提取

提取PCA特征矩阵中部分预定分量Ur，与所有非异常像素集合B1相乘。得到训练矢量主成分矩阵BU，记录训练矢量主成分与原训练矢量的对应关系索引IB。

(6)非异常像素的聚类

根据非异常像素矩阵B1、训练矢量主成分矩阵BU、主成分初始码书YU、记录训练矢量主成分与原训练矢量对应关系索引IB进行非异常像素的聚类，具体包括：依次将训练矢量主成分矩阵BU中的矢量利用高效排除不等式算法在初始码书主成分YU中快速搜索最佳匹配码字，依次将BU中的训练矢量划分到最佳匹配码字所对应的胞腔中，记录对应的码字索引，直到BU中所有的行矢量都训练完成，根据索引IB找到对应的原训练矢量，求原训练矢量的质心代替原码字，得到新的码书。

迭代结束的条件设置为：迭代次数达到预设值ite，或本次迭代的平均失真与前一次迭代的平均失真之间的相对误差小于失真阈值e。

判断迭代次数是否达到预设值ite，或者判断本次迭代平均失真与前一次迭代平均失真之间的相对误差是否小于预设失真阈值e，若满足迭代结束的条件，则终止训练，并保存最后一次迭代产生的各部分训练矢量码书和训练矢量码字索引作为最终压缩结果，进行存储和传输。否则，重新获取新的码书主成分YU进行训练矢量的聚类。

非异常像素的聚类和中断检查具体可采用如下方法：

定义B1＝{b0i,i＝1,2,…n}为原训练矢量集合，BU＝{bi,i＝1,2,...n}为主成分训练矢量集合，IB为BU和B1对应关系索引，YU＝{yi,i＝1,2,...,N}为初始码书主成分，P(yi)表示以yi为质心的胞腔，设定t的上限ite和失真阈值e。

初始化t＝0,D-1＝∞；

根据最小误差准则，将训练矢量主成分集合中的矢量划分给各由P(yi)表示的胞腔。即如果对于所有k＝1,2,…,N均有||bi-yj||≤||bi-yk||，则bi∈P(yj)，此时di＝||bi-yj||；

根据公式计算总的误差；

根据每个胞腔中的训练矢量主成分和索引IB找到原训练矢量b0k∈B1，得到原训练胞腔P0(yi)，i＝1,2,…,N，计算原训练矢量胞腔的质心Y＝{y0j,j＝1,2,…,N}，其中，|P0(yj)|表示胞腔P0(yj)中元素的数目。

如果t>ite，或相对误差则算法终止。否则，重新获取新的主成分码书YU，置t＝t+1，重新将矢量划分给相应的胞腔。

本发明在LBG算法基础上。将初始码书根据和值排序后再平均分成若干份求得，提升了初始码书的质量；将异常检测技术引入矢量量化中，压缩比略有降低的情况下，提高了图像恢复质量；将PCA引入LBG算法，让聚类过程可以在极低的维度下进行，在基本不影响图像恢复质量的前提下，大幅度降低了运算复杂度。本文提出的压缩方案实现过程的空间复杂度远远小于LBG算法，即使在处理器内存受限时也可以高效完成高光谱图像的较高质量的压缩。所以本文提出的压缩方案能在保证较好的图像恢复质量的前提下，高效快速地完成高光谱图像的压缩，而且操作简单，具有较高的应用价值。高光谱图像信息丰富、应用广泛，但具有波段数多和数据量庞大的特点，因此高光谱图像的压缩方法研究一直备受关注。本发明提出的压缩方案具有较高压缩比、较好图像恢复质量、较低计算复杂度的特点，可以作为高光谱图像的一种快速压缩方案。