基于Ganos GeoSOT地理网格模型：无人机路径规划能力实践_云原生数据库 PolarDB(PolarDB)-阿里云帮助中心

随着新能源技术的迅速发展，低空经济逐渐成为新的战略性新兴产业。然而，与传统的地面活动不同，低空活动具有立体性、区域性和融合性等特点，这些特点为安全引导低空活动的顺利开展带来了系列亟待解决的技术性问题。Ganos GeoSOT地理网格引擎提供了基于网格的路径规划能力，能够利用数字高程模型（DEM）、数字表面模型（DSM）、倾斜摄影等数据构建复杂环境中的无人机路径规划应用。

关于Ganos地理网格模型

地理网格是一种用于再现地球表面的多边形网格单元集合，能够有效表示地物在地理空间中的位置信息，并融合其他各类时空数据。地理网格计算通常采用由粗到细的逐级分割方式，对地球表面进行处理。通过将地球的曲面用一定大小的多边形网格进行近似模拟，实现地理空间定位与地理特征描述的有机结合，同时将误差范围控制在网格单元的可接受范围内。每个网格单元都会进行编码，网格与编码是一一对应的。三维地理网格不只考虑经纬度，还把高度维纳入剖分和编码范围。

Ganos地理网格引擎目前涵盖GeoSOT和H3两种地理网格。

GeoSOT是中国提出的一套地球空间剖分理论，并在此基础上发展出的一种离散化、多尺度区域位置标识体系。
H3是Uber研发的一种覆盖全球表面的二维地理网格，采用的基本网格是正六边形。关于H3网格的最佳实践可参考H3地理网格能力解析与最佳实践。

Ganos支持退化网格计算（如下图），充分利用网格的层级关系，通过更精简的网格组合对空间范围进行有效表达。此外，Ganos自研的地理网格索引，可用于高效查询网格编码并加速聚合计算。

业务背景

随着新能源技术的迅猛发展，低空经济已经逐步成为新的战略性新兴产业。低空经济是一种综合性经济形态，依托各种有人驾驶和无人驾驶航空器的低空飞行活动，推动相关领域的融合发展。其应用场景丰富，不仅涵盖了传统通用航空业态，还融合了以无人机为支撑的低空生产服务方式，广泛应用于工业、农业和服务业等多个领域。低空经济在构建现代产业体系方面具有重要作用，且发展空间极为广阔。

不同于传统的地表活动，低空活动具备立体性、区域性和融合性等特征。这些特征为安全引导低空活动的顺利开展带来了一系列亟待解决的技术问题。

立体性：主要体现为“活动空间上的立体性”。低空空域是指其垂直范围原则为真实高度1000米以下。根据不同地区的特点和实际需求，具体划定的高度范围构成了空域的活动空间。这一空域的划分使经济活动得以从地面向空中延伸，实现从“平面经济”向“立体经济”的转变，体现了一种三维空间的立体经济形态。
区域性：主要体现为“作用范围上的区域性”。与航空运输经济、高铁经济等大规模、大范围、一体化特征的经济形式不同，低空经济基于小飞机、小航线和小企业，具有小规模、小范围及分散性特点，且显著体现出地域性和区域性特征。
融合性：主要体现为“运行模式上的融合性”。与航空运输及铁路、公路、水运等交通运输经济相比，低空经济的主要作用不仅限于交通运输，更侧重于行业服务。其体现为“行业+航空”的模式，为相关行业提供空中解决方案或辅助手段，以提升工作效率、降低成本、增强获得感等。

通过分析低空活动的特性可以发现，低空区域并不存在路网。在三维空间中，如何确保低空设备的安全通行与到达，以及如何合理划分低空区域的活跃周期与飞行器的活动规划，都是低空经济顺利实施并保障区域活跃性与安全的重要工作。

功能介绍

Ganos地理网格引擎提供了基于网格的路径规划能力，能够利用数字高程模型（DEM）、数字表面模型（DSM）、倾斜摄影等数据构建复杂环境中的无人机路径规划应用场景。相关函数简介如下：

