关于Ganos

Ganos是阿里云数据库产品事业部联合阿里云飞天数据库与存储实验室联合共同研发的新一代云原生位置智能引擎，它将时空数据处理能力融入了云原生关系型数据库PolarDB、云原生多模数据库Lindorm、云原生数据仓库AnalyticDB和云数据库RDS PG等核心产品中。Ganos目前拥有几何、栅格、轨迹、表面网格、体网格、3D实景、点云、路径、地理网格、快显十大核心引擎，为数据库构建了面向新型多模多态时空数据的存储、查询、分析、服务等一体化能力。

本文介绍的地理网格引擎能力，依托阿里云云原生关系型数据库PolarDB、云原生数据仓库AnalyticDB和云数据库RDS建设输出。

关于Ganos地理网格引擎

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

Ganos地理网格引擎目前涵盖GeoSOT和H3两种地理网格。GeoSOT是我国提出的一套地球空间剖分理论，并在此基础上发展出的一种离散化、多尺度区域位置标识体系。H3是Uber研发的一种覆盖全球表面的二维地理网格，采用的基本网格是正六边形。利用Ganos地理网格的函数可以将不同的空间范围转换为网格编码，可以求出网格编码的空间范围、层级和父子网格。Ganos支持退化网格计算（如下图），即充分利用网格的层级关系，用更精简的网格组合对空间范围进行表达。此外，Ganos自研地理网格索引，可用于高效查询网格码以及加速聚合计算。

Ganos地理网格支持无人机路径规划能力实践

业务需求

随着新能源技术的迅猛发展，低空经济已经逐步成为新的战略性新兴产业。低空经济是以各种有人驾驶和无人驾驶航空器的各类低空飞行活动为牵引，辐射带动相关领域融合发展的综合性经济形态，它的应用场景丰富，既包括传统通用航空业态，又融合了以无人机为支撑的低空生产服务方式，在工业、农业、服务业等领域都有广泛应用，对构建现代产业体系具有重要作用，发展空间极为广阔。2024年，深圳率先发布《深圳经济特区低空经济产业促进条例》，为低空经济起飞铺就法治跑道；上海多年来积极打造华东无人机基地，目前华东无人机基地已先后成为国家级特色产业集群、上海市无人机特色园区、中国民航首批无人驾驶航空试验区、虹桥国际开放枢纽“华东无人机空港”等。

不同于传统的地表活动，低空活动具有立体性、区域性、融合性等特点，这些特点对于如何安全引导低空活动的顺利开展带来了一系列需要解决的技术问题。

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

通过低空活动的特点我们可以发现，低空区域是没有路网的，在三维空间下如何保障低空设备的安全通过与到达、如何合理的划分低空区域的活跃周期与飞行器的活动规划是低空经济可以顺利落地并保障区域性活跃与安全的重要工作；

功能介绍

Ganos地理网格引擎提供了基于网格的路径规划能力，可以使用DEM、DSM、倾斜摄影等数据构建复杂环境下的无人机路径规划应用场景，相关能力简介如下：

ST_SetCost

用于给障碍物设置通行开销，障碍物可以是实体比如地面、建筑物，也可以是空间范围比如雷达扫描范围、风场空间。使用时首先要将障碍物的空间范围转换为网格数组，再进行通行开销设置。

ST_3DGridPath

无人机路径规划的最核心函数，根据用户提供的起始点终止点以及通行的开销数据，计算出最优的三维网格路径。支持不同的路径算法比如A*和DIjkstra，不同的移动方式比如是否允许斜向飞行，以及不同的距离计算方法比如曼哈顿和欧几里得。

同时，Ganos地理网格引擎还提供了ST_MatchGridLevel（依据给定的DES/DSM计算合适的用于路径规划的网格级别）、ST_CostUnion（合并所有障碍物通行开销到一个统一的开销数组，作为最终规划路线的输入）等能力应对更为复杂的业务场景需求。

以上函数的详细使用请参见Ganos手册，手册计划2024年3月下旬上线。

最佳实践

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

数据概况

测试数据形如：

数据入库与处理

网格路径规划功能依赖的扩展需要在使用前创建：

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

该数据存储于 OSS 中，使用 ST_ImportOSGB 将其导入（导入函数的详细介绍可以参考这里）：

SELECT 
  ST_ImportOSGB(
    'test_oblique',
    'OSS://YOUR_OSS_AK:YOUR_OSS_SK@OSS_REGION/YOUR_BUCKET/YOUR_OSGB_PATH', 
    '{"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 轴向上。
倾斜数据往往以相对坐标记录点，因此需要将每个点都以锚点（ anchor ）做偏移（ ST_Translate ），使每个点的坐标都为绝对坐标，并通过 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_ZupToYup 函数将左手坐标系下的数据转换到 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：每一个网格中心点投影到 ECEF （ 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。

ST_3DGridPath 函数支持如下自定义参数：

algorithm

使用的算法。可选项为：

dij : DIjkstra算法，基础算法，较慢。
astar : A*算法。默认算法，一般速度最快。
nb_astar : 新双向A*算法，在数据量较大时可能会快于A*算法。

movement

移动方式。可选项为：

cross : 只允许向相邻格网移动，默认移动方式。
octothorpe : 允许向相邻格网及对角格网移动。
strict_octothorpe : 仅在相邻格网可通行时向对角格网移动。

distance

距离估算方式。可选项为：

euclidean：欧氏距离，默认距离估算方式。
manhattan：曼哈顿距离。
chebyshev : 切比雪夫距离。

三种距离计算方式代表着对于对角线移动距离的不同评估。

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

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

多成本区网格路径规划

在可通行区为不同区域设置不同通行成本后，算法将选择综合成本最小的路径。我们在上文场景中添加一个高成本通行区作为多成本网格路径规划的样例。

在 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内部函数完成，也提供了完备的路径规划的辅助函数；
无人机路径规划可以使用退化网格技术，节省存储开销，加快计算速度；
路径规划提供多种参数选择，包括规划算法、移动方式和距离计算方法，以满足不同需求；

总结

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

瑶池解决方案体验馆

欢迎大家通过瑶池解决方案体验馆免费体验阿里云瑶池数据库系列产品独特的能力与优秀的价格力，瑶池解决方案体验馆，真免费体验馆~

截屏2024-03-15 10.29.49.png

Ganos地理网格引擎支撑无人机路径规划能力实践