一种基于FPGA的二维平面阵数字波束形成方法

摘要

本发明属于阵列信号处理技术领域，具体涉及一种基于FPGA的二维平面阵数字波束形成方法。本发明针对大规模阵列情况下，采样矩阵协方差求逆(SMI)运算量过大情况，提出一种基于FPGA的二维平面阵数字波束形成方法，该方法基于共轭梯度算法(CG)的迭代求解权值，并进行了基于FPGA的硬件实现，以提高运算实时性。

权利要求

1.一种基于FPGA的二维平面阵数字波束形成方法，其特征在于，包括以下步骤：S1、设雷达天线阵列共有M×N个阵元，阵元分布于xoy平面，阵元沿x轴和y轴的阵元间距分别为dx和dy，接收信号为x(n)，M×N阵元的二维导向矢量为 θ为波束方位角，为波束俯仰角，将 用列向量表示为：为阵元m的相对相位差：

;

λ为信号波长；

S2、计算采样数据的协方差矩阵：将x(n)排成一维向量：

T

x(n)＝[x11(n) x12(n) … xMN(n)]根据采样数据x(n)和采样快拍数L，可得协方差矩阵估计S3、采用最小方差无失真响应进行波束形成，令最优权向量w满足：构造代价函数J(w)，代价函数取得最小值解即为所求权向量，J(w)表达式为：S4、采用共轭梯度算法迭代求解代价函数，具体为：将导向矢量 和协方差矩阵估计 作为FPGA计算模块的输入，并设置初始误差e，定义{p1,p2,p3,…}为一组共轭矢量，约束向量px与协方差矩阵估计 正交，满足x≠y，向量py是与向量px不同的向量，x，y为任意正整数，设置最大迭代次数2

imax；在进行迭代时，初始化迭代向量w0＝0，残差向量 迭代误差ρ0＝||r0|| ，迭代次数i＝0，p1＝r0；迭代过程为：S41、将导向矢量 以流水线方式输入，协方差矩阵估计 按行输入，每个时钟周期都并行MN路输入矩阵行元素，通过乘加运算得到辅助矢量S42、通过两路并行流水输入辅助矢量qi与共轭向量pi的每一个元素，先通过浮点数乘法IP核完成向量相乘操作，再通过浮点加法器与求和反馈完成复数累加计算，最后通过除法器进行除法运算，得到权值矢量步长因子S43、根据共轭向量pi、辅助矢量qi、权值矢量步长因子αi，通过乘法器和加法器完成乘法和累加运行，更新权值矢量wi+1＝wi+αipi和残差向量ri+1＝ri‑αiqi；

2

S44、设第i次误差ρi＝||ri+1|| ，判断终止条件是否ρi≤e或i＝imax是否成立，若是，则停止迭代，进入步骤S47，若否，则进入步骤S45；

S45、通过浮点除法器，计算共轭矢量迭代步长因子βi＝ρi/ρi‑1；

S46、计算第i+1次迭代的共轭向量pi+1＝ri+1+βipi，回到步骤S42；

S47、输出权向量w中的每个元素值；

S5、根据得到的权向量w，将其和接收信号x(n)输入多波束加权求和模块，具体为：将nH

时刻每个阵元权向量排成列向量w(n)，再将每个阵元通道的采样数据x(n)与加权矢量w(n)相乘求和，得到波束数据y(n)：

说明书

一种基于FPGA的二维平面阵数字波束形成方法

技术领域

[0001] 本发明属于阵列信号处理技术领域，具体涉及一种基于FPGA的二维平面阵数字波束形成方法。

背景技术

[0002] 数字波束形成技术应用于相控阵雷达能提高波束形成的灵活性，充分发挥其多波束和抗干扰能力优势。采样矩阵协方差求逆(SMI)的开环算法常被用于自适应波束形成。但
随着阵列规模的增大，采样矩阵协方差求逆运算量巨大，计算时延提高，难以保证波束形成
的低时延特性，影响雷达数据实时性。

发明内容

[0003] 本发明针对大规模阵列情况下，采样矩阵协方差求逆(SMI)运算量过大情况，提出一种基于FPGA的波束形成方法，该方法基于共轭梯度算法(CG)的迭代求解权值，并进行了
基于FPGA的硬件实现，以提高运算实时性。基于FPGA的CG迭代模块设计通用性，通过重复调
用或分时复用，可完成多波束权值计算。

