为什么SqlGeography STDistance有区别？

Question

我一直在追踪我的代码中的一个bug，我发现这是因为微软c# SqlGeography 2014库为STDistance返回的结果与计算点间距离的常规代码略有不同。

我编写了一个小型控制台exe来演示这个问题，但我不知道为什么会得到如此不同的结果？

static void Main(string[] args) {
        double earthRadius = 6378137; // meters => from both nad83 & wgs84
        var a = new { lat = 43.68151632, lng = -79.61162263 };
        var b = new { lat = 43.67575602, lng = -79.59586143 };

        // sql geography lib
        SqlGeographyBuilder sgb;
        sgb = new SqlGeographyBuilder();
        sgb.SetSrid(4326);
        sgb.BeginGeography(OpenGisGeographyType.Point);
        sgb.BeginFigure(a.lat, a.lng);
        sgb.EndFigure();
        sgb.EndGeography();
        SqlGeography geoA = sgb.ConstructedGeography;

        sgb = new SqlGeographyBuilder();
        sgb.SetSrid(4326);
        sgb.BeginGeography(OpenGisGeographyType.Point);
        sgb.BeginFigure(b.lat, b.lng);
        sgb.EndFigure();
        sgb.EndGeography();
        SqlGeography geoB = sgb.ConstructedGeography;

        // distance cast from SqlDouble
        double geoDistance = (double)geoA.STDistance(geoB);

        // math!
        double d2r = Math.PI / 180; // for converting degrees to radians
        double lat1 = a.lat * d2r,
            lat2 = b.lat * d2r,
            lng1 = a.lng * d2r,
            lng2 = b.lng * d2r,
            dLat = lat2 - lat1,
            dLng = lng2 - lng1,
            sin_dLat_half = Math.Pow(Math.Sin(dLat / 2), 2),
            sin_dLng_half = Math.Pow(Math.Sin(dLng / 2), 2),
            distance = sin_dLat_half + Math.Cos(lat1) * Math.Cos(lat2) * sin_dLng_half;

        // math distance
        double mathDistance = (2 * Math.Atan2(Math.Sqrt(distance), Math.Sqrt(1 - distance))) * earthRadius;

        // haversine
        double sLat1 = Math.Sin(a.lat * d2r),
                sLat2 = Math.Sin(b.lat * d2r),
                cLat1 = Math.Cos(a.lat * d2r),
                cLat2 = Math.Cos(b.lat * d2r),
                cLon = Math.Cos((a.lng * d2r) - (b.lng * d2r)),
                cosD = sLat1 * sLat2 + cLat1 * cLat2 * cLon,
                d = Math.Acos(cosD);

        // math distance
        double methDistance = d * earthRadius;


        // write the outputs
        Console.WriteLine("geo distance:\t" + geoDistance);    // 1422.99560435875
        Console.WriteLine("math distance:\t" + mathDistance);  // 1421.73656776243
        Console.WriteLine("meth distance:\t" + methDistance);  // 1421.73656680185
        Console.WriteLine("geo vs math:\t" + (geoDistance - mathDistance));     // 1.25903659632445
        Console.WriteLine("haversine vs math:\t" + (methDistance - methDistance)); // ~0.00000096058011
    }

微软是否使用了不同的计算方法？当计算距离小于1.5Km时，距离超过1米是的一个巨大的误差。

Answer 1

好的，在深入挖掘之后，我找到了答案，而微软是更正确的。

具体来说，他们使用的是文琴蒂公式。精度在0.5mm (不是米，半毫米)，而不是0.3%与Haversine公式。

原因是Haversine (我、Google和Bing地图显然也使用过)速度快，但依赖的是球形地球而不是椭球体。微软使用椭球地球来计算距离，而不是用球体来提供更精确的结果。

我在c#中像这样实现了Vincenty的方法，到目前为止它还在工作，但是还没有接近生产的地方。

    const double d2r = Math.PI / 180;  // degrees to radians
    const double EARTH_RADIUS = 6378137;  // meters
    const double EARTH_ELLIPSOID = 298.257223563; // wgs84
    const double EARTH_BESSEL = 1 / EARTH_ELLIPSOID;
    const double EARTH_RADIUS_MINOR = EARTH_RADIUS - (EARTH_RADIUS * EARTH_BESSEL); // 6356752.3142 meters => wgs84

    static double vincentyDistance(double lat1, double lng1, double lat2, double lng2) {
        double L = (lng2 - lng1) * d2r,
                U1 = Math.Atan((1 - EARTH_BESSEL) * Math.Tan(lat1 * d2r)),
                U2 = Math.Atan((1 - EARTH_BESSEL) * Math.Tan(lat2 * d2r)),
                sinU1 = Math.Sin(U1),
                cosU1 = Math.Cos(U1),
                sinU2 = Math.Sin(U2),
                cosU2 = Math.Cos(U2),
                lambda = L,
                lambdaP,
                iterLimit = 100,
                sinLambda,
                cosLambda,
                sinSigma,
                cosSigma,
                sigma,
                sinAlpha,
                cosSqAlpha,
                cos2SigmaM,
                C;
        do {
            sinLambda = Math.Sin(lambda);
            cosLambda = Math.Cos(lambda);
            sinSigma = Math.Sqrt((cosU2 * sinLambda) * (cosU2 * sinLambda) + (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) * (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda));
            if (0 == sinSigma) {
                return 0; // co-incident points
            };
            cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
            sigma = Math.Atan2(sinSigma, cosSigma);
            sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;
            cosSqAlpha = 1 - sinAlpha * sinAlpha;
            cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
            C = EARTH_BESSEL / 16 * cosSqAlpha * (4 + EARTH_BESSEL * (4 - 3 * cosSqAlpha));
            //  if (isNaN(cos2SigmaM)) {
            //      cos2SigmaM = 0; // equatorial line: cosSqAlpha = 0 (§6)
            //  };
            lambdaP = lambda;
            lambda = L + (1 - C) * EARTH_BESSEL * sinAlpha * (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)));
        } while (Math.Abs(lambda - lambdaP) > 1e-12 && --iterLimit > 0);

        if (iterLimit == 0) {
            return 0; // formula failed to converge
        };

        double uSq = cosSqAlpha * (EARTH_RADIUS * EARTH_RADIUS - EARTH_RADIUS_MINOR * EARTH_RADIUS_MINOR) / (EARTH_RADIUS_MINOR * EARTH_RADIUS_MINOR),
            A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq))),
            B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq))),
            deltaSigma = B * sinSigma * (cos2SigmaM + B / 4 * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) - B / 6 * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM * cos2SigmaM))),
            s = EARTH_RADIUS_MINOR * A * (sigma - deltaSigma);
        return s;
    }

这段代码是从我在这里找到的一个JavaScript实现转换而来的：https://gist.github.com/mathiasbynens/354587

Answer 2

微软可能更正确，但它是不正确的！

根据Kallay (2010年)，SqlGeography STDistance不返回测地线的距离，而是返回连接两点的大椭圆的距离。大椭圆距离不是最短距离，它不服从三角形不等式。大椭圆返回近地点测地线距离的合理近似。对于远距离点，误差可能高达33公里。

C#库NETGeographicLib (我编写了底层C++库)给出了测地线距离的精确计算。这比Vincenty (10纳米而不是0.5毫米)更精确，更重要的是，总是得到正确的结果。(Vincenty没有收敛到几乎相反的点。)