一种超谱信号的快速压缩编码方法及图像压缩方法

摘要

本发明针对超谱图像巨大的数据量，提出了一种快速的超谱信号压缩方案。压缩编码方案采用矢量量化聚类技术消除超谱图像的空间相关性和谱间相关性，结合传统的LBG聚类算法以及快速码字搜索算法生成优越的码书，并同时完成编码。对输入矢量进行哈达玛变换，然后在哈达玛域中，构造初始码书，同时利用3步排除不等式，实施快速聚类，高效生成最终码书和编码。本发明提供的超谱图像压缩方案，结合哈达玛变换性质和矢量量化的特点，能够以较小的图像质量损失来换取较大的压缩比，达到快速编码的目的。

说明书

技术领域

本发明属于图像压缩编码技术领域，具体涉及一种基于矢量量化数据压缩技术的超谱图像压缩。

背景技术

   随着遥感技术在地球资源管理、环境检测、军事侦测等应用领域的不断扩大，以往的多光谱图像数据由于仅仅在空间分辨率上有所提高，己经不能满足人们日益增长的科研和生产需求，而成像光谱仪的问世使能获取超谱分辨率的图像数据成为可能。超谱图像具有较高的波段数（其波段数一般都为几百个）和大数据量的特点，典型AVIRIS图像的大小为614×512，成像波段数为224，如果每像素灰度值以两个字节存储，一组图像的数据量则约为140M字节，如此海量的超谱数据给有限带宽的通信系统带来了传输、存储和管理带来了巨大的困难，寻求一种有效的压缩技术具有重要的意义。

   由成像光谱仪在获取地物空间分布信息的同时，对每一个像元形成一条地物光谱曲线。根据像元的光谱曲线可以确定像元处物体的特性。所谓像元是指采样单位，它是与分辨率相关的，假设分辨率为1千米，像元就代表地面1千米×1千米的面积，仅超谱图像上的一个像元就描述了一大块地面信息。对超谱图像压缩的同时，又必须考虑图像的质量，即压缩后的超谱图像必须能够得到清晰的地物光谱曲线以供分析人员使用。

   超谱图像具有两种相关性：空间相关性和谱间相关性。空间相关性是指每个谱段内某一像元与其相邻像元间的相似性。谱间相关性是指相邻波段间对应位置的像元具有较强的相关性，谱间相关性又分为谱间统计相关性和谱间结构相关性。

   矢量量化技术是一种有效的有损压缩技术，具有压缩比大、编解码简单、失真较小的优点，不仅广泛应用到图像和语音压缩编码等领域，也被广泛应用于其它众多领域，包括移动通信、图像检索、语音识别、遥感分类等方面。它的基本原理就是将一个包含多个矢量的对象分成若干个区域，然后将这个区域划分为若干个区域，通过寻找各个区域中最能匹配此区域的一个矢量来代替这个区域，编码完成后，发送这些矢量的索引，接收端通过查找码表来获得此区域的匹配点，最终以较低的失真恢复这个对象。在发送端只需发送匹配矢量的索引，相比数据的发送，可以有效的提高压缩比。将矢量量化这种聚类技术应用于超谱图像领域中，可以有效的消除超谱图像空间和谱间相关性，达到高压缩比的要求。

   码书设计和码字搜索是矢量量化的核心，优越的码书能够使图像获得较好的图像质量，快速的码字搜索则能够提高搜索最佳匹配码字的速度，从而提高搜索效率，降低运算复杂度。

   由Linde、Buzo和Gray提出的LBG算法（也叫做GLA算法），是矢量量化发展的基石，它是在每一轮迭代中将训练矢量放入各个胞腔，然后取胞腔的质心作为码书中的码字。

LBG算法是一种基于初始码书的迭代算法，一个不合理的初始码书容易使最终码书陷入局部最优，而且迭代多次后，效果微乎其微。而优越的码书可以使LBG算法的最终码书接近全局最优。在LBG迭代过程中始终满足最邻近条件和质心条件。其具体步骤如下：

