李西順, 吳元偉, 姚 當(dāng), 劉 佳, 馬浪明, 楊旭海, 張首剛

(1. 中國科學(xué)院國家授時中心, 陜西 西安 710600; 2. 中國科學(xué)院大學(xué)電子電氣與通信工程學(xué)院,北京 100049; 3. 中國科學(xué)院精密導(dǎo)航定位與定時技術(shù)重點實驗室, 陜西 西安 710600;4. 中國科學(xué)院時間頻率基準(zhǔn)重點實驗室, 陜西 西安 710600)

0 引 言

世界時(universal time, UT1)是以地球自轉(zhuǎn)定義的時間尺度,是構(gòu)成標(biāo)準(zhǔn)時間,即協(xié)調(diào)UT1(UT1 coordinated, UTC)的重要組成部分。UT1與極移、章動改正項組成了地球定向參數(shù)(earth orientation parameter, EOP),關(guān)于需要在地球坐標(biāo)和天球坐標(biāo)之間建立坐標(biāo)變換的應(yīng)用,UT1都必不可少,UT1數(shù)據(jù)服務(wù)是衛(wèi)星導(dǎo)航、深空探測等領(lǐng)域所必須的基礎(chǔ)服務(wù)。

當(dāng)前權(quán)威的UT1產(chǎn)品由國際地球自轉(zhuǎn)服務(wù)組織(International Earth Rotation Service, IERS)提供。該產(chǎn)品由多種空間大地測量技術(shù)綜合得到。在這些技術(shù)中甚長基線干涉測量(very long baseline interferometry, VLBI)技術(shù)和全球衛(wèi)星導(dǎo)航定位系統(tǒng)(global navigation satellite system, GNSS)技術(shù)發(fā)揮著重要作用。其中VLBI技術(shù)不用考慮動力學(xué)因素效應(yīng),是僅有的能夠做到監(jiān)測UT1變化的空間大地測量技術(shù),也是目前唯一能夠解算全部地球定向參數(shù)的空間技術(shù)。VLBI技術(shù)因具有極高的測角精度和長期的穩(wěn)定性,在解算EOP 和維持地球及天球參考框架方面發(fā)揮著舉足輕重的作用;但由于成本原因,尚未進行長期的連續(xù)觀測,導(dǎo)致其觀測頻度不足。GNSS 技術(shù)因觀測站遍布全球,可實現(xiàn)不間斷的連續(xù)觀測,通過解算可得到日長變化ΔLOD、極移,但GNSS解算的ΔLOD受制于軌道升交點進動問題,長期穩(wěn)定性有待提高。將上述兩種空間技術(shù)融合解算可以有效解決VLBI 技術(shù)觀測UT1頻率不足以及GNSS技術(shù)解算ΔLOD長期穩(wěn)定性上的不足,從而得到更連續(xù)、精度更高的UT1產(chǎn)品。

本文使用仿真與實測相結(jié)合的方法,首先通過馬爾可夫鏈蒙特卡羅(Markov chain Monte Carlo, MCMC)模擬仿真了UT1隨機誤差、時間間隔、系統(tǒng)偏差對于融合精度的影響以及權(quán)比的依賴關(guān)系。然后將國內(nèi)首套寬帶VLBI2010系統(tǒng)解算的UT1序列與國際GNSS監(jiān)測評估系統(tǒng)(international GNSS monitoring and assessment system, iGMAS)解算的ΔLOD序列,運用Vondrak方法融合在一起,并通過MCMC算法估計最優(yōu)的權(quán)重比。最后將該融合結(jié)果與IERS 14 C04 UT1序列進行比對,用以評估該融合方法的效能。

1 融合算法

1.1 理論支撐

UT1是由地球自轉(zhuǎn)角定義的時間尺度,是與地球固連的地球坐標(biāo)系(轉(zhuǎn)動的地球坐標(biāo)系,非慣性系)和天球坐標(biāo)系(無旋轉(zhuǎn)慣性系)間的重要的坐標(biāo)變換參數(shù)。ΔLOD是描述地球自轉(zhuǎn)運動的主要參量之一,代表了天文意義上的1 d與標(biāo)準(zhǔn)日長86 400 s之間的差別,體現(xiàn)了地球自轉(zhuǎn)角速度的快慢。UT1和ΔLOD間存在的關(guān)系如下:

d(UT1-TAI)/d=ΔLOD

(1)

式中:TAI為國際原子時。

對于連續(xù)的ΔLOD序列,積分后就是UT1 相對某一初值的變化量,利用GNSS衛(wèi)星數(shù)據(jù)估計ΔLOD已經(jīng)成為例行的工作。國外,美國海軍天文臺、噴氣推進實驗室、布拉格天文研究所已進行了相關(guān)融合算法的研究。

