朱夢(mèng)祎 吳元偉 姚 當(dāng) 弓劍軍 劉 佳李西順 馬浪明 韋 沛 雷 輝 楊旭海

(1 中國(guó)科學(xué)院國(guó)家授時(shí)中心西安 710600)

(2 中國(guó)科學(xué)院大學(xué)北京 100049)

(3 中國(guó)科學(xué)院精密導(dǎo)航定位與定時(shí)技術(shù)重點(diǎn)實(shí)驗(yàn)室西安 710600)

(4 中國(guó)科學(xué)院大學(xué)天文與空間科學(xué)學(xué)院北京 100049)

1 引言

甚長(zhǎng)基線干涉測(cè)量(VLBI)誕生于上世紀(jì)60年代,是一種重要的高精度射電干涉測(cè)量技術(shù)[1].通過(guò)極長(zhǎng)的基線長(zhǎng)度達(dá)到極高的角分辨率,在天體測(cè)量、空間大地測(cè)量、深空探測(cè)、衛(wèi)星定軌及導(dǎo)航等領(lǐng)域有著廣泛的應(yīng)用.傳統(tǒng)雷達(dá)測(cè)距和多普勒測(cè)速對(duì)飛行器僅在視向方向有約束,VLBI可以實(shí)現(xiàn)高精度測(cè)角觀測(cè),提供飛行器的角位置信息,具有對(duì)飛行器的橫向約束能力.美國(guó)國(guó)家航天局(National Aeronautics and Space Administration,NASA)最早將VLBI技術(shù)應(yīng)用于深空探測(cè),在20世紀(jì)70年代VLBI技術(shù)已被應(yīng)用于阿波羅登月工程,同時(shí)開(kāi)展了差分VLBI應(yīng)用于深空探測(cè)器導(dǎo)航的驗(yàn)證研究[2].20世紀(jì)90年代,NASA致力于發(fā)展同波束VLBI技術(shù),以更高的精度來(lái)測(cè)量?jī)蓚€(gè)飛行器之間的相對(duì)位置.干涉測(cè)量技術(shù)近些年開(kāi)始被用于中高軌衛(wèi)星如導(dǎo)航衛(wèi)星的觀測(cè)[3–5].由于VLBI技術(shù)對(duì)衛(wèi)星信號(hào)的依賴較小,干涉測(cè)量技術(shù)在無(wú)源目標(biāo)定軌領(lǐng)域也有重要應(yīng)用價(jià)值.

我國(guó)的上海天文臺(tái)佘山站和新疆天文臺(tái)南山站的兩個(gè)25 m射電望遠(yuǎn)鏡分別于1987年和1993年開(kāi)始進(jìn)行VLBI觀測(cè).在中國(guó)探月系列工程的背景下,中國(guó)的VLBI網(wǎng)(Chinese VLBI Network,CVN)已經(jīng)發(fā)展成為4個(gè)永久天線、一個(gè)移動(dòng)臺(tái)站和一個(gè)相關(guān)中心[6].2013年,中國(guó)科學(xué)院國(guó)家授時(shí)中心按照VLBI2010系統(tǒng)標(biāo)準(zhǔn)[7]設(shè)計(jì)了國(guó)內(nèi)首套寬帶VLBI系統(tǒng)[8],該系統(tǒng)主要用于世界時(shí)UT1的測(cè)定及導(dǎo)航衛(wèi)星、地球同步衛(wèi)星的軌道測(cè)定研究[9–10].

在人造地球衛(wèi)星VLBI跟蹤觀測(cè)方面,上世紀(jì)80年代美國(guó)和日本進(jìn)行過(guò)地球同步衛(wèi)星的VLBI觀測(cè)實(shí)驗(yàn)[11].2000年日本通信情報(bào)研究機(jī)構(gòu)(National Institute of information and Communications Techonlogy,NICT)進(jìn)行寬帶VLBI軟件處理機(jī)研究,并將其用于跟蹤衛(wèi)星觀測(cè)[12].2004年中國(guó)科學(xué)院上海天文臺(tái)舒逢春等[3]利用國(guó)內(nèi)VLBI網(wǎng)跟蹤大橢圓軌道衛(wèi)星,采用基于條紋幅度的加權(quán)最小二乘條紋擬合方法,獲得了衛(wèi)星VLBI觀測(cè)量及其精度估計(jì),完成了衛(wèi)星VLBI觀測(cè)量的閉合誤差檢驗(yàn),經(jīng)系統(tǒng)差改正后,衛(wèi)星VLBI觀測(cè)量序列被用于“探測(cè)1號(hào)”衛(wèi)星的軌道確定.2007年Lanyi等[13]通過(guò)無(wú)線電干涉法確定航天器的角位置,精確度可達(dá)到1–2 nrad.2013年P(guān)lank等[4]利用VLBI觀測(cè)導(dǎo)航衛(wèi)星的測(cè)試及研究表明,在單次觀測(cè)時(shí)延精度為30 ps的前提下,通過(guò)每周觀測(cè)一次導(dǎo)航衛(wèi)星的方法,可實(shí)現(xiàn)毫米級(jí)的臺(tái)站坐標(biāo)確定.2017年P(guān)lank等[5]利用澳大利亞位于Hobart和Ceduna的望遠(yuǎn)鏡,使用VLBI技術(shù)對(duì)GNSS (Global Navigation Satellite System)衛(wèi)星進(jìn)行了觀測(cè)實(shí)驗(yàn),充分利用測(cè)地VLBI已有的方法,開(kāi)發(fā)了適用于衛(wèi)星觀測(cè)的觀測(cè)流程以及數(shù)據(jù)分析方法,以便于支持更多的觀測(cè)實(shí)驗(yàn).2018年弓劍軍等[14]基于中國(guó)科學(xué)院國(guó)家授時(shí)中心的VLBI系統(tǒng),對(duì)中星12號(hào)GEO (Geosynchronous Earth Orbit)衛(wèi)星成功開(kāi)展了觀測(cè)實(shí)驗(yàn),使用DiFX (Distributed FX)和HOPS (Haystack Observatory Postprocessing System)軟件進(jìn)行數(shù)據(jù)相關(guān)及后處理,成功得到干涉條紋,時(shí)延測(cè)量精度約為10 ps.韋沛[15]于2020年進(jìn)行了GEO衛(wèi)星無(wú)源測(cè)定軌系統(tǒng)關(guān)鍵技術(shù)的研究,針對(duì)衛(wèi)星干涉測(cè)量開(kāi)發(fā)了采集系統(tǒng)和軟件,通過(guò)試驗(yàn)表明能有效地采集衛(wèi)星觀測(cè)數(shù)據(jù)并用于干涉時(shí)間測(cè)量,單次測(cè)量精度約為1 ns.

