CN110703284A - 一种基于稀疏核学习的单站gnss瞬时速度和加速度构建方法 - Google Patents

一种基于稀疏核学习的单站gnss瞬时速度和加速度构建方法 Download PDF

Info

Publication number
CN110703284A
CN110703284A CN201910874925.1A CN201910874925A CN110703284A CN 110703284 A CN110703284 A CN 110703284A CN 201910874925 A CN201910874925 A CN 201910874925A CN 110703284 A CN110703284 A CN 110703284A
Authority
CN
China
Prior art keywords
vector
observation
epoch
time
kernel function
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201910874925.1A
Other languages
English (en)
Other versions
CN110703284B (zh
Inventor
常国宾
钱妮佳
陈超
张书毕
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China University of Mining and Technology CUMT
Original Assignee
China University of Mining and Technology CUMT
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China University of Mining and Technology CUMT filed Critical China University of Mining and Technology CUMT
Priority to CN201910874925.1A priority Critical patent/CN110703284B/zh
Publication of CN110703284A publication Critical patent/CN110703284A/zh
Application granted granted Critical
Publication of CN110703284B publication Critical patent/CN110703284B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
    • G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/52—Determining velocity
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/24—Acquisition or tracking or demodulation of signals transmitted by the system
    • G01S19/246—Acquisition or tracking or demodulation of signals transmitted by the system involving long acquisition integration times, extended snapshots of signals or methods specifically directed towards weak signal acquisition
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/24—Acquisition or tracking or demodulation of signals transmitted by the system
    • G01S19/27—Acquisition or tracking or demodulation of signals transmitted by the system creating, predicting or correcting ephemeris or almanac data within the receiver
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computational Mathematics (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Computing Systems (AREA)
  • Power Engineering (AREA)
  • Complex Calculations (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

本发明公开了一种基于稀疏核学习的单站GNSS瞬时速度和加速度构建方法，步骤如下：S1：采集预设时段内的GNSS载波相位观测量；S2：构建时间差分载波相位观测量及第一观测方程；S3：逐历元处理时间差分载波相位观测量，获取历元间位移和接收机钟差的增量向量、协方差矩阵和互协方差矩阵；S4：获取不同方向的位移增量序列观测向量，及其协方差矩阵；S5：获取相应方向位移的稀疏核函数模型；S6：根据稀疏核函数模型，获取瞬时速度和加速度大小。本发明在采用稀疏核学习方法进行建模的同时，不失一般性地将初始历元处的位移设为零，并据此进行消元，从而不仅减少了计算量，还提高了建模精度。

Description

一种基于稀疏核学习的单站GNSS瞬时速度和加速度构建方法
技术领域
本发明涉及GNSS卫星导航数据处理技术领域，尤其涉及一种基于稀疏核学习的单站GNSS瞬时速度和加速度构建方法。
背景技术
采用GNSS的时间差分载波相位观测量可以在单站情况下进行精密的速度测量，然而这种GNSS测速方法存在如下两点不足：
第一，该方法是一种历元间平均速度测量方法，而不是瞬时速度测量方法；
第二，该方法只能得到观测历元处的速度值，而无法直接得到其他时刻处(不在观测历元处)的速度。
对上述速度进行中心差分而得到的加速度也将存在类似的缺陷。然而在很多科学和工程领域，人们感兴趣的是载体的瞬时速度和/或瞬时加速度，而不是历元间平均速度和平均加速度。
发明内容
发明目的：针对单站情况下，通过GNSS载波相位观测量载体的瞬时速度和瞬时加速度的过程中，存在的无法测量瞬时速度和除观测历元处之外其他时刻处的速度无法直接获取的问题，本发明提出一种基于稀疏核学习的单站GNSS瞬时速度和加速度构建方法。
技术方案：为实现本发明的目的，本发明所采用的技术方案是：
一种基于稀疏核学习的单站GNSS瞬时速度和加速度构建方法，所述方法具体包括有如下步骤：
S1：采集预设时间段内的GNSS载波相位观测量；
S2：将所述GNSS载波相位观测量进行历元间差分操作，构建时间差分载波相位观测量及对应的第一观测方程；
S3：通过最小二乘法逐历元处理所述时间差分载波相位观测量，获取历元间位移和接收机钟差的增量向量及其协方差矩阵，并计算得到位移增量和接收机钟差增量组成的增量向量的估计值在相邻历元间的互协方差矩阵；
S4：从所述位移增量和接收机钟差增量组成的增量参数向量的估计值中，提取出三个不同方向的位移增量，并分别组成三个序列，所述三个序列即为相应方向的位移增量序列，同时将所述位移增量序列称为位移增量序列观测向量，将所述位移增量序列观测向量表示为y，对应的观测误差表示为e，观测误差的协方差矩阵表示为Qee，所述观测误差的协方差矩阵通过协方差矩阵和互协方差矩阵进行构建；
S5：根据所述位移增量序列观测向量y和观测误差的协方差矩阵Qee，获取相应方向位移的稀疏核函数模型；
S6：根据所述稀疏核函数模型，确定出相应方向的速度模型和加速度模型，获取瞬时速度和加速度大小。
进一步地讲，在所述步骤S2中，构建所述时间差分载波相位观测量及对应的第一观测方程，具体如下：
S2.1：根据所述预设时间段内的GNSS载波相位观测量，将前后两个历元处的载波相位观测量相减，得到所述预设时间段内所有的时间差分载波相位观测量；
S2.2：将任一历元处所有可见卫星所对应的时间差分载波相位观测量组成观测向量，并将所述观测向量表示为
Figure BDA0002204020990000021
所述观测向量对应于第一观测方程，所述第一观测方程具体为：
Figure BDA0002204020990000022
其中：
Figure BDA0002204020990000023
为第k历元处的观测向量，Hk为第k历元处的观测矩阵，hk为第k历元处的位移增量和接收机钟差增量组成的增量参数向量，εk为第k历元处的载波相位观测量的观测误差，εk-1为第k-1历元处的载波相位观测量的观测误差。
进一步地讲，在所述步骤S3中，通过最小二乘法逐历元处理所述时间差分载波相位观测量，具体为：
Figure BDA0002204020990000024
其中：
Figure BDA0002204020990000025
Figure BDA0002204020990000026
为第k历元处的位移增量和接收机钟差增量组成的增量参数向量的估计值，Hk为第k历元处的观测矩阵，Qεε,k为εk的协方差矩阵，εk为第k历元处的载波相位观测量的观测误差，Qhh,k为
Figure BDA0002204020990000031
的协方差矩阵，Qεε,k-1为εk-1的协方差矩阵，εk-1为第k-1历元处的载波相位观测量的观测误差，
Figure BDA0002204020990000032
为第k历元处的观测向量。
进一步地讲，在所述步骤S3中，所述位移增量和接收机钟差增量组成的增量参数向量的估计值在相邻历元间的互协方差矩阵的计算公式，具体为：
Figure BDA0002204020990000033
其中：
Rhh,k为
Figure BDA0002204020990000035
和
Figure BDA0002204020990000036
之间的互协方差矩阵，
Figure BDA0002204020990000037
为第k历元处的位移增量和接收机钟差增量组成的增量参数向量的估计值，
Figure BDA0002204020990000038
为第k-1历元处的位移增量和接收机钟差增量组成的增量参数向量的估计值，
Figure BDA0002204020990000039
为第k历元处的观测向量，
Figure BDA00022040209900000310
为第k-1历元处的观测向量，εk-1为第k-1历元处的观测向量的观测误差，Qεε,k-1为εk-1的协方差矩阵，Qεε,k为εk的协方差矩阵，εk为第k历元处的载观测向量
Figure BDA00022040209900000312
的观测误差，Qεε,k-2为εk-2的协方差矩阵，εk-2为第k-2历元处的观测向量
Figure BDA00022040209900000313
的观测误差，Hk为第k历元处的观测矩阵，Hk-1为第k-1历元处的观测矩阵。
进一步地讲，在所述步骤S4中，所述观测误差的协方差矩阵Qee的主对角线元素，即为所述位移增量序列观测向量对应的协方差矩阵中的对角线元素，所述观测误差的协方差矩阵Qee的次对角线元素，即为所述位移增量序列观测向量对应的互协方差矩阵中的对角线元素，所述观测误差的协方差矩阵，具体为：
其中：Qee为观测误差向量的协方差矩阵，令建模方向对应于向量hk的第j个元素，hk为第k历元处的位移增量和接收机钟差增量组成的增量参数向量，则：q1为Qhh,1的第j个对角线元素，q2为Qhh,2的第j个对角线元素，q3为Qhh,3的第j个对角线元素，qn-2为Qhh,n-2的第j个对角线元素，qn-1为Qhh,n-1的第j个对角线元素，Qhh,k为
Figure BDA0002204020990000041
的协方差矩阵，
Figure BDA0002204020990000042
为第k历元处的位移增量和接收机钟差增量组成的增量向量的估计值，r1为Rhh,1的第j个对角线元素，r2为Rhh,2的第j个对角线元素，rn-2为Rhh,n-2的第j个对角线元素，Rhh,k为
Figure BDA0002204020990000043
和
Figure BDA0002204020990000044
之间的互协方差矩阵，
Figure BDA0002204020990000045
为第k-1历元处的位移增量和接收机钟差增量组成的增量向量的估计值。
进一步地讲，在所述步骤S5中，获取相应方向位移的稀疏核函数模型，具体如下：
S5.1：预设N组超参数，从中选出第i组超参数进行使用，所述超参数包括核函数宽度参数和正则化参数；
S5.2：根据所述选出的核函数宽度参数，采用Gauss RBF核函数，将任一方向位移的时间连续函数表示为如下核函数模型，具体为：
Figure BDA0002204020990000046
其中：xt为t时刻相应方向的位移，n为观测历元的总个数，αj为第j个核函数的系数，K为核函数，t为时间，tj为第j个核函数的中心时刻；
S5.3：通过所述不同方向的位移增量序列和核函数模型，获取各方向的位移增量序列和核函数模型之间的关系，具体为：
Figure BDA0002204020990000047
其中：κkj＝K(tk,tj)-K(tk-1,tj)
dk为相应方向位移增量序列在第k历元处的值，
Figure BDA0002204020990000048
为tk时刻处相应方向的位移，
Figure BDA0002204020990000049
为tk-1时刻处相应方向的位移，n为观测历元的总个数，αj为第j个核函数的系数，tk为第k个观测历元对应的时间，tk-1为第k-1个观测历元对应的时间，tj为第j个核函数的中心时刻，K为核函数；
S5.4：根据所述各方向的位移增量序列和核函数模型之间的关系，确定出所述第二观测方程，具体为：
y＝Aα+e
其中：矩阵A的第k行第j列元素为：κkj＝K(tk,tj)-K(tk-1,tj)
y为位移增量序列观测向量，A为κkj组成的矩阵，K为核函数，tk为第k个观测历元对应的时间，tk-1为第k-1个观测历元对应的时间，tj为第j个核函数的中心时刻，α为所有核函数系数组成的向量，e为观测误差向量；
S5.5：将初始历元处的位移设置为零，确定出第一个核函数系数的求解公式，具体为：
α1＝-bTβ
其中：uj＝K(t1,tj)，u＝[u1 u2…un]T
α1为第一个核函数的系数，b为向量u中除第一个元素之外的向量部分，K为核函数，t1为初始时间，tj为第j个核函数的中心时刻，β为待求解参数；
S5.6：将所述第一个核函数系数的求解公式代入第二观测方程中，同时令A＝[cC]，确定出待求解参数的观测方程，具体为：
Figure BDA0002204020990000051
其中：κkj＝K(tk,tj)-K(tk-1,tj)，uj＝K(t1,tj)，u＝[u1 u2…un]T，B＝C-cbT
y为位移增量序列观测向量，A为κkj组成的矩阵，K为核函数，tk为第k个观测历元对应的时间，tk-1为第k-1个观测历元对应的时间，tj为第j个核函数的中心时刻，α为所有核函数系数组成的向量，e为观测误差向量，c为矩阵A中的第一列向量，C为矩阵A中除第一列向量之外的矩阵，α1为第一个核函数的系数，β为待求解参数，b为向量u中除第一个元素之外的向量部分，t1为初始时间；
S5.7：根据所述选出的正则化参数，所述待求解参数通过稀疏正则化方法进行求解，确定出所述待求解参数的大小；
S5.8：将所述确定出的待求解参数代入AIC求解公式中，所述AIC求解公式具体为：
Figure BDA0002204020990000061
其中：
Figure BDA0002204020990000062
B＝C-cbT，κkj＝K(tk,tj)-K(tk-1,tj)，uj＝K(t1,tj)，u＝[u1 u2…un]T
AICi为第i组超参数对应的AIC值，n为观测历元的总个数，y为位移增量序列观测向量，c为矩阵A中的第一列向量，C为矩阵A中除第一列向量之外的矩阵，b为向量u中除第一个元素之外的向量部分，A为κkj组成的矩阵，K为核函数，tk为第k个观测历元对应的时间，tk-1为第k-1个观测历元对应的时间，tj为第j个核函数的中心时刻，t1为初始时间，
Figure BDA0002204020990000063
为最终确定出的待求解参数，Qee为观测误差向量的协方差矩阵，n0为
Figure BDA0002204020990000064
中非零元素的个数；
S5.9：从所述预设的N组超参数中选择一组未被选择过的超参数，并重复步骤S5.1-步骤S5.8，直至所述预设的N组超参数中所有的超参数均被选择，获取每组超参数对应的AIC值，并从所有AIC值中选出最小值，所述最小AIC值对应的一组超参数即为最优超参数，同时将所述最优超参数对应的核函数模型作为最终的模型。
进一步地讲，在所述步骤S5.7中，所述最终确定出的待求解参数，即为能使如下代价函数最小时对应的待求解参数，所述代价函数具体为：
Figure BDA0002204020990000065
其中：B＝C-cbT，κkj＝K(tk,tj)-K(tk-1,tj)，uj＝K(t1,tj)，u＝[u1 u2…un]T
φ为稀疏正则化代价函数，y为位移增量序列观测向量，c为矩阵A中的第一列向量，C为矩阵A中除第一列向量之外的矩阵，A为κkj组成的矩阵，K为核函数，tk为第k个观测历元对应的时间，tk-1为第k-1个观测历元对应的时间，tj为第j个核函数的中心时刻，b为向量u中除第一个元素之外的向量部分，t1为初始时间，β为待求解参数，Qee为观测误差向量的协方差矩阵，μi为第i组正则化参数；
所述待求解参数的观测方程通过稀疏正则化方法进行求解，具体如下：
S5.7.1：将待求解参数进行初始化，具体为：
Figure BDA0002204020990000071
其中：B＝C-cbT，κkj＝K(tk,tj)-K(tk-1,tj)，uj＝K(t1,tj)，u＝[u1 u2…un]T
β0为迭代前待求解参数的初步估计值，θ1和s1均为第一次迭代运算引入的辅助变量，c为矩阵A中的第一列向量，C为矩阵A中除第一列向量之外的矩阵，A为κkj组成的矩阵，K为核函数，tk为第k历元对应的时间，tk-1为第k-1历元对应的时间，tj为第j个核函数的中心时刻，b为向量u中除第一个元素之外的向量部分，t1为初始时间，Qee为观测误差向量的协方差矩阵，μi为第i组正则化参数，I为单位矩阵，y为位移增量序列观测向量；
S5.7.2：将所述迭代前待求解参数的初步估计值进行迭代计算，直至算法收敛，确定出每次迭代过程中所述待求解参数的大小，具体为：
Figure BDA0002204020990000073
其中：
Figure BDA0002204020990000074
B＝C-cbT，κkj＝K(tk,tj)-K(tk-1,tj)，uj＝K(t1,tj)，u＝[u1 u2…un]T
βm为第m次迭代时待求解参数的估计值，sm+1和θm+1为第m+1次迭代运算引入的辅助变量，
Figure BDA0002204020990000076
为θm-2λBTQee -1(Bθ-y)的软阈值函数，θm和sm为第m次迭代运算引入的辅助变量，λmax为最大特征值，c为矩阵A中的第一列向量，C为矩阵A中除第一列向量之外的矩阵，A为κkj组成的矩阵，K为核函数，tk为第k历元对应的时间，tk-1为第k-1历元对应的时间，tj为第j个核函数的中心时刻，b为向量u中除第一个元素之外的向量部分，t1为初始时间，Qee为观测误差向量的协方差矩阵，y为位移增量序列观测向量，βm-1为第m-1次迭代时待求解参数的估计值。
进一步地讲，在所述步骤S5.7.2中，所述软阈值函数的求解过程，具体如下：
S5.7.2.1：给定向量w，对所述给定的向量进行软阈值操作，得到一个同维数的向量，所述同维数向量的第j个元素，具体为：
Figure BDA0002204020990000081
其中：
Figure BDA0002204020990000082
B＝C-cbT，κkj＝K(tk,tj)-K(tk-1,tj)，uj＝K(t1,tj)，u＝[u1 u2…un]T
w表示软阈值函数中的自变量，(|wj|-α)+为铰链函数，sgn(wj)为符号函数，θm为第m次迭代运算引入的辅助变量，λmax为最大特征值，c为矩阵A中的第一列向量，C为矩阵A中除第一列向量之外的矩阵，A为κkj组成的矩阵，K为核函数，tk为第k历元对应的时间，tk-1为第k-1历元对应的时间，tj为第j个核函数的中心时刻，b为向量u中除第一个元素之外的向量部分，t1为初始时间，Qee为观测误差向量的协方差矩阵，y为位移增量序列观测向量；
S5.7.2.2：对所述铰链函数进行求解，具体为：
Figure BDA0002204020990000084
其中：(|wj|-α)+为铰链函数，wj为w的第j个分量，w表示软阈值函数中的向量，α为所有核函数系数组成的向量；
S5.7.2.3：对所述符号函数进行求解，具体为：
Figure BDA0002204020990000085
其中：sgn(wj)为符号函数，wj为w的第j个分量，w表示软阈值函数中的向量。
进一步地讲，在所述步骤S5.7.2.1和步骤S5.8中，Qee -1(Bθm-y)和
Figure BDA0002204020990000092
均通过TMA算法进行求解，即依次执行如下公式，具体为：
Figure BDA0002204020990000094
Figure BDA0002204020990000095
其中：f＝Qee -1(Bθm-y)＝Qee -1g或
Figure BDA0002204020990000097
qk为Qhh,k的第k个对角线元素，Qhh,k为
Figure BDA0002204020990000098
的协方差矩阵，
Figure BDA0002204020990000099
为第k历元处的位移增量和接收机钟差增量组成的增量序列的估计值，rk为Rhh,k的第k个对角线元素，Rhh,k为
Figure BDA00022040209900000910
和
Figure BDA00022040209900000911
之间的互协方差矩阵，
Figure BDA00022040209900000912
为第k-1历元处的位移增量和接收机钟差增量组成的增量序列的估计值，n为历元的总个数，rk-1为Rhh,k的第k-1个对角线元素，fk为向量f的第k个元素，gk为向量g的第k个元素，和均为辅助变量。
进一步地讲，在所述步骤S6中，获取所述瞬时速度和加速度大小，具体如下：
S6.1：将所述确定出的待求解参数代入第一个核函数系数的求解公式中，获取第一个核函数系数的值；
S6.2：根据所述第一个核函数系数的值，通过所述最终的核函数模型，确定出相应方向的速度模型和加速度模型，具体为：
Figure BDA0002204020990000101
其中：
Figure BDA0002204020990000102
vt为方向分量中第t时刻对应的速度，at为方向分量中第t时刻对应的加速度，xt为方向分量中第t时刻对应的时间函数，t为时间，n为历元的总个数，αj为第j个核函数的系数，
Figure BDA0002204020990000103
为最优的核函数宽度参数，tj为第j个核函数的中心时刻；
S6.3：将任意时刻代入所述速度模型和加速度模型中，获取所述时刻对应的瞬时速度和加速度大小。
有益效果：与现有技术相比，本发明的技术方案具有以下有益技术效果：
(1)本发明在进行最小二乘估计时，除了得到向量估计值、估计值的协方差矩阵外，还将得到估计值序列在相邻历元处的互协方差，该互协方差可以反映出时间差分载波相位观测量的有色噪声特性，同时该互协方差在后续的建模过程中也会被使用，从而使得本发明所采用的随机模型更具有实用性；
(2)本发明的单站GNSS瞬时速度和加速度构建方法适用于单站场合，无需参考站等额外设施，也无需精密星厉等第三方产品，并且在采用稀疏核学习方法进行建模的同时，还不失一般性地将初始历元处的位移设为零，并据此进行消元，从而不仅减少了计算量，还提高了建模的精度；
(3)本发明的单站GNSS瞬时速度和加速度构建方法充分考虑了观测量历元间的相关性，并针对性地采用了高效的TMA算法，再进一步降低计算量的同时，还可以得到载体在任一时刻的瞬时速度和瞬时加速度，并且通过高阶微分，还可以得到除速度、加速度之外的其他量。
附图说明
图1是本发明单站GNSS瞬时速度和加速度构建方法的流程示意图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。其中，所描述的实施例是本发明一部分实施例，而不是全部的实施例。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。
实施例1
参考图1，本实施例提供了一种基于稀疏核学习的单站GNSS瞬时速度和加速度构建方法，该方法适用于单站GNSS数据，这里的单站是指单独的流动站，不需要参考站。同时该方法可以采用广播星厉，而不需要第三方提供的精密星厉，从而也不需要额外的设备和信息辅助。该方法具体包括如下步骤：
步骤S1：采集预设时间段内的GNSS载波相位观测量。在本实施例中，预设时间段可以是航空重力测量中的一个航线时间，也可以将一个航线时间划分为N个时间段，其中划分的任一时间段均可作为预设时间段。同时其中的GNSS可以为GPS、Beidou、Galileo等任一码分多址卫星导航系统。
步骤S2：将步骤S1采集的GNSS载波相位观测量进行历元间差分操作，从而构建时间差分载波相位观测量及其对应的第一观测方程。其中进行历元间差分操作，可以有效消除大部分历元间共模误差或时间慢变误差，譬如：电离层延迟误差、对流层延迟误差、卫星轨道误差、卫星钟差、多路径误差。在本实施例中，构建时间差分载波相位观测量，具体如下：
步骤S2.1：根据预设时间段内的GNSS载波相位观测量，将其中前后两个历元处的载波相位观测量都进行相减，从而可以得到预设时间段内所有的时间差分载波相位观测量。
步骤S2.2：将预设时间段内所有的时间差分载波相位观测量组成观测向量，并标记为
Figure BDA0002204020990000111
其中k表示历元的序号，n表示历元的总个数。
同时通过时间差分载波相位观测量组建观测向量，由于在时间差分的过程中，可以消去大部分历元间共模误差或时间慢变误差，从而可以通过第一观测方程获取观测向量，该第一观测方程具体为：
Figure BDA0002204020990000112
其中：
Figure BDA0002204020990000113
为第k历元处的观测向量，Hk为第k历元处的观测矩阵，hk为第k历元处的位移增量和接收机钟差增量组成的增量参数向量，εk为第k历元处的载波相位观测量的观测误差，εk-1为第k-1历元处的载波相位观测量的观测误差。
具体地讲，第k历元处的观测矩阵Hk是根据标准单点定位的结果计算得到的，同时第k历元处的观测向量
Figure BDA0002204020990000121
即为所有可见卫星所对应的时间差分载波相位观测量所组成的观测向量。
步骤S3：通过最小二乘法对步骤S2.2中获取得到的时间差分载波相位观测量逐历元进行处理，从而获取历元间位移增量向量，其中包括有历元间位移增量的估计值序列和估计值向量的协方差矩阵。同时在采用最小二乘进行数据处理时，除了计算常规的参数估计及其协方差矩阵外，还要计算相邻历元间的互协方差矩阵。
在本实施例中，通过最小二乘法逐历元处理时间差分载波相位观测量，具体为：
Figure BDA0002204020990000122
其中：
Figure BDA0002204020990000123
Figure BDA0002204020990000124
为第k历元处的位移增量和接收机钟差增量组成的增量参数向量的估计值，Hk为第k历元处的观测矩阵，Qεε,k为εk的协方差矩阵，εk为第k历元处的载波相位观测量的观测误差，Qhh,k为
Figure BDA0002204020990000125
的协方差矩阵，Qεε,k-1为εk-1的协方差矩阵，εk-1为第k-1历元处的载波相位观测量的观测误差，
Figure BDA0002204020990000126
为第k历元处的观测向量。
具体地讲，第k历元处的观测向量
Figure BDA0002204020990000127
即为所有可见卫星所对应的时间差分载波相位观测量所组成的观测向量。
同时还计算得到第k历元处的位移增量和接收机钟差增量组成的增量序列hk的估计值
Figure BDA0002204020990000128
第k-1历元处的位移增量和接收机钟差增量组成的增量序列hk-1的估计值
Figure BDA0002204020990000129
之间的互协方差矩阵，该互协方差矩阵来源与步骤S2.2中历元间差分操作中引入的历元间相关性。
在本实施例中，互协方差矩阵的计算公式具体为：
Figure BDA0002204020990000131
其中：
Figure BDA0002204020990000132
Rhh,k为
Figure BDA0002204020990000133
和之间的互协方差矩阵，
Figure BDA0002204020990000135
为第k历元处的位移增量和接收机钟差增量组成的增量参数向量的估计值，
Figure BDA0002204020990000136
为第k-1历元处的位移增量和接收机钟差增量组成的增量参数向量的估计值，为第k历元处的观测向量，
Figure BDA0002204020990000138
为第k-1历元处的观测向量，εk-1为第k-1历元处的观测向量
Figure BDA0002204020990000139
的观测误差，Qεε,k-1为εk-1的协方差矩阵，Qεε,k为εk的协方差矩阵，εk为第k历元处的载观测向量
Figure BDA00022040209900001310
的观测误差，Qεε,k-2为εk-2的协方差矩阵，εk-2为第k-2历元处的观测向量
Figure BDA00022040209900001311
的观测误差，Hk为第k历元处的观测矩阵，Hk-1为第k-1历元处的观测矩阵。
步骤S4：从步骤S2中的位移增量和接收机钟差增量组成的增量参数向量的估计值中，提取三个不同方向的位移增量，并分别组成三个序列，该三个序列分别为相应方向的位移增量序列，即东、北、天三个方向的位移增量序列，并分别标记为：其中
Figure BDA00022040209900001313
表示东方向的位移增量序列，
Figure BDA00022040209900001314
表示北方向的位移增量序列，表示天方向的位移增量序列。从而舍去位置增量或位移增量和接收机钟差增量组成的增量序列中的第四组变量，也就是接收机钟差参数。同时将每个方向的位移增量序列作为该方向的位移增量序列观测向量。
在本实施例中，将位移增量序列观测向量表示为y，对应的观测误差表示为e，观测误差的协方差矩阵表示为Qee。同时观测误差的协方差矩阵、位移增量序列观测向量对应的观测误差之间满足下列关系式，具体为：
Qee＝cov[e]
其中：Qee为观测误差向量的协方差矩阵，e为观测误差向量。
将每个方向的位移增量序列均代入步骤S3中，通过协方差矩阵计算公式获取每个方向的位移增量序列所对应的协方差矩阵，通过互协方差矩阵计算公式提取每个方向的位移增量序列中历元间的互协方差矩阵。其中，计算得到的协方差矩阵和互协方差矩阵均为3×3大小的矩阵。
将位移增量序列观测向量对应的协方差矩阵中的对角线元素，作为观测误差的协方差矩阵中的主对角线元素，将位移增量序列观测向量对应的互协方差矩阵中的对角线元素，作为观测误差的协方差矩阵中的次对角线元素。其中观测误差的协方差矩阵为对称三对角矩阵，同时该观测误差的协方差矩阵，具体为：
Figure BDA0002204020990000141
其中：Qee为观测误差向量的协方差矩阵，令建模方向对应于向量hk的第j个元素，hk为第k历元处的位移增量和接收机钟差增量组成的增量参数向量，则：q1为Qhh,1的第j个对角线元素，q2为Qhh,2的第j个对角线元素，q3为Qhh,3的第j个对角线元素，qn-2为Qhh,n-2的第j个对角线元素，qn-1为Qhh,n-1的第j个对角线元素，Qhh,k为
Figure BDA0002204020990000142
的协方差矩阵，
Figure BDA0002204020990000143
为第k历元处的位移增量和接收机钟差增量组成的增量向量的估计值，r1为Rhh,1的第j个对角线元素，r2为Rhh,2的第j个对角线元素，rn-2为Rhh,n-2的第j个对角线元素，Rhh,k为
Figure BDA0002204020990000144
和
Figure BDA0002204020990000145
之间的互协方差矩阵，
Figure BDA0002204020990000146
为第k-1历元处的位移增量和接收机钟差增量组成的增量向量的估计值。
步骤S5：根据位移增量序列观测向量y和观测误差的协方差矩阵Qee，获取相应方向位移的稀疏核函数模型。具体如下：
步骤S5.1：预设N组超参数，从中选出第i组超参数进行使用，其中超参数包括有核函数宽度参数σ和正则化参数μ。
步骤S5.2：根据步骤S5.1中选择的核函数宽度参数σ，采用Gauss RBF核函数，将任一方向位移的时间连续函数均表示为如下核函数模型，具体为：
Figure BDA0002204020990000151
其中：xt为t时刻相应方向的位移，n为观测历元的总个数，αj为第j个核函数的系数，K为核函数，t为时间，tj为第j个核函数的中心时刻。
在本实施例中，由于选择了东、北、天三个不同方向的位移增量序列，从而核函数模型也分为东、北、天三个不同方向的核函数模型，并分别标记为：
Figure BDA0002204020990000152
其中
Figure BDA0002204020990000153
表示东方向对应的核函数模型，
Figure BDA0002204020990000154
为北方向对应的核函数模型，
Figure BDA0002204020990000155
为天方向对应的核函数模型。同时第j个核函数的系数αj为待求解的参数，第j个核函数的中心时刻tj为第k个观测历元的时刻。
同时
K(t,tj)＝exp[-σi -2(t-tj)]
其中：K为核函数，σi为第i组核函数宽度参数，t为时间，tj为第j个核函数的中心时刻。
步骤S5.3：根据步骤S4中的三个不同方向的位移增量序列
Figure BDA0002204020990000156
步骤S5.2中三个不同方向的核函数模型
Figure BDA0002204020990000157
确定出各方向的位移增量序列和该方向的核函数模型之间的关系，具体为：
Figure BDA0002204020990000158
其中：κkj＝K(tk,tj)-K(tk-1,tj)
dk为相应方向位移增量序列在第k历元处的值，
Figure BDA0002204020990000159
为tk时刻处相应方向的位移，
Figure BDA00022040209900001510
为tk-1时刻处相应方向的位移，n为观测历元的总个数，αj为第j个核函数的系数，tk为第k个观测历元对应的时间，tk-1为第k-1个观测历元对应的时间，tj为第j个核函数的中心时刻，K为核函数。
步骤S5.4：根据各方向的位移增量序列和核函数模型之间的关系，确定出第二观测方程，具体为：
y＝Aα+e
其中：κkj＝K(tk,tj)-K(tk-1,tj)
y为位移增量序列观测向量，A为κkj组成的矩阵，K为核函数，tk为第k个观测历元对应的时间，tk-1为第k-1个观测历元对应的时间，tj为第j个核函数的中心时刻，α为所有核函数系数组成的向量，e为观测误差向量。
在本实施例中，所有核函数系数组成的向量α，具体为：
α＝[α1 α2…αn]T
其中：α1为第一个核函数的系数，α2为第二个核函数的系数，αn为第n个核函数的系数。
S5.5：在载波相位时间差分方法中缺乏初值信息，我们无法得到绝对位置，然而在速度和加速度测量中，绝对位置是不必要的。从而将初始历元处的位移设置为零，并根据该条件对要求解的核函数的系数进行转化，则待求解的核函数的系数的个数将从n个系数转化为n-1个系数。
将初始历元处的位移设置为零并不失一般性，因为在速度和加速度的估计中，初始位移的具体数值不会产生影响。且初始历元处位移为零是一个约束条件，根据该约束条件可以对问题进行消元处理，即将原来的n个待估参数转化为n-1个待估参数，这样做的好处在于：待估参数的减少意味着建模精度的提高。
在本实施例中，具体以东方向对应的核函数模型
Figure BDA0002204020990000161
为例进行具体描述，具体为：
Figure BDA0002204020990000162
其中：uj＝K(t1,tj)，u＝[u1 u2…un]T
为东向分量中初始时刻对应的时间函数，n为观测历元的总个数，αj为第j个核函数的系数，K为核函数，t1为初始时间，tj为第j个核函数的中心时刻。
令
Figure BDA0002204020990000164
则通过上式可以确定出第一个核函数系数α1的求解公式，具体为：
α1＝-bTβ
其中：uj＝K(t1,tj)，u＝[u1 u2…un]T
α1为第一个核函数的系数，b为向量u中除第一个元素之外的向量部分，K为核函数，t1为初始时间，tj为第j个核函数的中心时刻，β为待求解参数。
从而在确定了待求解参数β的大小之后，即可确定出第一个核函数的系数α1的大小，进而能够得到向量u。
在本实施例中，待求解参数β具体为：
β＝[α2 α3…αn]T
其中：β为待求解参数，α2为第二个核函数的系数，α3为第三个核函数的系数，αn为第n个核函数的系数。
步骤S5.6：将第一个核函数系数的求解公式代入第二观测方程中，并令A＝[c C]，从而确定出待求解参数的观测方程，具体为：
Figure BDA0002204020990000171
其中：κkj＝K(tk,tj)-K(tk-1,tj)，uj＝K(t1,tj)，u＝[u1 u2…un]T，B＝C-cbT
y为位移增量序列观测向量，A为κkj组成的矩阵，K为核函数，tk为第k个观测历元对应的时间，tk-1为第k-1个观测历元对应的时间，tj为第j个核函数的中心时刻，α为所有核函数系数组成的向量，e为观测误差向量，c为矩阵A中的第一列向量，C为矩阵A中除第一列向量之外的矩阵，α1为第一个核函数的系数，β为待求解参数，b为向量u中除第一个元素之外的向量部分，t1为初始时间。
步骤S5.7：根据步骤S5.1中选出的正则化参数，待求解参数的观测方程通过稀疏正则化方法进行求解，从而确定出待求解参数的大小。
在本实施例中，最终确定出的待求解参数，即为能使如下代价函数最小时对应的待求解参数，该代价函数具体为：
Figure BDA0002204020990000172
其中：B＝C-cbT，κkj＝K(tk,tj)-K(tk-1,tj)，uj＝K(t1,tj)，u＝[u1 u2…un]T
φ为稀疏正则化代价函数，y为位移增量序列观测向量，c为矩阵A中的第一列向量，C为矩阵A中除第一列向量之外的矩阵，A为κkj组成的矩阵，K为核函数，tk为第k个观测历元对应的时间，tk-1为第k-1个观测历元对应的时间，tj为第j个核函数的中心时刻，b为向量u中除第一个元素之外的向量部分，t1为初始时间，β为待求解参数，Qee为观测误差向量的协方差矩阵，μi为第i组正则化参数
具体地讲，待求解参数的观测方程通过稀疏正则化方法进行求解，具体如下：
步骤S5.7.1：将待求解参数进行初始化，具体为：
Figure BDA0002204020990000181
其中：B＝C-cbT，κkj＝K(tk,tj)-K(tk-1,tj)，uj＝K(t1,tj)，u＝[u1 u2…un]T
β0为迭代前待求解参数的初步估计值，θ1和s1均为第一次迭代运算引入的辅助变量，c为矩阵A中的第一列向量，C为矩阵A中除第一列向量之外的矩阵，A为κkj组成的矩阵，K为核函数，tk为第k历元对应的时间，tk-1为第k-1历元对应的时间，tj为第j个核函数的中心时刻，b为向量u中除第一个元素之外的向量部分，t1为初始时间，Qee为观测误差向量的协方差矩阵，μi为第i组正则化参数，I为单位矩阵，y为位移增量序列观测向量；
步骤S5.7.2：将迭代前待求解参数的初步估计值进行迭代计算，直至算法收敛，从而确定出第k次迭代时待求解参数的大小，具体为：
Figure BDA0002204020990000182
其中：
Figure BDA0002204020990000183
B＝C-cbT，κkj＝K(tk,tj)-K(tk-1,tj)，uj＝K(t1,tj)，u＝[u1 u2…un]T
βm为第m次迭代时待求解参数的估计值，sm+1和θm+1为第m+1次迭代运算引入的辅助变量，
Figure BDA0002204020990000192
为θm-2λBTQee -1(Bθm-y)的软阈值函数，θm和sm为第m次迭代运算引入的辅助变量，λmax为最大特征值，c为矩阵A中的第一列向量，C为矩阵A中除第一列向量之外的矩阵，A为κkj组成的矩阵，K为核函数，tk为第k历元对应的时间，tk-1为第k-1历元对应的时间，tj为第j个核函数的中心时刻，b为向量u中除第一个元素之外的向量部分，t1为初始时间，Qee为观测误差向量的协方差矩阵，y为位移增量序列观测向量，βm-1为第m-1次迭代时待求解参数的估计值。
具体地讲，软阈值函数
Figure BDA0002204020990000195
的求解过程，具体如下：
步骤S5.7.2.1：给定向量w，对该给定的向量进行软阈值操作，获取得到一个同维数的向量，该同维数向量的第j个元素，具体为：
Figure BDA0002204020990000196
其中：B＝C-cbT，κkj＝K(tk,tj)-K(tk-1,tj)，uj＝K(t1,tj)，u＝[u1 u2…un]T
w表示软阈值函数中的自变量，(|wj|-α)+为铰链函数，sgn(wj)为符号函数，θm为第m次迭代运算引入的辅助变量，λmax为最大特征值，c为矩阵A中的第一列向量，C为矩阵A中除第一列向量之外的矩阵，A为κkj组成的矩阵，K为核函数，tk为第k历元对应的时间，tk-1为第k-1历元对应的时间，tj为第j个核函数的中心时刻，b为向量u中除第一个元素之外的向量部分，t1为初始时间，Qee为观测误差向量的协方差矩阵，y为位移增量序列观测向量。
步骤S5.7.2.2：对铰链函数进行求解，具体为：
Figure BDA0002204020990000201
其中：(|wj|-α)+为铰链函数，wj为w的第j个分量，w表示软阈值函数中的向量，α为所有核函数系数组成的向量。
步骤S5.7.2.3：对符号函数进行求解，具体为：
Figure BDA0002204020990000202
其中：sgn(wj)为符号函数，wj为w的第j个分量，w表示软阈值函数中的向量。
步骤S5.8：将最终确定出的待求解参数代入AIC求解公式中，该AIC求解公式具体为：
Figure BDA0002204020990000203
其中：B＝C-cbT，κkj＝K(tk,tj)-K(tk-1,tj)，uj＝K(t1,tj)，u＝[u1 u2…un]T
AICi为第i组超参数对应的AIC值，n为观测历元的总个数，y为位移增量序列观测向量，c为矩阵A中的第一列向量，C为矩阵A中除第一列向量之外的矩阵，b为向量u中除第一个元素之外的向量部分，A为κkj组成的矩阵，K为核函数，tk为第k个观测历元对应的时间，tk-1为第k-1个观测历元对应的时间，tj为第j个核函数的中心时刻，t1为初始时间，
Figure BDA0002204020990000205
为最终确定出的待求解参数，Qee为观测误差向量的协方差矩阵，n0为中非零元素的个数。
在本实施例中，具体地讲，
Figure BDA0002204020990000207
中的
Figure BDA0002204020990000208
软阈值函数
Figure BDA0002204020990000209
中的Qee -1(Bθm-y)均通过TMA算法进行求解，即依次执行如下公式，具体为：
Figure BDA0002204020990000211
Figure BDA0002204020990000212
Figure BDA0002204020990000213
其中：f＝Qee -1(Bθm-y)＝Qee -1g或
Figure BDA0002204020990000215
qk为Qhh,k的第k个对角线元素，Qhh,k为
Figure BDA0002204020990000216
的协方差矩阵，
Figure BDA0002204020990000217
为第k历元处的位移增量和接收机钟差增量组成的增量序列的估计值，rk为Rhh,k的第k个对角线元素，Rhh,k为和
Figure BDA0002204020990000219
之间的互协方差矩阵，
Figure BDA00022040209900002110
为第k-1历元处的位移增量和接收机钟差增量组成的增量序列的估计值，n为历元的总个数，rk-1为Rhh,k的第k-1个对角线元素，fk为向量f的第k个元素，gk为向量g的第k个元素，
Figure BDA00022040209900002111
和
Figure BDA00022040209900002112
均为辅助变量。
步骤S5.9：从预设的N组超参数中选择一组未被选择过的超参数，并重复步骤S5.1-步骤S5.8，直至预设的N组超参数中所有的超参数均被选择，同时得到每组超参数对应的AIC值，并将所有的AIC值进行比较，从中选出最小值，该最小AIC值对应的一组超参数即为最优超参数，同时该最优超参数对应的核函数模型即为最终的模型。
步骤S6：通过稀疏核函数模型，确定出速度模型和加速度模型，获取瞬时速度和加速度大小。具体如下：
步骤S6.1：将最终的待求解参数
Figure BDA00022040209900002113
代入第一个核函数的系数α1的求解公式中，确定出第一个核函数的系数α1的大小。
步骤S6.2：根据第一个核函数的系数α1的大小，通过最终的核函数模型，确定出速度模型和加速度模型，在本实施例中，以东向做具体说明，具体为：
Figure BDA0002204020990000221
其中：
Figure BDA0002204020990000222
vt为方向分量中第t时刻对应的速度，at为方向分量中第t时刻对应的加速度，xt为方向分量中第t时刻对应的时间函数，t为时间，n为历元的总个数，αj为第j个核函数的系数，
Figure BDA0002204020990000223
为最优的核函数宽度参数，tj为第j个核函数的中心时刻。
步骤S6.3：将任一时刻t代入步骤S6.2中的速度模型和加速度模型中，从而即可获取得到该历元处载体在该时刻的瞬时速度和瞬时加速度。同时还可以进行更高阶地求导，得到加加速度等其他量。
以上示意性的对本发明及其实施方式进行了描述，该描述没有限制性，附图中所示的也只是本发明的实施方式之一，实际的结构和方法并不局限于此。所以，如果本领域的普通技术人员受其启示，在不脱离本发明创造宗旨的情况下，不经创造性的设计出与该技术方案相似的结构方式及实施例，均属于本发明的保护范围。

Claims (10)

1.一种基于稀疏核学习的单站GNSS瞬时速度和加速度构建方法，其特征在于，所述方法具体包括有如下步骤：
S1：采集预设时间段内的GNSS载波相位观测量；
S2：将所述GNSS载波相位观测量进行历元间差分操作，构建时间差分载波相位观测量及对应的第一观测方程；
S3：通过最小二乘法逐历元处理所述时间差分载波相位观测量，获取历元间位移和接收机钟差的增量向量及其协方差矩阵，并计算得到位移增量和接收机钟差增量组成的增量向量的估计值在相邻历元间的互协方差矩阵；
S4：从所述位移增量和接收机钟差增量组成的增量参数向量的估计值中，提取出三个不同方向的位移增量，并分别组成三个序列，所述三个序列即为相应方向的位移增量序列，同时将所述位移增量序列称为位移增量序列观测向量，将所述位移增量序列观测向量表示为y，对应的观测误差表示为e，观测误差的协方差矩阵表示为Qee，所述观测误差的协方差矩阵通过协方差矩阵和互协方差矩阵进行构建；
S5：根据所述位移增量序列观测向量y和观测误差的协方差矩阵Qee，获取相应方向位移的稀疏核函数模型；
S6：根据所述稀疏核函数模型，确定出相应方向的速度模型和加速度模型，获取瞬时速度和加速度大小。
2.根据权利要求1所述的一种基于稀疏核学习的单站GNSS瞬时速度和加速度构建方法，其特征在于，在所述步骤S2中，构建所述时间差分载波相位观测量及对应的第一观测方程，具体如下：
S2.1：根据所述预设时间段内的GNSS载波相位观测量，将前后两个历元处的载波相位观测量相减，得到所述预设时间段内所有的时间差分载波相位观测量；
S2.2：将任一历元处所有可见卫星所对应的时间差分载波相位观测量组成观测向量，并将所述观测向量表示为
Figure FDA0002204020980000011
所述观测向量对应于第一观测方程，所述第一观测方程具体为：
Figure FDA0002204020980000012
其中：为第k历元处的观测向量，Hk为第k历元处的观测矩阵，hk为第k历元处的位移增量和接收机钟差增量组成的增量参数向量，εk为第k历元处的载波相位观测量的观测误差，εk-1为第k-1历元处的载波相位观测量的观测误差。
3.根据权利要求1或2所述的一种基于稀疏核学习的单站GNSS瞬时速度和加速度构建方法，其特征在于，在所述步骤S3中，通过最小二乘法逐历元处理所述时间差分载波相位观测量，具体为：
其中：
Figure FDA0002204020980000022
为第k历元处的位移增量和接收机钟差增量组成的增量参数向量的估计值，Hk为第k历元处的观测矩阵，Qεε,k为εk的协方差矩阵，εk为第k历元处的载波相位观测量的观测误差，Qhh,k为
Figure FDA0002204020980000024
的协方差矩阵，Qεε,k-1为εk-1的协方差矩阵，εk-1为第k-1历元处的载波相位观测量的观测误差，
Figure FDA0002204020980000025
为第k历元处的观测向量。
4.根据权利要求1或2所述的一种基于稀疏核学习的单站GNSS瞬时速度和加速度构建方法，其特征在于，在所述步骤S3中，所述位移增量和接收机钟差增量组成的增量参数向量的估计值在相邻历元间的互协方差矩阵的计算公式，具体为：
Figure FDA0002204020980000026
其中：
Rhh,k为
Figure FDA0002204020980000028
和
Figure FDA0002204020980000029
之间的互协方差矩阵，
Figure FDA00022040209800000210
为第k历元处的位移增量和接收机钟差增量组成的增量参数向量的估计值，
Figure FDA00022040209800000211
为第k-1历元处的位移增量和接收机钟差增量组成的增量参数向量的估计值，
Figure FDA00022040209800000212
为第k历元处的观测向量，
Figure FDA00022040209800000213
为第k-1历元处的观测向量，εk-1为第k-1历元处的观测向量的观测误差，Qεε,k-1为εk-1的协方差矩阵，Qεε,k为εk的协方差矩阵，εk为第k历元处的载观测向量
Figure FDA00022040209800000215
的观测误差，Qεε,k-2为εk-2的协方差矩阵，εk-2为第k-2历元处的观测向量
Figure FDA0002204020980000031
的观测误差，Hk为第k历元处的观测矩阵，Hk-1为第k-1历元处的观测矩阵。
5.根据权利要求4所述的一种基于稀疏核学习的单站GNSS瞬时速度和加速度构建方法，其特征在于，在所述步骤S4中，所述观测误差的协方差矩阵Qee的主对角线元素，即为所述位移增量序列观测向量对应的协方差矩阵中的对角线元素，所述观测误差的协方差矩阵Qee的次对角线元素，即为所述位移增量序列观测向量对应的互协方差矩阵中的对角线元素，所述观测误差的协方差矩阵，具体为：
Figure FDA0002204020980000032
其中：Qee为观测误差向量的协方差矩阵，令建模方向对应于向量hk的第j个元素，hk为第k历元处的位移增量和接收机钟差增量组成的增量参数向量，则：q1为Qhh,1的第j个对角线元素，q2为Qhh,2的第j个对角线元素，q3为Qhh,3的第j个对角线元素，qn-2为Qhh,n-2的第j个对角线元素，qn-1为Qhh,n-1的第j个对角线元素，Qhh,k为
Figure FDA0002204020980000033
的协方差矩阵，
Figure FDA0002204020980000034
为第k历元处的位移增量和接收机钟差增量组成的增量向量的估计值，r1为Rhh,1的第j个对角线元素，r2为Rhh,2的第j个对角线元素，rn-2为Rhh,n-2的第j个对角线元素，Rhh,k为
Figure FDA0002204020980000035
和
Figure FDA0002204020980000036
之间的互协方差矩阵，
Figure FDA0002204020980000037
为第k-1历元处的位移增量和接收机钟差增量组成的增量向量的估计值。
6.根据权利要求5所述的一种基于稀疏核学习的单站GNSS瞬时速度和加速度构建方法，其特征在于，在所述步骤S5中，获取相应方向位移的稀疏核函数模型，具体如下：
S5.1：预设N组超参数，从中选出第i组超参数进行使用，所述超参数包括核函数宽度参数和正则化参数；
S5.2：根据所述选出的核函数宽度参数，采用Gauss RBF核函数，将任一方向位移的时间连续函数表示为如下核函数模型，具体为：
Figure FDA0002204020980000041
其中：xt为t时刻相应方向的位移，n为观测历元的总个数，αj为第j个核函数的系数，K为核函数，t为时间，tj为第j个核函数的中心时刻；
S5.3：通过所述不同方向的位移增量序列和核函数模型，获取各方向的位移增量序列和核函数模型之间的关系，具体为：
其中：κkj＝K(tk,tj)-K(tk-1,tj)
dk为相应方向位移增量序列在第k历元处的值，
Figure FDA0002204020980000043
为tk时刻处相应方向的位移，
Figure FDA0002204020980000044
为tk-1时刻处相应方向的位移，n为观测历元的总个数，αj为第j个核函数的系数，tk为第k个观测历元对应的时间，tk-1为第k-1个观测历元对应的时间，tj为第j个核函数的中心时刻，K为核函数；
S5.4：根据所述各方向的位移增量序列和核函数模型之间的关系，确定出所述第二观测方程，具体为：
y＝Aα+e
其中：矩阵A的第k行第j列元素为：κkj＝K(tk,tj)-K(tk-1,tj)
y为位移增量序列观测向量，A为κkj组成的矩阵，K为核函数，tk为第k个观测历元对应的时间，tk-1为第k-1个观测历元对应的时间，tj为第j个核函数的中心时刻，α为所有核函数系数组成的向量，e为观测误差向量；
S5.5：将初始历元处的位移设置为零，确定出第一个核函数系数的求解公式，具体为：
α1＝-bTβ
其中：uj＝K(t1,tj)，u＝[u1 u2 … un]T
α1为第一个核函数的系数，b为向量u中除第一个元素之外的向量部分，K为核函数，t1为初始时间，tj为第j个核函数的中心时刻，β为待求解参数；
S5.6：将所述第一个核函数系数的求解公式代入第二观测方程中，同时令A＝[c C]，确定出待求解参数的观测方程，具体为：
Figure FDA0002204020980000051
其中：κkj＝K(tk,tj)-K(tk-1,tj)，uj＝K(t1,tj)，u＝[u1 u2 … un]T，B＝C-cbT
y为位移增量序列观测向量，A为κkj组成的矩阵，K为核函数，tk为第k个观测历元对应的时间，tk-1为第k-1个观测历元对应的时间，tj为第j个核函数的中心时刻，α为所有核函数系数组成的向量，e为观测误差向量，c为矩阵A中的第一列向量，C为矩阵A中除第一列向量之外的矩阵，α1为第一个核函数的系数，β为待求解参数，b为向量u中除第一个元素之外的向量部分，t1为初始时间；
S5.7：根据所述选出的正则化参数，所述待求解参数通过稀疏正则化方法进行求解，确定出所述待求解参数的大小；
S5.8：将所述确定出的待求解参数代入AIC求解公式中，所述AIC求解公式具体为：
Figure FDA0002204020980000052
其中：
Figure FDA0002204020980000053
B＝C-cbT，κkj＝K(tk,tj)-K(tk-1,tj)，uj＝K(t1,tj)，u＝[u1 u2 … un]T
AICi为第i组超参数对应的AIC值，n为观测历元的总个数，y为位移增量序列观测向量，c为矩阵A中的第一列向量，C为矩阵A中除第一列向量之外的矩阵，b为向量u中除第一个元素之外的向量部分，A为κkj组成的矩阵，K为核函数，tk为第k个观测历元对应的时间，tk-1为第k-1个观测历元对应的时间，tj为第j个核函数的中心时刻，t1为初始时间，
Figure FDA0002204020980000054
为最终确定出的待求解参数，Qee为观测误差向量的协方差矩阵，n0为中非零元素的个数；
S5.9：从所述预设的N组超参数中选择一组未被选择过的超参数，并重复步骤S5.1-步骤S5.8，直至所述预设的N组超参数中所有的超参数均被选择，获取每组超参数对应的AIC值，并从所有AIC值中选出最小值，所述最小AIC值对应的一组超参数即为最优超参数，同时将所述最优超参数对应的核函数模型作为最终的模型。
7.根据权利要求6所述的一种基于稀疏核学习的单站GNSS瞬时速度和加速度构建方法，其特征在于，在所述步骤S5.7中，所述最终确定出的待求解参数，即为能使如下代价函数最小时对应的待求解参数，所述代价函数具体为：
Figure FDA0002204020980000062
其中：B＝C-cbT，κkj＝K(tk,tj)-K(tk-1,tj)，uj＝K(t1,tj)，u＝[u1 u2 … un]T
φ为稀疏正则化代价函数，y为位移增量序列观测向量，c为矩阵A中的第一列向量，C为矩阵A中除第一列向量之外的矩阵，A为κkj组成的矩阵，K为核函数，tk为第k个观测历元对应的时间，tk-1为第k-1个观测历元对应的时间，tj为第j个核函数的中心时刻，b为向量u中除第一个元素之外的向量部分，t1为初始时间，β为待求解参数，Qee为观测误差向量的协方差矩阵，μi为第i组正则化参数；
所述待求解参数的观测方程通过稀疏正则化方法进行求解，具体如下：
S5.7.1：将待求解参数进行初始化，具体为：
Figure FDA0002204020980000063
其中：B＝C-cbT，κkj＝K(tk,tj)-K(tk-1,tj)，uj＝K(t1,tj)，u＝[u1 u2 … un]T
β0为迭代前待求解参数的初步估计值，θ1和s1均为第一次迭代运算引入的辅助变量，c为矩阵A中的第一列向量，C为矩阵A中除第一列向量之外的矩阵，A为κkj组成的矩阵，K为核函数，tk为第k历元对应的时间，tk-1为第k-1历元对应的时间，tj为第j个核函数的中心时刻，b为向量u中除第一个元素之外的向量部分，t1为初始时间，Qee为观测误差向量的协方差矩阵，μi为第i组正则化参数，I为单位矩阵，y为位移增量序列观测向量；
S5.7.2：将所述迭代前待求解参数的初步估计值进行迭代计算，直至算法收敛，确定出每次迭代过程中所述待求解参数的大小，具体为：
Figure FDA0002204020980000071
其中：
Figure FDA0002204020980000072
B＝C-cbT，κkj＝K(tk,tj)-K(tk-1,tj)，uj＝K(t1,tj)，u＝[u1 u2… un]T
βm为第m次迭代时待求解参数的估计值，sm+1和θm+1为第m+1次迭代运算引入的辅助变量，
Figure FDA0002204020980000073
为θm-2λBTQee -1(Bθm-y)的软阈值函数，θm和sm为第m次迭代运算引入的辅助变量，λmax为最大特征值，c为矩阵A中的第一列向量，C为矩阵A中除第一列向量之外的矩阵，A为κkj组成的矩阵，K为核函数，tk为第k历元对应的时间，tk-1为第k-1历元对应的时间，tj为第j个核函数的中心时刻，b为向量u中除第一个元素之外的向量部分，t1为初始时间，Qee为观测误差向量的协方差矩阵，y为位移增量序列观测向量，βm-1为第m-1次迭代时待求解参数的估计值。
8.根据权利要求7所述的一种基于稀疏核学习的单站GNSS瞬时速度和加速度构建方法，其特征在于，在所述步骤S5.7.2中，所述软阈值函数的求解过程，具体如下：
S5.7.2.1：给定向量w，对所述给定的向量进行软阈值操作，得到一个同维数的向量，所述同维数向量的第j个元素，具体为：
Figure FDA0002204020980000074
其中：B＝C-cbT，κkj＝K(tk,tj)-K(tk-1,tj)，uj＝K(t1,tj)，u＝[u1 u2… un]T
w表示软阈值函数中的自变量，(|wj|-α)+为铰链函数，sgn(wj)为符号函数，θm为第m次迭代运算引入的辅助变量，λmax为最大特征值，c为矩阵A中的第一列向量，C为矩阵A中除第一列向量之外的矩阵，A为κkj组成的矩阵，K为核函数，tk为第k历元对应的时间，tk-1为第k-1历元对应的时间，tj为第j个核函数的中心时刻，b为向量u中除第一个元素之外的向量部分，t1为初始时间，Qee为观测误差向量的协方差矩阵，y为位移增量序列观测向量；
S5.7.2.2：对所述铰链函数进行求解，具体为：
Figure FDA0002204020980000081
其中：(|wj|-α)+为铰链函数，wj为w的第j个分量，w表示软阈值函数中的向量，α为所有核函数系数组成的向量；
S5.7.2.3：对所述符号函数进行求解，具体为：
Figure FDA0002204020980000082
其中：sgn(wj)为符号函数，wj为w的第j个分量，w表示软阈值函数中的向量。
9.根据权利要求8所述的一种基于稀疏核学习的单站GNSS瞬时速度和加速度构建方法，其特征在于，在所述步骤S5.7.2.1和步骤S5.8中，Qee -1(Bθm-y)和
Figure FDA0002204020980000083
均通过TMA算法进行求解，即依次执行如下公式，具体为：
Figure FDA0002204020980000084
Figure FDA0002204020980000085
Figure FDA0002204020980000091
其中：f＝Qee -1(Bθm-y)＝Qee -1g或
Figure FDA0002204020980000092
g＝Qeef
qk为Qhh,k的第k个对角线元素，Qhh,k为
Figure FDA0002204020980000093
的协方差矩阵，
Figure FDA0002204020980000094
为第k历元处的位移增量和接收机钟差增量组成的增量序列的估计值，rk为Rhh,k的第k个对角线元素，Rhh,k为
Figure FDA0002204020980000095
和
Figure FDA0002204020980000096
之间的互协方差矩阵，
Figure FDA0002204020980000097
为第k-1历元处的位移增量和接收机钟差增量组成的增量序列的估计值，n为历元的总个数，rk-1为Rhh,k的第k-1个对角线元素，fk为向量f的第k个元素，gk为向量g的第k个元素，
Figure FDA0002204020980000098
和
Figure FDA0002204020980000099
均为辅助变量。
10.根据权利要求6所述的一种基于稀疏核学习的单站GNSS瞬时速度和加速度构建方法，其特征在于，在所述步骤S6中，获取所述瞬时速度和加速度大小，具体如下：
S6.1：将所述确定出的待求解参数代入第一个核函数系数的求解公式中，获取第一个核函数系数的值；
S6.2：根据所述第一个核函数系数的值，通过所述最终的核函数模型，确定出相应方向的速度模型和加速度模型，具体为：
Figure FDA00022040209800000910
其中：
Figure FDA00022040209800000911
vt为方向分量中第t时刻对应的速度，at为方向分量中第t时刻对应的加速度，xt为方向分量中第t时刻对应的时间函数，t为时间，n为历元的总个数，αj为第j个核函数的系数，
Figure FDA00022040209800000912
为最优的核函数宽度参数，tj为第j个核函数的中心时刻；
S6.3：将任意时刻代入所述速度模型和加速度模型中，获取所述时刻对应的瞬时速度和加速度大小。
CN201910874925.1A 2019-09-17 2019-09-17 一种基于稀疏核学习的单站gnss瞬时速度和加速度构建方法 Active CN110703284B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910874925.1A CN110703284B (zh) 2019-09-17 2019-09-17 一种基于稀疏核学习的单站gnss瞬时速度和加速度构建方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910874925.1A CN110703284B (zh) 2019-09-17 2019-09-17 一种基于稀疏核学习的单站gnss瞬时速度和加速度构建方法

Publications (2)

Publication Number Publication Date
CN110703284A true CN110703284A (zh) 2020-01-17
CN110703284B CN110703284B (zh) 2022-11-29

Family

ID=69194426

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910874925.1A Active CN110703284B (zh) 2019-09-17 2019-09-17 一种基于稀疏核学习的单站gnss瞬时速度和加速度构建方法

Country Status (1)

Country Link
CN (1) CN110703284B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111669183A (zh) * 2020-06-30 2020-09-15 中南大学 一种压缩感知采样与重建方法、设备及存储介质
US20220221485A1 (en) * 2021-01-13 2022-07-14 PeiLiang Xu Speed and Acceleration Calculation and Measurement Method, Device, and Application Based on Regularization Algorithms

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106569241A (zh) * 2016-09-27 2017-04-19 北京航空航天大学 一种基于gnss的单频高精度定位方法
CN107917708A (zh) * 2017-11-10 2018-04-17 太原理工大学 基于贝叶斯压缩感知的gps/ins紧组合周跳探测和修复算法
CN109085629A (zh) * 2018-08-29 2018-12-25 广州市中海达测绘仪器有限公司 Gnss基线向量解算定位方法、装置和导航定位设备
CN110100190A (zh) * 2017-01-04 2019-08-06 高通股份有限公司 用于在视觉惯性测距中使用全球定位历元的滑动窗口的系统和方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106569241A (zh) * 2016-09-27 2017-04-19 北京航空航天大学 一种基于gnss的单频高精度定位方法
CN110100190A (zh) * 2017-01-04 2019-08-06 高通股份有限公司 用于在视觉惯性测距中使用全球定位历元的滑动窗口的系统和方法
CN107917708A (zh) * 2017-11-10 2018-04-17 太原理工大学 基于贝叶斯压缩感知的gps/ins紧组合周跳探测和修复算法
CN109085629A (zh) * 2018-08-29 2018-12-25 广州市中海达测绘仪器有限公司 Gnss基线向量解算定位方法、装置和导航定位设备

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
李昕等: "基于多普勒测速的GPS单历元动态定位算法", 《武汉大学学报·信息科学版》 *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111669183A (zh) * 2020-06-30 2020-09-15 中南大学 一种压缩感知采样与重建方法、设备及存储介质
CN111669183B (zh) * 2020-06-30 2022-04-19 中南大学 一种压缩感知采样与重建方法、设备及存储介质
US20220221485A1 (en) * 2021-01-13 2022-07-14 PeiLiang Xu Speed and Acceleration Calculation and Measurement Method, Device, and Application Based on Regularization Algorithms
WO2022151843A1 (zh) * 2021-01-13 2022-07-21 徐培亮 一种基于正则化算法的速度和加速度计算方法及测量装置
JP2022108737A (ja) * 2021-01-13 2022-07-26 培亮 徐 正則化アルゴリズムに基づく速度及び／又は加速度の計算方法、計算装置、測定装置及びその応用
JP7221562B2 (ja) 2021-01-13 2023-02-14 培亮 徐 正則化アルゴリズムに基づく速度及び／又は加速度の計算方法、計算装置、測定装置及びその応用

Also Published As

Publication number Publication date
CN110703284B (zh) 2022-11-29

Similar Documents

Publication Publication Date Title
CN111007557B (zh) 自适应运动学模型辅助的gnss载波相位与多普勒融合测速方法
CN107677272B (zh) 一种基于非线性信息滤波的auv协同导航方法
Lemoine et al. CNES/GRGS RL04 Earth gravity field models, from GRACE and SLR data
CN110703284B (zh) 一种基于稀疏核学习的单站gnss瞬时速度和加速度构建方法
CN111156987A (zh) 基于残差补偿多速率ckf的惯性/天文组合导航方法
CN110764127B (zh) 易于星载在轨实时处理的编队卫星相对定轨方法
CN107561565B (zh) 低轨航天器星载gnss差分相对导航换星的处理方法
CN110006427B (zh) 一种低动态高振动环境下的bds/ins紧组合导航方法
CN111024124B (zh) 一种多传感器信息融合的组合导航故障诊断方法
Yeomans et al. Radar astrometry of near-Earth asteroids
CN116774264B (zh) 基于低轨卫星机会信号多普勒的运动目标定位方法
CN110727002A (zh) 一种基于稀疏正则化的单频单站动态gnss载波相位信号周跳修复方法
CN115307643A (zh) 一种双应答器辅助的sins/usbl组合导航方法
Sarychev et al. Optimal regressors search subjected to vector autoregression of unevenly spaced TLE series
Biswas et al. Computationally efficient unscented Kalman filtering techniques for launch vehicle navigation using a space-borne GPS receiver
Kaniewski et al. Algorithms of position and velocity estimation in GPS receivers
CN115659196A (zh) 基于非线性偏差演化的天基光学观测短弧关联与聚类方法
US20230011501A1 (en) Particle filtering and navigation system using measurement correlation
US20100161284A1 (en) Nonlinear Variable Lag Smoother
CN113093092A (zh) 一种水下鲁棒自适应单信标定位方法
CN105953791B (zh) 单探测器x射线脉冲星导航分时观测方法及装置
Li et al. Exploring the Potential of Deep Learning Aided Kalman Filter for GNSS/INS Integration: A Study on 2D Simulation Datasets
CN117589178A (zh) 一种基于光学测量弧段的航天器机动检测方法
US20220358365A1 (en) Tightly coupled end-to-end multi-sensor fusion with integrated compensation
CN116148852A (zh) 基于空时连续的北斗InSAR三维高精度形变反演方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant