CN103501437A - 一种基于分形和h.264的高光谱图像压缩方法 - Google Patents

Abstract

本发明提出了一种基于分形和H.264的高光谱图像压缩方法。首先将高光谱图像立方数据体转换为YUV格式的视频，送入编码器。高光谱视频的I帧，使用基于宏块平坦度的H.264快速帧内预测模式判别方法进行帧内预测，以去除高光谱图像的空间相关性；高光谱视频的P帧进行块运动估计/补偿分形编码，以去除高光谱图像的谱间相关性。P帧编码宏块的各种树状结构分块首先由MSE准则在参考帧中寻找最佳匹配块，确定各块的迭代函数系统系数，然后比较所有块划分方式的率失真代价，将率失真代价最小的块划分方式作为最终的帧间编码模式，记录最终的分形参数。最后，I帧和P帧的残差数据经过DCT变换、量化之后由熵编码CABAC进行编码写入码流，P帧的分形参数也进行CABAC熵编码。

Description

一种基于分形和H.264的高光谱图像压缩方法

技术领域：

本发明属于高光谱图像压缩领域，针对高光谱图像中存在的空间相关性与谱间相关性，提出一种基于分形和H.264的高光谱图像压缩方法，在保证图像质量的前提下，大大加快了高光谱图像的压缩速度，提高了压缩比。

背景技术：

高光谱技术是21世纪遥感技术的发展前沿和当今遥感界关注的焦点之一，它把传统的二维成像遥感技术和光谱技术有机地结合在一起，在用成像系统获得被测物的空间信息的同时，通过光谱仪系统把被测物的辐射分解成不同波长的谱辐射，能在一个光谱区间内获得每个像素几十甚至几百个连续的载波段信息。高光谱遥感技术己经在地质学、生态学、大气研究、土壤研究很多领域得到了成功的应用，显示出很大的潜力和广阔的发展前景。

高光谱图像在地面图像二维信息的基础上，增加了第三维光谱信息。高光谱数据是一个光谱图像的立方体，其空间图像维描述地表二维空间特征，其光谱维揭示图像每一像元的光谱曲线特征。在高光谱图像数据的空间图像维存在空间冗余，光谱维存在谱间冗余。因此高光谱图像压缩与普通图像压缩的区别是既要考虑去除空间相关性，又要考虑去除谱间相关性。目前，常用的高光谱图像压缩方法分为三大类：基于预测的方法、基于变换的方法和基于矢量量化的方法。例如陈雨时与张晔等提出了基于线性模型最优预测的压缩方案，利用递归双向预测的思想，通过建立谱带间的线性模型，得出了信噪比下的最优预测（参见陈雨时,张晔,张钧平.基于线性模型最优预测的高光谱图像压缩[J].南京航空航天大学学报,2007,39(3):368-372.）。Penna B通过采用点数抽取策略降低KLT的计算量，并将这种改进的KLT用于高光谱图像谱间去相关，而重建图像质量并未受到明显影响（参见Penna B,Tillo T,MagliE,et al.Transform coding techniques for lossy hyperspectral data compression[J].IEEETransactions on Geoscience and Remote Sensing,2007,45(5):1408-1421.）。Shen-En Qian提出了一种快速的矢量量法算法来提高码书的产生的效率，在该算法无需进行全搜索，大大地降低了运算的复杂度（参见Shen-En Qian.Hyperspectral data compression using a fast vectorquantization algorithm[J].IEEE Transactions on Geoscience and Remote Sensing.2004,42(8):1791-1798.）。除了上述传统的编码方法外，近年来分形编码也被应用于高光谱图像的压缩，如夏丽丽设计了一个三维分形压缩编码算法，在去除图像空间相关性的同时去除了波段图像之间的相关性，得到了较好的压缩效果（参见夏丽丽.基于分形理论的高光谱图像压缩算法的研究[J].计算机与数字工程,2011,39(9):132-135.）。

本发明结合了预测编码、分形编码、DCT变换编码、CABAC熵编码，利用基于宏块平坦度的H.264快速帧内预测编码方法降低高光谱图像的空间相关性，利用分形编码降低其谱间相关性，并且对编码残差数据进行DCT变换、量化，用CABAC熵编码写入码流，分形参数也进行CABAC熵编码，实现了高光谱图像的有效压缩。

发明内容：

本发明提出了一种基于分形和H.264的高光谱图像压缩方法。首先将高光谱图像立方数据体转换为YUV格式的视频，送入编码器。高光谱视频的I帧，使用基于宏块平坦度的H.264快速帧内预测模式判别方法进行帧内预测，以去除高光谱图像的空间相关性；高光谱视频的P帧进行块运动估计/补偿分形编码，以去除高光谱图像的谱间相关性。P帧编码宏块的各种树状结构分块首先由MSE准则在参考帧中寻找最佳匹配块，确定各块的迭代函数系统系数，然后比较所有块划分方式的率失真代价，将率失真代价最小的块划分方式作为最终的帧间编码模式，记录最终的分形参数。最后，I帧和P帧的残差数据经过DCT变换、量化之后由熵编码CABAC进行编码写入码流，P帧的分形参数也进行CABAC熵编码。

