Voronoi图是地学计算与空间分析的一个基本结构^[1-2]，在平面空间和椭球面上^[3-4]，Voronoi图在诸如骨架提取、过程模拟和空间优化等众多方向上有着重要且成熟的应用^[5-7]。然而在不规则地形曲面上，由于测地距离映射这个基础性科学问题的挑战，Voronoi图还缺乏严密的基础^[8-9]。实践中经常使用欧氏距离近似方法或栅格聚类方法来逼近地形曲面的测地Voronoi图（geodesic Voronoi diagram，GVD） ^[10-11]。欧氏距离近似的一个直观误差是平坦化的Voronoi单元和缩小的表面积；栅格Voronoi图则具有锯齿状边界，会影响地表优化模型的效率甚至阻碍收敛^[12-13]。

空间曲面上的度量问题吸引了科学界的广泛关注^{[9, 14-15]}，离散曲面上的测地距离映射大体上分为快速的微分几何方法和精确的计算几何方法两大类^{[8, 16-17]}。微分几何的方法将测地距离映射看作光或热的光滑函数，因而只有一阶精度并且不能完全满足度量要求^{[16, 18]}。计算几何的方法将来自同一源点、同一方向且空间连续的测地线封装为边上的可视窗口，计算窗口在地形格网上传播时的几何关系，得到精确的测地距离场（geodesic distance field，GDF）^[15]。理论上，可以依据GDF将格网上的任一点归纳到距离最近的一个源点，即GDF蕴涵测地Voronoi图。但是，在不同源点的测地线相遇的目标面片上，微分几何的GDF没有提供对面片的平分线划分（平分线通常是双曲线结构），计算几何的GDF以线性的可视窗口回避了对面片的进一步划分^[12]。因而，GVD的边、顶点等关键信息都还隐藏在目标面片中。

文献[19-20]对真实地形上的Voronoi图进行了理论分析，指出目标面片的单条边在最坏情况下会被平分线划分n次（n为面片数），定性地给出了目标面片上双曲线划分的复杂度。文献[12]就三角化的2-流形表面模型给出GVD的数学描述。在计算几何的框架上，通过预先细分使目标面片至多包含3个源点，给出一个基于测地距离映射的GVD集成算法^{[12, 21]}。但可视窗口对格网的隐式划分严重影响了精确的计算几何方法的效率^{[8, 16]}；同时，基于挑战性测地距离映射框架的集成算法非常复杂，导致其很难被复制；集成算法缺乏Voronoi图核心的空间划分及双曲线边重建的清晰脉络，而这实际上却可以通过生长代表着源点间平分线的边上奇异来显式确定。

鉴于此，本文提出一个避免过度细分预处理的GVD直接算法，主要包括以下3方面的工作：（1）在离散地形曲面上引入计算几何的精确GDF；生长代表源点之间平分线的边上奇异，完成对目标面片的精确划分；（2）在目标面片划分的基础上，提出一个聚类-划分-再聚类的GVD直接算法；（3）以Voronoi单元的面积统计和优化构建不规则三角网（triangulated irregular network，TIN）为例，从定量和定性两方面对地形上的精确Voronoi图开展应用探索。

1 GDF的奇异拓扑

给定一个离散地形曲面M和M上的源点集S。定义M=（V，E，F）为一般意义上的三角网，V、E、F分别为顶点集、边集、面片集；定义S为限定在M上的点集，即源点不一定在顶点上。

对于M上任一点p，测地距离映射寻求一个距离函数D（p），它返回p与某一源点s∈S之间的最短距离。GDF以源点集为基准，这个浮动的场模型同时适合地形自身及地表分布的空间建模^{[12, 22]}。

以MMP（Mitchell，Mount，and Papadimitriou）算法^[23]为代表的计算几何方法主要从这样一个基本事实出发：测地线在某个面片上是直线；经过某条边时，如果将相邻面片沿着这条边展开到同一平面，则测地线也将是直线；同样，在某一时刻，来自同一源点的测地线前端会形成圆弧状的波前阵，除非它遇到了鞍点或其他的波前阵。

