CN108921035A - 基于空间引力和像元聚集度的亚像元定位方法和系统 - Google Patents

Abstract

本发明提供一种基于空间引力和像元聚集度的亚像元定位方法和系统，所述方法包括：S1，基于遥感影像混合像元中各类型地物的丰度值，得到中心混合像元中分属各类地物的亚像元与8邻域相邻混合像元中相同地物类型亚像元的总引力值；基于所计算的引力值，将初始化地物类型分配给中心混合像元中的亚像元；S2，对S1中的亚像元分布进行像元聚集度优化，包括：S21，计算像元聚集度目标函数；S22，交换像元以优化亚像元空间分布并更新S21中的目标函数；S23，遍历遥感影像中所有混合像元，对每个混合像元执行步骤S22；S24，使用迭代算法，重复步骤S22‑步骤S23。本发明方法和系统计算的亚像元定位结果既具有较高空间相关性，又具有有效地物空间分布，提高了亚像元定位模型对于不同尺寸地物的适用性。

Description

基于空间引力和像元聚集度的亚像元定位方法和系统

技术领域

本发明涉及遥感技术，更具体地的，涉及一种亚像元定位方法。

背景技术

遥感影像中普遍存在的混合像元现象已经成为进一步提高地物识别与遥感分类精度的主要障碍之一。针对混合像元分解，需要研究的是混合像元内像元地类属性、丰度及其空间分布，分别对应的技术是端元提取技术、混合像元分解技术以及亚像元定位技术。从现有研究来看，对于前两个技术的研究较多也较为深入。但是，只依靠端元提取技术和混合像元分解技术仅能获得混合像元内各端元比例信息，而无法获得各类地物分布，制图结果仍停留在较低空间分辨率尺度上，需要依靠亚像元定位技术才能获得所需要的更高空间分辨率的制图结果。亚像元定位技术作为混合像元分解的后续有效处理手段，是一种将软分类转换成更高空间尺度上的硬分类技术。可以使遥感制图结果达到亚像元级，这将有利于对混合像元问题的掌控，继而有利于提高地物识别与遥感分类精度。

从亚像元定位的理论来看，亚像元定位成功与否关键在于地物空间分布特征描述是否准确。而不同地物尺寸和遥感影像像元分辨率关系下，地物空间分布特征的描述也是不同的。地物尺寸和遥感影像像元分辨率之间存在以下两种情况：地物尺寸大于像元分辨率(高分辨率制图，H型)和地物尺寸小于像元分辨率(低分辨率制图，L型)。当地物尺寸大于像元分辨率时，目标地物由多个亚像元组成，地物空间分布特征的描述主要是基于空间相关性的，即像元内以及像元间不同地物的空间分布存在相关性，距离相距较近的像元/亚像元与距离相距较远的像元/亚像元相比更可能属于同一地类。一般情况下，基于空间相关性理论的亚像元定位技术对各地类的面积重建较准确且总体定位精度较高，较具代表性的基于空间相关性的亚像元定位模型包括空间引力模型、矢量边界模型、神经网络模型、马尔可夫随机场模型、模拟退火模型、元胞自动机模型以及遗传算法模型等。而当地物尺寸小于像元分辨率时，目标地物只占据混合像元内一部分，混合像元内斑块破碎程度较高，空间结构复杂，亚像元定位相对困难。此时地物空间分布特征的描述不再基于空间相关性理论，而是依据空间模式拟合，构建其理论模型的重点在于如何有效借助先验信息描述各种地物空间分布。较具代表性的基于空间模式拟合的亚像元定位模型包括两点直方图模型、神经网络模型、地统计插值模型以及景观异质性模型等。

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

发明内容

基于以上背景技术中存在的问题，本发明提供一种基于空间引力和像元聚集度的亚像元定位方法，包括：

