FPGA实时多普勒中心频率估计

甘翼，李向阳

(中国电子科技集团公司第十研究所，四川 成都 610036)

0 引言

合成孔径雷达(synthetic aperture radar,SAR)受气象条件的影响较小，具有全天候工作的能力，在测绘、军事等领域具有广范的应用[1-2]。弹载SAR能够实现对目标的快速高精度成像和持续跟踪，是导弹实现远距离精确打击的重要传感器之一[3]。弹载SAR相对其他平台体积小、质量轻、功耗低，并且对处理实时性和大场景处理能力提出了更高的要求。因此，在弹载SAR实现方面，基于纯FPGA的SAR成像实现方法，可以提高大场景处理能力、减小成像时间、提高系统集成度、降低系统功耗[4]，在并行处理能力和成像时间方面相对DSP处理器具备极大的优势[5]。

SAR成像算法包括距离向压缩、距离走动矫正、方位向压缩以及多普勒中心频率估计、多普勒调频率估计等模块。其中距离向压缩、距离走动矫正和方位向压缩主要完成FFT变换、点乘和IFFT变换，具有单一的数据处理结构，无迭代操作，因而适合采用FPGA(field-programmable gate array)进行流水化处理。多普勒中心频率(frequency of Doppler centroid,FDC)估计、多普勒调频率(frequency of Doppler rate,FDR)估计[6]，因为算法相对较为复杂，使用FPGA实现具有较高的难度[7]，该部分通常采用DSP处理器完成[8]。这种折中方案会引入的数据交换延迟和FPGA处理等待延迟，大幅增加系统成像处理时间。以常用的DSP实现多普勒参数估计流程为例：首先，为了在DSP上实现参数估计，需要FPGA向DSP通过SRIO发送脉压数据，而距离向4096、方位向2048的IQ单精度浮点数据高达64 MB，通过2Lane@5Gbps的SRIO单次传输需要64 ms；其次，在DSP进行FDC估计时，因为后续距离走动校正需要多普勒中心频率估计结果，所以FPGA需要暂停当前流程，等待DSP的FDC处理结果才能进行后续操作，通常会产生高达80 ms的处理等待延迟；因此，基于FPGA+DSP的SAR处理流程在弹载条件硬件资源受限的情况下，无法满足大场景实时成像的要求，只能通过减小成像区域，来满足系统对实时性的要求。而通过将FDC估计采用FPGA实现即可大幅降低系统处理时间，提高系统的实时性。文献[9]提出了一种基于Zynq系列FPGA和RD算法的SAR成像实现方法，其中运动参数估计、多普勒参数估计、多普勒调频率估计等放在ARM处理器上进行处理，而距离向压缩、方位向压缩以及距离走动校正则是放在FPGA上进行。文献[10-12]提出了3种基于FPGA的SAR实现方法，其中文献[10]没有对参数估计部分进行阐述，文献[11-12]分别提出了基于FPGA和CS算法的SAR成像实现方法，并且采用FPGA实现多普勒参数估计。但是在多普勒参数估计部分均采用一次曲线进行多普勒中心频率拟合。而在弹载大斜视系统中，为了提高多普勒中心频率估计的准确性，通常采用2次或者3次曲线拟合进行逼近[13]。

本文提出一种基于FPGA的FDC估计方法，通过调整FDC数据处理流程、实现了SAR信号处理机的“零等待”。在FDC估计算法实现过程中，提出了基于FPGA的流水线式相位解缠绕方法和基于最小二乘的参数曲线拟合算法，实现相位解缠绕的连续流水输出以及中心频率估计参数精确插值。

1 多普勒中心频率估计算法

多普勒中心频率的估计有2种途径：①通过惯导获得[14]；②通过算法从回波数据中提取。由于惯导本身存在测量误差，为了获得更好的聚焦性能，通常从回波数据中通过算法提取多普勒中心频率参数15]。利用回波数据对多普勒中心频率进行估计有多种方法[16]。其中，能量均衡法、匹配相关法和最大似然法需要把距离向功率谱与某一参考函数做相关，把相关函数的零点作为多普勒中心[17]；相关函数法和数据序列符号自相关法则利用回波数据中的相位信息通过自相关提取多普勒中心，不需要额外的参考函数[18]。基于FPGA的可实现性，本文采用相位函数法实现多普勒中心估计[19]。由维纳-辛钦定理可知，相关函数r(τ)可以表示为其功率谱函数E{S(f)}的IDFT变换。