[0004] 本发明的技术方案为：一种基于FPGA的波束形成方法，包括以下步骤：

[0005] S1、设雷达天线阵列共有M×N个阵元，阵元分布于xoy平面，阵元沿x轴和y轴的阵元间距分别为dx和dy，接收信号为x(n)，M×N阵元的二维导向矢量为 θ为波束方位
角，为波束俯仰角，将 用列向量表示为：

[0006]

[0007] 为阵元m的相对相位差：

[0008]

[0009] λ为信号波长；

[0010] S2、计算采样数据的协方差矩阵：

[0011] 将x(n)排成一维向量：

[0012] x(n)＝[x11(n)x12(n)…xMN(n)]T

[0013] 根据采样数据x(n)和采样快拍数L，可得协方差矩阵估计

[0014]

[0015] S3、采用最小方差无失真响应进行波束形成，令最优权向量w满足：

[0016]

[0017] 构造代价函数J(w)，代价函数取得最小值解即为所求权向量，J(w)表达式为：

[0018]

[0019] S4、采用共轭梯度算法迭代求解代价函数，具体为：

[0020] 将导向矢量 和协方差矩阵估计 作为FPGA计算模块的输入，并设置初始误差e，定义{p1,p2,p3,…}为一组共轭矢量，约束向量px与协方差矩阵估计 正交，满足
x≠y，x，y为任意正整数，设置最大迭代次数imax；在进行迭代时，初始化迭代向
2
量w0＝0，残差向量 迭代误差ρ0＝||r0||，迭代次数i＝0，p1＝r0；迭代过程为：

[0021] S41、将导向矢量 以流水线方式输入，协方差矩阵估计 按行输入，每个时钟周期都并行MN路输入矩阵行元素，通过乘加运算得到辅助矢量

[0022] S42、通过两路并行流水输入辅助矢量qi与共轭向量pi的每一个元素，先通过浮点数乘法IP核完成向量相乘操作，再通过浮点加法器与求和反馈完成复数累加计算，最后通
过除法器进行除法运算，得到权值矢量步长因子

[0023] S43、根据共轭向量pi、辅助矢量qi、权值矢量步长因子αi，通过乘法器和加法器完成乘法和累加运行，更新权值矢量wi+1＝wi+αipi和残差向量ri+1＝ri‑αiqi；

[0024] S44、判断终止条件是否ρi≤e或i＝imax是否成立，若是，则停止迭代，进入步骤S47，若否，则进入步骤S45；

[0025] S45、通过浮点除法器，计算共轭矢量迭代步长因子βi＝ρi/ρi‑1；

[0026] S46、计算第i+1次迭代的共轭向量pi+1＝ri+1+βipi，回到步骤S42；

[0027] S47、输出权向量w中的每个元素值；

[0028] S5、根据得到的权向量w，将其和接收信号x(n)输入多波束加权求和模块，具体为：将n时刻每个阵元权向量排成列向量w(n)，再将每个阵元通道的采样数据x(n)与加权矢量
H
w(n)相乘求和，得到波束数据y(n)：

[0029]

[0030] 本发明的有益效果是：

[0031] 1)模块通用性强

[0032] 子模块内部设计都以流水线方式运行，通过合理并行调用CG子模块可提高数据吞吐率，通过复用单个CG子模块可降低资源消耗，可根据FPGA资源情况结合使用，提高资源利
用率。

[0033] 2)运算实时性高

[0034] 采用基于FPGA的设计实现方案，充分利用并行和流水特点，极大提高了整体运算速，波束形成时延低。相较于使用CPU、DSP或“FPGA+DSP”处理模式，FPGA整体运算速度更快，
实时性更好，波束捷变性能好。

附图说明

[0035] 图1为本发明方法中迭代求解权值的程序流程图；

[0036] 图2为基于FPGA的共轭梯度算法迭代求解权值模块设计框图。

具体实施方式

[0037] 下面结合附图和实施例对本发明的技术方案进行详细描述。

[0038] 为便于理解，先对本发明方案中相关的计算公式进行说明：

[0039] 设阵列接收信号为x(n)，将二维导向矢量排成一维列向量形式表示为 M×N阵元的二维导向矢量 用列向量表示为