ST_SetCost
用于给障碍物设置通行开销。障碍物可以是实体，例如地面、建筑物，也可以是空间范围，如雷达扫描范围或风场空间。在使用该函数之前，必须首先将障碍物的空间范围转换为网格数组，然后进行通行开销的设置。
ST_3DGridPath
无人机路径规划的最核心函数，根据起始点、终止点以及通行的开销数据，计算出最优的三维网格路径。支持不同的路径算法，如A*和Dijkstra，不同的移动方式比如是否允许斜向飞行，以及不同的距离计算方法比如曼哈顿和欧几里得。
ST_MatchGridLevel
根据地形高程数据的分辨率确定最高的网格剖分层级，路径计算时网格计算的层级应低于或等于此层级。
ST_CostUnion
合并所有的障碍物网格开销数组，支持退化网格合并，重叠的网格采用最大的开销并只保留一个，作为最终规划路线的输入。

技术优势

基于Ganos GeoSOT地理网格引擎开展无人机路径规划能力的支撑，具备以下技术特点：

统一的标准数据模型，所有输入数据均在Ganos内部采用标准数据模型进行表达（包括栅格、表面网格等）。
所有相关的数据处理均可通过Ganos内部函数完成，同时还提供了完善的路径规划辅助函数。
无人机路径规划可以采用退化网格技术，以降低存储开销并提升计算速度。
路径规划提供多种参数选项，包括规划算法、移动方式和距离计算方法，以满足不同的需求。

最佳实践

本部分以某园区的倾斜摄影数据为例，演示三维路径规划在实际场景中的应用。

建议配置

为了得到良好的体验，PolarDB集群建议使用以下配置：

数据库引擎：PostgreSQL 14
内核版本：14.12.23.1
CPU：>= 4 核
内存：>= 16 GB
磁盘：>= 100 GB
Ganos版本：>=6.7

说明

在导入及计算OSGB数据的过程中，需要占用一定量的内存，建议内存容量大于16 GB。

测试数据

数据入库与处理

安装网格路径规划功能所需插件。

CREATE EXTENSION IF NOT EXISTS ganos_geomgrid CASCADE;
CREATE EXTENSION IF NOT EXISTS ganos_utility CASCADE;

测试数据存储于OSS中，使用ST_ImportOSGB函数将其导入，详细参数介绍可参考ST_ImportOSGB。

SELECT 
  ST_ImportOSGB(
    'test_oblique',
    'OSS://<access_id>:<secrect_key>@<Endpoint>/<bucket>/path_to/file', 
    '{"parallel":16,"project":"building"}');

将某一级别的瓦块整体转换为单一GLB格式，并存入临时表temp_glb中，采用SfMesh类型作为中间类型。
```
SELECT 
  ST_ImportGLB(
    'temp_glb', 
    ST_AsGLB(ST_Combine(tile)), 
    'temp_glb_1', 
    '{"ignore_texture":true,"ignore_normal":true}'
  ) 
FROM test_oblique_tile WHERE lod = 19 AND project_name = 'building';
```
说明
- 调用ST_ImportGLB函数时，忽略无用的纹理（ ignore_texture ）和法线（ ignore_normal ）。
- 此处选用lod = 19的瓦块作为采样瓦块。您可以自行选择适当级别的瓦块。请注意，更高LOD级别的瓦块虽然会显著延长采样时间，但其精度也相应提高。
- 测试倾斜数据为常见的倾斜数据命名，例如Tile_+006_+004_L14_0.osgb，可以从文件名中分析出其LOD级别。
- 对于命名没有规律的倾斜数据，因无法获取其LOD值，可以将WHERE lod = 19替换为WHERE children IS NULL，即选用最高精度的叶节点瓦块作为采样瓦块。
更新临时数据。将其投影至EPSG:4490坐标系下：
```
UPDATE temp_glb 
SET 
  gltf_data = ST_Transform(
    ST_Flatten(
      ST_SetSrid(
        ST_Translate(
          ST_YupToZup(gltf_data), 
          x_off, y_off,z_off), 
        srid), 
      FALSE), 
    4490) 
FROM 
  (SELECT 
      srid, 
      ST_X(anchor) x_off, 
      ST_Y(anchor) y_off, 
      ST_Z(anchor) z_off 
    FROM test_oblique WHERE project_name = 'building') metadata 
WHERE gltf_id = 'temp_glb_1';
```
上述语句执行了如下操作：
- OSGB入库后转换为GLB时是Y轴向上。为正常处理其坐标，需要调用ST_YupToZup函数将其转换为Z轴向上。
- 倾斜数据往往以相对坐标记录点，因此需要调用ST_Translate函数将每个点都以锚点（anchor）做偏移，使每个点的坐标都为绝对坐标，并通过ST_SetSrid函数设置正确的SRID。
  说明
  该操作要求原始倾斜数据具有完整的元数据描述文件（metadata.xml），否则无法正确地获取锚点和坐标系。
