【地理】地球上两点之间的最短路径的坐标（WGS84 椭圆体）附 matlab实现

原创于 2024-10-29 20:49:39 发布 · 1.1k 阅读

25 ·

CC 4.0 BY-SA版权

文章标签：

#matlab #数学建模 #算法

✅作者简介：热爱数据处理、数学建模、算法创新的Matlab仿真开发者。

🍎更多Matlab代码及仿真咨询内容点击 🔗：Matlab科研工作室

🍊个人信条：格物致知。

🔥 内容介绍

在地球表面寻找两点之间的最短路径，并非简单的直线距离计算。由于地球是一个近似于椭球体的旋转椭球体，而非完美的球体，因此最短路径并非欧几里得空间中的直线，而是沿着地球表面曲率变化而形成的测地线 (Geodesic)。WGS84椭圆体是目前全球广泛使用的地球椭球体模型，其参数精确地定义了地球形状和大小，因此在精确计算地球上两点间最短路径时，必须考虑WGS84椭圆体的特性。本文将深入探讨在WGS84椭圆体上计算两点之间最短路径坐标的方法，并提供相应的Matlab实现代码。

首先，需要明确的是，在WGS84椭圆体上计算两点之间最短路径，通常采用Vincenty公式或其改进算法。Vincenty公式是一个迭代算法，它通过一系列复杂的三角函数计算来逼近测地线距离和方位角。该公式考虑了地球椭球体的扁率，能够提供高精度的计算结果。相比之下，简单的球面距离计算公式（例如Haversine公式）忽略了地球椭球体的扁率，其精度较低，尤其是在高纬度地区误差较大，因此不适用于需要高精度的应用场景。

Vincenty公式涉及一系列中间变量的计算，包括：

经纬度坐标转换: 将地理坐标（经度λ和纬度φ）转换为大地坐标 (x, y, z)，以便进行椭球体上的计算。
正算 (Direct Problem): 已知一点的经纬度坐标和方位角以及距离，计算另一点的经纬度坐标。
反算 (Inverse Problem): 已知两点的经纬度坐标，计算两点之间的距离和方位角。

反算问题是本文的核心问题，即已知两点经纬度坐标，求解最短路径以及路径上的中间点坐标。Vincenty公式的反算过程是一个迭代过程，其迭代公式较为复杂，在此不赘述具体推导过程，但其核心思想是通过不断逼近来找到最短路径的方位角和距离。

lambdaP = lambda;
sinSigma = sqrt((cos(U2) * sin(lambda))^2 + (cos(U1) * sin(U2) - sin(U1) * cos(U2) * cos(lambda))^2);
cosSigma = sin(U1) * sin(U2) + cos(U1) * cos(U2) * cos(lambda);
sigma = atan2(sinSigma, cosSigma);
sinAlpha = cos(U1) * cos(U2) * sin(lambda) / sinSigma;
alpha = asin(sinAlpha);
cos2SigmaM = cos(sigma) - 2 * sin(U1) * sin(U2) / (cos(alpha)^2);
C = f / 16 * cos(alpha)^2 * (4 + f * (4 - 3 * cos(alpha)^2));
lambda = L + (1 - C) * f * sin(alpha) * (sigma + C * sin(sigma) * (cos2SigmaM + C * cos(sigma) * (-1 + 2 * cos2SigmaM^2)));
end

uSq = cos(alpha)^2 * (a^2 - b^2) / b^2;
A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)));
B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)));
deltaSigma = B * sin(sigma) * (cos2SigmaM + B / 4 * (cos(sigma) * (-1 + 2 * cos2SigmaM^2) - B / 6 * cos2SigmaM * (-3 + 4 * sinSigma^2) * (-3 + 4 * cos2SigmaM^2)));
distance = b * A * (sigma - deltaSigma);

azimuth1 = atan2(cos(U2) * sin(lambda), cos(U1) * sin(U2) - sin(U1) * cos(U2) * cos(lambda));
azimuth2 = atan2(cos(U1) * sin(lambda), -sin(U1) * cos(U2) + cos(U1) * sin(U2) * cos(lambda));

azimuth1 = rad2deg(azimuth1);
azimuth2 = rad2deg(azimuth2);
distance = distance;

end

该函数计算两点之间的距离 (distance, 米) 和方位角 (azimuth1, azimuth2, 度)。为了计算中间点坐标，需要利用Vincenty公式的正算部分，该部分需要根据距离和方位角迭代计算中间点的坐标。完整的中间点坐标计算需要更复杂的代码，但其核心思想是重复应用Vincenty正算公式。

总结而言，精确计算WGS84椭圆体上两点之间的最短路径需要采用Vincenty公式或其改进算法。本文提供了Matlab实现的Vincenty公式反算部分，能够计算两点之间的距离和方位角。完整的路径计算需要结合正算部分，并进行多次迭代计算，以得到路径上所有中间点的坐标。需要注意的是，Vincenty公式虽然精度较高，但在极端情况下（例如两点非常接近极点）仍然可能存在收敛性问题，需要进行相应的误差处理和算法优化。未来的研究可以关注更高效、更鲁棒的测地线计算算法，以及如何将这些算法应用于更复杂的地理空间分析任务中。