GPS轨迹压缩之Douglas-Peucker算法

创作时间:

作者:

@小白创作中心

GPS轨迹压缩之Douglas-Peucker算法

引用

来源

https://www.cnblogs.com/fangts/p/18128289

在处理大量GPS轨迹数据时，如何在保证精度的同时减少数据量是一个重要问题。Douglas-Peucker（DP）算法是一种常用的轨迹压缩算法，它通过递归地去除不重要的点来简化曲线，从而实现数据压缩。本文将详细介绍DP算法的原理、实现及效果展示。

轨迹压缩算法概述

轨迹压缩算法主要分为两大类：无损压缩和有损压缩。

无损压缩算法主要包括哈夫曼编码等方法，这类算法在压缩过程中不会丢失任何信息，但压缩率通常较低。
有损压缩算法则分为批处理方式和在线数据压缩方式。其中，批处理方式包括DP算法、TD-TR算法和Bellman算法等；在线数据压缩方式则包括滑动窗口、开放窗口、基于安全区域的方法等。

本次介绍的DP算法属于有损压缩算法中的批处理方式，它通过递归地去除不重要的点来简化曲线，从而实现数据压缩。

Douglas-Peucker算法原理

DP算法的基本步骤如下：

在轨迹曲线的首尾两点A和B之间连接一条直线AB，该直线为曲线的弦。
遍历曲线上其他所有点，求每个点到直线AB的距离，找到最大距离的点C，最大距离记为dmax。
比较该距离dmax与预先定义的阈值Dmax大小：

如果dmax < Dmax，则将该直线AB作为曲线段的近似，曲线段处理完毕。
如果dmax >= Dmax，则使C点将曲线AB分为AC和CB两段，并分别对这两段进行（1）~（3）步处理。

当所有曲线都处理完毕时，依次连接各个分割点形成的折线，即为原始曲线的路径。

显然，算法的抽稀精读是和阈值相关的，阈值越大，舍弃的点越多，提取的关键点越少。若绘制的轨迹较长，可在对地图进行缩放时动态调整阈值。

海伦公式在DP算法中的应用

DP算法中需要求点到直线的距离，该距离指的是垂直欧式距离，即直线AB外的点C到直线AB的距离d。我们采用三角形面积相等法求距离d，具体方法是：A、B、C三点构成三角形，该三角形的面积有两种求法，分别是普通方法（底x高/2）和海伦公式。

海伦公式如下：

其中p为半周长：

假设有一个三角形，边长分别为a、b、c，三角形的面积S可由以下公式求得：

我们通过海伦公式求得三角形面积，然后就可以求得高的大小，此处高即为距离d。要想用海伦公式，必须求出A、B、C三点两两之间的距离，该距离公式也是一个数学公式。

注意：求出距离后，要加上绝对值，以防止距离为负数。

平均误差与压缩率

平均误差：指的是压缩时忽略的那些点到对应线段的距离之和除以总点数得到的数值。
压缩率：压缩率的计算公式如下：

Java代码实现

为了DP算法压缩之后能匹配到本身数据库查询出的结果记录，（因为结果记录列表的每一条字段是可伸缩的KV对且不固定）准备一个点的实体。

@Data
@NoArgsConstructor
@AllArgsConstructor
@ToString
public class TYPoint {
    public String id;//点ID
    public double lng;//经度
    public double lat;//纬度
}

DP算法实现类：

public class DPUtil {
    /**
     * 默认的点到轨迹线最大距离误差阈值（单位：米）
     */
    private static double defaultDMax = 30.0;

    /**
     * DP算法入口
     * 传入压缩前的轨迹点集合
     * 输出压缩后的结果轨迹点集合
     * @param originPoints 压缩前的轨迹点集合
     * @param dMax 点到轨迹线最大距离误差阈值
     * @return
     */
    public static List<TYPoint> dpAlgorithm(List<TYPoint> originPoints, Double dMax){
        List<TYPoint> resultPoints = new ArrayList<>();
        resultPoints.add(originPoints.get(0));//获取第一个原始经纬度点坐标并添加到过滤后的数组中
        resultPoints.add(originPoints.get(originPoints.size()-1));//获取最后一个原始经纬度点坐标并添加到过滤后的数组中
        //最大距离误差阈值
        if(dMax == null){
            dMax = defaultDMax;
        }
        int start = 0;
        int end = originPoints.size()-1;
        compression(originPoints,resultPoints,start,end,dMax);
        return resultPoints;
    }

