一种多通道卫星图像无损压缩方法

摘要

本发明公开了一种多通道卫星图像无损压缩方法，属于卫星遥感数据处理和图像压缩技术领域。本发明针对多通道卫星图像数据，在现有无损压缩编码技术的基础上，充分利用不同通道之间卫星云图的相似性，先求图像的差，再进行无损压缩编码，减少了不同码字的个数，从而提高压缩比。本发明对多通道卫星图像的平均压缩比可以达到2以上，比现有方法节省一半以上的存储空间。

说明书

技术领域

本发明涉及一种多通道卫星图像无损压缩方法，包括编码和解码两部分，属于卫星遥感数据处理和图像压缩技术领域。

背景技术

随着气象卫星技术的发展，我们每天可以获得巨量的卫星遥感资料，而且这些资料覆盖面广、空间分辨率高、时间连续性强，不受时间、天气状况等限制，在气象观测和天气预报中所起的作用也越来越重要。

气象卫星一般带有多个光谱通道的辐射扫描仪，同时对地球表面不同光谱通道的辐射进行扫描，并生成多个记录文件。比如，我国的风云二号E星（FY-2E）搭载有5个光谱通道进行辐射扫描，分别称为4个红外通道（IR1、IR2、IR3、IR4）和1个可见光通道（VIS）。

FY-2E 卫星4个红外通道的数据量化等级为1024，即每个记录值的范围为0～1023，至少用10个二进制位表示；可见光通道的数据量化等级为64，相应的记录值范围为0～63，至少用6个二进制位表示。气象业务部门对于FY-2E卫星图像资料的存储一直采用比较简单的非压缩方式，比如GPF格式卫星云图数据，一个GPF格式的文件动辄几兆甚至几十兆，这给资料存储和传输带来很多困难。GPF格式文件保存卫星图像数据时规定红外通道的每个记录值用一个16位二进制数表示（2个字节），可见光通道的每个记录值用1个8为二进制数表示（1个字节）。因此，现有的卫星图像数据文件相对于卫星观测的实际量化等级来说是存在压缩的可能性。

对卫星图像进行压缩编码的目的是为了更方便的进行存储和传输，以便于更好地利用这些卫星遥感获得的珍贵资料，因此必须对卫星图像进行无损压缩。

通过对卫星云图的对比观察发现，卫星云图灰度分布极其不均匀，往往集中在某一部分灰度范围内，而且不同通道之间的云图之间存在较强的相关性，尤其是光谱范围相邻的的通道，比如IR1通道和IR2通道的卫星云图相关性达到99%以上。因此，也存在利用信息分布不均匀性和通道相关性对多通道卫星图像进行压缩的可能性。

发明内容

本发明所要解决的技术问题在于克服现有技术不足，提供一种多通道卫星图像无损压缩方法，利用信息分布不均匀性和通道之间的相关性对多通道卫星图像进行无损压缩，从而获得更高的压缩比。

本发明的多通道卫星图像无损压缩方法，包括卫星图像的编码和解码，所述编码包括以下步骤：

步骤A1、对卫星各通道的原始图像C_i分别进行无损压缩编码，得到编码序列记为H_i0，编码参数记为T_i0，平均码长记为L_i0，i=1…n，n为卫星的图像通道数；

步骤A2、对任意两个通道的原始图像做差运算，得到的差图像记为D_ij，i,j为两幅做差运算图像的通道编号；

步骤A3、对得到的各差图像分别进行无损压缩编码，得到的编码序列记为H_ij，编码参数记为T_ij，平均码长记为L_ij，i,j=1…n且i≠j；

步骤A4、求所有已编码图像中最小的平均码长，记为L_t0=min{L_i0| i=1…n}；

步骤A5、将H_t0，T_t0，放入存储队列S，即S={H_t0，T_t0}；

步骤A6、令K={1,2,…,n}-{t}，M={t}；

步骤A7、求剩余已编码图像中最小的平均码长，Lpq=min{ L_i0, L_ij|i∈K,j∈M }；

步骤A8、若 q=0，则将第p个通道的图像C_p的编码序列H_p0、编码参数T_p0放入存储队列S，即S=S∪{H_t0，T_t0}；然后转步骤A10；

步骤A9、若q≠0，则将第p,q两个通道的差图像的编码序列H_pq、编码参数T_pq放入存储队列S，即S=S∪{H_pq，T_pq}；然后转步骤A10；

步骤A10、K= K - {p}，M= M ∪ {p}

