基于降低映射预测残差的高光谱图像压缩算法*

李佳颖，朱文泉，孟繁蕴

(1.北京师范大学地理科学学部中药资源保护与利用北京市重点实验室，北京 100875；2.北京师范大学地理科学学部天然药物教育部工程研究中心，北京 100875；3.北京师范大学地理科学学部遥感科学与工程研究院，北京 100875)

1 引言

遥感技术在近20年得到了前所未有的发展，已成为对地观测的重要手段[1]。为满足各种分析和应用需要，遥感图像光谱分辨率不断提高，高光谱遥感应运而生[2]。高光谱数据兼具空间和光谱信息，可绘制地物光谱特征曲线，为定量遥感提供更丰富的光谱信息，为地物识别提供更高的识别能力及判别精度。海量的高光谱信息给数据存储和传输带来了挑战，因此研究高光谱图像压缩很有必要。考虑到高光谱图像获取难度大，具有完整保存的必要性[3]，无损压缩成为研究的重中之重[4]。

高光谱图像无损压缩技术可分为3类：基于预测的方法、基于变换的方法和基于矢量量化的方法[5]。基于预测的方法是最经典的图像压缩算法之一，2012年5月，国际空间数据系统咨询委员会CCSDS(Consultative Committee for Space Data Systems)针对多/高光谱图像提出了自适应的三维预测无损压缩标准CCSDS 123.0-B-1算法[6]。该算法是基于预测器的无损压缩，压缩系统由预测器和编码器组成(如图1所示)。预测器核心为FL(Fast Lossless)算法[7 - 9]，可减小数据冗余，得到预测残差；编码器对预测残差进行熵编码后输出。CCSDS 123.0-B-1标准提出的算法计算简单，复杂度低，对绝大多数多光谱及高光谱图像都能达到较好的预测效果[10]，是目前星载多/高光谱图像无损压缩的国际通用标准。但是，该算法对所有像素采用同样的预测模式，未针对像素位置进行区分；采用前向线性预测，未充分利用图像谱间相关性，压缩率仍有待提高。为此，何淑贤[11]提出了基于边缘检测的预测方法，首先对图像进行边缘检测，然后根据不同方向的边缘，对待预测像素周围像素点赋予不同权重，提高了预测准确性。王丽萍[12]分析了CCSDS 123.0-B-1标准中的参考波段数、权重分辨率、计数器初始化因子、计数器缩放区间等11个参数对压缩性能的影响，提出了一组优化参数，压缩性能明显提高。张雷等[13]提出了优化组合预测编码OCPC(Optimized Combination Processing Coding)算法，预测器根据像素点位置自适应选择预测模式，编码器采用Golomb-Rice编码，得到了较好的压缩效果。

由前述分析可知，在CCSDS 123.0-B-1算法中，如果映射预测残差越小，则其压缩后的码长就会越短，其压缩性能就会越好。基于此，本文提出一种降低映射预测残差的RMPR(Reduction of Mapped Prediction Residual)算法，通过提高预测值准确度和优化残差映射器来减小映射预测残差，从而提高压缩性能。

Figure 1 Compression system图1 压缩系统

2 数据与方法

2.1 数据

