一、前言
2024年暑假，武汉大学导航工程专业学生来到了湖北孝感天紫湖度假村进行了专业实习。实习内容主要包括GNSS静态测量和RTK测量。
外业数据采集完之后，需要对静态测量的数据进行数据质量分析，我们可以使用现有的软件，例如RTKlib和华测CGO。老师要求我们自编程序实现GNSS的数据质量分析，我选择使用Matlab的APP designer功能设计数据分析程序，名为WHUGNSS。
至于选择使用matlab的原因，主要有两个，第一是使用matlab处理数据非常方便，第二是matlab具有强大的数据可视化功能，画图十分方便，使用C++，C#来实现会比较困难。缺陷在于无法生成独立的可执行程序。
二、自编程序实现GNSS数据处理
1.Matlab App designer介绍
这是matlab提供的一个可以实现设计图形用户界面的功能。具有一些基本的控件，例如按钮，文本框，坐标框，条形框等等。同学们可以查阅官方文档学习，比较简单，易于上手（我就是现学的）。设计界面如下：
对不同的控件可以添加不同的回调函数，实现不同的功能。例如点击不同控件时，执行什么操作，勾选不同复选框时，执行什么操作。
2.rinex观测值文件读取
要数据处理，就必须要会读文件，该程序支持读取多种格式的文件，包括以O文件，txt文件等，我们外业采集的观测值文件是以24O为后缀的，文件格式如下：
头文件
数据块
文件较大，可以使用matlab直接打开查看，先读取头文件部分，主要包括接收机和基站的一些基本信息，然后读取数据块文件，将不同卫星系统的数据存储到不同的结构体里面，方便后续计算。我是将历元转换成了GPS时，方便存储数据。
文件读取代码如下：
clear;
filename1="D038184K.24O";
fop=fopen(filename1,'r');
linenumber=1;
base_station='0';
base_position=zeros(1,3);
base_h=0;
time_s=0;%采样间隔
i=1;n1=0;
k1=1;k2=1;k3=1;k4=1;k5=1;%用于计数
sat_GPS=struct('PRN','','Time',[],'Data',[],'GPST','');
sat_BDS=struct('PRN','','Time',[],'Data',[],'GPST','');
sat_GLONASS=struct('PRN','','Time',[],'Data',[],'GPST','');
sat_Galileo=struct('PRN','','Time',[],'Data',[],'GPST','');
sat_QZSS=struct('PRN','','Time',[],'Data',[],'GPST','');
GPS_Fre={};
GLONASS_Fre={};
Galileo_Fre={};
QZSS_Fre={};
BDS_Fre={};
while ~feof(fop)
tline=fgetl(fop);
if linenumber==3
base_station=sscanf(tline(1:10),'%s');
end%读取基站名称
if linenumber==8
base_position=sscanf(tline(1:end),"%f");
end%读取基站坐标
if linenumber==9
base_h=sscanf(tline(1:14),"%f");
end%读取天线相位中心高
if linenumber==11
n1=0;
for n=1:13
GPS_Fre{n}=tline(8+n1:8+n1+2);
n1=n1+4;
end
end
if linenumber==12
n1=0;
for n=14:24
GPS_Fre{n}=tline(8+n1:8+n1+2);
n1=n1+4;
end
end
if linenumber==13
n1=0;
for n=1:sscanf(tline(5:6),"%f")
GLONASS_Fre{n}=tline(8+n1:8+n1+2);
n1=n1+4;
end
end
if linenumber==14
n1=0;
for n=1:sscanf(tline(5:6),"%f")
Galileo_Fre{n}=tline(8+n1:8+n1+2);
n1=n1+4;
end
end
if linenumber==15
n1=0;
for n=1:sscanf(tline(5:6),"%f")
QZSS_Fre{n}=tline(8+n1:8+n1+2);
n1=n1+4;
end
end
if linenumber==16
n1=0;
for n=1:13
BDS_Fre{n}=tline(8+n1:8+n1+2);
n1=n1+4;
end
end
if linenumber==17
n1=0;
for n=14:26
BDS_Fre{n}=tline(8+n1:8+n1+2);
n1=n1+4;
end
end
if linenumber==18
n1=0;
for n=27:32
BDS_Fre{n}=tline(8+n1:8+n1+2);
n1=n1+4;
end
end
if linenumber==19
time_s=sscanf(tline(1:end),'%f');
end%读取采样间隔
if(strcmp(tline(1),'>'))
obshead(i,:)=sscanf(tline(2:end),'%d %d %d %d %d %f %d %d');
for k=1:obshead(i,8)
tline=fgetl(fop);
if(strcmp(tline(1),'G'))
sat_GPS(k1).PRN=tline(1:3);
sat_GPS(k1).Time=obshead(i,1:6);
sat_GPS(k1).Data=sscanf(tline(4:end),'%f');
sat_GPS(k1).GPST=(sat_GPS(k1).Time(4)+8)*3600+sat_GPS(k1).Time(5)*60+sat_GPS(k1).Time(6);
k1=k1+1;
end
if(strcmp(tline(1),'C'))
sat_BDS(k2).PRN=tline(1:3);
sat_BDS(k2).Time=obshead(i,1:6);
sat_BDS(k2).Data=sscanf(tline(4:end),'%f');
sat_BDS(k2).GPST=(sat_BDS(k2).Time(4)+8)*3600+sat_BDS(k2).Time(5)*60+sat_BDS(k2).Time(6);
k2=k2+1;
end
if(strcmp(tline(1),'E'))
sat_Galileo(k3).PRN=tline(1:3);
sat_Galileo(k3).Time=obshead(i,1:6);
sat_Galileo(k3).Data=sscanf(tline(4:end),'%f');
sat_Galileo(k3).GPST=(sat_Galileo(k3).Time(4)+8)*3600+sat_Galileo(k3).Time(5)*60+sat_Galileo(k3).Time(6);
k3=k3+1;
end
if(strcmp(tline(1),'R'))
sat_GLONASS(k4).PRN=tline(1:3);
sat_GLONASS(k4).Time=obshead(i,1:6);
sat_GLONASS(k4).Data=sscanf(tline(4:end),'%f');
sat_GLONASS(k4).GPST=(sat_GLONASS(k4).Time(4)+8)*3600+sat_GLONASS(k4).Time(5)*60+sat_GLONASS(k4).Time(6);
k4=k4+1;
end
if(strcmp(tline(1),'J'))
sat_QZSS(k5).PRN=tline(1:3);
sat_QZSS(k5).Time=obshead(i,1:6);
sat_QZSS(k5).Data=sscanf(tline(4:end),'%f');
sat_QZSS(k5).GPST=(sat_QZSS(k5).Time(4)+8)*3600+sat_QZSS(k5).Time(5)*60+sat_QZSS(k5).Time(6);
k5=k5+1;
end
end
i=i+1;
end
linenumber=linenumber+1;
end操作：点击文件，选择打开rinex文件，然后选中指定格式的文件，等待系统读取文件。（这个过程会有一点慢，因为我把一些计算的内容放在了这里，以便提高后续程序的执行效率。
操作图：
读取文件之后程序主界面：
左侧的一些文本框会显示基站的基本信息和采样时间间隔，右侧文本框显示原始数据。
点击不同的卫星系统，程序会筛选出不同卫星系统的原始数据。
如图：
3.信噪比（载噪比）分析
——卫星观测载噪比指标统计值；
n——观测卫星总数；
j——观测卫星号；
——卫星 j 的观测历元总数；
i——卫星 j 的观测历元序号；
——卫星 j 在历元ti时的载噪比观测值，单位dBHz。
利用以上公式计算信噪比的统计值。勾选不同的卫星系统，点击开始计算，下拉框内会显示该系统的可视卫星，同时会在文本框内显示每颗卫星每个历元每个频率下的信噪比。画图时，选取每颗卫星的第一个有观测值的频率，绘出其信噪比随时间的变化图。
选择整个系统
选择单颗卫星
4.电离层延迟
当电离层延迟变化率大于0.07m/s时，表示电离层延迟发生跳变。对所有卫星某个频点的电离层延迟跳变数求和，并作为各系统频点的电离层延迟跳变次数。
在电离层延迟的主界面里，又两个坐标区，第一个绘制电离层延迟时序，第二个绘制电离层延迟变化率时序。同时，在最下面的文本框会显示，频率数，和电离层跳变数，频率数默认为2 。
5.周跳探测
周跳探测采用如下算法：
界面如下：
图中，绘制了卫星的MW值和MW的均方差图。并记录了周跳数，标注在了图上，同时显示在了文本框里，并计算了周跳比，及周跳个数除以观测值个数。
6.多路径计算：
计算伪距多路径误差，必须依赖双频观测数据。采用伪距观测方程和载波相位观测方程组合， 消除对流层和电离层延迟影响，按公式(16)进行计算。
（4.13）
MPk1——k1 频率上包含多路径误差和整周模糊度信息的计算量；
rk1——k1 频率伪距观测量，单位m；
fk1——导航信号 k1 载波的频率，单位MHz；
fk2——导航信号 k2 载波的频率，单位MHz；
jk1——k1 频率载波相位观测量，单位m；
jk2——k2 频率载波相位观测量，单位m；
MPk2——k2 频率上包含多路径误差和整周模糊度信息的计算量；
rk2——k2 频率伪距观测量，单位m。
对于同一颗卫星在连续观测且无周跳的情况下组合的模糊度参数不会变化，在无周跳的多个历元间按公式(4.14)进行计算，得到多路径误差。
（4.14）
MPk
——接收机观测到卫星在 k 频率上多路径误差的评估值；
Nsw ——滑动窗口的历元个数，默认为50；
MPk(t1)
——在历元 t i 接收机观测到卫星在 k 频率上包含多路径误差和整周模糊度信息的计算量。
选取G18卫星和C05卫星为例，分析多路径误差。图1为Mp1指标的时间序列图，图2为Mp2指标的序列图，在两个文本框里分别输出了两个计算量的平均值。
如图4.17和图4.18所示，两个卫星的多路径误差都不超过20m，说明观测值的质量较好。而北斗卫星的多路径误差较GPS卫星来看，其多路径误差变化较为平稳，。说明在观测过程中，C05卫星的观测值质量要更好。D038点位虽然位于湖边，但是环境开阔，并未受到多路径效应较大的影响。而位于树旁的点位，受多路径效应影响反而更大。说明信号遮挡对于观测值的影响更大。
遍历了所有卫星，多路径效应最严重的卫星，多路径效应值可达到50m以上。
执行界面如下：
三、总结
偷偷说一句，没搞懂GLONASS的频率怎么算，所以忽略GLONASS。后续我会将实习报告，完整源代码和程序上传到资源里面，需要的自取。有问题欢迎大家一起讨论！本期的内容就分享到这了。
扫一扫
3215
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
