大尺度自然地理实体边界智能化提取方法 —

自然地理实体(physical geographical objects,PGO)指地名指称的自然形成的地理上实际存在的自成单元的客观个体。在现实生活中, 一些自然地理实体具有明确的语义, 如范围很大的华北平原、四川盆地和辽东半岛等, 范围较大的大别山和秦岭等, 范围较小的香山和永定河谷等。然而, 这些实体在地图上却没有准确的分布范围或边界, 有些甚至在地图上没有明确的标注, 或仅以点对象的方式加以标注, 主要原因在于地图所能表达的实体均需由测绘技术获取其位置与空间分布范围。然而, 许多地理实体并不具有明确的能够用测绘技术测定的明确边界[1]。因此, 除具有明确分界线的地理实体和人工建筑物外, 大量的 PGO 并没有通过测绘或其他方式给出明确的边界。

在测绘领域, 一个对象必须经过人工测定其边界, 才能在地图上进行准确的表达。但是, 通过测绘工程或摄影测量在户外实测获得的地形图, 其比例尺一般在 1:10000 以下, 更小比例尺的地图一般由大比例尺地图综合而来。所以, 就传统测绘技术而言, 并没有测绘 PGO 分布范围和边界的具体方法。涉及一些小范围分布的对象(如河流、草地和林地等边界)时, 一般采用人工界定的方法。

地理边界的界定涉及人类对所处空间的认知,是认识地理现象和规律的重要手段[2]。目前, 农田边界[3]、气候带边界[4]和综合区域边界(如陆地热带北界[5])均有多种界定方法, 但自然地理实体边界的界定方法尚欠缺。直观而言, 利用其地理特征(如地形起伏、气候、海拔或植被分布等)是界定 PGO的分布范围或边界的最好方法[6]。然而, 由于历史文化因素对 PGO 的命名有很大的影响, 所以 PGO的范围界定不能单纯地考虑地理特征, 而是应该遵从其名称所指含义。PGO 的命名背景较为复杂, 与当地文化有着密切的联系。以北京周边景区为例,香山、妙峰山和野三坡等风景区的命名具有人文意义和历史意义, 与人们的生活联系紧密。除行政区划规定外, 景区自然形态的山或坡在地图上并没有明确的标注, 其范围可能只是一个小山头, 也可能是一片起伏连绵的山地, 或可能是一小块区域, 只有当地常年生活的民众才比较熟悉命名所指的大致范围。因此, 采用何种方法来界定自然实体的地理空间分布范围是一个比较复杂的问题。

语义明确但地理范围不明确的情况, 在地理信息领域通常被归纳为地理定义的不确定性[7]或模糊不确定性中的一种。空间不确定性研究是地理信息系统领域的重要研究内容[8], 通常采用误差估计方法给出一些不确定性的度量数值供用户参考[9–10]。由于对象边界存在不确定性(如由于季节性或年度性的水位变化, 河床宽度存在自然的变化; 沼泽,荒漠, 林地, 或一些自然植被的分界线并不明确),在数据采集时, 存在人为的取舍与界定问题。

严格地说, 确定地理实体边界应归结为地理信息本体定义与分类的问题。早年对地理本体问题探讨较多[11–12], 其中 PGO 相关内容涉及的问题主要是界定PGO语义内涵所确定的地理范围, 比如回答“香山在哪里?” “大别山北面有多大面积?”之类问题。这些问题均涉及相关地理实体的本体定义, 而研究地理本体的主要目标旨在给出语义描述的规范, 对于如何界定地理实体的分布范围, 目前并没有明确的可操作的方法。

在地理学领域, 对于山川湖泊类的 PGO, 一般均有常识性的定义, 比如对山的定义是“具有一定高度(海拔高和相对高), 并具有一定形态(山顶, 山坡, 山麓)的地貌结构, 称为山” [13]。这种定义仅仅可以帮助人们认知不同种类的 PGO, 而没有具体范围与边界的界定方法。自然地理实体类型多样化,且不同类型实体涉及的文化与地理特征均有很大的差异, 能够较为清晰地确定其边界的自然实体可分为山、平原、盆地 3 种。

