CN101639934B

CN101639934B - 基于轮廓波域块隐马尔可夫模型sar图像去噪方法

Info

Publication number: CN101639934B
Application number: CN2009100237887A
Authority: CN
Inventors: 焦李成; 侯彪; 田福苓; 王爽; 张向荣; 马文萍
Original assignee: Xidian University
Current assignee: Xidian University
Priority date: 2009-09-04
Filing date: 2009-09-04
Publication date: 2012-07-04
Anticipated expiration: 2029-09-04
Also published as: CN101639934A

Abstract

本发明公开了一种基于轮廓波域块隐马尔可夫模型SAR图像去噪方法，它涉及图像处理领域，主要解决现有方法缺少空间适应性，细节信息丢失严重，边缘模糊的问题。其步骤为：(1)对数变换和Contourlet分解；(2)对Contourlet系数进行Block HMM建模并训练；(3)利用估计参数对Contourlet系数进行校正；(4)对校正后的Contourlet数依次进行Contourlet逆变换和反对数变换，得到一次去噪图像；(5)对差值图像进行去噪，得到二次去噪图像；(6)对两次去噪图像进行融合，并对融合后的图像进行旋转平移，获得最终去噪图像并输出。本发明用于SAR图像去噪时改善了空间适应性，降低了计算复杂度，提高了去噪后图像边缘的清晰度。

Description

基于轮廓波域块隐马尔可夫模型SAR图像去噪方法

技术领域

本发明属于图像处理技术领域，尤其涉及SAR图像斑点噪声抑制的方法，该方法可用于对SAR图像、医学图像和自然图像进行去噪处理。

背景技术

合成孔径雷达SAR具有全天候、全天时对地观测的能力，还可以透过地表和植被获取信息，在军事、遥感等领域发挥出越来越大的作用。SAR图像中斑点噪声的存在，严重影响了后续各种自动解译和场景理解任务完成的准确性程度。因此斑点噪声的去除对SAR图像的后续处理是非常重要的。

斑点噪声的去除通常包括空域滤波和变换域方法。空域滤波如Gamma-MAP对SAR图像进行去噪，在一定程度上可以有效的减弱噪声的影响，但此类滤波器都不同程度的对SAR图像产生了过平滑作用，使图像变模糊，细节信息丢失严重。变换域中的阈值滤波，包括软、硬阈值滤波，作用于整个图像或整个尺度将损失其中的大部分细节，信息丢失也很严重。Crouse等人将小波变换与隐马尔可夫模型联系起来，提出了小波域隐马尔可夫模型，开辟了变换域统计信号处理这一新的研究领域。随着多尺度几何分析的发展，众多学者提出了多尺度几何变换域下的HMT模型和各种改善的HMT模型。Po和Do提出Contourlet域用于图像去噪的HMT模型。这种模型可以应用到SAR图像去噪中，但这种模型只捕获了微小的尺度内的依赖性，对SAR图像去噪时，同质区域内的噪声去除不够彻底，去噪效果并不理想。

Po和Do对Contourlet系数的尺度内，尺度间和方向间的相关性进行了分析，并得出结论：尺度内Contourlet系数的相关性即聚集性比尺度间的相关性即延续性对系数的影响显著。因此把背景隐马尔可夫模型CHMM推广到Contourlet变换域，这种新的模型称之为Contourlet域背景隐马尔可夫模型CD-CHMM。然而应用CD-CHMM模型对SAR图像去噪缺少空间适应性，造成了边缘模糊现象。

发明内容

本发明的目的在于克服上述已有方法的不足，提出一种基于轮廓波域块隐马尔可夫模型的SAR图像去噪方法，以彻底去除同质区域内噪声，保持细节信息，改善空间适应性和提高边缘清晰度。

本发明的技术方案是：对小波域块隐马尔可夫模型Block HMM进行背景以及块的改进，把这个改进的统计模型推广到Contourlet域，建立了Contourlet域Block HMM模型，并用于SAR图像去噪。具体实现步骤包括：

A建模步骤

