CN103559506B - 一种基于矢量边界的亚像元制图方法 - Google Patents

Abstract

本发明公开了一种基于矢量边界的亚像元制图方法，包括利用混合像元自身及其邻域像元的组分（端元）面积矢量或称丰度，提取混合像元内地物的精确矢量边界，再将其转换成给定放大尺度的栅格（亚像元），进而完成影像的亚像元制图。本发明具有实用性广、计算速度快、模拟精度高等优点，适用于遥感数据进行目标提取、国土资源调查、全球气候变化、灾害监测与灾情评估等地学信息提取和地学数据挖掘工作。

Description

一种基于矢量边界的亚像元制图方法

技术领域

本发明涉及一种亚像元制图方法，属于地球空间信息技术领域。

技术背景

利用数字遥感影像制作专题地图（影像分类）是遥感技术研究中的热点领域。然而在遥感影像获取过程中，由于地表的复杂性、外部环境的影响以及传感器自身的局限性，使得影像中普遍存在混合像元(Mixed Pixel)。混合像元的光谱特征是所对应的不同土地覆被类型光谱响应的综合特征，在遥感分类过程中使传统的硬分类方法产生误差与不确定性。因此，软分类（SoftClassification）方法被提出，即以百分比（软信息、分数影像）的形式描述各地物类别在像元中的面积比例，但在像元内部的空间位置仍未知。因此，又提出了亚像元制图（Sub-pixel Mapping or Super-resolution Mapping），实现了比原始输入空间分辨率更精细的硬分类图（Atkinson，1997）。亚像元制图能够进一步细化和精确影像分类结果，使遥感技术应用从像元尺度进入到亚像元尺度，而亚像元制图也成为必要且有意义的研究内容。

亚像元制图方法可以分成两大类（Atkinson，2009）：一类是回归学习类算法，如地统计学的方法（Boucher和Kyriakidis，2006，2007）和基于BP神经网络的方法（Zhang等，2002）；另一类是空间相关性最大类算法，如基于Hopfield神经网络的方法（Tatem等，2001）和空间吸引力模型的方法（Mertens，2006）。回归学习类算法计算速度快，并且可以提供不确定性的信息，但是除了空间吸引力模型外都需要辅助数据，如矢量数据、先验知识等。而在空间相关性最大类算法中，亚像元可以和亚像元做比较，并在预测过程中可以适当的提供支持。

Ge等在2009年提出了一种用多边形来表示混合像元内各类别空间位置的方法，命名为GSPM，其中多边形的顶点位置由混合像元及其邻域像元的各组分比例确定的。GSPM与以上两类方法相比包括两大优点：（1）计算仅依赖于原始影像软分类结果，不需要辅助数据；（2）计算速度快，不需要反复迭代。但仍然有待改善：（1）当缩放因子为奇数时，在确定混合像元内多边形顶点位置时会引入较多的舍入误差。（2）对于存在多个土地覆被类型的混合像元，不同土地覆被类型的界限可能是不一致的，虽然允许微调，但仍可能影响输出精度。针对上述问题，本发明提出一种基于矢量边界的亚像元制图方法（VBSPM），延续了GSPM用多边形获取亚像元内几何划分的想法，并进一步提出了一种矢量边界提取模型以减少舍入误差和不同土地覆被类之间边界的不一致的问题。对于沿原始遥感影像边缘丢失的结果，VBSPM利用衰减函数以获取这些混合像元的补充边界信息，用于预测亚像元的分布。

发明内容

本发明解决的技术问题：提供一种基于矢量边界的亚像元制图方法，该方法利用混合像元自身及其邻域像元的软分类信息，先计算像元内各组分类别的精确矢量边界，再确定给定放大因子下的亚像元的类别信息，获得较像元尺度的分类方法更精细和准确的地物分布信息，从而提高遥感影像分类的分辨率和分类精度。

本发明的技术方案：一种基于矢量边界的亚像元制图方法包括如下步骤：

步骤1、预处理遥感数据，通过软分类获取影像上每个像元的软分类信息，即在像元内各组分所占的面积比例，也称为丰度；

步骤2、利用线性衰减函数填充边界像元在分数影像范围外的邻域像元的软分类信息；

步骤3、根据混合像元的软分类信息，提取地物的精确边界，主要过程为：首先估算每条线段的长度和位置，确定初始矢量边界多边形，然后对多边形的进行拓扑调整，以减少地物边界间的缝隙和重叠，获得最终的各类别地物的矢量多边形；

步骤4、设定放大尺度因子，将步骤3中提取出的矢量边界进行多边形转换，用射线法（ray-crossing）判断亚像元中心点与步骤3中矢量多边形的位置关系，进而判断每个亚像元的类别信息。

