球面两点距离计算公式

两点之间的距离怎么定义?距离,就是路径的长度。显然从一点到另一点有无数条路径,而每一条路径都有一个长度,那么我们用哪一条路径的长度来定义两点之间距离呢?毫无疑问,“最短长度”!为什么?因为它特殊,最短!

小时候我们就知道了,“平面上,两点之间直线最短”。平面上的两点距离可以用欧氏距离公式 √x2+y2 表示。不过我们现在关心的是地球上两点之间的距离。不幸的是地球近似一个球形(R ≈ 6378137m),表面近似一个球面,并不是一个平面。我们不能用欧式的直线距离,那样我们会是穿过地球内部从一点达到另外一点。(至于直线为嘛会穿过地球内部?因为球面是一个凸面!) 那球面上的最短距离是什么呢?

球面上的最短距离,也就是测地线距离,是两点与球心确定的大圆在球面上划过的弧的长度。

这里并不来论证为嘛是大圆的弧线长度,我们更关心如何计算。距离计算的第一步就是用坐标来表示点。地球上的点通常用经纬度来标识,譬如杭州阿里巴巴西溪园区的坐标:120.029993(经度),30.286352(维度)。现在为了计算一下阿里巴巴西溪园区到翠苑三区(坐标:120.129107,30.297394)的球面距离,我们先推导一下计算公式。

两点距离公式

空间中的点除了用笛卡尔坐标 (x, y, z) 描述之外,还可以使用球面坐标 (r,θ,φ) 表示。而地球上某个点的经纬度表示,正是球面坐标的变种,r 就是地球半径 R,φ 就是经度,θ 是 90° - 维度。根据球面坐标与笛卡尔坐标的转换关系,我们可以有:

假设地球上点 A 的经纬度为 (α,β),那么
xA = Rcosβcosα
yA = Rcosβsinα
zA = Rsinβ

假设球心为 O,那么 OA 和 OB 的夹角是:
夹角 γ = cos-1(OA·OB/(|OA|·|OB|)) = cos-1(cosβAcosβBcos(αAB)+sinβAsinβB)

那么 O、A、B 确定的大圆在球面上划过的弧的长度为:
dAB = Rγ = Rcos-1(cosβAcosβBcos(αAB)+sinβAsinβB)

计算距离的程序

知道了公式写程序其实是很简单的一件事。一般语言的库都包含了常见数学函数的实现,不过需要注意的是,日常中用的角度和科学中用的弧度是不一样的,而库函数譬如 sin(x) 的参数通常是弧度表示。下面是使用go实现的一个计算公式,仅供参考:

  1. package geo_ext
  2.  
  3. import (
  4.         "math"
  5. )
  6.  
  7. type Position struct {
  8.         Lon float64 // 经度
  9.         Lat float64 // 维度
  10. }
  11.  
  12. func Distence(p, q Position) float64 {
  13.         var α1, β1, α2, β2 = rad(p.Lon), rad(p.Lat), rad(q.Lon), rad(q.Lat)
  14.  
  15.         // L = R*arccos(cos(β1)cos(β2)cos(α1-α2)+sin(β1)sin(β2)) , (
  16.         const r = 6378137
  17.         return r * math.Acos(math.Cos(β1)*math.Cos(β2)*math.Cos(α1-α2)+math.Sin(β1)*math.Sin(β2))
  18. }
  19.  
  20. func rad(d float64) float64 {
  21.         return d * math.Pi / 180.0
  22. }

回归问题

现在就回归最初的计算阿里巴巴西溪园区到翠苑三区的距离。

  1. package main
  2.  
  3. import (
  4.         "fmt"
  5.  
  6.         "github.com/antark/golang-x/geo_ext"
  7. )
  8.  
  9. func main() {
  10.         p, q := geo_ext.Position{120.029993, 30.286352}, // 阿里巴巴西溪园区
  11.                 geo_ext.Position{120.129107, 30.297394} // 翠苑三区
  12.         fmt.Println(geo_ext.Distence(p, q)) // 9605.87811540484
  13. }
  14. </pre>
  15. </div>

Tags: 

Article type: