CN102565810B - 一种遥感影像上土地利用地物边界轮廓提取方法 - Google Patents

Abstract

一种遥感影像上土地利用地物边界提取方法，将遥感影像与GIS数据配准，以GIS中的土地利用对象边界作为影像上土地利用地物边界轮廓提取的初始值，建立轮廓节点和网格的对应关系；统计网格内影像的区域信息，获取网格内光谱灰度信息统计值；构建网格能量最小化函数，将得到的属性引入到网格能量函数中；移动任意节点，计算此局部节点位置变化是否能减少网格总能量，能则改变节点位置，否则不变；检查各节点间的距离并满足一定条件的节点间插入新的节点，进行轮廓的多级演化，在循环达到一定条件前，返回上步继续循环。本发明在GIS数据引导下，基于网格能量最小化的动态轮廓方法，实现整体最优的高分辨率遥感影像土地利用边界特征提取。

Description

一种遥感影像上土地利用地物边界轮廓提取方法

技术领域

本发明属于遥感地物提取技术领域，涉及一种遥感影像上土地利用地物边界提取方法。

背景技术

遥感影像地物自动/半自动提取是遥感技术应用的重要内容之一，可以辅助工作者进行地物目标的判别，从而避免人工提取地物目标中大量繁琐的操作；但是遥感影像有着高度的复杂性，这些复杂性如光照条件、大气透明度以及目标物纹理变化、遮挡等造成了遥感影像自动目标提取的困难。而引入先验知识，尤其是过往的地理信息系统（GIS）数据能够为地物提取提供预设的地物形状与位置信息，从而提高地物提取的效率与准确性。

动态轮廓方法是一种已有GIS数据为导引时地物目标提取的常用方法，但是通常需要划分前景和背景，也即只能对单个目标进行演化，限制了此方法的应用范围。

发明内容

本发明的目的在于提供一种效率较高，并可用于多目标、多连通和复杂拓扑关系影像目标提取的遥感影像上的土地利用地物边界提取方法。

本发明的技术方案为一种遥感影像上土地利用地物边界提取方法，包括以下步骤：

步骤1，将遥感影像与GIS数据进行配准；

步骤2，GIS数据提供的土地上每一个土地利用对象作为一个网格叠置于遥感影像之上，对所有的网格区域进行标号，编号的值从1到                                                ， 

为网格区域的个数；建立遥感影像中网格轮廓的节点和网格的对应关系，形成节点链表和网格号码链表，节点链表记录网格轮廓的节点，网格号码链表记录节点所属网格的号码；

步骤3，对步骤2所得每一个网格内影像的区域信息进行统计，获取网格内的光谱灰度信息统计值T，实现步骤如下，

步骤3.1，遍历网格轮廓上所有节点的位置，设网格轮廓上某个节点的位置为

，找到所有节点的

取值中最大值

与最小值

,以及所有节点的

取值中最大值

与最小值，根据

、、、

确定每个网格的多边形区域外包矩形；

步骤3.2，绘制网格轮廓上的直线段，在绘制完成后跟踪网格轮廓的多边形，建立用于描述多边形的Freeman码；

步骤3.3，依据步骤3.2中建立的Freeman码，基于二维查表法建立轮廓值链表，所述轮廓值链表中记录网格轮廓上坐标为

的节点所对应的链表值为d(x,y)；

步骤3.4，依据步骤3.1建立的多边形区域外包矩形建立多边形外接矩形区域，通过对多边形外接矩形区域进行由上至下、由左至右的扫描获取多边形区域；进行扫描获取多边形区域时，在开始将扫描线经过的区域设置标记c(x,y)为0，若扫描线经过的链表值d(x,y)为-1，则设置标记c(x,y)为1，若扫描线再扫描经过的链表值d(x,y)为1，则再设置标记c(x,y)为0，如此反复进行；

步骤3.5，计算步骤3.4所得多边形区域的轮廓上第

行像素的光谱灰度信息

阶统计值

，计算公式如下

其中，

为坐标为的像素的光谱灰度值，为阶数；

步骤3.6，通过分别计算步骤3.4所得多边形区域内外像素的个数得到多边形区域内面积、多边形区域外面积

；根据步骤3.4所得标记d(x,y)与步骤3.5所得

阶统计值

，计算网格内的光谱灰度信息统计值T，计算公式如下

其中，

表示网格，

为网格

