Python - Vincenty的反算公式不收敛（寻找地球上两点间的距离）

Question

Python - Vincenty的反算公式不收敛（寻找地球上两点间的距离）

pythonmathformulawiki

4

我试图按照维汶提反问题的方法实现维基上这里描述的内容。

问题在于lambda根本无法收敛。如果我尝试迭代公式序列，它的值就会保持不变，我真的不确定为什么。也许我只是对一个明显的问题盯得太久而看不到。

需要指出的是，我是Python的新手，还在学习语言，因此我不确定是否使用了错误的语言可能导致问题，或者我在执行某些计算时有错误。我只是似乎找不到公式中的任何错误。

基本上，我尽可能地按照维基文章的格式编写了代码，结果如下：

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）的值不会变小。我甚至尝试将我的代码与其他实现进行比较，但它们看起来几乎一样。

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

谢谢大家！

- CodingBeagle

1

在你的例子中，lamb 应该收敛到哪个值？ - Jivan

尝试使用 from __future__ import division。 - Jonathan Davies

1

然后，我认为您在链接页面上阅读的数学写作中存在误解。当他们说 sin α = blabla 时，它并不意味着 sinAlpha = blabla，而是意味着 α = arcsin blablabla。 - Jivan

1

@Bitious 在 L = (lo2 - lo1) * 0.0174532925 中，为什么要乘以 * 0.017...？维基百科的文章中没有提到这一点。 - Jivan

2

是的，但我在谈论的是前面两行中的la1和la2。 - Jim Lewis

显示剩余10条评论

2个回答

3

即使 Vincenty 算法被正确实现，也会在某些点上无法收敛。（Vincenty 注意到了这个问题。）我提供了一个可以保证收敛的算法，在 Algorithms for geodesics 中有介绍；这里有一个 Python 实现 here。最后，你可以在维基百科页面 Geodesics on an ellipsoid 上找到更多关于这个问题的信息。（讨论页面有一些例子，其中 Vincenty 由 NGS 实现时无法收敛。）

- cffk

我已经向geopy提交了一个补丁，以便在可用时使用准确的geographiclib算法。 - cffk

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Jivan · Accepted Answer

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

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。

注意：您可以在此处测试最终值s。