[0040]

[0041] 由公式1，采用最小方差无失真响应进行波束形成，最优权向量应满足

[0042]

[0043] μ为任意常数，计算出的权值w只存在倍数差异，对方向图形状没有影响，结合协方差矩阵采样快拍数L，可令μ＝1/L，也能省去协方差矩阵估计中求和后除法操作，协方差矩
记为 可得

[0044]

[0045] 构造代价函数J(w)，代价函数取得最小值解即为所求权向量，J(w)可表达式为

[0046]

[0047] 定义{p1,p2,p3,…}为一组共轭矢量，约束向量px与协方差矩阵估计 正交，满足关系式

[0048]

[0049] 设迭代步长因子为αi，延pi方向进行搜索可得权值w的迭代形式为

[0050] wi+1＝wi+αipi,i＝1,2,3…   (公式6)

[0051] 将wi+1代入代价函数得到J(wi+1)，之后令代价函数对αi求偏导值为0。再定义辅助矢量 可得αi表达式为

[0052]

[0053] 其中ri为第i次迭代的残差矢量，残差矢量迭代关系可表示为

[0054] ri+1＝ri‑αiqi,i＝1,2,3…   (公式8)

[0055] 令共轭矢量pi的迭代步长因子为βi，根据第i+1次迭代残差矢量ri+1可得到共轭矢量pi+1表达式为

[0056] pi+1＝ri+1+βipi,i＝1,2,3…   (公式9)

[0057] βi的值可由残差矢量ri+1和ri表示为

[0058]

[0059] 设第i次误差ρi＝||ri+1||2，若ρi小于设定的误差门限e，或迭代次数i大于设定的最大迭代次数imax则停止迭代，否则继续迭代。

[0060] CG算法迭代流程图如附图1所示，具体迭代步骤如下：

[0061] 输入：导向矢量 协方差矩阵估计 初始误差e。

[0062] 初始化：进行迭代向量初始化。初始迭代向量w0＝0，残差向量 迭代误2
差ρ0＝||r0||，迭代次数i＝0，p1＝r0。

[0063] 步骤1：计算第i次迭代的辅助矢量

[0064] 步骤2：由ρi、qi计算第i次迭代的权值矢量步长因子

[0065] 步骤3：根据(公式6)更新权值矢量wi+1＝wi+αipi；根据(公式8)更新残差矢量ri+1＝ri‑αiqi。

[0066] 步骤4：计算迭代误差ρi＝||ri+1||2。

[0067] 步骤5：判断终止条件是否ρi≤e或i＝imax，若成立，则停止迭代，跳转至输出。若不成立，继续进行步骤6。

[0068] 步骤6：根据(公式10)计算共轭矢量迭代步长因子βi＝ρi/ρi‑1。

[0069] 步骤7：根据(公式9)计算第i+1次迭代的共轭向量pi+1＝ri+1+βipi，并跳转步骤2继续执行迭代。

[0070] 输出：权矢量w＝wi+1。

[0071] 基于FPGA的CG算法实现方法如下：

[0072] 根据算法流程图进行每一步的迭代子模块设计。得到基于FPGA的CG算法实现设计框图如附图2所示。

[0073] 输入：导向矢量 协方差矩阵估计 初始误差e。

[0074] 导向矢量 以流水线方式输入，赋值到内部寄存器，协方差矩阵估计 按行输入加法树结构、输入初始误差门限值e。

[0075] 初始化：复位默认w0＝0，缓存 作为初始辅助向量p1，缓存RAM记为RAM‑p，开始如下迭代步骤。

[0076] 步骤1：计算辅助矢量

[0077] 此步骤功能模块记为子模块1。其中乘和加运算模块代表复数乘法和复数加法。每个时钟周期都并行MN路输入矩阵行元素，可保证数据流水输入输出特性，提高数据吞吐率。
输出端辅助矢量qi输出时同步缓存入RAM中，记为RAM‑q。由于是复数，需并行调用两个RAM
同步缓存qi实部虚部，待步骤3继续使用。

[0078] 步骤2：计算权值矢量步长因子

