一、准备工作

观测值数据读取是进行数据处理的前提，通常，观测值的数据格式有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双系统，单频数据，结果图如下：

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