CN102905137B - 超光谱信号的快速差值矢量量化压缩编码方法 - Google Patents

Abstract

本发明提出了一种基于差值矢量量化的超光谱信号快速压缩方案。在哈达玛域中，构造初次矢量量化的初始码书，快速聚类，生成初次矢量量化的最终码书和索引；将初次矢量量化恢复的图像与原图像获得差值图像数据作为差值矢量量化的输入矢量，对其进行哈达玛变换，构造差值矢量量化的初始码书，快速聚类，生成差值矢量量化的最终码书和索引。本发明提供的超光谱图像压缩方案，结合Hadamard变换性质和差值矢量量化减小量化误差的特点，与传统的LBG算法相比，在压缩比相同的情况下，不仅能够大幅减少计算量，而且可以提高图像恢复质量，达到快速压缩编码的目的。

Description

超光谱信号的快速差值矢量量化压缩编码方法

技术领域

本发明属于超光谱遥感图像处理领域，具体涉及一种基于差值矢量量化技术的超光谱图像数据压缩方法。

背景技术

随着遥感技术在地球资源管理、环境检测、军事侦测等领域的广泛应用，以往多光谱图像数据由于仅仅在图像的空间分辨率上有所改善，其数据特性已经满足不了人们在生产和科研上的需要，如SPOT和Landset TM图像，其图像仅包含4～7个离散的谱带，而人们感兴趣的是吸收特性宽度在20nm到40nm的地物目标，因此当谱带内的地物目标光谱特性相近时，多光谱图像的应用就受到限制，此时必须利用更多的谱带（通常需要包括几百个谱带），成像光谱仪的问世使获取包含几百个谱带的数据图像成为可能，这种图像就是超光谱图像。

与多光谱图像相比，由于成像光谱仪在获取地物空间分布信息的同时，对每一个像元形成一条地物光谱曲线，使得超光谱图像可以区分并识别更多种地物目标，但这是以较大的数据量和较高的数据维为代价来换取的。以典型的AVIRIS图像为例，其图像大小为614×512×224，如果对其以每个像素灰度值为16bit存储，那么总的数据存储量大约为140M bit，显然要比多光谱图像的数据量大得多,这给超光谱图像的传输和存储带来极大的困难。因此，寻求一种有效的压缩技术具有十分重要的意义。

超光谱图像具有两种相关性：空间相关性和谱间相关性。空间相关性是指每个谱段内某一像元与其相邻像元间的相似性。谱间相关性是指相邻波段间对应位置的像元具有的相似性，谱间相关性又分为谱间统计相关性和谱间结构相关性。使用多维信号处理技术，才能有效利用这些相关性。

将基于多维信号处理矢量量化技术，作为超光谱图像压缩的一种有效手段，具有压缩比大，编解码简单，失真较小的优点。矢量量化的核心内容是码书设计和码字快速搜索。码书设计就是寻找最优码书，使恢复图像与原始图像之间的计算失真达到最小，以保证重构图像的质量；码字快速搜索就是如何快速地找到与输入矢量失真最小的码字。

LBG算法作为矢量量化中码书设计的经典算法，基本原理是从某个初始码书开始，根据最近邻条件得到新的胞腔，然后根据质心条件由新的胞腔得到新的码字，形成新的码书；对于一个训练序列{X_j;j=1,…,L},其具体过程如下：首先初始化，设初始码书为Y⁽⁰⁾={Y_i;i=1,…,N},其中Y_i表示码字，码书大小为N，设置失真阈值为ε或最大迭代次数为ite，令初始迭代次数t=0,平均失真D_-1=∞，然后用码书中的各个码字作为聚类中心，根据最佳划分准则（即找到每个矢量与其失真最小的码字），把训练矢量集划分为N个胞腔C⁽ⁿ⁾={S_i;i=1,2,…,N}，其中,对任意l∈{1,2,….,N}成立；最后计算平均失真。如果相对误差满足(D_t-1-D_t)/D_t≤ε或迭代次数达到设定的最大值就停止运算，Y(t)即为最终码书；否则，计算每个最小失真划分的质心=，由这些新质心组成新的码书，置t=t+1,进行下一次迭代。然而从运算复杂度来看，采用LBG进行数据压缩存在码书自适应能力不强，运算量大的缺点，并且获得高压缩比时图像的恢复质量不高。

