[学习]RTKLib详解：ephemeris.c与rinex.c

RTKLib详解：ephemeris.c与rinex.c

本文是 RTKLlib详解 系列文章的一篇，目前该系列文章还在持续总结写作中，以发表的如下，有兴趣的可以翻阅。

[学习] RTKlib详解：功能、工具与源码结构解析
[学习]RTKLib详解：pntpos.c与postpos.c
[学习]RTKLib详解：rtkcmn.c与rtkpos.c
[学习]RTKLib详解：ppp.c与ppp_ar.c
[学习]RTKLib详解：ephemeris.c与rinex.c

---

PART A: ephemeris.c

---

一、代码整体作用与工作流程分析

1.1 整体作用

该代码实现卫星星历与钟差计算的核心功能，支持多种卫星导航系统（GPS/Galileo/QZSS/北斗/GLONASS/SBAS），提供从广播星历、精密星历到SSR（State Space Representation）校正的多级精度支持，该代码通过模块化设计实现了多系统兼容，采用迭代法保证计算精度，结合误差模型和校正机制满足不同应用场景需求。主要功能包括：

• 卫星位置与速度计算（ECEF坐标系）
• 卫星钟差与钟漂计算
• 星历数据选择与有效性验证
• 坐标系旋转与相对论效应修正
• 误差方差建模（URA/SSR）

1.2 工作流程

---

二、核心函数说明

2.1 alm2pos (Almanac to Position)

功能: 通过卫星历书计算概略位置
输入:

• time: 当前时间（GPST）
• alm: 阿尔曼历数据（alm_t结构体）
  输出:
• rs: 卫星位置 [x,y,z] (m)
• dts: 钟差 (s)
  数学原理:

1. 计算平近点角 M = M0 + √(μ/A³) * tk
2. 牛顿迭代法解开普勒方程 E = M + e·sinE
3. 真近点角转换: ν = atan2(√(1-e²)sinE, cosE - e)
4. 极坐标到直角坐标转换: r = A(1 - e·cosE)

2.2 eph2clk (Ephemeris to Clock)

功能: 广播星历钟差计算
输入:

• time: 卫星时钟时间（GPST）
• eph: 广播星历数据（eph_t结构体）
  输出:
• 返回钟差 Δt = f0 + f1·t + f2·t²
  特性: 双次迭代消除时间依赖性

2.3 eph2pos (Ephemeris to Position)

功能: 广播星历卫星位置计算
数学流程:

关键公式:

• 开普勒方程：E = M + e·sinE （迭代收敛精度1e-14）
• 相对论修正项：Δt_rel = -2√(μA)/c² * e·sinE
• 坐标旋转矩阵：涉及升交点赤经Ω的时变特性

---

2.4 geph2pos (GLONASS Ephemeris)

特殊性:

• 使用数值积分法解算轨道微分方程
• 包含地球扁率J2和自转角速度Ω的摄动项
• 微分方程形式：ẍ = [c+Ω²]x + 2Ωẏ + ax

2.5 satpos_ssr (SSR校正)

校正流程:

1. 获取轨道偏差Δr [径向/切向/横向]
2. 获取钟差校正Δt = dclk[0] + dclk[1]·t + dclk[2]·t²
3. 天线相位中心修正
4. 方差计算采用SSR URA模型

三、核心算法数学推导

3.1 开普勒方程求解

标准形式：
$$M=E-e\sin E$$
牛顿迭代公式：
$$E_{n+1}=E_n-\frac{E_n-e\sin E_n-M}{1-e\cos E_n}$$
代码实现中设置最大迭代次数 $30$ 次，相对容差 $10^{-14}$。

3.2 坐标系转换

从轨道坐标系到ECEF的转换矩阵：
$$\left[\begin{array}{ccc}\cos u& -\sin u& 0\\ \sin u& \cos u& 0\\ 0& 0& 1\end{array}\right]$$
结合升交点赤经 $ \Omega $ 的时变特性：
$$\Omega ={\Omega }_0+(\dot{\Omega }-{\omega }_e)t-{\omega }_e{t}_{\text{os}}$$

3.3 SSR轨道校正

采用径向-切向-横向（Radial-Across-Along）分解：
$$\Delta \vec{r}={\vec{e}}_r\cdot \Delta r_{\text{radial}}+{\vec{e}}_a\cdot \Delta r_{\text{along}}+{\vec{e}}_c\cdot \Delta r_{\text{cross}}$$
其中 { e ⃗ r , e ⃗ a , e ⃗ c } \{\vec{e}_r, \vec{e}_a, \vec{e}_c\} {e r​,e a​,e c​} 构成正交基矢量。

四、完整函数调用关系

五、误差模型

5.1 URA转换模型

