前前言最近在上海出差，有对数据挖掘和机器学习的实践感兴趣，想要面基的小伙伴可以联系我。

联系方式在 Resume 里面，也可以看一下我的介绍，到时候咱们可以面基面的更有针对性。

困死了，实在写不下去了，前两节定的内容公式太多，打的我脑壳子痛，急需一份睡眠来缓解一下……

前言

GPS数据是很重要的用户数据，但同时其处理也是特别困难的，在信号不好的时候容易产生噪声，很多时候用户的行为本身也自带噪声(相对于描述用户的模型来说)。GPS数据是分布在球面上的一组数据(如果忽略海拔和地球的偏心率)，而从我实际的工作中来看，能直接使用的代码也是特别的少。

GPS数据处理的好，是能挖出金矿的，本文力求从工程实践的角度出发来写GPS数据的处理，在公式的推导上也力求眼睛。

我以Java为语言(因为我司业务系统用Java，我只能用JVM系的语言，否则妥妥C+Python)，力求清楚的描述清楚GPS数据处理常用的算法和原理。

目前代码零零碎碎的都写完了，但是还没有好好的整理，我会在一个月内整理成一个完整的库，手册的话看情况写~

目录GPS数据简述GPS数据的采集

GPS数据的一般记录格式

球面上的测地线测地线的长度求取

点到测地线的最短测地线

困死爹了，实在写不下去了，以下内容作为下卷出现吧！GPS序列操作GPS序列的清洗

GPS轨迹长度求和

GPS序列的压缩

Geo Hash 计算概述简单的GeoHash模型

高精度GeoHash模型

GPS数据简述

GPS数据的采集

GPS数据主要由GPS模块采集，一般的GPS模块通过串口通信来报告自己的方位，一般都是遵循NMEA-0183协议的，这个协议可以百度到。 我们一般通过软件Filter(我还没见过哪家MCU支持这个协议的，是我孤陋寡闻了)配合模块自带的Filter功能来采集有效的GPS信息(也包括UTC时间的信息)，也可以通过串口设置采集频率和设置过滤器，各家有所不同，需要仔细阅读文档。 下图是一个随便走淘宝上找的GPS模块。

关于GPS的信号采集代码，我提供一份我写的STM32F407上的版本： gps.c

这里就不赘述了，这个不是重点。

GPS数据的一般记录格式

在这里我们不关心方向，速度，海拔，我们手动忽略这些东西，我们关心的其实是经度和纬度，一般的，我们使用负数表示南纬和西经。 经纬度的定义如下图所示

0度经线在格林尼治(领先就有定义的权利啊……)

球面上的测地线

在曲面上连接两点的最短路径我们称之为测地线，因为我们生活在一个球上，所以这个测地线还是比较好求的。

测地线的长度求取

本节的代码主要在类：GeometryPoint.java 里面可以找到。 仅仅通过经纬度来计算测地线长度是可以的，但是推导会麻烦一点，我们首先将球面坐标系(经纬度)映射到欧氏坐标：

然后根据求出的

求内积，最后得到这样的公式：

我们稍微合并一下同类项，做出来的代码是这样的：

/**** @param point1* @param point2* @return 3D-vector inner product of two normlized vectors*/

public double getInnerProduct(GeometryPoint point1, GeometryPoint point2){

return Math.sin(point1.latitude *DEG2RAD)*Math.sin(point2.latitude *DEG2RAD)+

Math.cos(point1.latitude *DEG2RAD)*Math.cos(point2.latitude *DEG2RAD)*

Math.cos(point1.longitude *DEG2RAD-point2.longitude *DEG2RAD);

}

求出了内积，求弧度并且计算出测地距离就是顺理成章的事情：

private double geodesicDistance(GeometryPoint point1, GeometryPoint point2){

double alpha = getInnerProduct(point1,point2);

if (alpha>=1.0D){

return 0.0D;

}else if (alpha<=-1.0D){

return Math.PI*EARTH_RADIUS;

}else{

return Math.acos(alpha)*EARTH_RADIUS;

}

}

点到测地线的最短距离

这部分没整理好，暂时不开源，会在一个月内整理到工程中 在此之前，我们先考虑一下我们如何描述一条测地线？ 球面上任意不同两点可以指定唯一一条测地线，其次，测地线上取3个不同点，所组成的平面一定过圆心(想一下球面上测地线的定义，这里就不赘述了，如果这里描述有误请指出)。

我们怎么描述呢？如果只用两点描述，求点到测地线的关系就难求，我们最好求测地线所在平面的法向量。 展开一下法向量其实就是：

