含无效云区域遥感图像在轨压缩方法

，，

1. 中国民航大学 计算机科学与技术学院，天津 300300 2. 中国科学院 国家空间科学中心,北京 100190

近年来，随着国防与国民经济的发展需求，航天光学遥感技术在理论与工程技术上发展迅速，新的在轨光学遥感载荷不断涌现，空间分辨率、时间分辨率、光谱分辨率及辐射分辨率等性能指标大幅提升，使得光学遥感图像的信息量呈指数增长[1-2]，对应的实时数据传输速率已高达数吉比特每秒。然而利用目前的通信技术进行星地图像传输，传输率仅为500 Mbit/s左右[3]，星地图像传输带宽与遥感图像数据量之间存在着巨大的差距和矛盾，而且差距呈现出扩大趋势，严重制约了高分辨率光学遥感技术的应用。遥感图像信源数据量增大，而传输带宽不变，这就需要增大对信源的压缩倍数。然而对于现有的图像有损压缩算法来讲，随着图像压缩比的增大，图像的质量会急剧下降。所以直接利用现有的图像压缩方法不能很好的解决这个矛盾。地球上平均有1/3～1/2的地区覆盖着云层。大部分遥感图像中含有一定的云。对非气象类遥感卫星而言，由于云的遮挡，这部分几乎不具备任何可用信息或可用信息较少，但占用了有限的星上存储空间与星地传输带宽。在压缩过程中充分考虑并合理设计云剔除策略，减少甚至消除云无效区域的编码消耗，可以提高压缩性能，缓解数传带宽压力。

对于包含无效区域的图像的压缩问题，处理策略主要分成3大类。第1类是基于预处理的填充方法，在压缩之前按照有利于压缩的方式填充图像中的无效区域，不改变后续的压缩算法。这类填充方法复杂度低，而且输出码流遵循原压缩标准。无效区域的填充值的选取对于压缩性能有很大影响，填充值的计算主要有Phagocyte[4]和ADR[5]两种方法。Phagecyte利用指定区域窗口内像素的算数平均迭代填充无效区域，算法简单，在一定程度上提到了压缩性能，并成功应用于商业GIS产品中，但最终的无效区域并非单一的灰度值，无效区域的压缩编码仍消耗较高。由Jorge Gonalez Conejero等人提出的ADR方法利用计算出的某个单一值填充整个无效区域，降低了无效区域的压缩编码消耗，但该方法以整个有效数据区域像素的算数平均作为平均值，这个值并不能精确反映填充部分的边缘信息，所以在小波变化后，填充的边缘部分会产生较多的高频信息，压缩编码消耗较高。第2类是基于感兴趣区域(Region Of Interest,ROI)的方法，此类方法通常应用于基于小波变换的图像压缩算法中，在小波变换之后通过提升感兴趣区域的位平面达到优先编码的目的，如JPEG2000 part1中的MaxShift方法[6-8]。在输出码率不变的情况下，可以有效提高感兴趣区域的压缩性能。由于提升后的感兴趣区域(可以看作有效区域)边缘处的小波系数在恢复时丢失了一部分信息，所以会造成解压缩后有效区域边缘区域的图像质量较差。第3类是针对此类问题设计专用的压缩算法。这类算法大多是通过改进主流压缩算法流程中的某些环节，达到降低无效区域编码消耗的目的。文献[9]改进SPIHT算法，采用集合划分策略，在编码阶段直接丢弃无效区域信息。文献[10]改进EBCOT编码方法，改进后算法可以编码图像内的任意区域。这类算法一般会在编码过程中直接跳过无效区域，这样无效区域的编码消耗几乎为0，编码性能较高。但此类算法都是通用压缩方法，并未针对遥感图像压缩设计，算法设计时并未考虑空间应用特殊性及在轨成像特点，而且复杂度较高。

针对遥感图像压缩，空间数据系统咨询委员会(Consultative Committee for Space Data System, CCSDS)推出了基于小波变换的图像压缩标准CCSDS-122.0-B-1(又称CCSDS-IDC)[11-13]。标准针对空间任务图像采集特点及应用的特殊性，在现有压缩方法的基础上，综合考虑算法的压缩性能与算法复杂度，具有复杂度较低、压缩性能较高、功能全面、性能优异、实用性强等特点。所以本文在CCSDS-IDC标准基础上，针对含云图像压缩问题，提出了形状自适应CCSDS-IDC算法，改进后，含云遥感图像压缩性能得到了较大提高。

1 CCSDS-IDC算法