VLBI數(shù)據(jù)的相關(guān)處理及時(shí)延觀測(cè)量提取是通過(guò)相關(guān)機(jī)及后處理程序?qū)崿F(xiàn)的.專業(yè)從事VLBI研究的機(jī)構(gòu)(如SHAO(Shanghai Astronomical Observatory)、NASA、USNO(United States Naval Observatory)、JIVE(Joint Institute for VLBI in Europe)、VERA(VLBI Exploration of Radio Astrometry)、KVN (Korean VLBI Network)等)均有自主研發(fā)的硬相關(guān)及軟相關(guān)處理軟件[16].早期相關(guān)機(jī)為硬件相關(guān)機(jī),硬相關(guān)系統(tǒng)由于硬件固化、穩(wěn)定性低、擴(kuò)展性差,無(wú)法滿足多變的任務(wù)和應(yīng)用場(chǎng)景,調(diào)試和升級(jí)困難.日本在20世紀(jì)80年代采用Fortran研制的XF型軟件相關(guān)處理機(jī)已經(jīng)應(yīng)用于條紋檢測(cè).自2010年前后軟件相關(guān)機(jī)得到飛速發(fā)展及應(yīng)用,逐步替代硬件相關(guān)機(jī).美國(guó)海軍天文臺(tái)、日本國(guó)立天文臺(tái)、韓國(guó)天文研究院、中國(guó)科學(xué)院上海天文臺(tái)等機(jī)構(gòu)也都先后開(kāi)發(fā)了各自的軟相關(guān)處理軟件,但大多為自用閉源軟件.歐洲的VLBI聯(lián)合研究所成功研制了SFXC (Super FX Correlator)軟相關(guān)處理器并將其開(kāi)源化.澳洲學(xué)者Adam Deller和WalterBrisen于2005年開(kāi)發(fā)了開(kāi)源軟相關(guān)處理軟件DiFX,其功能更全、適用性更強(qiáng),美國(guó)VLBA (Very Long Baseline Array)使用的軟相關(guān)處理機(jī)即為在DiFX基礎(chǔ)上二次開(kāi)發(fā)的.

本文介紹基于開(kāi)源軟件二次開(kāi)發(fā)的一套衛(wèi)星干涉測(cè)量數(shù)據(jù)處理系統(tǒng),用于實(shí)現(xiàn)衛(wèi)星干涉測(cè)量數(shù)據(jù)的相關(guān)及后處理.該系統(tǒng)基于開(kāi)源的天文軟件,有學(xué)術(shù)圈和軟件社區(qū)的支持,易于擴(kuò)展、可持續(xù)升級(jí),且不受制于商業(yè)軟件的知識(shí)產(chǎn)權(quán)限制.文中將對(duì)此工具的原理和性能作簡(jiǎn)要介紹.數(shù)據(jù)處理方法部分介紹了本文所采用的3種開(kāi)源軟件包,闡述了數(shù)據(jù)相關(guān)及校準(zhǔn)的方法,給出了實(shí)測(cè)總時(shí)延的計(jì)算方法.結(jié)果及精度評(píng)估部分分別展示了實(shí)測(cè)數(shù)據(jù)的時(shí)延測(cè)量精度及定軌結(jié)果.最后對(duì)本論文所有的研究?jī)?nèi)容和成果進(jìn)行總結(jié),并根據(jù)系統(tǒng)目前存在的問(wèn)題指出了后續(xù)工作的研究方向.

2 數(shù)據(jù)處理方法

2.1 開(kāi)源軟件包

2.1.1 DiFX

DiFX1DiFX參考網(wǎng)址:https://www.atnf.csiro.au/vlbi/dokuwiki/doku.php/difx/start是一種甚長(zhǎng)基線干涉測(cè)量軟件相關(guān)包,是目前軟相關(guān)的標(biāo)準(zhǔn),在天體測(cè)量與測(cè)地VLBI領(lǐng)域的干涉儀數(shù)據(jù)處理中廣泛使用,如用于EHT (Event Horizon Telescope)、LBA (Long Baseline Array)、VLBA、IVS (International VLBI Service for Geodesy and Astrometry)等觀測(cè)數(shù)據(jù)的相關(guān)處理.DiFX軟件相關(guān)機(jī)的開(kāi)發(fā)始于2005年,由學(xué)者Adam Deller與WalterBrisen開(kāi)發(fā),現(xiàn)由澳大利亞、美國(guó)及歐洲的DiFX Developer升級(jí)維護(hù).此軟件開(kāi)源、可擴(kuò)展性強(qiáng)、用戶群大、維護(hù)升級(jí)有保證.自2007年起,每年都會(huì)召開(kāi)User Development Meeting,持續(xù)對(duì)軟件進(jìn)行維護(hù)、升級(jí)和擴(kuò)展,截止2020年11月,軟件最新版本為2.6.1,文中所使用版本為2.5.2,目前版本暫不支持處理衛(wèi)星數(shù)據(jù)[17].

