CN104244018B - 快速压缩高光谱信号的矢量量化方法 - Google Patents

Abstract

本发明针对高光谱图像数据量巨大，提出了一种快速压缩高光谱信号的矢量量化方案。对矢量数据的各分量按离散度进行重新排序，提取矢量数据的前几维分量组成训练矢量集，余下分量组成余下矢量数据，对训练矢量集进行聚类训练；再利用训练矢量数据的编码索引对余下矢量数据直接进行归类；将训练矢量数据的码书和余下矢量数据的码书重组，得到矢量数据的最终码书，完成高光谱图像的压缩编码。本发明使矢量的大部分能量集中在低维部分，仅仅提取矢量数据低维部分进行训练，空间复杂度得到很大减小。不仅能保证图像恢复质量基本不变、大量缩小压缩过程的空间复杂度，还能达到大幅度降低各类计算量，快速高效完成高光谱图像压缩的目标。

Description

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

技术领域

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

背景技术

随着光学技术与传感技术的发展，在二十世纪八十年代，高光谱遥感技术应运而生。这一技术利用具有高光谱分辨率的成像光谱仪进行对地观测，产生一组将光谱信息和图像信息合二为一的三维高光谱图像。高光谱图像的每个光谱波段对应观测场景一幅完整的二维图像，用于把握物质空间结构，而场景中的每个像素点可以提取出一条连续的光谱曲线，用于确定地物性质。所以利用高光谱图像能精确地实现地物分类、目标检测和异常检测，在地质、环境、水文、大气、农业、军事等方面具有重要的应用价值。

但是在信息量丰富的同时，高光谱图像也具有数据量庞大的特点。高光谱图像的波段数往往达到数百个，光谱分辨率达到纳米级，如果每个像素灰度值用2字节存储，那么一幅大小为614×512×224的高光谱图像，大约需要140Mbit的存储空间。尤其，随着成像光谱技术的不断发展，光谱分辨率还在不断提高，高光谱图像的数据量还将继续增加，因此，研究出快速有效的高光谱图像压缩方法显得十分迫切。

高光谱图像的特征最为突出的是图像的相关性。高光谱图像的相关性包括两个方面：空间相关性和谱间相关性。空间相关性指每个光谱波段内某一像素与其相邻像素之间的相似性；谱间相关性是指图像同一空间位置像素与相邻波段像素有相似性。高光谱图像涉及的目标大，空间相关性较普通二维图像低，而高光谱图像不同波段涉及的地面目标相同，谱间相关性高，而且谱间相关性强于空间相关性。

矢量量化技术压缩比高、编解码简单，是一种快速有效的数据压缩手段。矢量量化的原理是把数据划分为数据块，直接对数据块进行量化，而不需要去相关处理。基于矢量量化的压缩方案理论上是以信息的高阶熵为下限，能在高压缩率和平均最小失真间获得最佳折衷。

矢量量化编码由两大关键技术组成：码书设计和码字搜索。码书设计指寻找最优码书，使恢复图像与原始图像之间的失真达到最小，以保证重构图像获得较好的质量；码字搜索指快速地找到与输入矢量最匹配的码字。矢量量化的主要问题是其较高的压缩复杂性度，尤其是码书产生的过程，其计算量随着矢量维数的增加呈指数增长。所以，矢量快速搜索算法这一领域一直是比较活跃的，不断有新的算法出现。

矢量量化中码书设计的经典算法是1980年由Linde、Buzo、Gray三人提出的LBG(Linde Bazo Gray)算法。它不需要事先知道输入矢量的概率分布，通过训练矢量集和一定的迭代算法来逼近最优的再生码本。LBG算法的思想是：①随意选取N个训练矢量作为码矢量；②由这N个码矢量对所有的训练矢量进行划分，即分成N个集合，使每个集合中的矢量，都是与各码矢量距离中，与对应的码矢量的距离最小的；③由这N个集合的重心，得到N个新的码矢量；④如果这些码矢量与原来的码矢量变化不大(收敛)，就完成码书的训练，否则重新进行②、③。图1为传统的LBG方法流程框图。LBG算法的思想比较简单，但是其对初始码书依赖性强，计算量大。

发明内容

本发明针对现有高光谱图像压缩方案，在获得较高的图像恢复质量时，计算复杂度高的问题，提出了一种快速压缩高光谱信号的矢量量化方案。该方案能在保证图像恢复质量几乎不变的情况下，大幅度降低计算复杂度，快速高效地实现高光谱图像的压缩。

一种快速压缩高光谱信号的矢量量化方法，包括步骤：读取高光谱图像数据，截取部分转化为2维的矢量数据并分割成M个小部分，其中，M是2到9之间的一个整数，对分割出的各部分矢量数据分别进行Hadamard变换；对Hadamard变换域的各部分矢量数据分别作离散度排序得到各部分离散度排序后的矢量数据(B_1D、B_2D，…，B_MD)；提取离散度排序后矢量数据行矢量的前n维分量组成训练矢量集(B_1T、B_2T，…，B_MT)，剩下的k-n维分量组成对应部分的余下矢量数据(B_1R、B_2R，…，B_MR)，其中，n＝log₂(k)，k＝2ⁿ为该部分矢量数据行矢量的总维数；对训练矢量集(B_1T、B_2T，…，B_MT)分别进行训练，得到矢量量化最后一次迭代产生的各部分训练矢量码书(Y_1die、Y_2die，…，Y_Mdie)和对应的训练矢量码字索引(I_1die、I_2die，…，I_Mdie)；根据最后一次迭代产生的训练矢量码书和对应的训练矢量码字索引计算余下矢量数据的码书，将训练矢量码书和余下矢量码书重组连接，通过离散度反排序和hadamard反变换恢复出空域完整码书；打包空域完整码书和对应的空域完整码字索引，存储或传输。其中，B_MD代表第M部分离散度排序后的矢量数据，B_MT代表第M部分训练矢量集，B_MR代表第M部分余下矢量数据，Y_Mdie、I_Mdie分别代表最后一次迭代产生的第M部分训练矢量码书和第M部分训练训练码码字索引。

