为什么boost :: geometry geographic Vincenty距离在赤道附近不准确？

Question

我需要一个函数来高精度地计算一对WGS 84位置之间的距离，我打算使用boost geometry中的geographic函数。

的boost geometry Design Rational状态：

还有就是Andoyer方法，快速和精确的，并且在Vincenty方法，更慢和更精确..

然而，测试boost::geometry::distance功能时与Andoyer和Vincenty战略，我得到以下结果：

WGS 84 values (metres) 
    Semimajor axis:   6378137.000000 
    Flattening:    0.003353 
    Semiminor axis:   6356752.314245 

    Semimajor distance:  20037508.342789 
    Semiminor distance:  19970326.371123 

Boost geometry near poles 
Andoyer function: 
    Semimajor distance:  20037508.151445 
    Semiminor distance:  20003917.164970 
Vincenty function: 
    Semimajor distance:  **19970326.180419** 
    Semiminor distance:  20003931.266635 

Boost geometry at poles 
Andoyer function: 
    Semimajor distance:  0.000000 
    Semiminor distance:  0.000000 
Vincenty function: 
    Semimajor distance:  **19970326.371122** 
    Semiminor distance:  20003931.458623 

沿着半长轴的距离（即，在赤道周围）小于北极和南极之间的半微米轴周围的距离。这不可能是正确的。

Semiminor和Andoyer距离看起来合理。除了点位于地球的另一侧时，当功能返回零时！

问题出在：Vincenty算法，boost geometry执行它，还是我的测试代码？

测试代码：

/// boost geometry WGS84 distance issue 

// Note: M_PI is not part of the C or C++ standards, _USE_MATH_DEFINES enables it 
#define _USE_MATH_DEFINES 
#include <boost/geometry.hpp> 
#include <cmath> 
#include <iostream> 
#include <ios> 

// WGS 84 parameters from: Eurocontrol WGS 84 Implementation Manual 
// Version 2.4 Chapter 3, page 14 

/// The Semimajor axis measured in metres. 
/// This is the radius at the equator. 
constexpr double a = 6378137.0; 

/// Flattening, a ratio. 
/// This is the flattening of the ellipse at the poles 
constexpr double f = 1.0/298.257223563; 

/// The Semiminor axis measured in metres. 
/// This is the radius at the poles. 
/// Note: this is derived from the Semimajor axis and the flattening. 
/// See WGS 84 Implementation Manual equation B-2, page 69. 
constexpr double b = a * (1.0 - f); 

