[算法系列]算法一 地理空间距离计算优化

Posted 2022-06-23 @SmartSi

1. 地理空间距离计算面临的挑战

打开美团app，不管是筛选团购还是筛选商家，默认的排序项都是“离我最近”或者“智能排序”（如下图所示）:

不管是“离我最近”还是“智能排序”，都涉及到计算用户位置与各个团购单子或者商家的距离（注：在智能排序中距离作为一个重要的参数参与排序打分）。以筛选商家为例，北京地区有5~6w个POI（本文将商家称之为POI），当用户进入商家页，请求北京全城+所有品类+离我最近/智能排序时，我们筛选服务需要计算一遍用户位置与北京全城所有POI的距离。

这种大量计算距离的场景十分消耗资源，从测试来看目前5w个点仅计算一遍距离就需要7ms，而到100w的时候就需要140多ms，随着数据的快速增长，筛选服务的性能将会非常堪忧。

如何优化距离的计算，进而提高计算速度、降低cpu使用率已经迫在眉睫。美团移动后台团购组在此方向上进行了些许探索，下文将分为4部分展开：

(1) 地理空间距离计算原理；

(2) Lucene使用的距离计算公式；

(3) 优化方案；

(4) 实际应用。

2. 地理空间距离计算原理

地理空间距离计算方法较多，目前我们使用的可以分为两类：

(1) 球面模型，这种模型将地球看成一个标准球体，球面上两点之间的最短距离即大圆弧长，这种方法使用较广，在我们服务端被广泛使用；

(2) 椭球模型，该模型最贴近真实地球，精度也最高，但计算较为复杂，目前客户端有在使用，但实际上我们的应用对精度的要求没有那么高。

下面着重介绍我们最常用的基于球面模型的地理空间距离计算公式，推导也只需要高中数学知识即可。

该模型将地球看成圆球，假设地球上有A(ja,wa)，B(jb,wb)两点

备注:

ja和jb分别是A和B的经度，wa和wb分别是A和B的纬度

弧长计算公式:

l = n（圆心角）× π（圆周率） × r（半径）/ 180
  = α(圆心角弧度数) × r（半径）

在半径是R的圆中，因为360°的圆心角所对的弧长就等于圆周长C=2πr，所以n°圆心角所对的弧长为

l = n° x 2πr / 360° = n°πr ÷ 180°

A和B两点的球面距离就是AB的弧长，AB弧长 = R * 角AOB圆心角弧度数。如何求出角AOB呢？可以先求AOB的最大边AB的长度，再根据余弦定律可以求夹角。

备注:

角AOB是A跟B的夹角，O是地球的球心，R是地球半径，约为6367000米

如何求出AB长度呢？

(1) 根据经纬度，以及地球半径R，将A、B两点的经纬度坐标转换成球体三维坐标:

(2) 根据A、B两点的三维坐标求AB长度:

(3) 根据余弦定理求出角AOB:

(4) AB弧长=R*角AOB:

3. Lucene使用的地理空间距离算法

目前美团团购app后台使用lucene来筛选团购单子和商家，lucene使用了spatial4j工具包来计算地理空间距离，而spatial4j提供了多种基于球面模型的地理空间距离的公式，其中一种就是上面我们推导的公式，称之为球面余弦公式。还有一种最常用的是Haversine公式，该公式是spatial4j计算距离的默认公式，本质上是球面余弦函数的一个变换，之前球面余弦公式中有cos(jb-ja)项，当系统的浮点运算精度不高时，在求算较近的两点间的距离时会有较大误差，Haversine方法进行了某种变换消除了cos(jb-ja)项，因此不存在短距离求算时对系统计算精度的过多顾虑的问题。

3.1 Haversine公式代码

public static double distHaversineRAD(double lat1, double lon1, double lat2, double lon2) 
  double hsinX = Math.sin((lon1 - lon2) * 0.5);
  double hsinY = Math.sin((lat1 - lat2) * 0.5);
  double h = hsinY * hsinY + (Math.cos(lat1) * Math.cos(lat2) * hsinX * hsinX);
  return 2 * Math.atan2(Math.sqrt(h), Math.sqrt(1 - h)) * 6367000;

3.2 Haversine公式性能

目前北京地区在线服务有5w个POI，计算一遍距离需要7ms。现在数据增长特别快，未来北京地区POI数目增大到100w时，我们筛选服务仅计算距离这一项就需要消耗144多ms，性能十分堪忧。（注：本文测试环境的处理器为2.9GHz Intel Core i7，内存为8GB 1600 MHz DDR3，操作系统为OS X10.8.3，实验在单线程环境下运行）:

POI数目	耗时（ms）
5w	7
10w	14
100w	144

4. 优化方案

通过抓栈我们发现消耗cpu较多的线程很多都在执行距离公式中的三角函数例如atan2等。因此我们的优化方向很直接---消减甚至消除三角函数。基于此方向我们尝试了多种方法，下文将一一介绍。

4.1 坐标转换方法

4.1.1 基本思路

之前POI保存的是经纬度(lon,lat)，我们的计算场景是计算用户位置与所有筛选出来的poi的距离，这里会涉及到大量三角函数计算。一种优化思路是POI数据不保存经纬度而保存球面模型下的三维坐标（x,y,z），映射方法如下：

x = Math.cos(lat) Math.cos(lon);
y = Math.cos(lat) Math.sin(lon);
z = Math.sin(lat);

