当前位置：首页 > news >正文

椭球面长度计算的两种公式及投影选择

news 2025/5/8 18:00:44

最近几年，涉及椭球面长度计算逐渐增多，它的主要目的是用来验算投影参数选取是否合理。

椭球面长度计算的公式有两种，分别是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作为投影是适宜的。

点击查看原图