关于求外积的代码，我是使用行列式生成的结果做的：

public static double[] outerProduct3(double[] v1, double[] v2) {

assert (v1.length == 3);

assert (v2.length == 3);

double[] returnValue = {

v1[1] * v2[2] - v1[2] * v2[1],

-v1[0] * v2[2] + v1[2] * v2[0],

+v1[0] * v2[1] - v1[1] * v2[0]

};

return returnValue;

}

这么写比较快，关于推导，我肯定不会用纸啦~Matlab大法好：

syms i j k x1 x2 x3 y1 y2 y3

z = det([i j k;x1 x2 x3;y1 y2 y3])

% 输出

z = i*x2*y3 - i*x3*y2 - j*x1*y3 + j*x3*y1 + k*x1*y2 - k*x2*y1

随后是求取法向量和球面上一点(目标向量v3)的夹角我们就可以间接得到平面和目标向量的夹角了~(用90°减一下就好)

(就这张图是自己画的，请叫我灵魂手绘师)

求完以后我们发现一个问题，有可能我们绘制的点到某个圆上的测地线并不落在v1和v2之间的测地线上，这个时候，如果我们尝试用v1,v2和n作为一组新的基去表征v3，如果描述v1和v2的两个量里有一个是负数，就说明落在外面了。这个是线性代数的基础，在此不赘述了。 那么如何用v1,v2和n作为一组新的基去表征v3呢？ 如果我们有

我们也能列出下面这个鬼：

这个时候不妨设矩阵：

那么：

求逆矩阵的函数：

/*** @param x 3×3 Matrix* @return inv(x)->3×3 Matrix*/

public static double[][] invO3(double[][] x) {

assert (x.length == 3);

double[][] returnValue = {

{

x[1][1] * x[2][2] - x[1][2] * x[2][1],

x[0][2] * x[2][1] - x[0][1] * x[2][2],

x[0][1] * x[1][2] - x[0][2] * x[1][1]

},

{

x[1][2] * x[2][0] - x[1][0] * x[2][2],

x[0][0] * x[2][2] - x[0][2] * x[2][0],

x[0][2] * x[1][0] - x[0][0] * x[1][2]

},

{

x[1][0] * x[2][1] - x[1][1] * x[2][0],

x[0][1] * x[2][0] - x[0][0] * x[2][1],

x[0][0] * x[1][1] - x[0][1] * x[1][0]

},

};

double factor = x[0][0] * x[1][1] * x[2][2]

- x[0][0] * x[1][2] * x[2][1]

- x[0][1] * x[1][0] * x[2][2]

+ x[0][1] * x[1][2] * x[2][0]

+ x[0][2] * x[1][0] * x[2][1]

- x[0][2] * x[1][1] * x[2][0];

for (double[] array : returnValue) {

assert (array.length == 3);

array[0] = array[0] / factor;

array[1] = array[1] / factor;

array[2] = array[2] / factor;

}

return returnValue;

}这样写快滴恨(手动陕西话)

不要BLAS

用Matlab配合正则之类的很轻松就生成代码了，你以为是我手写的么(手动抠鼻)

还有优化的余地

记得处理AssertionError

最后我们看一下关于测地线的类的一些业务代码： 1、在初始化的时候就求好逆矩阵，以后计算快得很。

private GeometryPoint point1;

private GeometryPoint point2;

private double[][] iM;

private double[] n;

public GeodesicLine(GeometryPoint p1, GeometryPoint p2) {

point1 = p1;

point2 = p2;

double[] v1 = point1.get3DPos();

double[] v2 = point2.get3DPos();

n = outerProduct(v1, v2);

norm2Vector(1.0, n);

double[][] M = {v1, v2, n};

iM = invO3(M);

}

2、先判断是否在范围内，如果不在范围内用端点的距离最小值作为到测地线段的距离：

public double distance(GeometryPoint point) {

double[] v3 = point.get3DPos();

double[] params = {0, 0, 0};

for (int i = 0; i < 3; i++) {

for (int j = 0; j < 3; j++) {

params[j] += iM[i][j] * v3[j];

}

}

if (params[0] > 0 && params[1] > 0) {

//返回和n的内积的角度的余角计算出的数值 return Math.asin(

n[0] * v3[0] + n[1] * v3[1] + n[2] * v3[2]

) * GeometryPoint.EARTH_RADIUS;

} else {

//两个端点中选一个近的 return Math.min(point.distance(point1), point.distance(point2));

}

}

这样，我们就求出了点到测地线的距离。