为什么助推:几何地理Vincenty绕赤道的距离不准确

Question

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

助推几何设计理性状态：

有Andoyer方法，快速而精确，还有Vincenty方法，较慢且更精确。。

然而，当用Andoyer和Vincenty策略测试boost::geometry::distance函数时，我得到了以下结果：

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

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

半小调和Andoyer的距离看起来是合理的。除了当点在地球的另一侧时，当boost 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:tt"   << a << "n";
  std::cout << "tFlattening:tt"       << f << "n";
  std::cout << "tSemiminor axis:tt"   << b << "nn";
  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 polesn";
  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 << "nn";
  // 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 polesn";
  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

我听从@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相同的距离（到5dp），并为赤道上的对点计算了相同的距离。

这是正确的，因为赤道上对极点之间的最短距离是通过其中一个极点，而不是默认的助推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:tt"   << a << "n";
  std::cout << "tFlattening:tt"       << f << "n";
  std::cout << "tSemiminor axis:tt"   << b << "nn";
  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 antipodaln";
  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 << "nn";
  std::pair<double, double> equator_east (0.0, 90.0);
  std::pair<double, double> equator_west (0.0, -90.0);
  std::cout << "GeographicLib antipodaln";
  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 << "nn";
  std::pair<double, double> JFK   (40.6, -73.8);
  std::pair<double, double> LHR   (51.6, -0.5);
  std::cout << "GeographicLib verify distancen";
  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;
}

在他的论文《测地线的算法》中，Charles F.F.Karney指出"Vincenty的方法几乎不能收敛于对足点"。这可能回答了我最初的问题，即Vincenty算法不适用于对跖点。

注：我提出了boost票证#111817，描述了以下问题CCD_ 14算法针对对足点返回零并且向CCD_

然而，对于不正确的距离，唯一正确的解决方案是使用正确的算法，即：geographiclib

非常感谢Charles F.F.Karney（@cfk）礼貌地指出了我愚蠢的错误！

Answer 2

作为替代方案，请查看Charles F.F.Karney的地理信息库。正如文件所说："重点是返回准确的结果，误差接近四舍五入（约5-15纳米）。"