CN104239892A - 基于ksvd训练字典的sar图像混合模型拟合方法 - Google Patents

Abstract

本发明公开了一种基于KSVD算法训练字典的SAR图像混合模型拟合方法，它属于图像处理技术领域，主要解决SAR图像杂波分布模型种类繁多且对实际地物的建模能力有限的问题。其拟合过程为：基于对数正态分布和威布尔分布的概率密度函数，建立初始化分布模型字典；对图像中的每一类提取训练数据；利用KSVD算法训练初始化分布模型字典，得到新的字典集合；对图像中的每一类提取测试数据；用测试数据和新的字典集合作为OMP算法的输入来重构测试数据的稀疏系数；用字典集合中的字典乘以与其对应的稀疏系数就可得到测试数据的概率密度，从而得到其拟合直方图。本发明具有拟合精度高的优点，可用于SAR图像统计模型的拟合。

Description

基于KSVD训练字典的SAR图像混合模型拟合方法

技术领域

本发明属于图像处理领域，特别是一种涉及SAR图像统计建模的方法，可应用于SAR图像统计建模。

背景技术

多年来，合成孔径雷达(Synthetic Aperture Radar，SAR)已经被广泛应用于洪水监测、城市规划、土地利用和医学等诸多领域，这样的成像系统较光学成像系统具有全天时、全天候的成像特点和较强的地表穿透能力。对SAR图像统计特性的研究已成为解译SAR图像的一部分。随着SAR图像的乘性相干斑噪声的存在和分辨率的提高，传统用于低分辨率SAR图像统计建模的方法已不能对SAR图像数据进行精确的统计建模，特别是高分辨率条件下的SAR图像。单一模型虽能描述SAR图像的强度实测数据，但其精度不高。

SAR图像统计模型可分为两大类：

参量模型。参量模型的基本思想是利用参数估计的方法从已知的统计分布中确定图像数据的统计模型。参量模型有经验分布模型、由乘积发展的统计模型、由广义中心极限定理发展的统计模型、联合分布模型四类。经验分布有对数正态(Log-normal)分布、威布尔(Weibull)分布和广义Gamma分布(Generalized GammaDistribution，GΓ)。由广义中心极限定理发展的统计模型包括称α稳定分布(Symmetricα-stable distribution，SαS)、广义伽马瑞利分布(Generalized GammaRayleigh Distribution，GΓR)和广义高斯瑞利分布(Generalized Gauss RayleighDistribution，GGR)。联合分布模型是对单一分布模型取不同的参数，再按照场景中不同地物所占比例将其加权，最终得到表征地物统计分布的联合分布模型。参量模型是通过几个简单的参数对SAR图像数据进行建模来描述数据的分布，可压缩大量的数据信息，因此它在SAR图像统计建模中得到了广泛的应用。

非参量模型。非参量模型有Parzen窗法、支持向量机法和神经网络等。与参量模型相比，非参量模型建模比较灵活，适合对未知的复杂的概率密度函数进行估计，一般需要大样本数据的支持，估计精度较高，但其操作复杂，计算耗时大。

发明内容

本发明的目的在于克服上述已有技术的不足，提出一种基于广义K均值(KSVD)算法训练字典的SAR图像混合模型拟合方法，利用混合模型能较好地拟合不同的地物类型，从而取得更精确的建模结果。

本发明具体实现过程如下：

(1)基于对数正态分布和威布尔分布的概率密度函数，通过扩大这两个分布参数的范围，建立图像的初始化分布模型字典T₀；

(2)对图像中的每一类提取训练数据{getdata1,getdata2,…,getdatac,…,getdataC}，其中C为图像中的类别数；

(3)用图像的初始化分布模型字典T₀和训练数据{getdata1,getdata2,…,getdatac,…,getdataC}作为广义K均值(KSVD)算法的输入得到字典集合{dictionary1,dictionary2,…,dictionaryc,…,dictionaryC}；

(4)对图像中的每一类提取测试数据{data1,data2,…,datac,…,dataC}；

