北斗IGSO/MEO卫星伪距码偏差一步建模法的制作方法

本发明属于全球导航系统领域的定位定轨模糊度固定技术，特别涉及一种北斗igso/meo卫星伪距码偏差一步建模法。

背景技术：

模糊度收敛时间是目前卫星导航定位定轨技术的研究热点和难点问题。北斗系统igso/meo卫星的伪距码观测值存在系统性偏差，见图1，影响载波相位模糊度解算的收敛速度和可靠性。andre′hauschild、olivermontenbruck等学者提出北斗igso/meo卫星的码-相组合观测值存在与高度角相关的系统性误差^[1-2]。

此外，lambertwanninger学者通过分析消除了电离层延迟、对流层延迟、距离相关的误差的mp(multi-path)组合观测值，采用第一步分弧段建模和第二步对相同类型卫星所有弧段增加mp组合观测值之和为零的约束，分别建立了igso/meo卫星与高度角及频率相关的两套伪距码偏差离散型改正模型^[3]。louyidong等学者按照相同的建模方法提出三阶多项式的连续型改正模型^[4]。由于每个弧段建模时对未知的模糊度参数采取的是时序平均值作为真值消去的策略，因此模糊度消除误差会被吸收到建模结果中。其次，不同测站不同弧段能追踪到的可视卫星不同，对同类型卫星所有参与建模的数据整体加上一个所有观测值之和为零的约束，也会影响北斗卫星的伪距码偏差建模精度，最后，每颗卫星的唯一性也会导致相同类型不同卫星的伪距码偏差存在差异。

文中涉及如下参考文献：

[1]andre′hauschild,olivermontenbruck,jean-mariesleewaegen,lennardhuisman,peterj.g.teunissen,characterizationofcompassm-1signals.gpssolut(2012)16:117–126,doi10.1007/s10291-011-0210-3.

[2]olivermontenbruck,andre′hauschild,petersteigenberger,urshugentobler,peterteunissen,shinichinakamura,initialassessmentofthecompass/beidou-2regionalnavigationsatellitesystem.gpssolut(2013)17:211–222,doi10.1007/s10291-012-0272-x.

[3]lambertwanninger,susannebeer,beidousatellite-inducedcodepseudorangevariations:diagnosisandtherapy.gpssolut(2015)19:639–648,doi10.1007/s10291-014-0423-3.

[4]louy,gongx,gus,etal.assessmentofcodebiasvariationsofbdstriple-frequencysignalsandtheirimpactsonambiguityresolutionforlongbaselines[j].gpssolutions,2016:1-10.

技术实现要素：

针对北斗系统长距离双差模糊度和精密单点定位非差模糊度难以快速有效固定的问题，本发明提出了一种可快速有效固定模糊度的北斗igso/meo卫星伪距码偏差一步建模法。

本发明思路为：

有效顾及了模糊度消除策略的正确性，及数据处理过程中各历元依据高度角定权、各弧段依据历元数定权、不同测站依据先验值定权的灵活性，并且充分利用法方程综合解算模型系数，便于各弧段有效衔接，合理规避了所有卫星观测值之和为零的约束条件下存在的忽略不同弧段可视卫星存在差异的客观问题。

本发明技术方案如下：

一种北斗igso/meo卫星伪距码偏差一步建模法，包括步骤：

步骤1，对指定单颗北斗igso/meo卫星指定测站，根据指定测站的观测数据进行法方程叠加，本步骤进一步包括：

1.1利用igs精密卫星轨道、精密卫星钟差产品和观测数据，采用精密单点定位法解算指定测站坐标；

1.2对观测数据逐历元进行周跳探测，将第一历元或相位观测值发生周跳的当前历元记为参考历元t0，计算t0对应的卫星高度角elv(t0)及各频率mp组合观测值mpi(t0)，将mpi(t0)表示为elv(t0)的三阶多项式，i表示频率；然后，将t0的下一历元记为当前历元t1，执行子步骤1.3；

1.3计算当前历元t1对应的卫星高度角elv(t1)及各频率mp组合观测值mpi(t1)，将mpi(t1)表示为elv(t1)的三阶多项式；

1.4将mpi(t0)和mpi(t1)的三阶多项式按权重做差得：

其中，a1、a2、a3为三阶多项式系数；和分别表示t1和t0观测数据的权重；

历元观测数据的权重根据历元对应的卫星高度角确定，若历元t对应的卫星高度角elv(t)不大于30°，该历元观测数据的权重为4sin²(elv(t))；否则，权重为1；

1.5构建矩阵x＝[a1a2a3]^t、x为待估系数矩阵，计算当前历元t1频率i的法方程要素和将和按频率分别对应累加至前一历元更新的法方程要素，得当前历元t1更新的法方程要素；

