基于卫星双向时间频率传递进行钟差预报的方法研究.doc

第 32 卷 第 1 期2007 年 1 月武 汉 大 学 学 报 信 息 科 学 版Geo matics a nd Info r matio n Science of Wuha n U niver sit yVol . 32 No . 1J a n. 2007文章编号 :167128860 (2007) 0120043204文献标志码 : A基于卫星双向时间频率传递进行钟差预报的方法研究郭海荣1生2杨元喜3何海波1 ,4杨( 1 信息工程大学测绘学院 ,郑州市陇海中路 66 号 ,450052)( 2 北京环球信息应用开发中心 ,北京市海淀区北清路 22 号 ,100094) ( 3 西安测绘研究所 ,西安市雁塔路中段 1 号 ,710054)( 4 武汉大学测绘学院 ,武汉市珞喻路 129 号 ,430079)摘 要 :用多项式拟合 、谱分析 、改进的 A R 模型三种方法对由卫星双向时间频率传递得出的钟差时间序列进行了拟合和预报分析 。为了抵制钟差时间序列中异常值的影响 ,引入了“抗差等价权”,利用卫星双向时间频 率传递得到的 1 d 的钟差 ,按不同采样率 、不同时间跨度进行计算分析 。结果表明 ,抗差估计的预报精度明显 高于最小二乘估计 ;平滑值的预报精度高于采样值 ;由于钟差时间序列中有明显的周期变化 ,多项式进行钟差预报的精度不可靠 ;用谱分析进行钟差预报的精度不高 ,但可以发现钟差时间序列中的主要周期变化 ;改进的A R 模型预报精度最高 ,预报 6 h 钟差的 RM S 在 1 ns 左右 。 关键词 :时间预报 ;时间比对 ;卫星双向时间频率传递 中图法分类号 : P228 . 41随着扩频技术和伪随机码技术的应用 ,时间比对技术得到了迅速提高 ,现在利用卫星双向时 间频率传递 ( t wo2way sat ellit e ti me a nd f reque n2cy t ra n sf er , T WS T F T) 进行时间比对 ,精度可达0 . 20 . 3 n s 123 。由 T WS T F T 得到的钟差时间 序列为随机时间序列 。在预报过程中 ,随机误差和粗差对模型的影响将会变成系统性影响 。对时 间序列中的随机误差 ,用多余观测减弱其影响 ;对 时间序列中的粗差 ,用抗差估计来抵制其对参数估值的影响 。本文采用多项式拟合 、频谱分析 、改进的 A R 模型三种方法对钟差时间序列进行拟合 和预报分析 。知参数向量 ; et 为模型误差 。式 ( 1) 表示成矩阵形式为 :X = Ha + e( 2)其中 ,X 为观测向量 ; H 为 n m 的设计矩阵 ; a为未知参数向量 ; e 为误差向量 。1 . 2谱分析模型pa0 + b0 ti + A k si n (2f k t i + k ) + ei( 3)x i =k = 1式中 , a0 、b0 分别为长期变化的常数项 、系数项 ; p为主要周期函数的个数 ; f k 为对应周期项的频率 ; A k 、k 分别为对应周期项的振幅和相位 ; ei 为 x i 的残差向量 ; p 、f k 可由能量2频率图确定 4 。式 ( 3 )线性化后的矩阵形式为 :函数模型及其估值1( 4)X = Ha + e式中 , X 为由 x i 组成的向量 ; H 为设计矩阵 ; a 为未知参数向量 ; e 为 x i 的残差向量 。1 . 1多项式模型假设相 对 于 时 间 t1 , t2 , x n , 则有 :, tn 的 钟 差 为 x 1 ,1 . 3改进的 AR 模型钟差时间序列用 A R 模型表示为 :x2 ,x t= a0 + a1t + a2t2+ amtm + et( 1)+t5( )x t = a1 x t - 1 + a2 x t - 2 + am x t - m + e式中 ,t = t -t0 ;x t = x t - x t ; a0 、a1 、 、am 为未0式中各符号的意义同式 ( 1) 。式 ( 5) 的矩阵形式同收稿日期 :2006210227 。式 ( 4) 。考虑到钟差时间序列中明显存在长期趋势的 变化规律 , 此外 , 还有似半天 、一天的周期变化 , 为了能采用 A R 序列模型来建模和预报 , 必须将钟差时间序列转化为平稳序列 , 为此 , 对钟差时间序 列作下述差分运算 5 ,依据 经 验 , A R 模 型 的 阶 次 应 当 小 于 n 或n/ 10 , n 为观测量的个数 7 。计算中 , A R 模型的 阶次取为 20 。3预报模型及其精度A AX =( X t - X t - )2 ( X t - 1X t - - 1 )-+( 6)多项式预报求得未知参数后 , 预报值 x t + l 为 :3 . 1( X t - 2X t - - 2 )-式中 , 为 时 间 序 列 中 最 长 的 周 期 项 。显 然 ,A AX 序列已消除了线性和二次项 ,及其倍频 周期项也基本消除而变换成一个平稳的序列 , 这样就可以根据 A R 模型的建模理论进行建模 。1 . 4参数估计式 ( 2) 、式 ( 4) 中未知参数向量的 L S 估计为 :+ HT( 12)x t + l=x tt + l a0式中 , HT l = tt2tm ,t= t- t。t +t + lt + lt + lt + lt + l0预报的残差均方根为 :l2RM S =( X t + i - X t + i )( 13)/ li = 13 . 2谱分析预报a = ( HT P H) - 1 HT PX( 7)d X t + l = HT a( 14)t + lTTT= V PV式 中 , a 1 tt + l= a0 b0 a1 b1si n ( 2f 1 tt + l )a p b p ;Ht + l =2( 8)0n - mco s ( 2f 1 tt + l )式中 , P 为观测向量的权矩阵 ;V 为残差向量 。为了抵制钟差观测序列中粗差的影响 , 引入 抗差估计 , 未知参数的抗 差解为 : si n (2f p t t + l co s (2f p t t + l ) T 。求出参数估值后 , 其预报值的残差均方根可 按式 ( 13) 计算 。a R = ( HT P H) - 1 HTPX( 9)3 . 3改进的 AR 模型预报TV PVX t + l - + 2 ( X t + t - 1 - X t + l - - 1 ) -mX t +1 =2( )10=0n - 3 - s- X t + l - - 2 ) + aj A A X t + l ( 15)式中 , P为等价权矩阵 6 ; s 为 p i = 0 的个数 。残差均方根为 :n( X t + l - 2j = 1求出 参 数 值 后 , 其 预 报 值 的 残 差 均 方 根 可 按式 ( 13) 计算 。2RM S =V i / n( 11)i = 14计算与分析模型阶次的确定24 . 1数据预处理为熟悉时间序列的误差特性和观测误差的近GPS 钟差预报都是采用一阶 、二阶多项式 ,考虑到计算中采用的原子钟的性能不如 GPS 系统的 原子钟 ,下面对钟差时间序列的一阶 、二阶 、三阶多 项式的拟合精度及其预报精度进行比较分析 。由能量谱分析可知 ,钟差时间序列除有明显的长期变化外 ,还有周期性变化规律 。15 s 采样(平滑) 的钟差时间序列有 6 . 8 h 、3 . 8 h 、2 . 84 h 、2 . 27 h的周期变化 ; 30 s 和 1 mi n 采样 ( 平滑) 的 钟差时间序列有 6 . 8 h 、4 . 27 h 、2 . 84 h 、2 . 27 h 的 周期变化 。由此可见 ,不同采样率的钟差时间序 列的主要周期项稍微有些差异 。高频采样率能发现一些高频的周期变化 ,但高频周期变化太多也 会干扰长周期变化的提取 ; 低频采样会掩盖一些 高频变化 ,有时有利于长周期变化的提取 。下面 将采用 15 s 、30 s 、1 mi n 采样和平滑的钟差时间 序列进行分析 。似分布 ,需要对钟差时间序列进行初步分析 。计算时 ,采用由 TWS T F T 方法得到的 1 d 的钟差时 间序列 ,断点用线性内差进行插值 。分析前 ,还要 进行由卫星运动引起的 Sa gnac 效应改正 、时钟的 零值改正 8 。图 1 到图 4 是钟差 、钟速 、钟加速度以及 Sa gnac 效应改正 1 d 内随时间的变化情况 。由图 1图 4 可知 , 钟差时间序列除有明显 的长期变化外 ,还有明显的周期变化 。另外 ,钟差时间序列还有明显的粗差点存在 。由卫星运动引 起的 Sa gnac 效应改正虽有周期性变化规律 ,但周 期变化的振幅很小 ,在 0 . 1 n s 量级 ,因此 , Sa gnac 效应对钟差周期性变化的贡献不大 ,鉴于此 ,应从 其他方面查找原因 。4 . 2 计算结果 、分析下面采用改正后的 0 . 5 h 、1 h 、6 h 钟差时间第 32 卷第 1 期郭海荣等 :基于卫星双向时间频率传递进行钟差预报的方法研究451 s 5 s 15 s 30 s 1 mi n) 的钟差时间序列进行分析 。然后用求得的模型系数预报钟差 ,并与序列拟合一阶 、二阶 、三阶多项式系数 、由能量谱分析确定的函数模型系数以及改进的 A R 模型系 数 ,并进行精度分析 。为了分析不同采样率的钟 差时间序列对预报精度的影响 ,计算时采用 0 . 2s 、1 s 、5 s 、15 s 、30 s 、1 mi n 采样和平滑 ( 0 . 2 s 相应的观测值进行比较 ,求出预报精度的 RM S 。在估计模型系数时 , 采用 L S 估计和抗差估计两种统计方法 。图 1 钟差时间序列Fig. 1 Clock Diff ere2nce Time Serie s图 2 钟速 (频差)时间序列Fig. 2 Clock Velocit yTime Series图 3 钟加速度 (频漂) 时间序列Fig. 3 Clock Acceleratio nTime Serie s图 4 Sagnac 效应改正Fig. 4 Sagnac Eff ectCo r rectio n跨度要足够长) ,因此结果比较理想 。又由于钟差时间序列中有明显的粗差存在 ,因此抗差估计的 预报精度明显高于 L S 估计 。下面仅给出改进的A R 模型抗差估计的预报精度 。表 1 改进的 A R 模型在不同采样率 、不同时间 跨度的预报精度/ nsTa b. 1 Predictio n Preci sio n of Diff erent Sa mpling Rate ,Diff erent Sp an U sing Mo dified A R Mo del由图 5 、图 6 可知 , 经过 Sa gnac 效应改正和零值改正后的钟差时间序列在去除一阶或二阶长 期变化 项 后 , 周期 变化 的 幅度 仍在 0 70 n s 之间 。由于周期变化的存在 ,用不同时间跨度的钟 差进行多项式预报 ,精度可靠性较差 ,计算结果也 证实了上述结论 。虽然 GPS 星历钟差预报采用 一阶 、二阶多项式 ,但这种方法对本文分析的原子 钟不适用 ,因为本文采用的原子钟的稳定性远远不如 GPS 。欲采用形式比较简单 的 多项 式进 行 预报 ,需要去除时间序列中周期变化的影响 。能量谱分析可以提取钟差时间序列中的主要 周期变化 ,但实际的钟差时间序列除主要周期变 化外 ,还有许多高频周期变化和对时间序列贡献不大的周期变化 。由于谱分析的函数模型不可能 包含所有的周期变化 ,因此 ,谱模型的预报精度不 会很高 。用 6 h 钟差 ( 15 s 、30 s 、1 mi n 采样值和 平滑值) 预报 6 h 钟差的 RM S 达到 60 n s 。由于改进的 A R 模型通过差分运算已去除了一阶 、二阶长期变化 ,周期变化也基本消除 (时间采样率1 h 预报 1 h RMS3 h 预报 3 h RMS6 h 预报 6 h RMS0 . 2 s-1 s 采样1 s 平滑5 s 采样5 s 平滑15 s 采样15 s 平滑30 s 采样30 s 平滑-0 . 46-1 . 00 . 651 . 00 . 31 . 60 . 31 . 90 . 6-0 . 34-1 . 00 . 81 . 240 . 91 . 40 . 71 . 430 . 76-0 . 45-0 . 940 . 71 . 10 . 651 . 20 . 51 . 340 . 56-由表 1 可知 ,改进的 A R 模型在不同采样率 、不同时间跨度的预报精度都在 1 n s 左右 ; 平滑值的预报精度明显高于采样值的预报精度 ,但平滑 值要对原始数据进行存储 ,而采样值只需对采样 点的数据进行存储 ,大大降低了数据的存储量 。 为消除钟差时间中似半天的周期变化 ,改进的 A R 模型实际上只用半天数据参与模型参数的拟合和预报 。要得到精度更高的预报精度 ,时间 序列的长度至少是周期最长的主要周期变化的两 倍 ,以便从钟差时间序列中消除比半天周期还长 的主要周期变化 。图 5 去除一阶长期变化项后的钟差Fig. 5 Clock Diff erenceWiped Off O ne2O r der Lo ng Ter m图 6 去除二阶长期变化项后的钟差Fig. 6 Clock Diff erenceWip ed Off Two2O r der Lo ng Ter mCalibratio n System fo r Two2Way Time and Fre2 quency Tra nsf er C . IE E E Intl . Freq. Co nt . Symp . , Pi scata Way , New J er sey , 1998郭海荣 , 杨元喜 , 焦文海 . 地心运动时间序列的抗差谱分析 J . 测绘学报 ,2003 ,32 (4) : 3082312丁月蓉 , 郑大伟 . 天文 测 量 数 据 的 处 理 方 法 M .南京 : 南京大学出版社 , 1990Ya ng Yuanxi . Ro bust Estimatio n of Geo detic Da2 t um Tra nsfo r matio n J . J o ur nal of Geo desy , 1999 (73) : 2682274Ya ng Yuanxi . So me N umerical Predictio n Met ho ds fo r t he Wind Sp eed in t he Sea Ba sed o n ERS21 Scat2 t ero meter Wind Data J . Survey Review , 2001 (36) : 1212131刘利 . 相对论时间比对理论与高精度时间同步技术 D . 郑州 :信息工程大学 , 2004结语5 4 本文仅用 1 d 的钟差时间序列进行了分析 ,得出了一些初步结论 ,重要的是给出了本文采用 的原子钟的钟差预报方法 。下一步希望能对更长 时间的钟差时间序列进行分析 ,以便给出可靠的预报结果 。另外 , A R 模型阶次的确定还有待进 一步研究 。 5 6 7 参考文献 1 A scar r unz F G , J eff ert s S R , Pa r ker T E. Enviro n2mental Eff ect s o n Er ro r s in Two2way Time a nd Fre2 8 quency Transf er C . IE E E CP EM 1996Dig. , Braunschweig , Ger ma ny , 1996A scar r unz F G , J eff ert s S R , Pa r ker T E.Co nf . 2 Ea rt hSt atio n Er ro r s in Two2Way Time a nd FrequencyTra nsf er J . IE E E Tra ns Inst r um Mea s , 1997 , 46 :2052208A scar r unz F G , J eff ert s S R , Pa r ker T E. A Delay第一作者简介 :郭海荣 ,博士生 。现主要从事原子钟时频特性和卫星导航等方面的研究工作 。E2mail : hai ro ngguo 263 . net 3 Numerical Predict ion Met hods f or Cl ock Diff erence Ba sed onTwo2Way Satell ite Time and Frequency Transf er DataGU O H ai ron g 1Y A N G S he n g 2Y A N G Y u a n x i 3H E H ai bo1 , 4( 1 In stit ut e of Sur veyi ng and Mappi ng , Info r matio n Engi neeri ng U ni ver sit y , 66 Mi ddle Lo nghai Road , Zhengzho u 450052 , Chi na)( 2 Beiji ng Huanqi u Co mpa ny of Info r matio n Applicatio n a nd Develop ment , 22 Beiqi ng Road , Haidia n Di st rict , Beiji ng 100094 , Chi na)( 3 Xia n Re sea rch In stit ut e of Sur veyi ng a nd Mappi ng , 1 Middle Ya nt a Road , Xia n 710054 , Chi na)( 4 School of Geo de sy a nd Geo matics , Wu han U ni ver sit y , 126 L uo yu Road , Wu ha n 430079 ,Chi na)Abstract : Three f unctio nal mo del s , pol yno mial , sp ect ral a nal ysi s , a nd mo difie d A R mo del ,a re fit t ed a nd co mp a re d i n fit ti ng a nd p re dicti ng t he clock diff ere nce s ba se d o n t he dat a serie s derived f ro m t wo2way sat ellit e ti me a nd f reque ncy t ra n sf e r . A ro bu st equivale nt wei ght i s ap2 p lie d , w hic h co nt rol s t he significa nt of o ut l yi ng o b se r vatio n s. So me co ncl u sio n s sho w t hatt he p re