CN108985154B - 基于影像聚集度的小尺寸地物亚像元定位方法和系统 - Google Patents

Abstract

一种基于影像聚集度的小尺寸地物亚像元定位方法和系统，所述方法包括：S1，根据混合像元分解模型计算的丰度值为试验遥感影像中的亚像元随机分配地类，以得到亚像元初始化结果；S2，基于先验信息和亚像元定位结果中亚像元的影像聚集度，建立影像聚集度目标函数；S3，优化S2中的目标函数；S4，遍历亚像元初始化结果中所有混合像元，对每个混合像元执行步骤S3；S5，使用迭代算法，重复步骤S3‑步骤S4。亚像元定位精度较K‑mean硬分类方法和像元交换亚像元定位模型制图精度和效果均有所提高，为更好地突破遥感影像空间分辨率的限制提供了新的思路。

Description

基于影像聚集度的小尺寸地物亚像元定位方法和系统

技术领域

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

背景技术

遥感影像中普遍存在着混合像元，它的存在已经成为进一步提高地物遥感分类与识别精度的主要障碍之一。针对混合像元分解问题，需要研究的是混合像元内像元地类属性、不同地类所占比例(丰度)以及其空间分布情况，分别对应的技术是端元提取技术、混合像元分解技术以及亚像元定位技术。从现有研究来看，对于前两个问题的研究较为深入，提出的理论模型也较多。但是，单纯依靠端元提取技术和混合像元分解技术仅能获得一个像元中各端元的比例信息，而无法获得混合像元内部各地类的分布情况，制图结果仍停留在低空间分辨率尺度上，需要依靠亚像元定位技术才能获得所需要的更高空间分辨率的制图结果。亚像元定位技术又称超分辨率制图技术，其理论于1997年由Atkinson提出，并将其定义为一种将软分类转换成更高空间尺度上的硬分类技术。此外，亚像元定位技术也不同于传统的硬分类技术。传统硬分类技术将占像元内比例最高的地物类型赋予像元属性，生成硬分类制图结果；而亚像元定位技术则将像元切割成更小单元—亚像元，并根据丰度信息将具体地物类型相应的分配到这些亚像元中，制图结果无论在精度还是效果上都优于硬分类制图结果。因此，亚像元定位技术作为混合像元分解的后续有效处理手段，可以使遥感制图结果达到亚像元级，同原始遥感影像硬分类结果相比，制图精度和效果也有较大提高，这将有利于对混合像元问题的掌控，继而有利于对于遥感影像的研究和后续应用。

目前，大多数亚像元定位模型主要针对大尺寸地物，即地物尺寸大于像元分辨率。这种情况下，目标地物由多个亚像元组成，地物空间分布特征的描述主要是基于空间相关性的，即像元内以及像元间不同地物的空间分布存在相关性，距离较近的像元/亚像元与距离较远的像元/亚像元相比更可能属于同一地类。一般情况下，基于空间相关性理论的亚像元定位模型对各地类的面积重建较准确且定位总体精度较高，但亚像元制图结果趋于聚集，混合像元内部独立的小斑块易被忽略掉，对于小尺寸地物的定位效果较差，从而影响亚像元定位的总体效果和制图质量。而在真实遥感影像中，尤其是中高空间分辨率遥感影像中，小尺寸地物大量存在，遥感影像中一些像元斑块破碎程度高，地物空间结构复杂，亚像元定位相对困难。此时地物空间分布特征的描述不能再基于空间相关性，而是依据空间模式拟合，构建其理论模型的重点在于如何有效借助先验信息描述各种地物空间分布模式。

除了对于亚像元定位理论的研究，如何对模型进行求解也是亚像元定位研究中另一个重要研究内容。现有求解算法大致可以分为直接求解算法和迭代求解算法两种类型。直接求解算法是指通过计算直接确定每个亚像元所属的地类，求解过程中不需要进行迭代。直接求解算法整体计算效率高，运行速度快，但通常来说一次计算难以获得混合像元内部准确的地物分布，因此直接求解算法在应用上具有一定局限性。迭代求解算法一般先随机赋予每个亚像元某种地类属性，然后通过一定的规则迭代改变其地类属性从而获得最终制图结果。迭代求解虽然需要迭代运行，运行速度偏慢，但通过迭代通常可获得更准确的亚像元空间分布，最终获得更高的亚像元定位精度，因此使用迭代求解算法可有效提高亚像元定位精度和制图效果。因此，现阶段亟需提出一种计算量较小、求解过程简单，针对小尺寸地物的亚像元定位模型。