中第

个像素的坐标，

取值为1,2,…

，

为网格

中像素的总数，

为网格

的轮廓，

为网格内的光谱灰度信息一次项统计值，

为网格

内的光谱灰度信息平方项统计值；       

步骤4，根据网格内的光谱灰度信息统计值T计算所有网格的网格总能量

，计算网格总能量采用网格能量函数，所述网格能量函数公式如下

其中，为第

个网格中像素的总数，

取值为1,2,…

，

为网格的总数；第

个网格的光谱灰度信息统计值

；表示

的平方；为网格中像素的灰度；为网格轮廓参数，表示每个网格轮廓上每个节点的位置；

、

是第

个网格中的光谱灰度信息一次项统计值和光谱灰度信息平方项统计值；

步骤5，通过迭代使网络能量函数的值最小，包括执行以下子步骤，

步骤5.1，等概率在步骤2所得节点链表中随机选取某个节点，将该节点由初始的位置(x,y)移动到(x+dx,y+dy)，其中dx,dy的变化范围为（-d,d），而dx，dy的取值也是等概率地在范围（-d，d）中选取； 

步骤5.2，进行移动节点的交叉测试，如果测试发现轮廓由于节点移动出现了自相交，则放弃当前的移动并返回步骤5.1，否则进入步骤5.4；

步骤5.3，对于移动的节点，通过步骤2所得网格号码链表找到所有与移动的节点对应要改变的网格，根据步骤3重新获取所有与移动的节点对应要改变的网格内的光谱灰度信息统计值T，根据步骤4重新计算所有网格的网格总能量

，如果计算出的网格总能量

小于上次的则移动有效，保留移动结果后进入步骤5.4，否则放弃移动并进入步骤5.4；

步骤5.4，返回步骤5.1，并对步骤5.1、5.2、5.3进行迭代，直至没有可以移动的节点为止，进入步骤6； 

步骤6，判断是否满足结束条件，结束条件为相邻节点之间的距离是否小于或等于预设阈值，不满足则在距离大于预设阈值的相邻节点间插入新的节点，通过步骤2所得网格号码链表找到所有与新的节点对应要改变的网格，根据步骤3重新获取所有与新的节点对应要改变的网格内的光谱灰度信息统计值T，根据步骤4重新计算所有网格的网格总能量，直到满足结束条件，得到最终轮廓提取结果。

而且，当被处理的遥感影像为多通道影像，且通道的数目表示为M时，将步骤4的网格能量最小化函数变为多通道形式如下，

。

本发明根据已有的GIS中的土地数据，自动对目标进行提取。本发明利用网格能量最小化方法限制地物边界轮廓的演化，同时利用快速算法对影像目标上的信息进行快速统计，提高了提取的效率并拓宽了自动提取算法的应用范围，提高了实用性。

附图说明

图1是本发明实施例的流程图。

图2是本发明实施例的节点链表示意图。

图3是本发明实施例的网格链表示意图。

图4是本发明实施例的Freeman码示意图。

图5是本发明实施例的使用Freeman码对节点链编码示意图。

具体实施方式   

本发明技术方案可采用计算机软件技术实现自动运行流程。以下结合附图和实施例的流程详细说明本发明技术方案。

参见图1，本发明实施例提供的遥感影像上的土地利用地物边界提取方法包括以下步骤：

步骤1，将遥感影像与GIS数据进行配准。

实施例实现步骤1包括以下步骤：

Ⅰ、在影像与GIS土地利用地图上选取控制点；

Ⅱ、根据控制点之间的关系依据几何模型建立坐标变换方程后进行几何纠正。

配准的几何模型可采用RST模型、多项式模型以及Delaunay三角网模型等，几何纠正使用的重采样方法包括最邻近法、双线性法以及三次卷积法等，都为现有技术，本发明不予赘述。

步骤2，GIS数据提供的土地上每一个土地利用对象作为一个网格叠置于遥感影像之上，对所有的网格区域进行标号，编号的值从1到， 

为网格区域的个数；建立遥感影像中网格轮廓的节点和网格的对应关系，形成节点链表和网格号码链表，节点链表记录网格轮廓的节点，网格号码链表记录节点所属网格的号码。