int main(int /*argc*/, char ** /*argv*/) 
{ 
    std::cout.setf(std::ios::fixed); 

    std::cout << "WGS 84 values (metres)\n"; 
    std::cout << "\tSemimajor axis:\t\t" << a << "\n"; 
    std::cout << "\tFlattening:\t\t"  << f << "\n"; 
    std::cout << "\tSemiminor axis:\t\t" << b << "\n\n"; 

    std::cout << "\tSemimajor distance:\t" << M_PI * a << "\n"; 
    std::cout << "\tSemiminor distance:\t" << M_PI * b << "\n"; 
    std::cout << std::endl; 

    // Min value for delta. 0.000000014 causes Andoyer to fail. 
    const double DELTA(0.000000015); 

    // For boost::geometry: 
    typedef boost::geometry::cs::geographic<boost::geometry::radian> Wgs84Coords; 
    typedef boost::geometry::model::point<double, 2, Wgs84Coords> GeographicPoint; 
    // Note boost points are Long & Lat NOT Lat & Long 
    GeographicPoint near_north_pole (0.0, M_PI_2 - DELTA); 
    GeographicPoint near_south_pole (0.0, -M_PI_2 + DELTA); 

    GeographicPoint near_equator_east (M_PI_2 - DELTA, 0.0); 
    GeographicPoint near_equator_west (-M_PI_2 + DELTA, 0.0); 

    // Note: the default boost geometry spheroid is WGS84 
    // #include <boost/geometry/core/srs.hpp> 
    typedef boost::geometry::srs::spheroid<double> SpheroidType; 
    SpheroidType spheriod; 

    //#include <boost/geometry/strategies/geographic/distance_andoyer.hpp> 
    typedef boost::geometry::strategy::distance::andoyer<SpheroidType> 
                   AndoyerStrategy; 
    AndoyerStrategy andoyer(spheriod); 

    std::cout << "Boost geometry near poles\n"; 
    std::cout << "Andoyer function:\n"; 
    double andoyer_major(boost::geometry::distance(near_equator_east, near_equator_west, andoyer)); 
    std::cout << "\tSemimajor distance:\t" << andoyer_major << "\n"; 
    double andoyer_minor(boost::geometry::distance(near_north_pole, near_south_pole, andoyer)); 
    std::cout << "\tSemiminor distance:\t" << andoyer_minor << "\n"; 

    //#include <boost/geometry/strategies/geographic/distance_vincenty.hpp> 
    typedef boost::geometry::strategy::distance::vincenty<SpheroidType> 
                   VincentyStrategy; 
    VincentyStrategy vincenty(spheriod); 

    std::cout << "Vincenty function:\n"; 
    double vincenty_major(boost::geometry::distance(near_equator_east, near_equator_west, vincenty)); 
    std::cout << "\tSemimajor distance:\t" << vincenty_major << "\n"; 
    double vincenty_minor(boost::geometry::distance(near_north_pole, near_south_pole, vincenty)); 
    std::cout << "\tSemiminor distance:\t" << vincenty_minor << "\n\n"; 

    // Note boost points are Long & Lat NOT Lat & Long 
    GeographicPoint north_pole (0.0, M_PI_2); 
    GeographicPoint south_pole (0.0, -M_PI_2); 

    GeographicPoint equator_east (M_PI_2, 0.0); 
    GeographicPoint equator_west (-M_PI_2, 0.0); 

    std::cout << "Boost geometry at poles\n"; 
    std::cout << "Andoyer function:\n"; 
    andoyer_major = boost::geometry::distance(equator_east, equator_west, andoyer); 
    std::cout << "\tSemimajor distance:\t" << andoyer_major << "\n"; 
    andoyer_minor = boost::geometry::distance(north_pole, south_pole, andoyer); 
    std::cout << "\tSemiminor distance:\t" << andoyer_minor << "\n"; 

    std::cout << "Vincenty function:\n"; 
    vincenty_major = boost::geometry::distance(equator_east, equator_west, vincenty); 
    std::cout << "\tSemimajor distance:\t" << vincenty_major << "\n"; 
    vincenty_minor = boost::geometry::distance(north_pole, south_pole, vincenty); 
    std::cout << "\tSemiminor distance:\t" << vincenty_minor << "\n"; 

    return 0; 
} 

Answer 1

作为替代退房查尔斯F. F. Karney的geographiclib。正如文件中所述：“重点在于返回准确的结果，并且误差接近四舍五入（约5-15纳米）。”

Answer 2

我按照@ jwd630的建议和检查了geographiclib。
下面是结果：

WGS 84 values (metres) 
    Semimajor distance: 20037508.342789 
    Semiminor distance: 19970326.371123 

GeographicLib near antipodal 
    Semimajor distance: 20003931.458625 
    Semiminor distance: 20003931.455275 

GeographicLib antipodal 
    Semimajor distance: 20003931.458625 
    Semiminor distance: 20003931.458625 

GeographicLib verify 
    JFK to LHR distance: 5551759.400319 

即它提供了与Vincenty相同的距离，用于极点之间的Semiminor距离（至5dp），并计算赤道对极点的相同距离。

这是正确的，因为在Equator上的对映点之间的最短距离是通过其中一个极点，而不是赤道周围的，因为默认提升算法会计算Andoyer。

所以@JWD630的回答是正确的，三种算法中，geographiclib是唯一一个计算整个WGS84大地水准面的正确距离。

下面是测试代码：

/// GeographicLib WGS84 distance 

// Note: M_PI is not part of the C or C++ standards, _USE_MATH_DEFINES enables it 
#define _USE_MATH_DEFINES 
#include <GeographicLib/Geodesic.hpp> 
#include <cmath> 
#include <iostream> 
#include <ios> 

