1、基于小波理论的遥感图像高保真压缩方法研究1 引 言随着航天遥感技术的迅速发展，遥感数据量日益庞大，有限的信道容量与传输大量遥感数据的需求之间的矛盾日益突出。数据压缩技术作为解决这一问题的有效途径，在遥感领域越来越受到重视。一般地，图像压缩技术可分为两大类：无失真压缩技术和限失真(率失真)压缩技术。遥感应用中还经常提到一种分类压缩方法，是指当用户只对图像中具有某种特征的部分感兴趣时，可在遥感平台上抽取特征并且只传输符合特征的图像子集。这实际上也应属于限失真压缩1 。在遥感信息处理领域，根据信息处理的阶段性，遥感图像压缩又可分为星上无损压缩、星上有损压缩和地面遥感数据压缩。为了最大限度地保持遥感传
2、感器所获取的目标信息，星上压缩一般采用无损压缩方法。但当信息量大到无损压缩难以满足要求时，也考虑失真量较小的有损压缩，即高保真压缩方法。同时，为了适应遥感数据采样率较高的特点，星上压缩的实时性要求较高，因而要求压缩方法计算简单，硬件复杂度低。本文在分析遥感图像及其小波分解系数特性的基础上，提出了基于小波分析理论的自适应标量、矢量混合量化压缩方法。该方法先对遥感图像进行小波分解，然后以纹理复杂程度作为区域重要性度量，通过对纹理复杂的重要区域进行标量编码来保证恢复图像的质量，通过对平坦区(即不重要区)进行矢量编码来提高压缩比。实验结果表明该方法具有压缩率较高，图像恢复质量好，速度快等优点，十分适合
3、遥感数据的高保真压缩。下面我们先对小波理论作简要介绍，然后着重讨论自适应标量、矢量混合量化压缩方法。2 小波分析基本理论及其在图像压缩中的应用与傅里叶变换相似，小波变换是一种同时具有时频二维分辨率的变换。其优于傅氏变换之处在于它具有时域和频域“变焦距”特性，十分有利于信号的精细分析。第一个正交小波基是 Harr 于1910 年构造的；但 Harr 小波基是不连续的。到 80 年代，Meyer, Daubechies 等人从尺度函数的角度出发构造出了连续正交小波基。1989 年，Mallat 等人在前人大量工作的基础上提出多尺度分析的概念和基于多尺度分析的小波基构造方法，将小波正交基的构造纳入统
4、一的框架之中，使小波分析成为一种实用的信号分析工具2 。以一维情况为例，信号的多尺度分析算法可表示如下：其中 h(n)和 g(n)分别是低通和高通滤波器，m 是多尺度分解的层数，Sm(n)是 Sm-1(n)在 2-m 上的近似，Cm(n)是信号 S 从 2-m+1到 2-m 的尺度近似时丢失的信息。上述分解算法有精确的重构算法：上述算法推广到二维情形，可对图像的每一行、然后对每一列分别进行小波分析，最后得到代表了原图不同频率特性和方向特性的 4 个子图。进行三级小波多尺度分解后的结果如图 1，其中 LL3为低频分量；HLi 为竖直边缘细节；LHi 为水平边缘细节；HHi 对应45，135方向上
5、的细节(i=1,2,3,)。关于小波变换和多尺度分解的更详细讨论可参考有关文献2，3 。小波变换在压缩中提供了如下优点：(1) 多尺度分解提供了不同尺度下图像的信息，并且变换后的能量大部分集中在低频部分，方便了我们对不同尺度下的小波系数分别设计量化编码方案，在提高图像压缩比的情况下保持好的视觉效果和较高的 PSNR。(2) 小波分解和重构算法是循环使用的，易于硬件实现。Mallat 阐述了小波分解的时间复杂度为 N，还给出了快速小波算法2 ，其时间复杂度为 logN。图 1 3 层塔式小波分解示意图Fig. 1 Three-layer pyrimid wavelet decomposition
6、3 遥感图像及其小波分解系数的特性分析与一般自然图像不同，遥感图像空间局部相关性较弱，信息熵较高，所以对其进行大压缩比的无损压缩是很困难的。而遥感图像的小波变换系数却具有很强的局部相关性。图 2(a)所示的遥感图像，采用 Antonini 的 9-7 双正交小波(biorthogonal wavelet)，其 3 层小波分解系数统计如表 1(我们将小波分解滤波器的因子放到了重构滤波器中，这不影响数据重构)。经实验分析，其所有细节子图近似遵循零均值 Gamma 分布：图 2 用于实验的 4 幅遥感图像Fig. 2 The four test image我们选择研究了图 2(a)的小波分解高频细节
7、子图的所有 44局部区域。选择块的大小为 44 是基于如下理由： (1) 基于 3 层小波分解的图像要求宽和高为 8 的倍数，因而块的宽和高最好也为2 的幂，如 2、4、8，使对图像分块时，块内数据都是有用数据，不必因为块的尺寸不合适而补 0；(2) 块的尺寸太小，如取 2 时，对下述的不重要区压缩比不高，而尺寸太大，如取 8 时，则会使非纹理重要区大幅度减少，难以实现对整幅图像的高效压缩。实验表明，选择 44 是合适的。对这些块进行统计分析还表明：除极少数纹理较强的局部区域之外，绝大部分纹理不太强的区域近似服从广义Gaussian 分布，特别是平坦区域，均值大多约等于 0。例如，图2(a)经
8、 3 层小波分解后除 LL3 外的细节子图中的 44 局部区域共16128 块，其中近 97%的非纹理区呈近似的 Gaussian 分布，全部平坦区的均值在0.2 之间。这从实验上为我们采用下节所述的自适应标量、矢量混合量化压缩方法提供了有力的理论依据。需要说明的是，尽管只分析了一幅遥感图像，但我们实现的遥感图像高保真压缩系统良好的性能反过来说明了我们的上述看法是基本正确的。遥感图像的另一个特点是纹理丰富复杂。JPEG 的 88 分块压缩方法压缩纹理复杂的块时恢复误差较大，具有比较明显的方块效应，而基于小波变换的图像压缩方法较好地克服了方块效应的影响。通过对不同区域采用不同编码方法，可以较好地
9、保持原图的纹理信息，并达到较高的压缩比(表 2)。表 1 图 2(a)的小波分解细节子图的统计数据 Table 1 The statistic data of wavelet high frequency of fig. 2(a)HL3 LH3 HH3 HL2 LH2 HH2 HL1 LH1 HH1 最大值 53.99 61.35 36.56 47.11 62.13 18.38 103.14 40.67 9.66 最小值 -40.76 -38.42 -23.06 -35.84 -52.44 -20.89 -31.28 -56.22 -27.24 平均值 0.06 0.04 0.05 0.2 -
10、0.03 0.01 0.35 -0.06 -0.14 Var 3.94 4.32 2.6 3.63 2.84 1.56 6.08 2.01 0.88 表 2 3 种压缩方法的性能比较 Table 2 The comparison of performances of three compression method压缩方法 图 2(a) 图 2(b) 图 2(c) 图 2(d) 压缩比 PSNR 压缩比 PSNR 压缩比 PSNR 压缩比 PSNR WINZIP 1.63 1.38 1.74 1.7 JPEG 9.06 40.18 7.02 36.67 12.79 39.21 7.58 40.
11、53 本方法 16.85 40.23 11.97 40.12 17.74 39.27 15.21 40.73 4 自适应标量、矢量混合量化编码方案 基于小波分解的图像压缩方法的一个重要因素是量化方案的选择。一般说来，量化方法分为标量量化和矢量量化两种。标量量化的关键是去相关和编码。目前主要的去相关技术是预测方法，如 DPCM 预测；而编码仍以熵编码为主。标量量化的特点是可保持较高的图像恢复质量，但压缩率一般较低。矢量量化有效地利用了矢量各分量间的相互关联的性质(线性依赖性、非线性依赖性、概率函数的形状以及矢量维数)来进行去相关处理，具有较高的压缩比4 。但矢量量化码书的训练和图像压缩时间长，图
12、像恢复质量相对较差，在时间或质量要求较高的图像压缩系统中不常采用。目前最有效的基于小波分解的矢量量化方法有法国 M. Barlaud 等人提出的 PLVQ 塔式格型矢量量化方法5和美国 J. M. Shapiro 提出的 EZW 方法6 。这两种方法编码效率较高，但计算非常复杂，不能适用于实时性要求较高的场合。近年来，人们开始研究将标量、矢量量化相结合的方法，以同时获取较高的压缩比、恢复质量和时间性能，这是图像压缩技术的一个重要发展方向。本文提出的编码方案对图像小波细节子图划分为 44 的块，采用块内的方差作为块的纹理复杂度和重要性度量，对纹理复杂的重要块用较多的位进行编码，而对于较平坦的区域
13、用较少的位进行编码。这实际上相当于将各块的元素组成一个矢量，对不重要的子块采用矢量编码方案，而对重要子块采用标量编码，使得各子块的恢复误差大致平衡。本压缩方法的具体步骤如下。(1) 对图像进行 3 层小波分解，对 LL3 子图进行熵编码，对HH1 不编码(解码时以 0 填充)。(2) 把小波分解图的其它各细节子图按 44 块划分，设定阈值0T0T1T2。将方差小于 T0 的块划分为平坦区；方差大于 T0 小于T1 的块划分为次平坦区；方差大于 T1 小于 T2 的块划分为次纹理区；方差大于 T2 的块划分为纹理区。(3) 按各块在图中的位置进行块的类别编码。由于块的类别数为 4，采用两位编码。
14、对不同类别的块采用下述编码方案： 对于平坦区，假设其均值为 0，可认为块中所有元素均为0； 对于次平坦区，用 1 位表示整个块的均值，对于块中每个元素再各用 1 位进行编码，即块中大于均值的元素对应码号为 1，否则对应码号为 0； 对于次纹理区，用 2 位表示整个块的均值，对于块中每个元素再各用 2 位进行编码，把块中各元素值对应到 03 这 4 个码号上去； 对于纹理区，用 6 位表示整个块的均值，用 6 位表示块内方差，对于块中每个元素再各用 5 位进行编码，把块中各元素值对应到 031 这 32 个码号上去。(4) 对上述结果进行算术编码。上述算法对于原图 44 块的 128 位数据，平
15、坦区只用 2 位编码，次平坦区用 19 位编码，次纹理区用 36 位编码，纹理区用2+165+6+6=94 位编码。5 算法关键问题5.1 小波基的选取多尺度分析中小波基的选择注意 5 个方面的因素。我们选择了4 组小波基对应的二次镜面滤波器(QMF)研究它们的性质： (1) Battle 和 Lemarie 的 27-系数滤波器(简称 B-L 小波)7 ；(2) I. Daubechies 的 4-系数滤波器(简称 D-4 小波)；(3) I. Daubechies的 20-系数滤波器(简称 D-20 小波)8 ；(4) Antonini 的一组双正交小波基对应的滤波器9 。(1) 正交性。
16、用正交小波基由多尺度分解得到的各子带数据分别落在相互正交的 L2(R2)的子空间中，使各子带数据相关性减小。但是能准确重建的、正交的、线性相位、有限冲击响应滤波器组是不存在的，此时一般放宽正交性条件为双正交。上述 4 组 QMF 中除第 4 组外，都是正交小波基。(2) 支撑集。为了得到有限长度的滤波器组 h(n)，g(n)；避免滤波过程中的截断误差，要求小波基是紧支集的。上述 2 组Daubechies 小波和双正交小波都是紧支集的。而 B-L 小波是非紧支集的。(3) 对称性。对称滤波器组具有两个优点： (1) 人类的视觉系统对边缘附近对称的量化误差较非对称误差更不敏感；(2) 对称滤波器
17、组具有线性相位特性，对图像边缘作对称边界扩展时，重构图像边缘部分失真较小，有利于复杂特性的分析(如序列目标检测和分类)9 。上述的 B-L 小波和双正交小波是对称的，两种Daubechies 正交小波是非对称的。由表 3 可见，虽然 B-L 小波基不是紧支集的，但由于其对称性，恢复图像的 PSNR 值仍高于两种Daubechies 正交小波。表 3 4 种小波基的性能比较 Table 3 The comparison of performances of 4 wavelet basis小波基 B-L 小波 D-4 小波 D-20 小波 双正交小波 PSNR 39.9693 32.5134 27
18、.08921 46.280867 (4) 规则性(Regularity)。文献7指出小波基的规则性对最小化量化误差是很重要的。规则度是度量小波函数光滑性的。这 4组小波基中，D-20 小波规则性最好，规则度为 2.9；双正交小波的规则度为 1.410；B-L 小波的规则度为 1；D-4 小波的规则度为0.5。(5) 消失矩阶数。可以证明消失矩阶数很大时，精细尺度下的高频部分数据值有许多是小至可以忽略的(奇异点例外)。上述 D-20小波的消失矩为 10，双正交小波的消失矩为 4，D-4 小波的消失矩为 2，B-L 小波的消失矩为 1。采用上述 4 种小波基，我们分别对图 2(a)作 3 级多尺度
19、分解后不作编码，立即进行重构，结果如表 3(其它图像也有相似结果)。可见，本系统采用的双正交小波基具有良好的性能。 5.2 阈值的选取本方法的一个关键因素是 3 个阈值 T0，T1 和 T2 的选取。直观地说，3 个阈值越大，压缩比越高，而图像恢复质量越差。另外，根据人类的视觉生理、心理特点以及实验结果，不同级别的小波分解系数所含的能量是不同的10 ，因而在图像重构时其重要性也有差异，应区别对待。级别越高，小波系数所含能量越大，量化应越精细。在矢量编码方案中，一般采取级别高的小波子图矢量维数低就是这个道理。通过对不同级别的子图采用不同的阈值可以实现对不同级别子图的区别对待，即级别越高，阈值越小
20、。进一步，上述 3 个阈值的选取有两条途径： 一是由用户根据需要交互地给出，而由系统给出一个较优的缺省值。通过对图 2 的 4幅遥感图像进行实验可以获得各阈值与压缩系统性能指标 PSNR 和CR(峰值信噪比和压缩比)的关系。由于 T0，T1，T2 与 PSNR 和 CR 是一个多对多的关系；不易表达，为简单起见，我们设T2=2T1=2T0。表 4 给出了图 2(a)的 HL2 子图的阈值设定与压缩性能关系，其它图像和其它子图与压缩性能具有类似的关系。由表 4 可见，当 T0=4 时，总体压缩性能 PSNR+CR 和 PSNRCR的变化趋于平缓。对图 2(b)、(c)、(d) 3 幅图，这种情况
21、分别发生在 T0 为 5，7，7；因此我们选 HL2 子图的 T0 值为 6。其它阈值选取方法类似。表 4 图 2(a)的 HL2 子图的阈值设定与压缩性能关系一览表 Table 4 The relation of HL2 threhold and compressionperformance of fig. 2(a)T0 2 3 4 5 6 7 8 PSNR 41.098 40.666 40.532 40.233 39.905 39.904 39.549 CR 13.611 16.866 18.660 18.852 19.593 19.762 20.408 PSNR+CR 54.709 57
22、.532 59.192 59.085 59.498 59.666 59.96 PSNRCR 559.38 685.87 756.33 758.47 781.86 788.58 807.12 另一种方法是通过对小波分解子图进行统计分析后自适应获得，由前面的讨论我们认为这是不必要的，理由有二： (1) 不同阈值的选取以及同一设定对不同图像造成的压缩性能影响不太大；(2) 自适应选取方法时间性能大大降低，不适合于实时性要求较高的场合。5.3 算法的实时性问题目前小波分解已经有快速算法，并可用硬件实现3 ，使研制基于小波分解的实时图像压缩技术成为可能，这对于星上数据压缩具有十分重要的意义。为了提高编码
23、过程的速度，我们没有采用一般用于度量数据能量的方差指标，而代之以 44 子块的块内数据变化范围(即最大最小值之差)，从而减少了一次对块中所有元素的扫描，且避免了求方差时的乘法运算，只需作 16 次浮点数比较(即减法)操作，缩短了编码时间，而图像恢复质量基本没有下降。具体数据见表 5。另外，对于多波段遥感数据，我们先对它们进行 K-L 变换，然后对各 K-L 变换子图使用上述方法进行压缩，效果良好，试验结果见表 6。我们采用 Antonini 的 9-7 双正交小波基。表 2 给出了图 2(a)、(b)使用 3 种压缩方法的实验结果。由表 2 可见，本方法在保持较高的保真度情况下压缩比远高于无损
24、压缩，而压缩比和 PSNR 值均优于JPEG 方法。表 5 分别给出了使用方差和变化范围作为子块重要性度量的比较结果。显然，在性能基本不变的情况下，使用变化范围的方法速度要快约 40%。表 6 给出了一幅 Landsat TM7 波段 512512 大小的遥感数据压缩性能，PCi(i=1，2，3，7)代表 17 主成分数据，总体压缩时我们只使用了前 5 个主成分。可见，该方法在多波段遥感数据的高保真压缩方面也具有良好的效果。表 5 4 幅遥感图像的压缩性能表 Table 5 The compression performance of four test image图像名 基于方差的块重要性判
25、别 基于变化范围的块重要性判别 压缩比 PSNR 压缩时间/s 压缩比 PSNR 压缩时间/s 图 2(a) 16.88 40.23 3.91 16.85 40.23 2.64 图 2(b) 11.94 40.20 4.77 11.97 40.12 3.02 图 2(c) 17.74 39.27 2.97 17.74 39.27 1.98 图 2(d) 15.23 40.70 3.85 15.21 40.73 2.68 注： 以上实验是在 Pentium 133 CPU, 32M 内存的硬件平台上进行的，其中小波变换未采用快速算法，为程序设计方便，还曾将中间结果存入硬盘，且算法未进行优化。考虑
26、这些因素，将压缩时间降到 1s 左右是可能的。 *本文受国家卫星应用重点项目支持.收稿日期： 1998-03-23； 收到修改稿日期： 1998-08-30参 考 文 献 (References)1 K. Sayood. Data Compression in Remote Sensing Application. IEEE Trans. on Geoscience and Remote Sensing Newsletter, 1992, 38(1): 7-15.2 I. Daubechies. Ten lecture on Wavelets. Philadelphies: Capital c
27、ity Press, 1992.3 S. G. Mallat. A Theory for Multiresolution signal Decomposition. IEEE Trans. on PAMI, 1989, 11(7): 674-693.4 Hu Zheng, Yang Youwei. The principle and application of vector quantization. Xian Electronical Technical University Press, 1988. (In Chinese) 胡征，杨有为. 矢量量化原理与应用. 西安： 西安电子科技大学
28、出版社，1988.5 M. Barlaud, M. Antonini, P. Mathieu. Pyramidal Lattice Vector Quantization for Multiscale Image Coding. IEEE Trans. on Image Processing, 1994, 3(4): 367-380.6 J. M. Shapiro. Embedding image coding using zerotrees of wavelet coefficients. IEEE Trans. on Signal Processing, 1993, 41(12): 344
29、5-3462.7 G. Battle. A block spin construction of ondelettes, Part I: Lemarie function. Commun. Math. Phys., 1987, 110(3): 601-6158 I. Daubechies. Ten lecture on Wavelets. Philadelphies: Capital city Press, 1992.9 M. Antonini, et al. Image coding using wavelet transforms. IEEE Trans. on Image Process
30、ing, 1992, 1(2): 205-220.10 A. Averbuch, et al. Image Compression Using Wavelet Transform and Multiresolution Decomposition. IEEE Trans. on Image Processing, 1996, 5(1): 4-15.作 者 简 介李 强，男，1972 年 10 月生，1993 年在中国人民解放军国防科技大学自动化系获学士学位，现为该校自动化系博士研究生，主要从事模式识别、遥感信息处理和地理信息系统方面的研究，已发表论文 4 篇。A Compression Met
31、hod for Remote Sensing ImageLI Qiang WANG Zheng-Zhi(University of Science and Technology for National Defense, Changsha 410073)Abstract With the increasing volumes of remote sensing data, data compression is receiving more and more attention. Adapting to the specialities of remote sensing datathe lo
32、w local correlation and the rich complex texture information, this paper presents an adaptive scalar-vector hybrid quantization method for compression based on wavelet transform. According to textural intensity of every block in wavelet-transformed high-frequency subimage, we classify them into four
33、 classes. Compressing the plain block is at high compression ratio, and the textural block at high fidelity, The method enable the balance of the restore error of every block. This method is time efficient by avoiding the codebook training and searching, while surpass the performence of JPEG for sin
34、gle image. By combining with K-L transformation, which is a kind of correlation-reduction methods, we apply the presented method to multi-band remote sensing image, and good results have been obtained.Key words Wavelet transform, Hybrid quantization, Image compression.- _ D - _,_ V - (*)/ (*) R专业硬盘录
35、像机，摄像机，楼宇对讲系统供应商我们一直在追求，构筑完美结晶。徐金13788580809http:/ 2004-5-24 16:25:28 fire 头衔：小姐你爱吃青椒嘛？ 等级：超级版主 威望：20文章：3617经验：112610魅力：111630现金：113799注册：2004-3-8第 2 楼 看不懂-“小朋友，你父母知道你出来灌水嘛？” “我让小白转告他们了”_E-mail:QQ:13225080 2004-5-24 16:31:59 GEORGE 头衔：流浪军团一等兵 等级：超级版主 威望：200文章：4523经验：266940魅力：264030现金：264099注册：2003-6
36、-3第 3 楼 TriMedia 和 PC 使用的视频压缩和解压缩库原始，1，2Mbyte小波变换，1:237，5KByteDR Wavelet Lib 是基于小波变换的高性能的库，用于图像和视频码流的压缩和解压缩。通过实际的应用例子，可以强有力的证明该库的特点。系统要求 TriMedia 或者 PC 机 Windows 32，至少是 Pentium 133输入 帧YUV 4:1:1 或者 greyscale 格式QCIF176x144CIF352x288仅奇数区域(CCIR_601)704x288奇偶区域(CCIR_601)704x576（在 TriMedia 面板上，奇数区域和奇偶区域至少
37、需要 16MB 的主内存） 输出 生成的位流的长度可以由用户定义。因此，压缩帧的大小和它的分辨率可以达到最理想的状态。算法 在相邻的帧之间使用空间的和时间的冗余，防止数据丢失 执行多级的小波变换 在视频序列中增加相关性 通过使用明确的、合适的算法对变换的帧进行编码 根据数据所归属的帧，对生成的位流数据进行排序 应用 有限带宽连接上的视频传输 视频序列和静止图像的密集存储 基于分块导向技术的小波变换压缩的优点（比如：JPEG） 高压缩速率下无分块损失 通过调整压缩速率很容易的得到恒定的帧速率 预览功能，这是小波算法的固有特点 小波和 JPEG 压缩的比较 JPEG，1:150，6Kbytes 小
38、波，1:150，6Kbytes小波，1:25，36Kbytes原始，900Kbytes，640x512压缩速率 低带宽的算法压缩帧的大小压缩速率Pentium III/450压缩速率TrimMedia TM1100600Byte70,0 QCIF 帧/秒或者22,5 CIF 帧/秒51,5 QCIF 帧/秒或者16,35 CIF 帧/秒1000Byte67,0 QCIF 帧/秒或者17,0 CIF 帧/秒46,6 QCIF 帧/秒或者15,84 CIF 帧/秒中等带宽的算法5KByte43,0 QCIF 帧/秒或者16,5 CIF 帧/秒或者10,2 单区域23,7 QCIF 帧/秒或者12,0 CIF 帧/秒或者7,24 单区域15KByte22,0 QCIF 帧/秒或者11,0 CIF 帧/秒或者8,3 单区域10,5 QCIF 帧/秒或者7,4 CIF 帧/秒或者5,2 单区域