一种基于分形和H.264的高光谱图像压缩方法，其特征在于实现步骤如下：

步骤一：将高光谱图像立方体数据转换成YUV格式的视频；

步骤二：若为高光谱视频的I帧（第一帧必须为I帧，其它帧可以设置是否为I帧），使用基于宏块平坦度的H.264快速帧内预测模式判别方法进行帧内预测模式选择，预测宏块，完成所有块的预测后可以得到I帧的预测帧。通过原始帧和预测帧之差得到编码端残差帧，转入步骤四编码残差；若为高光谱视频的P帧，转到步骤三；

步骤三：若为高光谱视频的P帧，依次对当前帧的所有宏块进行分形编码。在参考帧中的搜索窗内对当前宏块进行块匹配，子块的位置作为父块的起始搜索点，父块的大小与子块的大小相同。对每个宏块进行树状结构分块，即分块从大到小可分为16×16、16×8、8×16、8×8，8×8往下（亚宏块分割）可以分为8×4，4×8，4×4。首先对编码宏块的各种树状结构分块由MSE准则寻找最佳匹配块，确定各块的迭代函数系统系数即IFS系数；然后比较所有块划分方式的率失真代价；最后将率失真代价最小的块划分方式作为最终的帧间编码模式。记录最终的IFS系数，转入步骤五得到本块的重建块。如果当前帧所有的宏块都已编码完毕，所有的重建块组成重建图像（即下一帧的参考帧），通过原始图像与重建图像之差得到编码端残差图像，转到步骤四编码残差。所述搜索窗为在参考帧中的矩形搜索区域；所述IFS系数包括父块与子块的位置偏移，即运动矢量(x,y)和比例因子s、偏移因子o；

步骤四：残差图像的数据经过DCT、量化之后的系数一方面进行Zig-Zag扫描，然后用熵编码CABAC进行编码写入码流；另一方面经过反量化、反DCT变换后得到解码端残差帧。由预测帧和解码端残差帧之和得到重建帧（即下一帧的参考帧）。如果是P帧则还要对所有IFS系数进行CABAC熵编码。判断当前帧是否为最后一帧，如果是最后一帧结束编码；否则，返回步骤二继续处理下一帧图像；

步骤五：通过保存的IFS系数代入解码方程式

r_i=s·d_i+o                     (1)

计算得到预测值，由原始块和预测块之差得到编码端残差块，编码端残差块经过DCT变换、量化、反量化和反DCT变换得到解码端残差块，再由预测块和解码端残差块之和得到重建块。转入步骤三编码深度图序列P帧下一宏块。式中r_i为预测块的像素值，d_i为参考帧相应父块的像素值。

所述步骤一中将高光谱图像立方体数据转换成YUV格式的视频包括以下两个步骤：

1)高光谱图像数据归一化至[0,255]区间内的整数。高光谱图像数据表示地物的光谱反射率，其位深度一般为16比特/像素，需要转换成8比特/像素。转换方法为首先将光谱反射率的值归一化至[0,255]区间，然后将其转换为无符号8比特的整数，最终的数值为[0,255]区间的整数；

2)将高光谱图像数据的每个波段作为一帧，每个波段的归一化后的反射系数作为该帧的亮度元素（Y分量），色度分量Cb和Cr均设为中值128。循环执行步骤2)直到高光谱图像的所有波段处理完毕。所有帧即组成一个YUV格式的视频文件。

所述步骤二中基于宏块平坦度的H.264快速帧内预测模式判别方法包括以下五个步骤：

1)计算当前编码宏块二维灰度直方图，并记录二维灰度直方图中z轴的最大值z_max。

图像的二维灰度直方图是基于图像像素灰度及像素邻域灰度均值的联合分布构成的三维描述图。设图像f(x,y)的大小为M×N，灰度级为L。由f(x,y)采用k×k点阵平滑得到的邻域平均灰度图像为g(x,y)，记为

$g(x,y)=\frac{1}{k^2}\underset{m=-\frac{k}{2}}{\overset{\frac{k}{2}}{\mathrm{\Σ}}}\underset{n=-\frac{k}{2}}{\overset{\frac{k}{2}}{\mathrm{\Σ}}}f(x+m,y+n)---\left(2\right)$

式中，1≤x+m≤M，1≤y+n≤N，k取奇数。

则g(x,y)的图像大小和灰度级与f(x,y)相同。由f(x,y)和g(x,y)可以构成一个二元组(i,j)，每个二元组属于二维平面上的一个点。定义二维直方图为N(i,j)，表示图像f(x,y)像素灰度值为i，且像素邻域平均灰度图像的灰度值为j时，二维点(i,j)出现的次数（i,j=0,1,...,L-1）。二维灰度直方图中的z轴即为N(i,j)，代表二维点(i,j)出现的次数；

2)设定上限阈值Th_high和下限阈值Th_low，Th_high和Th_low均为[1,256]间的整数；

