CN108919313A - 利用最优数值导数的gnss多普勒观测值生成方法 - Google Patents

Abstract

本发明公开了一种利用最优数值导数的GNSS多普勒观测值生成方法。针对常规方法仅利用必要等间隔历元的GNSS载波观测值以及不具有稳健性、最优性等不足，本发明提出基于载波相位/伪距观测值的多历元连续性构建多项式基函数和数值导数方程，利用最小二乘估计确定最优数值导数算子并生成多普勒观测值。本发明相较于常规方法有以下优点：(1)同时适用于载波相位和伪距观测值；(2)利用多历元GNSS连续性，满足观测值存在非均匀采样时的多普勒生成任务；(3)顾及了观测值随卫星高度角的精度变化或方差异常，能够提高多普勒观测生成精度；(4)所用历元数根据GNSS采样率可灵活设置，通过最小二乘估计保证了数值导数算子最优性。

Description

利用最优数值导数的GNSS多普勒观测值生成方法

技术领域

本发明涉及一种利用最优数值导数的GNSS多普勒观测值生成方法，属于卫星导航定位技术和测绘科学技术领域。

背景技术

随着全球卫星导航系统(GNSS)的快速发展，利用GNSS进行位置、速度、姿态解算和授时等应用越来越广泛。其中，GNSS多普勒观测值被广泛用于载体位置速度测量、周跳探测和伪距平滑等研究与应用。然而，GNSS多普勒观测值并不是GNSS接收机原始观测值输出格式的标准要求，用户目前需要基于原始载波相位观测值进一步利用常规生成算子获得多普勒观测值。例如，专利号为CN106772472A的中国专利“一种顾及GNSS接收机钟跳的导出多普勒观测值构造方法”直接采用连续两个历元的载波相位观测值进行差分的方法导出多普勒观测值；专利号为CN102253399A的中国专利“一种利用载波相位中心值的多普勒差分补偿测速方法”利用载波相位一阶中心微分值对原始多普勒观测值进行补偿从而得到修正后的测速值。

现有导出多普勒观测值是利用连续两个或必要历元的载波相位观测值通过常规差分算子求得，该方法存在如下问题：一、常规差分算子仅适用于载波相位观测值，无法避免周跳难题；二、仅采用两个历元或等间隔历元的观测值难以利用多历元GNSS轨迹连续性，无法完成非均匀采样条件下多普勒生成任务；三、所用历元数是常规差分算子的必要历元数，不能根据采样率和移动时间窗口长度进行最优算子估计；四、不考虑观测值精度随历元的变化情况，多普勒观测值受非等精度观测或异常方差影响严重。

发明内容

本发明目的在于提供一种利用最优数值导数的GNSS多普勒观测值生成方法，该方法能同时适用于载波相位和伪距观测值、均匀和非均匀采样时间，且具有最优性和稳健性。

为实现上述目的，本发明提供一种利用最优数值导数的GNSS多普勒观测值生成方法，包括以下步骤：

第一步、定义第i个移动窗口的时间区间为t∈[t_i,t_i+K]，设置数值导数算子长度为K+1，数值导数算子记为d_i＝[d_i,0 … d_i,k … d_i,K]^T，k∈[0,K]。将第i个移动窗口的时间变量t_i+k∈[t_i,t_i+K]处理为标准化时间变量τ_i+k∈[-1,1]，则标准化时间变量及其导数分别为和

第二步、选取一般多项式或递推多项式作为时间序列基函数，并设置多项式阶数N，且满足N＜K；其中，第j阶时间序列基函数记为F_j(t_i+k)、基函数导数记为j∈[0,N]。基函数导数可采用如下任意一种计算方法。

①时序基函数采用N阶一般多项式，导数递推形式为：

②时序基函数采用N阶切比雪夫多项式T_N(τ_i+k)，则导数递推形式为：

式中，

③时序基函数采用N阶勒让德多项式P_N(τ_i+k)，则导数递推形式为：

式中，

④时序基函数采用N阶盖根鲍尔多项式G_N(τ_i+k)，则导数递推形式为：

式中，λ为任一非零实数。

