专利名称：基于并行稀疏谱聚类的sar图像分割方法
技术领域：
本发明属于图像处理技术领域，涉及SAR图像分割，可用于雷达目标检测和目标 识别。
背景技术：
合成孔径雷达SAR具有全天时、全天候的探测与侦察能力。它利用脉冲压缩技术 获得高的距离分辨率，利用合成孔径原理提高方位分辨率，从而相比真实孔径雷达在遥感 领域具有独特的优势。对SAR图像的理解和解译属于图像处理范畴，还涉及信号处理，模式 识别及机器学习等众多学科。SAR图像分割作为SAR图像处理的关键环节之一，在国防和民 用领域正受到越来越广泛的关注。现有的SAR图像分割方法大致可以分为基于区域的分割 方法和基于边界的分割方法，如阈值分割、形态学的方法、聚类的方法、及随机场的方法等。其中，基于聚类的SAR图像分割方法，是将SAR图像中相似的区域尽量划分成一 类，将不相似的区域划分在不同的类别中。已经有很多成熟的聚类算法被用到SAR图像分 割中。作为一门新兴学科，结合谱图理论的谱聚类算法充分发掘了数据的成对点相似特性， 通过图的拉普拉斯矩阵特征分解达到降维目的，使得聚类结果更加稳健。但是，当数据规模 n比较大的时候，谱聚类算法需要计算大小为nXn的相似矩阵，并计算相应拉普拉斯矩阵 的特征分解问题，其计算的时间复杂度和空间复杂度分别为0(n3)和0(n2)，不能有效地对 包含大规模数据的SAR图像进行分割。

发明内容
本发明的目的在于克服上述已有技术的不足，提出了一种基于并行稀疏谱聚类的 SAR图像分割方法，以减小谱聚类的时间和空间复杂度，从而能够有效地对包含大规模数据 的SAR图像进行分割。为实现上述目的，本发明包括以下步骤(1)对待分割的SAR图像进行3层平稳小波变换，图像像素点总个数为n，通过下 式对每个像素点提取出10维子带能量特征，构成大小为nX 10的输入数据样本E 其中，MXN为利用滑动窗口确定的子带大小，这里取值为16X16，coef(i，j)为平 稳小波子带中第i行第j列的系数值；(2)配置MATLAB7. 8R2009 (a)版本的并行计算环境；(3)进行并行任务划分，使用并行任务调度器将像素点总个数为n的输入数据样 本均衡划分到P个处理器上，遇到不能整除的情况则将剩余的数据样本存放到最后一个节 占上-(4)在各个处理器节点上分布式地计算n/p个数据样本到其他所有数据样本之间 的欧氏距离，得到距离矩阵；
(5)对距离矩阵的每一列按照从小到大顺序进行排序，并根据下面的公式计算分 块相似矩阵S 其中，Xl，...，xn表示n个数据样本点，| |X「Xj| |2表示排序后的距离矩阵，o表示 尺度参数，对相似矩阵只保留每个数据样本与周围离它最相近的t个数据样本之间的值， 其他位置都稀疏化成0值；(6)将各个处理器节点计算得到的相似矩阵通过并行任务调度器收集到主节点 上，并排成一行汇聚成一个大小为nXn的稀疏相似矩阵W ；(7)根据稀疏相似矩阵W计算对角线元素为< 的度矩阵D，进而计算出拉普
拉斯矩阵L = D_1/2WD_"2，利用MATLAB的eigs函数对拉普拉斯矩阵进行特征分解；(8)将特征分解后得到的特征值进行排序，得到前K个最大特征值1
其对应的特征向量表示成vW， ，/;(9)将得到的特征向量排成一行，构造出特征向量矩阵U = [v1，v2,…，/]，并对
其归一化，得到规范化特征向量矩阵 (10)将规范化特征向量矩阵Y的每一行看成是空间□ K中的一个点，用K-均值聚 类算法将其聚为K类，并使用正交初始化寻找初始聚类中心；(11)根据得到的聚类中心，使用欧氏距离测度把SAR图像中的所有像素点划分到 不同的类别中，每一个像素点获得一个聚类标签，将这些聚类标签重新组合成和原图像尺 寸相同的矩阵，将矩阵元素作为分割结果图像的灰度值并显示出来，从而得到最终分割结合.本发明与现有的技术相比具有以下优点1)本发明利用MATLAB并行计算环境，设计用于SAR图像分割的并行稀疏谱聚类方 法，能够有效解决现有谱聚类技术运算量过大的问题；2)本发明利用稀疏矩阵存储谱聚类的大规模相似矩阵，能够有效解决现有谱聚类 技术在存储空间上的局限性；3)本发明相对已有的SAR图像分割技术，设计的并行稀疏谱聚类方法能够充分发 掘SAR图像的内在信息，获得更好的分割效果仿真实验结果表明，本发明提出的并行稀疏谱聚类方法能够有效地应用于SAR图 像分割，并进一步应用于雷达目标检测和目标识别。


