从计算几何^[1]的观点来看，椭球面三角形外心就是与3个相邻顶点所对应的椭球面Voronoi图的公共顶点，而Voronoi图具有许多优秀的性质，研究椭球面Voronoi图对于全球数据的管理和空间关系的动态维护具有重要意义^[2]。根据Voronoi图和Delaunay三角网的对偶关系，采用间接Delaunay三角剖分算法由椭球面Voronoi图能够比较容易地构造椭球面Delaunay三角网，而Delaunay三角网被广泛应用于数字高程模型(DEM)构建、空中三角测量、大地控制网评估和连续运行参考系统(CORS)网布站优化^[3]。因此作为生成椭球面Voronoi图并构造相应的Delaunay三角网的关键步骤，研究椭球面三角形外心的求解具有重要的理论和应用价值。

在平面上由3个相邻顶点计算三角形外心的坐标可采用确定的公式和成熟的算法^[4]，然而大地线距离是由在椭球面某范围内足够精确且不产生二义的大地线尺度函数所定义的，直接在椭球面上计算三角形外心大地坐标异常复杂^[5]。因此在传统的地理信息系统中大多数采用球面或投影平面对椭球面进行高或低精度的近似以降低实现难度^{[6, 7]}，例如在球面上10³km跨度内误差达到了35m，而在投影平面上10³km跨度内误差竟达到了8km^[7]，可见两者的误差控制情况都不是很理想。

文献[4]以构建区域性椭球面数字地面模型(DTM)为目的，提出了一种在测地坐标系^[8]中由椭球面三角形3个相邻顶点求解其外心测地坐标的方法。对于平均边长在1～2km的椭球面三角形，计算所得外心到该三角形各顶点的大地线距离偏差小于0.1mm^[4]。因为测地坐标系是基于区域性椭球面建立起来的，适用范围仅限于较大的局部区域^[5]，所以该方法对于跨度远远超过2km甚至达到10³～10⁴km的很多实际问题，如无人机航迹规划^[9]、巡航导弹突防路径规划^[10]以及海域划界中间线及其三接点的求解^[11]等，无法保证其结果的高精度。

文献[12]以求解海域划界中间线上的转折点位置为目的，提出了一种在球面和椭球面上求解三角形外心大地坐标的初步思路，但是未见该方法的有效试验结果。对于椭球面三角形外心大地坐标的求解，首先要计算得到相应的球面三角形外心的经纬度坐标，然后采用求解线性方程组的一般迭代方法逼近外心的精确大地坐标，然而迭代过程中产生的向量序列并不一定逐步逼近线性方程组的解，且不同迭代方法的收敛效率差别较大^[13]。

为了能够在全球范围内高效计算椭球面三角形外心的精确大地坐标，解决上述方法各自存在的问题，本文基于地图代数^[14]中的测地变换^[15]，提出了一种求解椭球面三角形外心的矢栅结合方法。算法可以直接在椭球面上快速计算大量离散点集的高精度Voronoi顶点，对于椭球面Voronoi图的生成和椭球面Delaunay三角网的构造具有重要作用，从而为建设动态全球性地理信息系统^{[16, 17]}(Global GIS)并实现精确椭球几何计算^[18]与空间分析^[19]打下基础。

本文方法基于地图代数理论，利用“矢栅结合使用长处”^[20]的思想，通过栅格大地线尺度距离变换——测地变换，先精确求解椭球面空间中各栅格中心点到其最近地理实体上来源点的距离，后在临界相邻栅格范围内，以栅格中心点作为控制点，通过数值计算内插初始等距点大地坐标并最终逼近椭球面三角形外心的精确位置。算法的误差主要来自于在椭球面3～4个临界相邻栅格范围内由初始等距点逼近三角形外心所导致的精度损失，下面从椭球大地测量学^[21]和球面三角学^[22]的角度论证：在满足一定条件下上述逼近过程与大地主题解算的精度^[23]能够保持一致。

