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

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

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

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

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

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

(1) 地理空间距离计算原理&＃xff1b;

(2) Lucene使用的距离计算公式&＃xff1b;

(3) 优化方案&＃xff1b;

(4) 实际应用。

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

地理空间距离计算方法较多&＃xff0c;目前我们使用的可以分为两类&＃xff1a;

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

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

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

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

备注:

ja和jb分别是A和B的经度&＃xff0c;wa和wb分别是A和B的纬度

弧长计算公式:

l &＃61; n&＃xff08;圆心角&＃xff09;× π&＃xff08;圆周率&＃xff09; × r&＃xff08;半径&＃xff09;/ 180 &＃61; α(圆心角弧度数) × r&＃xff08;半径&＃xff09;

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

l &＃61; n° x 2πr / 360° &＃61; n°πr ÷ 180°

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

备注:

角AOB是A跟B的夹角&＃xff0c;O是地球的球心&＃xff0c;R是地球半径&＃xff0c;约为6367000米

如何求出AB长度呢&＃xff1f;

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

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

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

(4) AB弧长&＃61;R*角AOB:

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

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

3.1 Haversine公式代码

public static double distHaversineRAD(double lat1, double lon1, double lat2, double lon2) { double hsinX &＃61; Math.sin((lon1 - lon2) * 0.5); double hsinY &＃61; Math.sin((lat1 - lat2) * 0.5); double h &＃61; hsinY * hsinY &＃43; (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&＃xff0c;计算一遍距离需要7ms。现在数据增长特别快&＃xff0c;未来北京地区POI数目增大到100w时&＃xff0c;我们筛选服务仅计算距离这一项就需要消耗144多ms&＃xff0c;性能十分堪忧。&＃xff08;注&＃xff1a;本文测试环境的处理器为2.9GHz Intel Core i7&＃xff0c;内存为8GB 1600 MHz DDR3&＃xff0c;操作系统为OS X10.8.3&＃xff0c;实验在单线程环境下运行&＃xff09;:

POI数目	耗时&＃xff08;ms&＃xff09;
5w	7
10w	14
100w	144

4. 优化方案

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

4.1 坐标转换方法

4.1.1 基本思路

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

x &＃61; Math.cos(lat) Math.cos(lon); y &＃61; Math.cos(lat) Math.sin(lon); z &＃61; Math.sin(lat);

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

在得到夹角之后&＃xff0c;还需要执行arccos函数&＃xff0c;将其转换成角度&＃xff0c;AB弧长&＃61;角AOB*R&＃xff08;R是地球半径&＃xff09;。

此方法性能如下&＃xff1a;

POI数目	耗时&＃xff08;ms&＃xff09;
5w	3
10w	8
100w	88

4.1.2 进一步优化

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

POI数目	耗时&＃xff08;ms&＃xff09;
5w	0.2
10w	0.5
100w	4

4.2 简化距离计算公式方法

4.2.1 基本思路

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

南北方向AM &＃61; R纬度差Math.PI/180.0&＃xff1b; 东西方向BM &＃61; R经度差Cos<当地纬度数* Math.PI/180.0>

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

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

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

(1) 有效性验证

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

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

测试点对	distanceSimplify&＃xff08;米&＃xff09;	distHaversineRAD&＃xff08;米&＃xff09;	差别&＃xff08;米&＃xff09;
&＃xff08;39.941&＃xff0c;116.45&＃xff09;&＃xff08;39.94&＃xff0c; 116.451&＃xff09;	140.0285167225230	140.02851671981400	0.0
&＃xff08;39.96&＃xff0c; 116.45&＃xff09;&＃xff08;39.94&＃xff0c; 116.40&＃xff09;	4804.421262839180	4804.421153907680	0.0
&＃xff08;39.96&＃xff0c; 116.45&＃xff09;&＃xff08;39.94&＃xff0c; 117.30&＃xff09;	72444.81551882200	72444.54071519510	0.3
&＃xff08;39.26&＃xff0c; 115.25&＃xff09;&＃xff08;41.04&＃xff0c; 117.30&＃xff09;	263525.6167839860	263508.55921886700	17.1

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

(2) 性能验证

POI数目	耗时&＃xff08;ms&＃xff09;
5w	0.5
10w	1.1
100w	11

4.2.2 进一步优化

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

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

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

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

public static double[] trainPolyFit(int degree, int Length){ PolynomialCurveFitter polynomialCurveFitter &＃61; PolynomialCurveFitter.create(degree); double minLat &＃61; 10.0; //中国最低纬度 double maxLat &＃61; 60.0; //中国最高纬度 double interv &＃61; (maxLat - minLat) / (double)Length; List weightedObservedPoints &＃61; new ArrayList(); for(int i &＃61; 0; i WeightedObservedPoint weightedObservedPoint &＃61; new WeightedObservedPoint(1, minLat &＃43; (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 &＃61; lng1 - lng2; // 经度差值 double dy &＃61; lat1 - lat2; // 纬度差值 double b &＃61; (lat1 &＃43; lat2) / 2.0; // 平均纬度 //2) 计算东西方向距离和南北方向距离(单位&＃xff1a;米)&＃xff0c;东西距离采用三阶多项式 double Lx &＃61; (a[3] * b*b*b &＃43; a[2]* b*b &＃43;a[1] * b &＃43; a[0] ) * toRadians(dx) * 6367000.0; // 东西距离 double Ly &＃61; 6367000.0 * toRadians(dy); // 南北距离 //3) 用平面的矩形对角距离公式计算总距离 return Math.sqrt(Lx * Lx &＃43; Ly * Ly); }

我们对此优化方法进行有效性和性能验证。

(1) 有效性验证

我们的优化方法叫distanceSimplifyMore&＃xff0c;lucene的方法叫distHaversineRAD&＃xff0c;下表是在不同尺度下两个方法的相差情况。

测试点对	distanceSimplifyMore&＃xff08;米&＃xff09;	distHaversineRAD&＃xff08;米&＃xff09;	差别&＃xff08;米&＃xff09;
&＃xff08;39.941&＃xff0c;116.45&＃xff09;&＃xff08;39.94&＃xff0c; 116.451&＃xff09;	140.0242769266660	140.02851671981400	0.0
&＃xff08;39.96&＃xff0c; 116.45&＃xff09;&＃xff08;39.94&＃xff0c; 116.40&＃xff09;	4804.113098854450	4804.421153907680	0.3
&＃xff08;39.96&＃xff0c; 116.45&＃xff09;&＃xff08;39.94&＃xff0c; 117.30&＃xff09;	72438.90919479560	72444.54071519510	5.6
&＃xff08;39.26&＃xff0c; 115.25&＃xff09;&＃xff08;41.04&＃xff0c; 117.30&＃xff09;	263516.676171262	263508.55921886700	8.1

可以看到在百米尺度上两者几乎未有差别&＃xff0c;在千米尺度上仅有分米的区别&＃xff0c;在更高尺度上如72千米仅有5.6m米别&＃xff0c;在264千米也仅有8.1米区别&＃xff0c;因此该优化方法的精度能满足我们的应用需求。

(2) 性能验证

POI数目	耗时&＃xff08;ms&＃xff09;
5w	0.1
10w	0.3
100w	4

5. 实际应用

坐标转换方法和简化距离公式方法性能都非常高&＃xff0c;相比lucene使用的Haversine算法大大提高了计算效率&＃xff0c;然而坐标转换方法存在一些缺点&＃xff1a;

(1) 坐标转换后的数据不能被直接用于空间索引。lucene可以直接对经纬度进行geohash空间索引&＃xff0c;而通过空间转换变成三维数据后不能直接使用。我们的应用有附近范围筛选功能&＃xff08;例如附近5km的团购单子&＃xff09;&＃xff0c;通过geohash空间索引可以提高范围筛选的效率&＃xff1b;

(2) 坐标转换方法增大内存开销。我们会将坐标写入倒排索引中&＃xff0c;之前坐标是2列&＃xff08;经度和纬度&＃xff09;&＃xff0c;现在变成3列&＃xff08;x,y,z&＃xff09;&＃xff0c;在使用中我们往往会将这数据放入到cache中&＃xff0c;因此会增大内存开销&＃xff1b;

(3) 坐标转换方法增大建索引开销。此方法本质上是将计算从查询阶段放至到索引阶段&＃xff0c;因此提高了建索引的开销。

基于上述原因我们在实际应用中采用简化距离公式方法&＃xff08;通过三次多项式来拟合cos三角函数&＃xff09;&＃xff0c;此方法在团购筛选和商家筛选的距离排序、智能排序中已经开始使用&＃xff0c;与之前相比&＃xff0c;筛选团购时北京全城美食品类距离排序响应时间从40ms下降为20ms。