1.2 Vondrak算法

本文采用Vondrak平滑組合算法對中國科學(xué)院國家授時中心(National Time Service Center, NTSC)VLBI2010系統(tǒng)解算的UT1數(shù)據(jù)與iGMAS解算的ΔLOD數(shù)據(jù)進行加權(quán)融合。該組合平滑算法是由Vondrak和 Cepek在2000年針對 UT1與ΔLOD融合問題而提出的。該算法不但能夠?qū)Φ乳g距和等精度的測量數(shù)據(jù)進行平滑處理,而且能夠?qū)Σ坏乳g距和不等精度的數(shù)據(jù)進行處理,在天文學(xué)上使用廣泛,比如國際地球自轉(zhuǎn)服務(wù)、原子尺度平滑、恒星天文學(xué)分析。

關(guān)于一組測量數(shù)據(jù)(,) (=1,2,…,),其中是自變量,是對應(yīng)自變量的測量值,Vondrak算法的理論思想如下:

=+=min

(2)

(3)

(4)

Vondrak算法使用多項式形式來表示平滑函數(shù),多項式的一般形式如下:

(5)

對式(5)中求三次導(dǎo)并代入式(4)得到

(6)

當(dāng)式(2)取得最小值時，有

(7)

將式(3)與式(6)代入式(7),便可獲得Vondrak平滑方法的基礎(chǔ)方程組:

(8)

其中系數(shù)矩陣,和中的系數(shù)為

在原有的Vondrak方法上,增加一階導(dǎo)數(shù)擬合度的定義,即Vondrak組合平滑算法:

(9)

(10)

運用上述Vondrak求解的過程,即可求得具體的方程解。

1.3 MCMC方法及權(quán)重選取

由上述理論可得,觀測值越接近于平滑值,則權(quán)重應(yīng)該越大;產(chǎn)品的測量精度越高,則權(quán)重越大。實際使用中,權(quán)重直接影響最終融合產(chǎn)品的平滑效果和精度水平。由于實測UT1序列存在精度和時間間隔的不確定性,很難通過直接的理論推導(dǎo)方式得到最優(yōu)的權(quán)比,因此本文通過MCMC的方法結(jié)合與事后精密EOP數(shù)據(jù)的比對,嘗試得到對UT1觀測值及一階導(dǎo)數(shù)值ΔLOD的最優(yōu)權(quán)重比。

MCMC方法是處理多維、復(fù)雜問題的常用方法,在空間物理學(xué)、圖像分析、金融等領(lǐng)域發(fā)揮了重大的作用。該方法的基本思路是基于某種采樣的方法進行大樣本隨機采樣,來構(gòu)造一個平穩(wěn)分布的馬爾可夫鏈,再對這些樣本進行蒙特卡羅模擬,獲得有效的統(tǒng)計量。

本文針對UT1及ΔLOD的權(quán)重優(yōu)化問題,對UT1和ΔLOD的權(quán)重做大樣本隨機采樣,采樣方法為吉布斯采樣,通過將二維空間采樣變成一維空間,顯著提高了MCMC方法的收斂速度。

2 輸入數(shù)據(jù)源

2.1 UT1數(shù)據(jù)

采用的UT1數(shù)據(jù)是由NTSC的國內(nèi)寬帶VLBI2010系統(tǒng)自主測定的UT1序列。

為滿足國家UT1測量需求,NTSC和上海天文臺共同研制了我國首套寬帶VLBI2010系統(tǒng),由三亞觀測站、吉林觀測站、喀什觀測站,及西安數(shù)據(jù)處理中心組成。自2018年以來,NTSC利用VLBI2010系統(tǒng)已開展了上百次常規(guī)的UT1觀測試驗?；竟ぷ髁鞒淌羌?、三亞、喀什三站同時根據(jù)觀測綱要開始觀測,采集并記錄觀測數(shù)據(jù),數(shù)據(jù)上傳至西安數(shù)據(jù)處理中心;觀測數(shù)據(jù)進行相關(guān)處理及后相關(guān)處理得到經(jīng)過時延校準(zhǔn)的射電源時延序列;進而結(jié)合三站氣象站數(shù)據(jù),最終通過數(shù)據(jù)分析軟件對時延結(jié)果運用最小二乘算法進行解算,確定UT1。

本文使用了2018年和2020年由NTSC通過VLBI2010系統(tǒng)自主測量并解算的UT1序列。2018年選取了從2018年7月24號到2018年9月20號的時間段,共53次觀測數(shù)據(jù)。2020年選取了從2020年8月14號到2020年11月4號的時間段,共68次觀測數(shù)據(jù)。

