根据经纬度对球面距离的求解

作者: 阿杜S考特 | 来源:发表于2019-07-25 01:10 被阅读0次

摘要

本片文章主要是推导球面距离公式的过程，需要您了解一些高中数学的三角函数的基本公式，简单的向量运算，如果你还记得的话，或者有兴趣的话可以接着往下看；事情的起因是在项目的代码中看到有一段根据手机定位的经纬度数据计算两个地理位置的距离。

export function calcDistance(point1, point2) {
  if (!point1 || !point2) return null;
  const EARTH_RADIUS = 6378.137;
  const { longitude: lng1, latitude: lat1 } = point1;
  const { longitude: lng2, latitude: lat2 } = point2;
  const radLat1 = (lat1 * Math.PI) / 180.0;
  const radLat2 = (lat2 * Math.PI) / 180.0;
  const radDiff = radLat1 - radLat2;
  const b = (lng1 * Math.PI) / 180.0 - (lng2 * Math.PI) / 180.0;
  const distance =
    2 *
    Math.asin(
      Math.sqrt(
        Math.pow(Math.sin(radDiff / 2), 2) +
          Math.cos(radLat1) * Math.cos(radLat2) * Math.pow(Math.sin(b / 2), 2)
      )
    );
  return distance * EARTH_RADIUS;
}

下面让我来解释一下上面这段代码的意思，

球面坐标系

假设地球半径为 $R$ ，经度longitude用希腊字母 $\lambda$ 表示；纬度latitude用希腊字母 $\varphi$ 表示，那么P1的坐标就是 $(\lambda_1,\varphi_1)$ ；而P2的坐标就是 $(\lambda_2,\varphi_2)$ ，那么代码的核心就是下面这个公式：

$|P_1P_2|=2R\sqrt{\sin^2\frac{\Delta\varphi}{2}+\cos\varphi_1cos\varphi_2\sin^2\frac{\Delta\lambda}{2}}$

其中
$\Delta\varphi=\varphi_1-\varphi_2，\Delta\lambda=\lambda_1-\lambda_2$

这个公式就是传说中的的半正矢公式Haversine formula，然后我就很好奇为啥这个公式是这样的，纠结这个公式是怎么来的，光是根据球面坐标系和这个公式，并不能理解，它并不像欧几里得空间上的距离那么直观。

下面我们根据空间直角坐标系来定义一下球面坐标系，球面上的任意一点 $P(x,y,z)$ 用经纬度来转换可以用下面的公式来表示:

