allan方差 imu_MEMS IMU的噪聲和誤差源2(allan方差)

Allan 方差簡介

Allan方差分析主要用于分析陀螺量化噪聲，一般用于分析(以下内容針對于陀螺儀的allan方差分析，對于加速度計都是類似的)量化噪聲(一般不一定能從allan上看出來)

角度随機遊走

零偏不穩定性

角速率随機遊走

速率斜坡(一般不一定能從allan上看出來)

Allan方差分析的一個用途是分析陀螺的性能或者對比不同陀螺的性能。若是将幾個陀螺的Allan方差曲線同時繪制一個雙對數圖上，曲線在左下角的對應的陀螺性能一般比較好一些(雷射陀螺)，而右上差(MEMS)。

斜率意義：

matlab實踐環境：(matlab2019， 安裝sensor fusion toolbox)，主要參考matlab自帶的 AllanVarianceExample.m

其實allan方差計算并不複雜，下面的例子隻針對陀螺，把陀螺模組化為：

其中我們隻關心:N: angle random walk(-1/2斜率)

K: rate random walk(1/2斜率)

B: bias instability(0斜率)

理由引嚴老師的原文：Allan方差可用于分析五種典型誤差，但是某個陀螺中并不一定五種誤差都有同時表現出來，不畫Allan方差圖/不肉眼看圖，而盲目的采用最小二乘回歸方法強行地就想求解五個系數的人是典型的教條主義者，往往得不到好的結果。比如如下Allan方差圖中y曲線所示(摘自https://github.com/ethz-asl/kalibr/wiki/IMU-Noise-Model)，采用“三段論”大概讀圖計算角度随機遊走系數、零偏不穩定性和角速率随機遊走(或速率斜坡)即可，完全可滿足實際應用。過細分析想得出精确系數的人都是有潔癖的。

手撸一個allan曲線記錄靜态條件下陀螺資料，總長度為L, 采樣周期為

,

為所有采樣記錄的樣本， 這裡

,采樣率=100。 這個matlab自帶的資料集的陀螺的機關為rad/s

% Load logged data from one axis of a three-axis gyroscope. This recording

% was done over a six hour period with a 100 Hz sampling rate.

load('LoggedSingleAxisGyroscope', 'omega', 'Fs')

t0 = 1/Fs;

2. 計算每一時刻下的累積角度:

theta = cumsum(omega, 1)*t0;

3. 計算allan方差

maxNumM = 100;

L = size(theta, 1);

maxM = 2.^floor(log2(L/2));

m = logspace(log10(1), log10(maxM), maxNumM).';

m = ceil(m); % m must be an integer.

m = unique(m); % Remove duplicates.

tau = m*t0;

avar = zeros(numel(m), 1);

for i = 1:numel(m)

mi = m(i);

avar(i,:) = sum( ...

(theta(1+2*mi:L) - 2*theta(1+mi:L-mi) + theta(1:L-2*mi)).^2, 1);

end

avar = avar ./ (2*tau.^2 .* (L - 2*m));

4. 最後計算

并作圖：

adev = sqrt(avar);

figure

loglog(tau, adev)

title('Allan Deviation')

xlabel('\tau');

ylabel('\sigma(\tau)')

grid on

axis equal

這樣allan方差就畫出來了，一定要注意Y軸的機關，1.是方差，還是标準差？ 2.機關有的人用rad/s,有的人用deg/s, 更多的人用rad/h

使用matlab自帶的allan函數

自己手撸太麻煩，用起來也心驚膽戰的，不如拿來主義，直接用matlab自帶的allan函數allanvar。

[avarFromFunc, tauFromFunc] = allanvar(omega, m, Fs);

adevFromFunc = sqrt(avarFromFunc);

figure

loglog(tau, adev, tauFromFunc, adevFromFunc);

title('Allan Deviations')

xlabel('\tau')

ylabel('\sigma(\tau)')

legend('Manual Calculation', 'allanvar Function')

grid on

axis equal

這是自己撸和使用自帶allanvar的曲線，一樣得嘛！

使用allan曲線計算零篇穩定性，角度随機遊走和角速率随機遊走

這個地方就不展開了，我直接封裝了一個matlab函數ch_allan, 其中調用allanvar 并且自動拟合直線計算N,K,B值, 注意輸入機關必須是rad/s

%% ('NoiseDensity(角度随機遊走)', N, 'RandomWalk(角速率随機遊走)', K,'BiasInstability(零偏穩定性)', B);

% omega 必須為 rad/s

% see Freescale AN: Allan Variance: Noise Analysis for Gyroscopes

function[avar, tau, N, K, B] =ch_allan(omega, Fs, is_plot)[avar, tau] = allanvar(omega, 'octave', Fs);

adev = sqrt(avar);

%% Angle Random Walk

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 angle random walk coefficient from the line.

logN = slope*log(1) + b;

N = 10^logN;

% Plot the results.

tauN = 1;

lineN = N ./ sqrt(tau);

%% Rate Random Walk

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;

K = 10^logK;

% Plot the results.

tauK = 3;

lineK = K .* sqrt(tau/3);

%% Bias Instability

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);

