Allan方差：标准allan方差求取和重叠(Overlaping)方差求取 & 如何计算各个参数值！

Allan介绍：

       对于随机误差，利用常规的分析方法，例如计算样本均值和方差并不能揭示潜在的误差源，另一方面，虽然自相关函数和功率谱密度函数分别从时域和频域描述了随机误差的统计特性，但是在实际工作中通过这些函数加以分析将随机误差分离出来是很困难的；Allan方差法是20世纪60年代由美国国家标准局David Allan提出的，它是一种基于时域的分析方法，不仅可以用来分析光学陀螺的误差特性，而且还可以应用于其他任何精密测量仪器.Allan方差法的主要特点是能非常容易对各种误差源及其对整个噪声统计特性的贡献进行细致的表征和辨识，而且便于计算，易于分离。它提供了一种识别并量化存在于数据中的不同噪声项。

       Allan方差分析主要用于分析陀螺量化噪声、角度随机游走噪声（角速率白噪声）、零偏不稳定性、角速率随机游走（角加速度白噪声）和速率斜坡（角速率常值趋势项），等五种典型误差。

其中在组合导航中，主要有用的是

（1）、角度随机游走系数（角速率白噪声）：主要用于设置Q阵
（2）、零偏不稳定性系数（Bias Instability）:主要用于设置一阶马氏过程的方差

1、标准allan方差计算法

定义与计算：

（1）基于角速率测量值的ALLAN计算，此处以陀螺仪输出为例：

设以采样时间 τ0 对陀螺仪输出角速率进行采样，共采样N个点，把所获得的N个数据分成K组，每组包含M个采样点。

 K=N/M，M<（N-1）/2；如下：

每一组的持续时间 τM=Mτ0 ，称之为相关时间，每一组的平均值为:

Allan方差定义为:

备注：此处的方差计算与一般我们所熟知的方差（概率统计中）计算不同，它每个样本不是与均值相减，而是每个区间的均值与下一个区间的均值相减（这里的区间就是分组，即组与组之间相减）。

（2）基于角度增量测量值的ALLAN计算，此处以陀螺仪输出为例：

陀螺仪直接输出测量时刻的角度增量值：

设采样时间为 τ0 ，则角度增量测量值在离散时刻为tk=kτ0(k=1,2,……，N)上进行的；简记为θk=θ（kτ0 ，时刻tk与tk+τ之间的角速率为：

式中，τ=m*τ0，则有角度增量测量值定义的ALLAN方差为：

随机噪声的ALLAN方差分析结果图为：

ALLAN方差的估计精度：

在实际中，ALLAN方差估计是基于有限长度数据，估计精度依赖于独立数组（分组个数）的数量。设共有N个数据点，将其分成长度为M个数据组，则ALLAN方差估计的百分比误差为：

上式表明：

在短τ中估算误差：即分组数量很多（M比较大），因此估算误差较小；

在长τ中估算误差：即分组数量很少（M比较小），因此估算误差较大；

从ALLAN双对数图中，也可以看出，图中左部分τ较小，allan估计误差较小；而在图的右部分τ较大，allan估计误差较大！

现在将标准ALLAN方差求解的程序附上：

clear ; clc;close all;
% function [sigma, tau, Err] = avar(y0, tau0)
% 计算Allan方差
% 输入：y -- 数据(一行或列向量), tau0 -- 采样周期
% 输出：sigma -- Allan方差(量纲单位与输入y保持一致), tau -- 取样时间, Err -- 百分比估计误差
% 作者: Yan Gong-min, 2012-08-22
% example: 
    y0 = randn(100000,1) + 0.00001*[1:100000]';
    tau0=0.1;
%     [sigma, tau, Err] = avar(y, 0.1);
%  [sigma, tau, Err] = avar(y0, tau0)
    N = length(y0);
    y = y0; NL = N;
    for k = 1:100000
        sigma(k,1) = sqrt(1/(2*(NL-1))*sum([y(2:NL)-y(1:NL-1)].^2));
        tau(k,1) = 2^(k-1)*tau0;
        Err(k,1) = 1/sqrt(2*(NL-1));
        NL = floor(NL/2);
        if NL<3
            break;
        end
        y = 1/2*(y(1:2:2*NL) + y(2:2:2*NL));  % 分组长度加倍(数据长度减半)
    end
    subplot(211), plot(tau0*[1:N], y0); grid
    xlabel('\itt \rm/ s'); ylabel('\ity');
    subplot(212), 
    loglog(tau, sigma, '-+', tau, [sigma.*(1+Err),sigma.*(1-Err)], 'r--'); grid
    xlabel('\itt \rm/ s'); ylabel('\it\sigma_A\rm( \tau )');

运行结果：

参考链接：

Allan方差分析的使用要点

Analysis and Modeling of Inertial Sensors Using Allan Variance

IMU Noise Model

2、重叠方差计算法(Overlapping Allan)

重叠式和标准式的不同在于提高有限数据的利用率，具体如下图所示：

Cluster中为3*τ；每一个Cluster之间的距离为τ;

示例代码（基于角度增量计算）：