分割矢量数据具体包括：读取3维高光谱图像数据，截取部分相同空间位置的像素值转化为2维的矢量数据并分割为M个小部分，M最好取2到9之间的一个整数，分割原则为：每部分的矢量维数都为2ⁿ或最接近2ⁿ，维数不为2ⁿ的部分通过末尾补零使其维数为最接近2ⁿ，各部分矢量维数之和尽量接近原矢量维数，其中，n为不小于5的正整数。所述离散度排序具体为：将行矢量某一维分量的最大值减去最小值得到差值获取行矢量该维分量的离散度，将离散度值分别按照降序排序，得到对应的排序索引作为各部分矢量数据的离散度排序索引，按照离散度排序索引将对应的各部分矢量数据的行矢量分量重排序，得到离散度排序后的各部分矢量数据(B_1D、B_2D，…，B_MD)。提取训练矢量集具体为：提取离散度排序后的各部分矢量数据的前n维分量组成对应部分的训练矢量集(B_1T、B_2T，…，B_MT)，其中，n＝log₂(k)，k＝2ⁿ为该部分矢量数据行矢量的维数。对训练矢量集(B_1T、B_2T，…，B_MT)进行训练具体为：按照行矢量第一维分量值的大小对各部分训练矢量进行升序排序，排序索引作为各部分矢量排序索引，按照各部分矢量排序索引重新排序各部分训练矢量，将排序后的训练矢量按照码书大小进行平均分组，依次选取每组第一个矢量组成初始码书，分组号作为初始码字索引；根据初始码书和初始码字索引搜索训练矢量对应的最佳匹配码字，并将训练矢量划分到对应的胞腔中，记录对应的码字索引，直到所有矢量训练完成，更新码书，以各胞腔的质心代替原来胞腔对应的码字。计算余下矢量数据(B_1R、B_2R，…，B_MR)码书具体为：按照最终训练矢量码字索引(I_1T、I_2T，…，I_MT)分别把对应部分余下矢量数据(B_1R、B_2R，…，B_MR)中的行矢量依次分配到对应的胞腔中，以各胞腔质心作为对应胞腔余下矢量数据的最终余下矢量码字(Y_1R、Y_2R，…，Y_MR)。所述搜索最佳匹配码字具体包括如下步骤：①根据训练矢量X的第1维分量X₁、当前码字Y_j的第1维分量Y_j1，根据公式：D₁＝(X₁-Y_j1)²获取第一维分量差平方D₁，如果D₁≥D_min，且X₁≥Y_j1，排除码字Y_i(i＝1,…,j)，如果D₁≥D_min，且X₁≤Y_j1时，排除码字Y_i(i＝j,…,N)，否则执行步骤②；②根据公式：计算训练矢量X的方差V_X、码字Y_j的方差V_Y，根据公式D₂＝D₁+(V_X-V_Y)²确定D₂，D₂为第一维分量差平方与方差差平方之和，若不等式D₂≥D_min成立则排除码字Y_j，否则执行步骤③；③根据公式计算训练矢量X与当前判断码字Y_j的累积失真，如D_q(X,Y_j)≥D_min成立，排除码字Y_j，否则把训练矢量X划分到对应的胞腔中，其中，X_i为训练矢量X的第i维分量，Y_ji为码字Y_j的第i维分量，n为训练矢量的维数，D_min为当前最小失真，N为码书尺寸。训练矢量码书和余下矢量码书的重组连接方式具体为：按最终训练矢量码字(Y_1T、Y_2T，…，Y_MT)在前，最终余下矢量码字(Y_1R、Y_2R，…，Y_MR)在后的顺序，把相同序号的最终训练矢量码字和最终余下矢量码字连接起来，合并重组成2ⁿ维的最终完整码字，其中I_MT为第M部分最终训练矢量码字索引，Y_MT为第M部分最终训练矢量码书，Y_MR为第M部分最终余下矢量码书。

采用矢量量化对数据进行压缩编码，编码后只需要传输最终码书以及各个训练矢量对应的码字索引，压缩比大，解码过程只需要在接收端根据码字索引值在码书中查找到对应的码字代替训练矢量，解码过程简单。高光谱图像的数据量庞大，且有很强的谱间相关性，因此将相同空间位置像素的像素值组织成一个矢量数据，再用矢量量化进行压缩编码，可以有效减小数据量，又很好地保存光谱信息，且操作较简单，速度较快。

附图说明

图1传统的LBG方法流程框图；

图2矢量维数分割量化的超光谱图像压缩方法流程图；

图3本发明方案流程图。

具体实施方式

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

数据预处理阶段：读取高光谱图像三维数据，以同一空间位置对应像素各波段像素值作为一个矢量数据，即把高光谱图像每一个像元对应的光谱矢量对应一个矢量数据，将三维的原始数据转化为二维的矢量矩阵数据。然后按照2ⁿ(n为不小于5的正整数)准则将所有矢量数据分为较小的部分，使总维数最接近原维数，如果某部分矢量维数不满足2ⁿ，通过在矢量末尾补零进行维数扩展，使每个分割出的矢量部分中矢量的维数都要满足2ⁿ。对分割出的各部分矢量数据分别作Hadamard变换，得到Hadamard变换域的各部分矢量数据。对Hadamard变换域的各部分矢量数据依次进行如下操作：对该部分矢量数据的各维分量，首先计算该维分量对应的离散度并记录，该维分量的离散度指该维分量的最大值减去该维分量的最小值的差值；然后按照这些离散度值的大小进行降序排序，记录分量离散度排序索引；再按照离散度排序索引重新排序该部分矢量数据的各维分量；最后截取排序后的该部分矢量数据的前n＝log₂(k)维(该部分矢量的维数为k＝2ⁿ)组成该部分矢量数据对应的要参与训练的训练矢量并保存，且保存剩下的2ⁿ-n维作为该部分矢量数据对应的余下矢量数据。

初始化阶段：输入分割出的第一部分矢量数据对应的训练矢量，依次进行如下操作：按照训练矢量第一维分量值的大小进行升序排序，记录排序索引作为对应的矢量排序索引，再按照训练矢量排序索引对各个训练矢量进行重新排序，得到矢量排序后的训练矢量，然后对矢量排序后的训练矢量按照码书大小进行平均分组，依次选取每组第一个矢量组成第一部分训练矢量对应的初始码书并保存，且保存各个矢量对应的组号作为第一部分训练矢量对应的初始码字索引。依次输入分割出的第二部分直至最后一部分矢量数据对应的训练矢量，按照对第一部分训练矢量相同的操作，分别得到各部分训练矢量对应的初始码书和初始码字索引。最后设置迭代次数或失真阈值作为训练结束的条件。

训练阶段：输入排序后的第一部分训练矢量进行如下操作：利用高效三步排除不等式算法快速搜索训练矢量对应的最佳匹配码字，并将训练矢量划分到对应的胞腔中，记录对应的码字索引，直到所有当前训练的矢量都已经训练完成时，更新码书，以各胞腔的质心来代替原来胞腔对应的码字。依次输入分割出的第二部分直至最后一部分矢量数据对应的训练矢量进行训练，按照与第一部分训练矢量相同的操作逐步训练完所有训练矢量。

中断检查阶段：判断迭代次数是否已经达到预先设定的值，或者判断本次迭代产生的平均失真与前一次迭代产生的平均失真之间的相对误差是否小于设定的失真阈值。若满足迭代结束的条件，则终止训练，并保存最后一次迭代产生的各部分训练矢量码书和训练矢量码字索引；否则，返回到训练阶段继续训练各部分训练矢量。

码字重组阶段：输入分割出的第一部分矢量数据对应的训练矢量排序索引、余下矢量数据、以及最后一次迭代产生的训练矢量码书和训练矢量码字索引，进行如下操作：首先，按照训练矢量排序索引调整对应的最后一次迭代产生的训练矢量码字索引，使训练矢量码字索引与原来的训练矢量顺序匹配，把恢复顺序的训练矢量码字索引作为对应的最终训练矢量码字索引，把最后一次迭代产生的训练矢量码书直接作为对应的最终训练矢量码书；再按照最终训练矢量码字索引，把余下矢量数据分配到对应的胞腔中，以各胞腔质心作为胞腔中余下矢量数据对应的最终余下矢量码字；然后将最终训练矢量码字和最终余下矢量码字按码字序号合并重组，即成为Hadamard变换域的完整码字；再将Hadamard变换域的完整码字按照分量离散度排序索引重新排序，恢复到原来的分量顺序；最后将恢复顺序的完整码字全部作Hadamard反变换，如果本部分矢量数据在数据预处理阶段阶段作了添加零扩维处理，则Hadamard反变换后去除矢量末尾添加零的那几维，得到空域的最终完整码字，把最终训练矢量码字索引作为本部分矢量数据的最终码字索引，打包空域的最终完整码字和最终码字索引进行存储和传输，完成本部分矢量数据的压缩。依次输入分割出的第二至最后一部分矢量数据对应的余下矢量数据、最终训练矢量码书和最终训练矢量码字索引按照与第一部分相同的操作完成各部分数据的最终压缩。

图3本发明方案流程图。具体包括以下几个方面：

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

