CN103810680A - 基于方向波域hmt模型的sar图像降斑方法 - Google Patents

Abstract

本发明公开了一种基于方向波域HMT模型的SAR图像降斑方法，主要解决现有方法中不能有效地保持图像边缘结构等细节信息及不能够较好地保持SAR图像的辐射特性的问题。其实现步骤是：(1)输入图像；(2)建立非对数加性模型；(3)构造方向波域HMT模型参数集；(4)估计方向波域HMT模型参数集；(5)估计无斑图像的方向波高频子带系数；(6)方向波逆变换；(7)输出估计的无斑图像。本发明具有有效保持图像的边缘和细节信息、降斑效果更明显的优点，可用于SAR图像的降斑处理。

Description

基于方向波域HMT模型的SAR图像降斑方法

技术领域

本发明属于图像处理技术领域，更进一步涉及合成孔径雷达图像降斑技术领域中的一种基于方向波域隐马尔可夫树(Hidden Markov tree，HMT)模型的合成孔径雷达(Synthetic Aperture Radar，SAR)图像降斑方法。本发明可用于对SAR图像进行降斑处理。

背景技术

SAR图像中相干斑噪声的存在，严重影响了人类视觉对SAR图像的感知，制约了SAR图像在图像分割及变化检测方面的应用，因此对SAR图像进行有效的去噪是SAR图像处理中的一个迫切需求。

西安电子科技大学在其专利申请“基于SWT域改进粒子滤波的SAR图像降斑方法”（专利申请号：CN201110060825.9，公开号：CN102129672A）中提出了一种基于SWT域改进粒子滤波的SAR图像降斑方法。该方法将待降斑的空域图像变换到平稳小波域，并提取不同变换尺度水平、垂直和对角方向下的SWT域图像，对其应用改进粒子滤波降斑方法进行降斑处理，然后将得到的降斑后小波域图像进行逆变换得到最终的降斑图像。该方法虽然能够解决已有方法降斑后图像边缘和纹理模糊等问题，但仍然存在的不足是，平稳小波变换只能描述点奇异性，却不能有效地刻画图像中的二维线奇异性，会损失部分细节信息。

西安电子科技大学在其专利申请“基于隐马尔科夫树模型的sar图像去噪方法”（专利申请号：CN200910023171，公开号：CN101685158A）中提出了一种基于隐马尔科夫树模型的SAR图像去噪方法。该方法具体步骤是：首先，利用对数变换将SAR图像相干斑噪声转变为加性噪声，然后利用HMT模型和背景HMT模型对含噪SAR图像的轮廓波系数进行校正，得到第一幅去噪图像，之后利用各向异性扩散方法得到第二幅去噪图像，最后对两幅去噪图像进行融合得到最终去噪图像。该方法虽然利用了轮廓波系数尺度间、尺度内和方向内的相关性，克服轮廓波变换非平移不变性带来的伪吉布斯现象，提高了同质区域的等效视数，但是仍然存在的不足是，通过对数运算将相干斑噪声转换成加性噪声，不能够较好地保持SAR图像的辐射特性，导致后续处理无法将相干斑噪声彻底去除。

发明内容

本发明的目的在于针对上述现有技术的不足，提出一种基于方向波域HMT模型的SAR图像降斑方法。本发明能够在滤除相干斑噪声的同时，很好地保持图像的边缘等结构信息。

本发明实现的具体步骤包括如下：

1.一种基于方向波域HMT模型的SAR图像降斑方法，包括如下步骤：

(1)输入图像。

任选一幅SAR图像作为待降斑图像。

(2)建立非对数加性模型。

按照下式，对待降斑图像建立非对数加性模型，将待降斑图像转化为地物的真实后向散射强度和加性噪声的和。

I=RX=X+(R-1)X=X+N

其中，I表示待降斑图像，R表示相干斑噪声，R的均值为1，X表示地物的真实后向散射强度，N表示将要滤除的加性噪声。

