地理空间距离计算及优化（依据两个点经纬度计算距离）

时间：2017-08-18 14:41:21 阅读：5563 评论：0 收藏：0 [点我收藏+]

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

打开美团app。无论是筛选团购还是筛选商家，默认的排序项都是“离我近期”或者“智能排序”（例如以下图所看到的）。

技术分享

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

这样的大量计算距离的场景十分消耗资源，从測试来看眼下5w个点仅计算一遍距离就须要7ms，而到100w的时候就须要140多ms，随着数据的高速增长，筛选服务的性能将会很堪忧。

怎样优化距离的计算。进而提高计算速度、减少cpu使用率已经迫在眉睫。美团移动后台团购组在此方向上进行了些许探索，下文将分为4部分展开：1）地理空间距离计算原理；2）Lucene使用的距离计算公式；3）优化方案。4）实际应用。

2 地理空间距离计算原理

地理空间距离计算方法较多，眼下我们使用的能够分为两类：1）球面模型，这样的模型将地球看成一个标准球体，球面上两点之间的最短距离即大圆弧长。这样的方法使用较广，在我们服务端被广泛使用。2）椭球模型，该模型最贴近真实地球。精度也最高，但计算较为复杂，眼下client有在使用。但实际上我们的应用对精度的要求没有那么高。

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

技术分享

经纬度距离计算示意图

该模型将地球看成圆球，如果地球上有A(ja,wa)，B(jb,wb)两点（注：ja和jb各自是A和B的经度。wa和wb各自是A和B的纬度），A和B两点的球面距离就是AB的弧长，AB弧长=R*角AOB（注：角AOB是A跟B的夹角，O是地球的球心，R是地球半径，约为6367000米）。怎样求出角AOB呢？能够先求AOB的最大边AB的长度，再依据余弦定律能够求夹角。

怎样求出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)项，因此不存在短距离求算时对系统计算精度的过多顾虑的问题。

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;
    }

2）Haversine公式性能

依据2个经纬度点，计算这2个经纬度点之间的距离(通过经度纬度得到距离)

球面上随意两点之间的距离计算公式能够參考维基百科上的下述文章。

值得一提的是，维基百科推荐使用Haversine公式。理由是Great-circle distance公式用到了大量余弦函数，而两点间距离非常短时（比方地球表面上相距几百米的两点）。余弦函数会得出0.999...的结果。会导致较大的舍入误差。而Haversine公式採用了正弦函数，即使距离非常小，也能保持足够的有效数字。曾经採用三角函数表计算时的确会有这个问题，但经过实际验证。採用计算机来计算时。两个公式的差别不大。稳妥起见。这里还是採用Haversine公式。

技术分享

当中

技术分享

R为地球半径。可取平均值 6371km；
φ₁, φ₂ 表示两点的纬度；
Δλ 表示两点经度的差值。

依据2个经纬度坐标。距离计算函数

以下就是计算球面间两点(lat1, lon1) - (lat2, lon2)之间距离的函数。

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace HarvenSin
{
    class Program
    {
        /// <summary>
        /// 依据经纬度。计算2个点之间的距离。
        /// </summary>
        /// <param name="args"></param>
        static void Main(string[] args)
        {
            //39.94607,116.32793  31.24063,121.42575
            Console.WriteLine(Distance(39.94607, 116.32793, 31.24063, 121.42575));

        }


        public static double HaverSin(double theta)
        {
            var v = Math.Sin(theta / 2);
            return v * v;
        }


        static double EARTH_RADIUS = 6371.0;//km 地球半径 平均值。千米

        /// <summary>
        /// 给定的经度1。纬度1；经度2，纬度2. 计算2个经纬度之间的距离。
        /// </summary>
        /// <param name="lat1">经度1</param>
        /// <param name="lon1">纬度1</param>
        /// <param name="lat2">经度2</param>
        /// <param name="lon2">纬度2</param>
        /// <returns>距离（公里、千米）</returns>
        public static double Distance(double lat1,double lon1, double lat2,double lon2)
        {
            //用haversine公式计算球面两点间的距离。

            //经纬度转换成弧度
            lat1 = ConvertDegreesToRadians(lat1);
            lon1 = ConvertDegreesToRadians(lon1);
            lat2 = ConvertDegreesToRadians(lat2);
            lon2 = ConvertDegreesToRadians(lon2);

            //差值
            var vLon = Math.Abs(lon1 - lon2);
            var vLat = Math.Abs(lat1 - lat2);

            //h is the great circle distance in radians, great circle就是一个球体上的切面，它的圆心即是球心的一个周长最大的圆。
            var h = HaverSin(vLat) + Math.Cos(lat1) * Math.Cos(lat2) * HaverSin(vLon);

            var distance = 2 * EARTH_RADIUS * Math.Asin(Math.Sqrt(h));

            return distance;
        }

        /// <summary>
        /// 将角度换算为弧度。
        /// </summary>
        /// <param name="degrees">角度</param>
        /// <returns>弧度</returns>
        public static double ConvertDegreesToRadians(double degrees)
        {
            return degrees * Math.PI / 180;
        }

        public static double ConvertRadiansToDegrees(double radian)
        {
            return radian * 180.0 / Math.PI;
        }

    }
}

公式来历：

VERSINE(F)=1-cos(F)

技术分享

Haversine名字来历是Ha-VERSINE，即Half-Versine 。表示sin的一半的意思。

hav(A) = (1-cos(A))/2 = sin(A/2)* sin(A/2)

推倒过程：

例如以下一个半径为1 的圆，O是圆心。A、B是弦(chord)。

