RTKPPP定位算法流程.doc

RTKPPP定位算法流程.doc

《RTKPPP定位算法流程.doc》由会员分享，可在线阅读，更多相关《RTKPPP定位算法流程.doc（37页珍藏版）》请在冰豆网上搜索。

1基础知识

1.1GPS精密单点定位的基本原理

GPS精密单点定位一般采用单台双频GPS接收机，利用IGS提供的精密星历和卫星钟差，基于载波相位观测值进行的高精度定位。

观测值中的电离层延迟误差通过双频信号组合消除，对流层延迟误差通过引入未知参数进行估计。

1.2时间系统

RTKLIB内部使用GPST（GPST时间）用于GNSS的数据处理和定位算法。

数据在RTKLIB内部处理之前，需要转换成GPST时间。

使用GPST的原因是避免处理润秒。

RTKLIB使用以下结构体表示时间：

typedefstruct

time_ttime；/*time（s）expressedbystandardtime_t*/

doublesec;/*fractionofsecondunder1s＊/

}gtime_t；

1.2.1GPST和UTC（UniversalTimeCoordinated）

关系参考【图1】，参考【图2】:

图1转换关系公式

图2

通过使用GPS导航信息中的UTC参数，GPST到UTC或者UTC到GPST之前的转换可以用更准确的表达方式，如【图3】。

图3

这些参数是由GPS导航消息提供的.

1.2.2BDT（北斗导航卫星系统时间）

BDT（北斗导航卫星系统时间）是一个连续的时间系统，没有润秒。

开始历元的时间是【UTC2006年1月1号00:

00:

00】。

北斗时间计算公式【图4】：

图Error!

Bookmarknotdefined.

UTC和GPST时间转换同上面的GPS一样，只不过UTC参数来自与北斗导航信息中。

1.3坐标系统

接收机和卫星的位置在RTKLIB中表示为在ECEF（地心地固坐标系）坐标系统中的X,Y，Z组件.

1.3.1大地坐标到ECEF坐标的转换

转换公式如【图5】。

第三个公式最后一行有错，应该为:

（v（1–e2）+h）sin

图4

参数说明：

a：

地球参考椭球的长半径

f:

地球参考椭球的扁平率

h:

椭球高度

:

纬度

:

经度

当前版本的RTKLIB使用的值为【图6】：

图Error!

Bookmarknotdefined.

图Error!

Bookmarknotdefined.参考椭球体

1.3.2ECEF坐系到大地坐标的转换

转换公式如【图8】

图5

1.3.3本地坐标到ECEF坐标的转换

在接收机位置的本地坐标,也被称为ENU坐标，通常使用在GNSS导航处理。

ECEF坐标到本地坐标转换的旋转矩阵表示为【图9】。

图Error!

Bookmarknotdefined.Er旋转矩阵

参数说明：

：

接收机位置的纬度

：

接收机位置的经度

通过使用Er和接收机的坐标rr【ECEF】,坐标recef【ECEF】可以被转换到本地坐标的坐标rlocal，公式如【图10】.

图6

2RTKPPP定位算法

2.1单点定位（pntpos）

1:

satpossà2:

estposà3:

estvel

1.计算计算卫星位置、速度和时钟（satposs）

a）通过广播星历计算卫星钟差（ephclk）

Ø卫星编号到卫星系统的转换（satsys）

根据卫星的编号，获取到对应的卫星导航系统。

Ø选择星历（seleph）

1、传入信号传输时间，卫星编号,导航数据等参数。

2、遍历导航数据,遍历导航数据里面的星历数据，判断星历数据的卫星编号是否和传入的卫星编号相等.

3、如果星历数据的卫星编号和传入的卫星编号相等，就计算星历参考时间（toe）和信号传输时间的时间差。

如果不相等，继续处理下一条星历数据。

4、判断计算出来的时间差，如果时间差大于了允许的最大时间差，继续查找下一个星历数据。

否则,判断时间差最小的星历数据，记录星历数据的位置。

5、返回之前记录出来的星历数据。

Ø使用广播星历计算卫星时钟偏差（eph2clk）

1、传入信号发射时刻的时间和星历数据。

2、计算信号发射时刻的时间和本时段钟差参数参考时间（星历参数toc）的时间差。

3、通过下式计算钟差，这里还没有处理相对论校正项和tgd:

（代码中有个迭代过程，资料上没看到写）

卫星钟差计算出来之后，信号发射时刻的时间还要加上这个钟差。

b）计算卫星在信号发射时刻的位置、速度和时钟（satpos）

根据星历表选项来选择不同的处理，如下：

Ø广播星历（EPHOPT_BRDC）:

ephpos（广播星历到卫星位置和钟差）

1、根据公式计算出tk；

2、根据使用的卫星系统，选择使用的地球引力常数（mu）和地球的角速度（omge）

3、根据公式计算出平近点角M。

4、求解开普勒方程，按照以下公式迭代求解.

5、根据以下公式计算出u（改正后的纬度幅角）,r（改正后的径向）,i（改正后的轨道倾角）的值.

6、根据以下公式计算卫星在轨道平面内的坐标.

7、根据不同的卫星系统，做不同的计算。

GPS计算方式：

北斗计算方式:

其中

8、按照公式计算出时间tc。

9、按照以下公式计算出钟差和钟漂。

Ø精密星历（EPHOPT_PREC）:

peph2pos

Ø广播+SBAS（EPHOPT_SBAS）:

satpos_sbas

Ø广播+SSR_APC（EPHOPT_SSRAPC）：

satpos_ssr

Ø广播+SSR_COM（EPHOPT_SSRCOM）:

satpos_ssr

ØQZSSLEX星历（EPHOPT_LEX）：

lexeph2pos

2.使用伪距估算接收机的位置，返回估算状态结果（estpos）

a）伪距残差（rescode）

Ø把ecef坐标系转换成大地坐标系（ecef2pos）

1、按照以下公式做转换，暂时还没看懂。

Ø计算几何距离和接收机到卫星的单位矢量（geodist）

1、用卫星的坐标向量做欧几里德范数，返回值和地球长半轴（WGS84）比较。

小于地球长半轴（WGS84），返回—1;

2、计算卫星坐标和接收机坐标的差值向量。

3、用差值向量做欧几里德范数，再用上一步计算出来的差值向量和计算结果做除法，得到视线向量.按照如下公式：

4、使用以下公式计算几何距离。

Ø计算卫星方位角/仰角（satazel）

1、把接收机ecef坐标转换到大地坐标；

2、判断高度是否大于地球半长轴（WGS84）的负数值；

3、如果高度小于等于地球半长轴（WGS84）的负数值，方位角为0，仰角为PI/2;

4、如果高度大于地球半长轴（WGS84）的负数值，把ECEF向量转换到局部坐标.然后对转换出来的坐标做内积。

5、使用如下公式计算卫星方位角和仰角。

Ø伪距使用编码残差改正（prange）

暂时没找到对应的文档对应.

Ø电离层改正（ionocorr）

n通过广播电离层模型（klobuchar模型）计算出电离层延迟（ionmodel）

1、校验传入的电离层模型参数，校验失败,使用默认的电离层模型参数；

2、使用以下公式计算出地球为中心的角度（半圆）；

3、使用以下公式计算子的电离层的纬度/经度（半圆）；

4、使用以下公式计算地磁纬度。

5、计算本地时间，返回值按这个公式【tt-=floor（tt/86400。

0）＊86400.0；】处理,保证tt的范围（0<=tt<86400）.

6、计算倾斜因子。

7、电离层延迟计算,公式如下。

Ø对流层改正（tropcorr）

n通过标准大气压和saastamoinen模型计算对流层延迟（tropmodel）

1、使用以下公式计算总气压。

2、使用以下公式计算绝对温度。

3、使用以下公式计算水蒸汽的分压。

代码中没有使用。

4、使用以下公式计算【Saastamoinen模型】.

Ø伪距残差

残差值=伪距—（物理距离+dtr—光速＊时钟偏差+电离层误差+对流层误差）

Ø时钟系统和接收器的偏置补偿

Ø伪距测量误差方差（varerr）

b）方差权重值（weightbyvariance）

c）最小二乘估计（lsq）

Ø最小二乘估计通过求解正规方程（X=（A＊A'）^-1＊A＊Y）

1、计算矩阵A＊Y的结果，保存到矩阵Ay；

2、计算矩阵A＊A'的结果，保存到矩阵Q；

3、求的矩阵Q的逆矩阵,结果保存到矩阵Q；

4、最后x=Q＊Ay

d）欧几里德范数（norm）

公式:

1、通过最小二乘法估算出来的参数dx【长度为4】，分别加到位置向量上x【长度为4】。

2、对这个速度向量做欧几里德范数，返回值同1E—4比较，小于这个值时，就得到估算出的接收机的位置为向量x的值。

e）验证求解（valsol）

1、对伪距残差值做内积，然后同卡方分布的自由度分布值（alpha=0。

001）做比较，大于卡方分布的自由度分布值的数据，无效。

2、计算dops；暂时还没没找到资料。

3、用计算出来的dops和配置的最大dops阈值比较，大于配置的最大dops阈值的数据无效。

3.使用多普勒估算接收机速度（estvel）

a）多普勒残差（resdop）

Ø把ecef坐标系转换成大地坐标系（ecef2pos）

Ø把xyz坐标转换成enu坐标（xyz2enu）

Ø计算ecef中的瞄准线向量

Ø计算相对于接收机在ECEF中的卫星速度

Ø多普勒残差计算

1、速率公式:

b）最小二乘估计（lsq）

Ø最小二乘估计通过求解正规方程（X=（A＊A'）^-1*A＊Y）

1、计算矩阵A*Y的结果，保存到矩阵Ay；

2、计算矩阵A＊A'的结果,保存到矩阵Q;

3、求的矩阵Q的逆矩阵，结果保存到矩阵Q；

4、最后x=Q＊Ay

c）欧几里德范数（norm）

公式：

3、通过最小二乘法估算出来的参数dx【长度为4】，分别加到速度向量上x【长度为4】。

4、对这个速度向量做欧几里德范数,返回值同1E—6比较，小于这个值时，就得到估算出的接收机速度为向量x的值。

2.2精确定位（pppos）

1：

udstate_pppà2：

satpossà3：

testeclipseà4：

res_pppà5:

res_pppà6：

filterà7：

res_ppp

1.暂时更新状态（udstate_ppp）

a）位置更新

b）时钟更新

c）对流层参数更新

d）相位偏差更新

1、通过LLI检测周跳；

2、通过Geometry-Free相位跳变检测周跳（如果双频测量值可用）；

2.计算卫星的位置和速度和时钟（satposs）

c）通过广播星历计算卫星钟差（ephclk）

Ø卫星编号到卫星系统的转换（satsys）

根据卫星的编号，获取到对应的卫星导航系统。

Ø选择星历（seleph）

1、传入信号传输时间，卫星编号，导航数据等参数。

2、遍历导航数据，遍历导航数据里面的星历数据，判断星历数据的卫星编号是否和传入的卫星编号相等。

3、如果星历数据的卫星编号和传入的卫星编号相等，就计算星历参考时间（toe）和信号传输时间的时间差。

如果不相等，继续处理下一条星历数据。

4、判断计算出来的时间差，如果时间差大于了允许的最大时间差,继续查找下一个星历数据。

否则,判断时间差最小的星历数据，记录星历数据的位置。

5、返回之前记录出来的星历数据。

Ø使用广播星历计算卫星时钟偏差（eph2clk）

1、传入信号传输时间和星历数据。

2、传入信号发射时刻的时间和星历数据。

3、计算信号发射时刻的时间和本时段钟差参数参考时间（星历参数toc）的时间差.

4、通过下式计算钟差，这里还没有处理相对论校正项和tgd:

（代码中有个迭代过程，资料上没看到写）

d）计算卫星位置、速度和时钟（satpos）

根据星历表选项来选择不同的处理，如下:

Ø广播星历（EPHOPT_BRDC）：

ephpos（广播星历到卫星位置和钟差）

1、根据公式计算出tk；

2、根据使用的卫星系统，选择使用的地球引力常数（mu）和地球的角速度（omge）

3、根据公式计算出平近点角M。

4、求解开普勒方程,按照以下公式迭代求解.

5、根据以下公式计算出u（改正后的纬度幅角）,r（改正后的径向）,i（改正后的轨道倾角）的值。

6、根据以下公式计算卫星在轨道平面内的坐标。

7、根据不同的卫星系统，做不同的计算.

GPS计算方式：

北斗计算方式:

其中

8、按照公式计算出时间tc。

9、按照以下公式计算出钟差和钟漂。

Ø精密星历（EPHOPT_PREC）:

peph2pos

Ø广播+SBAS（EPHOPT_SBAS）:

satpos_sbas

Ø广播+SSR_APC（EPHOPT_SSRAPC）:

satpos_ssr

Ø广播+SSR_COM（EPHOPT_SSRCOM）：

satpos_ssr

ØQZSSLEX星历（EPHOPT_LEX）:

lexeph2pos

3.排除对被遮蔽卫星的测量（testeclipse）

a）获得太阳和月亮在ecef中的位置（sunmoonpos）

b）归一化三维矢量（normv3）

c）欧几里德范数的向量（norm）

d）计算太阳-地球-卫星角（dot）

4.相位和编码残差（res_ppp）

a）地球潮汐改正（tidedisp）

b）计算几何距离（geodist）

1、用卫星的坐标向量做欧几里德范数,返回值和地球长半轴（WGS84）比较。

小于地球长半轴（WGS84），返回—1；

2、计算卫星坐标和接收机坐标的差值向量.

3、用差值向量做欧几里德范数，再用上一步计算出来的差值向量和计算结果做除法，得到视线向量。

4、使用以下公式计算几何距离.

备注：

这个公式是地球自转改正公式。

地球自转改正产生地球自转改正的机制在于地固系随地球的自转而旋转,地固系是非惯性系统,用地固系中的坐标计算卫星到接收机几何距离时，由于卫星位置和接收机位置是两个不同时刻的位置矢量，而且这两个时刻的地固系相对于惯性系是变化的，因此要考虑地球自转引起的距离变化.其改正公式为：

其中,w为地球自转角速度;C为光速；xsi和xri分别代表卫星位置矢量和测站位置矢量的分量,i=1,2对应x，y分量。

地球自转对纬度影响甚小,对经度影响最大,其次是高度.同时地球的自转影响还跟测站的经纬度以及测站与卫星的几何关系有关，对两极的测站，影响为零;对赤道上的测站影响最大;当卫星在测站子午面内影响为零;卫星在测站东方时,影响为负,卫星在测站西方时,影响为正。

c）卫星方位角/仰角计算（satazel）

1、把接收机ecef坐标转换到大地坐标；

2、判断高度是否大于地球半长轴（WGS84）的负数值；

3、如果高度小于等于地球半长轴（WGS84）的负数值，方位角为0，仰角为PI/2;

4、如果高度大于地球半长轴（WGS84）的负数值,把ECEF向量转换到局部坐标。

然后对转换出来的坐标做内积。

5、使用如下公式计算卫星方位角和仰角.

d）对流层延迟修正：

Ø通过标准大气压和saastamoinen模型计算对流层延迟（tropmodel）

1、使用以下公式计算总气压。

2、使用以下公式计算绝对温度。

3、使用以下公式计算水蒸汽的分压.代码中没有使用。

4、使用以下公式计算【Saastamoinen模型】。

e）卫星天线模型（satantpcv）

f）接收机天线模型，通过天线相位中心参数计算天线偏移（antmodel）

g）相位缠绕校正（windupcorr）

h）电离层和天线相位校正测量（corrmeas）

i）计算卫星时钟和对流层延迟

5.卡尔曼滤波（filter）

卡尔曼滤波按照以下的公式更新状态：

K=P＊H*（H'＊P*H+R）^-1;

xp=x+K＊v；

Pp=（I–K＊H'）＊P；

参数解释：

x：

状态向量（nx1）

p:

状态的协方差矩阵（nxn）

H：

设计矩阵的转置矩阵（nxm）

v:

创新（测量—模型）（mx1）

R:

测量误差的协方差矩阵（mxm）

n,m:

状态和测量值的数量

xp：

状态更新后的向量（nx1）

Pp:

状态更新后的协方差矩阵（nxn）

6.组合残差（res_ppp）

7.解析PPP整周模糊度（pppamb）

a）averageLC（average_LC）

b）固定宽巷模糊（fix_amb_WL）

Ø计算LC波长度（m）（lam_LC）

Ø宽巷模糊

Ø宽巷模糊度的方差

c）fixnarrow—laneambiguity【ARmode:

PPP—AR】

d）fixnarrow-laneambiguitybyILS【ARmode：

PPP-ARILS】

3以下是上面一些算法的参考资料，信息来自于Rtklib的手册

3.1接收机和卫星天线之间的几何距离

3.2卫星方向的方位角和仰角

3.3GPS、Galileo和QNSS广播星历和时钟

3.4北斗广播星历和时钟

3.5对流层和电流层模型

3.5.1对流层模型

3.5.2电离层模型

3.6单点定位

3.6.1线性最小二乘估计

3.6.2估算接收机的位置和钟差

3.6.3估算接收机的速度和钟漂

3.6.4求解验证和RAIMFDE

3.7精密定位

3.7.1用于PPP的ZD测量模式

3.7.2接收机天线相位中心模型

3.7.3卫星天线相位中心模型

3.7.4由地球潮汐引起的位移

3.7.5相位缠绕校正

3.7.6通过PPP估算接收机的位置