CN103442236A - 一种多级和分维矢量量化的遥感信号压缩编码方法 - Google Patents

Abstract

本发明提出了一种基于多级分维矢量量化的高光谱信号编码方法。对第一级矢量量化的输入矢量首先进行去均值操作，构造第一级矢量量化的初始码书，并生成第一级矢量量化的最终码书和编码；由第一级矢量量化的重构图像和原去均值图像，构成矢量进行分维处理，分维后每一个部分的数据作为第二级矢量量化的输入矢量，构造第二级矢量量化的初始码书，同时利用快速搜索算法，实施快速聚类，生成每一个部分数据的第二级矢量量化的最终码书和索引。本发明结合去均值、矢量分维和多级矢量量化达到减小量化误差的目标，不仅能够提高压缩比和图像恢复质量，而且可以大幅降低算法的计算复杂度，达到对高光谱遥感信号快速压缩编码的目的。

Description

一种多级和分维矢量量化的遥感信号压缩编码方法

技术领域

本发明属于高光谱遥感图像处理领域，具体涉及一种结合多级和分维矢量量化技术的高光谱遥感图像数据压缩编码方法。

背景技术

一切物体均具有反射或者辐射不同波长电磁波的特性，通过识别电磁波而达到识别物体及其存在环境的技术称为遥感。它是20世纪60年代发展起来的一项综合对地观测技术，能够不直接接触目标、区域或现象而获得相关数据，并加以分析以获得所需要的信息。遥感技术是以电磁辐射理论为基础，集合了电磁波理论、光谱学和色度学、物理与几何光学、地理学、地质学、大气科学以及海洋科学等众多学科。遥感作为人类获取信息的重要手段，经过几十年的发展，已经在国民经济各个邻域有着广泛的应用，并且成为当今最为活跃的科技领域之一。

高光谱遥感（Hyperspectral remote sensing）又称超谱遥感，目前应用最多的是美国NASA喷气实验室研制的机载可见光成像光谱仪（AVIRIS）。高光谱遥感图像是指利用高光谱传感器对同一对象在多个窄光谱波段范围上获得的图像，它记录了成像地物在多个光谱波段的响应特性，提供了关于地物更细致的光谱信息，但同时其数据量急剧增加，给机载和星载数据传输和存储带来困难。典型AVIRIS图像的大小为614×512，成像波段数为224，如果每像素灰度值以两个字节存储，一组图像的数据量则约为140M字节。因此，随着高光谱遥感技术的发展，海量的高光谱遥感数据的传输、存储和管理带来的困难也日益加剧。为了提高高光谱传感器的工作效率，在有限的通信带宽实现高光谱遥感图像的高速传输，使各种机载/星载超谱传感器获得的高光谱遥感图像最大限度地得到利用，研究和开发高性能的高光谱遥感图像压缩技术具有非常重要的意义。对高光谱遥感图像进行压缩可以：（1）较快地传输各种信源，降低信道占用费用。（2）在现有的通信干线上开通更多的并行业务。（3）降低发射机的功率。（4）紧缩数据存储容量，降低存储费用。因此，高光谱遥感图像压缩的算法研究具有重要应用价值。

高光谱遥感图像的压缩要考虑去除两类相关性——谱间相关性和空间相关性。谱间相关性部分是由于传感器的频谱交叠，部分是由于地物的反射频谱覆盖了较大的频谱范围。谱间相关性又可以进一步分为两类，在统计意义上存在的图像像素数值上的相关性称为统计相关性。由于各个波段图像对应同一视场，图像内容具有相同的几何结构的相关性称为结构相关性。空间相关性主要是由于同一地物具有相同的反射特性。由于矢量量化可以通过聚类的方法去掉高光谱遥感图像的冗余度，它能有效地利用矢量中各分量间的4种相互关联的性质：线性依赖性、非线性依赖性、概率密度函数的形状以及矢量维数，并且矢量量化具有压缩比大、编解码简单和失真较小的特点。因此矢量量化作为一种对高光谱遥感图像压缩效果显著，图像恢复质量高的压缩算法，得到广泛的应用。申请号为201110272304.X的专利申请公开为一种超谱信号的快速压缩编码方法及图像压缩方法，但该方法在压缩比、图像恢复质量（峰值信噪比）和计算复杂度上仍有待改善。

发明内容

本发明针对现有高光谱遥感图像压缩采用无损压缩方法图像的压缩比低和有损压缩图像的恢复质量差的缺点，提出了一种基于多级分维矢量量化的高光谱遥感图像的快速压缩编码方法，在保证码书尺寸不变的情况下，提高了高光谱遥感图像的压缩比和恢复质量，同时也降低了算法的计算复杂度。

一种多级和分维矢量量化的高光谱遥感信号快速编码方法，包括：读取高光谱数据源，提取高光谱数据各行矢量的均值并保存，对高光谱的行矢量进行去均值操作获得去均值图像，构造第一级矢量量化的初始码书，获取第一级编码码书和编码；用第一级码书及编码索引重构的图像，与去均值图像相减得到差值图像，对差值图像矢量进行分维处理，对分维后各部分分别构造对应的第二级矢量量化的初始码书作为第二级矢量量化的输入数据，训练该初始码书得到各部分的第二级矢量量化的码书和编码索引；由行矢量均值、第一级码书和编码索引、以及分维的各部分第二级矢量量化码书和编码索引构成高光谱数据源的压缩数据。

其主要包括以下几个阶段：

第一级矢量量化初始化阶段：读取高光谱数据源，将高光谱数据源矩阵的每一个行矢量与其均值做差，对生成的去均值矩阵进行哈达玛变换，将得到的矢量数据按照行矢量均值的大小进行升序排列，获得排序后索引值和排序后的矢量矩阵，然后对排序后的矢量平均分组，依次选取每组首矢量得到第一级矢量量化的初始码书，设置第一级矢量量化过程的最大迭代次数或失真阈值。

