做完单帧匹配后,能得到局部的单帧重建车道线,这时以GPS位置为半径搜索高精度地图,横向定位的目的就是将单帧重建出来的车道线向高精度地图配准,之所以叫横向定位是因为车道线只能确定无人车处于那条车道,而不能确定其距离某个路标有多远。

由于GPS信号的误差大概在10米左右,一般都不会定位到当前车道,可能会错1-2个车道,如上图所示。并且由于偏航角bearing是通过与上一帧的向量确定的,并不太准,这时单帧重建结果与高精度地图还会存在一个角度上的偏差,当然这里的定位误差与偏航角误差可以通过卡尔曼滤波进行平滑与预测,这个后面会说到。总之,虽然它们存在一个偏差,但是这个偏差是符合刚性变换的,可以通过一个刚性变换RT纠正过来的。

要求RT,需要得到同名点,首先需要线线匹配,这里应用了一个相似度算法,步骤如下:

1、将单帧重建车道线编号,如s0、s1、s2、s3,同时将高精度地图车道线也进行编号,如h0、h1、h2、h3。这时可能匹配的组合就是(s0,h0)、(s0, h1)、(s0, h2)、(s0, h3) ... ... (s3, h3),共有16种组合,将其记为n0, n1, n2 ... n15,每一个组合都代表了一种匹配可能,我们就是要从中找到最符合实际情况的4对组合。

2、将n0, n1, ... n15也两两组合起来,如(n0, n0)、(n0, n1),..., (n15, n15),这样的组合共有225对,对每一对计算相似度

,记为score0-0,scpre0-1,....score15-15, 显然这里有些组合是不存在的,如(n0, n1)代表(s0, h0),(s0, h1)这样的匹配组合,s0不可能与h0,h1同时匹配,那就将(n0,n1)的相似度score0-1设为0,凡是不存在的组合相似度都设为0,对存在的组合计算相似度。

如上图所示,高精度地图上下最短距离分别为s2,s1,单帧重建上下最短距离分别为d2,d1,则定义几何相似度为:

geometry_score = exp(-std::fabs((d2/s2 + d1/ s1) / 2.0 - 1.0));

同时定义属性相似度:

double property_score = 0.0;if (l1.mark_type == r1.mark_type)property_score += 0.5;if (l2.mark_type == r2.mark_type)property_score += 0.5;

总得分:score = geometry_score + property_score

将(nk, np)与其他所有的组合进行相似度计算,并将相似度之和作为其总体的相似度scorek-p,对(n0, n0) , (n0, n1),...(n15, n15)每一对组合都如法炮制,最有可能的匹配组合就是他们之中的最大值。因此,最优匹配组合就选出来了。

3、对于次优组合的选择,要先更新组合相似度值s0-0,...s15-15,更新的原则是,如果与最优匹配组合有冲突,则相似度设置为0,否则不变,然后重新进行第二步计算,选出次优的。如果次优组合的score与最优组合的score相差不多,则代表这可能这种匹配方式是有歧义的,这发生在单帧重建结果与高精度地图不是1对1时,有歧义的话,很可能同名点会找错,遇到这种情况应该忽略掉这次匹配,靠后续的粒子滤波去大概估一个位置。

4、重复2~3过程,直到所有的组合都被选出来。

完成的代码如:

