使用「Cheap Ruler」进行快速测地线逼近

Posted on Mar 27, 2022


距离计算是个说容易也容易,说难也难的事情。

  • 在对精度不是很敏感的场景上,用地球是半径 6371 公里的球的假设也足够用;
  • 精度要求高,计算量会显著增加,比如文森特公式;
  • Mapbox 开源的 cheapruler 库在绝大多数城市级别(纬度相近)地理计算上比传统的计算方法效率高得多,并且误差在可接受的范围上。 实现上非常讨巧,在部分地理计算是瓶颈的场景下能显著提升处理速度;

测地线计算是许多地理空间应用的重要组成部分。例如查找两个位置之间的距离、道路的长度、道路段的方向、多边形的面积或最接近道路的点。

在 Mapbox,我们依靠这些计算来处理来自我们 SDK 的遥测数据,这样我们就可以推断出未标示的街道、道路速度、转弯线路、转弯限制以及卫星地图更新的优先区域。

我们通常使用 Turf 进行这些计算ーー它是一个用于地理空间分析的简单、快速、模块化的 JavaScript 库。然而,遥测数据的数量与日俱增,处理这些数据的计算难度也越来越大。在这个过程中,我们可以挤出的每一点性能都是值得的。

为了优化处理过程,我构建了 Cheap Ruler,这是一个用于在城市尺度上进行超快测地线计算的 JavaScript 库。它不能替代 Turf ー它只实现了其方法的一小部分,有一定的局限性,并且只适用于对性能敏感的应用。但是在适合的地方,它可以比使用传统方法快 100 倍。让我们来探索它是如何工作的。

测量距离

为了得到一个平面上两点的距离,我们使用勾股定理:

这对地球上的位置不起作用,因为经纬度不能直接映射到距离单位。如果你站在其中一个极点附近,经度 360 度相当于几英里,而在赤道,同样的 360 度相当于 24900 英里。

其次,地球不是一个平坦的地方ーー它是球形的!幸运的是,数学家们发现了一种测量球体上距离的方法—- 首先使用球面余弦定律,然后使用在 Turf 和大多数地理应用程序中使用的半正矢公式进行改进。

不幸的是,地球甚至不是一个球体ーー它是一个扁圆形的椭球体。因此,半正矢公式可能导致 0.5% 的误差。为了解决这个问题,Thaddeus Vincenty 开发了一个非常复杂的公式,精确度可达 0.5 毫米,使其成为所有严肃的科学目的的终极测地公式。

这两个公式都依赖于许多计算费用高昂的三角计算,我们怎样才能避免它们呢?

欧几里德测地线近似

我们需要的大多数性能敏感测量都是在城市街区范围内进行的。在遥测技术方面,我们通常测量几十英里(而不是几千英里)的距离。对于像这样的小距离,我们可以假装地球是平的,并且避免了大多数复杂的计算,而不会造成明显的精度损失。

左: 子午线(经线) ,右: 平行线(纬线)。

为了估计两个纬度之间的距离(距离的“垂直”部分) ,我们可以按比例缩小到子午线的长度(大约 12,430 英里) ,因为它对于任何经度都大致相同:

$$ \partial y=12430 \frac{\left|lat_{1}-lat_{2}\right|}{180} $$

对于经度差异,纬线的长度(360 度圆)取决于纬度,从赤道长度(24,901 英里)开始,再用纬度的余弦进行标度:

$$ \partial x =24901 \frac{\left|lng_{1}-lng_{2}\right|}{360} \cos(lat) $$

然后,我们可以将这两个公式与欧氏勾股定理结合起来,用中间的纬度进行经度校正,得到两个位置之间距离的欧氏近似值:

$$ \begin{aligned} \partial y &=12430 \frac{\left|lat_{1}-lat_{2}\right|}{180} \newline \partial x &=24901 \frac{\left|lng_{1}-lng_{2}\right|}{360} \cos \left(\frac{lat_{1}+lat_{2}}{2}\right) \newline d &=\sqrt{\partial x^{2}+\partial y^{2}} \end{aligned} $$

由此得到的公式只有一个三角函数调用,比重三角函数的半正矢公式要快得多。

但还有一个考虑。假设我们正在计算一个特定的城市区域(例如在缩放等级为 14 的瓦片中)。坐标纬度的余弦值在一个很小的区域内变化不大,因此我们可以使用相同的经度校正倍增器进行所有这些计算,而不会造成明显的精度损失。所以我们只计算一个区域的余弦,然后完全避免三角函数的计算,导致了巨大的加速比。