那么当我们求夹角AOB时，只需要做一次点乘操作。比如求(lon1,lat1）和 (lon2,lat2）的夹角，只需要计算x1x2 + y1y2 + z1*z2， 这样避免了大量三角函数的计算。

在得到夹角之后，还需要执行arccos函数，将其转换成角度，AB弧长=角AOB*R（R是地球半径）。

此方法性能如下：

POI数目	耗时（ms）
5w	3
10w	8
100w	88

4.1.2 进一步优化

美团是本地生活服务，我们的业务场景是在一个城市范围内进行距离计算，因此夹角AOB往往会比较小，这个时候sinAOB约等于AOB，因此我们可以将 Math.acos(cosAOB)R 改为Math.sqrt(1 - cosAOBcosAOB)*R，从而完全避免使用三角函数，优化后性能如下。

POI数目	耗时（ms）
5w	0.2
10w	0.5
100w	4

4.2 简化距离计算公式方法

4.2.1 基本思路

我们的业务场景仅仅是在一个城市范围内进行距离计算，也就是说两个点之间的距离一般不会超过200多千米。由于范围小，可以认为经线和纬线是垂直的，如图所示，要求A（116.8，39,78）和B（116.9，39.68）两点的距离，我们可以先求出南北方向距离AM，然后求出东西方向距离BM，最后求矩形对角线距离，即sqrt(AMAM + BMBM)。

南北方向AM = R纬度差Math.PI/180.0；
东西方向BM = R经度差Cos<当地纬度数* Math.PI/180.0>

这种方式仅仅需要计算一次cos函数:

public static double distanceSimplify(double lat1, double lng1, double lat2, double lng2, double[] a) 
     double dx = lng1 - lng2; // 经度差值
     double dy = lat1 - lat2; // 纬度差值
     double b = (lat1 + lat2) / 2.0; // 平均纬度
     double Lx = toRadians(dx) * 6367000.0* Math.cos(toRadians(b)); // 东西距离
     double Ly = 6367000.0 * toRadians(dy); // 南北距离
     return Math.sqrt(Lx * Lx + Ly * Ly);  // 用平面的矩形对角距离公式计算总距离
    

我们对这个方法的有效性和性能进行验证。

(1) 有效性验证

我们首先检验这种简化是否能满足我们应用的精度，如果精度较差将不能用于实际生产环境。

我们的方法叫distanceSimplify，lucene的方法叫distHaversineRAD。下表是在不同尺度下两个方法的相差情况。

测试点对	distanceSimplify（米）	distHaversineRAD（米）	差别（米）
（39.941，116.45）（39.94， 116.451）	140.0285167225230	140.02851671981400	0.0
（39.96， 116.45）（39.94， 116.40）	4804.421262839180	4804.421153907680	0.0
（39.96， 116.45）（39.94， 117.30）	72444.81551882200	72444.54071519510	0.3
（39.26， 115.25）（41.04， 117.30）	263525.6167839860	263508.55921886700	17.1

可以看到两者在百米、千米尺度上几乎没有差别，在万米尺度上也仅有分米的差别，此外由于我们的业务是在一个城市范围内进行筛选排序，所以我们选择了北京左下角和右上角两点进行比较，两点相距有260多千米，两个方法差别17m。从精度上看该优化方法能满足我们应用需求。

(2) 性能验证

POI数目	耗时（ms）
5w	0.5
10w	1.1
100w	11

4.2.2 进一步优化

我们看到这里计算了一次cos这一三角函数，如果我们能消除此三角函数，那么将进一步提高计算效率。 如何消除cos三角函数呢？

采用多项式来拟合cos三角函数，这样不就可以将三角函数转换为加减乘除了嘛！

首先决定多项式的最高次数，次数为1和2显然都无法很好拟合cos函数，那么我们选择3先尝试吧，注：最高次数不是越多越好，次数越高会产生过拟合问题。

使用org.apache.commons.math3这一数学工具包来进行拟合。中国的纬度范围在10~60之间，即我们将此区间离散成Length份作为我们的训练集。

public static double[] trainPolyFit(int degree, int Length)
    PolynomialCurveFitter polynomialCurveFitter = PolynomialCurveFitter.create(degree);
    double minLat = 10.0; //中国最低纬度
    double maxLat = 60.0; //中国最高纬度
    double interv = (maxLat - minLat) / (double)Length;
    List<WeightedObservedPoint> weightedObservedPoints = new ArrayList<WeightedObservedPoint>();
    for(int i = 0; i < Length; i++) 
        WeightedObservedPoint weightedObservedPoint = new WeightedObservedPoint(1,  minLat + (double)i*interv, Math.cos(toRadians(x[i])));
        weightedObservedPoints.add(weightedObservedPoint);
    
    return polynomialCurveFitter.fit(weightedObservedPoints);



public static double distanceSimplifyMore(double lat1, double lng1, double lat2, double lng2, double[] a) 
     //1) 计算三个参数
     double dx = lng1 - lng2; // 经度差值
算法进阶系列1 空间搜索 GeoHash 算法

 写代码如何合理使用和优化我们的机器资源（CPU内存网络磁盘）
 空间权重矩阵的那些事-球面距离权重矩阵
 计算球面上任意两点间的球面距离（C++实现）
 地理空间距离计算优化
 地理空间距离计算及优化（根据两个点经纬度计算距离）