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

摘要

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

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

最后编辑于：2019.07.25 14:24:03

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 199,711评论 5赞 468
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 83,932评论 2赞 376
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 146,770评论 0赞 330
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 53,799评论 1赞 271
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 62,697评论 5赞 359
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 48,069评论 1赞 276
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 37,535评论 3赞 390
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 36,200评论 0赞 254
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 40,353评论 1赞 294
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 35,290评论 2赞 317
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 37,331评论 1赞 329
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 33,020评论 3赞 315
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 38,610评论 3赞 303
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 29,694评论 0赞 19
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 30,927评论 1赞 255
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 42,330评论 2赞 346
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 41,904评论 2赞 341