初始化。设初始码书                                                , N为码字个数，失真阈值,一个训练序列,平均失真；用码书Y⁽ⁿ⁾中的各个码字作为聚类中心，根据最佳划分准则把训练矢量集划分为N个胞腔，其中，对任意成立。计算平均失真；如果相对误差满足，停止运算，Y⁽ⁿ⁾即为最终码书。否则转入下一步。计算每个最小失真划分的质心=，由这些新质心组成新的码字。然而从运算复杂度来看，采用LBG进行数据压缩主要存在码书自适应能力不强，运算量大的缺点。

发明内容

本发明针对现有技术数据压缩编码自适应能力不强，计算量大的缺陷，提出一种超谱信号的快速编码方法，在保证图像质量的前提下，得到较高的压缩比和高效率编码。

本发明解决上述技术问题的技术方案是，提出一种超谱信号的快速编码方法，主要包括以下几个阶段：

初始化阶段:按照输入矢量第一维大小排序，进行平均分组→选择每组第一个训练矢量生成初始码书，根据各个训练矢量和初始码书中码字的误差计算平均误差，并设置迭代次数。

训练阶段:利用不等式排除算法快速搜索当前最佳匹配码字，将当前训练矢量划分到当前最佳匹配码字所对应的胞腔，训练完所有的训练矢量，以胞腔的质心更新当前最佳匹配码字。

中断检查: 若迭代次数满足要求，迭代结束，打包码书和所有矢量匹配码字的索引。

具体包括步骤：读取超谱图像数据，选取每个波段中的一段图像块作为矢量量化的输入信源，根据图像尺寸和波段数构建3维矩阵，并转化为2维矩阵；将2维矩阵每一行的后面添加全0矢量，形成具有2的整数幂次的行矢量的矩阵B'；对矩阵B'的行进行哈达玛(Hadamard)变换，将变换后的矩阵按照第一列分量重新进行行排序，记录这时的排序索引I₁，获得行排序后的矩阵E，然后对矩阵E按行平均分组，依次选择每组第一个训练矢量生成初始码书；设定训练矢量的胞腔索引矩阵I，将初始码书中训练矢量归属到相应的胞腔中；利用不等式排除算法快速搜索当前最佳匹配码字，将当前训练矢量划分到最佳匹配码字所对应的胞腔，训练完所有的训练矢量，以胞腔的质心更新当前最佳匹配码字。

所述不等式排除算法具体为：根据编码索引I，选定当前匹配码字Yp，计算训练矢量X和匹配码字Yp的欧式距离，将欧式距离作为当前最小失真D_min，进一步搜索码字Yj，如果满足D₁=(X₁－ Yj1)²≥D_min，排除码字Yj；否则根据哈达玛方差公式计算D₂=(X₁－ Yj1)²+(V_X-Vj)² ，如果满足D₂≥D_min，则排除码字Yj；计算剩下的码字与训练矢量X的欧式距离，找出距离最近的索引I(i)，将训练矢量X划分到第I(i)个胞腔中。

在超谱信号快速编码基础上，本发明提出一种超谱图像压缩方法，包括步骤，读取超谱图像数据，选取每个波段中的一段图像块作为矢量量化的输入信源，根据图像尺寸和波段数构建3维矩阵，将3维矩阵的各列首尾相接转化为2维矩阵；将2维矩阵转换为具有2的整数幂次的行矢量的矩阵B'；对矩阵B'的行进行哈达玛变换，并按照第一列分量重新进行行排序，记录这时的排序索引I₁，然后对行排序后的矩阵E按行平均分组，依次选择每组第一个训练矢量生成初始码书；设定训练矢量的胞腔索引矩阵I，将初始码书中训练矢量归属到相应的胞腔中；利用排除算法快速搜索当前最佳匹配码字，将当前训练矢量划分到当前最佳匹配码字所对应的胞腔，训练完所有的训练矢量，以胞腔的质心更新作为当前最佳匹配码字；形成新的码书Y；以排序索引I₁重新排序编码索引I，获得新的编码索引I₂，该索引对应排序前矩阵C的每行待编码矢量编码后的对应码字的标号，经哈达玛反变换后，将码书Y以及编码索引I₂作为超谱数据块的压缩数据。

   在实际的超谱图像传输的时候，只需要传输码书以及各个训练矢量对应于码书中的索引，在接收端只需要利用索引查找码书中对应的码字来代替训练矢量，从而可以恢复超谱图像信号。值得注意的是，传输数据的大小相对于原始超谱图像的数据量来说是非常小的，因此，通过矢量量化编码可以达到压缩超谱图像的目的。采用本发明提出的超谱信号快速编码方法进行图像压缩，码书自适应能力强，运算复杂度低，在保证图像质量的前提下，得到较高的压缩比和高效率编码。

