CN101853491B - 基于并行稀疏谱聚类的sar图像分割方法 - Google Patents

Abstract

本发明公开了一种基于并行稀疏谱聚类的SAR图像分割方法，它涉及图像处理技术领域，主要解决已有谱聚类技术在解决大规模SAR图像分割应用上的局限性问题。其实现过程是：1)对待分割的SAR图像提取特征；2)配置MATLAB并行计算环境；3)将任务划分到各个处理器节点，计算分块稀疏相似矩阵；4)并行任务调度器收集计算结果并汇聚成完整的稀疏相似矩阵；5)求解拉普拉斯矩阵并进行特征分解；6)对经过归一化后的特征向量矩阵进行K-均值聚类；7)输出SAR图像分割结果。本发明能够有效克服已有谱聚类技术在运算量和存储空间上的瓶颈问题，对大规模SAR图像进行分割效果显著，适用于SAR图像目标检测和目标识别。

Description

基于并行稀疏谱聚类的SAR图像分割方法

技术领域

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

背景技术

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

其中，基于聚类的SAR图像分割方法，是将SAR图像中相似的区域尽量划分成一类，将不相似的区域划分在不同的类别中。已经有很多成熟的聚类算法被用到SAR图像分割中。作为一门新兴学科，结合谱图理论的谱聚类算法充分发掘了数据的成对点相似特性，通过图的拉普拉斯矩阵特征分解达到降维目的，使得聚类结果更加稳健。但是，当数据规模n比较大的时候，谱聚类算法需要计算大小为n×n的相似矩阵，并计算相应拉普拉斯矩阵的特征分解问题，其计算的时间复杂度和空间复杂度分别为O(n³)和O(n²)，不能有效地对包含大规模数据的SAR图像进行分割。

发明内容

本发明的目的在于克服上述已有技术的不足，提出了一种基于并行稀疏谱聚类的SAR图像分割方法，以减小谱聚类的时间和空间复杂度，从而能够有效地对包含大规模数据的SAR图像进行分割。

为实现上述目的，本发明包括以下步骤：

(1)对待分割的SAR图像进行3层平稳小波变换，图像像素点总个数为n，通过下式对每个像素点提取出10维子带能量特征，构成大小为n×10的输入数据样本E：

$E=\frac{1}{M\×N}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}+\left|\mathrm{coef}(i,j)\right|$

其中，M×N为利用滑动窗口确定的子带大小，这里取值为16×16，coef(i，j)为平稳小波子带中第i行第j列的系数值；

(2)配置MATLAB7.8R2009(a)版本的并行计算环境；

(3)进行并行任务划分，使用并行任务调度器将像素点总个数为n的输入数据样本均衡划分到p个处理器上，遇到不能整除的情况则将剩余的数据样本存放到最后一个节点上；

(4)在各个处理器节点上分布式地计算n/p个数据样本到其他所有数据样本之间的欧氏距离，得到距离矩阵；

(5)对距离矩阵的每一列按照从小到大顺序进行排序，并根据下面的公式计算分块相似矩阵S：

$S=\exp (-\frac{{\left|\right|x_i-x_j\left|\right|}^2}{{2\mathrm{\σ}}^2}),i,j=1,...,n$

其中，x₁，...，x_n表示n个数据样本点，||x_i-x_j||²表示排序后的距离矩阵，σ表示尺度参数，对相似矩阵只保留每个数据样本与周围离它最相近的t个数据样本之间的值，其他位置都稀疏化成0值；

(6)将各个处理器节点计算得到的相似矩阵通过并行任务调度器收集到主节点上，并排成一行汇聚成一个大小为n×n的稀疏相似矩阵W；

(7)根据稀疏相似矩阵W计算对角线元素为

的度矩阵D，进而计算出拉普拉斯矩阵L＝D^-1/2WD^-1/2，利用MATLAB的eigs函数对拉普拉斯矩阵进行特征分解；

(8)将特征分解后得到的特征值进行排序，得到前K个最大特征值1＝λ₁≥λ₂≥…≥λ_K，其对应的特征向量表示成v¹，v²，…，v^K；

(9)将得到的特征向量排成一行，构造出特征向量矩阵U＝[v¹，v²，…，v^K]，并对其归一化，得到规范化特征向量矩阵

(10)将规范化特征向量矩阵Y的每一行看成是空间□^K中的一个点，用K-均值聚类算法将其聚为K类，并使用正交初始化寻找初始聚类中心；

(11)根据得到的聚类中心，使用欧氏距离测度把SAR图像中的所有像素点划分到不同的类别中，每一个像素点获得一个聚类标签，将这些聚类标签重新组合成和原图像尺寸相同的矩阵，将矩阵元素作为分割结果图像的灰度值并显示出来，从而得到最终分割结果。

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

