一、MATLAB相关函数介绍
1、 uigetfile
[fname, fpath] = uigetfile('*.*N', 'File N');
功能：打开一个扩展名为“*N”的文件，获取文件名与其绝对路径； File N为弹出的对话框的名称。
2、strcat
filename = strcat(fpath, fname);
功能：将两个字符串连接在一起，这样就获取了N文件的路径。
3、fopen
fid = fopen(filename, 'rt');
功能：打开一个文件，rt表示打开文件用于只读，并且是以文本文件的形式打开，其他方式可参考相关文档；fid表示文件的句柄，或者说是文件的代号或标识，当其值为一个正整数，表示打开成功，若为-1，表示没有打开，因此进行文件读取应该添加是否成功打开文件的判断。
4、fgetl
strow = fgetl(fid);
功能：从文件中读取一行数据并丢弃其中的换行符。如果读取成功，strow容纳了读取到的文本字符串，如果遇到文件末尾的结束标志（EOF），则函数返回-1，即strow值为-1。
5、strfind
EndofHeader = strfind(strow, 'END OF HEADER');
功能：从字符串strow中查找是否含有字符串’END OF HEADER’，如果找到，返回被查找的字符串第一个字符的索引，即字母E的位置在strow中是第一个字符，如果没有查找到，返回一个空数组[]。
6、fclose
fclose(fid);
功能：关闭打开的文件。
7、str2num
X = str2num(strow);
功能：将字符串strow转换为数字。
8、strcmp
TF = strcmp(Type, 'L1')
功能：判断Type中的字符串是否等于’L1’，如果相等，返回1，否则返回0。至于什么时候用strcmp，什么时候用strfind，在编写NO文件读取的时候会有更深的体会。
二、N文件和O文件的读取
关于N文件和O文件的读取，最重要的是搞清楚文件中各个位置所代表的字符的含义，关于这方面的内容我在之前的文档中已经做过简要的介绍，还是建议自己对照N文件和O文件进行查验，加深自己的理解。感兴趣的参考
052RINEX中N文件示例说明
053RINEX中O文件示例说明
在此推荐一下vim编辑器，用其查看此类文件的比例之处就是可以显示索引。点此进入vim网站下载
三、时间系统的换算
N文件和O文件中的时间即为格里高利历的形式，即年月日时分秒的形式，但是文件中的时间系统皆属于GPS时间，所以在将该时间换算为GPS周及GPS秒的形式时，不需要考虑跳秒的问题。换算时引入儒略日，即先换算到儒略日，然后换算到GPS形式，下面简单罗列一下计算公式。
1、从格里高利到儒略日
有年月日时分秒：
y
,
m
,
d
,
h
,
m
,
s
y,m,d,h,m,s
y,m,d,h,m,s。
注意在rinex中，年份都只有后两位，所以首先进行判断所处具体年份，如果
y
y
y大于等于80，加1900，否则加2000，这样就求得完整的年份。
下面还有对月份的判断，如果月份大于2，不做任何处理；否则年份-1，月份+12。至于为什么，按照我的理解是跟公式中的30.6001有关，如果月份小于等于2，由于2月份的缘故就会使平均天数出问题，这时将新一年的两个月放置到前一年作为13月和14月就可以补偿，感兴趣的可以测试一下。
进行完上述调整，按照计算公式可得到：
J
D
=
f
i
x
(
365.25
∗
y
)
+
f
i
x
(
30.6001
∗
(
m
+
1
)
)
+
d
+
h
+
m
/
60
+
s
/
3600
24
+
1720981.5
JD = fix(365.25*y) + fix(30.6001*(m+1))+d+\frac{h+m/60+s/3600}{24}+1720981.5
JD=fix(365.25∗y)+fix(30.6001∗(m+1))+d+24h+m/60+s/3600+1720981.5
2、从儒略日到GPST
有儒略日
J
D
JD
JD，计算公式如下：
W
e
e
k
=
f
i
x
(
J
D
−
2444244.5
7
)
Week = fix(\frac{JD-2444244.5}{7})
Week=fix(7JD−2444244.5)
S e c o n d = ( J D − 2444244.5 7 − W e e k ) ∗ 7 ∗ 24 ∗ 3600 Second = (\frac{JD-2444244.5}{7}-Week)*7*24*3600 Second=(7JD−2444244.5−Week)∗7∗24∗3600
GPS标准时间定义为 J D = 2444244.5 JD=2444244.5 JD=2444244.5(1980年1月6日0时)
四、单点定位解算
在此之前要进行一些必要的说明。在O文件读取过程中，对不同的观测值类型要进行编号，在这里就需要提取出L1观测值(伪距)进行单点定位的解算。有多少个历元，就进行多少个历元的计算，得到多少个接收机的位置，在此就一个历元进行讨论。
1、卫星位置计算
首先要将之前读取的单历元数据与相应的导航电文进行链接，然后进行各个卫星位置的计算。在此不再进行码字介绍如何进行卫星的位置计算，截图显示一下。
卫星的位置应为卫星发射信号的时刻，而接收机中是接受信号的时刻，所以需要减去信号传播的时间获取卫星发射信号的时刻。因为信号传播时间我们无从知晓，所以首先设置一个大概的传播时间t1，求出卫星大概的发射信号时刻来推算卫星的位置。在此我们设置条件，求得卫星位置之后，可以根据其与测站的近似坐标得到近似站星距，从而得到卫星传播时间t2，比较t1和t2，获取差值，当差值小于设定的阈值时，可以认为得到了较为准确的卫星位置，可以用来计算接收机的位置。
当然，计算出卫星发射信号t1时刻的位置，还要知道收到卫星信号t2时刻卫星的位置，添加地球自转改正：
另外，卫星钟差与相对论效应也在此时顺便求取：
Δ t S V = a 0 + a 1 ( t − t 0 ) + a 2 ( t − t 0 ) 2 + Δ t r \Delta t_{SV} = a_0 + a_1(t-t_0) + a_2(t-t_0)^2+\Delta t_r ΔtSV=a0+a1(t−t0)+a2(t−t0)2+Δtr
其中， a 0 , a 1 , a 2 a_0, a_1, a_2 a0,a1,a2从导航电文中获取，分别为t0时刻钟差、钟速(频偏)和漂移项；
Δ t r = − 2 a μ c 2 e s i n E \Delta t_r = -\frac{2 \sqrt{a \mu} } {c^2} e sinE Δtr=−c22aμesinE
a
a
a：卫星轨道长半径；
e
e
e：卫星轨道偏心率；
E
E
E：卫星的偏近点角；
μ
\mu
μ：地心引力常数
39860047
×
1
0
7
m
3
/
s
2
39860047×10^7m^3/s^2
39860047×107m3/s2；
c
c
c: 光速
2.99792458
×
1
0
8
m
/
s
2.99792458×10^8 m/s
2.99792458×108m/s。
2、电离层延迟
电离层改正根据《GPS理论、算法与应用(第三版)》（许国昌）中的克罗布歇(Klobuchar)模型进行，可能是我的程序的原因，经过这个改正，我的解算结果变差了。手动尴尬！
3、对流层延迟
对流层改正依据RTKLIB中的萨斯塔莫宁(Saastamoinen)模型进行，如下：
4、平差计算
平差计算是单点定位最为重要的一部分，稍有不慎很容易出错，特别是不同的定义下得到的结果不同，很容易造成思维混乱，所以进行这一部分的内容一定要自己推导一下完整的公式！
伪距观测的观测方程
ρ ~ i = ( X i − X ) 2 + ( Y i − Y ) 2 + ( Z i − Z ) 2 − c V t R + c V t i S − ( V i o n ) i − ( V t r o p ) i \tilde \rho_i = \sqrt{(X^i-X)^2+(Y^i-Y)^2+(Z^i-Z)^2}-cV_{tR}+cV_{t_i^S}-(V_{ion})_i-(V_{trop})_i ρ~i=(Xi−X)2+(Yi−Y)2+(Zi−Z)2−cVtR+cVtiS−(Vion)i−(Vtrop)i
其中：
ρ
~
i
\tilde \rho_i
ρ~i : 伪距观测值；
X
i
,
Y
i
,
Z
i
X^i,Y^i,Z^i
Xi,Yi,Zi ： 第i颗卫星坐标；
X
,
Y
,
Z
X,Y,Z
X,Y,Z ： 测站坐标；
c
c
c ： 光速；
V
t
R
V_{tR}
VtR ： 接收机钟差；
V
t
i
S
V_{t_i^S}
VtiS ： 卫星钟差；
(
V
i
o
n
)
i
(V_{ion})_i
(Vion)i ： 第i颗卫星电离层延迟改正；
(
V
t
r
o
p
)
i
(V_{trop})_i
(Vtrop)i ： 第i颗卫星对流层延迟改正。
测站近似坐标设为
(
X
0
,
Y
0
,
Z
0
)
(X^0,Y^0,Z^0)
(X0,Y0,Z0)，在近似坐标点泰勒展开线性化：
ρ ~ i = r i 0 + X 0 − X i r i 0 V x + Y 0 − Y i r i 0 V y + Z 0 − Z i r i 0 V z + c V t R − c V t i S + ( V i o n ) i + ( V t r o p ) i \tilde \rho_i = r_i^0 + \frac{X^0-X^i}{r_i^0} V_x + \frac{Y^0-Y^i}{r_i^0} V_y + \frac{Z^0-Z^i}{r_i^0} V_z +cV_{tR}-cV_{t_i^S}+(V_{ion})_i+(V_{trop})_i ρ~i=ri0+ri0X0−XiVx+ri0Y0−YiVy+ri0Z0−ZiVz+cVtR−cVtiS+(Vion)i+(Vtrop)i
其中：
r i 0 = ( X i − X 0 ) 2 + ( Y i − Y 0 ) 2 + ( Z i − Z 0 ) 2 r_i^0 = \sqrt{(X^i-X^0)^2+(Y^i-Y^0)^2+(Z^i-Z^0)^2} ri0=(Xi−X0)2+(Yi−Y0)2+(Zi−Z0)2
令：
l i = X 0 − X i r i 0 m i = Y 0 − Y i r i 0 n i = Z 0 − Z i r i 0 L i = ρ ~ i − r i 0 + c V t i S − ( V i o n ) i − ( V t r o p ) i \begin{matrix} l_i = \frac{X^0-X^i}{r_i^0} \\ \\ m_i = \frac{Y^0-Y^i}{r_i^0} \\ \\ n_i = \frac{Z^0-Z^i}{r_i^0} \\ L_i = \tilde \rho_i - r_i^0 + cV_{t_i^S}-(V_{ion})_i - (V_{trop})_i \end{matrix} li=ri0X0−Ximi=ri0Y0−Yini=ri0Z0−ZiLi=ρ~i−ri0+cVtiS−(Vion)i−(Vtrop)i
得到误差方程形如：
v i = l i V x + m i V y + n i V z + c V t R − L i v_i = l_i V_x + m_i V_y + n_i V_z + cV_{tR} - L_i vi=liVx+miVy+niVz+cVtR−Li
卫星数量大于等于4颗时：
[ v 1 v 2 ⋅ ⋅ ⋅ v n ] = [ l 1 m 1 n 1 1 l 2 m 2 n 2 1 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ l n m n n n 1 ] [ V x V y V z c V t R ] − [ L 1 L 2 ⋅ ⋅ ⋅ L n ] \begin{bmatrix} v_1 \\ \\ v_2 \\ \\ ··· \\ \\ v_n\\ \end{bmatrix} = \begin{bmatrix} l_1 & m_1 & n_1 & 1 \\ \\ l_2 & m_2 & n_2 & 1 \\ \\ ··· & ··· & ··· & ···\\ \\ l_n & m_n & n_n & 1 \\ \end{bmatrix} \begin{bmatrix} V_x \\ \\ V_y \\ \\ V_z \\ \\ cV_{tR}\\ \end{bmatrix} - \begin{bmatrix} L_1 \\ \\ L_2 \\ \\ ··· \\ \\ L_n\\ \end{bmatrix} ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡v1v2⋅⋅⋅vn⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡l1l2⋅⋅⋅lnm1m2⋅⋅⋅mnn1n2⋅⋅⋅nn11⋅⋅⋅1⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡VxVyVzcVtR⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤−⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡L1L2⋅⋅⋅Ln⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤
也就是这样的形式：
V = B x ^ − L V = B \hat x-L V=Bx^−L
那么：
x ^ = ( B T P B ) − 1 B T P L \hat x = (B^T P B)^{-1}B^TPL x^=(BTPB)−1BTPL
坐标结果：
[ X Y Z ] = [ X 0 Y 0 Z 0 ] + [ V x V y V z ] \begin{bmatrix} X\\ \\ Y\\ \\ Z\\ \end{bmatrix} = \begin{bmatrix} X_0\\ \\ Y_0\\ \\ Z_0\\ \end{bmatrix} + \begin{bmatrix} V_x\\ \\ V_y\\ \\ V_z\\ \end{bmatrix} ⎣⎢⎢⎢⎢⎡XYZ⎦⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎡X0Y0Z0⎦⎥⎥⎥⎥⎤+⎣⎢⎢⎢⎢⎡VxVyVz⎦⎥⎥⎥⎥⎤
钟差：
V t R = c V t R c V_{tR} = \frac{cV_{tR}}{c} VtR=ccVtR
协因数阵：
Q x ^ x ^ = ( B T P B ) − 1 = [ q x x q x y q x z q x t q y x q y y q y z q y t q z x q z y q z z q z t q t x q t y q t z q t t ] Q_{\hat x \hat x} = (B^T P B)^{-1} = \begin{bmatrix} q_{xx} & q_{xy} & q_{xz} & q_{xt} \\ \\ q_{yx} & q_{yy} & q_{yz} & q_{yt} \\ \\ q_{zx} & q_{zy} & q_{zz} & q_{zt} \\ \\ q_{tx} & q_{ty} & q_{tz} & q_{tt} \\ \end{bmatrix} Qx^x^=(BTPB)−1=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡qxxqyxqzxqtxqxyqyyqzyqtyqxzqyzqzzqtzqxtqytqztqtt⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤
单位权中误差：
m 0 = V T P V n − 4 m_0 = \sqrt{\frac{V^T PV}{n-4} } m0=n−4VTPV
解算精度:
m x = m 0 q x x m y = m 0 q y y m z = m 0 q z z m t R = m 0 q t t m p o s = m 0 q x x + q y y + q z z \begin{matrix} m_x = m_0 \sqrt{q_{xx}}\\ \\ m_y = m_0 \sqrt{q_{yy}}\\ \\ m_z = m_0 \sqrt{q_{zz}}\\ \\ m_{tR} = m_0 \sqrt{q_{tt}}\\ \\ m_{pos} = m_0 \sqrt{q_{xx}+q_{yy}+q_{zz}} \end{matrix} mx=m0qxxmy=m0qyymz=m0qzzmtR=m0qttmpos=m0qxx+qyy+qzz
精度因子：
P D O P = q x x + q y y + q z z T D O P = q t t G D O P = q x x + q y y + q z z + q t t \begin{matrix} PDOP = \sqrt{q_{xx}+q_{yy}+q_{zz}} \\ \\ TDOP = \sqrt{q_{tt}} \\ \\ GDOP = \sqrt{q_{xx}+q_{yy}+q_{zz}+q_{tt}} \end{matrix} PDOP=qxx+qyy+qzzTDOP=qttGDOP=qxx+qyy+qzz+qtt
至此，解算完成。
参考：
[1] 李征航,黄劲松. GPS测量与数据处理(第三版)[M]. 武汉:武汉大学出版社,2016.
[2] 许国昌,许艳. GPS理论、算法与应用(第三版)[M]. 北京:科学出版社,2017.
[3] 许承东,李怀建,张鹏飞等. GNSS数学仿真原理及系统实现[M]. 北京:中国宇航出版社,2014.
[4] 周忠谟,易杰军,周琪. GPS卫星测量原理与应用(修订版)[M]. 北京:测绘出版社,1997.
[5] 张小红. GPS单点定位PPT
[6] RTKLIB参考手册
208
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