(5)用测试数据{data1,data2,…,datac,…,dataC}和字典集合{dictionary1,dictionary2,…,dictionaryc,…,dictionaryC}作为正交匹配追踪(OrthogonalMatching Pursuit,OMP)算法的输入来重构测试数据{data1,data2,…,datac,…,dataC}的稀疏系数{coeff1,coeff2,…,coeffc,…,coeffC}。

(6)用字典集合{dictionary1,dictionary2,…,dictionaryc,…,dictionaryC}中的第c类的字典dictionaryc乘以与其对应的重构的稀疏系数{coeff1,coeff2,…,coeffc,…,coeffC}中的第c类的稀疏系数coeffc，c取1至C，就可得到图像中每一类测试数据的概率密度。

实现本发明目的的技术方案：使用广义K均值(KSVD)算法训练由对数正态分布和威布尔分布构成的字典，得到的概率密度能够比较精确地反映图像统计信息，有效地克服了单一模型不能得到精确建模结果的缺陷。本发明与现有技术相比具有以下优点：

1、本发明由于使用广义K均值(KSVD)算法训练统计模型字典，克服了传统方法中对分布的参数估计速度慢、过程复杂的缺点，使其拟合过程简单；

2、本发明由于采用了混合模型，使其能够对均匀区域、不均匀区域和极不均匀区域都能得到较好的拟合结果；

3、仿真结果表明，本发明方法较单一模型能更有效地对SAR图像进行统计建模。

附图说明

图1是本发明的流程示意图；

图2是一幅含有均匀区域和不均匀区域的两类地物背景SAR图像；

图3是本发明对图2中不均匀区域仿真的结果图；

图4是本发明对图2中均匀区域仿真的结果图；

图5是一幅含有均匀区域和极不均匀区域的两类地物背景SAR图像；

图6是本发明对图5中极不均匀区域仿真的结果图。

具体实施方式

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

步骤一、基于对数正态分布和威布尔分布的概率密度函数，通过扩大这两个分布参数的范围，建立图像的初始化分布模型字典T₀。

对数正态分布的概率密度函数：

$f\left(x\right)=\frac{1}{\mathrm{\σ}\sqrt{2\mathrm{\π}}}\exp [-\frac{(\ln x-m)}{2{\mathrm{\σ}}^2}]$

其中，x为图像中像素点的灰度值，m为图像中像素点的灰度值的均值，σ为图像中像素点的灰度值的方差。

威布尔分布的概率密度函数：

$f\left(x\right)=\frac{\mathrm{\η}}{{\mathrm{\μ}}^{\mathrm{\η}}}{x}^{\mathrm{\η}-1}\exp [-{\left(\frac{x}{\mathrm{\μ}}\right)}^{\mathrm{\η}}]$

其中，x为图像中像素点的灰度值，μ为图像的尺度参数，η为图像的形状参数。

对分布的参数进行估计是分布的核心问题，是长期制约分布实际应用的主要技术瓶颈。在本发明中，优选地，选择基于Mellin变换的参数估计方法来估计上述对数正态分布和威布尔分布的参数。对数正态分布的参数估计公式如下：

