武汉大学卫星导航原理课程作业——MW，GF组合观测值

于 2024-08-27 16:10:59 发布
阅读量1.6k 收藏 51
点赞数 37
分类专栏： 卫星导航程序设计 文章标签： c++ matlab 算法
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qq_73763569/article/details/141573871
版权

一、作业要求

  • 作业内容
  1. 编程读取rinex观测文件，任选2颗GPS卫星和3颗北斗卫星（其中，北斗选1颗GEO+1颗IGSO+1颗MEO），提取一天中上述卫星的双频伪距和相位观测值（GPS使用C1C、C2W、L1C、L2W；BDS使用C2I、C6I、L2I、L6I）。
  2. 计算MW、GF组合观测值并画出其具体数值随时间变化的曲线图，分析其变化规律。
  3. 在上述基础上，进行历元间求差，得到GF和MW组合的差分序列，绘图展示并分析其变化规律（注: 选择某一连续观测时段，无观测值时不参与计算）
  • 给定文件说明
  1. 观测文件：WUH200CHN_R_20213650000_01D_30S_MO.rnx

（注：WUH2站，30s采样间隔）可在武汉大学IGS数据中心下载。

  1. Rinex 3.04格式说明文档。给出的观测文件格式为rinex 3.04，参数含义参照说明文档附录A2（文件头）、A3（文件体）。
  2. 北斗卫星类型：同学们在选择需要计算的北斗卫星时，可在以下网址查询北斗卫星是否属于MEO/IGSO/GEO

http://www.csno-tarc.cn/system/constellation

二、rinex观测值文件解析

头文件部分

  在计算MW和GF组合观测值时我们要用到的信息包括11-23行，这一段显示了每个卫星记录的数据，便于我们后续做双频组合。

数据块部分

最上面的一行分别是年、月、日、时、分、秒、历元状况（0表示正常）、可视卫星数

然后下面的每一行数据都是对应头文件所示的数据。

三、文件读取与数据赋值

  我选择使用C01和G05两颗卫星作为示例，首先我们定义历元，GPS卫星和BDS卫星的结构体

struct epoch
	{
		string icon;//标志<
		int year;
		int month;
		int day;
		int hour;
		int minute;
		double second;
		int k;
		int number;//卫星数
		double p;//钟偏
	};//存储历元信息的结构体
	vector<epoch> E;
	struct GPS_data
	{
		epoch E;
		double GPST;
		string GPS_name;
		double C1C;
		double C2W;
		double L1C;
		int n1;
		double L2W;
		int n2;
		string data_G[24];
		double MW;//MW组合观测值
		double GF;//GF组合观测值
	};//存储GPS卫星数据的结构体
	vector<GPS_data> G;
	struct BDS_data
	{
		epoch E;
		double GPST;
		string BDS_name;
		double C2I;
		double C6I;
		double L2I;
		int n1;
		double L6I;
		int n2;
		string data_B[32];
		double MW;
		double GF;
	};//存储BDS卫星数据的结构体
	vector<BDS_data> B;

  然后是遍历文件，查找历元行和C01和G05的数据，并与历元对应。由于历元行会给出观测卫星数，所以我们在读取历元行之后获得观测卫星数，在历元行的基础上往下读取相应的行数，这样就能实现数据与历元对应。

  同时在读取每一行数据时，我们不能选择直接用infile读取或者用stringstream解析之后读取，因为会有少量数据缺失，遇到数据缺失时，数据就会出现对应错误，导致MW和GF组合观测值有些算出来会很大。所以选择使用截取字符串的方法，读取了卫星号之后，每次读取16个字符，最后一个数据读取14个字符，就能实现数据的完美对应（RTKlib就是这样读的）。当数据为空时，我们不进行赋值。具体代码如下：