1.6对当前参考历元t0后的观测数据逐历元重复执行子步骤1.3～1.5，直至探测到相位观测值发生周跳的历元，将以当前参考历元t0为起点、该发生周跳历元的前一历元为终点间的历元弧段记为当前历元弧段，该发生周跳历元的前一历元更新的法方程要素即当前历元弧段的法方程要素，保存当前历元弧段法方程要素和历元数；然后，以该发生周跳的历元重置参考历元t0，采用子步骤1.3～1.5方法获得下一历元弧段法方程要素；

1.7将各历元弧段的法方程要素叠加得指定测站对应的法方程要素和ni和li表示指定测站频率i对应的法方程要素，ni,n和li,n表示频率i对应的第n个历元弧段的法方程要素，k表示历元弧段数，pi,n为频率i对应的第n个历元弧段的权值，表示频率i对应的第n个历元弧段的历元数；

步骤2，对指定单颗北斗igso/meo卫星变换指定测站，采用步骤1方法获得变换后指定测站各频率对应的法方程要素；

步骤3，根据需求以及指定测站所处环境或指定测站历史观测数据质量，自由设定不同指定测站观测数据的权重，对指定单颗卫星各指定测站对应的法方程要素按频率分别进行加权叠加，获得各指定单颗卫星的法方程要素；

步骤4，在各频率下，分别根据指定单颗北斗igso/meo卫星的法方程要素构建法方程，整体求解法方程得单颗北斗igso/meo卫星各频率对应的待估系数矩阵x＝[a1a2a3]^t，即三阶多项式改正模型系数。

上述周跳探测采用turboedit法进行。

步骤3中，各指定测站观测数据的权重采用等权策略设定。

上述方法还包括对步骤4所获得的三阶多项式改正模型系数进行循环迭代剔除粗差，具体为：

在各频率下分别进行：

5.1根据各指定北斗igso/meo卫星的三阶多项式改正模型系数，计算三阶多项式改正模型系数的均方差误差rms，根据rms设置阈值；

5.2逐指定北斗igso/meo卫星逐历元，利用步骤4估计的三阶多项式改正模型系数以及当前历元和当前参考历元的卫星高度角，根据mpi组合观测值的三阶多项式估算当前历元和当前参考历元的组合观测值之差，记为

5.3将利用原始观测值计算的当前历元和当前参考历元的实际mpi组合观测值做差，所得差值记为δmpi；

5.4将与δmpi做差，若差值绝对值小于阈值，则当前历元观测数据为有效数据；否则，剔除当前历元观测数据；

5.5根据保留的历元观测数据，重新执行步骤1～3，求解待估系数矩阵x＝[a1a2a3]^t。

上述阈值为3～10倍的rms。

与现有技术相比，本发明具有以下优点和有益效果：

(1)本发明利用相同弧段内各历元与参考历元做差，准确消除未知的整周模糊度参数，提高模型精度。

(2)本发明将所有建模数据加入特定的法方程进行整体解算，求取三阶多项式改正模型系数，对不同弧段依据历元数定权，不同测站利用历史观测数据质量的经验值定权，进一步提高建模精度。

(3)本发明依据一步建模法建立了各颗igso/meo卫星的高度角及频率相关的三阶多项式改正模型，提高各颗卫星码-相观测值的可靠性，有益于减少模糊度收敛时间。

附图说明

图1为igso/meo卫星伪距系统性偏差示意图；

图2为不同igso卫星及所有igso卫星统一的伪距码偏差随高度角变化的曲线图，其中，图(a)、(b)、(c)分别对应频率b1、b2、b3；

图3为不同meo卫星及所有meo卫星统一的伪距码偏差随高度角变化的曲线图，其中，图(a)、(b)、(c)分别对应频率b1、b2、b3；

图4为原始伪距码观测值经三阶多项式改正模型改正前后的mp组合观测值及高度角的时间序列，其中，图(a)为改正后的mp组合观测值及高度角的时间序列，图(b)为改正前的mp组合观测值及高度角的时间序列，历元采样间隔为30s；

图5为本发明的具体流程图。

具体实施方式

下面将结合具体实施方式进一步说明本发明，见图5，本发明包括建模和粗差剔除两部分。

建模部分的具体步骤如下：

步骤1，对指定单颗北斗igso/meo卫星指定测站，根据指定测站的观测数据进行法方程叠加。

本步骤进一步包括：

步骤1.1，利用igs精密卫星轨道、精密卫星钟差产品和指定测站的观测数据，采用精密单点定位(precisepointpositioning，简称ppp)法解算指定测站坐标。

