孫曉彤，童 力，鄭為民,4，余 赟

(1. 中國科學院上海天文臺，上海 200030；2. 中國科學院大學，北京 100049；3. 中國科學院射電天文重點實驗室，江蘇 南京 210008；4. 上海市導航定位重點實驗室，上海 200030)

甚長基線干涉測量是20世紀60年代后期發(fā)展起來的一種射電干涉技術(shù)，具有高分辨率、高精度的特點，在天文觀測、大地測量和深空探測等領(lǐng)域得到了廣泛應用[1]，它的基本工作方式是位于基線兩端的天線接收射電源輻射的電磁波，相關(guān)處理機利用干涉原理從實測數(shù)據(jù)中得到時延以及時延變化率，進一步獲取射電源和有關(guān)基線的位置等信息。因此，快速、準確地獲得射電源的干涉條紋，對甚長基線干涉測量的數(shù)據(jù)處理至關(guān)重要。

相關(guān)處理機是甚長基線干涉測量技術(shù)的核心數(shù)據(jù)處理設(shè)備，預報時延模型與實際時延的誤差直接影響干涉條紋的處理效果，因此，處理機對所需預報時延模型有很高的精度要求。但在實際數(shù)據(jù)處理時，地球模型誤差、大氣模型誤差、儀器誤差以及臺站坐標誤差等都可能導致時延模型不準確甚至無法引導處理機正常工作。在常規(guī)地面甚長基線干涉測量中，由于臺站坐標較穩(wěn)定，可以通過預先計算獲取預報時延模型與實際時延之間的誤差，然后經(jīng)過反復嘗試獲得準確的干涉條紋，在觀測過程中通常不需要再次進行搜索。但對于空間甚長基線干涉測量，觀測系統(tǒng)和軌道不穩(wěn)定，誤差隨時變化，與地面甚長基線干涉測量有較大區(qū)別，因此，無法通過預先計算對時延模型進行一次性修正。目前，針對探測器信號的無時延模型條紋搜索研究較為成熟，上海天文臺在這一領(lǐng)域取得了大量研究成果，解決的主要思路是利用探測器測控信號自身的特征設(shè)計條紋搜索算法[2]。但由于射電源信號是高斯白噪聲信號，無法直接利用信號特征進行條紋搜索。在實際數(shù)據(jù)處理過程中，無時延模型的射電源條紋搜索往往采用半人工的方式，主要依賴操作人員的經(jīng)驗，這種方式費時費力，難以滿足未來空間甚長基線干涉測量的需求。因此，解決不依賴時延模型對射電源條紋進行搜索這一難題，對推動空間甚長基線干涉測量的應用具有重要意義。

針對無時延模型的射電源條紋搜索，本文提出了基于多重網(wǎng)格的射電源條紋搜索算法，能夠在較大的時延和時延率平面內(nèi)進行搜索，利用滑窗技術(shù)以及多重網(wǎng)格技術(shù)快速自動定位到精確的時延、時延率。

1 甚長基線干涉測量相關(guān)處理原理

甚長基線干涉測量觀測臺站接收的射電源信號經(jīng)過變頻、濾波和二次變頻后，被劃分成多個頻率通道，經(jīng)過采樣、量化、編碼后記錄在磁盤上或通過網(wǎng)絡傳輸?shù)缴蹰L基線干涉測量中心[3]。設(shè)臺站1接收的信號為x1，臺站2接收的信號為x2，兩臺站信號之間滿足

x2=x1(t-τ)，

(1)

由于地球自轉(zhuǎn)和射電源的相對運動，τ實際上是一個隨時間變化的函數(shù)，在很短的時間Δt內(nèi)τ可以近似為一次函數(shù)：

(2)

經(jīng)過下變頻后兩臺站的輸出信號分別為

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

(3)

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

(4)

其中，f0為天空頻率。

將原始數(shù)據(jù)解碼后，固定臺站1信號，對臺站2信號做時延補償，設(shè)τs為搜索的補償時延，它與實際時延τr之差為Δτs，用公式表示為

τs=τr-Δτs.

(5)