void getdata(string filename)
	{
		ifstream infile(filename);
		string line;
		epoch E1;
		GPS_data G1;
		BDS_data B1;
		while (!infile.eof())
		{
			getline(infile, line);
			if (line.find(">")==0)
			{
				stringstream iss(line);
				iss >> E1.icon >> E1.year >> E1.month >> E1.day >> E1.hour >> E1.minute >> E1.second >> E1.k >> E1.number >> E1.p;
				E.push_back(E1);
			}
			//往下读取k行数据，k为卫星数
			for (int j = 0; j < E1.number; j++)
			{
				getline(infile, line);
				if (line.find("C01") == 0)
				{
					B1.BDS_name = line.substr(0, 3);
					for (int i = 0; i < (line.length() - 1) / 16; i++)
					{
						if (i == (line.length() - 1) / 16)
						{
							B1.data_B[i]=(line.substr(3 + 16 * i + 14));//最后一个数据只读取14个字节
						}
						else
						{
							B1.data_B[i]=(line.substr(3 + i * 16, 16));//每16个字节一个数据
						}
					}
					B1.E = E1;
					if (B1.data_B[3]== "                ")
					{
						continue;
					}
					else
					{
						B1.C2I = stod(B1.data_B[1]);
						B1.C6I = stod(B1.data_B[3]);
						B1.L2I = stod(B1.data_B[17].substr(0, 14));//截取前14个字符为载波相位观测值，最后一个字符为信号强度
						B1.n1 = stoi(B1.data_B[17].substr(14, 2));
						B1.L6I = stod(B1.data_B[19].substr(0, 14));
						B1.n2 = stoi(B1.data_B[19].substr(14, 2));
						B.push_back(B1);
					}
					
				}//选择你所需要的BDS卫星数据
				if (line.find("G05") == 0)
				{
					G1.GPS_name = line.substr(0, 3);
					for (int i = 0; i < (line.length() - 1) / 16; i++)
					{
						if (i == (line.length() - 1) / 16)
						{
							G1.data_G[i]=(line.substr(3 + 16 * i + 14));//最后一个数据只读取14个字节
						}
						else
						{
							G1.data_G[i]=(line.substr(3 + i * 16, 16));//每16个字节一个数据
						}
					}
					G1.E = E1;
					if (G1.data_G[3]== "                ")
					{
						continue;//如果没有数据，则跳出本次for循环
					}
					else
					{
						G1.C1C = stod(G1.data_G[0]);
						G1.C2W = stod(G1.data_G[3]);
						G1.L1C = stod(G1.data_G[12].substr(0, 14));//截取前14个字符为载波相位观测值，最后一个字符为信号强度
						G1.n1 = stoi(G1.data_G[12].substr(14, 2));
						G1.L2W = stod(G1.data_G[15].substr(0, 14));
						G1.n2 = stoi(G1.data_G[15].substr(14, 2));
						G.push_back(G1);
					}
					
				}//选择你所需要的GPS卫星数据
			}

		}
	}

三、MW和GF组合观测值计算

GF组合

编写代码时选的是第二个以m为单位的GF组合观测值

MW组合

也是以m为单位的MW组合

代码如下：

