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