发明内容

本发明针对现有超光谱图像压缩技术，在获得高压缩比时图像的恢复质量不高，计算量大的问题，提出了一种基于差值矢量量化的超光谱图像快速压缩编码方法，在保证压缩比不变的情况下，提高超光谱图像的恢复质量，

本发明解决上述技术问题的技术方案是：提出一种差值矢量量化的超光谱信号快速编码方法，包括步骤，读取超光谱数据源，得到初次矢量量化的初始码书，根据图像压缩质量设置最大迭代次数或失真阈值；计算训练矢量X和当前码字Y_b的欧氏距离，根据欧氏距离确定训练矢量的最佳匹配码字，依次将训练矢量划分到它所对应最佳匹配码字的胞腔中，直到划分完所有的训练矢量，最后用胞腔的质心来更新当前最佳匹配码字；用初次矢量量化的生成码书及训练矢量对应的最佳匹配码字的索引构造哈达玛域的恢复图像，经过哈达玛反变换及降维生成空域恢复图像，将原始图像与恢复图像作差生成差值图像，获得差值矢量量化的初始码书；搜索差值训练矢量的最佳匹配码字，依次将训练矢量划分到它所对应的胞腔中，直到划分完所有的训练矢量，最后用胞腔的质心来更新当前最佳匹配码字获得差值矢量量化的最终码书；当累计迭代次数达到设定的最大迭代次数或最近两次平均失真的相对误差值小于设置的失真阈值，打包差值矢量量化的最终码书以及所有矢量最佳匹配码字的索引值。

该方法主要包括以下几个阶段：

初次矢量量化初始化阶段：读取超光谱数据源，对超光谱数据源矩阵进行哈达玛变换得到的矢量数据按照第一维分量值的大小进行升序排列，获得排序后的矩阵E以及排序索引值，然后对排序后的矢量平均分组，依次选取每组首矢量得到初次矢量量化的初始码书，并根据压缩质量设置初次矢量量化过程的最大迭代次数或失真阈值。

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

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

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

生成初始码书具体为，将超光谱数据源矩阵进行扩展，形成列数（维数）为2的整数次幂的行矢量，对其进行Hadamard变换并进行升序排序，根据排序索引对行排序后的矩阵E按照行平均分N₁组（N₁为初次矢量量化码书的尺寸），依次选择每组第一个矢量作为该组的初始码字生成初始码书Y；将生成的差值矩阵进行扩展，形成列数（维数）为2的整数次幂的行矢量，对其进行Hadamard变换并进行升序排序，根据排序索引对行排序后的矩阵CZ_E按照行平均分N₂组（N₂为差值矢量量化码书的尺寸），依次选择每组第一个矢量作为该组的初始码字生成初始码书CZ_Y。

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

附图说明：

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

图2 本发明压缩编码方法流程框图；

图3 本发明实施方式具体流程图。

具体实施方式

以下使用具体事例和附图对本发明作进一步说明。本发明提出的基于差值矢量量化的超光谱信号快速编码方法，具体步骤如图1、图2所示。主要包括：

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

快速搜索训练矢量的最佳匹配码字，依次将训练矢量划分到它们所对应最佳匹配码字的胞腔中，直到划分完所有的训练矢量，最后用胞腔的质心来更新当前码字。

在每轮迭代之前预先判断累计迭代次数是否达到设定的最大迭代次数，若满足，就结束训练并打包初次矢量量化的最终码书和所有矢量最佳匹配码字的索引；否则，继续训练。

用初次矢量量化的最终码书、每个输入矢量对应匹配码字的索引、哈达玛反变换和变维（去除零矩阵）重构恢复图像，然后用原图像与恢复图像生成差值图像；对差值图像矢量进行哈达玛变换，对变换后矢量数据按照第一维分量值进行升序排列，获得排序索引，然后利用平均分组选取每组首矢量得到初始码字，形成差值矢量量化的初始码书，并设置差值矢量量化过程的迭代次数。

