求两条大地线经纬度交点的Java代码问题求助

问题概述 #

我编写的Java代码用于计算地球表面两条大地线的交点（考虑地球曲率），结果经度正确，但纬度偏差明显——正确值约为-26.5098，代码返回-26.373957。另外我怀疑当前使用的Haversine公式未考虑地球扁率，可能是问题根源。代码由网上片段整合而来，对背后数学原理不甚了解，希望明确修正方向。

核心错误分析 #

1. 平面相交法的本质错误 #

当前代码通过大地线端点的笛卡尔向量叉乘得到过地心平面的法向量，再求两个平面的交线来获取交点。但椭球上的大地线并非过地心平面与椭球的交线——只有球面上的大圆才符合这个特性，椭球大地线是测地线，轨迹不在过地心平面上，这是纬度计算错误的核心原因。

2. Haversine公式的局限性 #

Haversine公式基于球体模型计算距离，完全忽略地球扁率（WGS84是椭球），用它判断点是否在大地线上会产生系统误差。同时代码中用ah+bh == ch的精确相等判断，由于浮点数精度问题，几乎不可能命中，会误判交点是否在大地线上。

3. 未归一化的笛卡尔向量问题 #

代码注释掉了向量归一化步骤，直接用平面交线的原始向量转换大地坐标，该向量的模长不等于椭球半径，转换出的点根本不在椭球表面，必然导致纬度错误。

修正步骤 #

1. 替换核心计算逻辑：放弃平面相交法 #

椭球大地线的交点计算需要用到测地线的微分方程或专门的数值解法，不要自己从零实现，建议使用成熟的 geodesy 工具类，或基于Vincenty算法扩展实现交点计算。

2. 替换Haversine公式为椭球距离算法 #

用支持WGS84椭球的Vincenty逆公式计算大地距离，同时修改点在大地线上的判断逻辑：用误差阈值（如1毫米）比较ah+bh与ch的差值，而非精确相等。

3. 确保计算点在椭球表面 #

所有参与计算的点必须是椭球表面的点，笛卡尔向量需对应到椭球上的坐标，而非任意向量。

代码调整建议 #

（1）修复点在大地线上的判断逻辑 #

替换Haversine公式为Vincenty逆公式，并添加误差阈值判断：

private static boolean isPointOnGeodesic(double[] gp, double[] ls, double[] le) {
    double ah = vincentyDistance(ls[0], ls[1], gp[0], gp[1]);
    double bh = vincentyDistance(gp[0], gp[1], le[0], le[1]);
    double ch = vincentyDistance(ls[0], ls[1], le[0], le[1]);
    // 允许1毫米的误差范围
    return Math.abs(ah + bh - ch) < 1e-3;
}

// Vincenty逆公式：计算WGS84椭球上两点的大地距离
private static double vincentyDistance(double lat1, double lon1, double lat2, double lon2) {
    double lat1Rad = Math.toRadians(lat1);
    double lon1Rad = Math.toRadians(lon1);
    double lat2Rad = Math.toRadians(lat2);
    double lon2Rad = Math.toRadians(lon2);

    double a = SEMI_MAJOR_AXIS;
    double b = SEMI_MINOR_AXIS;
    double f = FLATTENING;

    double L = lon2Rad - lon1Rad;
    double U1 = Math.atan((1 - f) * Math.tan(lat1Rad));
    double U2 = Math.atan((1 - f) * Math.tan(lat2Rad));
    double sinU1 = Math.sin(U1), cosU1 = Math.cos(U1);
    double sinU2 = Math.sin(U2), cosU2 = Math.cos(U2);

    double lambda = L, lambdaP;
    int iterLimit = 100;
    double sinLambda, cosLambda;
    double sinSigma, cosSigma, sigma, sinAlpha, cosSqAlpha, cos2SigmaM;

    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 (sinSigma == 0) return 0; // 两点重合
        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;
        if (Double.isNaN(cos2SigmaM)) cos2SigmaM = 0; // 赤道附近特殊处理
        double C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha));
        lambdaP = lambda;
        lambda = L + (1 - C) * f * sinAlpha *
                (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)));
    } while (Math.abs(lambda - lambdaP) > 1e-12 && --iterLimit > 0);

    if (iterLimit == 0) return Double.NaN; // 迭代不收敛

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

    return s;
}

