课程设计报告-5818

班 级 学 号 本科课程设计报告 题 目 光学遥感图像地物分类方法与数据 处理实验 学 院 电子工程学院 专 业 遥感科学与技术 学生姓名 导师姓名 西安电子科技大学电 子 工 程 学 院课程设计（报告）任务书学生姓名 指导教师 职称副教授学生学号专业 遥感科学与技术题目 光学遥感图像地物分类方法与数据处理实验相关专业课程 MTLAB；图像处理与识别；模式识别 任务与要求1、 熟练掌握光学遥感图像中地物可分性的原理。2、 深入学习KNN与K-means光学遥感图像数据分类方法。3、 利用MATLAB软件编写KNN与K-means识别算法，并对微波遥感实现典型地物目标的分类。开始日期2016 年 12 月 26 日完成日期2017年 1 月 9 日课程设计所在单位电子工程学院 年 月 日本表格由电子工程学院网络信息中心 编辑录入 .光学遥感图像地物分类方法与数据处理实验（电子工程学院遥感科学与技术）摘要：遥感具有覆盖范围广、包含信息量大、获取信息快等优点,其所获得的图像已广泛应用于国防安全与国民经济的众多领域。遥感图像分类是遥感图像解译的关键技术,受到各国研究者的广泛关注,成为近年来的研究热点。而对遥感图像地物分类算法的进一步研究可以有效改善这些问题，所以我们进行光学遥感图像分类方法与数据处理实验是很有必要的。1本论文是通过Matlab，具体以KNN以及K-means算法对光学遥感图像进行地物分类仿真，并对仿真结果进行了数据分析。在最后，我们还总结了全篇所做工作并对未来光学遥感图像分类的发展趋势做了展望。关键词：遥感图像，图像分类，KNN， K-meansAbstract: Remote sensing has many advantages, such as wide coverage, large amount of information and fast access to information. The obtained image has been widely used in many fields of national defense security and national economy. Remote sensing image classification is the key technology of remote sensing image interpretation, which has been the focus of attention of researchers all over the world and has become a hotspot in recent years. Further research on remote sensing image feature classification algorithm can effectively improve these problems, so we carry out optical remote sensing image classification method and data processing experiment is necessary. In this paper, the classification of optical remote sensing image is simulated by using Matlab and KNN and K-means algorithm, and the simulation results are analyzed. In the end, we also summarize the work done in the whole article and prospect of the development trend of the future classification of optical remote sensing images.Key words: optical remote sensing image, feature classification, KNN, K-means1 引言1.1 目前遥感图像地物分类现状 先从传统的遥感分类方面说起，该分类方法是目前运用较多，算法比较成熟的方法。分为监督分类和非监督分类，他们的原理都是根据图像象元的光谱特征的相似度来进行的分类。监督分类用于用户对分类区比较熟悉，由用户自己控制，非监督分类则是将象元相似度大小进行归类合并。但是未充分利用遥感图像提供的多种信息，只考虑多光谱特征，没有利用到地物空间关系、空间位置形状、纹理等方面的信息。1、监督分类 监督分类可根据应用目标和区域，有选择地决定分类类别，可控制样本的选择，避免了非监督分类中对光谱集群组的重新归类。但个人认为其人为主观因素较强，操作者所选择的训练样本有可能不是很典型并且有可能不能反映图像的真实情况，所以图像中同一类别的光谱差异和人为因素，有可能造成样本没有代表性，并且训练样本的选取和评估需要花费较多的人力和时间。2、非监督分类非监督分类过程不需要任何的先验知识，仅凭遥感影像地物光谱特征的分布规律，随其自然地进行分类。但是看文献时看到，非监督分类还有一个前提，那就是：假定遥感影像上同类地物在同样条件下具有相同的光谱信息特征。如果产生的光谱万一不一定对应于操作者想要的类别，且操作者较难对产生的类别进行控制，比如图像中各类别的光谱特征会随时间、地形等变化，不同图像以及不同时段的图像之间的光谱无法保持其连续性，从而使不同图像之间的对比变得困难。2无论是监督分类还是非监督分类， 都是依据地物的光谱特性的点独立原则来进行分类的，且都是采用的统计方法。该方法只是根据各波段灰度数据的统计特征进行的，加上卫星遥感数据的分辨率的限制，一般图像的像元很多是混合像元，带有混合光谱信息的特点，致使计算机分类面临着诸多模糊对象，不能确定其究竟属于哪一类地物。而且，同物异谱和异物同谱的现象普遍存在，也会导致误分、漏分情况的出现，因此人们不断尝试新方法来加以改善和提高遥感图像分类的效率和质量。 这些新方法主要有决策树分类法、综合阈值法、专家系统分类法、多特征融合法、神经网络分类法以及基于频谱特征的分类法等。近年来的研究大多将传统方法与新方法加以结合。即在非监督分类和监督分类的基础上， 运用新方法来改进，减少错分和漏分情况，对遥感图像的分类精度有了一定程度的增强3。1.2 遥感发展历程1、无记录地面遥感阶段（1608-1838）1608年汉斯李波尔赛制造了世界第一架望远镜 1609年伽利略制作了放大三倍的科学望远镜并首次观测月球 1794年气球首次升空侦察为观测远距离目标开辟了先河，但望远镜观测不能把观测到的事物用图像的方式记录下来。 2、有记录地面遥感阶段（1839-1857）1839年达盖尔（Daguarre）发表了他和尼普斯（Niepce）拍摄的照片，第一次成功将拍摄事物记录在胶片上 1849年法国人艾米劳塞达特（Aime Laussedat）制定了摄影测量计划，成为有目的有记录的地面遥感发展阶段的标志。 3、 初期发展又可称为空中摄影遥感阶段（1858-1956） 1858年用系留气球拍摄了法国巴黎的鸟瞰像片 ，1903年飞机的发明 ，1909年第一张航空像片。 一战期间（1914-1918）形成独立的航空摄影测量学的学科体系。 二战期间（1931-1945）彩色摄影、红外摄影、雷达技术、多光谱摄影、扫描技术以及运载工具和判读成图设备。4、 现代遥感技术从1957年开始迅速发展1957年：前苏联发射了人类第一颗人造地球卫星。20世纪60年代：美国发射了TIROS、ATS、ESSA等气象卫星和载人宇宙飞船。1972年：发射了地球资源技术卫星ERTS-1（后改名为LandsatLandsat-1），装有MSS感器，分辨率79米。 1982年Landsat-4发射，装有TM传感器，分辨率提高到30米。1986年法国发射SPOT-1，装有PAN和XS遥感器，分辨率提10米。1999年美国发射IKNOS，空间分辨率提高到1米。【4】1.3 遥感图像分类原理概述通常我们所指的遥感图像是指卫星探测到的地物亮度特征, 它们构成了光谱空间。每种地物有其固有的光谱特征, 它们位于光谱空间中的某一点。但由于干扰的存在, 环境条件的不同, 例如: 阴影,地形上的变化, 扫描仪视角, 干湿条件, 不同时间拍摄及测量误差等, 使得测得的每类物质的光谱特征不尽相同, 同一类物质的各个样本在光谱空间是围绕某一点呈概率分布, 而不是集中到一点, 但这仍使我们可以划分边界来区分各类。因此, 我们就要对图像进行分类。图像分类的任务就是通过对各类地物波谱特征的分析选择特征参数, 将特征空间划分为不相重叠的子空间, 进而把影像内诸像元划分到各子间去, 从而实现分类。分类方法可以分为统计决策法( 判别理论识别法) 模式识别和句法模式识别。统计决策法模式识别指的是: 对研究对象进行大量的统计分析, 抽出反映模式的本质特点、特征而进行识别。主要的有监督分类中的最小距离法、逐次参数估计法、梯度法、最小均方误差法、费歇准则法和非监督分类中的按批修改的逐步聚类法、等混合距离法。此外还可以将两者结合起来, 互相补充以获得较好的效果。句法模式识别则需要了解图像结构信息, 从而对其进行分类5。2 基础原理2.1 KNN算法原理邻近算法，或者说K最近邻(kNN，k-NearestNeighbor)分类算法是数据挖掘分类技术中最简单的方法之一。所谓K最近邻，就是k个最近的邻居的意思，说的是每个样本都可以用它最接近的k个邻居来代表6。KNN算法的核心思想是如果一个样本在特征空间中的k个最相邻的样本中的大多数属于某一个类别，则该样本也属于这个类别，并具有这个类别上样本的特性。该方法在确定分类决策上只依据最邻近的一个或者几个样本的类别来决定待分样本所属的类别。 kNN方法在类别决策时，只与极少量的相邻样本有关。由于kNN方法主要靠周围有限的邻近的样本，而不是靠判别类域的方法来确定所属类别的，因此对于类域的交叉或重叠较多的待分样本集来说，kNN方法较其他方法更为适合。图2.1.1 KNN分类图 上图中，绿色圆要被决定赋予哪个类，是红色三角形还是蓝色四方形？如果K=3，由于红色三角形所占比例为2/3，绿色圆将被赋予红色三角形那个类，如果K=5，由于蓝色四方形比例为3/5，因此绿色圆被赋予蓝色四方形类。 K最近邻(k-Nearest Neighbor，KNN)分类算法，是一个理论上比较成熟的方法，也是最简单的机器学习算法之一。该方法的思路是：如果一个样本在特征空间中的k个最相似(即特征空间中最邻近)的样本中的大多数属于某一个类别，则该样本也属于这个类别。KNN算法中，所选择的邻居都是已经正确分类的对象。该方法在定类决策上只依据最邻近的一个或者几个样本的类别来决定待分样本所属的类别。 KNN方法虽然从原理上也依赖于极限定理，但在类别决策时，只与极少量的相邻样本有关。由于KNN方法主要靠周围有限的邻近的样本，而不是靠判别类域的方法来确定所属类别的，因此对于类域的交叉或重叠较多的待分样本集来说，KNN方法较其他方法更为适合。KNN算法不仅可以用于分类，还可以用于回归。通过找出一个样本的k个最近邻居，将这些邻居的属性的平均值赋给该样本，就可以得到该样本的属性。更有用的方法是将不同距离的邻居对该样本产生的影响给予不同的权值(weight)，如权值与距离成反比。算法流程：1. 准备数据，对数据进行预处理2. 选用合适的数据结构存储训练数据和测试元组3. 设定参数，如k4.维护一个大小为k的的按距离由大到小的优先级队列，用于存储最近邻训练元组。随机从训练元组中选取k个元组作为初始的最近邻元组，分别计算测试元组到这k个元组的距离，将训练元组标号和距离存入优先级队列5. 遍历训练元组集，计算当前训练元组与测试元组的距离，将所得距离L 与优先级队列中的最大距离Lmax6. 进行比较。若L=Lmax，则舍弃该元组，遍历下一个元组。若L Lmax，删除优先级队列中最大距离的元组，将当前训练元组存入优先级队列。7. 遍历完毕，计算优先级队列中k 个元组的多数类，并将其作为测试元组的类别。8. 测试元组集测试完毕后计算误差率，继续设定不同的k值重新进行训练，最后取误差率最小的k 值优点：1.简单，易于理解，易于实现，无需估计参数，无需训练；2. 适合对稀有事件进行分类；3.特别适合于多分类问题(multi-modal,对象具有多个类别标签)， kNN比SVM的表现要好。缺点：该算法在分类时有个主要的不足是，当样本不平衡时，如一个类的样本容量很大，而其他类样本容量很小时，有可能导致当输入一个新样本时，该样本的K个邻居中大容量类的样本占多数。 该算法只计算“最近的”邻居样本，某一类的样本数量很大，那么或者这类样本并不接近目标样本，或者这类样本很靠近目标样本。无论怎样，数量并不能影响运行结果。该方法的另一个不足之处是计算量较大，因为对每一个待分类的文本都要计算它到全体已知样本的距离，才能求得它的K个最近邻点。可理解性差，无法给出像决策树那样的规则。改进策略： kNN算法因其提出时间较早，随着其他技术的不断更新和完善，kNN算法的诸多不足之处也逐渐显露，因此许多kNN算法的改进算法也应运而生。针对以上算法的不足，算法的改进方向主要分成了分类效率和分类效果两方面。分类效率：事先对样本属性进行约简，删除对分类结果影响较小的属性，快速的得出待分类样本的类别。该算法比较适用于样本容量比较大的类域的自动分类，而那些样本容量较小的类域采用这种算法比较容易产生误分。分类效果：采用权值的方法（和该样本距离小的邻居权值大）来改进，Han等人于2002年尝试利用贪心法，针对文件分类实做可调整权重的k最近邻居法WAkNN (weighted adjusted k nearest neighbor)，以促进分类效果；而Li等人于2004年提出由于不同分类的文件本身有数量上有差异，因此也应该依照训练集合中各种分类的文件数量，选取不同数目的最近邻居，来参与分类。2.2 K-means算法原理K-means算法是硬聚类算法，是典型的基于原型的目标函数聚类方法的代表，它是数据点到原型的某种距离作为优化的目标函数，利用函数求极值的方法得到迭代运算的调整规则。K-means算法以欧式距离作为相似度测度，它是求对应某一初始聚类中心向量V最优分类，使得评价指标J最小。算法采用误差平方和准则函数作为聚类准则函数7。 K-means算法是很典型的基于距离的聚类算法，采用距离作为相似性的评价指标，即认为两个对象的距离越近，其相似度就越大。该算法认为簇是由距离靠近的对象组成的，因此把得到紧凑且独立的簇作为最终目标8。 (2.2.1) 个初始类聚类中心点的选取对聚类结果具有较大的公式影响，因为在该算法第一步中是随机的选取任意个对象作为初始聚类的中心，初始地代表一个簇。该算法在每次迭代中对数据集中剩余的每个对象，根据其与各个簇中心的距离将每个对象重新赋给最近的簇。当考察完所有数据对象后，一次迭代运算完成，新的聚类中心被计算出来。如果在一次迭代前后，J的值没有发生变化，说明算法已经收敛。 算法过程如下: 1)从N个文档随机选取个文档作为质心 2)对剩余的每个文档测量其到每个质心的距离，并把它归到最近的质心的类 3)重新计算已经得到的各个类的质心 4)迭代23步直至新的质心与原质心相等或小于指定阈值，算法结束 具体如下: 输入:k, datan; (1) 选择k个初始中心点，例如c0=data0,ck-1=datak-1; (2) 对于data0.datan，分别与c0ck-1比较，假定与ci差值最少，就标记为i; (3) 对于所有标记为i点，重新计算ci= 所有标记为i的dataj之和/标记为i的个数;(4) 重复(2)(3)，直到所有ci值的变化小于给定阈值。算法优点：1.算法快速、简单; 2.对大数据集有较高的效率并且是可伸缩性的; 3.时间复杂度近于线性，而且适合挖掘大规模数据集。K-Means聚类算法的时间复杂度是O(nkt) ,其中n代表数据集中对象的数量，t代表着算法迭代的次数，k代表着簇的数目。算法缺点： 在 K-means 算法中 K 是事先给定的，这个 K 值的选定是非常难以估计的。很多时候，事先并不知道给定的数据集应该分成多少个类别才最合适。这也是 K-means 算法的一个不足。有的算法是通过类的自动合并和分裂，得到较为合理的类型数目 K，例如 ISODATA 算法。关于 K-means 算法中聚类数目K 值的确定在文献中，是根据方差分析理论，应用混合 F统计量来确定最佳分类数，并应用了模糊划分熵来验证最佳分类数的正确性。在文献中，使用了一种结合全协方差矩阵的 RPCL 算法，并逐步删除那些只包含少量训练数据的类。而文献中使用的是一种称为次胜者受罚的竞争学习规则，来自动决定类的适当数目。它的思想是:对每个输入而言，不仅竞争获胜单元的权值被修正以适应输入值，而且对次胜单元采用惩罚的方法使之远离输入值。 在 K-means 算法中，首先需要根据初始聚类中心来确定一个初始划分，然后对初始划分进行优化。这个初始聚类中心的选择对聚类结果有较大的影响，一旦初始值选择的不好，可能无法得到有效的聚类结果，这也成为 K-means算法的一个主要问题。对于该问题的解决，许多算法采用遗传算法(GA)，例如文献 中采用遗传算法(GA)进行初始化，以内部聚类准则作为评价指标。 从 K-means 算法框架可以看出，该算法需要不断地进行样本分类调整，不断地计算调整后的新的聚类中心，因此当数据量非常大时，算法的时间开销是非常大的。所以需要对算法的时间复杂度进行分析、改进，提高算法应用范围。在文献中从该算法的时间复杂度进行分析考虑，通过一定的相似性准则来去掉聚类中心的侯选集。而在文献中，使用的 K-means 算法是对样本数据进行聚类，无论是初始点的选择还是一次迭代完成时对数据的调整，都是建立在随机选取的样本数据的基础之上，这样可以提高算法的收敛速度9。3 方案设计及仿真实现3.1 KNN部分主函数中包含两个子函数Data_Gen用于产生模拟数据，这里产生的数据一个1000*3的矩阵，其中前两列代表数据，即数据是一个包含两个元素的向量，有效数据为0-10之间的随机数，最后一列为分类号，这里假设一共分为4类，分类号为1,2,3,4 图3.1.1 数据分类图KNNdatgingTest_2Wei子函数用于实现KNN算法1参数设置，这里选取的K值为4，测试数据比例为0.1，也就是说1000组数据中有900组为样本，剩余的为测试数据，因为数据都是随机产生的，所以选定前100组数据位测试数据，其余的为样本数据；2导入数据后对数据做归一化处理；3按序遍历所有测试数据，每一个测试数据分别与所有测试样本求欧式距离；4对求得欧式距离按从小到大排序；5选取最小的K个，并选取分类号最多的结果作为测试数据的分类号。%该程序用于产生测试用的数据%function Data_Gen%num 代表数据个数clear;clc;row=3;%数据维数+1column=1000;%数据个数format=1;%数据格式选择，为0的时候产生类别明确的数，当为1时产生模糊的数Data=zeros(column,row);%产生一个1000*3的矩阵,前两列为数据，最后一列为类别if (format=0)for n= 1:1:column a=round(3*rand(1);%产生0-3的证书作为类别 if a=0 Data(n,:)=(rand(1)*3+1),(rand(1)*3+6),a+1;%当a=0时，在左上角产生数，并归为1类 elseif a=1 Data(n,:)=(rand(1)*3+6),(rand(1)*3+6),a+1;%当a=1时，在左上角产生数，并归为2类 elseif a=2 Data(n,:)=(rand(1)*3+1),(rand(1)*3+1),a+1;%当a=2时，在左上角产生数，并归为3类 else Data(n,:)=(rand(1)*3+6),(rand(1)*3+1),a+1;%当a=3时，在左上角产生数，并归为4类 endendelsefor n= 1:1:column a=round(3*rand(1);%产生0-3的证书作为类别 if a=0 Data(n,:)=(rand(1)*5),(rand(1)*5+5),a+1;%当a=0时，在左上角产生数，并归为1类 elseif a=1 Data(n,:)=(rand(1)*5+5),(rand(1)*5+5),a+1;%当a=1时，在左上角产生数，并归为2类 elseif a=2 Data(n,:)=(rand(1)*5),(rand(1)*5),a+1;%当a=2时，在左上角产生数，并归为3类 else Data(n,:)=(rand(1)*5+5),(rand(1)*5),a+1;%当a=3时，在左上角产生数，并归为4类 endendendfid=fopen(Data_Test.txt,wt);for n=1:1:column for k=1:1:(row-1) fprintf(fid,%d,Data(n,k); fprintf(fid,%s,); end fprintf(fid,%d,Data(n,row); fprintf(fid,n);endfclose(fid);function relustLabel = KNN(inx,data,labels,k)% % inx 为 输入测试数据，data为样本数据，labels为样本标签 % datarow , datacol = size(data);%求两点间距离，这里使用欧式距离diffMat = repmat(inx,datarow,1) - data ;distanceMat = sqrt(sum(diffMat.2,2);%按升序排列，IX为检索号B , IX = sort(distanceMat,ascend);len = min(k,length(B);%找出出现频率最高的分类号作为最后的分类结果relustLabel = mode(labels(IX(1:len);end% function KNNdatgingTest_2Wei% clcclearclose all% 导入数据 data = load(Data_Test.txt);dataMat = data(:,1:2);labels = data(:,3);len = size(dataMat,1);k = 4;%选取的点数error = 0; % 测试数据比例 Ratio = 0.1;numTest = Ratio * len; % 归一化处理maxV = max(dataMat);minV = min(dataMat);range = maxV-minV;newdataMat = (dataMat-repmat(minV,len,1)./(repmat(range,len,1); % 测试err_num=;%用于记录错误的数据检索号for i = 1:numTest % 前100个为测试数据，后900个为样本数据，并调用KNN分类器 classifyresult = KNN(newdataMat(i,:),newdataMat(numTest+1:len,:),labels(numTest+1:len,:),k); %如果发生分类错误，则记录错误数据索引号 if(classifyresult=labels(i) error = error+1; err_num=err_num,i; %更新分类号为判断的结果,并打印错误号% data(i,3)=classifyresult; fprintf(错误数据%d 测试结果为：%d 真实结果为：%dn,i classifyresult labels(i) endend % 计算准确率 fprintf(准确率为：%fn,1-error/(numTest); %画出分类后的结果图% figure(); %画出样本数据以.表示 for i=numTest+1:1:len if(data(i,3)=1) plot(data(i,1),data(i,2),-r);hold on%红色 elseif(data(i,3)=2) plot(data(i,1),data(i,2),-k);hold on%黑色 elseif(data(i,3)=3) plot(data(i,1),data(i,2),-b);hold on%蓝色 else plot(data(i,1),data(i,2),-g);hold on%绿色 end end%画出测试数据以。表示 for i=1:1:numTest if(data(i,3)=1) plot(data(i,1),data(i,2),-rs);hold on%红色 elseif(data(i,3)=2) plot(data(i,1),data(i,2),-ks);hold on%黑色 elseif(data(i,3)=3) plot(data(i,1),data(i,2),-bs);hold on%蓝色 else plot(data(i,1),data(i,2),-gs);hold on%绿色 end end %画出错误数据以x表示 for i=1:1:length(err_num) if(data(err_num(i),3)=1) plot(data(err_num(i),1),data(err_num(i),2),-rx);hold on%红色 elseif(data(err_num(i),3)=2) plot(data(err_num(i),1),data(err_num(i),2),-kx);hold on%黑色 elseif(data(err_num(i),3)=3) plot(data(err_num(i),1),data(err_num(i),2),-bx);hold on%蓝色 else plot(data(err_num(i),1),data(err_num(i),2),-gx);hold on%绿色 end endgrid on;运行结果：错误数据63 测试结果为：1 真实结果为：3准确率为：0. 图3.1.2 KNN仿真结果图 图中”.”代表样本数据，红色、黑色、蓝色、绿色分别代表第一、二、三、四类。图中方框代表测试数据，红色、黑色、蓝色、绿色分别代表判定为第一、二、三、四类。错号代表判别类别错误。此次判定只有一个错误，准确率达到99%。3.2 K-means部分10 图3.2.1 测试图像此测试图像是由小组成员拍摄的日常教学楼的光学图片。我们所做的实验是将这张照片中的事物分为5类，以求得最真实的分类效果。function kmeans_demo1()clear;close all;clc;% 读取测试图像im = imread(IMG.jpg);imshow(im), title(Imput image);% 转换图像的颜色空间得到样本cform = makecform(srgb2lab);lab = applycform(im,cform);ab = double(lab(:,:,2:3);nrows = size(lab,1); ncols = size(lab,2);X = reshape(ab,nrows*ncols,2);figure, scatter(X(1,:),X(2,:),3,filled); box on; %显示颜色空间转换后的二维样本空间分布图3.2.2 颜色空间转换后的二维样本空间分布图 在图像由RGB转换到LAB色彩模式，LAB的空间相互分量相关性比较小，图像分割效果会更好。% 对样本空间进行Kmeans聚类k = 5; % 聚类个数max_iter = 100; %最大迭代次数centroids, labels = run_kmeans(X, k, max_iter); % 显示聚类分割结果figure, scatter(X(1,:),X(2,:),3,labels,filled); %显示二维样本空间聚类效果hold on; scatter(centroids(1,:),centroids(2,:),60,r,filled)hold on; scatter(centroids(1,:),centroids(2,:),30,g,filled)box on; hold off;%print -dpdf 2D2.pdfpixel_labels = reshape(labels,nrows,ncols);rgb_labels = label2rgb(pixel_labels);figure, imshow(rgb_labels), title(Segmented Image);%print -dpdf Seg.pdfend图3.2.3二维样本空间聚类效果图function centroids, labels = run_kmeans(X, k, max_iter)% 该函数实现Kmeans聚类% 输入参数：% X为输入样本集，dxN% k为聚类中心个数% max_iter为kemans聚类的最大迭代的次数% 输出参数：% centroids为最终聚类中心 % labels为样本的类别标记% 采用K-means+算法初始化聚类中心 centroids = X(:,1+round(rand*(size(X,2)-1); labels = ones(1,size(X,2); for i = 2:k D = X-centroids(:,labels); D = cumsum(sqrt(dot(D,D,1); if D(end) = 0, centroids(:,i:k) = X(:,ones(1,k-i+1); return; end centroids(:,i) = X(:,find(rand D/D(end),1); ,labels = max(bsxfun(minus,2*real(centroids*X),dot(centroids,centroid