r(τ)=F-1{E[S(f)]}.

(1)

对于雷达发射系统，假设距离向功率谱函数E{S(f)}为天线功率方向图W(f)，则

r(τ)=F-1{W(f-fdc)}=r0(τ)exp{j2πfdcτ}.

(2)

从式(2)中可以看出自相关函数r(τ)是多普勒中心fdc的函数，当采用与k个方位向上的数据进行自相关时，也即τ=kT，其中T=1/PRF为重频周期，fdc可以表示为

(3)

不同方位向上回波S(nT)的自相关函数可以表示为

r(kT)=E{s*(nT)s(nT+kT)}∝

(4)

从式(3)和式(4)通过不同方位向上相同距离门的回波函数相关，即可得到fdc。因为非明显场景中无法提取正确的多普勒信息，需根据回波能量对距离门进行筛选，截取有效的场景数据进行多普勒相关估计。同时为了降低计算量，在进行方位向相关之前，会对距离门上的数据进行抽取，最终只能得到部分多普勒频率估计结果。为了得到完整的多普勒中心频率，提高多普勒中心频率估计的正确性，需要对多普勒频率估计结果进行二次曲线拟合和插值处理，通过求均值最终得到多普勒中心频率，用于距离走动校正。

由于FPGA内部存储资源的限制，无法保存所有的距离向数据，因而本文采用流水线的方位向相关累积方法，无需保存距离向数据，大幅降低对存储资源的消耗。在基于最小二乘的二次曲线拟合中，需要进行矩阵求逆操作，本文采用矩阵分解的形式，在FPGA上直接实现3×3矩阵的求逆操作。下面将详细介绍基于FPGA的多普勒中心估计处理流程和方法。

2 “零等待”的多普勒中心估计处理流程

多普勒中心频率估计包含距离向积累抽取、方位向相关积累、相位解缠绕和插值拟合等。FPGA的处理流程和DSP的处理流程[20]对比如图1所示。在DSP中距离向脉压块数据的读写和多普勒中心频率估计需顺序进行，相比FPGA实现，需等待一个孔径或子孔径的距离向脉压数据全部存储后，再完整读取进行多普勒参数估计，随后还须增加一次DDR数据读取时间，才能进行RCMC处理。而RCMC依赖多普勒中心频率的估计结果。而此时多普勒中心频率估计尚未完成，因此，一个成像孔径或子孔径的距离向脉压数据需要缓存。以孔径大小为4 096×2 048个处理单元为例，距离向脉压后的单精度IQ数据为64 MB，由于数据量较大，FPGA内部BRAM资源无法满足要求，需要转储到外部大容量DDR中，与前述DSP实现方法不同，FPGA可以发挥其并行处理的优势，实现边存储边处理。如图1所示，在对脉压结果进行存储的过程中，即可同步进行多普勒中心频率的估计。距离向压缩结果存储完成时，即可同步完成多普勒中心频率估计，随即开展RCMC处理，实现系统数据处理的“零等待”，降低了系统的处理时间，同时也简化了系统设计。

图1 FPGA与DSP对比处理流程Fig.1 Comparison of DSP and FPGA processes

如图2所示，基于FPGA的多普勒中心频率估计可以采用流水线的方式进行数据处理。传统非流水线处理方式需要的时间为t=N(tddr+tfdc)，而流水线的处理方式其时间消耗为t=Ntddr+tfdc，其中tddr为DDR读取一个距离向数据所需要的时间，tfdc为处理一个距离向所需要的时间，N为CPI周期内方位向积累点数。在相同tddr和tfdc的情况下，通过流水线的处理方式，可以“隐藏”N-1次FDC处理时间，对于4 096个距离向脉压数据，在数据接收完成后，只需额外的一个处理周期即可外出多普勒中心频率的估计，亦即“隐藏”了4 095次FDC处理时间。

图2 基于FPGA流水与非流水对比处理流程Fig.2 Comparison between FPGA-based pipelined and non-pipelined processes

以上是基于FPGA的多普勒中心估计采用的基本流程和优化方法，下面将从距离向积累抽取、相位解缠绕、二次曲线拟合和插值3个方面详细介绍实现过程。

2.1 距离向数据的积累抽取

