CN102129573A - 基于字典学习和稀疏表示的sar图像分割方法 - Google Patents

Abstract

本发明公开了一种基于字典学习和稀疏表示的SAR图像分割技术，主要解决现有特征提取需要花费大量时间以及距离测度具有一定缺陷的问题。其实现过程是：1)输入待分割图像，确定分割类数k；2)为待分割图像每个像素点提取p×p窗口得到测试样本集，从测试样本集中随机选择少量样本得到训练样本集；3)提取训练样本集小波特征；4)利用谱聚类算法对训练样本集进行划分；5)对每一类训练样本，应用K-SVD算法训练字典；6)求解测试样本在字典上稀疏表示向量；7)计算测试样本的重构误差函数；8)根据重构误差函数求测试样本标签，得到图像分割结果。本发明具有分割快速且效果好的优点，该技术可以进一步用于SAR图像自动目标识别。

Description

基于字典学习和稀疏表示的SAR图像分割方法

技术领域

本发明属于图像处理技术领域，涉及SAR图像分割，用于SAR图像目标识别。

背景技术

合成孔径雷达SAR具有全天时、全天候的探测与侦察能力。它利用脉冲压缩技术获得高的距离分辨率，利用合成孔径原理提高方位分辨率，从而相比真实孔径雷达在遥感领域具有独特的优势。对SAR图像的理解和解译属于图像处理范畴，还涉及信号处理，模式识别及机器学习等众多学科。由于SAR具有的独特作用，使得SAR图像的理解与解译在国防和民用领域正受到越来越广泛的关注，SAR图像分割作为SAR图像后续解译处理的关键环节之一，也就显得愈加重要。现有的SAR图像分割方法大致可以分为基于区域分割方法和基于边缘的分割方法，如阈值分割、聚类的方法、形态学的方法、及随机场的方法等。

目前已经有很多成熟的聚类算法被用到SAR图像分割中。但是大部分基于聚类的SAR图像分割方法都要对图像特征提取，特征提取虽然是线下过程，但是由于需要对图像的每一个像素点分别进行特征提取，所以这一过程十分耗时，同时在聚类过程中一般要使用某种距离测度，如欧氏距离、流形距离，来计算样本之间的相似度，使用某种距离则会受限于该种距离测度方法的固定缺陷的限制如欧氏距离仅对球形分布数据效果较好，而流形距离仅对具有流形分布的数据效果较好，在数据分布未知情况下，若强制使用某种距离测度方法，可能会导致最终SAR图像分割效果较差。

发明内容

本发明的目的在于克服上述已有技术的不足，提出了一种基于字典学习和稀疏表示的SAR图像分割方法，以降低图像特征提取花费的时间，并在图像数据分布未知情况下，提高SAR图像的分割结果。

实现本发明目的的技术思路是：从待分割图像随机选取的少量像素点提取纹理特征，得到纹理特征样本，然后对纹理特征样本利用谱聚类方法进行聚类，得到纹理特征样本的初始标签，根据纹理特征标签样本，进一步采用字典学习和稀疏表示方法得到图像中所有像素点的标签。具体步骤包括如下：

(1)输入待分割SAR图像，确定需要将图像划分的类数k；

(2)以待分割SAR图像中每个像素点为中心提取w×w的窗口样本，得到规模为m的测试样本集X，从提取的所有测试样本中随机选择n个作为训练样本集Y，其中n＜＜m；

(3)对规模为n的训练样本集Y，提取训练数据样本窗口中心所在像素点的小波能量特征，得到规模为n的小波特征数据集Y′，该小波特征数据集Y′中的样本与训练样本集Y中的样本一一对应；

(4)用谱聚类算法对小波特征数据集Y′进行聚类，得到小波特征数据集Y′的类别划分，同时得到训练样本集Y的类别划分Y_i，i＝1，2，3，...k；

(5)对划分后的每一类训练样本Y_i，i＝1，2，3，...k，利用K-SVD算法分别训练，得到相应的字典D_i，i＝1，2，3，...，k；

(6)利用正交匹配追踪OMP算法求解测试样本s在D_i上稀疏表示向量，得到稀疏向量αⁱ，i＝1，2，3，...，k，s∈X；

(7)令

