2016, Vol. 62
文章信息
- 殷彬, 裴加军, 蒙广元, 匡娇娇, 王慧芳, 吴敏渊
- YIN Bin, PEI Jiajun, MENG Guangyuan, KUANG Jiaojiao, WANG Huifang, WU Minyuan
- CCSDS多光谱图像无损压缩算法的FPGA实现
- The FPGA Implementation of CCSDS Lossless Multispectral Image Compression Algorithm
- 武汉大学学报(理学版), 2016, 62(6): 531-536
- Journal of Wuhan University(Natural Science Edition), 2016, 62(6): 531-536
- http://dx.doi.org/10.14188/j.1671-8836.2016.06.007
-
文章历史
- 收稿日期：2016-02-24
引用本文
|
2. 上海航天电子技术研究所，上海 201109
2. Shanghai Academy of Spaceflight Technology, Shanghai 201109, China
随着航天研究技术的不断发展，光谱图像的分辨率不断提高，卫星上产生的海量数据对有限的信道资源造成了巨大的压力，给数据的存储和传输都带来了很大的困难，因此进行有效的数据压缩就尤为重要.而多光谱图像信息量较大，如果对多光谱图像进行有损压缩，很多细节的信息将会丢失，导致图像模糊及光谱曲线的连续性变差[1]，为了保留图像的细节信息，我们更希望大容量的多光谱图像的压缩既要高保真又要保证实时性，因此研究多光谱图像的无损压缩具有十分重要的意义.
目前，多光谱图像的无损压缩技术大致可分为以下三类：基于预测的压缩方法、基于变换的压缩方法和基于矢量量化的压缩方法[2].其中，基于预测的压缩方法提出时间早，发展比较成熟，目前在无损压缩编码领域已经得到了广泛的应用.通常来说，预测方法可以分为三种，它们分别为一维预测、二维预测和三维预测.一维预测是需要与待测的像素点在同一扫描行中前面相邻的一个或者几个像素点；二维预测除了用同一扫描行的前几个相邻的像素点以外，还需要该像素点前一行或者前几行的若干个像素点来预测；三维预测除了需要二维预测所需的像素点外，还要前一帧或者前一个波段对应位置的像素点.CCSDS (空间数据系统咨询委员会)在2012年5月发布的第一个专门针对星载多光谱/高光谱图像无损压缩的标准CCSDS 123.0-B-1[3]，该标准是基于三维自适应预测的无损压缩算法[4].该算法复杂度低，压缩率高，目前已经成为星载多光谱/高光谱图像无损压缩的国际通用标准.本文采用该标准算法，基于FPGA搭建的硬件平台，实现了多光谱无损压缩算法.
1 算法分析CCSDS无损压缩算法与传统无损压缩算法的最大区别就是充分利用了多光谱图像谱段像素之间存在的谱间冗余，既使用谱段内，也使用谱间相邻的像素值去预测当前的像素值，通过对像素值和它的预测残差进行熵编码来实现数据的压缩.该算法由预测器和熵编码两部分组成，如图 1所示.预测器采用三维自适应线性预测的方法，对图像数据进行去相关处理，熵编码模块负责对去相关后的数据进行压缩编码.
|
| 图 1 压缩算法模型 Figure 1 Compression algorithm model |
预测就是根据之前时刻的样本值对当前的样本值进行一种模型化的运算[5].为了保证无损压缩，预测必须是可逆的.图 2表示三维预测模型.其中z表示波段，x表示图像宽度，y表示图像高度.预测器先计算待预测点sz, y, x的局部和，以此计算出局部差和方向局部差，得到当前谱段的局部差向量.再利用前p个谱段对应像素点的局部差向量与权重向量线性组合就得到了预测采样值和缩放采样预测值[6].每预测完成一个点，都要对当前权重向量进行更新以便计算下个点的预测值.最后通过一个映射函数把预测残差映射成无符号类型.
|
| 图 2 三维预测模型 Figure 2 Three-dimensional prediction model |
局部和、局部差主要有两种计算模式：全模式邻域导向和简化模式列导向.对帧成像仪产生的并校准过的图像，采用全模式邻域导向的局部和计算方法，此方法往往产生更少的压缩图像数据量.对推扫成像仪产生的未校准图像，采用简化模式列导向的局部和计算方法将获得更高的压缩比.根据项目需求，本文采用基于全模式邻域导向的预测器对图像进行去相关处理.
邻域导向的局部和σz, y, x是待预测像素点sz, y, x的西边(W)像素点sz, y, x-1、西北边(NW)像素点sz, y-1, x-1、北边(N)像素点sz, y-1, x和东北边(NE)像素点sz, y-1, x+1的加权和，具体见公式(1):
|
(1) |
全模式需要计算p个中央局部差和3个方向局部差，p表示预测所用的波段数量，可以自定义，只需保证p小于当前待预测像素所在波段之前的波段数即可.p个中央局部差分别是待预测像素的前p个波段与待预测像素同位置处像素的值与其局部和之差，例如待预测像素前一个波段与之同位置处像素sz-1, y, x的中央局部差dz-1, y, x的计算方法见公式(2)所示，其余波段的以此类推.
|
(2) |
3个方向局部差dz, y, xN，dz, y, xW，dz, y, xNW分别是北边、西边和西北边三个像素点与待预测像素点局部和之差，如图 2，相关的计算公式见(3)、(4)和(5).
|
(3) |
|
(4) |
|
(5) |
由此得到全模式下的差向量Uz(t).
|
(6) |
全模式对应的权重向量如下式，
|
(7) |
对于每个波段，当前波段的三个方向局部差对应的权重初始值如下，
|
(8) |
各个波段中央局部差的权重初始值如下，
|
(9) |
其中，Ω为用户可自定义的权重值分辨率，取值范围在[4,19], 从(8~9)式可以看出，距离待预测像素点越远的波段，计算局部差时所占的权重越小.
CCSDS算法实现自适应编码的重点就是权重向量的更新.对每个待预测像素点进行预测后，权重向量的各个权重值都会更新一次，从而保证对下一个待预测像素点的预测可以更加准确.更新的依据就是待预测像素点灰度值与预测值之差，即预测误差ez(t).对该值进行运算，若该值为正数，说明预测值偏小，适当增大权重值；若该值为负数，说明预测值偏大，适当减少权重值，从而使得预测值逐渐地接近待预测像素点的值.权重更新表达如(10)式所示:
|
(10) |
其中，裁剪函数clip的计算方法见(11)式
|
(11) |
权重更新因子ρ(t)决定了权重值每次调整的步长，该值越小，权重值调整的越快，对图像数据统计结果的适应性越好，但压缩效果的稳定性越差.
1.1.3 预测值的计算差向量的每个元素乘以相应的权重值并累加后得到加权和，称之为预测中央局部差
CCSDS 123.0-B-1标准算法中的熵编码指的是长度受限的Golomb编码，对符合几何对称性并且其分布特征是指数衰减概率分布的非负整数信号，Golomb编码具有较好的编码效果.此方法的使用有利于硬件的编解码.其中定义的码字长度不能超过最大码字长度Umax(比特)，Umax是整数，由用户输入.根据预测器输出的预测残差计算出kz(t)和uz(t)，其中kz(t)表示除数的指数，uz(t)表示相除后的商, δz(t)代表预测器输出的预测残差, 具体公式如下.
|
(12) |
当uz(t) < Umax时，码字由uz(t)个0、1个1和映射预测残差的低kz(t)位二进制数组成；否则，由Umax个0和D位的映射预测残差组成.由于该编码算法较为常见，具体计算方法不再赘述.
2 FPGA设计与实现本文采用的平台是Stratix Ⅱ[7]系列的FPGA芯片，考虑到在FPGA上实现CCSDS算法的复杂度和中间数据存储操作的便利，在本文中设定输入的多光谱图像的光谱波段数为4，单个波段的宽度是1 024个像素，高度也为1 024个像素.对于其他的不同波段数不同宽高的多光谱图像，我们只需要在上位机上更改相应的参数, 然后发送到FPGA配置相应的寄存器即可.
2.1 压缩系统的结构设计CCSDS无损压缩系统的功能流程图如图 3所示，为了保证原始图像数据的实时输入，未压缩的数据由PC上位机程序控制，先通过以太网发送到FPGA，FPGA转存到同步动态随机存储器(synchronous dynamic random access memory, SDRAM)，顶层控制模块读写SDRAM的速率为100 MHz，由于SDRAM的数据线位宽是16位，因此可以保证未压缩数据的输入速率达到200 MB/s.当原始图像数据发送完毕后，顶层控制模块将会读取前面写入的数据然后交给数据预处理模块处理，之后数据被送到压缩算法模块，经过预测器和熵编码的处理后输出不定长的码流数据送入缓冲区先入先出队列(first input first out put, FIFO).码流组包模块从码流缓冲区FIFO取出数据后进行组包，其输出的定长数据通过顶层控制模块写入SDRAM.最后顶层控制模块从SDRAM读出压缩后的定长码流数据, 并通过千兆网送回PC端进行解压验证.
|
| 图 3 CCSDS多光谱无损压缩系统功能流程图 Figure 3 Functional flow diagram of CCSDS multispectral image lossless compression system |
数据预处理模块实现的功能是根据输入的原始图像数据经过处理后输出符合预测模块输入要求的数据，输出的数据包括当前像素点与其西边、西北边、北边和东北边的像素点.该模块的原理流程图，如图 4所示.当使能信号为高电平时，原始图像数据连续进入BZ移位寄存器，其长度是4个字(由于当前用于预测的光谱波段数为4，同理对于BZ_Previous移位寄存器)，位宽为14位，经过移位后按照波段按像元交叉格式(band interleaved by line format, BIP)的顺序并行输出Zt、Zt_1、Zt_2、Zt_3作为预测波段局部差矢量的计算.Zt同时进入大容量FIFO和BZ_Previous移位寄存器，其中大容量FIFO存储当前处理像素点的上一行像素点.BZ_Previous移位寄存器经过4个字长的移位寄存器后的输出作为当前处理像素西边的像素点.大容量FIFO输出经过缓冲输出后得到Pixel_W、Pixel_NE、Pixel_N和Pixel_NW，其中B1、B2、B3移位寄存器的长度均为4，位宽均为14位.
|
| 图 4 数据预处理模块原理流程图 Figure 4 Block flow diagram of data preprocessing module |
压缩算法模块是本文的核心模块，包括预测器和熵编码器两个模块.其具体实现流程如图 5所示.
|
| 图 5 压缩算法实现流程图 Figure 5 Compression algorithm realization flow diagram |
预测模块包括预测值计算模块和映射预测残差计算模块.考虑到预测值计算模块和映射预测残差模块存在较多的数值运算，同时权重更新模块需要用到预测值计算模块的输出作为输入.为了最大程度提高系统的压缩速度，整个模块共包含两级流水线.在预测模块加入一级流水线后，预测值计算模块和映射预测残差模块分别需要一个时钟周期来完成.插入流水线后，预测模块的最大频率较未插入流水线提升将近一倍.
预测模块主要完成的功能包括计算局部和、局部差、矢量乘法和预测样值估算等.映射预测残差模块的主要功能是把预测后的差值全部转为非负整数，以适应后面的熵编码模块的输入要求.
2.3.2 熵编码模块熵编码模块的主要功能是实现数据的无损编码.该模块由编码采样和计数器累加器更新两个子模块组成，其中编码采样模块由乘法器、移位寄存器、加法器、除法器和求对数运算器组成.如果整个熵编码模块在一个时钟周期内完成，则这条路径的延时将会很长，很难满足实时性的要求.在这里编码采样子模块和计数器累加器子模块采用并行计算的方法，有效降低了关键路径的延时，提高了系统的频率.对于求对数则采用查表的方法计算，这是因为kz(t)是不大于D-2的最大非负整数，采用查表只需要使用D-2路的多路选择器，资源消耗没有明显增多，但是却能大大降低对数运算带来的路径延时和复杂度.计数器累加器更新模块由比较器、FIFO (长度为波段数目，位宽等于像素点位宽)和其他逻辑电路组成.熵编码实现的原理框图如图 6所示.
|
| 图 6 熵编码实现原理流程图 Figure 6 The flow diagram of entropy coding |
码流组包的作用是将熵编码模块输出的不定长的码流数据进行整合，最终按照规定的要求进行码流的输出[8].本文组包后输出的定长数据位为8位，位序为升序.码流组包功能框图如图 7所示，组包的过程如图 8所示.先从输入的不定长码流数据中取出有效数据，然后把有效数据送到寄存器的空余部分，当寄存器被填满时把数据送到发送端.
|
| 图 7 码流组包模块功能流程图 Figure 7 Functional flow diagram of data stream packet module |
|
| 图 8 不定长组包原理图 Figure 8 Process of variable length data stream packet |
本文采用Stratix Ⅱ系列的FPGA芯片作为设计验证平台，搭建了验证系统.验证采用的数据源是多幅不同宽高和波段数的多光谱图像.输入图像的宽度和高度都不超过4 095个像素，波段数不超过15个，位深度不超过14位.为了验证系统的正确性，本文测试了一系列的多光谱图像.图像来源包括用户给定的多光谱图像和Landsat卫星公开的多光谱图像.图 9是其中4幅测试用图.
|
| 图 9 多光谱图像测试图 (a)山脉；(b)农田；(c)城市；(d)海边 Figure 9 Multispectral image test charts (a) Mountains; (b) Farmland; (c) City; (d) Seaside |
其中前两幅图像的尺寸为1 024×1 024，波段数为8，位深度为10位.后两幅图像的尺寸为1 024×1 024，波段数为4，位深度为8位.这四幅多光谱图像的像素存储格式均为BIP格式.解压FPGA压缩后上传的数据并和原数据比较，结果完全一致.证明本系统可以准确无误地实现对多光谱图像进行无损压缩.表 1是CCSDS 3种压缩算法的性能对比，CCSDS 123.0的压缩性能明显好于CCSDS 121.0[9]，略好于CCSDS 122.0[10].对系统进行编译后查看Quartus II静态时序分析工具(Time Quest Timing Analyzer)后获知，压缩算法模块在布局布线后的数据吞吐率可达66 M sample/s，相比NASA JPL实验室实现的33 M sample/s快了一倍.在CPU主频2.53 GHz，内存2GB的PC机上，压缩一幅8MB的图像平均时间为1.7 s，满足项目的实时性要求.
| 测试图 | CCSDS121.0-B-1 | CCSDS122.0-B-1 | CCSDS123.0-B-1 |
| 山脉 | 1.526 93 | 1.706 00 | 1.952 60 |
| 农田 | 1.088 26 | 1.270 00 | 1.384 56 |
| 城市 | 1.421 08 | 1.208 25 | 1.501 18 |
| 海边 | 2.518 95 | 1.621 53 | 2.595 56 |
本文基于CCSDS 123.0-B-1的算法原理，设计了CCSDS多光谱无损压缩系统的整体结构.随后对预测模块、熵编码模块等进行分析并设计，最后在FPGA上实现了该算法.通过一系列的测试图验证了本系统的正确性与合理性，从而能够满足航天应用的需求.
但是本系统仍有需要完善的地方.一是系统的抗误码能力较差，当压缩后的码流中的某些位出现错误时后面的码流解压容易出错，造成恢复出来的图像存在较多错误，下一步可采用分段压缩的方法提高压缩算法模块的抗误码能力；二是本文实现的压缩算法模块部分所占用的FPGA逻辑资源较多，为了提高压缩速率，采用了面积换速度的方法，所以在资源占用上还可以进一步进行优化.
| [1] | 朱学伟, 李晓辉, 朱博. CCSDS压缩算法对高光谱数据质量的影响研究[J]. 遥感信息 , 2013, 28 (4) : 29–36 ZHU X W, LI X H, ZHU B. Influence of the hyperspectral data quality based on CCSDS compression algorithm[J]. Remote Sensing Information , 2013, 28 (4) : 29–36 |
| [2] | 万建伟, 粘永健, 苏令华, 等. 高光谱图像压缩技术研究进展[J]. 信号处理 , 2010, 26 (9) : 1397–1407 WAN J W, ZHAN Y J, SU L H, et al. Research progress on hyperspectral imagery compression technique[J]. Signal Processing , 2010, 26 (9) : 1397–1407 |
| [3] | CCSDS. CCSDS123.0-B-1. Lossless Multispectral & Hyperspectral Image Compression[S]. Washington D C: CCSDS Secretariat, 2012. |
| [4] | ARANKI N, BAKHSHI A, KEYMEULEN D, et al. Fast and adaptive lossless on-board hyperspectral data compression system for space applications[C]//2009 IEEE Aerospace conference, New York: IEEE Press, 2009:1-8. |
| [5] | 邓海涛.多光谱图像自适应无损压缩算法研究与FPGA实现[D].海口:海南大学, 2014. DENG H T. Research and Implementation of Adaptive Lossless Compression Algorithm for Multi spectral Images Based on FPGA[D]. Haikou: Hainan University, 2014(Ch). http://cdmd.cnki.com.cn/Article/CDMD-10589-1014240043.htm |
| [6] | 王丽萍.基于预测的高光谱图像无损压缩技术研究[D].西安:西安电子科技大学, 2014. WANG L P. Research on Prediction-based Lossless Compression for Hyperspectral Images[D]. Xi'an: Xidian University, 2014(Ch). http://cdmd.cnki.com.cn/Article/CDMD-10701-1014331212.htm |
| [7] | Altera.SII5V1-4.5. Stratix Ⅱ Device Handbook[S]. California Los Angeles: Altera, 2011. |
| [8] | 雷震霖.空间图像CCSDS压缩算法研究与FPGA实现[D].大连:大连理工大学, 2007. LEI Z L. CCSDS Space Image Compression Algorithm and FPGA Implementation[D]. Dalian: Dalian University of Technology, 2007(Ch). http://cdmd.cnki.com.cn/Article/CDMD-10141-2008034577.htm |
| [9] | CCSDS.CCSDS 121.0-B-1. Lossless multispectral & hyperspectral image compression[S]. Washington D C: CCSDS Secretariat, 1997. |
| [10] | CCSDS.CCSDS 122.0-B-1. Image Data Compression[S]. Washington D C: CCSDS Secretariat, 2005. |