Инерционный анализ шума датчика Используя отклонение Аллана

Открыть скрипт

В этом примере показано, как использовать отклонение Аллана, чтобы определить шумовые параметры гироскопа MEMS. Эти параметры могут использоваться, чтобы смоделировать гироскоп в симуляции. Измерение гироскопа моделируется как:

$$\Omega(t) = \Omega_{Ideal}(t) + Bias_N(t) + Bias_B(t) + Bias_K(t)$$

Три шумовых параметра N (угол случайный обход), K (случайное блуждание по скорости) и B (нестабильность смещения) оцениваются данные об использовании, регистрируемые от стационарного гироскопа.

Фон

Отклонение Аллана было первоначально разработано Дэвидом В. Алланом, чтобы измерить устойчивость частоты генераторов точности. Это может также использоваться, чтобы идентифицировать различные источники шума, существующие в стационарных измерениях гироскопа. Рассмотрите выборки L данных из гироскопа с шагом расчета$\tau_{0}$. Сформируйте кластеры данных длительности$\tau_{0}$$2\tau_{0}$...,$m\tau_{0}, (m < (L-1)/2)$ и получите средние значения суммы точек данных, содержавшихся в каждом кластере по длине кластера. Отклонение Аллана задано как 2D демонстрационное отклонение средних значений кластера данных как функция кластерного времени. Этот пример использует перекрывающееся средство оценки отклонения Аллана. Это означает, что расчетные кластеры перекрываются. Средство оценки выполняет лучше, чем неперекрывающиеся средства оценки для больших значений L.

Вычисление отклонения Аллана

Отклонение Аллана вычисляется можно следующим образом:

Регистрируйте стационарные выборки гироскопа L с периодом расчета$\tau_{0}$. Позвольте$\Omega$ быть регистрируемыми выборками.

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

Для каждой выборки вычислите выходной угол$\theta$:

$$\theta(t) = \int^{t}\Omega(t')dt'$$

Для дискретных выборок$\tau_{0}$ может использоваться совокупная сумма, умноженная на.

theta = cumsum(omega, 1)*t0;

Затем вычислите отклонение Аллана:

$$\sigma^2(\tau) =&#xA;\frac{1}{2\tau^2}<(\theta_{k+2m}-2\theta_{k+m}+\theta_{k})^2&#62;$$

где$\tau = m\tau_{0}$ и$<&#62;$ среднее значение ансамбля.

Среднее значение ансамбля может быть расширено до:

$$\sigma^2(\tau) =&#xA;\frac{1}{2\tau^2(L-2m)}\sum_{k=1}^{L-2m}(\theta_{k+2m} - 2\theta_{k+m}&#xA;+ \theta_{k})^2$$

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

Наконец, отклонение Аллана$\sigma(t) = \sqrt{\sigma^2(t)}$ используется, чтобы определить параметры шума гироскопа.

adev = sqrt(avar);

figure
loglog(tau, adev)
title('Allan Deviation')
xlabel('\tau');
ylabel('\sigma(\tau)')
grid on
axis equal

Отклонение Аллана может также быть вычислено с помощью 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

Шумовая идентификация параметра

Чтобы получить шумовые параметры для гироскопа, используйте следующее отношение между отклонением Аллана и двухсторонней степенью спектральной плотностью (PSD) шумовых параметров в исходном наборе данных$\Omega$. Отношение:

$$\sigma^2(\tau) = 4\int_{0}^{\infty}S_\Omega(f)&#xA;\frac{sin^4(\pi f\tau)}{(\pi f\tau)^2}df$$

От вышеупомянутого уравнения отклонение Аллана пропорционально общей шумовой степени гироскопа, когда прошли фильтр с передаточной функцией$sin^4(x)/(x)^2$. Эта передаточная функция является результатом операций, сделанных, чтобы создать и работать с кластерами.

Используя эту интерпретацию передаточной функции, полоса пропускания фильтра зависит от$\tau$. Это означает, что различные шумовые параметры могут быть идентифицированы путем изменения полосы пропускания фильтра или варьирования$\tau$.

Угол случайный обход

Угол случайный обход характеризуется белым шумовым спектром гироскопа выход. PSD представлен:

$$S_\Omega(f) = N^2$$

где

N = угол случайный коэффициент обхода

Замена в исходное уравнение PSD и выполнение выражений интегрирования:

$$\sigma^2(\tau) = \frac{N^2}{\tau}$$

Вышеупомянутое уравнение является линией с наклоном-1/2, когда построено на графике логарифмического журнала$\sigma(\tau)$ по сравнению с$\tau$. Значение N может быть считано непосредственно прочь этой линии в$\tau = 1$.

% Find the index where the slope of the log-scaled Allan deviation is equal
% to the slope specified.
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);
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

