WGS84坐标系下两点穿越圆形浮标的时间估算Python实现方法

核心逻辑 #

已知两点分别位于浮标内外，假设这段间隔内运动为大圆弧匀速运动，通过以下步骤估算穿越时间：

• 计算两点与浮标中心的球面距离，确认内外位置
• 计算两点间的大圆弧距离、运动速度与方位角
• 求解运动轨迹与浮标边界的交点，根据弧长比例插值得到穿越时间

1. 半正矢距离计算函数 #

实现haversine公式，计算WGS84坐标系下两点的球面距离（单位：米），同时用于判断点是否在浮标内部：

import math
from datetime import datetime

def haversine(lat1, lon1, lat2, lon2):
    # 经纬度转弧度
    lat1_rad = math.radians(lat1)
    lon1_rad = math.radians(lon1)
    lat2_rad = math.radians(lat2)
    lon2_rad = math.radians(lon2)
    
    d_lat = lat2_rad - lat1_rad
    d_lon = lon2_rad - lon1_rad
    
    # 半正矢公式计算球面距离
    a = math.sin(d_lat/2)**2 + math.cos(lat1_rad) * math.cos(lat2_rad) * math.sin(d_lon/2)**2
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))
    earth_radius = 6371000  # 地球平均半径，单位米
    return earth_radius * c

2. 判断点的内外位置 #

用上述函数计算两点到浮标中心的距离，确认哪个点在浮标内部、哪个在外部：

# 已知参数
lat1, lon1, t1 = 45.965467, 8.509283, datetime(2024, 6, 19, 11, 46, 39)  # P1
lat2, lon2, t2 = 45.968483, 8.5089, datetime(2024, 6, 19, 11, 57, 12)    # P2
lat_c, lon_c, r = 45.987168002, 8.452897189, 4500                        # 浮标中心坐标与半径

# 计算两点到浮标中心的距离
dist_p1_c = haversine(lat1, lon1, lat_c, lon_c)
dist_p2_c = haversine(lat2, lon2, lat_c, lon_c)

# 标记内外状态
is_p1_inside = dist_p1_c <= r
is_p2_inside = dist_p2_c <= r

3. 计算运动轨迹参数 #

计算两点间的大圆弧距离、平均运动速度，以及从起点到终点的方位角：

# 两点间总球面距离
dist_p1_p2 = haversine(lat1, lon1, lat2, lon2)
# 时间差（单位：秒）
time_diff = (t2 - t1).total_seconds()
# 平均运动速度（单位：米/秒）
speed = dist_p1_p2 / time_diff if time_diff != 0 else 0

# 计算两点间的方位角（0-360度）
def calculate_bearing(lat1, lon1, lat2, lon2):
    lat1_rad = math.radians(lat1)
    lon1_rad = math.radians(lon1)
    lat2_rad = math.radians(lat2)
    lon2_rad = math.radians(lon2)
    
    d_lon = lon2_rad - lon1_rad
    y = math.sin(d_lon) * math.cos(lat2_rad)
    x = math.cos(lat1_rad)*math.sin(lat2_rad) - math.sin(lat1_rad)*math.cos(lat2_rad)*math.cos(d_lon)
    bearing = math.atan2(y, x)
    bearing = math.degrees(bearing)
    return (bearing + 360) % 360  # 转为0-360度范围

bearing_p1_p2 = calculate_bearing(lat1, lon1, lat2, lon2)

4. 求解轨迹与浮标边界的交点并估算穿越时间 #

利用球面余弦定理，计算从起点到交点的弧长，再根据速度比例插值得到穿越时间：

def estimate_crossing_time(lat_start, lon_start, t_start, lat_end, lon_end, t_end, lat_c, lon_c, r):
    dist_start_c = haversine(lat_start, lon_start, lat_c, lon_c)
    dist_end_c = haversine(lat_end, lon_end, lat_c, lon_c)
    dist_start_end = haversine(lat_start, lon_start, lat_end, lon_end)
    
    # 过滤两点同在内/同在外的无效情况
    if (dist_start_c <= r and dist_end_c <= r) or (dist_start_c > r and dist_end_c > r):
        return None
    
    # 计算起点到终点的方位角与起点到浮标中心方位角的夹角
    def get_bearing_diff(lat_a, lon_a, lat_b, lon_b, lat_c, lon_c):
        bearing_a_b = calculate_bearing(lat_a, lon_a, lat_b, lon_b)
        bearing_a_c = calculate_bearing(lat_a, lon_a, lat_c, lon_c)
        diff = abs(bearing_a_b - bearing_a_c)
        return min(diff, 360 - diff)  # 取最小夹角
    
    theta_deg = get_bearing_diff(lat_start, lon_start, lat_end, lon_end, lat_c, lon_c)
    theta_rad = math.radians(theta_deg)
    
    # 球面三角形余弦定理计算起点到交点的弧长（转为弧度计算）
    r_rad = r / 6371000
    dist_start_c_rad = dist_start_c / 6371000
    cross_dist_rad = math.acos(
        math.cos(r_rad) - math.cos(dist_start_c_rad)*math.cos(r_rad) + 
        math.cos(dist_start_c_rad)*math.cos(r_rad)*math.cos(theta_rad)
    )
    dist_start_crossing = cross_dist_rad * 6371000
    
    # 根据弧长比例插值计算穿越时间
    time_ratio = dist_start_crossing / dist_start_end
    total_time_diff = (t_end - t_start).total_seconds()
    crossing_time = t_start + datetime.timedelta(seconds=time_ratio * total_time_diff)
    
    return crossing_time

# 根据内外点确定计算方向（从外部点指向内部点）
if is_p1_inside:
    crossing_time = estimate_crossing_time(lat2, lon2, t2, lat1, lon1, t1, lat_c, lon_c, r)
else:
    crossing_time = estimate_crossing_time(lat1, lon1, t1, lat2, lon2, t2, lat_c, lon_c, r)

print(f"估算的穿越时间：{crossing_time}")

关键说明 #

• 假设运动为匀速大圆弧运动，这是GPS信号丢失间隔内的合理近似
• 半正矢公式确保了WGS84坐标系下的球面距离计算精度，避免平面近似误差
• 球面几何求解交点的逻辑，适配大半径浮标的场景，计算结果更贴合实际运动轨迹

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