高光谱微波辐射计系统中2GHz带宽数字谱仪设计

许皓文，陆 浩，王振占

（1. 中国科学院国家空间科学中心微波遥感技术重点实验室，北京 100190；2. 中国科学院大学，北京 100049）

1 引言

微波临边探测是一种新型的观测地球大气的方式，主要通过天线从天底往上垂直扫描大气临边切点时测量大气毫米波、亚毫米波的热发射，得到各种大气变量的遥测值［1］. 在微波临边探测中，高光谱微波辐射计作为一种被动微波遥感器件，主要用来探测大气的痕量气体成分［2］，而数字谱仪是高光谱微波辐射计的核心部件，主要承担了后端的谱分析功能［3，4］.

在目前已经发射成功并获取到观测数据的临边探测卫星中，最早搭载微波临边探测器（Microwave Limb Sounder，MLS）的是美国于1991年发射的高层大气研究卫星（Upper Atmosphere Research Satellite，UARS），其搭载的是滤波器组谱仪（Filter Banks Spectrometer，FBS）［5，6］.2001 年，瑞典航天公司联合其他几个国家发射了高层大气物理研究和天文学联合任务卫星（Odin），其搭载的是数字自相关谱仪（Digital Auto-Correlator Spectrometer，DACS）和声光谱仪（Acousto-Optical Spectrometer，AOS）［7，8］. 美国对地观测系统（Earth Observing System，EOS）中的Aura 卫星于2004 年发射，是UARS的继任者，其搭载的是FBS和DACS［9，10］.日本在2009年底发射了超导亚毫米波临边辐射探测仪（JEM/SMILES），其后端是由AOS组成［11，12］.

随着数字信号处理技术和半导体技术的快速发展，数字谱仪也逐渐成熟并展示出了其特有的优势. 数字谱仪主要包括了DACS 和快速傅里叶变换谱仪（Fast Fourier Transform Spectrometer，FFTS），主要特点是利用超大规模集成电路实现宽频带信号的并行处理，而且还具有抗干扰能力强、稳定性高、体积小和更加灵活等优点［13］.

考虑到卫星上搭载时的限制因素，以及系统对输入带宽和谱分辨率的要求等，本文提出了一种新型宽带、实时快速傅里叶变换类型的数字谱仪设计方案. 与DACS相比，FFTS在实现方式上具有更大的优势，可以降低系统复杂度，节约硬件资源. 本文所设计的FFTS分析带宽2 GHz，通道数为1 024（1k），与各个卫星上搭载谱仪的参数指标对比如表1所示. 由于各类谱仪体制不同，AOS以像素点进行衡量，其余类型均以通道数进行衡量.

从表1可以看出，本文所设计的数字谱仪在分析带宽和通道数上均有所提高. 在实际设计中一般要求分析带宽尽可能大，通道数尽可能多，但是从数字谱仪研制的角度来看，高带宽和多通道是很难同时兼顾的，所以在微波临边探测中，数字谱仪要充分平衡分析带宽和通道数之间的关系，针对不同的探测需求对谱分辨率进行动态的调整. 例如在大气的中高层区域，降低谱分辨率会在一定程度上降低反演的精度，过低的谱分辨率无法辨别大气上层的气体吸收峰，所以在保持足够灵敏度的前提下，2 MHz的谱分辨率有助于提高中高层大气的反演精度［14］.

表1 谱仪参数指标对比

2 系统原理及核心算法设计

2.1 高光谱微波辐射计系统原理

高光谱微波辐射计是一种能够获取大量窄带连续频谱通道的微波遥感系统，其主要原理是测量天线的辐射温度，之后对射频前端接收机产生的中频模拟信号进行数字化处理，通过后端数字谱仪实现输入信号谱分析的功能.

高光谱微波辐射计系统框图如图1所示，系统中主要集成了天线、射频前端接收机以及后端数字谱仪. 接收机又包括了射频放大器、混频器和中频放大器等，其主要作用是对中心频率包含在带宽内的输入信号进行放大及带外抑制，输出同样带宽的中频信号［15］. 数字谱仪主要是利用集成的模拟数字转换器（Analog to Digital Converter，ADC）和现场可编程门阵列（Field Program⁃mable Gate Array，FPGA）等芯片，通过模数转换，快速傅里叶变换（Fast Fourier Transform，FFT），从而得到功率谱密度（Power Spectral Density，PSD）的结果. 通过对数字谱仪获取的大气中痕量气体的功率谱进行定标分析，可以得到大气成分的辐射亮温谱分布，从而反演大气参数和分布变化，最终得到大气成分的分布规律，以便对其进行合理的科学研究.