S1，基于遥感影像混合像元中各类型地物的丰度值，得到中心混合像元中分属各类地物的亚像元与8邻域相邻混合像元中相同地物类型亚像元的总引力值；基于所计算的引力值，将初始化地物类型分配给中心混合像元中的亚像元；

S2，对S1中的亚像元分布进行像元聚集度优化，包括：

S21，计算像元聚集度目标函数；

S22，交换像元以优化亚像元空间分布并更新S21中的目标函数；

S23，遍历遥感影像中所有混合像元，对每个混合像元执行步骤S22；

S24，使用迭代算法，重复步骤S22-步骤S23。

本发明提供一种基于空间引力和像元聚集度的亚像元定位系统，包括：处理器和存储器，所述存储器存储有处理器可执行的程序代码，所述程序代码被执行时，所述处理器实现如下步骤：

S1，基于遥感影像混合像元中各类型地物的丰度值，得到中心混合像元中分属各类地物的亚像元与8邻域相邻混合像元中相同地物类型亚像元的总引力值；基于所计算的引力值，将初始化地物类型分配给中心混合像元中的亚像元；

S2，对S1中的亚像元分布进行像元聚集度优化，包括：

S21，计算像元聚集度目标函数；

S22，交换像元以优化亚像元空间分布并更新S21中的目标函数；

S23，遍历遥感影像中所有混合像元，对每个混合像元执行步骤S22；

S24，使用迭代算法，重复步骤S22-步骤S23。

本发明的有益效果包括：

本发明在空间引力初始化的基础上，利用像元聚集度进行优化以得到最终亚像元定位结果。亚像元定位结果既具有较高空间相关性，又具有有效地物空间分布，提高了亚像元定位模型对于不同尺寸地物的适用性。

本发明的方法应用于我国重要农业区中高空间分辨率多光谱遥感影像中，与原始空间引力模型和像元聚集度模型进行对比，对所提新模型进行定位效果和制图精度的验证；本发明的方法提高了亚像元定位精度和亚像元级制图质量，为更好地突破遥感影像空间分辨率的限制提供了新的思路。

附图说明

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

图2为本发明的方法的一个实施方式的技术路线图。

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

图4-图8显示了亚像元定位实验结果，其中图4为实验的参考影像图，图5为实验的K-mean分类结果图，图6为空间引力模型亚像元定位结果图，图7为像元聚集度模型亚像元定位结果图，图8为本发明的方法的实验结果图。

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

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

具体实施方式

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

本发明方法包括两部分：空间引力初始化和像元聚集度优化。本发明的方法的原理图如图1-2所示。

步骤S1，基于遥感影像混合像元中各类型地物的丰度值，得到中心混合像元中分属各类地物的亚像元与8邻域相邻混合像元中相同地物类型亚像元的总引力值，通过归一化，以像元为单位计算遥感影像中每一类地物的归一化引力值。

具体地，由混合像元分解模型得到某中心混合像元p_m中含有G类地物且第d类地物丰度值为f(p)。将p_m分割为s×s个亚像元，则该混合像元中的亚像元与相邻混合像元之间的引力值为：

其中为中心混合像元p_m中d类地物的丰度值；为与p_m相邻混合像元p_n中d类地物的丰度值；i为p_m的亚像元x_i的脚标，i＝1,2,3,…,s²；k为p_n中d类地物的数目；j为p_n中属于d类地物的亚像元x_j的脚标，j＝1,2,3,…,k；为p_m中所有地物类型为d的亚像元与p_n中地物类型为d的亚像元间的引力值；R_ij为从亚像元x_i到亚像元x_j的欧式距离。

其中(m_i,n_i)为x_i的坐标，(m_j,n_j)为x_j的坐标。

其中为中心混合像元p_m中所有地物类型为d的亚像元与8邻域混合像元中地物类型为d的亚像元间的总引力值。

具体地，步骤S1包括步骤S11-S13：

