CN107341837B - 基于影像金字塔的栅格-矢量数据转换及连续尺度表达方法 - Google Patents

Abstract

本发明公开了一种基于影像金字塔的栅格‑矢量数据转换及连续尺度表达方法，首先将遥感影像进行预处理，对预处理后的遥感影像进行多尺度可视化表达，获取多尺度栅格遥感影像；然后对步骤2中获得的多尺度栅格遥感影像进行噪声处理，对经噪声处理后的多尺度栅格遥感影像进行自动矢量化操作；最后对自动矢量化操作后的影像进行简化处理，获得简化处理后的矢量影像的连续尺度表达。本发明将目前成熟的遥感图像处理技术和地图制图综合有机融合，实现遥感影像数据的栅格‑矢量数据转换的自动化及遥感栅格数据和地图矢量数据的连续尺度表达，大幅度提高地理空间数据处理的效率和效果。

Description

基于影像金字塔的栅格-矢量数据转换及连续尺度表达方法

技术领域

本发明属于遥感图像处理和地图制图综合的技术领域，涉及一种地理空间数据栅格-矢量自动化转换及连续尺度表达方法，具体涉及一种融合遥感数据处理和地图制图综合的地理空间数据栅格-矢量自动化转换及连续尺度表达方法。

背景技术

地理空间数据的多尺度表达是摄影测量与遥感领域和地图学与地理信息系统领域中的重要内容之一，在遥感领域，多尺度表达技术不仅仅应用与影像数据可视化，同时也为影像特征提取、影像分类和混合像元分解等提供了切入点；在地理信息系统领域中，地理要素的多尺度表达依据现有制图规则和专题地物的分布规律进行空间数据的尺度间联合并应用于多尺度空间分析和相关决策。栅格-矢量数据转换是地理信息系统中数据转换的基本内容，并在数据采集、检索查询、数据融合以及数据共享方面提供灵活方便的数据存储形式。

多尺度表达通常采用影像金字塔技术。但是，一方面传统的多尺度表达在遥感领域和地图制图综合领域是基于不同的数据源分别进行处理，这种数据处理方式大幅度增加了地理空间数据处理的时间复杂度；另一方面，矢量数据的多尺度表达一般仅表达出在有限尺度上的结果，难以获取有限尺度之外的连续尺度表达。传统的栅格-矢量数据转换采用手工矢量化手段，但随着地理空间数据的海量增长，手工矢量化的人力成本和人为误差限制了矢量化效率和精度的提高。

发明内容

为了解决上述技术问题，本发明提供了一种融合遥感数据处理和地图制图综合的地理空间数据栅格-矢量自动化转换及多尺度连续表达方法。

本发明所采用的技术方案是：一种基于影像金字塔的栅格-矢量数据转换及连续尺度表达方法，其特征在于，包括以下步骤：

步骤1：将遥感影像进行预处理；

具体实现包括以下子步骤：

步骤1.1：在原始遥感影像的上选择云体阴影区、山体阴影区或在近红外波段上选择灰度值小于8的水体区域或红色波段上选择灰度值小于5的植被区域作为感兴趣区并计算感兴趣区像元的平均灰度值，用每个波段减去该平均值获得大气校正后的影像；

步骤1.2：根据遥感影像数据来源分析影像各波段的特点，结合实验影像中的地物分布特征，选择能够最明显反应出影像地物特征的三个波段进行影像的彩色融合。

步骤2：对预处理后的遥感影像进行多尺度可视化表达，获取多尺度栅格遥感影像；

具体实现包括以下子步骤：

步骤2.1：选用高斯核函数对遥感影像进行高斯卷积运算；假设预处理后的遥感影像I₀(x,y)的大小为m×n，其中x＝1,2,…,m；y＝1,2,…,n；高斯核函数为

卷积运算后影像为L₀(x,y,σ)：L₀(x,y,σ)＝G(x,y,σ)*I₀(x,y)；

其中，σ为尺度空间因子，是呈高斯分布的标准差，σ的值越大表示图像对应的尺度越大，图像越模糊；

步骤2.2：将卷积运算后影像L₀(x,y,σ)进行降采样处理以获得尺寸及分辨率均不同的多尺度栅格遥感影像数据；

多尺度栅格遥感影像的具体层数K取决于原始影像大小，其具体计算如下：K＝[log₂(min(m,n))]-a；

其中，m,n分别为原始影像的行数和列数，a为位于0～log₂(min(m,n))之间的任意整数。

步骤3：对步骤2中获得的多尺度栅格遥感影像进行噪声处理；

具体实现包括以下子步骤：

步骤3.1：将多尺度栅格遥感影像进行相似块分组，分组时首先在多尺度栅格遥感影像中选择并定义k×k大小的参照块和r×r大小的搜索区域，再根据相似块定义标准G(P)在影像中寻找与参照块差异度最小的影像块并整合成三维矩阵Q(P)；相似块定义标准G(P)具体如下，G(P)＝{Q:d(P,Q)≤τ^step1}；

