本发明属于遥感数据分析领域，具体涉及一种高分辨率遥感影像聚类方法。
背景技术：
：高分辨率遥感影像是由高分辨率遥感卫星、高分辨率航拍设备拍摄的，地面分辨率在2米以上的地面影像数据。目前某些设备可以提供地面分辨率达到5厘米的遥感影像数据，可以十分精细的反应地物的特征信息。随着的进步，我们可以获得越来越多的高分辨率遥感影像。遥感影像数据只有经过人工解译、分类、聚类才能获得照片中地物的类型，才能真正有效利用遥感这些数据。影像的聚类可以在不输入样本的情况下对影像中的地物进行分类，可以较低的成本获得地表的地物区域、范围、边界情况，因此遥感影像聚类算法对于遥感数据的处理十分重要。对于遥感影像的聚类，目前业内广泛采用的算法是直接利用k-means(k均值算法)、em、分层次聚类和群智能算法对影像中的每个像元波段值进行分析，进而获得聚类结果。当这种方法面对高分辨率遥感影像时会遇到一个问题：由于高分辨率遥感影像展示了更加精细的地物特征，这种特征导致每一种地物的像元存在高度的类内异质性，进而导致同一地物内部像元有很大差异，此时利用传统的算法进行聚类，一方面获得的聚类结果会存在大量的杂点，导致结果影像中出现“胡椒盐”现象；另一方面，由于异质性的存在，很多明显隶属于同一类目的地物被划分到不同的类目中，导致无法获得期望的聚类结果，聚类准确率不高。所以需要针对高分辨率遥感影像的特点设计一种方法，可以应对单个像元与其周边邻域像元的异质特性，进而获得较好的高分辨率遥感影像聚类结果。技术实现要素：本发明为解决现有技术获得的聚类结果会存在大量的杂点、聚类准确率不高的问题，提供了一种基于范围能量差的高分辨率遥感影像聚类方法。本发明所述一种基于范围能量差的高分辨率遥感影像聚类方法，具体通过以下技术方案实现：步骤一、根据高分辨率遥感影像image的地面分辨率resolution建立第一范围变量scale1和第二范围变量scale2；对于高分辨率遥感影像中的每一个像元构造对应的能量差结构体enstruct；设置迭代次数计数器iteratorcounter＝0；高分辨率遥感影像image的宽度为width，高度为height，波段个数为bandnum，期望的聚类之后的类目数为categorynum；步骤二、对于image中的所有像元，利用k-means算法根据像元的波段值将所有像元聚类为categorynum个类目，获得categorynum个聚类中心；所有聚类中心通过聚类中心列表center来表示，center包含categorynum个元素，每个元素为一个bandnum维变量；利用center初始化image中每个像元对应enstruct的所属类目；步骤三、计算image中所有像元的enstruct的第一范围能量列表、第二范围能量列表、第一范围能量列表的能量差和第二范围能量列表的能量差，根据能量差设置像元的enstruct的最大能量差所在位置；步骤四、设置像元类目转换计数器transformcounter＝0；iteratorcounter＝iteratorcounter+1；步骤五、对image中的所有像元，根据像元的enstruct的最大能量差所在位置energypos的值设置所属类目category，统计category的值发生变化的像元的个数并存储到transformcounter中；步骤六、如果transformcounter<width×height×threshold或者iteratorcounter>maxiteartor则转到步骤七，否则转到步骤三；其中threshold为最小变化百分比，maxiteartor为最大迭代次数；步骤七、对image中的所有像元，将每个像元对应enstruct的category值作为像元最终的类目值，输出聚类结果。本发明与现有技术相比较，最为突出的特点和显著的有益效果是：本发明针对高分辨率遥感影像，通过本方法可以在聚类的过程中综合考虑一个像元本身的波段值和类目值、以及它周围邻域两个区域范围所有像元波段值、类目值，有效的应对高分辨率遥感影像过多地物细节信息带来的高度的类内异质性问题，通过本方法可以对高分辨率遥感影像进行聚类并获得较好的分类结果，聚类结果杂点少、聚类准确率高达94％，远超传统方法。附图说明图1是本发明方法的总流程图；图2是本发明步骤一的流程图；图3是本发明步骤二的流程图；图4是本发明步骤三的流程图；图5是本发明步骤五的流程图；图6是本发明实施例进行聚类的高分辨率遥感影像；图7是本发明实施例进行聚类处理后的结果。具体实施方式具体实施方式一：本实施方式给出的一种基于范围能量差的高分辨率遥感影像聚类方法，如图1所示，具体是按照以下步骤进行的：步骤一、根据高分辨率遥感影像image的地面分辨率resolution建立第一范围变量scale1和第二范围变量scale2；对于高分辨率遥感影像中的每一个像元构造对应的能量差结构体enstruct；设置迭代次数计数器iteratorcounter＝0；其中，round(·)是将数值的小数部分进行四舍五入的函数；resolution表示高分辨率遥感影像的地面分辨率；从上式可以看出，这里取了一大一小两个范围，scale1为小范围变量，scale2为大范围变量；高分辨率遥感影像image的宽度为width，高度为height，波段个数为bandnum，期望的聚类之后的类目数为categorynum；步骤二、步骤二如图3所示，对于image中的所有像元，利用k-means算法(k均值算法)根据像元的波段值将所有像元聚类为categorynum个类目，获得categorynum个聚类中心；所有聚类中心通过聚类中心列表center来表示，center包含categorynum个元素，每个元素为一个bandnum维变量；利用center初始化image中每个像元对应enstruct的所属类目；步骤三、计算image中所有像元的enstruct的第一范围能量列表、第二范围能量列表、第一范围能量列表的能量差和第二范围能量列表的能量差，根据能量差设置像元的enstruct的最大能量差所在位置；第一范围能量列表是指像元的第一范围内的像元的能量值列表，第二范围能量列表是指像元的第二范围内的像元的能量值列表。步骤四、设置像元类目转换计数器transformcounter＝0；iteratorcounter＝iteratorcounter+1；步骤五、对image中的所有像元，根据像元的enstruct的最大能量差所在位置energypos的值设置所属类目category，统计category的值发生变化的像元的个数并存储到transformcounter中；步骤六、如果transformcounter<width×height×threshold或者iteratorcounter>maxiteartor则转到步骤七，否则转到步骤三；其中threshold为最小变化百分比，默认值为0.05；maxiteartor为最大迭代次数，默认值为20；步骤七、对image中的所有像元，将每个像元对应enstruct的category值作为像元最终的类目值，输出聚类结果。具体实施方式二：本实施方式与具体实施方式一不同的是：如图2所示，步骤一具体包括以下步骤：步骤一一、地面的分辨率resolution对应的单位为米，其数值为拍摄高分辨率遥感影像的卫星的固定参数。按如下公式建立第一范围变量scale1和第二范围变量scale2：其中，round(·)是将数值的小数部分进行四舍五入的函数(即，将数值四舍五入进行取整)，resolution表示高分辨率遥感影像的地面分辨率；步骤一二、建立高分辨率遥感影像中的每一个像元p对应的能量差结构体p.enstruct；p.enstruct包含以下属性：(1)像元p所有波段的值构成的矢量pv，pv共有bandnum维，每一个维度对应像元p在对应波段上的值；(2)像元p在高分辨率遥感影像上处在的行数row；(3)像元p在高分辨率遥感影像上处在的列数col；(4)像元p的索引值index，index的计算公式为：index＝row*width+col(5)像元p所属类目category，类目取值范围是[1,categorynum]，初始值为1；(6)像元p的第一范围邻域enstruct列表enlist1，enlist1的初值为空；(7)像元p的第二范围邻域enstruct列表enlist2，enlist2的初值为空；(8)像元p的第一范围能量列表energy1，energy1包含categorynum个元素，每个元素的初值为0；(9)像元p的第二范围能量列表energy2，energy2包含categorynum个元素，每个元素的初值为0；(10)像元p的最大能量差所在位置energypos，初值为1；步骤一三、初始化image中的每一个像元p对应的enlist1和enlist2。其他步骤及参数与具体实施方式一相同。具体实施方式三：本实施方式与具体实施方式二不同的是：步骤一三中初始化image中的每一个像元p对应的enlist1和enlist2的具体过程为：enlist1的初始化方式为：以像元p所在的行数row和列数col为中心在image上取scale1×scale1为窗口范围内的所有像元，并将这些像元的enstruct加入到enlist1列表中；enlist1的元素个数为n1，n1的值为scale1×scale1；enlist2的初始化方式为：以像元p所在的行数row和列数col为中心在image上取scale2×scale2为窗口范围内的所有像元，并将这些像元的enstruct加入到enlist2列表中；enlist2的元素个数为n2，n2的值为scale2×scale2。其他步骤及参数与具体实施方式二相同。具体实施方式四：本实施方式与具体实施方式三不同的是：步骤二中利用center初始化image中每个像元对应enstruct的所属类目的方法为：计算image中的每一个像元p与center列表每一个元素的距离，获得距离最近的元素在center中所处的位置minpos，并将minpos写入到p.enstruct的category属性中；具体步骤包括以下步骤：步骤二一、计算像元p的暂存波段值矢量temppv＝p.enstruct.pv，p.enstruct.pv为像元p对应enstruct的pv属性；步骤二二、设置距离最近的值min＝0，最近的值在center的位置minpos＝0；步骤二三、设定中心点技术器的值centercounter＝1；步骤二四、计算temppv与center[centercounter]的距离dis，center[centercounter]为center的第centercounter个元素；步骤二五、如果dis小于min，或者minpos的值为0，则转到步骤二六，否则转到步骤二七；步骤二六、min＝dis；minpos＝centercounter；步骤二七、centercounter＝centercounter+1；步骤二八、如果centercounter小于等于categorynum则转到步骤二四，否则转到步骤二九；步骤二九、执行p.enstruct.category＝minpos，p.enstruct.category为p对应enstruct的category属性。其他步骤及参数与具体实施方式三相同。具体实施方式五：本实施方式与具体实施方式四不同的是：如图4所示，步骤三中计算像元的enstruct的第一范围能量列表、第二范围能量列表、第一范围能量列表的能量差和第二范围能量列表的能量差，根据能量差设置像元的enstruct的最大能量差所在位置的具体步骤包括：步骤三一、获得image上任意一个像元p的第一范围邻域enstruct列表变量tl1，tl1＝p.enstruct.enlist1；p.enstruct.enlist1表示p.enstruct的enlist1属性；步骤三二、获得像元p的第二范围邻域enstruct列表变量tl2，tl2＝p.enstruct.enlist2；p.enstruct.enlist2表示p.enstruct的enlist2属性；步骤三三、获得像元p在image上的行值变量trow，trow＝p.enstruct.row；获得像元p在image上的列值变量tcol，tcol＝p.enstruct.col；其中p.enstruct.row表示p.enstruct的row属性，p.enstruct.col表示p.enstruct的col属性；步骤三四、设置类目计数器的值lmcounter＝1；步骤三五、计算第一范围能量变量littleenergy的值；步骤三六、计算第二范围能量变量bigenergy的值；步骤三七、设置p.enstruct.energy1[lmcounter]＝littleenergy；p.enstruct.energy1表示p.enstruct的energy1属性，p.enstruct.energy1[lmcounter]指的是p.enstruct.energy1的第lmcounter个元素；步骤三八、设置p.enstruct.energy2[lmcounter]＝bigenergy；p.enstruct.energy2表示p.enstruct的energy2属性，p.enstruct.energy2[lmcounter]指的是p.enstruct.energy2的第lmcounter个元素；步骤三九、lmcounter＝lmcounter+1；步骤三十、如果lmcounter小于等于categorynum则转到步骤三五，否则转到步骤三十一；步骤三十一、计算p.enstruct.energy1中所有元素的均值mean1、最小值min1、最大值max1、最小值在energy1中的位置pos1；步骤三十二、计算第一范围能量列表的能量差littlediff＝(mean1-min1)/(max1-min1)；步骤三十三、计算p.enstruct.energy2中所有元素的均值mean2、最小值min2、最大值max2、最小值在energy2中的位置pos2；步骤三十四、计算第二范围能量列表的能量差bigdiff＝(mean2-min2)/(max2-min2)；步骤三十五、如果littlediff大于bigdiff则转到步骤三十六，否则转到步骤三十七；步骤三十六、p.enstruct.energypos＝pos1，转到步骤三十八；p.enstruct.energypos表示p.enstruct的energypos属性；步骤三十七、p.enstruct.energypos＝pos2；步骤三十八、计算过程结束。其他步骤及参数与具体实施方式四相同。具体实施方式六：本实施方式与具体实施方式五不同的是：如图5所示，步骤五中统计category的值发生变化的像元的个数并存储到transformcounter中具体为，对于image中的每一个像元，执行如下操作：步骤五一、获得像元p的p.enstruct.energypos的值；步骤五二、如果p.enstruct.energypos不等于p.enstruct.category则转到步骤五三，否则转到步骤五四；步骤五三、transformcounter＝transformcounter+1；步骤五四、p.enstruct.category＝p.enstruct.energypos。其他步骤及参数与具体实施方式五相同。具体实施方式七：本实施方式与具体实施方式四不同的是：步骤二四中，dis的计算公式为：dis＝norm(temppv-center[centercounter])其中norm为计算一个矢量的欧几里得范数(euclideannorm)。其他步骤及参数与具体实施方式四相同。具体实施方式八：本实施方式与具体实施方式五不同的是：步骤三五中，计算littleenergy的公式如下：littleenergy＝0.8×e(-m1)+0.2×e(-m2)其中e为自然对数，其中，tl1[i]是tl1中的第i个元素，i∈[1,n1]；因为tl1＝p.enstruct.enlist1；enlist1为第一范围邻域enstruct列表，可知tl1[i]表示像元p的第一范围邻域内某一像元的enstruct，因此tl1[i]同样包含enstruct的属性。由此可知，tl1[i].row表示tl1[i]的row属性，tl1[i].col表示tl1[i]的col属性，tl1[i].pv表示tl1[i]的pv属性，tl1[i].category表示tl1[i]的category属性。其他步骤及参数与具体实施方式五相同。具体实施方式九：本实施方式与具体实施方式五不同的是：步骤三六中，计算bigenergy的公式如下：bigenergy＝0.2×e(-m3)+0.8×e(-m4)其中e为自然对数，其中，tl2[j]表示tl2的第j个元素，j∈[1,n2]；因为tl2＝p.enstruct.enlist2；enlist2为第二范围邻域enstruct列表，可知tl2[j]表示像元p的第二范围邻域内某一像元的enstruct，因此tl2[j]同样包含enstruct的属性。由此可知，tl2[j].row表示tl2[j]的row属性，tl2[j].col表示tl2[j]的col属性，tl2[j].pv表示tl2[j]的pv属性，tl2[j].category表示tl2[j]的category属性。其他步骤及参数与具体实施方式五相同。实施例采用以下实施例验证本发明的有益效果：本实施例所述一种基于范围能量差的高分辨率遥感影像聚类方法按照以下步骤进行：将50张地面分辨率为2m～9cm的高分辨率遥感影像分别进行本发明的聚类方法、k均值算法、em、分层次聚类和群智能算法进行聚类(如图6所示为进行聚类的分辨率为9cm的高分辨率遥感影像)，得到不同的处理结果。如图7所示为采用本发明方法对图6进行聚类的结果，可以看出，本发明方法成功的将高分辨率遥感影像中的树木、道路、房屋区分出来了，聚类结果明显优于其他方法。分别对以上五种聚类方法得到的处理结果进行分析，比较结果如下表所示：本发明方法k均值算法em分层次聚类法群智能算法杂点较少较多较多较多较多准确率94％65％73％77％81％本发明还可有其它多种实施例，在不背离本发明精神及其实质的情况下，本领域技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。当前第1页12