基于对称初始化和聚集度优化的亚像元定位方法和系统与流程

本发明涉及遥感技术，更具体地的，涉及遥感亚像元定位和遥感影像分类技术。

背景技术

遥感影像中普遍存在的混合像元现象在一定程度上阻碍了地物遥感分类与识别精度的提高。对于混合像元分解，需要获取的参数是混合像元内的地物类别，地物所占比例以及地物空间分布情况，对应的求解技术分别为端元提取技术、混合像元分解技术以及亚像元定位技术(胡茂桂和王劲峰2014，凌峰等2011)。亚像元定位技术又称超分辨率制图技术，是一种特殊的遥感影像分类技术，其理论于1997年由atkinson提出(凌峰等2011)，并将其定义为一种将软分类转换成更高空间尺度上的硬分类技术。此外，亚像元定位技术也不同于传统的硬分类技术。硬分类技术将占像元内比例最高的地物类型分配给混合像元，生成硬分类结果；而亚像元定位技术则将混合像元切割成更小单元—亚像元，并根据地类比例信息将地物类型相应的分配到亚像元中，制图结果无论在精度还是效果上都优于硬分类结果。

从已有亚像元定位技术来看，亚像元定位效果的优劣主要在于地物空间分布特征描述是否准确。地物尺寸和遥感影像像元分辨率存在以下两种关系，且两种关系下地物空间分布特征的描述是不同的。当地物尺寸大于像元分辨率时，目标地物由多个亚像元组成，地物空间分布特征的描述主要是基于空间相关性的(陈跃红和葛咏2012)，即像元内以及像元间不同地物的空间分布存在相关性，距离相距较近的像元/亚像元与距离相距较远的像元/亚像元相比更可能属于同一地类。一般情况下，基于空间相关性理论的亚像元定位技术对各地类面积重建较准确且亚像元定位总体精度较高。而当地物尺寸小于像元分辨率时，目标地物只占混合像元的一部分，混合像元内斑块破碎程度较高，空间结构复杂，亚像元定位相对困难(林皓波等2011，任武和葛咏2011)。此时地物空间分布特征的描述不再基于空间相关性理论，而是依据空间模式拟合，构建其理论模型的重点在于如何有效借助先验信息描述各种地物空间分布模式。较具代表性的基于空间模式拟合的亚像元定位模型包括两点直方图模型(atkinson2008)、神经网络模型(tatemetal.2008)、地统计插值模型(张旸和陈沈良2005，boucher&kyriakidis2007，atkinsonetal.2008，王立国等2014，杜奕和张挺2016)以及景观异质性模型(frazier2015，huetal.2015)。其中，地统计插值模型、景观异质性模型以及一些神经网络模型求解计算量很大且求解过程繁琐，对于数据量和计算能力的要求很高。因此，对于小尺寸地物亟需提出一种计算量较小、求解过程简单的亚像元定位模型。

同时，上述大部分亚像元定位模型仅能针对两种地物尺寸和像元分辨率关系之一，而较少有模型能兼顾两者并存的情况。如基于空间相关性的亚像元定位模型主要针对地物尺寸大于像元分辨率的情况，亚像元制图结果趋于聚集，对于线形地物的表现能力较差，混合像元内部独立的小斑块也会被忽略掉，从而影响亚像元定位效果和制图质量。基于空间模式拟合的亚像元定位模型主要针对地物尺寸小于像元分辨率的情况，过多的依靠训练影像结构信息，而忽略像元间空间相关性，从而不能获得较高的亚像元定位精度。而在遥感影像中通常包括上述两种情况，因此如何兼顾两种尺寸地物进行亚像元定位以提高亚像元定位精度和制图质量是亚像元定位模型研究中亟待解决的技术问题之一。

