根据经纬度计算距离公式的推导和 C# 实现

我们这里介绍一段据说是从 Google Maps 里面扒出来的代码：

看得有点懵？小编根据代码还原了其中的距离计算公式：

\[ y_{1}=\frac{lat_{1}}{180}\pi, y_{2}=\frac{lat_{2}}{180}\pi\]

\[ \Delta A = y_{1}-y_{2} \]

\[ \Delta B = \frac{lng_{1}}{180}\pi - \frac{lng_{2}}{180}\pi \]

\[S=2R\arcsin{\sqrt{{\sin^2{\frac{\Delta A }{2}}+\cos{y_{1}}\cos{y_{2}}{\sin^2{\frac{\Delta B }{2}}}}}}\]

其中：

\(S\)——两点之间的距离，单位 \(km\)；

\(R\)——地球半径，单位 \(km\)，本文取值为 \(R=6378.137 km\)；

\( \left(lng_{1},lat_{1}\right) \)——点 A 的经纬度坐标，\(lng_{1}\) 是经度，\(lat_{1}\) 是纬度；

\(\left(lng_{2},lat_{2}\right) \)——点 B 的经纬度坐标，\(lng_{2}\) 是经度，\(lat_{2}\) 是纬度。

那该算法是否可信呢？小编使用百度地图的测距功能进行了简单验证。结果发现，通过上述算法计算出的距离（\(208.7m\)）与通过百度地图测量出来的长度（\(209m\)）基本相等，这从一定程度上说明了算法的可靠性。百度地图的坐标拾取系统地址请戳这里。需要注意的是，国内的地图都是加偏的，而且不同的地图服务商采用的坐标系统存在一定的差异，详见《分享一款坐标转换软件：万能坐标转换 By 未来交通实验室》。

那么这个公式具体是怎么推导出来的呢？下面小编来扒一扒：

假设地球是半径为 \(R=6378.137 km\) 的标准球体，O为球心。A \(\left(lng_{1},lat_{1}\right)\)、B \(\left(lng_{2},lat_{2}\right)\) 是需要计算距离的两个点，C \(\left(lng_{1},lat_{2}\right)\)、D \(\left(lng_{2},lat_{1}\right)\) 为辅助点，AB、AC、AD、BC、BD 均为弦长，显然四边形 ACBD 为等腰梯形。O'为 A、D 两点所在纬度圆的圆心。AH 垂直于 OM，OE 垂直于 AC，AF 垂直于 BC。

根据经纬度坐标的定义，\( \angle MON\) 是 A、B 两点的经度之差，\( \angle AOC\) 是 A、B 两点的纬度之差，即：

\[ \angle AO'D = \angle MON= x_{2} - x_{1} \]

\[ \angle AOC= y_{1} - y_{2} \]

图中 \( \Delta AOC \) 为等腰三角形，因此：

\[  BD=AC=2R \sin{\frac{\angle AOC}{2}}=2R \sin{\frac{y_{1} - y_{2}}{2}} \]

同理：

\[  AD=2 R_{O'} \sin{\frac{\angle AO'D}{2}} \]

在直角三角形 AHO 中，有 \( R_{O'}=AO'=HO=R \cos{y_{1}}\)，则：

\[  AD=2 R \cos{y_{1}} \sin{\frac{x_{2} - x_{1}}{2}} \]

同理：

\[  CB=2 R \cos{y_{2}} \sin{\frac{x_{2} - x_{1}}{2}} \]

在等腰梯形 ACMD 中，有：

\[ CF= \frac{CB-AD}{2}，BF= \frac{CB+AD}{2}\]

根据勾股定理，有 \( AB^2= AF^2+BF^2，AC^2=CF^2+AF^2\)，则：

\[ AB= \sqrt{AF^2+BF^2}\\=\sqrt{AC^2-CF^2+BF^2}\\=\sqrt{AC^2-\frac{(CB-AD)^2}{4}+\frac{(CB+AD)^2}{4}}\\=\sqrt{AC^2+CB \times AD}\]

在等腰三角形 AOB 中，可以计算出 \(\angle AOB\)，即：

\[\angle AOB =\arcsin{\frac{AB}{2R}}= \arcsin{\frac{\sqrt{AC^2+CB \times AD}}{2R}}\]

代入 \( AC=2R \sin{\frac{y_{1} - y_{2}}{2}} ， CB=2 R \cos{y_{2}} \sin{\frac{x_{2} - x_{1}}{2}}，AD=2 R \cos{y_{1}} \sin{\frac{x_{2} - x_{1}}{2}}\)，整理后有：

\[\angle AOB = 2\arcsin{\sqrt{\sin^2{\frac{y_{1} - y_{2}}{2}}+\cos{y_{2}} \sin^2{\frac{x_{2} - x_{1}}{2}}}}\]

根据弧长公式，有：

\[ \widehat{AB} = R \times\angle AOB =2R\arcsin{\sqrt{\sin^2{\frac{y_{1} - y_{2}}{2}}+\cos{y_{2}} \cos{y_{1}} \sin^2{\frac{x_{2} - x_{1}}{2}}}}\]

推导结束。