3)若z_max≥Th_high，认为宏块平坦，去除Intra4×4模式，选择Intra16×16模式，并将率失真开销最小的模式作为最优帧内预测模式；同时更新上限值

否则转入步骤4)；

4)若z_max≤Th_low，认为宏块细节丰富，去除Intra16×16模式，选择Intra4×4模式，并将率失真开销最小的模式作为最优帧内预测模式；同时更新下限值

否则转入步骤5)；

5)若Th_low<z_max<Th_high，认为宏块平坦度特征不显著，采用H.264标准帧内预测算法。

所述步骤三中P帧宏块分形编码包括以下十个步骤：

1)首先以16×16宏块为单位在参考帧中的搜索窗内对当前宏块进行块匹配。子块的位置作为父块的起始搜索点，父块的大小与子块的大小相同。按照全搜索方式遍历整个搜索窗，寻找匹配误差MSE最小的父块位置。记录该位置迭代函数系统系数即IFS系数（包括运动矢量(x,y)，比例因子s，偏移因子o），并计算该模式的率失真开销cost1；

2)将该16×16宏块划分为2个16×8的小块，对每个小块均按照全搜索方式遍历相应的整个搜索窗，寻找匹配误差MSE最小的相应父块位置。记录这两个位置的IFS系数，并计算两个块的率失真开销之和cost2；

3)将该16×16宏块划分为2个8×16的小块，对每个小块均按照全搜索方式遍历相应的整个搜索窗，寻找匹配误差MSE最小的相应父块位置。记录这两个位置的IFS系数，并计算两个块的率失真开销之和cost3；

4)将该16×16宏块划分为4个8×8的小块，对每个小块均按照全搜索方式遍历相应的整个搜索窗，寻找匹配误差MSE最小的相应父块位置。记录这4个位置的IFS系数，并分别计算这4个块的率失真开销cost4_1、cost4_2、cost4_3、cost4_4及四个率失真开销之和cost4；

5)比较步骤1)到步骤4)中的率失真开销cost1、cost2、cost3、cost4，若最小值为cost4，则需进一步将每个8×8块继续划分，转到步骤6)；否则率失真开销最小的划分方式即为该宏块最终的P帧编码模式，保留相应的IFS系数，结束；

6)将第1个8×8块首先划分为2个8×4小块，对每个小块均按照全搜索方式遍历相应的整个搜索窗，寻找匹配误差MSE最小的相应父块位置。记录这两个位置的IFS系数，并计算这两个块的率失真开销之和cost5；

7)将第1个8×8块划分为2个4×8小块，对每个小块均按照全搜索方式遍历相应的整个搜索窗，寻找匹配误差MSE最小的相应父块位置。记录这两个位置的IFS系数，并计算这两个块的率失真开销之和cost6；

8)将第1个8×8块划分为4个4×4小块，对每个小块均按照全搜索方式遍历相应的整个搜索窗，寻找匹配误差MSE最小的相应父块位置。记录这4个位置的IFS系数，并计算这4个块的率失真开销之和cost7；

9)比较cost4_1、cost5、cost6、cost7，将其中的最小值对应的块划分方式作为该8×8块最终的编码模式，保留相应的IFS系数；

10)依次对16×16宏块的每个8×8块按照步骤6)至步骤9)的方法选择最终的编码模式，保留相应的IFS系数。

其中匹配误差MSE的计算如下：

$\mathrm{MSE}(R,D)=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{[r_i-(s\&CenterDot;d_i+o)]}^2=\frac{1}{N}[\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{r}_i^2+s(s\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{d}_i^2-2\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{r}_i{d}_i+2o\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{d}_i^2)+o(N\&CenterDot;o-2\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{r}_i^2)]---\left(3\right)$

式中s、o的计算式分别为：

$s=\frac{[N\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{r}_i{d}_i-\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{r}_i\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{d}_i]}{[N\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{d}_i^2-{\left(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{d}_i\right)}^2]}---\left(4\right)$

$o=\frac{1}{N}[\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{r}_i-s\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{d}_i]---\left(5\right)$

其中，r_i是当前块的各点像素值，d_i是匹配块的各点像素值，N为当前块的像素个数，s为尺度因子，o为偏移因子。

率失真开销的计算如下：

cost=MSE(R,D)+λ_MODE·Rate(x,y,s,o,MODE|QP)(6)

其中，R代表当前块，D代表匹配块，MODE表示宏块划分方式，QP表示量化参数，λ_MODE是与QP有关的拉格朗日参数。

本发明所提出的基于分形和H.264的高光谱图像压缩方法的优点在于：

(1)本方法首先将高光谱图像立方体数据进行了格式转换，转换成YUV格式的视频，这样就可以像编码普通视频一样编码高光谱图像；

