CN102968796A - 基于采样学习的sar图像分割方法 - Google Patents

Abstract

本发明公开一种基于采样学习的SAR图像分割方法，主要解决现有算法计算量大，分割速度慢的问题。其实现步骤为：(1)输入待分割图像，提取特征；(2)对数据集进行M次随机采样；(3)采用谱聚类算法对M次样本数据集分别进行聚类；(4)将M次聚类后属于同一类的数据合在一起，每一类集合在一起的数据构成该类的新数据集，利用KSVD算法对新数据集训练一个字典；(5)求解测试样本在字典上的稀疏编码；(6)计算测试样本在字典上的重构误差；(7)利用重构误差确定测试样本的标签，得到最终的分割结果。本发明具有分割快速准确的优点，该方法进一步可用于SAR图像目标识别与分类。

Description

基于采样学习的SAR图像分割方法

技术领域

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

背景技术

随着科学技术的发展，人们越来越多的以图像的形式表示各种信息。图像分割也成为人们研究的热点。

合成孔径雷达SAR不受气候、昼夜等因素影响，具有全天候成像优点。它利用合成孔径原理提高方位分辨率，利用脉冲压缩技术获得高的距离分辨率，从而相比真实孔径雷达在遥感领域具有很大优势。由于SAR具有的独特作用，使得SAR图像的理解和解译，在国防和民用领域正受到越来越广泛的关注。

目前已经有很多成熟的聚类算法被用到SAR图像分割中。但是基于聚类的SAR图像分割方法例如谱聚类算法，进行计算时要以像素点为单位，其中用到的方阵的维数为图像中所有像素点个数，这就导致需要处理的数据维数非常大，使得分割时间很长，在一般的计算机上实现是比较困难的，从而也使该算法的应用受到限制。

发明内容

本发明的目的在于针对现有技术的不足，提出一种基于采样学习的SAR图像分割方法，以减小数据量，提高分割速度，易于在普通计算机上实现。

本发明的技术方案是：对待分割的图像提取特征，再进行M次随即采样，对每次采样的数据进行聚类，将M次聚类后属于一类的数据集合在一起训练得到一个字典，通过稀疏编码获得所有测试样本的标签。实现步骤包括如下：

（1）对待分割的SAR图像，提取三层小波特征和灰度共生特征，得到规模为26×N的数据集Y，N为图像中所有像素点的个数；

（2）对数据集Y进行M次随机采样，每次选择规模为26×n的样本数据集Y^P作为样本，n<<N，p=1,2,…,M；

（3）采用谱聚类算法对每次随机采样的样本数据集Y^P进行聚类，得到K个子数据集，并且再从每一个子数据集中选取一些靠近该子数据集聚类中心的数据，K为所要分割的类数；

（4）将M次样本聚类后属于同一类的数据合在一起，对每一类集合在一起的数据构成该类的新数据集Y_k，利用KSVD算法对新数据集Y_k训练一个字典D_k，k=1,2,3,…,K；

（5）利用正交匹配追踪OMP算法求解测试样本S在字典D=[D₁，...,D_K]上的稀疏编码，得到稀疏向量α=[(α¹)^T,...,(α^K)^T]^T，[]^T表示矩阵的转置；

（6）根据子字典D_k和稀疏子向量α^k计算测试样本S在字典D_k上的重构误差： $R_k\left(S\right)={\left|\right|S-D_k{\mathrm{\&alpha;}}^k\left|\right|}_2^2,$

表示l₂范数；

（7）利用重构误差R_k(S)，求解测试样本S的标签l(S)=argmin_{k＝1，...，K} R_k(S)，arg min为求最小值的函数，得到每个测试样本的标签，进而得到最终的分割结果。

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

1、本发明使用少量的数据来学习字典，有效的解决了现有算法运算量大，分割速度慢的问题。

2、本发明使用的多次采样数据集合在一起的方法，能获得更好的用于学习字典的数据。

3、本发明使用的字典学习和稀疏编码的方法，能充分发挥SAR图像内在的信息，获得更好的分割效果。