第三步、动态构建数值导数算子估计方程；记导出多普勒的历元为t′_I，t_i≤t′_I≤t_i+K，则第i个移动窗口的多普勒数值导数算子估计方程如下：

第四步、动态构建GNSS观测值方差-协方差矩阵；

读取GNSS广播星历并计算卫星高度角，进而采用高度角模型确定GNSS观测值方差阵Q_i：

Q_i＝diag[q(t_i),q(t_i+1),…,q(t_i+K)]

式中，E(t_i)表示历元t_i的卫星高度角，对于GNSS伪距或载波相位，q₁,q₂取不同的经验值。

第五步、利用最小二乘估计方法求解最优数值导数算子；

基于最小二乘准则求解第i个移动窗口的数值导数算子，最优数值导数算子的估计公式如下：

第六步、获取多普勒观测值时间序列；

读取移动时间窗口内的GNSS观测值y(t_i+k)，基于第i个移动窗口的数值导数算子计算移动窗口的多普勒观测值y⁽¹⁾(t′_I)，则最小二乘最优数值求导所得的多普勒观测值为：

式中，d_i,k为第i个数值导数算子的第k个值。若y(t_i+k)表示第i个移动窗口载波相位、伪距观测值，所得即为相位导出、伪距导出多普勒观测时序；

与现有多普勒生成方法相比，本发明不仅适用于GNSS载波相位导出多普勒观测值，也适用于GNSS伪距导出多普勒观测值。而且，充分利用了多历元GNSS观测值连续性、GNSS观测值随历元或卫星高度角的精度变化或方差异常影响，满足非均匀采样的多普勒观测值生成任务，并通过最小二乘估计保证了数值导数算子的最优性。

附图说明

图1为本发明方法流程图；

图2为实施例1中心数值导数算子的伪距导出多普勒观测值；

图3为实施例1伪距导出多普勒与原始多普勒观测值的偏差；

图4为实施例2后向数值导数算子的载波相位导出多普勒观测值；

图5为实施例2载波相位导出多普勒与原始多普勒观测值的偏差。

具体实施方式

下面结合附图对本发明进一步说明：

如图1所示，一种利用最优数值导数的GNSS多普勒观测值生成方法，包括以下步骤：

第一步、定义第i个移动窗口的时间区间为t∈[t_i,t_i+K]，设置数值导数算子长度为K+1，数值导数算子记为d_i＝[d_i,0 … d_i,k … d_i,K]^T，k∈[0,K]。将第i个移动窗口的时间变量t_i+k∈[t_i,t_i+K]处理为标准化时间变量τ_i+k∈[-1,1]，则标准化时间变量及其导数分别为和

第二步、选取一般多项式或递推多项式作为时间序列基函数，并设置多项式阶数N，且满足N＜K；其中，第j阶时间序列基函数记为F_j(t_i+k)、基函数导数记为j∈[0,N]。基函数导数可采用如下任意一种计算方法。

①时序基函数采用N阶一般多项式，导数递推形式为：

②时序基函数采用N阶切比雪夫多项式T_N(τ_i+k)，则导数递推形式为：

式中，

③时序基函数采用N阶勒让德多项式P_N(τ_i+k)，则导数递推形式为：

式中，

④时序基函数采用N阶盖根鲍尔多项式G_N(τ_i+k)，则导数递推形式为：

式中，λ为任一非零实数。

第三步、动态构建数值导数算子估计方程；记导出多普勒的历元为t′_I，t_i≤t′_I≤t_i+K，则第i个移动窗口的多普勒数值导数算子估计方程如下：

第四步、动态构建GNSS观测值方差-协方差矩阵；

读取GNSS广播星历并计算卫星高度角，进而采用高度角模型确定GNSS观测值方差阵Q_i：

Q_i＝diag[q(t_i),q(t_i+1),…,q(t_i+K)]

式中，E(t_i)表示历元t_i的卫星高度角，对于GNSS伪距或载波相位，q₁,q₂取不同的经验值。

第五步、利用最小二乘估计方法求解最优数值导数算子；

基于最小二乘准则求解第i个移动窗口的数值导数算子，最优数值导数算子的估计公式如下：