设K_A、K_B和K_C分别表示椭球面三角形ABC各顶点的高斯曲率，R_A、R_B和R_C分别表示对应顶点纬度处的曲率半径，是3个顶点的平均高斯曲率，是3个顶点平均纬度处的平均曲率半径，如图 1所示，以 为半径作密切球^[21]，在球面上做与三角形ABC对应边长度相等的三角形A₀B₀C₀。

图 1 椭球面三角形ABC与球面三角形A0B0C0 Fig. 1 Triangle ABC on ellipsoidal surface and triangle A0B0C0 on spherical surface

图选项

对于椭球面三角形ABC，由椭球面正弦定理^[21]可得

点击浏览公式

对于球面三角形A₀B₀C₀，由球面正弦定理^[28]可得

点击浏览公式

比较式(1)和式(2)可得

点击浏览公式

同理可得

点击浏览公式

因为式(3)、式(4)和式(5)中 和 都是微小量，所以 ，且 ，其中ε是球面角超，以同样的精度得出

点击浏览公式

式(6)为椭球面改正数公式，当椭球面三角形ABC的平均边长小于200 km时，椭球面改正数小于0.001″，而当200km线长的方向变化0.001″时，端点的位移小于0.001m。因此在平均边长小于200km的条件下可直接把椭球面三角形当作球面三角形，而对于平均边长在200～400km 之间的椭球面三角形需要经过适当改正转化为球面三角形^[21]。

如图 2所示，椭球面三角形DEF和DEG共用同一边DE，且大地线DF和DG长度相等，∠EDF和∠EDG也相等。对应到本文方法中，D表示初始等距点，F和G分别表示相邻两个栅格中心点的来源点，E表示∠FDG平分线(大地线)上的一点。当DF和DG的长度不超过10⁴km且DE的长度不大于1km时，如果能够论证EF和EG的长度偏差小于0.001m，由初始等距点逼近三角形外心与大地主题解算的精度^[23]相当的结论就是成立的。在密切球上做与椭球面三角形DEF和DEG对应边长度相等的三角形D₀E₀F₀和D₀E₀G₀，已知大圆弧D₀F₀和D₀G₀长度相等，假设∠E₀D₀F₀和∠E₀D₀G₀也相等，即两边及其夹角彼此相等。根据球面三角形全等定理^[22]，球面三角形D₀E₀F₀和D₀E₀G₀全等，由全等球面三角形的性质可知大圆弧E₀F₀和E₀G₀的长度也相等。实际上，将椭球面三角形转化到对应球面上以后，∠E₀D₀F₀和∠E₀D₀G₀并不相等，因此大圆弧E₀F₀和E₀G₀的长度不相等。

图 2 椭球面三角形DEF和DEG与球面三角形D0E0F0和D0E0G0 Fig. 2 Triangle DEF and DEG on ellipsoidal surface and triangle D0E0F0 and D0E0G0 on spherical surface

图选项

由式(3)可知，椭球面上的∠EDF与对应球面上的∠E₀D₀F₀仅相差一个椭球面改正数，其大小主要与椭球面上大地线DE和DF长度的乘积有关。因为大地线DF的长度不超过10⁴km，且大地线DE的长度不大于1km，所以椭球面改正数小于0.001″，由球面余弦公式^[22]可知对应球面上大圆弧E₀F₀端点的位移小于0.001m，同理可知大圆弧E₀G₀端点的位移也小于0.001m，即大地线EF与EG的长度偏差在0.001m以内可控，也即逼近过程与大地主题解算的精度^[23]保持一致。

一般情况下大地线DF、DG和DE的长度越短，大地线EF与EG的长度偏差就越小，逼近的精度也就越高。在本文方法中DF、DG的大地线长度对应于外心到椭球面三角形顶点的大地线距离，DE的大地线长度对应于栅格单元的跨度。因此理论上算法的精度与大地线距离的长短呈负相关，与栅格分辨率的高低呈正相关。