2.1.2 AIPS

美國(guó)國(guó)立射電天文臺(tái)(National Radio Astronomy Observatory,NRAO)開(kāi)發(fā)的天文圖像處理系統(tǒng)(Astronomical Image Processing System,AIPS2AIPS參考網(wǎng)址:http://www.aips.nrao.edu/index.shtml),是一個(gè)用于交互式校準(zhǔn)和編輯射電干涉數(shù)據(jù)以及用于校準(zhǔn)、構(gòu)建、顯示和分析天文圖像的軟件包,代碼開(kāi)源、說(shuō)明文檔完備,是甚長(zhǎng)基線干涉測(cè)量領(lǐng)域(尤其是天體測(cè)量領(lǐng)域)使用最廣泛的軟件包之一.

該軟件可在DEC Alpha (數(shù)字UNIX)、Sun (Solaris)和PC (Linux)等多種計(jì)算機(jī)體系架構(gòu)上運(yùn)行.受自由軟件基金會(huì)的通用公共許可證(GNU General Public License,GPL)保護(hù),可免費(fèi)提供給學(xué)術(shù)界使用.同時(shí)NRAO提供了所有平臺(tái)的源代碼和二進(jìn)制發(fā)行版[18].主要由NRAO的Eric W.Greisen進(jìn)行維護(hù).自2000年起,每年都會(huì)進(jìn)行版本更新,每個(gè)版本隨著debug增添補(bǔ)丁.隨著版本更新,31DEC17及更早版本的發(fā)行版不再維護(hù).目前最新版本為31DEC20,文中所使用版本為31DEC19.

2.1.3 ParselTongue

ParselTongue3ParselTongue參考網(wǎng)址:http://old.jive.nl/jivewiki/doku.php?id=parseltongue:parseltongue是經(jīng)典AIPS、Obit以及其他可能基于任務(wù)的數(shù)據(jù)縮減包的Python接口,允許運(yùn)行AIPS任務(wù),并通過(guò)Python訪問(wèn)AIPS數(shù)據(jù)表頭和擴(kuò)展表;支持運(yùn)行Obit任務(wù)和訪問(wèn)FITS文件中的數(shù)據(jù);還實(shí)現(xiàn)對(duì)AIPS UV數(shù)據(jù)的處理,允許使用現(xiàn)代編程語(yǔ)言編寫(xiě)AIPS腳本,從而實(shí)現(xiàn)復(fù)雜的自動(dòng)化數(shù)據(jù)處理任務(wù).只要滿足Python版本在2.2以上、工作環(huán)境為經(jīng)典AIPS、使用有效的Obit版本這3個(gè)條件,ParselTongue可在任何現(xiàn)代Linux或UNIX上安裝運(yùn)行[19],本文自用版本為最新版本2.3.

2.2 數(shù)據(jù)相關(guān)與校準(zhǔn)方法

數(shù)據(jù)處理主要包括相關(guān)和后處理兩個(gè)環(huán)節(jié),如圖1所示.相關(guān)環(huán)節(jié)使用DiFX,標(biāo)準(zhǔn)DiFX僅能夠處理射電源數(shù)據(jù),需要輸入測(cè)站數(shù)據(jù)、觀測(cè)綱要(VEX格式)以及初始鐘模型,而處理衛(wèi)星觀測(cè)數(shù)據(jù),除測(cè)站數(shù)據(jù)、觀測(cè)綱要及初始鐘模型外,還需要準(zhǔn)備衛(wèi)星的初始幾何時(shí)延模型.處理衛(wèi)星數(shù)據(jù)時(shí),v2d文件中積分時(shí)間和點(diǎn)數(shù)的設(shè)置均與處理射電源數(shù)據(jù)時(shí)不同,處理衛(wèi)星數(shù)據(jù)時(shí)單通道頻點(diǎn)數(shù)量更多.目前所設(shè)置的衛(wèi)星數(shù)據(jù)積分時(shí)間為0.128 s、通道頻點(diǎn)數(shù)量為1024,射電源數(shù)據(jù)積分時(shí)間為1.024 s、通道頻點(diǎn)數(shù)量為64.對(duì)于實(shí)時(shí)任務(wù),鐘模型一般采用前日鐘差數(shù)據(jù)外推的鐘差鐘速信息作為輸入量.初始的衛(wèi)星幾何時(shí)延模型則由衛(wèi)星初始軌道計(jì)算.相關(guān)后的數(shù)據(jù)被存儲(chǔ)為FITS文件格式.后處理環(huán)節(jié)主要基于ParselTongue和AIPS軟件包編寫(xiě),可實(shí)現(xiàn)數(shù)據(jù)的檢查和編輯,電離層、中性大氣、鐘模型、儀器時(shí)延的修正,條紋搜索、帶寬綜合以及總基線時(shí)延/時(shí)延率序列的計(jì)算.對(duì)時(shí)效性要求更高的實(shí)時(shí)任務(wù),因當(dāng)前系統(tǒng)無(wú)法獲取實(shí)時(shí)的大氣時(shí)延和鐘模型,時(shí)延修正環(huán)節(jié)則僅使用預(yù)報(bào)的電離層模型修正電離層時(shí)延,不修正中性大氣和鐘模型引起的時(shí)延.

圖2為經(jīng)過(guò)時(shí)延校準(zhǔn)和條紋搜索后的基線條紋信噪比(Signal-to-Noise Ratio,SNR)結(jié)果.圖3展示了校準(zhǔn)前后北斗C02衛(wèi)星的互相關(guān)功率譜,左圖為3條基線校準(zhǔn)前衛(wèi)星互相關(guān)功率譜,右圖為3條基線校準(zhǔn)后衛(wèi)星互相關(guān)功率譜.如圖所示,北斗C02衛(wèi)星2218 MHz信號(hào)的觀測(cè),衛(wèi)星信號(hào)有效帶寬8 MHz,條紋搜索的時(shí)長(zhǎng)為5 s,以JL站為參考站,JL-KS基線信噪比在100–600范圍內(nèi)變化,集中在300–400范圍內(nèi);JL-SY基線信噪比在100–800范圍內(nèi)變化,集中在300–400范圍內(nèi);3條基線校準(zhǔn)前相位均在0?–360?變化,校準(zhǔn)后相位均在0?±5?附近.

接下來(lái)的2.2.1–2.2.4節(jié)中,將介紹電離層時(shí)延、中性大氣時(shí)延、鐘模型和硬件時(shí)延校準(zhǔn)有關(guān)的技術(shù)細(xì)節(jié),與之對(duì)應(yīng)更詳盡的包含數(shù)據(jù)流、軟件命令的流程圖在附錄A中給出.

2.2.1 中性大氣時(shí)延

中性大氣由對(duì)流層、對(duì)流層及平流層組成,其中對(duì)流層時(shí)延是中性大氣時(shí)延的主要因素.對(duì)流層引起的時(shí)延可以分為兩個(gè)主要部分:干大氣時(shí)延和濕大氣時(shí)延.干大氣時(shí)延是由大氣中氣體的干燥部分引起的,而濕大氣時(shí)延僅由大氣中水蒸氣的變化引起.干大氣時(shí)延約占對(duì)流層總時(shí)延的90%,天頂方向的干大氣時(shí)延通常約為2.3 m[20].濕大氣時(shí)延是由于水蒸氣的存在而引起的,在干旱地區(qū)為幾厘米或更小,在潮濕地區(qū)則為35 cm左右.天頂方向的對(duì)流層時(shí)延使用Saasatamoinen模型計(jì)算,視向方向的對(duì)流層時(shí)延使用NMF映射函數(shù)模型計(jì)算,視向與天頂?shù)膶?duì)流層時(shí)延轉(zhuǎn)換由如下公式計(jì)算,其中c為光速,ε為觀測(cè)仰角,τt為視向?qū)α鲗涌倳r(shí)延,τdry、τwet分別為天頂干大氣時(shí)延和濕大氣時(shí)延,M為映射函數(shù)[21–22].

