劉麗佳，彭 勃

(中國(guó)科學(xué)院國(guó)家天文臺(tái)，北京 100012)

行星際閃爍(Interplanetary Scintillation，簡(jiǎn)稱IPS)是指來自遙遠(yuǎn)致密射電源的電磁波在穿過太陽系行星際空間時(shí)，被太陽風(fēng)等離子體散射形成的射電波強(qiáng)度和相位的一種隨機(jī)起伏現(xiàn)象。通過觀測(cè)這種閃爍現(xiàn)象，可以研究致密射電源的角徑，也可以測(cè)量太陽風(fēng)的速度和研究太陽風(fēng)等離子的不規(guī)則結(jié)構(gòu)[1-4]，是兼有理論研究和實(shí)際應(yīng)用價(jià)值的重要課題。

國(guó)際上英國(guó)、前蘇聯(lián)和日本都相繼開展了行星際閃爍測(cè)量用以研究太陽風(fēng)及瞬變現(xiàn)象[5-7]，印度上世紀(jì)70年代也開始了行星際閃爍觀測(cè)[8]。日本的三站觀測(cè)系統(tǒng)可直接測(cè)出投影太陽風(fēng)風(fēng)速，印度的Ooty射電望遠(yuǎn)鏡利用單站單頻模式觀測(cè)致密星體，通過譜擬合的方式得到的閃爍功率譜可測(cè)得太陽風(fēng)風(fēng)速及閃爍指數(shù)。

我國(guó)的行星際閃爍觀測(cè)研究開始于上世紀(jì)90年代末期[9]，北京天文臺(tái)利用密云綜合孔徑射電望遠(yuǎn)鏡的復(fù)合干涉儀對(duì)行星際閃爍觀測(cè)進(jìn)行了初步試驗(yàn)和研究，根據(jù)所得閃爍頻譜得到了太陽風(fēng)風(fēng)速及閃爍指數(shù)，并與同期印度Ooty觀測(cè)數(shù)據(jù)進(jìn)行了對(duì)比，符合較好。但由于采用單站單頻模式進(jìn)行觀測(cè)Fresnel膝不明顯，不好確定對(duì)應(yīng)的頻率，對(duì)估測(cè)速度造成困難。此種方式對(duì)模型選取依賴較大，存在一定誤差。

此后國(guó)際上發(fā)展了單站雙頻模式對(duì)行星際閃爍進(jìn)行觀測(cè)[10-11]。這種模式通過對(duì)歸一化互相關(guān)譜第一零點(diǎn)的測(cè)量來計(jì)算太陽風(fēng)風(fēng)速，同時(shí)也可以研究電子密度起伏的空間譜，比起單站單頻模式準(zhǔn)確性更高。

中國(guó)科學(xué)院國(guó)家天文臺(tái)正在密云基地50m射電望遠(yuǎn)鏡上建立一套單站雙頻行星際閃爍觀測(cè)系統(tǒng)。此系統(tǒng)服務(wù)于國(guó)家重大科學(xué)項(xiàng)目——子午工程。該系統(tǒng)建成投入使用后，將使國(guó)家天文臺(tái)密云基地成為國(guó)際上唯一的多組雙頻行星際閃爍觀測(cè)站。

由于執(zhí)行探月任務(wù)，目前密云基地只能在較高頻(S/X波段)進(jìn)行單站單頻的觀測(cè)。針對(duì)目前設(shè)備情況，作者研究編寫了開展行星際閃爍單站單頻觀測(cè)的數(shù)據(jù)處理軟件。本文的第一部分介紹行星際閃爍單站單頻觀測(cè)的相關(guān)理論，第二部分討論單站單頻數(shù)據(jù)處理，最后對(duì)數(shù)據(jù)處理軟件進(jìn)一步分析。

1 行星際閃爍單站單頻模式

