Python - 文森蒂的逆公式不收敛（寻找地球上点之间的距离）

Question

我正在尝试实现 Vincenty 的逆问题，如 wiki HERE 中所述

问题在于 lambda 根本没有收敛。如果我尝试迭代公式序列，该值保持不变，我真的不知道为什么。也许我只是对一个明显的问题视而不见。

需要注意的是，我是 Python 新手，还在学习这门语言，所以我不确定是否是语言的误用导致了问题，或者我在某些计算中确实存在一些错误我执行。我似乎在公式中找不到任何错误。

基本上，我以尽可能接近 wiki 文章的格式编写代码，结果是这样的：

import math

# Length of radius at equator of the ellipsoid
a = 6378137.0

# Flattening of the ellipsoid
f = 1/298.257223563

# Length of radius at the poles of the ellipsoid
b = (1 - f) * a

# Latitude points
la1, la2 = 10, 60

# Longitude points
lo1, lo2 = 5, 150

# For the inverse problem, we calculate U1, U2 and L.
# We set the initial value of lamb = L
u1 = math.atan( (1 - f) * math.tan(la1) )
u2 = math.atan( (1 - f) * math.tan(la2) )
L = (lo2 - lo1) * 0.0174532925

lamb = L

while True:
    sinArc = math.sqrt( math.pow(math.cos(u2) * math.sin(lamb),2) + math.pow(math.cos(u1) * math.sin(u2) - math.sin(u1) * math.cos(u2) * math.cos(lamb),2) )
    cosArc = math.sin(u1) * math.sin(u2) + math.cos(u1) * math.cos(u2) * math.cos(lamb)
    arc = math.atan2(sinArc, cosArc)
    sinAzimuth = ( math.cos(u1) * math.cos(u2) * math.sin(lamb) ) // ( sinArc )
    cosAzimuthSqr = 1 - math.pow(sinAzimuth, 2)
    cosProduct = cosArc - ((2 * math.sin(u1) * math.sin(u2) ) // (cosAzimuthSqr))
    C = (f//16) * cosAzimuthSqr  * (4 + f * (4 - 3 * cosAzimuthSqr))
    lamb = L + (1 - C) * f * sinAzimuth * ( arc + C * sinArc * ( cosProduct + C * cosArc * (-1 + 2 * math.pow(cosProduct, 2))))
    print(lamb)

如前所述，问题是“lamb”（lambda）的值不会变小。我什至尝试将我的代码与其他实现进行比较，但它们看起来几乎相同。

我在这里做错了什么？ :-)

谢谢大家！

Answer 1

首先，您也应该将纬度转换为弧度（您已经对经度进行了此操作）：

u1 = math.atan( (1 - f) * math.tan(math.radians(la1)) )
u2 = math.atan( (1 - f) * math.tan(math.radians(la2)) )
L = math.radians((lo2 - lo1)) # better than * 0.0174532925

一旦您这样做并摆脱 //（int 部门）并用 /（float 部门）替换它们，lambda 将停止在您的迭代中重复相同的值并开始遵循此路径（基于您的示例坐标）：

2.5325205864224847
2.5325167509030906
2.532516759118641
2.532516759101044
2.5325167591010813
2.5325167591010813
2.5325167591010813

您似乎期望 10^(−12) 的收敛精度，这似乎说明了这一点。

您现在可以退出循环（lambda 已收敛）并继续前进，直到您计算出所需的测地线距离 s。

注意：you can test your final value s here.

Answer 2

即使正确实现，Vincenty 的算法也会失败 收敛一些点。 （文森蒂注意到了这个问题。） 我给出了一个保证的算法 汇入Algorithms for geodesics；有一条蟒蛇 实现可用here。最后，您可以找到更多 维基百科页面上有关问题的信息， Geodesics on an ellipsoid。 （讨论页有例子 由 NGS 实施的 Vincenty 的成对点， 无法收敛。）