CN101483777B - 基于自适应多尺度Bandelet包的SAR图像去噪压缩方法 - Google Patents

Abstract

本文公开了一种基于自适应多尺度Bandelet包的SAR图像去噪压缩方法，它涉及遥感图像处理技术领域。其实现过程是：1)读入复值SAR图像的模值图像M进行对数变换，得到对数变换后的图像M′＝log(M+1)；2)将M′进行自适应多尺度Bandelet包变换，并进行软阈值法滤除斑点噪声，得到滤除噪声后的系数；3)将滤波后的系数进行逆自适应多尺度Bandelet包变换得到M″，再进行反对数变换得到去噪后的图像M_denoise＝e^M″-1；4)对去噪后的图像M_denoise再执行自适应多尺度Bandelet包变换，得到Bandelet系数；5)对Bandelet系数先进行量化和编码，再进行逆自适应多尺度Bandelet包变换，得到去噪压缩后的SAR图像。本发明具有压缩后图像峰值信噪比、等效视数和边缘保持系数高，有效保持了SAR图像高频纹理信息和目标几何结构，以及复杂度低、实现快速等优点；在保证图像质量不影响后续分析与解译的前提下，能获得较高的压缩比，可用于高分辨SAR图像的存储与传输处理。

Description

基于自适应多尺度Bandelet包的SAR图像去噪压缩方法

技术领域

本发明属于遥感图像处理技术领域，涉及合成孔径雷达(SAR)图像去噪压缩方法。该方法可用于高分辨SAR图像的存储和传输处理中。

背景技术

SAR是一种可全天候、全天时工作的遥感成像设备，可以对地面目标进行高分辨率、大面积成像，具有成像机理复杂、数据率高和图像数据量大的特点。例如一个实用的星载SAR系统的数据率大约为200MbPs，高速率、大数据量传输数据的下传相对于有限的带宽来说是非常困难的。尽管有一些方法如降低雷达PRF、降低A/D变换器的位数和系统带宽等能够降低数据率，但都会使系统性能下降。因此，必须对数据量巨大的SAR图像进行压缩。

相比于其它遥感图像，SAR图像具有如下优点：

1)高空间分辨率：包括距离分辨率和方位分辨率，可获得与光学图像相当的高分辨率图像。目前SAR的分辨能力已经可以达到0.3m(APG-76雷达)，但仍未达到其物理极限，在未来一段时间内，SAR的成像分辨率将会更高；

2)强穿透能力：除了能穿透云雾之外，对一些地物介质，如岩石、土壤、松散沉积物、植被、冰层等，有穿透一定深度的能力。因此，它不仅反映地球表面的信息，还可以在一定程度上反映地表以下物质的信息；

3)立体效应：雷达散射及雷达波束对地面倾斜照射，产生雷达阴影，即图像暗区，这种明暗效应能增强图像的立体感，这种明显的地形起伏感，对地形、地貌及地质构造等信息有较强的表现力和较好的探测效果；

4)几何特征：可用于进行地形、地物的测量和分析等；

5)对与水有关信息的识别能力更强，对松散沉积物的表面结构反映明显，对居民点及线性地物的表现尤为明显；

6)SAR图像与其它遥感图像，例如红外图像、可见光图像、多光谱图像等之间具有较好的互补性。

SAR图像反映的是被探测区域的雷达后向散射特性，所以具有一些不同于人物、景物和其它光学图像的特点。首先，由于SAR采用相干成像，SAR图像受到由成像机理所决定的乘性相干斑点噪声的污染；第二，由于不同的地表粗糙度对雷达回波的后向散射系数差别较大，SAR图像呈现出丰富的纹理特征；第三，由于地表的雷达后向散射系数的差异很大，SAR图像的动态范围很大。这些特点导致许多普通光学图像压缩算法应用到该领域时存在很大的局限性；另一方面，SAR图像的某些应用对于重建图像的准确性有更高的要求，从而对于图像压缩算法的设计提出更高要求。因此，SAR图像的处理不能简单采用光学图像的处理手段。