步骤S11，计算空间引力，根据式(1)依次计算中心混合像元p_m中分属各类地物的亚像元与相邻混合像元中相同地物类型亚像元的引力值z(ω_in)，然后依据式(3)计算p_m中分属各类地物的亚像元与8邻域相邻混合像元中相同地物类型亚像元的总引力值z(ω_i)。

步骤S12：引力归一化，设遥感影像中共有G类地物，则G类地物的加权引力值为：

对属于d类地物的总引力值z(ω_i)进行归一化，则归一化后的d类地物的引力值为：

步骤S13：初始化亚像元地物类型，依次以像元为单位计算遥感影像中每一类地物的归一化引力值并对比不同地类归一化引力值的大小，若对于亚像元x第e类地物的归一化引力值在所有地类归一化引力值中最大，即：

则将地物类型e分配给亚像元x。

步骤S14，遍历整幅遥感影像，将所得亚像元分布作为亚像元初始化分布。

步骤S2，对S1中的亚像元初始化分布进行像元聚集度优化，具体原理如下。

设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类地物的像元聚集度可以表示为：

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

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

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

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

设d_uv表示第u个点和第v个点之间的距离，δ_uv表示这两个亚像元的影像值之差，则：

当K很大时，若f(i+△i,j+△j)≠f(i,j)，则亚像元(i,j)处的像元聚集度可以简化为：

设X中含有k₁个C₁类地物亚像元，含有k₂个C₂类地物亚像元，……，k_n个C_n类地物亚像元。则k＝k₁+k₂+……+k_n，则该混合像元的像元聚集度可以表示为：

具体地，步骤S2包括步骤S21-S24。

步骤S21：计算像元聚集度目标函数。在空间引力初始化分布基础上，根据式(12)分别计算较高空间分辨率遥感影像和亚像元初始化结果(亚像元定位结果)中的像元聚集度，并建立目标函数：

其中，是较高空间分辨率遥感影像的像元聚集度，是亚像元初始化结果的像元聚集度；

步骤S22：交换像元以优化亚像元空间分布并更新S21中的目标函数；

步骤S221，选取亚像元初始化结果中某一混合像元p_m中的亚像元p_i和亚像元p_j，且亚像元p_i和亚像元p_j分属于不同类型地物，设p_i属性为C₁，p_j属性为C₂；

步骤S222，交换亚像元p_i与亚像元p_j的地物类型，若交换后目标函数值变小，则交换亚像元p_i与亚像元p_j的地物类型，即p_i属性为C₂，p_j属性为C₁，并更新像元聚集度目标函数；

步骤S223，对混合像元p_m中的每个亚像元执行步骤S221-步骤S222。

步骤S23：遍历遥感影像中所有混合像元，对每个混合像元执行步骤S22；

步骤S24：迭代算法，重复步骤S22-步骤S23、迭代H次后获得最终亚像元定位结果。更具体地，迭代至不再进行像元交换后获得最终亚像元定位结果。

为验证所提亚像元定位模型的可行性和有效性，本文以中国北方黄淮海平原区内的河北省石家庄市(图1灰色区域)为研究区进行亚像元定位实验，研究区覆盖面积为15km×15km。该研究区主要作物种植制度为冬小麦-夏玉米一年两熟制。其中，夏收作物主要以冬小麦为主(播种面积比例占夏收作物面积的90％以上)，秋收作物较为复杂主要为玉米、蔬菜、棉花以及豆类等。为减少实验结果的不确定性，本发明选取种植结构相对简单且以冬小麦为主的夏收作物为研究对象，开展研究区多光谱遥感影像亚像元定位实验。研究区范围如图3所示。

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

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

如图9所示，本发明的方法的验证过程如下：

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

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

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

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

kappa＝(p₀-p_e)/(1-p_e) (9)

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

