文章目录
目录
地图匹配技术概览
文章目录
前言
一、地图匹配概述
二、轨迹预处理
1.降噪
1.1 中值滤波
1.2 极值滤波
1.3 分段
1.4 插值算法
三、基于权重的地图匹配算法
3.1、计算方向相似度(ω1)
3.2、计算距离相似度(ω2)
3.3、计算形状相似度(ω3)
3.4、计算权重
总结
前言
随着生活质量的日益提高,地图相关技术已经渗透到人们生活的各个角落,本文对地图匹配技术,进行简单的概述和实现,为没接触过地图开放的朋友,又想要了解相关技术的提供一个参考信息,内容有误之处,还请大神指点。
一、地图匹配概述
由于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);