一種用于空間VLBI的射電源條紋搜索算法

朱曉娜，童 力，鄭為民，張 娟

(1. 中國科學(xué)院上海天文臺，上海 200030；2. 中國科學(xué)院大學(xué)，北京 100049；3. 國家基礎(chǔ)學(xué)科公共科學(xué)數(shù)據(jù)中心，北京 100190；4. 中國科學(xué)院射電天文重點實驗室，江蘇 南京 210033；5. 上海市導(dǎo)航定位重點實驗室，上海 200030)

甚長基線干涉測量是一種具有高分辨率、高精度特點的射電干涉技術(shù)[1]。它采用原子鐘控制的高穩(wěn)定獨立本振系統(tǒng)和數(shù)據(jù)采集裝置，通過多個望遠(yuǎn)鏡同時接收目標(biāo)源信號，記錄數(shù)據(jù)并發(fā)送到異地的數(shù)據(jù)處理中心，經(jīng)處理得到觀測結(jié)果。在甚長基線干涉測量觀測中，兩臺射電望遠(yuǎn)鏡間的距離稱為基線，它決定甚長基線干涉測量的最大分辨率[2]。若有空間望遠(yuǎn)鏡與地面望遠(yuǎn)鏡構(gòu)成空地甚長基線干涉測量系統(tǒng)，基線長度超過地球直徑，可實現(xiàn)更高的分辨率。我國探月工程四期任務(wù)將利用嫦娥七號中繼星上的4.1 m口徑拋物面天線，進(jìn)行月地空間甚長基線干涉測量試驗。不同于常規(guī)地基甚長基線干涉測量[3]，受中繼星軌道擾動和設(shè)備時延的影響，射電源預(yù)報時延模型可能和實際時延有較大差異，無法引導(dǎo)相關(guān)處理機正常工作。通過射電源條紋搜索，可確定實際相關(guān)處理機所需的高精度時延模型，用于空間甚長基線干涉測量空地基線相關(guān)處理。

目前，在探月工程中，探測器信號相關(guān)處理時，采用差分單向測距(Differential One-way Ranging, DOR)信號的點頻信號特征設(shè)計條紋搜索算法[4]。而射電源信號是白噪聲信號，無法利用現(xiàn)有的探測器信號搜索算法。針對射電源信號的條紋搜索有eFFT條紋搜索算法、CAF-W條紋搜索算法[5]和多重網(wǎng)格條紋搜索算法[6]。eFFT條紋搜索算法和CAF-W條紋搜索算法可以實現(xiàn)殘余時延和時延率的計算，但主要用于強源的時延模型搜索。多重網(wǎng)格條紋搜索算法的基本思想是建立多層次、不同尺寸的網(wǎng)格，通過這些網(wǎng)格，在不同范圍內(nèi)搜索需要的目標(biāo)。多重網(wǎng)格條紋搜索算法使用線性時延模型進(jìn)行相關(guān)處理，在時延-時延率平面內(nèi)搜索條紋，適用于信號較強、短時間積分就能得到條紋的應(yīng)用場景。由于空間小口徑天線接收信號較弱，獲得干涉條紋需要進(jìn)行較長積分時間，使用通常的線性時延模型無法保證較長時間范圍內(nèi)的相關(guān)處理機模型時延精度，需構(gòu)造精度更高的二次項時延模型。如果在常規(guī)多重網(wǎng)格中增加時延二次項，會使搜索空間從二維擴展到三維，計算復(fù)雜度增加一個量級。

本文基于多重網(wǎng)格條紋搜索算法，提出了一種用于空間甚長基線干涉測量的射電源條紋搜索算法：首先建立關(guān)于時延和時延率的二維搜索網(wǎng)格，使用每個網(wǎng)格點對應(yīng)的先驗時延和時延率信息構(gòu)建線性時延模型，對原始數(shù)據(jù)分時段相關(guān)處理；然后由二維傅里葉變換求解殘余時延和時延率后進(jìn)行補償，得到每個時段的時延和時延率，計算每個網(wǎng)格點對應(yīng)的相關(guān)幅度值；最后在最大相關(guān)幅度值對應(yīng)的網(wǎng)格點上重構(gòu)二次項時延模型。采用該算法，對俄羅斯RadioAstron空間甚長基線干涉測量的實際觀測數(shù)據(jù)進(jìn)行相關(guān)處理時取得了較好的效果。

