#include "stdafx.h"
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <algorithm>
#include <sstream>
#include <iomanip>
#include <E:/Study software/eigen
-3.3.9/Eigen/Dense>
using namespace Eigen;
using namespace std;
const
double c=299792458; //light speed
const
double GM=3.986005e14;
const
double we=7.292115e
-5; //地球自转速度
//
-------------------------------------------------------------------
typedef vector<vector<
double>> Mat;
//民用时转化为GPS时
void gpstime
(int year,int month,int day,int hour, int minute,
double second, int &week,
double &secofweek)
{
int ry[12]={31,29,31,30,31,30,31,31,30,31,30,31},fry[12]={31,28,31,30,31,30,31,31,30,31,30,31};
if
(year > 80)
year = 1900 + year;
else
year = 2000 + year;
int z=0;
for
(int i1=1980; i1 < year+1; i1++)
{
if
(i1%4 == 0 && i1%400 != 0)
z++; //记闰年数
}
if
(year%4 == 0 && year%100 != 0)
{
for
(int i2 = 0; i2 < month
-1; i2++)
{
day += ry[i2];
}
}
else
{
for
(int i3=0; i3 < month
-1; i3++)
{
day += fry[i3];
}
}
day +=
(year
-1980)*365 + z
- 5;
week = day/7;
secofweek= day%7*24*3600 + hour*3600 + minute*60 + second ;
}
//观测文件类
class Rea
dObsFile //观测文件类
{
public:
double X0,Y0,Z0; //估计坐标
int epochNum;
vector<int>svNum; //卫星个数
vector<vector<
double>>c1,p1,p2,l1,l2,d1,d2;
vector<vector<int>>prn;//卫星的prn号二维向量
int week;
vector<
double>secofweek;
void Rea
dObsHeader
(string ofile)
{
//用文件流读文件
ifstream in_o
(ofile);
if
(in_o)
cout<<"打开O文件成功"<<endl;
else
cerr<<"未找到O文件"<<endl;
string str;
//定位到初始坐标
while(str.find
("APPROX POSITION XYZ") == string::npos)
{
getline
(in_o,str);
}
stringstream ss
(str); //用字符串流提取出来初始坐标
ss>>X0>>Y0>>Z0; //初始坐标提取
//
--------------------------------------------------------------------
//定位到数据头
while(str.find
("END OF HEADER") == string::npos)
{
getline
(in_o,str);
}
//计算历元数
epochNum=0;
while(in_o)
{
getline
(in_o,str);
for
(int i = 0;
(size_t)i < str.length
(); i++)
if
(str[i] == 'G')
{
epochNum++;
break;
}
}
} //end of function Rea
dObsHeader
void Rea
dObsData
(string ofile)
{
//定位到文件头
ifstream in_o1
(ofile);
string str;
//定位到数据头
while(str.find
("END OF HEADER") == string::npos)
{
getline
(in_o1,str);
}
//用历元数定义向量的列数
c1.resize
(epochNum);
p1.resize
(epochNum);
p2.resize
(epochNum);
l1.resize
(epochNum);
l2.resize
(epochNum);
d1.resize
(epochNum);
d2.resize
(epochNum);
prn.resize
(epochNum);
svNum.resize
(epochNum);
secofweek.resize
(epochNum);
vector<
double> second
(epochNum);
int year = 0,month = 0,day = 0,hour = 0,minute = 0,mark = 0;
week=0;
char g;
for
(int i=0;i<epochNum;i++) //每个历元
{
//week
和secofweek提取
in_o1>>year>>month>>day>>hour>>minute>>second[i]>>mark;
gpstime
(year,month,day,hour,minute,second[i],week,secofweek[i]);//转化为周
和周内秒
//提取卫星数目
和prn号
getline
(in_o1,str);
stringstream sss
(str);
sss>>svNum[i];//卫星数目
//用卫星数定义向量的列数
c1[i].resize
(svNum[i]);
p1[i].resize
(svNum[i]);
p2[i].resize
(svNum[i]);
l1[i].resize
(svNum[i]);
l2[i].resize
(svNum[i]);
d1[i].resize
(svNum[i]);
d2[i].resize
(svNum[i]);
prn[i].resize
(svNum[i]);
for
(int k=0; k<svNum[i]; k++)
{
sss>>g>>prn[i][k];//prn号
}
//读取观测数据
for
(int j = 0; j < svNum[i]; j++)
{
in_o1>>c1[i][j]>>p1[i][j]>>p2[i][j]>>l1[i][j]>>l2[i][j]>>d1[i][j]>>d2[i][j];
}
getline
(in_o1,str);
}// end of for epochNum
}// end of Rea
dObsData
}; //end of class Rea
dObsFile
//星历文件类
class ReadNavFile //星历文件类
{
public:
int svNum;//星历
中的卫星个数
vector<vector<
double>>ephlist;
void ReadNavData
(string nfile)
{
//
--------------------------------------------------------------------
//用文件流读文件
ifstream in_n
(nfile);
if
(in_n)
cout<<"打开N文件成功"<<endl;
else
cerr<<"未找到N文件"<<endl;
string str;
//
--------------------------------------------------------------------
//定位到数据头
while(str.find
("END OF HEADER") == string::npos)
{
getline
(in_n,str);
}
//计算文件的行数
和段数
int count = 0;
while(in_n)
{
getline
(in_n, str);//从文件中读取一行
remove
(str.begin
(), str.end
(), ' ');//这个算法函数在algorithm头文件中，删除一行
中的空格
remove
(str.begin
(), str.end
(), '\t');//删除一行
中的制表符，因为制表符
和空格都是空的
if
(str.length
() > 0)
{
//如果删除制表符
和空格之后的一行数据还有其他字符就算有效行
count ++;
}
}
svNum =count/8;//卫星个数,每个卫星8行数据
ephlist.resize
(svNum);//用卫星个数定义星历的行数
//定位到文件头
ifstream in_n1
(nfile);
//定位到数据头
while(str.find
("END OF HEADER") == string::npos)
{
getline
(in_n1,str);
}
for
(int i=0; i < svNum; i++)
{
ephlist[i].resize
(38);//星历的列数=38
getline
(in_n1,str);
ephlist[i][0] = atoi
(str.substr
(0,2).c_str
());
ephlist[i][1] = atoi
(str.substr
(3,2).c_str
());
ephlist[i][2] = atoi
(str.substr
(7,1).c_str
());
ephlist[i][3] = atoi
(str.substr
(10,1).c_str
());
ephlist[i][4] = atoi
(str.substr
(13,1).c_str
());
ephlist[i][5] = atoi
(str.substr
(15,2).c_str
());
ephlist[i][6] = atof
(str.substr
(19,4).c_str
());
ephlist[i][7] = atof
(str.substr
(23,19).c_str
());
ephlist[i][8] = atof
(str.substr
(42,19).c_str
());
ephlist[i][9] = atof
(str.substr
(61,19).c_str
());
for
(int k = 0; k < 7; k++) //一共7行，一行行读取
switch
(k)
{
case 0:
{
getline
(in_n1,str);
ephlist[i][10] = atof
(str.substr
(3,19).c_str
());
ephlist[i][11] = atof
(str.substr
(23,19).c_str
());
ephlist[i][12] = atof
(str.substr
(42,19).c_str
());
ephlist[i][13] = atof
(str.substr
(61,19).c_str
());
break;
}
case 1:
{
getline
(in_n1,str);
ephlist[i][14] = atof
(str.substr
(3,19).c_str
());
ephlist[i][15] = atof
(str.substr
(23,19).c_str
());
ephlist[i][16] = atof
(str.substr
(42,19).c_str
());
ephlist[i][17] = atof
(str.substr
(61,19).c_str
());
break;
}
case 2:
{
getline
(in_n1,str);
ephlist[i][18] = atof
(str.substr
(3,19).c_str
());
ephlist[i][19] = atof
(str.substr
(23,19).c_str
());
ephlist[i][20] = atof
(str.substr
(42,19).c_str
());
ephlist[i][21] = atof
(str.substr
(61,19).c_str
());
break;
}
case 3:
{
getline
(in_n1,str);
ephlist[i][22] = atof
(str.substr
(3,19).c_str
());
ephlist[i][23] = atof
(str.substr
(23,19).c_str
());
ephlist[i][24] = atof
(str.substr
(42,19).c_str
());
ephlist[i][25] = atof
(str.substr
(61,19).c_str
());
break;
}
case 4:
{
getline
(in_n1,str);
ephlist[i][26] = atof
(str.substr
(3,19).c_str
());
ephlist[i][27] = atof
(str.substr
(23,19).c_str
());
ephlist[i][28] = atof
(str.substr
(42,19).c_str
());
ephlist[i][29] = atof
(str.substr
(61,19).c_str
());
break;
}
case 5:
{
getline
(in_n1,str);
ephlist[i][30] = atof
(str.substr
(3,19).c_str
());
ephlist[i][31] = atof
(str.substr
(23,19).c_str
());
ephlist[i][32] = atof
(str.substr
(42,19).c_str
());
ephlist[i][33] = atof
(str.substr
(61,19).c_str
());
break;
}
case 6:
{
getline
(in_n1,str);
ephlist[i][34] = atof
(str.substr
(3,19).c_str
());
break;
}
}
}
}//end of function ReadNavData
()
};// end of class ReadNavFile
//计算卫星位置函数
void satpos
(vector<vector<
double>>ephlist,int prn,
double t,
double *X,
double *Y,
double *Z,
double *dts)
{
int k;
for
(k=0;k<ephlist.size
();k++) //搜索星历，找到卫星prn所在的行数k
{
if
(ephlist[k][0]==prn)
break;
if
(k==ephlist.size
())
{
cout<<"未找到对应星历"<<endl;
break;
}
}
double af0=ephlist[k][7];
double af1=ephlist[k][8];
double af2=ephlist[k][9];
double IODE=ephlist[k][10];
double CRS=ephlist[k][11];
double DELTA_N=ephlist[k][12];
double M_0=ephlist[k][13];
double CUC=ephlist[k][14];
double e=ephlist[k][15];
double CUS=ephlist[k][16];
double SQRT_A=ephlist[k][17];
double toe=ephlist[k][18];
double CIC=ephlist[k][19];
double OMEGA_0=ephlist[k][20];
double CIS=ephlist[k][21];
double I_0=ephlist[k][22];
double CRC=ephlist[k][23];
double OMEGA=ephlist[k][24];
double OMEGA_
DOT=ephlist[k][25];
double I_
DOT=ephlist[k][26];
double dt=af0+af1*
(t
-toe)+af2*pow
((t
-toe),2);//卫星钟差改正
*dts=dt;
t=t
-dt;
double n0=sqrt
(GM*1.0)/pow
((SQRT_A),3);
double n=n0+DELTA_N;
double tk=t
-toe;
double M=M_0+n*tk;
double E=M;
double E0;
for
(int j=0;j<10;j++)
{
E0=E;
E=M+e*sin
(E);
if
(fabs
(E0
-E)<1e
-12)
break;
}
double f=atan
(sqrt
((1+e)/
(1
-e))*tan
(E/2))*2;
double u0=OMEGA+f;
double Tu=CUC*cos
(2*u0)+CUS*sin
(2*u0);
double Tr=CRC*cos
(2*u0)+CRS*sin
(2*u0);
double Ti=CIC*cos
(2*u0)+CIS*sin
(2*u0);
double u=u0+Tu;
double r=SQRT_A*SQRT_A*
(1
-e*cos
(E))+Tr;
double i=I_0+Ti+I_
DOT*tk;
double x=r*cos
(u);
double y=r*sin
(u);
double L=OMEGA_0+OMEGA_
DOT*
(t
-toe)
-we*t;
*X=x*cos
(L)
-y*cos
(i)*sin
(L);
*Y=x*sin
(L)+y*cos
(i)*cos
(L);
*Z=y*sin
(i);
}
//计算信号发射时刻函数
double emitTime
(double X0,
double Y0,
double Z0,int prn, vector<vector<
double>>ephlist,
double pseu
dorange,
double tr) //信号发射时刻计算
{
double t0;
double t1=0;
double dtr=0; //钟差初始化
t0=pseu
dorange/c
-dtr;//信号传播时间
//cout<<"初始信号传播时间"<<t0<<endl;
int iter=0;//迭代次数
double emitTime;
double Xs,Ys,Zs; //卫星位置
double dts;
double R;
double temp_t0=0;
while (fabs
(t0
-temp_t0)>1e
-7)
{
temp_t0=t0;
emitTime=tr
-t0;//信号发射时刻=接受时刻
-传播时间
satpos
(ephlist , prn , emitTime, &Xs, &Ys, &Zs, &dts);
//cout<<"卫星位置"<<Xs<<" "<<Ys<<" "<<Zs<<" "<<endl;
R=sqrt
(pow
(Xs
-X0,2)+pow
(Ys
-Y0,2)+pow
(Zs
-Z0,2)); //几何距离
t1=R/c;
t0=t1;
//cout<<"信号传播时间"<<t0<<endl;
iter++;
if
(iter>10)
{
cout<<"卫星发射时刻迭代未收敛"<<endl;
break;
}
}
return emitTime;
}
int _tmain
(int argc, _TCHAR* argv[])
{
string ofile="SITE247J.01O"; //O文件
string nfile="SITE247J.01N"; //N文件
//string ofile="ffff13050802.13O"; //O文件
//string nfile="ffff13050802.13N"; //N文件
//读取观测值文件
Rea
dObsFile Obs;
Obs.Rea
dObsHeader
(ofile);//读取观测文件头，获得接收机估计位置
和卫星个数
Obs.Rea
dObsData
(ofile);//读取所有历元的观测数据
Obs.X0=0;Obs.Y0=0;Obs.Z0=0;
cout<<"接收机估计位置： "<<Obs.X0<<" "<<Obs.Y0<<" "<<Obs.Z0<<endl;
cout<<"历元个数： "<<Obs.epochNum<<endl;
cout<<"起始观测时间： "<<Obs.week<<" "<<Obs.secofweek[0]<<endl;
//读取星历文件
ReadNavFile Eph;
Eph.ReadNavData
(nfile); //读取星历到数组Eph.ephlist[][]
/*
for
(int i=0;i<7;i++)
{
for
(int j=0;j<35;j++)
{
cout<<Eph.ephlist[i][j]<<" ";
}
cout<<endl;
}
*/
vector<vector<
double> > Pos
(Obs.epochNum,vector<
double>
(3)); //储存每个历元的接收机位置
double sumX=0,sumY=0,sumZ=0;
//逐个历元解算
for
(int epoch=0;epoch<Obs.epochNum;epoch++)
{
vector<vector<
double>>satPosition; //用来储存每颗卫星位置
和卫星钟差svNum*4
satPosition.resize
(Obs.svNum[epoch]);
//int epoch=1;
//计算每颗卫星的信号发射时刻
vector<
double>secOfEmit
(Obs.svNum[epoch]);
for
(int sv=0;sv<Obs.svNum[epoch];sv++)
{
satPosition[sv].resize
(4);//卫星位置+卫星钟差
double tr=Obs.secofweek[epoch]*1.0;
//使用p1 p2组成的双频改正
模型
secOfEmit[sv]=emitTime
(Obs.X0, Obs.Y0, Obs.Z0, Obs.prn[epoch][sv],Eph.ephlist,2.54573*Obs.p1[epoch][sv]
-1.54573*Obs.p2[epoch][sv],tr);
}
//计算每颗卫星的位置
for
(int ii=0;ii<Obs.svNum[epoch];ii++)
{
satpos
(Eph.ephlist,Obs.prn[epoch][ii],secOfEmit[ii],&satPosition[ii][0],&satPosition[ii][1],&satPosition[ii][2],&satPosition[ii][3]);
}
/*
cout<<"每颗卫星的位置"<<endl;
for
(int jj=0;jj<Obs.svNum;jj++)
{
//cout<<setprecision
(8)<<satPosition[jj][0]<<" "<<satPosition[jj][1]<<" "<<satPosition[jj][2]<<endl;
}
*/
//平差计算接收机位置
MatrixXd B
(Obs.svNum[epoch],4);
MatrixXd L
(Obs.svNum[epoch],1);
MatrixXd R
(Obs.svNum[epoch],1);
MatrixXd dX
(4,1);
MatrixXd BtB
(4,4);
MatrixXd BtL
(4,1);
dX
(3,0)=0;//初始化接收机钟差，单位米
int iteration=0;
do //迭代平差
{
for
(int i=0;i<Obs.svNum[epoch];i++) //给B L矩阵赋值
{
R
(i,0)=sqrt
(pow
(Obs.X0
-satPosition[i][0],2)+pow
(Obs.Y0
-satPosition[i][1],2)+pow
(Obs.Z0
-satPosition[i][2],2));
B
(i,0)=
(Obs.X0
-satPosition[i][0])/R
(i,0);
B
(i,1)=
(Obs.Y0
-satPosition[i][1])/R
(i,0);
B
(i,2)=
(Obs.Z0
-satPosition[i][2])/R
(i,0);
B
(i,3)=1;
L
(i,0)=2.54573*Obs.p1[epoch][i]
-1.54573*Obs.p2[epoch][i]
-R
(i,0)
-dX
(3,0)+satPosition[i][3]*c;//p1 p2双频改正
模型
}
BtB=B.transpose
()*B;
BtL=B.transpose
()*L;
dX=BtB.inverse
()*BtL;
//cout<<"dX"<<endl;
//cout<<dX<<endl;
//cout<<" "<<endl;
Obs.X0=Obs.X0+dX
(0,0);
Obs.Y0=Obs.Y0+dX
(1,0);
Obs.Z0=Obs.Z0+dX
(2,0);
iteration++;
if
(iteration>10)
{
cout<<"迭代未收敛"<<endl;
//break;
}
}
while(fabs
(dX
(0,0))>0.001&&fabs
(dX
(1,0))>0.001&&fabs
(dX
(2,0))>0.001); //迭代平差 end of
while
cout<<"第"<<epoch+1<<"个历元:"<<endl;
cout<<"iteration="<<iteration<<endl;
cout<<"定位结果："<<setprecision
(8)<<Obs.X0<<" "<<Obs.Y0<<" "<<Obs.Z0<<endl;
cout<<endl;
sumX=sumX+Obs.X0;
sumY=sumY+Obs.Y0;
sumZ=sumZ+Obs.Z0;
}//end of for
(epoch)
//计算平均位置
double X,Y,Z;
X=sumX/
(Obs.epochNum*1.0);
Y=sumY/
(Obs.epochNum*1.0);
Z=sumZ/
(Obs.epochNumx*1.0);
cout<<"平均接收机位置: "<<setprecision
(8)<<X<<" "<<Y<<" "<<Z<<endl;
system
("pause");
return 0;
}学习一下这段代码。然会交给我如何继续编写下面的这段代码。1、以下是reader.h的代码：#pragma once
#ifndef READER_H
#define READER_H
//头文件保护指令,避免多重包含导致的重复定义错误（如类重复定义、函数重声明等）
#define _CRT_SECURE_NO_DEPRECATE
#define _CRT_SECURE_NO_WARNINGS
#define _SCL_SECURE_NO_WARNINGS
//这三行是微软VC++编译器警告控制宏
#include <iostream>
#include <fstream>
#include <stdio.h>
#include <string>
#include <sstream>
#include <vector>
#include <cstddef>
#include <cstring>
#include <array>
#include <map>
#include <utility>
#include <iomanip> // set precision 设置精度
#include <math.h> /* pow */
#include <E:/Study software/eigen
-3.3.9/Eigen/Dense>
#include <E:/Study software/eigen
-3.3.9/Eigen/Core>
constexpr int MAX_STRLEN = 100; // 字符串长度
constexpr int MAX_BYTE = 1024; // 字节
constexpr
double PI = 3.14159265358979323846264;
using namespace std;
using namespace Eigen;
class reader
{
public:
reader
(const char* obsFilename, const char* navFilename);
//分别用于观测文件
和导航文件的文件指针
FILE* pFile; // File pointer for observation 用于观测的文件指针
FILE* NavpFile; // File pointer for navigation 用于导航的文件指针
//Flags
bool obsHeader = true;//表示程序正在读取观测文件（.obs）的头部信息
bool navHeader = true;//表示程序正在读取导航文件（.nav）的头部信息
int SENTINEL = 0;//定义一个哨兵值
struct Header
{
string version; //版本
string markername; //标记点名称
int nobstype; //观测类型数量 number
vector< string > obstype; // 存储观测类型字符串
double interval; //间隔
vector<
double > antpos; //存储天线位置数据
};
Header head; //声明了一个head的变量
struct RECORDS
{
vector< vector<
double> > record;//二维动态数组（向量的向量）
double time;
double YEAR;
double MONTH;
double DAY;
double HOUR;
double MIN;
double Indivi_Sec;// 当天累计秒数
double SEC;
vector<int> prn;// 卫星PRN编号列表
int NumOfSat;// 卫星数量
};
RECORDS obs;//声明了一个obs的变量
struct NAVRECORDS
{
double
PRN, //卫星的PRN号
(1,0)
SVClockBias, //卫星钟偏差
(1,2)
SVClockDrift, //卫星钟漂移
(1,3)
SVClockDriftRate,//卫星钟漂移速度
(1,4)
IODE, //星历数据发布号
Crs, //轨道半径改正项
(2,2)
DeltaN, //平均角速度改正项
(2,3)
Mo, //平近点角
(2,4)
Cuc, //升交点角距改正项
(3,1)
Eccentricity, //轨道偏心率
Cus, //升交点角距改正项
Sqrta, //轨道长半轴平方根
Toe, //星历的参考时刻
(4,1)
Cic, //轨道倾角的改正项
OMEGA, //升交点经度
Cis, //轨道倾角改正项
Io, //轨道倾角
(5,1)
Crc, //轨道半径的改正项
omega, //近地点角距
OMEGA
DOT, //升交点赤经变化率
I
DOT, //轨道倾角的变率
(6,1)
L2CodesChannel, // L2频道C/A码标识
GPSWeek, // GPS时间周
L2PDataFlag, //L2P码标识
SVAccuracy, //卫星精度
(7,1)
SVHealth, //卫星健康
TGD, //电离层延迟
IODC, //星钟的数据质量
TransmissionTime,//信息发射时间
(8,1)
FitInterval, //星历拟合区间
(8,2)
SEC, //GPS 周内秒
gpstime; //GPS 时间
//历元：TOC中卫星钟的参考时刻
int TOC_Y;//年
int TOC_M;//月
int TOC_D;//日
int TOC_H;//时
int TOC_Min;//分
int TOC_Sec;//秒
};
//函数声明
double GPSTime
(vector<
double> time); //时间转换为GPS秒数的函数
double scientificTo
double
(string str); //用于将科学计数法字符串转换为双精度浮点数的函数
void rinexOBSReader
(Header& head, RECORDS& obs_RECORDS, int& SENTINEL); //RINEX 格式的观测文件读取（头文件，记录文件，哨兵值）
void rinexNAVReader
(vector<NAVRECORDS>& nav_RECORDS); //RINEX 格式的导航文件读取（记录文件动态数组）
};
#define FILE_SIZE 200
#endif 2、以下是reader.cpp的代码：#define _CRT_SECURE_NO_WARNINGS
#include <iostream>
using namespace std;
#include <string>
#include "reader.h"
#define _CRT_SECURE_NO_DEPRECATE
#define _CRT_SECURE_NO_WARNINGS
#define _SCL_SECURE_NO_WARNINGS
// 此函数计算以秒为单位的GPS时间
double reader::GPSTime
(vector<
double> time)
{
double julian_date, JDholder1, JDholder2, gpsweekholder, y, m, gpsweek, gpstime;
//儒略日、计算儒略日时的中间变量、计算儒略日时的另一个中间变量、计算 GPS 周时的中间变量、年、月、gps周数、gps秒数
//time[0]= 年（后两位）、time[1]= 月、time[2]= 日、time[3]= 时、time[4]= 分、time[5]= 秒。
double UT = time[3] +
(time[4] / 60) +
(time[5] / 3600);//UT = 小时 + 分钟/60 + 秒/3600
//儒略日计算需要将 1
-2 月视为上一年的 13
-14 月（历史算法约定），因此需要调整：
if
(time[1] > 2) {//即3
-12月
y = time[0] + 2000;//这里的年是两位数，如time[0]=23 → 2023
m = time[1];//月份不变
}
else {//如果月份为1月
和2月
y = time[0] + 2000
- 1;//年份
-1
m = time[1] + 12;//月份＋12
}
JDholder1 = 365.25 * y; //计算年贡献（包含闰年因素的年累计天数）
JDholder2 = 30.6001 *
(m + 1);//计算月贡献
if
(JDholder1 > 0 && JDholder2 > 0) {//如果两个均为正数，floor为向下取整
julian_date = floor
(JDholder1) + floor
(JDholder2) + time[2] +
(UT / 24) + 1720981.5;
}
else {//否则向上取整
julian_date = ceil
(JDholder1) + ceil
(JDholder2) + time[2] +
(UT / 24) + 1720981.5;
}
gpsweekholder =
(julian_date
- 2444244.5) / 7;//计算gps总周数
if
(gpsweekholder > 0) {
gpsweek = floor
(gpsweekholder);//如果大于0就向下取整
}
else {
gpsweek = ceil
(gpsweekholder);//如果小于0就向上取整
}
gpstime = round
((((julian_date
- 2444244.5) / 7)
- gpsweek) * 7 * 24 * 3600);
return gpstime;
//返回gpstime，即当前时间在 GPS 周内的总秒数（范围 0~604799，因为 1 周 = 604800 秒）。
//整个过程本质是：本地时间 → 儒略日 → GPS 时间（周 + 秒） 的转换，用于 GPS 系统的时间统一。
}
double reader::scientificTo
double
(string str)//定义成员数返回类型为
double，输入字符串为采用科学计数法的数
//因为在导航文件中科学计数法采用的是用D代表指数，所以需要转换为c++中可以识别的e
{
istringstream os
(str.replace
(15, 1, "e")); //replace D with e in rinex Nav file //在RINEX导航文件中用e替换D
//在输入字符串str的第 15 个位置（索引从 0 开始），替换 1 个字符为"e"。
double str1;
os >> str1;
return str1;
}
reader::reader
(const char* obsFilename, const char* navFilename)//在reader类中构造了一个reader函数（接收o文件路径
和n文件路径，const保证不可修改）
{
pFile = fopen
(obsFilename, "rt"); //pFile用于保存打开的obs文件句柄
if
(pFile == NULL) perror
("Error opening file");
else cout << "观测文件打开成功！" << endl;
NavpFile = fopen
(navFilename, "rt");
if
(NavpFile == NULL) perror
("Error opening file");
else cout << "导航文件打开成功！" << endl;
}
//constexpr int MAX_STRLEN = 100;
//constexpr int MAX_BYTE = 1024;
void reader::rinexOBSReader
(Header& head, RECORDS& obs_RECORDS, int& SENTINEL) //RINEX 格式的obs文件读取（头文件，记录文件，哨兵值）
{
char line[MAX_STRLEN];
if
(obsHeader == true)
{
fgets
(line, MAX_STRLEN, pFile);//打开文件读取第一行
while ((!feof
(pFile)) &&
(strstr(line, "END OF HEADER") == NULL))
//// 循环：直到文件结束 或 读到"END OF HEADER"（头部结束）
{
if
(sizeof
(line) >= 61)
{
if
(strstr(line, "RINEX VERSION"))//解析 RINEX 版本号
{
string sline
(line);
string token = sline.substr
(5, 3); // from 5th character and five spaces // 从第5个字符开始，取3个字符（版本号位置）
head.version = token; // 存储到Header结构体的version成员
}
if
(strstr(line, "MARKER NAME"))//站点标记名
{
string sline
(line); // convert the line
(char) to string class type sLine
string token = sline.substr
(0, 10); // from 5th character and five spaces// 从第0个字符开始，取10个字符（站点名位置）
head.markername = token;// 存储到markername成员
}
if
(strstr(line, "INTERVAL"))//采样间隔
{
string sline
(line); // convert the line
(char) to string class type sLine
string token = sline.substr
(4, 6); // from 5th character and five spaces
head.interval = stod
(token);
}
if
(strstr(line, "APPROX POSITION XYZ")) //近似坐标
{
string sline
(line); // convert the line
(char) to string class type sLine
string token = sline.substr
(2, 56); // from 3rd character upto 56 提取XYZ坐标的字符串（长度56）
char cline[MAX_BYTE]; // cline is a char line
strncpy
(cline, token.c_str
(), sizeof
(cline));
//strncpy ：复制最多 sizeof
(cline) 个字符到 cline 。
cline[sizeof
(cline)
- 1] = 0;//确保字符数组 cline 以空字符（\0）结尾
for
(char* r = strtok
(cline, " "); r; r = strtok
(NULL, " "))//strtok是分割符，在cline中找空格进行分割，当r=null即r为空时就停止分割
{
double rr = atof
(r);//atof ：把字符串 "4916147.2345" 变成数值 4916147.2345 。
head.antpos.push_back
(rr);//尾插法插入rr
}
}
if
(line[0] == 'G' &&
strstr(line, "# / OBS TYPES"))//解析观测类型
{
string sline_new
(line); // convert the line
(char) to string class type sLine
string token_new = sline_new.substr
(4, 2); // from 5th character and two spaces
head.nobstype = stoi
(token_new);//利用stoi把 token_new 里的数字文本转成真正的整数，再保存到 nobstype 。
// 若观测类型数量≤13，当前行即可容纳所有类型
if
(head.nobstype <= 13)
{
string sline
(line); // convert the line
(char) to string class type sLine
string token = sline.substr
(7, 51); // from 10th character upto
char cline[MAX_BYTE]; // cline is a char line
strncpy
(cline, token.c_str
(), sizeof
(cline));
cline[sizeof
(cline)
- 1] = 0;
for
(char* p = strtok
(cline, " "); p; p = strtok
(NULL, " "))
{
head.obstype.push_back
(p);
}
}
}
}
//......... Get a new line.......................
fgets
(line, MAX_STRLEN, pFile);//读取下一行
}
}
obsHeader = false;//标记 “RINEX obs文件的头部已经读取完毕”
// 输出头部信息
cout << "RINEX版本号: " << head.version << endl;
cout << "站点标记名: " << head.markername << endl;
cout << "采样间隔: " << head.interval << endl;
cout << "近似坐标XYZ: ";
for
(double coord : head.antpos) {
cout << coord << " ";
}
cout << endl;
cout << "观测类型数量: " << head.nobstype << endl;
cout << "观测类型: ";
for
(const std::string& type : head.obstype) {
cout << type << " ";
}
//Read the observation file......读取观测文件的实际观测数据
//一行一行读文件，找到包含卫星观测数据的行，解析这个时间点（历元）的时间、卫星信息、每个卫星的观测值，最后把这些数据存起来。
// 读取观测文件的实际观测数据
// 循环读取每个包含>且历元标志为0的行，提取时间信息，
// 并获取该历元行到下一个包含>的行之间所有首字符为'G'的行数据
int Nobs = head.nobstype; // 观测类型数量
int Num_Sat; // 卫星数量
//char line[MAX_STRLEN]; // 存储读取的行
while (!feof
(pFile)) // 外层循环：处理所有历元直到文件结束
{
bool foundEpoch = false;
double timeStamp = 0.0, timeSEC = 0.0; // 时间相关变量
// 第一步：寻找历元行时，修改读取逻辑（记录每行起始位置）
while (true) { // 用无限循环+break替代原循环，更清晰控制指针
long epochStartPos = ftell
(pFile); // 记录当前行的起始位置（关键！）
if
(fgets
(line, MAX_STRLEN, pFile) == nullptr) break; // 读取失败则退出
if
(strstr(line, ">") != nullptr) { // 找到历元行
std::string sline
(line);
std::vector<std::string> tempEpoch;
std::stringstream ss
(sline);
std::string token;
ss >> token; // 跳过>
while (ss >> token && tempEpoch.size
() < 7) {
tempEpoch.push_back
(token);
}
// 检查是否有足够的字段且历元标志为0
if
(tempEpoch.size
() >= 7 && tempEpoch[6] == "0")
{
// 提取时间信息（年、月、日、时、分、秒）
try
{
double Y = stod
(tempEpoch[0]);
double M = stod
(tempEpoch[1]);
double D = stod
(tempEpoch[2]);
double H = stod
(tempEpoch[3]);
double Min = stod
(tempEpoch[4]);
double Sec = stod
(tempEpoch[5]);
// 计算时间戳
和当天总秒数
vector<
double> timeHolder;
timeHolder.push_back
(Y);
timeHolder.push_back
(M);
timeHolder.push_back
(D);
timeHolder.push_back
(H);
timeHolder.push_back
(Min);
timeHolder.push_back
(Sec);
timeStamp = GPSTime
(timeHolder);
timeSEC = H * 3600 + Min * 60 + Sec;
foundEpoch = true; // 标记找到有效历元行
/*
// 输出历元时间信息
cout << endl << "
-------------------------" << endl;
cout << "找到有效历元行" << endl;
cout << "历元时间: " << Y << "
-" << tempEpoch[1] << "
-" << tempEpoch[2]
<< " " << tempEpoch[3] << ":" << tempEpoch[4] << ":" << tempEpoch[5] << endl;
cout << "时间戳: " << timeStamp << endl;
cout << "当天总秒数: " << timeSEC << endl;
*/
break;
}
catch
(...)
{
// 处理时间转换异常
continue;
}
}
}
}
// 如果未找到有效历元行，结束循环
if
(!foundEpoch)
break;
// 第二步：提取当前历元的G行（在找到有效历元后）
bool reachNextEpoch = false;
long nextEpochStartPos =
-1; // 下一个历元行的起始位置
while (!reachNextEpoch) {
long currentPos = ftell
(pFile); // 记录当前读取位置
if
(fgets
(line, MAX_STRLEN, pFile) == nullptr) break;
if
(strstr(line, ">") != nullptr) { // 找到下一个历元行
nextEpochStartPos = currentPos; // 记录其起始位置（关键！）
reachNextEpoch = true;
break;
}
// 处理首字符为'G'的行
if
(line[0] == 'G')
{
std::string gLine
(line);
double gpsType = 0.0; // 第2
-3个字符组成的数值
double dataField = 0.0; // 第6
-19个字符组成的数值
// 提取第1
-3个字符（索引1
和2）
if
(gLine.length
() >= 20)
{
try
{
std::string typeStr = gLine.substr
(1, 2);
gpsType = stod
(typeStr);
}
catch
(...)
{
// 处理类型转换异常
}
}
int gpsnub = static_cast<int>
(gpsType);
// 提取第6
-19个字符（索引5到18，共14个字符）
if
(gLine.length
() >= 20) // 确保有足够长度
{
try
{
std::string dataStr = gLine.substr
(5, 14);
dataField = stod
(dataStr);
}
catch
(...)
{
// 处理数据转换异常
}
}
int gpsC1C = static_cast<int>
(dataField);
// TO
DO: 在这里将提取的卫星数据
(gpsnub, gpsC1C等)存储起来
// 例如：添加到数据结构中，与当前timeStamp关联
//cout << "G" << gpsnub << " " << "C1C=" << gpsC1C << endl;
}
}
// 回退到下一个>行的开头
if
(nextEpochStartPos !=
-1) {
fseek
(pFile, nextEpochStartPos, SEEK_SET);
}
}
// TO
DO: 在这里可以处理当前历元的所有数据（如统计卫星数量等）
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 导航文件读取器：按“G行+7行参数=8行/组”读取RINEX NAV文件
void reader::rinexNAVReader
(vector<NAVRECORDS>& nav_RECORDS) {
char aline[MAX_STRLEN]; // 行缓存（建议MAX_STRLEN定义为256，适配长参数行）
vector<string> groupData; // 存储每组8行数据（1个G行 + 7个参数行）
//
-------------------------- 1. 跳过文件头部（优化原逻辑）
--------------------------
if
(navHeader == true) {
// 读取首行（若头部未处理）
if
(fgets
(aline, MAX_STRLEN, NavpFile) == NULL) {
cerr << "Error: 导航文件为空或无法读取！" << endl;
return;
}
// 循环跳过头部，直到找到"END OF HEADER"或文件结束
while (!feof
(NavpFile) &&
strstr(aline, "END OF HEADER") == NULL) {
if
(fgets
(aline, MAX_STRLEN, NavpFile) == NULL) {
cerr << "Error: 未找到头部结束标记（END OF HEADER）！" << endl;
return;
}
}
navHeader = false; // 标记头部已处理，后续不再重复执行
}
//
-------------------------- 2. 核心：按“G行+7行”分组读取数据
--------------------------
while (true) {
// 2.1 清空上一组缓存，寻找下一个G开头的行（组起始标志）
groupData.clear
();
bool foundGLine = false;
// 循环读取行，直到找到G开头的行或文件结束
while (fgets
(aline, MAX_STRLEN, NavpFile) != NULL) {
string currentLine
(aline);
// 去除行尾换行符（避免解析干扰）
currentLine.erase
(currentLine.find_last_not_of
("\r\n") + 1);
// 跳过空行，判断首字符是否为'G'
if
(!currentLine.empty
() && currentLine[0] == 'G') {
groupData.push_back
(currentLine);
foundGLine = true;
break;
}
}
// 若未找到G行，文件读取结束，退出循环
if
(!foundGLine) {
break;
}
// 2.2 读取该G行后续的7行参数（组成8行/组的完整数据）
bool groupComplete = true;
for
(int i = 0; i < 7; i++) { // 已存1行G行，再读7行
if
(fgets
(aline, MAX_STRLEN, NavpFile) == NULL) {
groupComplete = false;
break; // 文件提前结束，组不完整
}
// 去除行尾换行符，存入缓存
string paramLine
(aline);
paramLine.erase
(paramLine.find_last_not_of
("\r\n") + 1);
groupData.push_back
(paramLine);
}
// 若组不完整（不足8行），丢弃该组并继续下一轮
if
(!groupComplete || groupData.size
() != 8) {
groupData.clear
();
continue;
}
//
-------------------------- 3. 解析当前完整组（8行数据）
--------------------------
NAVRECORDS navrecord_temp;
double prn, Y, M, D, H, Min, Sec, timeStamp, timeSEC;
// 3.1 解析第1行（G开头行：卫星编号+时间+卫星钟参数）
string gLine = groupData[0];
// 卫星编号（G01 → 截取索引1
-2："01"）
string prn_s = gLine.substr
(1, 2);
prn = stod
(prn_s);
// 时间字段（适配你的示例格式：G01 2022 09 03 00 00 00）
Y = stod
(gLine.substr
(4, 4)); // 年：索引4
-7（4位，如2022）
M = stod
(gLine.substr
(9, 2)); // 月：索引9
-10（2位，如09）
D = stod
(gLine.substr
(12, 2)); // 日：索引12
-13（2位，如03）
H = stod
(gLine.substr
(15, 2)); // 时：索引15
-16（2位，如00）
Min = stod
(gLine.substr
(18, 2));// 分：索引18
-19（2位，如00）
Sec = stod
(gLine.substr
(21, 2));// 秒：索引21
-22（2位，如00）
// 卫星钟参数（适配示例中字段位置：秒后紧跟钟差、钟速、钟速变化率）
string SVCB_s = gLine.substr
(23, 19); // 钟差：索引24
-41（如"2.907151356340E
-04"）
string SVCD_s = gLine.substr
(42, 19); // 钟偏：索引43
-60（如"
-6.707523425580E
-12"）
string SVCDR_s = gLine.substr
(61, 19); // 钟偏移：索引62
-79（如"0.000000000000E+00"）
// 科学计数法转换（RINEX的D→标准E）
double SVCB = scientificTo
double
(SVCB_s);
double SVCD = scientificTo
double
(SVCD_s);
double SVCDR = scientificTo
double
(SVCDR_s);
// 3.2 计算GPS时间
和当天总秒数
vector<
double> timeHolder = {Y, M, D, H, Min, Sec};
timeStamp = GPSTime
(timeHolder); // 调用你的GPS时间转换函数
timeSEC = H * 3600 + Min * 60 + Sec; // 当天总秒数（时→秒+分→秒+秒）
// 3.3 解析轨道参数（第2
-7行，共6行，每行4个参数 → 24个，取前24个适配结构体）
vector<
double> navRecHolder;
for
(int i = 1; i <= 6; i++) { // 组内第2
-8行（索引1
-6，共6行参数）
string paramLine = groupData[i];
// 关键修改1：固定跳过前4个字符（而非跳过空格），从第5个字符开始读取（索引4）
size_t startIdx = 4;
// 每行读取4个参数，每个参数固定19个字符宽度
for
(int j = 0; j < 4; j++) {
// 关键判断：确保当前起始索引 + 19字符不超出行长度（避免越界）
if
(startIdx + 19 > paramLine.length
()) {
cerr << "Warning: 参数行长度不足！行索引：" << i
<< "，当前起始位置：" << startIdx << endl;
break;
}
// 提取19个字符的参数（科学计数法字符串）
string p_s = paramLine.substr
(startIdx, 19);
// （可选）调试：打印提取的原始参数字符串，验证是否正确
// cout << "调试：提取参数[" << j << "]：" << p_s << endl;
// 科学计数法转换（D→E），存入临时容器
double p = scientificTo
double
(p_s);
navRecHolder.push_back
(p);
// 关键修改2：下一个参数起始位置 = 当前起始位置 + 19（无需跳过空格，按固定宽度偏移）
startIdx += 19;
}
}
// 校验轨道参数数量（至少24个，否则丢弃该组）
if
(navRecHolder.size
() < 24) {
groupData.clear
();
continue;
}
// 3.4 解析传输时间
和拟合区间（原代码逻辑修正：从第8行提取，而非额外行）
// 说明：你的示例中8行/组已包含所有参数，传输时间
和拟合区间从第8行末尾提取
string lastParamLine = groupData[7]; // 组内第8行（最后一行参数）
size_t fitStartIdx = lastParamLine.length
()
- 19; // 拟合区间：最后19字符
size_t transStartIdx = fitStartIdx
- 19
- 1; // 传输时间：拟合区间前19字符（跳过1个空格）
// 处理边界：确保起始索引合法
if
(transStartIdx < 0 || fitStartIdx < 0) {
navrecord_temp.TransmissionTime = NAN;
navrecord_temp.FitInterval = NAN;
} else {
// 提取传输时间
string TransmissionTime_s = lastParamLine.substr
(transStartIdx, 19);
navrecord_temp.TransmissionTime = scientificTo
double
(TransmissionTime_s);
// 提取拟合区间
string FitInterval_s = lastParamLine.substr
(fitStartIdx, 19);
navrecord_temp.FitInterval = scientificTo
double
(FitInterval_s);
}
//
-------------------------- 4. 数据存入结构体（与原代码完全一致）
--------------------------
navrecord_temp.PRN = prn;
navrecord_temp.SEC = timeSEC;
navrecord_temp.gpstime = timeStamp;
navrecord_temp.SVClockBias = SVCB;
navrecord_temp.SVClockDrift = SVCD;
navrecord_temp.SVClockDriftRate = SVCDR;
// 轨道参数赋值（24个字段，与navRecHolder索引一一对应）
navrecord_temp.IODE = navRecHolder[0];
navrecord_temp.Crs = navRecHolder[1];
navrecord_temp.DeltaN = navRecHolder[2];
navrecord_temp.Mo = navRecHolder[3];
navrecord_temp.Cuc = navRecHolder[4];
navrecord_temp.Eccentricity = navRecHolder[5];
navrecord_temp.Cus = navRecHolder[6];
navrecord_temp.Sqrta = navRecHolder[7];
navrecord_temp.Toe = navRecHolder[8];
navrecord_temp.Cic = navRecHolder[9];
navrecord_temp.OMEGA = navRecHolder[10];
navrecord_temp.Cis = navRecHolder[11];
navrecord_temp.Io = navRecHolder[12];
navrecord_temp.Crc = navRecHolder[13];
navrecord_temp.omega = navRecHolder[14];
navrecord_temp.OMEGA
DOT = navRecHolder[15];
navrecord_temp.I
DOT = navRecHolder[16];
navrecord_temp.L2CodesChannel = navRecHolder[17];
navrecord_temp.GPSWeek = navRecHolder[18];
navrecord_temp.L2PDataFlag = navRecHolder[19];
navrecord_temp.SVAccuracy = navRecHolder[20];
navrecord_temp.SVHealth = navRecHolder[21];
navrecord_temp.TGD = navRecHolder[22];
navrecord_temp.IODC = navRecHolder[23];
//
-------------------------- 5. 存入结果向量
--------------------------
nav_RECORDS.push_back
(navrecord_temp);
// 可选：打印组解析成功提示（便于调试）
std::cout << "Success: 解析第" << nav_RECORDS.size
() << "组数据（PRN:" << prn << "，时间:" <<
(int)Y << "
-" <<
(int)M << "
-" <<
(int)D << " " <<
(int)H << ":" <<
(int)Min << "）" << endl;
}
//
-------------------------- 6. 读取完成统计
--------------------------
std::cout << "\nRINEX NAV文件读取完成！共解析有效数据组：" << nav_RECORDS.size
() << " 组" << endl;
}
3、以下是源.cpp的代码：#include <iostream>
using namespace std;
#include <string>
#include<time.h>
#include <vector>
#include "reader.h" // 包含reader类的头文件
int main
() {
// 定义O文件
和N文件的路径
const char* obsFilename = "F:/学习/1.研究生科研文件夹/2025长安大学暑期培训/3.伪距单点定位原理与实现相关资料/spp数据参考/test/brst2460.22o";
const char* navFilename = "F:/学习/1.研究生科研文件夹/2025长安大学暑期培训/3.伪距单点定位原理与实现相关资料/spp数据参考/test/brdm2460.22p";
// 创建reader对象，传入文件路径
reader rinexReader
(obsFilename, navFilename);
// 声明需要的变量
reader::Header head;
reader::RECORDS obs_RECORDS;
int SENTINEL = 0; // 初始值设为0，根据需要调整
// 调用读取函数
rinexReader.rinexOBSReader
(head, obs_RECORDS, SENTINEL);
vector<reader::NAVRECORDS> navData; // 存储导航数据的向量
rinexReader.rinexNAVReader
(navData); // 调用导航文件读取函数
return 0;
}
本文通过一个具体的示例详细介绍了如何使用strstr函数来查找一个字符串在另一个字符串中首次出现的位置，并展示了如何利用循环来查找所有匹配项的方法。
扫一扫
5万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