由于數(shù)據(jù)是以數(shù)字信號的形式存儲，采樣間隔ΔT，則τs由整數(shù)部分nΔT和小數(shù)部分τsf兩部分組成：

τs=nΔT+τsf.

(6)

時延的補償分為兩部分：(1)在時域完成的整數(shù)比特補償；(2)在頻域完成的小數(shù)比特補償。對臺站2信號進行整數(shù)比特補償后為[3]

s2(t-nΔT)=x2(t-nΔT)e-j2πf0(t-nΔT)

=x2(t-τr-τsf+Δτs)e-j2πf0(t-τr+τsf+Δτs).

(7)

(8)

進行傅里葉變換后得到其頻譜S2(f)，完成小數(shù)比特補償后得到

(9)

對臺站1做傅里葉變換，得到兩臺站之間的互功率譜為

(10)

其中，φ1和φ2分別為臺站1和臺站2信號的原始相位。

通過上述分析可以看出，兩臺站互相關(guān)功率譜和補償時延τs與實際時延的差值Δτs、補償時延率

2 多重網(wǎng)格的條紋搜索

基線兩端的天線接收由射電源輻射的電磁波，地球和射電源的相對運動使得電磁波的波前到達兩個天線的時間差不斷改變。干涉條紋只存在于同一波前信號的互相關(guān)功率譜中[4]。因此，只有將接收的兩臺站信號進行準確的時延和時延率補償才可能得到射電源的干涉條紋。本文提出一種基于多重網(wǎng)格技術(shù)的算法，可以在無預報時延模型的情況下對射電源條紋進行搜索。

2.1 多重網(wǎng)格技術(shù)基本原理

多重網(wǎng)格的基本思想是建立多層次、不同尺寸的網(wǎng)格，通過這些網(wǎng)格在不同范圍內(nèi)搜索需要的目標值[5]。如圖1，搜索空間劃分為L層，每層網(wǎng)格由邊長相同的單元格構(gòu)成，單元格的邊長稱為搜索步長，每個單元格的結(jié)點表示不同的搜索值。在劃分過程中，隨著網(wǎng)層號的降低，網(wǎng)格尺寸減小，搜索精度提高。

通過對甚長基線干涉測量相關(guān)原理的分析可知，每一對補償時延和補償時延率與其互相關(guān)功率譜相對應。條紋搜索需要在時延-時延率平面中找出存在條紋的互相關(guān)功率譜對應的時延和時延率。將時延-時延率搜索平面劃分網(wǎng)格，網(wǎng)格尺寸的大小決定了時延和時延率的搜索精度，最底層的網(wǎng)格尺寸由相關(guān)處理機的精度需求決定。

圖1 多重網(wǎng)格示意圖Fig.1 Schematic diagram of multiple grids

甚長基線干涉測量處理機對時延和時延率的補償精度有很高的要求。時延補償?shù)木缺仨殱M足[6]

(11)

其中，W為頻率通道帶寬；N為譜通道個數(shù)(FX 型處理機)或時延通道個數(shù)(XF型處理機)。

時延率的誤差Δτ˙需要滿足

(12)

其中，T為積分時間；f為觀測頻率。

若直接選取滿足相關(guān)處理機精度要求的單一網(wǎng)格，相應的計算量極大，條紋搜索效率很低。為了提高搜索效率，采用粗網(wǎng)格定范圍、細網(wǎng)格定終值的策略。首先選取較大的搜索步長，在粗網(wǎng)格的搜索平面內(nèi)確定正確的時延和時延率補償所在的范圍，如圖2(a)，通過計算確定正確的值落在圖中紅色區(qū)域，用于第2步精搜索，這一過程稱為粗搜索。縮小搜索范圍后，同時縮小網(wǎng)格以滿足相關(guān)處理機的精度要求，如圖2(b)，進而可以在細網(wǎng)格內(nèi)搜索到更高精度的時延和時延率補償，完成精搜索。

圖2 網(wǎng)格操作示意圖。(a)粗網(wǎng)格；(b)細網(wǎng)格Fig.2 Schematic diagram of grid operation. (a) Rough grid; (b) Fine grid