var_uraeph:
$$\text{ura\_value}=\left[2.4,3.4,\dots ,6144.0\right]$$
$$\text{var}=\{\begin{array}{ll}(6144.0{)}^2& \text{if }\text{ura}<0\text{ or }\text{ura}>15\\ (\text{ura\_value}[\text{ura}]{)}^2& \text{otherwise}\end{array}$$

var_urassr:
$$\text{std}=\{\begin{array}{ll}(3^{(\lfloor \text{ura}/8\rfloor )}\cdot (1+(\text{ura}\mathrm{mod}7)/4)-1)\times 10^{-3}& \text{if }0<\text{ura}<63\\ (5.4665{)}^2& \text{if }\text{ura}\ge 63\\ (0.15{)}^2& \text{if }\text{ura}\le 0\end{array}$$

5.2 SSR校正限幅

• 最大轨道修正：$10\text{m}$
• 最大钟差修正：$c\cdot 10^{-6}\text{m}$
• 最大龄期：轨道 $90\text{s}$ / 钟差 $10\text{s}$

5.3 相对论修正项

$$\Delta t_{\text{rel}}=-\frac{2}{c^2}\sqrt{\mu A}\cdot e\cdot \sin E$$

5.4 微分方程形式（GLONASS）

$$\ddot{x}=[c+{\Omega }^2]x+2\Omega \dot{y}+a_x$$

5.5 SSR钟差修正

$$\Delta t_{\text{corr}}=\Delta t_{\text{clk}}+\frac{\Delta {\vec{r}}_{\text{radial}}}{c}$$

---

六、系统特性处理

系统	引力常数 $ \mu $	地球自转角速度 $ \Omega $	特殊处理
GPS	$3.9860050\times 10^{14}$	$7.292115\times 10^{-5}$	标准开普勒解算
GLONASS	$3.9860044\times 10^{14}$	$7.292115\times 10^{-5}$	数值积分法
Galileo	$3.986004418\times 10^{14}$	$7.2921151467\times 10^{-5}$	精密摄动模型
北斗GEO卫星	$3.986004418\times 10^{14}$	$7.292115\times 10^{-5}$	五度倾角坐标变换

---

PART B：rinex.c

---

一、整体作用与工作流程分析

该代码实现了一个RINEX（Receiver Independent Exchange Format）文件处理模块，主要用于GNSS（全球导航卫星系统）数据的读取、解析和输出，该代码通过模块化设计实现了复杂多系统的GNSS数据处理，其核心在于通过统一接口处理不同卫星系统的观测数据与导航数据，同时保持RINEX标准的兼容性。。其核心功能包括：

• RINEX文件头解析与生成
• 观测数据与导航电文的读取
• 卫星星历解码
• 相位偏移处理
• 多系统（GPS/GLO/GAL/QZS/SBS/CMP）支持
• 数据格式转换与输出

工作流程：

---

二、函数详细说明

2.1 init_sta

功能：初始化测站参数结构体
输入参数：

• sta_t *sta：测站参数结构体指针
  输出参数：
• sta：初始化后的测站参数

2.2 sat2code

功能：将卫星编号转换为RINEX格式的卫星代码
输入参数：

• int sat：卫星编号（含系统标识）
• char *code：输出缓冲区
  输出参数：
• code：格式如"G01" (GPS), “R24” (GLONASS)
  数学原理：
  通过satsys()函数分离系统标识符（SYS_GPS/SYS_GLO等）和PRN编号，使用ASCII字符拼接。

2.3 adjweek

功能：调整时间以处理周跳问题
输入参数：

• gtime_t t：当前时间
• gtime_t t0：参考时间
  输出参数：
• gtime_t：调整后的时间
  数学原理：
  $$\begin{gathered} \text{if }\Delta t<-302400s\Rightarrow t=t+604800s \\ \text{if }\Delta t>302400s\Rightarrow t=t-604800s \end{gathered}$$
  其中$\Delta t=\text{timediff}(t,t_0)$，确保时间连续性。

2.4 decode_eph

功能：解码GPS/北斗星历数据
输入参数：

• double ver：RINEX版本号
• int sat：卫星编号
• gtime_t toc：时钟参考时间
• double *data：原始星历数据数组
  输出参数：
• eph_t *eph：填充后的星历结构体
  关键字段处理：
• 计算卫星健康状态（svh）、精度（sva）
• 转换时间参数（toe, ttr）

2.5 uravalue & uraindex

功能：URA值与索引转换
数学关系：
$$\text{ura\_value}=\{\begin{array}{ll}\text{ura\_eph}[i]& \text{if }0\le i<15\\ 32767& \text{otherwise}\end{array}$$
$$\text{ura\_index}(x)=\min \left\{i\mid \text{ura\_eph}[i]\ge x\right\}$$

---

三、函数调用关系图

---

四、关键算法解析

4.1 星历解码（decode_eph）

北斗星历解码示例：

eph->iode = (int)data[3]; // AODE
eph->toes = data[11];     // Toe (s)
eph->week = (int)data[21]; 
eph->toe = bdt2gpst(bdt2time(eph->week, data[11])); 

将北斗时间（BDT）转换为GPS时间，并调整周跳。

4.2 相位偏移校正（set_index）

通过配置选项（如-CL1C=1.5）设置相位偏移：

ind->shift[i] = shift; // 以cycle为单位

输出时应用偏移：

val[i] = str2num(...) + ind->shift[i];

4.3 多系统信号优先级

ind->pri[i] = getcodepri(sys, ind->code[i], opt);

根据系统（sys）和信号类型（code）确定观测值优先级，用于选择主观测值。

---

五、输出流程

---

研究学习不易，点赞易。
工作生活不易，收藏易，点收藏不迷茫 ：）

---