与地图匹配相关的几个技术点

与地图匹配相关的几个技术点
- GPS数据采集
- - Matlab mobile
  - GPS独立硬件
  - 二者对比：
- GPS坐标转换
- GPS数据的可视化
- 交通路网
- 路网格式转换
- - 准备文件
  - 执行生成sql的命令。
- 路网搜索优化
- 地图匹配中的参数
- 分享：无人系统行业交流群

与地图匹配相关的几个技术点

自从发布若谷：Matlab算法案例— 基于隐马尔科夫模型(HMM)的地图匹配(Map-Matching) 之后，不少乎友对此很感兴趣，并联系我，需要付费索求源代码，以供学习、毕业设计使用。
奈何在提供源代码后，依旧对一些细节，不甚清楚。本文试图总结，与地图匹配相关的几个个问题，供大家借鉴。

GPS数据采集

两种方式：

Matlab mobile

Matlab有个手机版，可以获取移动轨迹（经纬度）信息，并且格式以.mat格式存储，且匹配Matlab相关函数。

GPS独立硬件

本文采用带串口的GPS硬件，用来读取GPS实时采集的信息。精度可能更高，取决于你的GPS硬件。比如，本文所示案例，其GPS硬件，精确到经纬度的秒，因此精度在15米以内。

关于GPS硬件数据采集，特别注意gPS的协议。详细解析方式，详见之前的知乎博文。

二者对比：

便捷性：Matlab mobile当然很方便，硬件成本为零。但是数据格式，主要支持Maltab。若需要转换数据格式，需要在Matlab中，进行数据读取。
数据精度：GPS独立硬件，数据精度更高，尤其是时间与位置的匹配上，由此产生的实时移动速度，也是精确的。在无遮挡的空旷地段，精度可以达到3米以内。比手机精度，要高。但是，在隧道等地段，手机端能获取数据，尽管数据不是精密测量的，很可能是通过手机自身的惯导和隧道内的通信技术，推算的结果。

GPS坐标转换

GPS坐标转换，主要是从经纬度转xyz坐标。不用多说，直接上程序。见下方：

function  [x,y,utmzone] = deg2utm(Lat,Lon)
% -------------------------------------------------------------------------
% [x,y,utmzone] = deg2utm(Lat,Lon)
%
% Description: Function to convert lat/lon vectors into UTM coordinates (WGS84).
% Some code has been extracted from UTM.m function by Gabriel Ruiz Martinez.
%
% Inputs:
%    Lat: Latitude vector.   Degrees.  +ddd.ddddd  WGS84
%    Lon: Longitude vector.  Degrees.  +ddd.ddddd  WGS84
%
% Outputs:
%    x, y , utmzone.   See example
%
% Example 1:
%    Lat=[40.3154333; 46.283900; 37.577833; 28.645650; 38.855550; 25.061783];
%    Lon=[-3.4857166; 7.8012333; -119.95525; -17.759533; -94.7990166; 121.640266];
%    [x,y,utmzone] = deg2utm(Lat,Lon);
%    fprintf('%7.0f ',x)
%       458731  407653  239027  230253  343898  362850
%    fprintf('%7.0f ',y)
%      4462881 5126290 4163083 3171843 4302285 2772478
%    utmzone =
%       30 T
%       32 T
%       11 S
%       28 R
%       15 S
%       51 R
%
% Example 2: If you have Lat/Lon coordinates in Degrees, Minutes and Seconds
%    LatDMS=[40 18 55.56; 46 17 2.04];
%    LonDMS=[-3 29  8.58;  7 48 4.44];
%    Lat=dms2deg(mat2dms(LatDMS)); %convert into degrees
%    Lon=dms2deg(mat2dms(LonDMS)); %convert into degrees
%    [x,y,utmzone] = deg2utm(Lat,Lon)
%
% Author:
%   Rafael Palacios
%   Universidad Pontificia Comillas
%   Madrid, Spain
% Version: Apr/06, Jun/06, Aug/06, Aug/06
% Aug/06: fixed a problem (found by Rodolphe Dewarrat) related to southern
%    hemisphere coordinates.
% Aug/06: corrected m-Lint warnings
%-------------------------------------------------------------------------% Argument checking
%
error(nargchk(2, 2, nargin));  %2 arguments required
n1=length(Lat);
n2=length(Lon);
if (n1~=n2)error('Lat and Lon vectors should have the same length');
end% Memory pre-allocation
%
x=zeros(n1,1);
y=zeros(n1,1);
utmzone(n1,:)='60 X';% Main Loop
%
for i=1:n1la=Lat(i);lo=Lon(i);sa = 6378137.000000 ; sb = 6356752.314245;%e = ( ( ( sa ^ 2 ) - ( sb ^ 2 ) ) ^ 0.5 ) / sa;e2 = ( ( ( sa ^ 2 ) - ( sb ^ 2 ) ) ^ 0.5 ) / sb;e2cuadrada = e2 ^ 2;c = ( sa ^ 2 ) / sb;%alpha = ( sa - sb ) / sa;             %f%ablandamiento = 1 / alpha;   % 1/flat = la * ( pi / 180 );lon = lo * ( pi / 180 );Huso = fix( ( lo / 6 ) + 31);S = ( ( Huso * 6 ) - 183 );deltaS = lon - ( S * ( pi / 180 ) );if (la<-72), Letra='C';elseif (la<-64), Letra='D';elseif (la<-56), Letra='E';elseif (la<-48), Letra='F';elseif (la<-40), Letra='G';elseif (la<-32), Letra='H';elseif (la<-24), Letra='J';elseif (la<-16), Letra='K';elseif (la<-8), Letra='L';elseif (la<0), Letra='M';elseif (la<8), Letra='N';elseif (la<16), Letra='P';elseif (la<24), Letra='Q';elseif (la<32), Letra='R';elseif (la<40), Letra='S';elseif (la<48), Letra='T';elseif (la<56), Letra='U';elseif (la<64), Letra='V';elseif (la<72), Letra='W';else Letra='X';enda = cos(lat) * sin(deltaS);epsilon = 0.5 * log( ( 1 +  a) / ( 1 - a ) );nu = atan( tan(lat) / cos(deltaS) ) - lat;v = ( c / ( ( 1 + ( e2cuadrada * ( cos(lat) ) ^ 2 ) ) ) ^ 0.5 ) * 0.9996;ta = ( e2cuadrada / 2 ) * epsilon ^ 2 * ( cos(lat) ) ^ 2;a1 = sin( 2 * lat );a2 = a1 * ( cos(lat) ) ^ 2;j2 = lat + ( a1 / 2 );j4 = ( ( 3 * j2 ) + a2 ) / 4;j6 = ( ( 5 * j4 ) + ( a2 * ( cos(lat) ) ^ 2) ) / 3;alfa = ( 3 / 4 ) * e2cuadrada;beta = ( 5 / 3 ) * alfa ^ 2;gama = ( 35 / 27 ) * alfa ^ 3;Bm = 0.9996 * c * ( lat - alfa * j2 + beta * j4 - gama * j6 );xx = epsilon * v * ( 1 + ( ta / 3 ) ) + 500000;yy = nu * v * ( 1 + ta ) + Bm;if (yy<0)yy=9999999+yy;endx(i)=xx;y(i)=yy;utmzone(i,:)=sprintf('%02d %c',Huso,Letra);
end

