0 原理解释
http://blog.sciencenet.cn/blog-3392538-1123730.html
在Allan方差分析中，共有5个噪声参数：量化噪声、角度随机游走、零偏不稳定性、速度随机游走和速度爬升。不同段的Allan方差曲线代表了不同的误差参数，我们要求解的 零偏噪声(Bias) 对应的曲线段的 斜率就是0。
Allan方差分析方法可以估计由频率漂移、温度误差、过程噪声等引起的系统误差或缺陷的时域信号的震荡稳定性。
1. 概述
在上一篇文章，我们对器件误差的成分进行了分析，并且介绍了最主要的五种误差成分在“方差-时间间隔”的双对数曲线上的直观表现及产生这种表现的原理。同时也提到，在两种误差分析方法中，Allan方差方法比直接方差统计要更有优势，所以本篇文章我们就一起讨论下Allan方差的具体实现方式。
2. 实现方式
Allan方差识别误差的方式是画“方差-时间间隔”双对数曲线，这包含三方面：方差、时间间隔、双对数(好像是一句废话，哈哈)。在测试中，我们按照一定的周期T采集了一段时间数据，根据这些离散的数据，便能计算其方差。这个周期T和方差在二维坐标系中便构成了一个点的横纵坐标。为了得到曲线，我们需要得到更多这样的点，把周期缩短是不可能了，因为不能凭空提高已有数据频率，所以只能把周期拉长。
如果我们对每两个相邻的数取平均，合并成一个，便构成了周期为2T的离散数据，同样也可以得到它的方差，也即又多得到了一个点。按照这样的方法，循环下去，可以得到周期为4T、8T、16T对应的点，直到数据少到无法再合并为止，这样就得到系列的点，也就是曲线了。这样做出来的曲线，其横纵坐标都会很长，各个点之间的间隔也极其不均匀，这也就是为什么用双对数曲线的原因了。
3. 参数拟合
我们先给出一张Allan方差曲线图，对着图来说吧。
在理论上，图中的曲线是由五个直线组成的，五个直线的参数分别包含了我们想要提取的五项误差。所以，只要把直线公式拟合出来，那么参数自然就知道了。我们把上一篇文章的公式再放过来看一下。
从公式上可以看出，我们只要找到直线在任一横坐标处的值，便可提取参数，方便起见，一般取横坐标为1处的值，因为这样，公式右边和时间相关的项就是零。
总结：
Allan方差分析的一个用途是分析陀螺的性能或对比不同陀螺的性能， 相比于其它分析法Allan法还是很好用的，比较全面。另一个用途是获得噪声参数，用于组合导航的Kalman滤波噪声参数设置。不是所有的Allan方差噪声系数都有用，主要有用的是角度随机游走系数（用于设置Q阵）和零偏不稳定性系数（用于设置一阶马氏过程的方差），其实这两个系数量级大小差不多就行了，太精细也没用。毕竟Allan方差分析的只是一个样本，你再采一个样本试试，肯定会有差异的；此外采集的是陀螺的静态性能数据，鬼知道实际动态应用的时候陀螺误差会变化多大，存在数量级差别都很有可能。这里需要的是理论指导下的工程和经验，而不是教条主义
4. 代码实践
我们使用开源系统imu_utils为例，来看提取参数的具体实现方式。
开源代码地址：https://github.com/gaowenliang/imu_utils
1) 代码思路
以上我们拟合直线，再求交点的方式，适合人工操作，按这种方式，给出一个曲线图，很快就可以计算出各误差参数。而对于机器来讲，则似乎不太适合，首先让它把曲线自适应地分成五段就要费一番功夫，此时就需要换一种思路。既然曲线是五条直线组成，而五条直线的公式已知，那么五条直线加一起不就是曲线的公式了吗，有了曲线公式，又有了曲线实际数据，那么这就变成了一个曲线拟合问题，用优化的方式就可以直接解决。比较来讲，这两种方式无好坏之分，只能说一个更适合人工，一个更适合机器，为各自量身打造。
2) 代码细节
首先是把原始数据划分成不同的时间间隔，思路就是上面提到的倍增法
int mode = numData / 2;
unsigned int maxStride = 1;
int shft =0;
while(mode){
mode = mode >> 1;
maxStride = 1 << shft;
shft++;
}然后每一个间隔下都计算一个方差，并把方差和真实曲线建立残差模型。
for ( int i = 0; i < num_samples; ++i ) {
ceres::CostFunction* f = new ceres::AutoDiffCostFunction< AllanSigmaError, 1, 5 > (new AllanSigmaError( sigma2s_tmp[i], m_taus[i] ) );
problem.AddResidualBlock( f, NULL, param );
}其中残差的计算如下:
// 计算残差
template< typename T >
bool operator( )( const T* const _paramt, T* residuals ) const {
T _Q = T( _paramt[0] );
T _N = T( _paramt[1] );
T _B = T( _paramt[2] );
T _K = T( _paramt[3] );
T _R = T( _paramt[4] );
T _tau = T( tau );
T _sigma2 = calcSigma2( _Q, _N, _B, _K, _R, _tau );
T _dsigma2 = T( calcLog10( _sigma2 ) ) - T( calcLog10( sigma2 ) );
residuals[0] = _dsigma2;
return true;
}
// 叠加五个直线，得到曲线方程
template< typename T >
T calcSigma2( T _Q, T _N, T _B, T _K, T _R, T _tau ) const {
return _Q * _Q / ( _tau * _tau )
+ _N * _N / _tau
+ _B * _B
+ _K * _K * _tau
+ _R * _R * _tau * _tau;
}最后，通过迭代优化，就可以直接计算出那五个参数了。
5. MATLAB程序 （转载于：https://blog.csdn.net/u012325601/article/details/60882949）
博客所对应的ALLAN方差的matlab源码地址为：https://github.com/XinLiGH/GyroAllan
allan.m文件
function [T,sigma] = allan(omega,fs,pts)
[N,M] = size(omega); % figure out how big the output data set is
n = 2.^(0:floor(log2(N/2)))'; % determine largest bin size
maxN = n(end);
endLogInc = log10(maxN);
m = unique(ceil(logspace(0,endLogInc,pts)))'; % create log spaced vector average factor
t0 = 1/fs; % t0 = sample interval
T = m*t0; % T = length of time for each cluster
theta = cumsum(omega)/fs; % integration of samples over time to obtain output angle θ
sigma2 = zeros(length(T),M); % array of dimensions (cluster periods) X (#variables)
for i=1:length(m) % loop over the various cluster sizes
for k=1:N-2*m(i) % implements the summation in the AV equation
sigma2(i,:) = sigma2(i,:) + (theta(k+2*m(i),:) - 2*theta(k+m(i),:) + theta(k,:)).^2;
end
end
sigma2 = sigma2./repmat((2*T.^2.*(N-2*m)),1,M);
sigma = sqrt(sigma2);nihe.m文件
function C=nihe(tau,sig,M)
X=tau';Y=sig';
B=zeros(1,2*M+1);
F=zeros(length(X),2*M+1);
for i=1:2*M+1
kk=i-M-1;
F(:,i)=X.^kk;
end
A=F'*F;
B=F'*Y;
C=A\B;gyro_data.m文件
clear
clc
close all
% --------------------.bin文件得读取方式 ---------------------%
fip=fopen('201214204200_imu.bin','rb');
[data,num]=fread(fip,[7 624300],'double');%inf表示读取文件中的所有数据，[M,N]表
fclose(fip)
data = data';
% --------------------.txt文件得读取方式 ---------------------%
% data = importdata("imu.txt"); %读取名字为 imu.txt 的文件
datanumber = 624300; % 注释后面的，先看一下自己的数据大小
IMU_raw_data = zeros(datanumber,7);
IMU_raw_gyro_data = zeros(datanumber,3);
for i = 1:1:datanumber
IMU_raw_data(i,1:7) = data(i,1:7)*200; %（1）IMU增量变成角速度; 否则不用*200
IMU_raw_gyro_data(i,1:3) = IMU_raw_data(i,2:4)*3600; %（2）把 deg/s 转为 deg/h
end
[A, B] = allan(IMU_raw_gyro_data, 100, 100); %求Allan标准差，用100个点来描述
loglog(A, B, 'o'); %画双对数坐标图
xlabel('time:sec'); %添加x轴标签
ylabel('Sigma:deg/h'); %添加y轴标签
legend('X axis','Y axis','Z axis'); %添加标注
grid on; %添加网格线
hold on; %使图像不被覆盖
C(1, :) = nihe(A', (B(:,1)').^2, 2)'; %拟合
C(2, :) = nihe(A', (B(:,2)').^2, 2)';
C(3, :) = nihe(A', (B(:,3)').^2, 2)';
Q = sqrt(abs(C(:, 1) / 3)); %量化噪声，单位：deg/h
N = sqrt(abs(C(:, 2) / 1)) / 60; %角度随机游走，单位：deg/h^0.5
Bs = sqrt(abs(C(:, 3))) / 0.6643; %零偏不稳定性，单位：deg/h
K = sqrt(abs(C(:, 4) * 3)) * 60; %角速率游走，单位：deg/h/h^0.5
R = sqrt(abs(C(:, 5) * 2)) * 3600; %速率斜坡，单位：deg/h/h
fprintf('量化噪声 X轴：%f Y轴：%f Z轴：%f 单位：deg/h\n', Q(1), Q(2), Q(3));
fprintf('角度随机游走 X轴：%f Y轴：%f Z轴：%f 单位：deg/h^0.5\n', N(1), N(2), N(3));
fprintf('零偏不稳定性 X轴：%f Y轴：%f Z轴：%f 单位：deg/h\n', Bs(1), Bs(2), Bs(3));
fprintf('角速率游走 X轴：%f Y轴：%f Z轴：%f 单位：deg/h/h^0.5\n', K(1), K(2), K(3));
fprintf('速率斜坡 X轴：%f Y轴：%f Z轴：%f 单位：deg/h/h\n', R(1), R(2), R(3));
D(:, 1) = sqrt(C(1,1)*A.^(-2) + C(1,2)*A.^(-1) + C(1,3)*A.^(0) + C(1,4)*A.^(1) + C(1,5)*A.^(2)); %生成拟合函数
D(:, 2) = sqrt(C(2,1)*A.^(-2) + C(2,2)*A.^(-1) + C(2,3)*A.^(0) + C(2,4)*A.^(1) + C(2,5)*A.^(2));
D(:, 3) = sqrt(C(3,1)*A.^(-2) + C(3,2)*A.^(-1) + C(3,3)*A.^(0) + C(3,4)*A.^(1) + C(3,5)*A.^(2));
loglog(A, D); %画双对数坐标图运行结果 （放置10个小时左右）
量化噪声 X轴：0.003226 Y轴：0.004128 Z轴：0.001843 单位：deg/h
角度随机游走 X轴：0.003048 Y轴：0.004053 Z轴：0.004220 单位：deg/h^0.5
零偏不稳定性 X轴：0.273572 Y轴：0.353873 Z轴：0.277789 单位：deg/h
角速率游走 X轴：2.661679 Y轴：5.058015 Z轴：3.378206 单位：deg/h/h^0.5
速率斜坡 X轴：1.636452 Y轴：4.985344 Z轴：7.571720 单位：deg/h/h
注意：当Allan标准差的拟合多项式中的拟合系数是负值时，所得误差项的拟合结果随着相关时间的微小改变变化很大，拟合误差很大，可信度差。
参考文献下载地址如下
链接：https://pan.baidu.com/s/1TTzZywJb1tdaVWuI-XMclw 密码：l7p4
6. 一些实践经验
至此为止，我们应该已经弄清了器件误差提取方法，然而，实践总是会给人惊喜，出现一些意料之外，而且我们也不必完全拘泥于理论，要看我们的目的是什么。
1) 有的时候不一定五个斜率全都能体现出来，我们知道，每个斜率对应一项误差，如果某项误差相对于其他项的量级比较小，那么曲线中就没有他什么事了，直接被掩盖，极端情况下，一个曲线只有一个斜率也是有可能的。
2) Allan方差对误差项提取的比较多，但是很多误差项的分析都是为了提高器件生产工艺用的，每种误差成分产生的原理不同，通过分析哪种误差成分为主，我们就知道为了提高精度，应该重点改进哪项工艺。如果只是为传感器融合找个方差的参数，大可不必计算那么细致，在图上看一下大致量级，就可以直接用了，因为参数早晚都是要调的。严格来讲，滤波的方差其实应该用白噪声强度，它不属于这五项误差的任何一项，但是很多时候手册上不给出这个强度值，而方差参数在用的时候都会结合系统再一次进行调整，所以有些时候直接用零偏不稳定性做方差初值也不是不可以。
3) 随机游走这种东西理论上是可以建模，但是实际滤波模型中却很少使用，甚至不建议使用。一是、因为它并不是构成零偏趋势项的主要成分，往往都被温变等的趋势掩盖，即使做了温补也是这样。二是、因为Allan方差识别出的参数并不完全准确，这会带来模型误差，不准确的模型还不如没有模型，因为它会使滤波变得不稳定，效果反而变差。
758
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
