066-松组合备忘之单点定位

一、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。
  注意在rinex中，年份都只有后两位，所以首先进行判断所处具体年份，如果$y$y大于等于80，加1900，否则加2000，这样就求得完整的年份。
  下面还有对月份的判断，如果月份大于2，不做任何处理；否则年份-1，月份+12。至于为什么，按照我的理解是跟公式中的30.6001有关，如果月份小于等于2，由于2月份的缘故就会使平均天数出问题，这时将新一年的两个月放置到前一年作为13月和14月就可以补偿，感兴趣的可以测试一下。
  进行完上述调整，按照计算公式可得到：
$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

  有儒略日$JD$JD，计算公式如下：
$Week = fix(\frac{JD-2444244.5}{7})$Week=fix(7JD−2444244.5​)

$Second = (\frac{JD-2444244.5}{7}-Week)*7*24*3600$Second=(7JD−2444244.5​−Week)∗7∗24∗3600

  GPS标准时间定义为$JD=2444244.5$JD=2444244.5(1980年1月6日0时)

四、单点定位解算

  在此之前要进行一些必要的说明。在O文件读取过程中，对不同的观测值类型要进行编号，在这里就需要提取出L1观测值(伪距)进行单点定位的解算。有多少个历元，就进行多少个历元的计算，得到多少个接收机的位置，在此就一个历元进行讨论。

1、卫星位置计算

  首先要将之前读取的单历元数据与相应的导航电文进行链接，然后进行各个卫星位置的计算。在此不再进行码字介绍如何进行卫星的位置计算，截图显示一下。

  卫星的位置应为卫星发射信号的时刻，而接收机中是接受信号的时刻，所以需要减去信号传播的时间获取卫星发射信号的时刻。因为信号传播时间我们无从知晓，所以首先设置一个大概的传播时间t1，求出卫星大概的发射信号时刻来推算卫星的位置。在此我们设置条件，求得卫星位置之后，可以根据其与测站的近似坐标得到近似站星距，从而得到卫星传播时间t2，比较t1和t2，获取差值，当差值小于设定的阈值时，可以认为得到了较为准确的卫星位置，可以用来计算接收机的位置。
  当然，计算出卫星发射信号t1时刻的位置，还要知道收到卫星信号t2时刻卫星的位置，添加地球自转改正：

  另外，卫星钟差与相对论效应也在此时顺便求取：

$\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$a0​,a1​,a2​从导航电文中获取，分别为t0时刻钟差、钟速(频偏)和漂移项；

$\Delta t_r = -\frac{2 \sqrt{a \mu} } {c^2} e sinE$Δtr​=−c22aμ​​esinE

$a$a：卫星轨道长半径；
$e$e：卫星轨道偏心率；
$E$E：卫星的偏近点角；
$\mu$μ：地心引力常数$39860047×10^7m^3/s^2$39860047×107m3/s2；
$c$c: 光速 $2.99792458×10^8 m/s$2.99792458×108m/s。

2、电离层延迟

  电离层改正根据《GPS理论、算法与应用(第三版)》（许国昌）中的克罗布歇(Klobuchar)模型进行，可能是我的程序的原因，经过这个改正，我的解算结果变差了。手动尴尬！

3、对流层延迟

  对流层改正依据RTKLIB中的萨斯塔莫宁(Saastamoinen)模型进行，如下：

4、平差计算

  平差计算是单点定位最为重要的一部分，稍有不慎很容易出错，特别是不同的定义下得到的结果不同，很容易造成思维混乱，所以进行这一部分的内容一定要自己推导一下完整的公式！

  伪距观测的观测方程

$\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​

其中：
$\tilde \rho_i$ρ~​i​ : 伪距观测值；
$X^i,Y^i,Z^i$Xi,Yi,Zi ： 第i颗卫星坐标；
$X,Y,Z$X,Y,Z ： 测站坐标；
$c$c ： 光速；
$V_{tR}$VtR​ ： 接收机钟差；
$V_{t_i^S}$VtiS​​ ： 卫星钟差；
$(V_{ion})_i$(Vion​)i​ ： 第i颗卫星电离层延迟改正；
$(V_{trop})_i$(Vtrop​)i​ ： 第i颗卫星对流层延迟改正。
  测站近似坐标设为$(X^0,Y^0,Z^0)$(X0,Y0,Z0)，在近似坐标点泰勒展开线性化：

$\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​+ri0​X0−Xi​Vx​+ri0​Y0−Yi​Vy​+ri0​Z0−Zi​Vz​+cVtR​−cVtiS​​+(Vion​)i​+(Vtrop​)i​

其中：

$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​

令：

$\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​=ri0​X0−Xi​mi​=ri0​Y0−Yi​ni​=ri0​Z0−Zi​Li​=ρ~​i​−ri0​+cVtiS​​−(Vion​)i​−(Vtrop​)i​​

  得到误差方程形如：

$v_i = l_i V_x + m_i V_y + n_i V_z + cV_{tR} - L_i$vi​=li​Vx​+mi​Vy​+ni​Vz​+cVtR​−Li​

  卫星数量大于等于4颗时：

$\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}$⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡​v1​v2​⋅⋅⋅vn​​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤​=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡​l1​l2​⋅⋅⋅ln​​m1​m2​⋅⋅⋅mn​​n1​n2​⋅⋅⋅nn​​11⋅⋅⋅1​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤​⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡​Vx​Vy​Vz​cVtR​​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤​−⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡​L1​L2​⋅⋅⋅Ln​​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤​

  也就是这样的形式：

$V = B \hat x-L$V=Bx^−L

那么：

$\hat x = (B^T P B)^{-1}B^TPL$x^=(BTPB)−1BTPL

坐标结果：

$\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​⎦⎥⎥⎥⎥⎤​=⎣⎢⎢⎢⎢⎡​X0​Y0​Z0​​⎦⎥⎥⎥⎥⎤​+⎣⎢⎢⎢⎢⎡​Vx​Vy​Vz​​⎦⎥⎥⎥⎥⎤​

钟差：

$V_{tR} = \frac{cV_{tR}}{c}$VtR​=ccVtR​​

协因数阵：

$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=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡​qxx​qyx​qzx​qtx​​qxy​qyy​qzy​qty​​qxz​qyz​qzz​qtz​​qxt​qyt​qzt​qtt​​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤​

单位权中误差：

$m_0 = \sqrt{\frac{V^T PV}{n-4} }$m0​=n−4VTPV​​

解算精度:

$\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​=m0​qxx​​my​=m0​qyy​​mz​=m0​qzz​​mtR​=m0​qtt​​mpos​=m0​qxx​+qyy​+qzz​​​

精度因子：

$\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​+qzz​​TDOP=qtt​​GDOP=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参考手册