文章目录
目录
前言
随着生活质量的日益提高，地图相关技术已经渗透到人们生活的各个角落，本文对地图匹配技术，进行简单的概述和实现，为没接触过地图开放的朋友，又想要了解相关技术的提供一个参考信息，内容有误之处，还请大神指点。
一、地图匹配概述
由于GPS技术特性原因，定位结果受客观存在的环境因素影响，我们拿到的原始轨迹数据并不是我们在轨迹结果上看到的那么整洁、规范和有序，如果直接把这些点打到地图上，那么我们的行动轨迹，可能看起来就是这样：
所以信号也有颗不安分的心，这份轨迹数据让人看起来很灰心，因为勺子里面装了一些不该有的东西。那么我们需要通过一系列的技术手段，来保留我们想要的、看起来合理的轨迹，去掉看起来不合理的数据。
这一系列的技术手段，包括降噪、分段、抽稀、路网匹配、轨迹补偿等，如下图
降噪：降噪技术在数据处理方面应用得非常广泛，例如BI、机器学习等领域，数据几乎都需要进行噪点处理，降噪后的数据通常会存在多多少少的失真，这个缺陷可以通过后期的路网匹配来进行修正；
抽稀：浮点运算是个比较消耗CPU性能的计算，我们通常需要对原始数据进行抽稀处理，将原始数据精简为足以还原轨迹的数量即可；
路网匹配：由于信号偏移、降噪处理等因素，导致有些轨迹点是不在道路网上的，当我们把这些点打到地图上以后，会发现我们可能行走在山顶、湖泊、大海，显然这并不合理，所以我们需要结合路网数据，将轨迹拟合到道路网。
二、轨迹预处理
1.降噪
剔除数据中的冗余点和噪音点，通常使用，噪音滤波技术，主要算法包括中值滤波、均值滤波、卡尔曼滤波等；
轨迹压缩技术，通常使用道格拉斯压缩算法；
轨迹分段技术，通常使用网格划分算法，以及角度、密度、速度相似度匹配算法
后续小节将对各个处理节点算法进行一一详解。
1.1 中值滤波
当时由于项目时间原因，中值滤波算法是本次主要降噪算法，该算法在当前版本中，实现成本、学习成本、处理效果是最优选择。
算法原理
1、对坐标点Pi，获取前后n（n为大于2的偶数）个节点的集合R；
2、将R的经度、纬度，分别进行正序排序；
3、分别取经度和纬度的中间值，作为该节点的最终经度和纬度；
图例：
关键代码：
/ 起始索引，前面range个坐标不处理
int startIndex = STEP_SIZE / 2 + 1;
// 结束索引，后range个坐标不处理
int endIndex = tracks.size() - STEP_SIZE / 2;
// STEP_SIZE每次处理的坐标个数，除以二即为前后处理的范围
int range = STEP_SIZE / 2;
for (int i = startIndex; i < endIndex; i++) {
UserTrackDTO track = tracks.get(i);
List<Double> lngs = new ArrayList<>();
List<Double> lats = new ArrayList<>();
// 获取预定范围间的所有经度和纬度
for (int j = i - range; j <= i + range; j++) {
lngs.add(tracks.get(j).getLngDouble());
lats.add(tracks.get(j).getLatDouble());
}
// 排序
Collections.sort(lngs);
Collections.sort(lats);
// 取中间
track.setLng(String.valueOf(lngs.get(3)));
track.setLat(String.valueOf(lats.get(3)));
}中值滤波的缺点：对于连续多个噪点的数据处理得并不理想，增加步长会增加时间复杂度和严重失真，后续会考虑引入卡尔曼滤波算法看看效果；
1.2 极值滤波
算法原理
1、计算各个坐标点距离上一个点的距离、速度；
2、按比例移除距离、速度较大的点；
关键代码：
tracks.get(0).setDistance(0);
for (int i = 1; i < tracks.size(); i++) {
UserTrackDTO lastTrack = tracks.get(i - 1);
UserTrackDTO track = tracks.get(i);
Double distance = DistanceHepler.distance(track, lastTrack);
track.setDistance(distance.intValue());
}
Double maxSize = tracks.size() * 0.80D;
tracks = tracks.stream()
.sorted(Comparator.comparing(UserTrackDTO::getDistance))
.collect(Collectors.toList());
tracks = tracks.stream()
.limit(maxSize.intValue())
.collect(Collectors.toList());
tracks = tracks.stream()
.sorted(Comparator.comparing(UserTrackDTO::getIndex))
.collect(Collectors.toList());极值滤波缺点：只能按预设定的阈值进行截取，未考虑当前路况等综合情况进行阈值的设定，会导致精确度不足；
1.3 压缩
算法原理
（1）在曲线首尾两点A，B之间连接一条直线AB，该直线为曲线的弦；
（2）得到曲线上离该直线段距离最大的点C，计算其与AB的距离d；
（3）比较该距离与预先给定的阈值threshold的大小，如果小于threshold，则该直线段作为曲线的近似，该段曲线处理完毕。
（4）如果距离大于阈值，则用C将曲线分为两段AC和BC，并分别对两段曲线进行1~3的处理。
（5）当所有曲线都处理完毕时，依次连接各个分割点形成的折线，即可以作为曲线的近似。
图例：
关键代码：
/**
*
* @param coordinate
* @param dMax
* @return
*/
public static List<UserTrackDTO> douglasPeucker(List<UserTrackDTO> coordinate, int dMax) {
// 抽稀点数量需要大于2
if (coordinate == null || coordinate.size() <= 2) {
return coordinate;
}
coordinate = checkCircle(coordinate);
// 起点和终点重合时，分段处理
int endIndex = coordinate.size() - 1;
List<UserTrackDTO> result = new ArrayList<>();
compressLine(coordinate, result, 0, endIndex, dMax);
result.add(coordinate.get(0));
result.add(coordinate.get(endIndex));
return result.stream()
.sorted(Comparator.comparing(UserTrackDTO::getIndex))
.collect(Collectors.toList());
}
/**
* 递归方式压缩轨迹
*
* @param coordinate
* @param result
* @param start
* @param end
* @param dMax
* @return
*/
private static void compressLine(List<UserTrackDTO> coordinate, List<UserTrackDTO> result, int start, int end, int dMax) {
if (start < end) {
double maxDist = 0;
int currentIndex = 0;
UserTrackDTO startPoint = coordinate.get(start);
UserTrackDTO endPoint = coordinate.get(end);
for (int i = start + 1; i < end; i++) {
double currentDist = distToSegment(startPoint, endPoint, coordinate.get(i));
if (currentDist > maxDist) {
maxDist = currentDist;
currentIndex = i;
}
}
if (maxDist >= dMax) {
// 将当前点加入到过滤数组中
result.add(coordinate.get(currentIndex));
// 将原来的线段以当前点为中心拆成两段，分别进行递归处理
compressLine(coordinate, result, start, currentIndex, dMax);
compressLine(coordinate, result, currentIndex, end, dMax);
}
}
}
这个算法处理环形轨迹的时候，会将整条轨迹全部精简掉，所以这里用到分段技术。
1.3 分段
之前一个朋友建议折半分段，但是折半分段只解决了单环形轨迹，如果来个奥运五环的轨迹，又吃不消了，所以考虑网格分段和基于方向角的分段，其实网格分段连环形轨迹的问题都解决不了，所以终极方案还是基于方向角的分段。当然，关于分段还有其他的方法，这里主要用到这两种，而且网格划分不仅仅是用在这里，其他地方也用得到，比如为提取停留事件的聚类做数据准备。
网格分段
算法原理
1、对一组轨迹进行边界检索，获取最小经纬度坐标和最大经纬度坐标；
2、计算需要划分的网格数量和增长因子；
3、创建网格；
4、分配坐标点到网格；
5、通过网格索引查询网格；
原理图
类图
网格数量=Math.sqr(轨迹数量)
边数量= Math.ceil(Math.sqr(网格数量))
MeshIndex索引检索算法：二分法(BinarySearch)
关键代码
/**
* 二分查找
*
* @param startIndex
* 开始索引
* @param endIndex
* 结束索引
* @param value
* 经纬度值
* @param arr
* 值数组
* @return 返回所在范围的索引
*/
private Integer binarySearch(int startIndex, int endIndex, Long value) {
if (startIndex > endIndex) {
return -1;
}
count++;
if (startIndex == 4 && endIndex == 4) {
System.out.println("");
}
int midIndex = startIndex + (endIndex - startIndex) / 2;
midIndex = midIndex < arr.length ? midIndex : arr.length - 1;
if (isInBoundary(midIndex, value)) {
return midIndex;
} else if (isInLeft(midIndex, value)) {
return binarySearch(startIndex, midIndex - 1, value);
} else {
return binarySearch(midIndex + 1, endIndex, value);
}
}
基于方向角的分段
算法原理
有坐标集合 P∈{1,2,3…i}
遍历集合P
初始化方向角度Angle_P1_P2，进入下一次循环
计算Pi-1到Pi的方向角Angle_Pi-1_Pi
计算Angle_Pi-1_Pi与Angle_P1_P2的方向角度差
If（小于15°）
进入下一次循环
else
将初始方向角度小于15的集合P1到Pi-1记录为一个分段
初始化方向角度Angle_Pi_Pi+1，进入下一次循环
关键代码
Map<Integer, List<T>> segMap = new HashMap<>();
int startIndex = 0;
double startAngleSegment = -1;
for (int i = 5; i < points.size(); i++) {
int tmpStartIndex = -1;
if (i > 4 && i - startIndex > 4) {
tmpStartIndex = i - 5;
} else {
tmpStartIndex = startIndex;
}
T start = points.get(tmpStartIndex);
T cur = points.get(i);
double angleSegment = AngleUtil.getAngle(start, cur);
// 初始化方向角
if (startAngleSegment < 0) {
startAngleSegment = angleSegment;
continue;
}
// 计算当前方向角与初始方向角的角度差
double δAngle = Math.abs(startAngleSegment - angleSegment);
if (δAngle > 15) {
segMap.put(i, loadSegments(startIndex, i, points));
startIndex = i;
startAngleSegment = -1;
}
if (i == points.size() - 1 && startIndex < points.size() - 1) {
segMap.put(i, loadSegments(startIndex, i, points));
}
}
正北方向角计算公式代码
/**
* 获取AB连线与正北方向的角度
*
* @param A
* A点的经纬度
* @param B
* B点的经纬度
* @return AB连线与正北方向的角度（0~360）
*/
public static <T extends AbstractPoint> double getAngle(T point1, T point2) {
MyLatLng A = new MyLatLng(point1.getLongitude(), point1.getLatitude());
MyLatLng B = new MyLatLng(point2.getLongitude(), point2.getLatitude());
double dx = (B.m_RadLo - A.m_RadLo) * A.Ed;
double dy = (B.m_RadLa - A.m_RadLa) * A.Ec;
double angle = 0.0;
angle = Math.atan(Math.abs(dx / dy)) * 180. / Math.PI;
double dLo = B.m_Longitude - A.m_Longitude;
double dLa = B.m_Latitude - A.m_Latitude;
if (dLo > 0 && dLa <= 0) {
angle = (90. - angle) + 90;
} else if (dLo <= 0 && dLa < 0) {
angle = angle + 180.;
} else if (dLo < 0 && dLa >= 0) {
angle = (90. - angle) + 270;
}
return angle;
}
1.4 插值算法
插值算法这里用的是flanagan包（Dr Michael Thomas Flanagan）中的三次样条插值算法，代码如下：
public static <T extends AbstractPoint> Map<Long, Double[]> interpolate(List<T> list) {
if (list == null || list.size() < 2) {
return null;
}
double[] timespans = new double[list.size()];
double[] lngs = new double[list.size()];
double[] lats = new double[list.size()];
for (int i = 0; i < list.size(); i++) {
T t = list.get(i);
timespans[i] = t.getTimeLong().doubleValue();
lngs[i] = t.getLongitude();
lats[i] = t.getLatitude();
}
CubicSplineFast csLng = new CubicSplineFast(timespans, lngs);
CubicSplineFast csLat = new CubicSplineFast(timespans, lats);
// 计算需要插值的时间
List<Double> needInterpolationTimespans = searchExceptionDuration(list);
Map<Long, Double[]> map = new HashMap<>();
for (int i = 0; i < needInterpolationTimespans.size(); i++) {
Double timespan = needInterpolationTimespans.get(i);
double lng = csLng.interpolate(timespan);
double lat = csLat.interpolate(timespan);
map.put(timespan.longValue(), new Double[] { lng, lat });
}
return map;
}
/**
* 查找间隔大于30秒，小于15分钟的节点，以30秒为间隔进行插值
*
* @param list
* @return 需要插值的时刻点
*/
private static <T extends AbstractPoint> List<Double> searchExceptionDuration(List<T> list) {
List<Double> needInterpolation = new ArrayList<>();
for (int i = 1; i < list.size(); i++) {
T last = list.get(i - 1);
T current = list.get(i);
Long duration = current.getTimeLong() - last.getTimeLong();
if (duration.doubleValue() > MAX_DURATION_THRESHOLD && duration.doubleValue() < EXCEPTION_THRESHOLD) {
Double ceil = Math.ceil(duration.doubleValue() / MAX_DURATION_THRESHOLD);
// 计算增长因子
Double growthFactor = duration.doubleValue() / ceil;
for (int j = 1; j < ceil; j++) {
double e = last.getTimeLong().doubleValue() + growthFactor * j;
System.out.println(e);
needInterpolation.add(e);
}
}
}
return needInterpolation;
}
说明：
1、将时间作为纵轴，经度和纬度分别作为横轴进行系数计算，得到曲线方程；
2、查询并计算缺失坐标的时刻；
3、根据曲线方程计算缺失时刻的坐标；
三、基于权重的地图匹配算法
匹配算法过程：
降噪->分段->压缩->附近道路检索->匹配->组合分段
降噪算法：中值滤波、极值滤波
分段算法：网格分段、基于方向角的分段
压缩算法：道格拉斯
匹配算法：方向、距离、形状相似度结果结合权重的综合计算，理论依据请参考《基于权重的地图匹配算法》，苏海滨，陈永利，刘强
3.1、计算方向相似度（ω1）
算法原理
通过计算候选路段和轨迹的夹角余弦得到相似度结果，该结果越接近1，相似度越高。
步骤如下：
1、计算候选路段方向与正北方向的夹角θ1；
2、计算轨迹方向与正北方向夹角为θ2；
3、计算θ1和θ2角度差Δθ；
4、计算cosΔθ得到夹角余弦结果ω1；
关键代码
/**
* 计算两条线段夹角余弦
* @param start1 线段1的开始坐标
* @param end1 线段1的结束坐标
* @param start2 线段2的开始坐标
* @param end2 线段2的结束坐标
* @return
*/
public static <T extends AbstractPoint> double cosθ(T start1, T end1, T start2, T end2) {
// 把两条线段平移到原点（线段开始、结束均减去开始点的坐标）
double x1 = end1.getLatitude() - start1.getLatitude();
double y1 = end1.getLongitude() - start1.getLongitude();
double x2 = end2.getLatitude() - start2.getLatitude();
double y2 = end2.getLongitude() - start2.getLongitude();
double xi = Math.sqrt(Math.pow(x1, 2) + Math.pow(y1, 2));
double yi = Math.sqrt(Math.pow(x2, 2) + Math.pow(y2, 2));
if (xi == 0 || yi == 0) {
return 0;
}
double cos = (x1 * x2 + y1 * y2) / (xi * yi);
return cos;
}3.2、计算距离相似度（ω2）
算法原理
计算轨迹点到候选路段的投影总距离totalDistance， totalDistance最小的路段则为最近的路段，因为计算规则是距离越小权重越大，所以ω2=1/(1+totalDistance)；
转存失败重新上传取消投影距离计算公式（参考《段宗涛，等:一种改进的轨迹地图匹配算法》）
关键代码：
/**
* 计算投影dest点在start和end两点之间的投影 《段宗涛，等:一种改进的轨迹地图匹配算法》
*
* @param start
* @param end
* @param dest
* @return
*/
public static <T extends AbstractPoint> Double[] calcProjectionPosition(T start, T end, T dest) {
double k = calcK(start, end);
double doubleK = k * k;
Double[] destPosition = new Double[2];
destPosition[0] = (doubleK * dest.getLongitude() + end.getLongitude() - k * end.getLatitude() + k * dest.getLatitude()) / (doubleK + 1);
destPosition[1] = (doubleK * end.getLatitude() + dest.getLatitude() - k * end.getLongitude() + k * dest.getLongitude()) / (doubleK + 1);
return destPosition;
}3.3、计算形状相似度（ω3）
算法原理
- 计算轨迹到候选路段的投影坐标，再计算出这投影坐标的距离da；
- 计算轨迹的两点间距离dp；
- 计算距离差值Δd=dp-da，Δd越小，形状越相似
- 因为计算规则是距离越小权重越大，所以ω3=1/(1+∑Δd)；
关键代码：
public static <T extends AbstractPoint> double caculatePerpen(List<T> ts, List<RoadPoint> road) {
double total = 0.0D;
for (int i = 1; i < ts.size(); i++) {
T start = ts.get(0);
T end = ts.get(ts.size() - 1);
double distance = DistanceHepler.distance(start, end);
double perpenDistance = caculatePerpen(start, end, road.get(0), road.get(road.size() - 1));
double Δdistance = 1 / (1 + (distance - perpenDistance));
total += Δdistance;
}
return total;
}
public static <T extends AbstractPoint> double caculatePerpen(T si, T ei, T sj, T ej) {
double x1 = sj.getLatitude() - si.getLatitude();
double y1 = sj.getLongitude() - si.getLongitude();
double x2 = ei.getLatitude() - si.getLatitude();
double y2 = ei.getLongitude() - si.getLongitude();
double x3 = ej.getLatitude() - si.getLatitude();
double y3 = ej.getLongitude() - si.getLongitude();
double x = Math.pow(x2, 2) + Math.pow(y2, 2);
double u1 = (x1 * x2 + y1 * y2) / x;
double u2 = (x3 * x2 + y3 * y2) / x;
double psx = si.getLatitude() + u1 * x2;
double psy = si.getLongitude() + u1 * y2;
double pex = si.getLatitude() + u2 * x2;
double pey = si.getLongitude() + u2 * y2;
double Lper1 = Math.sqrt(Math.pow(psx - sj.getLatitude(), 2) + Math.pow(psy - sj.getLongitude(), 2));
double Lper2 = Math.sqrt(Math.pow(pex - ej.getLatitude(), 2) + Math.pow(pey - ej.getLongitude(), 2));
double d_perpen;
if (Lper1 == 0 && Lper2 == 0) {
d_perpen = 0;
} else {
d_perpen = (Math.pow(Lper1, 2) + Math.pow(Lper2, 2)) / (Lper1 + Lper2);
}
return d_perpen;
}3.4、计算权重
权重总值为1，方向、距离、形状三个相似度结果的权重分配如下：
方向：0.5
距离：0.3
形状：0.2
最终权重计算公式如下：
ω = 0.5*ω1 + 0.3*ω2 + 0.2*ω3
计算完成后，对所有候选路段的权重进行倒排，取权重结果最大的候选路段；
关键代码：
Double[] WEIGHTING_COEFFICIENT = { 0.5D, 0.3D, 0.2D };
WeightResult weightResult = new WeightResult(i);
List<RoadPoint> road = entity.getRoads().get(i);
// 1、计算方向权重
double cosθ = Calculator.cosθ(tracks, road);
weightResult.setTrackRoadθ(cosθ);
double ω1 = WEIGHTING_COEFFICIENT[0] * cosθ;
weightResult.setΩ1(ω1);
// 2、计算距离权重
double totalProjectionDistance = Calculator.calcLineDistance(tracks, road);
weightResult.setTotalProjectionDistance(totalProjectionDistance);
double ω2 = WEIGHTING_COEFFICIENT[1] / (1 + totalProjectionDistance);
weightResult.setΩ2(ω2);
// 3、计算形状权重
double shapeSimilarity = Calculator.caculatePerpen(tracks, road);
weightResult.setShapeSimilarity(shapeSimilarity);
double ω3 = WEIGHTING_COEFFICIENT[2] * shapeSimilarity;
weightResult.setΩ3(ω3);
weightResult.setWeight(ω1 + ω2 + ω3);
1678
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