(2)本方法采用基于宏块平坦度的H.264快速帧内预测编码降低高光谱图像中存在的空间相关性。采用二维灰度直方图表征宏块平坦度特征，并设置了双阈值，对平坦宏块去除Intra4×4模式；对细节丰富的宏块去除Intra16×16模式；对平坦度特征不明显的宏块采用H.264标准算法。阈值动态迭代更新，使判断阈值可以根据已重建宏块平坦度的波动而实时调整。既保证了高光谱图像解压缩的质量，又加快了压缩速度；

(3)本方法采用分形编码降低高光谱图像中存在的谱间相关性。块匹配支持16×16、16×8、8×16、8×8、8×4、4×8、4×4共7种块划分方式，每种块划分方式首先寻找匹配误差MSE最小的匹配块，然后比较所有模式的率失真代价，最后将率失真代价最小的模式作为最优的分形编码模式。提高了分形编码的精度。

附图说明：

图1为本发明一种基于分形和H.264的高光谱图像压缩方法流程图；

图2为高光谱图像“DC Mall”的第17波段灰度子图像；

图3为本发明一种基于分形和H.264的高光谱图像压缩方法步骤二中基于宏块平坦度的H.264快速帧内预测模式判别方法流程图

图4为本发明一种基于分形和H.264的高光谱图像压缩方法步骤三中P帧宏块分形编码模式选择流程图；

图5为经本发明一种基于分形和H.264的高光谱图像压缩方法压缩再重建后的“DC Mall”高光谱图像第17波段灰度图；

图6(a)为本发明一种基于分形和H.264的高光谱图像压缩方法与Lucana等人方法对“DCMall”高光谱图像进行压缩编码的峰值信噪比的对比图；

图6(b)为本发明一种基于分形和H.264的高光谱图像压缩方法与Lucana等人方法对“DCMall”高光谱图像进行压缩的比特数的对比图；

图6(c)为本发明一种基于分形和H.264的高光谱图像压缩方法与Lucana等人方法对“DCMall”高光谱图像进行压缩的时间的对比图。

具体实施方式：

下面将结合附图对本发明方法作进一步的详细说明。

本发明提出了一种基于分形和H.264的高光谱图像压缩方法。首先将高光谱图像立方数据体转换为YUV格式的视频，送入编码器。高光谱视频的I帧，使用基于宏块平坦度的H.264快速帧内预测模式判别方法进行帧内预测，以去除高光谱图像的空间相关性；高光谱视频的P帧进行块运动估计/补偿分形编码，以去除高光谱图像的谱间相关性。P帧编码宏块的各种树状结构分块首先由MSE准则在参考帧中寻找最佳匹配块，确定各块的迭代函数系统系数，然后比较所有块划分方式的率失真代价，将率失真代价最小的块划分方式作为最终的帧间编码模式，记录最终的分形参数。最后，I帧和P帧的残差数据经过DCT变换、量化之后由熵编码CABAC进行编码写入码流，P帧的分形参数也进行CABAC熵编码。

如附图1所示，一种基于分形和H.264的高光谱图像压缩方法流程图。以弗吉尼亚州光谱信息技术应用中心提供的华盛顿哥伦比亚特区购物中心（DC Mall）高光谱图像为例，截取256×1024像素子图像。去除其中云覆盖率较大的0.9到1.4μm范围的波段，最终保留的波段数目为191。附图2是高光谱图像“DC Mall”的第17波段灰度子图像。

步骤一：将高光谱图像立方体数据转换成YUV格式的视频。具体步骤如下：

第1步、高光谱图像数据归一化至[0,255]区间内的整数。高光谱图像数据表示地物的光谱反射率，其位深度一般为16比特/像素，需要转换成8比特/像素。转换方法为首先将光谱反射率的值归一化至[0,255]区间，然后将其转换为无符号8比特的整数，最终的数值为[0,255]区间的整数；

第2步、将高光谱图像数据的每个波段作为一帧，每个波段的归一化后的反射系数作为该帧的亮度元素（Y分量），色度分量Cb和Cr均设为中值128。循环执行第2步直到高光谱图像的所有波段处理完毕。所有帧即组成一个YUV格式的视频文件。

步骤二：判断当前编码高光谱视频帧是否为I帧（第一帧必须为I帧，其它帧可以设置是否为I帧），如果是I帧，使用基于宏块平坦度的H.264快速帧内预测模式判别方法进行帧内预测模式选择，预测宏块。附图3为基于宏块平坦度的H.264快速帧内预测模式判别方法流程图。具体包括以下5个步骤：

第1步、计算当前编码宏块二维灰度直方图，并记录二维灰度直方图中z轴的最大值z_max。

图像的二维灰度直方图是基于图像像素灰度及像素邻域灰度均值的联合分布构成的三维描述图。设图像f(x,y)的大小为M×N，灰度级为L。由f(x,y)采用k×k点阵平滑得到的邻域平均灰度图像为g(x,y)，记为

$g(x,y)=\frac{1}{k^2}\underset{m=-\frac{k}{2}}{\overset{\frac{k}{2}}{\mathrm{\&Sigma;}}}\underset{n=-\frac{k}{2}}{\overset{\frac{k}{2}}{\mathrm{\&Sigma;}}}f(x+m,y+n)---\left(7\right)$