图1是本发明的实现流程图；图2是本发明的MATLAB并行计算环境示意图；图3是本发明仿真使用的集群环境示意图；图4是本发明与现有K-均值聚类算法和Nystrtoi近似谱聚类算法应用于 ChinaLake机场SAR图像分割的对比实验结果图；图5是本发明与现有K-均值聚类算法和Nystr6m近似谱聚类算法应用于厦门鼓浪屿SAR图像分割的对比实验结果图； 图6是本发明与现有K-均值聚类算法和Nystrdm近似谱聚类算法应用于瑞士湖区 SAR图像分割的对比实验结果图；图7是本发明对瑞士湖区SAR图像分割的并行加速比和效率示意图。
具体实施例方式参照图1，本发明的具体实施过程如下步骤1.提取待分割SAR图像的特征。SAR图像不仅数据量大，在成像过程中不同的地物有不同的后向发射和散射特性， 从而具有丰富的幅度、相位、极化和纹理信息，并且，图像固有的相干斑噪声对于分割性能 直接产生影响，因此，有必要在图像分割前对SAR图像进行纹理分析，提取有效的纹理特征 进行聚类。在以上分析基础上，先对待分割的SAR图像进行3层平稳小波变换，图像像素点总 个数为n，通过下式对每个像素点提取出10维子带能量特征，构成大小为nX 10的输入数据 样本E: 其中，MXN为利用滑动窗口确定的子带大小，这里取值为16X16，coef(i, j)为 平稳小波子带中第i行第j列的系数值，这样，对于大小为256X256的SAR图像，就构成 65536X10的矩阵作为并行稀疏谱聚类方法的输入数据。步骤2.配置MATLAB并行计算环境。本步骤的具体实施参照图2如下2a)在如图3所示的高性能64刀片集群环境中安装MATLAB 7. 8R2009 (a)版本的 并行计算工具箱和分布式计算服务器；2b)在各个集群节点上用mdce start命令启动分布式计算服务器；2c)在集群的主节点上用startjobmanager命令启动并行任务调度器；2d)在各个处理器节点上用startworker命令启动工作进程。步骤3.并行任务划分。3a)根据配置好的MATLAB并行计算环境，使用f indResource命令指定配置好的任 务调度器；3b)使用该任务调度器，用createjob命令创建作业；3c)在所创建的作业中，将规模为n的输入数据样本用createTask命令均衡划分 到P个处理器上，遇到不能整除的情况则将剩余的数据样本存放到最后一个节点上；3d)使用submit命令提交作业，然后用waitForState命令将各个任务安排到相应 的任务队列里等待各个处理器节点的执行；3e)各个处理器节点分布式地计算n/p个数据样本到其他所有数据样本之间的欧 氏距离，获得距离矩阵，对得到的距离矩阵的每一列按照从小到大顺序进行排序，并根据下 面的公式将其转换成分块相似矩阵S
6
其中，Xl，...，xn表示n个数据样本点，| |Xi_Xj| |2表示排序后的距离矩阵，o表示 尺度参数，对相似矩阵只保留每个数据样本与周围离它最相近的t个数据样本之间的值， 其他位置都稀疏化成0值。步骤4.主节点任务收集。用getAllOutputArguments命令将p个处理器节点的计算结果收集到主节点上， 并汇聚成一个大小为nXn的完整的稀疏相似矩阵W = [Sp Sp]，其中Sji = 1，...， P)表示P个处理器节点计算得到的分块相似矩阵。步骤5.计算拉普拉斯矩阵特征分解问题。根据这个稀疏相似矩阵W计算对角线元素为
度矩阵D，进而计算出拉普拉 斯矩阵L = D_1/2WD_"2，利用MATLAB的eigs函数对这个大规模稀疏拉普拉斯矩阵进行特征 分解，得到前K个最大特征值1 =入≥入2≥…≥，以及这些特征值对应的特征向量步骤6.对归一化的特征向量矩阵进行K-均值聚类。将得到的特征向量排成一行，构造出特征向量矩阵U= [vSv2,…，/]，对其归一
化得到规范化特征向量矩阵
，将规范化特征向量矩阵Y的每一行看成是空
间口 K中的一个点，用K-均值聚类算法将其聚为K类，并使用正交初始化寻找初始聚类中 心。步骤7.输出分割结果。根据得到的聚类中心，使用欧氏距离测度把SAR图像中的所有像素点划分到不同 的类别中，每一个像素点获得一个聚类标签，将这些聚类标签重新组合成和原图像尺寸相 同的矩阵，将矩阵元素作为分割结果图像的灰度值并显示出来，从而得到最终分割结果。
本发明的效果可以通过如下仿真实验进行验证。对比实验中给出本发明的并行稀疏谱聚类方法，以及现有的K-均值聚类算法和 NyStr6m近似谱聚类算法的分割结果。其中，现有的K-均值聚类和NystrSm近似谱聚类选取 的是10次实验视觉效果最好的分割结果，本发明的并行稀疏谱聚类方法由于结果比较稳 定选取的是1次实验的分割结果。NystrSm近似谱聚类选取的样本点个数为100，尺度参数 s = 0. 2。并行稀疏谱聚类选取的最近邻数为100，尺度参数s = 0. 05。使用以上三种分割方法对3幅SAR图像进行分割的对比实验结果如下1)对Ku波段SAR图像的分割结果如图4所示，其中图4 (a)为美国加利福尼亚州的China Lake机场3米分辨率的Ku波段SAR图像， 图像尺寸为400 X 400，包括跑道、空地和机场建筑物三类地物。图4(b)为现有的K-均值聚类方法的分割结果，可见跑道上有很多细小的斑点被 误分成其他两类。图4(c)为现有的Nystriim近似方法的分割结果，尽管机场跑道的区域一致性得到 了改善，但是右上方的机场建筑物则被分割得支离破碎，在实际应用中不利于准确进行目 标识别。
图4(d)为本发明的并行稀疏谱聚类方法的分割结果，它准确给出了完整的跑道 和建筑物的轮廓，并且，对于三岔形跑道中央密布的路网，相比其他两种算法获得了更加完 整的区域。2)对C-波段SAR图像的分割结果如图所示，其中图5 (a)为一幅RadarSAT-2卫星拍摄的中国厦门鼓浪屿地区C_波段SAR图像，极 化方式为HV极化，截取图像块尺寸为400X400。图5(b)为现有的K-均值聚类方法的分割结果，它受岛屿上地形分布干扰太大，造 成陆地部分很多误分，延伸进水域的一些陆地与岛屿之间被零散地分隔开来。图5 (c)为现有的Nystr6m近似方法的分割结果，尽管比K-均值的结果有些改善， 但还是未能获得完整的一块区域。图5(d)是本发明的并行稀疏谱聚类方法的分割结果，可以清晰地看出鼓浪屿和 厦门岛一样被完整划分成陆地部分，并且对于一些与水域对比度差异较小的海岸线也能和 岛屿部分保持较好的连接度。3)对X-波段SAR图像的分割结果如图6所示，其中图 6 (a)为 Space Radar Laboratory Missions 的一幅 1994 年瑞士一湖区的 3 视 X-波段SAR图像，图像尺寸为512X512，包含三类地物湖泊、城区和山地。图6 (b)为现有的K-均值聚类方法的分割结果，它很好地将湖泊划分出来，但是， 在左右两边的山地地区同时也存在了严重的错分成水域的情况。图6 (c)为现有的Nystr6m近似方法的分割结果，它对山地区域的分割情况有所改 善，但是城区和山地有很多的混淆。图6(d)为本发明的并行稀疏谱聚类方法的分割结果，它在三种算法中最佳，它不 仅将湖泊沿岸准确提取出来，在山地区域没有丝毫误分成湖泊的情况，并且，城区的分割结 果保持了较好的区域完整性。因此，本发明的并行稀疏谱聚类方法被用来分割这幅SAR图 像是十分合适的。在本发明的并行计算中，对于p个处理器，Ts表示在串行计算机上算法从开始到结 束所需要的时间，TP表示在并行计算机上算法从开始到最后一个处理器完成任务所需要的 时间，加速比SP和效率EP分别定义为SP = TS/TP, EP = Ss/P。加速比衡量的是并行算法的 执行速度相对于串行算法的执行速度加快了多少倍。效率衡量的是单个处理器的计算能力 被有效利用的比率。用本发明的并行稀疏谱聚类方法对图6(a)所示的瑞士湖区SAR图像进行分割，其 并行加速比曲线和并行效率曲线分别如图7(a)和图7(b)所示。其中，所用到的处理器个 数最多达到了 42个。从图7可以看出，由于本发明的并行稀疏谱聚类方法采用MATLAB分布式计算，除 了任务的调度和数据共享的传输之外，各个处理器之间并不需要进行通信，因此在处理器 个数小于24时可以得到近乎线性的加速比。而当处理器个数上升到42个时，由于大规模 稀疏矩阵传输以及多处理器任务调度耗时增加，使得加速比有所下降，但是还是保持着上 升的趋势。在效率方面，小于32个处理器的并行情况下，效率都可以达到90%以上，表明各 个处理器节点得到了相对充分的利用。处理器较多时，消息阻塞和消息延时降低了单个处 理器的利用率。
权利要求
一种基于并行稀疏谱聚类的SAR图像分割方法，包括以下步骤(1)对待分割的SAR图像进行3层平稳小波变换，图像像素点总个数为n，通过下式对每个像素点提取出10维子带能量特征，构成大小为n×10的输入数据样本E <mrow><mi>E</mi><mo>=</mo><mfrac> <mn>1</mn> <mrow><mi>M</mi><mo>&times;</mo><mi>N</mi> </mrow></mfrac><munderover> <mi>&Sigma;</mi> <mrow><mi>i</mi><mo>=</mo><mn>1</mn> </mrow> <mi>M</mi></munderover><munderover> <mi>&Sigma;</mi> <mrow><mi>j</mi><mo>=</mo><mn>1</mn> </mrow> <mi>N</mi></munderover><mo>|</mo><mi>coef</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo></mrow><mo>|</mo> </mrow>其中，M×N为利用滑动窗口确定的子带大小，这里取值为16×16，coef(i，j)为平稳小波子带中第i行第j列的系数值；(2)配置MATLAB7.8R2009(a)版本的并行计算环境；(3)进行并行任务划分，使用并行任务调度器将像素点总个数为n的输入数据样本均衡划分到p个处理器上，遇到不能整除的情况则将剩余的数据样本存放到最后一个节点上；(4)在各个处理器节点上分布式地计算n/p个数据样本到其他所有数据样本之间的欧氏距离，得到距离矩阵；(5)对距离矩阵的每一列按照从小到大顺序进行排序，并根据下面的公式计算分块相似矩阵S <mrow><mi>S</mi><mo>=</mo><mi>exp</mi><mrow> <mo>(</mo> <mo>-</mo> <mfrac><msup> <mrow><mo>|</mo><mo>|</mo><msub> <mi>x</mi> <mi>i</mi></msub><mo>-</mo><msub> <mi>x</mi> <mi>j</mi></msub><mo>|</mo><mo>|</mo> </mrow> <mn>2</mn></msup><msup> <mrow><mn>2</mn><mi>&sigma;</mi> </mrow> <mn>2</mn></msup> </mfrac> <mo>)</mo></mrow><mo>,</mo> </mrow>i，j＝1，...，n其中，x1，...，xn表示n个数据样本点，||xi-xj||2表示排序后的距离矩阵，σ表示尺度参数，对相似矩阵只保留每个数据样本与周围离它最相近的t个数据样本之间的值，其他位置都稀疏化成0值；(6)将各个处理器节点计算得到的相似矩阵通过并行任务调度器收集到主节点上，并排成一行汇聚成一个大小为n×n的稀疏相似矩阵W；(7)根据稀疏相似矩阵W计算对角线元素为的度矩阵D，进而计算出拉普拉斯矩阵L＝D-1/2WD-1/2，利用MATLAB的eigs函数对该拉普拉斯矩阵进行特征分解；(8)对特征分解后得到的特征值进行排序，得到前K个最大特征值1＝λ1≥λ2≥…≥λK，其对应的特征向量表示成v1，v2…，vK；(9)将得到的特征向量排成一行，构造出特征向量矩阵U＝[v1，v2，…，vK]，并对其归一化，得到规范化特征向量矩阵(10)将规范化特征向量矩阵Y的每一行看成是空间口K中的一个点，用K-均值聚类算法将其聚为K类，并使用正交初始化寻找初始聚类中心；(11)根据得到的聚类中心，使用欧氏距离测度把SAR图像中的所有像素点划分到不同的类别中，每一个像素点获得一个聚类标签，将这些聚类标签重新组合成和原图像尺寸相同的矩阵，将矩阵元素作为分割结果图像的灰度值并显示出来，从而得到最终分割结果。FSA00000086656300013.tif,FSA00000086656300021.tif
2.根据权利要求1所述的一种基于并行稀疏谱聚类的SAR图像分割方法，其中步骤 (2)所述的配置MATLAB7. 8R2009 (a)版本的并行计算环境，按照如下步骤进行2a)在高性能集群环境上安装MATLAB7. 8R2009 (a)版本的并行计算工具箱和分布式计 算服务器；2b)在各个集群节点上启动分布式计算服务器；2c)在集群的主节点上启动并行任务调度器； 2d)在各个处理器节点上启动工作进程。
3.根据权利要求1所述的一种基于并行稀疏谱聚类的SAR图像分割方法，其中步骤 (3)所述的并行任务划分，按照如下步骤进行3a)根据步骤(2)配置好的MATLAB并行计算环境，指定所使用的并行任务调度器； 3b)使用该并行任务调度器创建作业；3c)在所创建的作业中，将像素点总个数为η的输入数据样本均衡划分到P个处理器 上，遇到不能整除的情况则将剩余的数据样本存放到最后一个节点上；3d)提交作业，并将各个任务安排到相应的任务队列里等待各个处理器节点的执行。
全文摘要
本发明公开了一种基于并行稀疏谱聚类的SAR图像分割方法，它涉及图像处理技术领域，主要解决已有谱聚类技术在解决大规模SAR图像分割应用上的局限性问题。其实现过程是1)对待分割的SAR图像提取特征；2)配置MATLAB并行计算环境；3)将任务划分到各个处理器节点，计算分块稀疏相似矩阵；4)并行任务调度器收集计算结果并汇聚成完整的稀疏相似矩阵；5)求解拉普拉斯矩阵并进行特征分解；6)对经过归一化后的特征向量矩阵进行K-均值聚类；7)输出SAR图像分割结果。本发明能够有效克服已有谱聚类技术在运算量和存储空间上的瓶颈问题，对大规模SAR图像进行分割效果显著，适用于SAR图像目标检测和目标识别。
文档编号G01S13/90GK101853491SQ20101016149
公开日2010年10月6日 申请日期2010年4月30日 优先权日2010年4月30日
发明者庄雄, 张佳, 朱虎明, 李阳阳, 焦李成, 王爽, 缑水平, 钟桦 申请人:西安电子科技大学