基于步行速度，在两个GPS位置之间进行插值

Question

问题：

---

给定两个位置：

L₁ = (纬度₁, 经度₁, 时间戳₁), L₂ = (纬度₂, 经度₂, 时间戳₂),

以及可配置但恒定的移动速度：

v = 每秒1.39米（例如）。

在用户从L₁到L₂行驶的过程中，我们如何在这两个位置之间进行插值，以估算出用户的位置？

---

我一直在寻找解决这个问题的方法，到目前为止，我发现对于较小距离（远离极地），可以使用线性插值。因此，我在维基百科上查找了线性插值，并找到了以下内容：

// Imprecise method which does not guarantee v = v1 when t = 1,
// due to floating-point arithmetic error.
float lerp(float v0, float v1, float t) {
    return v0 + t*(v1-v0);
}

我想使用这个lerp函数在L₁和L₂之间插值纬度和经度。那是比较简单的部分。如何计算t？我猜我必须计算一些时间差，但是如何考虑移动速度呢？

---

编辑：

我正在测试各种收集GPS位置的方法。为此，我在整个步行过程中记录航点位置。我需要使用移动速度在这些航点之间进行插值，以估计沿途的位置。然后，我可以将我的结果与估计结果进行比较，查看它们的准确程度。

例如：

- transporter_room_3

---

所以您的意思是设备在固定时间内从A点移动到B点，但它们的总速度可能完全不相关？ - njzk2

我已经实施了不同的GPS位置收集策略，例如更新之间的最小距离、更新之间的最短时间、启用加速计的策略等等。我需要建立一个“基准真相”来比较结果。我通过在实验过程中收集单独的航点来实现这一点，并且需要根据这些航点之间的插值（使用移动速度），以便我可以将它们与其他结果进行比较。我希望这讲得清楚，否则我会尝试更好地解释。 - transporter_room_3

在你的问题中存在一个问题：你说你想要一个恒定可配置的速度v。然而，你测量了时间戳和位置，这隐含地代表了一个速度。如果一个人以10公里/小时的速度跑步，你不能用一个给定的错误速度来插值。我不明白为什么你需要速度。你有时间戳作为插值的参考。你不能同时进行两者的插值。决定你的参考是什么；要么是(时间和位置)对，要么是(速度和位置)对，后者更复杂一些。 - AlexWien

好观点，AlexWien。我刚刚也意识到了这一点。我得想出采用哪种方法。我认为前者最符合我的需求。 - transporter_room_3

前一种方法在我的回答中可以看到。你是为自己、应用程序、大学还是商业目的做这个任务的？（获得报酬） - AlexWien

我正在为大学课程开发这个应用程序。航点插值是该应用程序中的一个小但重要的部分。感谢您提供的帮助。我实现了自己的版本，并且它运行得非常好。谢谢！ - transporter_room_3

Answer 1

仔细查看计算纬度/经度点之间距离、方位和更多内容，其中包含几个公式和JavaScript示例，这可能会对您有所帮助。我知道它不是Java，但应该足够简单，可以将代码移植过来。尤其是给出了公式的详细说明。

编辑：

虽然在较短距离内使用线性插值似乎没问题，但实际上可能相当偏差，特别是当您靠近极地时。从您在汉堡的例子中可以看出，这将在数百米内就会有明显的影响。请参见此答案以获得良好的解释。

问题：经度1度之间的距离取决于纬度。

这是因为地球不是平面而是一个球体 - 实际上是一个椭球体。因此，在二维地图上的直线在地球上 不是 直线 - 反之亦然。

为了解决这个问题，我们可以采用以下方法：

1. 从起始坐标（L1）指向结束坐标（L2）获取方位角
2. 沿着大圆路径从起始坐标（L1）计算一个新的坐标，给定计算出的方位角和指定的距离
3. 重复此过程，但使用新计算出的坐标作为起始坐标

我们可以创建一些简单的函数来完成这个技巧：

double radius = 6371; // earth's mean radius in km

