2. 中国科学院大学, 北京 100049;
3. 中国科学院天体结构与演化重点实验室, 云南 昆明 650216;
4. 云南省太阳物理与空间目标监测重点实验室, 云南 昆明 650216
2. University of Chinese Academy of Sciences, Beijing 100049, China;
3. Key Laboratory for the Structure and Evolution of Celestial Objects, Chinese Academy of Sciences, Kunming 650216, China;
4. Yunnan Key Laboratory of the Solar Physics and Space Science, Kunming 650216, China
太阳耀斑、日冕物质抛射(Coronal Mass Ejection, CME) 等太阳活动是太阳系中最剧烈的能量释放过程,与空间天气密切相关。该过程中产生的射电辐射会携带着日冕等离子体和磁场的信息[1],产生的高能粒子和宇宙射线可能对空间环境和地磁场产生扰动,直接影响到卫星、通讯、电力等现代化技术系统等,进而影响人类的生活环境[2-4]。近些年,随着我国航天事业的快速发展,空间天气逐渐成为太阳物理领域的研究热点。
太阳射电动态频谱是观测太阳射电爆发的重要手段,文[5] 根据太阳动态谱上的爆发形态,将其分为5类(Ⅰ, Ⅱ, Ⅲ, Ⅳ和Ⅴ型暴)。随着谱分辨率和时间分辨率的提高,研究人员发现了各种各样的精细结构,比如“斑马纹”结构、新型“曲棍球棒”结构等[6-9]。高分辨率动态频谱对太阳射电爆发精细结构的研究具有重要作用,有助于我们理解太阳射电爆发过程中的物理机制[10]。太阳射电动态频谱不仅可以提供关于日冕磁场和等离子体的信息[11],如等离子体辐射可以提供等离子体密度信息,回旋同步辐射可以提供磁场信息,而且可以用于空间天气的预警预报,如Ⅱ型暴是日冕物质抛射激波的最佳指示器,根据Ⅱ型暴的漂移率可以估算激波的速度。因此,太阳射电动态频谱不仅对太阳物理的研究具有重要作用,而且对空间天气的预警预报具有重要作用[12-13]。
子午工程二期太阳-行星际监测链分系统即将建设完成的太阳射电频谱仪,其频率覆盖30 MHz~15 GHz,站址位于中国东部(来源于https://www.sjdz.org.cn/cn/article/doi/10.19975/j.dqyxx.2023-016)。本项目中的太阳射电频谱仪放置于中国西部(四川稻城),两站配合观测,可以延长对太阳射电爆发活动连续监测的时段,加强我国对空间天气的实时预报能力。
本文设计的超宽带超高分辨率的太阳射电频谱仪,共由4套设备组成。数字频谱仪是基于现场可编程门阵列研发的数字终端,是设备的关键部分,充分利用数字信号处理的灵活性,采用快速傅里叶变换算法完成频谱分析。基于超高谱分辨率的需求,需要的快速傅里叶变换点数为262144,在现场可编程门阵列上通过一个FFT IP不能实现如此高点数的快速傅里叶变换运算。为此,我们根据Cooley-Tukey FFT算法,首先对数据分通道,然后做并行快速傅里叶变换处理,最后通过数据重组完成所需点数的快速傅里叶变换运算[14]。
1 太阳射电频谱仪的组成太阳射电频谱仪的组成如图 1,主要包括天线系统、接收机和数字频谱仪。天线系统由反射面和馈源组成,用于接收太阳辐射的射电信号并将其传输到接收机;接收机对信号进行放大、变频和滤波处理;数字频谱仪实现信号的采样和频谱分析,将频谱数据传输到计算机并保存。
![]() |
图 1 太阳射电频谱仪的组成框图 Fig. 1 The block diagram of solar radio spectrometer |
从厘米波段到十米波段,天线系统分别为6 m抛物面和四脊喇叭馈源,8 m抛物面和四脊喇叭馈源,12 m抛物面和对数周期馈源,以及对数周期天线,均采用赤道式架构,可以自动计算太阳的轨道运动参数并精确跟踪。对于6 m口径的抛物面天线,根据波束宽度的经验公式,半功率波束宽度(Half Power Beam Width, HPBW)= 1.22 λ/D (λ为射电信号的波长,D为天线口径)。7 GHz以上的频率范围,波束宽度小于0.5°,无法覆盖全部日面。通过采用波束赋形技术,实现7 GHz以上频率的半功率波束宽度大于0.5°,仿真结果如图 2。图中(a), (b), (c)和(d) 对应的频点依次为6 GHz, 9 GHz, 12 GHz和15 GHz,从图中可见半功率波束宽度依次为0.78°,0.54°,0.50°和0.58°,满足半功率波束宽度覆盖全部日面的需求。
![]() |
图 2 6 m抛物面天线的辐射方向图。(a), (b), (c)和(d) 对应的频点依次为6 GHz, 9 GHz, 12 GHz和15 GHz Fig. 2 The radiation pattern of parabolic antenna with diameter of 6 m. (a), (b), (c) and (d) correspond to the frequency of 6 GHz, 9 GHz, 12 GHz and 15 GHz, respectively |
接收机分别包括厘米波接收机、分米波接收机、米波接收机和十米波接收机。接收机对信号进行低噪声放大、滤波、功率分配、变频、增益控制等处理。分米波段和厘米波段接收机的频率范围很宽,采用滤波器组划分频率范围再分段变频的形式,输出多路中频通道。如果中频通道的带宽太窄,则中频通道数量会很多,增加接收机的复杂度;如果中频通道太宽,不仅对接收机的器件要求苛刻,也会对数字频谱仪提出更高的要求。综合考虑数字频谱仪的采样频率以及现场可编程门阵列时钟速度等因素,将中频通道的带宽设计为1 GHz,输出频率范围为0.1~1.1 GHz。
分米波接收机和厘米波接收机分别有4个和11个中频通道,米波接收机和十米波接收机分别有1个通道。太阳射电频谱仪需要同时观测两个极化信号,共需要34个采样通道。数字频谱仪的硬件成本很高,常规的做法是在各个中频通道之间加入微波开关,通过微波开关的切换实现各中频通道的轮流采样。此方法在降低数字频谱仪硬件成本的同时降低了信号的积分时间,导致系统灵敏度降低。因此,我们设计的厘米波段和分米波段接收机采用左右旋的切换方式,可以减少一半的数字频谱仪采样通道,但没有明显降低积分时间。对米波段和十米波段,带宽比较窄,可以一次完成采样,因此左右旋通道同时采样,分别需要2个采样通道,共需要19个采样通道。各个波段接收机的设计框图如图 3。
![]() |
图 3 各个波段接收机的设计结构。(a), (b), (c)和(d) 依次对应厘米波段、分米波段、米波段和十米波段接收机 Fig. 3 The design structure of analog receivers for each waveband. (a), (b), (c) and (d) correspond to centimetric receiver, decimetric receiver, metric receiver and decametric receiver |
接收机主要包括定标系统、低噪声放大器(Low Noise Amplifier, LNA)、极化合成器和模拟接收机。(1) 定标系统包括微波开关、噪声源和50 Ω负载,通过微波开关在天线信号、噪声源和50 Ω负载之间切换,完成对接收机系统的定标;(2) 低噪声放大器实现信号的低噪声放大;(3) 圆极化转换器实现线极化信号到圆极化信号的转换;(4) 模拟接收机包括放大器、光收发组件、后级放大、数控衰减器、功分器、滤波器和混频组件,实现信号放大、室外到室内的传输以及信号变频等,最终输出所需的中频信号。除上述设计外,在混频组件中配备数控衰减器,可以实现增益在30 dB范围内,步进为1 dB的调节,从而增加接收机的动态范围。4套接收机研发完成后,测试指标如表 1。
Gain/dB | Gain flatness/dB | Noise figure/dB | Input VSWR | Output VSWR | Harmonic suppression/dBc | In-band emission/dBc | Out-band emission/dBc | Image rejection/dBc | Pout, 1 dB/dBm | Phase noise of LO | |
Centimetric receiver | 62.98 | 0.84 | 2.46/ 3.61 | 1.47 | 1.49 | 82.47 | 56.12 | 64.10 | 56.49 | 19.18 | -91.0 dBc/Hz @10 kHz |
Decimetric receiver | 66.32 | 1.20 | 1.96 | 1.73 | 1.49 | 62.06 | 58.87 | 65.72 | 66.07 | 16.32 | -93.8 dBc/Hz @10 kHz |
Metric receiver | 65.69 | 1.43 | 1.92 | 1.85 | 1.41 | 53.01 | 54.28 | — | — | 16.44 | — |
Decametric receiver | 65.00 | 0.70 | 2.15 | 1.16 | 1.42 | 48.70 | 63.07 | — | — | 17.26 | — |
注释:厘米波段的噪声系数在4~8.5 GHz范围内最大值为2.46 dB,在8.5~15 GHz范围内最大值为3.61 dB。 |
数字频谱仪是通过快速傅里叶变换将信号从时域转换为频域再求功率谱的仪器。信号处理流程如图 4:(1) 输入的模拟信号由模数转换器(Analog-to-Digital Converter, ADC) 采样,实现信号的数字化;(2) 快速傅里叶变换运算实现信号从时域到频域的转换,设计中的快速傅里叶变换点数很大,需要使用并行快速傅里叶变换;(3) 快速傅里叶变换运算结果的模值作为信号能量的估计;(4) 通过能量的叠加完成积分,可以提高信噪比;(5) 除以积分时间作为功率谱密度的估计;(6) 将频谱数据通过光纤传输到计算机并保存。
![]() |
图 4 数字频谱仪的信号处理流程 Fig. 4 The signal processing flow of digital spectrometer |
根据奈奎斯特采样定律,采样率必须大于输入信号最高频率的2倍,所以将厘米波段、分米波段、米波段数字频谱仪的采样率设计为2.5 Gsps,十米波段数字频谱仪的采样率为250 Msps。根据Δf=Fs/N,其中Δf为频率分辨率(谱分辨率),Fs为采样频率,N为快速傅里叶变换运算的点数,表 2给出各个波段数字频谱仪的主要技术指标。
Sample | Spectral resolution/kHz | Temporal resolution/ms | Word length of FFT | |
Centimetric digital spectrometer | 2.5 Gsps | 76 | 10 | 32 768 |
Decimetric digital spectrometer | 2.5 Gsps | 95 | 10 | 262 144 |
Metric digital spectrometer | 2.5 Gsps | 95 | 10 | 262 144 |
Decametric digital spectrometer | 250 Msps | 7.6 | 1 | 32 768 |
数字频谱仪是基于标准VPX-6U架构的数据采集处理平台,由VPX-6U机箱、载板(母板)、采集子板和通信子板组成的高速信号处理系统。采用VPX高速背板设计,可以完成载板与载板的高速互联,载板和子板之间采用FMC HPC模式链接。子板可以在任意载板上互相交换,从而保证板卡的一致性和可交换性。数字频谱仪实物结构如图 5。
![]() |
图 5 数字频谱仪的实物结构图 Fig. 5 The physical structure of digital spectrometer |
载板是基于现场可编程门阵列设计开发的数字信号处理板,搭载两块现场可编程门阵列芯片用于处理两个通道的数据,8片DDR3 SDRAM存储模块用于数据的缓存。现场可编程门阵列易于实现流水和并行结构,并且具备灵活性等特点。FGPA型号为JFM7Vx690T36,通过对算法所需资源的仿真,其资源(存储资源、计算资源等) 满足需求。ADC的型号为ADI AD9625,采样率为2.5 Gsps,采样位宽为12 bit,无杂散动态范围为71 dBfs,有效位优于8.5 bit,满足动态范围的需求。通信子卡不仅支持10/100 M以太网口,还支持光纤通信,传输速率可以达到10 Gbps,满足频谱数据的传输速度,保证数据的实时传输。
2.2 并行快速傅里叶变换算法数字频谱仪的核心算法是快速傅里叶变换的实时运算,在现场可编程门阵列内部可以直接调用快速傅里叶变换IP核完成变换运算,但是快速傅里叶变换IP核能完成变换点数的上限无法满足项目中大点数快速傅里叶变换运算的需求。同时,由于项目中的模数转换器采样率为2.5 Gsps,现场可编程门阵列的时钟速度达不到如此高的速率,因此需要采用并行快速傅里叶变换算法以便满足高速且大点数的快速傅里叶变换运算。比如,数字频谱仪所需的262 144点的快速傅里叶变换,可以分成8路并行快速傅里叶变换进行计算,每路快速傅里叶变换点数为32 768,每路的工作时钟为312.5 MHz,现场可编程门阵列的时钟满足需求。
根据Cooley-Tukey FFT算法,对N点序列 x[n] 做傅里叶变换,定义为
$ \begin{equation} X[k]=\sum\limits_{n=0}^{N-1} x[n] W_{\mathrm{N}}^{n k}, \end{equation} $ | (1) |
其中,k=0, 1, 2, …, N-1。为了实现并行运算,可以将 x[n] 的N点数据排列成M行L列的矩阵,其中N=M×L,则 x[n] 可以转换成一个二维的快速傅里叶变换实现,
$ \begin{equation} X[k]=\sum\limits_{m=0}^{M-1}\left\{\left[\sum\limits_{l=0}^{L-1} x[l, m] W_{\mathrm{L}}^{k_0 l}\right] W_{\mathrm{N}}^{m k_0}\right\} W_{\mathrm{M}}^{m k_1}, \end{equation} $ | (2) |
令
$ \begin{equation} Y\left[m, k_0\right]=\sum\limits_{l=0}^{L-1} x[l, m] W_{\mathrm{L}}^{k_0 l}, \end{equation} $ | (3) |
$ \begin{equation} Z\left[m, k_0\right]=Y\left[m, k_0\right] W_{\mathrm{N}}^{m k_0}, \end{equation} $ | (4) |
则
$ \begin{equation} X[k]=\sum\limits_{m=0}^{M-1} Z\left[m, k_0\right] W_{\mathrm{M}}^{m k_1}, \end{equation} $ | (5) |
其中,0≤l,k0≤L-1,0≤m,k1≤M-1,WNmk0为旋转因子。(3) 式是L点的快速傅里叶变换运算;(4) 式是运算结果乘以旋转因子;(5) 式是M点快速傅里叶变换运算。一维N点的快速傅里叶变换分解成二维M×L点的快速傅里叶变换,信号处理流程如图 6。首先对矩阵的每一行进行快速傅里叶变换,然后乘以旋转因子,再对每一列进行快速傅里叶变换,最后经过数据重组完成快速傅里叶变换。
![]() |
图 6 并行快速傅里叶变换算法 Fig. 6 Parallel FFT algorithm |
根据上述信号处理流程,对并行快速傅里叶变换算法进行仿真,输入一个频率为762.94 MHz的正弦信号(包含一个高斯噪声) 的采样数据流,采样频率为2.5 Gsps,快速傅里叶变换点数选为32 768。使用并行快速傅里叶变换算法(M=8,L=4 096) 的仿真结果如图 7 (a),并行算法在现场可编程门阵列上实现的结果如图 7 (b)。对比两图中的信号都是762.94 MHz,没有显著不同,表明并行快速傅里叶变换算法在现场可编程门阵列上可以成功运用。对262 144点的快速傅里叶变换也做同样的处理,得到相似的结果。
![]() |
图 7 32 768点快速傅里叶变换的仿真结果。(a) 在Matlab中使用并行算法的32 768点快速傅里叶变换;(b) 从现场可编程门阵列中实施并行算法的32 768点快速傅里叶变换的结果 Fig. 7 The simulation results of FFT with word length of 32 768. (a) The result of parallel FFT with word length of 32 768 in Matlab; (b) the result of parallel FFT with word length of 32 768 implemented in FPGA |
此4套太阳射电频谱仪可以覆盖15 MHz~15 GHz的频谱范围,用于监测从日冕底部到两个太阳半径范围内的太阳射电爆发,这里是太阳射电爆发初始能量释放、日冕物质抛射、激波形成和加速的主要区域。同时,此太阳射电频谱仪具备超高谱分辨率和时间分辨率,达到甚至超过世界同类设备的水平,不仅为太阳射电爆发的研究提供有效的数据,而且为空间天气的预报预警提供可靠的数据。目前除天线正在加工,接收机和数字频谱仪已经研发完成。数字频谱仪的测试数据表明,采用并行快速傅里叶变换方式可以实现大点数快速傅里叶变换,获得高分辨率的频谱。结合国内的日像仪,比如明安图厘米-分米波射电频谱日像仪和稻城圆环阵太阳射电成像望远镜等设备,可以成功完成谱-像结合方式对太阳射电的观测研究。
[1] | BASTIAN T, BAIN H, CHEN B, et al. Diagnostics of space weather drivers enabled by radio observations[J]. Bulletin of the American Astronomical Society, 2019, 51(3): 323–330. |
[2] |
施硕彪, 董亮, 高冠男, 等. 米波太阳射电频谱仪的科学目标和技术方案[J]. 天文研究与技术, 2011, 8(3): 229–235 SHI S B, DONG L, GAO G N, et al. Scientific objectives and technical design of a meter-wave spectrometer for solar radio observation[J]. Astronomical Research & Technology, 2011, 8(3): 229–235. |
[3] |
颜毅华. 中国科学院国家天文台太阳物理研究20年[J]. 科学通报, 2021, 66(11): 1363–1384 YAN Y H. Research advances in solar physics at National Astronomical Observatories of Chinese Academy of Sciences[J]. Chinese Science Bulletin, 2021, 66(11): 1363–1384. |
[4] | LE G M, CAI Z, WANG H N, et al. Solar cycle distribution of great geomagnetic storms[J]. Astrophysics and Space Science, 2012, 339(1): 151–156. DOI: 10.1007/s10509-011-0960-y |
[5] | WILD J P. The radio emission of the sun[M]//PALMER H P, LARGE M I. Radio astronomy today. Manchester: Manchester University Press. 1963: 1-23. |
[6] | BASTIAN T S, BENZ A O, GRAY D E. Radio emission from solar flare[J]. Annual Review of Astronomy and Astrophysics, 1998, 36: 131–188. DOI: 10.1146/annurev.astro.36.1.131 |
[7] | FU Q J, YAN Y H, LIU Y Y, et al. A new catalogue of fine structures superimposed on solar microwave bursts[J]. Chinese Journal of Astronomy and Astrophysics, 2004, 4(2): 176–188. DOI: 10.1088/1009-9271/4/2/176 |
[8] | CHEN B, YAN Y H. On the origin of the zebra pattern with pulsating superfine structures on 21 April 2002[J]. Solar Physics, 2007, 246(2): 431–443. DOI: 10.1007/s11207-007-9063-x |
[9] | HUANG J, YAN Y H, LIU Y Y. An analysis of solar radio burst events on December 1, 2004[J]. Advances in Space Research, 2007, 39(9): 1439–1444. DOI: 10.1016/j.asr.2007.03.061 |
[10] |
谭宝林, 程俊, 谭程明, 等. 尖峰爆发标度律及其对新一代太阳射电望远镜参数的约束[J]. 天文学报, 2018, 59(4): 64–76 TAN B L, CHENG J, TAN C M, et al. Scaling-laws of radio spike bursts and their constraints on new solar radio telescopes[J]. Acta Astronomica Sinica, 2018, 59(4): 64–76. |
[11] | ZHANG Y Y, ZHANG L, SHANG Z Q, et al. A new multichannel parallel real-time FFT algorithm for a solar radio observation system based on FPGA[J]. Publications of the Astronomical Society of the Pacific, 2022, 134(1033): 034502. DOI: 10.1088/1538-3873/ac5212 |
[12] |
高冠男, 汪敏, 董亮, 等. 空间甚低频太阳Ⅱ型射电暴研究进展[J]. 深空探测学报, 2021, 8(4): 423–432 GAO G N, WANG M, DONG L, et al. Advances in space VLF type Ⅱ solar radio bursts[J]. Journal of Deep Space Exploration, 2021, 8(4): 423–432. |
[13] | WINTER L M, LEDBETTER K. Type Ⅱ and Type Ⅲ radio bursts and their correlation with solar energetic proton events[J]. The Astrophysical Journal, 2015, 809(1): 105. DOI: 10.1088/0004-637X/809/1/105 |
[14] | COOLEY J W, TUKEY J W. An algorithm for the machine calculation of complex Fourier series[J]. Mathematics of Computation, 1965, 19: 297–301. DOI: 10.1090/S0025-5718-1965-0178586-1 |