圖4中左側(cè)為衛(wèi)星視線方向的對(duì)流層時(shí)延,紅色為使用GPS數(shù)據(jù)和PPP (Precise Point Positioning)技術(shù)得到的結(jié)果,藍(lán)色為使用氣象站數(shù)據(jù)計(jì)算的結(jié)果.JL站GPS數(shù)據(jù)與氣象站數(shù)據(jù)計(jì)算視線時(shí)延的差值峰峰值為0.07 ns,KS站GPS數(shù)據(jù)與氣象數(shù)據(jù)的差值峰峰值為0.17 ns,SY站GPS數(shù)據(jù)與氣象數(shù)據(jù)的差值最大可達(dá)0.7 ns.利用PPP技術(shù)計(jì)算的天頂大氣時(shí)延的精度在1–2 cm,而使用氣象數(shù)據(jù)計(jì)算的精度則在10 cm.觀測(cè)站同時(shí)配備氣象站和測(cè)地型GPS接收機(jī),對(duì)于大氣時(shí)延的校準(zhǔn)優(yōu)先使用PPP方法得到的大氣時(shí)延,當(dāng)GPS數(shù)據(jù)無(wú)法獲取時(shí)采用氣象站數(shù)據(jù)替代.

2.2.2 電離層時(shí)延

