WGS84/ПЗ-90椭球面上测地线交点计算及球面模型误差分析

问题背景 #

我找到了一套球面大圆交点的计算公式，但无法将其迁移至WGS84或ПЗ-90椭球面计算。具体疑问如下：

• 在纬度≤80°时，使用地球球面模型计算大圆交点的误差有多大？
• 若误差超过1米，如何结合地球椭球参数计算测地线交点？
• 尝试将ПЗ-90手册中的笛卡尔坐标转换公式与该球面算法结合，结果不正确，求解决思路？

附原球面算法代码：

LatLonSpherical LatLon_intersection(LatLonSpherical p1, LatLonSpherical p2, LatLonSpherical p3, LatLonSpherical p4) {
    //math from here: http://www.boeing-727.com/Data/fly%20odds/distance.html#arc2
    coords P1(p1), P2(p2), P3(p3), P4(p4);//转换为笛卡尔坐标系
    coords V1 = cross(P1, P2); //向量叉乘
    coords V2 = cross(P3, P4);
    //计算V1、V2的单位向量
    double length = V1.length();
    coords U1(V1.x / length, V1.y / length, V1.z / length);
    length = V2.length();
    coords U2(V2.x / length, V2.y / length, V2.z / length);
    //方向向量
    //coords D = U1 * U2;
    coords D = cross(V1, V2);
    length = D.length();
    //向量S1及其反方向S2即为球面上两个不重叠大圆的交点坐标
    coords S1(D.x / length, D.y / length, D.z / length);
    coords S2(-D.x / length, -D.y / length, -D.z / length);
    //将笛卡尔坐标系转换为地理坐标    
    LatLonSpherical point1(D.x, D.y, D.z);
    //检查交点是否在大圆弧段上
    LatLonSpherical s1(S1.x, S1.y, S1.z), s2(S2.x, S2.y, S2.z);
    if (abs(p1.distanceTo(p2) - p1.distanceTo(s1) - p2.distanceTo(s1)) < 1 &&   
        abs(p3.distanceTo(p4) - p3.distanceTo(s1) - p4.distanceTo(s1)) < 1) {
        return s1; //s1在弧段上
    }
    if (abs(p1.distanceTo(p2) - p1.distanceTo(s2) - p2.distanceTo(s2)) < 1 &&
        abs(p3.distanceTo(p4) - p3.distanceTo(s2) - p4.distanceTo(s2)) < 1) {
        return s2; //s2在弧段上
    }
    s2.error = -1;//返回错误
    return s2;
}

问题解答 #

1. 球面模型在纬度≤80°时的误差 #

地球是椭球体，赤道半径（WGS84为6378137m）比极半径（6356752m）约大21km。在纬度≤80°区域：

• 大圆交点的位置误差主要源于球面与椭球面的曲率差异。短距离弧段（几百公里内）误差通常在1-5米量级；跨半球长弧段误差可能达到几十米。
• 误差随纬度升高而增大，但80°以内大部分场景下误差不超过10米；若精度要求在1米内，球面模型无法满足。

2. 椭球面测地线交点的计算方法 #

椭球面上的测地线（大地线）交点无法像球面那样用简单向量叉乘直接求解，需采用迭代或数值方法，核心思路如下：

• 迭代逼近法：先使用球面模型计算初始交点作为迭代起点，将该点投影到椭球面上，再计算两条大地线在该点的切线方向，判断交点是否在弧段上，若不满足则调整坐标重新计算，直到精度达标（如误差<0.1米）。
• 大地线方程求解：利用椭球面大地线的微分方程（克莱罗方程），联立两条大地线的参数方程，通过牛顿迭代等数值解法求解交点坐标。
• 借助成熟库实现：直接使用专业地理计算库，如PROJ的geod_intersect函数、GeographicLib的Geodesic类，这些库已封装高精度的椭球面交点计算逻辑，支持WGS84、ПЗ-90等椭球参数。

3. 球面算法结合椭球坐标转换失败的原因及解决思路 #

失败核心原因是球面算法的笛卡尔坐标基于单位球/平均半径球，而椭球笛卡尔坐标基于椭球参数，两者坐标空间不兼容。解决思路：

• 若要复用球面算法，需先将椭球面上的点转换为球面投影点（用平均地球半径将经纬度转球面笛卡尔坐标），计算球面交点后再投影回椭球面；但这种方法精度提升有限，仍无法满足1米级要求。
• 追求高精度时，需放弃直接复用球面算法的思路，改用专门的椭球面大地线交点计算方法（如上述迭代法、方程求解法或专业库）。

内容的提问来源于stack exchange，提问作者Remerd