$\left\{\begin{array}{c}m=k_1\\ \mathrm{\σ}=\sqrt{K_2}\end{array}\right.$

其中，k_i是第i阶对数累计量，计算公式如下：

$\left\{\begin{array}{c}{k}_1=E\left\{\ln x\right\}\\ k_{2}E\left\{{(\ln x-k_1)}^2\right\}\end{array}\right.$

其中，x为图像中像素点的灰度值。

威布尔分布的参数估计公式如下：

$\left\{\begin{array}{c}{k}_1=\ln u+{\mathrm{\η}}^{-1}\mathrm{\ψ}\left(1\right)\\ k_2={\mathrm{\η}}^{-2}\mathrm{\ψ}\left(\mathrm{1,1}\right)\end{array}\right.$

其中，ψ(1)是普西(digamma)函数，ψ(1,1)是普西(digamma)函数的一阶导数，k₁和k₂的计算公式与对数正态分布中的计算公式一样。

利用上面的公式可得到对数正态分布和威布尔分布的参数m、σ、μ、η，对这四个参数分别取0.1、0.1、1、1的步长，扩大参数的范围。m分别取m-0.5、m-0.4、m-0.3、m-0.2、m-0.1、m、m+0.1、m+0.2、m+0.3、m+0.4，σ分别取σ-0.6、σ-0.5、σ-0.4、σ-0.3、σ-0.2、σ-0.1、σ、σ+0.1、σ+0.2、σ+0.3、σ+0.4、σ+0.5、σ+0.6、σ+0.7、σ+0.8，μ分别取μ-5、μ-4、μ-3、μ-2、μ-1、μ、μ+1、μ+2、μ+3、μ+4，η分别取η-6、η-5、η-4、η-3、η-2、η-1、η、η+1、η+2、η+3、η+4、η+5、η+6、η+7、η+8。把上面的10个含有符号m的参数和15个含有符号σ的参数进行笛卡尔积，得到150组既含有符号m、又含有符号σ的参数。把这150组参数代入对数正态分布的概率密度函数中，得到一个b行、n1列的矩阵，其中b为灰度的级数，优选地，b为256，n1为150。同理，把上面的10个含有符号μ的参数和15个含有符号η的参数进行笛卡尔积，得到150组既含有符号μ、又含有符号η的参数。把这150组参数代入威布尔分布的概率密度函数中，得到一个b行、n2列的矩阵，其中b为灰度的级数，优选地，b为256，n2为150。把以上两个矩阵合并就可得到一个b行(优选地，256行)300列的矩阵，这个矩阵即初始化分布模型字典T₀。

步骤二、对图像的每一类提取训练数据Y＝{getdata1,getdata2,…,getdatac,…,getdataC}。

设图像有C类，对图像C类中的每一类各截取一个训练样本，得到训练样本S1，S2,…,Sc,…,SC，在每类训练样本中提取n*n的块，计算第c块的概率密度f_c(x)，得到的概率密度即为训练数据{getdata1,getdata2,…,getdatac,…,getdataC}，getdatac＝f_c(x)。在本发明中，优选地，n＝20。

步骤三、用图像的初始化分布模型字典T₀和训练数据{getdata1,getdata2,…,getdatac,…,getdataC}作为广义K均值(KSVD)算法的输入得到字典集合{dictionary1,dictionary2,…,dictionaryc,…,dictionaryC}。

在广义K均值(KSVD)算法中，对于要被表示的训练数据Y，可利用分布模型字典T与稀疏系数X乘积得到。衡量字典好坏的标准是满足一定稀疏度的前提下的误差的大小，即：

$\underset{T,X}{\min }\left\{{\left|\right|Y-\mathrm{TX}\left|\right|}_F^2\right\}\mathrm{subject\; to}\∀i,{\left|\right|x_i\left|\right|}_0\≤K$

其中，Y在本发明中是要表示的SAR图像像素点的概率密度，T是分布模型字典，X是稀疏系数，||·||_F是弗罗贝尼乌斯矩阵范数,上面的2表示平方，||·||₀是零范数，x_i为稀疏系数X的第i个行矢量，i的取值范围为1至(n1+n2)，(n1+n2)是初始化分布模型字典T₀的列数，K为稀疏度，即稀疏系数X的第i个行矢量x_i中非零分量数目的上限，K的取值范围为1至C，C是图像的类别数。

实现该步骤的具体实现过程如下：

3a)令当前的分布模型字典T等于初始化分布模型字典T₀。

3b)控制稀疏度为K的情况下，用步骤二得到的训练数据Y＝{getdata1,getdata2,…,getdatac,…,getdataC}和当前的分布模型字典T作为正交匹配追踪(Orthogonal Matching Pursuit,OMP)算法的输入得到训练数据的的稀疏系数X＝{coef11,coef22,…,coefcc,…,coefCC}。

3c)更新分布模型字典T：

逐列更新字典T，t_i为分布模型字典T的第i个列矢量，x_i为稀疏系数X的第i个行矢量。其中，i的取值范围为1至(n1+n2)，(n1+n2)是初始化分布模型字典T₀的列数，X共有C个列矢量，C为图像的类别数。

1)将T的第i个列矢量和X的第i个行矢量置零，i取1至(n1+n2)，(n1+n2)是初始化分布模型字典T₀的列数，得到T'和X'，计算此时的误差矩阵E'＝Y-T'X'。