步骤1.2，对指定测站的观测数据进行周跳探测，若当前历元为第一历元或当前历元的相位观测值发生了周跳，将当前历元设置为参考历元t0，计算参考历元t0对应的卫星高度角elv(t0)及各频率mp组合观测值mpi(t0)，mpi(t0)表示参考历元t0频率i的mp组合观测值；若当前历元不为第一历元且当前历元的相位观测值未发生周跳，计算当前历元对应的卫星高度角elv(t1)及各频率mp组合观测值mpi(t1)，t1表示当前历元，i表示频率。本实施例中，采用已有的turboedit法进行周跳探测。

记i(i＝1,2,3)分别表示北斗系统的三个频率b1、b2、b3，b1、b2、b3分别表示1561.098mhz、1207.14mhz、1268.52mhz。j、k为与i相应的另外两个频率，即当i＝1时，j＝2，k＝3；当i＝2时，j＝1，k＝3；当i＝3时，j＝1，k＝2。λi、λj、λk分别为频率i、j、k对应的波长，分别为频率i、j、k对应的相位观测值；pi表示频率i对应的伪距码观测值，bi,r、mmp、δmp分别表示频率i对应的卫星端码偏差、接收机端码偏差、mp组合观测值的多路径、mp组合观测值的观测噪声。

当前历元t1频率i的mp组合观测值mpi(t1)计算公式如下：

式(1)～(2)中：

pi表示当前历元t1频率i对应的伪距码观测值；

mijk为比例系数，由式(2)计算获得。

实际应用中为计算方便，通常令j＝i或k＝i，即mp组合观测值中一个相位观测值与伪距码观测值频率相同。

步骤1.3，将子步骤1.2所得mp组合观测值表示为卫星高度角的三阶多项式，见式(3)～(4)，三阶多项式系数记为am，am即三阶多项式中m次项的系数，m＝0,1,2,3。

参考历元t0的mp组合观测值mpi(t0)所表示的三阶多项式如下：

mpi(t0)＝a0+a1·elv(t0)+a2·elv(t0)²+a3·elv(t0)³(3)

当前历元t1的mp组合观测值mpi(t1)所表示的三阶多项式如下：

mpi(t1)＝a0+a1·elv(t1)+a2·elv(t1)²+a3·elv(t1)³(4)

步骤1.4，根据各历元的卫星高度角确定各历元观测数据的权重，当历元t的卫星高度角elv(t)≤30°时，历元t观测数据的权重为4sin²(elv(t))；否则，历元t观测数据的权重为1；t表示任意历元。

记参考历元t0的mp组合观测值mpi(t0)依据卫星高度角确定的权重为记当前历元t1的mp组合观测值mpi(t1)依据卫星高度角确定的权重为则将mpi(t0)和mpi(t1)的三阶多项式按照卫星高度角定权后做差，得式(5)，如下：

步骤1.5，以矩阵形式表示公式(5)并构建法方程。

构建矩阵构建系数矩阵x＝[a1a2a3]^t，构建矩阵则式(5)等效于：

因待估系数矩阵x在不同历元是公共的，因此本发明采取可进行方程整体求解的法方程叠加原理进行系数矩阵估计，将式(6)两侧均乘以矩阵的转置矩阵，有：

记和即法方程要素，i表示北斗系统的三个不同频率。

以下一历元为当前历元，对下一历元重复执行子步骤1.2～1.5，并将所得历元的法方程要素累加到前一历元更新的法方程要素，直至探测到相位观测值发生周跳的当前历元，此时，以相位观测值发生周跳的当前历元重置参考历元t0，然后，基于重置的参考历元重新执行子步骤1.2～15。

重置参考历元时，保存当前历元弧段最后更新的法方程要素和历元数，然后对下一历元弧段按照上述方法计算法方程要素和历元数。

步骤1.6，重置参考历元时，则表示当前历元弧段结束，保存该当前历元弧段各历元法方程要素叠加的结果，即矩阵ni,n和li,n，将ni,n和li,n记为当前历元弧段的法方程要素；同时，保存当前历元弧段历元数numi,n，n＝1,2,3......k，n代表历元弧段编号，k表示历元弧段数量。

对当前历元弧段所有历元按照步骤1.2～1.6处理至最后历元，逐历元按照相位观测值的周跳情况更新法方程叠加要素。

步骤1.7采用子步骤1.2～1.5计算各历元弧段的法方程要素及历元数，求取各历元弧段权值pi，n为频率i对应的第n个历元弧段权值。

步骤1.8，叠加各历元弧段的法方程要素，见公式(8)～(9)：

其中，ni和li表示指定测站频率i对应的法方程要素，ni,n和li,n表示频率i对应的第n个历元弧段法方程要素，k表示历元弧段数，pi,n为频率i对应的第n个历元弧段权值，numi,n表示频率i对应的第n个历元弧段历元数。

