GPS广播星历计算卫星位置和速度

“GPS广播星历计算卫星位置和速度”
及“GPS伪距定位”计算试验
1．试验内容及上交成果
1.1试验内容
应用C语言按预定格式（函数、输入输出变量之名称、类型）编写“GPS广播星历计算卫星位置和速度”函数SatPos_Vel()、“GPS伪距定位”函数Positioning()。

将此两个函数组成文件F2.cpp，并包含于文件GPS_Positioning.cpp中。

编译、连接并运行文件GPS_Positioning.cpp，逐一时刻读取广播星历（Ephemeris.dat）、观测时间及伪距、卫星号（Observation.dat）信息，计算WGS84坐标系中观测时刻相应的卫星位置、速度以及载体位置，结果保存于文件Position.dat中。

1.2上交成果
磁盘文件F2.cpp、Position.dat，并存于“学号作者中文姓名”目录中。

2．函数说明
2.1星历文件读取函数
void EFileReading(Efile)
功能：读取星历文件，给星历数据结构体Ephemeris赋值。

输入变量：EFile字符串，文件名。

2.2观测数据读取函数
int ObsReading(fp_Obs,Time,Rho,Mark)
功能：从文件Observation.dat中读取某一时刻的伪距、卫星号。

读取成功函数值返回“1”，失败返回“－1”（读错，或至文件尾）。

输入变量：fp_Obs文件指针；
输出变量：Time double，时间（秒）；
Rho double[12]，伪距（米）；
Mark int[12]，卫星号，“－1”表示此通道无卫星、无伪距。

2.3最小二乘估计函数
int LeastSquareEstimation(Y,A,P,m,n,X)
功能：最小二乘方法求解观测方程Y=AX+ε，其中观测值方差阵的逆阵为P（也称为权阵），得未知参数X。

成功返回“1”，失败返回“－1”（亏秩）。

输入变量：Y double[m]，观测方程自由项（米）；
A double[m×n]，系数阵（无量纲），按第1行第1、2……n元素，
第2行……顺序存放；
P double[m]，权矩阵对角线元素，0或1（无量纲）；
m int，观测值个数；
n int，未知参数个数；
输出变量：X double[n]，未知参数（米）。

原理：X=(A T PA)-1A T PY
例：Y＝[1,2,3,6],A=[1,0,0,0,1,0,0,0,1,1,1,1],P=[1,1,1,1],m=4,n=3。

调用LeastSquareEstimation(Y,A,P,m,n,X)，求得X=[1,2,3]。

2.4计算卫星位置、速度函数
int SatPos_Vel(Num,Time,SatPos,SatVel)
功能：计算WGS84坐标系中卫星位置、速度。

成功返回“1”，失败返回“－1”（无星历）。

输入变量：Num int，卫星号；
Time double，时间（秒）。

（隐含引用全局变量Ephemeris[32]，星历数据结构体）输出变量：SatPos double[3]，位置（米，WGS84）；
SatVel double[3]，速度（米/秒，WGS84）。

2.5伪距测量定位
int Positioning(Time,Rho,Mark,Pos,dt)
功能：由单一时刻伪距计算载体位置，并保存接收机钟面时、卫星位置、速度及载体位置于定位结果文件Position.dat中。

成功返回“1”，失败返回“－1”（观测
值数量不够，亏秩）。

输入变量：Time double，伪距观测时刻接收机钟面时（秒）；
Rho double[12]，伪距（米）；
Mark int[12]，卫星号，“－1”表示此通道无信号；
输出变量：Pos double[3]，载体位置（米，WGS84）；
dt double，钟差（秒）。

3．全局变量说明
LightVelocity字符常量，光速（米/秒），299792458；
PI字符常量，圆周率，3.141592653589793;
MU字符常量，地球引力常数（米3秒－2），3.986005e14；
OMEGAE字符常量，地球自转角速度（弧度·秒－1），7.292115e-5；
fp_Obs文件指针，观测数据文件；
fp_Pos文件指针，位置结果文件。

Structure Eph_Str
{
Mark整型，标识，“1”有，“－1”无
wn整型，星期数
toc卫星钟差改正参考时刻（秒）
af2卫星钟差改正系数（秒/秒2）
af1（秒/秒）
af0（秒）
aode数据龄期
deltan平均角速度改正（半周/秒）
m0参考时刻平近点角（半周）
e偏心率；
roota a1/2（米1/2）
toe星历参考时刻（秒）
cic（弧度）
crc（米）
cis（弧度）
crs（米）
cuc（弧度）
cus（弧度）
omega0准升交点赤经（半周）
omega近地点角距（半周）
i0轨道倾角（半周）
omegadot升交点赤经变化率（半周/秒）
idot轨道倾角变化率（半周/秒）
}Ephemeris[32]星历数据结构体。

注：数据类型除注明整型外，其余为双精度型。

4．文件说明
4.1观测数据文件Observation.dat
Time,n伪距观测时间，卫星颗数
Rho_1,n_1伪距（米），卫星号
……
Rho_n,n_n伪距（米），卫星号
Time,n伪距观测时间，卫星颗数
……
4.2卫星星历数据文件Ephemeris.dat
svprn整数，卫星号
wn整数，星期数
……
idot轨道倾角变化率（半周/秒）
（详细说明请参见数据结构体Structure Eph_Str Ephemeris[32]）svprn整数，卫星号
……
4.3位置计算结果文件Position.dat
Time伪距观测时刻（秒）
n_1,x_1,y_1,z_1,xdot_1,ydot_1,zdot_1卫星号，卫星坐标（米，WGS84），速度（米/秒）……
n_n,x_n,y_n,z_n,xdot_n,ydot_n,zdot_n卫星号，卫星坐标（米，WGS84），速度（米/秒）x,y,z载体坐标（米，WGS84）
Time伪距观测时刻（秒）
……
4.4GPS_Positioning.cpp
作为主函数，调用Parameter.h、Function.h、F1.cpp、F2.cpp。

逐一时刻读取广播星历、伪距等信息，计算卫星位置、速度及载体位置。

4.5Parameter.h
字符常量定义、全局变量说明。

4.6Function.h
函数说明。

4.7F1.cpp
包含辅助函数如下；
广播星历读取函数EFileReading()；
观测数据（时间、卫星个数、伪距、卫星号）读取函数ObsReading()；
最小二乘法求解函数LeastSquareEstimation()；
矩阵相乘计算函数MatrixMultiply()；
高斯选主元消去法矩阵求逆函数MatrixInverse()；
三维向量长度计算函数VectorLength()。

4.8F2.cpp
包含函数如下：
卫星位置、速度计算函数SatPos_Vel()；
伪距定位函数Positioning()。

5．首历元观测时刻参考答案（Position.dat）
193645.0
29-15567798.7945675.621321624.5249.3-2743.6303.3 014665363.69394841.820224091.6-2310.31043.61194.9 20-10127606.021492092.3-10899994.3-90.1-1466.0-2807.9 21-1578334.526219741.9-284477.9-409.9-14.93151.8 28-20375527.015214719.77639734.1-975.7140.2-2976.3 24-2456690.116712706.920740220.7-2203.1-1370.0828.0 5-24924627.93359562.99069955.6937.9-673.82850.3
1311881411.212148251.420386065.1-509.32539.5-1217.6 -2196881.65177555.22998042.2。