2)对于行矢量x_i，x_ic是行矢量x_i中的第c个元素，其中c＝1,2,3,…,C，C为X列元素的个数。设w为行矢量x_i中非零元素的下标的集合，即w＝{c|1≤c≤C,x_ic≠0}，根据w中的c选择误差矩阵E'对应的列，即E_r＝{E_c'}，其中c∈w，E_c'表示E'的第c列，E_r是根据w中的c选出E'中对应的列组成的新矩阵。在本发明中，优选地，利用奇异值分解的方法对E_r进行分解，即E_r＝UΣV^T，U、V是E_r进行奇异值分解得到的矩阵，Σ.是一个对角矩阵，对角线上的值是E_r的奇异值，然后利用U中的第一列替换t_i，实现t_i的更新。然后选择V的第一列和Σ.对角线上的第一个值，也就是Σ(1,1)，两者相乘得到列向量，依次用这个列向量中的元素替换xi中的非零元素，完成x_i的更新。

3)得到更新的字典T'和稀疏系数X'之后，计算误差矩阵E'＝Y-T'X'，判断||E'||₂＜ε是否成立，ε是一个接近于零的常数，在本发明中，优选地，设置ε为0.001。如果成立，训练过程结束，否则返回步骤3b)。

步骤四、对图像中的每一类提取测试数据{data1,data2,…,datac,…,dataC}。

设图像有C类，对图像C类中的每一类各截取一个测试样本，得到测试样本T1,T2,…,Tc,…,TC，在每类测试样本中提取d*d的块，计算第c块的概率密度g_c(x)，得到的概率密度即为测试数据{data1,data2,…,datac,…,dataC}，datac＝g_c(x)。在本发明中，优选地，d＝80。

步骤五、用测试数据{data1,data2,…,datac,…,dataC}和步骤三得到的字典集合{dictionary1,dictionary2,…,dictionaryc,…,dictionaryC}作为正交匹配追踪(OrthogonalMatching Pursuit,OMP)算法的输入来重构测试数据{data1,data2,…,datac,…,dataC}的稀疏系数{coeff1,coeff2,…,coeffc,…,coeffC}。

步骤六、用字典集合{dictionary1,dictionary2,…,dictionaryc,…,dictionaryC}中的第c类的字典dictionaryc乘以与其对应的重构的测试数据的稀疏系数{coeff1,coeff2,…,coeffc,…,coeffC}中的第c类的稀疏系数coeffc，c取1至C，就可得到测试数据{data1,data2,…,datac,…,dataC}中每一类的概率密度，所求的概率密度即为SAR图像中每一类测试数据的混合模型的拟合。

本发明的效果可通过以下仿真进一步说明：

1仿真内容：应用本发明方法、对数正态单模型、威布尔单模型分别对两幅SAR图像的地物进行直方图拟合实验。

2仿真实验结果

A含有均匀区域和不均匀区域的两类简单地物背景SAR图像的实验结果

用本发明方法、对数正态单模型、威布尔单模型对两类简单地物背景SAR图像的不均匀区域和均匀区域分别进行直方图拟合，其效果比较分别如图3、图4所示。图3(a)为从图2中截取的不均匀区域的SAR图像，该SAR图像纹理信息不是很复杂；图3(b)为用对数正态单模型对图3(a)进行直方图拟合得到的结果；图3(c)为用威布尔单模型对图3(a)进行直方图拟合得到的结果；图3(d)为用本发明方法对图3(a)进行直方图拟合得到的结果。由图3可见，对数正态单模型和威布尔单模型对不均匀区域部分的拟合不太精确，而本发明的拟合精度较高。图4(a)为从图2中截取的均匀区域的SAR图像，该SAR图像纹理信息较为单调，地物背景组成简单；图4(b)为用对数正态单模型对图4(a)进行直方图拟合得到的结果；图4(c)为用威布尔单模型对图4(a)进行直方图拟合得到的结果；图4(d)为用本发明方法对图4(a)进行直方图拟合得到的结果。由图4可见，本发明方法、对数正态单模型、威布尔单模型都能较好地拟合均匀区域，但本发明方法的拟合精度相对较高。