射電源信號的信噪比低，需要增加積分時間解決這一問題，但計算機對相關(guān)計算長度和傅里葉變換計算的長度都有限制[7]，通過對信號的分段處理可平衡兩者之間的關(guān)系。在兩路信號中，需要固定一路信號作為參照信號，對另一路信號進行分段處理，第1路信號與第2路信號每段進行相關(guān)處理即實現(xiàn)滑窗相關(guān)。信號滑窗操作方法如圖3，將臺站信號在搜索范圍內(nèi)以網(wǎng)格要求的搜索步長滑動，并在相關(guān)計算后記錄互相關(guān)功率譜，提取相位后計算方差。

2.2 評估時延和時延率補償結(jié)果的依據(jù)和方法

由(10)式可知，若補償準確，互相關(guān)功率譜的相位與頻率呈線性關(guān)系，斜率表示殘余時延。如圖4(a)表示時延和時延率補償與實際時延和時延率相近的相位結(jié)果，補償不準確的時延和時延率的相頻圖如圖4(b)。

圖3 滑窗操作示意圖Fig.3 Schematic diagram of sliding window operation

從圖4可以看出，互相關(guān)功率譜的相位彌散程度可以衡量時延、時延率補償?shù)臏蚀_性。補償時延和時延率與實際時延和時延率誤差小，相位彌散程度??；補償誤差大，相位彌散程度大。要度量相位的彌散程度需要先對相頻圖進行線性擬合，確定基準點的位置。采用最小二乘法進行線性擬合，如圖4(a)。由于最小二乘的曲線擬合方法是通過最小誤差的平方找到擬合結(jié)果，雖然補償不準確的相頻譜不具備線性關(guān)系，但仍然可以擬合出一條與真實數(shù)據(jù)距離最近的直線，如圖4(b)。

圖4 在相頻圖上的擬合結(jié)果。(a)補償準確的相頻圖擬合結(jié)果；(b)補償不準確的相頻圖擬合結(jié)果

在相頻圖中，有一種特殊情況是相位卷繞。由(10)式可知，互相關(guān)功率譜是共軛相乘，可以得到相位是相加的關(guān)系，表示為

φ(f)=φ1(f)+φ2(f) .

(13)

已知φ1(f)，φ2(f)取值范圍在(-π, π)之間，但兩者相加后取值范圍可能不在(-π, π)之內(nèi)。由于計算機在處理相位時只能用(-π, π)表示，此時相位圖被截斷而發(fā)生跳變[8]，如圖5(a)。從圖5可以看出條紋是存在的，但相位發(fā)生了截斷，出現(xiàn)條紋平行的現(xiàn)象，導致無法對相位數(shù)據(jù)進行準確的線性擬合，過大的擬合方差也影響后續(xù)的判斷。

通過上面的分析可知，因為兩臺站信號互相關(guān)功率譜的相位超出(-π, π)而導致相位卷繞現(xiàn)象。為了準確表示互相關(guān)功率譜的相位，需要對相位進行解卷繞處理。解卷繞方法為規(guī)定跳變閾值，當相位跳變超過閾值時，將相位進行遷移，使截斷的相位連接在一起[9]，解卷繞后的相位如圖5(b)。

圖5 卷繞與解卷繞。(a)相位卷繞；(b)相位解卷繞Fig.5 Phase wrapping and phase unwapping. (a) Phase wrapping; (a) Phase unwrapping

2.3 多重網(wǎng)格條紋搜索流程

圖6是多重網(wǎng)格條紋搜索的具體流程。首先需要計算相關(guān)處理機的精度，確定粗搜索的時延和時延率補償?shù)乃阉鞑介L，讀取臺站記錄的信號，解碼后在絕對時間上對齊，固定臺站1信號，對臺站2信號進行粗搜索的滑窗處理，完成時延和時延率的補償后計算互相關(guān)功率譜，搜索到互相關(guān)功率譜彌散程度最小的位置，縮小搜索范圍，以滿足精度要求的搜索步長進行相關(guān)處理，在信號相關(guān)后得到互相關(guān)功率譜，提取相位。計算每組時延、時延率對應相位的方差并記錄，最后將記錄的方差結(jié)果繪制輸出，選取峰值對應的時延、時延率，最終得到精確的時延、時延率并輸出。