鞍点是地形曲面上的一类特殊的特征点，因为具有超过$ 2{\pi } $的锥顶角，鞍点附近的空间不能被同时展开到同一平面上。从几何上看，经过鞍点的测地线在下游形成对源点不可见的区域^{[17, 24]}。在这片区域，需要以鞍点为新的源点中继传播，形成分叉的波前阵和新的可视窗口。鞍点的作用与源点相似，又称为伪源点。鞍点对波前阵的延缓作用是平分线呈双曲线形状的起因^[12]。

在经典组合拓扑中，点、线、面是空间连接关系发生突变的位置。GDF上可以类比这一概念来建立奇异拓扑：GDF上同时具有两个以上源点的点为奇点；具有3个以上源点的奇点为GVD顶点；空间上连续的奇点形成奇异轨迹。因此，源点之间的平分线是奇异轨迹，GVD的边是由顶点截断的平分线。

测地线经过的顶点除了边界点外就是鞍点（伪源点），这导致了一种特殊的直线平分线。如图 1（a）所示，假设有伪源点P₀、P₁，则由波前阵的对称性可知，经过它们的测地线是直线。在这条直线上，每一点都具有P₀和P₁两个不同的伪源点。射线P₀P₁是“平分”来自同一源点的不同方向波前阵的奇异轨迹，称为侧向奇异轨迹。

图  1  侧向奇异轨迹与前向奇异轨迹

Figure  1.  Lateral Singular Loci and Front Singular Loci

下载: 全尺寸图片 幻灯片

在一般情况下，奇异轨迹是双曲线。如图 1（b）所示，来自不同源点P₀、P₁的波前阵首先相遇于边上的$ {S}_{0} $点，$ {S}_{0} $显然是一个奇点，则有：

式中，$ \left|{P}_{1}{S}_{0}\right| $是欧氏距离；D（P）是伪源点的测地距离。式（1）定义了一条双曲线：与两个定点（即焦点P₀、P₁）相差固定距离的点的轨迹（即弧段$\overparen{S_0 I_0}$），被称为前向奇异轨迹。

记$ D\left({P}_{1}\right)\mathrm{―}D\left({P}_{0}\right)=k $，$ {P}_{1.x}=a $，$ {P}_{1.y}=b $，$ {P}_{0.x}=c $，$ {P}_{0.y}=d $，式（1）可以写成双曲线的标准形式：

式（2）为从边上生长代表平分线的奇异轨迹提供手段。同时，这种解析的空间划分与隐式的拓扑信息是计算几何方法优于微分几何方法的关键所在。

2 奇异生长的GVD直接生成算法

2.1   算法概述

假设地形格网上有n个面片和m个源点，根据欧拉公式，GVD边的数目的复杂度为$ O\left(m\right) $。一般地$ m\ll n $，所以大部分面片不涉及GVD边的非目标面片。算法的第一步是过滤这些只含有一个源点的面片，并以这些源点为标签进行聚类。

算法的第二步以双曲线完成对目标面片的划分。多个伪源点在目标面片中竞争空间，每一个源点都要同其他源点达成空间均衡。但边上奇点只能生成空间相邻源点间的均衡平分线，因此从边上开始由近及远进行奇异生长，确定非直接相邻源点间的新奇点是整个算法的关键。

文献[12, 23]证明，由双曲线划分的子面片必定和边上的某个可视窗口直接连通（即不存在“飞地”）。于是算法第二步的划分子面片可以由边上窗口的源点标识，算法第三步以这个标签再次聚类，最后提取GVD的拓扑结构。

2.2   奇点的奇异生长

测地线波前阵首先在目标面片的边上相遇，即可视窗口的两个端点奇点。基于一个奇点的左右源点与这个奇点自身，可以由式（2）构造一条双曲线，这条双曲线将构成GVD边的一部分。

由式（2）构造3点双曲线需要定位。首先，可视窗口以边为定位定向，算法需要把所有窗口的定位定向变换到统一的局部坐标系；其次，奇异轨迹也需要定向。如图 1（b）所示，只有经过当前奇点、限定在目标面片上的那一部分双曲线才是需要的，即弧段$\overparen{S_0 I_0}$。这里$ {S}_{0}{I}_{T} $是过$ {S}_{0} $的切线，$\overparen{S_0 I_0}$的定向由$ \stackrel{⃑}{{S}_{0}{I}_{T}}\times \stackrel{⃑}{{S}_{0}{I}_{0}} $给出。