1)本发明利用MATLAB并行计算环境，设计用于SAR图像分割的并行稀疏谱聚类方法，能够有效解决现有谱聚类技术运算量过大的问题；

2)本发明利用稀疏矩阵存储谱聚类的大规模相似矩阵，能够有效解决现有谱聚类技术在存储空间上的局限性；

3)本发明相对已有的SAR图像分割技术，设计的并行稀疏谱聚类方法能够充分发掘SAR图像的内在信息，获得更好的分割效果：

仿真实验结果表明，本发明提出的并行稀疏谱聚类方法能够有效地应用于SAR图像分割，并进一步应用于雷达目标检测和目标识别。

附图说明

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

图2是本发明的MATLAB并行计算环境示意图；

图3是本发明仿真使用的集群环境示意图；

图4是本发明与现有K-均值聚类算法和近似谱聚类算法应用于ChinaLake机场SAR图像分割的对比实验结果图；

图5是本发明与现有K-均值聚类算法和

近似谱聚类算法应用于厦门鼓浪屿SAR图像分割的对比实验结果图；

图6是本发明与现有K-均值聚类算法和

近似谱聚类算法应用于瑞士湖区SAR图像分割的对比实验结果图；

图7是本发明对瑞士湖区SAR图像分割的并行加速比和效率示意图。

具体实施方式

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

步骤1.提取待分割SAR图像的特征。

SAR图像不仅数据量大，在成像过程中不同的地物有不同的后向发射和散射特性，从而具有丰富的幅度、相位、极化和纹理信息，并且，图像固有的相干斑噪声对于分割性能直接产生影响，因此，有必要在图像分割前对SAR图像进行纹理分析，提取有效的纹理特征进行聚类。

在以上分析基础上，先对待分割的SAR图像进行3层平稳小波变换，图像像素点总个数为n，通过下式对每个像素点提取出10维子带能量特征，构成大小为n×10的输入数据样本E：

$E=\frac{1}{M\×N}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\left|\mathrm{coef}(i,j)\right|$

其中，M×N为利用滑动窗口确定的子带大小，这里取值为16×16，coef(i，j)为平稳小波子带中第i行第j列的系数值，这样，对于大小为256×256的SAR图像，就构成65536×10的矩阵作为并行稀疏谱聚类方法的输入数据。

步骤2.配置MATLAB并行计算环境。

本步骤的具体实施参照图2如下：

2a)在如图3所示的高性能64刀片集群环境中安装MATLAB 7.8R2009(a)版本的并行计算工具箱和分布式计算服务器；

2b)在各个集群节点上用mdce start命令启动分布式计算服务器；

2c)在集群的主节点上用startjobmanager命令启动并行任务调度器；

2d)在各个处理器节点上用startworker命令启动工作进程。

步骤3.并行任务划分。

3a)根据配置好的MATLAB并行计算环境，使用findResource命令指定配置好的任务调度器；

3b)使用该任务调度器，用createJob命令创建作业；

3c)在所创建的作业中，将规模为n的输入数据样本用createTask命令均衡划分到p个处理器上，遇到不能整除的情况则将剩余的数据样本存放到最后一个节点上；

3d)使用submit命令提交作业，然后用waitForState命令将各个任务安排到相应的任务队列里等待各个处理器节点的执行；

3e)各个处理器节点分布式地计算n/p个数据样本到其他所有数据样本之间的欧氏距离，获得距离矩阵，对得到的距离矩阵的每一列按照从小到大顺序进行排序，并根据下面的公式将其转换成分块相似矩阵S：

$S=\exp (-\frac{{\left|\right|x_i-x_j\left|\right|}^2}{{2\mathrm{\σ}}^2}),$ i，j＝1， ...，n

其中，x₁，...，x_n表示n个数据样本点，||x_i-x_j||²表示排序后的距离矩阵，σ表示尺度参数，对相似矩阵只保留每个数据样本与周围离它最相近的t个数据样本之间的值，其他位置都稀疏化成0值。

步骤4.主节点任务收集。

用getAllOutputArguments命令将p个处理器节点的计算结果收集到主节点上，并汇聚成一个大小为n×n的完整的稀疏相似矩阵W＝[S₁，...，S_p]，其中S_i(i＝1，...，p)表示p个处理器节点计算得到的分块相似矩阵。

步骤5.计算拉普拉斯矩阵特征分解问题。

根据这个稀疏相似矩阵W计算对角线元素为度矩阵D，进而计算出拉普拉斯矩阵L＝D^-1/2WD^-1/2，利用MATLAB的eigs函数对这个大规模稀疏拉普拉斯矩阵进行特征分解，得到前K个最大特征值1＝λ₁≥λ₂≥…≥λ_K，以及这些特征值对应的特征向量v¹，v²，…，v^K。