为了提高多普勒中心频率估计的鲁棒性，降低系统噪声和场景突变对参数估计的影响，采用邻域积累实现距离向数据的平滑滤波。虽然会降低多普勒中心频率估计的精度，但是在弹载SAR面临的复杂电磁环境中可以提高算法的鲁棒性。同时为了降低方位向相关的数据处理量，对积累后的数据进行1/4抽取。如图3所示，在FPGA的实现过程中，采用定点累加的方式完成四邻域数据的积累。通过2 bit计数器对数据进行循环处理，当循环计数值为‘0’时，使用当前距离门数据对累加器进行初始化，其他情况下进行距离门数据累加，在累加完成4组数据后取走处理结果。虽然理论上应该采用均值来表示累积的结果，但是每个距离门上的邻域累积具有相同的长度，所以直接采用累加结果表示平滑滤波结果，从而可以减少逻辑资源的消耗。

图3 距离向平滑抽取Fig.3 Smooth extraction of range dimension

根据FPGA数据处理的特点，以及距离向脉冲压缩的实现过程，脉压后的数据与回波数据一样是按照固定重复周期采用流水线的方式依次接收和处理。在多普勒中心频率估计中，距离向累积抽取可以直接在数据流上直接完成。而方位向相关累积数据则需要缓存。

数据缓存在FPGA中可以通过LE(logic elements)资源实现，也可以通过内部BRAM实现。LE寄存器具有操作简单、读取延时低等特点。但是对于4 096点的64 bit复数据，其资源消耗是巨大的，会占用xc7v690t芯片20%的LE资源。通过BRAM进行缓存可以使用FPGA内部的专用存储资源，存储容量相对较多，但是BRAM会产生1～2个时钟周期的读写延迟。为了实现数据的流水线处理，需要进行预读取。对于距离向脉压数据，在4邻域平滑累加的4个时钟周期内，需要完成一次读取和一次写入操作。对于采用的2时钟延时的BRAM，数据读取操作需要在4邻域求和计算输出完成2个时钟周期之前进行，写入操作则在求和计算结果之后，为了避免读写操作的冲突，通常安排在求和结束后的第2个时钟周期写入。而对于方位向相关积累来说，则在5个时钟周期内需完成读取、累加和写入操作。其处理时间为3个时钟周期，预读取需要2个时钟周期。

如图4所示，方位向相关积累需要在CPI内的所有PRT的每个距离门上进行，抽取后的每个距离门需要独立的寄存器用于缓存方位向相关积累结果，对于4 096的距离向长度，需要2组1 024点深度的寄存器组用于分别存放IQ数据的相关累加结果。

图4 方位向相关积累Fig.4 Correlation accumulation of azimuth range

方位向相关累积完成后，对存放在BRAM中的各个距离门累加结果依次读取，作为cordic模块的输入，进行arctan操作。cordic生成的只是相关累积结果的瞬时相位信息，范围在[-π,π)之间，因而瞬时相位信息和多普勒频率相位之间需要进行解缠绕操作。

2.2 解缠绕

多普勒中心频率反映了脉冲发射单元与目标之间的径向速度差。不同PRT之间的相关性反映了不同距离门上的多普勒变化率，这是一种相位累积效应，这种累积结果通过反正切变换为角度信息时，映射在[-π,π)的范围内。解缠绕操作可以恢复相位对应于角度上的累加信息。当角度信息在-π和-π附近变化时，通过判断邻域数据的变化率，增加或者减去2π，来达到相位连续累加的目的。在得到多普勒相位值后，根据其对应圆周角度的比例，通过PRT换算即可得到多普勒频率。由于PRT的离散数据采集造成的多普勒模糊[21]，则通过惯导信息来提取多普勒模糊数进行修正，从而可以得到实际的多普勒频率值。解缠绕的FPGA实现伪代码如表1所示。

表1 解缠绕FPGA实现伪代码Table 1 Unwrapping pseudo code

2.3 基于最小二乘的曲线拟合

如前所述，为了降低计算量，多普勒中心频率估计是1/4抽取后进行处理的，但是在方位向压缩中，需要计算每个距离门上的多普勒频率信息，为了得到完整的多普勒频率估计值，需要对多普勒频率的估计值进行插值。本文采用了二次曲线拟合的方法来实现缺失部分多普勒中心频率的插值。基于最小二乘法的曲线拟合方法提供基于均方误差最小的最优解，可以有效地降低多普勒中心频率估计误差，提高估计精度。也能在一定程度上抑制由于平台抖动等带来的多普勒中心频率剧烈变化，对多普勒中心频率估计结果具有滤波的作用。多普勒中心频率在标准扫描SAR场景中可以用二次曲线对多普勒中心频率进行近似描述。

