基于相时延测量的皮秒级精度的卫星干涉时间传递方法

1.本发明属于卫星授时技术领域，具体涉及一种基于相时延测量的皮秒级精度的卫星干涉时间传递方法。背景技术：2.目前常用的时间传递技术主要可分为星基系统和陆基系统两大类。星基系统使用卫星开展时间传递，主要使用的有国际导航卫星系统(gnss)、区域导航卫星系统(rnss)，其中我国的系统包括有北斗导航卫星系统 (bds)和中国区域导航试验系统(caps)。陆基授时系统又可分为无线电授时和有线授时两种方式，其中无线电授时系统主要包括bpl长波授时系统、bpm短波授时系统以及bpc低频时码授时系统等，有线授时系统主要有ntp互联网授时系统、pstn公网电话授时系统等，另外新研发的基于光纤链路开展的光纤时间传递也在不断扩展应用。各授时技术的授时精度如下表1所示。3.表1各类授时技术的精度对比[0004][0005][0006]对于长基线(百公里-千公里)应用，卫星共视和卫星授时可提供纳秒级别的时间传递精度，这也是国际权度局(bipm)组织国际时间传递所使用的主要方法。[0007]目前亚纳秒级的时间传递技术主要有基于gnss的载波共视技术、精密单点定位(ppp)时间传递技术和基于光纤的时间传递技术。基于gnss 的ppp时间传递和载波共视技术精度可达亚纳秒级别，无法满足皮秒、亚皮秒级别的高精度的时间传递需求。光纤时间传递精度高，但作为一种路基有线技术，必须依赖铺设的光缆，以光缆作为传输介质，应用范围受到局限。目前尚缺少一种适用于短基线(十公里级以内)的亚皮秒级无线时间传递方法。技术实现要素：[0008]为了解决现有技术中存在的上述问题，本发明提供了一种基于相时延测量的皮秒级精度的卫星干涉时间传递方法。本发明要解决的技术问题通过以下技术方案实现：[0009]第一方面，本发明提供的一种基于相时延测量的皮秒级精度的卫星干涉时间传递方法包括：[0010]获取两个测站共同观测同一卫星的观测数据；[0011]根据gnss数据，使用ppp解算出钟差初值；[0012]根据卫星坐标，测站坐标、各项修正量和钟差初值，计算理论时延；[0013]对观测数据依次进行相关处理和相关后处理，得到干涉群时延结果以及干涉相时延结果；[0014]根据干涉群时延结果，解算出模糊度；[0015]根据解算出的模糊度、干涉相时延结果以及理论时延，计算两个测站间的最终钟差，实现时间传递。[0016]可选的，理论时延表示为：[0017][0018]其中：δt′ab指ab两站的钟差初值，是通过ppp时间传递获得的较低精度的纳秒级钟差产品；c是真空中的光速，das、dbs指卫星s分别到a、 b站传播距离，指a、b站观测卫星s的对流层时延，指a、 b站观测卫星s的电离层时延，指a、b站观测卫星s的引力时延，指a、b站的系统硬件时延，oth指其它改正。[0019]可选的，对观测数据进行相关处理包括：[0020]对两个测站的观测数据进行理论时延补偿以及相位补偿；[0021]对理论时延补偿和相位补偿后的观测数据进行傅里叶变换；[0022]将傅里叶变换后的两测站的观测数据进行共轭相乘，得到互相关函数；[0023]对互相关函数进行转化，得到关于频率的相位函数。[0024]可选的，对观测数据进行相关后处理包括：[0025]根据相位函数，计算干涉群时延结果以及干涉相时延结果。[0026]可选的，理论时延补偿以及相位补偿后的两个测站的观测数据表示为：[0027]wa(t)＝va(t+δ)e-iθ[0028]wb(t)＝vb(t)[0029]θ＝2πfloδ[0030]其中，wa(t)、wb(t)分别表示a、b两个测站补偿后的信号，flo为变频器的本振频率；[0031]傅里叶变换后的信号表示为：[0032][0033][0034]互相关函数表示为：[0035][0036]相位函数表示为：[0037]可选的，根据解算出的模糊度、干涉相时延结果以及理论时延，计算两个测站间的最终钟差，实现时间传递包括：[0038]根据解算出的模糊度、干涉相时延结果以及理论时延，计算出干涉测量的时差值；[0039]根据干涉测量的时差值，引入观测方程，计算测站间的最终钟差，实现时间传递。[0040]其中，干涉群时延结果表示为：[0041][0042]干涉相时延结果表示为：[0043][0044]模糊度表示为：[0045][0046]干涉测量的时差值表示为：[0047]δt＝δ+τphase；[0048]其中，δt指干涉测量的时差值。[0049]其中，观测方程表示为：[0050][0051]其中：δtab指ab两站的最终钟差，是用于时间传递的修正量。[0052]第二方面，本发明提供的一种基于相时延测量的皮秒级精度的卫星干涉时间传递系统包括：两个测站，两个测站中包括一个主站；主站，用于：[0053]获取两个测站共同观测同一卫星的观测数据；[0054]根据gnss数据，使用ppp解算出钟差初值；[0055]将卫星坐标，测站坐标、各项修正量和钟差初值，计算理论时延；[0056]对观测数据依次进行相关处理和相关后处理，得到干涉群时延结果以及干涉相时延结果；[0057]根据干涉群时延结果，解算出模糊度；[0058]根据解算出的模糊度、干涉相时延结果以及理论时延，计算两个测站间的最终钟差，实现时间传递。[0059]其中，理论时延表示为：[0060][0061]其中：δt′ab指ab两站的钟差初值，是通过ppp时间传递获得的较低精度的纳秒级钟差产品；c是真空中的光速，das、dbs指卫星s分别到a、 b站传播距离，指a、b站观测卫星s的对流层时延，指a、b站观测卫星s的电离层时延，指a、b站观测卫星s的引力时延，指a、b站的系统硬件时延，oth指其它改正。[0062]本发明提供的一种基于相时延测量的皮秒级精度的卫星干涉时间传递方法，通过获取两个测站共同观测同一卫星的观测数据；基于gnss数据解算出钟差初值；根据卫星坐标，测站坐标、各项修正量和钟差初值，计算理论时延；对观测数据依次进行相关处理和相关后处理，得到干涉群时延结果以及干涉相时延结果；根据干涉群时延结果，解算出模糊度；根据理论时延、干涉相时延结果以及解算出的模糊度，计算两个测站间的钟差，实现时间传递。本发明中地面测站不主动发射信号，被动接收轨道已知的卫星的任意下行信号，无需解调信号也无需在卫星上搭载原子钟载荷，通过群时延辅助进行干涉相时延测量，利用已知的轨道信息和测站坐标解算站间钟差，实现时间传递，精度可达亚皮秒级。[0063]以下将结合附图及实施例对本发明做进一步详细说明。附图说明[0064]图1是本发明实施例提供的一种基于相时延测量的皮秒级精度的卫星干涉时间传递方法的流程示意图；[0065]图2是本发明实施例提供的基于相时延测量的皮秒级精度的卫星干涉时间传递原理图；[0066]图3是本发明实施例提供的干涉相时延时间传递结果图；[0067]图4a是本发明实施例提供的干涉相时延时间传递结果的频率稳定度示意图；[0068]图4b是本发明实施例提供的干涉相时延时间传递结果的时间稳定度示意图。具体实施方式[0069]下面结合具体实施例对本发明做进一步详细的描述，但本发明的实施方式不限于此。[0070]如图1所示，本发明提供的一种基于相时延测量的皮秒级精度的卫星干涉时间传递方法包括：[0071]s1，获取两个测站共同观测同一卫星的观测数据；[0072]s2，根据gnss数据，使用ppp解算出钟差初值；[0073]本步骤可以根据两个测站的观测数据，使用ppp时间传递解算获得钟差初值。[0074]传播时延表示为：[0075][0076]s3，根据卫星坐标，测站坐标、各项修正量和钟差初值，计算理论时延；[0077]理论时延表示为：[0078][0079]其中：δt′ab指ab两站的钟差初值，是通过ppp时间传递获得的较低精度的纳秒级钟差产品；c是真空中的光速，das、dbs指卫星s分别到a、 b站传播距离，指a、b站观测卫星s的对流层时延，指a、 b站观测卫星s的电离层时延，指a、b站观测卫星s的引力时延，指a、b站的系统硬件时延，oth指其它改正。[0080]参考图2，测站a、b同时观测卫星s，由于测站坐标和卫星轨道已知，可以计算信号从卫星到达测定轨主站的理论传播时延。[0081]值得说明的是：引力时延、质心修正、相位中心修正等均可精确修正，这里不做详细描述。其他影响较大的误差项分析如下。[0082](1)传播路径误差[0083]卫星到测站的传播几何距离可通过卫星和站的坐标计算，以站a为例，令a站坐标为(xa,ya,za),卫星s坐标为(xs,ys,zs)：[0084][0085]其中，站坐标和卫星轨道均可能存在误差。[0086](2)站坐标误差[0087]在静态接收机授时过程中，一般通过已知的精确坐标代入求解计算。实际上，站坐标会受到固体引力潮汐、大陆板块迁移、极移、海潮等各种因素的影响。陆板块迁移的年变化为厘米级，尚不需考虑。可是，地球固体潮汐和海洋潮汐较大，一般通过经验模型精确改正。[0088](3)轨道误差[0089]本发明以地球同步轨道(geo)卫星为例。目前可提供精密轨道的geo 卫星主要为北斗的geo卫星，轨道精度可达米级。国家授时中心的转发测定轨技术可为任意通信geo卫星提供米级轨道，其观测的卫星也可作为本专利观测的卫星。[0090]由观测方程可知，卫星与两测站的距离差经过一次差分，可部分消除卫星轨道引入的误差。卫星轨道引入的误差可由下式计算。[0091](eb-ea)·δs/c(4)[0092]其中ea、eb分别是a站、b站指向卫星的单位向量，δs为轨道误差。由于geo卫星轨道高度约为三万六千公里，卫星轨道高度远大于基线长度时，为便于计算可粗略认为轨道引起的误差为：[0093](dab/轨道高度)·δs/c(5)[0094]其中，dab为a、b之间的距离。[0095]要求100公里的基线轨道引入的误差优于1ps，则粗略估计geo卫星轨道精度应优于0.1米，meo卫星轨道精度应优于6.6厘米，1000公里leo 卫星轨道精度应优于3毫米，可暂不考虑。[0096](4)对流层时延[0097]对流层模型可通过模型修正，以国际gnss服务(igs)使用的gmf 模型为例。[0098]测站观测卫星的对流层时延与卫星高度角有关。[0099][0100]其中τtro为对流层延迟，即公式1中的或es为卫星高度角，表示干燥大气引起的对流层偏差中的干分量；表示水汽分子的偶极距所引起的对流层偏差的湿分量；映射函数也相应的分为干、湿映射函数，记作mdry(es)和mwet(es)。[0101][0102]其中m(es)指mdry(es)或mwet(es)，a、b、c是远小于1的常数，且分为干气部分和湿气部分，分别对应于干分量投影函数和湿分量投影函数。[0103]当基线较短时，两站观测卫星的路径基本一致，卫星高度角十分接近，使得两站观测卫星的对流层时延有很大的相关性，可通过差分消除部分影响。若基线南北分布，则其观测卫星的高度角差异最大。[0104]最差情况下令测站与卫星位于同一经度，易知南北分布的100公里基线的纬度差为0.9°，观测geo卫星的高度角之差在赤道处达到最大约1.06°，由于ppp解算的对流层产品精度优于1cm，则高度角大于41°时经过差分后，对流层误差的残余影响优于1ps。若测站为东西分布，则易知其观测geo 卫星的高度角之差在赤道处最大为0.017°，高度角大于17°时经过差分后，对流层误差的残余影响优于0.1ps。[0105](5)电离层时延[0106]电离层模型可通过双频观测改正，也可通过模型改正。下文以mslm (modified single-layer model)模型为例。[0107]电离层延迟可通过下式计算，[0108][0109]stec＝mf(z)·vtec(9)[0110]mf(z)＝1/cos(arcsin(r/(r+h)·cos(a(90°‑z))))(10)[0111]其中τion为电离层延迟，即公式1中的或freq为观测频点， stec指倾斜电子浓度，vtec指垂直电子浓度，mf为投影函数，z为卫星高度角，r为地球半径，h为电离层高度，a为常数0.9782。[0112]电离层延迟可通过双频观测修正，在没有双频观测的时候可采用使用 igs或国际gnss检测评估系统(igmas)提供的电离层网格产品修正。电离层网格产品的精度为2-8tecu(总电子数单位)，假设单频信号为c 波段(4-8ghz)，电离层延迟比在gps l1状况下缩小6至25倍，单站电离层延迟的修正精度为0.17ns至0.67ns。[0113]与对流层类似，短基线情况下可通过站间差分消除部分电离层影响。测站南北分布时，100公里基线下电离层的映射函数差异小于0.04，则电离层的修正的残余影响8ps至34ps。测站东西分布时，100公里基线下电离层的映射函数差异小于0.0008，则电离层的修正的残余影响0.1ps至0.5ps。测站南北分布时，建议通过双频观测的方法，利用电离层时延的平方反比特性消除电离层影响。[0114](6)系统时延[0115]该系统存在设备引起的硬件时延，可通过移动站进行连线干涉的方式标定。即使用一个移动站c，在a站使用其时频信号开展相时延干涉时差测量；再将该移动站c移至b站，使用其时频信号开展相时延干涉时差测量；两次测量结果互差后可消除移动站c引入的误差，得到a-b基线的系统硬件时延，具体解算过程不再详述。[0116]令两测站a、b接收到的信号为va(t)、vb(t)，且均为同一信号源下发的信号，不妨令b站为基准，则存在时差δt，令va(t+δt)＝vb(t)。相关处理即是求解两测站钟差δt的过程。由于δt较大时相关处理难以获得准确的结果，所以要先用理论时延δ进行修正，使va(t+δ)与vb(t)较为接近，则可转为测量较小的残余时延τ＝δt-δ。[0117]相关处理首先需要计算理论时延进行补偿，补偿的理论时延δ是根据卫星坐标、各项修正量和精度较差的钟差初值，即初步钟差计算的。[0118]s4，对观测数据依次进行相关处理和相关后处理，得到干涉群时延结果以及干涉相时延结果；[0119]本步骤可以对两个测站的观测数据进行理论时延补偿以及相位补偿；对理论时延补偿和相位补偿后的观测数据进行傅里叶变换；将傅里叶变换后的两测站的观测数据进行共轭相乘，得到互相关函数；对互相关函数进行转化，得到关于频率的相位函数；根据相位函数，计算干涉群时延结果以及干涉相时延结果。[0120]值得说明的是：对观测数据进行理论时延补偿和相位补偿(条纹旋转)，令补偿后的信号为wa(t)、wb(t)，其中flo为变频器的本振频率。wa(t)、wb(t) 分别表示a、b两个测站补偿后的信号；[0121]wa(t)＝va(t+δ)e-iθ(11)[0122]wb(t)＝vb(t)(12)[0123]θ＝2πfloδ(13)[0124]对补偿后的信号进傅里叶变换后的信号表示为：[0125][0126][0127]互相关函数为：[0128][0129]据互相关函数获得相位值即相位函数，是关于频率的函数。[0130]相位函数表示为：[0131][0132]则可以获得残余时延τ的群时延结果τgroup和相时延结果τphase。[0133][0134][0135]s5，根据干涉群时延结果，解算出模糊度；[0136]其中，n为模糊度，在一定时间内为一个常整数，可根据群时延计算。[0137][0138]s6，根据解算出的模糊度、干涉相时延结果以及理论时延，计算两个测站间的钟差，实现时间传递。[0139]本发明可以根据解算出的模糊度、干涉相时延结果以及理论时延，计算出干涉测量的时差值；根据干涉测量的时差值，引入观测方程，计算测站间的最终钟差，实现时间传递。[0140]其中，观测方程表示为：[0141][0142]其中：δtab指ab两站的最终钟差，是用于时间传递的修正量，δt指干涉测量的时差值，是通过相关处理获得的测量值，c是真空中的光速，das指卫星s分别到a、b站传播距离，指a、b站观测卫星s的对流层时延，指a、b站观测卫星s的电离层时延，指a、 b站观测卫星s的引力时延，指a、b站的系统硬件时延，oth指其它改正，包括天线相位中心改正等。[0143]则两站的时差值δt则可根据相时延和理论时延获得，两个测站间的时差值表示为：[0144]δt＝δ+τphase(22)[0145]则两站的最终钟差值δtab则可根据观测方程获得，两个测站间的钟差表示为：[0146][0147]值得说明的是：两个测站通过抛物面天线观测同一颗卫星，测站基于本地时间采集卫星信号，通过相关处理可获得卫星信号到达两个测站天线的时间差，通过卫星轨道和测站坐标可计算卫星到测站的实际几何距离差，二者差分后可获得两个测站的钟差，从而实现时间传递。[0148]本发明提供的一种基于相时延测量的皮秒级精度的卫星干涉时间传递系统包括：两个测站，两个测站中包括一个主站；主站，用于：[0149]获取两个测站共同观测同一卫星的观测数据；[0150]根据gnss数据，使用ppp解算出钟差初值；[0151]将卫星坐标，测站坐标、各项修正量和钟差初值，计算理论时延；[0152]对观测数据依次进行相关处理和相关后处理，得到干涉群时延结果以及干涉相时延结果；[0153]根据干涉群时延结果，解算出模糊度；[0154]根据解算出的模糊度、干涉相时延结果以及理论时延，计算两个测站间的最终钟差，实现时间传递。[0155]下面结合实验验证本发明的最终效果。[0156]为验证本方法，2021年11月17日-18日组织了验证试验。由于该方法现处于预研阶段，试验系统中使用了同一个原子钟，因此未使用ppp时间传递结果作为初值。[0157]试验使用了两面3.7米口径的抛物面天线，基线长度为19.93米；观测卫星为北斗c01星，星下点位于144.5°e；观测频段为c波段；使用difx 相关处理软件系统和aips相关后处理软件进行相关处理和相关后处理。[0158]试验中的轨道数据为igmas提供的事后精密轨道。试验没有进行大气延迟修正和设备时延修正。试验结果如下图3所示，存在的起伏可能是设备和大气时延引入的，在实际使用中可进一步修正。[0159]时差结果的修正阿伦方差和时间方差如图4a和图4b所示，频率稳定度达到了7.35e-13(1s)，时间稳定度达到了4.25e-13秒(1s)，达到了亚皮秒级水平。[0160]以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。