发明内容

基于以上讨论，本发明首次引入影像聚集度概念模拟小尺寸地物空间分布，开展基于影像聚集度的亚像元定位新模型研究。并将所提新模型应用于我国重要农业区中高空间分辨率多光谱遥感影像中，与经典像元交换亚像元定位模型进行对比，验证所提新模型定位效果和制图精度，以期为更好地突破遥感影像空间分辨率的限制提供了新的思路。

本发明提出一种基于影像聚集度的小尺寸地物亚像元定位方法，包括：

S1，根据混合像元分解模型计算的丰度值为试验遥感影像中的亚像元随机分配地类，以得到亚像元初始化结果；

S2，基于先验信息和亚像元定位结果中亚像元的影像聚集度，建立影像聚集度的目标函数；

S3，优化S2中的目标函数；

S4，遍历亚像元初始化结果中所有混合像元，对每个混合像元执行步骤S3；

S5，使用迭代算法，重复步骤S3-步骤S4。

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

S1，根据混合像元分解模型计算的丰度值为试验遥感影像中的亚像元随机分配地类，以得到亚像元初始化结果；

S2，基于先验信息和亚像元定位结果中亚像元的影像聚集度，建立影像聚集度目标函数；

S3，优化S2中的目标函数；

S4，遍历亚像元初始化结果中所有混合像元，对每个混合像元执行步骤S3；

S5，使用迭代算法，重复步骤S3-步骤S4。

本发明的有益效果包括：

亚像元定位精度较K-mean硬分类方法和像元交换亚像元定位模型制图精度和效果均有所提高。具有一定可行性和适用性，为更好地突破遥感影像空间分辨率的限制提供了新的思路。

附图说明

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

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

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

图4-图7显示了亚像元定位试验结果，其中图4为试验的参考影像图，图5为试验的K-mean分类结果图，图6为像元交换模型亚像元定位结果图，图7为本发明的方法的实验结果图。

图8为亚像元定位试验精度图。

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

具体实施方式

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

影像聚集度既能反映混合像元中亚像元的数量，又能反映亚像元的空间分布情况，是有效借助先验信息描述地物空间分布模式的手段之一。本发明首次引入影像聚集度概念模拟小尺寸地物空间分布，提出了基于影像聚集度的亚像元定位方法。本发明的方法的原理图如图1-2所示。

S1，亚像元空间分布初始化：根据混合像元分解模型计算的丰度值为试验遥感影像中的亚像元随机分配地类，以得到亚像元初始化结果。

S2，计算影像聚集度目标函数：基于先验信息(较高空间分辨率遥感影像)和亚像元定位结果(初始化结果)中亚像元的影像聚集度，建立影像聚集度目标函数。

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

式(1)对所有满足-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，则该混合像元的影像聚集度可以表示为

根据式(6)分别计算先验信息(较高空间分辨率遥感影像)和亚像元定位结果(初始化结果)中亚像元的影像聚集度，并建立基于影像聚集度的亚像元定位模型目标函数：

其中，

是先验信息图像的影像聚集度，

是试验遥感影像的影像聚集度；

S3，优化S2中的目标函数。

S31，选取试验遥感影像中某一混合像元p_m中的亚像元p_i和亚像元p_j，且亚像元p_i和亚像元p_j分属于不同类型地物，则设p_i属性为C₁，p_j属性为C₂；

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

S33，对混合像元p_m中的每个亚像元执行步骤S31-步骤S32。

S4：遍历混合像元，对试验遥感影像中的每个混合像元执行步骤S3。

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

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

试验数据为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分类结果(见图8)总体精度和kappa系数分别为95.90％和0.939，参考影像的质量满足亚像元定位研究中的精度要求。

试验的技术路线如图9所示。本发明的方法的验证过程如下：

首先，使用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所示，从图8可以看出，两种亚像元定位模型的亚像元定位精度均优于K-mean硬分类结果。对比两种亚像元定位模型亚像元定位结果，影像聚集度模型亚像元定位精度均优于像元交换模型亚像元定位精度。当重建尺度为4时，本发明的结果整体精度、kappa系数分别为68.84％和0.581，分别高于K-mean硬分类精度17.26％和0.210，分别高于像元交换模型亚像元定位精度3.66％和0.026。(见图8)。

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