表示测试样本s在D_i上的重构误差，对该误差进行高斯函数平滑，得到新的误差函数

(8)根据新的误差函数

求解测试样本S的标签

得到每个测试样本的类别标签l(S)∈{1，2，3，...，k}，即得到图像分割结果。

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

(1)本发明不需要对图像中的每一个像素点提取纹理特征，而是直接对图像中的每个像素点提取窗口作为样本，所以大大节省了图像线下提取特征所要花费的时间。

(2)在本发明中使用的字典学习稀疏表示方法能够深层挖掘表示图像的信息，根据字典对样本的重构误差最小化原则来确定样本的标签将会避免使用一般距离测度方法带来的缺陷，得到更好的图像分割结果。

附图说明

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

图2是本发明与现有方法对二类SAR图像的分割结果对比图；

图3是本发明与现有方法对三类SAR图像分割结果对比图。

具体实施方式

参照图1，本发明的具体实施过程如下：

步骤1.输入待分割图像，根据图像内容，判断需要识别的主要目标及背景，确定分割类数k，本实例中的k取值为2和3。

步骤2.对待分割图像提取训练样本和测试样本。

以图像中每个像素点为中心提取大小为p×p的窗口，得到规模为m的测试样本集F，m为待分割图像的像素点总数目，然后从测试样本集中随机选择n个样本，得到训练样本集Y，n远远小于m。

步骤3.对训练样本集Y提取小波特征。

SAR图像具有丰富的幅度、相位、极化和纹理信息，为了使每一个训练样本获得更加准确标记，在对训练样本集划分前需要对其进行纹理分析，提每个取训练样本的小波纹理特征，具体步骤如下：

(3a)取训练样本集Y中的训练样本T；

(3b)对训练样本T进行三层小波分解；

(3c)计算小波分解后的训练样本T的小波特征e：

$e=\frac{1}{w\×w}\underset{i=1}{\overset{w}{\mathrm{\Σ}}}\underset{j=1}{\overset{w}{\mathrm{\Σ}}}\left|\mathrm{coef}(i,j)\right|;$

其中，w×w为训练样本T对应窗口大小，本实例取值为16×16，coef(i，j)为小波分解后得到的小波子带中第i行第j列的系数值；

(3d)重复步骤(3a)-(3c)，计算下一训练样本的小波特征，直到计算完成训练样本集Y中所有训练样本的小波特征，得到小波特征数据集Y′。

步骤4.用谱聚类算法对小波特征数据集Y′进行聚类。

(4a)根据下式计算小波特征数据集Y′中的两个不同小波特征样本y′_i和y′_j之间的相似度值W_ij，构建相似度矩阵W：

$W_{\mathrm{ij}}=\exp (-\frac{{\left|\right|{y^{\′}}_i-{y^{\′}}_j\left|\right|}^2}{{2\mathrm{\σ}}^2}),i,j=\mathrm{1,2,3},...,n;$

W＝[W_ij]i，j＝1，2，3，...，n；

其中y′_i∈Y′，y′_i为Y′中的第i个样本，y′_j∈Y′，y′_j为Y′中的第j个样本，σ为尺度参数；

(4b)利用相似度矩阵W计算拉普拉斯矩阵：L＝Λ^-1/2WΛ^-1/2，

其中：Λ为对角矩阵， $\mathrm{\Λ}=\left[\begin{array}{ccccccc}{\mathrm{\Λ}}_{11}& 0& 0& ...& 0& ...& 0\\ 0& {\mathrm{\Λ}}_{22}& 0& ...& 0& ...& 0\\ 0& 0& {\mathrm{\Λ}}_{\mathrm{ii}}& ...& 0& ...& 0\\ 0& 0& 0& ...& 0& ...& {\mathrm{\Λ}}_{\mathrm{nn}}\end{array}\right]$

Λ_ii为对角线元素， ${\mathrm{\Λ}}_{\mathrm{ii}}=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{W}_{\mathrm{ij}},i=\mathrm{1,2},...n;$

(4c)对拉普拉斯矩阵L进行矩阵特征值分解，得到拉普拉斯矩阵L最大的前k个特征值λ₁≥λ₂≥...≥λ_k和与这些特征值对应的k个特征向量v¹，v²，...，v^k，利用k个特征向量构建特征向量矩阵U＝[v¹，v²，...，v^k]；