（2）替换大地线交点计算逻辑 #

放弃平面相交法，改用基于测地线参数化的数值迭代法，或直接使用成熟库（如Apache Commons Math的GeodesicLine）：

• 用Geodesic类计算每条大地线的方位角和总距离
• 参数化每条大地线的位置（比如用0到1的参数表示从起点到终点的位置）
• 通过二分法或牛顿迭代法找到两个参数，使得两条大地线的位置重合

原始问题代码 #

import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;

public class Intersection {

    // WGS84 ellipsoid constants
    private static final double SEMI_MAJOR_AXIS = 6378137.0; // meters
    private static final double FLATTENING = 1 / 298.257223563;
    private static final double SEMI_MINOR_AXIS = SEMI_MAJOR_AXIS * (1 - FLATTENING);
    private static final double ECCENTRICITY_SQUARED = 2 * FLATTENING - FLATTENING * FLATTENING;

    public static void main(String[] args) {
        // Example geodesic lines (latitude and longitude in degrees)
        double[] line1Start = {-25.0, 10.0}; 
        double[] line1End = {-28.0, 13.0};   
        double[] line2Start = {-25.0, 13.0}; 
        double[] line2End = {-28.0, 10.0};   
        int option = 0;

        double[] intersection = findIntersection(line1Start, line1End, line2Start, line2End, option);

        if (intersection != null) {
            System.out.printf("Intersection Point: Latitude = %.6f, Longitude = %.6f%n",
                    intersection[0], intersection[1]);
        } else {
            System.out.println("No intersection found or lines are coincident.");
        }
    }

    public static double[] findIntersection(double[] line1Start, double[] line1End,
            double[] line2Start, double[] line2End, int option) {
        // Convert geodetic coordinates to 3D Cartesian coordinates
        Vector3D p1 = geodeticToCartesian2(line1Start[0], line1Start[1]);
        Vector3D p2 = geodeticToCartesian2(line1End[0], line1End[1]);
        Vector3D p3 = geodeticToCartesian2(line2Start[0], line2Start[1]);
        Vector3D p4 = geodeticToCartesian2(line2End[0], line2End[1]);

        double check[] = cartesianToGeodetic(p1);
        System.out.println (String.format("Check 1: %.6f", check[0]) + String.format(" %.6f", check[1]));

        // Compute normal vectors for the planes containing the geodesics
        Vector3D n1 = p1.crossProduct(p2);
        Vector3D n2 = p3.crossProduct(p4);

        // Find the intersection line of the two planes
        Vector3D intersectionLine = n1.crossProduct(n2);
        check = cartesianToGeodetic(intersectionLine);
        System.out.println (String.format("Check 2: %.6f", check[0]) + String.format(" %.6f", check[1]));

        // Normalize the intersection line to find the intersection points on the spheroid

        // These two steps are commented out as were giving very wrong results

        // Vector3D intersectionPoint1 = intersectionLine.normalize();
        // Vector3D intersectionPoint2 = intersectionLine.negate().normalize();
        Vector3D intersectionPoint1 = intersectionLine;
        Vector3D intersectionPoint2 = intersectionLine.negate();

        // Convert back to geodetic coordinates
        double[] geodeticPoint1 = cartesianToGeodetic(intersectionPoint1);
        double[] geodeticPoint2 = cartesianToGeodetic(intersectionPoint2);
        System.out.println (String.format("GP 1: %.6f", geodeticPoint1[0]) + String.format(" %.6f", geodeticPoint1[1]));
        System.out.println (String.format("GP 2: %.6f", geodeticPoint2[0]) + String.format(" %.6f", geodeticPoint2[1]));

        // Check which intersection point lies on both geodesics
        if (isPointOnGeodesic(geodeticPoint1, line1Start, line1End) &&
                isPointOnGeodesic(geodeticPoint1, line2Start, line2End)) {
            return geodeticPoint1;
        } 
        else {
            System.out.println("Checking point 2");
            if (isPointOnGeodesic(geodeticPoint2, line1Start, line1End) &&
                    isPointOnGeodesic(geodeticPoint2, line2Start, line2End)) {
                return geodeticPoint2;
            }
        } 
        return null;
    }

    private static boolean isPointOnGeodesic(double[] gp, double[] ls, double[] le) {
        // TODO Auto-generated method stub
        double ah = haversine(ls[0], ls[1], gp[0], gp[1]);
        double bh = haversine(gp[0], gp[1], le[0], le[1]);
        double ch = haversine(ls[0], ls[1], le[0], le[1]);

        System.out.println (String.format("ah : %.2f", ah) +String.format(" bh : %.2f", bh) + String.format(" ch : %.6f", ch));
        System.out.println (String.format("ah + bh : %.6f", ah + bh) + String.format(" ch : %.6f", ch));
        if ( (ah+bh) == ch ) return true;

        return false;
    }

    private static double haversine(double lat1, double lon1, double lat2, double lon2) {

        // Convert latitude and longitude from degrees to radians
        double lat1Rad = Math.toRadians(lat1);
        double lon1Rad = Math.toRadians(lon1);
        double lat2Rad = Math.toRadians(lat2);
        double lon2Rad = Math.toRadians(lon2);

        // Haversine formula
        // Why is this not using FLATTENING? Does it matter?
        double dLat = lat2Rad - lat1Rad;
        double dLon = lon2Rad - lon1Rad;
        double a = Math.pow(Math.sin(dLat / 2), 2) +
                Math.cos(lat1Rad) * Math.cos(lat2Rad) * Math.pow(Math.sin(dLon / 2), 2);
        double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));

        // Distance in metres
        return SEMI_MAJOR_AXIS * c;
    }

    private static double[] cartesianToGeodetic(Vector3D point) {
        double x = point.getX();
        double y = point.getY();
        double z = point.getZ();
        double longitude = Math.atan2(y, x); // Longitude in radians

        double p = Math.sqrt(x * x + y * y);
        double theta = Math.atan2(z * SEMI_MAJOR_AXIS, p * SEMI_MINOR_AXIS);

        double sinTheta = Math.sin(theta);
        double cosTheta = Math.cos(theta);

        double latitude = Math.atan2(
                z + ECCENTRICITY_SQUARED * SEMI_MINOR_AXIS * sinTheta * sinTheta * sinTheta,
                p - ECCENTRICITY_SQUARED * SEMI_MAJOR_AXIS * cosTheta * cosTheta * cosTheta
                );

        double sinLatitude = Math.sin(latitude);
        double N = SEMI_MAJOR_AXIS / Math.sqrt(1 - ECCENTRICITY_SQUARED * sinLatitude * sinLatitude);

        double altitude = p / Math.cos(latitude) - N;

        // Convert radians to degrees for latitude and longitude
        double latitudeDegrees = Math.toDegrees(latitude);
        double longitudeDegrees = Math.toDegrees(longitude);

        return new double[]{latitudeDegrees, longitudeDegrees, altitude};
    }

    public static Vector3D geodeticToCartesian2(double latitude, double longitude) {
        // Convert latitude and longitude from degrees to radians
        double altitude = 0;
        double latRad = Math.toRadians(latitude);
        double lonRad = Math.toRadians(longitude);

        // Calculate the radius of curvature in the prime vertical
        double N = SEMI_MAJOR_AXIS / Math.sqrt(1 - ECCENTRICITY_SQUARED * Math.pow(Math.sin(latRad), 2));

        // Calculate Cartesian coordinates
        double x = (N + altitude) * Math.cos(latRad) * Math.cos(lonRad);
        double y = (N + altitude) * Math.cos(latRad) * Math.sin(lonRad);
        double z = ((1 - ECCENTRICITY_SQUARED) * N + altitude) * Math.sin(latRad);

        return new Vector3D(x, y, z);
    }
}

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