2.2 ΔLOD數(shù)據(jù)

采用的ΔLOD數(shù)據(jù)有兩種來源,分別是來自于iGMAS解算的ΔLOD產(chǎn)品,此產(chǎn)品用于實測驗證。以及國際GNSS服務(wù)(International GNSS Service, IGS)提供的ΔLOD產(chǎn)品,此產(chǎn)品用于仿真驗證。

2.3 事后精密EOP 數(shù)據(jù)

事后精密C04序列,從巴黎天文臺服務(wù)器下載,其中eopc04_IAU2000.62-now和eopc04.62-now文件格式相同,主要的差異是天極軸改正項歲差-章動與IAU模型的偏差修正項,兩個文件對應(yīng)的歲差-章動模型分別為IAU2000A歲差-章動模型和IAU1980歲差-章動模型。由于C04序列是事后精密產(chǎn)品,故而時效性不高,主要用于精度比對和研究。

3 仿真分析

為了評估及優(yōu)化融合算法,同時優(yōu)化和指導(dǎo)NTSC 的VLBI2010系統(tǒng)UT1測定方案,根據(jù)文獻[40]中提到的NTSC的VLBI2010系統(tǒng)實測UT1精度為58.8 μs,且不能連續(xù)觀測的現(xiàn)狀,仿真環(huán)節(jié)進行權(quán)重的仿真、時間間隔仿真、隨機誤差仿真及系統(tǒng)偏差的仿真,并研究輸出序列與這些變量間的依賴關(guān)系。

3.1 時間間隔采樣仿真

天文觀測會受到惡劣天氣、設(shè)備維護、儀器故障乃至網(wǎng)絡(luò)中斷等突發(fā)狀況導(dǎo)致觀測的不連續(xù)和中斷。針對這種現(xiàn)象,需研究不同間隔的UT1數(shù)據(jù)對融合精度產(chǎn)生的影響。故而仿真了不同間隔的UT1實驗,用于指導(dǎo)實際的UT1觀測。

仿真分析的流程圖如圖1所示。

圖1 時間間隔仿真Fig.1 Time interval simulation

利用不同時間間隔采樣的精密C04序列替代NTSC解算的UT1序列,并與IGS的ΔLOD序列進行融合。分別進行間隔為1 d、2 d、3 d、4 d、5 d、6 d采樣的仿真實驗。對每種仿真實驗通過MCMC算法,確定最優(yōu)權(quán)重比和最優(yōu)的融合結(jié)果,最后將融合的結(jié)果與C04序列做差,統(tǒng)計其標(biāo)準(zhǔn)差和權(quán)重比值(UT1/ΔLOD)。如表1所示,列代表對應(yīng)的采樣間隔,依次為1 d、2 d、3 d、4 d、5 d、6 d,行代表對應(yīng)的權(quán)比和標(biāo)準(zhǔn)差,其中標(biāo)準(zhǔn)差單位為μs。

表1 不同采樣間隔權(quán)重比值及標(biāo)準(zhǔn)差

可以得出采樣的時間間隔越長,融合的精度就越低。同時,隨著采樣間隔的增加,UT1與ΔLOD的權(quán)重比呈現(xiàn)下降的趨勢,表明權(quán)重與數(shù)據(jù)量有很強的相關(guān)性。這是因為采樣間隔越長,UT1有效的數(shù)據(jù)量就越少,進而UT1對于融合的貢獻也就越少,權(quán)重也會相應(yīng)地降低,所以隨著采樣間隔的增加,融合精度出現(xiàn)下降的趨勢。

3.2 隨機誤差仿真

研究隨機誤差對融合精度的影響,分別在精密C04序列的基礎(chǔ)上加0 μs、10 μs、20 μs、30 μs、40 μs、50 μs、60 μs的高斯白噪聲,再對加噪聲之后的C04序列進行間隔為1 d、2 d、3 d、4 d、5 d、6 d的采樣,共計42種UT1仿真序列,再與ΔLOD序列進行融合,將融合后的UT1與精密C04序列的進行做差,最終統(tǒng)計標(biāo)準(zhǔn)差,42種隨機誤差仿真實驗的結(jié)果如表2所示,行代表對應(yīng)的采樣間隔,列代表對應(yīng)的隨機誤差值,表中數(shù)值即為融合得到的UT1精度值。

表2 誤差仿真