// Helper function to convert degrees to radians
double DegToRad(double deg) {
    return (deg * Math.PI / 180);
}

// Helper function to convert radians to degrees
double RadToDeg(double rad) {
    return (rad * 180 / Math.PI);
}

// Calculate the (initial) bearing between two points, in degrees
double CalculateBearing(Location startPoint, Location endPoint) {
    double lat1 = DegToRad(startPoint.latitude);
    double lat2 = DegToRad(endPoint.latitude);
    double deltaLon = DegToRad(endPoint.longitude - startPoint.longitude);

    double y = Math.sin(deltaLon) * Math.cos(lat2);
    double x = Math.cos(lat1) * Math.sin(lat2) - Math.sin(lat1) * Math.cos(lat2) * Math.cos(deltaLon);
    double bearing = Math.atan2(y, x);

    // since atan2 returns a value between -180 and +180, we need to convert it to 0 - 360 degrees
    return (RadToDeg(bearing) + 360) % 360;
}

// Calculate the destination point from given point having travelled the given distance (in km), on the given initial bearing (bearing may vary before destination is reached)
Location CalculateDestinationLocation(Location point, double bearing, double distance) {

    distance = distance / radius; // convert to angular distance in radians
    bearing = DegToRad(bearing); // convert bearing in degrees to radians

    double lat1 = DegToRad(point.latitude);
    double lon1 = DegToRad(point.longitude);

    double lat2 = Math.asin(Math.sin(lat1) * Math.cos(distance) + Math.cos(lat1) * Math.sin(distance) * Math.cos(bearing));
    double lon2 = lon1 + Math.atan2(Math.sin(bearing) * Math.sin(distance) * Math.cos(lat1), Math.cos(distance) - Math.sin(lat1) * Math.sin(lat2));
    lon2 = (lon2 + 3 * Math.PI) % (2 * Math.PI) - Math.PI; // normalize to -180 - + 180 degrees

    return new Location(RadToDeg(lat2), RadToDeg(lon2));
}

// Calculate the distance between two points in km
double CalculateDistanceBetweenLocations(Location startPoint, Location endPoint) {

    double lat1 = DegToRad(startPoint.latitude);
    double lon1 = DegToRad(startPoint.longitude);

    double lat2 = DegToRad(endPoint.latitude);
    double lon2 = DegToRad(endPoint.longitude);

    double deltaLat = lat2 - lat1;
    double deltaLon = lon2 - lon1;

    double a = Math.sin(deltaLat / 2) * Math.sin(deltaLat / 2) + Math.cos(lat1) * Math.cos(lat2) * Math.sin(deltaLon / 2) * Math.sin(deltaLon / 2);
    double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));

    return (radius * c);
}

使用平均地球半径为6371公里。有关这个数字及其准确性的解释，请参见维基百科。

现在，可以通过给定的行程距离（以公里为单位）计算两点之间的新中间位置：

double bearing = CalculateBearing(startLocation, endLocation);

Location intermediaryLocation = CalculateDestinationLocation(startLocation, bearing, distanceTravelled);

假设速度为v（例如1.39）米/秒，现在可以使用一个简单的for循环来获取相隔1秒的点：

List<Location> locations = new ArrayList<Location>();

// assuming duration in full seconds
for (int i = 0; i < duration; i++){
    double bearing = CalculateBearing(startLocation, endLocation);
    double distanceInKm = v / 1000;
    Location intermediaryLocation = CalculateDestinationLocation(startLocation, bearing, distanceInKm);

    // add intermediary location to list
    locations.add(intermediaryLocation);

    // set intermediary location as new starting location
    startLocation = intermediaryLocation;
}

作为额外的奖励，您甚至可以确定在任意两点之间旅行所需的时间：

double distanceBetweenPoints = CalculateDistanceBetweenLocations(startPoint, endPoint) * 1000; // multiply by 1000 to get meters instead of km

double timeRequired = distanceBetweenPoints / v;

