0 引　言

目标方位（Direction of Arrival，DOA）估计在水下目标定位、信号识别及干扰抑制上非常重要，其中主要影响DOA估计性能的包括水下阵列阵型设计、信号处理算法和硬件实现平台等因素^[1]。

互质阵列由Palghat P. Vaidyanathan与Piya Pal提出^[2]，是一种稀疏阵列，能够通过特殊的阵型排布抑制栅瓣的产生，避免角度模糊现象，提高角度估计的可靠性；在相同阵列物理孔径约束下可以用少量的水听器实现 DOA估计^[3]，适合无人平台对阵列轻量化、低成本、低功耗的应用需求。典型的互质阵列由两阵元数互为质数的子阵列嵌套而成，但存在阵列自由度低、接收信号利用不充分等问题。Qin等^[4]提出了一种基于具有位移子阵列的新型互质阵列（Coprime Array With Displaced Subarrays，CADis），该互质阵型具有自由度高、接收信号利用充分、阵元间耦合干扰小的优点，本文所讨论互质阵列阵型均为CADis阵型。

压缩感知算法在21世纪初被提出^[5]，在图像噪声剔除、画质优化等方面应用广泛。在DOA估计应用中，与多重信号分类算法（Multiple Signal Classification，MUSIC）和最小方差无畸变响应方法（Minimum Variance Distortionless Response，MVDR）不同，压缩感知算法利用目标在空间中为稀疏分布的特性，突破了Nyquist采样定理的限制，能够利用较少的样本信号估计出目标的真实方位。研究表明，基于压缩感知的DOA估计方法的目标检测能力、抗干扰能力均优于传统方法^[6]。在均匀线阵进行压缩感知DOA估计时，由于信号的稀疏性难以保证导致重构效果不佳。而在互质阵中虚拟阵列的构建使得对信号参数进行估计时角度分辨率的提高，进而能够对解凸优化过程进行一定的化简，使之更适合硬件平台的实现，取得更优的DOA估计结果。

在DOA估计的硬件实现中，由于精确性、实时性的要求，对硬件平台的计算能力要求较高。现阶段接收信号大多为宽带信号，接收阵列进行宽带DOA估计时多采用DSP平台进行处理，而FPGA相较于DSP拥有更好的数据处理能力和处理速度，在进行多通道低时延波束形成时具有更大优势。进入新世纪之后，GPU器件的兴起给了DOA估计实现的另一种可能，但GPU器件的高功耗和运算延时的不确定性具有一定限制。因此本文选择FPGA器件作为主控芯片，在进行硬件实现时，阵列参数数据的储存和调用如果放置在FPGA芯片上将会消耗大量空间，需要外置储存器件的参与，本文引入DDR和SD卡进行协同配合实现存储调用。综上所述，本文针对水下DOA估计需求，设计实现了基于FPGA的互质阵列压缩感知DOA估计，提高了估计性能，满足实时性要求，为后续目标定位、信号识别提供有利信息。

1 互质阵信号模型及压缩感知DOA估计方法 1.1 互质阵信号模型

CADis互质阵型如图1所示，该阵型由两子阵并排放置而成，两子阵列间距为$ Ld $，其中$ L=\hat{M}+N $，子阵1阵元间距为$ \hat{M}d $，其中$ \hat{M}=1 $，子阵2阵元间距为$ Nd $。

一般地，假设$ K $个宽带信号入射在互质阵列上，则第$ m $个阵元的接收信号为：

$ {x}_{m}(t)=\sum \limits_{k=1}^{K}{s}_{k}(t-{\tau }_{m}({\theta }_{k}))+{n}_{m}(t)。$

(1)

式中：$ {s}_{k}(t) $为第$ k $个信号；$ {\tau }_{m}({\theta }_{k}) $为阵元$ m $接收第$ k $个信号相对于参考阵元的延时量；$ {n}_{m}(t) $为噪声。由于接收信号大多为宽带信号，不能简单将各阵元接收信号的时延分量等效为相移分量，因此需转换至频域分析。将单位时间内接收信号经过$ J $点傅里叶变换，得到频域快排模型。

$ {{\boldsymbol{X}}}({f}_{j})={{\boldsymbol{A}}}({f}_{j}){{\boldsymbol{S}}}({f}_{j})+{{\boldsymbol{N}}}({f}_{j}) 。$

(2)

