发明名称 一种基于星间双向测距的时间同步方法
摘要 本发明公开了一种基于星间双向测距的时间同步方法,具体步骤为采用二次多项式对时钟钟差进行建模;考虑频率漂移变化的影响,星载原子钟模型采用3个状态分量,建立滤波的状态方程;针对星座中建立星间链路的两颗卫星,采用星间双向测距方法进行伪距测量,建立观测方程;根据系统的状态方程和测量方程,进行Kalman滤波;利用集中式卡尔曼滤波或分散式卡尔曼滤波对时间同步进行滤波;分析星间链路数对时间同步精度的影响,通过不同星间链路数下仿真结果的对比选择出性价比最高的方案,通过本发明中的集中滤波和分散滤波的对比,以及不同星间链路数目下的滤波精度的对比从而给出参考实验方案并提高了测量同步测度。
申请公布号 CN103957095A 申请公布日期 2014.07.30
申请号 CN201410205482.4 申请日期 2014.05.15
申请人 北京航空航天大学;中国空间技术研究院 发明人 王艳东;陈惠勇;赵欢;王虎妹;石俊;王世涛
分类号 H04L7/00(2006.01)I;G04F5/14(2006.01)I 主分类号 H04L7/00(2006.01)I
代理机构 代理人
主权项 一种基于星间双向测距时间同步方法,其特征在于,具体包括如下步骤:步骤一:在某一瞬时,时钟给出的时刻与标准时刻之间的差值称为钟差x(t),采用二次多项式对时钟钟差进行建模;步骤二:考虑频率漂移变化的影响,星载原子钟模型采用3个状态分量,建立滤波的状态方程;步骤三:针对星座中建立星间链路的两颗卫星,采用星间双向测距方法进行伪距测量,建立观测方程;步骤四:根据系统的状态方程和测量方程,进行Kalman滤波;步骤五:利用集中式卡尔曼滤波或分散式卡尔曼滤波对时间同步进行滤波:其中,分散式卡尔曼滤波方法:以各个节点的时钟参数作为状态变量,利用N个滤波器分别对钟差进行估计,其中N为系统的节点总数,具体步骤如下:1)采用双向卫星时间频率传递方法TWSTFT进行地面系统时间同步,将某一配备有高精度原子钟的地面站作为系统时钟基准站,其他地面站与之同步;2)采用星间双向测距系统进行系统时间同步,将某一卫星作为系统时钟的基准站,其他卫星与之同步;3)分别完成地面系统和卫星系统间的时间同步,采用无线电双向测距法进行地面系统和卫星系统间的时间同步,其中,只要对于地面系统和卫星系统中的基准站进行时间同步即可;集中式卡尔曼滤波方法:由于系统内各节点绝对钟差的不可估计性,选取卫星1为基准节点,状态变量取为各节点与基准节点的时钟参数之差,维数为(N‑1)×3,N为系统节点数,进行卡尔曼滤波,估计系统内各节点的相对钟差,其中时间同步方法的选取与分散滤波方法一致;步骤六:分析星间链路数对时间同步精度的影响,通过不同星间链路数下仿真结果的对比选择出性价比最高的方案。其中上述步骤一中,在某一瞬时,时钟给出的时刻与标准时刻之间的差值称为钟差x(t),采用二次多项式对时钟钟差进行建模的具体步骤为:<maths num="0001" id="cmaths0001"><math><![CDATA[<mrow><mi>x</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>=</mo><msub><mi>a</mi><mn>0</mn></msub><mo>+</mo><msub><mi>a</mi><mn>1</mn></msub><mrow><mo>(</mo><mi>t</mi><mo>-</mo><msub><mi>t</mi><mn>0</mn></msub><mo>)</mo></mrow><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><msub><mi>a</mi><mn>2</mn></msub><msup><mrow><mo>(</mo><mi>t</mi><mo>-</mo><msub><mi>t</mi><mn>0</mn></msub><mo>)</mo></mrow><mn>2</mn></msup><mo>+</mo><msub><mi>&epsiv;</mi><mi>x</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000505734720000011.GIF" wi="737" he="111" /></maths>式中前三项表示时钟的时间系统误差参数,a<sub>0</sub>为时钟初始相位偏差,a<sub>1</sub>为原子钟的初始频率偏差,a<sub>2</sub>为原子钟的线性频率漂移率,ε<sub>x</sub>(t)为时钟受噪声影响引起的时钟偏差的随机变化分量,t<sub>0</sub>为原子钟的起始时刻,t原子钟的读数;其中上述步骤二中,考虑频率漂移变化的影响,星载原子钟模型采用3个状态分量,建立滤波的状态方程具体步骤为:X<sub>k</sub>=Φ<sub>k‑1,k</sub>X<sub>k‑1</sub>+W<sub>k‑1</sub>其中:<img file="FDA0000505734720000021.GIF" wi="871" he="320" />a<sub>0,k</sub>、a<sub>1,k</sub>、a<sub>2,k</sub>分别为t<sub>k</sub>历元的卫星时钟偏差、频率误差和频率漂移率;τ为滤波周期;<img file="FDA0000505734720000022.GIF" wi="152" he="66" />w<sub>f,k‑1</sub>、w<sub>a,k‑1</sub>分别为t<sub>k‑1</sub>历元的时钟相位噪声、频率噪声和频率漂移率噪声,具体统计特性由选用的时钟的Allan方差确定,满足如下关系:<maths num="0002" id="cmaths0002"><math><![CDATA[<mrow><mi>E</mi><mo>{</mo><mi>W</mi><mrow><mo>(</mo><mi>k</mi><mo>)</mo></mrow><msup><mi>W</mi><mi>T</mi></msup><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow><mo>}</mo><mo>=</mo><msub><mi>Q</mi><mi>k</mi></msub><msub><mi>&delta;</mi><mi>kj</mi></msub><mo>=</mo><munderover><mo>&Integral;</mo><mn>0</mn><mi>&tau;</mi></munderover><mi>&Phi;</mi><mrow><mo>(</mo><mi>&tau;</mi><mo>,</mo><mi>T</mi><mo>)</mo></mrow><mi>Q</mi><msup><mi>&Phi;</mi><mi>T</mi></msup><mrow><mo>(</mo><mi>&tau;</mi><mo>,</mo><mi>T</mi><mo>)</mo></mrow><mi>dT</mi><mo>=</mo><mfenced open='[' close=']'><mtable><mtr><mtd><msub><mi>q</mi><mn>1</mn></msub><mi>T</mi><mo>+</mo><msub><mi>q</mi><mn>2</mn></msub><mfrac><msup><mi>T</mi><mn>3</mn></msup><mn>3</mn></mfrac><mo>+</mo><msub><mi>q</mi><mn>3</mn></msub><mfrac><msup><mi>T</mi><mn>5</mn></msup><mn>20</mn></mfrac></mtd><mtd><msub><mi>q</mi><mn>2</mn></msub><mfrac><msup><mi>T</mi><mn>2</mn></msup><mn>2</mn></mfrac><mo>+</mo><msub><mi>q</mi><mn>3</mn></msub><mfrac><msup><mi>T</mi><mn>4</mn></msup><mn>3</mn></mfrac></mtd><mtd><msub><mi>q</mi><mn>3</mn></msub><mfrac><msup><mi>T</mi><mn>3</mn></msup><mn>6</mn></mfrac></mtd></mtr><mtr><mtd><msub><mi>q</mi><mn>2</mn></msub><mfrac><msup><mi>T</mi><mn>2</mn></msup><mn>2</mn></mfrac><mo>+</mo><msub><mi>q</mi><mn>3</mn></msub><mfrac><msup><mi>T</mi><mn>4</mn></msup><mn>3</mn></mfrac></mtd><mtd><msub><mi>q</mi><mn>2</mn></msub><mi>T</mi><mo>+</mo><msub><mi>q</mi><mn>3</mn></msub><mfrac><msup><mi>T</mi><mn>3</mn></msup><mn>3</mn></mfrac></mtd><mtd><msub><mi>q</mi><mn>3</mn></msub><mfrac><msup><mi>T</mi><mn>2</mn></msup><mn>2</mn></mfrac></mtd></mtr><mtr><mtd><msub><mi>q</mi><mn>3</mn></msub><mfrac><msup><mi>T</mi><mn>3</mn></msup><mn>6</mn></mfrac></mtd><mtd><msub><mi>q</mi><mn>3</mn></msub><mfrac><msup><mi>T</mi><mn>2</mn></msup><mn>2</mn></mfrac></mtd><mtd><msub><mi>q</mi><mn>3</mn></msub><mi>T</mi></mtd></mtr></mtable></mfenced></mrow>]]></math><img file="FDA0000505734720000023.GIF" wi="1552" he="353" /></maths>其中,Q(t)δ(τ)=E{w(t+τ)w(t)}Q(t)=diag[q<sub>1</sub>q<sub>2</sub>q<sub>3</sub>]是系统过程噪声向量,w(t)的方差强度阵,q<sub>1</sub>、q<sub>2</sub>、q<sub>3</sub>为原子钟的功率谱密度系数;T为滤波间隔,Q<sub>k</sub>为状态噪声向量的协方差矩阵,Φ(τ,T)为状态转移矩阵;其中上述步骤三中,针对星座中建立星间链路的两颗卫星,采用星间双向测距方法进行伪距测量,建立观测方程具体步骤为:某一参考历元时刻t<sub>k</sub>的星间测量方程为:ρ<sub>ji</sub>=D+c·δt<sub>i</sub>‑c·δt<sub>j</sub>+n<sub>ji</sub>  (1.7)ρ<sub>ij</sub>=D+c·δt<sub>j</sub>‑c·δt<sub>i</sub>+n<sub>ij</sub>  (1.8)式中:ρ<sub>ji</sub>表示卫星S<sub>j</sub>到S<sub>i</sub>的改化伪距,n<sub>ji</sub>表示测量噪声,δt<sub>j</sub>与δt<sub>i</sub>分别表示卫星S<sub>j</sub>、S<sub>i</sub>的绝对钟差,D为卫星间的真实距离,c为光速;由(1.7)和(1.8)式,进行伪距差分可得:z<sub>k</sub>=(ρ<sub>ij</sub>‑ρ<sub>ji</sub>)/2c=δt<sub>j</sub>‑δt<sub>i</sub>+n<sub>k</sub>/2c    (1.9)其中,n<sub>k</sub>=n<sub>ij</sub>‑n<sub>ji</sub>;其中上述步骤四中,根据系统的状态方程和测量方程,进行Kalman滤波具体步骤为,采用下列方程进行Kalman滤波解算:<maths num="0003" id="cmaths0003"><math><![CDATA[<mfenced open='{' close=''><mtable><mtr><mtd><msub><mi>X</mi><mrow><mi>k</mi><mo>/</mo><mi>k</mi><mo>-</mo><mn>1</mn></mrow></msub><mo>=</mo><msub><mi>&phi;</mi><mrow><mi>k</mi><mo>-</mo><mn>1</mn><mo>,</mo><mi>k</mi></mrow></msub><msub><mi>X</mi><mrow><mi>k</mi><mo>-</mo><mn>1</mn></mrow></msub></mtd></mtr><mtr><mtd><msub><mi>P</mi><mrow><mi>k</mi><mo>/</mo><mi>k</mi><mo>-</mo><mn>1</mn></mrow></msub><mo>=</mo><msub><mi>&phi;</mi><mrow><mi>k</mi><mo>-</mo><mn>1</mn><mo>,</mo><mi>k</mi></mrow></msub><msub><mi>P</mi><mrow><mi>k</mi><mo>-</mo><mn>1</mn></mrow></msub><msubsup><mi>&phi;</mi><mrow><mi>k</mi><mo>-</mo><mn>1</mn><mo>,</mo><mi>k</mi></mrow><mi>T</mi></msubsup><mo>+</mo><msub><mi>Q</mi><mrow><mi>k</mi><mo>-</mo><mn>1</mn></mrow></msub></mtd></mtr><mtr><mtd><msub><mi>K</mi><mi>k</mi></msub><mo>=</mo><msub><mi>P</mi><mrow><mi>k</mi><mo>/</mo><mi>k</mi><mo>-</mo><mn>1</mn></mrow></msub><msubsup><mi>H</mi><mi>k</mi><mi>T</mi></msubsup><msup><mrow><mo>(</mo><msub><mi>H</mi><mi>k</mi></msub><msub><mi>P</mi><mrow><mi>k</mi><mo>/</mo><mi>k</mi><mo>-</mo><mn>1</mn></mrow></msub><msubsup><mi>H</mi><mi>k</mi><mi>T</mi></msubsup><mo>+</mo><msub><mi>R</mi><mi>k</mi></msub><mo>)</mo></mrow><mrow><mo>-</mo><mn>1</mn></mrow></msup></mtd></mtr><mtr><mtd><msub><mi>X</mi><mi>k</mi></msub><mo>=</mo><msub><mi>X</mi><mrow><mi>k</mi><mo>/</mo><mi>k</mi><mo>-</mo><mn>1</mn></mrow></msub><mo>+</mo><msub><mi>K</mi><mi>k</mi></msub><mrow><mo>(</mo><msub><mi>Z</mi><mi>k</mi></msub><mo>-</mo><msub><mi>H</mi><mi>k</mi></msub><msub><mi>X</mi><mrow><mi>k</mi><mo>/</mo><mi>k</mi><mo>-</mo><mn>1</mn></mrow></msub><mo>)</mo></mrow></mtd></mtr><mtr><mtd><msub><mi>P</mi><mi>k</mi></msub><mo>=</mo><mrow><mo>(</mo><mi>I</mi><mo>-</mo><msub><mi>K</mi><mi>k</mi></msub><msub><mi>H</mi><mi>k</mi></msub><mo>)</mo></mrow><msub><mi>P</mi><mrow><mi>k</mi><mo>/</mo><mi>k</mi><mo>-</mo><mn>1</mn></mrow></msub><msup><mrow><mo>(</mo><mi>I</mi><mo>-</mo><msub><mi>K</mi><mi>k</mi></msub><msub><mi>H</mi><mi>k</mi></msub><mo>)</mo></mrow><mi>T</mi></msup><mo>+</mo><msub><mi>K</mi><mi>k</mi></msub><msub><mi>R</mi><mi>k</mi></msub><msubsup><mi>K</mi><mi>k</mi><mi>T</mi></msubsup></mtd></mtr></mtable></mfenced>]]></math><img file="FDA0000505734720000031.GIF" wi="821" he="382" /></maths>式中,Q为系统观测噪声阵,R为系统噪声阵,其大小与测距精度相关,H为观测矩阵,Z为观测向量,K为增益矩阵,Φ为状态转移矩阵,状态变量取为各节点与基准节点的时钟参数之差,其中取卫星1为基准。
地址 100191 北京市海淀区学院路37号