地理空间距离计算

1、球面半正矢公式

来源：维基百科 http://en.wikipedia.org/wiki/Haversine_formula

1.1、代码

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

1.2、优化

将经纬度坐标转三维坐标减少三角函数计算

经纬度坐标转三维直角坐标系坐标（将地球看成圆）：

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

sinAOB≈AOB完全避免三角函数计算

2、其他计算方法

业务场景仅仅是在一个城市范围内进行距离计算，也就是说两个点之间的距离一般不会超过200多千米。由于范围小，可以认为经线和纬线是垂直的。由勾股定理，点间距离可由两点水平距离和垂直距离得到

2.1 代码

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);  // 用平面的矩形对角距离公式计算总距离
    }
}

2.2 优化

采用多项式来拟合cos三角函数来去掉上面的三角函数cos。

使用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; // 经度差值
     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);
}

3、参考资料

1、维基百科 球面距离公式

2、美团点评技术团队 地理空间距离计算优化