    /**
     * 函数功能：根据最大距离限制，采用DP方法递归的对原始轨迹进行采样，得到压缩后的轨迹
     * @param originPoints：原始经纬度坐标点数组
     * @param resultPoints：保持过滤后的点坐标数组
     * @param start：起始下标
     * @param end：终点下标
     * @param dMax：预先指定好的最大距离误差 计算
     */
    public static void compression(List<TYPoint> originPoints, List<TYPoint> resultPoints,
                                   int start, int end, double dMax){
        if(start < end){//递归进行的条件
            double maxDist = 0;//最大距离
            int cur_pt = 0;//当前下标
            for(int i=start+1;i<end;i++){
                //当前点到对应线段的距离
                double curDist = distToSegment(originPoints.get(start),originPoints.get(end),originPoints.get(i));
                if(curDist > maxDist){
                    maxDist = curDist;
                    cur_pt = i;
                }//求出最大距离及最大距离对应点的下标
            }
            //若当前最大距离大于最大距离误差
            if(maxDist >= dMax){
                resultPoints.add(originPoints.get(cur_pt));//将当前点加入到过滤数组中
                //将原来的线段以当前点为中心拆成两段，分别进行递归处理
                compression(originPoints,resultPoints,start,cur_pt,dMax);
                compression(originPoints,resultPoints,cur_pt,end,dMax);
            }
        }
    }

    /**
     * 函数功能：使用三角形面积（使用海伦公式求得）相等方法计算点pX到点pA和pB所确定的直线的距离
     * @param pA：起始点
     * @param pB：结束点
     * @param pX：第三个点
     * @return distance：点pX到pA和pB所在直线的距离
     */
    public static double distToSegment(TYPoint pA,TYPoint pB,TYPoint pX){
        double a = Math.abs(geoDist(pA, pB));
        double b = Math.abs(geoDist(pA, pX));
        double c = Math.abs(geoDist(pB, pX));
        double p = (a+b+c)/2.0;
        double s = Math.sqrt(Math.abs(p*(p-a)*(p-b)*(p-c)));
        double d = s*2.0/a;
        return d;
    }

    /**
     * 函数功能：根据数学公式求两个经纬度点之间的距离
     * @param pA：起始点
     * @param pB：结束点
     * @return distance：距离
     */
    public static double geoDist(TYPoint pA,TYPoint pB){
        double radLat1 = Rad(pA.lat);
        double radLat2 = Rad(pB.lat);
        double delta_lon = Rad(pB.lng - pA.lng);
        double top_1 = Math.cos(radLat2) * Math.sin(delta_lon);
        double top_2 = Math.cos(radLat1) * Math.sin(radLat2) - Math.sin(radLat1) * Math.cos(radLat2) * Math.cos(delta_lon);
        double top = Math.sqrt(top_1 * top_1 + top_2 * top_2);
        double bottom = Math.sin(radLat1) * Math.sin(radLat2) + Math.cos(radLat1) * Math.cos(radLat2) * Math.cos(delta_lon);
        double delta_sigma = Math.atan2(top, bottom);
        double distance = delta_sigma * 6378137.0;
        return distance;
    }

    /**
     * 函数功能：角度转弧度
     * @param d：角度
     * @return 返回的是弧度
     */
    public static double Rad(double d){
        return d * Math.PI / 180.0;
    }
}

从入口代码可知dMax是可以传进来的，也就是点到轨迹直线的偏移量阈值是可以设置的。这里我设置默认dMax为30.0测试了一下，4135条数据查询出来有500条，当设置dMax为100时查询出只有289条记录了。具体看业务方需要以多大的压缩限度去压缩。