N =

    0.0126

Случайное блуждание по скорости

Случайное блуждание по скорости характеризуется красным шумом (Броуновский шум) спектр гироскопа выход. PSD представлен:

$$S_\Omega(f) = (\frac{K}{2\pi})^2\frac{1}{f^2}$$

где

K = коэффициент случайного блуждания по скорости

Замена в исходное уравнение PSD и выполнение выражений интегрирования:

$$\sigma^2(\tau) = \frac{K^2\tau}{3}$$

Вышеупомянутое уравнение является линией с наклоном 1/2, когда построено на графике логарифмического журнала$\sigma(\tau)$ по сравнению с$\tau$. Значение K может быть считано непосредственно прочь этой линии в$\tau = 3$.

% Find the index where the slope of the log-scaled Allan deviation is equal
% to the slope specified.
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);
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

K =

   9.0679e-05

Сместите нестабильность

Нестабильность смещения характеризуется розовым шумом (шум мерцания) спектр гироскопа выход. PSD представлен:

$$S_{\Omega}(f) = \left\{\begin{array}{lr}(\frac{B^2}{2\pi})\frac{1}{f}&#xA;&#38; : f \leq f_0\\0 &#38; : f &#62; f_0\end{array}\right.$$

где

B = сместите коэффициент нестабильности

$f_0$ = частота среза

Замена в исходное уравнение PSD и выполнение выражений интегрирования:

$$\sigma^2(\tau) = \frac{2B^2}{\pi}[\ln{2} + \\&#xA;-\frac{sin^3x}{2x^2}(sinx + 4xcosx) + Ci(2x) - Ci(4x)]$$

где

$x = \pi f_0\tau$

Ci = функция интегрального косинуса

Когда$\tau$ намного более длинно, чем инверсия частоты среза, уравнение PSD:

$$\sigma^2(\tau) = \frac{2B^2}{\pi}\ln{2}$$

Вышеупомянутое уравнение является линией с наклоном 0, когда построено на графике логарифмического журнала$\sigma(\tau)$ по сравнению с$\tau$. Значение B может быть считано непосредственно прочь этой линии с масштабированием$\sqrt{\frac{2\ln{2}}{\pi}} \approx 0.664$.

% Find the index where the slope of the log-scaled Allan deviation is equal
% to the slope specified.
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));
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

B =

    0.0020

Теперь, когда все шумовые параметры были вычислены, постройте отклонение Аллана со всеми линиями, используемыми для определения количества параметров.

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

Симуляция гироскопа

Используйте imuSensor объект симулировать измерения гироскопа шумовыми параметрами, идентифицированными выше.

% Simulating the gyroscope measurements takes some time. To avoid this, the
% measurements were generated and saved to a MAT-file. By default, this
% example uses the MAT-file. To generate the measurements instead, change
% this logical variable to true.
generateSimulatedData = false;

if generateSimulatedData
    % Set the gyroscope parameters to the noise parameters determined
    % above.
    gyro = gyroparams('NoiseDensity', N, 'RandomWalk', K, ...
        'BiasInstability', B);
    omegaSim = helperAllanVarianceExample(L, Fs, gyro);
else
    load('SimulatedSingleAxisGyroscope', 'omegaSim')
end

Вычислите симулированное отклонение Аллана и сравните его с записанными данными.

[avarSim, tauSim] = allanvar(omegaSim, 'octave', Fs);
adevSim = sqrt(avarSim);
adevSim = mean(adevSim, 2); % Use the mean of the simulations.

figure
loglog(tau, adev, tauSim, adevSim, '--')
title('Allan Deviation of HW and Simulation')
xlabel('\tau');
ylabel('\sigma(\tau)')
legend('HW', 'SIM')
grid on
axis equal

График показывает, что модель гироскопа создала из imuSensor генерирует измерения с подобным отклонением Аллана к записанным данным. Измерения модели содержат немного меньше шума начиная с квантования, и связанные с температурой параметры не устанавливаются с помощью gyroparams. Модель гироскопа может использоваться, чтобы сгенерировать измерения с помощью перемещений, которые легко не получены с оборудованием.

Ссылки

  • Станд. IEEE 647-2006 руководств формата спецификации стандарта IEEE и процедура тестирования для гироскопов лазера Одно Оси

Документация Sensor Fusion and Tracking Toolbox

Поддержка

Для просмотра документации необходимо авторизоваться на сайте