基于物理IMU数据，如何为陀螺仪模拟集成ARW、BI与谱密度？

问题背景 #

我正在开展一个依赖IMU数据的状态估计项目，可获取物理IMU数据，但希望模拟具备特定参数（包括Angle Random Walk (ARW)、Bias Instability (BI)以及功率谱密度）的陀螺仪数据。已有集成ARW和BI参数的伪代码，但在将谱密度集成到模拟中遇到困难，需要修改现有代码加入谱密度参数。

现有代码片段：

# Pseudocode for modeling a noisy gyro sensor
# Created By: MicWro Engr
# 'w' variables are supposed to be the Greek letter omega, used to denote the angular rate
# Remember: variance = (std_dev)^2
# Gyro model:
#   w_meas = w_true + bias + noise
import numpy as np

dt = 0.01 # [sec] sample rate

# Gyro bias parameters
w_bias = -0.01  # This is the constant offset in rate measurements. Also acts as the integrator's initial condition
bias_coef = 1.31e-4  # Bias instability coefficient
bias_var = ((2 * bias_coef**2) / np.pi) * np.log(2)  # bias instability variance

# Angle random walk (AWR) parameters
arw_coef = 4.54e-4  # ARW coefficient
arw_var = arw_coef**2 / dt  # [rad/sqrt(sec)] ARW variance

# simulation loop
while True:
    # obtain an actual rate from ground truth simulation
    w_actual = GroundTruth(...)
    
    # simulate bias instability noise
    w_bias += dt * np.random.normal(0, np.sqrt(bias_var))  # (mean, std_dev)
    
    # simulate angle random walk noise
    w_arw = np.random.normal(0, np.sqrt(arw_var))
    
    # simulated gyro measurement
    w_meas = w_actual + w_bias + w_arw

解决方案 #

核心概念：功率谱密度（PSD）与现有噪声分量的关系 #

陀螺仪的典型功率谱密度可分解为多个核心分量：

• 角度随机走（ARW）：对应PSD中的白噪声分量（1/f⁰），PSD值为 ( S_{ARW} = ARW_{\text{coef}}^2 )（单位：( \text{rad}^2/\text{s} )）
• 零偏不稳定性（BI）：对应PSD中的闪烁噪声分量（1/f¹），PSD值为 ( S_{BI} = BI_{\text{coef}}^2 )（单位：( \text{rad}^2/\text{s}3 )）

现有代码已基于ARW/BI系数构建噪声，我们可以直接将PSD参数作为输入，反向推导对应噪声分量的方差，或扩展代码支持自定义谱密度分量。

修改后的代码：支持PSD参数输入 #

下面的代码直接接受PSD参数，同时保留ARW/BI系数的兼容性，还可扩展其他谱分量：

import numpy as np

dt = 0.01  # [sec] sample rate

# --------------------------
# 功率谱密度（PSD）参数配置
# --------------------------
# ARW对应的白噪声PSD (rad²/s)
psd_arw = (4.54e-4)**2  # 等价于arw_coef²
# BI对应的闪烁噪声PSD (rad²/s³)
psd_bi = (1.31e-4)**2   # 等价于bias_coef²

# 从PSD推导各噪声分量的标准差
# ARW噪声标准差：sqrt(psd_arw / dt)
arw_std = np.sqrt(psd_arw / dt)
# BI噪声（随机游走过程）的驱动噪声标准差：sqrt( (2 * psd_bi / np.pi) * np.log(2) )
bi_drive_std = np.sqrt( (2 * psd_bi / np.pi) * np.log(2) )

# 初始零偏
w_bias = -0.01

# 模拟循环
while True:
    # 获取真实角速度
    w_actual = GroundTruth(...)  # 替换为你的真实值生成逻辑
    
    # 更新零偏不稳定性（一阶高斯马尔可夫过程）
    w_bias += dt * np.random.normal(0, bi_drive_std)
    
    # 生成角度随机走噪声
    w_arw = np.random.normal(0, arw_std)
    
    # 生成最终测量值
    w_meas = w_actual + w_bias + w_arw

扩展：支持更复杂的谱密度分量 #

如果需要模拟完整的陀螺仪谱密度（比如加入速率随机走、量化噪声等），可继续添加对应分量：

• 速率随机走（RRW）：对应PSD中的1/f⁻²分量，需用二阶积分过程模拟
• 量化噪声：对应高频白噪声，可通过固定步长的均匀分布模拟

示例：添加速率随机走（基于PSD参数）

# 速率随机走PSD (rad²/s^5)
psd_rrw = 1e-8
rrw_drive_std = np.sqrt(psd_rrw * dt)  # 驱动噪声标准差
w_rrw = 0  # 速率随机走初始值

# 在循环中更新：
w_rrw += dt * np.random.normal(0, rrw_drive_std)
# 测量值加入该分量
w_meas = w_actual + w_bias + w_arw + w_rrw

关键说明 #

1. 所有噪声分量的模拟都是基于功率谱密度的频域特性转换为时域随机过程：
   • 白噪声（ARW）直接用正态分布采样
   • 闪烁噪声（BI）用一阶高斯马尔可夫过程（随机游走）模拟
   • 更低频的分量（RRW）用高阶积分过程模拟
2. 如果你的谱密度是通过Allan方差分析得到的，需确保PSD参数和Allan方差系数的转换关系正确：
   • ( ARW_{\text{coef}} = \sqrt{S_{ARW}} )
   • ( BI_{\text{coef}} = \sqrt{S_{BI}} )

内容的提问来源于stack exchange，提问作者user3106378