對(duì)于單頻VLBI觀測(cè),電離層的校準(zhǔn)通常通過(guò)使用全球單層電離層模型計(jì)算得到衛(wèi)星視線方向的電離層時(shí)延[23].本文使用的試驗(yàn)數(shù)據(jù)來(lái)自北斗C02衛(wèi)星的S頻段,因該頻段主要的時(shí)延誤差來(lái)自電離層.為評(píng)估電離層修正的水平,本節(jié)中我們采用了4種不同類型單層電離層模型來(lái)計(jì)算視向方向電離層時(shí)延,對(duì)比不同電離層模型后處理結(jié)果的差異.采用歐洲軌道確定中心(Center for Orbit Determination in Europe,CODE4CODE電離層產(chǎn)品:https://www.aiub.unibe.ch/forschung/code--analysezentrum/global _ionosphere_maps_ produced-by code/index _ger.html)的3種單層電離層模型產(chǎn)品C1P (1日預(yù)報(bào))、C2P (2日預(yù)報(bào))、COD(事后精密)以及IGS (International GNSS Service)多分析中心加權(quán)平均的產(chǎn)品,選取2019年12月3日至5日C02衛(wèi)星的觀測(cè),以IGS產(chǎn)品[24]為基準(zhǔn),將其余3種產(chǎn)品電離層時(shí)延計(jì)算結(jié)果與之比對(duì).通常電離層時(shí)延計(jì)算使用如下公式[25],其中τi為視向電離層時(shí)延,τzenith為天頂電離層時(shí)延,VTEC (Vertical Total Electron Content)為垂直電離層電子濃度,v為觀測(cè)頻率.

由公式可知,電離層計(jì)算結(jié)果與特定頻點(diǎn)頻率相關(guān),但在實(shí)際的數(shù)據(jù)校準(zhǔn)中,并未計(jì)算特定頻點(diǎn)電離層時(shí)延,而是計(jì)算不依賴頻率的色散時(shí)延(Dispersive Delay)[26–27],其定義為波長(zhǎng)為1 m的電磁波在電離層傳播中的電離層時(shí)延,以秒為單位,disp-delay=40.28×1016STEC/c3,其中STEC (Sight Total Electron Content)為以TECU (Total Electron Content Unit)為單位的視線方向總的電子柱密度[25],通過(guò)修正色散時(shí)延,AIPS可實(shí)現(xiàn)對(duì)不同頻點(diǎn)電離層時(shí)延的修正.

圖5中從左至右分別為JL站、KS站、SY站,第1行為各站基于4種產(chǎn)品的電離層視向時(shí)延計(jì)算結(jié)果,第2行為CODE的3種產(chǎn)品與IGS的差值結(jié)果.紅色、黃色、綠色、藍(lán)色分別為IGS、C1P、C2P、COD產(chǎn)品計(jì)算結(jié)果.基于4種產(chǎn)品的電離層時(shí)延JL站變化范圍為1.5–5 ns,KS站為2–4 ns,SY站為1–7.5 ns.C1P、C2P、COD 3種產(chǎn)品與IGS的差值結(jié)果從圖中可直觀看出,3站均為藍(lán)色COD-IGS的差值結(jié)果波動(dòng)最小,即COD與IGS的計(jì)算結(jié)果吻合度最高.

2019年12月3日3站差值的數(shù)值統(tǒng)計(jì)結(jié)果如圖6所示,其中綠色為C1P-IGS、藍(lán)色為C2P-IGS、紅色為COD-IGS.JL、KS、SY 3站COD-IGS差值的標(biāo)準(zhǔn)偏差(Standard Deviation,STD)值分別為8.25×10?2ns、1.23×10?1ns、1.03×10?1ns,差異集中分布在±0.4 ns內(nèi).通過(guò)比較我們發(fā)現(xiàn)C1P、C2P的預(yù)報(bào)模型與事后的實(shí)測(cè)模型間差異大多在1 ns以內(nèi),KS站1日預(yù)報(bào)和2日預(yù)報(bào)電離層結(jié)果與事后電離層結(jié)果吻合度更好,JL和SY站C1P、C2P結(jié)果與事后結(jié)果的偏差稍大.對(duì)于S頻段的實(shí)時(shí)任務(wù),使用C1P、C2P電離層模型可實(shí)現(xiàn)0.5–1 ns精度水平的電離層修正.對(duì)于非實(shí)時(shí)任務(wù),則可使用事后的全球電離層產(chǎn)品修正電離層時(shí)延.

2.2.3 鐘模型修正

在鐘模型的修正環(huán)節(jié),我們比較了兩種鐘模型修正方法:一種是基于PPP鐘差數(shù)據(jù)建立的鐘模型,另一種是基于DBBC (Digital Base Band Converter)的鐘差數(shù)據(jù)建立的鐘模型.PPP鐘模型使用測(cè)站GPS接收機(jī)的GPS數(shù)據(jù)和PPP方法得到的鐘差序列建立,基于DBBC數(shù)據(jù)的鐘模型通過(guò)將DBBC的鐘差數(shù)據(jù)一階線性擬合得到.相關(guān)環(huán)節(jié)使用的鐘模型由前日的鐘差數(shù)據(jù)擬合外推得到,基于兩種數(shù)據(jù)得到的鐘模型扣除相關(guān)環(huán)節(jié)使用的先驗(yàn)鐘模型后得到鐘差改正量.

選取2019年12月4日至5日的數(shù)據(jù),使用上述方法進(jìn)行比對(duì)分析,結(jié)果如圖7所示,左側(cè)黑色、藍(lán)色、紅色分別為初始相關(guān)鐘模型、PPP鐘模型、DBBC鐘模型,右側(cè)紅色與藍(lán)色為DBBC和PPP鐘差數(shù)據(jù)相對(duì)先驗(yàn)鐘模型的差值,黑線為對(duì)藍(lán)色數(shù)據(jù)點(diǎn)的連續(xù)分段線性擬合,綠色為紅色數(shù)據(jù)點(diǎn)的一階線性擬合.由于FITS和DBBC的鐘差在10?6量級(jí),而PPP鐘差數(shù)據(jù)在10?4量級(jí),圖中PPP數(shù)據(jù)JL、KS、SY各站分別補(bǔ)償了?2.504×10?4s、?9.832×10?5s、1.204×10?4s.由圖可知,基于兩種數(shù)據(jù)計(jì)算的鐘模型基本趨勢(shì)一致,使用兩個(gè)鐘模型修正基線殘余時(shí)延的結(jié)果比對(duì)見(jiàn)3.1節(jié).

2.2.4 硬件設(shè)備時(shí)延

信號(hào)在接收機(jī)、數(shù)采設(shè)備、記錄設(shè)備傳輸過(guò)程,不同臺(tái)站使用不同的設(shè)備,帶來(lái)不同硬件設(shè)備時(shí)延.該項(xiàng)時(shí)延的修正通過(guò)觀測(cè)射電源的方式確定.觀測(cè)期間每小時(shí)觀測(cè)一次射電源.射電源數(shù)據(jù)的校準(zhǔn)流程與衛(wèi)星流程基本相同,依次對(duì)電離層時(shí)延、中性大氣時(shí)延、鐘模型進(jìn)行校準(zhǔn),條紋搜索后,可分別得到衛(wèi)星和射電源的基線殘余時(shí)延.假設(shè)臺(tái)站坐標(biāo)準(zhǔn)確、射電源坐標(biāo)精確,因電離層時(shí)延、中性大氣時(shí)延和鐘模型已經(jīng)修正,則射電源的基線殘余時(shí)延即是接收機(jī)、采集、記錄設(shè)備引起的硬件時(shí)延.將射電源殘余時(shí)延從衛(wèi)星的殘余時(shí)延中扣除,即可得到校準(zhǔn)了電離層時(shí)延、中性大氣時(shí)延、鐘和硬件設(shè)備時(shí)延的衛(wèi)星基線殘余幾何時(shí)延.利用該序列與衛(wèi)星初始幾何時(shí)延模型,計(jì)算得到衛(wèi)星基線總的幾何時(shí)延/時(shí)延率序列,公式見(jiàn)2.3節(jié).

2.3 總時(shí)延計(jì)算方法

DiFX的初始幾何時(shí)延為測(cè)站相對(duì)地心的幾何時(shí)延,而最終計(jì)算得到的幾何時(shí)延為兩個(gè)測(cè)站間的基線幾何時(shí)延,由地心-測(cè)站幾何時(shí)延到測(cè)站-測(cè)站幾何時(shí)延轉(zhuǎn)換的過(guò)程,公式如下:

其中,T為協(xié)調(diào)世界時(shí)(UTC)的時(shí)刻,為A-B基線幾何時(shí)延,對(duì)于A-B基線,其中A站為參考站(以A站站鐘為基線時(shí)延的時(shí)間戳).為A、B站地心時(shí)延之差;為某時(shí)刻經(jīng)過(guò)地心的波前,該波前經(jīng)過(guò)A站和B站的模型時(shí)延(包括幾何時(shí)延、大氣時(shí)延,但不包括時(shí)鐘)(通常);分別為對(duì)時(shí)間的1階、2階導(dǎo)數(shù);分別為對(duì)時(shí)間的1階、2階導(dǎo)數(shù).公式(4)、(5)詳細(xì)推導(dǎo)方式見(jiàn)附錄B.

3 實(shí)測(cè)數(shù)據(jù)的結(jié)果展示及精度評(píng)估

3.1 基線殘余時(shí)延及精度評(píng)估

圖8為3基線校準(zhǔn)后的殘余時(shí)延.黑色為實(shí)時(shí)校準(zhǔn)結(jié)果,藍(lán)色為事后使用PPP數(shù)據(jù)校準(zhǔn)鐘差的后處理結(jié)果,紅色為使用DBBC數(shù)據(jù)校準(zhǔn)鐘差的后處理結(jié)果,由圖可看出3種處理結(jié)果的殘余時(shí)延趨勢(shì)一致,紅色DBBC和實(shí)時(shí)結(jié)果有10–20 ns的明顯偏差,該偏差部分源于實(shí)時(shí)校準(zhǔn)模式未校準(zhǔn)大氣時(shí)延,部分來(lái)自事后校準(zhǔn)的鐘模型的修正.從圖中可看出,JL-SY基線Real-time峰峰值約為20 ns、PPP峰峰值約為15 ns、DBBC峰峰值約為5 ns,KS-SY基線Real-time峰峰值約為15 ns、PPP峰峰值約為8 ns、DBBC峰峰值約為7 ns,使用DBBC數(shù)據(jù)進(jìn)行修正的結(jié)果比使用PPP修正的結(jié)果要好.JL-SY基線和KS-SY基線對(duì)于未修正的實(shí)時(shí)結(jié)果,可以明顯看到一個(gè)殘余的鐘速項(xiàng),此項(xiàng)是因?yàn)閷?shí)時(shí)任務(wù)中使用的鐘模型不準(zhǔn)引起的,通過(guò)DBBC修正后,該項(xiàng)的起伏明顯降低有所改善.而PPP鐘模型無(wú)法修正來(lái)自儀器和硬件的“鐘飄”,由于DBBC鐘模型同時(shí)還包含有部分信號(hào)鏈路(包括接收機(jī)、數(shù)采、記錄設(shè)備)的硬件時(shí)延,因此整體上采用DBBC鐘差數(shù)據(jù)擬合的鐘模型比采用PPP鐘差數(shù)據(jù)擬合的鐘模型修正的基線殘余時(shí)延要更平穩(wěn),起伏更小.