对于旅游地区的名山大川, 当地或周边的企业或商业网点往往基于其文化特点来冠名, 越是靠近目标的区域, 此类名称越多。随着网络大数据的发展, 众源地理信息也日益增长, 我们可以通过导航地图获取此类 POI 信息。为此, 本文提出利用网络地图具有关联大数据的这一特点, 充分利用相关的语义信息, 再结合 PGO 的地理特征, 建立一种智能化的 PGO 分布范围确定方法, 解决地图中此类实体的明确表达问题。

本文提出一种智能化方法: 基于现有网络地图中的大数据信息, 获得与 PGO 的人文、历史方面关联紧密的 POI 数据, 通过将语义搜索与空间位置融合, 并结合 PGO 的空间几何特征, 运用相关算法提取其精确的分布范围与边界。本文仅以山地为例,阐述边界界定方法的具体内容。

1 实验数据的获取

1.1 地形数据的获取

简化起见, 本研究仅以山区为例探讨相关方法,使用的数据以数字高程模型(digital elevation models,DEM)为主。同时, 本研究以较大尺度范围的 PGO为研究重点, 故选用 90 m 分辨率的 SRTM (shuttle radar topography mission)数据, 相关数据通过 NASA(美国航空航天局)官网下载(对于较小尺度范围的实体, 可以选用 30 m, 甚至 5 m 分辨率的 DEM 数据,相关技术方法与下面所述基本上一致)。

1.2 语义信息的获取

对于百度和高德等导航地图, 对象位置的查询通常是其主要功能, 使用特定的关键词, 可以搜索到大量与之相关的内容, 我们称其为语义信息。如以关键词“大别山”为例, 在高德网络地图平台进行搜索, 可得到所有与“大别山”关联的语义信息, 搜索结果为上千条带有位置信息的地名, 包含所有大别山位置的标注信息, 同时也包含一些与大别山位置无任何关系, 而以大别山为冠名元素的餐馆、商店、产业园或其他类型的信息。一般而言, 离真正的大别山越近, 这类名称出现的概率越高。下面阐述如何基于这些语义信息, 确定 PGO 的概略范围。

2 PGO边界提取方法

如图1 所示, 首先获取 PGO 地名相关的互联网POI 信息, 接着基于搜索到的语义数据, 估计 PGO的概略范围, 然后判别概略范围内语义点分布的合理性, 并根据空间几何特征进行 PGO 分布范围的精确计算, 最后进行矢量转换, 提取其边界。

2.1 基于语义数据的PGO概略范围估计

前面搜索的点位, 无论表达的是餐馆还是产业园或公园, 甚至是商店, 当顾及其命名的文化因素时, 其位于 PGO 周边的概率远高于其位于一些无关位置的概率。同时, 考虑到我们所要确定的山区实体一般具有面状的覆盖区域, 一个理想的搜索结果应该是所有语义相关点(简称语义点)均落入同一个连通的区域, 而不是落入多个不连通的区域中。

虽然可以使用简单的聚类算法(如 K 均值法)对相对离散的语义点进行删减处理, 但聚类算法一般仅考虑空间距离, 且不同的阈值会产生不同的聚类结果, 从而较难判别一个语义点的有效性。鉴于搜索结果带有地理语义信息, 我们直接使用网络地图中已有的连通区域(如省、县、乡或镇的行政区域),与语义点集进行叠置分析, 选取涵盖最多语义点的一个连通区域作为对象空间分布的概略范围。判别过程如下。

1) 设搜索所得语义点集为 S, 点位数量为 n, 查询 S 中任一点 Pi (i∈n)所在县行政区域, 得到的结果为 Ai (i∈n), 同时建立映射关系 Ai→Pi。

2) 对 S 中的每一点执行此查询操作, 将所得结果 Ai (i∈n)加入同一个集合 T 中。

