基于國內寬帶VLBI2010系統及iGMAS的UT1/ΔLOD數據融合方法

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

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

0 引 言

世界時(universal time, UT1)是以地球自轉定義的時間尺度,是構成標準時間,即協調UT1(UT1 coordinated, UTC)的重要組成部分。UT1與極移、章動改正項組成了地球定向參數(earth orientation parameter, EOP),關于需要在地球坐標和天球坐標之間建立坐標變換的應用,UT1都必不可少,UT1數據服務是衛星導航、深空探測等領域所必須的基礎服務。

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

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

1 融合算法

1.1 理論支撐

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

d(UT1-TAI)/d=ΔLOD

(1)

式中:TAI為國際原子時。

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

1.2 Vondrak算法

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

關于一組測量數據(,) (=1,2,…,),其中是自變量,是對應自變量的測量值,Vondrak算法的理論思想如下:

=+=min

(2)

(3)

(4)

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

(5)

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

(6)

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

(7)

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

(8)

其中系數矩陣,和中的系數為

在原有的Vondrak方法上,增加一階導數擬合度的定義,即Vondrak組合平滑算法:

(9)

(10)

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

1.3 MCMC方法及權重選取

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

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

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

2 輸入數據源

2.1 UT1數據

采用的UT1數據是由NTSC的國內寬帶VLBI2010系統自主測定的UT1序列。

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

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

2.2 ΔLOD數據

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

2.3 事后精密EOP 數據

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

3 仿真分析

為了評估及優化融合算法,同時優化和指導NTSC 的VLBI2010系統UT1測定方案,根據文獻[40]中提到的NTSC的VLBI2010系統實測UT1精度為58.8 μs,且不能連續觀測的現狀,仿真環節進行權重的仿真、時間間隔仿真、隨機誤差仿真及系統偏差的仿真,并研究輸出序列與這些變量間的依賴關系。

3.1 時間間隔采樣仿真

天文觀測會受到惡劣天氣、設備維護、儀器故障乃至網絡中斷等突發狀況導致觀測的不連續和中斷。針對這種現象,需研究不同間隔的UT1數據對融合精度產生的影響。故而仿真了不同間隔的UT1實驗,用于指導實際的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算法,確定最優權重比和最優的融合結果,最后將融合的結果與C04序列做差,統計其標準差和權重比值(UT1/ΔLOD)。如表1所示,列代表對應的采樣間隔,依次為1 d、2 d、3 d、4 d、5 d、6 d,行代表對應的權比和標準差,其中標準差單位為μs。

表1 不同采樣間隔權重比值及標準差

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

3.2 隨機誤差仿真

研究隨機誤差對融合精度的影響,分別在精密C04序列的基礎上加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序列的進行做差,最終統計標準差,42種隨機誤差仿真實驗的結果如表2所示,行代表對應的采樣間隔,列代表對應的隨機誤差值,表中數值即為融合得到的UT1精度值。

表2 誤差仿真

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

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

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

3.3 系統偏差仿真

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

表3 偏差仿真

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

4 實測數據集融合

4.1 2018年實測數據集融合

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

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

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

圖4 2018年UT1數據比對分析Fig.4 Comparison and analysis of UT1 data in 2018

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

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

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

為了更具體地展示融合精度,統計了差值的均值、標準差、中位數、最大值及最小值,如表4所示。

表4 2018年精度統計分析

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

4.2 2020年實測數據集融合

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

圖7 2020年UT1數據比對分析Fig.7 Comparison and analysis of UT1 data in 2020

均值、標準差、中位數、最大值及最小值統計結果如表5所示。

表5 2020年精度統計分析

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

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

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

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

5 結 論

針對國內UT1序列精度提升問題,本文研究了一種基于國內數據的UT1/ΔLOD融合方法,該方法通過Vondrak方法及MCMC算法進行最優權重的融合。并結合UT1采樣時間間隔、隨機誤差、系統偏差的仿真實驗來進一步指導NTSC的VLBI2010系統UT1觀測,主要結論如下。

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

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