(4d)根据下式规范化特征向量矩阵U，得到规范化后的特征向量矩阵U′：

U′＝[U′_ij]，i＝1，2，3，...，n，j＝1，2，3，...，k；

其中U_ij为特征向量矩阵U中第i行第j列的元素；

(4e)用k-means算法对规范化后的特征向量矩阵U′按行聚类，得到对应小波特征数据集Y′的划分，将小波特征数据集Y′的划分结果作为训练样本集Y的划分结果，得到Y_i，i＝1，2，3，...k。

步骤5.对划分后的每一类训练数据利用K-SVD算法训练一个字典得到D_i，i＝1，...，k：

(5a)取训练样本的第j类数据Y_i，j∈{1，2，3，...，k}，j初值为1，令Y_T＝Y_j；

(5b)根据下式利用正交匹配追踪OMP方法，完成对Y_T的稀疏编码：

Y_T＝DX s.t.|x_i|₀≤t₀

其中，Y_T为需要稀疏表示的样本集合，D为用来稀疏表示的字典，X为样本经字典稀疏表示后的得到稀疏矩阵；x_i为稀疏矩阵X的第i列，t₀为预设定稀疏向量中的非零元素个数，称为稀疏度；

(5c)对字典D中所有原子d_i i＝1，2，3...L进行如下更新，原子d_i指字典中的第i个列向量，L为字典中原子个数：

首先，从样本集合Y_T中找出使用原子d_i的所有样本，取样本序号，得到具有n_ω个元素的序号集合ω：ω＝{k|1≤k≤n，xⁱ(k)≠0}，其中xⁱ为对样本集合Y_T稀疏编码后得到的稀疏矩阵X的第i行；

接着，令稀疏矩阵X的第i行元素全部为零，得到新矩阵X′，计算字典D与新矩阵X′对样本集合Y_T的重构误差矩阵E_i：

E_i＝Y_T-DX′；

然后，依据序号集合ω中所有序号，取所有n_ω个序号对应的误差矩阵E_i中的n_ω个列向量，构成特异误差矩阵

最后，应用公式

对特异误差矩阵

进行矩阵奇异值分解，并用分解行矩阵U的第一列替换所要更新的原子d_i，并用对角矩阵Δ中的第一行第一列的元素Δ(1，1)乘以分解列矩阵V的第一个列向量得到一个新的向量，利用该新的向量替换xⁱ，xⁱ为稀疏矩阵X的第i行；

(5d)判断是否满足j≥k，若满足，则得到字典D_i，i＝1，...，k，否则j＝j+1，返回步骤(5a)计算下一类训练数据对应字典；

(5e)重复步骤步骤(5a)-(5d)，直到得到所有类别的训练数据的对应字典D_i，i＝1，...，k。

步骤6.按照下式利用OMP算法求解每一测试样本S的稀疏向量αⁱ，i＝1，2，3，...，k，S∈F，S为测试样本集F中的样本；

S＝D_iαⁱ，s.t.|αⁱ|₀≤t₀ i＝1，2，3，...，k；

t₀为预设定稀疏向量中的非零元素个数，称为稀疏度；

步骤7.计算误差函数

令

表示测试样本S在对应字典D_i上的重构误差函数，对该误差函数进行高斯函数平滑，得到新误差函数G_σ是一个均方差参数为σ的高斯函数；

步骤8.根据新误差函数

求解测试样本S的标签l(S)：

$l\left(S\right)=\arg \min \hat{R}(S,D_i),i=\mathrm{1,2,3},...,k$

根据上式得到每个测试样本的类别标签l(S)∈{1，2，3，...，k}，即得到图像分割结果。

本发明的效果可以通过如下仿真实验进行验证。

(一)实验条件设置

本发明中SAR图像划分类别数k取值为2和3，分别在两幅SAR图像上进行了验证。

k＝2时，二类待分割SAR图像为X-波段SAR图像，大小为256×256，它包括山地和平坦地区两个部分，由于有部分山脊从平坦地区向山地延伸，一般分割方法很难对山脊分叉细节精确勾勒，如图2(a)所示。

k＝3时，三类待分割SAR图像为X-波段SAR图像，大小为300×400，它包括植被，河流，以及河中央陆地等目标物，由于对比度差异较小，河道轮廓很难被提取出来，如图3(a)所示。