所述步骤3中利用混合像元自身及其邻域像元的软分类信息，计算每个像元内部各地物类别的亚像元空间分布：假如混合像元内部有两类地物C₁和C₂，首先分别计算各类别的邻域8个像元以及中心混合像元在垂直和水平方向上对中心混合像元的贡献，即每像元在中心混合像元A中对应线段的长度和位置，连接线段端点形成初始多边形。然后调整类别间不完全重合的顶点，使得各类别的线段端点尽量重合，以减少不同类别多边形间的重叠和缝隙。然后连接调整后的由线段确定的边界点，得到中心混合像元内各类别的矢量多边形，即地物C₁和C₂的模拟分布。所述步骤4中，设定放大因子，利用射线法判断亚像元中心点与各地物类别矢量多边形的拓扑关系，根据中心点所在的矢量多边形的类别属性确定亚像元的类别属性。

本发明与现有技术相比的优点在于：本发明利用空间相关性提出一种混合像元内部地物制图方法，从而在亚像元的级别上提供了地物类型的边界信息。与传统方法相比，此算法计算简单、解释性较强、算法效率高，且不需要其他辅助数据；与GSPM相比，此方法利用矢量边界提取模型，不受放大尺度因子的影响，也减少了不同土地覆被类别边界间的缝隙和重叠，从而获取更为精确的矢量边界信息，进而减少亚像元制图的误差。

附图说明

图1为本发明的主流程图；

图2为中心像元及其邻域像元示意图；

图3为矢量边界提取实例图，其中a是中心像元及其邻域像元的组分比例图，b是初始多边形边界图，c是最终的地物边界图。

具体实施方式

如图1所示，本发明的具体实施步骤如下：

步骤1、预处理遥感数据，通过软分类获取影像上每个像元的软分类信息，即在像元内各组分所占的面积比例，也称为丰度。

假设在一个3×3像元的区域包含两种土地覆被类型C₁和C₂。通过软分类获得中心像元P_m及其邻域像元P_i的两种土地覆被类型所占的面积比例分别为P_m{0.7,0.3},P₀{0.35,0.65},P₁{0.95,0.05}P₂{1,0}P₃{1,0}P₄{1,0}P₅{1,0}P₆{0.65,0.35}和P₇{0,1}，如图3中的(a)所示。

步骤2、利用线性衰减函数填充边界像元在分数影像范围外的邻域像元的软分类信息。

为了处理遥感影像边界地区的混合像元，引入了衰减函数进行增强，补充在分数影像范围外的邻域像元。衰减函数是可以有效获取边界地区像元的8个邻域像元的值，补充基于空间相关性的亚像元制图方法的边界信息。在此发明中，用线性衰减函数来填充边缘像元在分数影像之外的邻域像元丰度值，从而避免信息丢失。

步骤3、根据混合像元的软分类信息，提取地物的精确边界，主要过程为：首先估算每条线段的长度和位置，确定初始矢量边界多边形，然后对多边形的进行拓扑调整，以减少地物边界间的缝隙和重叠，获得最终的各类别地物的矢量多边形。

（1）估计每条边界线段的长度和位置

混合像元内每条线段的长度和位置的估算公式如下：

L_ijk＝S_ik,    (1)