- 对于包含旋转矩阵的SfMesh对象，需预先进行压平处理（ST_Flatten），以将旋转矩阵应用于实际坐标，便于进行坐标投影。
- 由于网格自身的坐标系规定为CGC 2000，因此需要利用ST_Transform函数将每个点投影至EPSG:4490坐标系。

将中间数据转换为几何网格（GeomGrid）。

-- 创建网格表
CREATE TABLE IF NOT EXISTS building_grid
(
   id   SERIAL,
   grid GEOMGRID,
   grid_type TEXT
);
-- 生成GeomGrid并插入网格表
INSERT INTO building_grid (grid, grid_type)
SELECT grid, 'border'
FROM
  (SELECT unnest(ST_As3DGrid(gltf_data, 24, TRUE)) grid
   FROM temp_glb WHERE gltf_id = 'temp_glb_1' ) tmp;
-- 移除临时表
DROP TABLE temp_glb;

其中，ST_As3DGrid函数能够根据实际需求指定地理网格的精度级别。精度越高，所得到的网格越为细致，然而，采样所需的时间也将显著增加。不同级别下每个网格的大致尺度可以参考下表：

级别	大致尺度	级别	大致尺度	级别	大致尺度
0	全球	11	29.6 km	22	15.5 m
1	1/4 地球	12	14.8 km	23	7.7 m
2	-	13	7.4 km	24	3.9 m
3	-	14	3.7 km	25	1.9 m
4	-	15	1.8 km	26	1.0 m
5	-	16	989.5 m	27	0.5 m
6	890.5 km	17	494.7 m	28	24.2 cm
7	445.3 km	18	247.4 m	29	12.0 cm
8	222.6 km	19	123.7 m	30	6.0 cm
9	111.3 km	20	61.8 m	31	3.0 cm
10	59.2 km	21	30.9 m	32	1.5 cm

网格展示

生成的网格可以导出为GLB格式，以便于浏览大致形状。
```
SELECT 
  ST_AsGLB(
    ST_ZupToYup(
      ST_3DRemoveDuplicateVertex(
        ST_MergeGeomByMaterial(
          ST_Triangulate(
            ST_Translate(
              ST_Transform(
                ST_Flatten(
                  ST_Collect(
                    ST_AsMeshgeom(array (SELECT grid FROM building_grid WHERE grid_type = 'border')):: sfmesh[]),
                  TRUE
                ), 
                srid
              ), 
              - ST_X(anchor), - ST_Y(anchor), - ST_Z(anchor)
            )
          )
        ), 
        0.1
      )
    )
  ) 
FROM test_oblique WHERE project_name = 'building';
```
上述语句执行了如下操作：
- 通过ST_AsMeshGeom函数将网格数组转换为MeshGeom数组。
- 将MeshGeom数组类型转换为SFMesh数组，并聚合（ST_Collect）成独立的SfMesh对象。
- 将聚合后的对象压平（ST_Flatten），并重新投影（ST_Transform）回原始坐标系。
- 将聚合后的对象整体偏移（ST_Translate），使全部坐标点都变回相对于原点的相对坐标。
- 将聚合后的对象三角化（ST_Triangulate），方便接下来的处理步骤。
- 联合调用按材质合并几何函数（ST_MergeGeomByMaterial）与按阈值移除点函数（ST_3DRemoveDuplicateVertex）可以减少数据量，降低渲染压力。
- 调用ST_YupToZup函数将左手坐标系下的数据转换到GLTF所使用的右手坐标系，最终调用ST_AsGLB函数导出GLB。
不同级别的效果对比可参见下表（以25级和24级为例）：
级别
效果展示
原始数据
25级
24级