使用这种方法比仅使用坐标差的简单线性插值在任何距离上都能得到更高的精度。虽然该方法并不完美，但通常会有0.3%或更低的误差，这是相当可接受的。如果您需要更好的解决方案，可以考虑使用Vincenty公式。

- davethebrave

---

我没有看到任何回答他问题的内容。 - AlexWien

仅包含链接的答案是不好的答案：http://meta.stackexchange.com/questions/8231/are-answers-that-just-contain-links-elsewhere-really-good-answers - Juan Lopes

这是一个非常有帮助的答案！我唯一的注释是在你的方法CalculateDistanceBetweenLocations中，你有double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - 1));。这应该改为：double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));。 - user3798780

我已经提交了一份编辑，以更正用户3798780报告的CalculateDistanceBetweenLocations错误。原始帖子现在应该正确显示Math.sqrt(1 - a)。 - hsoi

很好的回答，但是能否编辑原帖以正确拼写“longitude”？在CalculateDestinationLocation函数的第7行中，它被拼写为“logintude”。 - lolololol ol

Answer 2

如果你先将你的纬度/经度转换为n-向量 (https://en.wikipedia.org/wiki/N-vector)，那么像这样的计算实际上非常简单。 转换后，您可以使用标准插值，并且还将避免任何长距离、极点或日期变更线方面的问题。

如果你在维基百科页面上检查“外部链接”，会发现有一个网页 (http://www.navlab.net/nvector/) 解决了十个问题，该页面的第六个问题 (插值位置) 应该与您的问题相同。 如您所见，该解决方案对于任何距离都是精确的，而且也适用于任何地球位置，如极点。

- MathGuy

Answer 3

我猜我需要计算一些时间差，但是如何考虑移动速度？

在线性插值中，在您的情况下，您会在两个时间点之间迭代，使用迭代变量t，该变量从开始时间t1到结束时间t2运行，并具有预定义的步骤。 假设步长为1秒，这对于您的应用程序非常实用。

long t1 = location1.getTimeStamp(); // in milliseconds;
long t2 = location2.getTimeStamp();
double deltaLat = location2.latitude - location1.latitude;
doule deltaLon =  location2.longitude- location1.longtude;
// remove this line if you don't have measured speed:
double deltaSpeed =  location2.speed - location1.speed;

long step = 1 * 1000; // 1 second in millis 
for (long t = t1; t1 < t2; t+= step) {

   // t0_1 shall run from 0.0 to (nearly) 1.0 in that loop
  double t0_1 = (t - t1) / (t2 - t1);
  double latInter = lat1 + deltaLat  * t0_1;
  double lonInter = lon1 + deltaLon  * t0_1;
  // remove the line below if you dont have speed
  double speedInter = speed1 + deltaSpeed  * t0_1;
  Location interPolLocation = new Location(latInter, lonInter, speedInter);
  // add interPolLocation to list or plot.
}

- AlexWien

---

这段代码假设纬度和经度在一个平面直角坐标系上。但实际上地球并不是这样，地球是一个球体。 - Ali Ashraf

Answer 4

有一些比线性插值更好的插值策略，包括运动插值，它以锚点的初始和最终速度作为输入。例如，参见最近一篇论文中的比较（Long JA (2015) Kinematic interpolation of movement data. Int J Geogr Inf Sci 8816:1–15. doi: 10.1080/13658816.2015.1081909）：

对于运动插值，有R和Python的实现。编写Java版本应该很容易。

- Tales P. Nogueira

Answer 5

没错。线性插值。

L1 = (1, 2, 3)
L2 = (4, 5, 6)
desired_number_of_interpolation_points = 9
interpolation_points = []

lat_step = (L2[0] - L1[0]) / (desired_number_of_interpolation_points + 1)
lon_step = (L2[1] - L1[1]) / (desired_number_of_interpolation_points + 1)
time_step = (L2[2] - L1[2]) / (desired_number_of_interpolation_points + 1) 

for i in range(1, desired_number_of_interpolation_points + 1)
    interpolation_points.append((lat_step * i, lon_step * i, time_step * i))

- Chad Kennedy