CCSDS-IDC压缩算法模型与传统基于DWT的压缩算法类似，如图1所示，其压缩性能与当前效果最好的JPEG2000方法相当。主要包含离散小波变换(DWT)和比特平面编码(BPE)两个部分。离散小波变换将图像变换到小波域，去除图像像素间的相关性，比特平面编码根据重要程度扫描和编码小波系数。

1.1 离散小波变换

9/7浮点DWT的两组滤波器系数为：

{h-4,h-3,h-2,h-1,h0,h1,h2,h3,h4}

{g-3,g-2,g-1,g0,g1,g2,g3}

系数如表1所示。变换公式为：

(1)

(2)

式中：j=1,2…,N,Cj和Dj分别为第j个低频和高频系数。

表1 CCSDS-IDC小波变换系数

与此对应的小波反变换公式为：

(3)

(4)

式中：q和p分别为综合滤波器低通和高通系数。

CCSDS-IDC算法的小波变换为3级，由于小波系数的对称性，二维小波变换可以按照行、列分别进行。变换后，图像的低频子带包含了图像的大部分信息。编码时为低频子带和高频子带间分配不同比例的码字，起到压缩的效果。3级二维变换示意如图2所示。

1.2 位平面编码

位平面编码将低频子带和各高频子带中对应的64个系数组成一组，即一个块。如图3所示，包含1个直流系数(DC)，63个交流系数(AC)(3个父系数、12个儿子系数和48个孙子系数)。

CCSDS-IDC标准中，BPE编码包含DC编码和AC编码两部分。DC编码采用CCSDS-RICE算法进行；AC是BPE编码的核心所在，按照比特平面的重要程度，扫描编码一个块中的63个交流小波系数。

扫描编码过程中，BPE共包含Stage0～Stage4的5个阶段。Stage0对直流系数进行RICE编码，Stage1按比特平面扫描编码3个父系数，Stage2按比特平面扫描编码12个子系数，Stage3按比特平面扫描编码48个孙子系数，Stage4扫描量化编码比特平面中的非重要系数。编码顺序示意如图4所示。

2 形状自适应CCSDS-IDC算法

CCSDS-IDC算法未给出针对无效区域的编码方法。本文为了提高包含无效区域(云区域)遥感图像的压缩性能，对CCSDS-IDC的DWT和BPE进行改进，提出形状自适应CCSDS-IDC(SA-CCSDS-IDC)方法来尽量减少无效区域压缩成本。改进后得到的SA-CCSDS-IDC算法流程示意如图5所示。首先根据掩膜图像信息对原始图像中有效区域像素进行3级形状自适应小波变换(Shape Adaptive Discrete Wavelet Transform， SA-DWT)，为了标识小波变换后有效区域内各个小波系数的位置信息，对掩膜图像进行Lazy小波变换，并用变换后的结果指导形状自适应位平面编码(Shape Adaptive Bit Plane Encoding，SA-BPE)，在扫描和编码过程直接跳过无效区域。为了解码时正确恢复各个像素之间的位置关系，对二值的掩膜图像利用链码进行形状编码，并将编码后的结果与SA-BPE有效区域编码信息融合到一起，最终形成压缩码流并输出。图像解压缩过程中，首先通过解码形状信息恢复云掩膜，然后根据掩膜解码比特平面信息，将结果进行形状自适应小波反变换，最终恢复出原始图像的有效区域部分。

2.1 形状自适应小波变换

经过云剔除后的遥感图像中的有效区域不再是规则的矩形，为提高压缩性能，利用形状自适应小波变换，在小波变换阶段直接跳过无效的云区域像素点。SA-DWT的基本思想是：进行一维小波变换时，利用云掩膜信息，对图像一行中连续的任意长度像素段进行小波变换，变换时，主要考虑以下两个问题：1)经过SA-DWT后的小波系数个数不能增加，应该与原有有效区域的像素个数相同。2)SA-DWT应该保留原始DWT中影响编码性能的一些特性，如空间相关性，变换的局部特性及子带间的自相似性等。相比传统DWT，SA-DWT主要包含两个特性：边缘扩展和采样策略。

(1)边缘扩展

(9,7)小波是一种奇对称双正交小波，在进行变换时，为保证无损重建，图像的边缘处可采用镜像扩展方式，如图6所示。当连续的像素段比较短时，扩展时开始像素与结尾像素会相互影响，需要从开始像素到结束像素递归交替地扩展边缘处的像素点。所以边界扩展更为复杂，需要特殊处理。图7给了当连续像素段的长度为2时的扩展实例。

(2)下采样策略

