地理空间距离计算优化【转】

原理

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

球面模型，这种模型将地球看成一个标准球体
球面上两点之间的最短距离即大圆弧长，这种方法使用较广，在服务端被广泛使用。
椭球模型，该模型最贴近真实地球，精度也最高，但计算较为复杂
目前客户端有在使用，但实际上应用对精度的要求没有那么高。

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

该模型将地球看成圆球，假设地球上有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长度呢？

根据经纬度，以及地球半径R，将A、B两点的经纬度坐标转换成球体三维坐标：
$\begin{aligned} &A\left\{\begin{array}{l} X a=R \cos (w a) \cos (j a) \\ Y a=R \cos (w a) \sin (j a) \\ Z a=R \sin (w a) \end{array}\right\}\\ &B\left\{\begin{array}{l} X b=R \cos (w b) \cos (j b) \\ Y b=\operatorname{Rcos}(w b) \sin (j b) \\ Z b=\operatorname{Rsin}(w b) \end{array}\right\} \end{aligned}$
根据A、B两点的三维坐标求AB长度：
$\begin{array}{l} A B^{2}=(X a-X b)^{2}+(Y a-Y b)^{2}+(Z a-Z b)^{2} \\ =\ldots=2 R^{2}(1-\cos (w a) \cos (w b) \cos (j b-j a)-\sin (w a) \sin (w b)) \end{array}$
根据余弦定理求出角AOB：
$\begin{array}{l} A B^{2}=A O^{2}+B O^{2}-2 A O * B O * \cos (\triangle A O B) \\ \cos (\triangle A O B)=\frac{A B^{2}-A O^{2}-B O^{2}}{-2 A O * B O}=\frac{A B^{2}-2 R^{2}}{-2 R^{2}}=1-\frac{A B^{2}}{2 R^{2}} \end{array}$
AB弧长=R*角AOB：
$\hat{A B}=\operatorname{Rarccos}(\cos (w a) \cos (w b) \cos (j b-j a)+\sin (w a) \sin (w b))$

Lucene

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

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

优化

计算多消耗在三角函数中，因此优化目标是消除三角函数

坐标转换

一种优化思路是POI数据不保存经纬度而保存球面模型下的三维坐标（x, y, z）。

那么当求夹角AOB时，只需要做一次点乘操作。比如求(lon1,lat1）和 (lon2,lat2）的夹角，
只需要计算 $x1*x2 + y1*y2 + z1*z2$ ，这样避免了大量三角函数的计算。
在得到夹角之后，还需要执行arccos函数，将其转换成角度，AB弧长=角AOB*R（R是地球半径）。

进一步优化：
当业务场景是在一个城市范围内进行距离计算时，夹角AOB往往会比较小，这个时候sinAOB约等于AOB，
因此可以将 Math.acos(cosAOB) * R 改为 Math.sqrt(1 - cosAOB * cosAOB) * R，从而完全避免使用三角函数。

简化距离

基本思路

当业务场景仅仅是在一个城市范围内进行距离计算时，两个点之间的距离一般不会超过200多千米。
由于范围小，可以认为经线和纬线是垂直的，这种方式仅仅需要计算一次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);  // 用平面的矩形对角距离公式计算总距离
    }
}

有效性验证

测试点对	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。从精度上看该优化方法能满足应用需求。

进一步优化

采用多项式来拟合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++) {
        double x = minLat + (double) i * interv;
        WeightedObservedPoint weightedObservedPoint = new WeightedObservedPoint(1, x, Math.cos(toRadians(x)));
        weightedObservedPoints.add(weightedObservedPoint);
    }
    return polynomialCurveFitter.fit(weightedObservedPoints);
}

得到多项式系数
0.9972739367451966, 4.323852126263266E-4, -1.7522667541578619E-4, 4.966195349399909E-7

$y = 0.0000004966195 x^{3}-0.0001752266754 x^{2}+0.0004323852126 x+0.9972739367452$

函数图形：
fit

新方法加入多项式系数数组参数fit：

public static double distanceSimplifyMore(double lat1, double lng1, double lat2, double lng2, double[] fit) {
    // 1) 计算三个参数
    double dx = lng1 - lng2;            // 经度差值
    double dy = lat1 - lat2;            // 纬度差值
    double b = (lat1 + lat2) / 2.0;     // 平均纬度
    // 2) 计算东西方向距离和南北方向距离(单位：米)，东西距离采用三阶多项式
    double Lx = (fit[3] * b * b * b + fit[2] * b * b + fit[1] * b + fit[0]) * toRadians(dx) * 6367000.0; // 东西距离
    double Ly = 6367000.0 * toRadians(dy); // 南北距离
    // 3) 用平面的矩形对角距离公式计算总距离
    return Math.sqrt(Lx * Lx + Ly * Ly);
}

有效性验证

测试点对	distanceSimplifyMore（米）	distHaversineRAD（米）	差别（米）
(39.941,116.45) (39.94, 116.451)	140.0242769266660	140.02851671981400	0.0
(39.96, 116.45) (39.94, 116.40)	4804.113098854450	4804.421153907680	0.3
(39.96, 116.45) (39.94, 117.30)	72438.90919479560	72444.54071519510	5.6
(39.26, 115.25) (41.04, 117.30)	263516.676171262	263508.55921886700	8.1

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

Benchmark

benchmark

Intel(R) Core(TM)i5-6200 CPU @ 2.30GHz - 2.4GHz
16GB 1600MHz DDR3
Windows10 企业版 (2018)

distanceSimplifyMore 1000W 随机经纬测试耗时大约300ms
实际测试中，distanceSimplifyMore 耗时太短，为做参考在方法中加入生成一次随机数。

实际应用

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

坐标转换后的数据不能被直接用于空间索引。lucene可以直接对经纬度进行geohash空间索引，而通过空间转换变成三维数据后不能直接使用。
应用有附近范围筛选功能（例如附近5km的团购单子），通过geohash空间索引可以提高范围筛选的效率；
坐标转换方法增大内存开销。将坐标写入倒排索引中，之前坐标是2列（经度和纬度），
现在变成3列（x,y,z），在使用中往往会将这数据放入到cache中，因此会增大内存开销；
坐标转换方法增大建索引开销。此方法本质上是将计算从查询阶段放至到索引阶段，因此提高了建索引的开销。

基于上述原因在实际应用中采用简化距离公式方法（通过三次多项式来拟合cos三角函数），
与之前相比，筛选团购时北京全城美食品类距离排序响应时间从40ms下降为20ms。

地理空间距离计算优化