斑噪现象使得SAR图像不能正确反映地物目标的散射特性，严重影响了成像质量。因此在对SAR图像进行压缩之前，希望先对图像进行斑点噪声抑制；其次，在SAR图像中，地物目标的几何特征被削弱和淡化，它们的波谱特征主要以不同的亮度值，即灰度差异综合反映出来。有用信息与无用信息相对于不同专业和用途而言混杂在一起，相互干扰，彼此抑制。因此，希望能兼顾均匀区域噪声的滤除和各类目标信息的清晰；另外，较大的动态范围要求SAR图像压缩算法不能只考虑全图的统计特性，而要能自适应于动态范围大的局部图像统计特性；再次，SAR图像的纹理特征对于后续的解译较为关键，如地物判别、边缘检测和图像分类，因此SAR图像压缩时必须尽可能地保留纹理信息；最后，所获取的SAR图像对目标的几何特性，无论是微观的还是宏观的，都非常敏感，希望所设计的压缩算法能够较好的保持目标的这些几何特性。

分析SAR图像的统计特性可知，SAR图像原始数据在距离向和方位向上存在较弱的相关性。经过正交变换，强相关的空间采样值变换为不相关或弱相关的变换系数，这样就可以去除存在于相关性中的数据冗余度。因此，以去相关为基础的压缩算法可用于SAR图像压缩。本发明考虑SAR图像的特点，希望设计有效保持SAR图像各个层次纹理和目标几何特征的正交变换，对SAR图像进行有效的滤波与压缩，使无用信息与噪声等干扰因素得到抑制与排除，有用信息得到增强与突出，进而在保证图像质量不影响的前提下，得到尽可能高的压缩比，以利于图像的进一步处理、分析与解译，最终为SAR系统中卫星到地面的高数据传输瓶颈问题提供一条高效可行的解决方案。

发明内容

本发明的目的在于克服SAR图像已有压缩技术的不足，提供一种基于自适应多尺度Bandelet包的SAR图像去噪压缩方法，对原始的SAR图像复数据进行分析、滤波与压缩，抑制噪声与无用信息，保证图像质量的同时得到尽可能高的压缩比，为SAR系统高数据传输瓶颈问题提供有效的解决方案。

为实现上述目的，本发明的SAR图像去噪压缩方法，包括如下步骤：

1)读入复值SAR图像的模值图像M进行对数变换，得到对数变换后的幅值图像M＝log(M+1)；

2)将对数变换后的幅值图像M′进行自适应多尺度Bandelet包变换，并采用软阈值法滤除斑点噪声，得到滤除噪声后的系数；

3)将滤波之后的系数进行逆自适应多尺度Bandelet包变换，得到M″，对M″进行反对数变换，得到去噪后的图像M_denoise＝e^M″-1；

4)对去噪后的图像M_denoise再进行自适应多尺度Bandelet包变换，得到图像的Bandelet系数；

5)对Bandelet系数先进行量化和编码，再进行逆自适应多尺度Bandelet包变换，最终得到去噪压缩后的SAR图像。

本发明由于充分利用了自适应多尺度Bandelet包对图像几何正则信息的自适应稀疏表示，不仅可以用较少的比特数来表示图像，而且能更好保持图像中目标的几何结构信息，因而在相同的Bpp下能达到较好的视觉压缩效果和较高的峰值信噪比；同时，由于本发明充分利用了多尺度分解后各个高频子带的信息，故能够较好的保持SAR图像的纹理信息；此外，由于本发明采用固定大小的子块剖分，以及线性复杂度的最佳几何流方向搜索方法，因此算法的实现具有较低的时间和空间复杂度。

附图说明

图1是本发明的实现流程图；

图2是图像剖分子块的几何流方向与Lagrange函数值关系图；

图3是本发明与第二代Bandelet所确定几何流方向对比图；

图4是本发明在2级和3级分解时与现有两种方法得到的恢复图像PSNR对比图；

图5是本发明与现有两种方法得到的恢复图像的等效视数和边缘保持系数对比图。