式中：$ {f}_{j} $为频带中第$ j $个频点；$ {{\boldsymbol{X}}}({f}_{j}) $为阵列接收信号的频域表示；$ {{\boldsymbol{S}}}({f}_{j}) $和$ {{\boldsymbol{N}}}({f}_{j}) $分别为原信号和噪声的频域；$ {{\boldsymbol{A}}}({f}_{j}) $为互质阵列导向矢量。互质阵列导向矢量分别由图1中子阵1和子阵2的导向矢量组合而成，对于子阵1，其导向矢量$ {a}_{1}({\theta }_{k}) $写为：

$ \begin{split} {{{{\boldsymbol{a}}}}_1}({\theta _k}) =& [1,\exp ( - j2\text π Md{f_j}\sin {\theta _k}/c), \cdots , \\ & \exp ( - j2\text π (N - 1)Nd{f_j}\sin {\theta _k}/c){]^\text T} 。\end{split} $

(3)

对于子阵2，其导向矢量$ {a}_{2}({\theta }_{k}) $写为：

$ \begin{split} {{{{\boldsymbol{a}}}}_2}({\theta _k}) = &[1,\exp ( - j2\text π Nd{f_j}\sin {\theta _k}/c), \cdots , \\ & \exp ( - j2\text π (M - 2)Nd{f_j}\sin {\theta _k}/c){]^\text T} 。\end{split} $

(4)

式中：c为声速。构建协方差矩阵$ R $第$ i $行、第$ j $列的元素可以表示成：

