CN101303406B - 一种gps轨迹提取方法 - Google Patents

Abstract

本发明公开了属于轨迹测量和数据处理技术范围的一种GPS轨迹提取方法。该方法是利用精度为5m-15m的GPS测量仪对一段轨迹进行首次测量，划分简单轨迹，对每条简单轨迹动态多次测量，而对这段轨迹的两个端点进行静态长时间测量，得到端点坐标。两个端点连接的线段为初始解，然后逐步增加顶点来实现高精度轨迹的提取。根据点到线段、顶点和端点的距离将数据点划分到相应的最近邻域，利用非线性优化理论使得投影距离的平均值最小，找出小于设定误差的最少的顶点数及顶点的坐标，实现高精度描述这段轨迹。本发明与高精度的测量设备和系统相比较，有效地降低了成本；有效地节省了测量时间，提高了工作效率；有效地节省了数据存储空间。

Description

一种GPS轨迹提取方法

技术领域

本发明属于轨迹测量和数据处理技术范围，具体是从低成本的GPS接收机测量的多条低精度测量数据中自动提取高精度GPS轨迹的一种GPS轨迹提取方法。

背景技术

GPS即全球定位系统(Global Positioning System)。简单地说，这是一个由覆盖全球的24颗卫星组成的卫星系统。这个系统可以保证在任意时刻，地球上任意一点都可以同时观测到4颗卫星，以保证卫星可以采集到该观测点的经纬度和高度，以便实现导航、定位、授时等功能。这项技术可以用来引导飞机、船舶、车辆以及个人，安全、准确地沿着选定的路线，准时到达目的地。