本文方法基于地图代数测地变换，充分发挥了矢栅紧密结合数据结构的优势，在保证栅格计算效率的条件下，能够计算定位误差小于0.001m的大跨度(10³～10⁴km范围)内椭球面三角形外心大地坐标。按照算法流程的先后顺序，主要涉及以下几个关键技术：测地变换、边界跟踪配对和求解外心大地坐标,其中边界跟踪配对用来确定外心所在的临界相邻栅格范围。

所谓测地变换是指在椭球面上进行的栅格大地线尺度距离变换，其目的是建立椭球面空间的高精度距离场。地图代数测地变换的时间复杂度为O(m)，其中m为计算规模，是全空间的栅格单元总数^[15]，等于横向分辨率(地理网格阵列的列数)与纵向分辨率(地理网格阵列的行数)的乘积。与利用穷举算法进行距离变换相比，地图代数测地变换完全不受生成元个数和生成元复杂程度的影响，不但能够适用于复杂的自然图形，而且可以大大提高计算效率，适用于本文问题的解决。作为大地线尺度下的空间度量，测地变换的精度最终取决于大地主题解算公式的精度，在20000km的范围内误差小于数厘米^[23]，因此可以认为栅格中心间的大地线距离是没有偏差的，在全球范围内可用。

文献[15]规定测地变换必须满足3个基础定义，其中定义3明确指出两定点a、b间的栅格距离等同于a、b所属栅格中心间的距离，顾及整个栅格和舍入误差，有不大于0.5个栅格单元的长度^[15]。虽然该算法的精度可以由栅格分辨率来控制，理论上能够达到很高的精度，但是对于某些实际应用，例如要在10³～10⁴km范围内计算误差小于0.001m的椭球面三角形外心，各栅格单元的跨度都不能超过0.001m，相应的横向和纵向分辨率将达到10⁹～10¹⁰，全空间的栅格单元总数就增加到了10¹⁸～10²⁰，现阶段普通的计算机显然无法处理规模如此庞大的地理网格阵列。

文献[24]提出了一种基于八邻域边界跟踪的物体标号算法，一次扫描即可同时获得物体的边界点序列以及边界链码，为后续处理提供数据准备。在此基础上，本文通过以下方法对测地变换所得分配场进行边界跟踪配对：首先找到最左上角的一个边界栅格作为搜索起点，按顺时针方向，自上而下、从左至右，搜索其八邻域，并将搜索路径上的四邻域栅格作为上一个边界栅格的配对栅格，直到找到下一个边界栅格；然后以此边界栅格为新的起点继续搜索，这一搜索过程不断重复下去，直至回到搜索起点。经过边界跟踪配对以后，所有边界栅格的邻接匹配关系都可以被确定。遍历所得边界栅格对序列，按顺序从前往后每两对一组进行判别扫描，过滤出匹配关系发生变化的相邻两对边界栅格对，即为椭球面三角形外心所在的栅格范围。

经过边界跟踪配对确定椭球面三角形外心所在的栅格范围以后，最终求解足够精确的结果又可分为两个步骤：初始等距点的内插以及外心大地坐标的逼近，它们分别如图 3和图 4所示。其中O₁和O₂分别表示上述相邻两对边界栅格对中的配对栅格中心点，B₁和B₂分别表示O₁和O₂的来源点，B₁到O₁的大地线距离和B₂到O₂的大地线距离分别等于S₁和S₂，O₃和O₄分别表示与O₁和O₂相互对应的边界栅格中心点，B₃表示O₃和O₄的同一个来源点，B₃到O₃的大地线距离和B₃到O₄的大地线距离分别等于S₃和S₄。

图 3 内插初始等距点 Fig. 3 Interpolation of the initial equidistant point

图选项

图 4 逼近椭球面三角形外心 Fig. 4 Approximation of the circumcenter of triangle on ellipsoidal surface

图选项

经过大地线尺度下的地图代数测地距离变换，可以确保O₁距离B₁比B₂更近，O₂距离B₂比B₁更近。由大地线长度变化的连续性容易证明，在如图 3所示的O₁与O₂连接成的大地线段上必然存在一点O到B₁和B₂的大地线距离都等于S，即得到初始等距点O。

在10³～10⁴km范围内，当栅格单元的跨度不大于1km时，在如图 4所示的3～4个临界相邻栅格范围内，∠α的平分线(大地线)上任意一点P到B₁和B₂的大地线距离偏差都小于0.001m。由大地线长度变化的连续性同理可证，在该角平分线上必然存在一个到B₁和B₃的大地线距离相等或到B₂和B₃的大地线距离相等的点，不妨假设该点为P，则P到B₁、B₂和B₃的大地线距离S′、S″和S’’’的偏差都小于0.001m，即求解得到由B₁、B₂和B₃构成的椭球面三角形外心的精确位置。

充分利用测地变换所生成的高精度距离场，通过数值计算先内插初始等距点，再以同样的方式不断逼近椭球面三角形外心是本文方法的关键。原则上在10³～10⁴km范围内要使得各栅格单元的跨度都不超过1km，因此一般情况下将横向分辨率范围设定为10³～10⁴，纵向分辨率范围也随之确定，相应的计算规模(全空间的栅格单元总数)控制在10⁶～10⁸之间。本文方法的具体实现步骤如图 5所示，其中源点坐标场记录的是各控制点的最近来源点经纬度坐标信息，距离场记录的是每个控制点到对应来源点的大地线距离信息，分配场则记录了所有控制点与其最近地理实体之间的对应关系。

图 5 椭球面三角形外心的求解流程 Fig. 5 Flow chart of determination of circumcenter of triangle on ellipsoidal surface

图选项

为了保证方法在全球范围内的普遍适用性，本文采用了WGS-84坐标中的椭球基本参数：椭球长半轴a= 6378137.0m，椭球扁率f= 1.0/298.257223563，且采用了解算精度与距离无关的贝塞尔大地主题解算公式^[25]。

进行以下10组数据的计算，其中横向分辨率都设定在5000以内，各椭球面三角形3个顶点及其外心的位置如图 6所示，具体的大地坐标见表 1。为了验证所得结果的正确性，分别计算各外心到对应椭球面三角形3个顶点的大地线距离及其标准差，详细结果见表 2。在表 1和表 2中，λ都表示经度，Φ则表示纬度。

图 6 椭球面上各三角形及其外心与外接大地圆 Fig. 6 Triangles on ellipsoidal surface and their geodetic circumcenters and circumcircles

图选项

由表 2可知，随着外心到对应椭球面三角形各顶点的大地线距离变得越来越远，标准差呈现不断增大的趋势，说明本文方法的计算精度与大地线距离的长短呈负相关。1—9号椭球面三角形外心到对应各三角形3个顶点的大地线距离之差都小于0.001m，其中外心到对应三角形顶点的最短大地线距离为307648.8577m，标准差不超过0.00006m；最长大地线距离达到3827513.4529m，标准差也不超过0.0005m。10号椭球面三角形外心到对应各三角形3个顶点的大地线距离分别是5646875.7513m、5646875.7518m和5646875.7539m，虽然两两之差不都小于0.001m，但是标准差不超过0.002m。因此在中、长距离情况下，本文方法对于椭球面三角形外心的求解精度非常高。在超长距离情况下，计算精度虽有所降低，但仍然保持在0.001m左右。