附图说明

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

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

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

具体实施方式

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

步骤1，对待分割的SAR图像，提取三层小波特征和灰度共生特征，得到规模为26×N的数据集Y。

(1a)对原图进行三层平稳小波变换，得到系数矩阵coef_m1(i₁,j₁),m1=1,...,10，当m1=1时，coef_m1(i₁,j₁)代表低频系数，当m1>1时，coef_m1(i₁,j₁)代表高频系数，对每个像素点提取10维子带能量特征e(i,j)=[e₁(i,j),...,e₁₀(i,j)]^T，作为该像素点的小波特征：

$e_{m1}(i,j)=\frac{1}{w\&times;w}\underset{i_1=i-(w-1)/2}{\overset{i+(w-1)/j}{\mathrm{\&Sigma;}}}\underset{j_1=j-(w-1)/2}{\overset{j+(w-1)/2}{\mathrm{\&Sigma;}}}\left|{\mathrm{coef}}_{m1}(i_1,j_1)\right|,$

其中w×w为滑动窗口的大小，coef_m1(i₁,j₁)为平稳小波子带中第i₁行第j₁列对应的系数值；

(1b)计算特征窗口内的图像子块的灰度共生矩阵P_d,θ，u,v(i3,j3)，其中，u和v表示灰度共生矩阵计算中两个不同像素点的灰度值；d表示灰度共生矩阵计算中两个像素点的距离，取值为1；θ表示提取灰度共生矩阵的方向：

(1b1)选择的θ的方向为0°，45°，90°，135°；

(1b2)分别在θ的四个方向上计算灰度共生矩阵的能量，熵，同质性，对比度，总共16维特征：

g(i3,j3)=[g¹(i3,j3),g²(i3,j3),g³(i3,j3),g⁴(i3,j3)]^T，

能量： $g^1(i3,j3)=[g_0^1,g_{45}^1,g_{90}^1,g_{135}^1],$

$g_{\mathrm{\&theta;}}^1(i3,j3)={\mathrm{\&Sigma;}}_{i3=1}^B{\mathrm{\&Sigma;}}_{j3=1}^B{\left(P_{d,\mathrm{\&theta;},u,v}(i3,j3)\right)}^2;$

熵： $g^2(i3,j3)=[g_0^2,g_{45}^2,g_{90}^2,g_{135}^2],$

$g_{\mathrm{\&theta;}}^2(i3,j3)=-{\mathrm{\&Sigma;}}_{i3=1}^B{\mathrm{\&Sigma;}}_{j3=1}^B{P}_{d,\mathrm{\&theta;},u,v}(i3,j3)\log P_{d,\mathrm{\&theta;},u,v}(i3,j3);$

同质性： $g^3(i3,j3)=[g_0^3,g_{45}^3,g_{90}^3,g_{135}^3],$

$g_{\mathrm{\&theta;}}^3(i3,j3)={\mathrm{\&Sigma;}}_{i3}^B{\mathrm{\&Sigma;}}_{j3=1}^B{P}_{d,\mathrm{\&theta;},u,v}(i3,j3)/(1+|i3-j3\left|\right);$

对比度： $g^4(i3,j3)=[g_0^4,g_{45}^4,g_{90}^4,g_{135}^4],$

$g_{\mathrm{\&theta;}}^4(i3,j3)={\mathrm{\&Sigma;}}_{i3=1}^B{\mathrm{\&Sigma;}}_{j3=1}^B{(i3-j3)}^2{P}_{d,\mathrm{\&theta;},u,v}(i3,j3);$

其中，B为量化的灰度等级，log表示求对数函数；

将10维小波特征和16维灰度共生特征相结合Y_i，j=[e(i,j)^T,g(i,j)^T]^T，为了处理方便将所有像素点排成一行，从而就得到规模为26×N的数据集Y=[Y₁,...,Y_(i-1)×n2+j1,...,Y_N]，i=1,...,m2,j=1,...,n2，其中m2和n2为图像的行数和列数，N=m2×n2为图像中所有像素点的个数。

