一种基于熵权法的BDS卫星轨道位置插值方法与流程

本发明属于卫星轨道插值的技术领域，具体涉及一种基于熵权法的bds卫星轨道位置插值方法。

背景技术：

北斗导航系统(beidounavigationsatellitesystem，bds)是中国自主研制的全球卫星定位系统，是继美国全球定位系统(globalpositioningsystem，gps)和俄罗斯的格洛纳斯系统(globalnavigationsatellitesystem，glonass)之后第三个成熟的卫星导航系统。北斗导航系统是全天候，全天时提供卫星导航信息的区域导航系统，此系统由定位于赤道上空的两颗地球同步卫星，地面中心站，用户终端三部分组成。该系统能提供四大功能：快速定位、实时导航、简短通信、精密授时。

在北斗卫星定位中，卫星位置的计算是一个关键环节。在北斗卫星数据后处理中，德国地学研究中心(geoforschungszentrum，gfz)提供的精密星历数据时间间隔过长，是不同时间段等间隔的卫星位置坐标，而在实际的应用中，需要10s，5s甚至更密集的时间点卫星位置坐标。因此根据gfz提供的精密星历数据选择正确的插值方法解算任意观测时刻的卫星位置具有重大的意义。本发明使用了拉格朗日插值法、切比雪夫插值法计算不同轨道卫星任一时刻的卫星位置，再利用熵权法完成两种插值方法的组合得出最终插值结果。

技术实现要素：

针对现有技术的不足，本发明提出一种基于熵权法的bds卫星轨道位置插值方法。

一种基于熵权法的bds卫星轨道位置插值方法，包括以下步骤：

步骤1：获取bds卫星轨道系统精密星历中的等间隔的已知历元时刻ti和ti，根据已知历元时刻ti和tj构造任意插值时刻的历元基函数li(t)，其中，t为任意插值时刻，i，j∈(0，n)，n为待插值阶数；

步骤2：根据bds卫星轨道系统精密星历中的已知等间隔历元位置(xi，yi，zi)和任意插值时刻的历元基函数li(t)，确定在t时刻的卫星轨道位置的拉格朗日插值[xl(t)，yl(t)，zl(t)]，其中，xl(t)＝(xl1，xl2，…xln)，yl(t)＝(yl1，yl2，…yln)，zl(t)＝(zl1，zl2，…zln)；xli为x方向拉格朗日插值的第i个待加权值，yli为y方向拉格朗日插值的第i个待加权值，zli为z方向拉格朗日插值的第i个待加权值；

步骤3：获取观测时段[t0，t0+δt]，对观测时段[t0，t0+δt]进行归一化处理，得到归一化后的观测时段τ∈[-1，1]，其中，t0为观测历元首时刻，δt为观测时长；

步骤4：构造卫星轨道位置插值的切比雪夫误差方程，根据归一化后的观测时段τ计算切比雪夫误差方程的系数矩阵a；

步骤5：根据卫星轨道位置插值的切比雪夫误差方程，以及切比雪夫误差方程的系数矩阵a，得到等间隔历元时刻多项式系数c，从而确定在t时刻的卫星轨道位置的切比雪夫插值[xq(t)，yq(t)，zq(t)]，其中，xq(t)＝(xq1，xq2，…xqn)，yq(t)＝(yq1，yq2，…yqn)，zq(t)＝(zq1，zq2，…zqn)，xqi为x方向切比雪夫插值的第i个待加权值，yqi为y方向切比雪夫插值的第i个待加权值，zqi为z方向切比雪夫插值的第i个待加权值；

步骤6：将卫星轨道位置的拉格朗日插值和对应时刻的卫星轨道位置的切比雪夫插值组成待加权项目数据矩阵

步骤7：计算待加权项目数据矩阵中各个方同的各待加权项的拉格朗日加权比重lli＝[lxli，lyli，lzli]和切比雪夫加权比重lqi＝[lxqi，lyqi，lzqi]，其中，lxli为x方向第i个待加权值的拉格朗日加权比重，lyli为y方向第i个待加权值的拉格朗日加权比重，lzli为z方向第i个待加权值的拉格朗日加权比重，lxqi为x方向第i个待加权值的切比雪夫加权比重，lyqi为y方向第i个待加权值的切比雪夫加权比重，lzqi为z方向第i个待加权值的切比雪夫加权比重；

步骤8：根据待加权项目数据矩阵中各个方向的各待加权项的拉格朗日加权比重lli＝[lxli，lyli，lzli]和切比雪夫加权比重lqi＝[lxqi，lyqi，lzqi]，计算各个方向上待加权值的熵值si＝[sxi，syi，szi]；

所述各个方向上待加权值的熵值sxi，syi，szi的计算公式如下所示：

sxi＝-k(lxlilnlxli+lxqilnlxqi).

syi＝-k(lylilnlyli+lyqilnlyqi).

szi＝-k(lzlilnlzli+lzqilnlzqi).

其中，

步骤9：根据各个方向上待加权值的熵值si＝[sxi，syi，szi]，计算各个方向上待加权值的熵权ei＝[exi，eyi，ezi]；

所述各个方向上待加权值的熵权exi，eyi，ezi的计算公式如下所示：

步骤10：根据各个方向上待加权值的熵权ei＝[exi，eyi，ezi]，以及在t时刻的卫星轨道位置的拉格朗日插值[xl(t)，yl(t)，zl(t)]和切比雪夫插值[xq(t)，yq(t)，zq(t)]，确定各个方向基于熵权法的轨道位置插值(wxi，wyi，wzi)，其中，wxi为x方向上的基于熵权法的轨道位置插值，wyi为y方向上的基于熵权法的轨道位置插值，wzi为z方向上的基于熵权法的轨道位置插值。

所述各个方向基于熵权法的轨道位置插值(wxi，wyi，wzi)的计算公式如下所示：

wxi＝eixqi+(1-ei)xli；

wyi＝exiyqi+(1-exi)yli；

wzi＝ezizqi+(1-ezi)zli。

所述构造的任意插值时刻的历元基函数li(t)如下所示：

所述t时刻的卫星轨道位置的拉格朗日插值[xl(t)，yl(t)，zl(t)]的计算公式如下所示：

所述构造的卫星轨道位置插值的切比雪夫误差方程v＝[vx，vy，vz]的计算公式如下所示：

vx＝acx-x；

vy＝acy-y；

vz＝acz-z；

其中，c＝[cx，cy，cz]为等间隔历元时刻多项式系数，

为切比雪夫误差方程的系数矩阵，|τ≤1|，n≥2。本发明的有益效果：

本发明提出一种基于熵权法的bds卫星轨道位置插值方法，本发明使用切比雪夫和拉格朗日插值法完成bds轨道插值，再利用熵权法求得改进插值算法的轨道插值，本发明提高了单一差值方法的插值精度，解决了使用单一插值方法可能导致的误差较大情况。

附图说明

图1为本发明具体实施方式中基于熵权法的bds卫星轨道位置插值方法的流程图；

图2为本发明具体实施方式中得到的不同插值情况误差比较图；

其中，(a)为实验一得到的不同历元的插值误差值；

(b)为实验二得到的不同历元的插值误差值；

(c)为实验三得到的不同历元的插值误差值；

(d)为实验四得到的不同历元的插值误差值；

(e)为实验五得到的不同历元的插值误差值；

(f)为实验六得到的不同历元的插值误差值。

具体实施方式

下面结合附图对本发明具体实施方式加以详细的说明。

一种基于熵权法的bds卫星轨道位置插值方法，如图1所示，包括以下步骤：

步骤1：获取bds卫星轨道系统精密星历中的等间隔的已知历元时刻ti和tj，根据已知历元时刻ti和tj构造任意插值时刻的历元基函数li(t)。

本实施方式中，构造的任意插值时刻的历元基函数li(t)如式(1)所示：

其中，t为任意插值时刻，i，j∈(0，n)，n为待插值阶数。

步骤2：根据bds卫星轨道系统精密星历中的已知等间隔历元位置(xi，yi，zi)和任意插值时刻的历元基函数li(t)，确定在t时刻的卫星轨道位置的拉格朗日插值[xl(t)，yl(t)，zl(t)]。

本实施方式中，t时刻的卫星轨道位置的拉格朗日插值[xl(t)，yl(t)，zl(t)]的计算公式如式(2)、(3)和(4)所示：

其中，xl(t)＝(xl1，xl2，…xln)，yl(t)＝(yl1，yl2，…yln)，zl(t)＝(zl1，zl2，…zln)；xli为x方向拉格朗日插值的第i个待加权值，yli为y方向拉格朗日插值的第i个待加权值，zli为z方向拉格朗日插值的第i个待加权值。

步骤3：获取观测时段[t0，t0+δt]，对观测时段[t0，t0+δt]进行归一化处理，得到归一化后的观测时段τ∈[-1，1]。

本实施方式中，归一化后的观测时段τ的计算公式如式(5)所示：

其中，t0为观测历元首时刻，δt为观测时长。

步骤4：构造卫星轨道位置插值的切比雪夫误差方程，根据归一化后的观测时段τ计算切比雪夫误差方程的系数矩阵a。

本实施方式中，构造的卫星轨道位置插值的切比雪夫误差方程v＝[vx，vy，vz]的计算公式如式(6)、(7)和(8)所示：