第一级矢量量化训练阶段：根据欧式距离确定训练矢量的最佳匹配码字，依次将训练矢量划分到它所对应最佳匹配码字的胞腔中，直到划分完所有的训练矢量，最后用胞腔的质心来更新当前最佳匹配码字。直至累计迭代次数达到设定的最大迭代次数或最近两次平均失真的相对误差值小于设置的失真阈值，结束训练并打包第一级矢量量化的最终码书以及所有矢量最佳匹配码字的索引和行矢量均值；否则，继续训练。

第二级矢量量化初始化阶段：用第一级矢量量化生成的码书以及所有矢量对应的最佳匹配码字的索引构造哈达玛域的恢复图像，利用矢量去均值后的图像与所得到的恢复图像生成差值图像，对差值图像按照2的整数次幂进行分维处理，生成三个部分。以第一部分为例，对第一部分的矢量进行哈达玛变换，将变换后矢量数据按照第一维分量值进行升序排列，获得排序索引，然后平均分组，依次选取每组首矢量得到第一部分的第二级矢量量化的初始码书，并设置第一部分的第二级矢量量化过程的迭代次数或失真阈值。

第二级矢量量化训练阶段：采用快速算法搜索训练矢量的最佳匹配码字，依次将训练矢量划分到它所对应的胞腔中，直到划分完所有的训练矢量，最后用胞腔的质心来更新当前最佳匹配码字。当累计迭代次数达到设定的最大迭代次数或最近两次平均失真的相对误差值小于设置的失真阈值，结束训练并打包每一部分的第二级矢量量化的最终码书以及所有矢量最佳匹配码字的索引值；否则，继续训练。

生成初始码书具体为：第一级矢量量化初始码书生成，将高光谱数据源矩阵的每一个行矢量与其均值相减，生成去均值后的矢量，对去均值后的矩阵矢量矩阵进行扩展，形成列数（维数）为2的整数次幂的行矢量，然后对该变换矩阵进行Hadamard变换，并按照原来行矢量的均值升序进行排序，根据排序索引对行排序后的矩阵E按照行平均分N₁组（N₁为第一级矢量量化码书的尺寸），依次选择每组第一个矢量作为该组的初始码字生成初始码书Y；第二级矢量量化初始码书生成，对生成的差值图像按照2的整数次幂进行分维处理生成FC(i)（其中i表示分维的各部分的标记），对FC(i)矩阵进行Hadamard变换并按照变换后矢量的第一维按照升序进行排序列，根据排序索引对行排序后的矩阵FC(i)_E按照行平均分N₂组（N₂为第二级矢量量化码书的尺寸），依次选择每组第一个矢量作为该组的初始码字生成第二级矢量量化的初始码书FC(i)_Y。

由上述过程可见，矢量量化编码只需要传输最终码书以及各个训练矢量对应于最终码书中码字的索引，解码过程只需要在接收端利用接收到的索引值在码书中查找与其所对应的码字来代替训练矢量，就可以恢复出高光谱遥感图像。由于矢量量化编码过程传输的数据量远远小于原始高光谱遥感图像的数据量，因此，通过矢量量化编码的方法可以达到压缩高光谱遥感图像的目的。同时，由本发明产生的多级码书和索引，支持更加灵活的数据处理、传输和存储策略，如基于渐进传输的信道环境和数据的渐进式呈现。

附图说明：

图1是传统的LBG算法流程框图；

图2是本发明方案流程框图。

具体实施方式

下文中所使用变量的意义如下：x表示空域中的矢量；X表示哈达玛域中的矢量；y表是空域中的码字；Y表示哈达玛域中的码字；ite表示迭代次数；D_min表示当前最小失真；N表示码书尺寸；I表示索引矩阵；E表示哈达玛域矢量排序后的训练矩阵矢量；V表示胞腔矢量存储矩阵；CZ表示生成的差值图像；FC(i)表示第二级矢量量化分维后各个部分的前缀，其中i表示第i部分；PSNR表示峰值信噪比；CR表示压缩比；Complex表示计算复杂度。

以下使用具体事例和附图对本发明作进一步说明。本发明提出的基于多级和分维矢量量化的高光谱遥感信号快速编码方法，主要包括：

从高光谱遥感图像中获取输入矢量，将输入矩阵的行矢量与该矢量的均值进行差值操作，对差值后矩阵行矢量变维，并进行哈达玛变换，将变换后的矢量数据按照行矢量均值进行升序排列，获得排序索引，然后利用平均分组选取每组首矢量的方法得到第一级矢量量化的初始码书，并设置第一级矢量量化过程的迭代次数。

快速搜索训练矢量的最佳匹配码字，依次将训练矢量划分到它们所对应最佳匹配码字的胞腔中，直到划分完所有的训练矢量，最后用胞腔的质心来更新当前码字。每轮迭代之前预先判断累计迭代次数是否达到设定的最大迭代次数ite₁，若满足，就结束训练并打包第一级矢量量化的行矢量均值、最终码书和所有矢量最佳匹配码字的索引；否则，继续训练。

用第一级矢量量化的最终码书、每个输入矢量对应匹配码字的索引、哈达玛反变换和变维（去除零矩阵）重构恢复图像，然后用去均值图像与恢复图像生成差值图像；对差值图像矢量的维数按照2的整数次幂进行分维处理，然后将每一部分的矢量进行哈达玛变换。

以第一部分为例，对第一部分变换后矢量数据按照第一维分量值进行升序排列，获得排序索引，然后利用平均分组选取每组首矢量得到初始码字，形成第二级矢量量化第一部分矩阵的初始码书，并设置第一部分第二级矢量量化过程的迭代次数。搜索训练矢量的最佳匹配码字，依次将训练矢量划分到它所对应最佳匹配码字的胞腔中，直到划分完所有的训练矢量，最后用胞腔的质心来更新当前码字。判断累计迭代次数是否达到设定的最大迭代次数，若满足，结束训练并打包第一部分第二级矢量量化的最终码书以及所有矢量最佳匹配码字的索引；否则，继续训练。