搜索训练矢量的最佳匹配码字，依次将训练矢量划分到它所对应最佳匹配码字的胞腔中，直到划分完所有的训练矢量，最后用胞腔的质心来更新当前码字。判断累计迭代次数是否达到设定的最大迭代次数，若满足，结束训练并打包差值矢量量化的最终码书以及所有矢量最佳匹配码字的索引；否则，继续训练。

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

计算训练矢量X和当前码字Y_b的欧氏距离，将此欧氏距离作为当前最小失真D_min，对于当前码字Y_j1，如果满足D₁=(X₁-Y_j1)²≥D_min，排除码字Y_j1；如果不满足D₁≥D_min，计算哈达玛域矢量方差D₂=(X₁-Y_j1)²+(V_X-V_j)²，( ，M是矢量的维数），如果满足D₂≥D_min，则排除码字Y_j1 ；如果还不满足D₂≥D_min，计算部分失真，即判断1≤q≤M，q为计算的矢量的维数，M为矢量的维数，将q的取值从1逐渐增加到M，如果一旦满足D_q(X,Y_j)≥D_min，则排除码字，计算出排除后剩余的码字与训练矢量X的欧式距离，找出其距离最近的码字，记录相应码字的索引号I(j)，将训练矢量X划分到第I(j)个胞腔中。

其中， X₁、Y_j1分别表示训练矢量和码字的第一维分量，V_X、V_j分别表示X和Y_j的子矢量的方差，D_min为当前与匹配码字的子矢量的最小欧式距离。D(X, Y_j)表示训练矢量X与码字Y_j之间的子矢量的欧式距离。

以下以典型的AVIRIS超光谱图像（大小为512×614×224）为例具体描述本发明的实现方式：它含有224个波段，每个波段的图像尺寸是512×614。为了处理方便，在每个波段中截取大小为256×256的图像块，这样生成的图像块矩阵P为(256,256,224)。将此3维矩阵P转化为2维矩阵B，其转换过程是将矩阵P的前两维分量作为B矩阵的行向量，第3维数据作为B矩阵的列向量，即B为(256×256,224)。

在MATLAB7.9软件平台下结合附图对本发明方案的实例进行详细说明。使用512×614×224规格的超光谱图像数据，其每个波段中像元数据以两个字节的带符号整数表示，本方案的具体实施步骤如下：

初次矢量量化初始化阶段：

步骤1：利用MATLAB函数库中的fopen函数和fread函数读取超光谱数据源，得到超光谱数据源的3维矩阵A(包含512行，614列，224个波段，简记为

A(512,614,224）)；

步骤2：选取矩阵A的前256行、256列和所有波段，组成新的3维矩阵P(256,256,224)，将3维矩阵P转化为2维矩阵B，以便进行矢量量化处理。其变换方式是将3维矩阵P的前两维转换为2维矩阵B的第一维，第3维作为矩阵B的第二维，最终生成矩阵B(256×256)行和224列的矩阵，简记为（65536，224）)；

步骤3：在矩阵B的每一行矢量后添加32维的零矢量，形成列数（维数）为2的整数次幂（即256维）的行矢量，获得矩阵B'(65536,256)，为后面的Hadamard变换做准备；

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

步骤5：根据初次矢量量化过码书尺寸N₁，对于训练矢量集E（待编码矢量集），设定编码索引矩阵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：以码字Y_I(i)为中心，使用上下查找法，对于当前码字Y_j，如果满足D₁=(X₁-Y_j1)²≥D_min，即训练矢量X的第一维分量与码字Y_j的第一维分量差值的平方大于当前最小欧氏距离，就排除码字Y_j，当X₁≥Y_j1时，进入步骤a)，当X₁≤Y_j1时，进入步骤b)；否则，进入步骤5；

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

b)排除码字Y_i，i=j,…,N₁，转入步骤7；（j为当前搜索码字的下标，N₁表示码书尺寸。）

步骤 5：计算D₂=D₁+(V_X-V_j)²（其中，，M是矢量的维数），即计算训练矢量X与码字Y_j方差差值平方与D₁的和值，如果满足D₂≥D_min，则排除码字Y_j，转入步骤7；否则，转入下一步进行判断；

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

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

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

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

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

差值矢量量化初始化阶段：

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

步骤 2：将空域恢复矩阵HF与原图像的空域矩阵Ｂ(65536,224)做差值，生成差值矩阵CZ(65536,224)；