使用折線圖更直觀地展示采樣間隔與隨機誤差對融合UT1產(chǎn)品的影響,如圖2所示,橫坐標(biāo)代表采樣間隔,以d為單位,縱坐標(biāo)為標(biāo)準(zhǔn)差,以μs為單位,圖中折線從下往上依次為0 μs、10 μs、20 μs、30 μs、40 μs、50 μs、60 μs的隨機誤差。

圖2 誤差仿真折線圖Fig.2 Error simulation line chart

隨機誤差的仿真實驗,可以得出,原始數(shù)據(jù)隨機誤差越大,融合精度就越低。隨機誤差與采樣間隔共同影響著融合精度,對于采樣間隔1 d和2 d的UT1數(shù)據(jù)來說,隨機誤差的影響占主導(dǎo)地位。隨機誤差對于精度有很強的相關(guān)性。當(dāng)采樣間隔大于2 d之后,采樣間隔的作用逐漸明顯,此時融合精度與隨機誤差不再是簡單的正相關(guān)。當(dāng)隨機誤差為50～60 μs時,采樣間隔為3～5 d,融合效果較好。

3.3 系統(tǒng)偏差仿真

研究系統(tǒng)偏差對融合效果的影響,分別在精密C04序列的基礎(chǔ)上加0 μs、10 μs、20 μs、30 μs、40 μs、50 μs、60 μs的系統(tǒng)偏差,生成7種UT1仿真序列,用該仿真序列分別替代NTSC解算的UT1序列,再與ΔLOD序列進行融合,將融合后的UT1序列與精密C04序列進行做差,統(tǒng)計標(biāo)準(zhǔn)差和均值,如表3所示。

表3 偏差仿真

可以得出系統(tǒng)偏差對于融合精度沒有影響,改變系統(tǒng)偏差,融合結(jié)果的精度沒有變化,這是因為系統(tǒng)偏差是一個固定值,只會影響整體的平移,不會影響穩(wěn)定性,所以不會改變標(biāo)準(zhǔn)差,對于融合算法沒有明顯的作用。因此,UT1的系統(tǒng)偏差,如站坐標(biāo)、大氣模型不準(zhǔn)導(dǎo)致的系統(tǒng)偏差,只能通過其他方法改正。

4 實測數(shù)據(jù)集融合

4.1 2018年實測數(shù)據(jù)集融合

使用NTSC于2018年下半年通過VLBI2010系統(tǒng)自主觀測并解算的UT1實測數(shù)據(jù)集,選取簡化儒略日58 323～58 381共59天的數(shù)據(jù),采用3sigma原則對誤差較大的數(shù)據(jù)進行剔除,將處理后的數(shù)據(jù)與iGMAS解算的ΔLOD數(shù)據(jù)集進行Vondrak平滑算法融合,并通過MCMC進行權(quán)重的最優(yōu)化,在進行了1 000次的MCMC循環(huán)采樣后,得到UT1與ΔLOD的權(quán)重比值收斂于0.176(UT1/ΔLOD),此時融合之后的UT1序列收斂于一個穩(wěn)定值,將該穩(wěn)定值與精密C04數(shù)據(jù)進行比對分析,工作流程圖如圖3 所示。

圖3 工作原理流程圖Fig.3 Working principle flow chart

圖4(a)展示了3種UT1數(shù)據(jù),分別為精密C04序列的UT1數(shù)據(jù)、NTSC解算的初始UT1數(shù)據(jù)和使用融合算法得到的融合UT1數(shù)據(jù)。其中橫坐標(biāo)為簡化儒略日,縱坐標(biāo)單位為s,黑色、藍色、紅色的標(biāo)識分別代表精密C04的UT1、NTSC初始UT1、融合UT1。圖4(b)縱坐標(biāo)單位為μs,紅色、藍色標(biāo)識分別代表融合UT1與C04_UT1的差值、NTSC初始UT1與C04_UT1的差值。

圖4 2018年UT1數(shù)據(jù)比對分析Fig.4 Comparison and analysis of UT1 data in 2018

使用直方圖對圖4(b)差值進行統(tǒng)計分析,如圖5所示,橫坐標(biāo)為差值,單位為μs,縱坐標(biāo)為歸一化后的值,代表數(shù)量,因融合后的UT1更加密集,故而進行歸一化處理。藍色、紅色分別為NTSC初始UT1與精密C04_UT1的差值,融合UT1與精密C04_UT1的差值。圖6縱坐標(biāo)代表差值,單位為μs。

圖5 2018年融合前后UT1精度比對Fig.5 Comparison of UT1 accuracy before and after combination in 2018

圖6 2018年融合前后產(chǎn)品比對Fig.6 Comparison of products before and after combination in 2018