行星際閃爍單站單頻模式是指在一個(gè)觀測(cè)站使用單一頻率進(jìn)行IPS觀測(cè)，用這種模式測(cè)量太陽風(fēng)風(fēng)速有兩種方法，第一種稱為參數(shù)擬合法，第二種稱為特征頻率法。參數(shù)擬合法依賴于太陽風(fēng)的幾個(gè)主要參數(shù)：太陽風(fēng)電子密度功率譜的冪指數(shù)α，各向異性軸比AR，太陽風(fēng)風(fēng)速V，射電源角尺度θ0，選擇適當(dāng)?shù)奶栵L(fēng)模式，計(jì)算出模式譜，使其與實(shí)測(cè)譜吻合，從而確定太陽風(fēng)和射電源的上述參數(shù)。第二種方法是通過實(shí)測(cè)譜確定兩個(gè)特征頻率：Fresnel頻率fF和第一最小頻率fmin，根據(jù)兩個(gè)特征頻率中的任一個(gè)來計(jì)算太陽風(fēng)風(fēng)速，公式如下[5]：

(1)

(2)

其中λ是觀測(cè)頻率；Z是太陽與射電源視線方向之間的最短距離。圖1為一行星際閃爍單站單頻閃爍譜示例。

圖1 單站單頻閃爍譜示例Fig.1 IPS spectrum observed at a single station for single frequencies

2 單站單頻數(shù)據(jù)處理

2.1 數(shù)據(jù)處理介紹

圖2 數(shù)據(jù)處理流程圖Fig.2 Flow chart of data reduction

為了得到行星際閃爍譜，用Visual Fortran語言編寫了針對(duì)單站單頻觀測(cè)模式的數(shù)據(jù)處理軟件，具體流程如圖2。圖中所示的觀測(cè)數(shù)據(jù)采集和記錄系統(tǒng)是一個(gè)實(shí)時(shí)顯示系統(tǒng)。在觀測(cè)過程中數(shù)據(jù)可以被實(shí)時(shí)顯示在電腦顯示屏上，使觀測(cè)數(shù)據(jù)的質(zhì)量實(shí)時(shí)可見。在遇到強(qiáng)干擾時(shí)，可以通過調(diào)節(jié)望遠(yuǎn)鏡接收機(jī)終端上的步進(jìn)衰減或者輻射計(jì)的增益來控制數(shù)據(jù)的輸出數(shù)值范圍，將數(shù)據(jù)控制在線性范圍內(nèi)，以免數(shù)據(jù)因?yàn)楦蓴_而溢出產(chǎn)生失真。數(shù)據(jù)記錄軟件記錄范圍為0～65535，線性范圍為0～32000。為使數(shù)據(jù)不失真，當(dāng)天線不對(duì)準(zhǔn)目標(biāo)源(off-source)觀測(cè)時(shí)，調(diào)整步進(jìn)衰減和輻射計(jì)增益，使數(shù)據(jù)顯示在2000～5000的范圍內(nèi)；天線對(duì)準(zhǔn)目標(biāo)源(on-source)觀測(cè)時(shí)，用采集軟件觀察采集的數(shù)據(jù)，使數(shù)據(jù)不超過30000。

數(shù)據(jù)處理軟件濾波去干擾包括3個(gè)步驟：原始數(shù)據(jù)成圖、去掉緩變成分和去野點(diǎn)。圖3是一段密云50m射電望遠(yuǎn)鏡在S波段觀測(cè)3C48的原始數(shù)據(jù)圖示例。

選擇數(shù)據(jù)進(jìn)行去掉緩變成分處理方法是：根據(jù)采樣頻率、行星際閃爍現(xiàn)象的性質(zhì)和觀測(cè)數(shù)據(jù)情況，選取一定長(zhǎng)度數(shù)據(jù)，用滑動(dòng)平均方法計(jì)算出每段數(shù)據(jù)的平均值，然后用原始數(shù)據(jù)減去每段數(shù)據(jù)的平均值。在本圖中首先選取每段數(shù)據(jù)的長(zhǎng)度為30s，進(jìn)行一次試去緩變成分處理。密云觀測(cè)時(shí)積分時(shí)間為0.1s，因此每段數(shù)據(jù)的長(zhǎng)度為300個(gè)點(diǎn)。

