实用型模块化成像光谱仪多光谱图像的信噪比估算及压缩方法研究

实用型模块化成像光谱仪多光谱图像的信噪比估算及压缩方法研究
摘要： 采用局部标准差法和去相关法对实用型模块化成像光谱仪多光谱图像的信噪比进行估算。这两种方法已 将地物变化的影响降低到很低的程度。这样在大气订正后，图像的信噪比性能充分反映出遥感仪器的信噪比性 能。针对图像压缩，提出控制各波段恢复图像的峰值信噪比刚好大于原始图像的信噪比，使由压缩算法本身所带 来的噪声限制在原始图像的噪声范围之内。结合这种压缩思想，用基于离散余弦变换和基于离散小波变换的压缩 算法，对实用型模块化成像光谱仪多光谱图像进行压缩。实验表明，利用这种方法，对于高信噪比的波段，图像信 息得到了保真；对低信噪比的波段，压缩倍数提高迅速且恢复图像视觉无失真，对整幅成像光谱图像，压缩性能提 升显著— — 当压缩比等于37．95倍时，峰值信噪比等于45．86 dB。

1 引 言 实用型模块化成像光谱仪一 (Operational modular imaging spectrometer，OMIS)由中国科学 院上海技术物理研究所研制，是一种带有稳定平台 的新型机载遥感器，具有128个波段，覆盖了从可见 光到热红外的光谱范围，采用12 bit编码。通过对多 光谱图像在空间维和光谱维方向的数据特性进行分 析 ]，可知多光谱图像内部存在很大的相关性，为多 光谱图像的压缩奠定了理论基础。为了方便数据的 广泛应用和仪器的优化改进，转速计| 水份计| 水份仪| 分析仪| 溶氧计| 电导度计| PH计| 酸碱计| 糖度计| 盐度计| 酸碱度计| 电导计| 水分测定仪|对遥感器获取的成像 光谱数据要有定量的质量评价。各通道数据的信噪 比是一个非常重要的评价量，该度量对用户来说是 必需的，因为噪声的大小决定了识别地物光谱中吸收 特征所能达到的精度，以及地物识别的精度。而要识 别一个位置的吸收特征，要求噪声的幅度必须小于该 吸收峰的深度。所以噪声大小对于判断某个吸收特 征的识别是否可能是必需的。仅仅只有噪声信息还 不够，因为同样水平的噪声在信号弱时对数据质量的 影响要比信号强时大，所以必须考虑信噪比。 本文采用局部标准差法和去相关法两种信噪比 分析方法，对实用型模块成像光谱仪图像的各波段 数据的信噪比进行了估算，指出基于空间和光谱维 去相关的方法比较可靠，并给出了实际的计算结果。 根据实验结果，本文提出在图像压缩当中，应控制各 波段恢复图像的峰值信噪比(Peak signal—to-noise ratio，PSNR)大于原始图像的信噪比(Signal—to— noise ratio，SNR)，使由压缩算法本身所带来的噪声 限制在原始图像的噪声范围之内，对图像整体质量 的影响减至最小。结合这种压缩思想，用基于离散 余弦变换(Discrete cosine transform，DCT)的压缩 算法和基于离散小波变换(Discrete wavelet transform，DWT)的压缩算法(这两种算法结构简 单，计算量小，便于在数字信号处理器(DSP)上实现 对图像的实时压缩[3 )，对实用型模块成像光谱仪多 光谱图像进行压缩。实验结果表明，利用这种方法， 对于高信噪比的波段，图像质量得到了保真；对低信 噪比的波段，压缩得到很大提高，对整幅成像光谱图 像，它的压缩倍数就上去了，而图像的有用信息也几 乎完整地被保留下来了。 基于离散余弦变换的压缩算法参考了国际标准 化组织(ISO)发布的基于离散余弦变换的JPEG[4 (Joint photographic experts group)压缩算法思想， 先进行离散余弦变换正变换，再对离散余弦变换系 数进行量化，并对量化后的直流系数和交流系数分 别进行差分编码或游程编码，最后再进行熵编码，原 理框图如图1所示。 Fig．1 DCT—based JPEG compression algorithm 基于离散小波变换的压缩算法先将图像分解为 低频分量sJ及高频分量{w ，w；，w；) + ～一 。 一级小波变换将原始图像分解成一个低频信号和三 个方向的高频分量信号(水平方向、垂直方向及对角 方向)，即每一层分解为四个子带信号，低频信号又 可分解成四个子带，故总的子带数为3K+1，其中 K 为分解层数，如图2所示。对所得到的低频分量 及高频细节分量{w ，w ，w；)，一J+卜一。，根据人类 的视觉生理特性分别作不同策略的量化与编码处 理。例如，对于低频分量采用快速余弦变换(离散余 弦变换)结合“之”字形扫描、熵编码方法(如 Huffman编码、算术编码、矢量量化等)进行压缩。 对于高频细节分量可以采用量化，去掉人眼不敏感 的高频成分并结合熵编码方法的压缩方法。 s I II ． Fig．2 Coefficient distrubtion of 2DW T 2 多光谱图像信噪比分析 噪声主要有两大类：加性噪声和乘性噪声。加 性噪声的大小和信号大小无关，而乘性噪声则正好 相反。跟大多数成像光谱图像一样，实用型模块成 像光谱仪图像的局部方差和局部均值之间没有明显 的趋势，从而验证了实用型模块成像光谱仪信号中 的噪声也主要是加性噪声，而不是乘性噪声(限于篇 幅，不再展开论证，详见参考文献[1])；另外，实用型 模块成像光谱仪信号中的噪声自身是不相关的，主 要是随机噪声。实用型模块成像光谱仪，作为一种 具有128波段、瞬时视场为3 mrad、光谱分辨力为 10 nm的新型遥感器，系统非常复杂，噪声因素很 多，各因素的特性也多种多样。为了计算数据的信 噪比，不能直接从噪声特性入手，而应从信号的相关 性着手，运用去相关的方法可靠地求出各波段的噪 声大小，进而求出信噪比信噪比。下面就噪声分析 方法逐一介绍。

2．1 局部标准差法 该方法使用了局部均值(Local mean，LM)和局 部标准差 5](Local standard deviation，LSD)的概 念，其基本思想是：由于选择一定大小的均匀区域比 较困难，那么就把图像分割成一个一个的小区域，这 些小区域内基本上可以认为是均匀的；分别计算这 些小区域内的局部标准差作为局部噪声大小，并选 择众数最多的那个区间的局部标准差作为整个图像 的平均噪声值。这种方法可普遍适用于各种遥感图 像的信噪比估算，但仍在一定程度上受地物目标均 匀程度的影响。具体的操作步骤如下： 1)将图像分割成4×4，或5×5，. . ，或8×8 的小块，对于每一个图像子块，信号的局部均值ML 由下式得到 1 i 一 = N ML一 ：Sf， (1) 』 这里S 是图像子块中第i个像素的灰度值；N 是图 像子块中所有像素的总数。局部标准差由下式得到 1 = N I／2 D 一『L ‘’ 【_Σ(s —ML) ]． (2) t=1 j 对于均匀的图像子块，局部标准差较小，而对 不均匀的图像子块，如包含图像边缘或纹理特征的 子块，局部标准差则较大。计算出整幅图像的局部均 值(记为 )、所有图像子块的局部标准差，并找出 所有图像子块中最大和最小的局部标准差。 2)在最小和最大的局部标准差之间，建立若干 个等值间隔的区间。将所有子块的局部标准差按照值 的大小依次排入相应的区间。对每个区间的局部标准 差的个数进行计数，计数值最大的那个区间的局部标 准差的平均值即为整幅图像的噪声，记为跣。 3)由下式可求得整幅图像的信噪比R ： ^ RsN一201g ． (3) i Ls 以海水、荒漠、城市、农田4幅不同地物特征的 图像作为测试目标，进行实验，并计算出信噪比，如 图3所示：横轴表示波段号，纵轴表示信噪比。 这种分析方法的原理决定了其计算结果会受地 物目标、计算区域的子区域的均匀程度等因素的影 响，因而不能精确反映图像的真实信噪比，跟图像的 真实信噪比会有少许偏差。 11期 蒋青松等： 实用型模块化成像光谱仪多光谱图像的信噪比估算及压缩方法研究