本发明以10m空间分辨率Sentinel-2多光谱遥感影像为原始遥感影像，分别基于K-mean分类方法，空间引力亚像元定位模型，像元聚集度亚像元定位模型以及所提空间引力-像元聚集度亚像元定位模型开展了亚像元定位实验。实验遥感影像空间分辨率分别为20m、30m、40m、50m和60m，目标影像空间分辨率为10m，即亚像元重建尺度分别为2，3，4，5和6。在像元聚集度模型和空间引力-像元聚集度模型中，迭代次数设置为10。具体研究结果如图4至图8以及图10所示，从图10可以看出，各亚像元定位精度均优于K-mean分类精度。对比各亚像元定位结果，本发明的方法定位精度均优于空间引力模型和像元聚集度模型亚像元定位精度。当重建尺度为4时，所提新亚像元定位模型亚像元定位整体精度、kappa系数分别为69.50％和0.590，分别高于K-mean分类结果17.92％和0.219，分别高于空间引力模型亚像元定位结果2.94％和0.023，分别高于像元聚集度模型亚像元定位结果0.66％和0.009。

参考影像图4中不同类别地物边界较为清晰，能够较好反映不同地类空间结构和细节信息。退化后影像K-mean分类结果(图5)丢失了部分空间结构和大部细节信息，因此分类结果变得较模糊，这种不利影响对线宽较小的线型地物作用更为突出。通过空间引力模型，像元聚集度模型以及所提空间引力-像元聚集度模型进行亚像元定位后，可重建因影像退化所丢失的地物空间结构和部分细节信息，且所提空间引力-像元聚集度模型的亚像元定位效果更接近于参考影像。通过对不同亚像元定位模型的比较，证明了所提空间引力-像元聚集度模型的可行性和有效性。

为提高亚像元定位模型对于不同尺寸地物适用性，本发明提出一种基于空间引力和像元聚集度的亚像元定位模型，以期使新模型兼顾不同尺寸地物并在一定程度上提高模型定位精度和制图效果。利用Sentinel-2多光谱遥感影像开展了基于K-mean分类方法，空间引力模型，像元聚集度模型以及所提空间引力-像元聚集度模型的亚像元定位实验及其精度验证研究。结果表明，当重建尺度为4时，本发明的亚像元定位整体精度、kappa系数分别为69.50％和0.590，分别高于K-mean分类结果17.92％和0.219，分别高于空间引力模型亚像元定位结果2.94％和0.023，分别高于像元聚集度模型亚像元定位结果0.66％和0.009。可见，本发明亚像元定位结果均优于K-mean分类方法，空间引力模型以及像元聚集度模型，证明本发明进行遥感影像亚像元定位具有一定可行性和有效性，这为更好地突破遥感影像空间分辨率的限制提供了新的思路。虽然本发明仅以Sentinel-2多光谱遥感影像为例，阐述和实施了如何利用所提基于空间引力和像元聚集度的亚像元定位模型进行亚像元定位的思路和方法，但该方法对于一般的多光谱遥感影像的亚像元定位具有普遍的适用性。

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

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

Claims (10)

1.一种基于空间引力和像元聚集度的亚像元定位方法，其特征在于，包括：

S1，基于遥感影像混合像元中各类型地物的丰度值，得到中心混合像元中分属各类地物的亚像元与8邻域相邻混合像元中相同地物类型亚像元的总引力值；基于所计算的引力值，将初始化地物类型分配给中心混合像元中的亚像元；S2，对S1中的亚像元分布进行像元聚集度优化，包括：

S21，计算像元聚集度目标函数；

S22，交换像元以优化亚像元空间分布并更新S21中的目标函数；

S23，遍历遥感影像中所有混合像元，对每个混合像元执行步骤S22；

S24，使用迭代算法，重复步骤S22-步骤S23。

2.根据权利要求1所述的亚像元定位方法，其特征在于，步骤S1包括：

S11，依次计算中心混合像元p_m中分属各类地物的亚像元与相邻混合像元中相同地物类型亚像元的引力值z(ω_in)，然后计算p_m中分属各类地物的亚像元与8邻域相邻混合像元中相同地物类型亚像元的总引力值z(ω_i)；