式中，1≤x+m≤M，1≤_y+n≤N，k取奇数。

则g(x,y)的图像大小和灰度级与f(x,y)相同。由f(x,y)和g(x,y)可以构成一个二元组(i,j)，每个二元组属于二维平面上的一个点。定义二维直方图为N(i,j)，表示图像f(x,y)像素灰度值为i，且像素邻域平均灰度图像的灰度值为j时，二维点(i,j)出现的次数（i,j=0,1,...,L-1）。二维灰度直方图中的z轴即为N(i,j)，代表二维点(i,j)出现的次数；

第2步、设定上限阈值Th_high（经验值设为175）和下限阈值Th_low（经验值设为100），Th_high和Th_low均为[1,256]间的整数；

第3步、若z_max≥Th_high，认为宏块平坦，去除Intra4×4模式，选择Intra16×16模式，并将率失真开销最小的模式作为最优帧内预测模式；同时更新上限值

否则转入第4步；

第4步、若z_max≤Th_low，认为宏块细节丰富，去除Intra16×16模式，选择Intra4×4模式，并将率失真开销最小的模式作为最优帧内预测模式；同时更新下限值否则转入第5步；

第5步、若Th_low<z_max<Th_high，认为宏块平坦度特征不显著，采用H.264标准帧内预测算法。

完成所有块的预测后可以得到I帧的预测帧。通过原始帧和预测帧之差得到编码端残差帧，转入步骤四编码残差；如果为高光谱视频的P帧，转到步骤三。

步骤三：如果为高光谱视频的P帧，依次对当前帧的所有宏块进行分形编码。在参考帧中的搜索窗内对当前宏块进行块匹配，子块的位置作为父块的起始搜索点，父块的大小与子块的大小相同。对每个宏块进行树状结构分块，即分块从大到小可分为16×16、16×8、8×16、8×8，8×8往下（亚宏块分割）可以分为8×4，4×8，4×4。首先对编码宏块的各种树状结构分块由MSE准则寻找最佳匹配块，确定各块的迭代函数系统系数即IFS系数；然后比较所有块划分方式的率失真代价；最后将率失真代价最小的块划分方式作为最终的帧间编码模式。附图4为P帧宏块分形编码模式选择流程图，具体包括以下10个步骤：

第一步、首先以16×16宏块为单位在参考帧中的搜索窗内对当前宏块进行块匹配。子块的位置作为父块的起始搜索点，父块的大小与子块的大小相同。按照全搜索方式遍历整个搜索窗，寻找匹配误差MSE最小的父块位置。记录该位置迭代函数系统系数即IFS系数（包括运动矢量(x,y)，比例因子s，偏移因子o），并计算该模式的率失真开销cost1；

第二步、将该16×16宏块划分为2个16×8的小块，对每个小块均按照全搜索方式遍历相应的整个搜索窗，寻找匹配误差MSE最小的相应父块位置。记录这两个位置的IFS系数，并计算两个块的率失真开销之和cost2；

第三步、将该16×16宏块划分为2个8×16的小块，对每个小块均按照全搜索方式遍历相应的整个搜索窗，寻找匹配误差MSE最小的相应父块位置。记录这两个位置的IFS系数，并计算两个块的率失真开销之和cost3；

第四步、将该16×16宏块划分为4个8×8的小块，对每个小块均按照全搜索方式遍历相应的整个搜索窗，寻找匹配误差MSE最小的相应父块位置。记录这4个位置的IFS系数，并分别计算这4个块的率失真开销cost4_1、cost4_2、cost4_3、cost4_4及四个率失真开销之和cost4；

第五步、比较第一步到第四步中的率失真开销cost1、cost2、cost3、cost4，若最小值为cost4，则需进一步将每个8×8块继续划分，转到第六步；否则率失真开销最小的划分方式即为该宏块最终的P帧编码模式，保留相应的IFS系数，结束；

第六步、将第1个8×8块首先划分为2个8×4小块，对每个小块均按照全搜索方式遍历相应的整个搜索窗，寻找匹配误差MSE最小的相应父块位置。记录这两个位置的IFS系数，并计算这两个块的率失真开销之和cost5；

第七步、将第1个8×8块划分为2个4×8小块，对每个小块均按照全搜索方式遍历相应的整个搜索窗，寻找匹配误差MSE最小的相应父块位置。记录这两个位置的IFS系数，并计算这两个块的率失真开销之和cost6；

第八步、将第1个8×8块划分为4个4×4小块，对每个小块均按照全搜索方式遍历相应的整个搜索窗，寻找匹配误差MSE最小的相应父块位置。记录这4个位置的IFS系数，并计算这4个块的率失真开销之和cost7；

第九步、比较cost4_1、cost5、cost6、cost7，将其中的最小值对应的块划分方式作为该8×8块最终的编码模式，保留相应的IFS系数；

第十步、依次对16×16宏块的每个8×8块按照第六步至第九步的方法选择最终的编码模式，保留相应的IFS系数，结束。

