大地主题解算 Vincenty公式
椭球
起点纬度 起点经度
终点纬度 终点经度
距离
方位角

注:

©联系我们 E-mail:18801790654@163.com

Vincenty 的公式

数学基础可参考椭球大地测量(http://en.wikipedia.org/wiki/Geographical_distance)。

Vincenty公式的两种迭代方法用于大地测量中计算两点之间的球面距离,是 Thaddeus Vincenty (1975a)推导出来的,其假设地球表面是扁球,因此比认为地球是圆球要精确。

第一种方法(正算)是已知从一个点出发的距离和方位角(方向),计算端点的位置。第二种方法(反算)是已知两个点, 计算距离和方位角。公式被广泛应用于大地测量学,因为他们在椭球体上可准确至0.5毫米。

1. 背景

Vincenty的目标是将现有椭球大地测量算法用最小长度表达,他未发表的报告( 1975b)提到使用Wang720台式计算器,只有几千字节的内存。为了对长距离得到高精度, 采用基于辅助球的经典Legendre1806)、 Bessel1825)和Helmert1880) 解。Vincenty采用这种方法的公式是1955 Rainford给出的。Legendre证明椭球上的大地线可以完全投影至辅助圆上,地理纬度对应为 归化纬度,大圆的方位角等于大地线的方位角,椭球上沿大地线的经度和距离可在球上表示,沿大圆的弧长可容易积分出来。 BesselHelmert给出了收敛很快的积分的级数展开式,使得大地线可计算 至任意精度。

为了减小程序,Vincenty采用这些级数展开式,采用这些级数的第一项作为小量重新 展开,截断至 ,得到了经度和距离积分的简洁表达式,公式为嵌套形式,因为多项式的计算只用到一个临时变量。最后,正算反算都需要用简单的 迭代来解隐含公式,尽管计算速度慢(有时反算时不收敛),但代码很少。

2. 符号

符号定义如下:

a 椭球长半轴 WGS-84:6378137.0 m)
ƒ 椭球扁率 WGS-84:1/298.257223563)
b = (1 - ƒ) a 椭球短半轴
φ1, φ2 点的维度
U1 = arctan[(1 − ƒ) tan φ1]
U2 = arctan[(1 − ƒ) tan φ2]
归化纬度(即辅助圆上的纬度)sphere)
L = L2 - L1 两点的经度差
λ1, λ2 辅助球上的经度
α1, α2 点上的方位角
α 赤道处的方位角
s 椭球上的长度
σ 辅助球上的弧长

3. 反算

给定两点的坐标(φ1, L1)和(φ2, L2) ,反算问题是求方位角α1α2和椭球上的距离s

取初值λ = L,计算U1U2L,由下式计算,迭代至λ收敛:

\sin \sigma = \sqrt{ (\cos U_2 \sin \lambda)^2 + (\cos U_1 \sin U_2 - \sin U_1 \cos U_2 \cos \lambda)^2}

\cos \sigma = \sin U_1 \sin U_2 + \cos U_1 \cos U_2 \cos \lambda \,

\sigma = \arctan\frac{\sin\sigma}{\cos\sigma}\,

\sin \alpha = \frac{\cos U_1 \cos U_2 \sin \lambda}{\sin \sigma} \,

\cos^2 \alpha = 1 - \sin^2 \alpha \,

\cos (2 \sigma_m) = \cos \sigma - \frac{2 \sin U_1\sin U_2}{\cos^2 \alpha} \,

C = \frac{f}{16} \cos^2 \alpha \big[4 + f(4-3 \cos^2 \alpha) \big] \,

\lambda = L + (1-C) f \sin \alpha \left\{ \sigma + C \sin \sigma \left[\cos (2 \sigma_m) + C \cos \sigma (-1 + 2 \cos^2 (2 \sigma_m)) \right]\right\} \,

λ收敛至 (相当于0.06mm),计算以下公式:

u^2 = \cos^2 \alpha \frac{a^2 - b^2}{b^2} \,

A = 1 + \frac{u^2}{16384} \left\{ 4096 + u^2 \left[ -768 +u^2 (320 - 175u^2) \right] \right\}