图1 高光谱微波辐射计系统结构

2.2 核心算法设计

基16实时复数FFT 算法是实现数字谱仪功能的核心算法，其是在离散傅里叶变换（Discrete Fourier Trans⁃form，DFT）的基础上经过并行处理［16］和复数处理［17］的改进和结合而得到的. 并行处理可以使FFT 的计算更加快速，同时降低计算复杂度，提高运算效率；而复数处理更加适用于硬件的实现过程，可以节约硬件资源，提高运算速度.

2.2.1 并行处理

对于输入的离散序列，可以利用DFT 的快速运算将原本的串行计算转化为并行计算，通过数据分解、矩阵映射、一次计算、对应相位因子乘积和二次计算就可以完成并行FFT的处理.

假定离散序列x（n）包含N点数据，如果N可以分解为2 个整数R和T的乘积，如式（1）所示，且R和T均为2的整数次幂的形式，那么对x（n）做DFT可以得到长度为N的复序列X（k），如式（2）所示.

其中

如果将x（n）映射为矩形数组x（r，t），将DFT的结果X（k）映射为相应的矩形数组X（p，q），在此过程中，选择对x（n）采用式（4）的按列映射的方式，而对X（k）采用式（5）的按行映射的方式.

因此可以得到

其中

化简后可得

可以将式（8）分解为以下3 步，如式（9）～（11）所示.

（1）计算T点的DFT，即

对矩形数组x（r，t）按行计算有r=0，1，…，R-1.

（2）定义新的矩形数组G（r，q），即

（3）计算R点的DFT，即

对矩形数组G（r，q）按列计算有q=0，1，…，T-1.

2.2.2 复数处理

对于输入的实数序列，可以利用DFT 的对称性将2N点实数序列的计算转化为N点复数序列的计算，之后将输出的N点复数序列进行适当的排序，就可以得到原始的2N点实数序列FFT的结果.

假定x（n）为一个2N点的实数序列，令x1（n）代表该序列中的偶数分量，x2（n）代表该序列中的奇数分量，则复序列y（n）可被定义为

由于DFT是线性的，所以y（n）的DFT可以表示为

其中，YRe（k）和YIm（k）分别代表了Y（k）的实部和虚部，序列x1（n）和x2（n）可以用y（n）表示为

其中，y*为y的共轭对称，由x1（n）和x2（n）的DFT结果可得

将式（13）分别代入式（16）和式（17），可得

根据式（18）和式（19），可以定义4 组中间因子YRe1（k），YRe2（k），YIm1（k），YIm2（k），如式（20）所示：

对于2N点的实数序列x（n），对其进行FFT 还可以表示为

其中

综合以上公式，可以得到

根据式（23）可以得到2N点实数序列经过FFT运算后前一半的结果，后一半与前一半对称分布.

3 数字谱仪设计及实现

3.1 整体设计流程

数字谱仪主要通过ADC 和FPGA 等芯片，利用核心的基16 实时复数FFT 算法完成相应的宽带信号处理，累加后从而得到功率谱的结果.2 GHz 带宽数字谱仪整体数据流程如图2 所示，4.8 GHz 高速数据在ADC内部经过1∶8的多路复用（DMUX）之后［18］，输出的数据和时钟为低电压差分信号（Low Voltage Differential Sig⁃nal，LVDS），输入至FPGA 的LVDS 数据和时钟首先通过Buffer 转换成单端数据和时钟，随后数据进入FPGA的1∶2 DMUX 模块，之后通过核心算法进行处理，经过一系列的拆分、整合和计算后可以得到每一包数据的结果，累加后通过串口上传至上位机.

图2 数字谱仪整体数据流程

3.2 核心算法实现