B = 10^logB;

% Plot the results.

tauB = tau(i);

lineB = B * scfB * ones(size(tau));

if is_plot == true

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

end

end

使用執行個體，這裡用gnss-ins-sim 這個仿真庫産生靜态陀螺資料，然後用ch_allan來分析，并對比allan出來的參數是否和生成仿真資料時設定的參數一緻

clear;

clc;

close all;

%% https://github.com/Aceinna/gnss-ins-sim

% VRW 機關: m/s/sqrt(hr) 'accel_vrw': np.array([0.3119, 0.6009, 0.9779]) * 1.0,

% 加速度零偏不穩定性: m/s^(2) 'accel_b_stability': np.array([1e-3, 3e-3, 5e-3]) * 1.0,

% ARW 機關: deg/sqrt(hr) 'gyro_arw': np.array([0.25, 0.25, 0.25]) * 1.0,

% 角速度零偏穩定性: deg/h 'gyro_b_stability': np.array([3.5, 3.5, 3.5]) * 1.0,

M = csvread('gyro-0.csv',1 ,0);

gyroReading = M(:,1);

gyroReading = deg2rad(gyroReading);

M = csvread('accel-0.csv',1 ,0);

Fs = 100;

accelReading = M(:,1);

%% 運作陀螺 allan

[avar1, tau1 , N, K, B] = ch_allan(gyroReading , Fs, true);

fprintf('GYR: 零偏不穩定性 %frad/s 或 %fdeg/h \n', B, rad2deg(B)*3600);

fprintf('GYR: 噪聲密度(角度随機遊走, ARW, Noise density) %f(rad/s)/sqrt(Hz) 或 %f deg/sqrt(h)\n', N, rad2deg(N)*3600^(0.5));

fprintf('GYR: 角速率随機遊走, bias variations ("random walks") %f(rad/s)sqrt(Hz) 或 %f deg/h/sqrt(h)\n', K, rad2deg(K) * (3600^(1.5)));

%% 運作加速度計 allan

[avar1, tau1 , N, K, B] = ch_allan(accelReading, Fs, true);

fprintf('ACC: 零偏不穩定性 %fm/s^(2) 或 %fmg 或 %fug\n', B, B / 9.80665 *1000, B / 9.80665 *1000*1000);

fprintf('ACC: 噪聲密度(速率随機遊走,VRW, Noise Density, Rate Noise Density) %f(m/s^(2))/sqrt(Hz) 或 %f m/s/sqrt(hr)\n', N, N * 3600^(0.5));

fprintf('ACC: 加速度随機遊走，bias variations ("random walks") %f(m/s^(2)sqrt(Hz))\n', K);

運作結果：(和仿真時設定的參數基本相同)

參考