轨迹压缩之Douglas-Peucker算法之Java实现

第一部分 问题描述

1.1 具体任务

本次作业任务是轨迹压缩，给定一个GPS数据记录文件，每条记录包含经度和维度两个坐标字段，所有记录的经纬度坐标构成一条轨迹，要求采用合适的压缩算法，使得压缩后轨迹的距离误差小于30m。

1.2 程序输入

本程序输入是一个GPS数据记录文件。

1.3 数据输出

输出形式是文件，包括三部分，压缩后点的ID序列及坐标、点的个数、平均距离误差、压缩率

第二部分 问题解答

根据问题描述，我们对问题进行求解，问题求解分为以下几步：

2.1 数据预处理

本次程序输入为GPS数据记录文件，共有3150行记录，每行记录又分为若干个字段，根据题意，我们只需关注经度和纬度坐标字段即可，原始数据文件部分记录如图2.1所示：

图2.1 原始数据文件部分记录示意图

如图2.1所示，原始数据文件每条记录中经纬度坐标字段数据的保存格式是典型的GPS坐标表达方式，即度分格式，形式为dddmm.mmmm，其中ddd表示度，mm.mmmm表示分，小数点前面表示分的整数部分，小数点后表示分的小数部分；本次数据预处理，为方便后面两个坐标点之间距离的计算，我们需要将度分格式的经纬度坐标数据换算成度的形式，换算方法是ddd+mm.mmmm/60，此处我们保留小数点后6位数字，换算后的形式为ddd.xxxxxx。

我们以第一条记录中经纬度坐标（11628.2491，3955.6535）为例，换算后的结果为（116.470818，39.927558），所有记录中经纬度坐标都使用方法进行，并且可以为每一个转换后的坐标点生成一个ID，进行唯一标识，压缩后，我们只需输出所有被保留点的ID即可。

2.2 Douglas-Peucker轨迹压缩算法

轨迹压缩算法分为两大类，分别是无损压缩和有损压缩，无损压缩算法主要包括哈夫曼编码，有损压缩算法又分为批处理方式和在线数据压缩方式，其中批处理方式又包括DP（Douglas-Peucker）算法、TD-TR（Top-Down Time-Ratio）算法和Bellman算法，在线数据压缩方式又包括滑动窗口、开放窗口、基于安全区域的方法等。

由于时间有限，本次轨迹压缩，我们决定采用相对简单的DP算法。

DP算法步骤如下：

（1）在轨迹曲线在曲线首尾两点A，B之间连接一条直线AB，该直线为曲线的弦；

（2）遍历曲线上其他所有点，求每个点到直线AB的距离，找到最大距离的点C，最大距离记为dmax；

（3）比较该距离dmax与预先定义的阈值Dmax大小，如果dmax<Dmax，则将该直线AB作为曲线段的近似，曲线段处理完毕；

（4）若dmax>=Dmax，则使C点将曲线AB分为AC和CB两段，并分别对这两段进行（1）~（3）步处理；

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

2.3 点到直线的距离

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

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

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

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

2.4 平均误差求解

平均误差指的是压缩时忽略的那些点到对应线段的距离之和除以总点数得到的数值。

2.5 压缩率求解

压缩率的计算公式如下：

2.6 数据结果文件的生成

经过上面的处理和计算，我们将压缩后剩余点的ID和点的个数、平均距离误差、压缩率等参数都写入最终的结果文件中，问题解答完成。

第三部分 代码实现

本程序采用Java语言编写，开发环境为IntelliJ IDEA 14.0.2，代码共分为两个类，一个是ENPoint类，用于保存经纬度点信息，一个是TrajectoryCompressionMain类，用于编写数据处理、DP算法、点到直线距离、求平均误差等函数。

3.1 程序总流程

整个程序流程主要包括以下几个步骤：

（1）定义相关ArrayList数组和File对象，其中ArrayList数组对象有三个，分别是原始经纬度坐标数组pGPSArryInit、过滤后的点坐标数组pGPSArrayFilter、过滤并排序后的点坐标数组pGPSArrayFilterSort；File文件对象共有五个，分别是原始数据文件对象fGPS、压缩后的结果数据文件对象oGPS、保持转换后的原始经纬度坐标点的数据文件fInitGPSPoint、仿真测试文件fTestInitPoint和fTestFilterPoint。

（2）获取原始点坐标并将其写入到文件中，主要包括读文件和写文件两种操作；

（3）进行轨迹压缩；

（4）对压缩后的经纬度点坐标进行排序；