步骤A11、若K不为空集，转步骤A7；否则转步骤A12；

步骤A12、将当前存储队列S中的数据作为压缩后的数据进行存储或传输；

所述解码包括以下步骤：

步骤B1、从压缩后的数据中读取出所有的编码序列和编码参数；

步骤B2、找出所有由原始图像直接无损压缩编码得到的编码序列H_i0，根据所对应的编码参数T_i0对其直接进行无损压缩解码，得到原始图像C_i；

步骤B3、令M={i|编码序列H_i0存在}为已经得到原始图像的通道组成的集合，K={1,2,…,n}-M为尚未解码的通道的集合；

步骤B4、找第一个H_pq，使得H_pq存在且p∈K,q∈M；

步骤B5、根据所对应的编码参数T_pq对H_pq进行无损压缩解码，得到D_pq；

步骤B6、将已解码的原始图像C_q与D_pq作和运算，得到第p个通道的原始图像C_p；

步骤B7、K= K - {p}，M= M ∪ {p}

步骤B8、若K不为空集，转步骤B4；否则转步骤B9；

步骤B9、得到所有n个通道的原始图像。

上述技术方案中，所述无损压缩编码/解码可采用各种现有的无损压缩编码方法，例如香农-范诺编码、算数编码、RLE编码等，本发明优选哈夫曼编码。哈夫曼编码是Huffman于1952年提出一种编码方法，该方法完全依据字符出现概率来构造异字头的平均长度最短的二进制码字，对信号进行无损压缩，且压缩率接近信号的熵，因此有时称之为最佳编码。

相比现有技术，本发明具有以下有益效果：

本发明方法充分利用了不同通道之间卫星云图的相似性，先求图像的差，再进行压缩编码，减少了不同码字的个数，从而提高压缩比，本发明的平均压缩比可以达到2以上，比现有方法可以节省一半以上的存储空间；本发明利用哈夫曼编码方法作为主要的编码方法，既保证了编码的高压缩比，同时也保证了数据的无损压缩，有利于以后完整地利用卫星云图提供的信息。

附图说明

图1为本发明方法的编码流程示意图；

图2为本发明方法的解码流程示意图；

图3为哈夫曼编码的原理示意图

图4为采用本发明方法对FY-2C卫星图像进行编码的过程示意图；

图5为采用本发明方法对FY-2C卫星图像进行解码的过程示意图。

具体实施方式

下面结合附图对本发明的技术方案进行详细说明：

本发明的思路是针对多通道卫星图像数据，在现有无损压缩编码技术的基础上，充分利用不同通道之间卫星云图的相似性，先求图像的差，再进行无损压缩编码，减少了不同码字的个数，从而提高压缩比。

下面以哈夫曼编码为例，来对本发明技术方案进行详细说明。

本发明分成编码和解码两个部分。对具有n个通道的卫星图像进行编码的时候，首先根据哈夫曼编码的平均码长决定是直接对卫星云图本身进行哈夫曼编码，还是将该通道的图像与其他通道做差运算，然后再对差图像进行哈夫曼编码，其编码过程如图1所示，包括以下步骤：

(1.1)        对卫星每个通道的图像，记为C_i，进行哈夫曼编码，得到哈夫曼编码序列记为H_i0，哈夫曼码表记为T_i0，平均码长记为L_i0，i=1…n；

(1.2)        对任意两个通道的卫星图像做差运算，得到的差图像记为D_ij，D_ij(x,y)=C_i(x,y) - C_j(x,y)，i,j为通道编号，x,y为图像坐标；

(1.3)        对差图像进行哈夫曼编码，得到的哈夫曼编码序列记为H_ij，哈夫曼码表记为T_ij，平均码长记为L_ij，i,j=1…n且i≠j；

(1.4)        求所有哈夫曼编码中最小的平均码长，记为L_t0=min{L_i0| i=1…n}；

(1.5)        将H_t0，T_t0，放入存储队列S，即S={H_t0，T_t0}

(1.6)        令K={1,2,…,n}-{t}，M={t}；

(1.7)        求剩余通道图像哈夫曼编码以及差图像的哈夫曼编码中平均码长最小的结果，Lpq=min{ L_i0, L_ij|i∈K,j∈M }；

(1.8)        若 q=0，说明对第p个通道的图像直接进行哈夫曼编码的平均码长比较短，应直接存储该图像的哈夫曼编码，则将H_p0，T_p0，放入存储队列S，即S=S∪{H_t0，T_t0}；转(1.10)步；