3) 对T中所有元素 Ai (i∈n)进行合并处理, 如果任两个元素 Ax (x∈n)和 Ay (y∈n)相等, 则保留其中一个, 同时将其映射点 Px 和 Py 加入点集 Sx (x∈m)中。

经过上述处理得到语义点的 m 个分组, 即所有语义点所在的 m 个行政区域。考虑到前面的一般性假设, 即原则上 PGO 空间范围应该是一个连续分布的面状区域, 故其范围应位于相邻的行政区域内,对步骤 3 所得结果继续进行步骤 4 的处理。

4) 判断 T 中所有处理后的元素 Ai (i∈m)的相邻关系, 对存在任意邻接关系的元素进行合并, 同时合并其对应的点集 Sx (x∈m), 最终得到以是否相邻为判断条件的聚类结果。图2 为由此过程得到的最终分组示例。顾及上述“涵盖最多语义点的连通区域为对象空间的概略范围”的约定, 所求区域即为max{Sx}(x∈m)。

2.2 PGO 概略范围内语义点集分布的合理性判别

虽然我们对语义点进行了分组, 初步将包含语义点数量最多的连通行政区域确定为地理实体所在的空间范围, 但此语义点集是否符合地理学认为的聚集模式, 还需要做进一步的判断。为了对此判断给出一个定量的描述, 我们引入最近邻距离分析法。该方法通常用于研究居民点布局、景观格局以及其他 POI 点的空间分布模式和特征。

设 n 为语义数据集中点的数量, di 为第 i 个语义点与最近邻点间的平面距离, A 为研究区的面积; 又设测量每个点与最近邻点之间的平均距离为 robs,其对应的随机分布的期望平均距离为 rexp, 两者的表达式分别为

则最近邻距离分析的一个常用统计量 R 的计算式为

判断此统计显著性的标准化值 Z 由下式计算得到:

其中, SE 为最近邻点平均距离的标准误差, 由下式计算得到:

R 统计量用于反映空间分布是否集聚, 如果 R值小于 1, 说明语义点的分布呈聚集模式; 如果 R 值大于 1, 说明语义点的分布为离散或竞争模式。因此, R 值越小, 语义点聚类程度越显著, 理想的 R 值应该小于 1。标准化值 Z 用来判断统计的显著性,当 Z 值大于 1.96 或小于−1.96 时, 说明该分布与随机模式之间具有显著差异, 且其绝对值越大, 点位分布模式的确定性越大, 与语义的相关性越强。因此, 可以利用 R 值和 Z 值评估搜索到的语义信息是否涵盖所搜索实体的地理位置和范围。

2.3 基于空间几何特征的PGO分布范围精确计算

山地在地形上有明显的起伏特征, 按 1:1 万比例尺地形图中地形类别的划分标准(GB/T13977—92《1:5000, 1:10000 地形图航空摄影测量外业规范》, 简称《规范》), 将地形分为 4 类: 平地、丘陵地、山地和高山地。规定山地应同时满足地面倾角介于 6°～25°之间以及高差介于 150～500 m 之间,高山地应同时满足地面倾角大于 25°以及高差大于500 m。本研究中, 我们使用坡度和高差两个因子进行山地范围的判别。

基于 DEM 数据, 计算 2.1 节所得概略区域内的地形坡度 SL 以及高差 ∆H; 基于高差 ∆H, 对区域内的地形进行分类, 假设所得分类指数为 T; 将 SL 与T 组合为 M, 作为判断是否为山地的因子:

其中, k 为当前窗口编号, N 为总的窗口数; 为方便可视化表达, 将 λ 和 μ 作为归一化参数, 当 Tk 取值范围为[1, 4]时, λ=16; 当 SLk 取值范围为[0, 90]时,μ=0.71。坡度 SLk 的计算公式为

式中, ∆H 表示单元格两边的高差, D 表示单元格的宽度。

基于式(5)的计算, 结合《规范》的判定条件,不难得到整个山地的精确分布范围。但是, 由于实际山地环境中偶尔也存在一些起伏度较小或高差较小的区域, 所以理论计算结果中难免存在一些不连续的区域。此外, 上述计算过程是基于 DEM 网格单元进行的, 这些不连续的区域在可视化的结果中会呈现为一些离散的像素集。

