1.一种北斗导航卫星广播星历参数拟合方法，该方法步骤如下：步骤一：准备包含星历参考时刻的一段时间内的卫星位置和速度及对应的时间信息；

步骤二：根据星历参考时刻的位置和速度求出星历参考时刻的轨道半长轴相对于参考值的偏差 平近点角 偏心率 近地点幅角 倾角 以及周历元零时刻升交点经度步骤三：设定17个星历参数的求解范围，将求解星历参数的最优化问题转为等价的拉格朗日对偶问题，构建方程组；

步骤四：设定带求解参数初始值，使用Levenberg‑Marquardt算法对方程组迭代求解，完成对卫星广播星历参数的拟合。

2.根据权利要求1所述的一种北斗导航卫星广播星历参数拟合方法，其特征在于：所述步骤一具体过程如下：准备包含星历参考时刻的一段时间内的卫星位置和速度及对应的时间信息，星历参考时刻记为toe，第i个时刻的卫星位置和速度记为xi＝(ri,vi)，其中，ri＝(rx,i,ry,i,rz,i)，vi＝(vx,i,vy,i,vz,i)，为了能够求解出17个轨道参数xi的数据点数应不小于6。

3.根据权利要求1所述的一种北斗导航卫星广播星历参数拟合方法，其特征在于：所述步骤二具体过程如下：S21、将星历参数的求解简化为二体问题，则根据星历参考时刻的位置和速度求出星历参考时刻的轨道半长轴相对于参考值的偏差 平近点角 偏心率 近地点幅角 倾角 升交点经度S22、求周历元零时刻的格林尼治恒星时；

周历元零时刻的儒略日整数部分记为jd，单位为天，小数部分对应的时间记为ut，单位转为秒，周历元零时刻格林尼治恒星时GMS计算方法为：

2 ‑6 3

GMS＝24110.54841+8640184.812866·T+0.093104·T ‑6.2×10 ·T +

1.002737909350795·ut

S23、计算周历元零时刻的升交点经度；

式中，mod(GMS,86400)为GMS对86400的余数。

4.根据权利要求1所述的一种北斗导航卫星广播星历参数拟合方法，其特征在于：所述步骤三具体过程如下：设定17个星历参数的求解范围，将求解星历参数的最优化问题转为等价的拉格朗日对偶问题，构建方程组；

待求解的轨道参数记为向量形式:

p中元素含义依次为：参考时刻长半轴相对于参考值的偏差、长半轴变化率、参考时刻卫星平均角速度与计算值之差、参考时刻卫星平均角速度与计算值之差的变化率、参考时刻的平近点角、偏心率、近地点幅角、周历元零时刻计算的升交点经度、参考时刻的轨道倾角、升交点赤经变化率、轨道倾角变化率、轨道倾角的正弦调和改正项的振幅、轨道半径的正弦调和改正项的振幅、轨道半径的正弦调和改正项的振幅、轨道半径的余弦调和改正项的振幅、纬度幅角的正弦调和改正项的振幅、纬度幅角的余弦调和改正项的振幅；

时刻tk，由p到地心地固坐标系下卫星位置rk的函数关系记为：rk＝f(p,tk)，由x求p的最优化问题：min||F(p)‑b||

s.t.ci≤pi≤di i＝1,2,…,17

可转换为相应的拉格朗日对偶问题：

s.t.ci≤pi≤di i＝1,2,…,17

式中，b＝[r1,r2,...,rn]，b为1×3n的矩阵，n为数据点数，F(p)＝[f(p,t1),f(p,t2),...,f(p,tn)]，F(p)为1×3n的矩阵，pi为p中第i个元素，λi为拉格朗日乘子，ci和di为pi的上下界。

5.根据权利要求1所述的一种北斗导航卫星广播星历参数拟合方法，其特征在于：所述步骤四具体过程如下：S41、记 取初始点p0，取λi＝10

‑6 ‑3

,i＝1,2,...,17，终止控制常数ε，计算 k＝0,m0＝10 ，v＝10；，其中，S42、计算雅克比矩阵Jk，计算 构建增量正规方程

S43、求解增量正规方程得到δk；

S44、如果 则令pk+1＝pk+δk，若||δk||<εk，停止迭代，输出结果；否则令mk+1＝mk/v，转到步骤S42；

S45、如果 则令mk+1＝v·mk，重新解正规方程得到δk，返回步骤S42。

6.一种计算机可读存储介质，存储有计算机程序，所述计算机程序被处理器执行时，使得所述处理器执行如权利要求1至5中任一项所述方法的步骤。