(1.9)        若q≠0，说明第p,q两个通道的差图像的平均码长比较短，存储该差图像的哈夫曼编码，则将H_pq，T_pq，放入存储队列S，即S=S∪{H_pq，T_pq}；转(1.10)步；

(1.10)    K= K - {p}，M= M ∪ {p}

(1.11)    若K不为空集，转(1.7)步；否则转(1.12)步；

（1.12）将S中所有的哈夫曼编码H_ij及码表T_ij存储到云图数据文件中。

对压缩云图的解压缩过程首先要根据存储的数据判断该哈夫曼编码序列是由原始图像直接编码得来的，还是由差图像编码得来的。若是原始图像编码得来的，则直接解码得到原始图像，若是由差图像编码来的，解码后还需与相应的图像做和运算才能得到原始的图像。其解码过程如图2所示，具体包括以下步骤：

(2.1)        从云图数据文件中读取所有的H_ij，T_ij；

(2.2)        找出所有由云图原始图像直接编码得到的编码序列H_i0，直接解码得到C_i；

(2.3)        令M={i|编码序列H_i0存在}为已经得到原始图像的通道组成的集合，K={1,2,…,n}-M为尚未解码的通道的集合；

(2.4)        找第一个H_pq，使得H_pq存在且p∈K,q∈M；

(2.5)        根据T_pq对H_pq进行哈夫曼解码，得到D_pq；

(2.6)        计算原始图像C_p，C_p(x,y) = D_pq(x,y) + C_q(x,y)；

(2.7)        K= K - {p}，M= M ∪ {p}

(2.8)        若K不为空集，转(2.4)步；否则转(2.9)步；

（2.9）  所有通道的图像均已得到，算法结束。

其中，对两个图像进行差/和运算，就是对两个图像中对应像素点的值做差/和运算，可表示为：D_ij(x,y)=C_i(x,y) - C_j(x,y)，D_ij(x,y)=C_i(x,y) +C_j(x,y)，i,j为通道编号，（x,y）为图像中像素点的坐标。

哈夫曼编码算法为现有技术，主要包括以下步骤：1)统计图像中每个灰度值出现的概率，按灰度出现的概率从大到小排列；2)把最小的两个概率相加合并成新的概率，与剩余的概率组成新的概率集合；3)对新的概率集合重新排序，再次把其中最小的两个概率相加，组成新的概率集合，如此重复进行，直到最后两个概率的和为1；4）分配码字：码字的分配从最后一步开始反向进行，对于每次相加的两个概率一个赋予“0”一个赋予“1”，读出时由每个灰度值的概率值开始一直走到最后的概率和“1”，将路线上遇到的“0”和“1”按低位到高位的顺序排列，即可得到每个灰度值的哈夫曼编码。例如对一个包含6个不同信号的信源进行哈夫曼编码，每个信号的出现概率和最后编码结果如表1所示，编码过程如图3所示。

表1

字母abcdef频率0.10.50.120.180.070.03哈夫曼编码000101001100100011

将图像中出现的每一个灰度值的哈夫曼编码按照编码所对应的二进制数的大小排列成该图像哈夫曼编码的码表T，码表的大小与灰度值的个数和每个灰度值的码长有关。例如：FY-2E卫星的红外通道的数据量化等级为1024，即最多有1024个灰度级，而根据实验统计，每个灰度值的码长最多不超过24位=3字节，因此每个红外通道的码表长度最多为24×1024位=3K字节；可见光通的数据量化等级为64，每个灰度的码长不超过20位<3字节，因此可见光通道的码表最多为3×64=192字节。

将每个通道的图像中的灰度值按照相应的码表转化成二进制序列，即可得到该图像的哈夫曼编码序列H，哈夫曼序列的长度取决该图像哈夫曼编码的平均码长以及图像中像素点的个数，Hsize=L×n×m，其中，L为平均码长，n、m分别为图像的长度和宽度。

哈夫曼解码算法具体为：根据每个通道图像的码表T和编码序列H字节进行解码，因为哈夫曼编码为异字头码，因此可以快速、唯一地译出原始数据。

为了验证本发明方法，对FY-2C卫星的GPF格式的卫星云图文件用本发明压缩方法进行压缩编码，并给出编码结果。