1 射電源條紋搜索算法

目前，中國甚長基線干涉測量網(wǎng)相關(guān)機使用五次項時延模型進(jìn)行相關(guān)處理。在條紋搜索過程中，為減少計算復(fù)雜度，通常使用由時延常數(shù)項和時延一次項系數(shù)構(gòu)成的線性時延模型。但隨著時間范圍的擴大，線性時延模型的誤差與五次多項式模型的誤差也會增加。綜合考慮條紋搜索計算量與模型精度因素，我們構(gòu)建由常數(shù)項、一次項系數(shù)和二次項系數(shù)構(gòu)成的二次項時延模型。

1.1 相關(guān)處理算法

假設(shè)兩臺站接收的目標(biāo)源射頻信號分別為s1(t)和s2(t)，同一波前信號到達(dá)兩臺站時延差為τ，則

s2(t)=s1(t+τ).

(1)

由射頻信號得到的基頻信號

x1(t)=s1(t)e-j2πf0t，

(2)

x2(t)=s2(t)e-j2πf0t=s1(t+τ)e-j2πf0t，

(3)

其中，f0為通道的天空頻率。相關(guān)處理以臺站1為參考，使用預(yù)報模型τc對臺站2進(jìn)行時延補償，即

x2(t-τc)=s2(t-τc)e-j2πf0(t-τc)=s1(t+τ-τc)e-j2πf0(t-τc)，

(4)

經(jīng)過條紋旋轉(zhuǎn)后，

x2(t-τc)e-j2πf0τc=s2(t-τc)e-j2πf0t=s1(t+τ-τc)e-j2πf0t.

(5)

對臺站1和臺站2的信號進(jìn)行傅里葉變換，

(6)

S1(f+f0)為臺站1信號經(jīng)過頻率轉(zhuǎn)換后的功率譜，共軛相乘得到兩臺站的互相關(guān)功率譜為

(7)

由(7)式可知，當(dāng)殘余時延τc-τ達(dá)到最小值時，兩臺站的互相關(guān)功率譜達(dá)到最大峰值。

1.2 殘余時延搜索算法

根據(jù)二維傅里葉變換的性質(zhì)可知，相關(guān)處理得到的互相關(guān)譜二維傅里葉變換后，可以得到一個關(guān)于殘余時延和時延率的二維互相關(guān)函數(shù)[7]，如圖1，其公式為

(8)

其中，Δτ為殘余時延；Δτ·為殘余時延率；B為通道帶寬；T為一組二維傅里葉變換的數(shù)據(jù)時長；f0

為天空頻率；F(Δτ, Δτ·)為不同殘余時延和時延率下的相關(guān)幅度值；P12(f,t)為相關(guān)處理后的互相關(guān)功率譜。將(8)式離散化，即

(9)

其中，F(xiàn)(ar,br)為經(jīng)過二維傅里葉變換后得到的相關(guān)幅度值[8]；ar和br分別為離散化后的殘余時延和時延率；N為頻域的傅里葉變換點數(shù)；M為時域的傅里葉變換點數(shù)。

圖1 二維傅里葉變換殘余時延搜索圖Fig.1 Residual delay search of 2-D Fourier transform

1.3 條紋搜索與二次項時延模型重構(gòu)流程

(10)

時延率步長為

(11)

條紋搜索的步長與時延模型精度和計算量密切相關(guān)，如果搜索步長設(shè)置過大，得到的模型精度太低，無法用于后續(xù)數(shù)據(jù)處理；如果步長設(shè)置太小，由于目前的算法采用全網(wǎng)格搜索策略，依次計算每個網(wǎng)格點數(shù)據(jù)會極大影響整體的計算效率。我們注意到步長設(shè)置較小時可能有多組網(wǎng)格點搜索到滿足精度要求的模型，后續(xù)會研究基于二維隨機游走的網(wǎng)格搜索策略，找到合適的時延模型就停止計算，這樣可以有效提高效率。