(3)构造方向波域HMT模型参数集。

(3a)对待降斑图像进行尺度为2的方向波变换，得到待降斑图像的方向波低频系数子带和3个方向上的方向波高频系数子带。

(3b)利用待降斑图像的3个方向上的方向波高频系数子带，分别建立3个方向上的方向波域HMT模型。

(3c)利用待降斑图像的3个方向上的方向波高频系数子带，分别构造3个方向上的方向波域HMT模型参数集。

(4)估计方向波域HMT模型参数集。

采用期望最大化方法，分别对3个方向上的方向波域HMT模型参数集进行训练，得到3个方向上的方向波域HMT模型参数集的最优估计值。

(5)估计无斑图像的方向波高频子带系数。

利用3个方向上的方向波域HMT模型参数集的最优估计值，按照贝叶斯最小均方误差准则，估计无斑图像的方向波高频子带系数。

(6)方向波逆变换。

对待降斑图像的方向波低频子带系数和估计的无斑图像所有方向波高频子带系数，进行方向波逆变换，得到估计的无斑图像。

(7)输出估计的无斑图像。

本发明方法与现有的技术相比具有以下优点：

第一，本发明采用方向波变换对待降斑图像进行分解，克服了现有技术中采用平稳小波变换进行分解导致的不能有效地保持图像边缘结构等细节信息的缺点，方向波变换具有多尺度、多方向的特性，能够更好地稀疏表示具有直线和曲线奇异性的边缘，使得本发明具有有效保持图像的边缘和细节信息的优点。

第二，本发明采用非对数加性模型将SAR图像的相干斑噪声转化为加性噪声，克服了现有技术中采用对数变换方法不能够较好地保持SAR图像的辐射特性的缺点，使得本发明具有对SAR图像降斑效果更明显的优点。

附图说明

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

图2是本发明及现有技术对合成孔径雷达SAR图像field的降斑效果对比图；

图3是本发明及现有技术对合成孔径雷达SAR图像river的降斑效果对比图。

具体实施方式

下面结合附图对本发明做进一步的描述。

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

步骤1，输入图像。

任选一幅SAR图像作为待降斑图像。本发明实施例中使用的合成孔径雷达SAR图像分别如图2(a)及图3(a)所示。其中，图2(a)为合成孔径雷达SAR图像田野图filed，大小为256×256，图3(a)为合成孔径雷达SAR图像河流图river，大小为256×256。

步骤2，建立非对数加性模型。

按照下式，对待降斑图像建立非对数加性模型，将待降斑图像转化为地物的真实后向散射强度和加性噪声的和：

I=RX=X+(R-1)X=X+N

其中，I表示待降斑图像，R表示相干斑噪声，R的均值为1，X表示地物的真实后向散射强度，N表示将要滤除的加性噪声。

步骤3，建立方向波域HMT模型，构造方向波域HMT模型参数集。

对待降斑图像进行尺度为2的方向波变换，得到待降斑图像的方向波低频系数子带和3个方向上的方向波高频系数子带。方向波变换的具体步骤如下：

第一步，从方向集{0°，30°，45°，60°，90°，120°，135°，150°}中，任意选择两个不平行的方向，分别作为方向波变换的变换方向和队列方向，构造采样矩阵：

$M=\left[\begin{array}{c}{d}_1\\ d_2\end{array}\right]$

其中，M表示采样矩阵，d₁表示沿变换方向的整数行向量，d₂表示沿队列方向的整数行向量。

第二步，将整数行向量d₁和d₂组成的网格在待降斑图像上进行L次平移，每次平移过程中，保留待降斑图像中和网格结点重合的像素点，得到待降斑图像的L个陪集，其中L表示采样矩阵M的行列式的绝对值。

第三步，沿变换方向和队列方向对待降斑图像的每个陪集分别进行2次一维小波变换，得到待降斑图像的方向波低频系数子带和3个方向上不同尺度的方向波高频系数子带。

