CN101685158B - 基于隐马尔科夫树模型的sar图像去噪方法 - Google Patents

Abstract

本发明公开一种基于隐马尔科夫树模型的SAR图像去噪方法，它涉及图像处理领域，主要解决图像细节、纹理信息被平滑及同质区域等效视数低的问题。其步骤为：(1)对数变换和Contourlet分解；(2)对Contourlet系数进行HMT建模，并训练；(3)利用估计参数对Contourlet系数进行校正；(4)建立背景隐马尔科夫训练模型，用该估计参数重新校正系数，并进行Contourlet逆变换和指数变换，得到一次去噪图像；(5)对差值图像进行去噪，得到二次去噪图像；(6)对两次去噪图像进行融合，并对融合后的图像进行旋转平移，获得最终去噪图像并输出。本发明较好地保持了SAR图像的细节和纹理信息，减少SAR图像同质区域内斑点噪声，可用于对SAR图像进行去噪处理。

Description

基于隐马尔科夫树模型的SAR图像去噪方法

技术领域

本发明属于图像处理技术领域，特别是一种涉及SAR图像去噪的方法，可用于对SAR图像、自然图像及医学图像的去噪。 

背景技术

合成孔径雷达SAR作为主动雷达，具有不受光照、天气条件影响的特点，可以全天候、全天时对地观测，还可以透过地表和植被获取信息，在农业、林业、地质、环境、水文、海洋、灾害、测绘与军事领域得到了广泛的应用。SAR图像中由于成像散射体散射回波的相干作用造成的斑点噪声的存在不利于图像中场景的自动分析和SAR图像的理解，使得图像的解译工作变得困难，尤其是对SAR图像中的点目标和边缘的影响非常明显。因此斑点噪声的去除对SAR图像的后续处理例如边缘检测是非常重要的。 

近年来，SAR图像相干斑噪声抑制技术飞速发展，可分为成像前多视平滑技术和成像后滤波技术两大类，后处理技术中空域滤波和变换域方法得到了广泛地研究。空域滤波如Lee，Gamma-MAP等对SAR图像进行去噪实现过程简单易于理解，一定程度上可以有效的减弱噪声的影响，但此类滤波器都不同程度的对图像产生了过平滑作用，使图像变模糊，细节信息丢失严重。近年来，变换域的方法也得到了发展，其中基于小波变换的方法成为研究的焦点，如Crouse等人将小波变换与隐马尔科夫模型联系起来，提出了小波域隐马尔科夫模型，开辟了多尺度变换域统计信号处理这一新的研究领域。但对于具有线状奇异的目标函数， 例如图像中的边缘，小波系数不再稀疏，因此小波在SAR图像去噪中并不能够很好地保持图像中的细节信息。 

为了解决小波基函数各向同性带来的问题，Donoho等人提出了多尺度几何分析，其能够有效表示和处理高维奇异函数，并在图像处理得到了广泛的应用，Contourlet便是其中一种分析工具。随着多尺度几何分析的发展，众多学者提出了多尺度几何变换域下的HMT模型和各种改善的HMT模型。Contourlet变换作为新的多尺度信号表示方法，除了具有小波变换的优点外，还具有多方向和各向异性的特点。以Contourlet系数的统计特征为基础，Po和Do提出了Contourlet域用于图像去噪的HMT模型，Contourlet域HMT模型和小波域HMT模型相比，不但可以描述尺度间的相关性，而且可以描述不同方向间系数的相关性。但这种模型只捕获了微小的尺度内的依赖性，用于SAR图像去噪时，一些细节、纹理信息被平滑，同质区域内的噪声去除不够彻底导致同质区域等效视数低，去噪效果并不理想。 

发明内容

本发明的目的在于克服上述已有方法存在的细节、纹理信息被平滑及同质区域等效视数低的缺点，提出一种基于隐马尔可夫树模型的SAR图像去噪方法， 