为了处理这种理论计算结果中的离散性与实际山地分布连续性不符的问题, 同时去除一些计算过程的噪声, 我们引入数学形态学方法, 对上述计算结果进行闭运算处理。

形态学闭运算包含一次膨胀运算和一次腐蚀运算, 假设 A 是待处理点集, B 是结构元素点集, 则使用 B 对 A 进行膨胀计算和腐蚀计算的公式如下:

其中, x 表示 A 中的元素。闭运算的优点在于能够弥补细小的不连续问题, 保持总形态和位置不变。因此, 对式(5)的计算结果进行闭运算操作, 能够在解决离散化问题的同时, 达到尽可能保持原始计算形态的目的。

2.4 PGO边界的提取

从基于 DEM 的计算到 PGO 范围的计算, 所处理的数据是栅格形式的, 但我们期望在地图上表达PGO, 故需要进一步提取表达其边界的矢量数据。

提取一个栅格对象的边界, 可使用二值化图像边界提取算法。此类算法较多, 如 Matlab 和 OpenCV函数包中均包含基于图像梯度的经典算子, 如 Canny 算子和高斯滤波算子等。

虽然使用形态学的闭运算能够对 PGO 分布的一些离散结果进行处理, 但对于一些较大的孔洞,闭运算处理的效果有限。主要原因在于, 为了保持结果的可靠性, 不能使用多次闭运算, 否则会造成原始计算结果的形态失真。为了处理较大形态的孔洞, 我们使用最大区域面积判断法, 过滤掉所有离散的孔洞区。此方法依赖的前提条件是, 假设内部具有最大面积的孔洞边界就是 PGO 的边界, 其余孔洞边界则仅仅包含其内部不连通区域或外部的离散区域。从经验的角度看, 这一假设对于大尺度 PGO而言是成立的。

3 实验与分析

考虑到名胜地区更受人们关注, 本研究选择大别山作为实验区, 从 NASA 官网下载的尽可能覆盖大别山区的 DEM 数据如图3 所示。

使用关键词“大别山”, 在高德地图中利用爬虫技术进行搜索, 对搜索结果中名称不包含搜索词的语义点进行过滤处理, 实际使用条目为 860 条, 表1列出其中前 10 条和后 10 条。

3.1 大别山概略范围判断

3.1.1 使用语义点所在行政区域估计大别山的概略范围

同时考虑地理与文化两方面的因素, 我们使用行政区域作为初步判断 PGO 所在范围的参考。使用 2.1 节所述方法, 将点数最大的行政区域选为主体区域, 去除主体区域以外离散的区划范围(图4(a)中粉红色区域), 得到大别山的概略范围(图4(a)中绿色区域)。我们对绿色区域与原始 DEM 数据进行叠加处理, 得到裁剪后反映大别山概略范围的 DEM数据(图4(b))。同时, 为了从视觉上比较该估计结果的正确性, 将语义点与 DEM 数据进行叠加。可以看出, 概略区域的语义点较为密集, 概略范围符合实际情况。

3.1.2 对所估计概略范围合理性的判断

对位于概略范围内的语义点进行聚集模式的判断, 使用 2.2 节所述判别方法进行计算, 所得结果如表2 所示。可以看出, Z 值远小于−1.96 的显著性判别标准, 说明语义点集具有很强的相关性; 同时R=0.366020, 说明概略范围内的语义分布呈显著聚集模式, 意味着这些语义点结合行政区划估计的概略范围与地理文化层面的大别山具有强关联性。

3.2 基于几何特征的大别山精确范围计算与边界提取

3.2.1 大别山精确范围计算

利用经过概略范围裁剪的 DEM 数据(图4(b)),通过式(5), 基于坡度和高差因子计算 M 值。根据《规范》中的地形分类, 计算得到判断是否为山地的 M 阈值为 36.26。使用此阈值对初步计算的结果进行处理, 得到大别山的精确分布范围(图5(a))。