為定量評(píng)估最終得到的時(shí)延序列的測(cè)量精度,我們對(duì)圖8所示的紅色數(shù)據(jù)點(diǎn),以DBBC鐘差數(shù)據(jù)修正鐘模型得到的衛(wèi)星殘余時(shí)延序列做了樣條函數(shù)擬合.通過(guò)調(diào)整擬合參數(shù),使擬合曲線反映殘余時(shí)延序列的趨勢(shì)項(xiàng),殘差項(xiàng)用于評(píng)估總時(shí)延序列中由測(cè)量噪聲和校準(zhǔn)誤差共同貢獻(xiàn)的“測(cè)量”熱噪聲水平.結(jié)果如圖9所示,右側(cè)為3基線殘余時(shí)延扣除樣條函數(shù)擬合趨勢(shì)項(xiàng)后的擬合殘差.該殘差均在±1.5 ns內(nèi),標(biāo)準(zhǔn)偏差為0.4 ns.這里需要注意的是,該殘差水平和樣條函數(shù)擬合的參數(shù)是存在相關(guān)性的,理論上當(dāng)樣條曲線擬合的控制點(diǎn)足夠多時(shí),殘差將會(huì)逼近0.因此在擬合過(guò)程中,需要合理選擇和調(diào)整擬合參數(shù),使得左圖中的紅色趨勢(shì)項(xiàng)不出現(xiàn)過(guò)擬合.左圖中紅色趨勢(shì)項(xiàng),包含了衛(wèi)星初始軌道偏差導(dǎo)致的真實(shí)幾何時(shí)延偏差以及未完全修正的系統(tǒng)性偏差,如大氣、電離層、儀器乃至測(cè)站坐標(biāo)不準(zhǔn)引起的系統(tǒng)性偏差.圖9基本反映了當(dāng)前我們得到的衛(wèi)星干涉測(cè)量時(shí)延序列精度水平.