对图像数据矢量分维后的其余部分重复上面对第一部分的操作，相应获得各部分的编码码书和编码索引。

具体步骤如下：⑴读取高光谱遥感图像，根据图像尺寸和波段数构建3维矩阵，并将其转化为2维矩阵B；⑵将B的每一个行矢量与其该行矢量的均值相减，生成去均值后的矩阵C；⑶判断2维矩阵C的行维数是否为2的整数次幂，如果不是，则在C的行向量后添加零矩阵，使C的行向量的维数转化成与其最接近的2的整数次幂的维数，这样产生新的增维后的矩阵C’；⑷对2维矩阵C’做行矢量的哈达玛变换，将变换后的矩阵的行按照原行矢量均值进行升序排列，获得矢量排序后的矩阵E（第一级矢量量化训练矢量集）并记录此时排序索引I₁；⑸对矩阵E按行平均分为N₁（第一级矢量量化的码书尺寸）组，依次选择每组首矢量生成初始码书，根据变换质量要求设置迭代次数ite₁或失真阈值ε₁；⑹设定第一级矢量量化码书的训练矢量的胞腔索引矩阵I，将初始码书中训练矢量归属到相应的胞腔中；⑺利用快速算法搜索最佳匹配码字，依次将训练矢量划分到其最佳匹配码字所对应的胞腔中，直到训练完所有的矢量，用胞腔的质心来更新码字。在达到设定的最大迭代数或失真阈值时，输出第一级矢量量化过程的行矢量均值、最终码书Y和训练矢量的聚类索引（胞腔号）I；⑻利用第一级矢量量化生成的最终码书Y和聚类索引I重构图像，然后对其进行哈达玛反变换和变维（去除零矩阵）生成恢复图像；⑼利用图像C和恢复图像生成差值图像CZ；⑽按照2的整数次幂对差值图像的矢量进行分维处理，不满足的进行补零操作，生成矩阵FC(i)，其中i表示分维后各部分的标记；⑾对2维矩阵FC(i)进行哈达玛变换，将变换后的行矢量按照第一维分量值进行升序排列，获得矢量排序后的矩阵FC(i)_E（第二级矢量量化训练矢量集）并记录此时排序索引FC(i)_I₁；⑿对矩阵FC(i)_E按行平均分为N₂（第二级矢量量化的码书尺寸）组，依次选择每组首矢量生成初始码书，设置迭代次数ite₂或失真阈值ε₂；⒀设定第二级矢量量化的训练矢量的胞腔索引矩阵FC(i)_I，将初始码书中训练矢量归属到相应的胞腔中；⒁利用快速算法搜索最佳匹配码字，然后将训练矢量划分到最佳匹配码字所对应的胞腔，直到训练完所有的矢量，本轮迭代结束时以胞腔的质心更新码字。在达到设定的最大迭代数或失真阈值时，输出第二级矢量量化设计最终码书FC(i)_Y以及训练矢量对应的聚类索引（胞腔号）FC(i)_I。

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

在MATLAB7.9软件平台下结合附图对本发明方案的实例进行详细说明。图2为本发明方法流程框图，使用典型的AVIRIS高光谱遥感图像（大小为512×614×224），为了处理方便，在每个波段中截取大小为256×256的图像块，其每个波段中像元数据以两个字节的带符号整数表示，本方案的具体实施步骤如下：

第一级矢量量化初始化阶段：

步骤1：读取高光谱数据，并将其转化为二维矩阵B(65536,224)

步骤2：对矩阵B的行矢量都进行去均值操作，即将每一个行矢量都与其行矢量均值相减，生成新的矢量，记为矩阵C；

步骤3：在矩阵C的每一行矢量后添加32维的零矢量，获得矩阵C’(65536,256)；

步骤4：对C’进行Hadamard变换得到变换矩阵D(65536,256)，将矩阵D按照原行矢量均值的大小重新进行升序排序，获得行排序后的矩阵E，记录这时的排序索引I₁，然后对矩阵E按照行平均分N₁组（N₁为第一级矢量量化码书的尺寸），依次选择每组第一个矢量作为该组的初始码字，生成初始码书Y；

步骤5：置初始编码索引矩阵为I(1,65536)，根据第一级对于训练矢量集E的划分，将训练矢量按照初始划分分配到相应的胞腔中，获得训练矢量的初始编码（胞腔号）；

步骤6：设置初始化迭代次数ite₁，根据编码质量要求由用户确定ite₁的值；第一级矢量量化训练阶段：

步骤1：计算E中所有训练矢量和Y中码字矢量的方差，并保存在初始码书中；

步骤2：初始化存放训练矢量的胞腔V以及存放此胞腔训练矢量的个数U；