⑵按照2ⁿ(n为不小于5的正整数)准则将矩阵B的列矢量分成最接近原维数的M个部分，M是2到9之间的一个整数，如果某部分矢量维数不满足2ⁿ，通过在矢量末尾补零把维数扩展到2ⁿ，得到分割出的M个部分的矢量数据B₁、B₂，…，B_M；

⑶对分割出的各部分矢量数据分别作Hadamard变换，得到Hadamard变换域的各部分矢量数据B_1H、B_2H，…，B_MH；

⑷输入第一部分Hadamard变换域的矢量数据B_1H，依次进行如下操作：对B_1H的行矢量的各维分量，首先计算该维分量对应的离散度dis₁保存，该维分量的离散度指该维分量的最大值减去该维分量的最小值的差值；然后按照这些离散度值的大小进行降序排序，记录对应的分量离散度排序索引I_1d；再按照离散度排序索引重新排序该部分矢量数据的各维分量；最后截取排序后的该部分矢量数据的前n维(该部分矢量的维数为2ⁿ)组成该部分矢量数据对应的要参与训练的训练矢量B_1T并保存，且保存剩下的2ⁿ-n维分量组成的该部分矢量数据对应的余下矢量数据B_1R。依次输入第二至最后一部分部分Hadamard变换域的矢量数据B_2H、B_3H，…，B_MH，按照与B_1H相同的操作进行处理，分别得到各部分矢量数据对应的分量离散度排序索引I_2d、I_3d，…，I_Md、训练矢量B_2T、B_3T，…，B_MT和余下矢量数据B_2R、B_3R，…，B_MR。

⑸输入分割出的第一部分矢量数据对应的训练矢量B_1T，依次进行如下操作：按照B_1T第一列值的大小进行升序排序，记录排序索引作为对应的矢量排序索引I_1V，再按照I_1V重新排序B_1T中的各个行矢量，得到矢量矩阵B_1V，并保存对应的训练矢量排序索引I_1V；然后对B_1V的行矢量按照码书大小N进行平均分组，依次选取每组第一个矢量组成第一部分训练矢量对应的初始码书Y_1F并保存，且保存对应的初始码字索引I_1F。依次输入分割出的第二部分直至最后一部分矢量数据对应的训练矢量B_2T、B_3T，…，B_MT，按照对第一部分训练矢量B_1T相同的操作，分别得到各部分排序后的训练矢量B_2V、B_3V，…，B_MV以及对应的训练矢量排序索引I_2V、I_3V，…，I_MV，初始码书和初始码字索引。最后设置迭代次数ite或失真阈值e作为训练结束的条件。

⑹输入第一部分排序后的训练矢量B_1V，和对应的初始码书、初始码字索引，依次进行如下操作：快速搜索最佳匹配码字，依次将B_1V中的训练矢量划分到其最佳匹配码字所对应的胞腔中，记录对应的码字索引，直到B_1V中所有的行矢量都已经训练完成时，更新码书，以各胞腔的质心来代替原来胞腔对应的码字，得到新的码书。依次输入排序后的第二部分直至最后一部分训练矢量B_2V、B_3V，…，B_MV，和对应的初始码书和初始码字索引，按照与第一部分训练矢量相同的操作逐步训练完各部分排序后的训练矢量，得到对应的码字索引和码书。

⑺判断迭代次数是否已经达到预先设定的值ite，或者判断本次迭代产生的平均失真与前一次迭代产生的平均失真之间的相对误差是否小于设定的失真阈值e。若满足迭代结束的条件，则终止训练，并保存最后一次迭代产生的各部分训练矢量码书和训练矢量码字索引；否则，返回到训练阶段继续训练各部分训练矢量。

⑻输入分割出的第一部分矢量数据对应的训练矢量排序索引I_1V、余下矢量数据B_1R、以及最后一次迭代产生的训练矢量码书和训练矢量码字索引，进行如下操作：首先，按照训练矢量排序索引I_1V调整对应的最后一次迭代产生的第一部分训练矢量码字索引，使训练矢量码字索引与原来的训练矢量顺序匹配，把恢复顺序的训练矢量码字索引作为对应的最终训练矢量码字索引I_1T，把最后一次迭代产生的训练矢量码书直接作为对应的最终训练矢量码书Y_1T；再按照最终训练矢量码字索引I_1T，把余下矢量数据分配到对应的胞腔中，以各胞腔质心作为胞腔中余下矢量数据对应的最终余下矢量码字Y_1R；然后将n维的最终训练矢量码字Y_1T和2ⁿ-n维的最终余下矢量码字Y_1R按码字序号合并重组成2ⁿ维的码字，再将合并成的码字的各分量按照对应的分量离散度排序索引I_1d重新排序，使码字各分量恢复到原来的顺序，得到Hadamard域的最终完整码字Y_1H；最后将恢复分量顺序的Hadamard变换域最终完整码字Y_1H中的码字全部作Hadamard反变换，如果本部分矢量数据在步骤⑵作了添加零扩维处理，则Hadamard反变换后去除矢量末尾添加零的那几维，得到空域的最终完整码字Y_1N，把最终训练矢量码字索引I_1T作为本部分矢量数据的最终码字索引I_1N，打包该部分矢量数据对应的空域最终完整码字Y_1N和最终码字索引I_1N进行存储和传输，完成本部分矢量数据的压缩。依次输入分割出的第二至最后一部分矢量数据对应的训练矢量排序索引I_2V、I_3V，…，I_MV、余下矢量数据B_2R、B_3R、…，B_MR、以及最后一次迭代产生的训练矢量码书和训练矢量码字索引，按照与第一部分相同的操作得到各部分矢量数据对应的空域最终完整码字Y_2N、Y_3N，…，Y_MN和最终码字索引I_2N、I_3N，…，I_MN，进行存储和传输，完成各部分矢量数据的压缩。

高效三步排除不等式算法具体可为：

第一步：输入要编码的当前训练矢量X，按照上一次迭代后X对应的码字索引找到相应的码字Y_b，如果是首次迭代则按照初始码字索引找到初始码书中对应的码字Y_b，计算X和Y_b之间的欧式距离作为当前最小失真D_min；按照以Y_b为中心上下搜索的顺序，输入下一个要判断的码字Y_j，计算X第一维分量X₁与码字Y_j第一维分量Y₁的差的平方D₁＝(X₁-Y_j1)²，与当前最小失真进行比较，即判断不等式D₁≥D_min是否成立，若成立则进一步判断X₁与Y_j1的关系，若X₁≥Y_j1，进入步骤a)，若X₁≤Y_j1时，进入步骤b)；若不等式D₁≥D_min不成立，则进入第二步判断；

1排除码字Y_i，i＝1,…,j；

2排除码字Y_i，i＝j,…,N(j为当前搜索码字的下标，N表示码书尺寸)；

第二步：分别计算Hadamard变化域训练矢量X的方差和码字Y_j的方差其中变量i表示第i维分量，变量n表示训练矢量的维数，判断不等式D₂＝D₁+(V_X-V_Y)²≥D_min是否成立，若成立则排除码字Y_j，否则进入第三步判断。

第三步：采用部分失真搜索算法(PDS算法)第三步搜索，即计算训练矢量X各维分量的累积失真n是训练矢量的维数，q从2逐渐变化到n的过程可以构成n-1个部分失真不等式，若其中有某一步不等式成立，则排除码字Y_j，否则将判断继续进行，直到q＝n，若还未排除码字Y_j，则当前码字Y_j为训练矢量X对应的新的失真最小的码字，即最匹配码字，并把X对应的码字索引更新为j，把X对应的最小失真D_min更新为当前码字Y_j和X之间的欧式距离

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

⑴构造矢量数据