// WGS 84 parameters from: Eurocontrol WGS 84 Implementation Manual 
// Version 2.4 Chapter 3, page 14 

/// The Semimajor axis measured in metres. 
/// This is the radius at the equator. 
constexpr double a = 6378137.0; 

/// Flattening, a ratio. 
/// This is the flattening of the ellipse at the poles 
constexpr double f = 1.0/298.257223563; 

/// The Semiminor axis measured in metres. 
/// This is the radius at the poles. 
/// Note: this is derived from the Semimajor axis and the flattening. 
/// See WGS 84 Implementation Manual equation B-2, page 69. 
constexpr double b = a * (1.0 - f); 

int main(int /*argc*/, char ** /*argv*/) 
{ 
    const GeographicLib::Geodesic& geod(GeographicLib::Geodesic::WGS84()); 

    std::cout.setf(std::ios::fixed); 

    std::cout << "WGS 84 values (metres)\n"; 
    std::cout << "\tSemimajor axis:\t\t" << a << "\n"; 
    std::cout << "\tFlattening:\t\t"  << f << "\n"; 
    std::cout << "\tSemiminor axis:\t\t" << b << "\n\n"; 

    std::cout << "\tSemimajor distance:\t" << M_PI * a << "\n"; 
    std::cout << "\tSemiminor distance:\t" << M_PI * b << "\n"; 
    std::cout << std::endl; 

    // Min value for delta. 0.000000014 causes boost Andoyer to fail. 
    const double DELTA(0.000000015); 

    std::pair<double, double> near_equator_east (0.0, 90.0 - DELTA); 
    std::pair<double, double> near_equator_west (0.0, -90.0 + DELTA); 

    std::cout << "GeographicLib near antipodal\n"; 
    double distance_metres(0.0); 
    geod.Inverse(near_equator_east.first, near_equator_east.second, 
       near_equator_west.first, near_equator_west.second, distance_metres); 
    std::cout << "\tSemimajor distance:\t" << distance_metres << "\n"; 

    std::pair<double, double> near_north_pole (90.0 - DELTA, 0.0); 
    std::pair<double, double> near_south_pole (-90.0 + DELTA, 0.0); 

    geod.Inverse(near_north_pole.first, near_north_pole.second, 
       near_south_pole.first, near_south_pole.second, distance_metres); 
    std::cout << "\tSemiminor distance:\t" << distance_metres << "\n\n"; 

    std::pair<double, double> equator_east (0.0, 90.0); 
    std::pair<double, double> equator_west (0.0, -90.0); 

    std::cout << "GeographicLib antipodal\n"; 
    geod.Inverse(equator_east.first, equator_east.second, 
       equator_west.first, equator_west.second, distance_metres); 
    std::cout << "\tSemimajor distance:\t" << distance_metres << "\n"; 

    std::pair<double, double> north_pole (90.0, 0.0); 
    std::pair<double, double> south_pole (-90.0, 0.0); 

    geod.Inverse(north_pole.first, north_pole.second, 
       south_pole.first, south_pole.second, distance_metres); 
    std::cout << "\tSemiminor distance:\t" << distance_metres << "\n\n"; 

    std::pair<double, double> JFK (40.6, -73.8); 
    std::pair<double, double> LHR (51.6, -0.5); 

    std::cout << "GeographicLib verify distance\n"; 
    geod.Inverse(JFK.first, JFK.second, 
       LHR.first, LHR.second, distance_metres); 
    std::cout << "\tJFK to LHR distance:\t" << distance_metres << std::endl; 

    return 0; 
} 

在他的论文Algorithms for geodesics， 查尔斯F. F. Karney指出：“Vincenty的方法收敛失败了将近径点”。 这可能会回答我的原始问题，即Vincenty算法不适用于反对点。

注：我已经提出boost票#11817描述，其中 的Andoyer算法返回零径点和发送拉请求boost与它一个解决这个问题。

但是，不正确的距离，唯一正确的解决方法是使用正确的算法，即：geographiclib

非常感谢查尔斯F. F. Karney（@cffk）的礼貌指出我的愚蠢的错误！