除了对亚像元定位理论的研究外，如何对模型进行求解也是亚像元定位研究中另一个重要研究内容。现有求解算法大致可以分为直接求解算法和迭代求解算法两种类型。直接求解算法是指通过计算直接确定每个亚像元所属地物类型，求解过程中不需要进行迭代运算；而迭代求解算法一般先随机赋予每个亚像元某种地物属性，然后通过一定的规则迭代改变其地物属性从而获得最终制图结果。两类求解算法的区别主要在于是否需要进行迭代。直接求解算法不需要进行迭代，算法整体计算效率高，运行速度快，如在空间引力模型中可直接利用混合像元分解结果基于像元间引力排序对亚像元地类直接进行分配，模型运算速度快效率高(mertensetal.2006)。直接求解算法虽在运算效率上具有明显优势，但不是所有亚像元定位模型中都可使用直接求解算法，使得直接求解算法的应用具有一定局限性。迭代求解算法虽然运行速度偏慢，但通过迭代通常可获得更好的空间相关性，最终获得更高的亚像元定位精度，因此使用迭代求解算法可有效提高亚像元定位精度和制图效果。

但在迭代求解算法中，随机亚像元初始化分布的空间相关性较差，且迭代过程中混合像元内所有不同地类的亚像元都需要进行像元交换，这将导致亚像元定位模型进行过多的迭代运算并耗费较长运算时间。所以，如何改进迭代求解算法，也是需要解决的技术问题之一。

技术实现要素：

基于以上讨论，针对现有模型及算法难以兼顾两种尺寸地物进行亚像元定位的不足，以及进一步提高亚像元定位精度和制图效率的需求，本发明提出一种基于对称初始化和聚集度优化的亚像元定位方法，本方法引入像元聚集度概念模拟小尺寸地物空间分布，结合基于空间引力计算的对称亚像元初始化分布，开展基于对称初始化和聚集度优化的亚像元定位方法研究。

本发明提出一种基于对称初始化和聚集度优化的亚像元定位方法，包括：

s1，按照地物类型依次计算实验遥感影像中的混合像元与其8邻域中相邻混合像元之间的空间引力值；

s2，初始化亚像元空间分布，基于步骤s1所计算的空间引力值，将初始化地物类型分配给中心混合像元中的亚像元；

s3，基于混合像元内亚像元的空间分布，计算某一地物类型亚像元的像元聚集度，然后基于不同类型地物亚像元的像元聚集度，得到混合像元的像元聚集度；同样方法得到先验信息数据的像元聚集度；

s4，基于s3所计算的像元聚集度对对称区域的不同地类亚像元进行像元交换。

本发明还提出一种基于对称初始化和聚集度优化的亚像元定位系统，包括：处理器和存储器，所述存储器存储有处理器可执行的程序代码，所述程序代码被执行时，所述处理器实现如下步骤：

s1，按照地物类型依次计算实验遥感影像中的混合像元与其8邻域中相邻混合像元之间的空间引力值；

s2，初始化亚像元空间分布，基于步骤s1所计算的空间引力值，将初始化地物类型分配给中心混合像元中的亚像元；

s3，基于混合像元内亚像元的空间分布，计算某一地物类型亚像元的像元聚集度，然后基于不同类型地物亚像元的像元聚集度，得到混合像元的像元聚集度；同样方法得到先验信息数据的像元聚集度；

s4，基于s3所计算的像元聚集度对对称区域的不同地类亚像元进行像元交换。

本发明的有益效果包括：

1、本发明的方法在理论上兼顾两种尺寸地物，并进一步提高亚像元定位精度和运算效率。

2、本发明方法中使用对称初始化亚像元分布，只有对称区域的不同地类亚像元需要进行像元交换，有效降低亚像元定位所需时间，大大提高亚像元定位效率。同时，像元聚集度既能反映混合像元中亚像元的数量，又能反映亚像元的空间分布情况，是有效借助先验信息描述小尺寸地物空间分布模式的手段之一，经过像元聚集度优化有效提高了亚像元定位结果中小尺寸地物的定位效果和制图质量。

3、本发明的方法应用于我国重要农业区中高空间分辨率多光谱遥感影像中，与原始基于聚集度优化的亚像元定位模型进行对比，进行定位效果和制图精度的验证，以期进一步提高亚像元定位精度和亚像元级制图质量，为亚像元定位以及遥感影像分类提供模型理论依据和方法技术借鉴。

附图说明

图1为本发明的方法的一个实施方式的流程图。