3.2 定軌結(jié)果

圖10為時(shí)延數(shù)據(jù)定軌結(jié)果與iGMAS (International GNSS Monitoring and Assessment System)給出的事后精密軌道間的軌道偏差,橫坐標(biāo)為簡(jiǎn)化儒略日(Modified Julian Day,MJD).上、中、下3幅圖分別為日期MJD=[58830.5,58832]、MJD=[58831.5,58833]、MJD=[58832.5,58833.5] 3個(gè)弧段.對(duì)應(yīng)弧段徑向(R)、軌道切向(T)、軌道法向(N)和3維位置(pos)的RMS分別為32.48 m、29.68 m、52.86 m和68.77 m.該結(jié)果引自韋沛博士論文[16]中北斗衛(wèi)星的定軌結(jié)果.

4 總結(jié)與展望

本文簡(jiǎn)要回顧了衛(wèi)星干涉測(cè)量及數(shù)據(jù)處理領(lǐng)域的現(xiàn)狀,并基于已有的天文開(kāi)源軟件二次開(kāi)發(fā)了一套衛(wèi)星干涉測(cè)量數(shù)據(jù)處理系統(tǒng).使用該系統(tǒng)對(duì)國(guó)家授時(shí)中心的北斗衛(wèi)星干涉測(cè)量數(shù)據(jù)進(jìn)行處理和定軌,測(cè)試表明,該系統(tǒng)對(duì)于人造地球衛(wèi)星的干涉測(cè)量數(shù)據(jù)具有良好的處理效果,但在海量數(shù)據(jù)的自動(dòng)批處理上仍有改進(jìn)的空間,同時(shí)需要開(kāi)展更多類型的觀測(cè)試驗(yàn),以精準(zhǔn)定位各種系統(tǒng)性偏差的來(lái)源,并予以改進(jìn)和改正.后續(xù)將進(jìn)一步考慮優(yōu)化算法和代碼效率,提高校準(zhǔn)精度.

附錄

A 衛(wèi)星干涉測(cè)量數(shù)據(jù)相關(guān)及后處理細(xì)節(jié)流程

如圖A.1所示,衛(wèi)星干涉測(cè)量數(shù)據(jù)相關(guān)及后處理流程,依次包含以下7個(gè)步驟:

(1)DiFX軟相關(guān)處理機(jī)對(duì)于衛(wèi)星數(shù)據(jù)的處理,需要預(yù)先準(zhǔn)備衛(wèi)星軌道文件、鐘模型以及觀測(cè)綱要文件.輸出結(jié)果為標(biāo)準(zhǔn)的FITS格式干涉儀數(shù)據(jù);

(2)FITS數(shù)據(jù)經(jīng)由AIPS的FITLD命令讀入程序,TASAV命令保存原始數(shù)據(jù)表格文件,POSSM命令用于生成原始數(shù)據(jù)的互相關(guān)功率譜(用于數(shù)據(jù)質(zhì)量檢查).同時(shí)將初始CL1表格輸出;

(3)GPS數(shù)據(jù)經(jīng)由Bernese軟件的PPP模塊計(jì)算生成各臺(tái)站天頂大氣時(shí)延序列.使用天頂大氣時(shí)延序列和源的俯仰方位信息來(lái)計(jì)算視線方向的大氣時(shí)延,并將視線方向大氣時(shí)延寫(xiě)入CL2表的ATMOS列.同時(shí)利用全球電離層模型及衛(wèi)星軌道,計(jì)算衛(wèi)星視線方向的電離層時(shí)延,并將該時(shí)延序列寫(xiě)入CL2表的DDISP列.此處CL2表的生成使用自己編寫(xiě)的代碼實(shí)現(xiàn).通過(guò)TBIN命令載入CL2表,用于校準(zhǔn)原始觀測(cè)數(shù)據(jù),其次通過(guò)TECOR命令校準(zhǔn)射電源的電離層時(shí)延,也寫(xiě)入CL2表.至此電離層和中性大氣時(shí)延校準(zhǔn)完成,CL2表包含對(duì)衛(wèi)星和射電源的中性大氣時(shí)延及電離層時(shí)延的校準(zhǔn)信息;

(4)使用正文2.2.3節(jié)中描述的方法計(jì)算鐘模型改正文件,將該改正文件作為CLCOR命令的輸入,實(shí)現(xiàn)對(duì)干涉數(shù)據(jù)的鐘模型修正;

(5)對(duì)修正了介質(zhì)傳播時(shí)延和鐘模型時(shí)延的數(shù)據(jù)做條紋搜索,搜索結(jié)果存入SN1表.使用SN1表升級(jí)CL3表得到CL4表,補(bǔ)償殘余時(shí)延后干涉數(shù)據(jù)的相位被拉平.使用POSSM命令生成校準(zhǔn)后的衛(wèi)星和射電源互相關(guān)功率譜(用于校準(zhǔn)質(zhì)量的評(píng)判);

(6)SN1表包含衛(wèi)星和射電源的殘余時(shí)延/時(shí)延率信息.im-file則是相關(guān)環(huán)節(jié)使用的衛(wèi)星和射電源初始幾何時(shí)延模型.因電離層時(shí)延、中性大氣時(shí)延、鐘模型已修正,射電源殘余時(shí)延可以當(dāng)作是未修正的系統(tǒng)性時(shí)延,主要成分為儀器硬件時(shí)延.將射電源的殘余時(shí)延序列插值至衛(wèi)星殘余時(shí)延序列的時(shí)刻,得到的序列包含了衛(wèi)星數(shù)據(jù)的硬件儀器時(shí)延以及未完全修正的其他系統(tǒng)性時(shí)延.將該修正項(xiàng)從衛(wèi)星殘余時(shí)延序列中扣除,得到衛(wèi)星的純幾何殘余時(shí)延;