附图说明

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

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

图3是美国喷气推进实验室提供的Lunar Lake超谱图像第一个场景的第十波段部分数据以及根据本发明方法处理后的显示，从上到下分别是原图，以码书尺寸为2048时矢量量化后的图，以及这两个图的数据的误差图。

图4是美国喷气推进实验室提供的Low Altitude超谱图像第一个场景的第十波段部分数据以及根据本发明方法处理后的显示，从上到下分别是原图，以码书尺寸为2048时矢量量化后的图，以及这两个图的数据的误差图。

具体实施方式

以下针对附图和具体实例对本发明的实施作进一步说明。本发明提出的超谱信号的快速编码方法，主要包括：

1．初始化阶段: 按照输入矢量第一维大小排序，按照平均分组→选择每组第一个训练矢量的步骤生成初始码书，计算各个训练矢量和初始码书中码字的误差，计算平均误差，并设置迭代次数。

2．训练阶段: 利用不等式排除算法快速搜索当前最佳匹配码字，将当前训练矢量划分到最佳匹配码字所对应的胞腔，训练完所有的训练矢量，以胞腔的质心更新当前最佳匹配码字。

3．中断检查: 若迭代次数满足要求，迭代结束，存储码书和所有矢量最匹配码字的索引，否则返回步骤2继续训练。

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

 构造矢量：读取超谱图像数据，选取每个波段中的一段图像块作为矢量量化的输入信源，根据图像尺寸和波段数构建2维矩阵，可采用美国喷气实验室提供的免费超谱源数据作为实验数据。

以512×614×224（512行，614列，224个波段）尺寸超谱图像为例，其中512×614代表一个波段中的图像尺寸（行列数），即512×614个像元，224是波段数，假设在每个波段中截取大小为256行和256列的图像块，这样得到图像块的3维矩阵P为256行、256列和224页（记为(256,256,224)）。将此3维矩阵P的各列首尾相接，转化为2维矩阵B，便于后续的矢量量化处理。这样P矩阵前两维的乘积数就是B矩阵的行维数，第3维数就是B矩阵的列维数，B是(256×256,224)的矩阵。

 对构造矢量进行哈达玛（Hadamard）变换。

对矩阵B进行哈达玛变换。矢量的哈达玛变换，是将多个行矢量组成的矩阵(假设矩阵维数为M×N)右乘哈达玛(维数为N×N)矩阵，完成对于行矢量的哈达玛变换，获得M×N维变换后的矢量。哈达玛矩阵中只有1和-1，具有计算简单，各行向量和列向量之间满足正交性的特点。

哈达玛矩阵是一组只含有{1,-1}的矩阵，对数据进行哈达玛变换不需要乘法运算，而只需简单的加减法即可完成。假设H_n为2ⁿ×2ⁿ的哈达玛矩阵，空域中输入矢量x为k维矢量，k=2ⁿ(n>0)，变换后矢量为X，有如下性质：

1）X₁=sx，其中，X₁是矢量X的第一维分量，sx为空域中输入矢量x的和值，即哈达玛域中矢量X的第一维分量是空域中矢量x的和值；

2），其中，Lx为哈达玛域中X范数的开方，即哈达玛域中矢量X的范数是空域中矢量x范数的k倍；

3）D(X,Yj)=kd(x,yj), Yj为码字yj的哈达玛变换，这里d(x,yj)是空域中矢量x和码字yj的欧式距离，D(X,Yj)是哈达玛变换域中矢量X和码字Yj的欧式距离，即哈达玛变换前后的欧式距离成倍数关系。

因此，在哈达玛域和空域中搜索最近邻码字是等价的，而且可以利用性质1来代替和值运算，达到减少计算复杂度的目的。哈达玛变换域算法要求矢量维数必须是2ⁿ。对于B矩阵，其列的维数是224，必需将B矩阵列维数扩充为256=2⁸维，可在B矩阵后面添加32列0矢量，将矩阵列维数变为256，扩充后的矩阵的每一行就可以进行哈达玛变换了。

 结合快速搜索算法和LBG算法生成最终码书。