其中匹配误差MSE的计算如下：

$\mathrm{MSE}(R,D)=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{[r_i-(s\&CenterDot;d_i+o)]}^2=\frac{1}{N}[\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{r}_i^2+s(s\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{d}_i^2-2\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{r}_i{d}_i+2o\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{d}_i^2)+o(N\&CenterDot;o-2\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{r}_i^2)]---\left(8\right)$

式中s、o的计算式分别为：

$s=\frac{[N\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{r}_i{d}_i-\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{r}_i\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{d}_i]}{[N\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{d}_i^2-{\left(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{d}_i\right)}^2]}---\left(9\right)$

$o=\frac{1}{N}[\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{r}_i-s\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{d}_i]---\left(10\right)$

其中，r_i是当前块的各点像素值，d_i是匹配块的各点像素值，N为当前块的像素个数，s为尺度因子，o为偏移因子。

率失真开销的计算如下：

cost=MSE(R,D)+λ_MODE·Rate(x,y,s,o,MODE|QP)                   (11)

其中，R代表当前块，D代表匹配块，MODE表示宏块划分方式，QP表示量化参数，λ_MODE是与QP有关的拉格朗日参数。

由上述步骤得到P帧宏块最终的编码模式，记录最终的IFS系数，转入步骤五得到本块的重建块。如果当前帧所有的宏块都已编码完毕，所有的重建块组成重建图像（即下一帧的参考帧），通过原始图像与重建图像之差得到编码端残差图像，转到步骤四编码残差。所述搜索窗为在参考帧中的矩形搜索区域；所述IFS系数包括父块与子块的位置偏移，即运动矢量(x,y)和比例因子s、偏移因子o。

步骤四：残差图像的数据经过DCT、量化之后的系数一方面进行Zig-Zag扫描，然后用熵编码CABAC进行编码写入码流；另一方面经过反量化、反DCT变换后得到解码端残差帧。由预测帧和解码端残差帧之和得到重建帧（即下一帧的参考帧）。如果是P帧则还要对所有IFS系数进行CABAC熵编码。判断当前帧是否为最后一帧，如果是最后一帧结束编码；否则，返回步骤二继续处理下一帧图像。

步骤五：通过保存的IFS系数代入解码方程式

r_i=s·d_i+o                     (12)

计算得到预测值，由原始块和预测块之差得到编码端残差块，编码端残差块经过DCT变换、量化、反量化和反DCT变换得到解码端残差块，再由预测块和解码端残差块之和得到重建块。转入步骤三编码高光谱视频P帧下一宏块。式中r_i为预测块的像素值，d_i为参考帧相应父块的像素值。

本方法选择Visual C++6.0作为所述方法的实现平台，CPU为Intel Core^TM2Duo T8300，2.4GHz主频，内存大小为2G，对高光谱图像“DC Mall”191个波段立方体数据进行了基于分形和H.264的高光谱图像压缩实验，转换成YUV格式的视频文件后，总帧数等于总波段数为191帧，每帧图像大小为256×1024像素，编码图像组结构为IPPP…，I帧编码周期为50，即第一帧编码为I帧，之后每隔49个P帧编码1个I帧，搜索范围为±7，量化参数QP为20。附图5为经本发明一种基于分形和H.264的高光谱图像压缩方法压缩再重建后的“DC Mall”高光谱图像第17波段灰度图。

分别采用Lucana等人方法（参见Lucana Santos,Sebastián López,Gustavo M.Callicó,et al.Performance evaluation of the H.264/AVC video coding standard for lossy hyperspectral imagecompression[J].IEEE Journal of Selected Topics in Applied Earth Observations and RemoteSensing,2012,5(2):451-461.）和本发明方法对“DC Mall”高光谱图像进行压缩的峰值信噪比的对比图如附图6(a)所示；分别采用Lucana等人方法和本发明方法对“DC Mall”高光谱图像进行压缩的比特数的对比图如附图6(b)所示；分别采用Lucana等人方法和本发明方法对“DC Mall”高光谱图像进行压缩的时间的对比图如附图6(c)所示。

分别采用Lucana等人方法和本发明方法对“DC Mall”高光谱图像191个波段压缩性能均值对比结果如表1所示。其中△PSNR、△比特率、△压缩时间的定义如下：

△PSNR=PSNR_OURS-PSNR_Lucana                     (13)

△比特率=（比特率_OURS-比特率_Lucana）/比特率_Lucana           (14)

△压缩时间=（压缩时间_OURS-压缩时间_Lucana）/压缩时间_Lucana        (15)

深度图序列	△PSNR/dB	△比特率	△压缩时间
				DC Mall	0.67	-90.50％	-90.00％

表1Lucana等人方法和本发明方法对“DC Mall”高光谱图像191个波段压缩性能均值对比

