精华内容
下载资源
问答
  • GPS从入门到放弃(二十六) --- RTKLIB函数解析

    千次阅读 多人点赞 2020-06-07 21:57:21
    GPS从入门到放弃(二十六) — RTKLIB函数解析 为了贴合这个系列的标题“从入门到放弃”,在入门之后现在就要放弃此方向了。虽然感觉遗憾,暂时也没有办法。在此附上此系列最后一篇,希望能给大家一些帮助。 此文中...

    GPS从入门到放弃(二十六) — RTKLIB函数解析

    为了贴合这个系列的标题“从入门到放弃”,在入门之后现在就要放弃此方向了。虽然感觉遗憾,暂时也没有办法。在此附上此系列最后一篇,希望能给大家一些帮助。

    此文中一些函数解析参考了 https://www.cnblogs.com/taqikema/p/8819798.html,在此表示感谢!

    rtksvrthread

    void *rtksvrthread(void *arg)
    
    • 所在文件:rtksvr.c
    • 功能说明:rtk服务线程。
    • 参数说明:无
    • 处理过程:
    1. 检查svr->state,若为0,线程结束。
    2. 从input stream 读取数据。
    3. 写入到 log stream。
    4. 若为精密星历,调用 decodefile 解码;否则调用 decoderaw 解码。
    5. 对每一个观测数据,调用 rtkpos 进行定位计算。
    6. 若定位成功,调整时间,写solution;若没成功,写solution。
    7. 若有需要,发送 nmea request 给 base/nrtk input stream。
    8. 休眠等下一个cycle。

    rtkpos

    int rtkpos(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)
    
    • 所在文件:rtkpos.c
    • 功能说明:根据观测数据和导航信息,计算接收机的位置、速度和钟差。
    • 参数说明:
    函数参数,4个:
    rtk_t    *rtk      IO  rtk控制结构体
    obsd_t   *obs      I   观测数据
    int      n         I   观测数据的数量
    nav_t    *nav      I   导航数据
    返回类型:
    int                O   (0:no solution,1:valid solution)
    
    • 处理过程:
    1. 设置基准站位置,记录观测值数量。
    2. 调用 pntpos 进行接收机单点定位。若为单点定位模式,输出,返回。
    3. 若为 PPP 模式,调用 pppos 进行精密单点定位,输出,返回。
    4. 若无基准站观测数据,输出,返回。
    5. 若为移动基站模式,调用 pntpos 进行基站单点定位,并加以时间同步;否则只计算一下差分时间。
    6. 调用 relpos 进行相对基站的接收机定位,输出,返回。

    pntpos

    int pntpos (const obsd_t *obs, int n, const nav_t *nav, const prcopt_t *opt,
                sol_t *sol, double *azel, ssat_t *ssat, char *msg)
    
    • 所在文件:pntpos.c
    • 功能说明:依靠伪距和多普勒频移测量值来进行单点定位,给出接收机的位置、速度和钟差。
    • 参数说明:
    函数参数,8个:
    obsd_t   *obs      I   观测数据
    int      n         I   观测数据的数量
    nav_t    *nav      I   导航数据
    prcopt_t *opt      I   处理过程选项
    sol_t    *sol      IO  solution
    double   *azel     IO  方位角和俯仰角 (rad) (NULL: no output)
    ssat_t   *ssat     IO  卫星状态            (NULL: no output)
    char     *msg      O   错误消息
    返回类型:
    int                O   (1:ok,0:error)
    
    • 处理过程:
    1. 当处理选项 opt 中的模式不是单点模式时,电离层校正采用 broadcast 模型,即Klobuchar模型,对流层校正则采用 Saastamoinen 模型;相反,当其为单点模式时,对输入参数 opt 不做修改。
    2. 调用 satposs 计算卫星们位置、速度、时钟
    3. 调用 estpos 根据伪距估计接收机位置,其中会调用 valsol 进行 χ2\chi^2 检验和 GDOP 检验。
    4. 若3中的检验没通过,调用 raim_fde 进行接收机自主完好性监测,判决定位结果的有效性,并进行错误排除。
    5. 调用 estvel 根据多普勒频移测量值计算接收机的速度。
    6. 赋值给卫星状态结构体ssat。

    satposs

    void satposs(gtime_t teph, const obsd_t *obs, int n, const nav_t *nav,
                 int ephopt, double *rs, double *dts, double *var, int *svh)
    
    • 所在文件:ephemeris.c
    • 功能说明:按照所观测到的卫星顺序计算出每颗卫星的位置、速度、{钟差、频漂}
    • 参数说明:
    函数参数,9个:
    gtime_t  teph      I   time to select ephemeris (gpst)
    obsd_t   *obs      I   观测量数据
    int      n         I   观测量数据的数量
    nav_t    *nav      I   导航数据
    int      ephopt    I   星历选项 (EPHOPT_???)
    double   *rs       O   卫星位置和速度,长度为6*n,{x,y,z,vx,vy,vz}(ecef)(m,m/s)
    double   *dts      O   卫星钟差,长度为2*n, {bias,drift} (s|s/s)
    double   *var      O   卫星位置和钟差的协方差 (m^2)
    int      *svh      O   卫星健康标志 (-1:correction not available)
    返回类型:
    • 处理过程:
    1. 大循环针对每一条观测数据,按顺序处理。
    2. 首先初始化,将对当前观测数据的 rs、dts、var和svh数组的元素置 0。
    3. 通过判断某一频率下信号的伪距是否为 0,来得到此时所用的频率个数。注意,频率个数不能大于 NFREQ(默认为 3)。
    4. 用数据接收时刻减去伪距信号传播时间,得到卫星信号的发射时刻。
    5. 调用 ephclk 函数,由广播星历计算出当前观测卫星与 GPS 时间的钟差 dt。注意,此时的钟差是没有考虑相对论效应和 TGD 的。
    6. 用 4 中的信号发射时刻减去 5 中的钟差 dt,得到 GPS 时间下的卫星信号发射时刻。
    7. 调用 satpos 函数,计算信号发射时刻卫星的位置(ecef,m)、速度(ecef,m/s)、钟差((s|s/s))。注意,这里计算出的钟差是考虑了相对论效应的了,只是还没有考虑 TGD。
    8. 如果没有精密星历,则用广播星历的钟差替代。
    • 注意:
    1. 4中的公式原理:由 ρ=c(tuts)\rho = c(t_u-t_s) 可得 ts=tuρ/ct_s=t_u-\rho/c

    ephclk

    int ephclk(gtime_t time, gtime_t teph, int sat, const nav_t *nav, double *dts)
    
    • 所在文件:ephemeris.c
    • 功能说明:通过广播星历来确定卫星钟差
    • 参数说明:
    函数参数,5个:
    gtime_t  time      I   信号发射时刻
    gtime_t  teph      I   用于选择星历的时刻 (gpst)
    int      sat       I   卫星号 (1-MAXSAT)
    nav_t    *nav      I   导航数据
    double   *dts      O   卫星钟差,长度为2*n, {bias,drift} (s|s/s)
    返回类型:
    int                O     (1:ok,0:error)
    
    • 处理过程:
    1. 调用 satsys 函数,根据卫星编号确定该卫星所属的导航系统和该卫星在该系统中的 PRN编号。
    2. 对于 GPS导航系统,调用 seleph 函数来选择最接近 teph 的那个星历。
    3. 调用 eph2clk 函数,通过广播星历和信号发射时间计算出卫星钟差。
    • 注意:
    1. 此时计算出的卫星钟差是没有考虑相对论效应和 TGD的。

    eph2clk

    int eph2clk (gtime_t time, const eph_t *eph)
    
    • 所在文件:ephemeris.c
    • 功能说明:根据信号发射时间和广播星历,计算卫星钟差
    • 参数说明:
    函数参数,2个
    gtime_t  time    I   time by satellite clock (gpst)
    eph_t    *eph    I   broadcast ephemeris
    返回类型:
    double    satellite clock bias (s) without relativeity correction
    
    • 处理过程:
      1. 计算与星历参考时间的偏差 dt = t-toc。
      2. 利用二项式校正计算出卫星钟差,从 dt中减去这部分,然后再进行一次上述操作,得到最终的 dt。(这一部分不知道是为什么?)
      3. 使用二项式校正得到最终的钟差。

    satpos

    int satpos(gtime_t time, gtime_t teph, int sat, int ephopt, const nav_t *nav,
               double *rs, double *dts, double *var, int *svh)
    
    • 所在文件:ephemeris.c
    • 功能说明:计算信号发射时刻卫星的 P(ecef,m)、V(ecef,m/s)、C((s|s/s))
    • 参数说明:
    函数参数,9个:
    gtime_t time      I   time (gpst)
    gtime_t teph      I   用于选择星历的时刻 (gpst)
    int     sat       I   卫星号
    nav_t   *nav      I   导航数据
    int     ephopt    I   星历选项 (EPHOPT_???)
    double  *rs       O   卫星位置和速度,长度为6*n,{x,y,z,vx,vy,vz}(ecef)(m,m/s)
    double  *dts      O   卫星钟差,长度为2*n, {bias,drift} (s|s/s)
    double  *var      O   卫星位置和钟差的协方差 (m^2)
    int     *svh      O   卫星健康标志 (-1:correction not available)
    返回类型:
    int               O     (1:ok,0:error)
    
    • 处理过程:
    1. 根据不同的星历选项的值调用不同的处理函数,如果星历选项是 EPHOPT_BRDC,调用 ephpos 函数,根据广播星历计算出算信号发射时刻卫星的 P、V、C。如果星历选项是 EPHOPT_PREC,调用 pephpos 函数,根据精密星历和时钟计算信号发射时刻卫星的 P、V、C。
    • 注意:
    1. 此时计算出的卫星钟差考虑了相对论,还没有考虑 TGD。

    ephpos

    int ephpos(gtime_t time, gtime_t teph, int sat, const nav_t *nav,
               int iode, double *rs, double *dts, double *var, int *svh)
    
    • 所在文件:ephemeris.c
    • 功能说明:根据广播星历计算出算信号发射时刻卫星的 P、V、C
    • 参数说明:
    函数参数,9个:
    gtime_t  time      I   transmission time by satellite clock
    gtime_t  teph      I   time to select ephemeris (gpst)
    int      sat       I   卫星号 (1-MAXSAT)
    nav_t    *nav      I   导航数据
    int      iode      I   星历数据期号
    double   *rs       O   卫星位置和速度,长度为6*n,{x,y,z,vx,vy,vz}(ecef)(m,m/s)
    double   *dts      O   卫星钟差,长度为2*n, {bias,drift} (s|s/s)
    double   *var      O   卫星位置和钟差的协方差 (m^2)
    int      *svh      O   卫星健康标志 (-1:correction not available)
    返回类型:
    int                O    (1:ok,0:error)
    
    • 处理过程:
    1. 调用 satsys 函数,确定该卫星所属的导航系统。
    2. 调用 seleph 函数来选择广播星历。
    3. 根据选中的广播星历,调用 eph2pos 函数来计算信号发射时刻卫星的 位置、钟差和相应结果的误差。
    4. 在信号发射时刻的基础上给定一个微小的时间间隔,再次计算新时刻的 P、V、C。与3结合,通过扰动法计算出卫星的速度和频漂。
    • 注意:
    1. 这里是使用扰动法计算卫星的速度和频漂,并没有使用那些位置和钟差公式对时间求导的结果。
    2. 由于是调用的 eph2pos 函数,计算得到的钟差考虑了相对论效应,没有考虑 TGD。

    eph2pos

    void eph2pos(gtime_t time, const eph_t *eph, double *rs, double *dts, double *var)
    
    • 所在文件:ephemeris.c
    • 功能说明:根据广播星历计算出算信号发射时刻卫星的位置和钟差
    • 参数说明:
    函数参数,5个
    gtime_t  time      I   transmission time by satellite clock
    eph_t    *eph      I   广播星历
    double   *rs       O   卫星位置和速度,长度为6*n,{x,y,z,vx,vy,vz}(ecef)(m,m/s)
    double   *dts      O   卫星钟差,长度为2*n, {bias,drift} (s|s/s)
    double   *var      O   卫星位置和钟差的协方差 (m^2)
    返回类型:无
    
    • 处理过程:
    1. 通过卫星轨道半长轴 A 判断星历是否有效,无效则返回。
    2. 计算规化时间 tk。
    3. 根据不同卫星系统设置相应的地球引力常数 mu 和 地球自转角速度 omge。
    4. 计算平近点角 M。
    5. 用牛顿迭代法来计算偏近点角 E。参考 RTKLIB manual P145 (E.4.19)。
    6. 计算升交点角距 u。
    7. 计算摄动校正后的升交点角距 u、卫星矢径长度 r、轨道倾角 i。
    8. 计算升交点赤经 O。
    9. 计算卫星位置存入 rs 中。
    10. 计算卫星钟差,此处考虑了相对论效应,没有考虑 TGD,也没有计算钟漂。
    11. 用 URA 值来标定误差方差,具体对应关系可在 ICD-GPS-200H 20.3.3.3.1.3 SV Accuracy 中找到。

    estpos

    int estpos(const obsd_t *obs, int n, const double *rs, const double *dts,
               const double *vare, const int *svh, const nav_t *nav,
               const prcopt_t *opt, sol_t *sol, double *azel, int *vsat,
               double *resp, char *msg)
    
    • 所在文件:pntpos.c
    • 功能说明:通过伪距实现绝对定位,计算出接收机的位置和钟差,顺带返回实现定位后每颗卫星的{方位角、仰角}、定位时有效性、定位后伪距残差。
    • 参数说明:
    函数参数,13个:
    obsd_t   *obs      I   观测量数据
    int      n         I   观测量数据的数量
    double   *rs       I   卫星位置和速度,长度为6*n,{x,y,z,vx,vy,vz}(ecef)(m,m/s)
    double   *dts      I   卫星钟差,长度为2*n, {bias,drift} (s|s/s)
    double   *vare     I   卫星位置和钟差的协方差 (m^2)
    int      *svh      I   卫星健康标志 (-1:correction not available)
    nav_t    *nav      I   导航数据
    prcopt_t *opt      I   处理过程选项
    sol_t    *sol      IO  solution
    double   *azel     IO  方位角和俯仰角 (rad)
    int      *vsat     IO  卫星在定位时是否有效
    double   *resp     IO  定位后伪距残差 (P-(r+c*dtr-c*dts+I+T))
    char     *msg      O   错误消息
    返回类型:
    int                O   1表示成功,0表示出错
    
    • 处理过程:
    1. 初始化:将 sol->rr 的前 3 项位置信息(ECEF)赋值给 x 数组。
    2. 开始迭代定位计算,首先调用 rescode 函数,计算当前迭代的伪距残差 vv、几何矩阵 HH、伪距残差的方差 var、所有观测卫星的方位角和仰角 azel、定位时有效性 vsat、定位后伪距残差 resp、参与定位的卫星个数 ns 和方程个数 nv。
    3. 确定方程组中方程的个数要大于未知数的个数。
    4. 以伪距残差的标准差的倒数作为权重,对 HHvv 分别左乘权重对角阵,得到加权之后的 HHvv
    5. 调用 lsq 函数,根据 dx=(HHT)1Hvdx=(HH^T)^{-1}HvQ=(HHT)1Q=(HH^T)^{-1},得到当前 x 的修改量 dx 和定位误差协方差矩阵中的权系数阵 QQ
    6. 将 5 中求得的 dx 加入到当前 x 值中,得到更新之后的 x 值。
    7. 如果 5 中求得的修改量 dx 小于截断因子(目前是10410^{-4}),则将 6 中得到的 x 值作为最终的定位结果,对 sol 的相应参数赋值,之后再调用 valsol 函数确认当前解是否符合要求(伪距残差小于某个 χ2\chi^2值和 GDOP 小于某个门限值,参考 RTKLIB Manual P162, E.6.33, E.6.34)。否则,进行下一次循环。
    8. 如果超过了规定的循环次数,则输出发散信息后,返回 0。
    • 注意:
    1. 关于第 1步,如果是第一次定位,即输入的 sol 为空,则 x 初值为 0;如果之前有过定位,则通过 1 中操作可以将上一历元的定位值作为该历元定位的初始值。
    2. 关于定位方程,RTKLIB中的 HH 相当于定位方程解算中的 GTG^T,即 H=GTH=G^T
    3. 关于加权最小二乘,这里的权重值是对角阵,这是建立在假设不同测量值的误差之间是彼此独立的基础上的。大部分资料上这里都是把权重矩阵 W 保留到方程的解的表达式当中,而这里是直接对 H 和 v 分别左乘权重对角阵,得到加权之后的 H 和 v,其表示形式像是没有加权一样。
    4. 解方程时的 dtr 单位是 m,是乘以了光速之后的,解出结果后赋给 sol->dtr 时再除以光速。
    5. sol->time 中存储的是减去接收机钟差后的信号观测时间。

    rescode

    int rescode(int iter, const obsd_t *obs, int n, const double *rs,
                const double *dts, const double *vare, const int *svh,
                const nav_t *nav, const double *x, const prcopt_t *opt,
                double *v, double *H, double *var, double *azel, int *vsat,
                double *resp, int *ns)
    
    • 所在文件:pntpos.c
    • 功能说明:计算当前迭代的伪距残差 v、几何矩阵 H、伪距残差的方差 var、所有观测卫星的方位角和仰角 azel、定位时有效性 vsat、定位后伪距残差 resp、参与定位的卫星个数 ns 和方程个数 nv。
    • 参数说明:
    函数参数,17int      iter      I   迭代次数
    obsd_t   *obs      I   观测量数据
    int      n         I   观测量数据的数量
    double   *rs       I   卫星位置和速度,长度为6*n,{x,y,z,vx,vy,vz}(ecef)(m,m/s)
    double   *dts      I   卫星钟差,长度为2*n, {bias,drift} (s|s/s)
    double   *vare     I   卫星位置和钟差的协方差 (m^2)
    int      *svh      I   卫星健康标志 (-1:correction not available)
    nav_t    *nav      I   导航数据
    double   *x        I   本次迭代开始之前的定位值
    prcopt_t *opt      I   处理过程选项
    double   *v        O   定位方程的右端部分,伪距残差
    double   *H        O   定位方程中的几何矩阵
    double   *var      O   参与定位的伪距残差的方差
    double   *azel     O   对于当前定位值,所有观测卫星的 {方位角、高度角} (2*n)
    int      *vsat     O   所有观测卫星在当前定位时是否有效 (1*n)
    double   *resp     O   所有观测卫星的伪距残差,(P-(r+c*dtr-c*dts+I+T)) (1*n)
    int      *ns       O   参与定位的卫星的个数
    返回类型:
    int                O   定位方程组的方程个数
    
    • 处理过程:
    1. 将之前得到的定位解信息赋值给 rr 和 dtr 数组,以进行关于当前解的伪距残差的相关计算。
    2. 调用 ecef2pos 函数,将上一步中得到的位置信息由 ECEF 转化为大地坐标系。
    3. 将 vsat、azel 和 resp 数组置 0,因为在前后两次定位结果中,每颗卫星的上述信息都会发生变化。
    4. 调用 satsys 函数,验证卫星编号是否合理及其所属的导航系统。
    5. 检测当前观测卫星是否和下一个相邻数据重复;重复则不处理这一条,去处理下一条。
    6. 调用 geodist 函数,计算卫星和当前接收机位置之间的几何距离 rr 和接收机到卫星方向的观测矢量。然后检验几何距离是否 >0。此函数中会进行地球自转影响的校正(Sagnac效应)。
    7. 调用 satazel 函数,计算在接收机位置处的站心坐标系中卫星的方位角和仰角;若仰角低于截断值,不处理此数据。
    8. 调用 prange 函数,得到经过DCB校正后的伪距值 ρ\rho
    9. 可以在处理选项中事先指定定位时排除哪些导航系统或卫星,这是通过调用 satexclude 函数完成的。
    10. 调用 ionocorr 函数,计算电离层延时 II (m)。所得的电离层延时是建立在 L1 信号上的,当使用其它频率信号时,依据所用信号频组中第一个频率的波长与 L1 波长的关系,对上一步得到的电离层延时进行修正。
    11. 调用 tropcorr 函数,计算对流层延时 TT (m)。
    12. ρ(r+dtrcdts+I+T)\rho-(r+dt_r-c·dt_s+I+T),计算出此时的伪距残差。
    13. 组装几何矩阵 HH,前 3 行为 6 中计算得到的视线单位向量的反向,第 4 行为 1,其它行为 0。
    14. 处理不同系统(GPS、GLO、GAL、CMP)之间的时间偏差,修改矩阵 HH
    15. 将参与定位的卫星的定位有效性标志设为 1,给当前卫星的伪距残差赋值,参与定位的卫星个数 ns 加 1。
    16. 调用 varerr 函数,计算此时的导航系统误差,然后累加计算用户测距误差(URE)。
    17. 为了防止不满秩的情况,把矩阵 HH 补满秩了。
    • 注意:
    1. 输入参数x的size为7*1,前3个是本次迭代开始之前的定位值,第4个是钟差,后三个分别是gps系统与glonass、galileo、北斗系统的钟差。但是从代码看,后三个没有用上。
    2. 返回值 v和 resp的主要区别在于长度不一致, v是需要参与定位方程组的解算的,维度为 nv*1;而 resp仅表示所有观测卫星的伪距残余,维度为 n*1,对于没有参与定位的卫星,该值为 0。

    raim_fde

    int raim_fde(const obsd_t *obs, int n, const double *rs,
                 const double *dts, const double *vare, const int *svh,
                 const nav_t *nav, const prcopt_t *opt, sol_t *sol,
                 double *azel, int *vsat, double *resp, char *msg)
    
    • 所在文件:pntpos.c
    • 功能说明:使用伪距残差判决法对计算得到的定位结果进行接收机自主正直性检测(RAIM),每次舍弃一颗卫星测量值,用剩余的值组成一组进行定位运算,选择定位后伪距残差最小的一组作为最终结果。这样如果只有一个异常观测值的话,这个错误可以被排除掉;有两个或以上错误则排除不了。注意这里只会在对定位结果有贡献的卫星数据进行检测。
    • 参数说明:
    函数参数,13个:
    obsd_t   *obs      I   观测数据
    int      n         I   观测数据的数量
    double   *rs       I   卫星位置和速度,长度为6*n,{x,y,z,vx,vy,vz}(ecef)(m,m/s)
    double   *dts      I   卫星钟差,长度为2*n, {bias,drift} (s|s/s)
    double   *vare     I   卫星位置和钟差的协方差 (m^2)
    int      *svh      I   卫星健康标志 (-1:correction not available)
    nav_t    *nav      I   导航数据
    prcopt_t *opt      I   处理过程选项
    sol_t    *sol      IO  solution
    double   *azel     IO  方位角和俯仰角 (rad)
    int      *vsat     IO  卫星在定位时是否有效
    double   *resp     IO  定位后伪距残差 (P-(r+c*dtr-c*dts+I+T))
    char     *msg      O   错误消息
    返回类型:
    int                O   (1:ok,0:error)
    
    • 处理过程:
    1. 大循环是每次舍弃第 i 颗卫星。
    2. 舍弃第 i 颗卫星后,将剩下卫星的数据复制到一起,调用 estpos 函数计算使用剩下卫星进行定位的定位值。
    3. 累加使用当前卫星实现定位后的伪距残差平方和与可用卫星数目,如果 nvsat<5,则说明当前卫星数目过少,无法进行 RAIM_FDE 操作。
    4. 计算伪距残差平方和的标准差,如果小于 rms,则说明当前定位结果更合理,将 stat 置为 1,重新更新 sol、azel、vsat(当前被舍弃的卫星,此值置为0)、resp等值,并将当前的 rms_e更新到 rms 中。
    5. 继续弃用下一颗卫星,重复 2-4 操作。总而言之,将同样是弃用一颗卫星条件下,伪距残差标准平均值最小的组合所得的结果作为最终的结果输出。
    6. 如果 stat不为 0,则说明在弃用卫星的前提下有更好的解出现,输出信息,指出弃用了哪颗卫星。
    • 注意:
    1. 源码中有很多关于 i、j、k的循环。其中,i表示最外面的大循环,每次将将第 i颗卫星舍弃不用,这是通过 if (j==i) continue实现的;j表示剩余使用的卫星的循环,每次进行相应数据的赋值;k表示参与定位的卫星的循环,与 j一起使用。

    estvel

    void estvel(const obsd_t *obs, int n, const double *rs, const double *dts,
                const nav_t *nav, const prcopt_t *opt, sol_t *sol,
                const double *azel, const int *vsat)
    
    • 所在文件:pntpos.c
    • 功能说明:依靠多普勒频移测量值计算接收机的速度,用牛顿迭代法。
    • 参数说明:
    函数参数,9个:
    obsd_t   *obs      I   观测数据
    int      n         I   观测数据的数量
    double   *rs       I   卫星位置和速度,长度为6*n,{x,y,z,vx,vy,vz}(ecef)(m,m/s)
    double   *dts      I   卫星钟差,长度为2*n, {bias,drift} (s|s/s)
    nav_t    *nav      I   导航数据
    prcopt_t *opt      I   处理过程选项
    sol_t    *sol      IO  solution
    double   *azel     IO  方位角和俯仰角 (rad)
    int      *vsat     IO  卫星在定位时是否有效
    返回类型:
    int                O     (1:ok,0:error)
    
    • 处理过程:
    1. 在最大迭代次数限制内,调用 resdop,计算定速方程组左边的几何矩阵和右端的速度残余,返回定速时所使用的卫星数目。
    2. 调用最小二乘法 lsq 函数,解出{速度、频漂}的步长 dx,累加到 x 中。
    3. 检查当前计算出的步长的绝对值是否小于 1E-6。是,则说明当前解已经很接近真实值了,将接收机三个方向上的速度存入到 sol->rr 中;否,则进行下一次循环。
    • 注意:
    1. 最终向 sol_t 类型存储定速解时,并没有存储所计算出的接收器时钟频漂。
    2. 这里不像定位时,初始值可能为上一历元的位置(从 sol 中读取初始值),这里定速的初始值直接给定为 0.

    resdop

    int resdop(const obsd_t *obs, int n, const double *rs, const double *dts,
               const nav_t *nav, const double *rr, const double *x,
               const double *azel, const int *vsat, double *v, double *H)
    
    • 所在文件:pntpos.c
    • 功能说明:计算定速方程组左边的几何矩阵和右端的速度残余,返回定速时所使用的卫星数目
    • 参数说明:
    函数参数,11个:
    obsd_t   *obs      I   观测数据
    int      n         I   观测数据的数量
    double   *rs       I   卫星位置和速度,长度为6*n,{x,y,z,vx,vy,vz}(ecef)(m,m/s)
    double   *dts      I   卫星钟差,长度为2*n, {bias,drift} (s|s/s)
    nav_t    *nav      I   导航数据
    double   *rr       I   接收机位置和速度,长度为6{x,y,z,vx,vy,vz}(ecef)(m,m/s)
    double   *x        I   本次迭代开始之前的定速值,长度为4{vx,vy,vz,drift}
    double   *azel     IO  方位角和俯仰角 (rad)
    int      *vsat     I   卫星在定速时是否有效
    double   *v        O   定速方程的右端部分,速度残差
    double   *H        O   定速方程中的几何矩阵
    返回类型:
    int                O   定速时所使用的卫星数目
    
    • 处理过程:
    1. 调用 ecef2pos 函数,将接收机位置由 ECEF 转换为大地坐标系。
    2. 调用 xyz2enu 函数,计算此时的坐标转换矩阵。
    3. 去除在定速时不可用的卫星。
    4. 计算当前接收机位置下 ENU中的视向量,然后转换得到 ECEF 中视向量的值。
    5. 计算 ECEF 中卫星相对于接收机的速度。
    6. 计算考虑了地球自转的用户和卫星之间的几何距离变化率,校正公式见 RTKLIB manual P159 (F.6.29),此公式可由 P140 (E.3.8b) 对时间求导得到。
    7. 根据公式计算出定速方程组右端项的多普勒残差,构建左端项的几何矩阵,最后再将观测方程数增 1.
    • 注意:
    1. 这里与定位不同,构建几何矩阵时,就只有 4个未知数,而定位时是有 NX个。
    2. 多普勒定速方程中几何矩阵 G 与定位方程中的一样。
    3. 第7步中计算多普勒残差 b 与很多资料不同,因为 estvel 中用的是牛顿迭代法,其最小二乘法并不是直接求解x,而是求解dx,再加到x上。下面式子可以参考,其中 xryrx_r、y_r 分别为接收机位置的x,y分量,xsysx_s、y_s 分别为卫星s位置的x,y分量,rsr_s 为接收机到卫星s的距离。注意在计算 G 的时候忽略了地球自转校正项。
      x=[v,cδt˙]T=[vx,vy,vz,cδt˙]T \boldsymbol{x} = [\boldsymbol{v},\:c\cdot\dot{\delta_t}]^T=[v_x,\:v_y,\:v_z,\:c\cdot\dot{\delta_t}]^T

    rs=(vsv)es+ωec(vs,yxr+ysvxvs,xyrxsvy) r_{s} = (\boldsymbol{v_s}-\boldsymbol{v})\cdot\boldsymbol{e_s} + \frac{\omega_e}{c}(v_{s,y}x_r+y_sv_x-v_{s,x}y_r-x_sv_y)

    y=λfd=h(x)=rs+cδt˙cδ˙t,s y = -\lambda f_d = h(x) = r_{s} + c\cdot\dot{\delta_t} - c\cdot\dot{\delta}_{t,s}

    G=hx=[e1,k1e2,k1e3,k1e4,k1] \boldsymbol{G} = \frac{\partial h}{\partial x} = \left[ \begin{array}{cc} -\boldsymbol{e}_{1,k} & 1\\ -\boldsymbol{e}_{2,k} & 1\\ -\boldsymbol{e}_{3,k} & 1\\ -\boldsymbol{e}_{4,k} & 1 \end{array} \right]

    b=yh(x)=λfd(rs+cδt˙cδ˙t,s) b = y-h(x) = -\lambda f_d - (r_{s} + c\cdot\dot{\delta_t} - c\cdot\dot{\delta}_{t,s})

    pppos

    void pppos(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)
    
    • 所在文件:ppp.c
    • 功能说明:精确点定位
    • 参数说明:
    函数参数,4个:
    rtk_t    *rtk      IO  rtk控制结构体
    obsd_t   *obs      I   观测数据
    int      n         I   观测数据的数量
    nav_t    *nav      I   导航数据
    返回类型:
    • 处理过程:
    1. 调用 udstate_ppp 更新状态值 rtk->x 及其误差协方差 rtk->P。
    2. 调用 satposs 计算卫星位置和钟差。
    3. 若有需要排除观测的卫星,调用 testeclipse 进行排除。
    4. 根据设置的滤波迭代次数,循环调用 res_ppp 计算载波相位和伪距残差 v 和观测矩阵 H,以及测量误差的协方差 R,并将残差值存入 rtk->ssat[sat-1].resc 和 rtk->ssat[sat-1].resp 中。然后调用 filter 进行 Kalman 滤波运算。
    5. 4 中循环结束后,再次调用 res_ppp,更新残差值,并将残差值存入 rtk->ssat[sat-1].resc 和 rtk->ssat[sat-1].resp 中。
    6. 对某些周整模糊度计算模式,调用 pppamb 进行周整模糊度解算。
    7. 更新 solution 状态。
    • 注意:
    1. 状态变量包含接收机位置、接收机速度、接收机钟差、[对流层参数]、每颗卫星的载波偏移。(参考 RTKLIB Manual P177 E.8.16)。其中载波偏移包含周整模糊度以及小数部分,可参考 RTKLIB Manual P139 E.3.5。
    2. 过程5中的目的是在过程6中的某些周整模糊度计算中会用上,以及一些输出中可能会用上。

    udstate_ppp

    void udstate_ppp(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)
    
    • 所在文件:ppp.c
    • 功能说明:更新状态值 rtk->x。
    • 参数说明:
    函数参数,4个:
    rtk_t    *rtk      IO  rtk控制结构体
    obsd_t   *obs      I   观测数据
    int      n         I   观测数据的数量
    nav_t    *nav      I   导航数据
    返回类型:
    • 处理过程:
    1. 调用 udpos_ppp 根据不同模式初始化状态 rtk->x 中的位置值。
    2. 调用 udclk_ppp 初始化状态 rtk->x 中的钟差值(6个,因有6个系统)。
    3. 对某些对流层模型,初始化状态 rtk->x 中的对流层参数。
    4. 调用 udbias_ppp 更新载波相位偏移状态值以及其误差协方差。
    • 注意:
    1. 状态变量包含接收机位置、接收机速度、接收机钟差、[对流层参数]、每颗卫星的载波偏移。(参考 RTKLIB Manual P177 E.8.16)。其中载波偏移包含周整模糊度以及小数部分,可参考 RTKLIB Manual P139 E.3.5。
    2. 处理过程 1 的 udpos_ppp 中,PPP-static 模式只用 rtk->sol.rr 初始化一次状态 rtk->x,而 PPP-Kinematic 模式每次都会用 rtk->sol.rr 重新初始化 rtk->x。而 rtk->sol.rr 的值来源于单点定位。
    3. 处理过程 1 的 udpos_ppp 中,PPP-Fixed 模式只用于残差分析,不用于定位。

    udbias_ppp

    void udbias_ppp(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)
    
    • 所在文件:ppp.c
    • 功能说明:从观测值 obs 经过计算得到偏移值存入状态 rtk->x 中,若有周跳需要更新状态 rtk->x。更新偏移值的误差协方差到 rtk->P。
    • 参数说明:
    函数参数,4个:
    rtk_t    *rtk      IO  rtk控制结构体
    obsd_t   *obs      I   观测数据
    int      n         I   观测数据的数量
    nav_t    *nav      I   导航数据
    返回类型:
    • 处理过程:
    1. 调用 detslp_ll 通过 LLI 检查是否有周跳。
    2. 调用 detslp_gf 通过 geometry-free phase jump 检查是否有周跳。
    3. 若观测中断需要 reset 载波偏移值。
    4. 接收机位置转换到大地坐标系。
    5. 对每一组观测数据,调用 corrmeas 计算观测值 meas,是否有周跳,然后计算出偏移值 bias。
    6. 若载波和伪距跳变太大,为了保持一致性,需要进行校正。
    7. 对每一组观测值更新偏移值的误差协方差到 rtk->P。若有周跳,或者 rtk->x 状态还未初始化,则用偏移值 bias 重置 rtk->x。

    res_ppp

    int res_ppp(int iter, const obsd_t *obs, int n, const double *rs,
                const double *dts, const double *vare, const int *svh,
                const nav_t *nav, const double *x, rtk_t *rtk, double *v,
                double *H, double *R, double *azel)
    
    • 所在文件:ppp.c
    • 功能说明:计算载波相位和伪距残差 v 和观测矩阵 H,以及测量误差的协方差 R。并将残差值存入 rtk->ssat[sat-1].resc 和 rtk->ssat[sat-1].resp 中。
    • 参数说明:
    函数参数:14int     iter    I   迭代次数(该函数中没有用到)   
    obsd_t  *obs    I   观测数据
    int      n      I   观测数据的数量
    double  *rs     I   卫星位置和速度,长度为6*n,{x,y,z,vx,vy,vz}(ecef)(m,m/s)
    double  *dts    I   卫星钟差,长度为2*n, {bias,drift} (s|s/s)
    double  *vare   I   卫星位置和钟差的协方差 (m^2)
    int     *svh    I   卫星健康标志 (-1:correction not available)
    nav_t   *nav    I   导航数据
    double  *x      I   状态变量
    rtk_t   *rtk    IO  rtk控制结构体
    double  *v      O   实际观测量与预测观测量的残差 (2n个,因每个观测数据有phase和code两个观测值)
    double  *H      O   观测矩阵
    double  *R      O   测量误差的协方差
    double  *azel   O   方位角和俯仰角 (rad)
    返回类型:
    int             O   残差个数,<=0表示失败
    
    • 处理过程
    1. 接收机位置传给 rr。
    2. 若需要地球潮校正,调用 tidedisp 对 rr 进行校正。地球潮包含固体潮、极潮和海潮负荷。
    3. rr 转换坐标系到大地坐标系 pos。
    4. 大循环,对每一个观测量
    5. 计算卫星到接收机的几何距离 r,计算仰角,若仰角低于阈值,排除此观测量。
    6. 若设置有卫星需要排除的,排除掉。
    7. 根据对流层模型设置,计算对流层延时校正值 dtrp。
    8. 调用 satantpcv 根据卫星天线模型计算校正值 dants。
    9. 调用 antmodel 根据接收机天线模型计算校正值 dantr。
    10. 调用 windupcorr 计算相位缠绕校正值存入 rtk->ssat[sat-1].phw 中。
    11. 调用 corrmeas 计算经过电离层和天线相位校正后的测量值 meas(包含相位和伪距两个值),把 8~10 步中的计算值都合并到 meas 中了。
    12. 对几何距离 r 进行卫星钟差和对流层校正。
    13. 构造残差 v 和观测矩阵 H。v 由经过电离层和天线相位校正后的测量值 meas 减去经过卫星钟差和对流层校正后的 r,再减去接收机钟差(若 meas 为载波,再减去载波偏移)得到。 H 参考 RTKLIB Manual P177 E.8.21。并将残差值存入 rtk->ssat[sat-1].resc 和 rtk->ssat[sat-1].resp 中。
    14. 重复 5~13 直到大循环结束。
    15. 计算测量误差的协方差 R。
    • 注意:
    1. 状态变量包含接收机位置、接收机速度、接收机钟差、[对流层参数]、载波偏移。(参考 RTKLIB Manual P177 E.8.16)。其中载波偏移包含周整模糊度以及小数部分,可参考 RTKLIB Manual P139 E.3.5。

    satantpcv

    void satantpcv(const double *rs, const double *rr, const pcv_t *pcv,
                    double *dant)
    
    • 所在文件:ppp.c
    • 功能说明:根据模型计算卫星天线偏移校正值dant。
    • 参数说明:
    函数参数:4double   *rs     I   卫星位置和速度,长度为6*n,{x,y,z,vx,vy,vz}(ecef)(m,m/s)
    double   *rr     I   接收机位置和速度,长度为6*n,{x,y,z,vx,vy,vz}(ecef)(m,m/s)
    pcv_t    *pcv    I   天线相位中心参数结构体
    double   *dant   O   卫星天线校正值
    返回类型:无
    
    • 处理过程:
    1. 根据卫星位置和接收机位置计算天底角。
    2. 调用 antmodel_s 计算卫星天线偏移校正值 dant。
    • 注意:
    1. antmodel_s 中的计算只是对 phase center variation 进行了简单的插值,并没有考虑 phase center offset。这里存疑!?
    2. 类似的对接收机天线偏移进行校正计算的 antmodel 函数中则不仅考虑了phase center variation,也同时考虑了phase center offset。

    antmodel

    void antmodel(const pcv_t *pcv, const double *del, const double *azel,
                  int opt, double *dant)
    
    • 所在文件:rtkcmn.c
    • 功能说明:根据模型计算接收机天线偏移校正值dant。
    • 参数说明:
    函数参数:5个
    pcv_t    *pcv    I   天线相位中心参数结构体
    double   *del    I   相对天线参考点偏移值
    double   *azel   I   方位角和俯仰角
    int       opt    I   选项(非0则需要考虑pcv)
    double   *dant   O   接收机天线校正值
    返回类型:无
    
    • 处理过程:
    1. 根据方位角和俯仰角计算单位观测矢量。
    2. 对每一个频率,先计算偏移,再根据opt值看是否要加上pcv的影响,最后得出dant。

    windupcorr

    void windupcorr(gtime_t time, const double *rs, const double *rr,
                    double *phw)
    
    • 所在文件:rtkcmn.c
    • 功能说明:根据模型计算相位缠绕校正值 phw。
    • 参数说明:
    函数参数:4个
    gtime_t time     I   time (GPST)
    double  *rs      I   卫星位置 (ecef) {x,y,z} (m)
    double  *rr      I   接收机位置 (ecef) {x,y,z} (m)
    double  *phw     IO  相位缠绕校正值 (cycle)
    返回类型:无
    
    • 处理过程:
    1. 调用 sunmoonpos 获取太阳的位置。
    2. 计算卫星到接收机的单位矢量 ek。
    3. 计算卫星天线坐标系的三个轴的单位矢量 exs, eys, ezs。
    4. 计算站点坐标系的三个轴的单位矢量 exr, eyr。
    5. 计算有效偶极 ds 和 dr。
    6. 计算相位缠绕校正值 phw。
    • 注意:
    1. 使用这个函数时需要传入之前的相位缠绕校正值,因为此函数假设校正值不会跳变超过0.5个载波周期。

    corrmeas

    int corrmeas(const obsd_t *obs, const nav_t *nav, const double *pos,
                const double *azel, const prcopt_t *opt,
                const double *dantr, const double *dants, double phw,
                double *meas, double *var, int *brk)
    
    • 所在文件:ppp.c
    • 功能说明:计算经过电离层、DCB、卫星天线、接收机天线、相位缠绕校正后的测量值 meas(包含相位和伪距两个值)。
    • 参数说明:
    函数参数:11个
    obsd_t   *obs    I   观测数据
    nav_t    *nav    I   导航数据
    double   *pos    I   接收机位置 (lat,lon,h)(rad,m)
    double   *azel   I   方位角和俯仰角 (rad)
    prcopt_t *opt    I   处理过程选项
    double   *dantr  I   接收机天线校正值
    double   *dants  I   卫星天线校正值
    double   phw     I   相位缠绕校正值
    double   *meas   O   校正后的测量值
    double   *var    O   校正后的测量值的误差协方差
    int      *brk    O   用与判断是否有周跳
    返回类型:
    int              O   (>0:ok,0:sth wrong)
    
    • 处理过程:
    1. 若电离层校正模式为 iono-free LC,调用 ifmeas 计算 meas 值,然后返回。否则进行下面的步骤。
    2. 调用 testsnr 看观测值的信噪比是否过低,过低则忽略此观测值。
    3. 进行 DCB 校正。
    4. 调用 corr_ion 计算 slant ionospheric delay 值 ion。
    5. 合并各种校正值得到 meas 值。

    ifmeas

    int ifmeas(const obsd_t *obs, const nav_t *nav, const double *azel,
              const prcopt_t *opt, const double *dantr, const double *dants,
              double phw, double *meas, double *var)
    
    • 所在文件:ppp.c
    • 功能说明:计算经过电离层、DCB、卫星天线、接收机天线、相位缠绕校正后的测量值 meas(包含相位和伪距两个值)。电离层延时计算需要双频点数据。
    • 参数说明:
    函数参数:9个
    obsd_t   *obs    I   观测数据
    nav_t    *nav    I   导航数据
    double   *azel   I   方位角和俯仰角 (rad)
    prcopt_t *opt    I   处理过程选项
    double   *dantr  I   接收机天线校正值
    double   *dants  I   卫星天线校正值
    double   phw     I   相位缠绕校正值
    double   *meas   O   校正后的测量值
    double   *var    O   校正后的测量值的误差协方差
    返回类型:
    int              O   (>0:ok,0:sth wrong)
    
    • 处理过程:
    1. 选择所用频段,没有足够可用的频段则返回。
    2. 调用 testsnr 看观测值的信噪比是否过低,过低则返回。
    3. 按公式计算 γ,c1,c2\gamma, c_1, c_2, 获取L1, L2, P1, P2, P1_C1, P1_P2的值。
    4. 按公式计算电离层校正和相位缠绕校正后的载波测量值。
    5. 按公式计算电离层校正和DCB校正后的伪距测量值。
    6. 若有SBAS,加上SBAS钟差校正到伪距测量值。
    7. 若有GLONASS,加上GPS和GLONASS的硬件偏差校正到伪距测量值。
    8. 加上卫星天线、接收机天线校正值到载波测量值和伪距测量值。

    filter

    int filter(double *x, double *P, const double *H, const double *v,
              const double *R, int n, int m)
    
    • 所在文件:rtkcmn.c
    • 功能说明:kalman滤波运算
    • 参数说明:
    函数参数,7个:
    double   *x        IO 状态变量 (n x 1)
    double   *P        IO 状态变量的误差协方差阵 (n x n)
    double   *H        I  观测矩阵的转置 (n x m)
    double   *v        I  实际观测量与预测观测量的残差 (measurement - model) (m x 1)
    double   *R        I  测量误差的协方差 (m x m)
    int      n         I  状态变量个数
    int      m         I  观测值个数
    返回类型:
    int                O (0:ok,<0:error)
    
    • 处理过程:
    1. 选择需要更新的状态 x 和对应的 P、H 到 x_、p_、H_ 中。
    2. 调用 filter_ 进行kalman滤波更新。
    3. 将更新值存到 x、P中。
    • 注意:
    1. 若状态 x[i]==0.0, 则不会更新 x[i] 和 P[i+i*n]

    filter_

    int filter_(const double *x, const double *P, const double *H,
               const double *v, const double *R, int n, int m,
               double *xp, double *Pp)
    
    • 所在文件:rtkcmn.c
    • 功能说明:kalman滤波运算,K=PH(HTPH+R)1,xp=x+Kv,Pp=(IKHT)PK=PH(H^TPH+R)^{-1}, xp=x+Kv, Pp=(I-KH^T)P
    • 参数说明:
    函数参数,9个:
    double   *x        I  状态变量 (n x 1)
    double   *P        I  状态变量的误差协方差阵 (n x n)
    double   *H        I  观测矩阵的转置 (n x m)
    double   *v        I  实际观测量与预测观测量的残差 (measurement - model) (m x 1)
    double   *R        I  测量误差的协方差 (m x m)
    int      n         I  状态变量个数
    int      m         I  观测值个数
    double   *xp       O  更新后的状态变量 (n x 1)
    double   *Pp       O  更新后的状态变量的误差协方差阵 (n x n)
    返回类型:
    int                O (0:ok,<0:error)
    
    • 处理过程:
    1. 调用矩阵运算函数按照公式进行矩阵运算
    • 注意:
    1. 矩阵是按列优先的顺序存的 (fortran convention)

    relpos

    int relpos(rtk_t *rtk, const obsd_t *obs, int nu, int nr, const nav_t *nav)
    
    • 所在文件:rtkpos.c
    • 功能说明:相对定位
    • 参数说明:
    函数参数,5个:
    rtk_t    *rtk      IO  rtk控制结构体
    obsd_t   *obs      I   观测数据
    int      nu        I   接收机观测数据的数量
    int      nr        I   基站观测数据的数量
    nav_t    *nav      I   导航数据
    返回类型:
    int                O   (1:ok,0:error)
    
    • 处理过程:
    1. 一系列值、状态、中间变量的初始化。
    2. 调用 satposs 计算卫星们的位置、速度和钟差。
    3. 调用 zdres 计算基站的没有差分的相位/码残差,若出错则返回0。
    4. 若为后处理,需要插值的,调用 intpres 进行插值。
    5. 调用 selsat 选择接收机与基站共同观测的卫星,返回共同观测的卫星个数,输出卫星号列表sat、在接收机观测值中的index值列表 iu 和在基站观测值中的index值列表 ir。
    6. 调用 udstate 更新状态值 rtk->x 及其误差协方差 rtk->P。
    7. 按迭代次数循环第8~10步。
    8. 调用 zdres 计算接收机的没有差分的相位/码残差。
    9. 调用 ddres 计算双差相位/码残差。
    10. 调用 filter 进行 Kalman 滤波运算。
    11. 再次调用 zdresddres 计算双差相位/码残差,调用 valpos 进行验证,若通过则更新 rtk->x 以及 rtk->P,并更新模糊度控制结构体。
    12. 根据选项模式调用不同的 resamb_* 函数解析周整模糊度。
    13. 保存solution状态。

    zdres

    int zdres(int base, const obsd_t *obs, int n, const double *rs,
              const double *dts, const int *svh, const nav_t *nav,
              const double *rr, const prcopt_t *opt, int index, double *y,
              double *e, double *azel)
    
    • 所在文件:rtkpos.c
    • 功能说明:计算接收机或基站的没有差分的相位/码残差(Zero-Difference Residuals)
    • 参数说明:
    函数参数,13个:
    int      base      I   0表示接收机,1表示基站
    obsd_t   *obs      I   观测数据
    int      n         I   观测数据的数量
    double   *rs       I   卫星位置和速度,长度为6*n,{x,y,z,vx,vy,vz}(ecef)(m,m/s)   
    double   *dts      I   卫星钟差,长度为2*n, {bias,drift} (s|s/s)
    int      *svh      I   卫星健康标志 (-1:correction not available)   
    nav_t    *nav      I   导航数据
    double   *rr       I   接收机/基站的位置和速度,长度为6*n,{x,y,z,vx,vy,vz}(ecef)(m,m/s)
    prcopt_t *opt      I   处理过程选项
    int      index     I   0表示接收机,1表示基站,与参数 base 重复了
    double   *y        O   相位/码残差
    double   *e        O   观测矢量 (ecef)
    double   *azel     O   方位角和俯仰角 (rad)
    返回类型:
    int                O   (1:ok,0:error)
    
    • 处理过程:
    1. 若没有接收机位置,返回0。
    2. 接收机位置传给 rr_。
    3. 若需要地球潮校正,调用 tidedisp 对 rr_ 进行校正。地球潮包含固体潮、极潮和海潮负荷。
    4. rr_ 转换坐标系到大地坐标系 pos。
    5. 大循环,对每一个观测量
    6. 调用 geodist 计算卫星到接收机的几何距离 r,调用 satazel 计算仰角,若仰角低于阈值,排除此观测量。
    7. 若设置有卫星需要排除的,排除掉。
    8. 根据卫星钟差校正 r。
    9. 根据对流层模型设置,调用 tropmodel 和 tropmapf 计算对流层延时校正值并校正 r。
    10. 根据接收机天线模型调用 antmodel 计算校正值 dant(对每一个频率都有一个值)。
    11. 调用 zdres_sat 计算没有差分的相位/码残差 y。
    12. 重复6~11直到大循环结束。

    zdres_sat

    void zdres_sat(int base, double r, const obsd_t *obs, const nav_t *nav,
                  const double *azel, const double *dant,
                  const prcopt_t *opt, double *y)
    
    • 所在文件:rtkpos.c
    • 功能说明:计算接收机或基站对某一颗卫星的没有差分的相位/码残差(Zero-Difference Residuals)。y = 观测值 - r - dant。
    • 参数说明:
    函数参数,13个:
    int      base      I   0表示接收机,1表示基站
    double   r         I   经过钟差和对流层校正后的几何距离。
    obsd_t   *obs      I   观测数据
    nav_t    *nav      I   导航数据
    double   *azel     I   方位角和俯仰角 (rad)
    double   *dant     I   接收机天线校正值
    prcopt_t *opt      I   处理过程选项
    double   *y        O   相位/码残差
    返回类型:
    • 处理过程:
    1. 按电离层校正模式是否为 IONOOPT_IFLC 分两种情况。
    2. 若是:检测信噪比,计算天线校正值 dant_if,然后计算残差。
    3. 若否:检测信噪比,然后计算残差。

    tidedisp

    void tidedisp(gtime_t tutc, const double *rr, int opt, const erp_t *erp,
                  const double *odisp, double *dr)
    
    • 所在文件:ppp.c
    • 功能说明:计算因地球潮汐而引起的站点位移校正值 dr。
    • 参数说明:
    函数参数,6个:
    gtime_t  tutc     I   time in utc
    double  *rr       I   站点位置 (ecef) (m)
    int      opt      I   选项(指定包含哪些潮的影响)
                          1: solid earth tide
                          2: ocean tide loading
                          4: pole tide
                          8: elimate permanent deformation
    double  *erp      I   地球自转参数
    double  *odisp    I   海潮负荷参数
                          odisp[0+i*6]: consituent i amplitude radial(m)
                          odisp[1+i*6]: consituent i amplitude west  (m)
                          odisp[2+i*6]: consituent i amplitude south (m)
                          odisp[3+i*6]: consituent i phase radial  (deg)
                          odisp[4+i*6]: consituent i phase west    (deg)
                          odisp[5+i*6]: consituent i phase south   (deg)
                          (i=0:M2,1:S2,2:N2,3:K2,4:K1,5:O1,6:P1,7:Q1,8:Mf,9:Mm,10:Ssa)
    double  *dr       O   因地球潮汐而引起的站点位移校正值 (ecef) (m)
    返回类型:
    • 处理过程:
    1. 若有erp,调用 geterp 获取erp参数。
    2. 若有选项,调用 tide_solid 计算固体潮。
    3. 若有选项,调用 tide_oload 计算海潮负荷。
    4. 若有选项,调用 tide_pole 计算极潮。

    geterp

    int geterp(const erp_t *erp, gtime_t time, double *erpv)
    
    • 所在文件:rtkcmn.c
    • 功能说明:获取ERP参数值。
    • 参数说明:
    函数参数,3个:
    erp_t  *erp        I   earth rotation parameters
    gtime_t time       I   time (gpst)
    double *erpv       O   erp values {xp,yp,ut1_utc,lod} (rad,rad,s,s/d)
    返回类型:
    int                O   (1:ok,0:error)
    
    • 处理过程:
    1. 计算当前时间与ERP参数中时间的插值。
    2. 若当前时间早于ERP参数中最早的时间,采用最早的时间来计算。
    3. 若当前时间晚于ERP参数中最晚的时间,采用最晚的时间来计算。
    4. 若当前时间在ERP参数中最早与最晚的时间之间,则先找到最接近的两个时间,然后用插值。

    udstate

    void udstate(rtk_t *rtk, const obsd_t *obs, const int *sat,
                const int *iu, const int *ir, int ns, const nav_t *nav)
    
    • 所在文件:rtkpos.c
    • 功能说明:更新状态值 rtk->x 及其误差协方差 rtk->P。
    • 参数说明:
    函数参数,7个:
    rtk_t    *rtk      IO  rtk控制结构体
    obsd_t   *obs      I   观测数据
    int      sat       I   接收机和基站共同观测的卫星号列表
    int      *iu       I   接收机和基站共同观测的卫星在接收机观测值中的index值列表
    int      *ir       I   接收机和基站共同观测的卫星在基站观测值中的index值列表
    int      ns        I   接收机和基站共同观测的卫星个数
    nav_t    *nav      I   导航数据
    返回类型:
    • 处理过程:
    1. 调用 udpos 根据不同模式更新rtk中的位置、速度、加速度值和协方差。
    2. 若 电离层模式>=IONOOPT_EST,调用 udion 更新状态 rtk->x 中的电离层参数(MAXSAT个)及其协方差。
    3. 若 对流层模式>=TROPOPT_EST,调用 udtrop 更新状态 rtk->x 中的对流层参数(2或6个)及其协方差。TROPOPT_EST初始化状态 rtk->x 中的对流层参数。
    4. 若为 GLONASS AR模式,调用 udrcvbias 更新接收机硬件偏移。
    5. 若 模式>PMODE_DGPS,调用 udbias 更新载波相位偏移状态值以及其误差协方差。
    • 注意:
    1. 状态变量包含接收机位置、速度、加速度值、[每颗卫星的电离层参数]、[对流层参数]、[接收机硬件偏移]、每颗卫星的载波偏移。其中载波偏移包含周整模糊度以及小数部分,可参考 RTKLIB Manual P139 E.3.5。
    2. 过程2、3中更新状态x,只根据需要作初始化,给初值;更新时只更新协方差。

    udpos

    void udpos(rtk_t *rtk, double tt)
    
    • 所在文件:rtkpos.c
    • 功能说明:更新rtk中的位置、速度、加速度值和协方差。
    • 参数说明:
    函数参数,2个:
    rtk_t    *rtk      IO  rtk控制结构体
    double   tt        I   本次更新与上次更新的时间差
    返回类型:
    • 处理过程:
    1. 若为 PMODE_FIXED 模式,直接从选项中取得位置值给rtk->x,然后返回。
    2. 若为第一个历元,用rtk->sol中的位置值初始化rtk->x。若为dynamics模式(即需要估计速度和加速度),一并初始化。
    3. 若为 PMODE_STATIC 模式,返回。
    4. 若非dynamics模式,用rtk->sol中的位置值初始化rtk->x,然后返回。
    5. 检查位置协方差,若大于阈值VAR_POS则用rtk->sol中的位置值重置rtk->x
    6. 根据Kalman滤波的预测方程 x=Fx 和 P=FP*F+Q 更新(参考RTKLIB Manual P161 E.7.4, E.7.5)。其中更新Q时需要坐标转换。

    udbias

    void udbias(rtk_t *rtk, double tt, const obsd_t *obs, const int *sat,
                const int *iu, const int *ir, int ns, const nav_t *nav)
    
    • 所在文件:rtkpos.c
    • 功能说明:更新载波相位偏移状态值到 rtk->x 以及更新其误差协方差到 rtk->P。
    • 参数说明:
    函数参数,8个:
    rtk_t    *rtk      IO  rtk控制结构体
    double   tt        I   本次更新与上次更新的时间差
    obsd_t   *obs      I   观测数据
    int      sat       I   接收机和基站共同观测的卫星号列表
    int      *iu       I   接收机和基站共同观测的卫星在接收机观测值中的index值列表
    int      *ir       I   接收机和基站共同观测的卫星在基站观测值中的index值列表
    int      ns        I   接收机和基站共同观测的卫星个数
    nav_t    *nav      I   导航数据
    返回类型:
    • 处理过程:
    1. 对每一颗卫星,循环2~4步。
    2. 调用 detslp_ll 通过 LLI 检查接收机和基站观测数据是否有周跳。
    3. 调用 detslp_gf_L1L2 和 detslp_gf_L1L5 通过 geometry-free phase jump 检查是否有周跳。
    4. 调用 detslp_dop 通过多普勒和相位差检查接收机和基站观测数据是否有周跳。
    5. 对每一个频率,循环6~10步。
    6. 若为instantaneous AR 模式或者超出 obs outage counter,重置载波偏移值。
    7. 若检测到周跳,重置载波偏移值。
    8. 用 phase - code 的值 (ϕρλ\phi - \frac{\rho}{\lambda}) 来估计载波偏移值。
    9. 校正载波偏移值以保持载波与伪距的一致性。
    10. 更新载波偏移值及其误差协方差。

    ddres

    int ddres(rtk_t *rtk, const nav_t *nav, double dt, const double *x,
              const double *P, const int *sat, double *y, double *e,
              double *azel, const int *iu, const int *ir, int ns, double *v,
              double *H, double *R, int *vflg)
    
    • 所在文件:rtkpos.c
    • 功能说明:计算接收机或基站的双差相位/码残差(Double-Difference Residuals)
    • 参数说明:
    函数参数,16个:
    rtk_t    *rtk      IO  rtk控制结构体
    nav_t    *nav      I   导航数据
    double   dt        I   接收机和基站的时间差
    double   *x        IO  状态变量
    double   *P        IO  状态变量的误差协方差阵
    int      sat       I   接收机和基站共同观测的卫星号列表
    double   *y        IO  相位/码残差
    double   *e        IO  观测矢量 (ecef)
    double   *azel     O   方位角和俯仰角 (rad)
    int      *iu       I   接收机和基站共同观测的卫星在接收机观测值中的index值列表
    int      *ir       I   接收机和基站共同观测的卫星在基站观测值中的index值列表
    int      ns        I   接收机和基站共同观测的卫星个数
    double   *v        O   实际观测量与预测观测量的残差
    double   *H        O   观测矩阵
    double   *R        O   测量误差的协方差
    int      *vflg     O   数据有效标志
    返回类型:
    int                O   (>0:ok,0:error)
    
    • 处理过程:
    1. 一些初始化和坐标转换。
    2. 若 电离层模式>=IONOOPT_EST,调用 ionmapf 计算电离层延迟因子。
    3. 若 对流层模式>=TROPOPT_EST,调用 prectrop 计算对流层延迟因子。
    4. 大循环,对每一个频率,循环5~
    5. 寻找仰角最高的参考卫星。
    6. 小循环,对每一种导航系统,对每一颗卫星,循环7~14步,计算双差。
    7. 用传入的没有差分的相位/码残差y计算双差残差v,并计算对应的H。
    8. 若要估计电离层参数,模式IONOOPT_EST,用电离层延迟因子修正v和H。
    9. 若要估计对流层参数,模式TROPOPT_EST,用对流层延迟因子修正v和H。
    10. 用相位偏移修正v和H。
    11. 若是GLONASS系统观测值,做相关修正。
    12. 根据选项maxinno的值检测是否要排除此观测数据。
    13. 计算单差的测量误差协方差Ri、Rj。
    14. 设置数据有效标志。
    15. 若为移动基站模式PMODE_MOVEB,计算移动基站限制并设置相应的数据有效标志。
    16. 用Ri、Rj计算双差的测量误差协方差R。
    展开全文
  • matlab cov 函数解析

    万次阅读 2012-03-09 10:50:42
    最近在用matlab ,一直搞不懂cov()函数怎么算出来了。从网上查了一下,结合一些程序例子总结如下: x = 6 9 3 4 5 4 2 1 6 7 7 8 7 8 9 10 >>cov(x) ans = 0.6667 1.3333

                   最近在用matlab ,一直搞不懂cov()函数怎么算出来了。从网上查了一下,结合一些程序例子总结如下:

    x =
     
         6    9     3     4
         5    4     2     1
         6    7     7     8
         7    8     9    10
     
    >>cov(x)
     
    ans =
     
        0.6667   1.3333    2.3333    3.0000
        1.3333   4.6667    3.0000    5.0000
        2.3333   3.0000   10.9167   13.0833
        3.0000   5.0000   13.0833   16.2500


    归纳起来为:cov对角线是相应列的方差,非对角线列是相应列的协方差,你是4*4的原始方阵,所以就是4*4的矩阵。

    1.先来验证cov 对角线是相应列方差:
    >>var(x(:,1))
     
    ans =
     
        0.6667
     
    >>var(x(:,2))
     
    ans =
     
        4.6667
     
    >>var(x(:,3))
     
    ans =
     
       10.9167
     
    >>var(x(:,4))
     
    ans =
     
       16.2500


    >> diag(cov(x))
     
    ans =
     
        0.6667
        4.6667
       10.9167
       16.2500
    >> 
    从上面结果可以看出cov对角线就是每一列方差。

    2.下面来验证非对角线列是相应列的协方差
         

    D(X+Y)=D(X)+D(Y)+2COV(X,Y)
    因此 COC(X,Y)=(D(X+Y)-D(X)-D(Y))/2  (a)

    我们来验证cov(X,Y)(4,3)
    据a式
    >> (var(x(:,3)+x(:,4))-var(x(:,3))-var(x(:,4)))/2
    ans =
     
       13.0833
    >> 
    值正好等于cov(X,Y)(4,3), cov(X,Y)(3,4)。

    3.疑问

    cov(X,Y)=EXY - EX*EY

    但是我按这种方法算,不对。如下:
    >> mean(x(:,4).*x(:,3))-mean(x(:,4)).*mean(x(:,3))
    
    ans =
    
        9.8125
    ≠ 13.0833

    有大牛知道为什么不? 谢谢了。





    展开全文
  • Android 开发中基本都要使用到文件的保存和读取操作,我们一般遇见的文件读写问题有几个:文件保存在哪?以及如何考虑相关函数,目录权限的问题?...接下来就来解析一下这几个相关目录的操作函数解析

      Android 开发中基本都要使用到文件的保存和读取操作,我们一般遇见的文件读写问题有几个:文件保存在哪?以及如何使用相关函数,目录权限的问题?以及删除应用之后目录是否会随之删除的问题?接下来就来解析一下这几个相关目录的操作函数。
      转载请注明出处:http://blog.csdn.net/self_study/article/details/58587412
      对技术感兴趣的同鞋加群 544645972 一起交流。
      我们这里假设应用的名字叫做 com.android.framework:

    /data/data/package_name/ 目录

      对应的目录名字为 /data/data/com.android.framework/。
      该目录是只对应用可见的,而且如果手机没有 root,用普通权限的 adb 也看不了这个目录,该目录用来存储和应用周期相关的文件,会随着应用的卸载一起删除,相关的子目录如下所示:

    • /data/data/com.android.framework/shared_prefs/用来存储 SharedPreference,对应函数为:getSharedPreferences(String fileName, int mode);
    • /data/data/com.android.framework/databases/用来存储数据库 DB,相关函数还有 getDatabasePath();
    • /data/data/com.android.framework/app_webview 和 /data/data/com.android.framework/xxxwebviewcachexxx来存储应用内置 webview 所产生的 cache 和 cookies 等,该目录由于 android 版本不同名字和位置也可能不同;
    • /data/data/com.android.framework/lib 用来存储该应用的 .so 静态库文件;
    • /data/data/com.android.framework/cache该目录可以使用函数 getCacheDir() 获取;
    • /data/data/com.android.framework/files该目录可以使用函数 getFilesDir() 获取,openFileInput() 和 openFileOutput() 函数也是在该目录下操作文件, fileList() 函数是用来列出该 files 目录下的所有文件,deleteFile(String name) 用来删除该 files 目录下的文件;
    • /data/data/com.android.framework/XXXX这个目录下面当然也能够创建子集的目录,使用的方法就是 getDir(String name, int mode),参数中的 name 就是需要在该目录下创建的子目录名字。
    如果能够打开应用的该目录,一般会在该目录下看到很多子目录。
      PS:还有一个特别奇怪的函数,无意中看见的 getCodeCacheDir() 函数,该函数的解释为 This location is optimal for storing compiled or optimized code generated by your application at runtime,该目录适合在运行时存放应用产生的编译或者优化的代码,但是我调用就报了 NoSuchMethodError,不知道为什么,知道的可以告诉我。
      另外还有其他两个相关函数:
    getPackageCodePath() = /data/app/com.android.framework-1.apk
    getPackageResourcePath() = /data/app/com.android.framework-1.apk
      注:以上函数没有明确指定的都是使用 Context 调用。

    SD 卡下的目录

      SD 卡下的目录,顾名思义就是需要插入 SD 卡,当 SD 卡不可用时这两个目录都是无效的,SD 卡下面也分为应用的私有目录和共有目录,私有目录的生命周期也是和应用挂钩的,卸载之后就会被删除,共有目录不会随着应用的卸载而删除。
      可以通过 Environment.getExternalStorageState() 函数来获取 SD 卡的挂载状态,当该函数返回 mounted 的时候,代表 SD 卡可用。
      注意 SD 卡使用时需要注册相关权限:<uses-permission android:name=”android.permission.WRITE_EXTERNAL_STORAGE”/>。

    SD 卡私有目录

      该目录下的文件卸载应用之后会自动删除。

    /sdcard/Android/data/package_name/

      对应为 /sdcard/Android/data/com.android.framework/。

    • Android/data/com.android.framework/files/该目录可用 getExternalFilesDir(String type) 和 getExternalFilesDirs(String type) 获取,参数 type 为子目录名字,null 则为根目录,后者调用之后会自动生成该目录,并且后者返回的是一个数组,如果插入外置存储卡,外置存储卡目录也会一并返回,具体区别可以查看后面的源码和结果;
    • Android/data/com.android.framework/cache/该目录可用 getExternalCacheDir() 和 getExternalCacheDirs() 获取,后者调用之后自动生成该目录,和上面一样,后者返回的是一个数组,如果插入外置存储卡,外置存储卡目录也会一并返回,具体区别可以查看后面的源码和结果。

    /sdcard/Android/obb/package_name/

      对应为 /sdcard/Android/obb/com.android.framework/。
      需要注意的是,obb 目录也可能不存在,原文:Note if the application does not have any OBB files, this directory may not exist,一般游戏 APP 会将游戏相关的数据包放到这个目录下。
      该目录的的相关函数简单只有两个:getObbDirs()getObbDir(),前者调用之后会自动生成该目录,后者在插入外置存储卡之后会在前者的基础上另外返回外置存储卡的 obb 目录,具体的结果可以查看后面的源码和结果。

    SD 卡共有目录

      该目录下的文件卸载应用之后还会留存,所以为了 SD 卡的整洁度,不要随便在 SD 卡的根目录下面创建文件,最好以应用名字创建一个目录,所有的需要卸载之后留存或者需要给其他应用共享的文件都放到该目录下,不要在根目录下创建文件,还有一个公司的多个应用最好共享一个目录,特别讨厌 360 ,恨不得把 360XXX 式的目录全部创建完才行,极其讨厌和反对。

    /sdcard/(any_folder_name)

      对应为 /sdcard/XXX。

    • Environment.getExternalStorageState()这个函数用来获取 SD 卡的挂载状态,如果传入参数 path 则是获取该路径的的挂载状态,比如这个目录被用户的 PC 挂载,或者从设备中移除,或者其他问题发生,状态的返回是不一样的;
    • Environment.getExternalStorageDirectory()该函数用来返回 SD 卡的根目录,即 /storage/emulated/0,注意不要在根目录下创建文件,强烈建立创建一个子目录去操作,要不然会污染 SD 卡的主目录,该目录所有应用都可操作,为共享目录;
    • Environment.getDownloadCacheDirectory()该函数用来返回 SD 卡下面的下载缓存目录;
    • Environment.getDataDirectory()该函数用来获取用户的数据目录;
    • Environment.getExternalStoragePublicDirectory(String type)该函数用来根据类型返回相关目录,类型为 Environment 的一些变量,传入的类型参数不能是 null,返回的目录路径有可能不存在,所以必须在使用之前确认一下,没有就创建该目录;
    • Environment.getRootDirectory()该函数用来返回根 System 目录,只挂载为只读;

    源码及结果

      注意:手机如果还能够SD卡扩展,就相当于能挂载两张SD卡,下面的测试结果也是两张SD卡的结果:

    L.e("getDatabasePath():"+getDatabasePath("student.db"));
    L.e("getCacheDir():" + getCacheDir());
    L.e("getFilesDir():" + getFilesDir());
    String[] strings = fileList();
    for (String path : strings){//为空
       L.e("fileList():---" + path);
    }
    L.e("getDir(\"zhao\"):" + getDir("zhao", MODE_PRIVATE));
    //        L.e("getCodeCacheDir():" + getCodeCacheDir()); //java.lang.NoSuchMethodError
    L.e("getPackageCodePath():" + getPackageCodePath());
    L.e("getPackageResourcePath():" + getPackageResourcePath());
    L.e("getExternalFilesDir():" + getExternalFilesDir(null));
    File[] paths = getExternalFilesDirs(null);
    for (File path : paths){
        L.e("getExternalFilesDirs():---" + path.getPath());
    }
    L.e("getExternalCacheDir():" + getExternalCacheDir());
    paths = getExternalCacheDirs();
    for (File path : paths){
       L.e("getExternalCacheDirs():---" + path.getPath());
    }
    L.e("getObbDir():" + getObbDir());
    paths = getObbDirs();
    for (File path : paths){
       L.e("getObbDirs():---" + path.getPath());
    }
    L.e("Environment.getExternalStorageState():"+ Environment.getExternalStorageState());
    L.e("Environment.getExternalStorageDirectory():"+Environment.getExternalStorageDirectory());
    L.e("Environment.getDownloadCacheDirectory():"+Environment.getDownloadCacheDirectory());
    L.e("Environment.getExternalStoragePublicDirectory(Environment.DIRECTORY_MUSIC):"+Environment.getExternalStoragePublicDirectory(Environment.DIRECTORY_MUSIC));
    L.e("Environment.getRootDirectory():"+Environment.getRootDirectory());

    对应的结果为,注意看外部存储卡的路径:

    E/com.android.framework﹕ getDatabasePath():/data/data/com.android.framework/databases/student.db
    E/com.android.framework﹕ getCacheDir():/data/data/com.android.framework/cache
    E/com.android.framework﹕ getFilesDir():/data/data/com.android.framework/files
    E/com.android.framework﹕ getDir("zhao"):/data/data/com.android.framework/app_zhao
    E/com.android.framework﹕ getPackageCodePath():/data/app/com.android.framework-1.apk
    E/com.android.framework﹕ getPackageResourcePath():/data/app/com.android.framework-1.apk
    E/com.android.framework﹕ getExternalFilesDir():/storage/emulated/0/Android/data/com.android.framework/files
    E/com.android.framework﹕ getExternalFilesDirs():---/storage/emulated/0/Android/data/com.android.framework/files
    E/com.android.framework﹕ getExternalFilesDirs():---/storage/ext_sd/Android/data/com.android.framework/files*******
    E/com.android.framework﹕ getExternalCacheDir():/storage/emulated/0/Android/data/com.android.framework/cache
    E/com.android.framework﹕ getExternalCacheDirs():---/storage/emulated/0/Android/data/com.android.framework/cache
    E/com.android.framework﹕ getExternalCacheDirs():---/storage/ext_sd/Android/data/com.android.framework/cache*******
    E/com.android.framework﹕ getObbDir():/storage/emulated/0/Android/obb/com.android.framework
    E/com.android.framework﹕ getObbDirs():---/storage/emulated/0/Android/obb/com.android.framework
    E/com.android.framework﹕ getObbDirs():---/storage/ext_sd/Android/obb/com.android.framework*******
    E/com.android.framework﹕ Environment.getExternalStorageState():mounted
    E/com.android.framework﹕ Environment.getExternalStorageDirectory():/storage/emulated/0
    E/com.android.framework﹕ Environment.getDownloadCacheDirectory():/cache
    E/com.android.framework﹕ Environment.getExternalStoragePublicDirectory(Environment.DIRECTORY_MUSIC):/storage/emulated/0/Music
    E/com.android.framework﹕ Environment.getRootDirectory():/system
    展开全文
  • 函数式编程思想概论

    千次阅读 2019-06-26 10:57:16
    函数式编程思想概论前言函数λ 演算λ项绑定变量和自由变量约简α 变换β 约简η 变换纯函数、副作用和引用透明性函数式编程与并发编程总结 原文地址 前言 在讨论函数式编程(Functional Programming)的具体内容...


    原文地址

    前言

    在讨论函数式编程(Functional Programming)的具体内容之前,我们首先看一下函数式编程的含义。在维基百科上,函数式编程的定义如下:“函数式编程是一种编程范式。它把计算当成是数学函数的求值,从而避免改变状态和使用可变数据。它是一种声明式的编程范式,通过表达式和声明而不是语句来编程。” (见 Functional Programming

    函数式编程的思想在软件开发领域由来已久。在众多的编程范式中,函数式编程虽然出现的时间很长,但是在编程范式领域和整个开发社区中的流行度一直不温不火。函数式编程有一部分狂热的支持者,在他们眼中,函数式编程思想是解决各种软件开发问题的终极方案;而另外的一部分人,则觉得函数式编程的思想并不容易理解,学习曲线较陡,上手起来也有一定的难度。大多数人更倾向于接受面向对象或是面向过程这样的编程范式。这也是造成函数式编程范式一直停留在小众阶段的原因。

    这样两极化的反应,与函数式编程本身的特性是分不开的。函数式编程的思想脱胎于数学理论,也就是我们通常所说的λ演算( λ-calculus)。一听到数学理论,可能很多人就感觉头都大了。这的确是造成函数式编程的学习曲线较陡的一个原因。如同数学中的函数一样,函数式编程范式中的函数有独特的特性,也就是通常说的无状态或引用透明性(referential transparency)。一个函数的输出由且仅由其输入决定,同样的输入永远会产生同样的输出。这使得函数式编程在处理很多与状态相关的问题时,有着天然的优势。函数式编程的代码通常更加简洁,但是不一定易懂。函数式编程的解决方案中透露出优雅的美。

    函数式编程所涵盖的内容非常广泛,从其背后的数学理论,到其中包含的基本概念,再到诸如 Haskell 这样的函数式编程语言,以及主流编程语言中对函数式编程方式的支持,相关的专有第三方库等。通过本系列的学习,你可以了解到很多函数式编程相关的概念。你会发现很多概念都可以在日常的开发中找到相应的映射。比如做前端的开发人员一定听说过高阶组件(high-order component),它就与函数式编程中的高阶函数有着异曲同工之妙。流行的前端状态管理方案 Redux 的核心是 reduce 函数。库 reselect 则是记忆化( memoization)的精妙应用。很多 Java 开发人员已经切实的体会到了 Java 8 中的 Lambda 表达式如何让对流(Stream)的操作变得简洁又自然。

    近年来,随着多核平台和并发计算的发展,函数式编程的无状态特性,在处理这些问题时有着其他编程范式不可比拟的天然优势。不管是前端还是后端开发人员,学习一些函数式编程的思想和概念,对于手头的开发工作和以后的职业发展,都是大有裨益的。本系列虽然侧重的是 Java 平台上的函数式编程,但是对于其他背景的开发人员同样有借鉴意义。

    下面介绍函数的基本概念。

    函数

    我们先从数学中的函数开始谈起。数学中的函数是输入元素的集合到可能的输出元素的集合之间的映射关系,并且每个输入元素只能映射到一个输出元素。比如典型的函数 f(x)=xx 把所有实数的集合映射到其平方值的集合,如 f(2)=4 和 f(-2)=4。函数允许不同的输入元素映射到同一个输出元素,但是每个输入元素只能映射到一个输出元素。比如上述函数 f(x)=xx 中,2 和-2 都映射到同一个输出元素 4。这也限定了每个输入元素所对应的输出元素是固定的。每个输入元素都必须被映射到某个输出元素,也就是说函数可以应用到输入元素集合中的每个元素。

    用专业的术语来说,输入元素称为函数的参数(argument)。输出元素称为函数的值(value)。输入元素的集合称为函数的定义域(domain)。输出元素和其他附加元素的集合称为函数的到达域(codomain)。存在映射关系的输入和输出元素对的集合,称为函数的图形(graph)。输出元素的集合称为像(image)。这里需要注意像和到达域的区别。到达域还可能包含除了像中元素之外的其他元素,也就是没有输入元素与之对应的元素。

    图 1 表示了一个函数对应的映射关系(图片来源于维基百科上的 Function 条目)。输入集合 X 中的每个元素都映射到了输出集合 Y 中某个元素,即 f(1)=D、f(2)=C 和 f(3)=C。X 中的元素 2 和 3 都映射到了 Y 中的 C,这是合法的。Y 中的元素 B 和 A 没有被映射到,这也是合法的。该函数的定义域是 X,到达域是 Y,图形是 (1, D)、(2, C) 和 (3, C) 的集合,像是 C 和 D 的集合。

    图 1. 函数的映射关系
    在这里插入图片描述
    我们通常可以把函数看成是一个黑盒子,对其内部的实现一无所知。只需要清楚其映射关系,就可以正确的使用它。函数的图形是描述函数的一种方式,也就是列出来函数对应的映射中所有可能的元素对。这种描述方式用到了集合相关的理论,对于图 1 中这样的简单函数比较容易进行描述。对于包含输入变量的函数,如 f(x)=x+1,需要用到更加复杂的集合逻辑。因为输入元素 x 可以是任何数,定义域是一个无穷集合,对应的输出元素集合也是无穷的。要描述这样的函数,用我们下面介绍的 λ 演算会更加直观。

    λ 演算

    λ 演算是数理逻辑中的一个形式系统,在函数抽象和应用的基础上,使用变量绑定和替换来表达计算。讨论 λ 演算离不开形式化的表达。在本文中,我们尽量集中在与编程相关的基本概念上,而不拘泥于数学上的形式化表示。λ 演算实际上是对前面提到的函数概念的简化,方便以系统的方式来研究函数。λ 演算的函数有两个重要特征:

    • λ 演算中的函数都是匿名的,没有显式的名称。比如函数 sum(x, y) = x + y 可以写成 (x, y)|-> x + y。由于函数本身仅由其映射关系来确定,函数名称实际上并没有意义。因此使用匿名函数是合理的。
    • λ演算中的函数都只有一个输入。有多个输入的函数可以转换成多个只包含一个输入的函数的嵌套调用。这个过程就是通常所说的柯里化(currying)。如 (x, y)|-> x + y 可以转换成 x |-> (y |-> x + y)。右边的函数的返回值是另外一个函数。这一限定简化了λ演算的定义。

    对函数简化之后,就可以开始定义 λ 演算。λ 演算是基于 λ 项(λ-term)的语言。λ 项是 λ 演算的基本单元。λ 演算在 λ 项上定义了各种转换规则。

    λ项

    λ 项由下面 3 个规则来定义:

    • 一个变量 x 本身就是一个 λ 项。
    • 如果 M 是 λ 项,x 是一个变量,那么 (λx.M) 也是一个 λ 项。这样的 λ 项称为 λ 抽象(abstraction)。x 和 M 中间的点(.)用来分隔函数参数和内容。
    • 如果 M 和 N 都是 λ 项,那么 (MN) 也是一个 λ 项。这样的λ项称为应用(application)。

    所有的合法 λ 项,都只能通过重复应用上面的 3 个规则得来。需要注意的是,λ 项最外围的括号是可以省略的,也就是可以直接写为 λx.M 和 MN。当多个 λ 项连接在一起时,需要用括号来进行分隔,以确定 λ 项的解析顺序。默认的顺序是左向关联的。所以 MNO 相当于 ((MN)O)。在不出现歧义的情况下,可以省略括号。

    重复应用上述 3 个规则就可以得到所有的λ项。把变量作为λ项是重复应用规则的起点。λ 项 λx.M 定义的是匿名函数,把输入变量 x 的值替换到表达式 M 中。比如,λx.x+1 就是函数 f(x)=x+1 的 λ 抽象,其中 x 是变量,M 是 x+1。λ 项 MN 表示的是把表达式 N 应用到函数 M 上,也就是调用函数。N 可以是类似 x 这样的简单变量,也可以是 λ 抽象表示的项。当使用λ抽象时,就是我们通常所说的高阶函数的概念。

    绑定变量和自由变量

    在 λ 抽象中,如果变量 x 出现在表达式中,那么该变量被绑定。表达式中绑定变量之外的其他变量称为自由变量。我们可以用函数的方式来分别定义绑定变量(bound variable,BV)和自由变量(free variable,FV)。

    对绑定变量来说:

    • 对变量 x 来说,BV(x) = ∅。也就是说,一个单独的变量是自由的。
    • 对 λ 项 M 和变量 x 来说,BV(λx.M) = BV(M) ∪ { x }。也就是说,λ 抽象在 M 中已有的绑定变量的基础上,额外绑定了变量 x。
    • 对 λ 项 M 和λ项 N 来说,BV(MN) = BV(M) ∪ BV(N)。也就是说,λ 项的应用结果中的绑定变量的集合是各自 λ 项的绑定变量集合的并集。

    对自由变量来说,相应的定义和绑定变量是相反的:

    • 对变量 x 来说,FV(x) = { x }。
    • 对 λ M 和变量 x 来说,FV(λx.M) = FV(M) − { x }。
    • 对 λ 项 M 和 λ 项 N 来说,FV(MN) = FV(M) ∪ FV(N)。

    在 λ 项 λx.x+1 中,x 是绑定变量,没有自由变量。在 λ 项 λx.x+y 中,x 是绑定变量,y 是自由变量。

    在λ抽象中,绑定变量的名称在某些情况下是无关紧要的。如 λx.x+1 和 λy.y+1 实际上表示的是同样的函数,都是把输入值加 1。变量名称 x 或 y,并不影响函数的语义。类似 λx.x+1 和 λy.y+1 这样的 λ 项在λ演算中被认为是相等的,称为 α 等价(alpha equivalence)。

    约简

    在 λ 项上可以进行不同的约简(reduction)操作,主要有如下 3 种。

    α 变换

    α 变换(α-conversion)的目的是改变绑定变量的名称,避免名称冲突。比如,我们可以通过 α 变换把 λx.x+1 转换成 λy.y+1。如果两个λ项可以通过α变换来进行转换,则这两个 λ 项是 α 等价的。这也是我们上一节中提到的 α 等价的形式化定义。

    对 λ 抽象进行 α 变换时,只能替换那些绑定到当前 λ 抽象上的变量。如 λ 抽象 λx.λx.x 可以 α 变换为 λx.λy.y 或 λy.λx.x,但是不能变换为 λy.λx.y,因为两者的语义是不同的。λx.x 表示的是恒等函数。λx.λx.x 和 λy.λx.x 都是表示返回恒等函数的 λ 抽象,因此它们是 α 等价的。而 λx.y 表示的不再是恒等函数,因此 λy.λx.y 与 λx.λy.y 和 λy.λx.x 都不是 α 等价的。

    β 约简

    β 约简(β-reduction)与函数应用相关。在讨论 β 约简之前,需要先介绍替换的概念。对于 λ 项 M 来说,M[x := N] 表示把 λ 项 M 中变量 x 的自由出现替换成 N。具体的替换规则如下所示。A、B 和 M 是 λ 项,而 x 和 y 是变量。A ≡ B 表示两个 λ 项是相等的。

    • x[x := M] ≡ M:直接替换一个变量 x 的结果是用来进行替换的 λ 项 M。
    • y[x := M] ≡ y(x ≠ y):y 是与 x 不同的变量,因此替换 x 并不会影响 y,替换结果仍然为 y。
    • (AB)[x := M] ≡ (A[x := M]B[x := M]):A 和 B 都是 λ 项,(AB) 是 λ 项的应用。对 λ 项的应用进行替换,相当于替换之后再进行应用。
    • (λx.A)[x := M] ≡ λx.A:这条规则针对 λ 抽象。如果 x 是 λ 抽象的绑定变量,那么不需要对 x 进行替换,得到的结果与之前的 λ 抽象相同。这是因为替换只是针对 M 中 x 的自由出现,如果 x 在 M 中是不自由的,那么替换就不需要进行。
    • (λy.A)[x := M] ≡ λy.A[x := M](x ≠ y 并且 y ∉ FV(M)):这条规则也是针对λ抽象。λ 项 A 的绑定变量是 y,不同于要替换的 x,因此可以在 A 中进行替换动作。
      在进行替换之前,可能需要先使用 α 变换来改变绑定变量的名称。比如,在进行替换 (λx.y)[y := x] 时,不能直接把出现的 y 替换成 x。这样就改变了之前的 λ 抽象的语义。正确的做法是先进行 α 变换,把 λx.y 替换成 λz.y,再进行替换,得到的结果是 λz.x。

    替换的基本原则是要求在替换完成之后,原来的自由变量仍然是自由的。如果替换变量可能导致一个变量从自由变成绑定,需要首先进行 α 变换。在之前的例子中,λx.y 中的 x 是自由变量,而直接替换的结果 λx.x 把 x 变成了绑定变量,因此 α 变换是必须的。在正确的替换结果 λz.x 中,z 仍然是自由的。

    β 约简用替换来表示函数应用。对 ((λV.E) E′) 进行 β 约简的结果就是 E[V := E′]。如 ((λx.x+1)y) 进行 β 约简的结果是 (x+1)[x := y],也就是 y+1。

    η 变换

    η 变换(η-conversion)描述函数的外延性(extensionality)。外延性指的是如果两个函数当且仅当对所有参数的结果相同时,才被认为是相等的。比如一个函数 F,当参数为 x 时,它的返回值是 Fx。那么考虑声明为 λy.Fy 的函数 G。函数 G 对于输入参数 x,同样返回结果 Fx。F 和 G 可能由不同的 λ 项组成,但是只要 Fx=Gx 对所有的 x 都成立,那么 F 和 G 是相等的。

    以 F=λx.x 和 G=λx.(λy.y)x 来说,F 是恒等函数,而 G 则是在输入参数 x 上应用恒等函数。F 和 G 虽然由不同的 λ 项组成,但是它们的行为是一样,本质上都是恒等函数。我们称之为 F 和 G 是 η 等价的,F 是 G 的 η 约简,而 G 是 F 的 η 扩展。F 和 G 互为对方的 η 变换。

    纯函数、副作用和引用透明性

    了解函数式编程的人可能听说过纯函数和副作用等名称。这两个概念与引用透明性紧密相关。纯函数需要具备两个特征:

    • 对于相同的输入参数,总是返回相同的值。
    • 求值过程中不产生副作用,也就是不会对运行环境产生影响。

    对于第一个特征,如果是从数学概念上抽象出来的函数,则很容易理解。比如 f(x)=x+1 和 g(x)=x*x 这样的函数,都是典型的纯函数。如果考虑到一般编程语言中出现的方法,则函数中不能使用静态局部变量、非局部变量,可变对象的引用或 I/O 流。这是因为这些变量的值可能在不同的函数执行中发生变化,导致产生不一样的输出。第二个特征,要求函数体中不能对静态局部变量、非局部变量,可变对象的引用或 I/O 流进行修改。这就保证了函数的执行过程中不会对外部环境造成影响。纯函数的这两个特征缺一不可。下面通过几个 Java 方法来具体说明纯函数。

    在清单 1 中,方法 f1 是纯函数;方法 f2 不是纯函数,因为引用了外部变量 y;方法 f3 不是纯函数,因为使用了调用了产生副作用的 Counter 对象的 inc 方法;方法 f4 不是纯函数,因为调用 writeFile 方法会写入文件,从而对外部环境造成影响。

    清单 1. 纯函数和非纯函数示例

    int f1(int x) {
      return x + 1;
    }
     
    int f2(int x) {
      return x + y;
    }
     
    int f3(Counter c) {
      c.inc();
      return 0;
    }
     
    int f4(int x) {
      writeFile();
      return 1;
    }
    

    如果一个表达式是纯的,那么它在代码中的所有出现,都可以用它的值来代替。对于一个函数调用来说,这就意味着,对于同一个函数的输入参数相同的调用,都可以用其值来代替。这就是函数的引用透明性。引用透明性的重要性在于使得编译器可以用各种措施来对代码进行自动优化。

    函数式编程与并发编程

    随着计算机硬件多核的普及,为了尽可能地利用硬件平台的能力,并发编程显得尤为重要。与传统的命令式编程范式相比,函数式编程范式由于其天然的无状态特性,在并发编程中有着独特的优势。以 Java 平台来说,相信很多开发人员都对 Java 的多线程和并发编程有所了解。可能最直观的感受是,Java 平台的多线程和并发编程并不容易掌握。这主要是因为其中所涉及的概念太多,从 Java 内存模型,到底层原语 synchronized 和 wait/notify,再到 java.util.concurrent 包中的高级同步对象。由于并发编程的复杂性,即使是经验丰富的开发人员,也很难保证多线程代码不出现错误。很多错误只在运行时的特定情况下出现,很难排错和修复。在学习如何更好的进行并发编程的同时,我们可以从另外一个角度来看待这个问题。多线程编程的问题根源在于对共享变量的并发访问。如果这样的访问并不需要存在,那么自然就不存在多线程相关的问题。在函数式编程范式中,函数中并不存在可变的状态,也就不需要对它们的访问进行控制。这就从根本上避免了多线程的问题。

    总结

    作为 Java 函数式编程系列的第一篇文章,本文对函数式编程做了简要的概述。由于函数式编程与数学中的函数密不可分,本文首先介绍了函数的基本概念。接着对作为函数式编程理论基础的λ演算进行了详细的介绍,包括λ项、自由变量和绑定变量、α变换、β约简和η变换等重要概念。最后介绍了编程中可能会遇到的纯函数、副作用和引用透明性等概念。本系列的下一篇文章将对函数式编程中的重要概念进行介绍,包括高阶函数、闭包、递归、记忆化和柯里化等。

    展开全文
  • 解析LOOKUP函数

    千次阅读 2014-02-20 13:54:29
    返回向量或数组中的数值。...函数 LOOKUP 的数组形式在数组的第一行或第一列查找指定的数值,然后返回数组的最后一行或最后一中相同位置的数值。 函数 LOOKUP 有两种语法形式:向量和数组。   注意: LOOKUP
  • EXCEL公式获取幂函数系数解析

    千次阅读 2019-12-12 19:57:49
    创建时间:2019-07-17 工具:Excel 乘幂函数:y = a * x^(-b) 衰减系数 b = -INDEX(LINEST(LN(y_value_array),LN(x_value_array)),1) 详见:...解析: 1.函数说明 1.1 LN(y_value_array...
  • java8函数式编程

    千次阅读 2018-03-29 17:56:43
    什么是函数式编程函数式编程是java8的一大特色,也就是将函数作为一个参数传递给指定方法。别人传的要么是基本数据类型,要么就是地址引用 ,我们要穿一个“动作”。Stream说到函数式编程,就不得不提及Stream,...
  • Python函数式编程指南

    千次阅读 2014-04-14 12:38:45
    1. 函数式编程概述 1.1. 什么是函数式编程? 函数式编程使用一系列的函数解决问题。函数仅接受输入并产生输出,不包含任何能影响产生输出的内部状态。任何情况下,使用相同的参数调用函数始终能产生同样的结果。...
  • java8函数式编程实例

    千次阅读 2018-07-19 11:54:02
    什么是函数式编程 函数式编程是java8的一大特色,也就是将函数作为一个参数传递给指定方法。别人传的要么是基本数据类型,要么就是地址引用 ,我们要穿一个“动作”。 Stream 说到函数式编程,就不得不提及Stream...
  • Tensorflow操作与函数全面解析

    千次阅读 2018-09-14 16:23:46
    其中tf.mul(a, b)函数便是tf的一个基本的算数运算,接下来介绍跟多的相关函数。 2、tf函数 TensorFlow 将图形定义转换成分布式执行的操作, 以充分利用可用的计算资源(如 CPU 或 GPU。一般你不需要显式指定使用 CPU ...
  • Python函数式编程指南(二):函数

    千次阅读 2013-07-19 20:44:47
    这是此系列的第二篇,试图说明在Python中如何更好地使用函数并引导诸位使用函数式的思维进行思考。掌握并应用这些内容,就已经是至少形似的函数式风格的代码了,至于思维么,这个真靠自己。 作者水平有限,如有错漏...
  • 《Python核心编程》第十一章:函数函数式编程——介绍函数的创建、调用方式,内部函数函数装饰器、函数参数的定义和传递、函数式编程、变量作用域、闭包。
  • INDEX+SMALL+IF+ROW函数组合使用解析

    万次阅读 2014-05-24 23:48:23
    [转载]INDEX+SMALL+IF+ROW函数组合使用解析   (2014-03-08 13:03:09) 转载▼ 标签:  转载 分类: excel 原文地址:INDEX+SMALL+IF+ROW函数组合使用解析作者:EluneXY ...
  • 支持百亿数据场景,海量高性能列式数据库HiStore技术架构解析 HiStore介绍 HiStore是阿里中间件团队研发的数据库产品,是一款基于独特的知识网格技术的列式数据库,定位于海量数据高压缩比列式存储,是低存储...
  • 文章目录Python 如何修改列表中的每一个值说明写法:一、最二的写法二、标准的写法三、高级的写法四、列表解析式扩展 说明 本人Python小菜鸡,新手一枚,分享自己平时学习的笔记内容 写的不好还请见谅,欢迎大佬...
  • 上一篇博客我简单介绍了下基于金字塔的图像修复 我clone了该项目,并逐步分析,本篇文章主要讲解一下这个项目的损失函数的定义 传统基于均方误差损失 从一幅缺失的图像转换到一幅修复的图像 我们最常想到的就是MSE...
  • 导语 上一篇文章中解释了最小二乘损失函数的由来,本篇将继续向下推导,即系数WW的推导。前置知识 里面用到了几个常见的与矩阵相关的求导公式 ∂Xθ∂X=XT\frac{\partial X\theta}{\partial X}=X^T ∂θTX∂θT=...
  • 最近华人数学家黄皓,用了3页论文,云淡风轻般的证明了布尔函数敏感度的猜想。布尔函数敏感度是一个困扰学界多年的问题,甚至有很多科学家在听说这个猜想被证明后,都以为证明会非常繁锁,甚至直接打算用几周来学习...
  • 列表解析式实现筛选出大于5的数[1,2,3,4,5,6,7,8,9 list(filter(lambda x:x&gt;5,[1,2,3,4,5,6,7,8,9])) #filter函数 python 中一个高阶函数,过滤器 filter 函数接受一个函数func和一个列表,这个函数func的...
  • 列表生成(list comprehension) 生成器(generator) 迭代器(iterator) 可迭代对象(iterable) Iterable、Iterator与Generator之间的关系 语法糖 语法糖(Syntactic sugar),是由Peter J. Landin(和图灵一样的...
  • 什么是函数式编程?实用指南

    千次阅读 2021-05-15 11:17:27
    从最早的时候开始,函数式编程就已经成为软件开发中的一种潮流,但是在现代时代中,它已经具有了新的重要性。本文着眼于函数式编程背后的概念,并通过JavaScript和Java中的示例提供了实用的理解。 定义功能性程序...
  • Python内置函数

    千次阅读 多人点赞 2019-05-10 15:10:36
    Python内置函数Python abs()函数Python all() 函数Python any() 函数Python basestring() 函数Python bin() 函数Python bool() 函数Python bytearray() 函数Python callable() 函数Python chr() 函数Python ...
  • 总结了这67个pandas函数,完美解决数据处理,拿来即用!
  • Oracle SQL提供了用于执行特定操作的专用函数。这些函数大大增强了SQL语言的功能。 SQL函数的分类: 单行函数 对每一个函数应用在表的记录中时,只能输入一行结果,返回一个结果,可以出现在 SELECT 子句中和 ...
  • Word中怎么打分段函数

    千次阅读 2017-05-22 13:33:00
    分段函数也是中学时代需要学习的一种函数,它的特有特征就是自带一个大括号,然后包含了几个函数解析式,其实对于分段函数,大家应该都不陌生,那么如何在Word中打分段函数呢?下面就给大家分享具体的编辑技巧。 ...
  • caffe中最典型且常用的卷积运算,是通过将卷积操作转化成矩阵乘法来实现的,因此,卷积层的一系列程序实际上就是在为矩阵的卷积展开和矩阵乘法函数做准备,caffe_cpu_gemm也就是在调用矩阵乘法函数cblas_sgemm。...
  • SRGAN损失函数(目标函数)详解

    千次阅读 2020-02-02 21:06:05
    SRGAN损失函数解析及其Keras代码实现
  • 损失函数loss

    万次阅读 2014-04-12 16:43:39
    损失函数(loss function)是用来估量你模型的预测值f(x)与真实值Y的不一致程度,它是一个非负实值函数,通常使用L(Y, f(x))来表示,损失函数越小,模型的鲁棒性就越好。 损失函数是经验风险函数的核心部分,也是...
  • numpy函数

    千次阅读 2019-02-23 15:49:27
    1、np.c_和np.r_的用法解析 2、shape函数 3、ones函数 4、eyes函数 6、numpy.ceil() 7、linalg模块 7.1、创建矩阵 7.2、使用inv函数计算逆矩阵 7.3、求解线性方程组 7.4、np.linalg.det() 7.5、linalg.eig...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 75,700
精华内容 30,280
关键字:

如何列函数解析式