double utils::Score(const DlLane3D &l1,const DlLane3D &l2,const hdmap_lane &r1,const hdmap_lane &r2)
{if (!geometry::IsOverlap(l1.geometry, l2.geometry))return 0.0;if (!geometry::IsOverlap(r1.geometry, r2.geometry))return 0.0;double head_disL, tail_disL, head_disR, tail_disR;geometry::GetHeadTailDistance(l1.geometry, l2.geometry, head_disL, tail_disL);geometry::GetHeadTailDistance(r1.geometry, r2.geometry, head_disR, tail_disR);double ratioL = std::fabs(head_disR) > EPS ? head_disL / head_disR : head_disL / EPS;if (std::fabs(head_disL - head_disR) > 2.0)return 0.0;double ratioR = std::fabs(tail_disR) > EPS ? tail_disL / tail_disR : tail_disL / EPS;if (std::fabs(tail_disL - tail_disR) > 2.0)return 0.0;double ratio = (ratioL + ratioR) / 2.0;double geometry_score = std::exp(-std::fabs(ratio - 1.0));double property_score = 0.0;if (l1.mark_type == r1.mark_type)property_score += 0.5;if (l2.mark_type == r2.mark_type)property_score += 0.5;return geometry_score + property_score;
}bool utils::Matching(const std::vector<DlLane3D> &dlLanes, const std::vector<hdmap_lane> &hdLanes, std::vector<match_pair> &match_results)
{std::vector<match_pair> candidate_match_pairs;for (std::size_t i = 0; i < dlLanes.size(); i++){for (std::size_t j = 0; j < hdLanes.size(); j++){if (geometry::IsOverlap(dlLanes[i].geometry, hdLanes[j].geometry)){candidate_match_pairs.emplace_back(i, j);}}}if (candidate_match_pairs.empty()){LOG("candidate_match_pairs is empty, match failed!") << std::endl;return false;}unsigned int dim = candidate_match_pairs.size();Eigen::MatrixXd affineMtx = Eigen::MatrixXd::Zero(dim, dim);for (std::size_t i = 0; i < dim; i++)  {unsigned int idL1 = candidate_match_pairs[i].source;unsigned int idR1 = candidate_match_pairs[i].target;for (std::size_t j = i + 1; j < dim; j++){unsigned int idL2 = candidate_match_pairs[j].source;unsigned int idR2 = candidate_match_pairs[j].target;if ((idL1 == idL2) || (idR1 == idR2)){continue; //not satisfy the one to one match condition}DlLane3D l1 = dlLanes[idL1];DlLane3D l2 = dlLanes[idL2];hdmap_lane r1 = hdLanes[idR1];hdmap_lane r2 = hdLanes[idR2];affineMtx(i, j) = affineMtx(j, i) = Score(l1, l2, r1, r2);}}std::vector<unsigned int> match_array;auto is_exist = [&](unsigned int id) -> bool {for (size_t i = 0; i < match_array.size(); i++){if (match_array[i] == id)return true;}return false;};auto is_conflict = [&](unsigned int id) -> bool {for (size_t i = 0; i < match_array.size(); i++){if (candidate_match_pairs[match_array[i]].source == candidate_match_pairs[id].source || candidate_match_pairs[match_array[i]].target == candidate_match_pairs[id].target){return true;}}return false;};auto is_ambiguous = [&](const std::vector<std::pair<unsigned int, double>> &sort_candidates, double simility_delta) -> bool {unsigned int sId = sort_candidates[0].first;for (std::size_t i = 1; i < sort_candidates.size(); i++){match_pair match_point = candidate_match_pairs[sort_candidates[i].first];if (match_point.source == candidate_match_pairs[sId].source || match_point.target == candidate_match_pairs[sId].target){double score0 = sort_candidates[0].second;double score1 = sort_candidates[i].second;double score_delta = score0 - score1;if (score_delta < simility_delta){return true;}break;}}return false;};unsigned int num = 0;while (num < dim){std::vector<std::pair<unsigned int, double>> sort_candidates;for (unsigned int i = 0; i < dim; i++){sort_candidates.emplace_back(i, 0.0);}for (unsigned int i = 0; i < dim; i++){if (!is_exist(i) && !is_conflict(i)){for (int j = 0; j < dim; j++){if (!is_conflict(j) || is_exist(j)){sort_candidates[i].second += affineMtx(i, j);}}}}std::sort(sort_candidates.begin(), sort_candidates.end(),[=](const std::pair<unsigned int, double> &lhs, const std::pair<unsigned int, double> &rhs) {return lhs.second > rhs.second;});if (sort_candidates[0].second < EPS){break;}if (num == 0 && sort_candidates.size() > 1 && hdLanes.size() != dlLanes.size()){if (is_ambiguous(sort_candidates, 0.7)){LOG("the matching is ambiguous, matching failed!") << std::endl;return false;}}unsigned int sId = sort_candidates[0].first;match_array.push_back(sId);num++;}for (unsigned sId : match_array){match_results.push_back(candidate_match_pairs[sId]);}return true;
}

5、接下来就很简单了,让相互匹配的线线之间做垂足,就可得到同名点集合matching_points。

得到同名点后,下一步就是计算RT,由于同名点的精度并不是很高(收到语义分割、已经pitch、yaw、roll等误差的影响),我们并不想优化RT全部的自由度,而是有选择性的优化,即不优化roll,只优化pitch、yaw、t1、t2、t3这5个自由度。传统的ICP方法没法固定某一个自由度,这里采用的是非线性优化的方法,效率略有下降,但好处是可以只对部分自由度进行优化。

X'表示X经过一次平移旋转后的结果,,,

表示绕X轴旋转角的旋转矩阵:

表示绕Y轴旋转角的旋转矩阵:

表示绕Z轴旋转角的旋转矩阵:

假设将R设为矩阵模式,如下:

则有:

写成完整形式为:

根据泰勒公式展开:

                         (式1)

而:

旋转矩阵稍复杂一些,推导过程如下:

而:

代入前式得:

同理

而:

代入前式得:

k角微分稍有不同:

而:

代入前式得:

这里有个性质:正交矩阵各元素等于其代数余子式,因此

将(式1)变形,同时将X‘’写成向量形式可得:

等式左边就是雅克比矩阵,右边是残差,可采用高斯牛顿算法,为了鲁棒性及效率,可采用阻尼法及M估计。本来这是6个自由度的问题,但是因为精度的原因,并不会优化roll,因为我认为它是不会变化的。这也是为什么采用欧拉角作为状态的原因,至于欧拉角固有的万向锁问题在这里是不存在的,因为万向锁只会在pitch=90度程序,显然这在现实世界中是不存在的:

bool utils::Translate3d_3d_LM(const std::vector<accordance_pair> &match_points, Eigen::Matrix3d &R, Eigen::Vector3d &T, double& error)
{if (match_points.empty()){LOG("match_points is empty, translate failed!") << std::endl;return false;}int N = match_points.size();std::vector<Eigen::Vector3d> q1(N), q2(N);for (int i = 0; i < N; i++){q1[i] = match_points[i].source;q2[i] = match_points[i].target;}double omga = .0, kappa = .0, t1 = .0, t2 = .0, t3 = .0;double stopThresholdLM = 0.0;double currentEro = 0.0;double currentLambda = 0.0;double ni = 2.0;//Generate init lamda{Eigen::MatrixXd H = Eigen::MatrixXd::Zero(5, 5);Eigen::MatrixXd g = Eigen::MatrixXd::Zero(5, 1);for(int i = 0; i < N; i++){Eigen::Vector3d trans = geometry::RotationAroundY(omga) * geometry::RotationAroundZ(kappa) * Eigen::Vector3d(q1[i].x(), q1[i].y(), q1[i].z()) + Eigen::Vector3d(t1, t2, t3);Eigen::Vector3d residual = trans - q2[i];currentEro += residual.transpose() * residual;double x = trans.x();double y = trans.y();double z = trans.z();Eigen::Matrix<double, 3, 5> J_aux;J_aux << (t3 - z), -std::cos(omga) * (y - t2), 1.0, 0.0, 0.0,0.0, (std::cos(omga) * (x - t1) + std::sin(omga) * (z - t3)), 0.0, 1.0, 0.0,(x - t1), -std::sin(omga) * (y - t2), 0.0, 0.0, 1.0;H += J_aux.transpose() * J_aux;g -= J_aux.transpose() * residual;}stopThresholdLM = 1e-6 * currentEro;double maxDiagonal = 0;for (int i = 0; i < 5; ++i){maxDiagonal = std::max(std::fabs(H(i, i)), maxDiagonal);}double tau = 1e-5;currentLambda = tau * maxDiagonal;}auto robustWeightCauchy = [=](double norm_res) -> double {return std::sqrt(1.0 / (1.0 + norm_res * norm_res));};int iter = 0;bool stop = false;while (!stop && iter < 30){bool oneStepSuccess = false;int false_cnt = 0;while (!oneStepSuccess){Eigen::MatrixXd H = Eigen::MatrixXd::Zero(5, 5);Eigen::MatrixXd g = Eigen::MatrixXd::Zero(5, 1);for (int i = 0; i < N; i++){Eigen::Vector3d trans = geometry::RotationAroundY(omga) * geometry::RotationAroundZ(kappa) * Eigen::Vector3d(q1[i].x(), q1[i].y(), q1[i].z()) + Eigen::Vector3d(t1, t2, t3);double x = trans.x();double y = trans.y();double z = trans.z();double err_i_norm = (trans - q2[i]).norm();double weight = robustWeightCauchy(err_i_norm);Eigen::Matrix<double, 3, 5> J_aux;J_aux << (t3 - z), -std::cos(omga) * (y - t2), 1.0, 0.0, 0.0,0.0, (std::cos(omga) * (x - t1) + std::sin(omga) * (z - t3)), 0.0, 1.0, 0.0,(x - t1), -std::sin(omga) * (y - t2), 0.0, 0.0, 1.0;H += J_aux.transpose() * J_aux * weight;g -= J_aux.transpose() * (trans - q2[i]) * weight;}for (int i = 0; i < 5; i++){H(i, i) += currentLambda;}Eigen::Matrix<double, 5, 1> delta = H.colPivHouseholderQr().solve(g);if ((delta.head(2).squaredNorm() < EPS && delta.tail(3).squaredNorm() < EPS) || false_cnt > 10){stop = true;break;}omga += delta[0];kappa += delta[1];t1 += delta[2];t2 += delta[3];t3 += delta[4];double scale = 0;scale = delta.transpose() * (currentLambda * delta + g);scale += 1e-3;double temp_err = 0.0;for (int i = 0; i < N; i++){Eigen::Vector3d trans = geometry::RotationAroundY(omga) * geometry::RotationAroundZ(kappa) * Eigen::Vector3d(q1[i].x(), q1[i].y(), q1[i].z()) + Eigen::Vector3d(t1, t2, t3);Eigen::Vector3d residual = trans - q2[i];temp_err += residual.transpose() * residual;}double rho = (currentEro - temp_err) / scale;if (rho > 0 && std::isfinite(temp_err)) // last step was good{double alpha = 1.0 - pow((2 * rho - 1), 3);alpha = std::min(alpha, 2. / 3.);double scaleFactor = (std::max)(1. / 3., alpha);currentLambda *= scaleFactor;ni = 2;currentEro = temp_err;oneStepSuccess = true;}else{currentLambda *= ni;ni *= 2;oneStepSuccess = false;}if(!oneStepSuccess){omga -= delta[0];kappa -= delta[1];t1 -= delta[2];t2 -= delta[3];t3 -= delta[4];false_cnt++;}else{false_cnt = 0;}}iter++;if (sqrt(currentEro) <= stopThresholdLM)stop = true;}R = geometry::RotationAroundY(omga) * geometry::RotationAroundZ(kappa);T = Eigen::Vector3d(t1, t2, t3);error = std::sqrt(currentEro / N);return true;
}

至此,RT就可以求出,将其应用与gps坐标和航向,完成横向定位:

高精地图应用(四)横向定位相关推荐

  1. 基于Autoware制作高精地图(四)

    基于Autoware制作高精地图(四) 来了来了!它来了!肯定有小伙伴遇到过,当用Autoware导入自己制作的高精地图(也就是.csv文件)的时候会出现带有方向的车道线lane不显示的情况,或者显示 ...

  2. LaneLoc:基于高精地图的车道线定位

    文章:LaneLoc: Lane Marking based Localization using Highly Accurate Maps 作者:Markus Schreiber, Carsten ...

  3. 智能网联汽车高精地图白皮书(2020)

    1. 前言 高精地图的发展与智慧交通.智能网联汽车紧密相关,从智能网联汽车上路伊始,高精地图产业就应势而生并飞速发展.相对于以往的导航地图,高精地图是智能网联汽车交通的共性基础技术,其服务的对象并非仅 ...

  4. 高精地图数据应用分发引擎建设实践

    1. 什么是高精数据分发引擎 1.1 高精地图概述 高精地图(High Definitation Map,HD MAP),和普通导航电子地图的主要区别是精度更高.信息更丰富.精度更高主要体现在高精地图 ...

  5. 从零学习自动驾驶—百度Apollo高精地图

    作者|晓畅Auto 来源|知乎.深蓝AI 本篇文章主要介绍主流自动驾驶实现方案中举足轻重的一个部分--高精地图.之所以称其为主流方案,是因为总有一个奇葩,那就是特斯拉,偏偏不走寻常路.而除此之外,包括 ...

  6. Apollo进阶课程 ⑥ | 高精地图与自动驾驶的关系

    目录 1)高精地图与自动驾驶 2)什么是高精地图 3)高精地图与导航地图 4)高精地图---无人驾驶的核心基础模块 5)高精地图与定位模块的关系 6)高精地图与感知模块的关系 7)高精地图与规划.预测 ...

  7. Apollo自动驾驶入门课程第②讲 — 高精地图

    目录 1. 高精地图与传统地图 2. 高精地图与定位.感知规划的关系 2.1 高精地图用于定位 2.2 高精地图用于感知 2.3 高精地图用于规划 3. Apollo高精度地图与构建 3.1 Apol ...

  8. 自动驾驶系统入门(二) - 车辆定位与高精地图

    1. 车辆定位 自动驾驶中按照不同的定位实现技术,高精度定位可以分为三类:1)基于信号的定位,典型代表就是GNSS定位,即全球卫星导航系统:2)航迹推算 -IMU(惯性测量单元),其根据上一时刻的位置 ...

  9. 自动驾驶入门技术(2) —— 车辆定位与高精地图

    1. 车辆定位 自动驾驶中按照不同的定位实现技术,高精度定位可以分为三类:1)基于信号的定位,典型代表就是GNSS定位,即全球卫星导航系统:2)航迹推算 -IMU(惯性测量单元),其根据上一时刻的位置 ...

  10. Apollo星火计划学习笔记第四讲2——高精地图定位模块

    Apollo学习笔记 零.目录 一.定位的作用 二.定位用到的算法 2.1 GPS 2.2 IMU 2.3 GNSS(GPS+IMU) 2.4 先验地图定位 2.5 实时定位和建图 2.6 小结 三. ...