本文的高光谱测试数据为CCSDS测试图像，来源于CCSDS官网(https://public.ccsds.org/default.aspx)。表1所示为10幅测试图像的具体信息，图2所示为测试图像的部分区域。

Table 1 Image information

2.2 方法

2.2.1 CCSDS 123.0-B-1算法简介

CCSDS 123.0-B-1算法的压缩流程如图3所示，下面重点介绍该算法的预测器和编码器。

(1)预测器。

Figure 2 Test images图2 测试图像

预测器通过当前像元周围像元和相邻谱带像元产生像元预测值，在此过程中去除了谱间相关和空间相关，从而减小冗余；像元预测值和真实像元值之间产生残差，由于残差有正负之分，通过映射器区分正负，且将正残差和负残差均映射为较小正整数输入编码器。

首先计算局部和σz,y,x，其中，x,y,z分别表示高光谱图像的列数、行数、波段数，用户可选择计算同列局部和或周边局部和，两者的预测点不同：同列局部和用同列像元Sz,y-1,x作预测点进行计算(如图4a所示)；周边局部和用周围4个像元Sz,y-1,x-1、Sz,y-1,x、Sz,y-1,x+1、Sz,y,x-1作预测点进行计算(如图4b所示)。

Figure 3 Compression flow chart of CCSDS 123.0-B-1 algorithm图3 CCSDS 123.0-B-1算法的压缩流程

针对局部差向量Uz(t)的计算，预测器中采用了2种模式：简化模式和完整模式。2种模式都采用向前线性单向预测。简化模式中，局部差向量取决于当前所预测谱带的先前谱带；完整模式中，局部差向量取决于当前预测谱带和先前谱带，当前预测谱带的局部差由当前像元周围3个像元决定(如图5所示)。每1个局部差向量对应1个权重向量Wz(t)，权值的更新自适应于每一像元，当前像元的权值依赖于上一像元。

Figure 4 Selection of samples for calculating local sums图4 计算局部和的预测点选取

Figure 5 Selection of current band samples for calculating local differences图5 计算局部差的当前谱带像元选取

残差映射器是阈值选择映射[14]，表达式如式(1)所示，其中θz(t)是阈值。当|Δx(t)|>θz(t)时，阈值选择映射可以减小映射预测残差δz(t)，有利于图像压缩[15]。

(1)

(2)编码器。

编码器对映射预测残差δz(t)进行了熵编码[16]，目的是消除残差统计的冗余[17]。CCSDS 123.0-B-1算法采用Golomb-Rice编码[18,19]，其基本思想是将映射预测残差除以已知整数，该整数是2的次幂，分别存储商和余数。编码器还设置了累加器和计数器以确定每个映射预测残差的编码长度，两者在编码中自适应迭代。这种编码方式适用于小数字出现概率高于大数字出现概率的情况，因为它使用较短码长编码小数字，用较长码长编码大数字。

2.2.2 改进的RMPR算法

由于映射预测残差δz(t)越小，压缩后码长越短，压缩性能越好，基于此，本文提出RMPR算法，通过提高预测值准确度和优化残差映射器来减小映射预测残差，从而提高压缩性能。其中，预测值准确度的提高主要通过预测点改进和预测方向改进2种途径。

(1)预测点改进。

If (x==1,2,3‖y==1,2,3‖x==Nx-1,Nx-2,Nx-3)

return getColumnSum(z,y,x)

else

return getNeighborSum(z,y,x)

(2)预测方向改进。

考虑到高光谱图像较强的谱间相关性，RMPR算法采用双向线性预测提高预测准确度，即将待预测波段B个前向波段和B个后向波段作为参考波段。过少的参考波段数会带来较大误差，甚至导致预测发散；过多的参考波段数会导致计算量增加，收敛条件更严格[21]，综合考虑，B值推荐为3。

(3)残差映射器改进。

预测残差有正负之分，为了便于后续编码，一般将负值变为正值，同时要将正负数映射为奇偶数，在解码时通过奇偶区分正负。CCSDS 123.0-B-1算法采用阈值选择映射，其中θz(t)是阈值，当预测较为准确，预测残差的数值|Δz(t)|较小，即|Δz(t)|≤θz(t)时，正残差被映射为正奇数：2|Δz(t)|-1，负残差被映射为正偶数：2|Δz(t)|；当预测不准导致预测残差的数值|Δz(t)|过大，即 |Δz(t)|>θz(t)时，预测残差被映射为|Δz(t)|与阈值θz(t)的和，由于阈值θz(t)小于|Δz(t)|。因此，|Δz(t)|与阈值θz(t)的和小于2倍的|Δz(t)|。相比于直接映射为2倍的|Δz(t)|，采用这种映射方式得到的映射预测残差δz(t)有所减小。 此种情况下映射预测残差δz(t)减小的具体数值取决于θz(t)，即要想使δz(t)有所降低，可减小与之相加的阈值θz(t)。RMPR算法采用式(2)作为残差映射器，m取值为0.1，可进一步减小映射预测残差δz(t)。

(2)

(4)RMPR算法流程。

综上，RMPR算法的处理流程如下所示：

①输入待压缩图像并初始化权重向量。

⑤e用式(2)计算映射预测残差δz(t)，输入编码器后进行Golomb-Rice编码。

⑥循环至下一像元继续进行②～⑤步，直至最后一个像元编码完成，补0后输出码流。

(5)算法测试环境及统计检验。

本文用Java语言实现CCSDS 123.0-B-1、对CCSDS 123.0-B-1仅改进预测点、对CCSDS 123.0-B-1仅改进预测方向、对CCSDS 123.0-B-1仅改进残差映射器及RMPR算法，在Eclipse平台上对高光谱图像进行测试。测试版本为Oxygen.2 Release (4.7.2)，设置虚拟机参数为-Xms2048m-Xmx16384m，测试环境为Windows 8.1 64位操作系统、4 GB内存、Intel Core i5-4210M CPU。压缩性能用比特每像素BPP(Bits Per Pixel)统计，压缩时间单位为毫秒，各算法之间差异的显著性用配对样本T检验进行比较，P值小于0.05表示差异显著。

2.2.3 其它用于比较的经典压缩算法

为了将本文RMPR算法与经典的高光谱图像无损压缩算法进行压缩性能方面的比较，本文还参考并比较了以下经典和新型的压缩算法。其中，RMPR、CCSDS 123.0-B-1、基于查找表的LUT(Look-Up Tables)算法[22]、联合图像专家小组发布的JPEG-LS(Joint Photographic Experts Group-Lossless or near Lossless)标准[23]和基于上下文的自适应无损压缩M-CALIC(Context-based Adaptive Lossless Image Coding)[24]均为基于预测的压缩算法；可逆的卡胡南-洛维变换RKLT(Reversible Karhunen-Loève Transform)[25]、主成分分析PCA(Principal Component Analysis)[26]、5/3离散小波变换DWT(Discrete Wavelet Transform)[27]为基于变换的压缩算法。本文采用Spectral Transform软件[28]实现RKLT，去除测试图像谱间冗余，再用Kakadu软件[29]实现JPEG2000图像压缩，得到RKLT+JPEG2000算法的压缩结果，其余算法结果来源于文献[8,30,31]。

Figure 6 Partial test images before and after compression图6 压缩前后的部分测试图像

3 实验结果

在压缩时间均无显著差异的情况下，改进预测点算法的压缩性能显著优于CCSDS 123.0-B-1算法的(如表2所示)，BPP降低了0.007 0(如表3所示，其中*表示各改进算法相较于CCSDS 123.0-B-1算法在压缩性能或压缩时间上有显著差异)，这是因为图像边缘映射预测残差降低明显；改进预测方向的压缩算法的性能也显著优于CCSDS 123.0-B-1算法的(如表2所示)，BPP降低0.044 0(如表3所示)，因为改进预测方向后大部分像素映射预测残差都有所下降；改进残差映射器的算法的压缩性能无显著提高(如表2所示)，Yellowstone Calibrated Sc00和Mt.St.Helens Uncalibrated 2幅测试图像压缩效果甚至无变化(如表3所示)，原因是预测器带来了较好的预测效果，预测差值绝对值较小，无法落入|Δz(t)|>θz(t)的区间；RMPR算法显著优于原始算法(如表2所示)，BPP降低0.065 3(如表3所示)，主要是预测精度有所提高、像素映射预测残差降低、压缩码长缩短。

Table 2 P values of paired-samples T test

为进一步验证RMPR算法的压缩性能，对RMPR算法压缩后的图像进行解压缩。压缩前后的图像在视觉效果上没有差别(如图6所示)，进一步在Matlab软件中对比原图像及解压缩后图像可知，两者像素值完全一致，RMPR算法可实现高光谱图像的无损压缩。

Table 3 Compression test results for different algorithms

4 讨论

4.1 RMPR算法的压缩性能与压缩效率

10幅测试图像压缩结果表明，优化前CCSDS 123.0-B-1算法可实现4.449 5BPP的压缩效果，RMPR算法可实现4.3842BPP的压缩效果，压缩比提升幅度为0.5%～3.8%，压缩性能显著优于 CCSDS 123.0-B-1算法的，且算法的改进并未显著降低压缩效率。当前针对CCSDS 123.0-B-1的改进算法有优化参数的MHDC(the Multispectral & Typerspectral Data Compression)[12]、OCPC[13]、基于V型扫描的CCSDS[32]和基于边缘检测的预测算法[11]等。优化参数下的MHDC算法在压缩效率不变的情况下，压缩比提高了0.7%～3.6%[12]；OCPC算法压缩比提升幅度为2.4%～5.2%，压缩效率未讨论[13]；基于V型扫描的CCSDS算法在压缩复杂度下降的基础上，压缩比提升了1.1%～4.2%[32]；基于边缘检测的预测算法压缩比提高了1.5%～11.2%，压缩时间无显著变化[11]。

将本文算法与经典的高光谱图像无损压缩算法进行压缩性能方面的比较，结果如表4所示，其中，“—”表示相关文献未针对该数据进行实验。对于AVIRIS高光谱数据，RMPR算法、PCA算法能够达到相对较优的压缩性能,两者均具有良好的去相关性能；DWT 5/3算法针对HYPERION数据表现出较好的压缩性能；从平均BPP来看，RKLT+JPEG2000算法的压缩性能最优，其次是RMPR算法和PCA算法。RKLT+JPEG2000算法综合表现良好，但并未针对单个图像实现最优压缩；在基于变换的压缩算法中，PCA算法能够实现良好的光谱去相关，但是协方差矩阵的计算、特征向量的提取大大增加了运算复杂度。综上，与经典高光谱压缩算法相比，RMPR算法并不逊色。

4.2 RMPR算法的参数选择

RMPR算法涉及2个参数的选取：预测波段数B和残差映射器参数m。

预测波段数B决定了参与预测的波段信息，进而影响压缩性能。将B数值分别定为1,3,5,…,13,15对10幅测试图像进行压缩。随着预测波段数B的不断增加，Cuprite Uncalibrated、Lake Monona Uncalibrated、Mt.St.Helens Uncalibrated 3幅图像的压缩性能越来越好，BPP减小，且B取值为1～9时，BPP下降速率较大(如图7所示)；Hawaii Uncalibrated Sc01、Yellowstone Uncalibrated Sc00、Yellowstone Calibrated Sc03、Yellowstone Uncalibrated Sc03 5幅图像压缩后BPP呈现先下降再上升的趋势，压缩性能由高到低，由于图像复杂度不同，BPP上升的临界点不同，但B取值为1～5时BPP均呈现下降趋势(如图7所示)；Maine Uncalibrated Sc10、Yellowstone Calibrated Sc00 2幅图像的压缩性能越来越差，BPP不断增加(如图7所示)。从实验结果来看，B越大，能利用的谱间相关信息越多，越利于压缩性能的提高，但是B的数值过大会导致待预测波段与参考波段相隔较远，两者相关性下降，压缩性能降低，且算法复杂度增加(如图8所示)。综合考虑压缩性能和压缩时间，预测波段数B取1～5为宜。

Table 4 Comparison of compression performance between RMPR algorithm and classical algorithms

Figure 7 Compression performance vs. the number of bands for prediction B图7 压缩性能随预测波段数B的变化图

在|Δz(t)|>θz(t)的情况下，残差映射器参数m决定了映射预测残差的数值，进而影响压缩性能。将m数值分别定为0.1,0.2,…,0.9,1对10幅测试图像进行压缩。随着残差映射参数m的不断减小，大部分测试图像的BPP呈下降趋势，且在m为0.1时下降速率达到最大，此时压缩性能最优(如图9所示)。因此,本文将残差映射器参数m取值为0.1。

Figure 8 Compression time vs. the number of bands for prediction B图8 压缩时间随预测波段数B的变化图

Figure 9 Image of compression performance with residual mapper parameter m图9 压缩性能随残差映射器参数m的变化图

4.3 RMPR算法的应用潜力及改进方向

RMPR算法在不影响图像解压缩且不影响压缩效率的前提下，达到了较好的预测效果，压缩性能显著优于原算法，且算法简单，易于硬件实现，适用于星载多/高光谱的数据压缩应用。

但是本算法仍有需要进一步完善的地方：一是未对地物进行分类，仅利用当前像元邻近像元作为预测点，改进方法是对图像进行预处理，用同类像元进行预测，以提高预测精度；二是无论双向线性预测还是单向预测，都是用相邻波谱段作为参考波段，但是有些待预测波段的最佳参考波段可能不与之相邻或相隔较远[33]，会使预测产生偏差，解决方法是为每个波段寻找最佳参考波段[34]。

5 结束语

本文对CCSDS 123.0-B-1算法的预测器进行了优化，提出了RMPR算法。RMPR算法从预测点、预测方向、残差映射器3个方面进行了优化，通过减小映射预测残差来缩短压缩码长，提高压缩性能。实验结果表明，在压缩时间无显著差异的情况下，RMPR算法压缩性能显著优于CCSDS 123.0-B-1算法的，可应用于资源有限的星上无损压缩，下一步可从地物分类和寻找最佳参考波段两方面进一步完善该算法。