步骤3：对于上述E矩阵的每一个行矢量，即每一个训练矢量（记为X），从上一次迭代后其归属的胞腔开始搜索，计算当前训练矢量X与码字Y_I(i)(I(i)是X上一次的编码索引值，对应的Y_I(i)就是其上一次X的编码码字）的欧式距离，并将此欧式距离作为当前最小失真D_min；

步骤4：计算D₁=(V_X-V_j)²（其中，

M是矢量的维数），即计算训练矢量X与码字Y_j方差差值平方和，如果满足D₁≥D_min，则排除码字Y_j，转入步骤6；否则，转入下一步进行判断；

步骤5：根据部分失真公式

其中q为不大于M(维数）的整数。计算训练矢量X与码字Y_j的部分失真，如果满足D_q≥D_min，则排除码字Y_j，转入步骤6；否则，继续累计部分失真，直到q的值为M（矢量的维数），此时如果D(X,Y_j)=D_q<D_min，则令D_min=D_q，更新训练矢量的编码索引值I(i)为当前最近的码字的下标值j，即I(i)=j，进入下一步；

步骤6：返回步骤4搜索下一个码字的最佳匹配，直到搜索完所有的剩余码字，获得当前训练矢量X的最佳匹配码字以及其所对应编码索引I(i)。此时，将训练矢量X归为第I(i)个胞腔中，置第I(i)个胞腔矢量个数计数器U(I(i))→U(I(i))+1和矢量存储器V(I(i))→V(I(i))+X，转入下一步；

步骤7：输入下一个训练矢量，返回步骤2继续训练，直到所有的训练矢量被训练完为止，求出这N₁个胞腔的质心（质心计算公式∑V(I(i))/U(I(i))，即胞腔矢量的和与胞腔矢量个数之商），并将这些质心替换原有码书中的码字，形成新的码书Y并置迭代次数ite₁→ite₁+1；

步骤8：根据设定的迭代次数ite₁，判断是否满足停止迭代的要求，如果满足，将所得Y作为最终码书，否则转到步骤2进行下一轮训练；

步骤9：以排序索引I₁重新排序编码索引I，获得新的编码索引I₂，这个索引正好对应排序前矩阵C的每行待编码矢量编码后的对应码字的标号，可经哈达玛反变换后，重建原数据块的编码数据块，码书Y以及编码索引I₂作为第一级矢量量化的最终压缩数据，可存储或传输；

第二级矢量量化矢量分维处理阶段：

步骤1：利用初级码书设计的最终码书Y和编码索引I₂生成哈达玛域的重构图像数据矩阵，然后对哈达玛域的恢复矩阵进行Hadamard反变换，并进行去除32维零矩阵操作，生成空域恢复矩阵HF(65536,224)；

步骤2：将空域恢复矩阵HF与原图像的去均值矩阵C(65536,224)相减，得到差值矩阵CZ(65536,224)；

步骤3：对差值矩阵按照2的整数次幂进行矢量分维，不满足的进行补零操作，本发明通过矢量分维将差值矩阵分为3层，即FC1(65536,32)、FC2(65536,64)和FC3(65536,128)，无需补零操作；

⑴对分维后第一部分FC1初始化阶段：

步骤1：对FC1进行哈达玛变换得到变换矩阵FC1_C(65536,32)，将变换矩阵FC1_C按照第一列分量的大小重新进行升序排列，记录这时的排序索引FC1_I₁，获得行排序后的矩阵FC1_E，然后对矩阵FC1_E按行平均分为N₂组（N₂为第二级矢量量化的码书尺寸），依次选择每组第一个矢量作为该组的初始码字，生成初始码书FC1_Y；

步骤2：置编码索引矩阵FC1_I(1,65536)，依据划分，为每个训练矢量分配索引值（胞腔号），对应填入索引矩阵FC1_I(1,65536)，得初始编码索引；

步骤3：设置初始化迭代次数ite₂，根据编码质量要求由用户确定ite₂的值；

对分维后第一部分FC1化训练阶段：

步骤1：计算所有训练矢量与其初始码书中码字的方差，并保存在第二级矢量量化的第一部分设计的初始码书中；

步骤2：初始化存放差值训练矢量的胞腔FC1_V以及此胞腔存放的训练矢量的个数FC1_U；

步骤3：对于上述FC1_E矩阵的每一个行矢量，即每一个训练矢量（记为FC1_X），从上一次迭代后其归属的胞腔开始搜索，计算当前训练矢量FC1_X与码字FC1_Y_{FC1_I(i)}(FC1_I(i)是FC1_X上一次的编码索引值，对应的FC1_Y_{FC1_I(i)}就是其上一次FC1_X所对应的编码码字）的欧式距离，并将此欧式距离作为当前最小失真FC1_D_min；

步骤4：以码字FC1_Y_{FC1_I(i)}为中心，采用上下查找法，对于当前码字FC1_Y_j，如果满足FC1_D₁=(FC1_X₁-FC1_Y_j1)²≥FC1_D_min，即矢量FC1_X的第一维分量FC1_X₁与码字FC1_Y_j的第一维分量FC1_Y_j1差值的平方大于当前最小欧式距离，就排除码字FC1_Y_j，当FC1_X₁≥FC1_Y_j1时，进入步骤a)，当FC1_X₁≤FC1_Y_j1时进入步骤b)；否则，进入步骤5；

a)排除码字FC1_Y_i,i=1,…,j，转入步骤7；

b)排除码字FC1_Y_i,i=1,…,N₂，转入步骤7；（其中，j为当前搜索

码字的标号，N₂表示第二级矢量量化码书的大小）

步骤5：调用公式FC1_D₂=FC1_D₁+(FC1_V_{FC1_X}-FC1_V_j)²其中，

${\mathrm{FC}1\_V}_{\mathrm{FC}1\_X}=\sqrt{\underset{i=2}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{FC}1\_X}_i^2},$ ${\mathrm{FC}1\_V}_j=\sqrt{\underset{i=2}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{FC}1\_Y}_{\mathrm{ji}}^2},$ M是矢量的维数，计算矢量FC1_X的方差和当前码字FC1_Y_j的方差差值平方与FC1_D₁的和值，如果满足FC1_D₂≥FC1_D_min，则排除码字FC1_Y_j，转入步骤7；否则，转入下一步进行判断；

步骤6：结合部分失真算法，根据公式

计算训练矢量FC1_X与码字FC1_Y_j的部分失真FC1_D_q，其中q为不大于M（维数）的整数。逐渐增加q的数值，一旦存在q满足FC1_D_q≥FC1_D_min，则排除码字FC1_Y_j，转入步骤7；否则，继续累计部分失真，直到q=M，此时如果存在FC1_D(FC1_X,FC1_Y_j)=FC1_D_q<FC1_D_min，则令FC1_D_min=FC1_D_q，更新训练矢量的编码索引值FC1_I(i)为当前最近的码字的标号j，即FC1_I(i)=j；