clear ; clc;close all;
omega = randn(100000,1) + 0.00001*[1:100000]';
fs=10; %频率
pts=100; %组
% [T,sigma]=allan(omega,fs,pts)
[N ,M]=size(omega); %计算数据集大小
n=2.^(0:floor(log2(N/2)))'; %最大数据集
maxN=n(end);
endLogInc=log10(maxN);
m=unique(ceil(logspace(0,endLogInc,pts)))';
t0=1/fs;
T=m*t0;
theta=cumsum(omega)/fs;  %速率值变为角度增量
sigma2=zeros(length(T),M);
for i=1:length(m)
    for k=1:N-2*m(i)
        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);
%plot
loglog(T, sigma, '-+'); grid
xlabel('\itt \rm/ s'); ylabel('\it\sigma_A\rm( \tau )');

运行结果：

然后根据ALLAN方差分析结果图求取在标准单位下（SI）的噪声值：

量化噪声：slope=-1

角度随机游走：slope=-1/2

零偏不稳定性：slope=0(曲线的平直部分，最下端)

角速率随机游走：slope=1/2

速率斜坡：slope=1;

以MEMS陀螺仪为例，其噪声主要是由： N (angle random walk), K (rate random walk), B (bias instability)

示例代码：

% Noise Parameter Identification
% Find the index where the slope of the log-scaled Allan deviation is equal to the slope specified.
% Angle Random Walk
% The above equation is a line with a slope of -1/2 when plotted on a log-log plot of σ(τ)versusτ. 
% The value of N can be read directly off of this line at τ=1.
slope = -0.5;
logtau = log10(tau);
logadev = log10(adev);
dlogadev = diff(logadev) ./ diff(logtau);
[~, i] = min(abs(dlogadev - slope)); %找到slope=-0.5 最接近的点；

% Find the y-intercept of the line.
b = logadev(i) - slope*logtau(i);
% Determine the angle random walk coefficient from the line.
logN = slope*log(1) + b;
disp('random walks')
N = 10^logN %角度随机游走值，单位为：(rad/s/root-Hz)

% Plot the results.
tauN = 1;
lineN = N ./ sqrt(tau);
figure
loglog(tau, adev, tau, lineN, '--', tauN, N, 'o')
title('Allan Deviation with Angle Random Walk')
xlabel('\tau')
ylabel('\sigma(\tau)')
legend('\sigma', '\sigma_N')
text(tauN, N, 'N')
grid on
axis equal

%% Rate Random Walk
% The above equation is a line with a slope of 1/2 when plotted on a log-log plot of σ(τ)versusτ. 
% The value of N can be read directly off of this line at τ=3.
slope = 0.5;
logtau = log10(tau);
logadev = log10(adev);
dlogadev = diff(logadev) ./ diff(logtau);
[~, i] = min(abs(dlogadev - slope));

% Find the y-intercept of the line.
b = logadev(i) - slope*logtau(i);

% Determine the rate random walk coefficient from the line.
logK = slope*log10(3) + b;
disp('rate random walks')
K = 10^logK %速率随机游走 angle rate random walks (rad/s^2/root-Hz)

% Plot the results.
tauK = 3;
lineK = K .* sqrt(tau/3);
figure
loglog(tau, adev, tau, lineK, '--', tauK, K, 'o')
title('Allan Deviation with Rate Random Walk')
xlabel('\tau')
ylabel('\sigma(\tau)')
legend('\sigma', '\sigma_K')
text(tauK, K, 'K')
grid on
axis equal

%% Bias Instability
%  The above equation is a line with a slope of 0 when plotted on a log-log plot of σ(τ)versusτ. 
% The value of B can be read directly off of this line with a scaling of sqrt(2*log(2)/pi) .
slope = 0;
logtau = log10(tau);
logadev = log10(adev);
dlogadev = diff(logadev) ./ diff(logtau);
[~, i] = min(abs(dlogadev - slope));

% Find the y-intercept of the line.
b = logadev(i) - slope*logtau(i);

% Determine the bias instability coefficient from the line.
scfB = sqrt(2*log(2)/pi);
logB = b - log10(scfB);
disp('Bias Instability')
B = 10^logB %零偏不稳定性 gyros dynamic biases or bias instabilities (radians/s)
disp('BI, correlation time, in seconds\n')
tau(i)

% Plot the results.
tauB = tau(i);
lineB = B * scfB * ones(size(tau));
figure
loglog(tau, adev, tau, lineB, '--', tauB, scfB*B, 'o')
title('Allan Deviation with Bias Instability')
xlabel('\tau')
ylabel('\sigma(\tau)')
legend('\sigma', '\sigma_B')
text(tauB, scfB*B, '0.664B')
grid on
axis equal


%% Now that all the noise parameters have been calculated, 
% plot the Allan deviation with all of the lines used for quantifying the parameters.
tauParams = [tauN, tauK, tauB];
params = [N, K, scfB*B];
figure
loglog(tau, adev, tau, [lineN, lineK, lineB], '--', ...
    tauParams, params, 'o')
title('Allan Deviation with Noise Parameters')
xlabel('\tau')
ylabel('\sigma(\tau)')
legend('\sigma', '\sigma_N', '\sigma_K', '\sigma_B')
text(tauParams, params, {'N', 'K', '0.664B'})
grid on
axis equal

结果输出：

如有疑问，请和我沟通！