2．2 去相关法 成像光谱仪是个复杂的系统，在信号中引入噪 声的因素很多。直接从噪声特性着手要估计噪声的 大小有一定的难度。另一方面，成像光谱仪高空间 分辨率和高光谱分辨率的特性，加上地物的连续性， 使得成像光谱数据在空间和光谱维上均具有很高程 度的相关性。去相关的方法则是利用了这一点，将 原始数据中，这种具有高相关性的信号去掉，余下的 就是噪声。具体去相关的方法是采用多元线性回 归[6j，即 2 ， ， — a~ri， ， 1 4-bxl， ，蚪1 4- c p，^4- d， (4) 其中2为图像灰度的线性拟合值，z为图像灰度值，i， 为图像二维空间坐标，分别代表行和列方向的坐标 值，五为波段序号。口、6、c、d为线性回归系数。并且有 rz l，j，^ i> 1 I z 一_ 1， (5) 【无意义， i一1， 一1 去相关后的残差图像为 ，．“j，l— zf，J，l一2“j，^， (6) 0 口 其中r为残差值。最后，噪声方差为 h 一 ，j． ， N 一硼× h一1， (7) 其中， 、h分别为图像的宽和高，r为图像去相关后 的残差值。 记第走波段图像均值为 h 一— w X L h zllj，¨ (8) 则第走波段图像的信噪比为 RsN一201g ． (9) tT 仍以海水、荒漠、城市、农田4幅不同地物特征 的图像作为测试目标进行实验，并计算出信噪比，如 图4所示。考虑到实际的实用型模块成像光谱仪原 始图像在光谱维的不够平稳性，原始的实用型模块 成像光谱仪图像数据在光谱维的相关性并不十分理 想，如果对原始图像数据进行增益校正，提高光谱维 的相关性，则所求信噪比将会更接近实际值。 ，．、，、 ～ 一 二 苎． 一 一 ’ 一一 一一 J ’ ’I ． J 王 -一-Field ] ____-City ·---一DCSCrt 。。。·一Sea I 1 FIg．4 Signal—to-noise ratio calculated by de-correlation method 这种方法基本不受地物目标、区域均匀程度的 3 基于信噪比分析结果的实用型模块 的!实 际信噪比煮情 况好，因而可以较真实地反映图像 成像光谱仪图像压缩思想 。 . . ’ 。 从计算数据结果看，我们知道原始图像本身质量 并不是太好，第1，65，81等波段的图像，因为处于探 测器探测波段的初始波段，探测响应率低，其信噪比 要比相邻波段图像的信噪比低很多；在可见近红外波 段(Ⅵ )，中红外波段(MIR)及热红外波段(TIR)， 图像的信噪比较高，也大多都在50 dB以下；而在短 波红外(s、Ⅳ1R)的一些波段，如65～96波段的信噪比 在35 dB以下。因此在用数据压缩算法对图像进行压 缩时，如果充分利用图像本身的这一特点，采用适当 的压缩策略，就可以既保证图像有用信息基本无损 失，又可以使图像得到较大程度的压缩。 与一般的静止图像相比，遥感图像数据的纹理 更丰富，具有高熵值和低冗余度的特点；另外，由于 遥感图像应用领域的特殊性，对图像数据的压缩不 能仅仅满足于视觉上的清晰，还要求图像信息得到 最大程度的保真，这就产生了压缩倍数和信息保真 的矛盾。为解决这一矛盾，结合图像本身的噪声污 染，我们提出如下压缩思想： 在图像压缩当中，对于信噪比较高的波段，由于 噪声干扰因素较小，图像质量较好，为尽量保真图 像，需要恢复图像的峰值信噪比较高；而对信噪比较 低的波段，由于噪声干扰严重，图像本身质量差，就 可让恢复图像的峰值信噪比较低，总之，是要控制各 波段恢复图像的峰值信噪比大于原始图像的信噪 比，这样就可以将由压缩算法(比如由量化所带来的 噪声)本身所带来的噪声限制在原始图像的噪声范 围之内，不会污染到高比特域的有用图像信息，而低 比特域的信息因为已淹没在噪声当中无法辨别，有 损失也只能是图像信息和噪声同时损失，对图像整 体质量基本不会产生影响。这样，对于高信噪比的 波段，压缩就少了一些；对低信噪比的波段，压缩就 多了一些，对整幅成像光谱图像，它的压缩倍数就上 去了，而图像的有用信息也几乎完整地被保留下来 了，是一种比较理想的压缩策略。