实施例使用链表来表示节点与所属网格之间的联系，共建立两个链表：其一是节点链表，记录每个节点的信息。参见图2，节点链表中记录有节点0、节点1…节点n，每个节点都记录节点坐标、节点属性以及所属网格号码链表。其二是网格链表，记录每个网格的信息。参见图3，网格链表中记录有网格0、网格1…网格n，每个网格包含所有此网格内所有节点编号链表，以及网格所对应的累加量：像素数量及网格区域内的光谱灰度信息统计值T。

步骤3，对步骤2所得每一个网格内影像的区域信息进行统计，获取网格内的光谱灰度信息统计值T.

实施例对网格内影像的区域信息采用依据Green-Ostrogradsky公式采用以下快速算法进行统计，实现步骤如下：

步骤3.1，遍历网格轮廓上所有节点的位置，设网格轮廓上某个节点的位置为

，找到所有节点的

取值中最大值

与最小值

,以及所有节点的

取值中最大值

与最小值

，根据

、、

、确定每个网格的多边形区域外包矩形；具体确定方式为设多边形区域外包矩形的左上角位置坐标为,多边形区域外包矩形的右下角坐标为；

步骤3.2，绘制网格轮廓上的直线段，在绘制完成后跟踪网格轮廓的多边形，建立用于描述多边形的Freeman码；绘制网格轮廓上的直线段可采用现有的数值微分法。

步骤3.3，依据步骤3.2中建立的Freeman码，基于二维查表法建立轮廓值链表，所述轮廓值链表中记录网格轮廓上坐标为

的节点所对应的链表值为d(x,y)；如图5上轮廓的链表值d(x,y)为-1、-1、-1、0、1、1。

步骤3.4，依据步骤3.1建立的多边形区域外包矩形建立多边形外接矩形区域，通过对多边形外接矩形区域进行由上至下、由左至右的扫描获取多边形区域；具体实施时，步骤3.4中，进行扫描获取多边形区域时，可以在开始将扫描线经过的区域标记c(x,y)为0，若扫描线经过的链表值d(x,y)为-1，则设置标记c(x,y)为1，若扫描线再扫描经过的链表值d(x,y)为1，则再设置标记c(x,y)为0，如此反复进行。标记c(x,y)在3.6步骤中用于光谱灰度信息的统计，c（x,y）为1表示像素在轮廓围成的多边形区域内部，为0表示在多边形区域区域外部。

步骤3.5，计算步骤3.4所得多边形区域的轮廓上第

行像素的光谱灰度信息

阶统计值

，计算公式如下

其中，

为坐标为的像素的光谱灰度值，

为阶数，

的取值为1和2；

步骤3.6，通过分别计算步骤3.4所得多边形区域内外像素的个数得到多边形区域内面积、多边形区域外面积

；根据步骤3.4所得标记d(x,y)与步骤3.5所得阶统计值

，计算网格内的光谱灰度信息统计值T，计算公式如下

其中，

表示网格，

为网格

中第

个像素的坐标，

取值为1,2,…

，为网格

中像素的总数，

为网格

的轮廓，

为网格

内的光谱灰度信息一次项统计值，

为网格

内的光谱灰度信息平方项统计值。       

Freeman码为依线段的方向进行编码的一种现有链码编码技术，如图4按照0、1、2、3、4、5、6、7的顺序编码。依据步骤Ⅱ中建立的Freeman码，基于二维查表法建立轮廓值链表，下表为查表法所利用的Freeman码与轮廓值对应表，记轮廓上坐标为（x,y）的点的链表值为d(x,y)： 

 其中，out和in分别表示轮廓上（x,y）位置进入和发出的边缘线。

步骤4，根据网格内的光谱灰度信息统计值T计算所有网格的网格总能量

，计算网格总能量采用网格能量函数。

实施例构建的网格能量函数公式如下，

其中，

为第

个网格中像素的总数，

取值为1,2,…

，

为网格的总数；第

个网格的光谱灰度信息统计值

；

表示

的平方；

为网格中像素的灰度；

为网格轮廓参数，表示每个网格轮廓上每个节点的位置；

、

是第

个网格中的光谱灰度信息一次项统计值和光谱灰度信息平方项统计值。

步骤5，通过迭代使网络能量函数的值最小，包括执行以下子步骤，

步骤5.1，等概率在步骤2所得节点链表中随机选取某个节点，将该节点由初始的位置(x,y)移动到(x+dx,y+dy),其中dx,dy的变化范围为（-d,d），而dx，dy的取值也是等概率地在范围（-d，d）中选取； d的取值可由人工给定，通常设置为3个像素位置；

