利用Allan方差分析惯性传感器噪声

打开脚本

本例展示了如何使用Allan方差来确定MEMS陀螺仪的噪声参数。这些参数可用于陀螺仪的仿真建模。陀螺仪测量模型为:

$$ \ω(t) = \ Omega_{理想}(t) + Bias_N (t) + Bias_B (t) + Bias_K (t) $$

三个噪声参数N(角度随机游走)，K(概率随机游走)，和B(偏置不稳定性)估计使用来自固定陀螺仪的记录数据。

背景

Allan方差最初由David W. Allan开发，用于测量精密振荡器的频率稳定性。它也可以用来识别各种噪声源存在于固定陀螺仪测量。考虑l来自陀螺仪的采样数据具有采样时间 $美元\ tau_ {0} $$ ．形成持续时间的数据集群 $美元\ tau_ {0} $$ ， $2美元\ tau_ {0} $$ 、…… $m\tau_{0}， (m < (L-1)/2)$ 并获得每个聚类中包含的数据点的和在聚类长度上的平均值。Allan方差被定义为数据簇平均的两样本方差作为簇时间的函数。本例使用重叠Allan方差估计器。这意味着计算的聚类是重叠的。的较大值时，估计器的性能优于非重叠估计器l．

艾伦方差计算

Allan方差计算如下:

日志l固定式陀螺仪具有采样周期 $美元\ tau_ {0} $$ ．让 $ω\美元$ 记录样本。

从三轴陀螺仪的一个轴加载记录数据。这段录音以100 Hz采样率在6小时内完成。负载(“LoggedSingleAxisGyroscope”，“ω”，“Fs”) t0 = 1/Fs;

对于每个样本，计算输出角度 $\θ美元$ ：

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

对于离散样本，累积和乘以 $美元\ tau_ {0} $$ 可以使用。

= cumsum(， 1)*t0;

接下来，计算Allan方差:

$$ \σ^ 2(\τ)= & # xA; \压裂{1}{2 \τ^ 2}< (\ theta_ {k + 2 m} 2 \ theta_ {k + m} + \ theta_ {k}) ^ 2 & # 62; $$

在哪里 $\tau = m\tau_{0}$ 而且 < & # 62;美元是整体平均值。

集合平均值可以扩展为:

$$ \σ^ 2(\τ)= & # xA; \压裂{1}{2 \τ^ 2 (L-2m)} \ sum_ {k = 1} ^ {L-2m} (\ theta_ {k + 2 m} - 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 = cell (M);% m必须为整数。M =唯一的(M);删除重复项。Tau = m*t0;Avar = 0(数字(m)， 1);为I = 1: number (m) mi = m(I);Avar (i，:) = sum(.．.(θ(1 + 2 * mi: L) - 2θ* (1 + mi: L-mi) +θ(1:l2 * mi))。^ 2, 1);结束Avar = Avar ./ (2*tau.;)^2 .* (L - 2*m));

最后，艾伦偏差 $\sigma(t) = \根号{\sigma²(t)}$ 用于确定陀螺仪的噪声参数。