步骤2，对指定单颗北斗igso/meo卫星变换指定测站，采用步骤1方法获得变换后指定测站对应的法方程要素。

步骤3，根据需求以及指定测站所处环境或指定测站历史观测数据质量自由设定不同指定测站观测数据的权重。本发明测站均为igs测站，因此选取等权策略。

根据各指定测站观测数据的权重，对指定单颗卫星的各指定测站对应的法方程要素进行加权叠加，获得各指定单颗卫星对应的法方程要素。

步骤4，根据指定单颗北斗igso/meo卫星对应的法方程要素，逐卫星构建与频率对应的法方程，整体求解法方程得各频率对应的待估系数矩阵x＝[a1a2a3]^t，即三阶多项式改正模型系数，从而获得各单颗北斗igso/meo卫星各频率对应的三阶多项式改正模型。

本步骤进一步包括：

步骤4.1，各指定单颗北斗igso/meo卫星各频率均对应一组三阶多项式改正模型系数。

步骤4.2，各指定单颗北斗igso/meo卫星各频率按照步骤1～2处理。

步骤4.2，依据已有的综合解方法，整体求解各组三阶多项式改正模型系数。

由上述步骤有nixi＝li，xi表示频率i对应的一组三阶多项式改正模型系数，两边同时左乘n的逆矩阵，即可求解出三阶多项式改正模型系数。求解公式如下：

xi^prn＝ni^-1prnli^prn(10)

式(10)中：

i表示北斗系统的频率；

prn表示各颗北斗igso/meo卫星，即prn依次为c06、c07、c08、c09、c10、c11、c12、c13、c14；

ni^prn和li^prn为北斗igso/meo卫星prn频率i对应的法方程要素；

xi^prn为北斗igso/meo卫星prn频率i对应的三阶多项式改正模型系数。

粗差剔除部分的具体步骤如下：

步骤5，采用已有的循环迭代剔除粗差法，确定最终改正模型。

在各频率下分别进行如下步骤：

步骤5.1，采用步骤1～4的方法整体估计各指定北斗igso/meo卫星与卫星高度角相关的伪距码偏差三阶多项式改正模型系数，计算三阶多项式改正模型系数的均方差误差rms。

步骤5.2，设置10倍rms的限制，即阈值，针对各卫星逐历元迭代剔除部分粗差数据，具体方法如下：

首先，逐卫星逐历元，利用估计的三阶多项式改正模型系数以及当前历元和当前参考历元的卫星高度角，根据公式(3)和(4)估计当前历元和当前参考历元的组合观测值之差，记为接着，将利用原始观测值按照公式(1)计算的当前历元和当前参考历元的实际mpi组合观测值做差，所得差值记为δmpi；最后，将与δmpi做差，若该差值绝对值小于阈值，则视该当前历元观测数据为有效数据；反之，剔除该当前历元观测数据，即该当前历元观测数据不参与重新估计三阶多项式改正模型系数的估计过程。

步骤5.3，由步骤5.2逐卫星逐历元处理后，得到新的法方程要素和按照公式(10)解算法方程，重新估计待估系数，即更新各卫星各频率的三阶多项式改正模型系数。

实施例中，依次设置10、7、5、3倍rms为阈值，并重复子步骤5.1～5.3，更新伪距码偏差三阶多项式改正模型系数，并最终确定以3倍rms为阈值的粗差剔除数据策略的解算结果为最终改正模型系数。

作为优选，所述的igso/meo卫星端伪距码偏差改正模型采用三阶多项式函数，为精密定位用户提供对伪距毫米级的改正量。

本发明基于多路径观测值的特性，通过动态设定参考历元，利用与相同弧段内的参考历元做历元间差分，准确地消除整周模糊度参数。同时，采取各历元依据高度角定权、不同弧段依据历元个数定权、不同测站依据先验值定权的定权策略，灵活合理地确定数据处理过程中的权重。此外，在建模求解改正模型系数时还利用了法方程综合解算，便于各颗卫星在不同测站不同弧段的有效衔接，成功规避了所有卫星观测值之和为零的约束条件中存在的忽略不同弧段可视卫星存在差异的客观问题，从而精确地得到各颗卫星伪距码偏差改正模型，保证其在模糊度固定中的实用性。

本文中所描述的具体实施例中的三阶多项式改正模型仅仅是对本发明思想做举例说明。本发明所述技术领域的技术人员可以对所描述的具体实施案例做各种各样的修改或补充或采用其他的函数模型替代，但并不会偏离本发明的思想或者超越所附权利要求书所定义的范围。