对任意长度的像素段进行小波变换后，产生低通和高通两个子带系数。不同的采样系数方式，会产生不同的效果。SA-DWT的采样方式主要有两种：局部下采样和全局下采样，如图8所示。对于局部下采样来讲，同一行中的每一个单独的连续像素段，从起始像素开始进行偶采样或奇采样，而对于全局下采样, 总是从当前行的由于连续像素段的长起始像素开始进行采样。度是不固定的，有可能是奇数，也有可能是偶数。当长度为奇数时，如果进行全局下采样, 有可能造成高通系数比低通系数多，这样更多小波系数信息位于高频系数中，这非常不利于压缩编码，但这种方式可以保持小波系数间的相位关系不变。根据文献[14]，利用局部下采样策略可以得到更优的编码压缩性能。故本文为提高剔除云后的有效区域的压缩性能，采用局部下采样策略。

2.2 形状位平面编码

除了对DWT进行形状自适应的改进外，BPE阶段同样需要改进。在编码过程中，改进后的SA-BPE直接跳过无效区域的小波系数，可以大大提高编码的效率。CCSDS-IDC编码过程中，将DWT后的小波系数按照子带分成多个块，并在各个块之间进行独立扫描和处理。改进后的BPE，跳过无效区域，直接编码有效区域，如图9所示。

原CCSDS-IDC算法对位平面扫描后得到扫描字和转义字，这些扫描结果和小波系数位置关系有关，而且长度是固定的(3 bit或者4 bit)。但是经过SA-DWT后，高频系数个数不再和低频系数个数相同，这样每个家族块中的小波系数有可能缺失。本文SA-BPE在块内扫描时，对缺失的小波填充0来处理。在对这些扫描结果进行熵编码时，由于熵编码本身算法复杂度较高，CCSDS为折中算法复杂度和算法压缩性能，统计了大量遥感图像的位平面扫描数据，依据统计编码结果给出了固定的编码表，及利用查找编码表代替熵编码计算进行扫描结果的编码，大大降低了算法复杂度。但本文SA-BPE扫描过程中，人为地在缺失像素的位置添加了0，使扫描结果的统计特性发成了变化，这样原CCSDS-IDC的编码表不在是统计上最优的。所以本文对123幅含云遥感图像的SA-BPE的扫描结果进行统计，利用Huffman编码重新设计了转义字及扫描字的编码查找表。

2.3 形状编码

对于SA-CCSDS-IDC算法，经过云检测后的地物区域的形状往往是不规则的，需要在编码过程中编码地物区域边缘所围成的区域形状信息。轮廓像素编码后，就可以在解码时定位该区域，精确恢复地物的形状信息。

链码是一种简单有效的轮廓信息编码方法[15]。链码的编码流程如下：以轮廓上的任意一点(x,y)为起点开始扫描编码，将相邻轮廓点的相对位置关系用0～7的8个数字，每种位置关系用3比特表示，分别是000～111，如图10所示。然后顺序扫描编码其他轮廓点。这种编码轮廓的方法可以精确的表示轮廓，并起到压缩的作用。若需要进一步提高压缩性能，可以首先将轮廓进行修正，增加拐点位置的轮廓像素，使轮廓的8邻域转化成4邻域，这样位置关系就可以用2比特表示。这在一定程度上造成了原轮廓的失真。本文结合以上两种方式，采用一个改进的链码来编码轮廓信息。相邻的3个轮廓点之间有一定的约束关系，已知前2个轮廓点，第3个轮廓点的可能邻域信息只有4个不同的链码。这样就可以在保证无失真的编码轮廓的同时，只采用2比特来编码图像边缘轮廓。

3 试验结果及分析

3.1 测试图像集

目前高分辨率光学遥感卫星的主要成像载荷是推扫式被动光学传感器，探测器阵列横向排列于光学焦平面上即垂直于卫星运动方向，在卫星完成向前飞行的过程中纵向对地面扫描，从而得到地面目标的一个二维条带状图像，在轨压缩处理也是按照条带完成的。以当前在轨商业光学遥感卫星QuickBird及OrbView遥感全色图像为基础，通过插值、采样等变换后，得到条带状的模拟测试图像。图像测试集包括不同场景、不同云量的共计123张测试图。图11列出了部分典型测试图像的原图及其云区域掩膜图像，包括Cloud_1、Cloud_2、Cloud_3、Cloud_4及Cloud_5，其中云区域掩膜中黑色代表地物区域，白色代表云无效区域。表2给出了部分测试集中图像的一些参数，包括图像分辨率、无效区域占比、灰度均值、灰度方差。灰度方差在一定程度上可以反应图像的复杂程度，灰度方差越大，图像场景越复杂。