对于使用 OpenStreetMap 公路的典型缩放 14 瓦片的简单基准,使用这种近似值计算公路长度比使用 Turf line-distance 计算公路长度要快 30 倍。

提高近似精度

从定义上看,这种近似方法不如传统方法精确,但精确度又有多少呢?

我写了一个小脚本来测量这一点,这里有一个表格,显示了上面描述的欧几里德近似和不同距离(单位是英里)和纬度的精确 Vincenty 公式之间的误差幅度:

正如你所看到的,用几千公里近似距离并不是一个好主意。但是如果你测量的范围在几百个之内,那么误差很小ーー不到 0.5% 。

现在让我们看看半正矢公式的误差(不依赖于距离) :

如果你只测量几百英里以内的距离,平坦地球近似值的误差只比适用于大多数应用场合的半正矢公式差略大。

但是更重要的是,有一种方法可以使它更加精确,同时保持同样的性能,即使用 FCC 规定的椭球投影到平面的公式:

$$ \begin{aligned} K_{1}&=111.13209-0.56605 \cos(2lat)+0.0012 \cos(4lat) \newline K_{2}&=111.41513 \cos (lat)-0.09455 \cos (3lat)+0.00012 \cos (5lat) \newline d&=\sqrt{\left(K_{1}\left(lat_{1}-lat_{2}\right)\right)^{2}+\left(K_{2}\left(lng_{1}-lng_{2}\right)\right)^{2}} \end{aligned} $$

这与之前描述的方法非常相似,只是我们使用了两个经纬度的系数。此外,我们可以使用切比雪夫方法来计算角度倍数的余弦值,将余弦值的计算从五个减少到一个。让我们来看看这些错误:

对于小于几百英里的数值,这种方法比半正矢公式法快得多,也更精确。

介绍 Cheap Ruler

我已经采用了上面的近似方法,并将其扩展到大多数在 Turf 中可用的测量方法。其结果是一个名为 cheap-ruler 的 JavaScript 模块。

这里有一个简短的列表,列出了它可以做什么,以及你在 Turf 上获得的加速(在 Node v6 上运行这些基准测试) :

  • distance: 快了 31 倍
  • bearing: 快了 3.6 倍
  • destination: 快了 7.2 倍
  • lineDistance: 快了 31 倍
  • area: 快了 3.4 倍
  • along: 快了 31 倍
  • pointOnLine: 快了 78 倍
  • lineSlice: 快了 60 倍

此外,我实现了几个实用函数,它们不直接映射 Turf 模块,但是可以非常快速地模拟通常由 Turf 调用组合完成的常见任务:

  • lineSliceAlong: 在给定的距离之间沿线划一条线,比 turf.lineSlice(turf.along(... 快 285 倍
  • bufferPoint: 给定一个缓冲区距离,在一个点周围创建一个边界框; 比用 turf.destination 创建两个对角线的边界框快 260 倍
  • bufferBBox: 快 260 倍
  • insideBBox: 比 turf.inside(turf.point(p), turf.bboxPolygon(bbox)) 快 19 倍

这个用法是这样的ー首先创建一个带有特定纬度(围绕计算所在的区域)和单位(例如英里)的标尺对象:

var ruler = cheapRuler(35.05, "miles");

如果你之前不知道纬度,你可以为每个特征创建一个单独的标尺(例如使用道路上第一个点的纬度)。

为了方便起见,您还可以根据瓦片坐标(yz)创建一个标尺:

var ruler = cheapRuler.fromTile(1567, 12, "kilometers");

然后你可以像使用 Turf 一样使用 ruler 对象,只不过你可以直接传递坐标,而不是把它们包装成 GeoJSON 特性:

var distance = ruler.distance([30.51, 50.32], [30.52, 50.312]);

避免特性包装器和输入验证使我们能够进一步压缩更多的性能。这个库用于性能至关重要的高级用途,而且我们并不针对普通用户,所以这样很好。

廉价尺子的另一个小优点是,它的代码行数不超过 200 行——非常适合在网页上使用,你希望它尽可能轻便,而且很容易移植到其他编程语言。

享受图书馆,如果你有任何建议或意见,请不要犹豫在 Twitter 上与我联系。谢谢!