从附图6和表1可以看出，本发明方法与传统的国际视频编码标准H.264标准测试模型JM18.1方法相比，峰值信噪比PSNR平均提高0.67dB，编码码流比特率平均降低17.73%，压缩时间平均减少90.00%，表现出优越的压缩性能。这是因为本发明方法在深度图序列I帧编码上，采用了改进的基于宏块平坦度的H.264快速帧内预测编码方法，去除高光谱图像的空间相关性；在深度图序列P帧编码上，采用了块运动估计/补偿分形编码算法，去除高光图像的谱间相关性。

Claims (4)

1.一种基于分形和H.264的高光谱图像压缩方法，其特征在于如下步骤：

步骤一：将高光谱图像立方体数据转换成YUV格式的视频；

步骤二：若为高光谱视频的I帧（第一帧必须为I帧，其它帧可以设置是否为I帧），使用基于宏块平坦度的H.264快速帧内预测模式判别方法进行帧内预测模式选择，预测宏块，完成所有块的预测后可以得到I帧的预测帧。通过原始帧和预测帧之差得到编码端残差帧，转入步骤四编码残差；若为高光谱视频的P帧，转到步骤三；

步骤三：若为高光谱视频的P帧，依次对当前帧的所有宏块进行分形编码。在参考帧中的搜索窗内对当前宏块进行块匹配，子块的位置作为父块的起始搜索点，父块的大小与子块的大小相同。对每个宏块进行树状结构分块，即分块从大到小可分为16×16、16×8、8×16、8×8，8×8往下（亚宏块分割）可以分为8×4，4×8，4×4。首先对编码宏块的各种树状结构分块由MSE准则寻找最佳匹配块，确定各块的迭代函数系统系数即IFS系数；然后比较所有块划分方式的率失真代价；最后将率失真代价最小的块划分方式作为最终的帧间编码模式。记录最终的IFS系数，转入步骤五得到本块的重建块。如果当前帧所有的宏块都已编码完毕，所有的重建块组成重建图像（即下一帧的参考帧），通过原始图像与重建图像之差得到编码端残差图像，转到步骤四编码残差。所述搜索窗为在参考帧中的矩形搜索区域；所述IFS系数包括父块与子块的位置偏移，即运动矢量(x,y)和比例因子s、偏移因子o；

步骤四：残差图像的数据经过DCT、量化之后的系数一方面进行Zig-Zag扫描，然后用熵编码CABAC进行编码写入码流；另一方面经过反量化、反DCT变换后得到解码端残差帧。由预测帧和解码端残差帧之和得到重建帧（即下一帧的参考帧）。如果是P帧则还要对所有IFS系数进行CABAC熵编码。判断当前帧是否为最后一帧，如果是最后一帧结束编码；否则，返回步骤二继续处理下一帧图像；

步骤五：通过保存的IFS系数代入解码方程式

r_i=s·d_i+o                     (1)

计算得到预测值，由原始块和预测块之差得到编码端残差块，编码端残差块经过DCT变换、量化、反量化和反DCT变换得到解码端残差块，再由预测块和解码端残差块之和得到重建块。转入步骤三编码深度图序列P帧下一宏块。式中r_i为预测块的像素值，d_i为参考帧相应父块的像素值。

2.根据权利要求1所述一种基于分形和H.264的高光谱图像压缩方法，其特征在于：所述步骤一中将高光谱图像立方体数据转换成YUV格式的视频包括以下两个步骤：

1)高光谱图像数据归一化至[0,255]区间内的整数。高光谱图像数据表示地物的光谱反射率，其位深度一般为16比特/像素，需要转换成8比特/像素。转换方法为首先将光谱反射率的值归一化至[0,255]区间，然后将其转换为无符号8比特的整数，最终的数值为[0,255]区间的整数；

2)将高光谱图像数据的每个波段作为一帧，每个波段的归一化后的反射系数作为该帧的亮度元素（Y分量），色度分量Cb和Cr均设为中值128。循环执行步骤2)直到高光谱图像的所有波段处理完毕。所有帧即组成一个YUV格式的视频文件。

3.根据权利要求1所述一种基于分形和H.264的高光谱图像压缩方法，其特征在于：所述步骤二中基于宏块平坦度的H.264快速帧内预测模式判别方法包括以下五个步骤：

1)计算当前编码宏块二维灰度直方图，并记录二维灰度直方图中z轴的最大值z_max。

图像的二维灰度直方图是基于图像像素灰度及像素邻域灰度均值的联合分布构成的三维描述图。设图像f(x,y)的大小为M×N，灰度级为L。由f(x,y)采用k×k点阵平滑得到的邻域平均灰度图像为g(x,y)，记为

$g(x,y)=\frac{1}{k^2}\underset{m=-\frac{k}{2}}{\overset{\frac{k}{2}}{\mathrm{\&Sigma;}}}\underset{n=-\frac{k}{2}}{\overset{\frac{k}{2}}{\mathrm{\&Sigma;}}}f(x+m,y+n)---\left(2\right)$

式中，1≤x+m≤M，1≤y+n≤N，k取奇数。