2.3   可视窗口的奇异生长

双曲线划分的主要挑战是因为从边上生长的奇异轨迹以不同的次序进入目标面片，并相互缠绕。同时，可视窗口是不同源点波前阵的线性截断，不同源点在目标面片上的空间竞争需要可视窗口之间达成空间均衡。

因而这里以边上的奇点生长来定义可视窗口的生长：可视窗口的奇异生长以它两端的奇点生长为开始；若从两端生长的双曲线在目标面片上存在有效的近端交点，则可视窗口的生长得到满足。有效近端交点的定义是：如果左曲线的第一个交点同时也是右曲线的第一个交点，则称为存在有效近端交点；如果这个存在的交点同时距离当前源点更近，则确定为有效近端交点。交点和一次交点可由x单调的双曲线排列（如CGAL）得到。

从算法实现来看，可视窗口的生长流程描述如下：若当前的奇点轨迹与左右窗口的奇点轨迹都存在一次交点，则距离近的那个点是有效近端交点，新奇点生成，新窗口代替对应的旧窗口，新弧段补全隐藏的GVD边的一部分（示例见下节）。

2.4   GVD的直接生成算法

通过窗口的奇异生长能逐步确定有效近端交点，建立新奇点和新窗口，由近及远地解缠双曲线的复杂划分，最终完成可视窗口所代表的波前阵生长和目标窗口的精确划分。

算法唯一使用的复杂数据结构是可视窗口，参照文献[12, 17]对可视窗口的定义，围绕可视窗口生长的直接GVD算法描述如下：

1）对所有面片检查其源点数，若只含有1个源点，则以该源点ID标记并聚类；若含有2个以上源点，则将该目标面片压入队列。

2）对目标面片队列中的每一个面片执行循环：①对于所有窗口，由端点生长成对的双曲线；②在双曲线形成的边界面片排列中，为每一个窗口寻找有效近端交点；③对于每一个有效近端交点，连接其原窗口两个端点和此交点形成两个新窗口，代替原窗口，进入下一轮循环；④检查是否还有新交点。

3）以源点ID聚类所有面片，抽取GVD矢量结构。

直接GVD算法的关键是以可视窗口的生长完成GVD边的显式构造，本质上是以一个队列循环取代文献[12, 18]的递归过程，避免过度的细分预处理。

图 2展示了一个相对简单的目标面片上的双曲线划分。面片$ ABC $上共有17个可视窗口，2对边上的奇点（S₀₁、S₁₁和S₁₄、S₂₅）具有不同源点，是GVD的边经过的奇点。

图  2  目标面片上的双曲线划分与边上的奇点

Figure  2.  Conic Tessellation and Singularities on Edge of a Target Facet

下载: 全尺寸图片 幻灯片

图 3展示了自奇点S₀₁生长的双曲线，它与直接相邻的奇点S₁₁/$ {B}_{0} $的双曲线各有一个交点I₀₁/I₀₂。I₀₁/I₀₂都是一次交点，但因为I₀₁离S₀₁更近，所以I₀₁为有效近端交点。这里顶点B同时是鞍点和奇点，其中继的波前阵在目标面片上形成分叉$ {B}_{0}/{B}_{1} $。有效近端交点$ {I}_{01} $将替代原奇点$ {S}_{01} $成为新奇点，其右窗口更新为$ B{I}_{01} $。

图  3  有效近端交点计算与可视窗口奇异生长

Figure  3.  Effective Near-End Intersection and Singular Growth of Visibility Window

下载: 全尺寸图片 幻灯片

图 4展示了整个算法的流程。非目标面片以聚类类别标识，将含有2个源点的目标面片渲染为土红色，含有3个以上源点的目标面片渲染为绿色，GVD边渲染为白色。

图  4  算法流程

Figure  4.  Flowchart of the Algorithm

下载: 全尺寸图片 幻灯片

3 算法验证与应用

3.1   算法对比与验证