以达到更彻底的去噪效果，并很好的保持图像的细节和纹理信息。 

本发明的技术方案是：以Contourlet系数邻域间和尺度间都具有很大的相关性为基础，利用隐马尔可夫树模型HMT模型来捕获系数尺度间的依赖性，利用背景隐马尔科夫模型CHMM来捕获系数邻域内的依赖性，把HMT和CHMM结合起来建立Contourlet域改进的HMT统计模型，模型结构如图3所示，用该统计模型对SAR图像进行去噪处理。 具体实现步骤如下： 

(1)对输入SAR图像依次进行对数变换和Contourlet变换，获得不同方向的Contourlet变换系数； 

(2)采用单方向传递HMT模型和多方向传递HMT模型分别对不同方向的Contourlet系数进行尺度间HMT建模； 

(3)用EM算法对建立的HMT模型进行训练，得到均值、方差、状态转移概率和状态概率最优估计参数，并利用该最优估计参数按照贝叶斯最小均方误差准则对Contourlet系数进行萎缩，得到初始无斑系数； 

(4)对初始无斑系数建立尺度内的背景隐马尔科夫模型CHMM，利用EM算法对尺度内的CHMM进行训练，得到尺度内的估计参数，利用尺度内的估计参数对初始无斑系数进行萎缩，得到最终无斑系数； 

(5)对最终无斑系数依次进行Contourlet逆变换和指数变换，得到初步的去噪图像I1； 

(6)采用各向异性扩散方法对原始SAR图像和去噪图像I1的差值图像进行去噪处理，得到二次去噪图像I2； 

(7)用模极大值方法对初步去噪图像I1和二次去噪图像I2进行融合； 

(8)对融合后的图像用CycleSpinning方法进行旋转平移处理，获得最终去噪图像，并输出。 

本发明与现有技术相比具有如下优点： 

1)本发明充分利用了Contourlet系数之间的相关性 

传统的HMT去噪方法，只考虑了Contourlet系数微小的尺度内相关性，没有充分利用系数之间的相关性，导致去噪图像的一些细节、纹 理信息被平滑，图像噪声去除的不彻底，等效视数较低，本发明提出的改进HMT算法克服了系数相关性没有充分利用的缺点，利用了Contourlet系数尺度间、尺度内和方向内的相关性，得到了较好的去噪结果。 

2、本发明提供了更好的相干斑抑制结果 

传统的空域滤波不同程度的对图像产生了过平滑作用，细节信息丢失严重，传统的Contourlet域去噪算法，由于Contourlet变换非平移不变性，会出现划痕即伪吉布斯现象；本发明采用各向异性扩散方法对差值图像进行去噪，补充了去噪中丢失的图像细节信息，使得融合后图像的细节特征也得到了较好的保持；引入了CycleSpinning，克服Contourlet变换非平移不变性带来的伪吉布斯现象，使得同质区域内噪声去除的彻底，用于评价去噪后同质区域内噪声去除多少的等效视数得到提高。 

仿真实验的结果表明，本发明方法更好的保持了图像的细节和纹理信息，同质区域内噪声去除的较彻底。 

附图说明

图1是本发明的主要操作过程示意图； 

图2是本发明Contourlet变换域尺度间HMT模型图； 

图3是本发明改进的Contourlet变换域HMT模型图； 

图4是用本发明与已有方法对SAR1图像去噪效果对比图； 

图5是用本发明与已有方法对SAR2图像去噪效果对比图。 

具体实施方式

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

步骤1：对SAR图像作对数变换，将乘性噪声转换成为加性高斯白 噪声来处理：logy＝log z+log x，其中y表示输入SAR图像，z表示噪声图像，x表示不含噪声的图像，采用方向不相同的两组方向滤波器，该两组方向滤波器的方向分别为4，4，4和4，8，8，对经对数变换后的数据进行Contourlet分解，得到的Contourlet变换系数分别为y₁和y₂。 