步骤5.2，进行移动节点的交叉测试，如果测试发现轮廓由于节点移动出现了自相交，则放弃当前的移动并返回步骤5.1，否则进入步骤5.4；测试可以通过判断以被移动节点为其中一个端点的轮廓上的直线段与轮廓上其他直线段是否相交的方法进行，直线段的相交算法为计算几何中常用方法；

步骤5.3，对于移动的节点，通过步骤2所得网格号码链表找到所有与移动的节点对应要改变的网格，根据步骤3重新获取所有与移动的节点对应要改变的网格内的光谱灰度信息统计值T，根据步骤4重新计算所有网格的网格总能量

，如果当前计算出的网格总能量

小于上次的则移动有效，保留此节点的移动后进入步骤5.4，否则放弃移动并进入步骤5.4；

步骤5.4，返回步骤5.1，并对步骤5.1、5.2、5.3进行迭代直至没有可以移动的节点为止，进入步骤6。步骤5通过不断的迭代使得网格能量最小化，在没有节点可以移动时就认为网络能量达到了最小。

步骤6，判断是否满足结束条件，结束条件为相邻节点之间的距离是否小于或等于预设阈值，不满足则在距离大于预设阈值的相邻节点间插入新的节点，通过步骤2所得网格号码链表找到所有与新的节点对应要改变的网格，根据步骤3重新获取所有与新的节点对应要改变的网格内的光谱灰度信息统计值T，直到满足结束条件，得到最终轮廓提取结果。预设阈值可以预先人工设定，一般可设置为5个像素距离。新的节点在距离大于预设阈值的相邻节点间等距插入，使节点之间距离小于阈值即可。

本发明方案也适用于当被处理的遥感影像为多通道影像，且通道的数目表示为M时，将步骤4的网格能量最小化函数变为多通道形式如下，

其中，m表示通道的标号，取值为1,2,…M。

即、

是第

个网格中第m通道的光谱灰度信息一次项统计值和光谱灰度信息平方项统计值

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (2)

1.一种遥感影像上土地利用地物边界提取方法，其特征在于，包括以下步骤：

步骤1，将遥感影像与GIS数据进行配准；

步骤2，GIS数据提供的土地上每一个土地利用对象作为一个网格叠置于遥感影像之上，对所有的网格区域进行标号，编号的值从1到R，R为网格区域的个数；建立遥感影像中网格轮廓的节点和网格的对应关系，形成节点链表和网格号码链表，节点链表记录网格轮廓的节点，网格号码链表记录节点所属网格的号码；

步骤3，对步骤2所得每一个网格内影像的区域信息进行统计，获取网格内的光谱灰度信息统计值T，实现步骤如下，

步骤3.1，遍历网格轮廓上所有节点的位置，设网格轮廓上某个节点的位置为(x,y)，找到所有节点的x取值中最大值x_max与最小值x_min,以及所有节点的y取值中最大值y_max与最小值y_min，根据x_max、x_min、y_max、y_min确定每个网格的多边形区域外包矩形；

步骤3.2，绘制网格轮廓上的直线段，在绘制完成后跟踪网格轮廓的多边形，建立用于描述多边形的Freeman码；

步骤3.3，依据步骤3.2中建立的Freeman码，基于二维查表法建立轮廓值链表，所述轮廓值链表中记录网格轮廓上坐标为(x,y)的节点所对应的链表值为d(x,y)；

步骤3.4，依据步骤3.1建立的多边形区域外包矩形建立多边形外接矩形区域，通过对多边形外接矩形区域进行由上至下、由左至右的扫描获取多边形区域；进行扫描获取多边形区域时，在开始将扫描线经过的区域设置标记c(x,y)为0，若扫描线经过的链表值d(x,y)为-1，则设置标记c(x,y)为1，若扫描线再扫描经过的链表值d(x,y)为1，则再设置标记c(x,y)为0，如此反复进行；

步骤3.5，计算步骤3.4所得多边形区域的轮廓上第y行像素的光谱灰度信息k阶统计值F_k[x,y]，计算公式如下

$F_k[x,y]=\underset{t=0}{\overset{x}{\mathrm{\Σ}}}{I}^k(t,y)$

其中，I为坐标为(t,y)的像素的光谱灰度值，k为阶数；

步骤3.6，通过分别计算步骤3.4所得多边形区域内外像素的个数得到多边形区域内面积N_a、多边形区域外面积N_b；根据步骤3.4所得标记d(x,y)与步骤3.5所得k阶统计值F_k[x,y]，计算网格内的光谱灰度信息统计值T，计算公式如下

$T=\left(\begin{array}{c}{T}_1\\ T_2\end{array}\right)=\left(\begin{array}{c}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}I(x_i,y_i)\\ \underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\left(I(x_i,y_i)\right)}^2\end{array}\right)$

$\underset{(x,y)\∈a}{\mathrm{\Σ}}{I}^k(x,y)=\underset{(x,y)\∈\∂C}{\mathrm{\Σ}}c(x,y)F_k[x,y]$

其中，a表示网格，(x_i,y_i)为网格a中第i个像素的坐标，i取值为1,2,…N，N为网格a中像素的总数，

为网格a的轮廓，T₁为网格a内的光谱灰度信息一次项统计值，T₂为网格a内的光谱灰度信息平方项统计值；

步骤4，根据网格内的光谱灰度信息统计值T计算所有网格的网格总能量J(I(x,y),θ)，计算网格总能量采用网格能量函数，所述网格能量函数公式如下

$J(I(x,y),\mathrm{\θ})=\underset{\mathrm{ar}=1}{\overset{R}{\mathrm{\Σ}}}{N}_{\mathrm{ar}}\log (T_{2\mathrm{ar}}/N_{\mathrm{ar}}-{(T_{1\mathrm{ar}}/N_{\mathrm{ar}})}^2)$

其中，N_ar为第ar个网格中像素的总数，ar取值为1,2,…R，R为网格的总数；第ar个网格的光谱灰度信息统计值I²(x,y)表示I(x,y)的平方；I(x,y)为网格中像素(x,y)的灰度；θ为网格轮廓参数，表示每个网格轮廓上每个节点的位置；T_1ar、T_2ar是第ar个网格中的光谱灰度信息一次项统计值和光谱灰度信息平方项统计值；

步骤5，通过迭代使网络能量函数的值最小，包括执行以下子步骤，

步骤5.1，等概率在步骤2所得节点链表中随机选取某个节点，将该节点由初始的位置(x,y)移动到(x+dx,y+dy)，其中dx,dy的变化范围为（-d,d），而dx，dy的取值也是等概率地在范围（-d，d）中选取，d的取值设置为3个像素位置；

步骤5.2，进行移动节点的交叉测试，如果测试发现轮廓由于节点移动出现了自相交，则放弃当前的移动并返回步骤5.1，否则进入步骤5.4；

步骤5.3，对于移动的节点，通过步骤2所得网格号码链表找到所有与移动的节点对应要改变的网格，根据步骤3重新获取所有与移动的节点对应要改变的网格内的光谱灰度信息统计值T，根据步骤4重新计算所有网格的网格总能量J(I(x,y),θ)，如果计算出的网格总能量J(I(x,y),θ)小于上次的则移动有效，保留移动结果后进入步骤5.4，否则放弃移动并进入步骤5.4；

步骤5.4，返回步骤5.1，并对步骤5.1、5.2、5.3进行迭代，直至没有可以移动的节点为止，进入步骤6；

步骤6，判断是否满足结束条件，结束条件为相邻节点之间的距离是否小于或等于预设阈值，不满足则在距离大于预设阈值的相邻节点间插入新的节点，通过步骤2所得网格号码链表找到所有与新的节点对应要改变的网格，根据步骤3重新获取所有与新的节点对应要改变的网格内的光谱灰度信息统计值T，根据步骤4重新计算所有网格的网格总能量J(I(x,y),θ)，直到满足结束条件，得到最终轮廓提取结果。

2.根据权利要求1所述遥感影像上土地利用地物边界提取方法，其特征是：当被处理的遥感影像为多通道影像，且通道的数目表示为M时，将步骤4的网格能量最小化函数变为多通道形式如下，

$J(I(x,y),\mathrm{\θ})=\underset{\mathrm{ar}=1}{\overset{R}{\mathrm{\Σ}}}{N}_{\mathrm{ar}}\left(\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\log (T_{2\mathrm{ar}}^m/N_{\mathrm{ar}}-{(T_{1\mathrm{ar}}^m/N_{\mathrm{ar}})}^2)\right).$

Images