其中，d(P,Q)表示两个影像块之间的欧式距离，τ^step1为距离阈值；

对分组后的三维矩阵进行三维变换并设定合适的阈值进行阈值收缩处理，将噪声大于阈值和不大于阈值的像元进行区分，并根据各影像块中低于阈值的像元数目及其噪声强度进行影像块权重的确定；对矩阵进行三维逆变换并将影像块恢复到原来位置得到带权重的影像，完成影像的基础估计；

步骤3.2：在基础估计后影像和噪声影像中分别选择并定义一定大小的影像参照块和搜索区域，分别在基础估计后影像和噪声影像中寻找与参照块差异度最小的影像块并形成三维矩阵，将三维矩阵进行三维变换；然后采用维纳滤波进行系数放缩处理，该系数取决于基础估计的三维矩阵及噪声强度；影像块内的像素值由维纳滤波系数和噪声强度作为权重进行加权更新获得，最后将影像块融合到原来的位置获得去噪后的影像I₀'(x,y)。

步骤4：对经噪声处理后的多尺度栅格遥感影像进行自动矢量化操作；

具体实现过程为：

将经噪声处理后的多尺度栅格遥感影像I₀'(x,y)分为待矢量化目标地物I_0S'(x,y)和影像背景I_0B'(x,y)两个部分，假设待矢量化地物的边界为C₀，C为演化追踪曲线，C将影像分为曲线内部C_in和曲线外部C_out两部分；构建能量泛函E(C)，将检测待矢量化地物边界问题转化为解能量泛函最值问题，采用欧拉-拉格朗日方程求解能量泛函获得矢量化影像。

步骤5：对自动矢量化操作后的影像进行简化处理；

具体实现过程为：

遍历自动矢量化操作后的影像，记录影像中的节点位置及对应的斑块并对斑块进行编号i，计算矢量化地物面积S_i；对不同尺度的矢量化影像上设定不同的面积阈值t，大于面积阈值的地物保留，其斑块属性值r_i设为1；小于面积阈值的地物舍去，其斑块属性r_i设为0：

步骤6：对简化处理后的矢量影像的连续尺度表达。

具体实现过程为：

首先选择在各不同比例尺图像中均存在的特征点，并记录特征点坐标；根据选择的特征点进行Delaunay三角剖分处理，获得大小两尺度上相对应的像素点坐标；读取两个不同尺度上相对应像元上的像素值，进行加权处理获得大小两尺度中间任意尺度图像以实现矢量图像的连续尺度表达。

本发明在现有遥感数据处理技术和地理信息系统地图综合方法的基础上，将地理信息数据的自动矢量化和栅格-矢量连续尺度表达结合在一起，实现地理信息数据的快速处理以及遥感技术与地理信息系统技术的结合。本发明主要采用高斯金字塔技术和块匹配与三维滤波算法进行遥感栅格数据多尺度表达和影像噪声控制，利用Chan-Vese模型实现栅格影像的自动矢量化，同时采用Morphing技术进行矢量化图像的连续尺度表达。

附图说明

图1是本发明实施例的原理图。

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

图3是本发明实施例的噪声处理流程示意图。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。

参见图1和图2，本发明实施例提供的一种基于影像金字塔的栅格-矢量数据自动转换及连续表达方法，包括以下步骤：

步骤1：将遥感影像进行预处理；

步骤1.1，在原始遥感影像的上选择云体阴影区、山体阴影区或在近红外波段上选择灰度值小于8的水体区域或红色波段上选择灰度值小于5的植被区域作为感兴趣区并计算其平均灰度值，并用每个波段减去这个平均值获得大气校正后影像；

步骤1.2，根据所用影像的特征，在ENVI平台上选用波段5(NIR)、波段6(SWIR1)和波段4(Red)进行RGB假彩色合成以区分陆地和水体。

步骤2：对预处理后遥感影像I₀(x,y)进行多尺度可视化表达，获取多尺度遥感栅格影像。

遥感影像多尺度可视化表达的方法主要有影像金字塔(含高斯金字塔和拉普拉斯金字塔)、小波变换等，本实施例中遥感影像的多尺度表达实现采用高斯金字塔方法，利用高斯核函数进行平滑处理获得平滑图像L₀(x,y,σ)。高斯核函数的具体形式及其卷积运算如下，

L₀(x,y,σ)＝G(x,y,σ)*I₀(x,y) (式2)

其中σ为尺度空间因子，是呈高斯分布的标准差，其取值为0.5.