步骤2，对数据集Y，进行M次随机采样，每次选择规模为26×n的样本数据集Y^P作为样本，n<<N，p=1,2,…,M，即随机生成一个大小为N的向量，该向量中每一个数对应于数据集Y中的一个数据，取该向量前n个数对应于数据集Y中的数据，得到样本数据集Y^P。

步骤3，采用谱聚类算法对每次随机采样的样本数据集Y^P进行聚类，得到K个子数据集，并且再从每一个子数据集中选取一些靠近该子数据集聚类中心的数据。

(3a)计算任意两点之间的相似度W_i2j2，构建相似度矩阵W，W=[W_i2j2]：

$W_{i2j2}=\exp (-\frac{{\left|\right|y_{i2}-y_{j2}\left|\right|}^2}{2{\mathrm{\&sigma;}}^2}),i2,j2=\mathrm{1,2,3},...,n;$

其中y_i2和y_j2为样本数据集Y^P中的数据，σ为尺度参数，exp为指数函数；

(3b)利用相似度矩阵W计算拉普拉斯矩阵L：L=I-A^-1/2WA^-1/2；

其中I为单位矩阵，A为度矩阵且只有主对角线上的元素非零，A=[a_i2i2]， $a_{i2i2}={\mathrm{\&Sigma;}}_{j2=1}^n{W}_{i2j2};$

(3c)对拉普拉斯矩阵L进行特征值分解，得到特征值和其对应的特征向量，利用前K个最小的特征值对应的特征向量构造特征向量矩阵T，即从分解得到的特征向量中选出K个组成新的特征向量矩阵T，用K-means算法对特征向量矩阵T进行聚类，得到样本的聚类结果；

(3d)从每一个子数据集中选取一些靠近该子数据集聚类中心的数据。

步骤4，将M次样本聚类后属于同一类的数据合在一起，对每一类集合在一起的数据构成该类的新数据集Y_k，利用KSVD算法对新数据集Y_k训练一个字典D_k。

(4a)将M次样本聚类后属于同一类的数据合在一起，构成该类的新数据集Y_k；

(4b)随机初始化一个字典矩阵

(4c)利用正交匹配追踪OMP算法，完成对每一类集合在一起的新数据集Y_k的稀疏编码：

$X=\underset{X^{\&prime;\&prime;}}{\min }\left|\right|Y_k-D_k^0{X}^{\&prime;\&prime;}\left|\right|,s.t.{\left|x_i\right|}_0\&le;t_0;$

其中，X为样本经稀疏编码后得到的稀疏矩阵，X″为求解过程中对X的近似解，s.t.表示约束条件，x_i为稀疏矩阵的第i列，t₀为稀疏向量的非零元素的个数，称为稀疏度，||₀表示l₀范数；

(4d)对字典矩阵

中所有列d_i进行如下更新：

从每一类集合在一起的数据Y_k中找出使用到d_i的所有样本序号，得到大小为n_ω的序号集合ω：

ω={h｜xⁱ(h)≠0}，1≤h≤n；

其中，h表示每一类集合在一起的数据Y_k中使用到d_i的样本序号，xⁱ为稀疏矩阵X的第i行；

令稀疏矩阵X的第i行元素全部为零，得到新矩阵X′，利用字典矩阵

与新矩阵X′得到每一类集合在一起的数据Y_k的重构误差矩阵E_i：

$E_i=Y_k-D_k^0{X}^{\&prime;};$

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

再对特异误差矩阵

进行奇异值分解，得到它的三个分量U，Δ，V；

$E_i^R={\mathrm{U\&Delta;V}}^T;$

其中，U为对特异误差矩阵

进行奇异值分解后得到的分解行矩阵，Δ为奇异值分解后的对角矩阵，V为奇异值分解后的分解列矩阵；

用分解行矩阵U的第一列替换所要更新的列d_i，用对角矩阵Δ中的第一个的元素乘以分解列矩阵V的第一个列向量得到一个新的向量，用该向量替换稀疏矩阵X的第i行xⁱ；

