1.使用数据:
M40采集的imu观测值、gps观测值
2.数据处理程序
NovAtelConvert、GINS、RTKLIB、ins_gnss_losse(matlab版本)
默认这四个程序已经安装到电脑上。
3. 数据预处理
1.数据解码、获取导航文件
M40采集的数据文件如下:
1.1解码
- step1.解码.gnss文件:选中文件,右键convertto,选择rinex3.04
- step2,打开gins,选择文件,解码原始文件,导入数据,设置高级里边的IMU轴向调整如下图,解码即。
因为安装m40时,前向并未对准行进方向因此进行调整,将坐标系转为前右下。
注:使用gins也可以解码gnss数据但是,观测值卫星数较少,缺少鬼子和阿三的卫星。
- step3.下载当天的导航文件
(1)计算年纪日,计算链接:http://www.igs.gnsswhu.cn/index.php/home/data_product/igs.html
(2)ftp下载
https://cddis.nasa.gov/archive/gnss/data
文件路径说明 https://cddis.nasa.gov/archive/gnss/00readme
/daily/YYYY/DDD/YYd Observation files for year YYYY and day DDD
(Hatanaka format)
/YYf Beidou navigation files for year YYYY and day DDD
/YYg GLONASS navigation files for year YYYY and day
DDD
/YYh SBAS navigation files for year YYYY and day DDD
/YYi IRNSS navigation files for year YYYY and day DDD
/YYl Galileo navigation files for year YYYY and day
DDD
/YYm Met files for year YYYY and day DDD
/YYn GPS navigation files for year YYYY and day DDD
/YYo Observation files for year YYYY and day DDD
多系统星历文件 /YYp Mixed navigation files for year YYYY and day DDD
/YYq QZSS navigation files for year YYYY and day DDD
/YYs Summary files for year YYYY and day DDD
本例年纪日为选111故进入路径:https://cddis.nasa.gov/archive/gnss/data/daily/2022/111/22p/
下载一个大点的文件,然后解压。
注意rtklib暂时不支持rinix4.0版本,下载广播星历文件后需要查看版本号。
- step4.保留导航文件、观测文件、解码后的imu.bin文件
2.使用RTKLIB处理gps原始观测
操作步骤
(1)打开rtkpose.exe 可执行文件
(2)设置输出历元间隔、观测值文件路径、导航文件路径、输出文件路径、选用的卫星系统
(3)option中设置
定位模式选择单点定位、高度角阈值设置为15
输出设置,结果用经纬度高程表示,设置输出速度。
解算:
(4)点击plot查看轨迹和可见卫星数,可以初步查看定位效果。
(5)结果文件,pos
GPST latitude(deg) longitude(deg) height(m) Q ns sdn(m) sde(m) sdu(m) sdne(m) sdeu(m) sdun(m) age(s) ratio vn(m/s) ve(m/s) vu(m/s) sdvn sdve sdvu sdvne sdveu sdvun
2022/04/21 06:47:34.000 39.988723818 116.346130987 53.1772 5 9 8.5511 6.6179 21.2390 -3.5887 -9.0447 7.8971 0.00 0.0 -0.10206 0.10557 -0.21679 0.24766 0.19705 0.61898 -0.11752 -0.26352 0.26150
4. ins-gps-losse松组合
4.1 所需数据:
- rtklib解算的定位结果文件
- imu观测文件
4.2数据预处理
rtklib输出保留经纬度高程和北东地的速度信息,时间改成gps周内秒。
//matlab
%将rtklib解算单点定位结果时间转换成gps周内秒
%保留时间、纬度、精度、高程、enu速度
%%
clear;
clc;
data=importdata('../resource/220421071334.pos');
size=length(data)-14;%rtklib前14行为解释说明
enddata=zeros(size,7);%处理后的数据
%%提取有用的数据
if size> 0 %文件不为空
for j = 1:size
dataj = regexp(data{j+14,1},'\d*\.?\d*','match');%提取单个历元待处理数据
dataj=str2double(dataj);%数据格式转换
enddata(j,2:4)=dataj(7:9);%纬度、经度、高程
enddata(j,5:7)=dataj(20:22);%neu速度
utc=[dataj(1) dataj(2) dataj(3) dataj(4) dataj(5) dataj(6)];
gpst=calculationgpstime(utc)+18;
enddata(j,1)=gpst(2);%gps周内秒
end
end
%保存
Save(enddata,"../resource/220421071334.txt");
utc时间转换gps时间函数:calculationgpstime()
function t = calculationgpstime(t1)
%Calculating GPS time (gpsweek,seconds);
% t1 is the observation time
% TYPE : [YYYY MM DD HH MIN SEC]
% function to return gpsweek and the seconds in week
t0 = [1980 1 6 0 0 0.0];
sec = etime(t1,t0);
t = [fix(sec/3600/24/7),sec - fix(sec/3600/24/7)*3600*24*7];
4.3 ins_gps解算
更改文件输入路径:
imu_inpath = "../resource/220421071334_imu.bin"; % IMU文件路径, 严格以.txt/.dat/.bin结尾,推荐使用二进制, 文件数据要求见ImuFile
gnss_inpath = "../resource/220421071334.txt"; % GNSS文件路径,文件严格为文本文件
其他参数已经标定好了。
%%
%M40 data 主函数
%2022.3.17
%%
clear
clear global var
%% 全局变量
tic;
run arg
t = toc;
fprintf('设置全局变量用时:%f\n', t);
%% 用户设置
global GRAVITY RE %重力和地球半径
imu_inpath = "../resource/220421071334_imu.bin"; % IMU文件路径, 严格以.txt/.dat/.bin结尾,推荐使用二进制, 文件数据要求见ImuFile
gnss_inpath = "../resource/220421071334.txt"; % GNSS文件路径,文件严格为文本文件
pos = [ 39.991478610; 116.344457418; 56.1510]; % 初始位置; % 初始位置
vel = [0; 0; 0];% 初始速度
att = [0;0;0];% 初始姿态[roll, pitch, head]
fs = 200; % IMU采样频率
gyrobias = [0.01;0.01; 0.01]; % 陀螺仪零偏°/h
accebias = [0.05 * GRAVITY * 0.001; 0.05 * GRAVITY * 0.001; 0.05 * GRAVITY * 0.001];% 加速度计零偏 m/s
ARW = [0.002492; 0.001993; 0.002423]; % 陀螺仪白噪声(随机游走结果) °/sqrt(h)
VRW = [0.049549; 0.083993; 0.051888]; % 加速度计白噪声(随机游走结果) m /sqrt(s)
obnosie = [4 / RE; 4 / RE; 1; 0.1; 0.1; 0.1]; % GNSS观测标准差 [°; °; °; m/s; m/s; m/s]
mod = [1, 1, 1, 1, 1, 1]; % 马尔科夫模型相关时间 s
%%
tic;%启动秒表计时器
err = ImuErr(gyrobias, accebias, ARW, VRW);%读取imu误差参数
t = toc;%读取秒表计时器时间
fprintf('ImuErr用时:%f\n', t);
tic;%启动秒表计时器
imudata = ImuFile(imu_inpath);%读取imu数据
t = toc;%读取秒表计时器时间
fprintf('读取IMU文件用时:%f\n', t);
tic;%启动秒表计时器
gnssdata = Gnss(gnss_inpath);%读取gnss数据
t = toc;%读取秒表计时器时间
fprintf('读取GNSS文件用时:%f\n', t);
tic;%启动秒表计时器
inipva = Pva(pos, vel, att);%设置初始pva
t = toc;%读取秒表计时器时间
fprintf('Pva:%f\n', t);
tic;%启动秒表计时器
inertial = Ins(inipva, fs, err);%ins更行(机械编排)
t = toc;%读取秒表计时器时间
fprintf('Ins用时:%f\n', t);
%% 卡尔曼滤波
tic;%启动秒表计时器
kf_debug = Kf(err, obnosie, mod);
t = toc;%读取时间
fprintf('Kf用时:%f\n', t);
tic;%启动秒表计时器
[kf_debug, inertial, imudata] = kf_debug.Update(inertial, gnssdata, imudata);
t = toc;%读取时间
fprintf('卡尔曼滤波更新用时:%f\n', t);
%% 保存结果
tic;
imudata.Save("../resource/m40.dat");
t = toc;
fprintf('数据保存用时:%f\n', t);