[0079] 此步骤功能模块记为子模块2。两路并行流水输入辅助矢量qi与共轭向量pi的每一个元素，调用浮点数乘法IP核Floating point 7.1完成向量相乘操作。乘法运算后调用浮
点加法器与求和反馈完成复数累加计算，除法运算调用除法器，输入端读入上次迭代的误
差相除得到权值迭代步长αi，输出至步骤3中子模块3、4，开始步骤3。

[0080] 步骤3：权值向量w和残差向量r更新。

[0081] 根据(公式6)更新权值矢量wi+1＝wi+αipi；根据(公式8)更新残差矢量ri+1＝ri‑αiqi。权值向量更新模块记为子模块3，残差向量更新模块记为子模块4。流水读出RAM‑p中的
辅助向量p的每个元素至子模块3，流水读出RAM‑q种元素至子模块4，保持两个子模块并行
处理，由乘法器和加法器完成输入乘法和累加运算，输出向量w、r的元素分别缓存至RAM‑w、
RAM‑r中待下一次迭代使用。

[0082] 步骤4：计算迭代误差ρi＝||ri+1||2。

[0083] 此步骤运算模块记为子模块5，完成残差向量r的求模运算。输入端调用浮点数乘法器流水输入每个残差向量元素ri完成共轭相乘，输出至浮点加法器，通过循环反馈完成
累加操作。

[0084] 步骤5：判断终止条件是否ρi≤e或i＝imax，若成立，则停止迭代，跳转至输出。若不成立，继续进行步骤6。

[0085] 此步骤负责逻辑判断和权值输出控制，逻辑判断部分记为子模块6，将子模块5输出的误差值ρi与初始误差门限e相比较，低于门限值则将权矢量输出使能拉高，跳转至输
出，读出RAM‑w中缓存的权值w，停止迭代。否则继续步骤6。

[0086] 步骤6：计算共轭矢量迭代步长因子βi＝ρi/ρi‑1。

[0087] 此步骤完成共轭矢量pi的迭代步长βi计算，功能模块记为子模块7,根据(公式10)，通过调用浮点除法器完成除法计算，输入端被除数为步骤5输出的误差值ρi，除数为上一次
缓存的迭代误差ρi‑1，输出迭代步长βi。开始步骤7。

[0088] 步骤7：计算第i+1次迭代的共轭向量pi+1＝ri+1+βipi，之后继续步骤2。

[0089] 此步骤完成共轭矢量pi更新，功能模块记为子模块8。根据(公式9)，输入迭代步长βi，流水读出上一次迭代缓存的的共轭向量pi‑1。跟新方式与步骤3中权值向量更新类似，由
浮点乘法器和加法器完成累加操作。同样调用RAM缓存共轭矢量pi，等待下一次迭代使用，
继续步骤2。

[0090] 输出：连续读出RAM‑w中缓存的权值向量w的每个元素值wmn。

[0091] 若有K个方向导向矢量，可调用多个CG模块完成多路并行迭代，或一个CG模块复用K次分时迭代。当CG模块迭代计算完一组权值则流水输出权值w至加权求和模块。

[0092] 为进一步说明本发明的实用性，下面给出本发明方案下的时延和FPGA资源消耗说明。

[0093] 为保证计算精度，FPGA采用位浮点数计算，数据为复数。浮点复乘时延设为T1,浮点加减法时延设为T2,除法时延记为T3。设迭代次数为K，FPGA主时钟为Fm，CG算法迭代模块
的运行时延TCG可按(公式11)、(公式12)计算。

[0094]

[0095] TCG＝k(TS1+TS2+TS3+TS4)+TS5+(k‑1)(TS6+TS7)   (公式12)

[0096] 取16×16规模协方差阵为例计算，迭代次数k＝6，典型时延T1＝T2＝4、T3＝28，判断时延T4＝2，可计算出TCG＝6×(20+44+8+26)+2+5×(28+8)＝770。FPGA主时钟Fm设为
100MHZ，CG算法迭代模块求出权值时间为770×10ns＝7.7μs。

[0097] 由于算法进行乘法加法操作主要消耗DSP资源，这里计算DSP资源消耗。

[0098] 浮点数复乘消耗DSP slice设为D1，复加消耗DSP slice设为D2。算法每一步消耗DSP slice资源及总消耗可按(公式13)、(公式14)计算。

[0099]

[0100] DCG＝DS1+DS2+DS3+DS4+DS7   (公式14)