则g(x,y)的图像大小和灰度级与f(x,y)相同。由f(x,y)和g(x,y)可以构成一个二元组(i,j)，每个二元组属于二维平面上的一个点。定义二维直方图为N(i,j)，表示图像f(x,y)像素灰度值为i，且像素邻域平均灰度图像的灰度值为j时，二维点(i,j)出现的次数（i,j=0,1,...,L-1）。二维灰度直方图中的z轴即为N(i,j)，代表二维点(i,j)出现的次数；

2)设定上限阈值Th_high和下限阈值Th_low，Th_high和Th_low均为[1,256]间的整数；

3)若z_max≥Th_high，认为宏块平坦，去除Intra4×4模式，选择Intra16×16模式，并将率失真开销最小的模式作为最优帧内预测模式；同时更新上限值

否则转入步骤4)；

4)若z_max≤Th_low，认为宏块细节丰富，去除Intra16×16模式，选择Intra4×4模式，并将率失真开销最小的模式作为最优帧内预测模式；同时更新下限值

否则转入步骤5)；

5)若Th_low<z_max<Th_high，认为宏块平坦度特征不显著，采用H.264标准帧内预测算法。

4.根据权利要求1所述一种基于分形和H.264的高光谱图像压缩方法，其特征在于：所述步骤三中P帧宏块分形编码包括以下十个步骤：

1)首先以16×16宏块为单位在参考帧中的搜索窗内对当前宏块进行块匹配。子块的位置作为父块的起始搜索点，父块的大小与子块的大小相同。按照全搜索方式遍历整个搜索窗，寻找匹配误差MSE最小的父块位置。记录该位置迭代函数系统系数即IFS系数（包括运动矢量(x,y)，比例因子s，偏移因子o），并计算该模式的率失真开销cost1；

2)将该16×16宏块划分为2个16×8的小块，对每个小块均按照全搜索方式遍历相应的整个搜索窗，寻找匹配误差MSE最小的相应父块位置。记录这两个位置的IFS系数，并计算两个块的率失真开销之和cost2；

3)将该16×16宏块划分为2个8×16的小块，对每个小块均按照全搜索方式遍历相应的整个搜索窗，寻找匹配误差MSE最小的相应父块位置。记录这两个位置的IFS系数，并计算两个块的率失真开销之和cost3；

4)将该16×16宏块划分为4个8×8的小块，对每个小块均按照全搜索方式遍历相应的整个搜索窗，寻找匹配误差MSE最小的相应父块位置。记录这4个位置的IFS系数，并分别计算这4个块的率失真开销cost4_1、cost4_2、cost4_3、cost4_4及四个率失真开销之和cost4；

5)比较步骤1)到步骤4)中的率失真开销cost1、cost2、cost3、cost4，若最小值为cost4，则需进一步将每个8×8块继续划分，转到步骤6)；否则率失真开销最小的划分方式即为该宏块最终的P帧编码模式，保留相应的IFS系数，结束；

6)将第1个8×8块首先划分为2个8×4小块，对每个小块均按照全搜索方式遍历相应的整个搜索窗，寻找匹配误差MSE最小的相应父块位置。记录这两个位置的IFS系数，并计算这两个块的率失真开销之和cost5；

7)将第1个8×8块划分为2个4×8小块，对每个小块均按照全搜索方式遍历相应的整个搜索窗，寻找匹配误差MSE最小的相应父块位置。记录这两个位置的IFS系数，并计算这两个块的率失真开销之和cost6；

8)将第1个8×8块划分为4个4×4小块，对每个小块均按照全搜索方式遍历相应的整个搜索窗，寻找匹配误差MSE最小的相应父块位置。记录这4个位置的IFS系数，并计算这4个块的率失真开销之和cost7；

9)比较cost4_1、cost5、cost6、cost7，将其中的最小值对应的块划分方式作为该8×8块最终的编码模式，保留相应的IFS系数；

10)依次对16×16宏块的每个8×8块按照步骤6)至步骤9)的方法选择最终的编码模式，保留相应的IFS系数。

其中匹配误差MSE的计算如下：

$\mathrm{MSE}(R,D)=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{[r_i-(s\&CenterDot;d_i+o)]}^2=\frac{1}{N}[\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{r}_i^2+s(s\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{d}_i^2-2\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{r}_i{d}_i+2o\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{d}_i^2)+o(N\&CenterDot;o-2\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{r}_i^2)]---\left(3\right)$

式中s、o的计算式分别为：

$s=\frac{[N\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{r}_i{d}_i-\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{r}_i\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{d}_i]}{[N\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{d}_i^2-{\left(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{d}_i\right)}^2]}---\left(4\right)$

$o=\frac{1}{N}[\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{r}_i-s\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{d}_i]---\left(5\right)$

其中，r_i是当前块的各点像素值，d_i是匹配块的各点像素值，N为当前块的像素个数，s为尺度因子，o为偏移因子。

率失真开销的计算如下：

cost=MSE(R,D)+λ_MODE·Rate(x,y,s,o,MODE|QP)                 (6)

其中，R代表当前块，D代表匹配块，MODE表示宏块划分方式，QP表示量化参数，λ_MODE是与QP有关的拉格朗日参数。

Images