读取高光谱图像，截取每个波段中的相同空间位置的一段图像块作为输入信源，根据截取的图像块尺寸和光谱波段数构建三维数据矩阵，以同一空间位置对应像素各波段像素值作为一个矢量数据，即把截取的图像块每一个像元对应的光谱矢量对应一个矢量数据，把三维数据转化为二维矩阵数据。本发明中实验使用的高光谱图像是Lunar Lake和Jasper的第一场景，来源于美国喷气实验室提供的免费高光谱数据集AVIRIS，这些数据源被绝大多数研究者使用，研究结果具有可比性。AVIRIS中的图像大小为512×614×224，表示图像含有224个光谱波段，每个波段对应的二维图像大小为512×614。为了处理方便，分别截取Lunar Lake和Jasper第一场景中每个波段相同空间位置256×256大小的图像，组成256×256×224的三维图像块数据，提取此三维图像块每一相同空间位置对应的224个波段的像素值作为一个矢量数据，将此三维图像块数据转化为大小为65536×224的二维矩阵数据B，B的每一行代表一个224维矢量数据，对应一个空间像元，B的每一列表示一个光谱波段，对应相同空间位置一幅完整的256×256的二维图像。

⑵分割矢量数据

按照2ⁿ(n为不小于5的正整数)准则将数据矩阵B的列矢量分成最接近原维数的M个部分，M是2到9之间的一个整数，如果某部分矢量维数不满足2ⁿ，通过在矢量末尾补零把维数扩展到2ⁿ，得到分割出的M个部分的矢量数据B₁、B₂，…，B_M。由于使用的具体数据B中的矢量维数为224，所以根据实际恰好把B分割为矢量维数为32、64、128的三个部分，即M＝3，得到大小分别为65536×32、65536×64、65536×128的三个数据矩阵B₁、B₂、B₃，使每一部分数据矩阵中的矢量维数都满足2ⁿ条件，B₁中的矢量维数为2⁵、B₂中的矢量维数为2⁶、B₃中的矢量维数为2⁷。

⑶对分割出的各部分矢量数据分别作Hadamard变换

对分割出的各部分矢量数据B₁、B₂，…，B_M分别作Hadamard变换，得到Hadamard变换域的各部分矢量数据B_1H、B_2H，…，B_MH。对一个维数为2ⁿ的行矢量进行Hadamard变换即将矢量右乘一个2ⁿ×2ⁿ的Hadamard方阵，所以对本文中分割出的三部分矢量数据B₁、B₂、B₃分别作Hadamard变换，具体操作步骤为：将B₁中每个行矢量右乘一个2⁵×2⁵的Hadamard方阵，得到对应的Hadamard变换域的数据矩阵B_1H；将B₂中每个行矢量右乘一个2⁶×2⁶的Hadamard方阵，得到对应的Hadamard变换域的数据矩阵B_2H；将B₃中每个行矢量右乘一个2⁷×2⁷的Hadamard方阵，得到对应的Hadamard变换域的数据矩阵B_3H。

由于Hadamard方阵中只含有1和-1两种元素，所以Hadamard变换仅仅需要简单的加减运算不需要乘法运算，且Hadamard变换有快速算法，所以变换所需要的计算量小，速度快。不仅如此，矢量的Hadamard变换还具有一些优越性质，下面作具体介绍：

令H_n为2ⁿ×2ⁿ的Hadamard矩阵，记x为空域中维数为k＝2ⁿ的行矢量，X为矢量x经Hadamard变换后的矢量，则x的Hadamard变换定义为X＝xH_n，矢量Hadamard变换的性质可以表述为：

性质1X₁＝s_x，s_x是空域中矢量x各维分量之和，X₁是x经Hadamard变换后所得矢量X的第一维分量，即Hadamard域中矢量的第一维分量是对应空域中矢量各分量的和值；

性质2l_x ²表示空域中矢量x的范数，L_x ²表示x经Hadamard变换后所得矢量X的范数，即Hadamard域中矢量范数是对应空域中矢量范数的k倍；

性质3D(X,Y_j)＝kd(x,y_j)，d(x,y_j)为空域中矢量x和码字y_j的欧氏距离，D(X,Y_j)是对应的Hadamard域中矢量X和码字Y_j的欧式距离，即Hadamard域中的矢量失真是空域中失真的k倍，由此看出，在Hadamard域搜索最佳码字和在空域搜索是等价的。

⑷对各部分Hadamard变换域的矢量数据分别作分量离散度排序并构造出训练矢量集和余下矢量集

定义一组矢量数据某一维分量的离散度为该维分量的最大值减去该维分量的最小值的差值，例如B_1H中所有行矢量第一维分量的离散度为B_1H的第一列的最大值减去第一列的最小值的差值。则可以分为以下几步进行：

第一步，计算各部分Hadamard变换域的矢量数据B_1H，B_2H，…，B_MH中所有行矢量的各维分量的离散度dis₁、dis₂，…，dis_M，然后按照这些离散度值的大小进行降序排序，记录对应的分量离散度排序索引I_1d，I_2d，…，I_Md；

第二步，按照分量离散度排序索引I_1d，I_2d，…，I_Md重新排序各部分Hadamard变换域的矢量数据B_1H，B_2H，…，B_MH中行矢量的各维分量，得到对应的分量离散度排序后的矢量数据B_1D，B_2D，…，B_MD；

第三步，对各部分分量离散度排序后的矢量数据B_1D，B_2D，…，B_MD分别截取前n＝log₂(k)维(该部分矢量的维数为k＝2ⁿ)组成对应部分的训练矢量B_1T，B_2T，…，B_MT并保存，且保存剩下的2ⁿ-n维分量组成对应部分的余下矢量数据B_1R，B_2R，…，B_MR。

由于本文中M＝3，分割出的三部分矢量数据中，第一部分矢量维数为2⁵，所以截取第一部分分量离散度排序后的矢量数据B_1D的前log₂(2⁵)＝5维组成第一部分训练矢量B_1T并保存，且保存剩下的2⁵-5维分量组成第一部分余下矢量数据B_1R；第二部分矢量维数为2⁶，所以截取第二部分分量离散度排序后的矢量数据B_2D的前log₂(2⁶)＝6维组成第二部分训练矢量B_2T并保存，且保存剩下的2⁶-6维分量组成第二部分余下矢量数据B_2R；第三部分矢量维数为2⁷，所以截取第三部分分量离散度排序后的矢量数据B_3D的前log₂(2⁷)＝7维组成第三部分训练矢量B_3T并保存，且保存剩下的2⁷-7维分量组成第三部分余下矢量数据B_3R。

实验表明哈达码变换域的矢量数据能量主要集中在离散度大的分量，各矢量分量经过离散度排序后，能量主要集中在低维部分，且约90％以上的能量集中在前n＝log₂(k)维(k为矢量维数)，因此，本方案只选取离散度排序后矢量的前n＝log₂(k)维作为训练数据，即用集中了矢量绝大部分能量的低维少数数据代表整个矢量参加聚类训练，然后用训练数据的聚类类别代表整个矢量的类别。这样不仅能达到比较准确地实现聚类，还能达到大幅度降低计算量的目的。

⑸训练初始化

第一步，按照各部分训练矢量B_1T，B_2T，…，B_MT第一列值的大小分别进行升序排序，记录各个排序索引作为对应的各部分矢量排序索引I_1V，I_2V，…，I_MV，再按照矢量排序索引I_1V，I_2V，…，I_MV分别重新排序B_1T，B_2T，…，B_MT中的行矢量，得到矢量排序后的各部分训练矢量B_1V，B_2V，…，B_MV；

第二步，分别对各部分矢量排序后的训练矢量B_1V，B_2V，…，B_MV的行矢量按照码书大小N进行平均分组，依次选取每组第一个矢量组成各部分训练矢量对应的初始码书，且记录各个矢量对应的组号作为各部分训练矢量对应的初始码字索引；