$ {\begin{split}{{{\boldsymbol{R}}}}(i,j) =\, & {r_{i,j}} = E[{{{{\boldsymbol{X}}}}_i}({f_j}){{{\boldsymbol{X}}}}_j^*({f_j})] = \\ &\left\{ {\begin{array}{*{20}{l}} \displaystyle\sum\limits_{k = 1}^K {\sigma _k^2} \exp \left( - j\displaystyle \frac{{2\text π c}}{{{f_j}}} ({d_{\text{i}}} - {d_j})\sin {\theta _k}\right),i \ne j ，\\ \displaystyle \sum\limits_{k = 1}^K {\sigma _k^2} \exp \left( - j\displaystyle \frac{{2\text π c}}{{{f_j}}} ({d_{\text{i}}} - {d_j})\sin {\theta _k}\right) + \sigma _n^2,i = j 。\end{array}} \right. \end{split}}$

(5)

可以看出，协方差矩阵中的元素从形式上可以看成位置为$ {d}_{i}-{d}_{j} $的均匀阵列阵元所接收的单快拍信号，而该均匀线阵称之为虚拟阵列，位置差集可具体表示为：

$ {{\boldsymbol{S}}}=\left\{{d}_{i}-{d}_{j}\left| i,j=1,2,\cdots ,M+N-1\right.\right\}。$

(6)

将式(5)所构建的协方差矩阵$ R $矢量化，将其作为以位置差集所构建的虚拟阵列的接收信号，该信号具有单快拍、相干的性质。

$ {{\boldsymbol{v}}}=vec({{\boldsymbol{R}}}) 。$

(7)

式中：$ vec $为矩阵的矢量化操作，需要注意的是，$ {{\boldsymbol{v}}} $中可能存在部分冗余，需去除冗余。

1.2 压缩感知DOA估计方法

压缩感知理论是利用信号的稀疏性对低维信号进行重构。对于$ N $维实信号$ {{\boldsymbol{X}}}\in {{{\boldsymbol{R}}}}^{N} $，当信号含有$ K $个非零值，则称信号$ {{\boldsymbol{X}}} $是稀疏的($ K \lt \lt N $)。观测信号为$ {{\boldsymbol{Y}}}\in {{{\boldsymbol{R}}}}^{M} $时，重构过程可表示为：

$ \min ||{{\boldsymbol{X}}}|{|}_{0}s.t.{{\boldsymbol{Y}}}={\boldsymbol\varPhi {\boldsymbol{X}}} 。$

(8)

式中：$ \boldsymbol{\varPhi }\in {{{\boldsymbol{R}}}}^{M\times N} $为观测矩阵。对于大多数信号很难满足稀疏性，但经过某个合适的稀疏基变换后，可表现出稀疏特性，即$ \boldsymbol{X}{=\boldsymbol{\varPsi\varTheta}\mathit{\mathit{ }}} $，其中$ \boldsymbol{\varPsi }\in {\boldsymbol{R}}^{N\times N} $为稀疏矩阵。需要注意的是，选用的稀疏基$ {{{{\boldsymbol{\varPsi }}}}}\in {{{\boldsymbol{R}}}}^{N\times N} $需与观测矩阵$ {{\boldsymbol{\varPhi }}}\in {{{\boldsymbol{R}}}}^{M\times N} $满足RIP性质，即$ {{\boldsymbol{\varPhi }}}\in {{{\boldsymbol{R}}}}^{M\times N} $与$ \boldsymbol{\varPsi }\in {\boldsymbol{R}}^{N\times N} $不相关。经过稀疏基变换后，式(8)等价为：

$ \min \left\|{\boldsymbol{\varPsi \varTheta}} \right\|_{0}s.t.{{\boldsymbol{Y}}}=\boldsymbol{\varPhi X}={{\boldsymbol{\varPhi \varPsi \varTheta }}}=\boldsymbol{A\varTheta }。$

(9)

式中：$ {{\boldsymbol{A}}}={{\boldsymbol{\varPhi \varPsi }}}\in {{{\boldsymbol{R}}}}^{M\times N} $称为感知矩阵。式中求解$ {l}_{0} $范数是个NP-hard问题，数学上通常用求解$ {l}_{1} $范数代替，可得：

$ \min ||{{\boldsymbol{\varTheta }}}|{|}_{1}s.t.{{\boldsymbol{Y}}}={{\boldsymbol{A\varTheta}} }。$

(10)

这使得问题转化为凸优化求解问题，通过求解方程可到到稀疏向量$ {{\boldsymbol{\varTheta }}}\in {{{\boldsymbol{R}}}}^{N} $。

进行DOA估计时，需将互质阵列信号模型对应到压缩感知原理式中。考虑到虚拟阵列信号为宽带信号，计算时需对每个频域子带分别求解，则对单一子带$ {f}_{j} $而言，式(2)变为：

$ {{{\boldsymbol{X}}}}_{a}({f}_{j})={{{\boldsymbol{A}}}}_{a}({f}_{j},\hat{\theta }){{{\boldsymbol{S}}}}_{a}({f}_{j})+{{{\boldsymbol{N}}}}_{a}({f}_{j})。$

(11)

式中：$ {{{\boldsymbol{X}}}}_{a}({f}_{j}) $表示式(7)中协方差向量化后的$ {{\boldsymbol{v}}} $，$ {{{\boldsymbol{A}}}}_{a}({f}_{j},\hat{\theta }) $表示虚拟阵列的导向矢量；$ \hat{\theta } $表示搜索方向；$ {{{\boldsymbol{S}}}}_{a}({f}_{j}) $和$ {{{\boldsymbol{N}}}}_{a}({f}_{j}) $为接收信号与噪声在频域中的稀疏表示。对应到式(10)中，将$ {{{\boldsymbol{X}}}}_{a}({f}_{j}) $作为观测序列$ Y $，$ {{{\boldsymbol{A}}}}_{a}({f}_{j},\hat{\theta }) $作为感知矩阵$ {{\boldsymbol{A}}} $，$ {{{\boldsymbol{S}}}}_{a}({f}_{j}) $作为待求解稀疏分量$ { {\boldsymbol{\varTheta}} } $，$ {\boldsymbol{N}}_{a}({f}_{j}) $为测量噪声^[7]。将式(11)代入式(10)，则可求解信号：

$ \min ||{{{\boldsymbol{S}}}}_{a}({f}_{j})|{|}_{1}s.t.{{{\boldsymbol{X}}}}_{a}({f}_{j})={{{\boldsymbol{A}}}}_{a}({f}_{j},\hat{\theta }){\boldsymbol{S}}_{a}({f}_{j})。$

(12)

通过求解对应搜索角度$ \hat{\theta } $的信号功率值$ {P}_{a}(\theta ) $来判断信号的波达方向，即$ {P}_{a}({f}_{j},\theta )=|{{{\boldsymbol{S}}}}_{a}({f}_{j}){|}^{2} $。则$ J $个子带平均意义上的功率谱为：

$ P(\theta )=\frac{1}{J}\sum \limits_{j=1}^{J}{{\boldsymbol{P}}}({f}_{j},\theta )=\frac{1}{J}\sum \limits_{j=1}^{J}|{{{\boldsymbol{S}}}}_{a}({f}_{j}){|}^{2}。$

(13)

求解式(12)的凸优化问题，本文采用正交匹配追踪（Orthogonal Matching Pursuit，OMP）算法，表1展示了阵列信号处理中OMP算法的详细步骤。

由于CADis互质阵列构建的虚拟阵元数远多于待估目标数，稀疏性可以保证，为使信号重建速度更快，稀疏度选择为1即可，即$ K=1 $。

2 开发环境及系统整体方案 2.1 开发环境

在硬件系统开发过程中，不再仅仅局限于算法本身的实现过程，同时应考虑的还有数据的传输调用过程、储存器件的选择、传输协议的利用、测试流程等方面，这是一个系统工程。因此在开发硬件系统时，需针对实际使用需求选择合适的开发平台、编程语言和验证方法。

FPAG最早的历史可追溯到20世纪80年代初，其出现标志着电子电路设计领域的一大突破，它允许工程师通过软件配置硬件电路的实现逻辑，给了数字电路涉及极大的自主性。尽管FPGA有着强大的可编程能力，但其仍在复杂逻辑、串行场景下存在较大的局限性，使得应用受到了极大的限制，在大多数串行应用时不得以重复等待或自行设计软核，这严重浪费了硬件处理资源。为了弥补这一缺点，FPGA逐渐加入了ARM端，将SoC处理器、存储器、接口等组件集成在同一芯片上，使得FPGA向着系统级芯片不断发展。

本文使用的开发板搭载一片XCZU3EG-SFVC784-2I芯片，它集成了ARM处理器和FPGA数字逻辑，FPGA中内置了大量IP核供用户调用，并提供了丰富的硬件逻辑资源，包括约70000个逻辑单元，能够满足本文的需要。该芯片从架构上分为两部分：系统端（PS端）和逻辑端（PL端）：PS端集成了高性能SoC处理器，能够运行操作系统和高级软件应用；PL端则由传统FPGA数字电路逻辑组成，支持硬件加速和定制化的编程。两者之间的传输通过高速AXI总线协议进行链接，使得硬件加速和软件控制能够更加高效。

本设计在PL与PS分别完成各自的开发部分，最后整合完成系统设计。根据选用芯片所支持的软件版本，PL部分编程需在Vivado2020.2中完成，编程语言为Verilog HDL；PS部分须在Vitis 2020.2中完成，编程语言为C/C++。其中PL部分开发为纯硬件开发，使用Verilog HDL语言，在Vivado软件中进行模块化编程；而PS部分编程除开Vivado软件外还需要使用到Vitis软件，首先在Vivado中的Block Design界面搭建嵌入式系统，完成时钟、使用接口、传输速率等设置后生成包含平台硬件信息的XSA（Xilinx Shell Archive）文件，随后导入Vitis进行进一步编程。在Vitis中完成软件部分设计，分别包括主程序编写、库函数调用、运行地址设置等工作，最后联合FPGA部分上板调试。PS端开发流程如图2所示。

2.2 系统整体方案设计

系统整体框图如图3所示，包含PL与PS两部分，其中PL部分为硬件电路逻辑端，PS部分为ARM系统端。PL端中主要负责对阵列接收信号的预处理和波束形成的实现，PS端主要负责对阵元扫描向量$ {{\boldsymbol{A}}}({f}_{j},\theta ) $的存储及发送。本文实现基于8元互质阵，CADis阵型，阵元位置为$ \{\text{0,1},\text{2},\text{3},\text{4},\text{10},\text{15,20}\} $，数据采用32位定点数（Q18.14，即1位符号位、17位整数、14位小数），能满足运算所需的精度要求。下面对算法实现过程进行介绍。

PL端主要包括时频转换模块、数据预处理模块、虚拟矩阵构建模块和波束形成模块。时频转换模块负责将时域信号通过傅里叶变换（fast Fourier transform，FFT）至频域，FFT算法的实现采用自带IP核实现，其具有方便、准确的特点。数据预处理模块主要完成对频域信号异常值的剔除，以免影响DOA估计的效果。虚拟矩阵构建模块通过对频域信号$ {{{\boldsymbol{X}}}}_{*}(k) $分子带求取协方差矩阵$ {{\boldsymbol{R}}}(i,j) $并将其向量化。协方差矩阵计算部分通过式(5)进行计算，矩阵的向量化部分依据式(7)将协方差矩阵中元素存入ROM后按时钟进行对应输出完成排列。

波束形成模块主要完成OMP算法的实现，式(12)及表1描述了这一过程。该模块中需完成对PS端所发送数据的接收，由于$ K=1 $时，表1步骤4中$ {({{{{\boldsymbol{E}}}}_{t}}{{}^{\text T}}{{\boldsymbol{E}}_{t}})}^{-1} {{{\boldsymbol{E}}}}_{t}{}^{\text T}{{\boldsymbol{y}}}\propto {{{\boldsymbol{E}}}}_{t}{}^{\text T}{{\boldsymbol{y}}} $，其中矩阵求逆部分$ {({{{{\boldsymbol{E}}}}_{t}}{{}^{\text T}}{{{{\boldsymbol{E}}}}_{t}})}^{-1} $为固定常数，因此仅计算$ {{{\boldsymbol{E}}}}_{t}^{\text{T}}y $即可。根据表1中步骤2，要完成$ {\mathbf{\lambda }}_{t}=\underset{j=1,\cdots ,N}{\arg \max }| \lt {{{\boldsymbol{r}}}}_{t-1},{\phi }_{j} \gt | $，其中$ {{{\boldsymbol{r}}}}_{t-1} $为式(7)所得去除冗余后的$ {{\boldsymbol{v}}} $，$ {\phi }_{j} $即为该频点处虚拟阵元扫描向量矩阵$ {{\boldsymbol{\varPhi }}} $的第$ j $列，按照矩阵乘法结果保存最大值索引序列，从$ {{\boldsymbol{\varPhi }}} $中取出对应向量作为$ {{{\boldsymbol{E}}}}_{t} $，随后计算$ {{{\boldsymbol{E}}}}_{t}{}^{\text T}\mathbf{y} $求得该频点处重构信号$ {\hat{{\boldsymbol{x}}}} $，随后根据式(13)求取方位功率谱。

在进行波束形成时，需动态调用阵元扫描向量$ {{\boldsymbol{A}}}({f}_{j},\theta ) $，但是将大量的数据存储在FPGA芯片上不仅会造成硬件存储资源的浪费，也会影响信号处理速度，因此可将其存储在DDR中，既方便数据的随时调用同时节约片上资源。受限于DDR的断电不保存特性，为避免数据的人工重复性写入，本文引入了SD卡，将阵元扫描向量以二进制格式存储在SD卡中，系统上电时自动将SD卡中的数据写入DDR，PS端在接收到阵元扫描向量数据调用信号后以DMA的方式通过AXI-stream协议发送至PL端。PS端与PL端数据交互过程如图4所示。

可以看出，在PS端与PL端的数据交互数据交互过程中主要分为数据流和控制流。具体交互过程为：在上电之初将存储在SD卡中的阵元扫描向量数据通过ARM内核写入到DDR4中，减少系统初始化的时间；同时通过HPM接口以AXI-lite协议对AXI-DMA传输中的突发传输长度、传输起始地址、传输速率等进行配置，完成对传输过程的控制；当接收到PL部分传来的传输申请时，ARM通过HP接口以AXI-steam协议通过DMA的方式传输至FIFO模块中，然后由波束形成模块顺序读出，完成数据流传输，其中FIFO模块起到数据在跨时钟域传输过程中的缓存作用，MM2S（Memory Mapped to Streaming Interface）为储存器到外设映射。

3 实验验证与结果分析 3.1 仿真数据验证

在数据仿真中设定信号频率范围为500～3000 Hz的线性调频信号，来向为80°，采样率为15 kHz，叠加信噪比为−10 dB高斯噪声，单拍处理信号时间长度为1 s，信号时域波形及功率谱如图5所示。

为验证FPGA的计算准确性，FPGA与Matlab分别进行运算，DOA估计结果如图6所示。可以看到，两结果一致，对两结果求取互相关系数来表征相似性，相关性为0.99，表明本设计进行DOA估计的有效性。

3.2 湖试数据验证

对某次湖上实验目标进行DOA估计。8元互质阵列阵元位置为$ \{\text{0,1},\text{2},\text{3},\text{4},\text{10},\text{15,20}\} $，与前序理论分析及仿真一致。同样取1 s数据进行分析，结果如图7所示。图7(a)为在线逻辑分析仪（Integrated Logic Analyzer，ILA）中对运行结果的抓取，可以看到系统正常运行。从图7(b)中可以看出FPGA与Matlab运算结果基本一致，误差主要原因是FPGA定点数计算与Matlab双精度浮点数计算之间的差异，并叠加各模块计算误差的累计，对其求取互相关系数为0.98，表明FPGA进行互质阵列压缩感知DOA估计的有效性。

对系统实现实时性进行分析，本设计系统时钟为100 MHz，处理1 s湖试数据所需时间为0.032 s，满足实时性要求。

4 结　语

本文采用FPGA平台来实现互质阵列压缩感知算法的DOA估计，分析了互质阵列下压缩感知的算法原理和硬件实现流程，介绍了系统整体开发环境及实现方案，并作出了相应改进：软件层面选择了正交匹配追踪算法来解决凸优化问题，同时由于互质阵列对信号稀疏性可以保证进而对算法进行了一定的简化，提高了系统运算能力；硬件层面PL端与PS端协同配合，针对DDR内存断电不保存的特性，引入了SD卡作为数据的最终存储端，简化了数据导入流程。最终分析了系统运算延时验证处理实时性，并通过Matlab计算仿真与ILA抓取信号的数据对比来协同验证此设计运算结果的正确性。