B含有均匀区域和极不均匀区域的两类复杂地物背景SAR图像的实验结果

用本发明方法、对数正态单模型、威布尔单模型对两类复杂地物背景SAR图像的极不均匀区域进行直方图拟合，其效果比较如图6所示。图6(a)为图5中截取的极不均匀区域的SAR图像；图6(b)为用对数正态单模型对图6(a)进行直方图拟合得到的结果；图6(c)为用威布尔单模型对图6(a)进行直方图拟合得到的结果；图6(d)为用本发明方法对图6(a)进行直方图拟合得到的结果。由图6可见，本发明的拟合精度较高，本发明拟合的直方图与原始直方图几乎重合。相比之下，对数正态单模型和威布尔单模型方法的直方图拟合不够精确。

Claims (10)

1.一种基于KSVD训练字典的SAR图像混合模型拟合方法，包括如下步骤：

步骤1，基于对数正态分布和威布尔分布的概率密度函数，通过扩大这两个分布参数的范围，建立图像的初始化分布模型字典T₀；

步骤2，对图像中的每一类提取训练数据{getdata1,getdata2,…,getdatac,…,getdataC}，其中C为图像的类别数；

步骤3，用图像的初始化分布模型字典T₀和训练数据{getdata1,getdata2,…,getdatac,…,getdataC}作为广义K均值(KSVD)算法的输入得到字典集合{dictionary1,dictionary2,…,dictionaryc,…,dictionaryC}；

步骤4，对图像中的每一类提取测试数据{data1,data2,…,datac,…,dataC}；

步骤5，用测试数据{data1,data2,…,datac,…,dataC}和步骤三得到的字典集合{dictionary1,dictionary2,…,dictionaryc,…,dictionaryC}作为正交匹配追踪(OrthogonalMatching Pursuit,OMP)算法的输入来重构测试数据{data1,data2,…,datac,…,dataC}的稀疏系数{coeff1,coeff2,…,coeffc,…,coeffC}；

步骤6，用步骤3得到的字典集合{dictionary1,dictionary2,…,dictionaryc,…,dictionaryC}中的第c类的字典dictionaryc乘以与其对应的重构的测试数据的稀疏系数{coeff1,coeff2,…,coeffc,…,coeffC}中的第c类的稀疏系数coeffc，c取1至C，就可得到测试数据{data1,data2,…,datac,…,dataC}中每一类的概率密度，所求的概率密度即为SAR图像中每一类测试数据的混合模型的拟合。

2.根据权利要求1所述的混合模型拟合方法，其中，步骤1具体包括：

(1a)对数正态分布的概率密度函数：

$f\left(x\right)=\frac{1}{\mathrm{\σ}\sqrt{2\mathrm{\π}}}\exp [-\frac{(\ln x-m)}{2{\mathrm{\σ}}^2}]$

其中，x为图像中像素点的灰度值，m为图像中像素点的灰度值的均值，σ为图像中像素点的灰度值的方差；

威布尔分布的概率密度函数：

$f\left(x\right)=\frac{\mathrm{\η}}{{\mathrm{\μ}}^{\mathrm{\η}}}{x}^{\mathrm{\η}-1}\exp [-{\left(\frac{x}{\mathrm{\μ}}\right)}^{\mathrm{\η}}]$

其中，x为图像中像素点的灰度值，μ为图像的尺度参数，η为图像的形状参数。

3.根据权利要求2所述的混合模型拟合方法，其中，选择基于Mellin变换的参数估计方法来估计上述对数正态分布和威布尔分布的参数。

4.根据权利要求3所述的混合模型拟合方法，其中，对数正态分布的参数估计公式如下：