网格对象也可以通过以下语句生成适合Cesium读取的数据：

-- 创建函数方便后期调用
CREATE OR REPLACE FUNCTION GridToJson(geomgrid[])
RETURNS json AS $$
SELECT to_json(
        array_agg(
            json_build_object('center', array_to_json(array[ST_X(center), ST_Y(center), ST_Z(center)]), 
                           'size', array_to_json(array[ST_DistanceSpheroid(min_point, ST_Point(ST_X(max_point), ST_Y(min_point), 4490)), ST_DistanceSpheroid(min_point, ST_Point(ST_X(min_point), ST_Y(max_point), 4490)), z]))))
FROM
  (SELECT ST_Transform(ST_PointZ(ST_X(st_centroid(geom)), ST_Y(st_centroid(geom)), (ST_ZMax(geom) + ST_ZMin(geom)) / 2, 4490), 4479) center,
          ST_SetSrid(BOX[0]::geometry, 4490) min_point,
          ST_SetSrid(BOX[1]::geometry, 4490) max_point,
          ST_ZMax(geom) - ST_ZMin(geom) z
   FROM
     (SELECT ST_ASBox(grid) BOX, ST_AsGeometry(grid) geom
      FROM (SELECT unnest($1) grid) a)b)c $$
LANGUAGE 'sql';
-- 调用GridToJson函数生成结果
SELECT GridToJson(array (SELECT grid
                     FROM building_grid WHERE grid_type = 'border')) ;

返回结果为JSON格式的数组，如下：

[
  {
    "center": [-1583397.2647956165,5313841.088458569,3142388.7651142543],
    "size": [3.3600142470488144,3.848872238575258,3.3680849811062217]
  },
  ...
]

其中：

center：每一个网格中心点投影至EPSG:4479坐标系下的坐标。
size：每个网格在纬度方向/经度方向/高度方向的大致尺寸。

在Cesium中，可以通过如下方式加载上述数据：

<!DOCTYPE html>
<html lang="en">
<head>
  <meta charset="utf-8">
  <script src="https://cesium.com/downloads/cesiumjs/releases/1.89/Build/Cesium/Cesium.js"></script>
  <link href="https://cesium.com/downloads/cesiumjs/releases/1.89/Build/Cesium/Widgets/widgets.css" rel="stylesheet">
  <style>
    html,
    body,
    #cesium_container {
      width: 100%;
      height: 100%;
      margin: 0;
      padding: 0;
      overflow: hidden;
      }
  </style>