(7)將衛(wèi)星的純幾何殘余時(shí)延與衛(wèi)星初始幾何時(shí)延模型相加,同時(shí)補(bǔ)償時(shí)延率等修正項(xiàng),計(jì)算得到衛(wèi)星的基線幾何時(shí)延/時(shí)延率序列(干涉測(cè)量的觀測(cè)量),見(jiàn)附錄B.

B 衛(wèi)星基線純幾何時(shí)延計(jì)算方法

初始幾何時(shí)延模型為臺(tái)站相對(duì)于地心的純幾何時(shí)延,而定軌需要得到基于觀測(cè)殘余時(shí)延修正后的基線純幾何時(shí)延.本節(jié)介紹由初始的地心幾何時(shí)延模型、殘余時(shí)延/時(shí)延率序列,計(jì)算基線純幾何時(shí)延/時(shí)延率模型的方法.

B.1 定義

tgc:波前經(jīng)過(guò)地心的時(shí)刻,假設(shè)該時(shí)刻為理想的UTC時(shí);

tA:波前經(jīng)過(guò)臺(tái)站A時(shí)參考站A上的時(shí)鐘時(shí)刻;

tB:波前經(jīng)過(guò)臺(tái)站B時(shí)參考站B上的時(shí)鐘時(shí)刻;

cA(tgc):tgc時(shí)刻A站的時(shí)鐘偏移量,使tA?cA為理想的UTC時(shí);

cB(tgc):tgc時(shí)刻B站的時(shí)鐘偏移量,使tB?cB為理想的UTC時(shí);

(tgc):tgc時(shí)刻經(jīng)過(guò)地心的波前,該波前經(jīng)過(guò)站A時(shí)的模型時(shí)延(包括幾何時(shí)延、大氣時(shí)延,但不包括時(shí)鐘),通常<0;

τ′B(tgc):tgc時(shí)刻經(jīng)過(guò)地心的波前,該波前經(jīng)過(guò)站B時(shí)的模型時(shí)延(包括幾何時(shí)延、大氣時(shí)延,但不包括時(shí)鐘),通常<0;

τAB:臺(tái)站A到臺(tái)站B的時(shí)延(包括A、B的鐘差),τAB≡tB?tA.

B.2 以地心為基準(zhǔn)的先驗(yàn)時(shí)延

在UTC時(shí)的T時(shí)刻,通過(guò)地心的波前,該波前通過(guò)測(cè)站A、B的時(shí)刻如表B.1所示.

表B.1 以經(jīng)過(guò)地心為基準(zhǔn)的波前經(jīng)過(guò)各點(diǎn)的UTC時(shí)間Table B.1 UTC time at each point of the wavefront passing through the center of the earth

因此,UTC時(shí)T時(shí)刻,經(jīng)過(guò)地心的波前先后經(jīng)過(guò)A、B兩站,B站相對(duì)于A站的基線時(shí)延tB?tA可由下式表示.

B.3 以臺(tái)站為基準(zhǔn)的先驗(yàn)時(shí)延

在UTC時(shí)的T時(shí)刻,波前經(jīng)過(guò)臺(tái)站的模型中,波前通過(guò)臺(tái)站和地心的時(shí)刻如表B.2所示:

表B.2 以經(jīng)過(guò)臺(tái)站為基準(zhǔn)的波前經(jīng)過(guò)各點(diǎn)的UTC時(shí)間Table B.2 UTC time at each point of the wavefront passing through the center of the station

定義?t ≡tgc?T≡?.?t以波前經(jīng)過(guò)臺(tái)站A為參考時(shí)刻,該波前再經(jīng)過(guò)地心時(shí),相對(duì)于經(jīng)過(guò)臺(tái)站A的時(shí)延;因參考點(diǎn)從地心變化至臺(tái)站A,?t近似但不完全等于?是T時(shí)刻經(jīng)過(guò)地心的波前,該波前經(jīng)過(guò)A站的時(shí)刻相對(duì)T的時(shí)延.由此,可得A-B基線的先驗(yàn)基線時(shí)延為:

其中所有變量在T時(shí)刻進(jìn)行估計(jì),對(duì)于地球上的固定天線,結(jié)果精度小于1 fs.注意?t不取決于臺(tái)站的時(shí)鐘,它的值僅由傳播時(shí)延及其時(shí)延率決定.

結(jié)果可達(dá)到亞飛秒精度.

B.4 以臺(tái)站為基準(zhǔn)的先驗(yàn)時(shí)延率

可以簡(jiǎn)單地在通過(guò)將(B2)式對(duì)時(shí)間求導(dǎo),來(lái)推導(dǎo)以臺(tái)站為參考和以地心為參考的先驗(yàn)時(shí)延率之間的關(guān)系.

該公式可達(dá)到1×10?22精度.

致謝感謝上海天文臺(tái)舒逢春老師關(guān)于基線時(shí)延計(jì)算方法的討論,感謝國(guó)家授時(shí)中心王源昕博士協(xié)助GPS數(shù)據(jù)的大氣時(shí)延計(jì)算.本文的數(shù)據(jù)處理系統(tǒng)基于DiFX、AIPS、ParselTongue 3種開(kāi)源天文軟件包,GPS數(shù)據(jù)的處理使用了Bernese軟件的PPP功能.文中使用了CODE和IGS的事后和預(yù)報(bào)的全球單層電離層網(wǎng)格數(shù)據(jù),使用了iGMAS的精密衛(wèi)星軌道產(chǎn)品,在此一并致謝.