adv =根号(avar);图logog (tau, adev) title(“艾伦偏差”)包含(“\τ”）;ylabel (“\σ(\τ)”网格)在轴平等的

艾伦方差也可以用allanvar函数。

[avarFromFunc, tauFromFunc] = allanvar(ω， m, Fs);adevFromFunc = sqrt(avarFromFunc);figure loglog(tau, adev, tauFromFunc, adevFromFunc);标题(“艾伦偏差”)包含(“\τ”) ylabel (“\σ(\τ)”)传说(手工计算的，“allanvar函数”网格)在轴平等的

噪声参数识别

为了获得陀螺仪的噪声参数，可以使用原始数据集中噪声参数的Allan方差与双侧功率谱密度(PSD)之间的关系 $ω\美元$ ．关系是:

$$ \σ^ 2(\τ)= 4 \ int_ {0} ^ {\ infty} S_ \ω(f) & # xA; \压裂{罪^ 4 f(\π\τ)}{(f \ \πτ)^ 2}df $$

由上式可知，当陀螺仪通过传递函数为的滤波器时，Allan方差与陀螺仪的总噪声功率成正比罪^ 4美元(x) / (x) ^ 2美元．这个传递函数产生于创建和操作集群所做的操作。

使用此传递函数解释，滤波器的带通取决于 $\τ美元$ ．这意味着不同的噪声参数可以通过改变滤波器带通或变化来识别 $\τ美元$ ．

角度随机游走

利用陀螺仪输出的白噪声谱来表征角度随机游走。PSD表示为:

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

在哪里

N=角度随机游走系数

代入原PSD方程进行积分得到:

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

上面的方程是一条斜率为-1/2的直线 $\σ(\τ)美元$ 与 $\τ美元$ ．的价值N可以直接从这一行读到 $\tau = 1$ ．的单位N是 $(rad / s) / \√{赫兹}$ ．

求log- scaling Allan deviation的斜率相等的指数%到指定的斜率。斜率= -0.5;Logtau = log10(tau);Logadev = log10(adev);Dlogadev = diff(logadev) ./ diff(logtau);[~， i] = min(abs(dlogadev - slope));求直线的y轴截距。B = logadev(i) - slope*logtau(i);从直线上确定角度随机游走系数。logN = slope*log(1) + b;N = 10^logN绘制结果。tauN = 1;亚麻= N ./平方根(tau);figure loglog(tau, adev, tau, lineN，“——”，褐色，N，“o”)标题(“带角度随机游走的艾伦偏差”)包含(“\τ”) ylabel (“\σ(\τ)”)传说(‘\σ’，“\ sigma_N”)文本(牛角，N，“N”网格)在轴平等的

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的直线 $\σ(\τ)美元$ 与 $\τ美元$ ．的价值K可以直接从这一行读到 $\tau = 3$ ．的单位K是 $(rad / s \√{赫兹}$ ．

求log- scaling Allan deviation的斜率相等的指数%到指定的斜率。斜率= 0.5;Logtau = log10(tau);Logadev = log10(adev);Dlogadev = diff(logadev) ./ diff(logtau);[~， i] = min(abs(dlogadev - slope));求直线的y轴截距。B = logadev(i) - slope*logtau(i);从该线确定随机游走系数。logK = slope*log10(3) + b;K = 10^logK绘制结果。tauK = 3;lineK = K .*平方根(tau/3);figure loglog(tau, adev, tau, lineK，“——”， tauK, K，“o”)标题(“速率随机游走的艾伦偏差”)包含(“\τ”) ylabel (“\σ(\τ)”)传说(‘\σ’，“\ sigma_K”text(tauK, K，“K”网格)在轴平等的

K = 9.0679 -05

偏见不稳定

陀螺仪输出的粉色噪声(闪烁噪声)谱表征了偏置不稳定性。PSD表示为:

$$ S_{\ω}左(f) = \ \{\{数组}{lr}开始(\压裂{B ^ 2}{2 \π})\压裂{1}{f} & # xA; & # 38;: f \leq f_0\\0 &: f >f_0 \结束数组{}\ $$

在哪里

B偏差不稳定系数

f_0美元截止频率

代入原PSD方程进行积分得到:

$$ \σ^ 2(\τ)= \压裂{2 b ^ 2}{\π}[\ ln {2} + \ \ & # xA; - \压裂{罪^ 3 x} {2 x ^ 2} (sinx + 4 xcosx) + Ci (2 x) - Ci (4 x)) $$

在哪里

$x = \pi f_0\tau$

Ci=余弦积分函数

当 $\τ美元$ 远长于截止频率的倒数，则PSD方程为:

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

上面的方程是一条斜率为0的直线 $\σ(\τ)美元$ 与 $\τ美元$ ．的价值B可以直接从这一行读取 $\sqrt{\frac{2\ln{2}}{\pi}} \约0.664$ ．的单位B是 rad / s美元．

求log- scaling Allan deviation的斜率相等的指数%到指定的斜率。斜率= 0;Logtau = log10(tau);Logadev = log10(adev);Dlogadev = diff(logadev) ./ diff(logtau);[~， i] = min(abs(dlogadev - slope));求直线的y轴截距。B = logadev(i) - slope*logtau(i);从直线上确定偏置不稳定系数。scfB =√(2*log(2)/pi);logB = b - log10(scfB);B = 10^logB绘制结果。tauB = tau(i);lineB = B * scfB * ones(size(tau));figure loglog(tau, adev, tau, lineB，“——”， tauB, scfB*B，“o”)标题(“带偏差不稳定性的艾伦偏差”)包含(“\τ”) ylabel (“\σ(\τ)”)传说(‘\σ’，“\ sigma_B”(tauB, scfB*B，“0.664”网格)在轴平等的

B = 0.0020

现在已经计算了所有的噪声参数，用所有用于量化参数的线绘制艾伦偏差。

tauParams = [tauN, tauK, tauB];params = [N, K, scfB*B];figure loglog(tau, adev, tau，[亚麻，lineK, lineB]，“——”，.．.tauParams参数,“o”)标题(“带噪声参数的艾伦偏差”)包含(“\τ”) ylabel (“\σ(\τ)”)传说(“\σ(rad / s)的美元，“美元\ sigma_N ((rad / s) / \√{赫兹})的美元，.．.“美元\ sigma_K (rad / s \√{赫兹})的美元，“美元\ sigma_B (rad / s)的美元，“翻译”，“乳胶”) text(tauParams, params， {“N”，“K”，“0.664”})网格在轴平等的

陀螺仪模拟

使用imuSensor对象用上面确定的噪声参数模拟陀螺仪测量。

模拟陀螺仪测量需要一些时间。为了避免这种情况，%测量值生成并保存到mat文件。默认情况下，% example使用mat文件。要生成度量值，请更改%此逻辑变量为true。generateSimulatedData = false;如果generateSimulatedData将陀螺仪参数设置为确定的噪声参数%以上。陀螺=陀螺参数(“NoiseDensity”N“随机散步”、钾、.．.“BiasInstability”, B);omegaSim = helperAllanVarianceExample(L, Fs，陀螺);其他的负载(“SimulatedSingleAxisGyroscope”，“omegaSim”）结束

计算模拟阿伦偏差，并将其与测井数据进行比较。

[avarSim, tauSim] = allanvar(omegaSim，“八度”Fs);adevSim = sqrt(avarSim);adevSim = mean(adevSim, 2);使用模拟的平均值。figure loglog(tau, adev, tauSim, adevSim，“——”)标题(“HW和模拟的艾伦偏差”)包含(“\τ”）;ylabel (“\σ(\τ)”)传说(“HW”，“SIM卡”网格)在轴平等的

图中显示的陀螺仪模型创建自imuSensor生成与记录数据具有相似艾伦偏差的测量值。由于量化和温度相关参数不使用，模型测量的噪声略小gyroparams．陀螺仪模型可用于使用不容易被硬件捕获的运动来生成测量值。

参考文献

单轴激光陀螺的标准规范格式指南和测试程序

利用Allan方差分析惯性传感器噪声

背景

艾伦方差计算

噪声参数识别

角度随机游走

随机游走率

偏见不稳定

陀螺仪模拟

参考文献

传感器融合和跟踪工具箱文档

金宝app

自主系统的传感器融合与跟踪