$\begin{equation} \left\{ \begin{array}{lr} x=R\cos{\varphi}\cos{\lambda}\\ y=R\cos{\varphi}\sin{\lambda}\\ z=R\sin{\varphi} \end{array} \right. \end{equation}$

方法一：从弦长求大圆距离

用欧几里得空间距离求解，对于球面上的两点 $P_1(x_1,y_1,z_1)$ 和 $P_2(x_2,y_2,z_2)$ 两点之间的距离可以表示为：

$|P_1P_2|=\sqrt{(x_1-x_2)^2+(y_1-y_2)^2+(z_1-z_2)^2}$

那么转换成球面坐标以后就是将球面公式(1)代入到欧式距离公式中，则就有：

$|P_1P_2|^2=R^2[(\cos{\varphi_1}cos{\lambda_1}-\cos{\varphi_2}\cos{\lambda_2})^2+(\cos{\varphi_1}\sin{\lambda_1}-\cos{\varphi_2}\sin{\lambda_2})^2\\+(\sin{\varphi_1}-\sin{\varphi_2})^2]$

等价于

$|P_1P_2|^2=R^2[\cos^2{\varphi_1}\cos^2{\lambda_1}+\cos^2{\varphi_2}\cos^2{\lambda_2}-2\cos{\varphi_1}\cos{\varphi_2}\cos{\lambda_1}\cos{\lambda_2}\\ +\cos^2{\varphi_1}\sin^2{\lambda_1}+\cos^2{\varphi_2}\sin^2{\lambda_2}-2\cos{\varphi_1}\cos{\varphi_2}\sin{\lambda_1}\sin{\lambda_2}\\ +\sin^2{\varphi_1}+\sin^2{\varphi_2}-2\sin{\varphi_1}\sin{\varphi_2}]$

等价于

$|P_1P_2|^2=R^2[\cos^2{\varphi_1}+\cos^2{\varphi_2}+\sin^2{\varphi_1}+\sin^2{\varphi_2}\\ -2\cos{\varphi_1}\cos{\varphi_2}(\cos{\lambda_1}\cos{\lambda_2}+\sin{\lambda_1}\sin{\lambda_2})\\ -2\sin{\varphi_1}\sin{\varphi_2}]$

等价于

$|P_1P_2|^2=R^2[2-2\cos{\varphi_1}\cos{\varphi_2}\cos{(\lambda_1-\lambda_2)}-2\sin{\varphi_1}\sin{\varphi_2}]$

等价于

$|P_1P_2|^2=2R^2(1-\sin{\varphi_1}\sin{\varphi_2}-\cos{\varphi_1}\cos{\varphi_2}\cos{\Delta\lambda})$

即

$|P_1P_2|=\sqrt{2}R\sqrt{1-\sin{\varphi_1}\sin{\varphi_2}-\cos{\varphi_1}\cos{\varphi_2}\cos{\Delta\lambda}}$

接下来我们把问题转化到平面几何上，参照下图， $|P_1P_2|$ 的长度我们可以通过上门的公式求出，现在我们需要求出 $\stackrel{\frown}{|P_1P_2|}$ ，现在已知 $|OP_1|=|OP_2|=R$ ，假设 $\angle P_1OP_2=\delta$ ，那么 $\stackrel{\frown}{|P_1P_2|}=\delta R$

球体在圆上的切片

由上图显然:

$\sin{\frac{\delta}{2}}=\frac{|P_1Q|}{R}=\frac{|P_2Q|}{R}=\frac{|P_1P_2|}{2R}$

这就等价于

$\delta=2\arcsin{\frac{|P_1P_2|}{2R}}$

所以

$\delta=2\arcsin{\sqrt{\frac{{1-\sin{\varphi_1}\sin{\varphi_2}-\cos{\varphi_1}\cos{\varphi_2}\cos{\Delta\lambda}}}{2}}}$

最终得到

$\stackrel{\frown}{|P_1P_2|}=2R\arcsin{\sqrt{\frac{{1-\sin{\varphi_1}\sin{\varphi_2}-\cos{\varphi_1}\cos{\varphi_2}\cos{\Delta\lambda}}}{2}}}$

以上是最大圆距离(Creat-circle distance)的一种表现形式。

方法二：用向量的余弦夹角公式求解大圆距离

用球面坐标表示 $P_1(x_1,y_1,z_1)$ 和 $P_2(x_2,y_2,z_2)$ 的坐标，可以得到

$P_1(R\cos{\varphi_1}\cos{\lambda_1},R\cos{\varphi_1}\sin{\lambda_1},R\sin{\varphi_1})$

$P_2(R\cos{\varphi_2}\cos{\lambda_2},R\cos{\varphi_2}\sin{\lambda_2},R\sin{\varphi_2})$

那么向量 $\overrightarrow{OP_1}$ 和向量 $\overrightarrow{OP_2}$ 的夹角 $\delta$ 的余弦：

$\cos{\angle P_1OP_2}=\cos{\delta}=\frac{\overrightarrow{OP_1}\cdot\overrightarrow{OP_2}}{|\overrightarrow{OP_1}||\overrightarrow{OP_2}|}$

那么

$\cos{\delta}=\frac{x_1*x_2+y_1*y_2+z_1*z_2}{R^2}$

$\cos{\delta}=\cos{\varphi_1}\cos{\varphi_2}\cos{\lambda_1}\cos{\lambda_2}+\cos{\varphi_1}\cos{\varphi_2}\sin{\lambda_1}\sin{\lambda_2}+\sin{\varphi_1}\sin{\varphi_2}$

$\cos{\delta}=\cos{\varphi_1}\cos{\varphi_2}[\cos{\lambda_1}\cos{\lambda_2}+\sin{\lambda_1}\sin{\lambda_2}]+\sin{\varphi_1}\sin{\varphi_2}$

$\cos{\delta}=\cos{\varphi_1}\cos{\varphi_2}\cos{(\lambda_1-\lambda_2)}+\sin{\varphi_1}\sin{\varphi_2}$

$\cos{\delta}=\cos{\varphi_1}\cos{\varphi_2}\cos{\Delta \lambda}+\sin{\varphi_1}\sin{\varphi_2}$

由此可知

$\delta=\arccos{(\cos{\varphi_1}\cos{\varphi_2}\cos{\Delta \lambda}+\sin{\varphi_1}\sin{\varphi_2})}$

最终得到最大圆距离(Creat-circle distance)的另一种表现形式，一般网上找到的Creat-circle distance指的就是这个公式

$\stackrel{\frown}{|P_1P_2|}=R\arccos{(\cos{\varphi_1}\cos{\varphi_2}\cos{\Delta \lambda}+\sin{\varphi_1}\sin{\varphi_2})}$

注：这个公式的余弦函数在地球上距离较近的两点，有计算精度丢失的问题，因此在此基础上我们可以利用余弦定理等价变换成另一种易于计算机较精确计算的公式，半正矢函数公式Haversine公式

方法三、用半正矢公式定义的距离

球体在圆上的切片

由余弦定理可知

$|P_1P_2|^2=|OP_1|^2+|OP_2|^2-2|OP_1|\cdot|OP_2|\cdot \cos{\delta}$

将 $\cos{\delta}=\cos{\varphi_1}\cos{\varphi_2}\cos{(\lambda_1-\lambda_2)}+\sin{\varphi_1}\sin{\varphi_2}$ 代入上式得到:

$|P_1P_2|^2=2R^2-2R^2(\cos{\varphi_1}\cos{\varphi_2}\cos{\Delta\lambda}+\sin{\varphi_1}\sin{\varphi_2})$

$|P_1P_2|^2=2R^2-2R^2[\cos{\varphi_1}\cos{\varphi_2}(1-2\sin^2{\frac{\Delta\lambda}{2}})+\sin{\varphi_1}\sin{\varphi_2}]$

$|P_1P_2|^2=2R^2-2R^2[\cos{\varphi_1}\cos{\varphi_2}-2\cos{\varphi_1}\cos{\varphi_2}\sin^2{\frac{\Delta\lambda}{2}}+\sin{\varphi_1}\sin{\varphi_2}]$

$|P_1P_2|^2=2R^2-2R^2[\cos{\varphi_1}\cos{\varphi_2}+\sin{\varphi_1}\sin{\varphi_2}-2\cos{\varphi_1}\cos{\varphi_2}\sin^2{\frac{\Delta\lambda}{2}}]$

$|P_1P_2|^2=2R^2-2R^2[\cos{(\varphi_1-\varphi_2)}-2\cos{\varphi_1}\cos{\varphi_2}\sin^2{\frac{\Delta\lambda}{2}}]$

$|P_1P_2|^2=2R^2-2R^2[\cos{\Delta\varphi}-2\cos{\varphi_1}\cos{\varphi_2}\sin^2{\frac{\Delta\lambda}{2}}]$

$|P_1P_2|^2=2R^2-2R^2\cos{\Delta\varphi}+4R^2\cos{\varphi_1}\cos{\varphi_2}\sin^2{\frac{\Delta\lambda}{2}}$

$|P_1P_2|^2=2R^2(1-\cos{\Delta\varphi})+4R^2\cos{\varphi_1}\cos{\varphi_2}\sin^2{\frac{\Delta\lambda}{2}}$

$|P_1P_2|^2=4R^2\sin^2{\frac{\Delta\varphi}{2}}+4R^2\cos{\varphi_1}\cos{\varphi_2}\sin^2{\frac{\Delta\lambda}{2}}$

由正弦定义可知：

$\sin{\frac{\delta}{2}}=\frac{|P_1P_2|}{2R}$

这等价于

$\sin^2{\frac{\delta}{2}}=\frac{|P_1P_2|^2}{4R^2}$

即

$\sin^2{\frac{\delta}{2}}=\frac{|P_1P_2|^2}{4R^2}=\sin^2{\frac{\Delta\varphi}{2}}+\cos{\varphi_1}\cos{\varphi_2}\sin^2{\frac{\Delta\lambda}{2}}$

$\sin{\frac{\delta}{2}}=\sqrt{\sin^2{\frac{\Delta\varphi}{2}}+\cos{\varphi_1}\cos{\varphi_2}\sin^2{\frac{\Delta\lambda}{2}}}$

所以

$\delta=2\arcsin{\sqrt{\sin^2{\frac{\Delta\varphi}{2}}+\cos{\varphi_1}\cos{\varphi_2}\sin^2{\frac{\Delta\lambda}{2}}}}$

最终得到

$\stackrel{\frown}{|P_1P_2|}=2R\arcsin{\sqrt{\sin^2{\frac{\Delta\varphi}{2}}+\cos{\varphi_1}\cos{\varphi_2}\sin^2{\frac{\Delta\lambda}{2}}}}$

这便是开篇代码中用的Haversine公式推导的距离。

总结三个公式

对于球面上的两个点 $P_1$ 和 $P_2$ ，已知他们的经度 $\lambda$ 和维度 $\varphi$ ，就会有以下三种距离

$\stackrel{\frown}{|P_1P_2|}=2R\arcsin{\sqrt{\frac{{1-\sin{\varphi_1}\sin{\varphi_2}-\cos{\varphi_1}\cos{\varphi_2}\cos{\Delta\lambda}}}{2}}}$

$\stackrel{\frown}{|P_1P_2|}=R\arccos{(\cos{\varphi_1}\cos{\varphi_2}\cos{\Delta \lambda}+\sin{\varphi_1}\sin{\varphi_2})}$

$\stackrel{\frown}{|P_1P_2|}=2R\arcsin{\sqrt{\sin^2{\frac{\Delta\varphi}{2}}+\cos{\varphi_1}\cos{\varphi_2}\sin^2{\frac{\Delta\lambda}{2}}}}$

其实正半矢公式是余弦公式的一个转换，只是在相同纬度或者比较接近的两个位置距离求解的情况下更适用于计算机的精确计算，因为不容易丢失精度，但是由于摩尔定律，目前的64位的机器应该表现的越来越好，具体是否有所改善有待于实验。说到底求解的过程其实就是，当知道了两点的经度和纬度，然后求出球心到两点形成的两个向量的夹角，用夹角和地球半径就求出了球面距离；还有一种思路就是，直接在直角坐标系下求出两点的欧式距离，然后利用毕达哥拉斯定理(装逼说法，其实就是勾股定理)求出两个向量的夹角，然后同理求出球面距离。

参考文献

https://zh.wikipedia.org/wiki/%E5%A4%A7%E5%9C%86%E8%B7%9D%E7%A6%BB

根据经纬度对球面距离的求解

摘要

方法一：从弦长求大圆距离

方法二：用向量的余弦夹角公式求解大圆距离

方法三、用半正矢公式定义的距离

总结三个公式

参考文献

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读