第三步，打包保存各部分矢量排序后训练矢量B_1V，B_2V，…，B_MV、对应的各部分训练矢量初始码书和训练矢量初始码字索引。

第四步，根据编码质量要求设置迭代次数ite或失真阈值e作为训练结束的条件(迭代次数越多，编码质量越好；失真阀值越小，编码质量越好)。

对于本文的数据来说，由于M＝3，所以仅需按以上第一步至第四步操作，得到对应的三部分矢量排序后的矩阵B_1V，B_2V，B_3V、训练矢量排序索引I_1V，I_2V，I_3V、训练矢量初始码书和训练矢量初始码字索引，完成初始化。

⑹训练阶段

引入变量t(t∈{1,2,…,M})，初始化t＝1，t代表第t部分，则本部分具体操作步骤可描述如下：

第一步，输入当前要处理的训练矢量数据B_tV(B_tV∈{B_1V,B_2V,…,B_MV})、上一次迭代产生的第t部分训练矢量码书和对应的训练矢量码字索引(如果是首次迭代则输入第t部分训练矢量初始码书和对应的训练矢量初始码字索引)。

第二步，利用高效三步排除不等式算法快速搜索当前处理的训练矢量矩阵B_tV中的训练矢量对应的失真最小的码字，记录对应的训练矢量码字索引，并将训练矢量划分到对应的胞腔中。

第三步，检查当前处理中的矢量矩阵B_tV中的训练矢量是否已经训练完，如果没有，则按顺序输入B_tV中下一个训练矢量作为当前训练矢量，利用高效三步排除不等式算法继续训练B_tV中的其它训练矢量，如果B_tV中的训练矢量已经训练完，则更新码书，以各胞腔的质心来代替原来胞腔对应的码字，得到本轮迭代对应的第t部分新的训练矢量码书。

第四步，判断各部分训练矢量数据B_1V，B_2V，…，B_MV是否都已经训练完成，如果已经训练完，则本轮迭代结束；如果没有，则顺序输入下一个训练矢量数据矩阵，继续训练。即判断t>M是否成立，如果成立，则保存本轮迭代产生的各部分训练矢量码书和对应的各部分训练矢量码字索引，并转入步骤⑺；如果不成立则t＝t+1，转入⑹训练阶段中的第一步。

上文中所述的高效三步排除不等式算法具体操作步骤为：

①输入要编码的当前训练矢量X，按照上一次迭代后X对应的码字索引找到相应的码字Y_b，如果是首次迭代则按照初始码字索引找到初始码书中相应的码字Y_b，计算X和Y_b之间的欧式距离作为当前最小失真D_min，按照以Y_b为中心上下搜索的顺序，输入下一个要判断的码字Y_j，计算X第一维分量X₁与码字Y_j第一维分量Y₁的差的平方D₁＝(X₁-Y_j1)²，与当前最小失真D_min进行比较，即判断不等式D₁≥D_min是否成立，若成立则进一步判断Y_j1与X₁的关系，若X₁≥Y_j1，进入步骤a)，若X₁≤Y_j1时，进入步骤b)；若不等式D₁≥D_min不成立，则进入步骤②进行第二步判断；

a)排除码字Y_i，i＝1,…,j；

b)排除码字Y_i，i＝j,…,N(j为当前搜索码字的下标，N表示码书尺寸)；

②分别计算Hadamard变换域当前训练矢量X的方差和码字Y_j的方差其中变量i表示矢量的第i维分量，变量n表示训练矢量的维数，判断不等式D₂＝D₁+(V_X-V_Y)²≥D_min是否成立，若成立则排除码字Y_j，否则进入下一步判断；

③采用部分失真搜索算法(PDS算法)进行判断，即计算当前训练矢量X与当前判断码字Y_j的累积失真n是训练矢量的维数，q从2逐渐变化到n的过程可以构成n-1个部分失真不等式，若其中有某一步不等式成立，则排除码字Y_j，否则q增加1继续判断，直到q＝n，若还未排除码字Y_j，则当前码字Y_j为训练矢量X对应的新的失真最小的码字，即最匹配码字，并把X对应的码字索引更新为j，把X对应的最小失真D_min更新为当前码字Y_j和X之间的欧式距离

⑺中断检查

判断迭代次数是否已经达到预先设定的值ite，或者判断本次迭代产生的平均失真与前一次迭代产生的平均失真之间的相对误差是否小于设定的失真阈值e。若满足迭代结束的条件，则终止训练，并保存最后一次迭代产生的各部分训练矢量码书Y_1die、Y_2die，…，Y_Mdie和对应的训练矢量码字索引I_1die、I_2die，…，I_Mdie；否则，返回到步骤⑹训练阶段进入下一次迭代训练。

⑻码字重组

引入变量c(c∈{1,2,…,M})，初始化c＝1，c代表第c部分，则本部分的具体操作步骤可描述如下：

第一步，输入数据：输入当前要处理的那部分矢量数据对应的训练矢量排序索引I_cV、余下矢量数据B_cR、以及对应的最后一次迭代产生的训练矢量码书Y_cdie和训练矢量码字索引I_cdie。

第二步，调整训练矢量码字索引：由于最后一次迭代产生的训练矢量码书Y_1die、Y_2die，…，Y_Mdie和训练矢量码字索引I_1die、I_2die，…，I_Mdie是分别和B_1V、B_2V，…，B_MV对应的，而不是和矢量排序前的数据B_1T、B_2T，…，B_MT对应，所以需要调整码字索引。首先按照本部分训练矢量排序索引I_cV调整最后一次迭代产生的本部分训练矢量码字索引I_cdie，得到新的索引，使新的索引与最后一次迭代产生的码书Y_cdie和原来的训练矢量B_cT对应，然后把调整后的训练矢量码字索引作为本部分最终训练矢量码字索引I_cT，把最后一次迭代产生的训练矢量码书Y_cdie直接作为第本部分最终训练矢量码书Y_cT。

第三步，余下矢量数据码书设计：按照最终训练矢量码字索引I_cT，把各部分余下矢量数据B_cR中的矢量依次分配到对应的胞腔中，以各胞腔质心作为对应胞腔的余下矢量数据的余下矢量码字Y_cR。

第四步，码字重组：将n维的最终训练矢量码字Y_cT和2ⁿ-n维的最终余下矢量码字Y_cR按码字序号合并重组成2ⁿ维的最终完整码字Y_cZ，即按最终训练矢量码字在前最终余下矢量码字在后的顺序，把相同序号的最终训练矢量码字和最终余下矢量码字连接起来。

对本文数据来说，M＝3，c∈{1,2,3}。对应c＝1时，就是将5维的第一部分训练矢量码字Y_1T和2⁵-5维的最终余下矢量码字Y_1R按码字序号合并重组成2⁵维的最终完整码字Y_1Z，Y_1Z的前5维为Y_1T，Y_1Z的后2⁵-5维为Y_1R；c＝2就是将6维的第二部分训练矢量码字Y_2T和2⁶-6维的最终余下矢量码字Y_2R按码字序号合并重组成2⁶维的最终完整码字Y_2Z，Y_2Z的前6维为Y_2T，Y_2Z的后2⁶-6维为Y_2R；c＝3就是将7维的第三部分训练矢量码字Y_3T和2⁷-7维的最终余下矢量码字Y_3R按码字序号合并重组成2⁷维的最终完整码字Y_3Z，Y_3Z的前7维为Y_3T，Y_3Z的后2⁷-7维为Y_3R。

第五步，恢复码字分量顺序：将重组成的最终完整码字Y_cZ的各分量按照对应的分量离散度排序索引I_cd重新排序，使最终完整码字各分量恢复到原来的顺序，得到Hadamard变换域的最终完整码字Y_cH。