3.2 结果及性能分析

利用以上测试图像集，测试和比较了本文给出的针对含云遥感图像压缩的SA-CCSDS-IDC方法相比未经云剔除的原始CCSDS-IDC算法、JPEG2000算法和经过ADR填充后的CCSDS-IDC算法的压缩性能，表3～表7给出了Cloud_1 ～ Cloud_5典型测试图像的压缩性能测试结果。

表2 典型测试图像参数

表3 典型测试图像参数Cloud_1不同方案不同压缩倍数下的压缩性能PSNR

通过对比分析以上结果，可以得出以下结论：

1)在所有压缩倍数下，本文SA-CCSDS-IDC方法相比未经云剔除处理的CCSDS-IDC、JPEG2000方法及经过ADR填充后的CCSDS-IDC方法解压后的图像质量均有提高。对于Cloud_1图像，在压缩倍数为4、8、16、32、64、128、256下，相比原始CCSDS-IDC算法，解压缩后图像PSNR分别提高了5.87 dB、5.88 dB、4.16 dB、3.50 dB、2.73 dB、2.75 dB、2.50 dB，图像压缩性能得到了很大的提升。

2)在压缩质量相当的情况下，SA-CCSDS-IDC方法可以在一定程度上提高图像的压缩倍数。对于Cloud_2图像，JPEG2000在压缩32倍时的PSNR为38.81 dB，而SA-CCSDS-IDC方法在压缩64倍时的PSNR为39.40 dB，这说明图像质量相当的情况下，压缩倍数可以极大地提高，图像压缩后的数据量为之前的1/2。经过整个图像测试集进行测算，在图像质量固定的条件下，对于地物探测场景，压缩比平均提高75%，数传数据量减少42%(云占比平均43%)，可以有效减轻数传带宽压力。卫星每轨下行数传总量固定，随着数据量的减小，单轨成像时间延长75%，卫星将可获得更多有价值数据，有效提升遥感卫星工作效能。

表4 典型测试图像参数Cloud_2不同方案不同压缩倍数下的压缩性能PSNR

表5 典型测试图像参数Cloud_3不同方案不同压缩倍数下的压缩性能PSNR

3)随着压缩倍数的提高，本文SA-CCSDS-IDC算法图像压缩性能的提高程度呈下降趋势。

表6 典型测试图像参数Cloud_4不同方案不同压缩倍数下的压缩性能PSNR

表7 典型测试图像参数Cloud_5不同方案不同压缩倍数下的压缩性能PSNR

图12给出了不同测试图像组(不同颜色曲线)在不同码率下的本文方法的PSNR提升情况，可以看出，在较低码速率下，本文算法在压缩性能上并不能获得较大的提高。另外，Cloud_1到Cloud_5这5幅测试图像的云区域占比是逐渐下降的。通过下图可以看出，云区域占比的越大，在各个不同压缩倍数情况下，本文SA-CCSDS-IDC算法的压缩性能提升程度越大，即压缩性能提升程度与云区域占比正相关。另外，图像场景越复杂，本文改进方法性能提升程度越高，由于Cloud_4图像灰度方差较大，场景复杂，所以经过云剔除处理后压缩性能提高更多，与Cloud_3接近。

为进一步说明本文两种改进算法对于压缩后图像质量的提升情况，图13给出了Cloud_2图像压缩倍数为4和64情况下，CCSDS-IDC算法、经过ADR填充的CCSDS-IDC算法与本文SA-CCSDS-IDC算法解压后图像细节展示。可以看出，在压缩4倍的情况下，由于3种方法的PSNR均比较高，主观视觉上并没有明显的区别；在压缩64倍的情况，CCSDS-IDC的细节丢失明显，ADR+CCSDS-IDC次之，本文方法图像质量最高，保留了原图像中的大部分细节。

4 结束语

遥感图像中含有大量的云无效区域，占用了大量的数传带宽。在压缩过程中剔除云无效区域，可以进一步提高压缩性能。本文在CCSDS-IDC算法基础上，针对云剔除后的遥感图像形状不固定的特点，改进了小波变换和位平面编码，使其能够根据云掩膜图像自适应地压缩编码图像中的有效区域。利用图像测试集，对方法进行了比较分析，试验结果表明，本文SA-CCSDS-IDC方法压缩性能明显优于原CCSDS-IDC方法，适用于含云遥感图像在轨压缩应用。下一步拟对算法进行硬件FPGA实现，并在实际空间任务中得到应用。