3 實測信號驗證

為了驗證多重網(wǎng)格條紋搜索的正確性，我們以嫦娥四號探月工程中一段實測數(shù)據(jù)進行驗證，其觀測代碼為s8b06a，觀測源為3C 273B。采用FX型相關(guān)處理機。本次實驗選擇天馬站和烏魯木齊站的信號，采樣率為4 MHz, 采用2 bit量化，劃分為16個通道，通道帶寬為2 MHz，積分時間2 s。實驗選取第1通道的信號，中心頻率f0為2 233.5 MHz。由(11)式和(12)式計算可知，時延、時延率補償?shù)木刃枰獫M足：

(14)

(15)

算法本身對時延和時延率的搜索范圍在理論上沒有硬性約束，為了方便，在實驗操作時，以先驗模型為初值確定較大的搜索范圍。在工程應用中，若有更大的搜索范圍需求而帶來更大的計算量，可通過信息傳遞接口(Message Passing Interface, MPI)、統(tǒng)一計算設(shè)備架構(gòu)(Compute Unified Device Architecture, CUDA)等并行計算技術(shù)實現(xiàn)。

首先選擇粗網(wǎng)格進行搜索，時延搜索步長為1 × 10-5s，時延率搜索步長為1 × 10-9s/s。在搜索范圍內(nèi)，分別計算不同的時延和時延率對應的互相關(guān)功率譜的相位方差σ2。為了使結(jié)果更加直觀，取σ2的倒數(shù)Q繪制成圖7(a)，得到峰值位置為補償時延2.88 × 10-3s，補償時延率7.15 × 10-7s/s。因此，時延補償?shù)姆秶鸀?.87 × 10-3～2.89 × 10-3s，時延率補償?shù)姆秶鸀?.14 × 10-7～7.16 × 10-7s/s，在此范圍內(nèi)進行更為精確的條紋搜索，得到時延、時延率補償?shù)木_值。

圖6 多重網(wǎng)格條紋搜索流程圖Fig.6 Multi-grid fringe search flow chart

圖7 多重網(wǎng)格條紋搜索。(a)粗網(wǎng)格搜索圖；(b)細網(wǎng)格搜索圖Fig.7 Multi-grid fringe search. (a) Rough grid fringe search; (b) Fine grid fringe search

經(jīng)過上述操作，確定時延的精搜索范圍為2.87 × 10-3～2.89 × 10-3s，時延率的精搜索范圍為7.14 × 10-7～7.16 × 10-7s/s；精搜索時延步長為1 × 10-6s，時延率步長為1 × 10-10s/s。分別計算不同時延、時延率對應的Q值，繪制成圖7(b)，求得峰值對應的時延為2.881 × 10-3s，時延率為7.151 × 10-7s/s。將搜索到的時延、時延率代入相關(guān)處理機進行驗證，得到射電源的條紋如圖8。

本文實驗使用的中央處理器型號是Intel(R)Xeon(R)CPU X5650@2.67GHz。實測數(shù)據(jù)庫應用中多重網(wǎng)格條紋搜索的窗口大小、程序運行時間等相關(guān)參數(shù)如表1。

表1 多重網(wǎng)格條紋搜索的相關(guān)參數(shù)Table 1 Parameters of model-free fringe searches

圖8 實測數(shù)據(jù)條紋Fig.8 Fringe phases of actual data

4 結(jié) 論

本文針對射電源甚長基線干涉測量自動條紋搜索這一難題，基于多重網(wǎng)格的方法提出了一種射電源信號的條紋搜索算法，并通過嫦娥四號實測信號進行驗證，該方法解決了射電源條紋搜索依賴于時延模型這一難題?；谠撍惴ǖ膽帽尘埃?jīng)過后續(xù)深入的研究，有望應用于空間甚長基線干涉測量的數(shù)據(jù)處理。