步骤6.对归一化的特征向量矩阵进行K-均值聚类。

将得到的特征向量排成一行，构造出特征向量矩阵U＝[v¹，v²，…，v^K]，对其归一化得到规范化特征向量矩阵

将规范化特征向量矩阵Y的每一行看成是空间□^K中的一个点，用K-均值聚类算法将其聚为K类，并使用正交初始化寻找初始聚类中心。

步骤7.输出分割结果。

根据得到的聚类中心，使用欧氏距离测度把SAR图像中的所有像素点划分到不同的类别中，每一个像素点获得一个聚类标签，将这些聚类标签重新组合成和原图像尺寸相同的矩阵，将矩阵元素作为分割结果图像的灰度值并显示出来，从而得到最终分割结果。

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

对比实验中给出本发明的并行稀疏谱聚类方法，以及现有的K-均值聚类算法和近似谱聚类算法的分割结果。其中，现有的K-均值聚类和

近似谱聚类选取的是10次实验视觉效果最好的分割结果，本发明的并行稀疏谱聚类方法由于结果比较稳定选取的是1次实验的分割结果。

近似谱聚类选取的样本点个数为100，尺度参数s＝0.2。并行稀疏谱聚类选取的最近邻数为100，尺度参数s＝0.05。

使用以上三种分割方法对3幅SAR图像进行分割的对比实验结果如下：

1)对Ku波段SAR图像的分割结果如图4所示，其中：

图4(a)为美国加利福尼亚州的China Lake机场3米分辨率的Ku波段SAR图像，图像尺寸为400×400，包括跑道、空地和机场建筑物三类地物。

图4(b)为现有的K-均值聚类方法的分割结果，可见跑道上有很多细小的斑点被误分成其他两类。

图4(c)为现有的

近似方法的分割结果，尽管机场跑道的区域一致性得到了改善，但是右上方的机场建筑物则被分割得支离破碎，在实际应用中不利于准确进行目标识别。

图4(d)为本发明的并行稀疏谱聚类方法的分割结果，它准确给出了完整的跑道和建筑物的轮廓，并且，对于三岔形跑道中央密布的路网，相比其他两种算法获得了更加完整的区域。

2)对C-波段SAR图像的分割结果如图所示，其中：

图5(a)为一幅RadarSAT-2卫星拍摄的中国厦门鼓浪屿地区C-波段SAR图像，极化方式为HV极化，截取图像块尺寸为400×400。

图5(b)为现有的K-均值聚类方法的分割结果，它受岛屿上地形分布干扰太大，造成陆地部分很多误分，延伸进水域的一些陆地与岛屿之间被零散地分隔开来。

图5(c)为现有的

近似方法的分割结果，尽管比K-均值的结果有些改善，但还是未能获得完整的一块区域。

图5(d)是本发明的并行稀疏谱聚类方法的分割结果，可以清晰地看出鼓浪屿和厦门岛一样被完整划分成陆地部分，并且对于一些与水域对比度差异较小的海岸线也能和岛屿部分保持较好的连接度。

3)对X-波段SAR图像的分割结果如图6所示，其中：

图6(a)为Space Radar Laboratory Missions的一幅1994年瑞士一湖区的3视X-波段SAR图像，图像尺寸为512×512，包含三类地物：湖泊、城区和山地。

图6(b)为现有的K-均值聚类方法的分割结果，它很好地将湖泊划分出来，但是，在左右两边的山地地区同时也存在了严重的错分成水域的情况。

图6(c)为现有的近似方法的分割结果，它对山地区域的分割情况有所改善，但是城区和山地有很多的混淆。

图6(d)为本发明的并行稀疏谱聚类方法的分割结果，它在三种算法中最佳，它不仅将湖泊沿岸准确提取出来，在山地区域没有丝毫误分成湖泊的情况，并且，城区的分割结果保持了较好的区域完整性。因此，本发明的并行稀疏谱聚类方法被用来分割这幅SAR图像是十分合适的。

在本发明的并行计算中，对于p个处理器，T_S表示在串行计算机上算法从开始到结束所需要的时间，T_P表示在并行计算机上算法从开始到最后一个处理器完成任务所需要的时间，加速比S_P和效率E_P分别定义为S_P＝T_S/T_P，E_P＝S_S/P。加速比衡量的是并行算法的执行速度相对于串行算法的执行速度加快了多少倍。效率衡量的是单个处理器的计算能力被有效利用的比率。

用本发明的并行稀疏谱聚类方法对图6(a)所示的瑞士湖区SAR图像进行分割，其并行加速比曲线和并行效率曲线分别如图7(a)和图7(b)所示。其中，所用到的处理器个数最多达到了42个。