步骤2：对变换系数y₁建立单方向传递HMT模型，如图2(a)所示，黑色方形为父节点，四个空方形是其子节点，该模型的参数集为Θ₁；对变换系数y₂建立多方向传递HMT模型，如图2(b)所示，黑色方形为父节点，四个空方形是其子节点，每一个父节点对应的四个子节点分布在两个不同的方向子带内，该模型的参数集为Θ₂。 

步骤3：对包括均值、方差、状态转移概率及根节点的初始状态概率的两个参数集Θ₁和Θ₂进行初始化，再用EM算法对单方向传递HMT模型和多方向传递HMT模型进行训练，得到Θ₁和Θ₂的最优估计参数；利用该最优估计参数按照贝叶斯最小均方误差准则对Contourlet系数进行萎缩，得到初始无斑系数，萎缩公式如下： 

${\hat{c}}_{j,k,i}=\underset{m=\left\{\mathrm{1,2}\right\}}{\mathrm{\Σ}}p(S_{j,k,i}=m)\×\frac{{\mathrm{\σ}}_{j,k,i,m}^2}{{\mathrm{\σ}}_{j,\mathrm{speckle}}^2+{\mathrm{\σ}}_{j,k,i,m}^2}{c}_{s_{j,k,i}}$

其中m表示节点的状态，取值为1或2，p(S_j，k，i＝m)是j尺度，k方向，i位置的节点状态概率，σ_{j，k，i，m} ²是状态为m，j尺度，k方向，i位置的节点的方差，σ_j，speckle ²是尺度j的噪声方差， 

是j尺度，k方向，i位置的Contourlet系数， 

是萎缩后的初始无斑系数。 

步骤4：对初始无斑系数建立尺度内的背景隐马尔科夫模型CHMM，如图3所示，黑色的节点是Contourlet变换系数，白色的圆形节点是Contourlet的系数的隐状态变量，状态数为2，菱形节点为Contourlet的系数的背景变量，一个父节点的状态对应4个子节点的状态，j尺度，k方向，i位置的节点的背景v_j，k，i是由邻域内八个节点的局部能量平均值决定的，如下式： 