</head>
<body>
  <div id="cesium_container"></div>
  <script>
    Cesium.Ion.defaultAccessToken = YOUR_CESIUM_TOKEN
    const viewer = new Cesium.Viewer('cesium_container', ['timeline', 'animation', 'infoBox', 'navigationHelpButton'].reduce((opts, opt) => (opts[opt] = false) || opts, {}))
    
    // 指定需要显示的网格
    const grids = [{ "center": [-1583315.44750805, 5313895.98214751, 3142303.525767127], "size":... 
    const add_grids = (grids, material) =>
      grids.forEach(({ center, size }) => viewer.entities.add({
          position: new Cesium.Cartesian3(...center),
          box: {dimensions: new Cesium.Cartesian3(...size),
            material, outline: true, outlineColor: Cesium.Color.BLACK}}))
      
    // 使用半透明的红色展示
	add_grids(grids, Cesium.Color.RED.withAlpha(0.2))
      
    // 跳转到网格位置
    viewer.zoomTo(viewer.entities)
  </script>
</body>
</html>

将24级网格与倾斜数据叠加显示，叠加倾斜数据代码请参考附录效果如下：

网格路径规划

构建网格后，可以在上述场景中以某两点作为起止点，使用ST_3DGridPath函数规划基于网格的三维路径：

-- 定义起点/终点/可搜索范围
WITH vals (start_p, end_p, box_range) AS (
  VALUES 
    (
      ST_PointZ(106.59182291666669, 29.70644097140958, 355.322855208642, 4490), 
      ST_PointZ(106.59244791666666, 29.707135415854022, 348.58669235736755, 4490), 
      'BOX3D(106.591788194444 29.7062326380763 341.85053662386713, 106.592899305556 29.7071701380762 362.0590251699541)' :: box3d
    )), 
border AS (
  SELECT grid, grid_type 
  FROM building_grid, vals 
  WHERE 
    -- 仅计算范围内的网格
    -- 所有参与计算的几何类型均要求空间参考为EPSG:4490
    ST_3DIntersects(ST_SetSrid(box_range :: meshgeom, 4490), grid) 
    OR ST_3DContains(ST_SetSrid(box_range :: meshgeom, 4490), grid)
) 
SELECT ST_3DGridPath(start_p, end_p, box_range, array[ST_SetCost( array (SELECT grid FROM border WHERE grid_type = 'border' ), -1)], '{"algorithm":"astar","movement":"cross","distance":"euclidean"}') result 
FROM vals;

其中：

ST_SetCost函数为网格设定同行成本。

不可通行区固定为-1。
若无其他设置，任意非不可通行区的通行成本固定为1。

使用上文创建的GridToJson函数，可以在Cesium中浏览效果：

tiu6hs7lgzu44_9abd608ffe3b40a0b61e50c8eb68f3ad

tiu6hs7lgzu44_2bd2c03c8d3d414e9a4d465d96f33e83 tiu6hs7lgzu44_d496f4033f4748e1877bbf968e656824

若将移动方式修改为strict_octothorpe，则其效果为：

tiu6hs7lgzu44_e4c154f88f4c4018969a2492b97f49f2

tiu6hs7lgzu44_64f6ef2d064d4e11ba56915ccbf0f1d7 tiu6hs7lgzu44_a3d48bf8de084ffb9c5b0d7a1d82bfe3

多成本区网格路径规划

在可通行区域为不同区域设置不同的通行成本后，算法将选择综合成本最小的路径。我们在前述场景中添加了一个高成本通行区域，以示范多成本网格路径规划的应用。

在building_grid表中插入代表高成本通行区的网格：

INSERT INTO building_grid (grid, grid_type)
SELECT grid, 'risk_zone'
FROM
  (SELECT unnest(ST_As3DGrid(ST_SetSrid('BOX3D(106.59192708333332 29.706614582520686 341.8505366235368, 106.59220486111111 29.706892360298472 365.4271128197744)' :: box3d :: meshgeom, 4490), 24)) grid) tmp;

该区域形如（橘色部分）：

说明

与非通行网格重合的高成本通行区网格也将视为非通行网格。

执行路径规划。

WITH vals (start_p, end_p, box_range) AS (
  VALUES 
    (
      ST_PointZ(106.59182291666669, 29.70644097140958, 355.322855208642, 4490), 
      ST_PointZ(106.59244791666666, 29.707135415854022, 348.58669235736755, 4490), 
      'BOX3D(106.591788194444 29.7062326380763 341.85053662386713, 106.592899305556 29.7071701380762 362.0590251699541)' :: box3d
    )), 
border AS (
  SELECT grid, grid_type 
  FROM building_grid, vals 
  WHERE 
    ST_3DIntersects(ST_SetSrid(box_range :: meshgeom, 4490), grid) 
    OR ST_3DContains(ST_SetSrid(box_range :: meshgeom, 4490), grid)
) 
-- 设定 risk_zone 的通行成本为2
SELECT ST_3DGridPath(start_p, end_p, box_range, array[ST_SetCost( array (SELECT grid FROM border WHERE grid_type = 'border' ), -1), ST_SetCost( array (SELECT grid FROM border WHERE grid_type = 'risk_zone' ), 2)], '{"algorithm":"astar","movement":"cross","distance":"euclidean"}') result 
FROM vals;

此时效果为：

对比没有高成本通行区域的路径（紫色路径）：

总结

利用Ganos地理网格引擎的相关方法，通过少量代码实现无人机路径规划。Ganos作为全球首个专业级空间数据库，已经将狭义的空间数据拓展至“空天地、室内外、地上下、动静态”等全空间范畴，从数据库系统最底层为物理世界数字化提供时空处理框架，未来Ganos还将提供更多高效的库内空间分析与可视化能力，推动各行业的空间信息应用真正走向“视算一体”。

附录

使用Node.js实现简易的后端服务，代理3DTiles请求，并以Cesium作为前端框架展示结果。

在自定义目录下编写依赖文件package.json：

{
  "dependencies": {
    "koa": "^2.14.2",
    "koa-router": "^12.0.0",
    "koa-send": "^5.0.1",
    "pg": "^8.10.0"
  }
}

在该目录执行npm install安装依赖库。在同一目录下创建Node.js脚本文件index.js：

const Koa = require('koa');
const Send = require('koa-send');
const Router = require('koa-router');
const router = new Router()
const { Pool } = require('pg');
const pool = new Pool({ user: <YOUR_USER>, host: <YOUR_HOST>, database: <YOUR_DB_NAME>, port: <YOUR_PORT> });

const TABLE_NAME = 'test_osgb'
const PROJECT_NAME = 'prj1'

router.get('/', (_) => Send(_, '/index.html'))

/* get metadata of current project */
router.get('/project', async (ctx) => {
  const sql = `SELECT (REGEXP_MATCHES(ST_ASTEXT(ST_TRANSFORM(ANCHOR,4326)), 'POINT Z \\((.*?)\\)'))[1] _ANCHOR, PROJECT_ID FROM ${TABLE_NAME} WHERE PROJECT_NAME='${PROJECT_NAME}';`
  const { rows: [{ _anchor, project_id }] } = await pool.query(sql)
  const anchor = _anchor.split(' ').map(x => parseFloat(x))
  const url = `/tileset/${project_id}/${project_id}`
  ctx.body = { url, anchor }
})

router.get('/tileset/:project_id/:uid', async (ctx) => {
  const { params: { project_id, uid } } = ctx
  const sql = `SELECT TILESET FROM ${TABLE_NAME}_TILESET WHERE PROJECT_ID=$1::UUID AND UID=$2::UUID;`
  const { rows: [{ tileset }] } = await pool.query(sql, [project_id, uid])
  ctx.body = tileset
})

router.get('/b3dm/:project_id/:uid', async (ctx) => {
  const { params: { project_id, uid } } = ctx
  const sql = `SELECT ST_ASB3DM(TILE) TILE FROM ${TABLE_NAME}_TILE WHERE PROJECT_ID=$1::UUID AND UID=$2::UUID;`
  const { rows: [{ tile }] } = await pool.query(sql, [project_id, uid])
  ctx.body = tile
})

new Koa().use(router.routes()).listen(5500, '0.0.0.0');

在同一目录下创建HTML文件index.html。

<!DOCTYPE html>
<html lang="en">

<head>
  <meta charset="utf-8">
  <script src="https://cesium.com/downloads/cesiumjs/releases/1.89/Build/Cesium/Cesium.js"></script>
  <link href="https://cesium.com/downloads/cesiumjs/releases/1.89/Build/Cesium/Widgets/widgets.css" rel="stylesheet">
  <style>
    html,
    body,
    #cesium_container {
      width: 100%;
      height: 100%;
      margin: 0;
      padding: 0;
      overflow: hidden;
    }
  </style>
</head>

<body>
  <div id="cesium_container"></div>
  <script>
    Cesium.Ion.defaultAccessToken = <SET_YOUR_OWN_TOKEN_HERE>
    const disable_opt = ['timeline', 'animation', 'infoBox', 'navigationHelpButton']
      .reduce((opts, opt) => (opts[opt] = false) || opts, {})
    const viewer = new Cesium.Viewer('cesium_container', disable_opt)
    fetch('/project')
      .then(res => res.json())
      .then(({ anchor, url }) => {
        const position = Cesium.Cartesian3.fromDegrees(...anchor)
        const modelMatrix = Cesium.Transforms.eastNorthUpToFixedFrame(position)
        const tileset = new Cesium.Cesium3DTileset({ url, modelMatrix })
        tileset.readyPromise.then(() => viewer.zoomTo(tileset))
        viewer.scene.primitives.add(tileset)
      })
  </script>
  </div>
</body>

</html>

在同一目录下执行如下命令启动服务后，浏览器访问localhost:5500即可浏览数据。
```
node index.js
```

级别	效果展示
原始数据
25级
24级