最新文章

  1. 美团面试题:JVM 堆内存溢出后,其他线程是否可继续工作?
  2. jython在MyEclipse控制台出现Failed to install
  3. SFB 项目经验-12-为某上市企业的Skype for Business购买Godday证书
  4. Java Math的 floor,round和ceil的总结 ,java基础知识
  5. dev 报表设计器 怎么设置每页10行_可嵌入您系统的.NET 报表控件ActiveReports:带状列表组件...
  6. mongodb上一篇下一篇_如何使用Microsoft office word—上一篇
  7. Python数据挖掘与机器学习,快速掌握聚类算法和关联分析
  8. 20190913:(leetcode习题)罗马数字转整数
  9. 工程师追查线上问题(或运维)常用的shell命令
  10. python中分支语句elif与else的区别_浅谈对python中if、elif、else的误解
  11. 医院管理数据库课程设计
  12. Windows Server 2012/2012R2 中配置 MSDTC,令其使用特定端口
  13. r7 6800u核显相当于什么显卡
  14. 从键盘输入一个整数,判断它是正数,负数,0
  15. Windows——卸载MinGW的方法
  16. AI智能视频批量剪辑软件开发-云罗企客-视频一键批量处理
  17. ffmpeg源码分析与应用示例(一)——H.264解码与QP提取
  18. Hive 基础知识(二)
  19. html网页弹幕特效,jquery仿哔哩哔哩弹幕文字动画特效
  20. 代码大全2 --- 33章 个人性格

热门文章

  1. linux系统底层,干货|七点,用计算机底层知识教你安装Linux系统!
  2. 计算机软件存储位置,微信电脑版存储位置在什么地方?查看微信电脑版存储路径的方法...
  3. [转载]Core animation简介
  4. Oracle技巧查询,很香
  5. CSS常用定位方法(绝对定位、相对定位、固定定位)
  6. android root工具排行榜,可root任何机?史上最强安卓root工具出炉
  7. 二叉搜索树的模拟及其实现(c++)
  8. Oracle使用函数达到drop table if exists
  9. Linux系统并搭建Sip server平台
  10. linux系统amd驱动怎么安装教程,ubuntu amd显卡驱动安装教程