进行以下5组数据的计算，其中构成椭球面三角形ABC各顶点的大地坐标分别是A(88.1314°,55.2526°)、B(100.3333°,40.5152°)和C(111.6666 °,66.3738 °)，分别计算不同分辨率下外心的大地坐标和相应外心到各顶点的大地线距离及其标准差，详细结果见表 3，在表 3中λ和φ分别表示经度和纬度。

由表 3可知，随着分辨率变得越来越高，标准差呈现不断减小的趋势，说明本文方法的计算精度与分辨率的高低呈正相关。当横向分辨率从1000逐渐增大到4000时，标准差从0.0015567060m逐渐减小到0.0002081666m，计算结果的精度在不断提高。当横向分辨率继续增大到5000时，标准差反而增加到0.0003214551m，说明当分辨率增加到一定程度后，对于算法精度的提高作用是有限的，由于大地主题解算的精度限制^[23]，此时大地主题解算的误差对结果的影响超过逼近误差，精度发生异常波动属于正常现象。因此本文方法能够在可计算的栅格分辨率范围内(横向分辨率在10³～10⁴之间，全空间的栅格单元总数在10⁶～10⁸之间)高效计算定位误差小于0.001m的大跨度(10³～10⁴km范围内)椭球面三角形外心大地坐标。

由于本文方法在划分网格单元的步骤中采用的是基于大地坐标系统的传统平面格网系统构建方法，不可避免地存在面积变形、形状变形以及由赤道向南北极递增的问题^[26]，因此分析方法精度与椭球面三角形各顶点所处纬度的关系也很有必要。观察表 2中的6号和7号椭球面三角形，两者的横向分辨率都设定为4000，前者的外心到各顶点的大地线距离分别是1680030.9983m、1680030.9983m和1680030.9991m，标准差是0.0004618801m，后者的外心到各顶点的大地线距离分别是2092296.4951m、2092296.4942m和2092296.4946m，标准差是0.0004509249m，同属中长距离情况，计算精度基本相当。由表 1可知，6号椭球面三角形各顶点处于低纬度地区，而7号椭球面三角形各顶点处于中高纬度地区。因此在正确划分网格单元的条件下，本文方法的精度与所处纬度无关。

选取如图 7所示的位于大西洋两岸的100个城市结点进行试验。在横向分辨率设定为4000的条件下，采用配置为Intel Core 2 Duo E7500处理器和2GB内存的台式电脑，本文方法计算100个城市的所有167个Voronoi顶点的整个过程耗时不超过20s(包括测地变换过程在内)，与相同条件下计算单个椭球面三角形外心速度相当，说明算法效率基本不受对象个数的影响。作为一种面向椭球面空间且内容无关的算法，本文方法更适合于海量空间数据的高精度快速计算。

图 7 大西洋两岸100个城市及其Voronoi顶点 Fig. 7 The hundred cities on both sides of the Atlantic and their Voronoi vertexes

图选项

综上所述，本文方法精度与外心到对应椭球面三角形各顶点的大地线距离的长短呈负相关，与栅格分辨率的高低呈正相关，与椭球面三角形各顶点所处纬度无关，能够满足在全球范围内以可计算的栅格分辨率(全空间的栅格单元总数在10⁶～10⁸之间)快速求解外心精确大地坐标(定位误差小于0.001m)的需求。

本文方法采用矢栅紧密结合的数据结构，基于地图代数测地变换建立了椭球面空间的高精度距离场，在对分配场进行边界跟踪、配对、遍历和过滤操作之后，通过数值计算内插逼近的方法，最终计算出精确的大跨度椭球面三角形外心大地坐标。算法充分发挥了矢栅结合的优势，具有面向椭球面空间且内容无关的特性，能够同时兼顾计算精度和效率。该方法可以用来生成椭球面Voronoi图和构造椭球面Delaunay三角网，且已经成功应用于海域划界三接点的高精度快速求解。后续工作将研究如何将本文方法拓展到椭球面缓冲区和等比例线的求解问题^[27]，从而提高海域划界相关计算的精度和效率。