利用待降斑图像的3个方向上的方向波高频系数子带，分别建立3个方向上的方向波域HMT模型。建立方向波域HMT模型的具体步骤如下：

第一步，在3个方向上的方向波高频系数子带中，任选一个方向上的尺度为2的方向波高频子带H，从方向波高频系数子带H中任选一个坐标为(i,j)的高频系数C，作为方向波域HMT模型中的一个父节点，其中，i表示高频系数的行序号，j表示高频系数的列序号。

第二步，在3个方向上的方向波高频系数子带中，选取和方向波高频子带H方向相同、尺度为1的方向波高频系数子带Q，从方向波高频系数子带Q中选出坐标为(2i-1,2j-1)、(2i-1,2j)、(2i,2j-1)、(2i,2j)的四个高频系数，作为高频系数C的四个子节点，其中，i表示高频系数的行序号，j表示高频系数的列序号。

第三步，对其他两个方向上的方向波高频系数子带，按照第一步和第二步的步骤，分别建立方向波域HMT模型，得到3个方向上的方向波域HMT模型。

利用待降斑图像的3个方向上的方向波高频系数子带，按照下式，分别构造3个方向上的方向波域HMT模型参数集：

$\mathrm{\θ}=\{{\mathrm{\μ}}_{j,m},{\mathrm{\σ}}_{j,m}^2,P_j\left(m\right),{\mathrm{\ξ}}_{j,k}^{m,n}\}$

其中，θ表示方向波域HMT模型参数集，μ_j,m表示节点j的状态变量取值为m时方向波高频子带的均值，

表示节点j的状态变量取值为m时方向波高频子带的方差，P_j(m)表示节点j的状态变量取值为m时的概率，

表示节点j的状态变量取值为m时节点k的状态变量取值为n的转移概率，节点k是节点j的父节点，m、n表示状态值，m,n∈{1,2}，j、k表示节点序号。

步骤4，估计方向波域HMT模型参数集。

采用期望最大化方法，分别对3个方向上的方向波域HMT模型参数集进行训练，得到3个方向上的方向波域HMT模型参数集的最优估计值。期望最大化方法的具体步骤如下：

第一步，将方向波域HMT模型参数集θ中的所有参数初始化为0，设置迭代计数器l=0，完成初始化。

第二步，对每个方向波变换高频系数，计算其节点状态变量联合概率密度函数p(S|w,θ^l)，其中，S表示节点状态，w表示该高频系数的值，θ^l表示方向波域HMT模型参数集θ迭代l次后的值。

第三步，按照下式，计算方向波域HMT模型参数集θ迭代l+1次后的值：

θ^l+1=arg max_θE[ln f(w,S|θ^l)|w,θ^l]

其中，θ^l+1表示方向波域HMT模型参数集θ迭代l+1次后的值，arg max(·)表示取最大函数值对应参数的操作，E(·)表示取期望操作，ln(·)表示取对数操作，f(w,S|θ^l)表示方向波高频子带系数的后验概率，S表示节点状态，w表示该高频系数的值，θ^l表示方向波域HMT模型参数集θ迭代l次后的值。

第四步，计算方向波域HMT模型参数集两次迭代结果之间的误差值，若该误差值小于预设阈值1×10^-4，则循环终止，得到方向波域HMT模型参数集的最优估计值

否则迭代计数器l=l+1，方向波域HMT模型参数集θ^l=θ^l+1，返回第二步继续循环，其中，θ^l表示方向波域HMT模型参数集θ迭代l次后的值，θ^l+1表示方向波域HMT模型参数集θ迭代l+1次后的值。

步骤5，估计无斑图像的方向波高频子带系数。

利用3个方向上的方向波域HMT模型参数集的最优估计值，按照贝叶斯最小均方误差准则，估计无斑图像的方向波高频子带系数。降斑公式如下：