GVD算法目前还没有公开的源码，但文献[12]提供了一个演示程序和实验数据。在这个实验数据上，目标面片20上的奇异生长结果和文献[12]的结果对比如图 5所示。注意文献[12]对面片20的定位与本文不一样，但两者定向是一致的（顶点对应关系：A-4/B-5/C-13）。图 5（a）中的前向奇异被渲染为蓝色，侧向奇异被渲染为绿色。

图  5  验证数据目标面片20上的奇异生长对比

Figure  5.  Singular Growth on Facet 20 of the Validation Data

下载: 全尺寸图片 幻灯片

对比图 5（b）与图 5（c）可以看出，本文目标面片的划分结果与文献[12]的划分结果一致。这种解析结果的吻合证明了奇异生长方法的可靠性。

3.2   算法效率分析

如§2.1所述，可将GVD直接生成算法分为3个步骤：（1）非目标面片的聚类，算法效率很高；（2）使用窗口的奇异生长直接得到部分GVD边并精确划分目标面片，避免了文献[12, 25]过度的递归细分和GVD结构的动态维护；（3）再次采用成熟的聚类算法并提取GVD结构，整体算法直观易行。

GVD的生成效率与底层的GDF映射算法密切相关。本文与文献[12, 25-26]基于经典MMP算法^[17]，窗口分支的开销较大；改进的CH（improved CH algorithm，ICH）算法^[27]通过过滤无效窗口，大大节省内存，并使算法逼近理论时间复杂度$ O\left({n}^{2}\right) $；鞍点图（saddle vertex graph，SVG）算法^{[8, 24]}充分发展计算局部化技术，能够在多核并行架构上处理千万数量级点云的GDF，但是其预处理的基础仍然是MMP类的精确算法。本文主要从原理上说明并验证GVD，提出的直接算法受MMP算法的内存局限性影响。

算法时间主要消耗在3个方面：（1）MMP算法的可视窗口以边定向定位，对目标面片上的可视窗口及其源点需要进行刚性平移与旋转，统一变换到面片的局部坐标系；（2）双曲线划分与有效近端交点计算，其高精度的双曲线数字运算消耗最大；（3）双曲线对目标面片进行非常破碎的细分，对§2.4中步骤3）的聚类过程影响不大，但对GVD的拓扑抽取过程影响较大。

本文在Windows10/VC++ 2019/I7 8650U CPU/8 GB RAM的普通电脑上生成10万三角面片的GVD，耗时不超过2 min。

3.3   Voronoi单元的面积误差

Voronoi图的单元面积经常作为基础参数出现在分层地理统计^[28]（如地表生物质分类统计^[29]）、自然邻近插值（如通过激光雷达点云生产数字高程模型^[30]）和科学可视化^[31]等应用中。Voronoi图单元的面积误差会从空间结构和空间模式等不同层面影响分析应用的结果^[5]。

在地学分析中，高程经常如同顶点的属性一样使用^[32]，这时候使用的是平面Voronoi图和投影面积，高程的影响则以间接的参数化形式提供。地形Voronoi图有时采用Delaunay三角网对偶的近似方法，如著名的QHull软件包。由于Delaunay三角化实际上是将高程点投影到平面再优化，所以这是一种以欧氏距离为条件的近似优化，所得到的对偶Voronoi图是一种近似地形Voronoi图。栅格Voronoi图假设面片总是可以细分，是一种更谨慎的近似方法。

以图 4中的地形格网（具有1 092个顶点和1 990个三角面片）为例，按相同的5个源点产生5个Voronoi单元；以GVD的曲面面积为基准，统计GVD的投影面积、Delaunay三角网对偶的面积和投影面积、栅格Voronoi图的曲面面积，以及均方根误差（root mean square error，RMSE）如表 1所示。

表  1  GVD与近似Voronoi单元的面积比较/m²

Table  1.  Comparison of Areas Between GVD and Approximated Voronoi Cells/m²