(4e)重复(4c)到(4d)的步骤，得到字典D_k。

步骤5，利用正交匹配追踪OMP算法求解测试样本S在字典D=[D₁，...，D_K]上的稀疏编码，得到稀疏向量α=[(α₁)^T，...,(α^K)^T]^T。

步骤6，根据子字典D_k和稀疏子向量α^k计算测试样本S在字典D_k上的重构误差： $R_k\left(S\right)={\left|\right|S-D_k{\mathrm{\&alpha;}}^k\left|\right|}_2^2,$

表示l₂范数。

步骤7，利用重构误差R_k(S)，求解测试样本S的标签l(S)=arg min_{k＝1，...，K} R_k(S)，arg min为求最小值的函数，依次得到每个测试样本的标签，进而得到最终的分割结果。

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

1．实验条件设置

本发明中所用SAR图像为二类图和三类图，并且分别在这2副SAR图像上进行了实验。

实验中采样的样本数n为200，采样次数M取5，从每一个子数据集中选取靠近该子数据集聚类中心的20个数据，每一类字典的个数为50，小波特征窗口为15×15，灰度共生特征窗口为16×16，谱聚类中的尺度参数σ=2，稀疏度为3，谱聚类选取的样本点数为100。

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

三类待分割SAR图像的大小为256×256，它包括山地河流，森林和平原地区三个部分，如图3(a)所示。

2．仿真内容及结果

仿真1：实验中用K-means聚类算法和

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

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

从图2(b)可以看出，K-means方法分割的结果虽然分割出了一些纵向延伸的山脊，但边缘部分错分很多。

从图2(c)可以看出，

谱聚类方法分割的结果比K-means方法错分的部分要少一些，但对于一些地方还是有很多错分。

从图2(d)可以看出，本方法的分割的结果能分出更多的岔道，而且在边缘上，也要比以上2种方法好。

仿真2：实验中用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图像，提取三层小波特征和灰度共生特征，得到规模为26×N的数据集Y，N为图像中所有像素点的个数；

（2）对数据集Y进行M次随机采样，每次选择规模为26×n的样本数据集Y^P作为样本，n<<N，p=1,2,…,M；

（3）采用谱聚类算法对每次随机采样的样本数据集Y^P进行聚类，得到K个子数据集，并且再从每一个子数据集中选取一些靠近该子数据集聚类中心的数据，K为所要分割的类数；

（4）将M次样本聚类后属于同一类的数据合在一起，对每一类集合在一起的数据构成该类的新数据集Y_k，利用KSVD算法对新数据集Y_k训练一个字典D_k，k=1,2,3,…,K；

