GPS 數據處理簡述（下）

01-29

前言

事實上，工作中用到的都是GPS時間序列的處理，硬體採集的也大多是帶UTC時間的時間序列，而且採樣間隔並不一定均勻，往往還有雜訊，對於這樣的時間序列，我們第一要做的是清洗。清洗結束之後是壓縮，一輛商用車，每天大約能產生10k~100k左右的GPS數據點，這樣的量對於一些計算來說數據量還不是很大，但是對於展示來說就顯得比較大。這個時候我們考慮將軌跡進行壓縮，壓縮的時候，本壓縮演算法思路基於道格拉斯-普客抽稀演算法。但是並不一樣。最後，我們通過計算GeoHash可以做到快速找距離最近的點。通過找最近的點，可以做到對路徑進行匹配。當然，我們還會使用一些小Tricks，來加速我們的計算。

代碼在：

TsingJyujing/GeoScala

GPS序列操作

GPS序列的清洗

一般來說，沒有人能在地面以高速移動，尤其是超過光速。正常的軌跡應該是醬紫的：

所以，看到這樣的GPS軌跡，我們可以肯定伊是一段雜訊：

那麼這樣的雜訊有什麼特點呢？速度他娘的快了！有的時候能超過光速！！一般我們的重卡不會超過100km/h，超過了一般交警會請你喝茶，所以我把速度限制定在150km/h是很合理的，如果是小車一般定200左右沒問題，但是德國例外，德國你得定250。定好了，我們就可以開始我們的工作了，針對每兩個點：

$v=frac{d_s}{d_t}$

其中ds的單位是km，dt的單位是小時。關於序列差分的實現，我特地請教了一下網友，Java的我當然知道，Scala里還有很優雅的實現：https://www.zhihu.com/question/65742957 所以最後清理GPS軌跡的代碼就是：

def cleanGPS[T <: Tickable](n gpsArray: Iterable[GeometryPoint[T]],n velocityLimit: Doublen ): immutable.IndexedSeq[GeometryPoint[T]] = gpsArray.iterator.sliding(2, 1).filter(data => {n (data.head.distance(data.last) / math.abs(data.head.userInfo.getTick - data.last.userInfo.getTick)) < velocityLimitn}).map(_.head).toIndexedSeqn

多一句嘴，寫代碼的時候，要麼無量綱，要麼就標清楚量綱。

GPS軌跡長度求和

有了上文的方法，那麼求和甚至可以寫成一行：

def GPSMile[T](gpsIter: Iterator[GeometryPoint[T]]): Double = gpsIter.sliding(2, 1).map {n case Seq(a, b) => a.distance(b)n}.sumn

當然，記得先清洗

GPS序列的壓縮

如果用本文的演算法，在我的破i5的機子上單核壓縮就可以達到300ms/10k點的速度，應該說速度是可以的。壓縮的問題，其實可以理解為這樣一個問題：我在軌跡中找到這樣的一個點，這個點和起始點連接以後，中途所有的點到這條直線（其實是測地線）的（測地）距離均小於某一定值（比如20m）。那麼我們首先要理解，一般的情況是什麼樣的：

對於某一個GPS序列，如果取第二個點作為終結點，那麼距離一定是0，從第三個點開始，中間存在大於0的數值，如果我們取離直線的最遠距離作為Loss，我們可以得到一個近似單調遞增的函數：

一般來說，在我們的抽稀精度和重卡的運營距離來說，你作為單調遞增序列沒毛病。那麼問題就變為一個searchInSorted的問題，一般用二分法查找：

def sparsifyGPS[T <: Tickable](n gpsArray: IndexedSeq[GeometryPoint[T]],n sparsityParam: Double,n sparsitySearchParam: Intn ) = sparsifyGPSIndexed(n gpsArray, sparsityParam, sparsitySearchParamn).map(n i => gpsArray(i)n)nndef sparsifyGPSIndexed(n gpsArray: IndexedSeq[GeometryPoint[_ <: Tickable]],n sparsityParam: Double,n sparsitySearchParam: Intn ): IndexedSeq[Int] = {n if (gpsArray.size < 10) return gpsArray.indicesn val returnList = new mutable.MutableList[Int]n returnList += 0n var nowIndex = 0n val getDistance = (startIndex: Int, endIndex: Int) => {n val line = new GeodesicLine(n gpsArray(startIndex),n gpsArray(endIndex)n )n gpsArray.slice(startIndex + 1, endIndex).map((point) => line.distance(point)).maxn }n n val loop = Breaksn loop.breakable(n while (true) {n val indexFound = SeqUtil.searchInSorted(n (i) => getDistance(nowIndex, i),n sparsityParam,n nowIndex + 2,n math.min(n nowIndex + sparsitySearchParam,n gpsArray.size - 2n )n )._1n if (indexFound >= gpsArray.size - 4) {n returnList += gpsArray.size - 1n loop.breakn } else {n returnList += indexFoundn nowIndex = indexFoundn }n }n )n returnList.toIndexedSeqn}n

只要距離不超出我們的允許範圍，我們就認可這一條折線，也不論其是不是最優。

Geo Hash 計算概述

為了快速查找地圖上的點（例如最近點），我們需要一個快速的索引演算法，這個索引演算法的原理很粗暴，就是將墨卡托坐標系劃分成M*N個網格，但是為了達到精確查找的效果我們需要對邊界進行求導，對求導點也納入計算，防止出現漏點的情況。

簡單的GeoHash模型