4 实验方法、结果和结论 本文选用了基于离散余弦变换的压缩算法 和 基于离散小波变换的压缩算法来实现图像压缩。基 于离散余弦变换编码的过程先进行离散余弦变换正 变换，再对离散余弦变换系数进行量化，并对量化后 的直流系数和交流系数分别进行差分编码或游程编 码，最后再进行熵编码。其中，量化器量化步长的大 小将决定图像的压缩倍数(压缩比)及恢复图像的质 量(峰值信噪比)，量化步长越大，则压缩比越大，峰 值信噪比越小，图像质量越差。 基于离散小波变换的压缩算法选用了8阶 Daubechies小波基，对实用型模块成像光谱仪成像 光谱图像进行三层小波分解(三级小波变换后的图 像系数分布如图5所示)，并用零树编码的方法实现 对图像的压缩。跟基于离散余弦变换的JPEG算法 一样，可通过调整量化步长控制图像的压缩倍数(压 缩比)及恢复图像质量(峰值信噪比)。 Fig．5 Im age coefficient distribution after 3 layers of 2DW T algorithm 用基于离散余弦变换的压缩算法和基于离散小 波变换的压缩算法，并结合第三节中提出的压缩思 想，分别对实用型模块成像光谱仪拍摄的石家庄城 市和农田图像进行压缩实验，结果如图6所示。图 中：粗实线代表实用型模块成像光谱仪原始图像的 信噪比，细实线代表恢复图像的峰值信噪比，虚线代 表图像压缩倍数；横轴代表波段号，左纵轴代表峰值 信噪比、信噪比，右纵轴代表压缩比。 从图6可以看出，在短波红外的65～96波段， 原始图像信噪比较小，图像质量较差，本文对这些波 段进行了高倍压缩，而仍然控制峰值信噪比大于信 噪比，使图像的整体压缩倍数得到很大提高，且图像 质量并无多大损失；而在短波红外之外的可见近红 外、中红外、热红外波段，原始图像信噪比较大，图像 质量好，本文对这些波段只进行较小倍数的压缩，使 恢复图像高度保真。整体上看，这种压缩方法使图 像的压缩倍数上去了，图像的有用信息也基本完整 地保留下来了，是一种比较理想的压缩方法。总体 压缩效果如表1所示。 ∞ Z ∽ Z ∽ BalldN0 E ’二 U BandNO E ’二 U Fig．6 128 bands OM IS image compression results． (a)DCT-based compression for city image： (b)DW T-based compression for field image 图7是用以上压缩方法压缩的部分波段图像， 其中，石家庄城市图像用了基于离散余弦变换的压 缩算法，石家庄农田图像用了基于离散小波变换的 压缩算法。图中，a、b、C、d分别代表石家庄城市第 32、75、105、118波段图像，e、f、g、h分别代表石家庄 农田第32、75、105、118波段图像，z一1代表原始图 像，z一2代表恢复图像(z表示a～h)。