步骤7：返回步骤4进行下一个码字的匹配，直到搜索完所有的剩余码字，获得当前训练矢量FC1_X的最佳匹配码字和对应编码索引FC1_I(i)。将训练矢量FC1_X划分到第FC1_I(i)个胞腔中，置第FC1_I(i)个胞腔矢量个数计数器FC1_U(FC1_I(i))→FC1_U(FC1_I(i))+1和矢量存储器FC1_V(FC1_I(i))→FC1_V(FC1_V(FC1_I(i))+FC1_X，进入下一步；

步骤8：输入下一个训练矢量，返回步骤2继续训练，直到所有的第二级矢量量化过程的训练矢量被训练完为止，此时求出这N₂个胞腔的质心，即可采用每个胞腔内所有训练矢量的和值∑FC1_V(FC1_I(i))与这个胞腔内训练矢量的数量FC1_U(FC1_I(i))之商∑(FC1_V(FC1_I(i))/FC1_U(FC1_I(i)))，并且将这些质心替换原有码书中的码字，形成新的码书FC1_Y，然后计算更新码字后的平均量化误差，置迭代次数ite₂→ite₂+1；

步骤9：根据设定的迭代次数ite₂，判断是否满足停止迭代的要求，如果满足，所得FC1_Y作为最后的码书，否则转到步骤2进行下一轮训练；

步骤10：以排序索引FC1_I₁重新排序编码索引FC1_I，获得新的编码索引FC1_I₂，这个索引正好对应排序前矩阵FC1_C的每行待编码矢量编码后的对应码字的标号，可经哈达玛反变换后，重建原数据块的编码数据块。码书FC1_Y以及编码索引FC1_I₂作为第一部分第二级矢量量化的压缩数据，可存储或传输；

⑵对分维后第二部分FC2初始化阶段：

步骤1：对FC2进行哈达玛变换得到变换矩阵FC2_C(65536,64)，将变换矩阵FC2_C按照第一列分量的大小重新进行升序排列，记录这时的排序索引FC2_I₁，获得行排序后的矩阵FC2_E，然后对矩阵FC2_E按行平均分为N₂组（N₂为第二级矢量量化的码书尺寸），依次选择每组第一个矢量作为该组的初始码字，生成初始码书FC2_Y；

步骤2：置编码索引矩阵FC2_I(1,65536)，依据对于训练矢量的划分，分配相应索引（胞腔号），填入FC2_I(1,65536)，获得初始编码索引；

步骤3：设置初始化迭代次数ite₂，根据编码质量要求由用户确定ite₂的值；

对分维后第二部分FC2化训练阶段：

步骤1：计算所有训练矢量与其初始码书中码字的方差，并保存在第二级矢量量化第二部分设计的初始码书中；

步骤2：初始化存放差值训练矢量的胞腔FC2_V以及此胞腔存放的训练矢量的个数FC2_U；

步骤3：对于上述FC2_E矩阵的每一个行矢量，即每一个训练矢量（记为FC2_X），从上一次迭代后其归属的胞腔开始搜索，计算当前训练矢量FC2_X与码字FC2_Y_{FC2_I(i)}(FC2_I(i)是FC2_X上一次的编码索引值，对应的FC2_Y_{FC2_I(i)}就是其上一次FC2_X所对应的编码码字）的欧式距离，并将此欧式距离作为当前最小失真FC2_D_min；

步骤4：以码字FC2_Y_{FC2_I(i)}为中心，采用上下查找法，对于当前码字FC2_Y_j，如果满足FC2_D₁=(FC2_X₁-FC2_Y_j1)²≥FC2_D_min，即矢量FC2_X的第一维分量FC2_X₁与码字FC2_Y_j的第一维分量FC2_Y_j1差值的平方大于当前最小欧式距离，就排除码字FC2_Y_j，当FC2_X₁≥FC2_Y_j1时，进入步骤c)，当FC2_X₁≤FC2_Y_j1时进入步骤d)；否则，进入步骤5；

c)排除码字FC2_Y_i,i=1,…,j，转入步骤7；

d)排除码字FC2_Y_i,i=1,…,N₂，转入步骤7；（其中，j为当前搜索

码字的标号，N₂表示第二级矢量量化码书的大小）

步骤5：调用公式FC2_D₂=FC2_D₁+(FC2_V_{FC1_X}-FC2_V_j)²其中，

$\mathrm{FC}2\_V_{\mathrm{FC}2\_X}=\sqrt{\underset{i=2}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{FC}2\_X}_i^2},$ ${\mathrm{FC}2\_V}_j=\sqrt{\underset{i=2}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{FC}2\_Y}_{\mathrm{ji}}^2},$ M是矢量的维数，计算矢量FC2_X的方差和当前码字FC2_Y_j的方差差值平方与FC2_D₁的和值，如果满足FC2_D₂≥FC2_D_min，则排除码字FC2_Y_j，转入步骤7；否则，转入下一步进行判断；

步骤6：结合部分失真算法，根据公式

计算训练矢量FC2_X与码字FC2_Y_j的部分失真FC2_D_q，其中q为不大于M（维数）的整数。逐渐增加q的数值，一旦存在q满足FC2_D_q≥FC2_D_min，则排除码字FC2_Y_j，转入步骤7；否则，继续累计部分失真，直到q=M，此时如果存在FC2_D(FC2_X,FC2_Y_j)=FC2_D_q<FC2_D_min，则令FC2_D_min=FC2_D_q，更新训练矢量的编码索引值FC2_I(i)为当前最近的码字的标号j，即FC2_I(i)=j；