（5）生成仿真测试文件，并用R语言工具进行图形绘制，得到最终的结果；

（6）求平均误差和压缩率，平均误差通过函数求得，压缩率直接计算获得；

（7）将最终结果写入结果文件中，包括过滤后的点的ID，点的个数、平均误差和压缩率；

3.2 具体实现代码

（1）ENPoint.java

 1 package cc.xidian.main;   2    3 import java.text.DecimalFormat;   4    5 /**   6  * Created by hadoop on 2015/12/20.   7  */   8 public class ENPoint implements Comparable<ENPoint>{   9     public int id;//点ID  10     public double pe;//经度  11     public double pn;//维度  12   13     public ENPoint(){}//空构造函数  14     public String toString(){  15         //DecimalFormat df = new DecimalFormat("0.000000");  16         return this.id+"#"+this.pn+","+this.pe;  17     }  18     public String getTestString(){  19         DecimalFormat df = new DecimalFormat("0.000000");  20         return df.format(this.pn)+","+df.format(this.pe);  21     }  22     public String getResultString(){  23         DecimalFormat df = new DecimalFormat("0.000000");  24         return this.id+"#"+df.format(this.pn)+","+df.format(this.pe);  25     }  26     @Override  27     public int compareTo(ENPoint other) {  28         if(this.id<other.id)  return -1;  29         else if(this.id>other.id)  return 1;  30         else  31             return 0;  32     }  33 }

（2）TrajectoryCompressionMain.java

  1 package cc.xidian.main;    2     3 import java.io.*;    4 import java.text.DecimalFormat;    5 import java.util.*;    6 import java.util.List;    7     8 /**    9  * Created by hadoop on 2015/12/19.   10  */   11 public class TrajectoryCompressionMain{   12     public static void main(String[] args)throws Exception{   13         //-----------------------1、相关ArrayList数组和File对象的声明和定义-------------------------------------------------//   14         ArrayList<ENPoint> pGPSArrayInit = new ArrayList<ENPoint>();//原纪录经纬度坐标数组   15         ArrayList<ENPoint> pGPSArrayFilter = new ArrayList<ENPoint>();//过滤后的经纬度坐标数组   16         ArrayList<ENPoint> pGPSArrayFilterSort = new ArrayList<ENPoint>();//过滤并排序后的经纬度坐标数组   17         File fGPS = new File("2007-10-14-GPS.log");//原始数据文件对象   18         File oGPS = new File("2015-12-25-GPS-Result.log");//过滤后的结果数据文件对象   19         //保持转换成度后的原始经纬度数据文件，保持格式为“ID#经纬值，纬度值”，其中经度和维度单位为度，并保留小数点后6位数字   20         File fInitGPSPoint = new File("2007-10-14-GPS-ENPoint.log");//保持转换后的原始经纬度坐标点的数据文件   21         File fTestInitPoint = new File("2007-10-14-GPS-InitTestPoint.log");//用于仿真的原始经纬度坐标点数据文件   22         File fTestFilterPoint = new File("2015-12-25-GPS-FilterTestPoint.log");//用于仿真的过滤后的经纬度坐标点数据文件   23         //-------------------------2、获取原始点坐标并将其写入到文件中-------------------------------------------------------//   24         pGPSArrayInit = getENPointFromFile(fGPS);//从原始数据文件中获取转换后的经纬度坐标点数据，存放到ArrayList数组中   25         writeInitPointToFile(fInitGPSPoint, pGPSArrayInit);//将转换后的原始经纬度点数据写入文件中   26         System.out.println(pGPSArrayInit.size());//输出原始经纬度点坐标的个数   27         //-------------------------3、进行轨迹压缩-----------------------------------------------------------------------//   28         double DMax = 30.0;//设定最大距离误差阈值   29         pGPSArrayFilter.add(pGPSArrayInit.get(0));//获取第一个原始经纬度点坐标并添加到过滤后的数组中   30         pGPSArrayFilter.add(pGPSArrayInit.get(pGPSArrayInit.size()-1));//获取最后一个原始经纬度点坐标并添加到过滤后的数组中   31         ENPoint[] enpInit = new ENPoint[pGPSArrayInit.size()];//使用一个点数组接收所有的点坐标，用于后面的压缩   32         Iterator<ENPoint> iInit = pGPSArrayInit.iterator();   33         int jj=0;   34         while(iInit.hasNext()){   35             enpInit[jj] = iInit.next();   36             jj++;   37         }//将ArrayList中的点坐标拷贝到点数组中   38         int start = 0;//起始下标   39         int end = pGPSArrayInit.size()-1;//结束下标   40         TrajCompressC(enpInit,pGPSArrayFilter,start,end,DMax);//DP压缩算法   41         System.out.println(pGPSArrayFilter.size());//输出压缩后的点数   42         //-------------------------4、对压缩后的经纬度点坐标数据按照ID从小到大排序---------------------------------------------//   43         ENPoint[] enpFilter = new ENPoint[pGPSArrayFilter.size()];//使用一个点数组接收过滤后的点坐标，用于后面的排序   44         Iterator<ENPoint> iF = pGPSArrayFilter.iterator();   45         int i = 0;   46         while(iF.hasNext()){   47             enpFilter[i] = iF.next();   48             i++;   49         }//将ArrayList中的点坐标拷贝到点数组中   50         Arrays.sort(enpFilter);//进行排序   51         for(int j=0;j<enpFilter.length;j++){   52             pGPSArrayFilterSort.add(enpFilter[j]);//将排序后的点坐标写到一个新的ArrayList数组中   53         }   54         //-------------------------5、生成仿真测试文件--------------------------------------------------------------------//   55         writeTestPointToFile(fTestInitPoint,pGPSArrayInit);//将原始经纬度数据点写入仿真文件中，格式为“经度，维度”   56         writeTestPointToFile(fTestFilterPoint, pGPSArrayFilterSort);//将过滤后的经纬度数据点写入仿真文件中，格式为“经度，维度”   57         //-------------------------6、求平均误差-------------------------------------------------------------------------//   58         double mDError = getMeanDistError(pGPSArrayInit,pGPSArrayFilterSort);//求平均误差   59         System.out.println(mDError);   60         //-------------------------7、求压缩率--------------------------------------------------------------------------//   61         double cRate = (double)pGPSArrayFilter.size()/pGPSArrayInit.size()*100;//求压缩率   62         System.out.println(cRate);   63         //-------------------------8、生成最终结果文件--------------------------------------------------------------------//   64         //将最终结果写入结果文件中，包括过滤后的点的ID，点的个数、平均误差和压缩率   65         writeFilterPointToFile(oGPS,pGPSArrayFilterSort,mDError,cRate);   66         //------------------------------------------------------------------------------------------------------------//   67     }   68    69     /**   70      *函数功能：从源文件中读出所以记录中的经纬度坐标，并存入到ArrayList数组中，并将其返回   71      * @param fGPS：源数据文件   72      * @return pGPSArrayInit：返回保存所有点坐标的ArrayList数组   73      * @throws Exception   74      */   75     public static ArrayList<ENPoint> getENPointFromFile(File fGPS)throws Exception{   76         ArrayList<ENPoint> pGPSArray = new ArrayList<ENPoint>();   77         if(fGPS.exists()&&fGPS.isFile()){   78             InputStreamReader read = new InputStreamReader(new FileInputStream(fGPS));   79             BufferedReader bReader = new BufferedReader(read);   80             String str;   81             String[] strGPS;   82             int i = 0;   83             while((str = bReader.readLine())!=null){   84                 strGPS = str.split(" ");   85                 ENPoint p = new ENPoint();   86                 p.id = i;   87                 i++;   88                 p.pe = (dfTodu(strGPS[3]));   89                 p.pn = (dfTodu(strGPS[5]));   90                 pGPSArray.add(p);   91             }   92             bReader.close();   93         }   94         return pGPSArray;   95     }   96    97     /**   98      * 函数功能：将过滤后的点的经纬度坐标、平均距离误差、压缩率写到结果文件中   99      * @param outGPSFile：结果文件  100      * @param pGPSPointFilter：过滤后的点  101      * @param mDerror：平均距离误差  102      * @param cRate：压缩率  103      * @throws Exception  104      */  105     public static void writeFilterPointToFile(File outGPSFile,ArrayList<ENPoint> pGPSPointFilter,  106                                               double mDerror,double cRate)throws Exception{  107         Iterator<ENPoint> iFilter = pGPSPointFilter.iterator();  108         RandomAccessFile rFilter = new RandomAccessFile(outGPSFile,"rw");  109         while(iFilter.hasNext()){  110             ENPoint p = iFilter.next();  111             String sFilter = p.getResultString()+"\n";  112             byte[] bFilter = sFilter.getBytes();  113             rFilter.write(bFilter);  114         }  115         String strmc = "#"+Integer.toString(pGPSPointFilter.size())+","+  116                 Double.toString(mDerror)+","+Double.toString(cRate)+"%"+"#"+"\n";  117         byte[] bmc = strmc.getBytes();  118         rFilter.write(bmc);  119   120         rFilter.close();  121     }  122     /**  123      * 函数功能：将转换后的原始经纬度数据点存到文件中  124      * @param outGPSFile  125      * @param pGPSPointFilter  126      * @throws Exception  127      */  128     public static void writeInitPointToFile(File outGPSFile,ArrayList<ENPoint> pGPSPointFilter)throws Exception{  129         Iterator<ENPoint> iFilter = pGPSPointFilter.iterator();  130         RandomAccessFile rFilter = new RandomAccessFile(outGPSFile,"rw");  131         while(iFilter.hasNext()){  132             ENPoint p = iFilter.next();  133             String sFilter = p.toString()+"\n";  134             byte[] bFilter = sFilter.getBytes();  135             rFilter.write(bFilter);  136         }  137         rFilter.close();  138     }  139     /**  140      * 函数功能：将数组中的经纬度点坐标数据写入测试文件中，用于可视化测试  141      * @param outGPSFile：文件对象  142      * @param pGPSPointFilter：点数组  143      * @throws Exception  144      */  145     public static void writeTestPointToFile(File outGPSFile,ArrayList<ENPoint> pGPSPointFilter)throws Exception{  146         Iterator<ENPoint> iFilter = pGPSPointFilter.iterator();  147         RandomAccessFile rFilter = new RandomAccessFile(outGPSFile,"rw");  148         while(iFilter.hasNext()){  149             ENPoint p = iFilter.next();  150             String sFilter = p.getTestString()+"\n";  151             byte[] bFilter = sFilter.getBytes();  152             rFilter.write(bFilter);  153         }  154         rFilter.close();  155     }  156   157     /**  158      * 函数功能：将原始经纬度坐标数据转换成度  159      * @param str：原始经纬度坐标  160      * @return ：返回对于的度数据  161      */  162     public static double dfTodu(String str){  163         int indexD = str.indexOf('.');  164         String strM = str.substring(0,indexD-2);  165         String strN = str.substring(indexD-2);  166         double d = Double.parseDouble(strM)+Double.parseDouble(strN)/60;  167         return d;  168     }  169     /**  170      * 函数功能：保留一个double数的小数点后六位  171      * @param d：原始double数  172      * @return 返回转换后的double数  173      */  174     public static double getPointSix(double d){  175         DecimalFormat df = new DecimalFormat("0.000000");  176         return Double.parseDouble(df.format(d));  177     }  178     /**  179      * 函数功能：使用三角形面积（使用海伦公式求得）相等方法计算点pX到点pA和pB所确定的直线的距离  180      * @param pA：起始点  181      * @param pB：结束点  182      * @param pX：第三个点  183      * @return distance：点pX到pA和pB所在直线的距离  184      */  185     public static double distToSegment(ENPoint pA,ENPoint pB,ENPoint pX){  186         double a = Math.abs(geoDist(pA, pB));  187         double b = Math.abs(geoDist(pA, pX));  188         double c = Math.abs(geoDist(pB, pX));  189         double p = (a+b+c)/2.0;  190         double s = Math.sqrt(Math.abs(p*(p-a)*(p-b)*(p-c)));  191         double d = s*2.0/a;  192         return d;  193     }  194   195     /**  196      * 函数功能：用老师给的看不懂的方法求两个经纬度点之间的距离  197      * @param pA：起始点  198      * @param pB：结束点  199      * @return distance：距离  200      */  201     public static double geoDist(ENPoint pA,ENPoint pB)  202     {  203         double radLat1 = Rad(pA.pn);  204         double radLat2 = Rad(pB.pn);  205         double delta_lon = Rad(pB.pe - pA.pe);  206         double top_1 = Math.cos(radLat2) * Math.sin(delta_lon);  207         double top_2 = Math.cos(radLat1) * Math.sin(radLat2) - Math.sin(radLat1) * Math.cos(radLat2) * Math.cos(delta_lon);  208         double top = Math.sqrt(top_1 * top_1 + top_2 * top_2);  209         double bottom = Math.sin(radLat1) * Math.sin(radLat2) + Math.cos(radLat1) * Math.cos(radLat2) * Math.cos(delta_lon);  210         double delta_sigma = Math.atan2(top, bottom);  211         double distance = delta_sigma * 6378137.0;  212         return distance;  213     }  214     /**  215      * 函数功能：角度转弧度  216      * @param d：角度  217      * @return 返回的是弧度  218      */  219     public static double Rad(double d)  220     {  221         return d * Math.PI / 180.0;  222     }  223       224     /**  225      * 函数功能：根据最大距离限制，采用DP方法递归的对原始轨迹进行采样，得到压缩后的轨迹  226      * @param enpInit：原始经纬度坐标点数组  227      * @param enpArrayFilter：保持过滤后的点坐标数组  228      * @param start：起始下标  229      * @param end：终点下标  230      * @param DMax：预先指定好的最大距离误差  231      */  232     public static void TrajCompressC(ENPoint[] enpInit,ArrayList<ENPoint> enpArrayFilter,  233                                      int start,int end,double DMax){  234         if(start < end){//递归进行的条件  235             double maxDist = 0;//最大距离  236             int cur_pt = 0;//当前下标  237             for(int i=start+1;i<end;i++){  238                 double curDist = distToSegment(enpInit[start],enpInit[end],enpInit[i]);//当前点到对应线段的距离  239                 if(curDist > maxDist){  240                     maxDist = curDist;  241                     cur_pt = i;  242                 }//求出最大距离及最大距离对应点的下标  243             }  244             //若当前最大距离大于最大距离误差  245             if(maxDist >= DMax){  246                 enpArrayFilter.add(enpInit[cur_pt]);//将当前点加入到过滤数组中  247                 //将原来的线段以当前点为中心拆成两段，分别进行递归处理  248                 TrajCompressC(enpInit,enpArrayFilter,start,cur_pt,DMax);  249                 TrajCompressC(enpInit,enpArrayFilter,cur_pt,end,DMax);  250             }  251         }  252     }  253     /**  254      * 函数功能：求平均距离误差  255      * @param pGPSArrayInit：原始数据点坐标  256      * @param pGPSArrayFilterSort：过滤后的数据点坐标  257      * @return ：返回平均距离  258      */  259     public static double getMeanDistError(  260             ArrayList<ENPoint> pGPSArrayInit,ArrayList<ENPoint> pGPSArrayFilterSort){  261         double sumDist = 0.0;  262         for(int i=1;i<pGPSArrayFilterSort.size();i++){  263             int start = pGPSArrayFilterSort.get(i-1).id;  264             int end = pGPSArrayFilterSort.get(i).id;  265             for(int j=start+1;j<end;j++){  266                 sumDist += distToSegment(  267                         pGPSArrayInit.get(start),pGPSArrayInit.get(end),pGPSArrayInit.get(j));  268             }  269         }  270         double meanDist = sumDist/(pGPSArrayInit.size());  271         return meanDist;  272     }  273 }

第四部分 程序结果

4.1 程序输出结果

压缩后的结果：

（1）总点数：140个点；（2）平均距离误差：7.943786；（3）压缩率：4.4444%

4.2 仿真结果

经过轨迹压缩，我们将原始经纬度坐标点转换为压缩过滤后的经纬度坐标点，我们将这两种点坐标数据分别写入两个文件中，然后根据这两个文件使用R语言进行图形绘制，分别画出压缩前和压缩后的轨迹 ，进行对比，根据对比结果就可以看出我们轨迹压缩算法是否有效，最终对比结果如图4.1所示：

第五部分 总结

本程序编写过程中，遇到了各种各样的问题，也学到了很多编程经验，下面就遇到的问题及解决方案做一个总结，最后对程序存在的不足提出改进建议。

5.1 遇到的问题及解决方案

问题1：经纬度坐标顺序问题

解决：距离公式中的参数是纬度在前经度在后，需要调整下经纬度坐标点的顺序。

问题2：距离不能为负值

解决：保证求出的距离不能为负值，加绝对值函数即可。

问题3：DP算法实现细节

解决：开始使用ArrayList数组解决下标问题，递归求解时出现巨大误差，后来改用普通数组下标进行递归，结果好多了。

5.2 存在的不足与展望

（1）进行轨迹压缩时，DP算法是最简单的一种算法，并不是最优的，可选用一些效果好的算法再次进行轨迹压缩；

（2）本次实验数据的记录为3150条，数据量不算大，如果有10亿条数据，该怎么办呢？我们可以从硬件、分布式、数据预处理、数据切分、性能好的数据仓库等方面考虑。

注：原始数据文件见360云盘链接：https://yunpan.cn/cu8BiNsvH66SD  访问密码 7be0