imu数据仿真程序地址: https://github.com/Aceinna/gnss-ins-sim. Allan方差分析代码地址： https://github.com/gaowenliang/imu_utils.

如表所示：

import os import math import numpy as np from gnss_ins_sim.sim import imu_model from gnss_ins_sim.sim import ins_sim  # globals D2R = math.pi/180  motion_def_path = os.path.abspath('.//demo_motion_def_files//') fs = 100.0          # IMU sample frequency  def test_allan():     ''' An Allan analysis demo for Sim. '''     #### Customized IMU model parameters, typical for IMU381     ''' 关于imu建立模型： 1.accuracy精度 该系统定义了可直接使用的三种精度 ： 'low-accuracy', 'mid-accuracy' , 'high accuracy'. 或以下自定义精度模型： 'gyro_b'：陀螺仪bias，度/小时 'gyro_arw随机游走陀螺角度，度/ RT-小时 'gyro_b_stability'：陀螺仪bias不稳定，度/小时 'gyro_b_corr陀螺仪偏置不稳定相关性(秒)。 将其设置为‘’inf使用随机游走模型 使用一阶高斯-马尔可夫模型将其设置为正实数 'accel_b加速度计bias，米/秒^ 2 'accel_vrw加速度计速度随机游走，米/秒/ RT-小时 'accel_b_stability加速度计bias不稳定性，米/秒^ 2 'accel_b_corr加速度计偏置不稳定性相关性(秒)。 'mag_std磁噪声STD，UT 2.axis 6:只生成陀螺仪和加速度计数据。 'mag_std磁噪声STD，UT 2.axis 6：仅生成陀螺仪和加速度计数据。 9:除陀螺仪和加速度计数据外，生成磁力计数据。 3.gps True生成GPS数据，False不生成。 '''     imu_err = 
       
        { 
        
        'gyro_b'
        : np
        .array
        (
        [
        0.0
        , 
        0.0
        , 
        0.0
        ]
        )
        , 
        'gyro_arw'
        : np
        .array
        (
        [
        0.25
        , 
        0.25
        , 
        0.25
        ]
        ) 
        * 
        1.0
        , 
        'gyro_b_stability'
        : np
        .array
        (
        [
        3.5
        , 
        3.5
        , 
        3.5
        ]
        ) 
        * 
        1.0
        , 
        'gyro_b_corr'
        : np
        .array
        (
        [
        100.0
        , 
        100.0
        , 
        100.0
        ]
        )
        , 
        'accel_b'
        : np
        .array
        (
        [
        0.0e-3
        , 
        0.0e-3
        , 
        0.0e-3
        ]
        )
        , 
        'accel_vrw'
        : np
        .array
        (
        [
        0.03119
        , 
        0.03009
        , 
        0.04779
        ]
        ) 
        * 
        1.0
        , 
        'accel_b_stability'
        : np
        .array
        (
        [
        4.29e-5
        , 
        5.72e-5
        , 
        8.02e-5
        ]
        ) 
        * 
        1.0
        , 
        'accel_b_corr'
        : np
        .array
        (
        [
        200.0
        , 
        200.0
        , 
        200.0
        ]
        )
        , 
        'mag_std'
        : np
        .array
        (
        [
        0.2
        , 
        0.2
        , 
        0.2
        ]
        ) 
        * 
        1.0 
        } 
        # do not generate GPS and magnetometer data imu 
        = imu_model
        .IMU
        (accuracy
        =imu_err
        , axis
        =
        6
        , gps
        =
        False
        ) 
        #### Allan analysis algorithm 
        from demo_algorithms 
        import allan_analysis algo 
        = allan_analysis
        .Allan
        (
        ) 
        #### start simulation 
        ''' 关于Sim的说明： [fs, 0.0, 0.0]：imu（陀螺仪和加速度计），GPS和磁力计的采样频率，此demo中没有后两者 motion_def_path+"//motion_def-Allan.csv"：初始条件和运动定义，需要指定算法和运动场景。参考ins_sim.py 对于Allan方差分析，需要采集的是IMU静止时的传感器数据，因此，运动场景非常简单：保持静止。 在完整的运动仿真中，运动场景定义需要包含一系列的运动模式或机动指令。 ref_frame：0代表NED坐标系，1代表惯性系 mode：车辆操纵能力，[最大加速度，最大角加速度，最大角速度] env： algorithm：算法对象 ''' sim 
        = ins_sim
        .Sim
        (
        [fs
        , 
        0.0
        , 
        0.0
        ]
        , motion_def_path
        +
        "//motion_def-Allan.csv"
        , ref_frame
        =
        1
        , imu
        =imu
        , mode
        =
        None
        , env
        =
        None
        , algorithm
        =algo
        ) 
        ''' run():参数为执行次数，默认1 results()：保存数据，参数为路径，空则不保存 plot()：画图 ''' sim
        .run
        (
        ) 
        # generate simulation results, summary, and save data to files sim
        .results
        (
        ) 
        # save data files 
        # plot data sim
        .plot
        (
        [
        'ad_accel'
        , 
        'ad_gyro'
        ]
        ) 
        if __name__ 
        == 
        '__main__'
        : test_allan
        (
        ) 
       