void caculate()
	{
		for (int i = 0; i < B.size(); i++)
		{
			B[i].MW= (c*B[i].L2I - c * B[i].L6I) / (Bf1 - Bf2) - (B[i].C2I*Bf1 + B[i].C6I*Bf2) / (Bf1 + Bf2);//以米为单位的MW组合观测值
			B[i].GF= B[i].C2I - B[i].C6I;//以m为单位的GF组合
			B[i].GPST = CT2GPST(B[i].E);//通用时转换为GPST
		}
		for (int i = 0; i < G.size(); i++)
		{
			G[i].MW = (c*G[i].L1C - c * G[i].L2W) / (Gf1 - Gf2) - (G[i].C1C*Gf1 + G[i].C2W*Gf2) / (Gf1 + Gf2);//以米为单位的MW组合观测值
			G[i].GF = G[i].C1C - G[i].C2W;//以m为单位的GF组合
			G[i].GPST = CT2GPST(G[i].E);//通用时转换为GPST
		}

这一步顺便把时间转换成了以s为单位，方便后续画时序图。

四、计算结果

计算完成之后需要有一个数据的保存，我写了一个保存为文本文件的函数

void savefile(string filename1,string filename2)
	{
		ofstream outfile1(filename1); 
		for (int i = 0; i < G.size(); i++)
		{
			outfile1 << G[i].GPS_name << " " << G[i].GPST << " " << G[i].MW << " " << G[i].GF << endl;
		}
		outfile1.close();
		ofstream outfile2(filename2);
		for (int i = 0; i < B.size(); i++)
		{
			outfile2 << B[i].BDS_name << " " << B[i].GPST << " " << B[i].MW << " " << B[i].GF << endl;
		}
		outfile2.close();
	}

C01卫星计算结果

G05卫星计算结果

五、结果可视化

使用matlab绘图进行计算结果的可视化。

C01卫星MW值和GF值时间序列图

C01卫星MW和GF值差分时间序列图

G05卫星MW值和GF值时间序列图

G05卫星MW和GF值差分时间序列图

Matlab源代码如下：

clear;
data_G=readtable("G05.txt");
data_B=readtable("C01.txt");
fig=figure;
set(fig, 'Position', [200, 100, 800, 350]);
segmentStart = 1;
time=data_G.Var2;
values1=data_G.Var3;
values2=data_G.Var4;
subplot(2,1,1);
hold on;
for i = 2:length(time)
    if (time(i) - time(i-1)) > 300 % 设置时间间隔阈值
        % 绘制当前段数据
        plot(time(segmentStart:i-1), values1(segmentStart:i-1),  'g-');
        segmentStart = i; % 更新分段起点
    end
end
% 绘制最后一段数据
plot(time(segmentStart:end), values1(segmentStart:end),  'g-');
xlabel('时间');
xticks([0,14400,14400*2,14400*3,14400*4,14400*5,14400*6]);
xticklabels({"0:00","04:00","08:00","12:00","16:00","20:00","24:00"})
ylabel('MW组合观测值');
title('G05卫星MW值时间序列图');
hold off;
subplot(2,1,2);
hold on;
segmentStart = 1;
for i = 2:length(time)
    if (time(i) - time(i-1)) > 300 % 设置时间间隔阈值
        % 绘制当前段数据
        plot(time(segmentStart:i-1), values2(segmentStart:i-1),  'r-');
        segmentStart = i; % 更新分段起点
    end
end
% 绘制最后一段数据
plot(time(segmentStart:end), values2(segmentStart:end),  'r-');
xlabel('时间');
xticks([0,14400,14400*2,14400*3,14400*4,14400*5,14400*6]);
xticklabels({"0:00","04:00","08:00","12:00","16:00","20:00","24:00"})
ylabel('GF组合观测值');
title('G05卫星GF值时间序列图');
hold off
%%
values1_diff=diff(values1);
values2_diff=diff(values2);
time_diff=time(2:end);
set(figure,'Position',[200,100,800,350]);
subplot(2,1,1);
segmentStart = 1;
hold on;
for i = 2:length(time_diff)
    if (time_diff(i) - time_diff(i-1)) > 300 % 设置时间间隔阈值
        % 绘制当前段数据
        plot(time_diff(segmentStart:i-1), values1_diff(segmentStart:i-1),  'g-');
        segmentStart = i; % 更新分段起点
    end
end
% 绘制最后一段数据
plot(time_diff(segmentStart:end), values1_diff(segmentStart:end),  'g-');
xlabel('时间');
xticks([0,14400,14400*2,14400*3,14400*4,14400*5,14400*6]);
xticklabels({"0:00","04:00","08:00","12:00","16:00","20:00","24:00"})
ylabel('MW组合观测值差分');
title('G05卫星MW值差分时间序列图');
hold off;
subplot(2,1,2);
hold on;
segmentStart = 1;
for i = 2:length(time_diff)
    if (time_diff(i) - time_diff(i-1)) > 300 % 设置时间间隔阈值
        % 绘制当前段数据
        plot(time_diff(segmentStart:i-1), values2_diff(segmentStart:i-1),  'r-');
        segmentStart = i; % 更新分段起点
    end
end
% 绘制最后一段数据
plot(time_diff(segmentStart:end), values2_diff(segmentStart:end),  'r-');
xlabel('时间');
xticks([0,14400,14400*2,14400*3,14400*4,14400*5,14400*6]);
xticklabels({"0:00","04:00","08:00","12:00","16:00","20:00","24:00"})
ylabel('GF组合观测值差分');
title('G05卫星GF值差分时间序列图');
hold off

六、总结

  本篇文章所有的源代码和结果文件我都打包放在了资源里，可以免费下载，希望大家点点下载，支持一下，谢谢大家。

【免费】武汉大学卫星导航原理课程作业-MW，GF组合观测值.zip资源-CSDN文库

两眼一睁就是学
关注
专栏目录
武大IGS中心下载 导航数据和精密星历数据
努力努力努力
06-14 1107
目录0、各种数据命名一、打开官网连接二、下载导航数据2.1 下载2024年6月12天的导航数据（.n），打开连接后，按照步骤来即可：2.2 下载2024年6月12天的精密轨道文件（.sp3），打开连接后，按照步骤来即可：2.3 下载2024年6月12天的精密钟差文件（.clk），打开连接后，按照步骤来即可：连接如下：四、 RTKLIB代码中的GNSS数据格式、基本概念_rinex-CSDN博客连接如下：武汉大学IGS数据中心 (gnsswhu.cn)
北斗导航 | 观测值的线性组合（公式）
尘世冰封的专栏
09-14 500
基于**不同频率或者不同类型观测值之间的线性组合**，可以形成具有某些特定性能的虚拟观测值。常用的基于双频的线性观测值组合包括宽巷组合、无电离层组合以及MW组合。
RTKLIB中解算时为何使用L2W，不使用L2S
qq_43398531的博客
10-22 637
使用RTKLIB跑数据时，观测值文件里面可能含有不同得信号，如下 下图为Rinex 3.04具体说明有兴趣得可以看一下 使用RTKLIB进行解算时对于观测值得选取有一定得优先级，不同得信号有不同得精度，RTKLIB里面对于优先级得定义如下： 0代表L1，1代表L2.因此对于GPS中L2信号L2W优先级高于L2S因此使用L2W进行解算。 ...
武汉大学IGS数据中心网站的下载方法
weixin_42717467的博客
05-27 4916
目的 进行GNSS解算学习时，从书本中学到了算法。然而网上有不少现成的开源GNSS代码可以使用，这时需要数据样本，导入开源代码如RTKlib@Copyright © 2007-2013, T. Takasu, All rights reserved.来进行学习。数据来源有很多（需要后期整理），武汉大学IGS数据中心是其一。 过程 欲下载能进行双差解算的数据样本。 因此需要下载两个基站的obs文件，同一时间的星历。 打开网站 选择标题栏的下载中心–IGS，不要从右下角的蓝色按钮进去，因为进不去。 下载obs
GPS术语 -- 词汇与概念解释(三)
红茶绿茶和奶茶的博客
06-29 1316
1基频 (Fundamental Frequency) GPS中使用的基频 F为10.23MHz 。L1、 L2载波频率是基频的整数倍。 L1＝154 F＝1575.42 M Hz L2＝120 F＝1227.60MHz 2地心 (Geocenter) 地球质量中心 3大地基准点（ Geodetic Datum） 设计用来最佳拟合一部分或全部大地水准面的一种数学模式。它由大地椭球体及该椭球体与由大地基准原点所决定的地形表面的关系来定义的。这种关系一般（但不是必须）由六个要素来确定：大地纬...
Gps应用知识
weixin_30721077的博客
04-20 670
1． GPS系统组成GPS gloabal Positioning System,这玩意是美国人搞的。主要分三大块，地面的控制站、天上飞的卫星、咱们手里拿的接收机。简单唠叨唠叨先说说设备， 当然大个的都是老美给咱准备好的，地上，有一个主控制站，当然在老美的本土了，在科罗拉多。三个地面天线，五个监测站，分布在全球。主要是收集数据，计算导航信息，诊断系统状态，调度卫星这些杂事。天上，有27颗卫星，距离...
武汉大学卫星导航原理课程作业-MW，GF组合观测值.zip
最新发布
08-27
源代码包括rinex格式的观测值文件读取，MW，GF组合观测值的计算，和结果的可视化源代码，内容齐全。
MW组合观测值序列_lipsg96_MW组合_mW序列_
10-02
"MW组合观测值序列"这一标题暗示了我们正在处理的是一个时间序列的数据集，可能包含了不同时间段、不同波段的观测值。这些数据可能来自多个卫星或雷达系统，经过整合后形成了一种综合观测序列。这种序列可以帮助科学...
GPS卫星组合观测值计算.zip
06-14
《GPS卫星组合观测值计算——基于武汉大学测绘学院GNSS作业》 GPS（全球定位系统）作为现代导航与定位技术的核心，其应用广泛且深入。在武汉大学测绘学院的GNSS（全球导航卫星系统）作业中，学生们通过处理2.11版本...
过程控制课程设计（300MW火电机组协调控制系统的设计）
07-21
通过本课程设计，使学生能较好的运用过程控制的基本概念、基础理论与方法，根据大型火电机组的生产实际，对火电机组的过程控制系统进行分析，设计出原理正确，功能较为全面的300MW火电机组协调控制系统。
GNSS数据下载网址与格式转换
qq_41531197的博客
03-16 4551
数据下载网址： 天线改正文件(atx)下载 ftp://ftp.igs.org/pub/station/general 通用广播星历(brdc/brdm)： ftp://cddis.gsfc.nasa.gov/pub/gps/data/daily/YYYY/brdc ftp://cddis.gsfc.nasa.gov/pub/gps/data/campaign/mgex/daily/rinex3/...
GPS载频信号简述（L1、L2、L3、L5）
热门推荐
Cail466210445的博客
06-22 7万+
L1:1575.42±1.023MHzL2:1227.60±1.023MHz L3:用途和频率暂未公开。资料显示其用于发现核爆炸或其他高能量红外辐射事件的核爆炸侦察系统(NDS)平台提供通讯联系。L5:1176.45±1.023MHz（2009年测试，2010年正式播发）单频接收机只能接收L1载波信号，测定载波相位观测值进行定位。由于不能有效消除电离层延迟影响，只适用于短基线（&lt;15km）的...
武汉大学GNSS中心给IGS提供的数据产品下载地址，包括精密GPS轨道，钟差，EOP，以及实时轨道和钟差
p312011150的博客
09-26 4932
中心对外提供数据服务     中心igs服务器 地址：ftp://ics.gnsslab.cn PANDA软件tables目录，每天更新；  包含我们自己的精密轨道和卫星钟差：whu目录；  igs数据和产品；  IERS数据和产品。   IGS实时数据流 主页：http://ntrip.gnsslab.cn/home caster: ntrip.gnsslab.cn port...
RTKLIB软件介绍
qq_32201015的博客
03-23 1万+
转载自南京信息工程大学遥感与测绘工程学院的赵乐文老师的网课。（来源：B站） 开源GNSS数据处理软件介绍 RTKLIB下载：http://www.rtklib.com CSDN里的资源：快速入门 RTKLIB增强版：https://rtklibexplorer.wordpress.com/ 针对低成本接收机进行算法优化，附有大量的测试过程和数据 RTKLIB二次开发(GAMP): https://geodesy.noaa.gov/gps-toolbox/GAMP.htm 针对PPP数据处理进行算法优化
GPS 信号L2C介绍L2C(L) L2C(M)
技术、思维
03-28 8140
http://en.wikipedia.org/wiki/GPS_signals L2C One of the first announcements was the addition of a new civilian-use signal, to be transmitted on a frequency other than the L1 frequency used
最新IGS MGEX测站数据选择及下载详细教程（附IGS登录方法）
SmartTiger_CSL的博客
10-24 1万+
MGEX数据的下载 MGEX测站选择 1 首先登陆IGS的MGEX官网（ http://mgex.igs.org/index.php），选择“MultiGNSS Stations”查看所有测站分布图，选定大致区域这里选择澳大利亚的测站。 2 在下面的option中可以按照多个字段信息排序，查看各个测站的详细配置信息，包括：测站名，所属的网络、城市、国家和机构，接收机及时钟和天线类型，标定方式，接收的卫星新系统和数据格式及最新日期等： 3 选择其中的3个测站：STR1、STR2和TID1。 打开IGS的
武汉大学GNSS中心给IGS提供的数据产品下载地址
whysnlc的专栏
04-05 1万+
中心对外提供数据服务 中心igs服务器 地址：ftp://ics.gnsslab.cn PANDA软件tables目录，每天更新； 包含我们自己的精密轨道和卫星钟差：whu目录； igs数据和产品； IERS数据和产品。 IGS实时数据流 主页：http://ntrip.gnsslab.cn/home caster: ntri
用matlab写一个MW-GF组合周跳探测代码
06-07
以下是一个用 Matlab 实现的 MW-GF（Moving Window-Geometry-Free）组合周跳探测代码： ```matlab function [L1_valid, L2_valid] = MW_GF_cycle_slip_detection(L1, L2, window_size, threshold) % MW-GF（Moving Window-Geometry-Free）组合周跳探测 % % 参数： % L1: L1 信号的相位观测值 % L2: L2 信号的相位观测值 % window_size: 滑动窗口大小 % threshold: 周跳探测门限 % 返回值： % L1_valid: 经过周跳探测后的 L1 信号相位观测值 % L2_valid: 经过周跳探测后的 L2 信号相位观测值 % 初始化变量 L1_valid = L1; L2_valid = L2; n = length(L1); L1_diff = zeros(n, 1); L2_diff = zeros(n, 1); for i = window_size+1:n % 计算 L1 和 L2 的差分值 L1_diff(i) = L1(i) - L1(i-window_size); L2_diff(i) = L2(i) - L2(i-window_size); % 计算 MW-GF 残差 MW = (L1_diff(i) + L2_diff(i)) / 2; GF = (L1_diff(i) - L2_diff(i)) / 2; % 判断是否发生周跳 if abs(MW) > threshold && abs(GF) < threshold % 修复周跳 L1_valid(i:end) = L1_valid(i:end) - L1_diff(i); L2_valid(i:end) = L2_valid(i:end) - L2_diff(i); end end end ``` 这个算法是基于 MW-GF 组合算法实现的，使用一个滑动窗口来计算 L1 和 L2 的差分值，然后计算 MW-GF 残差，根据门限判断是否发生周跳，并将周跳修复。其中，`window_size` 参数表示滑动窗口的大小，`threshold` 参数表示周跳探测的门限，可以根据具体情况进行调整。
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

两眼一睁就是学

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付：¥1
获取中
扫码支付

您的余额不足，请更换扫码支付或充值

打赏作者

实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值