爲什麼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）的禮貌指出我的愚蠢的錯誤！