$\left\{\begin{array}{c}m=k_1\\ \mathrm{\σ}=\sqrt{K_2}\end{array}\right.$

其中，k_i是第i阶对数累计量，计算公式如下：

$\left\{\begin{array}{c}{k}_1=E\left\{\ln x\right\}\\ k_{2}E\left\{{(\ln x-k_1)}^2\right\}\end{array}\right.$

其中，x为图像中像素点的灰度值；

威布尔分布的参数估计公式如下：

$\left\{\begin{array}{c}{k}_1=\ln u+{\mathrm{\η}}^{-1}\mathrm{\ψ}\left(1\right)\\ k_2={\mathrm{\η}}^{-2}\mathrm{\ψ}\left(\mathrm{1,1}\right)\end{array}\right.$

其中，ψ(1,1)是普西(digamma)函数，ψ(1,1)是普西(digamma)函数的一阶导数。

5.根据权利要求2所述的混合模型拟合方法，其中，得到对数正态分布和威布尔分布的参数m、σ、μ、η后，

(1b)对(1a)的四个参数分别取0.1、0.1、1、1的步长，扩大参数的范围；

(1c)将(1b)含有符号m的参数和含有符号σ的参数进行笛卡尔积，得到n1组既含有符号m、又含有符号σ的参数，把这n1组参数代入对数正态分布的概率密度函数中，得到一个b行、n1列的矩阵，其中b为灰度的级数；把含有符号μ的参数和含有符号η的参数进行笛卡尔积，得到n2组既含有符号μ、又含有符号η的参数，把这n2组参数代入威布尔分布的概率密度函数中，得到一个b行、n2列的矩阵，其中b为灰度的级数；把这两个矩阵合并得到一个b行(n1+n2)列的矩阵，这个矩阵称为初始化分布模型字典T₀。

6.根据权利要求5所述的混合模型拟合方法，其中，b为256。

7.根据权利要求1所述的混合模型拟合方法，其中，步骤3具体包括：

3a)令当前的分布模型字典T等于初始化分布模型字典T₀；

3b)控制稀疏度为K的情况下，用步骤2得到的训练数据Y＝{getdata1,getdata2,…,getdatac,…,getdataC}和当前的分布模型字典T作为正交匹配追踪(Orthogonal Matching Pursuit,OMP)算法的输入得到训练数据的的稀疏系数X＝{coef11,coef22,…,coefcc,…,coefCC}；

3c)逐列更新字典T：t_i为分布模型字典T的第i个列矢量，x_i为稀疏系数X的第i个行矢量，其中，i的取值范围为1至(n1+n2)，(n1+n2)是初始化分布模型字典T₀的列数，X共有C个列矢量，C为图像的类别数；

1)将T的第i个列矢量和X的第i个行矢量置零，i取1至(n1+n2)，(n1+n2)是初始化分布模型字典T₀的列数，得到T'和X'，计算此时的误差矩阵E'＝Y-T'X'；

2)对于行矢量x_i，x_ic是行矢量x_i中的第c个元素，其中c＝1,2,3,…,C，C为X列元素的个数，设w为行矢量x_i中非零元素的下标的集合，即w＝{c|1≤c≤C,x_ic≠0}，根据w中的c选择误差矩阵E'对应的列，即E_r＝{E_c'}，其中c∈w，E_c'表示E'的第c列，E_r是根据w中的c选出E'中对应的列组成的新矩阵；利用奇异值分解的方法对E_r进行分解，即E_r＝UΣV^T，U、V是E_r进行奇异值分解得到的矩阵，Σ.是一个对角矩阵，对角线上的值是E_r的奇异值，然后利用U中的第一列替换t_i，实现t_i的更新；然后选择V的第一列和Σ.对角线上的第一个值，也就是Σ(1,1)，两者相乘得到列向量，依次用这个列向量中的元素替换x_i中的非零元素，完成x_i的更新；

3)得到更新的字典T'和稀疏系数X'之后，计算误差矩阵E'＝Y-T'X'，判断||E'||₂＜ε是否成立，ε是一个接近于零的常数；如果成立，训练过程结束，否则返回步骤3b)。

8.根据权利要求1所述的混合模型拟合方法，其中，对图像中的每一类提取训练数据{getdata1,getdata2,…,getdatac,…,getdataC}时，取的块大小为20*20，每个块的概率密度就是一个训练数据。

9.根据权利要求1所述的混合模型拟合方法，其中，对图像中的每一类提取测试数据{data1,data2,…,datac,…,dataC}时，取的块大小为80*80，每个块的概率密度就是一个测试数据。

10.根据权利要求1所述的混合模型拟合方法，ε为0.001。