图2为本发明的方法的技术路线图。

图3为空间引力计算示意图。

图4为中心混合像元及其相邻混合像元的示意图。

图5为亚像元赋值顺序的示意图。

图6为初始化后的亚像元分布图。

图7为验证本发明的方法的研究区概况图。

图8-图11显示了亚像元定位实验结果，其中图8为实验的参考影像图，图9为实验的k-mean分类结果图，图10为实验的聚集度优化亚像元定位结果图，图11为本发明的方法的实验结果图。

图12为本发明和其他方法的亚像元定位精度对比表。

图13为实验的技术路线图。

具体实施方式

下面参照附图描述本发明的实施方式，其中相同的部件用相同的附图标记表示。在不冲突的情况下，下述的实施例及实施例中的技术特征可以相互组合。

如图1-2所示，本发明的亚像元定位方法包括：

s1，按照地物类型依次计算实验遥感影像中的混合像元与其8邻域中相邻混合像元之间的空间引力值。

像元间空间引力是描述像元空间相关性的重要指标，在空间引力模型中如何评定像元间空间引力大小是模型研究重点之一。本发明引入万有引力描述像元间空间引力，即通过万有引力描述像元空间相关性。由万有引力定律可知，空间内任何两个有质量的物体间都存在着引力关系，且引力大小与物体质量乘积成正比，与物体间距离的平方成反比。因此，在空间引力模型和万有引力定律的基础上，本发明将空间内存在引力关系的两个像元假设为两个质点，则万有引力定律中的质量乘积就可以描述像元的丰度乘积，物体间距离的平方就可用于表达亚像元的空间分布情况。像元间引力计算原理具体如下：

通过端元选择模型和混合像元分解模型得到某中心混合像元pm中含有ω类地物且第z类地物丰度值为z(p)，pn为pm8邻域内的相邻混合像元。将pm分割为s×s个亚像元，则pm中亚像元与pn之间的空间引力值为：

其中，z(pm)为中心混合像元pm中z类地物的丰度值；z(pn)为与pm8邻域内相邻混合像元pn中z类地物的丰度值；i为pm中的亚像元xi的脚标，i＝1,2,3,…,s²；k为pn中z类地物的数目；j为pn中属于z类地物的亚像元xj的脚标，j＝1,2,3,…,k；rij为亚像元xi与亚像元xj间的欧式距离。设两像元/亚像元间的引力存在于两像元/亚像元的中心之间，则rij的计算是从一个像元/亚像元的中心到另一个像元/亚像元的中心。

其中，(mi,ni)为xi的坐标，(mj,nj)为xj的坐标。

其中，z(ωi)为中心混合像元pm与其8邻域内所有地物类型为z的相邻混合像元亚像元间的总引力值。

在s1中，根据式(3)按照地物类型依次计算中心混合像元pm与其8邻域内所有相邻混合像元间的空间引力值。此时没有对亚像元分配地物类型，因此在计算空间引力时需将混合像元pm和pn看作为一个像元，则式(3)简化为式(4)。空间引力计算示意图如图2所示。

其中，rmn是混合像元pm与混合像元pn间的欧氏距离。设两混合像元间的引力存在于两混合像元的中心之间，则rmn的计算从一个混合像元的中心到另一个混合像元的中心。

s2，初始化亚像元空间分布，基于步骤s1所计算的空间引力值，将初始化地物类型分配给中心混合像元中的亚像元。

若亚像元随机分配地物类型，则混合像元内所有不同地类的亚像元都需要进行像元交换，这将导致过多的迭代并耗费较长运算时间。为有效降低亚像元定位所需时间，本发明使用对称模式初始化混合像元中的亚像元空间分布。在这种初始化方式下，只有对称区域的不同地类亚像元需要进行像元交换，亚像元定位效率将大大提高。