確定步長后，根據(jù)搜索范圍和先驗時延模型劃分網(wǎng)格點。每個網(wǎng)格點對應(yīng)不同的時延和時延率，代表經(jīng)過不同時延補償?shù)臄?shù)據(jù)處理。如圖3，以先驗時延模型中的時延和時延率為中心點，等步長計算每個網(wǎng)格點對應(yīng)的時延和時延率。

圖2 射電源條紋搜索與二次項時延模型重構(gòu)流程圖Fig.2 Flow chart of radio source fringe search and quadratic term delay model reconstruction

由于弱信號使用線性時延模型進(jìn)行長時間積分無法保證時延精度，獲得干涉條紋，所以將原始數(shù)據(jù)劃為I組較短的數(shù)據(jù)段進(jìn)行處理，每一組為一個二維傅里葉變換組，時間長度為T。在每個網(wǎng)格點上構(gòu)造線性時延對原始數(shù)據(jù)補償后進(jìn)行相關(guān)處理，相關(guān)處理的積分時間為T/M，一個二維傅里葉變換組的數(shù)據(jù)大小為MN。得到兩臺站的互相關(guān)功率譜后，利用(9)式對每一組的相關(guān)譜進(jìn)行二維傅里葉變換，得到每個網(wǎng)格點上的I組殘余時延、時延率和相關(guān)幅度值，利用最大相關(guān)幅度的平均值確定最優(yōu)解對應(yīng)的網(wǎng)格點在該網(wǎng)格點上進(jìn)行時延和時延率的擬合，得到最優(yōu)的二次項時延系數(shù)。

圖3 二維網(wǎng)格示意圖Fig.3 Schematic diagram of two-dimensional grids

(12)

(13)

甚長基線干涉測量觀測具有時間連續(xù)性，時延模型也應(yīng)該是連續(xù)的。但在求解殘余時延的過程中，數(shù)據(jù)分段導(dǎo)致每一組之間的時延具有間斷點，在擬合二次項時延系數(shù)時，需要根據(jù)連續(xù)的約束條件，即保持第1組的時延a0不變，修正其余組的時延ai，使之與前一組末的時延相等，即

(14)

經(jīng)過上述時延改正后，再進(jìn)行二次項時延系數(shù)擬合。設(shè)未知的二次項時延系數(shù)為[x0x1x2]，求解時延的多項式為

(15)

其中，ti=iT為每個二維傅里葉變換組開始的時間。對(15)式右端求導(dǎo)，得到時延率的多項式為

bi=x1+2x2ti.

(16)

聯(lián)立(15)式和(16)式建立超定方程組，再利用正交三角分解法求解[9]。

2 算法驗證

RadioAstron項目是前蘇聯(lián)在20世紀(jì)80年代中期提出的一項空間甚長基線干涉測量計劃。2011年，俄羅斯發(fā)射了10 m口徑的RadioAstron空間射電望遠(yuǎn)鏡(Spektr-R, Ra)衛(wèi)星。Ra與地面多個大型射電望遠(yuǎn)鏡構(gòu)成空間甚長基線干涉測量系統(tǒng)，接收頻段為0.3 GHz，1.6 GHz，5.0 GHz和22 GHz，軌道高度近地點600 km，遠(yuǎn)地點330 000 km，通過聯(lián)合地面望遠(yuǎn)鏡進(jìn)行干涉測量。

為驗證空間甚長基線干涉測量射電源條紋搜索算法，我們使用2014年空間射電望遠(yuǎn)鏡Ra與美國阿雷西博射電望遠(yuǎn)鏡(Arecibo Radio Telescope, Ar)聯(lián)合觀測射電源0823 + 033的數(shù)據(jù)進(jìn)行干涉處理。觀測起始時間為2014年4月11日23時07分00秒，結(jié)束時間為2014年4月11日23時07分59秒。其中，Ar站為2比特量化，4比特扇出的數(shù)據(jù)；空間站Ra為1比特量化，2比特扇出的數(shù)據(jù)；天空頻率4 836 MHz；帶寬16 MHz；積分時間為8.192 × 10-3s。相關(guān)處理后的數(shù)據(jù)規(guī)格為60 × 120 × 8 192，利用每120 × 8 192個點的數(shù)據(jù)為一組進(jìn)行二維傅里葉變換。由(10)式和(11)式確定時延和時延率方向的搜索步長為