第六步、获取多普勒观测值时间序列；

读取移动时间窗口内的GNSS观测值y(t_i+k)，基于第i个移动窗口的数值导数算子计算移动窗口的多普勒观测值y⁽¹⁾(t′_I)，则最小二乘最优数值求导所得的多普勒观测值为：

式中，d_i,k为第i个数值导数算子的第k个值。若y(t_i+k)表示第i个移动窗口载波相位、伪距观测值，所得即为相位导出、伪距导出多普勒观测时序；

所述的构建数值导数算子估计方程有K+1种方式，两种典型方式包括中心数值导数算子和后向数值导数算子。其中，中心数值导数算子的理论精度最高，后向数值导数算子的历史数据利用率最高。下面是选取2017年6月1日(年积日152，采样率5s)10时至14时的香港CORS网站(http://www.geodetic.gov.hk)HKKS的观测值文件和广播星历文件(GPS卫星PRN2)，分别采用后向数值导数算子和中心数值导数算子的实施例作为具体实施方式对本发明进行详细说明。

实施例1：中心数值导数算子的伪距导出多普勒观测值实例

基于伪距观测值C1利用中心数值导数算子获取伪距导出多普勒观测值，并以原始多普勒观测值D1作为参考进行比较，以偏差和标准中误差评定本发明效果。此例，所得伪距导出多普勒观测值时间序列分别见图2所示，与原始多普勒观测值的偏差见图3所示。从图3可以看出，中心数值导数算子所生成的伪距导出多普勒观测值与原始多普勒观测值之间的偏差不超过0.3周/秒，标准中误差约为0.03周/秒，从而验证了本发明方法的正确性和有效性。中心数值导数算子导出多普勒观测值精度高于其它数值导数算子，采用中心数值导数算子的伪距导出多普勒观测值时间序列生成的具体步骤如下：

第一步、定义第i个移动窗口的时间区间为t∈[t_i,t_i+K]，设置数值导数算子长度为K+1＝31。第i个移动窗口的标准化时间变量及其导数分别为和

第二步、选取切比雪夫递推多项式作为时间序列基函数，并设置多项式阶数N＝3；其中，基函数导数(j∈[0,N])递推计算公式如下：

式中，

第三步、动态构建中心数值导数算子估计方程；记导出多普勒的历元为t′_I＝(t_i+t_i+K)/2，则第i个移动窗口的多普勒观测值中心数值导数算子估计方程如下：

第四步、动态构建GNSS观测值方差-协方差矩阵。

读取GNSS广播星历并计算卫星高度角，进而采用高度角模型确定GNSS观测值方差阵Q_i：

Q_i＝diag[q(t_i),q(t_i+1),…,q(t_i+K)]

式中，E(t_i)为历元t_i卫星高度角，GNSS伪距取q1＝0.3,q2＝0.5m。

第五步、利用最小二乘估计方法求解最优中心数值导数算子。

基于最小二乘准则求解第i个移动窗口的中心数值导数算子，最优中心数值导数算子的估计公式如下：

第六步、获取伪距导出多普勒观测值时间序列。

读取移动时间窗口内的GNSS伪距观测值y(t_i+k)，基于第i个移动窗口的最优中心数值导数算子计算移动窗口的多普勒观测值则最优中心数值导数算子所得的伪距导出多普勒观测值为：

式中，为第i个最优中心数值导数算子的第k个值。

实施例2：后向数值导数算子的相位导出多普勒观测值实例

基于载波相位观测值L1利用后向数值导数算子获取相位导出多普勒观测值，并以原始多普勒观测值D1作为参考进行比较，以标准偏差和中误差评定本发明效果。此例，所得载波相位导出多普勒观测值时间序列分别见图4所示，与原始多普勒观测值的偏差见图5所示。从图5可以看出，后向数值导数算子所生成的相位导出多普勒观测值与原始多普勒观测值之间的偏差不超过0.3周/秒，标准中误差约为0.03周/秒，从而验证了本发明方法的正确性和有效性。后向数值导数算子导出多普勒观测值的历史数据利用率高于其它数值导数算子，采用后向数值导数算子的相位导出多普勒观测值时间序列生成的具体步骤如下：

第一步、定义第i个移动窗口的时间区间为t∈[t_i,t_i+K]，设置数值导数算子长度为K+1＝31。第i个移动窗口的标准化时间变量及其导数分别为和

第二步、选取切比雪夫递推多项式作为时间序列基函数，并设置多项式阶数N＝3；其中，基函数导数(j∈[0,N])递推计算公式如下：

式中，

第三步、动态构建后向数值导数算子估计方程；记导出多普勒的历元为t′_I＝t_i+K，则第i个移动窗口的多普勒观测值后向数值导数算子估计方程如下：

第四步、动态构建GNSS观测值方差-协方差矩阵。

读取GNSS广播星历并计算卫星高度角，进而采用高度角模型确定GNSS观测值方差阵Q_i：

Q_i＝diag[q(t_i),q(t_i+1),…,q(t_i+K)]

式中，E(t_i)为历元t_i的卫星高度角，载波相位取q₁＝3,q₂＝5mm。

第五步、利用最小二乘估计方法求解最优后向数值导数算子。

基于最小二乘准则求解第i个移动窗口的后向数值导数算子，最优后向数值导数算子的估计公式如下：

第六步、获取载波相位导出多普勒观测值时间序列。

读取移动时间窗口内的GNSS载波相位观测值y(t_i+k)，基于第i个移动窗口的最优后向数值导数算子计算移动窗口的多普勒观测值则最优后向数值导数算子所得的载波相位导出多普勒观测值为：

式中，为第i个最优后向数值导数算子的第k个值。

Claims (3)

1.利用最优数值导数的GNSS多普勒观测值生成方法，其特征在于，包括以下步骤：

第一步、定义第i个移动窗口的时间区间为t∈[t_i,t_i+K]，设置数值导数算子长度为K+1，数值导数算子记为d_i＝[d_i,0 … d_i,k … d_i,K]^T，k∈[0,K]；

第二步、选取一般多项式或递推多项式作为时间序列基函数，设置多项式阶数N，且满足N＜K；同时，第j阶时序基函数记为F_j(t_i+k)、基函数导数记为j∈[0,N]；

第三步、动态构建数值导数算子估计方程；记导出多普勒的历元为t′_I，t_i≤t′_I≤t_i+K，则第i个移动窗口的多普勒数值导数算子估计方程如下：

第四步、动态构建GNSS观测值方差-协方差矩阵；

读取GNSS广播星历并计算卫星高度角，进而采用高度角模型确定GNSS观测值方差阵Q_i：

Q_i＝diag[q(t_i),q(t_i+1),…,q(t_i+K)]

式中，E(t_i)表示历元t_i的卫星高度角，对于GNSS伪距或载波相位，q₁,q₂取不同的经验值；

第五步、利用最小二乘估计方法求解最优数值导数算子；

基于最小二乘准则求解第i个移动窗口的数值导数算子，最优数值导数算子的估计公式如下：

第六步、获取多普勒观测值时间序列；

读取移动时间窗口内的GNSS观测值y(t_i+k)，基于第i个移动窗口的数值导数算子计算移动窗口的多普勒观测值y⁽¹⁾(t′_I)，则最小二乘最优数值求导所得的多普勒观测值为：

式中，d_i,k为第i个数值导数算子的第k个值。若y(t_i+k)表示第i个移动窗口载波相位、伪距观测值，所得即为相位导出、伪距导出多普勒观测时序。

2.根据权利要求1所述的GNSS多普勒观测值生成方法，其特征在于，第三步、第五步和第六步所述的最优数值导数算子采用中心数值导数算子，其中t′_I＝(t_i+t_i+K)/2，所述中心数值导数算子估计公式为：

相应多普勒观测值生成方法为：

3.根据权利要求1所述的GNSS多普勒观测值生成方法，其特征在于，第三步、第五步和第六步所述的最优数值导数算子采用后向数值导数算子，其中t′_I＝t_i+K，所述后向数值导数算子估计公式为：

相应多普勒观测值生成方法为：