在s2中，基于步骤s1所计算的空间引力值，将初始化地物类型分配给中心混合像元中的亚像元。设中心混合像元pm的坐标为(i,j)，则从pm8邻域左上角混合像元开始，以顺时针方向，其相邻混合像元可依次设为a，b，c，d，e，f，g，h。如图4所示，左侧的图为遥感影像中的一部分，其中包括中心混合像元pm及其8邻域混合像元，右侧的图是混合像元中的亚像元。相邻混合像元a到h的坐标分别为(i-1,j-1),(i-1,j),(i-1,j+1),(i,j+1),(i+1,j+1),(i+1,j),(i+1,j-1),(i,j-1)。

首先，在s21中，比较混合像元pm与其8邻域中所有相邻混合像元间的空间引力值z(ωi)的大小，并得到最大空间引力值z(ωi)max的坐标。然后，在s22中，根据上述最大空间引力值的坐标和亚像元赋值顺序(见图4)，确定亚像元的赋值顺序得到亚像元的初始化空间分布。

优选地，为更好的形成较对称的亚像元空间分布，在初始化赋值时按照丰度值从小到大的顺序对混合像元中的地类进行赋值。经过初始化赋值后，混合像元pm中属于z类地物的亚像元将靠近含有z类地物数目最多的相邻混合像元集中分布。

步骤3：基于混合像元内亚像元的空间分布，计算某一地物类型亚像元的像元聚集度，然后基于不同类型地物亚像元的像元聚集度，得到混合像元的像元聚集度；同样方法得到先验信息数据的像元聚集度。进一步，基于先验信息数据中亚像元的像元聚集度和实验遥感影像中亚像元的像元聚集度，建立聚集度目标函数，以供后续迭代用。

下面进行详细描述，其中，像元聚集度计算方法如下：

设x＝(f(i,j))m×n是遥感影像中某一混合像元内的亚像元空间分布，其中以灰度值f(i,j)表示亚像元(i,j)处的地物类型。当地物类型为c，即f(i,j)＝c时，若f(i+△i,j+△j)＝c，则亚像元(i,j)处c类地物的像元聚集度可以表示为：

式(5)对所有满足-i+1≤△i≤m-i，-j+1≤△j≤n-j的△i，△j求和。设灰度值为c的亚像元共有k个，则c类地物的像元聚集度为：

显然，jc的大小不仅与地物类型为c的亚像元个数相关，也与这些亚像元的空间分布情况相关。设duv表示地物类型为c的第u个点和第v个点之间的距离，则

两个亚像元的像元聚集度不仅与它们之间的距离相关，还与它们的地物类型相关。亚像元(i,j)处的像元聚集度可以表示为：

其中，k为常数，为x的像元聚集度，k是x所含亚像元的个数。

设duv表示第u个点和第v个点之间的距离，δuv表示这两个亚像元的灰度值之差，则

当k很大时，若f(i+△i,j+△j)≠f(i,j)，则

亚像元(i,j)处的像元聚集度可以简化为：

设x中含有k1个c1类地物亚像元，含有k2个c2类地物亚像元，……，kn个cn类地物亚像元。则k＝k1+k2+……+kn，则该混合像元的像元聚集度可以表示为

在s31中，根据式(10)分别计算先验信息数据(较高空间分辨率遥感影像，地面样方数据)和实验遥感影像中亚像元的像元聚集度，并建立目标函数：

其中，是先验信息数据的像元聚集度，是实验遥感影像的像元聚集度。

步骤4：基于s3所计算的像元聚集度对对称区域的不同地类亚像元进行像元交换。

为有效节省亚像元定位所需时间，本发明仅需对对称区域的不同地类亚像元进行交换就能获得最终亚像元定位结果，大大提高了亚像元定位效率。

s41，根据亚像元初始化后地物空间分布的对称性，选取属于不同地类的亚像元pi及处于对称区域的亚像元pj。具体地，亚像元初始化后的地物空间分布是关于水平轴、垂直轴、主对角线以及副对角线对称的，对称轴如图6虚线所示。混合像元pm被这4条虚线划分为完全相等的8个部分，每个部分从左上角开始按顺时针方向分别标记为i到viii，相邻区域的亚像元以虚线对称。选取亚像元pi及其以虚线对称的亚像元pj，且pi和pj分别属于不同地类，设pi地物类型为c1，pj地物类型为c2；

