一种基于计时残差的综合脉冲星时确定方法及系统

1.本发明涉及脉冲星时技术领域，尤其涉及一种基于计时残差的综合脉冲星时确定方法及系统。背景技术：2.在定位导航授时(pnt)体系中，授时决定着定位和导航能否成功，而高精度授时更是以高精度、高稳定度守时为基础，保证通信、金融等各行各业的平稳、安全运转，其中，脉冲星属于中子星，其具有磁场强、体积小、密度大、自转高速的特点，脉冲星每自转一个周期，地面测站就会接收到脉冲星两极发射的高能射线，相比传统的原子钟守时方法，脉冲星守时技术在长期稳定度上具有可观的优势，尤其是毫秒脉冲星，有些自转周期变化率达到10-19～10-21,被誉为自然界最稳定的钟，而脉冲星计时残差的稳定度可以用于衡量脉冲星发射高能脉冲信号的稳定度。3.目前，在对综合脉冲星时进行稳定度评估时，传统的脉冲星时生成方法主要采用vondrak滤波预处理、三次样条插值方法对其脉冲星计时稳定度进行分析，但是在利用传统方法评估综合脉冲星时过程中，不仅无法保证脉冲星计时残差的所有点数对齐，使得计时总残差和总误差难以取信，而且仅采用vondrak滤波进行数据处理，对计时稳定度贡献有限。技术实现要素：4.本发明提供了一种基于计时残差的综合脉冲星时确定方法及系统，解决的技术问题是，传统的脉冲星时生成方法不仅无法保证脉冲星计时残差的所有点数对齐，使得计时总残差和总误差难以取信，而且仅采用vondrak滤波进行数据处理，对计时稳定度贡献有限。5.为解决以上技术问题，本发明提供了一种基于计时残差的综合脉冲星时确定方法及系统。6.第一方面，本发明提供了一种基于计时残差的综合脉冲星时确定方法，所述方法包括以下步骤：7.根据预设的最优筛选条件选取至少两颗最优脉冲星，并获取每颗所述最优脉冲星的初始计时残差和初始误差；8.对每颗所述最优脉冲星的初始计时残差依次进行预处理和滤波，得到对应的平滑计时残差；9.选取基准脉冲星，并根据所述基准脉冲星的观测点数，对每颗所述最优脉冲星的平滑计时残差和初始误差进行等间距降采样，得到对应的降采样后的计时残差和降采样后的初始误差；10.根据所述降采样后的计时残差和所述降采样后的初始误差，得到每颗最优脉冲星的稳定度；11.根据每颗最优脉冲星的稳定度，得到综合加权残差和综合加权误差；12.根据所述综合加权残差、所述综合加权误差和所述基准脉冲星的儒略日，得到综合脉冲星时稳定度。13.在进一步的实施方案中，所述预设的最优筛选条件包括：14.所述最优脉冲星为毫秒脉冲星，且所述最优脉冲星的观测时间跨度须满足预设观测时间跨度阈值，观测点数须满足预设观测点数阈值。15.在进一步的实施方案中，所述对每颗最优脉冲星的初始计时残差依次进行预处理和滤波，得到对应的平滑计时残差的步骤包括：16.利用绝对中位数偏差法对每颗最优脉冲星的初始计时残差进行粗差点提取；17.对每颗最优脉冲星的所述粗差点进行置零处理，并利用分段三次hermite插值的方法，对所述粗差点进行插值处理，得到预处理后的初始计时残差；18.利用vondrak滤波对预处理后的初始计时残差进行滤波降噪处理，得到每颗最优脉冲星的平滑计时残差。19.在进一步的实施方案中，所述综合加权残差和所述综合加权误差的计算公式分别为：[0020][0021][0022]其中，[0023][0024]式中，res总表示综合加权残差；resi′表示第i颗最优脉冲星的降采样后的计时残差；err总表示综合加权误差；erri′表示第i颗最优脉冲星的降采样后的初始误差；wi表示第i颗最优脉冲星的权值；σzi表示第i颗最优脉冲星的稳定度。[0025]在进一步的实施方案中，所述综合脉冲星时稳定度的计算公式为：[0026][0027]式中，σ总(τ)表示综合脉冲星时稳定度；τ表示时间尺度；c3表示三次多项式系数；res总表示综合加权残差；x(mjd3)表示三次多项式拟合函数；mjd3表示基准脉冲星的儒略日；err总表示综合加权误差。[0028]在进一步的实施方案中，所述选取基准脉冲星的步骤包括：根据预设的基准筛选条件对所有所述最优脉冲星进行筛选，得到基准脉冲星，并将基准脉冲星观测儒略日作为基准儒略日。[0029]在进一步的实施方案中，所述预设的基准筛选条件包括：选取观测点数最少的最优脉冲星作为基准脉冲星。[0030]第二方面，本发明提供了一种基于计时残差的综合脉冲星时确定系统，所述系统包括：[0031]最优脉冲星获取模块，用于根据预设的最优筛选条件选取至少两颗最优脉冲星，并获取每颗所述最优脉冲星的初始计时残差和初始误差；[0032]脉冲星数据处理模块，用于对每颗所述最优脉冲星的初始计时残差依次进行预处理和滤波，得到对应的平滑计时残差；还用于选取基准脉冲星，并根据所述基准脉冲星的观测点数，对每颗所述最优脉冲星的平滑计时残差和初始误差进行等间距降采样，得到对应的降采样后的计时残差和降采样后的初始误差；根据所述降采样后的计时残差和所述降采样后的初始误差，得到每颗最优脉冲星的稳定度；[0033]计时残差加权模块，用于根据每颗最优脉冲星的稳定度，得到综合加权残差和综合加权误差；[0034]综合脉冲星时确定模块，用于根据所述综合加权残差、所述综合加权误差和所述基准脉冲星的儒略日，得到综合脉冲星时稳定度。[0035]第三方面，本发明还提供了一种计算机设备，包括处理器和存储器，所述处理器与所述存储器相连，所述存储器用于存储计算机程序，所述处理器用于执行所述存储器中存储的计算机程序，以使得所述计算机设备执行实现上述方法的步骤。[0036]第四方面，本发明还提供了一种计算机可读存储介质，所述计算机可读存储介质中存储有计算机程序，所述计算机程序被处理器执行时实现上述方法的步骤。[0037]本发明提供了一种基于计时残差的综合脉冲星时确定方法及系统，所述方法通过对选取的最优脉冲星的初始计时残差进行预处理、vondrak滤波和降采样，同时对误差进行降采样，以根据降采样后的计时残差和误差获取各颗最优脉冲星的稳定度，根据各颗最优脉冲星的稳定度得到加权计时残差和加权误差，从而计算得到综合脉冲星时。与现有技术相比，该方法基于预处理-vondrak滤波-降采样的脉冲星计时残差处理，确定综合脉冲星时，同时利用加权平均脉冲星时算法，使得综合脉冲星时的稳定度比单颗脉冲星稳定度平均水平提高1个数量级以上，提高了综合脉冲星时的精度和可靠性，实用性强。附图说明[0038]图1是本发明实施例提供的基于计时残差的综合脉冲星时确定方法流程示意图；[0039]图2是本发明实施例提供的综合脉冲星时计算过程示意图；[0040]图3是本发明实施例提供的脉冲星j0437-4715通过mad法去除粗差点前后的计时残差对比示意图；[0041]图4是本发明实施例提供的脉冲星j0437-4715进行预处理前后的计时残差对比示意图；[0042]图5是本发明实施例提供的脉冲星j0437-4715计时残差在vondrak滤波前后的对比示意图；[0043]图6是本发明实施例提供的综合脉冲星时稳定度和最优脉冲星单独守时下的稳定度对比示意图；[0044]图7是本发明实施例提供的基于计时残差的综合脉冲星时确定系统框图；[0045]图8是本发明实施例提供的计算机设备的结构示意图。具体实施方式[0046]下面结合附图具体阐明本发明的实施方式，实施例的给出仅仅是为了说明目的，并不能理解为对本发明的限定，包括附图仅供参考和说明使用，不构成对本发明专利保护范围的限制，因为在不脱离本发明精神和范围基础上，可以对本发明进行许多改变。[0047]参考图1，本发明实施例提供了一种基于计时残差的综合脉冲星时确定方法，如图1所示，该方法包括以下步骤：[0048]s1.根据预设的最优筛选条件选取至少两颗最优脉冲星，并获取每颗所述最优脉冲星的初始计时残差和初始误差。[0049]具体地，本实施例根据预设的最优筛选条件，在ipta2019-release2中，从nanograv观测到的37颗脉冲星中筛选出若干最优脉冲星，并使用tempo2软件对每一颗所述脉冲星的星历文件和模型参数进行运算，获取每颗最优脉冲星的初始计时残差和初始误差，在本实施例中，所述预设的最优筛选条件包括：[0050](1)所述最优脉冲星为毫秒脉冲星；由于毫秒脉冲星自转稳定度高，计时残差小，因此，将其作为最优脉冲星的基本要求之一；[0051](2)所述最优脉冲星的观测时间跨度须满足预设观测时间跨度阈值，以便于后续计算该脉冲星的长期稳定度，并研究噪声特性，本实施例优先将预设观测时间跨度阈值设置为10yr(即10年)以上；[0052](3)所述最优脉冲星的观测点数(即toas)须满足预设观测点数阈值，以提高计时模型拟合精度，本实施例优先将预设观测点数阈值设置为8000以上。[0053]基于上述的最优筛选条件，本实施例从nanograv观测到的37颗脉冲星中筛选出如表1所示的8颗脉冲星，表1如下所示：[0054]表1[0055][0056]需要说明的是，虽然psr j0437-4715脉冲星的观测年份较长，但是其原始观测数据中的toas只有5302，远少于其他脉冲星，若将其一起选择，会对其他脉冲星降采样后的点数产生一定的不利影响，从而不利于综合脉冲星时的稳定度提高，因此，本实施例排除psr j0437-4715脉冲星。[0057]s2.对每颗所述最优脉冲星的初始计时残差依次进行预处理和滤波，得到对应的平滑计时残差。[0058]在一个实施例中，所述对每颗最优脉冲星的初始计时残差依次进行预处理和滤波，得到对应的平滑计时残差的步骤包括：[0059]利用绝对中位数偏差法对每颗最优脉冲星的初始计时残差进行粗差点提取；[0060]对每颗最优脉冲星的所述粗差点进行置零处理，并利用分段三次hermite插值的方法，对所述粗差点进行插值处理，得到预处理后的初始计时残差；[0061]利用vondrak滤波对预处理后的初始计时残差进行滤波降噪处理，得到每颗最优脉冲星的平滑计时残差。[0062]如图2所示，本实施例分别对这8颗最优脉冲星的初始计时残差进行预处理和vondrak滤波，其中，预处理步骤采用绝对中位数偏差法(mad)对脉冲星计时残差的粗差点进行提取，具体为：[0063][0064]式中，xj(j＝1,2,...,n)表示初始计时残差向量。[0065]若xj满足上式，则说明该初始计时残差超过粗差阈值，本实施例将该初始计时残差为粗差点，并将粗差点进行置零处理，然后利用分段三次hermite插值的方法，对所述粗差点进行插值处理，得到预处理后的初始计时残差，使得预处理后的初始计时残差更为合理且平滑，从而避免出现runge现象；在对初始计时残差进行预处理之后，本实施例还通过vondrak滤波对预处理后的初始计时残差进行滤波降噪处理，以滤除高频噪声。[0066]需要说明的是，虽然脉冲星j0437-4715不被纳入8颗最优脉冲星中，但它可以用于演示实验(其他8颗脉冲星的结果类似)，图3为脉冲星j0437-4715通过mad法去除粗差点前后的计时残差对比示意图；图4为脉冲星j0437-4715进行预处理前后的计时残差对比示意图，其中，在图4中，带有龙格现象的预处理曲线是在预处理中采用了三次样条插值替换了三次hermite插值生成的曲线；图5为脉冲星j0437-4715计时残差在vondrak滤波前后的对比示意图。[0067]s3.选取基准脉冲星，并根据所述基准脉冲星的观测点数，对每颗所述最优脉冲星的平滑计时残差和初始误差进行等间距降采样，得到对应的降采样后的计时残差和降采样后的初始误差。[0068]具体地，本实施例首先根据预设的基准筛选条件对所有所述最优脉冲星进行筛选，得到基准脉冲星，并将基准脉冲星观测儒略日作为基准儒略日；其中，所述预设的基准筛选条件包括：选取观测点数(toas)最少的最优脉冲星作为基准脉冲星，根据表1选出的最优脉冲星，本实施例筛选出toas最少的脉冲星psr j1643-1224作为基准脉冲星，该基准脉冲星的观测儒略日((modified julian date，mjd)范围，即mjd3作为基准儒略日。[0069]然后，本实施例将ni(i＝1,2,…8)记为第i颗最优脉冲星的toas，特别地，记n＝min{ni}＝n3，其中，n3表示表1中的第3颗基准脉冲星psr j1643-1224的toas，同时，本实施例记resi和erri分别为第i颗最优脉冲星的平滑计时残差和初始误差，分别表示为：[0070][0071][0072]根据所述基准脉冲星的观测点数，对每颗所述最优脉冲星的平滑计时残差和初始误差进行等间距降采样，得到对应的降采样后的计时残差和降采样后的初始误差，本实施例通过降采样将各颗最优脉冲星的计时残差和误差点数对齐，并选择基准脉冲星的儒略日作为数据处理时的统一条件，形成了一种可解释(可信)的综合脉冲星时，其中，降采样后的计时残差和降采样后的初始误差分别为：[0073][0074][0075]其中，resi′表示第i颗最优脉冲星的降采样后的计时残差；erri′表示第i颗最优脉冲星的降采样后的初始误差；《x》表示对x进行四舍五入。[0076]s4.根据所述降采样后的计时残差和所述降采样后的初始误差，得到每颗最优脉冲星的稳定度。[0077]为了消除线性频漂的影响，脉冲星的稳定度通常采用σz方差(又称hadamard方差)评估，与allan方差等统计量不同，σz计算方式无需要求残差以相等间隔分布，对分析脉冲星数据尤为方便，σz的计算方式是使用三次多项式函数拟合时间尺度为τ的(拟合后)计时残差，并利用三次多项式系数c3确定其数值，具体为：记数据总观测时间为t＝tn-t1，并划分为等间隔τ的子序列，设定t0为参考时间，对每个子序列用三次多项式函数x(t)进行最小二乘拟合，即：[0078]x(t)＝c0+c1(t-t0)+c2(t-t0)2+c3(t-t0)3[0079]其中，每颗最优脉冲星的稳定度σzi(τ)的计算公式具体为：[0080][0081]式中，ti(i＝1,2,…,n)表示第i颗最优脉冲星的观测时间(单位为儒略日)，由于根据上述方法选取了8颗最优脉冲星，因此，n＝8；xi表示第i颗最优脉冲星的计时残差；σi表示第i颗最优脉冲星的误差；尖括号表示在所有子序列上，以c3不确定度的平方倒数为权做加权平均。[0082]s5.根据每颗最优脉冲星的稳定度，得到综合加权残差和综合加权误差。[0083]在一个实施例中，所述综合加权残差和所述综合加权误差的计算公式分别为：[0084][0085][0086]其中，[0087][0088]式中，res总表示综合加权残差；resi′表示第i颗最优脉冲星的降采样后的计时残差；err总表示综合加权误差；erri′表示第i颗最优脉冲星的降采样后的初始误差；wi表示第i颗最优脉冲星的权值；σzi表示第i颗最优脉冲星的稳定度。[0089]s6.根据所述综合加权残差、所述综合加权误差和所述基准儒略日，得到综合脉冲星时稳定度。[0090]在一个实施例中，所述综合脉冲星时稳定度的计算公式为：[0091][0092]式中，σ总(τ)表示综合脉冲星时稳定度；τ表示时间尺度；c3表示误差；res总表示综合加权残差；x(mjd3)表示三次多项式拟合函数；mjd3表示基准脉冲星的儒略日；err总表示综合加权误差。[0093]本实施例通过基准脉冲星的儒略日mjd3、综合加权残差res总和综合加权误差err总计算得到8颗最优脉冲星组成的综合脉冲星时稳定度，图6为综合脉冲星时稳定度和8颗最优脉冲星单独守时下的稳定度对比示意图，在图6中，横坐标表示计算稳定度的时间尺度，数值上等同于每个子序列的长度(如横坐标取0时，即时间尺度τ取对数的值为0，故τ为1yr，则在σz稳定度计算时，应把观测数据以1yr为单位分组拟合)；纵坐标表示σz稳定度数值的对数，由图6可知，本实施例基于经典加权平均算法的综合脉冲星时稳定度σ总(τ)，在22yr的时间尺度下达到了1e-15.9，表2为8颗最优脉冲星各自计算出的长、短期稳定度，表2如下所示：[0094]表2[0095][0096]表2中的长、短期稳定度的时间尺度分别对应图6中8颗最优脉冲星单独守时曲线的左端点和右端点，它们均不是一个具体数值，而是一个涵盖两倍跨度的范围，这是因为每颗脉冲星的观测时间均不相同，而σz稳定度的计算是在时域以2为底的对数离散分布的；在图6中，由于在22yr下稳定度最高的单颗最优脉冲星psr j1713+0747只达到了1e-15.0(psr j1909-3744的观测时间只达22yr的一半，未被考虑在内)，而图6中综合脉冲星时的22yr稳定度为1e-15.9，由此可以得出结论，基于本实施例提供的计时残差数据处理算法获取的综合脉冲星时稳定度比单颗脉冲星普遍高出1个数量级。[0097]本发明实施例提供了一种基于计时残差的综合脉冲星时确定方法，所述方法对筛选出的每颗最优脉冲星的初始计时残差依次进行预处理和滤波，得到对应的平滑计时残差，并对每颗所述最优脉冲星的平滑计时残差和初始误差进行等间距降采样，以根据降采样后的计时残差、降采样后的初始误差和基准脉冲星的儒略日，确定综合脉冲星时稳定度。与现有技术相比，本实施例采用基于预处理-vondrak滤波-降采样的脉冲星计时残差处理方法，同时利用经典加权平均算法对8颗最优脉冲星处理后的计时残差建立了综合脉冲星时，实现了在22yr的稳定度比稳定度最高的单颗脉冲星提高了0.9个数量级，比该时间尺度下单颗脉冲星稳定度的平均水平提高1个数量级以上的技术效果，对于脉冲星守时具有重要的实际意义。[0098]需要说明的是，上述各过程的序号的大小并不意味着执行顺序的先后，各过程的执行顺序应以其功能和内在逻辑确定，而不应对本技术实施例的实施过程构成任何限定。[0099]在一个实施例中，如图7所示，本发明实施例提供了一种基于计时残差的综合脉冲星时确定系统，所述系统包括：[0100]最优脉冲星获取模块101，用于根据预设的最优筛选条件选取至少两颗最优脉冲星，并获取每颗所述最优脉冲星的初始计时残差和初始误差；[0101]脉冲星数据处理模块102，用于对每颗所述最优脉冲星的初始计时残差依次进行预处理和滤波，得到对应的平滑计时残差；根据预先选取的基准儒略日，对每颗所述最优脉冲星的平滑计时残差和初始误差进行等间距降采样，得到对应的降采样后的计时残差和降采样后的初始误差；还用于根据所述降采样后的计时残差和所述降采样后的初始误差，得到每颗最优脉冲星的稳定度；[0102]计时残差加权模块103，用于根据每颗最优脉冲星的稳定度，得到综合加权残差和综合加权误差；[0103]综合脉冲星时确定模块104，用于根据所述综合加权残差、所述综合加权误差和所述基准儒略日，得到综合脉冲星时稳定度。[0104]关于一种基于计时残差的综合脉冲星时确定系统的具体限定可以参见上述对于一种基于计时残差的综合脉冲星时确定方法的限定，此处不再赘述。本领域普通技术人员可以意识到，结合本技术所公开的实施例描述的各个模块和步骤，能够以硬件、软件或者两者结合来实现。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用使用不同方法来实现所描述的功能，但是这种实现不应认为超出本技术的范围。[0105]本发明实施例提供了一种基于计时残差的综合脉冲星时确定系统，所述系统的脉冲星数据处理模块采用基于预处理-vondrak滤波-降采样的脉冲星计时残差处理方式，计时残差加权模块通过加权平均算法对最优脉冲星处理后的计时残差和误差进行处理，从而确定综合脉冲星时稳定度。与现有技术相比，本技术对后续的定位导航授时可以进行有效的指导，提高了脉冲星计时的可靠性和可解释性，具有适用范围广、实用性强、导航效率高等优点。[0106]图8是本发明实施例提供的一种计算机设备，包括存储器、处理器和收发器，它们之间通过总线连接；存储器用于存储一组计算机程序指令和数据，并可以将存储的数据传输给处理器，处理器可以执行存储器存储的程序指令，以执行上述方法的步骤。[0107]其中，存储器可以包括易失性存储器或非易失性存储器，或可包括易失性和非易失性存储器两者；处理器可以是中央处理器、微处理器、特定应用集成电路、可编程逻辑器件或其组合。通过示例性但不是限制性说明，上述可编程逻辑器件可以是复杂可编程逻辑器件、现场可编程逻辑门阵列、通用阵列逻辑或其任意组合。[0108]另外，存储器可以是物理上独立的单元，也可以与处理器集成在一起。[0109]本领域普通技术人员可以理解，图8中示出的结构，仅仅是与本技术方案相关的部分结构的框图，并不构成对本技术方案所应用于其上的计算机设备的限定，具体的计算机设备可以包括比图中所示更多或更少的部件，或者组合某些部件，或者具有相同的部件布置。[0110]在一个实施例中，本发明实施例提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现上述方法的步骤。[0111]本发明实施例提供的一种基于计时残差的综合脉冲星时确定方法及系统，其一种基于计时残差的综合脉冲星时确定方法实现了基于预处理-vondrak滤波-降采样的脉冲星计时残差处理，同时利用经典加权平均脉冲星时算法，得到具有高可信度的综合脉冲星时稳定度，可用于为脉冲星导航提供数据支持。[0112]在上述实施例中，可以全部或部分地通过软件、硬件、固件或者其任意组合来实现。当使用软件实现时，可以全部或部分地以计算机程序产品的形式实现。所述计算机程序产品包括一个或多个计算机指令。在计算机上加载和执行所述计算机程序指令时，全部或部分地产生按照本发明实施例所述的流程或功能。所述计算机可以是通用计算机、专用计算机、计算机网络、或者其他可编程装置。所述计算机指令可以存储在计算机可读存储介质中，或者从一个计算机可读存储介质向另一个计算机可读存储介质传输，例如，所述计算机指令可以从一个网站站点、计算机、服务器或数据中心通过有线(例如同轴电缆、光纤、数字用户线)或无线(例如红外、无线、微波等)方式向另一个网站站点、计算机、服务器或数据中心进行传输。所述计算机可读存储介质可以是计算机能够存取的任何可用介质或者是包含一个或多个可用介质集成的服务器、数据中心等数据存储设备。所述可用介质可以是磁性介质，(例如，软盘、硬盘、磁带)、光介质(例如，dvd)、或者半导体介质(例如ssd)等。[0113]本领域技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。[0114]以上所述实施例仅表达了本技术的几种优选实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和替换，这些改进和替换也应视为本技术的保护范围。因此，本技术专利的保护范围应以所述权利要求的保护范围为准。