${\hat{g}}_m=\underset{m=1}{\overset{2}{\mathrm{\Σ}}}p(S=m)\×\frac{{\mathrm{\σ}}_m^2}{{\mathrm{\σ}}_N^2+{\mathrm{\σ}}_m^2}{g}_m$

其中，

表示估计得到的无斑图像的方向波高频子带系数，g_m表示待降斑图像的节点状态变量取值为m的方向波高频系数，∑表示求和操作，S表示节点状态，p(S=m)表示节点状态S取值为m时的节点状态概率，

表示方向波高频系数g_m所在的方向波高频子带的方差，

表示方向波高频系数g_m所在的方向波高频子带的噪声N的方差，m表示节点的状态值，m∈{1,2}。

步骤6，方向波逆变换。

对待降斑图像的方向波低频子带系数和估计的无斑图像所有方向波高频子带系数进行方向波逆变换，得到估计的无斑图像。

步骤6，输出估计的无斑图像。

下面结合附图2和附图3的仿真图对本发明的效果做进一步说明。

1.仿真实验条件：

本发明的硬件测试平台是：处理器为Inter Core2Duo CPU E8200，主频为2.67GHz，内存2GB，软件平台为：Windows7旗舰版32位操作系统和MatlabR2012b。本发明的输入图像分别为合成孔径雷达SAR图像田野图field和合成孔径雷达SAR图像机场图river，大小都为256×256，格式都为BMP。

2.仿真内容：

本发明用到的两个对比方法分别如下：

黄海燕等人在文献“一种改进的SAR图像LEE滤波算法”中提出的合成孔径雷达SAR图像降斑方法，简称Lee滤波方法。

武昕伟等人在文献“一种基于小波域隐马尔可夫模型的SAR相干斑抑制算法”中提出的合成孔径雷达SAR图像降斑方法，简称基于WHMT的滤波。

在两组实验中，我们采用以下指标来衡量本发明方法的性能：

1，等效视数（ENL），其作用为衡量降斑方法对相干斑噪声的抑制效果。

2，边缘保持指数（EPI），该指标用来衡量降斑图像中边缘的保持效果，该指数越大，说明降斑方法在边缘的保持方面越好。

图2为仿真实验中本发明与现有技术对合成孔径雷达SAR图像田野图field降斑效果对比图。其中，图2(a)为输入的合成孔径雷达SAR图像的田野图field，其大小为256×256，图2(b)为现有技术中采用Lee滤波方法对合成孔径雷达SAR图像田野图field进行降斑的结果图，图2(c)为现有技术中采用基于WHMT的方法对合成孔径雷达SAR图像田野图field进行降斑的结果图，图2(d)为本发明方法对合成孔径雷达SAR图像田野图field进行降斑的结果图。

图3为仿真实验中本发明与现有技术对合成孔径雷达SAR图像河流图river降斑效果对比图。其中，图3(a)为输入的合成孔径雷达SAR图像河流图river，其大小为256×256，图3(b)为现有技术中采用Lee滤波方法对合成孔径雷达SAR图像河流图river进行降斑的结果图，图3(c)为现有技术中采用基于WHMT的方法对合成孔径雷达SAR图像河流图river进行降斑的结果图，图3(d)为本发明方法对合成孔径雷达SAR图像河流图river进行降斑的结果图。

3.仿真结果分析：

表1为这两幅真实SAR图像三种滤波结果的等效视数（ENL）及边缘保持指数（EPI），其中等效视数（ENL）分别利用图2(a)和图3(a)中标示的A区域计算。

结合表1、附图2及附图3来看，Lee滤波方法滤波效果一般，主要是在平滑区域仍有大量相干斑噪声残留，而基于WHMT的滤波方法在平滑区域对相干斑噪声的抑制效果较为明显，但是其未能较好地保持图像的细节，尤其是在附图2(b)中田垄等地方出现了模糊现象，而本发明方法无论在平滑区域的滤波效果方面还是边缘结构保持方面，都明显优于前两种方法，表1中的数据也证明了这一点。