角度AOB=theta。

则角度AOC=theta/2。OC是垂直于AB的垂线(perpendicular)。AC长度是sin(theta/2)，AB长度是2*sin(theta/2)。

技术分享

(图1)

例如以下地球图所看到的，如果半径R为1。O是球心，A (lat1,lon1) 和 B (lat2,lon2) 是我们感兴趣的2个点。2跟经度线 lon1。lon2相交于北极(north pole)N。EF所在的线是赤道(equator)。ACBD是平面上的等腰梯形的四个顶点(vertice)。

AC和DB的弦（直线）在图上没有画出。CD的位置是：C (lat2,lon1) and D (lat1,lon2)。

角度AOC是A点与C点的纬度差 dlat。角度EOF是经度E点和经度F点的差dlon。

技术分享

(图2)

弦AC的长度。參照图1的方式，那么是AC=2*sin(dlat/2)，弦BD也是一样的长度。

E、F 2个点是赤道上的2个点，它们的纬度是0。

EF的距离是EF=2*sin(dlon/2)

A、D2个点所在的纬度是lat1。

AD所在纬度的圆平面的半径是cos(lat1)。从A作一条垂线(perpendicular)到OE为AG，AO是球半径，则OG=cos(lat1)，即A、D所在纬度圆圈的半径（AO`）。

这时候，AD的弦长AD= 2*sin(dlon/2)*cos(lat1)，类似的能够推出CB的长度= CB=2*sin(dlon/2)*cos(lat2)

以下看一下怎样求AB的长度。回到平面等腰梯形。例如以下图：

技术分享

(图3)

AH是到CB的垂线(perpendicular)，CH= (CB-AD)/2。

依据勾股定理(Pythagorean theorem): 【^2表示2的平方】

AH^2 = AC^2 - CH^2

= AC^2 - (CB-AD)^2/4

HB 的长度是HB=AD+CH = AD+(CB-AD)/2 = (CB+AD)/2。依据勾股定理得到：

AB^2 = AH^2 + HB^2

= AC^2 - (CB-AD)^2/4 + (CB+AD)^2/4

= AC^2 + CB*AD

依据前面球面上的求经纬距离的方式。我们已经得到 AC、AD和CB的长度。代入公式得到：

AB^2 = 4*(sin^2(dlat/2) + 4*cos(lat1)*cos(lat2)*sin^2(dlon/2))

如果中间值h 是AB长度一半的平方，例如以下

h = (AB/2)^2

= (sin^2(dlat/2)) + cos(lat1) * cos(lat2) * sin^2(dlon/2)

（请參看代码里的h）

最后一步。是求得代表AB长度的角度AOB。

參照图1的方式。我们能够知道

技术分享

(图4)

设AC= 技术分享，依据勾股定理(Pythagorean theorem)得到：

OC= 技术分享 = sqrt(OA^2 - AC^2)

= 技术分享 = sqrt(1-a) // sqrt表示开根号

假设设c是角AOB的度数值。

tan(<AOC) = tan(c)= AC/OC = sqrt(a)/sqrt(1-a)

则：

c = 2 * arctan(sqrt(a)/sqrt(1-a))，

最后的AB真实距离，把地球半径带上就能够了。

distance = 2 * EARTH_RADIUS * c。

2)第二种方法：

SQL Server本身是支持空间数据索引的(Spatial Indexing)，具有空间数据计算能力。

他是通过一个扩展DLL Microsoft.SqlServer.Types.dll 来实现这些功能的。

这是一个托管DLL，那意味着.NET C# asp.net 也能够使用些功能。

比如通过 reference 引用： Microsoft.SqlServer.Types.dll 这个dll。

var a = SqlGeography.Point(22.54587746 , 114.12873077, 4326); //上海的某个点
var b = SqlGeography.Point(23, 115, 4326); //上海的某个点,4236代表WGS84这样的坐标參照系统。
Console.WriteLine(a.STDistance(b)); //距离

这个算出来的距离，与上面使用haversine公式算出的距离，误差在几米之内。

參考资料：

http://mathforum.org/library/drmath/view/51879.html

http://blog.charlee.li/location-search/

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

技术分享

4.优化方案

通过抓栈我们发现消耗cpu较多的线程非常多都在运行距离公式中的三角函数比如atan2等。因此我们的优化方向非常直接---消减甚至消除三角函数。

基于此方向我们尝试了多种方法。下文将一一介绍。

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是地球半径）。

此方法性能例如以下：

技术分享

2）进一步优化

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

技术分享

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<当地纬度数* 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.1）有效性验证

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

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

技术分享

能够看到两者在百米、千米尺度上差点儿没有区别，在万米尺度上也仅有分米的区别，此外因为我们的业务是在一个城市范围内进行筛选排序，所以我们选择了北京左下角和右上角两点进行比較，两点相距有260多千米，两个方法区别17m。

从精度上看该优化方法能满足我们应用需求。

1.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 weightedObservedPoints = new ArrayList();
    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; // 经度差值
     double dy = lat1 - lat2; // 纬度差值
     double b = (lat1 + lat2) / 2.0; // 平均纬度
     //2) 计算东西方向距离和南北方向距离(单位：米)。东西距离採用三阶多项式
     double Lx = (a[3] * b*b*b  + a[2]* b*b  +a[1] * b + a[0] ) * toRadians(dx) * 6367000.0; // 东西距离
     double Ly = 6367000.0 * toRadians(dy); // 南北距离
     //3) 用平面的矩形对角距离公式计算总距离
     return Math.sqrt(Lx * Lx + Ly * Ly);
}