步骤7：返回步骤4进行下一个码字的匹配，直到搜索完所有的剩余码字，获得当前训练矢量FC2_X的最佳匹配码字和对应编码索引FC2_I(i)。将训练矢量FC2_X划分到第FC2_I(i)个胞腔中，置第FC2_I(i)个胞腔矢量个数计数器FC2_U(FC2_I(i))→FC2_U(FC2_I(i))+1和矢量存储器FC2_V(FC2_I(i))→FC2_V(FC2_V(FC2_I(i))+FC2_X，进入下一步；

步骤8：输入下一个训练矢量，返回步骤2继续训练，直到所有的第二级矢量量化过程的训练矢量被训练完为止，此时求出这N₂个胞腔的质心，即可采用每个胞腔内所有训练矢量的和值∑FC2_V(FC2_I(i))与这个胞腔内训练矢量的数量FC2_U(FC2_I(i))之商∑(FC2_V(FC2_I(i))/FC2_U(FC2_I(i)))，并且将这些质心替换原有码书中的码字，形成新的码书FC2_Y，然后计算更新码字后的平均量化误差，置迭代次数ite₂→ite₂+1；

步骤9：根据设定的迭代次数ite₂，判断是否满足停止迭代的要求，如果满足，所得FC2_Y作为最后的码书，否则转到步骤2进行下一轮训练；

步骤10：以排序索引FC2_I₁重新排序编码索引FC2_I，获得新的编码索引FC2_I₂，这个索引正好对应排序前矩阵FC2_C的每行待编码矢量编码后的对应码字的标号，可经哈达玛反变换后，重建原数据块的编码数据块。码书FC2_Y以及编码索引FC2_I₂是第二部分第二级矢量量化的最终压缩数据；

⑶对分维后第三部分FC3初始化阶段：

步骤1：对FC3进行哈达玛变换得到变换矩阵FC3_C(65536,128)，将变换矩阵FC3_C按照第一列分量的大小重新进行升序排列，记录这时的排序索引FC3_I₁，获得行排序后的矩阵FC3_E，然后对矩阵FC3_E按行平均分为N₂组（N₂为第二级矢量量化的码书尺寸），依次选择每组第一个矢量作为该组的初始码字，生成初始码书FC3_Y；

步骤2：置编码索引矩阵FC3_I(1,65536)，依据对于训练矢量的划分，分配相应索引（胞腔号），填入FC3_I(1,65536)，获得初始编码索引；

步骤3：设置初始化迭代次数ite₂，根据编码质量要求由用户确定ite₂的值；

对分维后第三部分FC3化训练阶段：

步骤1：计算所有训练矢量与其初始码书中码字的方差，并保存在第二级矢量量化第三部分设计的初始码书中；

步骤2：初始化存放差值训练矢量的胞腔FC3_V以及此胞腔存放的训练矢量的个数FC3_U；

步骤3：对于上述FC3_E矩阵的每一个行矢量，即每一个训练矢量（记为FC3_X），从上一次迭代后其归属的胞腔开始搜索，计算当前训练矢量FC3_X与码字FC3_Y_{FC3_I(i)}(FC3_I(i)是FC3_X上一次的编码索引值，对应的FC3_Y_{FC3_I(i)}就是其上一次FC3_X所对应的编码码字）的欧式距离，并将此欧式距离作为当前最小失真FC3_D_min；

步骤4：以码字FC3_Y_{FC3_I(i)}为中心，采用上下查找法，对于当前码字FC3_Y_j，如果满足FC3_D₁=(FC3_X₁-FC3_Y_j1)²≥FC3_D_min，即矢量FC3_X的第一维分量FC3_X₁与码字FC3_Y_j的第一维分量FC3_Y_j1差值的平方大于当前最小欧式距离，就排除码字FC3_Y_j，当FC3_X₁≥FC3_Y_j1时，进入步骤e)，当FC3_X₁≤FC3_Y_j1时进入步骤f)；否则，进入步骤5；

e)排除码字FC3_Y_i,i=1,…,j，转入步骤7；

f)排除码字FC3_Y_i,i=1,…,N₂，转入步骤7；（其中，j为当前搜索

码字的标号，N₂表示第二级矢量量化码书的大小）

步骤5：调用公式FC3_D₂=FC3_D₁+(FC3_V_{FC1_X}-FC3_V_j)²其中， ${\mathrm{FC}3\_V}_{\mathrm{FC}3\_X}=\sqrt{\underset{i=2}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{FC}3\_X}_i^2},$ ${\mathrm{FC}3\_V}_j=\sqrt{\underset{i=2}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{FC}3\_Y}_{\mathrm{ji}}^2},$ M是矢量的维数，计算矢量FC3_X的方差和当前码字FC3_Y_j的方差差值平方与FC1_D₁的和值，如果满足FC3_D₂≥FC3_D_min，则排除码字FC3_Y_j，转入步骤7；否则，转入下一步进行判断；

步骤6：结合部分失真算法，根据公式

计算训练矢量FC3_X与码字FC3_Y_j的部分失真FC3_D_q，其中q为不大于M（维数）的整数。逐渐增加q的数值，一旦存在q满足FC3_D_q≥FC3_D_min，则排除码字FC3_Y_j，转入步骤7；否则，继续累计部分失真，直到q=M，此时如果存在FC3_D(FC3_X,FC3_Y_j)=FC3_D_q<FC3_D_min，则令FC3_D_min=FC3_D_q，更新训练矢量的编码索引值FC3_I(i)为当前最近的码字的标号j，即FC3_I(i)=j；

