《概率机器人》里程计运动模型gmapping中代码解析

里程计运动模型（odometery motion model）用距离测量代替控制。实际经验表明虽然里程计虽然仍存在误差，但通常比速度运动模型更加的精确。相比于速度运动模型运动信息ut由

(x¯t−1x¯t)

为了提取相对的距离，ut被转变为三个步骤的序列：旋转，平移，另一个旋转，

上图就是里程计的测距模型，同样这些旋转和平移都是有噪声的。
首先由里程计算p(xt|ut,xt−1)的算法，算法的输入是机器人的初始位姿xt−1，从机器人里程计获得的一对位姿ut=(x¯t−1,x¯t),以及一个假定的最终姿态xt，输出的数值概率是p(xt|ut,xt−1)

算法的2-4行是从里程计的读数获取相对运动参数(δrot1 δtrans δrot2)T
第5-7行是相同的，就是计算位姿为xt−1 到x−t的相对运动参数，
第8-10行是计算误差概率，
第11行返回各自的误差概率P1 ，P2 ,P3相乘得到的组合误差概率，假定不同误差源之间是相互独立的，变量α1 α4是指定机器人运动噪声的机器人的特定参数。
基于里程计运动模型的采样算法

数学推导：根据上上面的图应该不难推导出算法中的(δrot1 δtrans δrot2)T，为了建立运动误差模型，假设旋转和平移的“真”值是测量值减去均值为0 方差为b2的独立噪声ϵb2得到也就是上述算法中的5-7行，为了计算旋转和平移的真值引入了误差参数α1 α4，所以实际位置xt，从xt−1经过初始旋转δ^rot1跟随平移距离δ^trans和另一个旋转δ^rot2得到

这里的时刻t的位姿用xt−1=(x y θ)T
那么根据给定的运动模型采样算法，，对于不同的误差参数也会有不同的概率分布：

第一个模型的采样参数是中等的，可以说是正常的，第二个和第三个扥别是比较大的平移和旋转的误差所造成的。

Gmapping中实现里程计运动模型的采样程序如下：

这个函数的输入是机器人的当前位姿pnew和上一时刻的机器人的位姿pold，

OrientedPoint MotionModel::drawFromMotion(const OrientedPoint& p, const OrientedPoint& pnew, const OrientedPoint& pold) const{    double sxy=0.3*srr;  //srr我理解为两轮子里程计的方差    OrientedPoint delta=absoluteDifference(pnew, pold); //具体如下面的介绍    OrientedPoint noisypoint(delta);//存储噪声估计    noisypoint.x+=sampleGaussian(srr*fabs(delta.x)+str*fabs(delta.theta)+sxy*fabs(delta.y));    noisypoint.y+=sampleGaussian(srr*fabs(delta.y)+str*fabs(delta.theta)+sxy*fabs(delta.x));    noisypoint.theta+=sampleGaussian(stt*fabs(delta.theta)+srt*sqrt(delta.x*delta.x+delta.y*delta.y));    noisypoint.theta=fmod(noisypoint.theta, 2*M_PI);if (noisypoint.theta>M_PI)        noisypoint.theta-=2*M_PI;    return absoluteSum(p,noisypoint);}

首先解释一下函数：

OrientedPoint delta=absoluteDifference(pnew,pold);    double 

具体的 内容如下

orientedpoint<T,A> absoluteDifference(const orientedpoint<T,A>& p1,const orientedpoint<T,A>& p2){    orientedpoint<T,A> delta=p1-p2;    delta.theta=atan2(sin(delta.theta), cos(delta.theta));    double s=sin(p2.theta), c=cos(p2.theta);    return orientedpoint<T,A>(c*delta.x+s*delta.y,                              -s*delta.x+c*delta.y, delta.theta);}

就是计算位姿的变化量，这个OrientedPoint delta的计算 结果对应的理论公式的结果就是

⎛⎝⎜(x′−x)cosθ+(y′−y)sinθ−(x′−x)sinθ+(y′−y)cosθΔθ⎞⎠⎟

计算新旧帧的绝对误差.具体就不再深入。
那么对于其中的三行代码是分别给位姿的三个两添加噪声进去，为什么要这样写呢？

noisypoint.x+=sampleGaussian(srr*fabs(delta.x)+str*fabs(delta.theta)+sxy*fabs(delta.y));    noisypoint.y+=sampleGaussian(srr*fabs(delta.y)+str*fabs(delta.theta)+sxy*fabs(delta.x));    noisypoint.theta+=sampleGaussian(stt*fabs(delta.theta)+srt*sqrt(delta.x*delta.x+delta.y*delta.y));

首先我们对sampleGaussian（b2）函数已经是有了解了，这个意思就是以均值为0 方差为b2的近似的正态分布的采样算法。所以这些括号里的都是计算方差的，为什么参数不同呢？仔细看看就能找到规律了，但是明白之前的一篇博客里是有对srr srt str stt这些高斯模型参数的设置的，分别代表delta的三个变量之间的方差而已。
接下来的程序

noisypoint.theta=fmod(noisypoint.theta, 2*M_PI);    if (noisypoint.theta>M_PI)        noisypoint.theta-=2*M_PI;

百度一下fmod()函数是对浮点型数据进行取模运算，就是计算
noisypoint.theta/2*M_PI的余数。 因为要把角度差限定在[−π π]之间。
最后是返回累计的噪声：
absoluteSum(p,noisypoint);
这个函数的实现代码是

template <class T, class A>orientedpoint<T,A> absoluteSum(const orientedpoint<T,A>& p1,const orientedpoint<T,A>& p2){    double s=sin(p1.theta), c=cos(p1.theta);    return orientedpoint<T,A>(c*p2.x-s*p2.y,                  s*p2.x+c*p2.y, p2.theta) + p1;

P1是当前机器人的位姿，前面的c*p2.x-s*p2.y, s*p2.x+c*p2.y, p2.theta是分别在P1的位姿上加上各自的噪声分量。

*这是我个人的理解，可能有一些偏差，或者错误，有错误还请指正，当然不喜勿喷