vx＝acx-x(6)

vy＝acy-y(7)

vz＝acz-z(8)

其中，c＝[cx，cy，cz]为等间隔历元时刻多项式系数，

根据公式(6)可以推出其展开公式(9)如下所示：

其中，

同理，可得公式(7)和公式(8)的展开公式。

切比雪夫误差方程的系数矩阵a如式(10)所示：

其中，|τ≤1|，n≥2。。

步骤5：根据卫星轨道位置插值的切比雪天误差方程，以及切比雪天误差方程的系数矩阵a，得到等间隔历元时刻多项式系数c，从而确定在t时刻的卫星轨道位置的切比雪夫插值[xq(t)，yq(t)，zq(t)]。

本实施方式中，将公式(6)采用最小二乘概念整理为公式(11)所示：

vx^tpvx＝min(11)

其中，p＝i为权值。

可以得到等间隔历元时刻多项式系数cx如式(12)所示：

cx＝(a^ta)^-1(a^tx)(12)

同理，可以得到公式等间隔历元时刻多项式系数cy和cz。

本实施方式中，根据公式等间隔历元时刻多项式系数cx、cy和cz确定在t时刻的卫星轨道位置的切比雪夫插值xq(t)，yq(t)和zq(t)如式(13)、(14)和(15)所示：

xq(t)＝(xq1，xq2，…xqn)(13)

xq(t)＝(yq1，yq2，…yqn)(14)

zq(t)＝(zq1，zq2，…zqn)(15)

其中，x方向切比雪夫插值的第i个待加权值xqi如式(16)所示：

y方向切比雪夫插值的第i个待加权值yqi如式(17)所示：

z方向切比雪夫插值的第i个待加权值zqi如式(18)所示：

步骤6：将卫星轨道位置的拉格朗日插值和对应时刻的卫星轨道位置的切比雪夫插值组成待加权项目数据矩阵

如式(19)、(20)和(21)所示：

步骤7：计算待加权项目数据矩阵中各个方向的各待加权项的拉格朗日加权比重lli＝[lxli，lyli，lzli]和切比雪夫加权比重lqi＝[lxqi，lyqi，lzqi]。

本实施方式中，x方向第i个待加权值的拉格朗日加权比重lxli的计算公式如式(22)所示：

同理，可得y方向第i个待加权值的拉格朗日加权比重lyli，z方向第i个待加权值的拉格朗日加权比重lzli。

x方向第i个待加权值的切比雪夫加权比重lxqi的计算公式如式(23)所示：

同理，可得y方向第i个待加权值的切比雪夫加权比重lyqi，z方向第i个待加权值的切比雪夫加权比重lzqi。

步骤8：根据待加权项目数据矩阵中各个方向的各待加权项的拉格朗日加权比重lli＝[lxli，lyli，lzli]和切比雪夫加权比重lqi＝[lxqi，lyqi，lzqi]，计算各个方向上待加权值的熵值si＝[sxi，syi，szi]。

本实施方式中，各个方向上待加权值的熵值sxi，syi，szi的计算公式如式(24)、(25)和(26)所示：

sxi＝-k(lxlilnlxli+lxqilnlxqi)(24)

syi＝-k(lylilnlyli+lyqilnlyqi)(25)

szi＝-k(lzlilnlzli+lzqilnlzqi)(26)

其中，

步骤9：根据各个方向上待加权值的熵值si＝[sxi，syi，szi]，计算各个方向上待加权值的熵权ei＝[exi，eyi，ezi]。

本实施方式中，各个方向上待加权值的熵权exi，eyi，ezi的计算公式如下如式(27)、(28)和(29)所示：

步骤10：根据各个方向上待加权值的熵权ei＝[exi，eyi，ezi]，以及在t时刻的卫星轨道位置的拉格朗日插值[xl(t)，yl(t)，zl(t)]和切比雪夫插值[xq(t)，yq(t)，zq(t)]，确定各个方向基于熵权法的轨道位置插值(wxi，wyi，wzi)，其中，wxi为x方向上的基于熵权法的轨道位置插值，wyi为y方向上的基于熵权法的轨道位置插值，wzi为z方向上的基于熵权法的轨道位置插值。

本实施方式中，各个方向基于熵权法的轨道位置插值(wxi，wyi，wzi)的计算公式如式(30)、(31)和(32)所示：

wxi＝eixqi+(1-ei)xli(30)

wyi＝exiyqi+(1-exi)yli(31)

wzi＝ezizqi+(1-ezi)zli(32)

本发明实施方式中，采用了六组实验得到的不同历元的插值误差值，如图2所示。不同实验的误差值平均值如表1所示：

表1不同实验的误差值平均值