GPS数据的可视化

不考虑背景地图的话，GPS数据的可视化，是特别简单的。本文主要考虑带有地图遥感影像的GPS可视化。本文，主要考虑将gps数据，制作成KML文件，在google Earth上，进行可视化显示。

经过地图辅助后，就变得漂亮多了。

以上地图，是基于高德导航地图，进行的数据可视化。高德地图，还提供了动态演示效果的展示。如下图所示：

交通路网

交通路网数据，有很多下载渠道。本文选择了最通用的渠道，就是从开源地图Open Street Map上，免费下载。如图，选择地图区域，并且导出。

不过，需要指出的是，OSM（Open Street Map）的地图，格式不一定满足数据处理的要求，往往需要进行格式转换，才能进一步使用。

路网格式转换

本文，就是把Open Street Map 转成为SHAPE file格式后，进行地图匹配的。毕竟，shapefile格式的软件基础ArcGIS是最广泛的。
其中一个工具是ArcMap，进行格式转换。

具体格式转换，详见下方图示中的按钮和命令。

转化成shapefile格式之后，还要转成编程语言能便捷操作的，比如存储至关系数据库中。
有一个比较常用的工具，shp2mysql。值得推荐，从2009年第一次使用，到现在，陆陆续续地使用着。之前的博客，也有介绍。
https://zhuanlan.zhihu.com/p/102120251

准备文件

一个shp文件有多个文件组成, .dbf, .sbn, .sbx, .shp, .shx.
将shp以及shp的相关文件和DOShere文件放在同一目录下。

执行生成sql的命令。

shp2MySQL shp文件名.shp 表名数据库名 > 生成的sql文件名.sql
进入该文件目录，执行下述命令：

shp2MySQL Line.shp line hmm4mmYuan > line4mm.sql

路网搜索优化

面对浩如烟海的路网数据，为了优化路网的搜索效率，GIS的空间分析中，常常采用的是四叉树的数据结构，对路网进行分块、分级存储。本文不进行细致说明。有意者，请关注之前的HMM一文，其包含的代码部分，有详细介绍。

地图匹配中的参数

基于隐马尔科夫链的地图匹配中，由于要建设5项参数，其中概率转移时，要考虑前后两个点，到线的映射关系，这其中的转移概率，受到线段长度与相邻GPS点的平均距离。

相邻点的最大、最小、平均距离；
道路线的最大、最小、平均长度；
道路线的总数量；

分享：无人系统行业交流群

若谷：自动驾驶行业交流群及公约

自动驾驶技术是实现更高程度的智能生活的技术基础，目前已经广泛应用于无人车/船/机。
欢迎大家交流无人系统行业的技术发展（工业数据采集，感知(雷达/毫米波/红外)，人工智能算法，图像识别，精密控制，建模与仿真）、产业全供应链融合、岗位招聘与人才求职、创新创业与融资扩产、行业羽毛球联谊赛。
我是无人驾驶辅助系统ADAS企业的开发者，撰写《无人驾驶–视觉·感知·控制》