本实例可采用如下所述排除法搜索匹配码字：

1）对于选定码字Yj，计算训练矢量X和码字Yj的欧式距离，根据公式：计算距离，如果距离满足条件，排除码字Yj。其中，D(X,Yj)为矢量X和码字Yj在哈达玛变换域的欧式距离，d_min代表空域中的矢量x和当前最匹配码字yp的欧式距离，D_min是空域中x的哈达玛变换X与空域中yp的哈达玛变换Yp的欧式距离， D_min为哈达玛变换域中的欧式距离，其中D_min=kd_min，k为矢量的维数。在哈达玛域中矢量量化等效于在空域中矢量量化。这样，只需算出D₁的值来代替复杂的D(X,Yj)计算，因为如果，那么肯定满足，顺利的排除码字Yj，因避免计算D(X,Yj)，可大大减少搜索的运算量。

2）对于采用步骤1）不能排除的码字，可以进一步使用不等式进行判断，即，如果满足，可以排除码字Yj，也避免计算D(X,Yj)。其中，X₁、Yj1分别表示矢量X和码字Yj的第一维分量，、分别表示矢量X的方差和码字Yj的方差。

3）经过上述不等式排除不必要的码字后，对于剩余码字，如果需要可以根据部分失真搜索算法(PDS)进一步提高排除效率。

在哈达玛变换域中，首先将所有的训练矢量按照第一维分量进行排序，并设定编码索引I，对于每一个矢量X，依据其索引值p，在码字中找出初始最佳匹配码字Yp，计算它与X之间的欧式距离，作为当前的最小距离D_min；然后，以Yp为中心，上下交替搜索其他码字Yj，利用上面所述的3步排除算法逐步排除不必要的码字；

最后，在剩余的码字中，计算它们与矢量X的欧式距离，选出最匹配码字作为X的编码码字，同时更新其索引值。本方法加快LBG算法的聚类速度，快速的生成最终码书。

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

假设码书大小为N，使用512×614×224规格的超谱图像数据，其每个波段中像元数据以两个字节的带符号整数记录，通过实验对本发明方案和基于空域的LBG超谱图像压缩算法进行比较。本发明方案的具体实施步骤如下：

初始化阶段：

步骤1：利用MATLAB函数库中的fopen函数和fread函数读取超谱数据源，得到超谱数据源的3维矩阵A(有512行，614列，224页，简记为（512,614,224）)；

步骤 2：截取矩阵A的前256行、256列和所有页，组成新的3维矩阵P(256,256,224)，将3维矩阵P转化为2维矩阵B，以便矢量量化处理。其变换方式是将3维矩阵P的各列首尾相接，形成256×256行、1列和224页的矩阵，等效为2维矩阵B (256×256行和224列的矩阵，简记为（256×256，224）)；

步骤 3：在2维矩阵B的每一行后添加32维的0矢量，形成维数为2的整数幂次（即256维）的行矢量，获得矩阵B'（256×256，256），为下一步的行变换做准备；

步骤 4：对B'进行哈达玛变换得到变换矩阵C(256×256,256)，将变换矩阵C按照第一列分量的大小重新进行排序，记录这时的排序索引I₁，获得行排序后的矩阵E，然后对矩阵E按行平均分组，依次选择每组第一行生成初始码书Y，这个初始码书自然是排了序的；

步骤 5：依据编码码书大小N，对于排序后的训练矢量（即待编码矢量），设定编码索引矩阵I(1,256×256)，并置训练矢量的初始编码索引为I=，为每个训练矢量分配索引值。将训练矢量按照初始编码索引的索引值归属到相应的胞腔中（具有相同索引值的矢量属于同一个胞腔），后续搜索时，将从这个胞腔开始，实施上下搜索，以减少搜索范围，降低运算复杂度；

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

训练阶段：

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

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

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

步骤 4：以码字YI (i)为中心，采用上下查找法，对于当前码字Yj，如果满足D₁=(X₁－ Yj1)²≥D_min，即矢量X的第一维分量X₁与码字Yj的第一维分量Yj1差值的平方大于当前最小欧式距离，就排除码字Yj，进入下列步骤a)或b)。否则进入步骤5；