挖个坑，还没尝试～～

对于其中的命令类型，官方有如下说明：

对于使用转台的IMU分立级标定而言，至少需要6个不同的转台位置：Z轴朝上，Z轴朝下，X轴朝上，X轴朝下，Y轴朝上，Y轴朝下（定义时多了两个复位）。 根据上述的说明定义转台的运动场景： 在demo_motion_def_files文件夹下定义motion_def_my_imu.csv： 接下来参考demo_allan.py，在gnss-ins-sim文件夹下使用我们刚刚定义的运动模式产生仿真的imu数据，在这里我们要将加速度计三个轴的零偏设置为[0.1,0.2,0.3]，具体实现Sim_IMU.py如下：

import numpy as np
import os
import math
from gnss_ins_sim.sim import imu_model
from gnss_ins_sim.sim import ins_sim

# globals
D2R = math.pi/180

motion_def_path = os.path.abspath('.//demo_motion_def_files//')
fs = 100.0          # IMU sample frequency

def create_sim_imu():     #生成 imu 方针imu数据
    ### Customized IMU model parameters, typical for IMU381
    imu_err = { 
        'gyro_b': np.array([0.0, 0.0, 0.0]),
               'gyro_arw': np.array([0.25, 0.25, 0.25]) * 1.0,
               'gyro_b_stability': np.array([3.5, 3.5, 3.5]) * 1.0,
               'gyro_b_corr': np.array([100.0, 100.0, 100.0]),
               'accel_b': np.array([0.1, 0.2, 0.3]),
               'accel_vrw': np.array([0.03119, 0.03009, 0.04779]) * 1.0,
               'accel_b_stability': np.array([4.29e-5, 5.72e-5, 8.02e-5]) * 1.0,
               'accel_b_corr': np.array([200.0, 200.0, 200.0]),
               'mag_std': np.array([0.2, 0.2, 0.2]) * 1.0
              }
    # do not generate GPS and magnetometer data
    imu = imu_model.IMU(accuracy=imu_err, axis=6, gps=False)

    #### start simulation
    sim = ins_sim.Sim([fs, 0.0, 0.0],
                      motion_def_path+"//motion_def_my_imu.csv",
                      ref_frame=1,
                      imu=imu,
                      mode=None,
                      env=None,
                      algorithm=None)
    sim.run(1)
    # generate simulation results, summary, and save data to files
    sim.results('.//my_sim_imu//')  # save data files

if __name__ == '__main__':
    create_sim_imu()

在进行imu仿真的时候，ins_sim.Sim函数内部会调用acc_gen（）函数，该函数在gnss-ins-sim/gnss_ins_sim/pathgen/pathgen.py中定义，在这里函数里可以修改刻度系数和安装误差，修改函数如下：

ks=np.array([[1.02,0.02,-0.03],[0.03,0.96,0.01],[0.04,0.05,1.05]])
ks=np.transpose(ks)
ref_a2=np.dot(ref_a,ks)
a_mea = ref_a2 + acc_bias + acc_bias_drift + acc_noise + acc_vib

修改后产生的加速度数据的内参矩阵为 [ 1.02 0.03 0.04 0.02 0.96 0.05 − 0.03 0.01 1.05 ] \begin{bmatrix} 1.02&0.03 &0.04 \\ 0.02& 0.96 & 0.05\\ -0.03& 0.01 &1.05 \end{bmatrix} ⎣⎡​1.020.02−0.03​0.030.960.01​0.040.051.05​⎦⎤​ 产生了如下仿真数据：

The following results are saved:
        time: sample time
        ref_pos: true position in the local NED frame
        ref_vel: true vel in the NED frame
        ref_att_euler: true attitude (Euler angles, ZYX)
        ref_accel: true accel in the body frame
        ref_gyro: true angular velocity in the body frame
        accel: accel measurements
        gyro: gyro measurements
        ref_att_quat: true attitude (quaternion)

代码采样频率为100hz，每个位置停留20s，即每个位置会产生2000个数据.

在读取的时候遇到了一些问题，我首先打开了文件查看数据格式（默认wps），看起来非常的规范（没有逗号之类的间隔），可是我单个数据读取的时候，莫名其妙的多了很多的间隔，一个数据有些被分成了两个，完全错乱，所以后来改成行读入，按逗号分隔，读入成功.

#include <vector>
#include <Eigen/Core>
#include <iostream>
#include <fstream>
#include <Eigen/Geometry>

#define D2R 0.017453292519943295 //角度转弧度