(A1)对输入SAR图像依次进行对数变换和Contourlet变换，获得Contourlet变换系数，得到高频Contourlet系数；

(A2)对得到的高频Contourlet系数，按如下步骤建立Block HMM模型：

①定义改进的背景为：

context = w_{0} \cdot Σ_{i = 1}^{4} | {NA}_{i} | + w_{1} \cdot Σ_{i = 1}^{4} | {NB}_{i} | + w_{2} \cdot | P | + w_{3} \cdot (| C_{1} | + | C_{2} |)

其中NA是水平和垂直的相邻系数，NB是对角线相邻系数，C₁，C₂是最近的两个兄弟节点，P是高频Contourlet系数的父节点，w₀＝1，w₁＝0.6，w₂＝0.2，w₃＝0.6；

②利用上述公式计算每一个高频Contourlet系数的改进背景；

③利用下面的公式对每一个高频Contourlet系数的改进背景进行二值化，得到二值化值v：

v = \{\begin{matrix} 1, context < 4 w_{0} \cdot M + 4 w_{1} \cdot M + w_{2} \cdot Mp + w_{3} \cdot (M c_{1} + M c_{2}) \\ 0, context &GreaterEqual; 4 w_{0} \cdot M + 4 w_{1} \cdot M + w_{2} \cdot Mp + w_{3} \cdot (M c_{1} + {Mc}_{2}) \end{matrix}

其中，w₀＝1，w₁＝0.6，w₂＝0.2，w₃＝0.6，v是改进背景二值化值，M，Mp，Mc₁，Mc₂分别表示高频Contourlet系数子节点所在子带的绝对平均值，父子带的绝对平均值，两个同尺度相邻方向子带的绝对平均值；

④将高频Contourlet系数分成4×4或8×8或16×16的多个小块，再根据高频Contourlet系数的改进背景的二值化值把一个小块划分为两个互不相关的块；

B模型训练步骤

(B1)初始化Block HMM模型的均值、方差和初始状态概率，并通过EM算法训练Block HMM模型，得到该模型的最优估计参数；

(B2)利用得到的最优估计参数按照贝叶斯最小均方误差准则，对高频Contourlet系数进行降斑处理，得到高频Contourlet无斑系数；

(B3)对得到的高频Contourlet无斑系数依次进行Contourlet逆变换和反对数变换，得到一次去噪图像I1；

C后处理步骤

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

(C2)用模极大值方法对一次去噪图像I1和二次去噪图像I2进行融合处理；

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

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

1)本发明更加充分地捕获了Contourlet系数之间的相关性

小波域Block HMM的背景是传统的由八邻域系数决定的背景，只反映了尺度内的相关性；传统的隐马尔可夫树模型用于SAR图像去噪时，只考虑了尺度间的相关性，导致同质区域内噪声去除的不彻底；传统的背景隐马尔可夫模型用于SAR图像去噪时，只考虑了尺度内的相关性，导致细节信息被平滑；本发明中改进的背景由尺度间的父节点、方向间的兄弟节点和尺度内的邻域节点共同决定，充分捕捉了Contourlet系数尺度间、方向间和尺度内的相关性，用于SAR图像去噪时较好的保持了SAR图像中的细节信息，同质区域内噪声去除的较彻底。

2)本发明提供了更好的空间适应性

小波域Block HMM的高频系数通过背景信息分为两个块，用于SAR图像去噪时，空间适应性差，计算复杂度较大；而本发明中把方向尺度内的高频Contourlet系数通过改进的背景分成较多的块，达到在可靠相关Contourlet系数数目和空间适应性的一个平衡，克服了SAR图像去噪缺少空间适应性、边缘模糊及计算复杂度大等缺点

仿真实验的结果表明，本发明方法更好的保持了图像的细节信息，改善了空间适应性，同质区域内噪声去除的较彻底，边缘比较清晰。

附图说明

图1是本发明的图像去噪流程图；

图2是本发明的Contourlet域块隐马尔可夫模型图；

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

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

具体实施方式

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