步骤7：返回步骤4进行下一个码字的匹配，直到搜索完所有的剩余码字，获得当前训练矢量FC3_X的最佳匹配码字和对应编码索引FC3_I(i)。将训练矢量FC3_X划分到第FC3_I(i)个胞腔中，置第FC3_I(i)个胞腔矢量个数计数器FC3_U(FC3_I(i))→FC3_U(FC3_I(i))+1和矢量存储器FC3_V(FC3_I(i))→FC3_V(FC3_V(FC3_I(i))+FC3_X，进入下一步；

步骤8：输入下一个训练矢量，返回步骤2继续训练，直到所有的第二级矢量量化过程的训练矢量被训练完为止，此时求出这N₂个胞腔的质心，即可采用每个胞腔内所有训练矢量的和值∑FC3_V(FC3_I(i))与这个胞腔内训练矢量的数量FC3_U(FC3_I(i))之商∑(FC3_V(FC3_I(i))/FC3_U(FC3_I(i)))，并且将这些质心替换原有码书中的码字，形成新的码书FC3_Y，然后计算更新码字后的平均量化误差，置迭代次数ite₂→ite₂+1；

步骤9：根据设定的迭代次数ite₂，判断是否满足停止迭代的要求，如果满足，所得FC3_Y作为最后的码书，否则转到步骤2进行下一轮训练；

步骤10：以排序索引FC3_I₁重新排序编码索引FC3_I，获得新的编码索引FC3_I₂，这个索引正好对应排序前矩阵FC3_C的每行待编码矢量编码后的对应码字的标号，可经哈达玛反变换后，重建原数据块的编码数据块。码书FC3_Y以及编码索引FC3_I₂是第三部分第二级矢量量化的压缩数据；

在接收端使用每一个行矢量的均值数据，第一级矢量量化过程以及设计的最终码书Y及其索引I₂，第二级矢量量化过程的分维后三部分的最终码书FC1_Y、FC2_Y和FC3_Y以及其相对应的索引值FC1_I₂、FC2_I₂和FC3_I₂，可重构矢量量化处理后的高光谱数据。

高光谱遥感图像压缩的性能主要用图像的峰值信噪比(PSNR)、图像的压缩比(CR)以及生成最终码书所需的计算复杂度来进行评价。

①PSNR计算公式：

$\mathrm{PSNR}=10\&times;\log 10\left(\frac{{\left(\mathrm{Peaksignal}\right)}^2}{\mathrm{MSE}}\right)$

其中，Peaksignal表示超光谱图像矩阵B中的最大像素值，MSE表示第二级矢量量化的均方量化误差；

$\mathrm{MSE}=\frac{\underset{k=1}{\overset{3}{\mathrm{\&Sigma;}}}{\mathrm{MSE}}_k\&times;N_k}{\underset{k=1}{\overset{3}{\mathrm{\&Sigma;}}}{N}_k}$

其中：MSE_k分别为分维后第二级矢量量化第k部分图像平均每一个像素点的误差；N_k表示分维后第k部分第二级矢量量化过程中图像的列数（维数）。

其中根据如下公式计算MSE_k：

${\mathrm{MSE}}_k=\underset{i=1}{\overset{65536}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N_k}{\mathrm{\&Sigma;}}}{(\mathrm{fc}\left(k\right)\_y(i,j)-\mathrm{fc}\left(k\right)\_b(i,j))}^2/(65536\&times;N_k),k=\mathrm{1,2,3}$

其中，fc(k)_y(i,j)表示分维后第k部分第二级矢量量化中编码后第i个矢量的第j维分量，fc(k)_b(i,j)表示分维后第k部分第二级矢量量化编码前的第i个矢量的第j维分量，65536表示差值图像中矢量个数（行数），N_k表示分维后第k部分第二级矢量量化过程中图像的列数（维数）。

②CR计算公式：

$\mathrm{CR}=\frac{65536\&times;16\&times;224}{(N_1+N_2)\&times;16\&times;224+65536\&times;{\log }_2\left(N_1\right)+3\&times;65536\&times;{\log }_2\left(N_2\right)+65536\&times;16}$

其中，N₁为第一级矢量量化的码书尺寸,N₂为第二级矢量量化的码书尺寸，224为矢量维数，16代表每个像元用2个字节表示，65536表示矢量个数，65536×16表示行矢量均值。

③计算复杂度计算公式：

$\mathrm{Complex}=\frac{{\mathrm{Num}1}_i\&times;{\mathrm{ite}}_1+{\mathrm{Num}2}_i\&times;{\mathrm{ite}}_2}{{\mathrm{ite}}_1+{\mathrm{ite}}_2},i=\mathrm{1,2,3,4}$

其中：Num1_i依次表示第一级矢量量化过程中平均每个矢量（长度为224）每次迭代的比较（com），乘法（×），加法（±）和开方（sqrt）的次数；Num2_i依次表示第二级矢量量化过程中平均每个矢量（长度为224）每次迭代的比较（com），乘法（×），加法（±）和开方（sqrt）的数；i=1,2,3,4依次表示比较（com），乘法（×），加法（±）和开方（sqrt）；ite₁,ite₂分别表示第一级和第二级矢量量化的迭代次数。

Num2_i的计算公式为：

${\mathrm{Num}2}_i=\frac{\underset{j=1}{\overset{3}{\mathrm{\&Sigma;}}}{\mathrm{Num}2\left(j\right)}_i\&times;N_j}{224},i=\mathrm{1,2,3,4}$

其中：Num2(j)_i表示分维后第j部分第二级矢量量化过程中平均每个矢量每次迭代的比较（com），乘法（×），加法（±）和开方（sqrt）的次数；N_j表示分维后第j部分矢量的维数（列数）。

以Lunar Lake和Low Altitude的第一场景的部分数据作为仿真测试信源，通过以上3个评价标准对LBG算法、文献（申请号：201110272304.X的专利申请）与本发明方案进行比较，表1列出了三种方法在码书相同条件下经过10次迭代的压缩比CR、峰值信噪比PSNR，其中LBG算法和文献（申请号：201110272304.X）的码书大小为N，迭代次数ite=10，而本发明采用多级量化方式，第一级矢量量化过程码书为N₁，迭代次数为ite₁=3，第二级矢量量化过程码书为N₂，迭代次数为ite₂=7，其中N=N₁×N₂；表2列出了Lunar Lake图像在码书尺寸相同的情况下经过三种算法处理后的运算复杂度，表3列出了LowAltitude图像在码书尺寸相同的情况下经过三种算法处理后的运算复杂度，其中表2和表3的参数设置与表1相同。

