最近几年，涉及椭球面长度计算逐渐增多，它的主要目的是用来验算投影参数选取是否合理。
椭球面长度计算的公式有两种，分别是Haversine(半正矢公式)和Vincenty公式。
Haversine(半正矢公式)是一种近似公式（近似为球），公式如下：
Haversine公式比较简单，可在Excel的单元格中编辑。
Vincenty公式是由Thaddeus Vincenty于1975年提出的一种用于计算地球表面两点间最短距离（即大圆距离）的迭代算法。它考虑了地球的椭球形状，因此比基于正球体模型的Haversine公式更为精确。Vincenty公式特别适用于需要高精度距离计算的场景，如航空导航、远洋航行和精密地图制作。
Python程序如下：
import math
def vincenty_formula(lat1, lon1, lat2, lon2):
# 地球椭球参数 (WGS-84)
a = 6378137.0 # 半长轴 (meters)
f = 1 / 298.257223563 # 扁率
b = (1 - f) * a# 将经纬度转换为弧度
lat1, lon1, lat2, lon2 = map(math.radians, [lat1, lon1, lat2, lon2])
# 计算U值
U1 = math.atan((1 - f) * math.tan(lat1))
U2 = math.atan((1 - f) * math.tan(lat2))
sinU1 = math.sin(U1)
cosU1 = math.cos(U1)
sinU2 = math.sin(U2)
cosU2 = math.cos(U2)
# 初始化迭代变量
L = lon2 - lon1
lambda_ = L
sin_lambda = cos_lambda = sin_sigma = cos_sigma = sigma = None
sin_alpha = cos_sq_alpha = cos2_sigma_m = C = None
iter_limit = 100
for _ in range(iter_limit):
sin_lambda = math.sin(lambda_)
cos_lambda = math.cos(lambda_)
sin_sigma = math.sqrt((cosU2 * sin_lambda) ** 2 +
(cosU1 * sinU2 - sinU1 * cosU2 * cos_lambda) ** 2)
cos_sigma = sinU1 * sinU2 + cosU1 * cosU2 * cos_lambda
sigma = math.atan2(sin_sigma, cos_sigma)
sin_alpha = cosU1 * cosU2 * sin_lambda / sin_sigma
cos_sq_alpha = 1 - sin_alpha ** 2
cos2_sigma_m = cos_sigma - 2 * sinU1 * sinU2 / cos_sq_alpha
C = f / 16 * cos_sq_alpha * (4 + f * (4 - 3 * cos_sq_alpha))
lambda_prev = lambda_
lambda_ = L + (1 - C) * f * sin_alpha * (
sigma + C * sin_sigma * (cos2_sigma_m + C * cos_sigma * (-1 + 2 * cos2_sigma_m ** 2))
)
if abs(lambda_ - lambda_prev) < 1e-12:
break
else:
raise ValueError("Vincenty formula failed to converge")
u_sq = cos_sq_alpha * (a ** 2 - b ** 2) / (b ** 2)
A = 1 + u_sq / 16384 * (4096 + u_sq * (-768 + u_sq * (320 - 175 * u_sq)))
B = u_sq / 1024 * (256 + u_sq * (-128 + u_sq * (74 - 47 * u_sq)))
delta_sigma = B * sin_sigma * (cos2_sigma_m + B / 4 * (
cos_sigma * (-1 + 2 * cos2_sigma_m ** 2) - B / 6 * cos2_sigma_m * (-3 + 4 * sin_sigma ** 2) * (-3 + 4 * cos2_sigma_m ** 2)
))
s = b * A * (sigma - delta_sigma)
将之改写为VBA程序代码如下：
Public Function distance(ByVal lon1 As Double, ByVal lat1 As Double, ByVal lon2 As Double, ByVal lat2 As Double) As Double
'地球椭球参数 (WGS-84)
Dim a As Double, f As Double, b As Double
a = 6378137# ' 半长轴 (meters)
f = 1 / 298.257223563 ' 扁率
b = (1 - f) * a
'将经纬度转换为弧度
lat1 = lat1 * Pi() / 180#
lon1 = lon1 * Pi() / 180#
lat2 = lat2 * Pi() / 180#
lon2 = lon2 * Pi() / 180#
' 计算U值
Dim U1 As Double, U2 As Double, sinU1 As Double, cosU1 As Double, sinU2 As Double, cosU2 As Double
U1 = Atn((1 - f) * Tan(lat1))
U2 = Atn((1 - f) * Tan(lat2))
sinU1 = Sin(U1)
cosU1 = Cos(U1)
sinU2 = Sin(U2)
cosU2 = Cos(U2)
'初始化迭代变量
Dim L As Double, lambda_ As Double, sin_lambda As Double, sin_alpha As Double
Dim sigma As Double, cos_sigma As Double, sin_sigma As Double, cos_lambda As Double
Dim C As Double, cos2_sigma_m As Double, cos_sq_alpha As Double
Dim iter_limit As Integer, i As Integer
L = lon2 - lon1
lambda_ = L
iter_limit = 100
For i = 1 To iter_limit
sin_lambda = Sin(lambda_)
cos_lambda = Cos(lambda_)
sin_sigma = Sqr((cosU2 * sin_lambda) ^ 2 + (cosU1 * sinU2 - sinU1 * cosU2 * cos_lambda) ^ 2)
cos_sigma = sinU1 * sinU2 + cosU1 * cosU2 * cos_lambda
sigma = Atn2(sin_sigma, cos_sigma)
sin_alpha = cosU1 * cosU2 * sin_lambda / sin_sigma
cos_sq_alpha = 1 - sin_alpha ^ 2
cos2_sigma_m = cos_sigma - 2 * sinU1 * sinU2 / cos_sq_alpha
C = f / 16 * cos_sq_alpha * (4 + f * (4 - 3 * cos_sq_alpha))
lambda_prev = lambda_
lambda_ = cos2_sigma_m + C * cos_sigma * (-1 + 2 * cos2_sigma_m ^ 2)
lambda_ = sigma + C * sin_sigma * (lambda_)
lambda_ = L + (1 - C) * f * sin_alpha * (lambda_)
If Abs(lambda_ - lambda_prev) < 0.000000000001 Then
Exit For
End If
Next
Dim u_sq As Double, A1 As Double, B1 As Double, delta_sigma As Double, s As Double
u_sq = cos_sq_alpha * (a ^ 2 - b ^ 2) / (b ^ 2)
A1 = 4096 + u_sq * (-768 + u_sq * (320 - 175 * u_sq))
A1 = 1 + u_sq / 16384 * (A1)
B1 = 256 + u_sq * (-128 + u_sq * (74 - 47 * u_sq))
B1 = u_sq / 1024 * (B1)
delta_sigma = cos_sigma * (-1 + 2 * cos2_sigma_m ^ 2) - B1 / 6 * cos2_sigma_m * (-3 + 4 * sin_sigma ^ 2) * (-3 + 4 * cos2_sigma_m ^ 2)
delta_sigma = B1 * sin_sigma * (cos2_sigma_m + B1 / 4 * (delta_sigma))
s = b * A1 * (sigma - delta_sigma)
distance = s
End Function
下表是某线段的投影公式（UTM123）、Haversine公式和Vincenty公式计算的长度，从表中可以看出，线段投影（UTM123）长度487787m，Haversine公式计算的长度为488200m，Vincenty公式计算的长度为487884m，后两者比前者分别长412m、96m。以Vincenty公式计算的长度为准，它与投影长度的差值相对于总长的比例为0.2‰<1‰，因此，选择UTM123作为投影是适宜的。
参考文献
https://www.chimaozy.com/vincenty%E5%85%AC%E5%BC%8F/
https://www.oryoy.com/news/python-shi-xian-vincenty-suan-fa-jing-que-ding-wei-liang-dian-jian-ju-li-ji-suan-xiang-jie.html
https://www.cnblogs.com/yzyeilin/archive/2013/03/16/2962604.html
https://blog.csdn.net/sinat_36912383/article/details/132561815
https://www.cnblogs.com/aoldman/p/4241117.html
https://zhuanlan.zhihu.com/p/373411796
发表评论：