表1两幅真实SAR图像三种滤波结果的等效视数（ENL）及边缘保持指数（EPI）

以上实验表明：本发明方法无论是在平滑区域的滤波效果方面还是边缘结构保持方面都表现出了良好的性能，能够解决现有技术中存在边缘结构等信息无法有效保持及相干斑噪声滤除不彻底的问题，是一种非常实用的SAR图像滤波方法。

Claims (6)

1.一种基于方向波域HMT模型的SAR图像降斑方法，包括如下步骤：

(1)输入图像：

任选一幅SAR图像作为待降斑图像；

(2)建立非对数加性模型：

按照下式，对待降斑图像建立非对数加性模型，将待降斑图像转化为地物的真实后向散射强度和加性噪声的和：

I=RX=X+(R-1)X=X+N

其中，I表示待降斑图像，R表示相干斑噪声，R的均值为1，X表示地物的真实后向散射强度，N表示将要滤除的加性噪声；

(3)构造方向波域HMT模型参数集：

(3a)对待降斑图像进行尺度为2的方向波变换，得到待降斑图像的方向波低频系数子带和3个方向上的方向波高频系数子带；

(3b)利用待降斑图像的3个方向上的方向波高频系数子带，分别建立3个方向上的方向波域HMT模型；

(3c)利用待降斑图像的3个方向上的方向波高频系数子带，分别构造3个方向上的方向波域HMT模型参数集；

(4)估计方向波域HMT模型参数集：

采用期望最大化方法，分别对3个方向上的方向波域HMT模型参数集进行训练，得到3个方向上的方向波域HMT模型参数集的最优估计值；

(5)估计无斑图像的方向波高频子带系数：

利用3个方向上的方向波域HMT模型参数集的最优估计值，按照贝叶斯最小均方误差准则，估计无斑图像的方向波高频子带系数；

(6)方向波逆变换：

对待降斑图像的方向波低频子带系数和估计的无斑图像所有方向波高频子带系数，进行方向波逆变换，得到估计的无斑图像；

(7)输出估计的无斑图像。

2.根据权利要求1所述基于方向波域HMT模型的SAR图像降斑方法，其特征在于：步骤(3a)所述方向波变换的具体步骤如下：

第一步，从方向集{0°，30°，45°，60°，90°，120°，135°，150°}中，任意选择两个不平行的方向，分别作为方向波变换的变换方向和队列方向，构造采样矩阵：

$M=\left[\begin{array}{c}{d}_1\\ d_2\end{array}\right]$

其中，M表示采样矩阵，d₁表示沿变换方向的整数行向量，d₂表示沿队列方向的整数行向量；

第二步，将整数行向量d₁和d₂组成的网格在待降斑图像上进行L次平移，每次平移过程中，保留待降斑图像中和网格结点重合的像素点，得到待降斑图像的L个陪集，其中L表示采样矩阵M的行列式的绝对值；

第三步，沿变换方向和队列方向对待降斑图像的每个陪集分别进行2次一维小波变换，得到待降斑图像的方向波低频系数子带和3个方向上不同尺度的方向波高频系数子带。

3.根据权利要求1所述基于方向波域HMT模型的SAR图像降斑方法，其特征在于：步骤(3b)所述的方向波域HMT模型，按如下步骤建立：

第一步，在3个方向上的方向波高频系数子带中，任选一个方向上的尺度为2的方向波高频子带H，从方向波高频系数子带H中任选一个坐标为(i,j)的高频系数C，作为方向波域HMT模型中的一个父节点，其中，i表示高频系数的行序号，j表示高频系数的列序号；