[0101] 一个复乘拆成实部虚部后包含4个乘法器和2个加法器。Xilinx IP中，选择全用DSP slice模式时，一个单精度浮点乘法器最多消耗2个DSP slice，单精度浮点加法器最多
消耗2个DSP slice。因此，可计算D1＝2×(4+2)＝12，D2＝2×2＝4，以16阵元为例，最大总消
耗DSP slice个数可计算DCG＝256+24+32+24+16＝352个。若使用选用一半DSP slice一半
logic cell模式，DSP资源消耗也可减半，即352/2＝176个。目前Xilinx 7系列支持最大Dsp 
slice个数S‑7系列为160，A‑7系列为740，K‑7系列为1920，V‑7系列为3600，一般7系列FPGA
便能满足设计所需DSP资源要求。

[0102] 实施例

[0103] 本例阵面有为4×4＝16个阵元，协方差矩阵 规模为16×16。流水输入16个阵元导向矢量 并赋值到内部寄存器。缓存 作为初始辅助向量p1，缓存RAM记为
RAM‑p，初始w0＝0。按行输入协方差矩阵估计 到子模块1的加法树结构。输入初始误差门限
值e。开始如下迭代步骤。

[0104] 步骤1：计算辅助矢量

[0105] 每个时钟周期都并行16路输入矩阵行元素，输出端辅助矢量qi输出时同步缓存入RAM中，记为RAM‑q。由于是复数，需并行调用两个RAM同步缓存qi实部虚部，待步骤c继续使
用。

[0106] 步骤2：计算权值矢量步长因子

[0107] 两路并行流水输入辅助矢量qi与共轭向量pi的每一个元素，调用浮点数乘法IP核Floating point 7.1完成向量相乘操作。乘法运算后调用浮点加法器与求和反馈完成复数
累加计算，除法运算调用除法器，输入端读入上次迭代的误差相除得到权值迭代步长αi，输
出至步骤3中子模块3、4，开始步骤c。

[0108] 步骤3：权值向量w和残差向量r更新。

[0109] 根据(公式6)更新权值矢量wi+1＝wi+αipi；根据(公式8)更新残差矢量ri+1＝ri‑αiqi。保持两个子模块并行处理，由乘法器和加法器完成输入乘法和累加运算，输出向量w、r
的元素分别缓存至RAM‑w、RAM‑r中待下一次迭代使用。

[0110] 步骤4：计算迭代误差ρi＝||ri+1||2。

[0111] 输入端调用浮点数乘法器流水输入每个残差向量元素ri完成共轭相乘，输出至浮点加法器，通过循环反馈完成累加操作。

[0112] 步骤5：判断终止条件是否ρi≤e或i＝imax，若成立，则停止迭代，跳转至输出。若不成立，继续进行步骤f。

[0113] 将步骤4中计算输出的误差值ρi与初始误差门限e相比较，低于门限值则将权矢量输出使能拉高，跳转至输出，读出RAM‑w中缓存的权值w，停止迭代。否则继续步骤f。

[0114] 步骤6：计算共轭矢量迭代步长因子βi＝ρi/ρi‑1。

[0115] 根据(公式10)，通过调用浮点除法器完成除法计算，输入端被除数为步骤5输出的误差值ρi，除数为上一次缓存的迭代误差ρi‑1，输出迭代步长βi。开始步骤g。

[0116] 步骤7：计算第i+1次迭代的共轭向量pi+1＝ri+1+βipi，之后继续步骤1。

[0117] 根据(公式9)，输入迭代步长βi，流水读出上一次迭代缓存的的共轭向量pi‑1。跟新方式与步骤c中权值向量更新类似，由浮点乘法器和加法器完成累加操作。同样调用RAM缓
存共轭矢量pi，等待下一次迭代使用，继续步骤1。

[0118] 输出：连续读出RAM‑w中缓存的权值向量w的每个元素值w11～w44。

[0119] 并行调用k个CG迭代模块，或分时复用k次CG迭代模块，以相同的输入方式和迭代方式完成a～g步骤得到k个波束权值向量数据流w1～wk。

[0120] 步骤8：得到权值向量w后，即可采用目前的多波束加权求和模块，结合接收信号x(n)得到波束数据，得到k个波束权值向量，即可得到k个方向的波束数据yk(n)。