表1.在相同码书尺寸下经10次迭代三种算法的峰值信噪比

表2.Lunar Lake图像在相同码书尺寸下三种算法处理后的计算复杂度

表3.Low Altitude图像在相同码书尺寸下三种算法处理后的计算复杂度

从实验结果可以看出：本发明提出的方法在码书尺寸相同的情况下，与LBG算法和文献（专利申请号:201110272304.X）相比，在压缩比、图像恢复质量（峰值信噪比）和计算复杂度上均有所改善，其中在码书尺寸为4096时，本发明相对于LBG算法，Lunar Lake的压缩比提高了一倍，图像峰值信噪比提高了1.23dB（表1），计算复杂度除了开方次数有些许提高，比较次数降低了81.2%，乘法次数和加法次数都不到LBG算法的1%（表2）。Low Altitude图像的压缩比提高了一倍，峰值信噪比提高了1.70dB（表1），计算复杂度除了开方次数有些许提高，比较次数降低了81.7%，乘法次数和加法次数都不到LBG算法的1%（表3）。在码书尺寸为4096时，本发明Lunar Lake的压缩比提高了一倍，图像峰值信噪比提高了0.80dB（表1），计算复杂度除了开方次数有些许提高，比较次数降低了81.7%，乘法次数降低了78.3%，加法次数降低了72.3%（表2）；LowAltitude图像的压缩比提高了一倍，峰值信噪比提高了0.42dB（表1），计算复杂度除了开方次数有些许提高，比较次数降低了80.2%，乘法次数和加法次数都不到只有LBG算法的1%。（表3）。同时在码书尺寸为4096时，对于Lunar Lake图像的部分数据的处理，量化编码的MSE=360.8，对于Low Altitude图像的部分数据的处理，量化编码的MSE=909.8，即每一个像元的平均量化误差在31以下，对于以两个字节表示一个像元数据的高光谱数据而言，如此小的量化误差，使得这时的矢量量化压缩近似于无损压缩，本发明的处理数据具有很高的工程应用价值。随着高光谱技术的发展和广泛应用，高光谱遥感图像的波段数不断增加，成像光谱仪获取的数据量越来越大，这给数据的存储和传输带来了巨大的压力，对高光谱数据进行有效压缩显得非常迫切。这里提出的高光谱遥感图像压缩编码方法，是面向工程应用的重要探索。

Claims (8)

1.一种多级和分维矢量量化的高光谱遥感信号快速编码方法，其特征在于，读取高光谱数据源，提取高光谱数据各行矢量的均值并保存，对高光谱数据的行矢量进行去均值操作获得去均值图像，构造第一级矢量量化的初始码书，获取第一级编码码书和编码；用第一级编码码书及编码索引重构的图像，与去均值图像相减得到差值图像，对差值图像矢量进行分维处理，对分维处理后各部分分别构造对应的第二级矢量量化的初始码书作为第二级矢量量化的输入数据，训练该初始码书得到各部分的第二级矢量量化的码书和编码索引；由行矢量均值、第一级码书和编码索引、以及分维的各部分第二级矢量量化码书和编码索引构成高光谱数据源的压缩数据。

2.根据权利要求1所述的方法，其特征在于，所述去均值操作具体为：对高光谱的行矢量减去该行矢量均值，生成差值数据矩阵。

3.根据权利要求1所述的方法，其特征在于，获得第一级矢量量化的初始码书具体包括：对高光谱数据源矩阵进行补零操作，进行哈达玛变换得到的矢量数据按照原行矢量均值的大小进行升序排列，获得排序后索引值构成矩阵，然后对排序后的矩阵平均分组，依次选取每组第一个矢量作为第一级矢量量化的初始码书。

4.根据权利要求1所述的方法，其特征在于，所述差值图像具体为：去均值后的原图像减掉由第一级矢量量化得到的重构图像，生成差值图像。

5.根据权利要求1所述的方法，其特征在于，所述对差值图像矢量进行分维处理具体为：按照2的整数次幂对差值图像矢量进行分维，将矢量数据分为三个部分，每一个部分都单独作为第二级矢量量化的一个输入数据。

6.根据权利要求1所述的方法，其特征在于，对于分维后的第二级矢量量化的初始码书构造：分别对分维后的图像矩阵进行哈达玛变换，将变换后矢量数据按照第一维分量值进行升序排列，获得排序索引，对排序索引平均分组，依次选取每组首矢量分别得到每一个部分第二级矢量量化的初始码书。

7.根据权利要求1所述的方法，其特征在于，构造第一级矢量量化的初始码书具体包括：⑴读取高光谱遥感图像矩阵，将其转化为2维矩阵B；⑵将矩阵B的每一行矢量都减去该行矢量所对应的均值，生成去均值矩阵C；⑶使C的行向量的维数转化成2的整数次幂的维数得到矩阵C’；⑷对矩阵C’进行哈达玛变换，将变换后的矩阵按照原行矢量均值对矢量进行升序排列，获得矢量排序后的矩阵E并记录此时排序索引I₁；⑸对矩阵E按行根据第一级矢量量化的初始码书尺寸N₁平均分组，依次选择每组首矢量生成初始码书。

8.根据权利要求1所述的方法，其特征在于，获得第二级矢量量化的初始码书具体包括：利用去均值图像和第一级矢量量化的重构图像，生成差值图像CZ；对差值图像矩阵的行矢量按照2的整数次幂进行分维，生成分维矩阵FC1、FC2和FC3；分维后各部分行矢量进哈达玛变换，变换后的矢量依其第一维分量值进行升序排列，获得相应排序后的矩阵作为第二级矢量量化训练矢量集；对此矩阵按行根据第二级矢量量化的初始码书尺寸N₂平均分组，依次选择每组首矢量生成对应的第二级各个分维部分的初始码书。

Images