第二步，在3个方向上的方向波高频系数子带中，选取和方向波高频子带H方向相同、尺度为1的方向波高频系数子带Q，从方向波高频系数子带Q中选出坐标为(2i-1,2j-1)、(2i-1,2j)、(2i,2j-1)、(2i,2j)的四个高频系数，作为高频系数C的四个子节点，其中，i表示高频系数的行序号，j表示高频系数的列序号；

第三步，对其他两个方向上的方向波高频系数子带，按照第一步和第二步的步骤，分别建立方向波域HMT模型，得到3个方向上的方向波域HMT模型。

4.根据权利要求1所述的基于方向波域HMT模型的SAR图像降斑方法，其特征在于：步骤(3c)所述的方向波域HMT模型参数集，按下式构造：

$\mathrm{\θ}=\{{\mathrm{\μ}}_{j,m},{\mathrm{\σ}}_{j,m}^2,P_j\left(m\right),{\mathrm{\ξ}}_{j,k}^{m,n}\}$

其中，θ表示方向波域HMT模型参数集，μ_j,m表示节点j的状态变量取值为m时方向波高频子带的均值，

表示节点j的状态变量取值为m时方向波高频子带的方差，P_j(m)表示节点j的状态变量取值为m时的概率，

表示节点j的状态变量取值为m时节点k的状态变量取值为n的转移概率，节点k是节点j的父节点，m、n表示状态值，m,n∈{1,2}，j、k表示节点序号。

5.根据权利要求1所述的基于方向波域HMT模型的SAR图像降斑方法，其特征在于：步骤(4)所述的期望最大化方法，按如下步骤进行：

第一步，将方向波域HMT模型参数集θ中的所有参数初始化为0，设置迭代计数器l=0，完成初始化；

第二步，对每个方向波变换高频系数，按照p(S|w,θ^l)公式，计算其节点状态变量联合概率密度，S表示节点状态变量，w表示该高频系数的值，θ^l表示方向波域HMT模型参数集θ迭代l次后的值；

第三步，按照下式，更新方向波域HMT模型参数集：

θ^l+1=arg max_θE[ln f(w,S|θ^l)|w,θ^l]

其中，θ^l+1表示方向波域HMT模型参数集θ迭代l+1次后的值，arg max(·)表示取最大函数值对应参数的操作，E(·)表示取期望操作，ln(·)表示取对数操作，f(w,S|θ^l)表示方向波高频子带系数的后验概率，S表示节点状态，w表示该高频系数的值，θ^l表示方向波域HMT模型参数集θ迭代l次后的值；

第四步，计算方向波域HMT模型参数集两次迭代结果之间的误差值，若该误差值小于预设阈值1×10^-4，则循环终止，得到方向波域HMT模型参数集的最优估计值

否则迭代计数器l=l+1，方向波域HMT模型参数集θ^l=θ^l+1，返回第二步继续循环，其中，θ^l表示方向波域HMT模型参数集θ迭代l次后的值，θ^l+1表示方向波域HMT模型参数集θ迭代l+1次后的值。

6.根据权利要求1所述的基于方向波域HMT模型的SAR图像降斑方法，其特征在于：步骤(5)所述的贝叶斯最小均方误差准则是指，按下式估计无斑图像方向波高频子带系数：

${\hat{g}}_m=\underset{m=1}{\overset{2}{\mathrm{\Σ}}}p(S=m)\×\frac{{\mathrm{\σ}}_m^2}{{\mathrm{\σ}}_N^2+{\mathrm{\σ}}_m^2}{g}_m$

其中，

表示估计得到的无斑图像的方向波高频子带系数，g_m表示待降斑图像的节点状态变量取值为m的方向波高频系数，∑表示求和操作，S表示节点状态，p(S=m)表示节点状态S取值为m时的节点状态概率，

表示方向波高频系数g_m所在的方向波高频子带的方差，

表示方向波高频系数g_m所在的方向波高频子带的噪声N的方差，m表示节点的状态值，m∈{1,2}。

Images