步骤 3：在2维矩阵CZ的每一行矢量后添加32维的零矢量，形成维数为2的整数次幂（即256维）的行矢量，获得矩阵CZ＇(65536,256)，为下一步的行变换做准备；

步骤 4：对CZ＇进行哈达玛变换得到变换矩阵CZ_C(65536,256)，将变换矩阵CZ_C按照第一列分量的大小重新进行升序排列，记录这时的排序索引CZ_I₁，获得行排序后的矩阵CZ_E，然后对矩阵CZ_E按行平均分为N₂组（N₂为差值矢量量化的码书尺寸），依次选择每组第一个矢量作为该组的初始码字生成初始码书CZ_Y；

步骤 5：依据差值矢量量化编码过程码书大小N₂，对于排序后的训练矢量（即待编码矢量），设定编码索引矩阵CZ_I(1,65536)，并置训练矢量的初始编码索引为，为每个训练矢量分配索引值。将训练矢量按照初始编码索引的索引值归属到相应的胞腔中（具有相同索引值的矢量属于同一个胞腔）；

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

差值矢量量化训练阶段：

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

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

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

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

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

d)排除码字CZ_Y_i,i=1,…,N₂，转入步骤7；（其中，j为当前搜索码字的标号，N₂表示差值矢量量化码书的大小）

步骤 5：调用公式CZ_D₂=CZ_D₁+(CZ_V_{CZ_X}-CZ_V_j)²其中， ，M是矢量的维数，计算矢量CZ_X的方差和当前码字CZ_Y_j的方差差值平方与CZ_D₁的和值，如果满足CZ_D₂≥CZ_D_min，则排除码字Y_j，转入步骤7；否则，转入下一步进行判断；

步骤 6：结合部分失真算法，根据公式计算训练矢量CZ_X与码字CZ_Y_j的部分失真CZ_D_q，其中q为不大于M（维数）的整数。逐渐增加q的数值，一旦存在q满足CZ_D_q≥CZ_D_min，则排除码字CZ_Y_j，转入步骤7；否则，继续累计部分失真，直到q=M，此时如果存在CZ_D(CZ_X,CZ_Y_j)=CZ_D_q＜CZ_D_min，则令CZ_D_min=CZ_D_q，更新训练矢量的编码索引值CZ_I(i)为当前最近的码字的标号j，即CZ_I(i)=j；

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

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

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

步骤 10：以排序索引CZ_I_i重新排序编码索引CZ_I，获得新的编码索引CZ_I₂，这个索引正好对应排序前矩阵CZ_C的每行待编码矢量编码后的对应码字的标号，可经哈达玛反变换后，重建原数据块的编码数据块，将码书CZ_Y以及编码索引CZ_I₂作为差值矢量量化过程的最终码书和索引，进行存储或传输；

在接收端按照初次矢量量化过程设计的最终码书Y以及其索引I₂和差值矢量量化过程的最终码书CZ_Y和其索引值CZ_I₂，即可恢复经过矢量量化处理后的信源信息。

文中所使用的参数和变量为： x表示空域中的矢量，X表示哈达玛域中的矢量，y表示空域中的码字，Y表示哈达玛域中的码字，ite表示迭代次数，D_min表示当前最小失真，N表示码书尺寸，I表示索引矩阵，E表示哈达玛域矢量排序后的训练矩阵矢量，V表示胞腔矢量存储矩阵，CZ表示差值矢量量化过程的前缀，PSNR表示峰值信噪比，CR表示压缩比。

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

以Lunar Lake和Low Altitude的第一场景的部分数据作为仿真测试信源，通过以上3个评价标准对LBG算法与本发明方案进行比较，表1列出了两种算法在压缩比相同条件下经过10次迭代的峰值信噪比PSNR，其中LBG算法的码书大小为N，迭代次数ite=10，而本发明采用差值量化方式，初次矢量量化过程码书为N₁，迭代次数为ite₁=5，差值矢量量化过程码书为N₂，迭代次数为ite₂=5；表2列出了Lunar Lake图像经过两种算法处理后的运算复杂度，表3列出了Low Altitude图像经过两种算法处理后的运算复杂度，其中表2和表3的参数设置与表1相同。

表1. 在相同压缩比时经10次迭代两种算法的峰值信噪比