${\mathrm{Pos}}_{\mathrm{ijk}}={\mathrm{Pos}}_{\mathrm{ijk}}^{\left(0\right)}+\left\{\begin{array}{cc}(1-S_{\mathrm{ik}})\·\frac{S_{(i+1)k}}{S_{(i+1)k}+S_{(i-1)k}},& S_{(i+1)k}\≠S_{(i-1)k}\\ (1-S_{\mathrm{ik}})\·\frac{1}{2},& S_{(i+1)k}=S_{(i-1)k}\end{array}\right.---\left(2\right)$

其中，L_ijk表示像元i(i=0,1,..,7,c)对应于线段j(j=AB,BCD,..,KLA,MN,PQ)的第k(k=1,2,..,K)类地物的长度，如图2所示；S_ik表示像元i上第k类地物的面积比（软分类结果）；S_(i+1)k，S_(i-1)k分别表示像元i左右两侧像元上第k类地物的面积比（软分类结果）；表示像元i对应于线段j的第k类地物的初始位置；Pos_(i+1)k表示用S_(i+1)k，S_(i-1)k计算出的平移线段的终止位置。

举例说明公式（2）如何计算代表类别C₁的多边形边界的长度和位置。在多边形边界AB上的线段是V₁V₂，其中|V₁V₂,|=P₀₁|AB|=0.35|AB|；在多边形BCD上的边界是V₃CV₄，|V₃CV₄|=0.95|BCD|；以此类推。中心像元内部的多边形顶点由线段T₁T₃和T₄T₂确定，其中|T₁T₃|=0.7|PQ|，|T₄T₂|=0.7|MN|。

（2）获取每个土地覆被类型对应的多边形

沿逆时针方向连接各顶点得到初始多边形，要求一定要通过PQ和MN。例如，C₁类别的初始多边形顶点是V₂、T₁、V₃、C、V₄、T₂、V₅、V₆、V₇、F、V₈、T₃、V₉、V₁₀、T₃、V₁₁、V₁₂、T₄、T₁和V₁，其中顶点T₁、T₂、T₃和T₄在像元内部。其中，部分顶点重合，如B和V₂，D和V4、V5，M和T₂等，如图3中的（b）所示。为了消除不同类别的顶点不完全重合可能带来的多边形的缝隙和重叠，对部分顶点进行删除，调整规则如下：

规则一：当向量V_jT_i和向量T_iV_j+1间的夹角小于设定的阈值θ，并且顶点V_j和V_j+1不在同一条线段上时，删除顶点T_i。

规则二：当向量V_jT_i和向量T_iV_j+1间的夹角小于设定的阈值θ，并且顶点V_j和V_j+1在同一条线段上时，删除顶点V_j和V_j+1。

然后，初始多边形中的一些顶点被删除。如图3中的（b）所示，顶点T2和T3被删除。然后，沿逆时针顺序连接删除后剩余的标记顶点，形成最终的边界多边形，如图3中的（c）所示。

步骤4、设定放大尺度因子，将步骤3中提取出的矢量边界进行多边形转换，用射线法（ray-crossing）判断亚像元中心点与步骤3中矢量多边形的位置关系，进而判断每个亚像元的类别信息。

在提取出混合像元内的多边形和设定放大尺寸因子后，选用射线法确定各亚像元的类别。以亚像元的中心为端点，任意画一条射线，计算射线与多边形边界的交点数。如果是奇数，表示中心点在多边形内，赋予亚像元与多边形相同的类别标记。如果是偶数，则表示中心点在多边形外，对其他多边形重复上述过程直至确定亚像元类别信息。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本邻域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的权利要求范围之内。

Claims (2)

1.一种基于矢量边界的亚像元制图方法，其特征在于步骤如下：

步骤1、获取基于矢量边界亚像元制图方法的输入数据：对遥感影像进行预处理，通过软分类获取影像上每个像元的软分类信息，即像元内各组分所占的面积比例，也称为丰度；

步骤2、补充边缘信息：为避免边缘像元的信息丢失，利用线性衰减函数填充边缘像元在影像范围外的邻域像元的软分类信息；

步骤3、提取各地物在混合像元内的矢量边界：首先根据公式利用混合像元的软分类信息估算矢量边界线段的长度和位置，确定初始矢量多边形，然后对矢量多边形进行拓扑调整，以减少地物边界间的缝隙和重叠，获得最终的各类别地物的矢量多边形；

步骤4、确定混合像元内各亚像元的类别属性：设定放大尺度因子，将步骤3中提取出的矢量边界进行矢量多边形转换，用射线法判断亚像元中心点与步骤3中矢量多边形的拓扑关系，进而判断每个亚像元的类别属性。

2.根据权利要求1所述的基于矢量边界的亚像元制图方法，其特征在于：所述步骤3中利用混合像元自身和邻域像元的软分类信息，计算对应的中心混合像元边界线段的长度L_ikj和位置Pos_ijk ⁽⁰⁾，具体公式如下：

L_ijk＝S_ik

${\mathrm{Pos}}_{\mathrm{ijk}}={\mathrm{Pos}}_{\mathrm{ijk}}^{\left(0\right)}+\left\{\begin{array}{cc}(1-S_{\mathrm{ik}})\·\frac{S_{(i+1)k}}{S_{(i+1)k}+S_{(i-1)k}},& S_{(i+1)k}\≠S_{(i-1)k}\\ (1-S_{\mathrm{ik}})\·\frac{1}{2},& S_{(i+1)k}=S_{(i-1)k}\end{array}\right.$

其中，L_ikj表示像元i对应于线段j的第k类地物的长度；S_ik表示像元i中第k类地物的软分类信息；S_(i+1)k和S_(i-1)k分别表示像元i左右两侧像元中第k类地物的软分类信息；Pos_ijk ⁽⁰⁾表示像元i中对应于线段j的第k类地物的初始位置；Pos_ijk表示用S_(i+1)k和S_(i-1)k计算出的平移线段的终止位置；

然后，根据线段确定各类别的初始矢量多边形，再通过拓扑调整减少各类别初始矢量多边形之间的缝隙和重叠，获得最终各类别地物的矢量多边形和矢量边界。