当前位置：首页 > news >正文

卫星导航算法程序设计——单点定位测速（SPV）

news 来源：原创 2025/5/11 11:08:59

本文使用C++语言独立编写SPV软件，可使用二进制文件数据流和实时数据流进行SPV解算。相关成果可通过资源免费下载。

一、程序总体设计

图1展示了整个个单点定位测速系统的架构。整个流程可以分为几个主要模块：

（1）底层运算模块

这一模块包括三个子模块：时间转换模块、坐标转换模块和矩阵运算模块。这些模块为上述数据处理和SPV解算提供底层的运算支持。

（2）数据获取与处理

数据解码模块从网络端口或者文件流中获取原始二进制数据，并将其解码为星历数据、观测数据和电离层参数。

粗差探测用于对获取的数据进行初步的错误检测，剔除粗差。

（3）SPV 解算

这一模块在流程中起到关键作用，它基于 “数据获取与处理” 模块提供的数据，进行信号发射时刻的卫星位置和速度的计算。然后计算地球自转改正、电离层延迟、对流层延迟、相对论改正，而后进行最小二乘单点定位，多普勒单点测速，最后对单点定位的结果进行精度评定。

图1 软件架构

二、算法设计

2.1 时间板块

图2 时间转换

2.2 坐标转换板块

图3 坐标转换

2.3 解码板块

（1）Binary数据头结构

表2.1 Binary数据头结构

1	AA
2	44
3	12
4	Header Length
5	Message ID
8	Message Length
12	Week
13	ms

其中，AA,44,12为同步符，标志着消息的开始，解码过程中匹配到同步符之后才开始数据的解码。注意第13个字节处的GPSsec解码出来的单位为ms，需要转换成s

（2）NonAtel CRC校验程序

数据通信领域中最常用的一种差错校验码，用于验证数据在通信过程中是否存在错误，可靠性很高。

CRC校验原理

被校验的原数据转换成二进制序列，假设共K位

以一定规则产生一个新的二进制序列，假设共R位，即CRC校验码

把新的二进制序列附加在原数据二进制序列后，共K+R位，发送出去

◼ 接收端接收数据后，把原数据的K位二进制序列按相同规则产生一个R位二进制序列与附加的R位二进制序列进行比较，相同表示传递的数据没有问题，不相同表示传递的数据出现错误与误差。

对应函数：

#define POLYCRC32 0xEDB88320u /* CRC32 polynomial */

unsigned int crc32(const unsigned char *buff, int len)