第六步，码字Hadamard反变换：将恢复分量顺序的Hadamard变换域最终完整码字Y_cH中的码字全部作Hadamard反变换，如果矢量数据在⑵分割矢量数据阶段作了添加零扩维处理，则去除矢量末尾添加零的那几维，得到空域最终完整码字Y_cN。

第七步，打包存储：把最终训练矢量码字索引I_cT直接作为空域最终完整码字索引I_cN，打包空域最终完整码字Y_cN和对应的最终完整码字索引I_cN，进行存储或传输，完成矢量数据的压缩。

第八步，判断各部分码字是否都已经重组完成，如果已经完成，则压缩结束；如果没有，则顺序输入下一部分相关数据，继续重组下一部分码字。即，判断c>M是否成立，如果成立则保存本轮迭代产生的各部分训练矢量码书和对应的各部分训练矢量码字索引，并转入步骤⑺，如果不成立则t＝t+1，转入⑹训练阶段中的第一步。

在MATLAB7.9软件平台下结合附图对本发明方案进行详细说明。

使用512×614×224规格的高光谱图像Lunar Lake和Jasper作为实验数据，其每个像素用两个字节的带符号整数表示。仿真实验对本发明方案、全搜索LBG方法和矢量维数分割量化的超光谱图像压缩方法3种方法进行了仿真对比。以下将其简称为分割算法。本发明方案的具体实施步骤如下：

数据预处理阶段：

步骤1：利用MATLAB函数库中的fopen函数和fread函数快速读取高光谱图像，得到3维的高光谱数据矩阵A，A含有512行，614列，224个波段，简记为A(512,614,224)；

步骤2：截取A的前256行、前256列和所有224个波段得到新的3维矩阵P(256,256,224)；

步骤3：将3维矩阵P转化为2维矩阵B，以便进行矢量量化处理。变换方式是将P的前2维转换为B的第1维，P的第3维作为矩阵B的第2维，最终生成65536行224列的B，简记为B(65536,224)；

步骤4：将B的所有列矢量分割为三部分，第一部分为B的1-32列，简记为B₁(65536,32)，第二部分为B的33-96列，简记为B₂(65536,64)，第三部分为B的97-128列，简记为B₃(65536,128)；

步骤5：将B₁中每个行矢量右乘一个32×32的Hadamard方阵，得到B_1H(65536,32)，将B₂中每个行矢量右乘一个64×64的Hadamard方阵，得到B_2H(65536,64)，将B₃中每个行矢量右乘一个128×128的Hadamard方阵，得到B_3H(65536,128)；

步骤6：计算B_1H中行矢量各维分量的离散度dis₁(1,32)，按照这些离散度值的大小进行降序排序，记录排序索引I_1d(1,32)，计算B_2H中行矢量各维分量的离散度dis₂(1,64)，按照这些离散度值的大小进行降序排序，记录排序索引I_2d(1,64)，计算B_3H中行矢量各维分量的离散度dis₃(1,128)，按照这些离散度值的大小进行降序排序，记录排序索引I_3d(1,128)；

步骤7：按照I_1d重新排序B_1H的各列，得到B_1D(65536,32)，按照I_2d重新排序B_2H的各列，得到B_2D(65536,64)，按照I_3d重新排序B_3H的各列，得到B_3D(65536,128)；

步骤8：截取B_1D的1-5列组成B_1T(65536,5)，且保存B_1D的6-32列组成B_1R(65536,27)，截取B_2D的1-6列组成B_2T(65536,6)，且保存B_2D的7-64列组成B_2R(65536,58)，截取B_3D的1-7列组成B_3T(65536,7)，且保存B_3D的8-128列组成B_3R(65536,121)；

初始化阶段：

步骤1：按B_1T第一列值的大小进行升序排序，记录排序索引I_1V(1,65536)，按照I_1V对B_1T的各个行矢量进行重新排序，得到B_1V(65536,5)，按B_2T第一列值的大小进行升序排序，记录排序索引I_2V(1,65536)，按照I_2V对B_2T的各个行矢量进行重新排序，得到B_2V(65536,6)，按B_3T第一列值的大小进行升序排序，记录排序索引I_3V(1,65536)，按照I_3V对B_3T的各个行矢量进行重新排序，得到B_3V(65536,7)；

步骤2：设置码书大小为N，设定第一部分初始码字索引矩阵I_1ini(1,65536)，把B_1V的行矢量平均分成N组，依次提取每组第一个矢量组成第一部分初始码书Y_1ini(N,5)，并记录各个矢量对应的组号作为第一部分初始码字索引然后将B_1V中的训练矢量按照初始索引I_1ini分配到相应的胞腔中，即具有相同索引值的矢量归属于相同的胞腔；设定第二部分初始码字索引矩阵I_2ini(1,65536)，把B_2V的行矢量平均分成N组，依次提取每组第一个矢量组成第二部分初始码书Y_2ini(N,6)，并记录各个矢量对应的组号作为第二部分初始码字索引然后将B_2V中的训练矢量按照初始索引I_2ini分配到相应的胞腔中，即具有相同索引值的矢量归属于相同的胞腔；设定第三部分初始码字索引矩阵I_3ini(1,65536)，把B_3V的行矢量平均分成N组，依次提取每组第一个矢量组成第三部分初始码书Y_3ini(N,7)，并记录各个矢量对应的组号作为第三部分初始码字索引然后将B_3V中的训练矢量按照初始索引I_3ini分配到相应的胞腔中，即具有相同索引值的矢量归属于相同的胞腔；

步骤3：设置迭代次数ite＝10作为训练结束的条件。

训练阶段和中断检查阶段：

步骤1：引入变量t(t∈{1,2,3})，初始化t＝1，t代表第t部分训练矢量；

步骤2：输入当前要处理的训练矢量数据B_tV(B_tV∈{B_1V,B_2V,B_3V})并计算B_tV中所有行矢量的方差并保存，对于B_tV中的某个行矢量X，其方差V_X的计算方式为其中变量i表示矢量的第i维分量，变量n表示B_tV中行矢量的维数，即B_tV的列数；

步骤3：引入变量num_ite(num_ite∈{1,2,…,10})，初始化num_ite＝1，num_ite代表当前正在进行第num_ite次迭代；

步骤4：初始化存放各胞腔中训练矢量之和的变量V(N,n)为全零矩阵，以及更新存放各胞腔内训练矢量个数的变量U(1,N)为全零矢量；

步骤5：输入上一次迭代产生的第t部分训练矢量码书和对应的训练矢量码字索引，并计算训练矢量码书中所有码字的方差并保存，对于码字Y_j，其方差V_Y的计算方式为其中变量i表示码字的第i维分量，变量n表示码书中行矢量的维数，即码书矩阵的列数(如果是首次迭代，则输入第t部分训练矢量初始码书和对应的训练矢量初始码字索引，并计算相应的初始码字的方差并保存)；

步骤6：引入变量h(h∈{1,2,…,65536})，初始化h＝1，h代表B_tV中第h个行矢量，即第h行；

步骤7：按顺序输入B_tV中第h个行矢量作为当前要编码的训练矢量X，按照上一次迭代后X对应的码字索引b＝I(1,h)找到相应的码字Y_b(如果是首次迭代则按照初始码字索引b＝I_tini(1,h)找到初始码书中相应的码字Y_b)，计算X和Y_b之间的欧式距离作为当前最小失真变量i代表矢量的第i维分量；

步骤8：判断码书中的码字是否都已经与当前训练矢量X比较完成，若已经比较完成，则进入训练阶段步骤12；否则，按照以Y_b为中心上下搜索的顺序，输入下一个要判断的码字Y_j；