s42，交换亚像元pi与亚像元pj的地物类型，若交换后目标函数值变小，则交换亚像元pi与亚像元pj的地物类型，即pi地物类型为c2，pj地物类型为c1，并更新像元聚集度目标函数。

s43，对混合像元pm中的每个亚像元执行步骤s41-步骤s42。

步骤s5，遍历实验遥感影像中每个混合像元，并执行步骤s4。

步骤s6：使用迭代算法，重复步骤s4-步骤s5，迭代至不再进行像元交换后获得最终亚像元定位结果。

为验证本发明的方法的有效性，以及针对真实遥感影像进行亚像元定位的适用性，以中国河北省石家庄市(图7灰色区域)为研究区进行亚像元定位实验，研究区覆盖面积为15km×15km。该研究区主要作物种植制度为冬小麦—夏玉米一年两熟制。其中，夏收作物主要以冬小麦为主(播种面积占夏收作物总播种面积的90％以上)，秋收作物较为复杂主要为玉米、蔬菜、棉花和豆类等。为减少亚像元定位结果的不确定性，本发明选取种植结构相对简单且以冬小麦为主的夏收作物为研究对象，开展研究区多光谱遥感影像亚像元定位实验。

实验数据为sentinel-2卫星10m空间分辨率多光谱遥感影像，影像成像时间为2017年4月18日，数据合成波段为2，3，4，8a。截取影像研究区覆盖区域1500像元×1500像元，通过对影像进行几何校正、影像配准和融合等处理，获得了研究区数字正射影像。为避免端元选择和混合像元分解过程中引入额外误差，且保证实验结果仅关注所提对称初始化—聚集度优化亚像元定位模型的有效性，本发明将k-mean分类结果直接作为参考影像来进行亚像元定位精度和制图效果的验证，参考影像如图8所示。

此外，为保证亚像元定位研究中参考影像的质量，开展了对sentinel-2遥感影像k-mean分类结果的验证工作。研究中，采用系统抽样与随机抽样相结合的方式，获取了sentinel-2卫星过境时30个地面样方的地物分布调查结果，每个样方面积不低于200m×200m。利用30个地面样方对实验影像k-mean分类结果进行验证，验证结果表明实验遥感影像k-mean分类结果总体精度和kappa系数分别为95.90％和0.939，参考影像的质量满足亚像元定位研究中的精度要求。

本发明的方法的验证过程如下：

首先，使用k-mean分类方法对原始遥感影像进行分类，并将分类结果作为评价亚像元定位精度和制图效果的参考影像；其次，对原始遥感影像进行退化处理并获得实验遥感影像。为便于精度评价，退化尺度与重建尺度相同。若退化尺度为s，则将原影像中每s×s个像元值加权平均一次，加权平均值作为退化影像中对应空间位置的新像元值；然后，使用各影像分类方法或亚像元定位模型对实验遥感影像进行亚像元定位，为客观评价亚像元定位模型的有效性，直接对参考影像中像元进行统计，获得混合像元中各类地物丰度值；最后，使用参考影像验证各亚像元定位精度和制图效果。实验具体技术路线如图13所示。

为验证亚像元定位模型定位精度，本发明选用总体精度和kappa系数作为精度评价指标。遥感影像中纯净像元可能会高估亚像元定位精度，特别是当纯净像元在遥感影像中占主导地位时。为了减少这种不利影响，本发明引入混合像元的总体精度和kappa系数对亚像元定位结果进行评价。其中，混合像元总体精度计算公式如下：

式中oamixed为混合像元总体精度，m0为被正确定位的混合像元总数，m为混合像元总数。总体精度oamixed越高，表示亚像元定位精度越高。

混合像元kappa系数计算公式如下

kappa＝(p0-pe)/(1-pe)(13)

式中为亚像元定位一致率，pii＝aii/n，aii为亚像元定位一致数，n为混合像元总数，为期望一致率，即亚像元定位结果中的混合像元和参考影像中的混合像元由于偶然机会所造成的一致率，pi＝ri/n，βi＝ci/n，ri，ci分别为第i个亚像元所对应的行合计和列合计。kappa系数计算结果通常在[0,1]区间内，kappa系数越高说明亚像元定位精度越高。