具体实施方式

参照附图1，本发明的具体实现步骤如下：

步骤1，读入原始的复值SAR图像的模值图像M，并对其进行对数变换。对输入的模值图像M按如下公式进行对数变换，得到幅值图像M′：

M＝log(M+1)。(1)

步骤2，将幅值图像M′进行自适应多尺度Bandelet包变换，得到原图的Bandelet系数。

(2.1)对幅值图像M′进行二维子波包变换，分解层数一般取2～3，得到一系列多尺度图像{C₁，C₂，…}。

(2.2)对于分解得到的若干多尺度子带C_k(k＝1，…)进行如下变换：

首先，将大小为P×P的子带C_k剖分为m×m的子块；

其次，在每个子块中根据Lagrange函数最小化的原则，采用折半查找法确定最佳几何流方向，最佳方向满足使得Lagrange函数

取值最小。式中f和f_R分别为原图和恢复后图像，R_jS表示编码四叉树的比特数目，R_jG表示编码每个剖分子块中的几何流的比特数，R_jB表示编码量化Bandelet系数的比特数目；

最后，在每个子块中搜索得到的最佳方向上计算弯曲子波变换，得到各个子带的Bandelet系数。

(2.3)将得到各个子带的Bandelet系数进行组合，得到原图的Bandelet系数B_e。步骤3，对Bandelet系数B_c进行软阈值法滤除斑点噪声。

(3.1)在Bandelet变换域内使用软阈值法滤除斑点噪声，得到滤除噪声后的系数；

滤除斑点噪声公式为：

D(B_c)＝sgn(B_c)×max(0，|B_c|＞T_d)(3)

其中B_c为Bandelet系数，T_d为设定的阈值，D为滤除噪声后的系数；

(3.2)对滤除噪声后的系数D进行逆自适应多尺度Bandelet包变换得到M″，再进行反对数变换，得到去噪后的图像：

M_denoise＝e^M″-1。(4)

步骤4，对去噪后的图像M_denoise执行自适应多尺度Bandelet包变换。

对去噪后的图像M_denoise进行如步骤2所述的自适应多尺度Bandelet包变换，得到去噪后图像的Bandelet系数。

步骤5，对去噪后图像的Bandelet系数进行量化和编码。

(5.1)对系数采用一致均匀量化，设x表示Bandelet系数，T为量化阈值，则量化值Q(x)为，