步骤9：计算X第一维分量X₁与码字Y_j第一维分量Y₁的差的平方D₁＝(X₁-Y_j1)²，判断不等式D₁＝(X₁-Y_j1)²≥D_min是否成立，若成立则进一步判断X₁与Y_j1的关系，若X₁≥Y_j1，进入步骤a)，若X₁≤Y_j1时，进入步骤b)；若不等式D₁＝(X₁-Y_j1)²≥D_min不成立，则进入步骤10进行下一步判断；

a)排除码字Y_i(i＝1,…,j)，转入步骤8；

b)排除码字Y_i(i＝j,…,N)，转入步骤8(j为当前搜索码字的下标，N表示码书尺寸)；

步骤10：计算D₂＝D₁+(V_X-V_Y)²，判断不等式D₂≥D_min是否成立，若成立则排除码字Y_j，转入步骤11，否则进入步骤12进行下一步判断；

步骤11：

2.1初始化变量q＝2，q∈{2,3,…,n}，n是当前处理的训练矢量矩阵B_tV的列数；

2.2计算当前训练矢量X与当前判断码字Y_j的累积失真

2.3判断是否成立，若成立则排除码字Y_j，转入步8，否则转入d)；

2.4q增加1即q＝q+1，判断q>n是否成立，若成立则转入步骤11中的e)；否则转入步骤11中的b)继续判断；

2.5把当前码字Y_j作为训练矢量X对应的新的失真最小的码字，即最匹配码字，并把X对应的码字索引I(1,h)更新为I(1,h)＝j，把X对应的最小失真D_min更新为当前码字Y_j和X之间的欧式距离然后转入步骤8；

步骤12：把当前训练矢量X划分到对应的胞腔中，即把X对应胞腔的训练矢量之和变量更新为V(I(1,h),1:n)＝V(I(1,h),1:n)+X，把X对应胞腔的训练矢量个数变量更新为U(1,I(1,h))＝U(1,I(1,h))+1；

步骤13：h增加1即h＝h+1，判断h>65536是否成立，若成立则转入步骤14；否则转入训练阶段步骤7，继续输入下一个训练矢量进行编码；

步骤14：本轮迭代结束，更新码书，以各胞腔的质心来代替原来胞腔对应的码字。

步骤15：num_ite增加1即num_ite＝num_ite+1，判断num_ite>ite是否成立，若成立则转入训练阶段步骤16，否则转入训练阶段步骤5，继续下一次迭代；

步骤16：t增加1，即t＝t+1，判断t>3是否成立，若成立则转入训练阶段步骤17，否则转入训练阶段步骤2继续训练下一部分训练矢量；

步骤17：打包各部分训练矢量码书和对应的训练矢量码字索引保存。

码字重组阶段：

步骤1：引入变量c(c∈{1,2,3})，初始化c＝1，c代表第c部分矢量数据；