在本文数字谱仪核心算法的实现过程中，根据式（9）～（11），将T和R的值分别设定为64和16. 数字谱仪核心算法架构如图3 所示，实际架构包括了16 路并行数据的处理，这里仅以其中2 路作为代表进行说明.核心算法架构主要包括了64 点FFT 模块、旋转因子存储模块、数据相乘模块、基16-FFT模块、数据排序模块、中间因子模块、数据计算模块、数据累加模块和并转串模块［17］.

图3 数字谱仪核心算法架构

64 点FFT 模块首先利用FFT IP 核同时完成16 路，每一路64 点的FFT，IP 核设置为流水线的方式可以保证数据的实时性. 旋转因子存储模块是通过FPGA ROM 存储复数旋转因子，为了节约资源，本设计中将复数旋转因子的实部和虚部合并为一个数据进行存储.数据相乘模块对前面的结果进行复数的乘法运算. 中间因子模块和数据计算模块分别按照式（20）和式（23）进行相应的计算，累加后通过并转串最终可以得到1 024个复数点的结果.

本文在算法实现中为了达到实时处理的目的，特别对基16-FFT 模块和数据排序模块进行了改进和优化. 在基16-FFT 模块，根据式（11），采取了最直接的矩阵相乘的运算方式. 在数据排序模块，采取了“乒乓”RAM 的数据处理方式，从而保证了整个设计过程中实时计算的持续性.

4 仿真结果分析

通过模拟正弦信号以及真实亮温信号可以对核心算法完成前期的设计仿真以及验证，同时根据真实亮温信号仿真得到的功率谱和亮温谱还可以对数字谱仪的性能进行评估.

利用MATLAB 产生2 GHz 正弦信号，采样率设定为4.8 GHz，将产生的信号通过核心算法进行运算，经过一定次数的累加后输出，并和原始信号直接进行FFT 运算，并经同样累加次数的结果进行对比. 图4（a）给出了同样累加次数的情况下，2 GHz 正弦信号直接进行FFT 运算和通过核心算法运算后的功率谱结果和峰值通道对比，图4（b）为2 种算法对应通道的功率结果差值.

从图4中可以看出，对于输入的正弦信号，2种算法的功率谱结果高度重合，且峰值通道完全一致，每个通道对应差值的量级为10-3，在峰值处因为信号能量较大，所以该处附近差值的散点相对聚集，而其他地方差值的散点相对分散.

图4 2GHz正弦信号仿真结果

除了模拟正弦信号外，还可以利用大气辐射传输模拟器（Atmospheric Radiative Transfer Simulator，ARTS）模拟临边探测时的真实大气亮温. 将由ARTS 产生的亮温值转化为进入辐射计的热噪声，从而可以根据谱辐射计的仿真模型得到功率谱的结果，对功率谱进行定标可以得到真实的亮温谱分布［19］，这样就可以从实际的目标场景角度对宽带数字谱仪进行仿真验证. 这里以240 GHz 辐射计的目标亮温谱为例进行说明，仿真条件设置为辐射计噪声温度1 000 K 以及积分时间50 ms.

图5 为240 GHz 辐射计中利用ARTS 模拟的天底角大气亮温，其中包含了所探测目标O3和HNO3的吸收峰. 辐射计接收的是目标的热噪声，根据谱辐射计的仿真模型，经过热辐射噪声信号模型、射频前端模型、数字后端模型和两点定标模型可以得到功率谱结果和亮温谱结果.

图5 240 GHz辐射计中模拟的天底角大气亮温

因为直接FFT 算法和核心算法对应的功率谱和亮温谱都高度重合，所以这里仅列出2 种算法对应通道的功率结果差值和对应通道的亮温结果差值，分别如图6、图7所示.

图6 临边探测中真实场景目标的功率偏差仿真结果

从图6 和图7 可以看出，2 种算法对应通道的功率偏差量级都在10-5左右，2 种算法对应通道的亮温偏差最大值为0.02 K 左右，通过模拟不同的信号从多个维度验证了核心算法的合理性和正确性.

图7 临边探测中真实场景目标的亮温偏差仿真结果

灵敏度是辐射计的一项重要指标，可以将它定义为接收机在频域上最小可分辨的亮温，通过标准差的分析方法可以获得辐射计的灵敏度. 此外，灵敏度的计算公式可以简化为