影像聚集度既能反映混合像元中亚像元的数量，又能反映亚像元的空间分布情况，是有效借助先验信息描述地物空间分布模式的手段之一。针对现有亚像元定位模型主要基于空间相关性理论难以模拟小尺寸地物空间分布情况的现状，为了进一步提高亚像元定位效果和制图精度，本发明在亚像元定位模型中首次引入影像聚集度概念，提出了基于影像聚集度的亚像元定位方法。本发明以Sentinel-2多光谱遥感影像为例开展了基于K-mean硬分类方法，像元交换亚像元定位模型和本文所提影像聚集度亚像元定位模型的亚像元定位试验以期对所提新模型的亚像元定位精度和制图效果进行验证。结果表明，所提新模型的亚像元定位结果优于其它两种制图方法的制图结果，更接近于参考影像；当重建尺度为4时，本发明方法的亚像元定位结果整体精度、kappa系数分别为68.84％和0.581，分别高于K-mean硬分类精度17.26％和0.210，分别高于像元交换模型亚像元定位精度3.66％和0.026。可见，本发明的方法均优于K-mean硬分类方法和像元交换亚像元定位模型，证明了本发明的方法进行遥感影像亚像元定位具有一定可行性和适用性。虽然本发明仅以Sentinel-2多光谱遥感影像为例，阐述和实施了如何利用所提影像聚集度亚像元定位模型进行亚像元定位的思路和方法，但该方法对于一般的多光谱遥感影像的亚像元定位具有普遍的适用性。

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

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

Claims (8)

1.一种基于影像聚集度的小尺寸地物亚像元定位方法，其特征在于，包括：

S1，根据混合像元分解模型计算的丰度值为试验遥感影像中的亚像元随机分配地类，以得到亚像元初始化结果；

S2，基于先验信息和亚像元定位结果中亚像元的影像聚集度，建立影像聚集度目标函数，

其中，影像聚集度的含义如下：

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

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

影像聚集度目标函数为：

其中，

是先验信息图像的影像聚集度，

是试验遥感影像的影像聚集度；

S3，优化S2中的目标函数，包括：

S31，选取试验遥感影像中某一混合像元p_m中的亚像元p_i和亚像元p_j，且亚像元p_i和亚像元p_j分属于不同类型地物，则设p_i属性为C₁，p_j属性为C₂；

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

S33，对混合像元p_m中的每个亚像元执行步骤S31-步骤S32；

S4，遍历亚像元初始化结果中所有混合像元，对每个混合像元执行步骤S3；

S5，使用迭代算法，重复步骤S3-步骤S4。

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

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

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

S31，选取试验遥感影像中某一混合像元p_m中的亚像元p_i和亚像元p_j，且亚像元p_i和亚像元p_j分属于不同类型地物，则设p_i属性为C₁，p_j属性为C₂；

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

S33，对混合像元p_m中的每个亚像元执行步骤S31-步骤S32。

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

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

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

S1，根据混合像元分解模型计算的丰度值为试验遥感影像中的亚像元随机分配地类，以得到亚像元初始化结果；

S2，基于先验信息和亚像元定位结果中亚像元的影像聚集度，建立影像聚集度目标函数

其中，影像聚集度的含义如下：

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

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

影像聚集度目标函数为：

其中，

是先验信息图像的影像聚集度，

是试验遥感影像的影像聚集度；

S3，优化S2中的目标函数，包括：

S31，选取试验遥感影像中某一混合像元p_m中的亚像元p_i和亚像元p_j，且亚像元p_i和亚像元p_j分属于不同类型地物，则设p_i属性为C₁，p_j属性为C₂；

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

S33，对混合像元p_m中的每个亚像元执行步骤S31-步骤S32；

S4，遍历亚像元初始化结果中所有混合像元，对每个混合像元执行步骤S3；

S5，使用迭代算法，重复步骤S3-步骤S4。

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

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

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

S31，选取试验遥感影像中某一混合像元p_m中的亚像元p_i和亚像元p_j，且亚像元p_i和亚像元p_j分属于不同类型地物，则设p_i属性为C₁，p_j属性为C₂；

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

S33，对混合像元p_m中的每个亚像元执行步骤S31-步骤S32。

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

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

Images