{

　int i,j；

　unsigned int crc=0；

　for (i=0；i<len；i++)

　{

　crc^=buff[i]；

　for (j=0；j<8；j++)

　{

　　　if (crc&1) crc=(crc>>1)^P0LYCRC32；

　　　else crc>>=1；

　　}

　　return crc；

}

（3）原始观测值数据解码

表2.2 原始观测值数据结构

Fielb　　　　　　　　　　　 Format Binary Binary 　　　　Type Bytes Offset
1	RANGE header	Log header.See Messages on page 30 formore information.		H	0
2	# obs	Number of observations with information to follow 1	Ulong	4	H
3	PRN/slot	Satellite PRN number of range measurement Refer to PRN Numbers on page 49	Ushort	2	H+4
4	glofreq	(GLONASS Frequency +7)(see GLONASS Slot and Frequency Numbers section ofthis manual)	Ushort	2	H+6
5	psr	Pseudorange measurement(m)	Double	8	H+8
6	Psr a	Pseudorange measurement standard deviation(m)	Float	4	H+16
7	adr	Carrierphase,in cyces (accumulated Dopplerrange)	Double	8	H+20
8	adr a	Estimated carrler phase standard deviation(cycles)	Float	4	H+28
9	dopp	Instantaneous carrler Doppler frequency (Hz)	Float	4	H+32
10	C/No	Carrler to noisedensity ratio C/No=10[log10(S/N₀)](dB-Hz)	Float	4	H+36
11	locktime	Number of seconds of continuous tracking(no cycle slipplng)	Float	4	H+40
12	ch-tr- status	Tracking status(see Table 149:Channel Tracking Status on the next pageand the example in Figure 15:Channel Tracking Example on the next page)	ulong	4	H+44
13	Next PRN offset =H+4+(#obs×44)
variable	XXXX	32-bit CRC(ASCII and Binary only)	ulong	4
variable	[CR][LF]	Sentence terminator(ASCII only)	-	-	-

原始观测值解码，就是 OEM 中的 Range 模块，解码内容其实就是对应 RINEX 文件中的“.o”文件中的内容。Range 的数据标识 Message ID 为43。Range 数据内容包括表2.2中的内容。其中尤其需要注意的是，由于载波相位的观测值单位是周，这里最好直接转成m单位，这样方面后续使用；另外，最终的相位测量值需要乘以“-”号，DOPP 测量值也是。当解码出上面的数据后，最后的信号跟踪状态数据“ch-tr-status”中就有卫星系统、信号频率，信号跟踪状态等信息。根据信号的标识，可以解码得到对应的频率，本次本次实习中，我们需要用到的卫星信号频率为：

对于 GPS 卫星

L1： 1575.42 MHz; 具体频点L1I

L2： 1227.6 MHz; 具体频点L2I

对于 BDS 卫星

B1：1561.098 MHz; 具体频点B1I

B3：1268.52 MHz；具体频点B3I

在解码具体数值，因为不同数值所占字节数不同，所以解码的时候我们需要定义截取不同字节的函数。在C++里面，可以定义一个类来提取不同字节的数据。

表2.3 提取不同数据类型的方法

uint8_t	无符号8位整数
int8_t	有符号8位整数
uint16_t	无符号16位整数（小端）
uint32_t	无符号32位整数（小端）
int32_t	有符号32位整数（小端）
float	单精度浮点数（小端）
double	双精度浮点数（小端）

（4）广播星历数据解码

1.GPSEPHEM 7

2.BDSEPHEMERIS 1696

表2.4 BDS广播星历参数

参数	定义	数据类型	单位	解码位置
toe	星历参考时间	Ulong	s	H+72
√A	长半轴的平方根	double	m1/2	H+76
e	偏心率	double	-	H+84
ω	近地点幅角	double	π	H+92
∆n	卫星平均运动速率与计算值之差	double	π/s	H+100
M0	参考时间的平近点角	double	π	H+108
Ω0	按参考时间计算的升交点经度	double	π	H+116
Ω˙	升交点赤经变化率	double	π/s	H+124
i0	参考时间的轨道倾角	double	π	H+132
IDOT	轨道倾角变化率	double	π/s	H+140
Cuc	纬度幅角的余弦调和改正项的振幅	double	rad	H+148
Cus	纬度幅角的正弦调和改正项的振幅	double	rad	H+156
Crc	轨道半径的余弦调和改正项的振幅	double	m	H+164
Crs	轨道半径的正弦调和改正项的振幅	double	m	H+172
Cic	轨道倾角的余弦调和改正项的振幅	double	rad	H+180
Crs	轨道倾角的正弦调和改正项的振幅	double	rad	H+188

GPS卫星星历和 BDS 卫星星历存在部分差异，但是绝大部分数据内容都是相同的，因此只需要定义其中一个星历即可，不需要为每一个卫星系统都定义星历结构体，后面计算卫星位置的模块也如此，这样可以大大减少冗余代码，提高代码的可复用性。

除了两类必需的数据之外，还解码了接收机参考定位结果数据，电离层参数数据，后者主要用于单频定位

图4 解码板块

2.4 卫星位置、速度、钟差计算

卫星位置计算参见文章卫星导航算法程序设计——广播星历与卫星位置计算。

GPS卫星速度及钟差计算：

$\Delta t_{sv} = ClkBias + ClkDrift \cdot (t - t_{oc}) + ClkDriftRate \cdot (t - t_{oc})^2 + \Delta t_r$ $\Delta t_r = Fe \sqrt{A} \sin E_k$

$F = \frac{-2\sqrt{\mu}}{c^2} = -4.442807633 \cdot 10^{-10} \, (\text{sec} \cdot \text{m}^{-1/2})$

Δtr项为相对论效应改正项，ClkBias为钟偏，ClkDrift为钟速，ClkDriftRate为钟漂。按照图5进行迭代解算，迭代出卫星信号的发射时刻，并计算信号发射时刻的卫星位置。

BDS的卫星位置和速度计算可参见北斗接口文档。

图5 钟差迭代计算

2.5 误差改正

地球自转改正

对流层改正

电离层改正

1 双频消电离层模型

2 klobuchar模型

粗差探测

图6 误差改正

粗差探测之后，若探测到该观测值为粗差，则标记该观测值对应的卫星不可用。不进入后续的最小二乘解算。先地球自转改正再对流层改正，这两个顺序不能变，对流层和电离层的改正顺序可以改变。地球自转改正时不要重复改正，迭代计算改正值。在进行电离层改正时，根据选择的定位模式调用不同的函数。单频则使用Klobuchar模型，并作接收机时延的改正，双频使用双频改正。第一次迭代不计算对流层和电离层改正，因为接收机的初始位置设置为（0，0，0），无法计算高度角和方位角。同时注意对流层和电离层改正值的单位，应该以m为单位参与解算。