方法	单元1	单元2	单元3	单元4	单元5	面积和	RMSE
GVD曲面	7 379.99	7 116.12	5 463.38	6 991.84	6 337.97	33 289.30	-
GVD投影	5 743.07	5 424.64	3 348.78	5 348.23	4 935.36	24 800.08	1 713.53
Delaunay三角网对偶	6 479.36	5 809.39	2 918.87	6 312.48	6 259.35	27 779.45	1 375.57
Delaunay三角网对偶投影	5 772.80	5 120.10	2 744.92	5 584.25	5 578.00	24 800.07	1 817.47
栅格Voronoi图	7 407.72	7 131.24	5 447.30	6 980.24	6 322.83	33 289.31	18.00

下载: 导出CSV 

| 显示表格

从表 1中可以看出，Delaunay三角网对偶及Delaunay三角网对偶投影的单元面积误差数量级相当（RMSE都在1 300 m²以上），这是因为它们本质上都是欧氏距离的近似度量。直观上，欧氏距离近似导致曲面平坦化和Voronoi单元面积缩小；数值上，本例中平坦化的Voronoi单元面积的偏差达到19.85%，几乎达到地形粗糙度指数（34.23%）^[33]的2/3。

本例中栅格Voronoi图具有比较小的RMSE，这主要是因为源点数（5）远小于面片数（1 990）。这说明如果误差要求不高且距离度量得当，栅格Voronoi图通常是一个好的近似。但也必须注意到，栅格Voronoi图单元具有锯齿状边缘^[34]，这种“好的近似”在更精巧的空间优化场合不仅会影响迭代效率，还会影响空间优化的结果^{[12, 13]}。

3.4   优化构建TIN

以一组采样点为源点构建GDF和GVD，则GVD单元所对应的是采样点沿地形表面的外接空圆。这样，从采样点优化构建TIN的一个计算困难的判断采样点最小邻接空圆问题已经得到解决^{[10, 35-37]}（在非欧空间中搜索最邻近的k个点需要超线性时间），这是GVD最直接的应用。

实验2数据是由美国St. Helens火山10 m水平分辨率的DEM处理而来的TIN，具有17 222个面片和8 800个高程点。使用100个相对均匀的随机点作为源点，得到GVD图，如图 6所示（源点以红色点标示，地形以高程渲染），放大的GVD曲线格网（绿色）与地形格网的局部如图 6（b）所示。

图  6  实验2的地形格网、源点和GVD格网

Figure  6.  TIN Grids, Source Points, and GVD Grids of Experiment 2

下载: 全尺寸图片 幻灯片

以GVD单元为指示，产生GVD的对偶TIN，如图 7（a）所示（白色格网）。作为对比，以同样的100个源点为顶点，使用2D的Delaunay优化构建TIN，结果如图 7（b）所示（蓝色格网）。图 7（b）中未被蓝色格网覆盖的7处白色线段指示出传统的Delaunay优化可能不符合最小角最大化的错分情形，底图为被渲染成不同颜色的GVD单元。

图  7  GVD的对偶与传统2D Delaunay三角化的潜在错分

Figure  7.  GVD Duals and the Potential Wrong Triangulations of Traditional 2D Delaunay

下载: 全尺寸图片 幻灯片

使用不同数目（30、50、80）的源点产生均匀与不均匀分布的GVD，都能重现上述2D Delaunay的潜在错分。究其根源，是因为GVD能够提供沿地形的外接空圆，而2D Delaunay外接空圆是基于欧氏平面空间条件判定的。

4 结语

本文在离散地形曲面上引入GDF和发展奇异生长的方法来应对Voronoi图的精度问题和复杂算法的挑战。通过奇异生长直接得到部分GVD边结构并突破Voronoi图核心的空间划分的技术屏障；经过常规的面片聚类后再抽取GVD，避免了目标窗口过度的递归细分，降低了动态维护GVD双曲线结构的复杂性；整体算法直观易行。以GVD为基准，实验中发现近似Voronoi图的单元面积误差达到地表粗糙度指数的数量级，而传统的Delaunay优化构建TIN会出现不少不符合最小角最大化原则的潜在错分。这表明GVD可以为数字地形分析中的许多应用带来基础性改进。

在地形曲面上开展基于GVD的自然邻近插值、表面过程模拟与空间优化应用探索，将是下一步值得深入研究的内容。