将图5(a)的局部放大可以看出, 由于噪声以及微小平缓地貌的存在, 基于 DEM 栅格数据的计算结果中存在大量孔洞和离散区域。在进一步清晰地界定大别山的 PGO 范围之前, 使用形态学闭运算对原始计算结果进行处理。为了尽可能保持原始边界的形态, 我们将式(7)中的 B 元素设置为 6×6 的矩形结构, 并对图5(a)数据进行计算, 得到图5(b)所示结果。

3.2.2 大别山矢量边界提取

经过形态学算法处理后, 图5(b)中仍然存在很多不连续的区域和孔洞, 这是由于实际山区内部存在大量较为平缓的区域。从图4(b)中语义点与 DEM叠加的效果可以看出, 山区存在大量居民生活区。在宏观地理学意义上, 这些居民生活区属于大别山区, 应该将其划分在整个大别山的 PGO 之中。但是, 使用形态学算法会使图5(a)中 PGO 范围失真。对这些区域, 我们使用二值化边界提取方法, 得到尽可能接近真实形态的大别山边界。

我们对图5(a)所示结果使用 OpenCV 库的 Find-Contours 函数进行处理, 通过扫描整张二值图像,找到所有的轮廓, 并将其存储在一个轮廓的向量中,提取包括内部孔洞的所有边界, 结果如图6(a)所示。FindContours 函数提取的边界实际上是每个封闭区域的边界, 符合我们对山地实体表达的需求,故在得到所有区域的边界后, 使用 ContourArea 函数计算这些区域的面积, 选取面积最大的区域作为大别山 PGO 的边界和范围(图6(b))。

3.3 结果与分析

将图6(b)叠加到 Google 地图上, 得到图7(a)所示结果。如果对该边界采用多边形填充来表达, 则可以展示基于电子地图的 PGO 查询效果(图7(b))。图7(c)展示基于现有高德地图查找“大别山”的结果,仅能得到一些零散的标注点, 虽然用户可以基于遥感影像目视识别山地范围, 但高德地图缺乏大别山地区边界的数据, 因而无法对其做出完整的表达。

为了验证本文方法所得结果的正确性, 我们采用人工解译方法, 基于 Google 地图的地形数据, 提取大别山区的轮廓线(图7(d))。分别使用余弦相似度、平均结构相似性(structural similarity, SSIM)和汉明距离来度量自动提取边界与人工解译边界之间的相似性, 对仅包含实体边界的两张图片进行计算,结果如表3 所示。

余弦相似度和平均结构相似性度量指数越接近1, 说明图片越相似; 汉明距离为 0, 说明图片完全一致; 介于 0～5 之间, 说明图片有所不同, 但基本上相似; 汉明距离大于 10, 说明图片完全不相似。综合 3 种相似性判别因子, 能够说明基于本文提出的智能化分析和计算方法, 可以得到比较符合实际的 PGO 范围和边界, 为其在数字化地形图上的表达提供数据, 也为进一步基于电子地图分析这些实体提供支撑。

3.3.2 分析与讨论

1) 结果的精度。由于 PGO 的边界具有不确定性, 故本文方法提取的边界是 PGO 的参考边界。就总体效果而言, 使用本文算法得到的 PGO 边界符合日常地理认知。考虑到大别山是宏观自然地理实体,其表达应聚焦于小比例尺地形图, 同时考虑到数据处理效果, 本文选用 90 m 格网的 DEM 数据。但是,现实生活中对山地的认定一般遵从中等比例尺, 故选用 1:1 万地图对地形的划分标准。粗略地估计,本文方法提取的边界误差在 100 m 左右, 能够达到1:50 万的成图要求。需要说明的是, 自然地理实体范围的界定受 POI 数量及其分布影响, 尤其是人文信息方面 POI 语义的模糊性影响, 我们无法给出一个所有人都能接受的, 被认为是准确无误的大别山范围, 此处所言精度仅仅指所提取边界在图上的表达精度。