進(jìn)行去緩變成分處理之后是去野點(diǎn)處理。所謂野點(diǎn)是由于干擾等因素出現(xiàn)的高出正常值的數(shù)據(jù)。野點(diǎn)的存在會(huì)使閃爍譜結(jié)果不準(zhǔn)確。去野點(diǎn)的方法是：根據(jù)采樣頻率和實(shí)際觀測(cè)數(shù)據(jù)情況選擇一定長(zhǎng)度的數(shù)據(jù)，計(jì)算這段數(shù)據(jù)的均方根(rms)，然后使這段數(shù)據(jù)中的每一個(gè)值同一個(gè)設(shè)定的閾值相比較(通常是選擇3rms或5rms)，絕對(duì)值大于這個(gè)閾值的點(diǎn)被其前后兩個(gè)點(diǎn)的平均值所取代。圖4是對(duì)原始數(shù)據(jù)去緩變并同3rms比較后所得。

圖3 密云50m射電望遠(yuǎn)鏡IPS觀測(cè)3C48原始數(shù)據(jù)Fig.3 Raw IPS data toward 3C48 obtained with the Miyun 50m radio telescope；x axis:recorder moments;y axis:ralative values of voltages

圖4 IPS數(shù)據(jù)去緩變、野點(diǎn)示例(3rms比較)Fig.4 Plot for IPS after the subtraction of slow changing component and the removal of excessive changing values (“wild point”)

去野點(diǎn)過程中選擇半分鐘的數(shù)據(jù)來計(jì)算rms值，閾值選定為3rms，比較后，共去掉638個(gè)點(diǎn)，占數(shù)據(jù)總量的1.47%。

軟件去干擾后是濾波再采樣、減去平均值和三角加權(quán)。設(shè)軟件去干擾后的數(shù)據(jù)為DATA。在減去平均值這個(gè)步驟中選擇一定長(zhǎng)度的窗函數(shù)，計(jì)算這段數(shù)據(jù)的長(zhǎng)度，然后減去平均值。在減去平均值時(shí)采用滑動(dòng)平均的方法(同軟件去干擾中去掉緩變成分類似)?；瑒?dòng)平均值同樣選擇一定長(zhǎng)度的數(shù)據(jù)進(jìn)行平均。記減去平均值后的數(shù)據(jù)為DATA1。減去平均值后是三角加權(quán)，在這一步中選擇與上一步相同長(zhǎng)度的三角函數(shù)，三角函數(shù)的兩端為0，中間數(shù)值最大，兩側(cè)對(duì)稱，使三角函數(shù)與DATA1相乘，所得數(shù)據(jù)記為DATA2。在下一步中DATA2將進(jìn)行傅里葉變換，數(shù)據(jù)從時(shí)域變換到頻域，頻域中的數(shù)據(jù)記錄為DATA3，DATA3再進(jìn)行自相關(guān)后即為閃爍譜。圖5是密云50m望遠(yuǎn)鏡觀測(cè)得到的一個(gè)閃爍譜。

圖5 密云50米射電望遠(yuǎn)鏡觀測(cè)到3C48的閃爍譜;2008年3月27日；S波段Fig.5 IPS Spectrum observed with the Miyun 50m radio telescope through the S band on March 27,2008

圖6 2005年9月Ooty望遠(yuǎn)鏡觀測(cè)3C279數(shù)據(jù)原始圖Fig.6 Raw data toward 3C279 observed by the Ooty radio telescope in September 2005

密云現(xiàn)有的輻射計(jì)積分時(shí)間過長(zhǎng)且增益可調(diào)范圍小，且IPS現(xiàn)象一般在較低頻率觀測(cè)明顯，目前不適合進(jìn)行行星際閃爍觀測(cè)。在密云IPS專用輻射計(jì)研制期間，我們已開始申請(qǐng)烏魯木齊25m望遠(yuǎn)鏡開展試驗(yàn)觀測(cè)研究，取得了初步成功(在另文討論)。

2.2 數(shù)據(jù)處理驗(yàn)證

為了對(duì)自主編寫的數(shù)據(jù)處理軟件進(jìn)行驗(yàn)證，2009年3月，聯(lián)系了印度IPS組的Prof Monoharan P K,請(qǐng)他提供一些Ooty實(shí)測(cè)數(shù)據(jù)和處理結(jié)果，以檢驗(yàn)處理軟件。Prof Monoharan提供了2005～2007年間印度Ooty射電望遠(yuǎn)鏡觀測(cè)到的行星際閃爍原始數(shù)據(jù)及對(duì)應(yīng)的閃爍譜圖。