步骤1：对输入SAR图像作对数变换，将乘性噪声转换成为加性高斯白噪声来处理：log y＝log z+log x，其中y表示输入SAR图像，z表示噪声图像，x表示不含噪声的图像；对数变换后的图像进行Contourlet分解，其中方向滤波器的方向为4，4，4，获得Contourlet变换系数，得到高频Contourlet系数。

步骤2：对得到的高频Contourlet系数建立Contourlet域改进的BlockHMM模型如图2所示。

具体实施如下：

(2a)由尺度间的父节点、方向间的兄弟节点和尺度内的邻域节点共同决定改进的背景，其计算公式如下：

context = w_{0} \cdot Σ_{i = 1}^{4} | {NA}_{i} | + w_{1} \cdot Σ_{i = 1}^{4} | {NB}_{i} | + w_{2} \cdot | P | + w_{3} \cdot (| C_{1} | + | C_{2} |)

利用上面公式计算每一个高频Contourlet系数的改进背景，并利用下面的二值化公式对计算得到的高频Contourlet系数的改进背景进行二值化处理，得到高频Contourlet系数的改进背景二值化值，

v = \{\begin{matrix} 1, context < 4 w_{0} \cdot M + 4 w_{1} \cdot M + w_{2} \cdot Mp + w_{3} \cdot (M c_{1} + M c_{2}) \\ 0, context &GreaterEqual; 4 w_{0} \cdot M + 4 w_{1} \cdot M + w_{2} \cdot Mp + w_{3} \cdot (M c_{1} + {Mc}_{2}) \end{matrix}

其中，context是高频Contourlet系数的改进背景，NA是水平和垂直的相邻系数，NB是对角线相邻系数，C₁，C₂是最近的两个兄弟节点，P是高频Contourlet系数的父节点，w₀＝1，w₁＝0.6，w₂＝0.2，w₃＝0.6，v是改进背景二值化值，M，Mp，Mc₁，Mc₂分别表示高频Contourlet系数子节点所在子带的绝对平均值，父子带的绝对平均值，两个同尺度相邻方向子带的绝对平均值，图2中的Ω₁，Ω₂，Ω₃是上下文向量；

(2b)把高频Contourlet系数分成4×4、8×8或16×16的多个小块，再利用改进的背景二值化值对每一个小块中的高频Contourlet系数进行块标签，标签值为1或0，有相同块标签值的高频Contourlet系数就认为属于同一个块，把一个小块划分为两个互不相关的块；

(2c)对上面得到的块进行修订，把非常小的块删除掉，得到修订好的块。

步骤3：对Block HMM模型的均值、方差和初始状态概率进行初始化处理，其中均值均设为0，在同一个块中，因为背景都是相同的，不用再考虑背景信息，所以对每一个块的训练可以看成是一个高斯混合场GMF的训练，用EM算法进行训练，得到高频Contourlet系数的最优估计参数，训练过程如下：

(3a)初始化参数集∏_k，b ^p，并利用贝叶斯公式计算高频Contourlet系数的状态转移概率：

P (S_{k, b} = m | c_{i, k}, c_{i, k} = v, b_{i, k} = b, Π_{k, b}^{p}) = \frac{P_{S_{k, b}} (m) g (c_{i, k}; 0, σ_{k, b, m}^{2})}{Σ_{m = 0}^{1} P_{S_{k, b}} (m) g (c_{i, k}; 0, σ_{k, b, m}^{2})}

(3b)通过下面的公式计算高频Contourlet系数的参数集∏_k，b ^p+1：

P_{S_{k, b}} (m) = \frac{Σ_{| (i, k) | B_{i, k} = b |} P (S_{k, b} = m | c_{i, k}, v_{i, k} = v, B_{i, k} = b, Π_{k, b}^{p})}{N}

σ_{k, b, m}^{2} = \frac{Σ_{| (i, k) | B_{i, k} = b |} c_{i, k}^{2} P (S_{k, b} = m | c_{i, k}, v_{i, k} = v, B_{i, k} = b, Π_{k, b}^{p})}{Σ_{| (i, k) | B_{i, k} = b |} P (S_{k, b} = m | c_{i, k}, v_{i, k} = v, B_{i, k} = b, Π_{k, b}^{p}) P_{S_{k, b}} (m)}

其中，高频Contourlet系数的参数集∏_k，b ^p包括均值，方差和状态概率，m表示节点的状态，k是高频Contourlet系数子带的标签，i是高频Contourlet系数的位置，B_i，k是i位置，k子带块的标签，b是块标签的值，m是高频Contourlet系数节点的状态，N是在块中Contourlet系数的个数，p是循环次数，g(c_i，k；0，σ_k，b，m ²)是高斯函数，c_i，k是i位置，k子带高频Contourlet系数，v_i，k是i位置，k子带Contourlet系数的背景，v是背景值，

是k子带，块标签的值为b的节点状态概率，σ_k，b，m ²是k子带，块标签的值为b，节点状态为m的方差；

更新迭代：p＝p+1，直到收敛为止。

步骤4：利用得到的高频Contourlet系数的最优估计参数，按照贝叶斯最小均方误差准则，对高频Contourlet系数进行降斑处理，得到高频Contourlet无斑系数，高频Contourlet系数降斑公式如下：

{\hat{c}}_{i, k} = \underset{m = {1,2}}{Σ} p (S_{k, b} = m) \times \frac{σ_{k, b, m}^{2}}{σ_{k, speckle}^{2} + σ_{k, b, m}^{2}} c_{i, k}

其中m表示节点的状态，取值为1或2，p(S_k，b＝m)是k子带，块标签的值为b的节点状态概率，

是k子带，块标签的值为b，节点状态为m的方差，σ_k，speckle ²是k子带的噪声方差，c_i，k是i位置，k子带高频Contourlet系数，

是降斑后的高频无斑Contourlet系数。

步骤5：对低频Contourlet系数和得到的高频Contourlet无斑系数相结合，得到处理后的Contourlet无斑系数，对Contourlet无斑系数依次进行Contourlet逆变换和指数变换，得到一次去噪图像I1。

步骤6：用各向异性扩散方法对输入SAR图像和一次去噪图像I1的差值图像进行去噪处理，得到二次去噪图像I2。对差值图像的扩散去噪方程如下：

v_{i, j}^{k + 1} = v_{i, j}^{k} + τ Σ_{\underset{(r, s) &NotEqual; (0,0)}{r, s = - 1}}^{1} \frac{g ({\sqrt{2}}^{1 - | r | - | s |} (| v_{i + r, j + s} - v_{i, j} |)) (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是扩散率中的对比参数。

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

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

(7b)对一次去噪图像I1和二次去噪图像I2的高频Contourlet系数进行模极大值融合得到融合图像的高频Contourlet系数，两个去噪图像的低频Contourlet系数直接相加得到融合图像的低频Contourlet系数；

(7c)对融合图像的高频Contourlet系数和融合图像的低频Contourlet系数进行组合，得到融合图像的Contourlet系数，对融合图像的Contourlet系数进行Contourlet反变换得到融合后的图像。

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

(8a)通过对高频Contourlet系数在行和列方向上平移处理，得到一组不同的去噪结果I_i，j：

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

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

(8b)对所有去噪结果进行线性平均：

S_{e} = \frac{1}{K_{1} K_{2}} Σ_{i = 0, j = 0}^{K_{1}, K_{2}} I_{i, j}

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

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

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

1、仿真条件

本发明中的两幅测试SAR图像分别为一幅x波段的SARfield图像和一幅Ku波段的Horsetrack图像。采用Gamma_MAP滤波方法、Contourlet域阈值去噪方法、Contourlet域隐马尔可夫树HMT模型去噪方法、Contourlet域背景隐马尔可夫模型CHMM去噪方法和本发明方法对SARfield图像和Horsetrack图像进行去噪处理。

2、仿真结果分析

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

图3(a)为SARfield原图像，

图3(b)为SARfield图像经Gamma_MAP滤波后的去噪图像，

图3(c)为SARfield图像经轮廓波Contourlet域阈值去噪后的去噪图像，

图3(d)为SARfield图像经轮廓波Contourlet域隐马尔可夫树HMT模型去噪后的去噪图像，

图3(e)为SARfield图像经轮廓波Contourlet域背景隐马尔可夫模型CHMM去噪后的去噪图像，

图3(f)为SARfield图像经本发明方法去噪后的去噪图像。

图4(a)为Horsetrack原图像，

图4(b)为Horsetrack图像经G-MAP滤波后的去噪图像，

图4(c)为Horsetrack图像经轮廓波Contourlet域阈值去噪后的去噪图像，

图4(d)为Horsetrack图像经轮廓波Contourlet域隐马尔可夫树HMT去噪后的去噪图像，

图4(e)为Horsetrack图像经轮廓波Contourlet域背景隐马尔可夫模型CHMM去噪后的去噪图像，

图4(f)为Horsetrack图像经本发明方法去噪后的去噪图像。

从图3(b)可以看出，Gamma_MAP滤波对SARfield图像的斑点噪声进行了较好的平滑，但边界模糊，图像细节信息丢失严重；从图3(c)可以看出，由于过平滑的原因使得图像细节信息丢失严重，说明噪声滤除能力有限；从图3(d)可以看出，Contourlet域HMT模型去噪方法能较好的保持细节信息，但SARfield图像去噪后同质区域内噪声去除的不够彻底；从图3(e)可以看出，基于Contourlet域CHMM方法处理后的SARfield图像，其同质区域内噪声去除不够彻底，边缘模糊；从图3(f)可以看出，SARfield图像去噪后空间适应性好，同质区域内噪声去除彻底，细节信息保持的较好，边缘较清晰。

从图4可以看出，本发明在提高边缘的清晰度的同时较好的保持了Horsetrack图像的细节信息、改善了空间适应性；可见，本发明在细节保持，提高边缘清晰度和改善空间适应性方面有相当大的优势。

SARfield和Horsetrack图像去噪后的评价指标包括等效视数ENL、均值和均值比如表1和表2所示：

表1对SARfield图像进行去噪的去噪结果指标比较

SARfield	区域1ENL	区域2ENL	区域3ENL	均值	均值比
						Gamma_MAP滤波	28.1532	19.2887	38.4854	105.5474	0.9876
Contourlet域阈值	18.5148	12.7519	10.7738	101.7181	0.9518
						Contourlet域HMT	12.9199	9.8039	10.5558	102.3756	0.9580
Contourlet域CHMM	23.1863	16.4171	18.6424	100.4813	0.9402
						本发明方法	58.0173	33.3215	58.1845	106.6627	0.9981

表2对Horsetrack进行去噪的去噪结果指标比较

Horsetrack	区域1ENL	区域2ENL	区域3ENL	均值	均值比
						Gamma_MAP滤波	66.5108	222.8523	69.6148	72.7787	0.9805
Contourlet域阈值	69.8477	218.1453	46.3711	70.2888	0.9470
						Contourlet域HMT	39.2305	85.2133	30.4376	70.6914	0.9524
Contourlet域CHMM	64.0009	190.3675	50.4084	68.4158	0.9218
						本发明方法	91.2217	432.2075	101.9164	74.1334	0.9988

从表1可以看出，本发明的等效视数值均比其它几种经典方法的等效视数值大，表明本发明对SARfield图像相干斑噪声抑制程度较深，斑点噪声去除的较彻底，去噪效果较好；本发明的均值和SARfield原图像的均值基本相同，表明本发明对SARfield图像去噪后，图像的平均亮度基本保持不变；本发明的均值比与1的差距比其它经典方法小，表明本发明对SARfield图像的去噪过程中，辐射度失真很小，空间适应性好，细节和边缘保持的好。

从表2可以看出，本发明对Horsetrack图像去噪后，在等效视数、均值和均值比上都比其它方法好，表明本方法对Horsetrack图像相干斑抑制程度深，去噪过程中辐射度失真较小，去噪效果较好。

综上，本发明方法既较好地改善了空间适应性，降低了计算复杂度，提高了去噪后图像的边缘清晰度，也有效地保持了图像的细节信息和彻底的去除了同质区域内的斑点噪声。