簡單的GeoHash很簡單，我將地圖劃分為1616的棋盤，那麼可以用0~F表示我在哪一塊，這個時候可以用一個byte表示我們的GeoHash。例如第3行第11列就是0x2A。如果進一步對網格進行劃分，劃分為256256，我們就用1個uint16描述，例如0x2AFC 進一步思考，對於任意給定的邊界我們都可以得到這樣的劃分，我們可以開發出兩層的GeoHash，即將棋盤的一個格子劃分為更加細密的棋盤，一般來說，兩層夠用了，這也就是本文代碼使用的模型，每一層的精度都是可以控制的。在初始化GeometryHashConnectionLayer的時候的兩個參數就是每一層劃分的精度。

高精度GeoHash模型

GeoHash的邊界並不是方的，是球面上一個有邊界的區域。所以最遠點、最近點不一定落在四個角上，也有可能落在邊上。這個時候我們考慮對四條邊界進行求導：設點是

$latitude=theta$

$longitude=alpha$

對上下邊界進行求導：設緯線的緯度是

$theta_1$

經度取β，這個時候緯線上任意一點的距離

$x = sin(theta_1) times sin(theta) + cos(theta_1) times cos(theta) times cos(beta-alpha)$

這個時候如果求令導數為0：

$frac{partial{d}}{partial{beta}}=frac{-sin(beta-alpha)}{sqrt{1-x^2}}=0$

很輕鬆就可以得到經度相同或者差180°的整數倍的時候導數為0，和我們的幾何直覺相同。我們也很輕鬆的可以證明，在經度相同的時候取得最小值（因為是工程問題，不考慮在極點導數處處為0的情況）。

下面我們論證經線（東西）上的導數：這個時候經度我們取α1，緯度取β。

$x = sin(beta) times sin(theta) + cos(beta) times cos(theta) times cos(alpha_1-alpha)$

$d = cos^{-1} (x)$

$frac{partial{d}}{partial{beta}}=frac{ cos(beta) times sin(theta) - sin(beta) times cos(theta) times cos(alpha_1-alpha) }{sqrt{1-x^2}}=0$

最後我們得到，導數為0的時候：

$beta = tg^{-1} (frac{tg(theta)}{cos(alpha_1-alpha)})$

只要考慮所有導數為0的點，就可以得到所有的疑似最遠最近點（當然還要配合四個角），最後我們得到：

public List<GeometryPoint> geometryHashBlockBoundary(long blockHash, long accuracy, GeometryPoint centerPoint) {n List<GeometryPoint> returnPoints = new ArrayList<GeometryPoint>();n long latCode = blockHash % POW2E31;n long lngCode = (blockHash - latCode) / POW2E31;n double lngMin = (180.0D * lngCode - 180.0D * accuracy) / ((double) accuracy);n double latMin = (180.0D * latCode - 90.00D * accuracy) / ((double) accuracy);n double lngMax = (180.0D * (lngCode + 1) - 180.0D * accuracy) / ((double) accuracy);n double latMax = (180.0D * (latCode + 1) - 90.00D * accuracy) / ((double) accuracy);n returnPoints.add(new GeometryPoint(lngMin, latMin));n returnPoints.add(new GeometryPoint(lngMin, latMax));n returnPoints.add(new GeometryPoint(lngMax, latMin));n returnPoints.add(new GeometryPoint(lngMax, latMax));nn if (lngMin < centerPoint.longitude && lngMax > centerPoint.longitude) {n returnPoints.add(new GeometryPoint(centerPoint.longitude, latMin));n returnPoints.add(new GeometryPoint(centerPoint.longitude, latMax));n }nn double[] hashBoundary = new double[2];n hashBoundary[0] = Math.cos((lngMin - centerPoint.longitude) * DEG2RAD);n hashBoundary[1] = Math.cos((lngMax - centerPoint.longitude) * DEG2RAD);n double upBoundary = Math.max(hashBoundary[0], hashBoundary[1]);n double downBoundary = Math.min(hashBoundary[0], hashBoundary[1]);n double partialDiffVariable = Math.tan(centerPoint.latitude * DEG2RAD);n double judgeMinLatitude = partialDiffVariable / Math.tan(latMin * DEG2RAD);n double judgeMaxLatitude = partialDiffVariable / Math.tan(latMax * DEG2RAD);n if (judgeMinLatitude < upBoundary && judgeMinLatitude > downBoundary) {n returnPoints.add(n new GeometryPoint(n centerPoint.longitude + Math.acos(judgeMinLatitude) * RAD2DEG,n judgeMinLatitude));n }n if (judgeMaxLatitude < upBoundary && judgeMaxLatitude > downBoundary) {n returnPoints.add(n new GeometryPoint(n centerPoint.longitude + Math.acos(judgeMaxLatitude) * RAD2DEG,n judgeMinLatitude));n }n return returnPoints;n}n

與GeoHash和而不同的HeatMap

相比之下，HeatMap的思路就簡單多了，就將網格進行計數即可：

軌跡匹配的Fréchet演算法

詳細演算法參見：https://en.wikipedia.org/wiki/Fréchet_distance

Fréchet演算法是用來評價軌跡相似度的，但是其評價的思路可以借鑒到用於匹配路徑上。一般的，路徑匹配之後又三種情況：

車輛完整的走完了所有路程
車輛走完了部分的路程
車輛沒有走該路

那我們針對軌跡上每一點，都找給定路徑上最近的一點，隨後根據每個點距離路徑的距離，我們就可以得到哪些點是在路徑上的。

按道理，軌跡有M個點，路徑有N個點，複雜度就是O(M*N)

但是用了GeoHash，複雜度就是O(Mlog(N)) 而一般軌跡比較大，如果處理得當，採用反向適配的方法可以做到O(Nlog(M))

這一塊繪圖比較困難，展開了又是一片文章，先放出代碼，以後有空詳細寫一下吧！