bool ReadData(const std::string &path ,std::vector<Eigen::Vector3d ,Eigen::aligned_allocator<Eigen::Vector3d>> &dataVector)
{ 
        
    //清除原有数据
    dataVector.clear();

    std::ifstream fin(path);

    std::string line="";
    std::string temp="";
    int cnt=0;
    while (std::getline(fin,line))
    { 
        
        cnt++;
        if (cnt==1)
        { 
        
            line="";
            temp="";    
            continue;
        }
        std::vector<double> acc;
        for (int i=0;i<line.size();i++)
        { 
        
            if (line[i]==',')
            { 
        
                acc.push_back(std::stod(temp));
                temp="";
            }
            else
            { 
        
                temp+=line[i];
            }         
        }
        acc.push_back(std::stod(temp));
        line="";
        temp="";

        dataVector.push_back(Eigen::Vector3d(acc[0],acc[1],acc[2]));
        std::cout<<cnt<<": "<<acc[0]<<", "<<acc[1]<<", "<<acc[2]<<std::endl;
    }
    if (cnt==0)
    { 
        
        std::cerr<<"read fail"<<std::endl;
        return false;
    }
    return true;
}

在读取数据的时候我注意到其实数据总共只有15725个（仿真出来的数据也就这么多，不是没有读取到），但是按照推算，8个场景（算上两次复位动作），100hz，20s。那么应该产生16000个数据，显然这存在一些问题.但是我们的标定最少只需要6个不同位置的数据，所以，我们要对这些数据进行一些优化,然后进行标定，代码如下：

void Calibration(std::vector<Eigen::Vector3d ,Eigen::aligned_allocator<Eigen::Vector3d>> &accs)
{ 
        
    Eigen::Vector3d truth_z_down=CalculateMean(accs,1000,1500);
    Eigen::Vector3d truth_z_up=CalculateMean(accs,3000,3500);
    Eigen::Vector3d truth_y_down=CalculateMean(accs,7000,7500);
    Eigen::Vector3d truth_y_up=CalculateMean(accs,9000,9500);
    Eigen::Vector3d truth_x_up=CalculateMean(accs,13000,13500);
    Eigen::Vector3d truth_x_down=CalculateMean(accs,15000,15500);
/* std::cout<<"truth_z_up :"<<std::endl; std::cout<<truth_z_up<<std::endl; std::cout<<"truth_z_down :"<<std::endl; std::cout<<truth_z_down<<std::endl; std::cout<<"truth_y_up :"<<std::endl; std::cout<<truth_y_up<<std::endl; std::cout<<"truth_y_down :"<<std::endl; std::cout<<truth_y_down<<std::endl; std::cout<<"truth_x_up :"<<std::endl; std::cout<<truth_x_up<<std::endl; std::cout<<"truth_x_down :"<<std::endl; std::cout<<truth_x_down<<std::endl; */
    Eigen::Vector3d calib_bias;
    Eigen::Matrix3d calib_sacle_factor;
    Eigen::Matrix3d calib_installation_error;
    Eigen::Matrix3d calib_internal_matrix=Eigen::Matrix3d::Identity();

    //解析法求解 利用6组数据建立六个方程
    calib_bias=(truth_z_up+truth_z_down)/2;
    Eigen::Vector3d calib_internal_matrix_col3=(truth_z_up-calib_bias)/g;
    Eigen::Vector3d calib_internal_matrix_col2=(truth_y_up-calib_bias)/g;
    Eigen::Vector3d calib_internal_matrix_col1=(truth_x_up-calib_bias)/g;

    calib_internal_matrix<<calib_internal_matrix_col1[0],calib_internal_matrix_col1[1],calib_internal_matrix_col1[2],calib_internal_matrix_col2[0],calib_internal_matrix_col2[1],calib_internal_matrix_col2[2],calib_internal_matrix_col3[0],calib_internal_matrix_col3[1],calib_internal_matrix_col3[2];
    calib_internal_matrix.transpose();
    std::cout<<"calib_internal_matrix: "<<std::endl;
    std::cout<<calib_internal_matrix<<std::endl;
    std::cout<<"calib_bias: "<<std::endl;
    std::cout<<calib_bias<<std::endl;
}

int main(int argc ,char** argv)
{ 
        
    //加速度计仿真数据
    std::string accel_csv_path="./gnss-ins-sim/my_sim_imu/accel-0.csv";
    
    std::vector<Eigen::Vector3d ,Eigen::aligned_allocator<Eigen::Vector3d>> accs;
    ReadData(accel_csv_path,accs);
    Calibration(accs);
}

设定的真实内参矩阵： [ 1.02 0.03 0.04 0.02 0.96 0.05 − 0.03 0.01 1.05 ] \begin{bmatrix} 1.02&0.03 &0.04 \\ 0.02& 0.96 & 0.05\\ -0.03& 0.01 &1.05 \end{bmatrix} ⎣⎡​1.020.02−0.03​0.030.960.01​