从图7可以看出，由于本发明的并行稀疏谱聚类方法采用MATLAB分布式计算，除了任务的调度和数据共享的传输之外，各个处理器之间并不需要进行通信，因此在处理器个数小于24时可以得到近乎线性的加速比。而当处理器个数上升到42个时，由于大规模稀疏矩阵传输以及多处理器任务调度耗时增加，使得加速比有所下降，但是还是保持着上升的趋势。在效率方面，小于32个处理器的并行情况下，效率都可以达到90％以上，表明各个处理器节点得到了相对充分的利用。处理器较多时，消息阻塞和消息延时降低了单个处理器的利用率。

Claims (3)

1.一种基于并行稀疏谱聚类的SAR图像分割方法，包括以下步骤：

(1)对待分割的SAR图像进行3层平稳小波变换，图像像素点总个数为n，通过下式对每个像素点提取出10维子带能量特征，构成大小为n×10的输入数据样本中特征的计算公式E：

$E=\frac{1}{M\×N}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\left|\mathrm{coef}(i,j)\right|$

其中，M×N为利用滑动窗口确定的子带大小，这里取值为16×16，coef(i，j)为平稳小波子带中第i行第j列的系数值；

(2)配置MATLAB7.8 R2009(a)版本的并行计算环境；

(3)进行并行任务划分，使用并行任务调度器将像素点总个数为n的输入数据样本均衡划分到p个处理器上，遇到不能整除的情况则将剩余的数据样本存放到最后一个节点上；

(4)在各个处理器节点上分布式地计算n/p个数据样本到其他所有数据样本之间的欧氏距离，得到距离矩阵；

(5)对距离矩阵的每一列按照从小到大顺序进行排序，并根据下面的公式计算分块相似矩阵S：

$S=\exp (-\frac{{\left|\right|x_i-x_j\left|\right|}^2}{2{\mathrm{\σ}}^2}),i,j=1,...,n$

其中，x₁，…，x_n表示n个数据样本点，‖x_i-x_j‖²表示排序后的距离矩阵，σ表示尺度参数，对相似矩阵只保留每个数据样本与周围离它最相近的t个数据样本之间的值，其他位置都稀疏化成0值；

(6)将各个处理器节点计算得到的相似矩阵通过并行任务调度器收集到主节点上，并排成一行汇聚成一个大小为n×n的稀疏相似矩阵W；

(7)根据稀疏相似矩阵W计算对角线元素为

的度矩阵D，进而计算出拉普拉斯矩阵L＝D^-1/2WD^-1/2，利用MATLAB的eigs函数对该拉普拉斯矩阵进行特征分解；

(8)对特征分解后得到的特征值进行排序，得到前K个最大特征值1＝λ₁≥λ₂≥…≥λ_K，其对应的特征向量表示成v¹，v²，…，v^K；

(9)将得到的特征向量排成一行，构造出特征向量矩阵U＝[v¹，v²，…，v^K]，并对其归一化，得到规范化特征向量矩阵

(10)将规范化特征向量矩阵Y的每一行看成是空间中的一个点，用K-均值聚类算法将其聚为K类，并使用正交初始化寻找初始聚类中心；

(11)根据得到的聚类中心，使用欧氏距离测度把SAR图像中的所有像素点划分到不同的类别中，每一个像素点获得一个聚类标签，将这些聚类标签重新组合成和原图像尺寸相同的矩阵，将矩阵元素作为分割结果图像的灰度值并显示出来，从而得到最终分割结果。

2.根据权利要求1所述的一种基于并行稀疏谱聚类的SAR图像分割方法，其中步骤(2)所述的配置MATLAB7.8 R2009(a)版本的并行计算环境，按照如下步骤进行：

2a)在高性能集群环境上安装MATLAB7.8 R2009(a)版本的并行计算工具箱和分布式计算服务器；

2b)在各个集群节点上启动分布式计算服务器；

2c)在集群的主节点上启动并行任务调度器；

2d)在各个处理器节点上启动工作进程。

3.根据权利要求1所述的一种基于并行稀疏谱聚类的SAR图像分割方法，其中步骤(3)所述的并行任务划分，按照如下步骤进行：

3a)根据步骤(2)配置好的MATLAB并行计算环境，指定所使用的并行任务调度器；

3b)使用该并行任务调度器创建作业；

3c)在所创建的作业中，将像素点总个数为n的输入数据样本均衡划分到p个处理器上，遇到不能整除的情况则将剩余的数据样本存放到最后一个节点上；

3d)提交作业，并将各个任务安排到相应的任务队列里等待各个处理器节点的执行。

Images