步骤2：输入第c部分训练矢量数据B_cV(B_cV∈{B_1V,B_2V,B_3V}、排序索引I_cV、余下矢量数据B_cR、以及对应的最后一次迭代产生的训练矢量码书Y_cdie和训练矢量码字索引I_cdie；

步骤3：按照I_cV调整I_cdie，得到新的索引，使新的索引与最后一次迭代产生的码书Y_cdie和原来的训练矢量B_cT对应，并把新的索引作为B_cT对应的第c部分最终训练矢量码字索引I_cT，把最后一次迭代产生的训练矢量码书Y_cdie直接作为B_cT对应的第c部分最终训练矢量码书Y_cT；

步骤4：按照最终训练矢量码字索引I_cT，把对应的余下矢量数据B_cR中的行矢量依次分配到对应的胞腔中，以各胞腔质心作为对应的最终余下矢量码字Y_cR；

步骤5：将n维的最终训练矢量码字Y_cT和2ⁿ-n维的最终余下矢量码字Y_cR按码字序号合并重组成2ⁿ维的最终完整码字Y_cZ，即按最终训练矢量码字在前最终余下矢量码字在后的顺序，把相同序号的最终训练矢量码字和最终余下矢量码字连接起来:；

步骤6：将Y_cZ的各列按照对应的分量离散度排序索引I_cd重新排序，使最终完整码字各分量恢复到原来的顺序，得到Hadamard变换域的最终完整码字Y_cH；

步骤7：将Y_cH中的码字全部作Hadamard反变换，得到空域最终完整码字Y_cN；

步骤8：把最终训练矢量码字索引I_cT直接作为空域最终完整码字索引I_cN，打包空域最终完整码字Y_cN和对应的空域最终完整码字索引I_cN，进行存储或传输，完成第c部分矢量数据的压缩；

步骤9：c增加1，即c＝c+1，判断c>3是否成立，如果成立则整个高光谱图像的压缩完成，如果不成立则转入码字重组阶段步骤2，继续下一部分矢量的码字重组。

在接收端只需按照第一部分空域最终完整码字Y_1N和对应的空域最终完整码字索引I_2N恢复第一部分矢量数据，按照第二部分空域最终完整码字Y_2N和对应的空域最终完整码字索引I_3N恢复第二部分矢量数据，按照第三部分空域最终完整码字Y_3N和对应的空域最终完整码字索引I_3N恢复第三部分矢量数据，再将恢复的各部分矢量数据按顺序组合起来，即相同序号的矢量按第一部分、第二部分、第三部分的顺序连接起来，组成维数为三部分矢量维数之和的矢量，即可得到一个恢复的矢量数据矩阵B_recon，经过维数变换即可得到恢复的高光谱图像。

高光谱图像压缩的性能主要用峰值信噪比(PSNR，单位：dB)、压缩比(CR)、以及计算复杂度三个指标来进行评价，这3个指标的具体计算方式如下：

①PSNR计算公式：

其中，Peaksignal表示高光谱图像矩阵B中的最大像素值，MSE表示整个高光谱图像的平均量化误差；

MSE计算公式：

B是原始矢量矩阵，B_recon是恢复后的矢量矩阵，B和B_recon都是65536行224列的。

②CR计算公式：

其中，N为码书尺寸，224为原始矢量矩阵B中的矢量维数，32、64、128分别为第一部分、第二部分、第三部分矢量的维数，16代表每个像素值用2个字节即16bit表示，65536表示矢量个数。

③计算复杂度计算公式：

其中，Num₁表示比较计算总次数(com)，Num₂表示乘法计算总次数(×)，Num₃表示加法计算总次数(±)，Num₄表示开平方计算总次数(sqrt)，ite表示矢量量化过程的迭代次数，65536表示矢量的总的个数。

以Lunar Lake和Jasper的第一场景的部分数据作为仿真测试数据，设置迭代次数10次，其他参数设置相同，通过3个评价指标对LBG方法、矢量维数分割量化的超光谱图像压缩方法(图2为矢量维数分割量化的超光谱图像压缩方法流程图)和本发明方案进行了仿真对比，实验结果如下：表1列出了Lunar Lake图像经过3种方法处理后的峰值信噪比(PSNR)和压缩比(CR)比较；表2列出了Jasper图像经过3种方法处理后的峰值信噪比(PSNR)和压缩比(CR)比较；表3列出了Lunar Lake图像经过3种方法处理后的运算复杂度，表4列出了Jasper图像经过3种方法处理后的运算复杂度。

表1.Lunar Lake图像经过3种方法处理后的峰值信噪比(PSNR)和压缩比(CR)比较

表2.Jasper图像经过3种方法处理后的峰值信噪比(PSNR)和压缩比(CR)比较

表3.Lunar Lake图像经过3种方法处理后的运算复杂度

表4.Jasper图像经过3种方法处理后的运算复杂度

从实验结果可以看出：本文提出的压缩方案相对于LBG方法，压缩比稍有降低，但是峰值信噪比有较大提高，平均峰值信噪比提高约2.7，最高峰值信噪比提高2.98，而且计算量得到大量减少，集中主要计算量的乘法计算次数平均不到LBG方法的0.5％，加法次数平均不到LBG方法的0.1％；同时，本文提出的发明方案相比于分割方法，压缩比相同，峰值信噪比略微降低，平均降低仅为0.12，最小降低仅为0.06，但是计算复杂度大幅度降低，尤其是集中了主要计算量的乘法计算仅为分割方法的7.7％，加法计算仅为分割方法的11.4％。

此外，在空间复杂度方面，LBG方法与分割方法都是将高光谱图像的全部波段输入参与训练，训练矢量的维数等于图像的总波段数，而本文提出的压缩方案，仅仅提取少数的log₂(k)(k为原始波段数)个特征波段参与训练，训练矢量的维数仅为log₂(k)个。例如本文实验中使用的224个波段的高光谱图像，LBG方法与分割方法输入的参与训练的训练矢量维数为224，而本文提出的压缩方案将高光谱图像经变换分割成32维、64维、128维的三部分矢量后，仅需输入log₂(32)+log₂(64)+log₂(128)＝18维矢量分量作为训练矢量，压缩中所需的内存空间为LBG方法与分割方法的18/224≈0.08。所以本文提出的压缩方案空间复杂度远小于LBG方法与分割方法，即使在处理机内存受限的情况下，也可以高效率实现较高质量的高光谱图像压缩。

总之，本文提出的压缩方案，与LBG方法比较，在压缩比略有降低的情况下，提高了图像恢复质量，同时大量减少了计算量；与分割方法比较，在压缩比相同，基本不影响图像恢复质量的条件下，大幅度降低了计算复杂度。而且，本文提出的压缩方案实现过程的空间复杂度远远小于LBG方法和分割方法，即使在处理机内存受限时也可以高效完成高光谱图像较高质量的压缩。所以本文提出的压缩方案能在保证较好的图像恢复质量的前提下，高效快速地完成高光谱图像的压缩，且简单易操作，具有实际应用的价值。高光谱图像信息丰富、应用广泛，但具有波段数较多和数据量海量的特点，因此高光谱图像的压缩方法研究一直备受关注，本文提出的压缩方案具有较高压缩比、较好图像恢复质量、低计算复杂度的特点，可以作为高光谱图像的一种快速压缩方案。

Claims (7)

1.一种快速压缩高光谱信号的矢量量化方法，其特征在于，读取高光谱图像数据，截取每个波段中的相同空间位置的一段图像块作为输入信源，根据截取的图像块尺寸和光谱波段数构建三维数据矩阵，以同一空间位置对应像素各波段像素值作为一个矢量数据，即把截取的图像块每一个像元对应的光谱矢量对应一个矢量数据，把三维数据转化为二维矩阵数据；截取转化为2维的矢量数据并分割成M部分，即按照2ⁿ准则将所有矢量数据进行分割，使每个分割出的矢量部分中矢量的维数都要满足2ⁿ，对分割出的各部分矢量数据分别进行哈达玛Hadamard变换；对Hadamard变换域的各部分矢量数据分别作离散度排序得到矢量数据(B_1D、B_2D，…，B_MD)；提取矢量数据行矢量的前n维分量组成训练矢量集(B_1T、B_2T，…，B_MT)，剩下的k-n维分量组成对应部分的余下矢量数据(B_1R、B_2R，…，B_MR)，其中，n＝log₂(k)，k为该部分矢量数据行矢量的维数；对训练矢量集(B_1T、B_2T，…，B_MT)分别进行训练，得到矢量量化最后一次迭代产生的各部分训练矢量码书(Y_1die、Y_2die，…，Y_Mdie)和对应的训练矢量码字索引(I_1die、I_2die，…，I_Mdie)；计算余下矢量数据的码书，将训练矢量码书和余下矢量数据的码书重组连接，通过离散度反排序和hadamard反变换恢复出空域完整码书；打包空域完整码书和对应的空域完整码字索引，存储或传输，其中，B_MD代表第M部分离散度排序后的矢量数据，B_MT代表第M部分训练矢量集，B_MR代表第M部分余下矢量数据，Y_Mdie、I_Mdie分别代表最后一次迭代产生的第M部分训练矢量码书和第M部分训练矢量码字索引。

2.根据权利要求1所述的方法，其特征在于，分割矢量数据具体包括：读取3维高光谱图像数据，截取部分相同空间位置的像素值转化为2维矢量数据并按照2ⁿ的原则将矢量数据分割为M个小部分，分割原则为：每部分的矢量维数都为2ⁿ，维数不为2ⁿ的部分通过末尾补零使其维数为2ⁿ，各部分矢量维数之和等于原矢量维数。

3.根据权利要求1所述的方法，其特征在于，所述离散度排序具体为：将行矢量某一维分量的最大值减去最小值得到差值获取行矢量该维分量的离散度，将离散度值分别按照降序排序，得到对应的排序索引作为各部分矢量数据的离散度排序索引，按照离散度排序索引将对应的各部分矢量数据的行矢量分量重排序，得到离散度排序后的各部分矢量数据。

4.根据权利要求1所述的方法，其特征在于，对训练矢量集(B_1T、B_2T，…，B_MT)进行训练具体为：按照行矢量第一维分量值的大小对各部分训练矢量进行升序排序，排序索引作为各部分矢量排序索引，按照各部分矢量排序索引重新排序各部分训练矢量，将排序后的训练矢量按照码书大小进行平均分组，依次选取每组第一个矢量组成初始码书；根据初始码书和初始码字索引搜索训练矢量对应的最佳匹配码字，并将训练矢量划分到对应的胞腔中，记录对应的码字索引，直到所有矢量训练完成，更新码书，以各胞腔的质心代替原来胞腔对应的码字。

5.根据权利要求1所述的方法，其特征在于，计算余下矢量数据的码书具体为：按照训练矢量码字索引(I_1T、I_2T，…，I_MT)分别把对应部分余下矢量数据(B_1R、B_2R，…，B_MR)中的行矢量依次分配到对应的胞腔中，以各胞腔质心作为对应胞腔的余下矢量数据的余下矢量码字Y_cR，根据余下矢量码字构建余下矢量的数据码书。

6.根据权利要求4所述的方法，其特征在于，搜索最佳匹配码字具体包括如下步骤：①根据训练矢量X的第1维分量X₁、当前码字Y_j的第1维分量Y_j1，根据公式：D₁＝(X₁-Y_j1)²获取第一维分量差平方D₁，如果D₁≥D_min，且X₁≥Y_j1，排除码字Y_i(i＝1,…,j)，如果D₁≥D_min，且X₁<Y_j1，排除码字Y_i(i＝j,…,N)，否则执行步骤②；②根据公式：计算训练矢量X的方差V_X、码字Y_j的方差V_Y，根据公式D₂＝D₁+(V_X-V_Y)²确定变量D₂，若不等式D₂≥D_min成立则排除码字Y_j，③根据公式计算训练矢量X与当前码字Y_j的累积失真，如D_q(X,Y_j)≥D_min成立，排除码字Y_j，否则把训练矢量X划分到对应的胞腔中，其中，X_i为训练矢量X的第i维分量，Y_ji为码字Y_j的第i维分量，n为训练矢量的维数，D_min为当前最小失真，N为码书尺寸。

7.根据权利要求1-5其中之一所述的方法，其特征在于，将训练矢量码书和余下矢量数据的码书重组连接具体为：按最终训练矢量码字在前、最终余下矢量码字在后的顺序，把相同序号的最终训练矢量码字和最终余下矢量码字连接起来，合并重组成2ⁿ维的完整码书。