GPS导航系统的基本原理是测量出已知位置的卫星到用户接收机之间的距离，然后综合多颗卫星的数据就可知道接收机的具体位置。要达到这一目的，卫星的位置可以根据星载时钟所记录的时间在卫星星历中查出。而用户到卫星的距离则通过纪录卫星信号传播到用户所经历的时间，再将其乘以光速得到(由于大气层电离层的干扰，这一距离并不是用户与卫星之间的真实距离，而是伪距(PR)：当GPS卫星正常工作时，会不断地用1和0二进制码元组成的伪随机码(简称伪码)发射导航电文。GPS系统使用的伪码一共有两种，分别是民用的C/A码和军用的P(Y)码。C/A码频率1.023MHz，重复周期一毫秒，码间距1微秒，相当于300m；P码频率10.23MHz，重复周期266.4天，码间距0.1微秒，相当于30m。而Y码是在P码的基础上形成的，保密性能更佳。导航电文包括卫星星历、工作状况、时钟改正、电离层时延修正、大气折射修正等信息。它是从卫星信号中解调制出来，以50b/s调制在载频上发射的。导航电文每个主帧中包含5个子帧每帧长6s。前三帧各10个字码；每三十秒重复一次，每小时更新一次。后两帧共15000b。导航电文中的内容主要有遥测码、转换码、第1、2、3数据块，其中最重要的则为星历数据。当用户接受到导航电文时，提取出卫星时间并将其与自己的时钟作对比便可得知卫星与用户的距离，再利用导航电文中的卫星星历数据推算出卫星发射电文时所处位置，用户在WGS-84大地坐标系中的位置速度等信息便可得知。

目前，普通低成本的GPS接收机的误差一般在15m左右，价格在2000元左右；采用高精度GPS差分系统或GPS接收机，GPS精度可以达到0.1-1m，但是成本要增加到几万到几十万。利用低成本的GPS接收机采用静态长时间测量可以提高精度，但是要测量一条轨迹，逐点静态测量显然耗费时间太长。所以，利用低成本GPS设备，高效率的提取高精度的GPS轨迹数据，必须采用新的测量方法和数据处理方法。

发明内容

本发明的目的是提供一种GPS轨迹提取方法；该方法是从低成本的GPS接收机测量的多条低精度测量数据中自动提取高精度GPS轨迹，不需要进行长时间的静态测量，提高了工作效率，满足了在低成本设备条件下获取高质量数据的市场性价比要求。

为实现上述发明目的，本发明采取如下步骤：

步骤101，设定误差精度e，用标称数据精度为p米的GPS测量仪对待测轨迹进行首次测量，将测量所得的原始数据点输入到计算机上，连接得到一条轨迹，判断待测轨迹类型是否为简单轨迹，若不是将待测轨迹分为多条简单轨迹；

步骤102，用标称数据精度为p米的GPS测量仪测量得到各条简单轨迹的两个端点的坐标和轨迹上的原始数据点坐标，将两端点连接得到初始线段；

步骤103，计算初始线段是否小于设定误差，满足则跳到判断整条轨迹提取过程结束的步骤108，不满足则进行下一步骤104；

步骤104，顶点数量增加1，计算顶点初始位置坐标；

步骤105，计算顶点折线段误差是否小于设定误差的r倍，不满足则重复上一步骤104直到满足，满足则进入下一步骤106；

步骤106，优化顶点坐标；

步骤107，计算误差平均值，看是否小于设定误差，不满足则在误差最大的区域增加一个顶点，重复上一步骤106直到满足，满足则进入下一步骤108；

步骤108，判断整条轨迹提取过程是否完毕，察看待测轨迹所拆分成的每条简单轨迹是否都已经提取完毕，如果完毕则将各条简单轨迹连接得到完整的轨迹，未完毕则回到测量简单轨迹的两个端点坐标的步骤102，进行下一条简单轨迹的测量计算。

本发明的有益效果：

1.与高精度的测量设备和系统相比较，有效地降低了成本；

2.与采用低精度设备长时间静态测量方法相比，有效地节省了测量时间，提高了工作效率；

3.与逐点存储数据的方法比较，有效地节省了数据存储空间。

附图说明

图1是本发明的流程图。

图2是测量所得原始数据图。

图3是区域划分原理图。

图4是对特定轨迹首次测量结果示意图。

图5是对特定轨迹测量所得原始数据图。

图6是顶点个数为1的计算示意图。

图7是顶点个数为2的计算示意图；

图8是顶点个数为3的计算示意图。

图9是顶点个数为4的计算示意图。

图10是顶点个数为5的计算示意图。

图11是顶点个数为6的计算示意图。

图12是第一次优化前顶点坐标示意图。

图13是第一次优化后顶点坐标示意图。

图14是第二次优化前顶点坐标示意图。

图15是第二次优化后顶点坐标示意图。

图16是第三次优化前顶点坐标示意图。

图17是第三次优化后顶点坐标示意图。

具体实施方式

下面结合附图对本发明作进一步说明：

图1为本发明的流程图，

步骤101、设定一个误差精度e，0.5m≤e≤3m。对待测轨迹进行首次测量，判断待测轨迹类型是否为简单轨迹，若不是将待测轨迹分为多条简单轨迹。

步骤102，用标称数据精度为p米(5≤p≤15)的低精度GPS测量仪对待测轨迹进行首次测量，将测量所得的原始数据点输入到计算机上，连接得到一条轨迹，将轨迹两端点相连得到一条线段，判断轨迹是否满足：无环绕，无交叉，并且轨迹上每个点向此线段及其延长线做垂线，所得垂线与轨迹无此点以外的其他相交点。如果满足则此轨迹为简单轨迹不需要拆分，如果不满足则将轨迹分为多段满足上述条件的简单轨迹分别进行测量。测量得到各条简单轨迹的两个端点的坐标，轨迹上的原始数据点坐标，将两端点连接得到初始线段，如图2所示为测量所得原始数据图。

对简单轨迹的两个端点E₁和E₂分别进行k次静态测量，得到它们的坐标，端点坐标为k次静态测量的平均值，k的取值范围为1800-3600，每秒测1次，确定端点坐标的公式为：

将两个端点直接连接得到初始线段L。

以任一端点为原点，初始线段及其延长线为x轴，建立二维直角坐标系。

步骤103、计算初始线段是否小于设定误差。

计算所有测量的原始数据点到初始线段的距离的平均值，看所得结果是否小于e(预先设定的误差要求)，如果满足则轨迹提取过程结束，否则进行步骤104。

步骤104、顶点数量增加1，计算顶点初始位置坐标。

将初始线段按照横坐标值平均分为n个平行于y轴的区域，每个区域内的原始数据点的横坐标和纵坐标分别取算术平均值得到n个顶点的初始位置坐标，将两个端点和n个顶点连接起来得到顶点折线段。

步骤105：计算顶点折线段误差是否小于设定误差的r倍。

结合图3进行说明，X_di和Y_di是数据点d_i的x，y坐标；X_Pj+1和Y_Pj+1是顶点P_j+1的x，y坐标；X_Pj和Y_Pj是顶点P_j的x，y坐标，j是顶点号，顶点数为n；i是数据点号，数据点总数为m个，优化目标e为误差平均值。

计算数据点到顶点折线段上顶点和端点的距离；数据点向顶点折线段中各线段做垂线，如果垂足在线段上，则计算数据点到此线段的距离，比较数据点到各个顶点，端点和线段的距离，取最小值，将所有到同一线段，端点或顶点距离最小的数据点划分到同一区域。如图3以2个顶点为例，图中P₁和P₂为顶点，E₁、E₂为初始线段的两个端点；所有的原始数据点划分到P₁、P₂、S₁、S_12、S₂、P_E1、P_E2七个区域，其中S₁中的数据点到线段s₁、的距离最小，P₁中的数据点到顶点p₁的距离最小，S_1，2中的数据点到线段s_1，2的距离最小，P₂中的数据点到顶点p₂的距离最小，S₂中的数据点到线段s₂的距离最小，P_E1中的数据点到端点E₁的距离最小，P_E2中的数据点到端点E₂的距离最小。

利用公式

计算数据点到顶点折线段上顶点和端点的距离；数据点向顶点折线段中各线段做垂线，利用公式 ${\mathrm{\λ}}_1=\frac{(X_{\mathrm{di}}-X_{E1})*(X_{E2}-X_{E1})+(Y_{\mathrm{di}}-Y_{E1})*(Y_{E2}-Y_{E1})}{{(X_{E2}-X_{E1})}^2+{(Y_{E2}-Y_{E1})}^2}$ 判断垂足是否在线段上，如果λ₁∈[0，1]，则说明垂足在线段上，则利用公式 ${\mathrm{DistLine}1}_{i,j}=\frac{|(X_{\mathrm{Pj}+1}-X_{\mathrm{Pj}})\×Y_{\mathrm{di}}-(Y_{\mathrm{Pj}+1}-Y_{\mathrm{Pj}})\×X_{\mathrm{di}}+X_{\mathrm{Pj}}\×Y_{\mathrm{Pj}+1}-Y_{\mathrm{Pj}}\×X_{\mathrm{Pj}+1}|}{\sqrt{{(X_{\mathrm{Pj}+1}-X_{\mathrm{Pj}})}^2+{(Y_{\mathrm{Pj}+1}-Y_{\mathrm{Pj}})}^2}}$ 计算数据点到此线段的距离，比较数据点到各个顶点，端点和线段的距离，取最小值，则原始数据点到顶点折线段的平均距离为 $D1=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\left(\mathrm{MIN}\right\{\mathrm{DistLine}{1}_{i,j}(j=\mathrm{1,2},...,n),\mathrm{DistP}{1}_{i,j}(j=\mathrm{1,2},...,n)\left\}\right)/m,$ 看所得结果D1是否满足D1≤r*e，如果不能满足则重复步骤104直到满足，如果满足则进入步骤106。

步骤106、优化顶点坐标。

在原始数据点到折线段的平均距离最大的区域内增加一个顶点P_k，P_k的初始坐标为到此线段距离最大的数据点的坐标值，然后在对此区域内三个顶点的坐标(如果区域内包含端点，不需要考虑端点)进行局部优化，即调整原有顶点P_j、P_j+1和新加顶点P_k的位置。

在此三顶点所在的区域及相邻的区域内，利用公式。

${\mathrm{DistLine}2}_{i,j}=\frac{|(X_{\mathrm{Pk}}-X_{\mathrm{Pj}})\×Y_{\mathrm{di}}-(Y_{\mathrm{Pk}}-Y_{\mathrm{Pj}})\×X_{\mathrm{di}}+X_{\mathrm{Pj}}\×Y_{\mathrm{Pk}}-Y_{\mathrm{Pj}}\×X_{\mathrm{Pk}}|}{\sqrt{{(X_{\mathrm{Pk}}-X_{\mathrm{Pj}})}^2+{(Y_{\mathrm{Pk}}-Y_{\mathrm{Pj}})}^2}}$

$+\frac{|(X_{\mathrm{Pj}+1}-X_{\mathrm{Pk}})\×Y_{\mathrm{di}}-(Y_{\mathrm{Pj}+1}-Y_{\mathrm{Pk}})\×X_{\mathrm{di}}+X_{\mathrm{Pk}}\×Y_{\mathrm{Pj}+1}-Y_{\mathrm{Pk}}\×X_{\mathrm{Pj}+1}|}{\sqrt{{(X_{\mathrm{Pj}+1}-X_{\mathrm{Pk}})}^2+{(Y_{\mathrm{Pj}+1}-Y_{\mathrm{Pk}})}^2}}$

计算得到数据点d_i到相应线段的距离。

该公式包括2个部分，在第一个部分数据点di∈S_j，k，如图3所示，S_j，k是顶点p_j和顶点p_k所连接线段对应的区域；在公式的第二个部分di∈S_k，j+1，S_k，j+1是顶点p_k和顶点p_j+1所连接线段对应的区域。

利用公式 $\begin{array}{c}\mathrm{DistP}{2}_{i,j}=\sqrt{{(X_{\mathrm{di}}-X_{\mathrm{Pj}})}^2+{(Y_{\mathrm{di}}-Y_{\mathrm{Pj}})}^2}+\sqrt{{(X_{\mathrm{di}}-X_{\mathrm{Pk}})}^2+{(Y_{\mathrm{di}}-Y_{\mathrm{Pk}})}^2}\\ +\sqrt{{(X_{\mathrm{di}}-X_{\mathrm{Pj}+1})}^2+{(Y_{\mathrm{di}}-Y_{\mathrm{Pj}+1})}^2}\end{array}$

计算得到数据点d_i到相应顶点和端点的距离；

该公式包括三个部分：第一部分中d_i∈P_j，如图3所示，P_j是顶点p_j所对应的区域；第二部分中d_i∈P_k，P_k是顶点p_j所对应的区域；第三部分中d_i∈P_j+1，P_j+1是顶点p_j+1所对应的区域。

设此三顶点所在区域及相邻的区域内的原始数据点为K个，建立优化模型然后调用MATLAB ToolBox，用最优化工具箱进行编程求解和计算，其中调用的核心求解函数为：

[x，fval，exitflag，output]＝fmincon(fun，x0，lb，ub，options)fmincon是有约束的非线性最小化求解函数，上式中的输入参数fun为要求解的模型函数，x0(矩阵)为函数的初始解，即顶点的初始位置坐标(最大误差区域三顶点的初步坐标)，lb和ub是矩阵形式的约束参数，约束条件为：

lb≤x≤ub，lb和ub中元素的坐标通过下列公式计算得出：

$X_{\mathrm{Ubpi}}=X_{\mathrm{Pi}}+\min (\mathrm{abs}\left(\frac{X_{\mathrm{Pi}}-X_{\mathrm{Pi}-1}}{2}\right),\mathrm{abs}\left(\frac{X_{\mathrm{Pi}+1}-X_{\mathrm{Pi}}}{2}\right))$

$Y_{\mathrm{Ubpi}}=Y_{\mathrm{Pi}}+\min (\mathrm{abs}\left(\frac{Y_{\mathrm{Pi}}-Y_{\mathrm{Pi}-1}}{2}\right),\mathrm{abs}\left(\frac{Y_{\mathrm{Pi}+1}-Y_{\mathrm{Pi}}}{2}\right))$

将P_j、P_j+1和P_k的坐标分别代入P_i得到ub矩阵的三个元素。

$X_{\mathrm{Lbpi}}=X_{\mathrm{Pi}}+\min (\mathrm{abs}\left(\frac{X_{\mathrm{Pi}}-X_{\mathrm{Pi}-1}}{2}\right),\mathrm{abs}\left(\frac{X_{\mathrm{Pi}+1}-X_{\mathrm{Pi}}}{2}\right))$

$Y_{\mathrm{Lbpi}}=Y_{\mathrm{Pi}}-\min (\mathrm{abs}\left(\frac{Y_{\mathrm{Pi}}-Y_{\mathrm{Pi}-1}}{2}\right),\mathrm{abs}\left(\frac{Y_{\mathrm{Pi}+1}-Y_{\mathrm{Pi}}}{2}\right))$

将P_j、P_j+1和P_k的坐标分别代入P_i得到lb矩阵的三个元素。

options是函数内部参数，输出参数x(矩阵)为优化后的求解结果(三顶点坐标)，fval为目标函数的返回值(最小误差和)，exitflag描述退出条件(程序自动设定)，output为包含的优化信息(优化方法)。

步骤107：计算误差平均值，看是否小于设定误差。

利用公式

计算数据点到顶点折线段上顶点和端点的距离；数据点向顶点折线段中各线段做垂线，利用公式 ${\mathrm{\λ}}_1=\frac{(X_{\mathrm{di}}-X_{E1})*(X_{E2}-X_{E1})+(Y_{\mathrm{di}}-Y_{E1})*(Y_{E2}-Y_{E1})}{{(X_{E2}-X_{E1})}^2+{(Y_{E2}-Y_{E1})}^2}$ 判断垂足是否在线段上，如果λ₁∈[0，1]，则说明垂足在线段上，则利用公式 ${\mathrm{DistLine}1}_{i,j}=\frac{|(X_{\mathrm{Pj}+1}-X_{\mathrm{Pj}})\×Y_{\mathrm{di}}-(Y_{\mathrm{Pj}+1}-Y_{\mathrm{Pj}})\×X_{\mathrm{di}}+X_{\mathrm{Pj}}\×Y_{\mathrm{Pj}+1}-Y_{\mathrm{Pj}}\×X_{\mathrm{Pj}+1}|}{\sqrt{{(X_{\mathrm{Pj}+1}-X_{\mathrm{Pj}})}^2+{(Y_{\mathrm{Pj}+1}-Y_{\mathrm{Pj}})}^2}}$ 计算数据点到此线段的距离，比较数据点到各个顶点，端点和线段的距离，取最小值，则原始数据点到顶点折线段的平均距离为 $D2=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\left(\mathrm{MIN}\right\{\mathrm{DistLine}{1}_{i,j}(j=\mathrm{1,2},...,n),\mathrm{DistP}{1}_{i,j}(j=\mathrm{1,2},...,n)\left\}\right)/m,$ 看所得结果是否满足D2≤e，如果不能满足则重复步骤106直到满足，如果满足则本条简单轨迹提取过程结束，进入步骤108。

步骤108：判断整条轨迹提取过程是否完毕。

察看待测轨迹所拆分成的每条简单轨迹是否都已经提取完毕，如果完毕则将各条简单轨迹连接得到完整的轨迹，如果没有完成则选取下一条未提取的简单轨迹重复步骤102进行新的简单轨迹提取过程。

下面通过一个具体实施例来说明本发明的方法：

步骤101、设定误差精度，对待测轨迹进行首次测量，判断待测轨迹类型是否为简单轨迹，若不是将待测轨迹分为多条简单轨迹。

设定一个误差精度e＝0.5m。

用标称数据精度为p＝5米的低精度GPS测量仪对待测长度约为170米的轨迹进行首次测量，花费228s测量得到约228个数据点，将测量所得的原始数据点输入到计算机上，连接得到一条轨迹，将轨迹两端点相连得到一条线段，如图4为对特定轨迹首次测量结果示意图，此轨迹满足：无环绕，无交叉，并且轨迹上每个点向此线段及其延长线做垂线，所得垂线与轨迹无此点以外的其他相交点。因此，此轨迹为简单轨迹不需要拆分。

步骤102、测量得到简单轨迹的两个端点的坐标，轨迹上的原始数据点坐标，将两端点连接得到初始线段。

用标称数据精度为p＝5米的低精度GPS测量仪对简单轨迹进行动态多次测量，将测得的原始数据也输入到计算机上，如图5为对特定轨迹测量所得原始数据图。

对简单轨迹的两个端点E₁和E₂分别进行1800次静态测量，得到它们的坐标，端点坐标为1800次静态测量的平均值，每秒测1次，共耗时1800*2＝3600s。确定端点坐标的公式为：

然后动态对轨迹测量六次，每次约228个点，加上最初的一次动态测量，共七次，共约1600s。

测量端点和原始数据点共需要1800*2+1600＝5200s，约1.44小时。

将两个端点直接连接得到初始线段L。

以起点为原点，初始线段及其延长线为x轴，建立二维直角坐标系。

步骤103、计算初始线段是否小于设定误差。

计算所有测量的原始数据点到初始线段的距离的平均值，所得结果为33.701m＞e＝0.5m，进行步骤104。

步骤104、顶点数量增加1，计算顶点初始位置坐标。

步骤105：计算顶点折线段误差是否小于设定误差的r倍。

r设为5。

计算所有测量的原始数据点到顶点折线段上各顶点、端点和线段的距离的平均值，如图6，当顶点个数为1时，误差为D1＝17.468m＞r*e＝2.5m，未达到要求，顶点坐标为：28.3091、40.0339。

计算所有测量的原始数据点到顶点折线段上各顶点、端点和线段的距离的平均值，如图7，当顶点个数为2时，误差为D1＝8.6018m＞r*e＝2.5m，未达到要求，顶点坐标为：8.9777、34.7341；59.7227、48.6460。

计算所有测量的原始数据点到顶点折线段上各顶点、端点和线段的距离的平均值，如图8，当顶点个数为3时，误差为D1＝4.7355m＞r*e＝2.5m，未达到要求，顶点坐标为：5.5176、30.3690；36.8317、62.2494；64.6958、44.5779。

计算所有测量的原始数据点到顶点折线段上各顶点、端点和线段的距离的平均值，如图9，当顶点个数为4时，误差为D1＝3.3903m＞r*e＝2.5m，未达到要求，顶点坐标为：

3.9861、27.6297；27.3959、60.9486；46.3212、61.4811；67.2496、41.4373。

计算所有测量的原始数据点到顶点折线段上各顶点、端点和线段的距离的平均值，如图10，当顶点个数为5时，误差为D1＝2.6608m＞r*e＝2.5m，未达到要求，顶点坐标为：

3.1340、25.7073；21.8216、58.7939；36.7884、62.7028；51.8192、59.8703；68.5142、39.4141。

计算所有测量的原始数据点到顶点折线段上各顶点、端点和线段的距离的平均值，如图11，当顶点个数为6时，误差为D1＝2.2296m＜r*e＝2.5m，达到要求，顶点坐标为：

2.5550、24.1530；18.0722、56.7105；30.6540、62.0805；42.8906、62.4150；55.4589、58.1747；69.4735、37.5451。。

图6----图11为寻找初始解，图11为最终确定的初始解获得初始解的时间为0.281s。

步骤106、优化顶点坐标。

步骤107：计算误差平均值，看是否小于设定误差。

在原始数据点到折线段的距离最大的区域内增加一个顶点P₁₂，如图12，P₁₂的初始坐标为到此线段距离最大的数据点的坐标值(4.7583 43.7561)，然后在对此区域内三个顶点的坐标(如果区域内包含端点，不需要考虑端点)进行局部优化，即调整原有顶点P₁、P₂和新加顶点P₁₂的位置，如图13，优化后平均误差为D2＝1.4359m＞e＝0.5m，未达到要求。

优化前顶点的坐标为：

2.5550、24.1530；4.7583、43.7561；18.0722、56.7105；30.6540、62.0805；42.8906、62.4150；55.4589、58.17469.4735、37.5451。

优化后顶点的坐标为：

1.4533、33.7704；5.8600、44.9452；17.0723、57.2799；30.6540、62.0805；42.8906、62.4150；55.4589、58.1747；69.4735、37.5451。

在原始数据点到折线段的距离最大的区域内增加一个顶点P₆₂，如图14，P₆₂的初始坐标为到此线段距离最大的数据点的坐标值(75.2614 32.0848)，然后在对此区域内三个顶点的坐标(如果区域内包含端点，不需要考虑端点)进行局部优化，即调整原有顶点P₆和新加顶点P₆₂的位置，如图15，优化后平均误差为D2＝0.75028m＞e＝0.5m，未达到要求。

优化前的顶点坐标为：

1.4533、33.7704；5.8600、44.9452；17.0723、57.279930.6540、62.0805；42.8906、62.4150；55.4589、58.1747；69.4735、37.5451；75.2614、32.0848。

优化后顶点的坐标为：

1.4533、33.7704；5.8600、44.9452；17.0723、57.2799；30.6540、62.0805；42.8906、62.4150；55.4589、58.1747；72.3675、40.2753；75.4859、30.3162。

在原始数据点到折线段的距离最大的区域内增加一个顶点P₅₆，如图16，P₅₆的初始坐标为到此线段距离最大的数据点的坐标值(62.5410 54.7695)，然后在对此区域内三个顶点的坐标(如果区域内包含端点，不需要考虑端点)进行局部优化，即调整原有顶点P₅、P₆和新加顶点P₅₆的位置，如图17，优化后平均误差为D2＝0.45113m＜e＝0.5m，达到要求。

优化前顶点的坐标为：

1.4533 33.7704；5.8600 44.9452；17.0723 57.2799；30.6540 62.0805；42.8906 62.4150；55.4589 58.1747；62.5410 54.7695；72.3675 40.2753；75.4859 30.3162。

优化后顶点的坐标为：

1.4533 33.7704；5.8600 44.9452；17.0723 57.2799；30.6540 62.0805；42.8906 62.4150；53.3059 59.8772；63.7794 53.0669；72.0904 41.5639；75.4859 30.3162。

图12——图17为局部优化的过程，图12——图13为第一次优化，图14——图15为第二次优化，图16——图17为第三次优化。图17为最终优化后的结果

优化了三次，总的优化时间为2.293s

程序最终运行结果：

最终所有顶点的坐标为：

1.4533 33.7704；5.8600 44.9452；17.0723 57.2799；30.6540 62.0805；42.8906 62.4150；53.3059 59.8772；63.7794 53.0669 72.0904 41.5639；75.4859 30.3162。

端点E1的坐标为：0、0。

端点E2的坐标为：73.8090、0。

步骤108：判断整条轨迹提取过程已完毕。

利用本方法获得的整个轨迹与利用高精度的差分GPS系统获得的高精度的GPS轨迹相比较进行检验，发现本发明最终的数据精度约为1m.同数据精度为5m低精度GPS采集设备相比，精度大大提高。

利用本方法获得的整个轨迹所需时间约为5200s+63s+1748.6s＝2小时，同利用逐点精确测量获取轨迹需要228*0.5小时＝114个小时的方法相比较，大大节省了时间。

利用本方法获得的整个轨迹所需存储空间约为228*7+2*1800＝5196个点，同利用逐点精确测量获取轨迹需要存储空间约为228*1800＝41万个点的方法相比较，大大节省了存储空间。

以上所述的实施例，只是本发明较优选的具体实施方式，本领域的技术人员在本发明技术方案范围内进行的通常变化和替换都应包含在本发明的保护范围内。

Claims (2)

1.一种GPS轨迹提取方法，其特征在于包括下列步骤：

步骤101，设定目标误差e值的取值范围为：0.5m≤e≤3m，用标称数据精度p的取值范围为：5米≤p≤15米的GPS测量仪对待测轨迹进行首次测量，将测量所得的原始数据点输入到计算机上，连接得到一条轨迹，判断待测轨迹类型是否为简单轨迹，若不是将待测轨迹分为多条简单轨迹；其中简单轨迹是将轨迹两端点相连得到一条线段，此轨迹满足：无环绕，无交叉，并且轨迹上每个点向此线段及其延长线做垂线，所得垂线与轨迹无此点以外的其他相交点；

步骤102，用标称数据精度为p米的GPS测量仪测量得到各条简单轨迹的两个端点的坐标和轨迹上的原始数据点坐标，将两端点连接得到初始线段；

步骤103，计算初始线段是否小于设定目标误差e，满足则跳到判断整条轨迹提取过程结束的步骤108，不满足则进行下一步骤104；其中计算初始线段是否小于设定目标误差e的具体方法为：计算所有测量的原始数据点到初始线段的距离的平均值，看所得结果是否小于设定目标误差e，若小于设定目标误差e，则满足误差要求；若大于设定目标误差e，则不满足误差要求；

步骤104，顶点数量增加1，计算顶点初始位置坐标；所述计算顶点初始位置坐标的具体步骤为：将初始线段按照横坐标值平均分为n个平行于y轴的区域，每个区域内的原始数据点的横坐标和纵坐标分别取算术平均值得到n个顶点的初始位置坐标，将两个端点和n个顶点连接起来得到顶点折线段；

步骤105，计算顶点折线段误差是否小于设定目标误差e的r倍，不满足则重复上一步骤104直到满足，满足则进入下一步骤106；其中r的取值范围为3-5；

所述计算顶点折线段误差是否小于设定目标误差e的r倍的具体步骤为：

(1)对原始数据点进行区域划分，数据点向顶点折线段中各线段做垂线，如果垂足在线段上，则计算数据点到此线段的距离，比较数据点到各个顶点，端点和线段的距离，取最小值，将所有到同一线段，端点或顶点距离最小的数据点划分到同一区域；

(2)利用公式

计算数据点到顶点折线段上顶点和端点的距离；数据点向顶点折线段中各线段做垂线，利用公式 ${\mathrm{\λ}}_1=\frac{(X_{\mathrm{di}}-X_{E1})*(X_{E2}-X_{E1})+(Y_{\mathrm{di}}-Y_{E1})*(Y_{E2}-Y_{E1})}{{(X_{E2}-X_{E1})}^2+{(Y_{E2}-Y_{E1})}^2}$ 判断垂足是否在线段上，如果λ₁∈[0，1]，则说明垂足在线段上，则利用公式 ${\mathrm{DistLine}1}_{i,j}=\frac{|(X_{\mathrm{Pj}+1}-X_{\mathrm{Pj}})\×Y_{\mathrm{di}}-(Y_{\mathrm{Pj}+1}-Y_{\mathrm{Pj}})\×X_{\mathrm{di}}+X_{\mathrm{Pj}}\×Y_{\mathrm{Pj}+1}-Y_{\mathrm{Pj}}\×X_{\mathrm{Pj}+1}|}{\sqrt{{(X_{\mathrm{Pj}+1}-X_{\mathrm{Pj}})}^2+{(Y_{\mathrm{Pj}+1}-Y_{\mathrm{Pj}})}^2}}$ 计算数据点到此线段的距离，比较数据点到各个顶点，端点和线段的距离，取最小值，则原始数据点到顶点折线段的平均距离为 $D1=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\left(\mathrm{MIN}\right\{\mathrm{DistLine}{1}_{i,j}(j=\mathrm{1,2},...,n),\mathrm{DistP}{1}_{i,j}(j=\mathrm{1,2},...,n)\left\}\right)/m,$ 看所得结果D1是否满足D1≤r*e；

步骤106，优化顶点坐标；

步骤107，计算误差平均值，看是否小于设定目标误差e，不满足则在误差最大的区域增加一个顶点，重复上一步骤106直到满足，满足则进入下一步骤108；所述计算误差平均值，看是否小于设定目标误差e的具体方法为：

利用公式

计算数据点到顶点折线段上顶点和端点的距离；数据点向顶点折线段中各线段做垂线，利用公式 ${\mathrm{\λ}}_1=\frac{(X_{\mathrm{di}}-X_{E1})*(X_{E2}-X_{E1})+(Y_{\mathrm{di}}-Y_{E1})*(Y_{E2}-Y_{E1})}{{(X_{E2}-X_{E1})}^2+{(Y_{E2}-Y_{E1})}^2}$ 判断垂足是否在线段上，如果λ₁∈[0，1]，则说明垂足在线段上，则利用公式 ${\mathrm{DistLine}1}_{i,j}=\frac{|(X_{\mathrm{Pj}+1}-X_{\mathrm{Pj}})\×Y_{\mathrm{di}}-(Y_{\mathrm{Pj}+1}-Y_{\mathrm{Pj}})\×X_{\mathrm{di}}+X_{\mathrm{Pj}}\×Y_{\mathrm{Pj}+1}-Y_{\mathrm{Pj}}\×X_{\mathrm{Pj}+1}|}{\sqrt{{(X_{\mathrm{Pj}+1}-X_{\mathrm{Pj}})}^2+{(Y_{\mathrm{Pj}+1}-Y_{\mathrm{Pj}})}^2}}$ 计算数据点到此线段的距离，比较数据点到各个顶点，端点和线段的距离，取最小值，则原始数据点到顶点折线段的平均距离为

$D2=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\left(\mathrm{MIN}\right\{\mathrm{DistLine}{1}_{i,j}(j=\mathrm{1,2},...,n),{\mathrm{DistP}1}_{i,j}(j=\mathrm{1,2},...,n)\left\}\right)/m,$ 看所得结果D2是否满足D2≤e；

式中X_di和Y_di是数据点d_i的x，y坐标；X_Pj+1和Y_Pj+1是顶点P_j+1的x，y坐标；X_Pj和Y_Pj是顶点P_j的x，y坐标，j是顶点号，顶点数为n；i是数据点号，数据点总数为m个，优化设定目标误差e为误差平均值；

步骤108，判断整条轨迹提取过程是否完毕，察看待测轨迹所拆分成的每条简单轨迹是否都已经提取完毕，如果完毕则将各条简单轨迹连接得到完整的轨迹，未完毕则回到测量简单轨迹的两个端点坐标的步骤102，进行下一条简单轨迹的测量计算。

2.根据权利要求1所述的GPS轨迹提取方法，其特征在于所述用标称数据精度为p米的GPS测量仪测量得到简单轨迹的两个端点的坐标和轨迹上的原始数据点坐标的具体步骤包括：

(1)用标称数据精度为p米的低精度GPS测量仪对简单轨迹进行动态a次测量，将测得的原始数据也输入到计算机上；其中a的取值范围为5-20，

(2)对简单轨迹的两个端点E₁和E₂分别进行k次静态测量，得到它们的坐标，端点坐标为k次静态测量的平均值，每秒测1次，确定端点坐标的公式为：

其中k的取值范围为1800-3600；

(3)以任一端点为原点，初始线段及其延长线为x轴，建立二维直角坐标系。

Images