為了更具體地展示融合精度,統(tǒng)計了差值的均值、標(biāo)準(zhǔn)差、中位數(shù)、最大值及最小值,如表4所示。

表4 2018年精度統(tǒng)計分析

圖5顯示,融合后的數(shù)據(jù)集相對IERS事后精密序列的偏差集中在-50～50 μs之間,相比之下,初始UT1數(shù)據(jù)的偏差分布集中在-80～80 μs之間。圖6直觀地展示了融合后的數(shù)據(jù)要比原始數(shù)據(jù)更為穩(wěn)定。通過表4得出NTSC初始UT1標(biāo)準(zhǔn)差為53.18 μs,融合UT1標(biāo)準(zhǔn)差為31.96 μs,融合后的均值、最值均比原始序列更小,經(jīng)過融合算法,得到更連續(xù)、更穩(wěn)定的UT1序列,其精度提高了37%。

4.2 2020年實測數(shù)據(jù)集融合

根據(jù)2020年NTSC開展的UT1觀測實驗,選取了簡化儒略日59 075～59 157共83 d較為連續(xù)的觀測數(shù)據(jù)。將處理后的2020年NTSC解算的UT1數(shù)據(jù)與iGMAS解算的ΔLOD數(shù)據(jù)集進行Vondrak平滑算法融合,并通過MCMC進行權(quán)重的最優(yōu)化,在進行了1 000次的MCMC循環(huán)采樣,得到UT1與ΔLOD的權(quán)重比值收斂于2.774 6,此時融合之后的UT1序列收斂于一個穩(wěn)定值,將該穩(wěn)定值與精密C04數(shù)據(jù)進行比對分析(方法與2018年相同),結(jié)果如圖7所示。

圖7 2020年UT1數(shù)據(jù)比對分析Fig.7 Comparison and analysis of UT1 data in 2020

均值、標(biāo)準(zhǔn)差、中位數(shù)、最大值及最小值統(tǒng)計結(jié)果如表5所示。

表5 2020年精度統(tǒng)計分析

圖8顯示,融合后的數(shù)據(jù)集相對IERS事后精密序列的偏差集中在-60～60 μs,相比之下,初始UT1數(shù)據(jù)的偏差集中在-100～100 μs。圖9直觀地展示了融合后的數(shù)據(jù)要比原始數(shù)據(jù)更為穩(wěn)定。通過表5得出NTSC初始UT1標(biāo)準(zhǔn)差為82.51 μs,融合UT1標(biāo)準(zhǔn)差為51.04 μs,融合后的中值、最值均比原始序列更小,經(jīng)過融合算法,得到更連續(xù)、更穩(wěn)定的UT1序列,其精度提高了38%。

圖8 2020年融合前后UT1精度比對Fig.8 Comparison of UT1 accuracy before and after combination in 2020

圖9 2020年融合前后產(chǎn)品比對Fig.9 Comparison of products before and after combination in 2020

需要指出的是,通過比對與MCMC結(jié)合的方法,得到的最優(yōu)權(quán)重比,是一種事后評估行為,實際使用過程中,該最優(yōu)權(quán)比的有效性建立在短期內(nèi)測量精度及頻度基本穩(wěn)定的前提下。若實測精度與頻度發(fā)生改變,則權(quán)比也應(yīng)酌情調(diào)整。在絕大多數(shù)實際應(yīng)用中所采用的權(quán)比,并非理論的最優(yōu)值,而是通過比對和分析,努力從觀測和算法兩方面共同去維持和逼近理論上最優(yōu)的權(quán)比,以得到更穩(wěn)定的UT1序列。

5 結(jié) 論

針對國內(nèi)UT1序列精度提升問題,本文研究了一種基于國內(nèi)數(shù)據(jù)的UT1/ΔLOD融合方法,該方法通過Vondrak方法及MCMC算法進行最優(yōu)權(quán)重的融合。并結(jié)合UT1采樣時間間隔、隨機誤差、系統(tǒng)偏差的仿真實驗來進一步指導(dǎo)NTSC的VLBI2010系統(tǒng)UT1觀測,主要結(jié)論如下。

(1) 對于NTSC的 VLBI2010 系統(tǒng),由于其精度為58.8 μs,通過隨機誤差仿真實驗得出,實測VLBI觀測的采樣間隔為3～5 d時(平均每周2次觀測),融合后的UT1產(chǎn)品精度最高。

(2) 合理觀測頻度、持續(xù)的比對及融合權(quán)比的優(yōu)化,可有效提升最終UT1產(chǎn)品序列的精度、穩(wěn)定度及連續(xù)性。該融合方法對2018年和2020年自主測量的UT1序列精度的提升效果均超過30%。