$v_{j,k,i}=\left\{\begin{array}{c}1,{\mathrm{\λ}}_{j,k,i}^2>{\mathrm{\δ}}_j^2\\ 0,\mathrm{others}\end{array}\right.$

其中，λ_j，k，i ²为系数c_j，k，i的八领域内Contourlet系数的局部能量平均值，δ_j ²是尺度j内系数的平均能量；利用EM算法对尺度内的CHMM进行训练，先将步骤(3)得到的均值、方差和状态概率最优估计参数作为CHMM训练的初始参数，然后进行训练，得到尺度内的估计参数，利用尺度内的估计参数对初始无斑系数进行萎缩，得到最终无斑系数。 

步骤5：对第(4)步得到的最终无斑系数依次进行Contourlet逆变换和指数变换，得到初步的去噪图像I1。 

步骤6：用各向异性扩散方法对原始SAR图像和去噪图像I1的差值图像进行去噪处理，得到二次去噪图像I2。对差值图像利用梯度算子来区分边缘和噪声，然后对差值图像进行扩散去噪，扩散去噪方程如下： 

$v_{i,j}^{k+1}=v_{i,j}^k+\mathrm{\τ}\underset{\underset{(r,s)\≠\left(\mathrm{0,0}\right)}{r,s=-1}}{\overset{1}{\mathrm{\Σ}}}\frac{g\left({\sqrt{2}}^{1-\left|r\right|-\left|s\right|}\left(\right|v_{i+r,j+s}-v_{i,j}\left|\right)\right)(v_{i+r,j+s}-v_{i,j})}{r^2+s^2}$

其中v_i，j ^k是进行各向异性扩散之前的差值图像像素值，v_i，j ^k+1是进行一 次各向异性扩散后的差值图像像素值，τ是尺度参数，取值范围是：0＜τ＜1/6，r，s是目标点邻域系数坐标相对于目标点的位置，g是扩散率函数，在上式中使用的扩散率函数是Perona-Malik扩散率：g(x)＝1/(1+x²/rr²)，rr是扩散率中的对比参数。 

对于rr的选取，根据多次实验，通过比较发现，rr取50时效果最好。 

步骤7：用模极大值方法对初步去噪图像I1和二次去噪图像I2进行融合，具体步骤如下： 

(7a)对初步去噪图像I1和二次去噪图像I2分别进行Contourlet分解，选择的滤波器方向均为4，4，4，得到初步去噪图像I1和二次去噪图像I2的Contourlet变换系数； 

(7b)对两个去噪图像的Contourlet变换系数的细节项进行模极大值融合得到融合图像的Contourlet系数的细节项，两个去噪图像的Contourlet变换系数的模糊项直接相加得到融合图像Contourlet系数的模糊项； 

(7c)对融合图像的Contourlet系数进行Contourlet反变换得到融合后的图像。 

步骤8：对融合后的图像用CycleSpinning方法进行旋转平移处理，具体步骤如下： 

(8a)通过行和列方向上的每组平移量得到一个不同的去噪结果s_i，j： 

s_i，j＝S_-i，-j(T^-1(Λ[T(S_i，j(x))])) 

其中，S是循环平移算子，T是Contourlet变换算子，T^-1是Contourlet逆变换算子，Λ为Contourlet域HMT去噪算子，下标-i，-j，i和j分别为行和列方向上的平移量； 

(8b)对所有去噪结果进行线性平均，得到抑制伪吉布斯现象的去噪结果： 

$s_t=\frac{1}{K_1{K}_2}\underset{i=0,j=0}{\overset{K1,K2}{\mathrm{\Σ}}}{s}_{i,j}$

其中，K₁，K₂分别表示行和列方向上的最大平移量，均为4。 

步骤9：把步骤(8)中得到的去噪结果s_t作为最终去噪结果，并输出。 

以下给出仿真实验结果分析，以进一步说明本发明的效果： 

1、仿真条件 

本发明中的两幅测试SAR图像均是从美国Sandia Labs网站上获得的Horse track SAR图像上截取的，大小均为256×256，分别命名为SAR1和SAR2。采用Gamma_MAP滤波方法、Lee滤波方法、小波域HMT去噪方法、Contourlet域HMT去噪方法和本发明方法1改进单方向传递HMT模型及本发明方法2改进多方向传递HMT模型对SAR1和SAR2图像进行去噪处理。 

2、仿真结果分析 

仿真的结果如图4和图5，其中： 

图4(a)为SAR1原图像， 

图4(b)为SAR1图像经Gamma_MAP滤波后的去噪图像， 

图4(c)为SAR1图像经Lee滤波后的去噪图像， 

图4(d)为SAR1图像经小波域HMT去噪后的去噪图像， 

图4(e)为SAR1图像经Contourlet域HMT去噪后的去噪图像， 

图4(f)为SAR1图像经本发明改进单方向传递HMT模型方法去噪后的去噪图像， 

图4(g)为SAR1图像经本发明改进多方向传递HMT模型方法去噪后的去噪图像。 

图5(a)为SAR2原图像， 

图5(b)为SAR2图像经Gamma MAP滤波后的去噪图像， 

图5(c)为SAR2图像经Lee滤波后的去噪图像， 

图5(d)为SAR2图像经小波域HMT去噪后的去噪图像， 

图5(e)为SAR2图像经Contourlet域HMT去噪后的去噪图像， 

图5(f)为SAR2图像经本发明改进单方向传递HMT模型方法去噪后的去噪图像， 

图5(g)为SAR2图像经本发明改进多方向传递HMT模型方法去噪后的去噪图像。 

从图4(b)和图4(c)可以看出，Gamma MAP滤波和Lee滤波对SAR1图像的斑点噪声去除的较好，但SAR1图像的边缘被模糊，图像细节信息丢失较严重；由图4(d)和图4(e)，从同质区域内噪声的去除可以看出Contourlet域HMT方法要优于小波域HMT方法，但这两种方法对SAR1图像去噪后的同质区域内噪声去除的都不够彻底。从图4(f)和图4(g)，可以看出SAR1去噪后同质区域内噪声去除的很彻底，细节信息、边缘信息及纹理信息得到了很好的保持。 

从图5可以看出，本发明在彻底去除噪声的同时更多的保持了SAR2图像的边缘和纹理细节；可见，本发明在细节，纹理保持方面有相当大的优势。 

对SAR图像去噪性能的评价指标包括等效视数ENL、均值及均值比，这些评价标准的定义如下： 

(1)等效视数(ENL)： $\mathrm{ENL}=A\frac{u^2}{{\mathrm{\σ}}^2}$

其中u和σ分别为相干斑抑制后图像中某同质区域的均值和标准 差，对于强度图像A＝1，对于幅度图像A＝4/π-1。它是衡量图像中相干斑相对强度的一种指标，相干斑抑制程度越深，等效视数越大。对于SAR1图像选择尺寸为40×40，50×50，40×40的同质区域作为计算ENL的测试数据，如图4(a)中的方形框标示区域；对于SAR2图像选择尺寸为30×30，30×30，30×30的同质区域作为计算ENL的测试数据，如图5(a)中的方形框标示区域。 

(2)图像的均值：均值反映图像的平均亮度，相干斑抑制前后的图像均值要基本保持一致。 

(3)均值比：原图像I与相干斑抑制后图像D均值之比， 

$\mathrm{\γ}=\frac{\mathrm{mean}\left(D\right)}{\mathrm{mean}\left(I\right)}$

理想情况下比率γ为纯斑点噪声的统计幅值1。实际的γ值与1相差越大，表明辐射度失真越厉害。 

SAR1和SAR2图像去噪后的评价指标如表1和表2所示： 

表1 对SAR1图像进行去噪后去噪图像的评价指标比较 

SAR1图像	区域1   ENL	区域2   ENL	区域3   ENL	均值	均值比
						Gamma_MAP	217.1888	122.2114	107.9559	83.1247	0.9913
Lee滤波	217.1888	122.2114	107.9559	83.6654	0.9978
						小波HMT	78.8185	58.1684	36.5429	80.3143	0.9577
Contourlet HMT	120.2160	80.0639	52.2560	79.9457	0.9533
						本发明方法1	327.6217	158.3174	140.7742	83.7811	0.9992
本发明方法2	338.5326	153.0471	143.3387	83.7741	0.9991

[0084] 表2 对SAR2图像进行去噪后去噪图像的评价指标比较 

SAR2图像	区域1   ENL	区域2   ENL	区域3   ENL	均值	均值比
						Gamma_MAP	66.5108	222.8523	69.6148	72.7787	0.9805
Lee滤波	66.5108	222.8523	69.6148	73.7958	0.9942
						小波HMT	35.2699	71.6358	26.7679	71.1371	0.9584
Contourlet HMT	39.7765	89.1000	30.8420	70.7065	0.9526
						本发明方法1	83.4655	360.5115	71.6091	74.1435	0.9989
本发明方法2	82.3402	334.9756	73.9066	74.1274	0.9987

从表1及表2可以看出： 

(a)本发明的等效视数值均比其它几种经典方法的等效视数值大，表明本发明对SAR图像相干斑噪声抑制程度较深，斑点噪声去除的较彻底，去噪效果较好； 

(b)本发明的均值和相应原图像的均值基本相同，表明本发明去噪后图像的平均亮度基本保持不变； 

(c)本发明的均值比与1的差距比其它经典方法小，表明本发明去噪过程中辐射度失真很小，细节信息和纹理信息保持的较好。 

综上，本发明方法既较好地保持了SAR图像中的点目标及其边缘等细节特征和纹理信息，也大幅度的减少了SAR图像同质区域内的斑点噪声。 

Claims (6)

1.一种基于隐马尔科夫树模型HMT的SAR图像去噪方法，包括如下步骤：

(1)对输入SAR图像依次进行对数变换和Contourlet变换，获得不同方向的Contourlet变换系数；

(2)采用单方向传递模型和多方向传递模型分别对不同方向的Contourlet系数进行尺度间HMT建模；

(3)用EM算法对建立的HMT模型进行训练，得到均值、方差、状态转移概率和状态概率最优估计参数，并利用该最优估计参数按照贝叶斯最小均方误差准则对Contourlet系数进行萎缩，得到初始无斑系数；

(4)对初始无斑系数建立尺度内的背景隐马尔科夫模型CHMM，利用EM算法对尺度内的CHMM进行训练，得到尺度内的估计参数，利用尺度内的估计参数对初始无斑系数进行萎缩，得到最终无斑系数；

(5)对最终无斑系数依次进行Contourlet逆变换和指数变换，得到初步的去噪图像I1；

(6)采用各向异性扩散方法对原始SAR图像和去噪图像I1的差值图像进行去噪处理，得到二次去噪图像I2；

(7)用模极大值方法对初步去噪图像I1和二次去噪图像I2进行融合；

(8)对融合后的图像用CycleSpinning方法进行旋转平移处理，获得最终去噪图像，并输出。

2.根据权利要求1所述的的SAR图像去噪方法，其中步骤(1)中所述的Contourlet变换，采用方向不相同的两组方向滤波器，该两组方向滤波器的方向分别为4，4，4和4，8，8，得到的Contourlet变换系数分别为y₁和y₂。

3.根据权利要求1所述的SAR图像去噪方法，其中步骤(2)所述的采用单方向传递模型和多方向传递模型分别对不同方向的Contourlet系数进行尺度间HMT建模，是对变换系数y₁建立单方向传递模型，该模型的参数集为Θ₁；对变换系数y₂建立多方向传递模型，该模型的参数集为Θ₂，这两种参数集Θ₁和Θ₂内的参数包括均值、方差、状态转移概率及根节点的初始状态概率。

4.根据权利要求1或3所述的SAR图像去噪方法，其中步骤(3)中所述的用EM算法对建立的HMT模型进行训练，是先对包括均值和方差，状态转移概率，根节点的初始状态概率的两个参数集Θ₁和Θ₂进行初始化，再用EM算法对HMT模型进行训练，得到Θ₁和Θ₂的最优估计参数。

5.根据权利要求1所述的SAR图像去噪方法，其中步骤(4)中所述的利用EM算法对尺度内的CHMM进行训练，是将步骤(3)得到的均值、方差、状态转移概率和状态概率最优估计参数作为CHMM训练的初始参数。

6.根据权利要求1所述的SAR图像去噪方法，其中步骤(8)中所述的对融合后的图像用CycleSpinning方法进行旋转平移处理，按如下步骤进行：

(6a)通过行和列方向上的每组平移量得到一个不同的去噪结果s_i，j：

s_i，j＝S_-i，-j(T^-1(Λ[T(S_i，j(x))]))

其中，S是循环平移算子，T是Contourlet变换算子，T^-1是Contourlet逆变换算子，Λ为Contourlet域HMT去噪算子，下标-i，-j，i和j分别为行和列方向上的平移量；

(6b)对所有去噪结果进行线性平均，得到抑制伪吉布斯现象的去噪结果S_t：

$s_t=\frac{1}{K_1{K}_2}\underset{i=0,j=0}{\overset{K1,K2}{\mathrm{\Σ}}}{s}_{i,j}$

其中，K₁，K₂分别表示行和列方向上的最大平移量。

Images