其中，Tsys表示辐射计的系统温度，β表示等效带宽，dτ表示单次测量的积分时间. 将2 GHz 全带内以自相关方式计算得到的灵敏度认为是理想全带灵敏度，在系统温度、积分时间、通道数等条件保持不变的情况下，理想子带灵敏度和理想全带灵敏度应该存在C之间的关系，其中C表示通道数.

通过对功率谱的结果进行定标分析，建立定标方程后多次仿真得到的理想全带灵敏度的均值为0.07 K左右，经过核心算法后子带灵敏度的均值为2.62 K 左右. 多次仿真结果表明，核心算法子带灵敏度与理想子带灵敏度相比恶化程度低于20%. 图8为其中一次仿真后核心算法子带灵敏度的结果.

图8 核心算法子带灵敏度仿真结果

5 测试结果分析

数字谱仪的谱分辨率Δf、信噪比（Signal to Noise Ratio，SNR）、信纳比（Signal to Noise and Distortion Ra⁃tio，SINAD）和有效位数（Effective Number Of Bits，ENOB）是衡量数字谱仪工作性能的重要指标. Δf，SNR，SINAD 和ENOB 的定义分别如式（25）～（28）所示［20，21］. 式（25）中fs和K分别代表采样频率和采样点数. 式（26）中PS和PN分别代表了信号和噪声的功率.式（27）中PD代表失真功率，式（28）中根据SINAD 的结果可以计算ENOB.

因为高光谱微波辐射计本身测量的就是噪声信号，但是在实际测试中直接产生宽带的噪声信号是较为困难的，所以同前期的仿真保持一致，实测中同样产生2 GHz的正弦信号来进行数字谱仪的性能分析.2 GHz正弦信号实测谱线的频率-功率图如图9所示，图9的子图中是对应的通道-功率图.

图9 2GHz正弦信号实测的功率谱结果

从图9 的通道-功率图中可以看出，信号实测功率谱和仿真功率谱的峰值通道保持一致，而且在该处附近能量较大. 因为受到ADC 量化噪声、内核相位偏差以及谐波的影响，在信号之外会产生一些小的峰值，由于其功率都在-40 dB以下，所以对数字谱仪ENOB的影响也很小. 从图9 频率-功率图中可以看出，对于2 GHz点频信号，实测功率谱的峰值频率为1 999 MHz，和输入的正弦信号相比偏差在1 MHz 左右. 而实际实际情况中，气体分子的每条吸收谱线，并不是只在一个频率上存在吸收的直线形式，一些效应会使吸收谱线发生谱线展宽，从而具有一定的宽度. 所以在谱分辨率恒定的情况下，1 MHz 的频率偏差可以满足本文中2.3 MHz 谱分辨率的设计要求.

根据得到的谱线结果，按照式（26）～（28），可以分别计算SRN，SINAD 和ENOB. 本设计方案中，结合不同频率的测试结果，实测的ENOB 至少为3 bit，同时SNR优于20 dB，同等条件下和仿真结果相比SNR 降低了2%左右.

6 结论

高光谱微波辐射计主要是利用后端的数字谱仪来实现模拟信号的谱分析，从而得到各种遥测大气成分的功率谱结果. 本文提出了一种新型的基于高速ADC和FPGA-Kintex7 的宽带、实时数字谱仪的设计和实现方案，并通过改进、结合FFT 的并行处理和复数处理，形成了核心的基16 实时复数FFT 算法，可以完成并行信号的实时处理. 对正弦信号和ARTS 模拟的真实亮温信号进行仿真，结果表明，核心算法和直接FFT算法后对应通道的亮温偏差最大值为0.02 K 左右，而且与理想子带灵敏度相比，核心算法的子带灵敏度恶化程度低于20%，从多个维度验证了核心算法的合理性和正确性. 本文所设计的2 GHz 带宽、1k 通道的数字谱仪谱分辨率恒定为2.3 MHz，实测的ENOB 至少为3 bit，同时SNR 优于20 dB，可以满足中高层大气谱线探测的需求，同时其硬件设计与核心算法设计的方案为实现更高带宽的数字谱分析系统提供了参考，对于临边探测中高层大气的痕量气体观测具有重要的意义.