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

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

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

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

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

1 射電源條紋搜索算法

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

1.1 相關處理算法

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

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為通道的天空頻率。相關處理以臺站1為參考，使用預報模型τc對臺站2進行時延補償，即

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

(4)

經過條紋旋轉后，

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

(5)

對臺站1和臺站2的信號進行傅里葉變換，

(6)

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

(7)

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

1.2 殘余時延搜索算法

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

(8)

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

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

(9)

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

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

1.3 條紋搜索與二次項時延模型重構流程

(10)

時延率步長為

(11)

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

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

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

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

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

(12)

(13)

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

(14)

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

(15)

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

bi=x1+2x2ti.

(16)

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

2 算法驗證

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

為驗證空間甚長基線干涉測量射電源條紋搜索算法，我們使用2014年空間射電望遠鏡Ra與美國阿雷西博射電望遠鏡(Arecibo Radio Telescope, Ar)聯(lián)合觀測射電源0823 + 033的數(shù)據(jù)進行干涉處理。觀測起始時間為2014年4月11日23時07分00秒，結束時間為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。相關處理后的數(shù)據(jù)規(guī)格為60 × 120 × 8 192，利用每120 × 8 192個點的數(shù)據(jù)為一組進行二維傅里葉變換。由(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)格點相關處理，得到峰值點的干涉條紋圖；再使用本文的射電源條紋搜索算法構建二次項時延模型；最后對兩者的結果比較驗證。

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

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

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

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

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

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

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

3 結 論

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