将卷积运算后影像L₀(x,y,σ)进行降采样处理以获得尺寸及分辨率均不同的多尺度栅格影像数据。降采样处理一般舍去下层影像中的部分行和列，上一层影像由下层影像的部分像元数组成，多尺度栅格影像的具体层数K取决于原始影像大小；根据选取的Landsat-8影像分辨率及所选影像的尺寸特点，将多尺度层数K计算式如下，

K＝[log₂(min(m,n))]-a(式3)

其中a设为5，则获得多尺度层数K为4的多尺度栅格影像。

步骤3：对步骤2中获得的多尺度栅格遥感影像进行噪声控制，目前图像噪声处理的方法主要有基于空间域、变换域、偏微分、变分以及形态学这几大类的算法。本发明根据所选遥感影像特征采用块匹配与三维滤波技术。

如图3，该步骤包含以下子步骤：

步骤3.1，将多尺度栅格遥感影像进行相似块分组，分组时首先在多尺度栅格遥感影像中选择并定义k×k大小的参照块和r×r大小的搜索区域，再根据相似块定义标准G(P)在影像中寻找与参照块差异度最小的影像块并整合成三维矩阵Q(P)；相似块定义标准G(P)具体如下，

G(P)＝{Q:d(P,Q)≤τ^step1}(式4)

其中，d(P,Q)表示两个影像块之间的欧式距离，τ^step1为距离阈值。

对分组后的三维矩阵进行三维变换并设定合适的阈值进行阈值收缩处理，将噪声大于阈值和不大于阈值的像元进行区分，并根据各影像块中低于阈值的像元数目及其噪声强度进行影像块权重的确定。对矩阵进行三维逆变换并将影像块恢复到原来位置得到带权重的影像，完成影像的基础估计。

步骤3.2：与步骤3.1类似，在基础估计后影像和噪声影像中分别选择并定义一定大小的影像参照块和搜索区域，分别在基础估计后影像和噪声影像中寻找与参照块差异度最小的影像块并形成三维矩阵，将三维矩阵进行三维变换；然后采用维纳滤波进行系数放缩处理，该系数取决于基础估计的三维矩阵及噪声强度。影像块内的像素值由维纳滤波系数和噪声强度作为权重进行加权更新获得，最后将影像块融合到原来的位置获得去噪后的影像I₀'(x,y)。

步骤4：对经噪声处理后的多尺度栅格遥感影像实现自动矢量化操作。

对经噪声处理后的多尺度栅格遥感影像进行自动矢量化操作。自动矢量化操作尽管在一些软件上已经实现，如SuperMap,Raster 2Vector等。但本发明在多尺度的基础上主要采用Chan-Vese主动轮廓算法实现影像自动矢量化。将影像I₀'(x,y)分为待矢量化目标地物I_0S'(x,y)和影像背景I_0B'(x,y)两个部分，假设待矢量化地物的边界为C₀，C为演化追踪曲线,C将影像分为曲线内部C_in和曲线外部C_out两部分。构建能量泛函E(C)，将检测待矢量化地物边界问题转化为解能量泛函最值问题，采用欧拉-拉格朗日方程求解能量泛函获得矢量化影像。

步骤5，对自动矢量化图像进行制图综合处理。

遍历矢量化影像，记录影像中的节点位置及对应的斑块并对斑块进行编号i，计算矢量化地物面积S_i；对不同尺度的矢量化影像上设定不同的面积阈值t，大于面积阈值的地物保留,其斑块属性值r_i设为1；小于面积阈值的地物舍去，其斑块属性r_i设为0：

步骤6：矢量图像的连续尺度表达。

连续尺度表达技术主要有地图综合技术、Morphing模型以及渐进式传输技术。本发明主要采用Morphing模型实现栅格和矢量图像的大尺度和小尺度图像之间进行双线性插值处理，获得连续尺度表达。首先选择在各不同尺度图像中均存在的特征点，并记录特征点坐标；根据选择的特征点进行Delaunay三角剖分处理，获得大小两尺度上相对应的像素点坐标；读取两个不同尺度上相对应像元上的像素值，进行加权处理获得大小两尺度中间任意尺度图像以实现矢量图像的连续尺度表达。

本文中所描述的具体实施例仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施案例做各种修改或补充或采用类似的方式替代(例如：选用其它类型的金字塔或者分割算法等)，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (5)

1.一种基于影像金字塔的栅格-矢量数据转换及连续尺度表达方法，其特征在于，包括以下步骤：

步骤1：将遥感影像进行预处理；

步骤2：对预处理后的遥感影像进行多尺度可视化表达，获取多尺度栅格遥感影像；

步骤2的具体实现包括以下子步骤：

步骤2.1：选用高斯核函数对遥感影像进行高斯卷积运算；假设预处理后的遥感影像I₀(x,y)的大小为m×n，其中x＝1,2,…,m；y＝1,2,…,n；高斯核函数为