下面详细介绍基于最小二乘的二次曲线拟合方法中矩阵求逆的FPGA简化实现方法[22]。对于线性函数y=Ax，若A,y为已知量，x为待解量，则x可以表示为：x=inv(ATA)ATy。其中，由于矩阵A可能不是方阵，通过左乘AT得到方阵ATA。

对于二次曲线：

(5)

其矩阵A：

(6)

x=(a,b,c).

对于3×3矩阵，矩阵求逆操作可以通过代数余子式来实现[23]：

(7)

定义矩阵M=ATA，亦即

(8)

(9)

|A|=a11|A11|+a12|A12|+a13|A13|，

(10)

(11)

式中：

(12)

在FPGA实现过程中，矩阵M的行列式可以通过其代数余子式进行简化实现：

|M|=a11·detA11+a12·detA12+a13·detA13=

fx4·detA11+fx3·detA12+fx2·detA13.

(13)

在得到矩阵M的代数余子式和行列式后，即可以通过cordic模块实现除法操作，实现矩阵的求逆，从而可以求解方阵ATA的逆矩阵，从而可以实现二次曲线系数(a，b，c)。再得到二次曲线的系数后，可以通过动态生成x及其相关x2，从而完成二次曲线的插值。

综上所述，在FPGA实现基于最小二乘的二次曲线拟合中起主要关键点在于：①动态生成采用流水线方式生成x通过乘法器的流水线操作从x中得到x2，随后通过流水线方式直接输出fx,fx2,fx3,fx4。②通过代数余子式的形式求解方阵ATA的逆矩阵inv(ATA)，从而为使用FPGA实现最小二乘算法提供有力保证。

3 实验结果

基于FPGA的多普勒中心频率估计为弹载SAR纯FPGA实现的一部分，已经在某项目上实现，并完成了各项测试，所有指标均达到或优于相关要求。SAR成像模块采用一片Xilinx公司的Virtex 7系列xc7vx690t芯片外扩2片DDR3芯片，完成多普勒中心频率估计、RCMC、多普勒调频率估计以及方位向压缩，生成4 096×2 048大小的图像。系统信号处理模块主频为200 MHz，系统LE资源占用约占芯片总资源的58%。信号处理板数量由原理样机的4块减少为2块，处理功耗下降85 W以上。FDC模块资源使用情况如表2所示。其中arcTan函数为采用cordic完成相关数据的反正切计算；CorrComplex函数计算复数据相关；corrMemRe函数为使用LUT资源例化相关累积。lsqFDC函数使用最小二乘法实现二次曲线拟合系数估计。genFDC函数使用lsqFDC计算的系数，实时生成2 048点有效场景的FDC值；thetaWarp函数完成相位角解缠绕。FDC估计模块在RC数据接收过程中直接流水完成，不消耗额外的运行时间，接收完数据可以直接得到相关累积结果。为了使参数拟合和插值模块不产生等待，FDC相关累积模块抛弃最后一个方位向累积，直接转入多普勒参数拟合状态，因基于最小二乘的二次曲线拟合和插值模块以及后续的多普勒模糊数叠加模块总处理时间为4 124个时钟周期，在200 MHz处理时钟下，为2.1 μs，远小于重频周期。因而可以在RC数据接收完最后一包距离向脉压数据前直接输出FDC估计结果，直接转入RCMC处理，实现系统“0等待”，因此FDC估计模块系统消耗时间为0。

表2 FDC模块资源使用统计Table 2 FDC module resource usage statistics

4 结束语

本文根据弹载SAR对信号处理软硬件实现的需求，提出了一种基于FPGA的快速多普勒中心频率估计方法，采用流水线并行处理的模式，实现FDC估计时间的零消耗。为了实现基于FPGA的FDC估计，创新性地提出了基于FPGA的流水线式相位解缠绕方法和基于最小二乘的参数曲线拟合算法，实现相位解缠绕的连续流水输出以及中心频率估计的二次取下拟合和插值，并在实际工程中进行了有效验证。测试结果说明，该方法充分发挥了FPGA并行处理的优势，大大降低了SAR信号处理对硬件资源的需求，在弹载SAR和无人机载SAR等对结构、质量和功耗限制较严格的条件下，有着非常广阔的应用前景。