Matlab读取GNSS 观测值o文件代码示例

一、准备工作
观测值数据读取是进行数据处理的前提，通常，观测值的数据格式有rtcm、ubx、rinex，各家厂商还有自定义的格式。数据读取是最简单的工作，掌握了数据组织格式即掌握了数据读取策略，简言之，我们是对字符串进行一系列操作。我认为对o文件读取要做以下准备工作，分别是：
（1）安装Matlab的PC;
（2) 学习rinex.pdf资料；
（3）o文件准备；
（4）伪代码编写，理清数据流；
（5）代码测试
我们在使用定位方法SPP/PPP/RTK时，是按照历元来读取卫星的观测值（PR/CP），Doppler，SNR等信息，后续的处理是对所有卫星数进行循环遍历，进行不同的操作。
二、需解决的问题
问题是：要实现使用Matlab绘制不同卫星系统的卫星在观测弧段的伪距观测值时序图。
在上述问题分析中，易知卫星PRN号是唯一的，可作为key；而卫星对应着观测值（PC/CP/Doppler），可作为value。换言之，key和value可以使用map的映射来作为数据结构，其中，value可以使用结构体保存。最重要的一环已经确定了，即使用map和struct。
接下来，我们拿一颗卫星的观测弧段内的伪距值来说，可能的情况有以下几种：
（1） 在观测弧段内，卫星一直被捕获，即100%存在伪距值；
（2） 在观测弧段内，卫星在开始历元没有被捕获，而后才被捕获；
（3）在（2）中还存在未捕获，捕获，未捕获，捕获等情况；
在未捕获的历元时观测值为0，同时，我们还要注意，有伪距不一定有载波，没有载波肯定没有伪距。我们的伪代码主要是对上述三种情况进行判断即可。
三、代码实现
具体问题具体分析，解决上述提到的问题，具体代码如下：

% function: read the body of rinex obs file
clear all
close all
clc

fid=fopen(' example.obs');
if fid==-1
    disp('fail to open obs file!');
end

% 创建结构体obs，使用Map进行映射，prn为键(char)，sv为值(table)
obs=containers.Map;
record_counter=0;

% 读取body记录
while feof(fid)==0
    line=fgets(fid);
    if strcmp(line(1),'>')
        sv_sum=str2double(line(34:35));
        time=line(3:29);
        if sv_sum>0
            for i=1:sv_sum
                line=fgets(fid);
                % 键 -> prn
                key=line(1:3);
                % 值 -> sv(pr cp dop)
                sv.pr=str2double(line(6:17));
                cp=line(21:33);
                if isnan(cp)
                    sv.cp=0;
                else
                    sv.cp=str2double(line(21:33));
                end
                sv.dop=str2double(line(41:49));
                % 值逻辑判断
                if isKey(obs,key)
                    jump_sec=record_counter-length(obs(key));
                    if jump_sec~=0 % 解决掉星后重新锁定的问题
                        temp.pr=0;
                        temp.cp=0;
                        temp.dop=0;
                        tmp_strcut=repmat(temp,1,jump_sec);
                        value=[obs(key) tmp_strcut sv];
                    else % 卫星一直锁定，直接追加
                        if strcmp(key,'C18')==1
                            if length(obs(key))>5713
                                disp(time);
                            end
                        end
                        value=[repmat(obs(key),1) repmat(sv,1)]; % mark
                    end
                else % 卫星首次锁定
                    if record_counter~=0 % 第一秒未锁定
                        temp.pr=0;
                        temp.cp=0;
                        temp.dop=0;
                        tmp_strcut=repmat(temp,1,record_counter);
                        value=[tmp_strcut sv];
                    else % 第一秒锁定
                        value=sv;
                    end
                end
                obs(key)=value;
            end
        end
        record_counter=record_counter+1;
    end
end
disp('create struct obs(body) successfully!');

% 绘图
sv_sum=obs.Count; % 可视卫星数量
keySet=keys(obs); % 可视卫星名称
disp(keySet);

v=299792458.0;
bf1=1.561098E9; % bds L1 频率
gf1=1.57542E9; % gps L1 频率
gps_sum=0;bds_sum=0;

for i=1:sv_sum
    prn=char(keySet(i));
    if strcmp(prn(1),'G')==1 % GPS
        gps_sum=gps_sum+1;
    else if strcmp(prn(1),'C')==1 % BDS
            bds_sum=bds_sum+1;
        end
    end
end
subfig1=ceil(sqrt(bds_sum));
subfig2=ceil(sqrt(gps_sum));

x=1:record_counter;
for i=1:sv_sum
    prn=char(keySet(i));
    info=obs(prn);
    pr=zeros(1,record_counter);cp=pr;
    for j=1:length(info)
        pr(j)=info(j).pr;
        cp(j)=info(j).cp;
    end
    if i==1 || i==bds_sum+1
        figure;
    end
    subplot(n,n,pos);
    plot(x,pr,'r',x,cp*f,'g'); 
    title(char(prn));
end

四、示例结果
我只针对工作中遇到的rinex版本数据进行处理，gps/bds双系统，单频数据，结果图如下：

示例代码比较简单，也有相应的注释。