合肥网站建设技术外包,怎样给公司做一个网站做推广,韩国虚拟空间网站,网站企业管理培训课程一、参考资料
深入浅出理解卡尔曼滤波
二、Allan方差分析方法
1. 引言
传统的误差指标往往是采用均值误差#xff08;反映整个误差序列有无宏观偏置#xff09;、标准差#xff08;反映整个误差序列的波动情况#xff09;#xff0c;以及均方根#xff08;RMS#x…一、参考资料
深入浅出理解卡尔曼滤波
二、Allan方差分析方法
1. 引言
传统的误差指标往往是采用均值误差反映整个误差序列有无宏观偏置、标准差反映整个误差序列的波动情况以及均方根RMS可以认为是宏观偏置与波动情况的综合。如下图所示它们都是反映误差序列的整体情况的指标其中含有短时间快速变化和长时间缓变的各个成份无法细分出不同时间尺度上的误差波动情况。
Allan方差是一个相对而言直观形象、朴实无华的概念它当初是在实践中针对高精度时钟稳定性分析的需求而提出的。Allan方差适用于分析惯导和时钟看重长期稳定性的传感器Allan方差曲线对中长期时间尺度上的误差特性有很强的表现力。Allan方差相比于普通的方差具有更好的描述长时间误差的优势因此Allan方差被用来当作IMU噪声标定的标准方法。
2. Allan方差的概念
Allan方差又称为阿伦方差或阿兰方差Allan方差法是20世纪60年代由美国国家标准局的David Allan提出的基于时域的分析方法。
设计Allan方差分析方法的目的是为了实现这样一种分析对惯性传感器误差的关注主要体现在不确定性和稳定性具体来说是不同时间尺度的稳定性。例如战术武器领域关心的是多少秒到多少分钟时间尺度的稳定性潜水艇领域关心的是天、轴周、月时间尺度的稳定性。那么有没有什么方法将不同时间尺度上的不确定性分门别类的表现出来Allan教授针对这个问题提出了Allan方差分析方法。
Allan方差是时频分析和惯性导航领域常用的一种误差分析方法它有效地刻画了待研究的误差时间序列在不同时间尺度上的波动水平不稳定性并可根据不同时间尺度上的Allan方差值所构成的曲线的形状特征来辨识其中包含的随机过程模型。Allan方差分析方法对中长期的随机波动具有很强的表现力它完全可以作为一个通用的时间序列分析工具来推广到其它应用领域。
在我们对惯性导航器件陀螺和加速度计进行误差分析时常采用Allan方差分析方法。这是当初从时频领域高精度时钟的频率稳定性借鉴来的一种时间序列分析方法 非常适合于对中长期波动不稳定性进行定量描述和分析因此正是我们惯性器件误差分析所需要的。
3. Allan方差计算过程
Allan方差是将误差序列在某个指定的时间尺度上的波动情况进行了精确提取其具体计算步骤如下
Step 1 分块
按窗口长度 τ ( n − 1 ) t 0 \tau\left(\left.n-1\right.\right)t_0 τ(n−1)t0 将序列 y 分成 N c N_c Nc 块(clusters)每块包n个数据点块间无重叠。
Step 2 块平均
分别计算各块内n个数据点的均值记为 y ‾ k 1 n ∑ i k k n − 1 y i \begin{aligned}\overline{y}_k\frac1n\sum_{i\operatorname{}k}^{kn-1}y_i\end{aligned} ykn1ik∑kn−1yi。
Step 3 计算Allan方差
相邻块的平均值求差。 σ 2 ( τ ) 1 2 ( N c − 1 ) ∑ k 1 N c − 1 ( y ˉ k 1 − y ˉ k ) 2 \begin{aligned}\sigma^2(\tau)\frac{1}{2(N_c-1)}\sum_{k1}^{N_c-1}(\bar{y}_{k1}-\bar{y}_k)^2\end{aligned} σ2(τ)2(Nc−1)1k1∑Nc−1(yˉk1−yˉk)2
Step 4 统计RMS值
根据Step 3求得的差值统计RMS值。
可选Step 5 绘制Allan方差曲线
改变块长度重复1~3并画出Allan标准差随块长度变化的双对数曲线称为Allan方差曲线。从图中可知每个时间块对应一个Allan方差将所有的不同块时间和对应的Allan方差连起来即绘制成Allan标准曲线。那么Allan方差分析就是分析Allan标准曲线。
4. Allan方差图读取误差系数
Allan方差法可用于5种随机误差的标定
量化噪声误差系数为 Q Q QAllan方差双对数曲线上斜率为-1的线段延长线与t1的交点的纵坐标读数为 3 T Q \frac{\sqrt{3}}{T}Q T3 Q角度随机游走其误差系数 N N NAllan方差双对数曲线上斜率为 − 1 / 2 -1/2 −1/2的线段延长线与 t 1 t1 t1交点的纵坐标读数即为 N T \frac{N}{\sqrt{T}} T N零偏不稳定性其误差系数 B B BAllan方差双对数曲线上斜率为0的线段延长线与t1交点的纵坐标读数为 2 ln 2 π B \sqrt{\frac{2\ln2}{\pi}}B π2ln2 B 一般常取底部平坦区的最小值或取 t 1 0 1 t10^1 t101或 t 1 0 2 t10^2 t102 处的值角速率随机游走其误差系数 K K K斜率为1/2的线段延长线与 t 1 t 1 t1交点的纵坐标读数为 K T / 3 \frac{K}{\sqrt{T/3}} T/3 K角速率斜坡其误差系数 R R R斜率为1的线段延长线与 t 1 t1 t1交点的纵坐标读数为 R T 2 \frac{RT}{\sqrt{2}} 2 RT。
假设各种误差源统计独立那总的艾伦方差为各种误差源之和即将量化噪声的平方 σ Q σ_Q σQ、角度随机游走的平方 σ R A W σ_{RAW} σRAW、零偏不稳定性的平方 σ b i a s σ_{bias} σbias、角速率随机游走的平方 σ R R W σ_{RRW} σRRW、角速率斜坡的平方 σ R R σ_{RR} σRR的总和。
5. Allan方差分析
通过块内求平均去除短时间的不确定性通过相邻块间求差去除长时间的不确定性。
如果我们只对某一时间尺度上的误差即误差波动感兴趣的话那么比这个时间尺度更小的细节变化短时间快速跳动和比这个时间尺度更大的宏观变化长时间缓慢漂移就都不关心了希望在我们的误差指标中都被消除掉。而Allan方差是这样做到的 通过分块确定所要考察的时间块长度 利用块内求平均的办法把短于块长度的那些快速变化成份细节都抹掉 再利用相邻块求差的办法把长于两块长度的那些缓慢变化成份宏观都抹掉 最后对差值序列统计其方值这是处理任何随机样本的标准操作这样统计出来的就是介于1倍块长度和2倍块长度这样一个很窄的时间尺度范围内的误差波动情况。
6. Allan方差曲线
如果我们对误差序列的各个时间尺度上的成份都感兴趣的话可以将块长度由短到长“扫描”一遍得出一组Allan方差值然后画个“Allan方差 vs. 块长度”的曲线这样就能全面反映被研究的误差序列的特性了。具体的实际上是“Allan方差的开方Allan Deviation vs. 块长度”的双对数曲线以便在更大的范围上有更强的表现力。以下是一张经典的高精度陀螺的Allan方差示意图。
从Allan方差曲线上我们可以根据曲线的形状特征来识别出不同类型的随机误差即随机过程模型并提取其参数。例如斜率为-1/2的直线段代表白噪声。具体可参见IEEE制定的一个关于陀螺测试的技术标准中的附录C内容[1]。这个附录中还给出了一个从Allan方差曲线中分析陀螺的多种随机误差的案例如上图。但我想提醒大家实际的Allan方差曲线往往只能表现出少数两三个主要误差类型因为其它误差都被这几个主要误差给淹没了。而由于我们在工程上关心的恰恰也就是主要误差源因此我们并不受此困扰。更多的使用Allan方差分析过程中的注意事项可参考西工大严恭敏老师的博文《Allan方差分析的使用要点》写得非常精准务实。
7. Allan方差的应用
Allan方差分析方法并不局限于对惯性器件和时钟误差的分析也可推广到其它误差序列分析而且它也不局限于是对误差序列进行分析也可应用于任何时间序列的分析例如建筑形变、地质演化等。
8. Allan方差与ROS
下面是读取bag包并获取Allen方差的图片会生成对应的Allen方差图。这里主要参考了allan_variance_ros 对应的C文件。
sudo pip install allantools#!/usr/bin/env python
import rospy
import sys
import allantools
import rosbag
import numpy as np
import csv
import rospkg
import os
import matplotlib.pyplot as plt # only for plotting, not required for calculations
import mathdef getRandomWalkSegment(tau,sigma):m -0.5 # slope of random walk Find point where slope -0.5 randomWalk Nonei 1idx 1mindiff 999logTau -999while (logTau0):logTau math.log(tau[i],10) logSigma math.log(sigma[i],10)prevLogTau math.log(tau[i-1],10)prevLogSigma math.log(sigma[i-1],10)slope (logSigma-prevLogSigma)/(logTau-prevLogTau)# 随机漫步的斜率diff abs(slope-m)# 当前斜率与目标斜率的差值if (diffmindiff):mindiff diff# 更新最小差值idx ii i 1 Project line to tau 10^0 x1 math.log(tau[idx],10)# 当前点的横坐标y1 math.log(sigma[idx],10)# 当前点的纵坐标x2 0y2 m*(x2-x1)y1return (pow(10,x1),pow(10,y1),pow(10,x2),pow(10,y2))def getBiasInstabilityPoint(tau,sigma):i 1while (itau.size):if (tau[i]1) and ((sigma[i]-sigma[i-1])0): # only check for tau 10^0breaki i 1return (tau[i],sigma[i])def main(args):rospy.init_node(allan_variance_node)t0 rospy.get_time() Parameters bagfile rospy.get_param(~bagfile_path,~/data/static.bag)# 输入的bag文件路径topic rospy.get_param(~imu_topic_name,/imu)# 输入的imu topic名称axis rospy.get_param(~axis,0)# 输入的轴0为x轴1为y轴2为z轴sampleRate rospy.get_param(~sample_rate,100)# 输入的采样频率isDeltaType rospy.get_param(~delta_measurement,False)# 是否是delta measurementnumTau rospy.get_param(~number_of_lags,1000)# 输入的tau数目resultsPath rospy.get_param(~results_directory_path,None)# 输出的结果路径 Results Directory Path if resultsPath is None:paths rospkg.get_ros_paths()path paths[1] # path to workspaces develidx path.find(ws/)# 找到路径workspacePath path[0:(idx3)]# 取到workspace的路径resultsPath workspacePath av_results/# 结果输出路径if not os.path.isdir(resultsPath):os.mkdir(resultsPath)print(\nResults will be save in the following directory: \n\n\t %s\n%resultsPath) Form Tau Array taus [None]*numTau# 初始化tau数组cnt 0;for i in np.linspace(-2.0, 5.0, numnumTau): # lags will span from 10^-2 to 10^5, log spacedtaus[cnt] pow(10,i)# 将tau数组赋值,维度在10^-2 到 10^5cnt cnt 1 Parse Bagfile bag rosbag.Bag(bagfile)N bag.get_message_count(topic) # 在bag文件中找到该topic的消息数量data np.zeros( (6,N) ) # 初始化数据矩阵维度为6*Nif isDeltaType:scale sampleRateelse:scale 1.0cnt 0for topic, msg, t in bag.read_messages(topics[topic]):# 遍历bag文件中的消息data[0,cnt] msg.linear_acceleration.x * scaledata[1,cnt] msg.linear_acceleration.y * scaledata[2,cnt] msg.linear_acceleration.z * scaledata[3,cnt] msg.angular_velocity.x * scaledata[4,cnt] msg.angular_velocity.y * scaledata[5,cnt] msg.angular_velocity.z * scalecnt cnt 1bag.close()print ([%0.2f seconds] Bagfile parsed\n%(rospy.get_time()-t0)) Allan Variance if axis is 0:currentAxis 1 # 循环通过所有轴1-6else:currentAxis axis # 只需循环一次然后中断while (currentAxis 6):# taus_used 对应了是否可以使用taus的数组adev是allan deviation degree的缩写为allan偏差adev_err是allan偏差的误差adev_norm是allan偏差的标准化;(taus_used, adev, adev_err, adev_n) allantools.oadev(data[currentAxis-1], data_typefreq, ratefloat(sampleRate), tausnp.array(taus) )# 计算allan variancerandomWalkSegment getRandomWalkSegment(taus_used,adev)# 计算random walk segmentbiasInstabilityPoint getBiasInstabilityPoint(taus_used,adev)# 计算bias instability pointrandomWalk randomWalkSegment[3]# 获取random walk segment的纵坐标biasInstability biasInstabilityPoint[1]# 获取bias instability point的纵坐标 Save as CSV if (currentAxis1):fname allan_accel_xtitle Allan Deviation: Accelerometer Xelif (currentAxis2):fname allan_accel_ytitle Allan Deviation: Accelerometer Yelif (currentAxis3):fname allan_accel_ztitle Allan Deviation: Accelerometer Zelif (currentAxis4):fname allan_gyro_xtitle Allan Deviation: Gyroscope Xelif (currentAxis5):fname allan_gyro_ytitle Allan Deviation: Gyroscope Yelif (currentAxis6):fname allan_gyro_ztitle Allan Deviation: Gyroscope Zprint ([%0.2f seconds] Finished calculating allan variance - writing results to %s%(rospy.get_time()-t0,fname))f open(resultsPath fname .csv, wt)try:writer csv.writer(f)writer.writerow( (Random Walk, Bias Instability) )writer.writerow( (randomWalk, biasInstability) )writer.writerow( (Tau, AllanDev, AllanDevError, AllanDevN) )for i in range(taus_used.size):writer.writerow( (taus_used[i],adev[i],adev_err[i],adev_n[i]) )finally:f.close() Plot Result plt.figure(figsize(12,8))ax plt.gca()ax.set_yscale(log)ax.set_xscale(log)plt.plot(taus_used,adev)plt.plot([randomWalkSegment[0],randomWalkSegment[2]],[randomWalkSegment[1],randomWalkSegment[3]],k--)plt.plot(1,randomWalk,rx,markeredgewidth2.5,markersize14.0)plt.plot(biasInstabilityPoint[0],biasInstabilityPoint[1],ro)plt.grid(True, whichboth)plt.title(title)plt.xlabel(Tau (s))plt.ylabel(ADEV)for item in ([ax.title, ax.xaxis.label, ax.yaxis.label] ax.get_xticklabels() ax.get_yticklabels()):item.set_fontsize(20)plt.show(blockFalse)plt.savefig(resultsPath fname)currentAxis currentAxis 1 axis*6 # increment currentAxis also break if axis is not 0inpinput(Press Enter key to close figures and end program\n)if __name__ __main__:main(sys.argv)