B = \frac{u^2}{1024} \left\{ 256 + u^2 \left[ -128 + u^2 (74-47 u^2) \right] \right\}

 \Delta \sigma = B \sin \sigma \Big\{ \cos(2 \sigma_m) + \tfrac{1}{4} B \big[ \cos \sigma \big(-1+2 \cos^2(2 \sigma_m) \big) - \tfrac{1}{6} B \cos(2 \sigma_m)  (-3+4 \sin^2 \sigma) \big(-3+4 \cos^2 (2 \sigma_m)\big)  \big] \Big\}

 s = b A(\sigma - \Delta \sigma) \,

 \alpha_1 = \arctan \left( \frac{\cos U_2 \sin \lambda}{\cos U_1  \sin U_2 - \sin U_1 \cos U_2 \cos \lambda} \right)

 \alpha_2 = \arctan \left( \frac{\cos U_1 \sin \lambda}{-\sin U_1 \cos U_2 + \cos U_1 \sin U_2 \cos \lambda} \right)

当两个点为对应点时,公式可能不收敛,这是由于第一次计算的λ可能大于π

4. 正算

给定起点坐标(φ1, L1)和起始方位角α1、大地线长度s,求终点坐标 (φ2, L2)和反方位角α2

先计算以下量:

 \tan U_1 = (1 - f)\tan \phi_1 \,

 \sigma_1 = \arctan \left ( \frac{ \tan U_1}{ \cos \alpha_1} \right ) \,

 \sin \alpha = \cos U_1 \sin \alpha_1; \,\,\,\, \cos^2 \alpha = (1 - \sin \alpha)(1 + \sin \alpha)

 u^2 = \cos^2 \alpha \frac{a^2 - b^2}{b^2} \,

 A = 1 + \frac{u^2}{16384} \left\{ 4096 + u^2 \left[ -768 +u^2 (320 - 175u^2) \right] \right\}

 B = \frac{u^2}{1024} \left\{ 256 + u^2 \left[ -128 + u^2 (74-47 u^2) \right] \right\}

然后取初值  \sigma = \tfrac{s}{bA} ,迭代计算,直至σ不变:

 2 \sigma_m = 2 \sigma_1 + \sigma \,

 \Delta \sigma = B \sin \sigma \Big\{ \cos(2 \sigma_m) + \tfrac{1}{4} B \big[ \cos \sigma \big(-1+2 \cos^2(2 \sigma_m) \big) - \tfrac{1}{6} B \cos(2 \sigma_m) (-3+4 \sin^2 \sigma) \big(-3+4 \cos^2 (2 \sigma_m)\big) \big] \Big\}

 \sigma = \frac{s}{bA} + \Delta \sigma \,

得到足够准确的σ后,计算:

 \phi_2 = \arctan \left( \frac{\sin U_1 \cos \sigma + \cos U_1 \sin \sigma \cos \alpha_1}{(1 - f) \sqrt{\sin^2 \alpha + (\sin U_1 \sin \sigma - \cos U_1 \cos \sigma \cos \alpha_1 )^2 } } \right) \,

 \lambda = \arctan \left( \frac{\sin \sigma \sin \alpha_1}{\cos U_1 \cos \sigma - \sin U_1 \sin \sigma \cos \alpha_1} \right) \,

 C = \frac{f}{16} \cos^2 \alpha \big[4 + f(4-3 \cos^2 \alpha) \big] \,

 L = \lambda - (1-C) f \sin \alpha \left\{ \sigma + C \sin \sigma \left[\cos (2 \sigma_m) + C \cos \sigma (-1 + 2 \cos^2 (2 \sigma_m)) \right]\right\} \,

 \alpha_2 = \arctan \left( \frac{\sin \alpha}{-\sin U_1 \sin \sigma + \cos U_1 \cos \sigma \cos \alpha_1} \right) \,

若起点是南极或北极,第一个方程不能解,若起始方位角正东或正西,第二个方程不能解,但可采用 函数来回避 无解问题。

5. Vincenty的修正

Vincenty 1976年的信中,建议替换A和B的表达式,其中k1采用Helmert的展开参数:

A = \frac {1 + \frac {1}{4} (k_1)^2}{1 - k_1}

B = k_1(1 - \tfrac {3}{8}(k_1)^2)

其中:

 k_1 = \frac { \sqrt {(1 + u^2)} - 1}{ \sqrt {(1 + u^2)} + 1}

6. 近似对应点

如前所述,对反算问题,若两点近似为对应点,则不能收敛或收敛很慢,例如,若取WGS84椭球, (φ1, L1) = (0°, 0°) 、(φ2, L2) = (0.5°, 179.5°), 需要迭代130次才能得到1mm的精度。算法可能给出准确的结果(19936288.579 m),也可能得到错误结果或无解。