卷积运算后影像为L₀(x,y,σ)：

L₀(x,y,σ)＝G(x,y,σ)*I₀(x,y)；

其中，σ为尺度空间因子，是呈高斯分布的标准差，σ的值越大表示图像对应的尺度越大，图像越模糊；

步骤2.2：将卷积运算后影像L₀(x,y,σ)进行降采样处理以获得尺寸及分辨率均不同的多尺度栅格遥感影像数据；

多尺度栅格遥感影像的具体层数K取决于原始影像大小，其具体计算如下：K＝[log₂(min(m,n))]-a；

其中，m,n分别为原始影像的行数和列数，a为位于

之间的任意整数；

步骤3：对步骤2中获得的多尺度栅格遥感影像进行噪声处理；

步骤4：对经噪声处理后的多尺度栅格遥感影像进行自动矢量化操作；

步骤4的具体实现过程为：

将经噪声处理后的多尺度栅格遥感影像I₀'(x,y)分为待矢量化目标地物I_0S'(x,y)和影像背景I_0B'(x,y)两个部分，假设待矢量化地物的边界为C₀，C为演化追踪曲线，C将影像分为曲线内部C_in和曲线外部C_out两部分；构建能量泛函E(C)，将检测待矢量化地物边界问题转化为解能量泛函最值问题，采用欧拉-拉格朗日方程求解能量泛函获得矢量化影像；

步骤5：对自动矢量化操作后的影像进行简化处理；

步骤6：对简化处理后的矢量影像的连续尺度表达。

2.根据权利要求1所述的基于影像金字塔的栅格-矢量数据转换及连续尺度表达方法，其特征在于，步骤1的具体实现包括以下子步骤：

步骤1.1：在原始遥感影像的上选择云体阴影区、山体阴影区，或在近红外波段上选择灰度值小于8的水体区域，或在红色波段上选择灰度值小于5的植被区域作为感兴趣区并计算感兴趣区像元的平均灰度值，用每个波段减去该平均值获得大气校正后的影像；

步骤1.2：根据遥感影像数据来源分析影像各波段的特点，结合实验影像中的地物分布特征，选择能够最明显反应出影像地物特征的三个波段进行影像的彩色融合。

3.根据权利要求1所述的基于影像金字塔的栅格-矢量数据转换及连续尺度表达方法，其特征在于，步骤3的具体实现包括以下子步骤：

步骤3.1：将多尺度栅格遥感影像进行相似块分组，分组时首先在多尺度栅格遥感影像中选择并定义k×k大小的参照块和r×r大小的搜索区域，再根据相似块定义标准G(P)在影像中寻找与参照块差异度最小的影像块并整合成三维矩阵Q(P)；相似块定义标准G(P)具体如下，

G(P)＝{Q:d(P,Q)≤τ^step1}；

其中，d(P,Q)表示两个影像块之间的欧式距离，τ^step1为距离阈值；

对分组后的三维矩阵进行三维变换并设定合适的阈值进行阈值收缩处理，将噪声大于阈值和不大于阈值的像元进行区分，并根据各影像块中低于阈值的像元数目及其噪声强度进行影像块权重的确定；对矩阵进行三维逆变换并将影像块恢复到原来位置得到带权重的影像，完成影像的基础估计；

步骤3.2：在基础估计后影像和噪声影像中分别选择并定义一定大小的影像参照块和搜索区域，分别在基础估计后影像和噪声影像中寻找与参照块差异度最小的影像块并形成三维矩阵，将三维矩阵进行三维变换；然后采用维纳滤波进行系数放缩处理，该系数取决于基础估计的三维矩阵及噪声强度；影像块内的像素值由维纳滤波系数和噪声强度作为权重进行加权更新获得，最后将影像块融合到原来的位置获得去噪后的影像I₀'(x,y)。

4.根据权利要求1所述的基于影像金字塔的栅格-矢量数据转换及连续尺度表达方法，其特征在于，步骤5的具体实现过程为：

遍历自动矢量化操作后的影像，记录影像中的节点位置及对应的斑块并对斑块进行编号i，计算矢量化地物面积S_i；对不同尺度的矢量化影像上设定不同的面积阈值t，大于面积阈值的地物保留，其斑块属性值r_i设为1；小于面积阈值的地物舍去，其斑块属性r_i设为0：

5.根据权利要求1所述的基于影像金字塔的栅格-矢量数据转换及连续尺度表达方法，其特征在于，步骤6的具体实现过程为：

首先选择在各不同比例尺图像中均存在的特征点，并记录特征点坐标；根据选择的特征点进行Delaunay三角剖分处理，获得大小两尺度上相对应的像素点坐标；读取两个不同尺度上相对应像元上的像素值，进行加权处理获得大小两尺度中间任意尺度图像以实现矢量图像的连续尺度表达。

Images