使用我們的軟件對(duì)多組印度Ooty望遠(yuǎn)鏡的觀測(cè)數(shù)據(jù)進(jìn)行了處理。處理結(jié)果與Ooty提供的處理結(jié)果符合很好，證明我們的數(shù)據(jù)處理達(dá)到了預(yù)期目標(biāo)。

圖6～8為處理的Ooty觀測(cè)數(shù)據(jù)對(duì)比圖。圖中觀測(cè)源為3C279，觀測(cè)時(shí)間是2005年9月2日，觀測(cè)頻率為327MHz。其中圖6為Ooty望遠(yuǎn)鏡所觀測(cè)的原始數(shù)據(jù)；圖7為該數(shù)據(jù)印度提供的線性坐標(biāo)閃爍譜圖與利用自主編寫的數(shù)據(jù)處理軟件得到的線性坐標(biāo)閃爍譜圖對(duì)比；圖8為利用自主編寫的數(shù)據(jù)處理軟件得到的對(duì)數(shù)坐標(biāo)閃爍譜圖及擬合圖，圖中實(shí)線為觀測(cè)所得閃爍譜圖，虛線為擬合圖。

(a) (b)

圖7 2005年9月Ooty望遠(yuǎn)鏡觀測(cè)3C279的IPS閃爍譜
(a)印度提供的閃爍譜；(b)自編寫數(shù)據(jù)處理軟件閃爍譜
Fig.7 IPS spectrum toward 3C279 observed by the Ooty radio telescope，(a) provided by the Indian group,(b) reduced by our own software

圖8 自編數(shù)據(jù)處理軟件的Ooty望遠(yuǎn)鏡數(shù)據(jù)閃爍譜及參數(shù)擬合結(jié)果(橫坐標(biāo)取對(duì)數(shù))觀測(cè)時(shí)間：2005年9月。Fig.8 IPS Spectrum from the data obtained with the Ooty radio telescope in September,2005 and reduced by our own software;x axis:frequenuy;y axis：logarithm of spectral power

由圖7可見，印度和我們的兩種數(shù)據(jù)處理方法得出的結(jié)果符合得很好。

對(duì)2005年Ooty望遠(yuǎn)鏡的觀測(cè)結(jié)果進(jìn)行了參數(shù)擬合，模型模擬結(jié)果如圖8。模型采用的參數(shù)為：AR=2.5，α=3.5，V=600km/s。由fF=1.1Hz和公式(1)得到V為672.4km/s,兩者在誤差允許范圍(～10%)內(nèi)吻合很好。 spaceweather.com公布的當(dāng)天太陽風(fēng)速為672.4km/s,與觀測(cè)的結(jié)果基本一致。

通過對(duì)Ooty射電望遠(yuǎn)鏡IPS觀測(cè)的數(shù)據(jù)處理，以及與spaceweather.com相應(yīng)時(shí)間觀測(cè)結(jié)果比較，可以得出以下結(jié)論：數(shù)據(jù)處理方法是正確的，數(shù)據(jù)處理結(jié)果也達(dá)到預(yù)期。

3 結(jié) 論

編寫的數(shù)據(jù)處理軟件適用于行星際閃爍單站單頻模式的數(shù)據(jù)處理。當(dāng)進(jìn)行單站雙頻模式觀測(cè)時(shí)，可以對(duì)此軟件進(jìn)行相應(yīng)的拓展：觀測(cè)得到的數(shù)據(jù)先分別進(jìn)行以上處理流程，在傅里葉變換后，得到的兩組數(shù)據(jù)不僅要做自相關(guān)(記為DATA11，DATA22)，同時(shí)還要進(jìn)行互相關(guān)運(yùn)算(記為DATA12，DATA21)。得到的4組數(shù)據(jù)再進(jìn)行減平均值的計(jì)算就可以得到歸一化互相關(guān)譜。通過計(jì)算歸一化互相關(guān)譜第一零點(diǎn)頻率可估算太陽風(fēng)風(fēng)速。