经过对FY-2C卫星云图每个通道图像的哈夫曼编码码长分析，我们对红外3通道IR3和可见光通道VIS数据直接进行哈夫曼编码，得到的编码序列记为H₃₀、H_VIS，相应的码表分别记为T₃₀、T_VIS；将IR1和IR2、IR2和IR4、IR3和IR4分别做差运算，得到的差图像记为D₁₂、D₂₄、D₃₄；再对D₁₂、D₂₄、D₃₄进行哈夫曼编码，得到编码序列H₁₂、H₂₄、H₃₄，和码表T₁₂、T₂₄、T₃₄；最后对所有的编码序列和码表进行存储。采用本发明方法对FY-2C卫星图像进行编码、解码的过程分别如图4、图5所示。

第一个文件是我国中部地区的GPF格式卫星云图，经度范围从100oE到115oE，纬度范围从20oN到35oN，云图的像素分辨率为333×336，共有111888个像素点。按照原有GPF格式文件的存储方式，红外通道每个像素点用16位=2个字节表示，可见光通道每个像素用8位=1个字节表示，因此该文件的数据部分总共需要字节数为：B1=111888×(2*4+1)= 1006992字节≈0.96MB。

采用本发明方法进行压缩编码：首先对IR3和VIS通道的图像直接进行哈夫曼编码，其平均码长分别为：L₃₀=7.4477位、L_VIS=4.4296位，最大码长分别是Lmax₃=16位、Lmax_VIS=16位。然后求出通道之间的差图像D₁₂，D₂₄，D₃₄并进行哈夫曼编码，其平均码长分别为：L₁₂=5.7974位、L₂₄=8.0546位、L34=7.2105位，最大码长分别是Lmax₁₂=17位、Lmax₂₄=17、Lmax₃₄=17位。因此，保存哈夫曼编码序列时需要的字节数为：

B2=n×m×(L12+L24+L34+L3+LVIS)/8                 

=333×336×(5.7974+8.0546+7.2105+7.4477+4.4296)/8

=460696.0428字节≈0.4394MB。

另外，还需要保存每个哈夫曼编码的码表。每个编码的最长码字都小于24，即可以用3个字节来表示一个码字。 FY-2E卫星的红外通道的数据量化等级为1024，码字个数也最多为1024；可见光通道数据量化等级为64，因此码字最多为64；而三个差图像由于云图通道之间存在着明显的相关性，比如IR1和IR2的相似度达到0.99，因此差图像中不同差值的个数也远远小于1024。所以保存所有码表的字节数为：B3=3×(1024×4+64)= 12480字节≈0.0119MB。

因此，用本发明方法对原有GPF文件的数据进行压缩后，其占用的字节数为B2+B3=0.4394+0.0119= 0.4513MB，仅为原GPF文件数据量的一半不到。其压缩比为：C=B1/(B2+B3)=0.96/0.4513=2.1272。

第二个文件是覆盖我国所有地区的GPF格式卫星云图，经度范围从70oE到140oE，纬度范围从15oN到55oN，云图的像素分辨率为889×1560，共有1386840个像素点。按照原有GPF格式文件的存储方式，需要字节数为：B1=1386840×(2*4+1)= 12481560字节≈11.9MB。

采用本发明方法进行编码的话， IR3和VIS通道的图像平均码长分别为：L3=8.3316位、LVIS=4.5873位，最大码长分别是Lmax3=21位、LmaxVIS=16位。差图像D₁₂，D₂₄，D₃₄的平均码长分别为：L₁₂=5.0317位、L₂₄=8.1044位、L₃₄=7.9188位，最大码长分别是Lmax₁₂=21位、Lmax₂₄=20、Lmax₃₄=21位。因此，保存哈夫曼编码序列时需要的字节数为：

B2=n×m×(L12+L24+L34+L3+LVIS)/8                 

=889×1560×(5.0317+8.1044+7.9188+8.3316 4.5873)/8 

=5889528.1字节≈5.6167MB。

保存所有码表的字节数仍旧为：B3= 12480字节≈0.0119MB。

因此，第二个GPF文件字节数仅为B2+B3=5.6167+0.0119= 5.6286MB，其压缩比为：C=B1/(B2+B3)= 5.6286/11.9=2.1142。

下表2显示了两个文件的压缩编码结果：

表2

文件图像分辨率原GPF文件数据量哈夫曼编码序列数据量码表数据量本方法数据量压缩比文件1333×3360.96MB0.4394MB0.0119MB0.4513MB2.1272文件2889×158011.9MB5.6167MB0.0119MB5.6286MB2.1142

从表2的结果可以看出，两个文件的压缩比都达到了2以上，证明利用本发明方法可以将卫星云图文件大小压缩到原文件的一半不到，从而可以大大节省卫星资料的存储空间和传输所使用的信道资源。