(17)

(18)

以先驗時延模型作為搜索的中心點，時延和時延率搜索的網(wǎng)格大小為20 × 200，其中，時延的搜索范圍為6.336 00 × 10-3～8.281 60 × 10-3s，時延率的搜索范圍為1.627 00 × 10-6～4.138 58 × 10-6s/s。算法驗證時，首先采用常規(guī)多重網(wǎng)格條紋搜索算法，僅使用線性時延模型對每個網(wǎng)格點相關(guān)處理，得到峰值點的干涉條紋圖；再使用本文的射電源條紋搜索算法構(gòu)建二次項時延模型；最后對兩者的結(jié)果比較驗證。

圖4(a)為直接利用線性時延模型進(jìn)行條紋搜索的相關(guān)幅度圖，峰值點對應(yīng)的網(wǎng)格位置為(15, 158)，相關(guān)幅度為2.869 00 × 10-4；圖4(b)為峰值點對應(yīng)的條紋相位。圖中沒有明顯的干涉條紋。

圖4 使用條紋搜索和線性時延無法得到正確條紋。(a)相關(guān)幅度圖；(b)干涉條紋圖Fig.4 Fringe search using linear delay can not produce the fringe. (a) Correlation amplitude; (b) interference fringe

圖5為空間甚長基線干涉測量射電源條紋搜索算法得到的不同網(wǎng)格點對應(yīng)的相關(guān)幅度，峰值點對應(yīng)的網(wǎng)格位置為(11， 101)，相關(guān)幅度為1.96700 × 10-3，經(jīng)過擬合后得到二次項時延系數(shù)[7.368 75 × 10-3s， 2.887 84 × 10-6s/s， -4.299 89 × 10-11s/s2]。使用該二次項模型做相關(guān)處理，可得到圖6(a)的干涉條紋，圖6(b)為解卷繞之后的條紋結(jié)果?？梢?，本文提出的算法在處理RadioAstron數(shù)據(jù)時優(yōu)于多重網(wǎng)格條紋搜索算法。

圖5 二次項時延補償?shù)南嚓P(guān)幅度Fig.5 Correlation amplitude using quadratic delay

圖6 采用二次項時延補償后獲得清晰的干涉條紋。(a)解卷繞之前的干涉條紋；(b)解卷繞之后的干涉條紋

該算法原理樣機在Linux集群平臺，采用計算節(jié)點間MPI+節(jié)點內(nèi)POSIX多線程的兩級并行實現(xiàn)，計算節(jié)點的中央處理器型號為英特爾Xeon E5-2620，主頻為2 GHz，每個節(jié)點的中央處理器個數(shù)為12，運行耗時如表1。后續(xù)通過優(yōu)化，運行速度可以進(jìn)一步提高。

表1 空間甚長基線干涉測量射電源條紋搜索算法運行耗時Table 1 Time consumption of the radio source fringe search algorithm for space VLBI

3 結(jié) 論

在未來探月四期任務(wù)中，地月空間甚長基線干涉測量試驗系統(tǒng)受中繼星軌道擾動和設(shè)備時延影響，可能導(dǎo)致射電源預(yù)報時延模型和實際時延有較大差異。本文在常規(guī)多重網(wǎng)格條紋搜索的基礎(chǔ)上，根據(jù)空間甚長基線干涉測量信號弱的特點，提出了一種適用于空間甚長基線干涉測量的射電源條紋搜索算法。該算法先建立關(guān)于時延和時延率的二維搜索網(wǎng)格，使用線性時延相關(guān)處理，然后由二維傅里葉變換搜索殘余時延和時延率，最后重構(gòu)二次項時延模型。本文使用RadioAstron數(shù)據(jù)驗證算法，結(jié)果表明在使用常規(guī)線性時延模型搜索無法實現(xiàn)空地基線條紋的情況下，新算法可以得到清晰的空間甚長基線干涉測量射電源干涉條紋。經(jīng)進(jìn)一步測試和完善，該算法可應(yīng)用于我國后續(xù)空間甚長基線干涉測量。