行星際閃爍的單站單頻模式，無論從觀測(cè)、數(shù)據(jù)采集設(shè)備，還是數(shù)據(jù)流程來講都較簡(jiǎn)單，易于操作，被國(guó)際上多數(shù)觀測(cè)站采用。但是，這種觀測(cè)模式要求有較高信噪比(大于25dB)和靈敏度。在觀測(cè)中幾乎所有觀測(cè)到的閃爍都可以找到fF，這是Fresnel過濾因子與電子密度冪函數(shù)譜相互作用的結(jié)果。進(jìn)行模式擬合時(shí)，由fF可知選擇合適的太陽風(fēng)速度值。另外，由于路徑積分效應(yīng)，F(xiàn)resnel振蕩現(xiàn)象已不明顯，極小值退變成平臺(tái)，只能觀測(cè)到第一極小值fmin，而其它的極小值總是被淹沒。當(dāng)各向異性軸比AR增大時(shí)，fF就會(huì)變得不明顯，從而對(duì)估測(cè)速度造成困難。同時(shí)fmin也容易受到噪聲、AR和ε的影響，因此參數(shù)擬合的準(zhǔn)確性容易受到所選參數(shù)的影響而存在一定的不確定性[12]。

相對(duì)于單站單頻模式，單站雙頻模式具有以下幾個(gè)優(yōu)點(diǎn)：更高的靈敏度，更精確的測(cè)量太陽風(fēng)速，對(duì)于各種太陽風(fēng)參數(shù)具有更高的穩(wěn)定性[6]。但是這種模式對(duì)觀測(cè)系統(tǒng)和接收系統(tǒng)有較高要求，同時(shí)數(shù)據(jù)處理流程也較單站單頻模式復(fù)雜。

目前國(guó)家天文臺(tái)密云基地50m射電望遠(yuǎn)鏡開展的單站雙頻模式的接收機(jī)系統(tǒng)已處于研制過程中。

致謝：感謝印度Prof Monoharan P K 提供Ooty望遠(yuǎn)鏡觀測(cè)數(shù)據(jù)，特別感謝張喜鎮(zhèn)老師的有益討論和指導(dǎo)。

[1] Hewish A,Symonds M D.Radio investigation of the solar plasma.Planetary and Space Science[J].Nature,1969,17:313.

[2] Armstrong J W,Coles W A.Analysis of three-station interplanetary scintillation[J].Geophysics Res,1972,77:4602-4610.

[3] Watanabe T,Shibasaki K,Kakinuma T.Latitudinal distribution of solar wind velocity and its relation to solar EUV corona[J].Geophys Res,1974,79:3841-3843.

[4] Hewish A,Scott P F,Wills D.Radio investigation of the solar plasma[J].Nature,1964,203:1214.

[5] Purvis A,Tappin S J,Rees W G,et al.The Cambridge IPS survey at 81.5 MHz[J]. Royal Astronomical Society,Monthly Notices,1987,229:589-619.

[6] Vitkevich V V，Glushaev A A,Iliasov Iu P,et al.Antenna equipment of the Lebedev Institute BSA radio telescope facility[J].Radiofizika，1976，19:1594-1606.

[7] Kojima M，Asai，Kozuka Y，et al.Velocity observations at high latitude and the acceleration phenomena[J].Advances in Space Research，1995,16(9):101-110.

[8] Swarup G，Sarma N V G，Joshi M N,et al.Large Steerable Radio Telescope at Ootacamund,India[J].Nature Physical Science,1971,230:185-188.

[9] Wu J H,Zhang X Z,Zheng Y J.IPS Observations at Miyun Station,Bao[J].ApSS,2001，278:189-192.

[10] Scott L,Rickett B J,Armstrong J W.The velocity and the density spectrum of the solar wind from simultaneous three-frequency IPS observations[J].Astronomy and Astrophysics,1983,123:191-206.

[11] Zhang X Z.A Study on the Technique of Observing Interplanetary Scintillation with Simultaneous Dual-frequency Measurement[J].Chin J Astron & Astrophys,2007, 7(5):712-720.

[12] 葉品中，邱育海.行星際閃爍(IPS)單站觀測(cè)對(duì)太陽風(fēng)速度的診斷[J].天體物理學(xué)報(bào)1996,16(4):389-394.