（5）利用正交匹配追踪OMP算法求解测试样本S在字典D=[D₁，...,D_K]上的稀疏编码，得到稀疏向量α＝[(α¹）^T,...,(α^K)^T]^T，[]^T表示矩阵的转置；

（6）根据子字典D_k和稀疏子向量α^k计算测试样本S在字典D_k上的重构误差：

$R_k\left(S\right)={\left|\right|S-D_k{\mathrm{\&alpha;}}^k\left|\right|}_2^2,$

表示l₂范数；

（7）利用重构误差R_k(S)，求解测试样本S的标签l(S)＝argmin_{k＝1，...,K} R_k(S)，arg min为求最小值的函数，得到每个测试样本的标签，进而得到最终的分割结果。

2.根据权利要求1所述的基于采样学习的SAR图像分割方法，其中步骤（1）所述的提取三层小波特征和灰度共生特征，按如下步骤进行：

(2a)对原图进行三层平稳小波变换，得到系数矩阵coef_m1(i₁,j₁),m1＝1,...,10，当m1=1时，coef_m1(i₁,j₁)代表低频系数；当m1>1时，coef_m1(i₁,j₁)代表高频系数。对每个像素点提取10维子带能量特征e(i,j)＝[e₁(i,j),...,e₁₀(i,j)]^T，作为该像素点的小波特征：

$e_{m1}(i,j)=\frac{1}{w\&times;w}\underset{i_1=i-(w-1)/2}{\overset{i+(w-1)/j}{\mathrm{\&Sigma;}}}\underset{j_1=j-(w-1)/2}{\overset{j+(w-1)/2}{\mathrm{\&Sigma;}}}\left|{\mathrm{coef}}_{m1}(i_1,j_1)\right|,$

其中w×w为滑动窗口的大小，coef_m1(i₁,j₁)为平稳小波子带中第i₁行第j₁列对应的系数值；

(2b)通过灰度共生矩阵，计算角度为0°，45°，90°，135°四个方向的能量，熵，同质性，对比度，得到16维特征作为灰度共生特征。

3.根据权利要求1所述的基于采样学习的SAR图像分割方法，其中步骤（3）所述的采用谱聚类算法对每次随机采样的样本数据集YP进行聚类，按如下步骤进行：

(3a)构建相似度矩阵W；

(3b)利用相似度矩阵W计算拉普拉斯矩阵L：L＝I-A^-1/2WA^-1/2；

其中I为单位矩阵，A为度矩阵且只有主对角线上的元素非零，A=[a_i2i2]，

$a_{i2i2}={\mathrm{\&Sigma;}}_{j2=1}^n{W}_{i2j2}$

(3c)对拉普拉斯矩阵L进行特征值分解，得到特征值和其对应的特征向量，利用前K个最小的特征值对应的特征向量构造特征向量矩阵T，即从分解得到的特征向量中选出K个组成新的特征向量矩阵T，用K-means算法对特征向量矩阵T进行聚类，得到样本的聚类结果。

4.根据权利要求1所述的基于采样学习的SAR图像分割方法，其中步骤（4）所述的对每一类集合在一起的数据构成该类的新数据集Y_k，利用KSVD算法对新数据集Y_k训练一个字典D_k，按如下步骤进行：

(4a)随机初始化一个字典矩阵

(4b)利用正交匹配追踪OMP算法，对集合在一起的新数据集Y_k进行稀疏编码；

$X=\underset{X^{\&prime;\&prime;}}{\min }\left|\right|Y_k-D_k^0{X}^{\&prime;\&prime;}\left|\right|,s.t.{\left|x_i\right|}_0\&le;t_0;$

其中，X为样本经稀疏编码后得到的稀疏矩阵，X″为求解过程中对X的近似解，s.t.表示约束条件，x_i为稀疏矩阵的第i列，t₀为稀疏向量中非零元素的个数，称为稀疏度，‖₀表示l₀范数；

(4c)对字典矩阵

中所有列d_i进行如下更新：

从每一类集合在一起的数据Y_k中找出使用到d_i的所有样本序号，得到大小为n_ω的序号集合ω：

ω={h｜xⁱ(h)≠0}，1≤h≤n；

其中，h表示每一类集合在一起的数据Y_k中使用到d_i的样本序号，xⁱ为稀疏矩阵X的第i行；

令稀疏矩阵X第i行的元素全部为零，得到新矩阵X′，利用字典矩阵

与新矩阵X′得到每一类集合在一起的数据Y_k的重构误差矩阵E_i：

$E_i=Y_k-D_k^0{X}^{\&prime;};$

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

再对特异误差矩阵

进行奇异值分解，得到它的三个分量U，Δ，V；

$E_i^R={\mathrm{U\&Delta;V}}^T;$

其中，U为对特异误差矩阵进行奇异值分解后得到的分解行矩阵，Δ为奇异值分解后的对角矩阵，V为奇异值分解后的分解列矩阵；

用分解行矩阵U的第一列替换所要更新的列d_i，用对角矩阵Δ中的第一个的元素乘以分解列矩阵V的第一个列向量得到一个新的向量，用该向量替换稀疏矩阵X的第i行xⁱ；

(4d)重复(4b)到(4d)的步骤，得到字典D_k。

Images