CCSDS多光谱图像无损压缩算法的FPGA实现

CCSDS多光谱图像无损压缩算法的FPGA实现 The FPGA Implementation of CCSDS Lossless Multispectral Image Compression Algorithm 1 first-author 殷彬，男，硕士生，现从事计算机图像处理方面的研究. E-mail:595846525@qq.com 殷彬，男，硕士生，现从事计算机图像处理方面的研究. E-mail: 595846525@qq.com 595846525@qq.com 2 1 1 1 1 corresp 吴敏渊，E-mail:wmy@whu.edu.cn 吴敏渊，E-mail:wmy@whu.edu.cn wmy@whu.edu.cn 武汉大学 电子信息学院，湖北 武汉 430072 武汉大学 电子信息学院，湖北 武汉 430072 School of Electronic Information, Wuhan University, Wuhan 430072, Hubei, China School of Electronic Information, Wuhan University, Wuhan 430072, Hubei, China 上海航天电子技术研究所，上海 201109 上海航天电子技术研究所，上海 201109 Shanghai Academy of Spaceflight Technology, Shanghai 201109, China Shanghai Academy of Spaceflight Technology, Shanghai 201109, China 本文基于CCSDS 123.0-B-1多光谱图像无损压缩算法，采用Stratix II系列的现场可编程门阵列(FPGA)芯片，按算法的功能模块分别在FPGA上进行RTL级的设计，实现了多光谱图像无损压缩算法.经过测试，数据吞吐率达到66 M Sample/s，比NASA JPL实验室实现的33 M Sample/s速度提升一倍.对压缩后的码流数据进行解压，解压后的数据和原图完全相同，证明了设计的正确性和合理性. Based on CCSDS 123.0-B-1 lossless multispectral image compression algorithm, a field programmable gate array (FPGA) chip of Stratix II series is implemented in all functional modules of this algorithm in RTL level. The result of the experiment showed that the data throughput of this system was 66 M sample/s, which was two times faster than 33 M sample/s achieved by NASA JPL laboratory. After decompressing the compressed bit stream, the decompressed data was identical to the original, which has proved the correctness and rationality of this implementation. CCSDS 无损压缩 多光谱图像 现场可编程门阵列(FPGA) CCSDS lossless compression multispectral image field programmable gate array (FPGA) 引言 1 随着航天研究技术的不断发展，光谱图像的分辨率不断提高，卫星上产生的海量数据对有限的信道资源造成了巨大的压力，给数据的存储和传输都带来了很大的困难，因此进行有效的数据压缩就尤为重要.而多光谱图像信息量较大，如果对多光谱图像进行有损压缩，很多细节的信息将会丢失，导致图像模糊及光谱曲线的连续性变差 [ 1 1 ] ，为了保留图像的细节信息，我们更希望大容量的多光谱图像的压缩既要高保真又要保证实时性，因此研究多光谱图像的无损压缩具有十分重要的意义. 目前，多光谱图像的无损压缩技术大致可分为以下三类：基于预测的压缩方法、基于变换的压缩方法和基于矢量量化的压缩方法 [ 2 2 ] .其中，基于预测的压缩方法提出时间早，发展比较成熟，目前在无损压缩编码领域已经得到了广泛的应用.通常来说，预测方法可以分为三种，它们分别为一维预测、二维预测和三维预测.一维预测是需要与待测的像素点在同一扫描行中前面相邻的一个或者几个像素点；二维预测除了用同一扫描行的前几个相邻的像素点以外，还需要该像素点前一行或者前几行的若干个像素点来预测；三维预测除了需要二维预测所需的像素点外，还要前一帧或者前一个波段对应位置的像素点.CCSDS (空间数据系统咨询委员会)在2012年5月发布的第一个专门针对星载多光谱/高光谱图像无损压缩的标准CCSDS 123.0-B-1 [ 3 3 ] ，该标准是基于三维自适应预测的无损压缩算法 [ 4 4 ] .该算法复杂度低，压缩率高，目前已经成为星载多光谱/高光谱图像无损压缩的国际通用标准.本文采用该标准算法，基于FPGA搭建的硬件平台，实现了多光谱无损压缩算法. 算法分析 1 CCSDS无损压缩算法与传统无损压缩算法的最大区别就是充分利用了多光谱图像谱段像素之间存在的谱间冗余，既使用谱段内，也使用谱间相邻的像素值去预测当前的像素值，通过对像素值和它的预测残差进行熵编码来实现数据的压缩.该算法由预测器和熵编码两部分组成，如 图 1 图 1 所示.预测器采用三维自适应线性预测的方法，对图像数据进行去相关处理，熵编码模块负责对去相关后的数据进行压缩编码. 预测器 2 预测就是根据之前时刻的样本值对当前的样本值进行一种模型化的运算 [ 5 5 ] .为了保证无损压缩，预测必须是可逆的. 图 2 图 2 表示三维预测模型.其中 z 表示波段， x 表示图像宽度， y 表示图像高度.预测器先计算待预测点 s z, y, x 的局部和，以此计算出局部差和方向局部差，得到当前谱段的局部差向量.再利用前 p 个谱段对应像素点的局部差向量与权重向量线性组合就得到了预测采样值和缩放采样预测值 [ 6 6 ] .每预测完成一个点，都要对当前权重向量进行更新以便计算下个点的预测值.最后通过一个映射函数把预测残差映射成无符号类型. 局部和、局部差的计算 3 局部和、局部差主要有两种计算模式：全模式邻域导向和简化模式列导向.对帧成像仪产生的并校准过的图像，采用全模式邻域导向的局部和计算方法，此方法往往产生更少的压缩图像数据量.对推扫成像仪产生的未校准图像，采用简化模式列导向的局部和计算方法将获得更高的压缩比.根据项目需求，本文采用基于全模式邻域导向的预测器对图像进行去相关处理. 邻域导向的局部和 σ z, y, x 是待预测像素点 s z, y, x 的西边( W )像素点 s z, y, x -1 、西北边( NW )像素点 s z, y -1, x -1 、北边( N )像素点 s z, y -1, x 和东北边( NE )像素点 s z, y -1, x +1 的加权和，具体见公式(1): http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458602&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458602&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458602&type=middle 全模式需要计算 p 个中央局部差和3个方向局部差， p 表示预测所用的波段数量，可以自定义，只需保证 p 小于当前待预测像素所在波段之前的波段数即可. p 个中央局部差分别是待预测像素的前 p 个波段与待预测像素同位置处像素的值与其局部和之差，例如待预测像素前一个波段与之同位置处像素 s z -1, y, x 的中央局部差 d z -1, y, x 的计算方法见公式(2)所示，其余波段的以此类推. http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458606&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458606&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458606&type=middle 3个方向局部差 d z, y, x N ， d z, y, x W ， d z, y, x NW 分别是北边、西边和西北边三个像素点与待预测像素点局部和之差，如 图 2 图 2 ，相关的计算公式见(3)、(4)和(5). http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458610&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458610&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458610&type=middle http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458612&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458612&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458612&type=middle http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458615&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458615&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458615&type=middle 由此得到全模式下的差向量 U z ( t ). http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458618&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458618&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458618&type=middle 权重向量的计算 3 全模式对应的权重向量如下式， http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458621&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458621&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458621&type=middle 对于每个波段，当前波段的三个方向局部差对应的权重初始值如下， http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458624&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458624&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458624&type=middle 各个波段中央局部差的权重初始值如下， http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458628&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458628&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458628&type=middle 其中， Ω 为用户可自定义的权重值分辨率，取值范围在[4,19], 从(8~9)式可以看出，距离待预测像素点越远的波段，计算局部差时所占的权重越小. CCSDS算法实现自适应编码的重点就是权重向量的更新.对每个待预测像素点进行预测后，权重向量的各个权重值都会更新一次，从而保证对下一个待预测像素点的预测可以更加准确.更新的依据就是待预测像素点灰度值与预测值之差，即预测误差 e z ( t ).对该值进行运算，若该值为正数，说明预测值偏小，适当增大权重值；若该值为负数，说明预测值偏大，适当减少权重值，从而使得预测值逐渐地接近待预测像素点的值.权重更新表达如(10)式所示: http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458631&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458631&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458631&type=middle 其中，裁剪函数clip的计算方法见(11)式 http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458633&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458633&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458633&type=middle 权重更新因子 ρ ( t )决定了权重值每次调整的步长，该值越小，权重值调整的越快，对图像数据统计结果的适应性越好，但压缩效果的稳定性越差. 预测值的计算 3 差向量的每个元素乘以相应的权重值并累加后得到加权和，称之为预测中央局部差 http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458650&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458650&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458650&type=middle z ，再由 http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458650&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458650&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458650&type=middle z 计算出预测值和预测误差，最后通过映射函数把预测残差映射成无符号类型，满足熵编码的输入要求. 熵编码 2 CCSDS 123.0-B-1标准算法中的熵编码指的是长度受限的Golomb编码，对符合几何对称性并且其分布特征是指数衰减概率分布的非负整数信号，Golomb编码具有较好的编码效果.此方法的使用有利于硬件的编解码.其中定义的码字长度不能超过最大码字长度 U max (比特)， U max 是整数，由用户输入.根据预测器输出的预测残差计算出 k z ( t )和 u z ( t )，其中 k z ( t )表示除数的指数， u z ( t )表示相除后的商, δ z ( t )代表预测器输出的预测残差, 具体公式如下. http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458636&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458636&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458636&type=middle 当 u z ( t ) < U max 时，码字由 u z ( t )个0、1个1和映射预测残差的低 k z ( t )位二进制数组成；否则，由 U max 个0和D位的映射预测残差组成.由于该编码算法较为常见，具体计算方法不再赘述. FPGA设计与实现 1 本文采用的平台是Stratix Ⅱ [ 7 7 ] 系列的FPGA芯片，考虑到在FPGA上实现CCSDS算法的复杂度和中间数据存储操作的便利，在本文中设定输入的多光谱图像的光谱波段数为4，单个波段的宽度是1 024个像素，高度也为1 024个像素.对于其他的不同波段数不同宽高的多光谱图像，我们只需要在上位机上更改相应的参数, 然后发送到FPGA配置相应的寄存器即可. 压缩系统的结构设计 2 CCSDS无损压缩系统的功能流程图如 图 3 图 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端进行解压验证. 数据预处理模块 2 数据预处理模块实现的功能是根据输入的原始图像数据经过处理后输出符合预测模块输入要求的数据，输出的数据包括当前像素点与其西边、西北边、北边和东北边的像素点.该模块的原理流程图，如 图 4 图 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位. 压缩算法模块 2 压缩算法模块是本文的核心模块，包括预测器和熵编码器两个模块.其具体实现流程如 图 5 图 5 所示. 预测模块 3 预测模块包括预测值计算模块和映射预测残差计算模块.考虑到预测值计算模块和映射预测残差模块存在较多的数值运算，同时权重更新模块需要用到预测值计算模块的输出作为输入.为了最大程度提高系统的压缩速度，整个模块共包含两级流水线.在预测模块加入一级流水线后，预测值计算模块和映射预测残差模块分别需要一个时钟周期来完成.插入流水线后，预测模块的最大频率较未插入流水线提升将近一倍. 预测模块主要完成的功能包括计算局部和、局部差、矢量乘法和预测样值估算等.映射预测残差模块的主要功能是把预测后的差值全部转为非负整数，以适应后面的熵编码模块的输入要求. 熵编码模块 3 熵编码模块的主要功能是实现数据的无损编码.该模块由编码采样和计数器累加器更新两个子模块组成，其中编码采样模块由乘法器、移位寄存器、加法器、除法器和求对数运算器组成.如果整个熵编码模块在一个时钟周期内完成，则这条路径的延时将会很长，很难满足实时性的要求.在这里编码采样子模块和计数器累加器子模块采用并行计算的方法，有效降低了关键路径的延时，提高了系统的频率.对于求对数则采用查表的方法计算，这是因为 k z ( t )是不大于D-2的最大非负整数，采用查表只需要使用D-2路的多路选择器，资源消耗没有明显增多，但是却能大大降低对数运算带来的路径延时和复杂度.计数器累加器更新模块由比较器、FIFO (长度为波段数目，位宽等于像素点位宽)和其他逻辑电路组成.熵编码实现的原理框图如 图 6 图 6 所示. 码流组包模块 2 码流组包的作用是将熵编码模块输出的不定长的码流数据进行整合，最终按照规定的要求进行码流的输出 [ 8 8 ] .本文组包后输出的定长数据位为8位，位序为升序.码流组包功能框图如 图 7 图 7 所示，组包的过程如 图 8 图 8 所示.先从输入的不定长码流数据中取出有效数据，然后把有效数据送到寄存器的空余部分，当寄存器被填满时把数据送到发送端. 验证与性能分析 1 本文采用Stratix Ⅱ系列的FPGA芯片作为设计验证平台，搭建了验证系统.验证采用的数据源是多幅不同宽高和波段数的多光谱图像.输入图像的宽度和高度都不超过4 095个像素，波段数不超过15个，位深度不超过14位.为了验证系统的正确性，本文测试了一系列的多光谱图像.图像来源包括用户给定的多光谱图像和Landsat卫星公开的多光谱图像. 图 9 图 9 是其中4幅测试用图. (a)山脉；(b)农田；(c)城市；(d)海边 (a) Mountains; (b) Farmland; (c) City; (d) Seaside 其中前两幅图像的尺寸为1 024×1 024，波段数为8，位深度为10位.后两幅图像的尺寸为1 024×1 024，波段数为4，位深度为8位.这四幅多光谱图像的像素存储格式均为BIP格式.解压FPGA压缩后上传的数据并和原数据比较，结果完全一致.证明本系统可以准确无误地实现对多光谱图像进行无损压缩. 表 1 表 1 是CCSDS 3种压缩算法的性能对比，CCSDS 123.0的压缩性能明显好于CCSDS 121.0 [ 9 9 ] ，略好于CCSDS 122.0 [ 10 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 结论 1 本文基于CCSDS 123.0-B-1的算法原理，设计了CCSDS多光谱无损压缩系统的整体结构.随后对预测模块、熵编码模块等进行分析并设计，最后在FPGA上实现了该算法.通过一系列的测试图验证了本系统的正确性与合理性，从而能够满足航天应用的需求. 但是本系统仍有需要完善的地方.一是系统的抗误码能力较差，当压缩后的码流中的某些位出现错误时后面的码流解压容易出错，造成恢复出来的图像存在较多错误，下一步可采用分段压缩的方法提高压缩算法模块的抗误码能力；二是本文实现的压缩算法模块部分所占用的FPGA逻辑资源较多，为了提高压缩速率，采用了面积换速度的方法，所以在资源占用上还可以进一步进行优化. 参考文献 朱学伟, 李晓辉, 朱博. 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(Ch). 万建伟, 粘永健, 苏令华, 等.高光谱图像压缩技术研究进展[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(Ch). CCSDS. CCSDS123.0-B-1. Lossless Multispectral & Hyperspectral Image Compression[S]. Washington D C: CCSDS Secretariat, 2012. 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. 邓海涛.多光谱图像自适应无损压缩算法研究与FPGA实现[D].海口:海南大学, 2014. http://cdmd.cnki.com.cn/Article/CDMD-10589-1014240043.htm DENG H T. Research and Implementation of Adaptive Lossless Compression Algorithm for Multi spectral Images Based on FPGA[D]. Haikou: Hainan University, 2014(Ch). 王丽萍.基于预测的高光谱图像无损压缩技术研究[D].西安:西安电子科技大学, 2014. http://cdmd.cnki.com.cn/Article/CDMD-10701-1014331212.htm WANG L P. Research on Prediction-based Lossless Compression for Hyperspectral Images[D]. Xi'an: Xidian University, 2014(Ch). Altera.SII5V1-4.5. Stratix Ⅱ Device Handbook[S]. California Los Angeles: Altera, 2011. 雷震霖.空间图像CCSDS压缩算法研究与FPGA实现[D].大连:大连理工大学, 2007. http://cdmd.cnki.com.cn/Article/CDMD-10141-2008034577.htm LEI Z L. CCSDS Space Image Compression Algorithm and FPGA Implementation[D]. Dalian: Dalian University of Technology, 2007(Ch). CCSDS.CCSDS 121.0-B-1. Lossless multispectral & hyperspectral image compression[S]. Washington D C: CCSDS Secretariat, 1997. CCSDS.CCSDS 122.0-B-1. Image Data Compression[S]. Washington D C: CCSDS Secretariat, 2005. 10.14188/j.1671-8836.2016.06.007 TN919.8 whdxxblxb-62-6-531 CAST创新基金项目(CAST201240) 2016-02-24 2016-12-24 2020-12-09 版权所有© 2016 《武汉大学学报(理学版)》编辑部 Copyright © 2016 Journal of Wuhan University (Natural Science Edition). All rights reserved. 2016 武汉大学学报(理学版) 6 62