$Q\left(x\right)=\left\{\begin{array}{cc}0,& \left|x\right|\&le;T\\ \mathrm{sgn}\left(x\right)(q+\frac{1}{2})T,& \mathrm{qT}\&le;\left|x\right|<(q+1)T\end{array}\right.---\left(5\right)$

其中q∈N；

(5.2)采用熵编码对Bandelet系数进行编码。

步骤6，对编码后的图像进行逆自适应多尺度Bandelet包变换，得到去噪压缩后图像。

上述步骤(2.2)中所述的根据该Lagrange函数最小化的原则选取几何流方向，依图2所示的剖分子块的几何方向与Lagrange函数值关系图确定搜索策略，其中图2(a)和图2(b)是对于如图3所示的几何规则图像采用固定剖分8×8的子块后，任取两子块，其中定义的Lagrange函数在(0，π)范围内的函数值；图2(c)和图2(d)是图像采用固定剖分16×16的子块后，任取两子块，其中定义的Lagrange函数在(0，π)范围内的函数值。从图2中可以看出：在不同大小的子块中，两副图像Lagrange函数值与角度都可近似的认为服从二次函数关系。同时考虑到折半查找法能有效的求解二次函数的极值，并且具有线性的复杂度等特点，因此本发明采用折半查找法来实现低复杂度的最佳几何流方向的搜索。记

为Lagrange函数，初始化(min，max)在(0，π)的范围内，折半查找法实现步骤如下：

A.计算

B.若

①当

min←mid，

②当

max←mid；

C.若

则max←mid；

D.若

则min←mid；

E.若min＜max-2，则转到A；

否则有

对应的方向即为最佳几何方向。

本发明对于几何规则图像其最佳几何流的查找结果如附图3(c)所示，对比另外两种方法所查找的几何流的效果，如附图3(a)和附图3(b)所示，本发明方法查找几何流的准确程度有所改进。

本发明的效果可通过以下仿真实验和对比客观评价指标进行说明：

1.本发明方法与第二代Bandelet方法分别对几何规则图像计算得到的几何流结果。

图3给出了第二代Bandelet方法和本发明提出的方法对于一副几何规则图像计算得到的几何流，其中图3(a)为第二代Bandelet方法在最大和最小分块分别设置为4×4和32×32得到的结果；图3(b)为第二代Bandelet方法在最大和最小分块分别设置为4×4和16×16得到的结果；图3(c)是用本发明方法所计算出的几何流结果。从图3(a)和图3(b)可以看出：最大和最小分块的大小对于几何流的确定有着较大的影响。由于在第二代Bandelet中采用的是近似的直线几何流表示，当图像分块越小时，对于几何流的位置和方向的定位就越准确，但遍历各个尺度和子块合并所需的计算量越大，并且编码几何流所需的比特数也会增多；当图像分块变大时，得到的方向和位置的定位精度有所降低，编码数目减少，但子块经方向离散化后得到的方向数目增多，使得优化所需的计算量大。因为图像中的几何规则性是预先未知的，对于大小为n×n的图像通常需要遍历各个尺度j(j＝2，…，logn-1)，即最大图像分块设为n/2×n/2，最小图像分块设为2×2。在第j个尺度上每个2^j×2^j子块中确定最佳的几何流方向，然后再执行子块合并，进而确定最佳四叉树，算法实现的计算复杂度较高。比较来说，图3(b)中选择的几何流要比图3(a)相对恰当一些，但对于几何流位置和方向的定位仍不太准确。图3(c)给出了本发明的结果，可以看出：本发明方法通过进一步分解图像的各高频子带，在多个高频子带中搜索确定图像的规则几何结构，从而能比图3(a)和图3(b)采用的第二代Bandelet方法定位出图像更精确的几何结构。最为重要的是，本发明由于采用了线性复杂度的折半查找方法来搜索最佳几何流，查找的速度也有着明显的改进；此外，固定大小剖分的方式能够避免如图3(a)和图3(b)所示的第二代Bandelet中的最佳四叉树搜索，在运行速度上相比第二代Bandelet有了较明显的提高。

2.本发明在客观指标和视觉效果上的仿真结果。

采用图4中的图4(a)、图4(b)、图4(c)和图4(d)所示的来自于USC-SIPI Image Database，Univ.of Southern California Signal and Image Processing Institute的两幅大小为512×512的SAR图像，使用本发明的去噪压缩方法进行处理，并使用客观评价指标峰值信噪比PSNR评定结果。峰值信噪比越大，说明去噪压缩取得的效果越好。假定两幅图像的大小为X×Y，令f(x，y)表示原始图像，表示解码后的图像，客观评价指标峰值信噪比PSNR定义为：

$\mathrm{PSNR}=101g\left\{\frac{{255}^2}{\frac{1}{\mathrm{XY}}\underset{x=1}{\overset{X}{\mathrm{\&Sigma;}}}\underset{y=1}{\overset{Y}{\mathrm{\&Sigma;}}}{[f(x,y)-\hat{f}(x,y)]}^2}\right\}---\left(6\right)$

实验结果如图4所示，其中图4(a)是对SAR图像A进行2级分解，三种方法在不同压缩码率Bpp下得到的PSNR结果对比图；图4(b)是对SAR图像A进行3级分解，三种方法在不同压缩码率Bpp下得到的PSNR效果对比图；图4(c)是对SAR图像B进行2级分解，三种方法在不同压缩码率Bpp下得到的PSNR效果对比图；图4(d)是对SAR图像B进行3级分解，三种方法在不同压缩码率Bpp下得到的PSNR效果对比图。通过压缩的效果对比，可以清楚的看到本发明方法相对子波方法和第二代Bandelet方法，在客观评价指标PSNR上均具有优势，在分解级数较高时，本发明方法相对第二代Bandelet方法的优势更为明显。

除了客观评价指标峰值信噪比PSNR之外，通常还采用客观评价指标边缘保持系数EPI和客观评价指标等效视数ENL对去噪压缩后的图像进行评价。客观评价指标边缘保持系数EPI的定义为：

$\mathrm{EPI}=\frac{\mathrm{\&Sigma;}\left(\right|P_e(i,j)-P_e(i+1,j)|+{|P}_e(i,j)-P_e(i,j+1)|)}{\mathrm{\&Sigma;}\left(\right|P_b(i,j)-P_e(i+1,j)|+|P_b(i,j)-P_b(i,j+1)\left|\right)}---\left(7\right)$

其中P_e(i，j)是平滑图像像元，P_b(i，j)是原始图像像元，P_e(i，j)和P_b(i，j)都位于边缘区域。边缘保持系数EPI越大，说明恢复图像的边缘保持效果越好。客观评价指标等效视数ENL定义为：

ENL＝μ²/σ²(8)

式(8)中μ和σ分别表示图像的局部均值和标准差。抑制性滤波后，斑点噪声的强度减小，图像的方差减小，客观评价指标等效视数ENL值增大。噪声抑制效果越好，客观评价指标等效视数ENL值越大。

图5给出了对来自USC-SIPI Image Database，Univ.of Southern California Signal andImage Processing Institute的SAR图像A采用三种方法进行压缩，得到的两种客观评价指标等效视数ENL和边缘保持系数EPI的曲线对比图。其中图5(a)是在不同压缩码率Bpp下，三种方法得到的等效视数ENL对比曲线；图5(b)是在不同压缩码率Bpp下，三种方法得到的边缘保持系数EPI的对比曲线图；通过图5(a)可以看出，本发明的方法得到的去噪压缩图像，其客观评价指标等效视数ENL和第二代Bandelet方法相当，但是明显优于子波方法；通过图5(b)可以看出，本发明的方法得到的去噪压缩图像，其客观评价指标边缘保持系数EPI在低比特率时优于子波方法，略高于第二代Bandelet方法。

表1给出了选择不同的阈值T时三种方法去噪后图像的客观评价指标等效视数ENL、客观评价指标边缘保持系数EPI和客观评价指标峰值信噪比PSNR的具体结果。表中2-Bandelet为第二代Bandelet，黑体为三者中性能最佳。

表1三种方法去噪后SAR图像各评价指标比较

	ENL	EPI	PSNR(dB)	T
					本发明	1.8710	0.9671	43.5910	0.02
2-Bandelet	1.8719	0.9660	43.3051	0.02
					Wavelet	1.8677	0.9664	42.5219	0.02
本发明	1.9626	0.7443	28.2436	0.2
					2-Bandelet	1.9645	0.7331	28.0315	0.2
Wavelet	1.9199	0.7272	26.6627	0.2
					本发明	2.0986	0.2229	18.5484	1.6

	ENL	EPI	PSNR(dB)	T
					2-Bandelet	2.0986	0.2227	18.5466	1.6
Wavelet	1.9115	0.2872	16.3249	1.6
					本发明	2.0938	0.1871	17.8154	2
2-Bandelet	2.0938	0.1871	17.8148	2
					Wavelet	1.9027	0.2549	15.1308	2

从表1中可以看出：本发明的去噪压缩方法得到的PSNR均高于现有的第二代Bandelet方法和子波方法，在客观评价指标等效视数ENL、客观评价指标边缘保持系数EPI方面和第二代Bandelet基本相当，但和子波方法相比，本发明方法具有明显的优势。

3.本发明在运行速度上的对比结果。

相比第二代Bandelet，本发明除了能够更好的保持SAR图像高频的纹理信息和目标的几何结构特征之外，另外一个优势就是在于其快速的实现速度。比较使用子波方法、第二代Bandelet和本发明对附图4中所示的两幅SAR图像去噪之后进行压缩的速度，得到的实验结果如表2所示。实验在Intel Core2PC，3.06G/1G PC机的硬件环境和Matlab2006R的软件环境上独立运行10次，取各种方法运行所需时间的平均结果，如表2所示。

表2三种方法压缩时间对比

	子波	第二代Bandelet	本发明
				图4(a)(2级分解)	4.2361s	3258.9s	7.8669s
图4(c)(2级分解)	4.5024s	3577.0s	7.9056s
				图4(b)(3级分解)	4.1778s	4947.3s	24.6600s
图4(d)(3级分解)	4.3045s	5456.1s	25.3836s

从表2可以看出，子波方法所需的时间最短，本发明方法相比第二代Bandelet变换，耗时大为减少；比起第二代Bandelet方法，本发明的实现速度得到了质的提高。这种提高源于两个方面，一方面，由于在子波包分解得到的若干多尺度子带上采用固定大小的剖分方式，避免了最佳四叉树的构造；另一方面，采用具有线性复杂度的折半查找法代替了第二代Bandelet中的穷尽搜索，因此实现速度上有较显著的提高，这更有利于实际系统中的实时实现。

Claims (4)

1.一种基于自适应多尺度Bandelet包的SAR图像去噪压缩方法，包括如下步骤：

(1.1)读入复值SAR图像的模值图像M进行对数变换，得到对数变换后的幅值图像M′＝log(M+1)；

(1.2)将幅值图像M′进行自适应多尺度Bandelet包变换，并进行软阈值法滤除斑点噪声，得到滤除噪声后的系数；

(1.3)将滤除噪声后的系数进行逆自适应多尺度Bandelet包变换得到M″，接着再进行反对数变换得到去噪后的图像Mdenoise＝eM″-1；

(1.4)对去噪后的图像M_denoise再进行自适应多尺度Bandelet包变换，得到Bandelet系数；

(1.5)对Bandelet系数先进行量化和编码，再进行逆自适应多尺度Bandelet包变换，得到去噪压缩后的SAR图像。

2.根据权利要求1所述SAR图像去噪压缩方法，其中步骤(1.2)和步骤(1.4)中所述的自适应多尺度Bandelet包变换，实现步骤如下：

(2.1)对幅值图像M′或去噪后的图像M_denoise进行二维子波包变换，分解层数一般取2～3，得到一系列多尺度图像{C₁，C₂，...}；

(2.2)对于分解得到的若干多尺度图像C_k(k＝1，...)，执行固定子块的剖分，在每个剖分的子块上搜索最佳几何流方向，然后在该方向上执行弯曲子波变换，得到各子带的Bandeiet系数；

(2.3)将各个子带的Bandelet系数进行组合，得到原图的Bandelet系数。

3.根据权利要求2所述的基于自适应多尺度Bandelet包的SAR图像去噪压缩方法，其中步骤(2.2)的实现步骤如下：

(3.1)将大小为P×P的子带分为大小为m×m的子块；

(3.2)在每个m×m子块中进行方向的离散化，在所有可能的方向上使用折半查找法搜索使得Lagrange函数值最小的方向，该方向即为该子块内几何流的最佳方向；

(3.3)在每个子块中的最佳方向上计算弯曲子波变换，得到其Bandelet系数。

4.根据权利要求3所述的SAR图像去噪压缩方法，其中步骤(3.2)中所述的搜索几何流的最佳方向，是根据Lagrange函数最小化的原则，采用折半查找法搜索使得Lagrange函数值最小的方向，即为该子块中的最佳几何流方向，实现步骤如下：首先初始化(min，max)在(0，π)的范围内，依次执行如下步骤，

(4.1)计算

(4.2)若

当

min←mid，

当

max←mid；

(4.3)若

则max←mid；

(4.4)若则min←mid；

(4.5)若min＜max-2，则转到(4.1)；

否则有

对应的方向即为最佳几何方向。

Images