本发明中所用训练样本数n＝6000，每一类的训练字典个数L＝128，窗口大小w×w取为13×13，稀疏度t₀＝4，高斯平滑参数σ＝8，

近似谱聚类选取的样本点个数为100，尺度参数θ＝2。

(二)仿真内容及结果

A)用现有的k-means聚类算法和

谱聚类算法以及本发明方法对二类SAR图像进行仿真实验，其分割结果如图2所示。其中图2(b)为k-means分割结果，图2(c)为

谱聚类分割结果，图2(d)为本发明方法分割结果。

从图2(b)可以看出，k-means分割结果尽管分割出了一些向纵深延伸的山脊，但同时山地部分也存在很多误分。

从图2(c)可以看出，

谱聚类分割结果相对于k-means的分割结果在山地误分方面有所改善，但是结果仍然不够好，划分结果杂点较多。

从图2(d)可以看出，本发明方法由于存在高斯平滑处理，山地部分的细节基本被平滑掉了，分割结果明显地呈现山地和平坦地区两部分，这也是符合地质勘探上对地形整体把握需要的。

B)用现有的k-means聚类算法和

谱聚类算法以及本发明方法对三类SAR图像进行仿真实验，其分割结果如图3所示。其中图3(b)为k-means分割结果，图3(c)为谱聚类分割结果，图3(d)为本发明方法分割结果。

从图3(b)可以看出，k-means聚类的结果中不仅河中央陆地与河道轮廓有很多混杂，甚至上方的植被也有很多被错分成河中央陆地的。

从图3(c)可以看出，

谱聚类分割结果在一些斑点细节的划分上有所改善，但是还存在严重的错分现象。

从图3(d)可以看出，本发明方法的分割结果中河中央陆地和植被的一致性得到了很大改善，同时，河道的轮廓也得到了较为准确的定位。

Claims (4)

1.一种基于字典学习和稀疏表示的SAR图像分割方法，包括以下步骤：

(1)输入待分割SAR图像，确定需要将图像划分的类数k；

(2)以待分割SAR图像中每个像素点为中心提取p×p的窗口样本，得到规模为m的测试样本集F，从提取的所有测试样本中随机选择n个作为训练样本集Y，其中n＜＜m；

(3)提取训练样本集Y中的训练样本的小波特征，得到规模为n的小波特征数据集Y′，该小波特征数据集Y′中的样本与训练样本集Y中的样本一一对应；

(4)用谱聚类算法对小波特征数据集Y′进行聚类，得到小波特征数据集Y′的类别划分，同时得到训练样本集Y的类别划分Y_i，i＝1，2，3，...k；

(5)将划分后的每一类训练样本Y_i，i＝1，2，3，...k，利用K-SVD算法分别训练，得到相应的字典D_i，i＝1，2，3，...，k；

(6)利用正交匹配追踪OMP算法求解测试样本S在D_i上稀疏表示向量，得到稀疏向量αⁱ，i＝1，2，3，...，k，S∈F；

(7)令

表示测试样本S在D_i上的重构误差函数，对该误差函数进行高斯函数平滑，得到误差函数

(8)根据新的误差函数求解测试样本S的标签得到每个测试样本的类别标签l(S)∈{1，2，3，...，k}，即得到图像分割结果。

2.根据权利要求1所述的字典学习和稀疏表示的SAR图像分割方法，其中步骤(3)所述的提取训练样本集Y中的训练样本的小波特征，按如下步骤进行：

(3a)取训练样本集Y中的训练样本T；

(3b)对训练样本T进行三层小波分解；

(3c)计算小波分解后的训练样本T的小波特征e：

$e=\frac{1}{w\×w}\underset{i=1}{\overset{w}{\mathrm{\Σ}}}\underset{j=1}{\overset{w}{\mathrm{\Σ}}}\left|\mathrm{coef}(i,j)\right|;$

其中，w×w为提取的训练样本T对应窗口的大小，coef(i，j)为训练样本T经小波分解后得到的小波子带中第i行第j列的系数值；

(3d)重复步骤(3a)-(3c)，计算下一训练样本的小波特征，直到计算完成训练样本集Y中所有样本的小波特征，得到小波特征数据集Y′。