本发明以10m空间分辨率sentinel-2多光谱遥感影像为原始遥感影像，分别基于k-mean分类方法，原始聚集度优化模型以及本发明开展了20m、30m、40m、50m和60m空间分辨率影像亚像元定位研究，亚像元定位目标影像空间分辨率均为10m，即重建尺度分别为2，3，4，5和6。具体实验结果如图8-12所示。从图12可以看出，两种亚像元定位模型定位结果精度均优于k-mean分类结果精度。对比两种亚像元定位模型定位结果，本发明定位精度均优于原始聚集度优化亚像元定位精度。当重建尺度为4时，本发明亚像元定位总体精度、kappa系数分别为69.88％和0.595，分别高于k-mean分类结果18.30％和0.224，分别高于原始聚集度优化模型亚像元定位结果1.04％和0.014(见图12)。

如图8-11所示，参考影像(图8)中不同地类边界较为清晰，能够较好反映不同地类空间结构和细节信息。退化后影像k-mean分类结果(图9)丢失了部分空间结构和大部细节信息，因此分类结果变得较模糊，这种不利影响对线宽较小的线形地物作用更加突出。通过原始聚集度优化模型和本发明进行亚像元定位后，可重建因影像退化所丢失的地物空间结构和部分细节信息，且本发明的亚像元定位效果更接近于参考影像。通过对比硬分类结果和原始聚集度优化模型亚像元定位结果，证明了本发明的有效性和利用真实遥感影像进行亚像元定位的适用性。

为提高亚像元定位模型对于不同尺寸地物适用性，并在一定程度上提高亚像元定位精度，本发明首次引入像元聚集度概念模拟小尺寸地物空间分布；结合基于空间引力计算的对称亚像元初始化分布，开展基于对称初始化和聚集度优化的亚像元定位模型研究，以期使新模型在理论上兼顾两种尺寸地物，并进一步提高模型定位精度和运算效率。利用sentinel-2多光谱遥感影像开展了基于k-mean分类方法，原始聚集度优化模型以及本发明的亚像元定位实验及其精度验证研究。结果表明，当重建尺度为4时，本发明亚像元定位整体精度、kappa系数分别为69.88％和0.595，分别高于k-mean分类结果18.30％和0.224，分别高于原始聚集度优化模型亚像元定位结果1.04％和0.014。可见，本发明均优于k-mean分类方法和原始聚集度优化模型，证明本发明进行遥感影像亚像元定位具有一定可行性和适用性，为亚像元定位以及遥感影像分类提供模型理论依据和方法技术借鉴。虽然本发明仅以sentinel-2多光谱遥感影像为例，阐述和实施了如何利用所提基于对称初始化和聚集度优化的亚像元定位模型进行亚像元定位的思路和方法，但该方法对于一般的多光谱遥感影像的亚像元定位具有普遍的适用性。

本发明的方法的定位效果和制图精度会受到亚像元丰度值的影响，在本文试验中为避免引入额外误差，客观评价亚像元定位模型有效性，使用了理想像元丰度值模拟混合像元分解模型输出的像元丰度值。在实际应用中，本发明的方法可与端元提取模型、混合像元分解模型联合使用，从而解决实际问题中的遥感影像分类和制图问题，为更好地突破遥感影像空间分辨率的限制提供了新的途径。

将本发明的方法应用于我国重要农业区中高空间分辨率多光谱遥感影像进行模型验证，其中，中高空间分辨率多光谱遥感影像为覆盖河北省石家庄市研究区的sentinel-2遥感影像。结果表明，当退化尺度为2-6时，利用本发明进行亚像元定位精度较k-mean分类方法分类精度，原始聚集度优化模型亚像元定位精度均有所提高。可见，本发明具有一定可行性和适用性，为亚像元定位以及遥感影像分类提供模型理论依据和方法技术借鉴。

以上所述的实施例，只是本发明较优选的具体实施方式，本领域的技术人员在本发明技术方案范围内进行的通常变化和替换都应包含在本发明的保护范围内。