2) 语义信息完整性的影响。本文方法适用于相关 POI 名称信息较为丰富的地理实体, 尤其是具有一定历史的与文化相关的实体。本文选用大别山为研究对象, 由于其具有相当大的历史文化影响力,故搜索得到的网络语义信息比较充分, 能满足本文处理流程的要求。对于一些文化影响力较弱的PGO, 在网络上很难搜索到足够的语义信息, 除借助地理学与地质学的知识外, 还应将地表影像和气候等更多的地理信息纳入提取流程中, 增强提取结果的可靠性。此外, 鉴于目前遥感与地理信息获取的便捷性, 此处仅强调语义信息完整性的影响, 事实上, 一个令人信服的、占据宏观尺度范围的自然地理实体, 除受 POI 点的人文信息影响外, POI 的空间分布对识别其分布范围也非常关键, 因此除了语义信息, 地理空间信息也必不可少。

3) 不同类型 PGO 提取结果的差异。本文方法对不同类型 PGO 的提取均有效, 概略范围的计算方法完全一致。在计算精确范围时, 由于不同类型PGO 的地表形态完全不同, 故应采用不同的约束规则。例如, 平原和高原地区可以采用起伏度和海拔高度, 盆地可以采用起伏度及其变化曲线等作为主要计算因子。对于网络地图缺乏足够语义约束信息的情况, 还应纳入地理与地质学的判断准则, 如天山与昆仑山分界处采用土壤颜色加以区分, 大兴安岭与小兴安岭以嫩江为界进行划分, 等等。

4 结论

本文针对现有地图上较大尺度 PGO 的表达问题, 提出基于语义大数据的, 结合自然行政区划以及 PGO 的空间几何特征等信息, 自动提取符合认知的 PGO 范围或边界的矢量数据, 解决了此类实体在地图上作为一个面对象的表达问题。以大别山为例,有效地验证了所提方法的可行性。

我们目前仅探讨了具有强烈文化色彩且相对孤立的山地实体, 下一步将重点探讨诸如秦岭、大巴山之类更复杂的山地以及需要考虑更多地理信息和人文信息的 PGO 智能化提取方法。

[1] Burrough P A, Frank A U.Geographic objects with indeterminate boundaries.Population & Development Review, 1996, 23(2): 437–439

[2] Tang G, Li J, Xiong L, et al.Scientific attributes and expression methods of geographical boundary.Journal of Geographical Sciences, 2022, 32(6): 1119–1135

[3] 王劲峰, 徐成东.地理探测器: 原理与展望.地理学报, 2017, 72(1): 116–134

[4] 晏路明.基于GIS的自然地理界线的判别分析//2003中国现场统计研究会第十一届学术年会论文集(上).北京: 中国现场统计研究会, 2003: 292–297

[5] 董玉祥, 徐茜, 杨忍, 等.基于地理探测器的中国陆地热带北界探讨.地理学报, 2017, 72(1): 135–147

[6] 张伟, 李爱农, 江晓波.基于DEM的中国山地空间范围定量界定.地理与地理信息科学, 2013, 29(5):58–63

[7] Fusco G, Caglioni M, Emsellem K, et al.Questions of uncertainty in geography.Environment and Planning A, 2017, 49(10): 2261–2280

[8] 胡圣武, 余旭.空间数据不确定性研究进展.河南理工大学学报(自然科学版), 2016, 35(6): 815–822

[9] 陈鹏飞.空间大数据中的不确定性分析与数据质量控制[D].武汉: 武汉大学, 2019

[10] 胡圣武.GIS质量评价与可靠性分析.北京: 测绘出版社, 2006

[11] 孙敏, 陈秀万, 张飞舟.地理信息本体论.地理与地理信息科学, 2004, 20(3): 6–11

[12] 陈建军, 周成虎, 王敬贵.地理本体的研究进展与分析.地学前缘, 2006, 13(3): 81–90

[13] 王明业, 朱国金, 贺振东, 等.中国的山地与山地系统.山地研究, 1986, 4(1): 67–74

大尺度自然地理实体边界智能化提取方法——以大别山为例