基于平均速度与到达时间的GPS点坐标计算及功能实现问询

问题描述 #

已知起点GPS点（A）、终点GPS点（B），目标物体（X）的平均速度（S）及预计到达B点的时间（T），需实现两个功能：

1. 估算目标物体的当前GPS位置
2. 计算物体在A到B路径上的行程占比（如中点对应50%）
   要求需考虑地球曲面特性，解决此前算法出现的不合理负值问题。

核心逻辑（球面几何） #

地球为近似球体，两点间最短路径为大圆航线，必须用球面几何公式计算，而非平面直线插值，否则会出现距离误差、位置偏移甚至负值问题。核心步骤：

1. 计算A-B大圆总距离：使用Haversine公式计算球面两点间最短距离，保证单位统一（如千米）。
2. 推导剩余行程参数：
   • 剩余时间：T - 当前时间，结果为负则表示已到达终点。
   • 剩余距离：平均速度S * 剩余时间，需处理时间为负或正的边界情况。
3. 行程占比计算：(总距离 - 剩余距离) / 总距离 * 100%，强制限制占比范围在0%-100%之间，避免异常值。
4. 当前位置插值：基于大圆航线，从A点正向（或B点反向）移动对应距离，用球面插值公式计算当前GPS坐标。

Golang实现代码 #

基础数据结构 #

import (
    "math"
    "time"
)

// GPSPoint 表示GPS坐标，单位：度（纬度lat，经度lon）
type GPSPoint struct {
    Lat float64
    Lon float64
}

// 地球半径常量，单位：千米
const earthRadiusKm = 6371.0

1. 计算两点大圆距离（Haversine公式） #

// GreatCircleDistance 计算两个GPS点之间的大圆距离，返回千米
func GreatCircleDistance(a, b GPSPoint) float64 {
    // 将度转换为弧度
    lat1 := a.Lat * math.Pi / 180
    lon1 := a.Lon * math.Pi / 180
    lat2 := b.Lat * math.Pi / 180
    lon2 := b.Lon * math.Pi / 180

    // Haversine公式核心计算
    dLat := lat2 - lat1
    dLon := lon2 - lon1
    aVal := math.Sin(dLat/2)*math.Sin(dLat/2) +
        math.Cos(lat1)*math.Cos(lat2)*
            math.Sin(dLon/2)*math.Sin(dLon/2)
    c := 2 * math.Atan2(math.Sqrt(aVal), math.Sqrt(1-aVal))

    return earthRadiusKm * c
}

2. 计算行程占比 #

// CalculateProgress 计算行程占比（0-100），输入总距离（km）、平均速度（km/h）、预计到达时间、当前时间
func CalculateProgress(totalDistanceKm, speedKmH float64, arrivalTime, currentTime time.Time) float64 {
    if speedKmH <= 0 {
        // 速度为0时，根据当前时间判断是否到达
        if currentTime.After(arrivalTime) {
            return 100.0
        }
        return 0.0
    }

    // 计算剩余时间（小时）
    remainingDuration := arrivalTime.Sub(currentTime)
    remainingHours := remainingDuration.Hours()

    // 计算剩余距离与已行驶距离
    remainingDistance := speedKmH * remainingHours
    traveledDistance := totalDistanceKm - remainingDistance

    // 限制占比在0%-100%区间内
    progress := (traveledDistance / totalDistanceKm) * 100
    switch {
    case progress < 0:
        return 0.0
    case progress > 100:
        return 100.0
    default:
        return progress
    }
}

3. 估算当前GPS位置（球面插值） #

// IntermediatePoint 计算大圆航线上，从A点出发向B点移动distanceKm后的位置
func IntermediatePoint(a, b GPSPoint, distanceKm float64) GPSPoint {
    totalDistance := GreatCircleDistance(a, b)
    if totalDistance == 0 {
        return a // 两点重合，直接返回A点
    }

    // 转换为弧度
    lat1 := a.Lat * math.Pi / 180
    lon1 := a.Lon * math.Pi / 180
    lat2 := b.Lat * math.Pi / 180
    lon2 := b.Lon * math.Pi / 180

    // 计算球面插值参数
    angularDistance := distanceKm / earthRadiusKm
    totalAngularDistance := totalDistance / earthRadiusKm

    sinTotal := math.Sin(totalAngularDistance)
    if sinTotal == 0 {
        return a
    }

    sinAngular := math.Sin(angularDistance)
    cosAngular := math.Cos(angularDistance)
    cosTotal := math.Cos(totalAngularDistance)

    // 计算目标点纬度弧度
    lat := math.Asin(math.Sin(lat1)*cosAngular + math.Cos(lat1)*sinAngular*math.Cos(lon2-lon1))

    // 计算目标点经度弧度
    lon := lon1 + math.Atan2(math.Sin(lon2-lon1)*sinAngular*math.Cos(lat1),
        cosAngular-math.Sin(lat1)*math.Sin(lat))

    // 转换回度并返回
    return GPSPoint{
        Lat: lat * 180 / math.Pi,
        Lon: lon * 180 / math.Pi,
    }
}

// EstimateCurrentPosition 估算当前位置
func EstimateCurrentPosition(a, b GPSPoint, speedKmH float64, arrivalTime, currentTime time.Time) GPSPoint {
    totalDistance := GreatCircleDistance(a, b)
    progress := CalculateProgress(totalDistance, speedKmH, arrivalTime, currentTime)

    switch {
    case progress <= 0:
        return a
    case progress >= 100:
        return b
    default:
        traveledDistance := totalDistance * (progress / 100)
        return IntermediatePoint(a, b, traveledDistance)
    }
}

异常问题说明 #

此前生成的算法出现不合理负值，主要原因：

• 未采用球面几何计算，误用平面直线插值，当两点跨经度大或纬度较高时，平面模型误差会导致插值结果超出合理范围，甚至出现负坐标。
• 未处理边界情况（如当前时间早于出发时间、晚于到达时间、速度为0等），导致逻辑触发异常计算。

使用示例 #

func main() {
    // 示例GPS点：A（北京），B（上海）
    pointA := GPSPoint{Lat: 39.9042, Lon: 116.4074}
    pointB := GPSPoint{Lat: 31.2304, Lon: 121.4737}
    speed := 80.0 // km/h
    arrivalTime := time.Date(2024, 5, 20, 18, 0, 0, 0, time.UTC)
    currentTime := time.Date(2024, 5, 20, 12, 0, 0, 0, time.UTC)

    // 计算行程占比
    progress := CalculateProgress(GreatCircleDistance(pointA, pointB), speed, arrivalTime, currentTime)
    println("行程占比：", progress, "%")

    // 估算当前位置
    currentPos := EstimateCurrentPosition(pointA, pointB, speed, arrivalTime, currentTime)
    println("当前位置：纬度", currentPos.Lat, "，经度", currentPos.Lon)
}

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