S12，将S11中计算的d类地物的总引力值z(ω_i)进行归一化；

S13，初始化亚像元地物类型，依次以像元为单位计算遥感影像中每一类地物的归一化引力值并对比不同地类归一化引力值的大小，若对于亚像元x第e类地物的归一化引力值在所有地类归一化引力值中最大，则将地物类型e分配给亚像元x；

S14，遍历整幅遥感影像，将所得亚像元分布作为亚像元初始化分布。

3.根据权利要求2所述的亚像元定位方法，其特征在于，

在S21中，基于较高空间分辨率遥感影像和亚像元初始化结果中的像元聚集度，来建立像元聚集度目标函数。

4.根据权利要求3所述的亚像元定位方法，其特征在于，步骤S22包括：

S221，选取亚像元初始化结果中某一混合像元p_m中的亚像元p_i和亚像元p_j，且亚像元p_i和亚像元p_j分属于不同类型地物；

S222，交换亚像元p_i与亚像元p_j的地物类型，若交换后目标函数值变小，则交换亚像元p_i与亚像元p_j的地物类型；

S223，对混合像元p_m中的每个亚像元执行步骤S221-步骤S222。

5.根据权利要求1所述的亚像元定位方法，其特征在于，

在S24中，迭代至不再进行像元交换后获得最终亚像元定位结果。

6.一种基于空间引力和像元聚集度的亚像元定位系统，其特征在于，包括：处理器和存储器，所述存储器存储有处理器可执行的程序代码，所述程序代码被执行时，所述处理器实现如下步骤：

S1，基于遥感影像混合像元中各类型地物的丰度值，得到中心混合像元中分属各类地物的亚像元与8邻域相邻混合像元中相同地物类型亚像元的总引力值；基于所计算的引力值，将初始化地物类型分配给中心混合像元中的亚像元；

S2，对S1中的亚像元分布进行像元聚集度优化，包括：

S21，计算像元聚集度目标函数；

S22，交换像元以优化亚像元空间分布并更新S21中的目标函数；

S23，遍历遥感影像中所有混合像元，对每个混合像元执行步骤S22；

S24，使用迭代算法，重复步骤S22-步骤S23。

7.根据权利要求6所述的亚像元定位系统，其特征在于，步骤S1包括：

S11，依次计算中心混合像元p_m中分属各类地物的亚像元与相邻混合像元中相同地物类型亚像元的引力值z(ω_in)，然后计算p_m中分属各类地物的亚像元与8邻域相邻混合像元中相同地物类型亚像元的总引力值z(ω_i)；

S12，将S11中计算的d类地物的总引力值z(ω_i)进行归一化；

S13，初始化亚像元地物类型，依次以像元为单位计算遥感影像中每一类地物的归一化引力值并对比不同地类归一化引力值的大小，若对于亚像元x第e类地物的归一化引力值在所有地类归一化引力值中最大，则将地物类型e分配给亚像元x；

S14，遍历整幅遥感影像，将所得亚像元分布作为亚像元初始化分布。

8.根据权利要求7所述的亚像元定位系统，其特征在于，

在S21中，基于较高空间分辨率遥感影像和亚像元初始化结果中的像元聚集度，来建立像元聚集度目标函数。

9.根据权利要求8所述的亚像元定位系统，其特征在于，步骤S22包括：

S221，选取亚像元初始化结果中某一混合像元p_m中的亚像元p_i和亚像元p_j，且亚像元p_i和亚像元p_j分属于不同类型地物；

S222，交换亚像元p_i与亚像元p_j的地物类型，若交换后目标函数值变小，则交换亚像元p_i与亚像元p_j的地物类型；

S223，对混合像元p_m中的每个亚像元执行步骤S221-步骤S222。

10.根据权利要求6所述的亚像元定位系统，其特征在于，

在S24中，迭代至不再进行像元交换后获得最终亚像元定位结果。