a) 当X₁≥Yj1时，排除码字Yi，i=1,…,j，转入步骤7；

b)  当X₁≤Yj1时，排除码字Yi，i=j,…N，转入步骤7。（其中，j为当前搜索码字的标号，N表示码书的大小）

步骤 5：调用公式计算矢量X的方差与当前码字Yj的方差的差值平方与D₁的和值，如果满足D₂≥D_min，则排除码字Yj，转入步骤7，否则转入下一步进行判断；

步骤 6：采用PDS算法，算法模块根据公式计算训练矢量X与码字Yj的部分失真Dq，其中q为不大于k维的整数。并进一步判断，一旦满足Dq≥D_min，则排除码字Yj，转入步骤7。否则继续累计部分失真，直到q=k，此时若D(X,Yj)=Dq<D_min，令D_min=D，更新训练矢量的编码索引值I(i)为当前最近的码字的标号j，即I(i)= j；

步骤 7：返回步骤4进行下一个码字的匹配，直到搜索完所有的剩余码字，获得当前训练矢量X的最佳匹配码字和对应编码索引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))之商（V(I(i))/ U(I(i))），并且将这些质心替换原有码书中的码字，形成新的码书Y，然后计算更新码字后的平均量化误差，置迭代次数ite→ite+1。

步骤 9：根据平均量化误差的相对下降比例（如0.5%或0.1%）或直接设定的迭代次数ite（如实验中可设为10），判断是否满足停止迭代的要求，如果满足，所得Y作为最后的码书，否则转到步骤2进行下一轮训练。

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

 在接收端只需按照I₂中各个索引值查找码书Y，即可恢复经过矢量量化处理后的信源信息。

超谱图像压缩性能的好坏主要从图像的峰值信噪比(PSNR)，图像的压缩比(CR)，以及生成最终码书所需的运算复杂度进行评价，PSNR，CR的具体表达式如下：

PSNR=，而，其中Peaksignal表示超谱图像矩阵A中的最大值，MSE为均方误差，k为矢量维数224，y(i，j),b(i,j)表示编码后第i个矢量的第j维分量，b(i,j)表示编码前的第i个矢量的第j维分量，256×256表示矢量个数。

CR=，其中k为矢量维数224，N为码书尺寸，16代表每个像元用2个字节表示，256×256表示矢量个数。

以Lunar Lake和Low Altitude的第一场景的部分数据作为仿真测试信源，通过以上3个评价标准对LBG算法的超谱图像压缩和本发明方案进行了比较，N表示码书大小，表1列出了两种方法经过10次迭代的PSNR和压缩比，表2列出了平均每次迭代每个矢量距离计算的次数，表3 列出了Lunar Lake图像经过两种算法处理后的运算复杂度，表4列出了Low Altitude图像经过两种算法处理后的运算复杂度，

表1  经过10次迭代的PSNR和压缩比

表2  平均每次迭代每个矢量距离计算的次数

表3  Lunar Lake图像经过两种算法处理后的计算复杂度

表4  Low Altitude图像经过两种算法处理后的计算复杂度

 本发明结合排序初始化码书方案、基于LBG的聚类算法和高效率的快速码字搜索算法，快速和高质量的生成码书和编码，在提高处理数据质量的情况下，大大的降低运算复杂度，得到预期的压缩效果。从实验结果可以得到在码书尺寸为2048，对于Lunar Lake图像的部分数据的处理，量化编码的MSE=682，对于Low Altitude图像的部分数据的处理，量化编码的MSE=1542，这样每一个像元的平均量化误差基本在40以下，此时，压缩比为29：1，对于以两个字节表示一个像元数据的超谱数据，如此小的量化误差，使得这时的矢量量化有损压缩属于近似无损压缩，处理数据具有工程应用价值，从图3和图4也可以看到处理后的成像质量是相当好的。表格中给出了各种压缩比的情况，相应的处理数据可以满足从一般的浏览领域到实际的工程应用领域的各种需求。超谱图像所具有较高的波段数和大数据量的特点，目前针对超谱图像的压缩一直是关注的研究领域。不过由于种种原因，到目前为止一直没有形成一套成熟或标准的超谱图像压缩技术，本发明方案可以填补这个空缺，并作为进一步的工作方向。