表2. Lunar Lake图像在相同压缩比下经过两种算法处理后的计算复杂度

表3. Low Altitude图像在相同压缩比下经过两种算法处理后的计算复杂度

从实验结果可以看出：本发明提出的方法相对于LBG算法在保证压缩比不变、计算复杂度大幅降低的情况下，提高了图像恢复质量（峰值信噪比），其中在压缩比CR=29.1时，Low Altitude图像峰值信噪比提高了2.40dB，Lunar Lake图像峰值信噪比提高了2.46dB（表1）；同时在CR=29.1对于Low Altitude图像的部分数据的处理，量化编码的MSE=793.5，对于Lunar Lake图像的部分数据的处理，量化编码的MSE=379.6，即每一个像元的平均量化误差在30以下，此时，对于以两个字节表示一个像元数据的超光谱数据而言，如此小的量化误差，使得这时的有损矢量量化压缩近似于无损压缩，本发明的处理数据具有很高的工程应用价值。表1中给出了各种压缩比的情况，相应的处理数据可以满足从一般的浏览领域到实际的工程应用领域的各种需求。

Claims (6)

1.一种差值矢量量化的超光谱信号快速编码方法，其特征在于，读取超光谱数据源，得到初次矢量量化的初始码书，根据图像压缩质量设置最大迭代次数或失真阈值；计算训练矢量X和当前码字Y_b的欧氏距离，根据欧氏距离确定训练矢量的最佳匹配码字，依次将训练矢量划分到它所对应最佳匹配码字的胞腔中，直到划分完所有的训练矢量，最后用胞腔的质心来更新当前最佳匹配码字；用初次矢量量化的生成码书及训练矢量对应的最佳匹配码字的索引重构图像，并经过哈达玛反变换及降维生成空域恢复图像，将原始图像与空域恢复图像作差生成差值图像，获得差值矢量量化的初始码书；搜索差值训练矢量的最佳匹配码字，依次将训练矢量划分到它所对应的胞腔中，直到划分完所有的训练矢量，最后用胞腔的质心来更新当前最佳匹配码字获得差值矢量量化的最终码书；当累计迭代次数达到设定的最大迭代次数或最近两次平均失真的相对误差值小于设置的失真阈值，打包差值矢量量化的最终码书以及所有矢量最佳匹配码字的索引值。

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

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

4.根据权利要求1所述的方法，其特征在于，对于单次差值矢量量化，生成初始码书具体为：读取超光谱图像构建3维矩阵，将其转化为2维矩阵B；使B的行向量的维数转化成2的整数次幂的维数得到矩阵B’；对矩阵B’进行哈达玛变换，将变换后的矩阵按照第一维分量值对矢量进行升序排列，获得矢量排序后的矩阵E并记录此时排序索引I₁；对矩阵E按行根据码书尺寸N₁平均分组，依次选择每组首矢量生成初次矢量量化的初始码书。

5.根据权利要求1所述的方法，其特征在于，对于单次差值矢量量化，获得差值矢量量化的初始码书具体为：利用原始图像B和恢复图像生成2维差值图像CZ；判断差值矩阵CZ的维数是否为2的整数次幂，如果不是，将差值矩阵行向量转化为2的整数次幂的矩阵CZ’；对CZ’进行哈达玛变换，将变换后的差值矢量的第一维分量值进行升序排列，获得矢量排序后的矩阵CZ_E为差值矢量量化训练矢量集，并记录此时排序索引CZ_I₁；对矩阵CZ_E按行根据差值矢量量化的码书尺寸N₂平均分组，依次选择每组首矢量生成差值矢量量化的初始码书。

6.根据权利要求1所述的方法，其特征在于，将欧氏距离作为当前最小失真D_min，将码字Y_j与训练矢量X₁比较，如果满足D₁＝(X₁-Y_j1)²≥D_min，排除码字Y_j；否则，根据公式D₂＝(X₁-Y_j1)²+(V_X-V_j)²计算，如果满足D₂≥D_min，则排除码字Y_j，其中否则，根据公式 1≤q≤M计算部分失真，将q的取值从1逐渐增加到M，一旦D_q(X,Y_j)≥D_min，就排除码字Y_j，其中，M为矢量的维数。