3.根据权利要求1所述的字典学习和稀疏表示的SAR图像分割方法，其中中步骤(4)所述的用谱聚类算法对小波特征数据集Y′进行聚类，按如下步骤进行：

(4a)根据下式计算小波特征数据集Y′中的特征样本y′_i和y′_j之间的相似度值W_ij，构建相似度矩阵W：

$W_{\mathrm{ij}}=\exp (-\frac{{\left|\right|{y^{\′}}_i-{y^{\′}}_j\left|\right|}^2}{{2\mathrm{\σ}}^2}),i,j=\mathrm{1,2,3},...,n;$

W＝[W_ij]，i，j＝1，2，3，...n；

其中y′_i∈Y′，y′_i为Y′中的第i个样本，y′_j∈Y′，y′_j为Y′中的第j个样本，σ为尺度参数；

(4b)利用相似度矩阵W计算拉普拉斯矩阵L：

L＝Λ^-1/2WΛ^-1/2；

其中Λ为对角矩阵，其对角线元素

(4c)对拉普拉斯矩阵L进行矩阵特征值分解，得到拉普拉斯矩阵L最大的前k个特征值λ₁≥λ₂≥...≥λ_k和与这些特征值对应的k个特征向量v¹，v²，...，v^k，利用k个特征向量构建特征向量矩阵U＝[v¹，v²，...，v^k]；

(4d)根据下式规范化特征向量矩阵U，得到规范化后的特征向量矩阵U′：

U′＝[U′_ij]，i＝1，2，3，...n，j＝1，2，3，...，k；

其中

U_ij为特征向量矩阵U中第i行第j列的元素；

(4e)用k-means算法对规范化后的特征向量矩阵U′按行聚类，得到对应小波特征数据集Y′的划分，将小波特征数据集Y′的划分结果作为训练样本集Y的划分结果，得到Y_i，i＝1，2，3，...k。

4.根据权利要求1所述的字典学习和稀疏表示的SAR图像分割方法，其中中步骤(5)所述的对划分后的每一类训练数据利用K-SVD算法训练一个字典得到D_i，i＝1，...，k，按如下步骤进行：

(5a)取训练样本集合中的第j类数据Y_i，j初值为1，令Y_T＝Y_j；

(5b)根据下式利用正交匹配追踪OMP方法，完成对Y_T的稀疏编码：

Y_T＝DX s.t.|x_i|₀≤t₀

其中，Y_T为需要稀疏表示的样本集合，D为用来稀疏表示的字典，X为样本经字典稀疏表示后的得到稀疏矩阵；x_i为稀疏矩阵X的第i列，t₀为预设定稀疏向量中的非零元素个数，称为稀疏度；

(5c)对字典D中所有原子d_i i＝1，2，3...L进行如下更新，原子d_i指字典中的第i个列向量，L为字典中原子个数：

首先，从样本集合Y_T中找出使用原子d_i的所有样本，取样本序号，得到具有n_ω个元素的序号集合ω：ω＝{k1|≤k≤n，xⁱ(k)≠0}，其中xⁱ为对样本集合Y_T稀疏编码后得到的稀疏矩阵X的第i行；

接着，令稀疏矩阵X的第i行元素全部为零，得到新矩阵X′，计算字典D与新矩阵X′对样本集合Y_T的重构误差矩阵E_i：

E_i＝Y_T-DX′；

然后，依据序号集合ω中所有序号，取所有n_ω个序号对应的误差矩阵E_i中的n_ω个列向量，构成特异误差矩阵

最后，应用公式

对特异误差矩阵

进行矩阵奇异值分解，并用分解行矩阵U的第一列替换所要更新的原子d_i，并用对角矩阵Δ中的第一行第一列的元素Δ(1，1)乘以分解列矩阵V的第一个列向量得到一个新的向量，利用该新的向量替换xⁱ，xⁱ为稀疏矩阵X的第i行；

(5d)判断是否满足j≥k，若满足，则得到字典D_i，i＝1，...，k，否则j＝j+1，返回步骤(5a)计算下一类训练数据对应字典；

(5e)重复步骤步骤(5a)-(5d)，直到得到所有类别的训练数据的对应字典D_i，i＝1，...，k。

Images