于涤非， 张春华， 黄勇

(中国科学院大学, 北京 100190)



改进型圆面阵的快速波束形成算法设计

于涤非， 张春华， 黄勇

(中国科学院大学, 北京 100190)

为降低基于圆面阵的三维成像声纳波束形成算法的运算量，提出改进圆面阵布阵阵型设计及优化波束方位的方法。改进型圆面阵由多条通过圆心的线阵组成，其波束形成等效于所有线阵波束形成结果的叠加；通过对波束方位进行优化，使得每一条线阵的波束形成过程等效为等长度的一维快速傅里叶变换运算。仿真结果表明，改进型圆面阵的波束形成运算量较直接计算显著降低。运用改进型圆面阵可显著降低三维成像声纳信号处理系统的复杂度,提高了声纳的综合性能。

声学; 三维成像声纳; 圆面阵; 波束形成; 快速傅里叶变换

0 引言

近年来，利用平面换能器阵通过波束形成对水下目标进行三维成像技术发展迅速。实时三维成像声纳[1-6]常用的平面基阵为方阵[1-3]、同心圆面阵[2-6]等。

基于方阵的三维成像声纳的优点是存在大量快速波束形成算法，能够运用快速傅里叶变换(FFT)[2]或Chirp z变换(CZT)[7]等方法降低运算量。可以采用编程方便灵活的数字信号处理器(DSP)实现其信号处理系统[8]。但也存在一定的缺点，在基阵阵元数一定的情况下，声纳的视场范围与测向精度是一对矛盾。为提高测向精度，需要增加阵元间距以扩大基阵孔径，而阵元间距的增加会缩小视场范围。如果同时提高测向精度和视场范围，则必须大量增加基阵阵元数，从而增加系统复杂度。

基于同心圆面阵的三维成像声纳在一定程度上可以解决视场范围与测向精度间的矛盾。同心圆面阵阵元间距大于半波长时，其波束方向图上不会出现栅瓣[4]，因此在阵元数一定的条件下，声纳可以同时获得较大的视场范围与较高的测向精度。但基于同心圆面阵的三维成像声纳同样有一定的缺点，文献[9]指出圆面阵的波束形成可等效为广义傅里叶变换，运算量较大，且无法通过FFT、CZT等方法降低运算量，因此对于基于同心圆面阵的三维成像声纳，其信号处理系统较为复杂，一般只能采用图形处理器(GPU)[10]或现场可编程门阵列(FPGA)[11]实现。其实现难度较大，灵活性较差。

因此，三维成像声纳设计时很难兼顾视场范围、测向精度、运算量、阵元数这4方面的性能，一般都需要牺牲某些方面的性能以保证主要指标满足设计要求。为解决上述问题，本文提出一种改进型圆面阵设计，基于该圆面阵的三维成像声纳能够兼顾这4方面的性能，从而提高声纳整体性能。

1 基于同心圆面阵三维成像声纳的特点

假设布阵间距为d，声纳发射单频脉冲信号，波长为λ，使用如图1所示的俯仰角θ及方位角φ定义。

图1 声纳方位角及俯仰角定义Fig.1 Definition of azimuth and elevatioan angle of sonar

文献[1]中采用半波长布阵的正方形栅格圆面阵，视场范围为θ∈[0°，90°]，φ∈[0°，360°]，阵元数为3 136，角分辨率为1.6°. 文献[2]中采用波长布阵的方阵，角分辨率为1°，阵元数为2 304，视场范围为θ∈[0°,30°],φ∈[0°,360°]. 波长布阵相比于半波长布阵可以获得更高的测向精度，但视场范围有所降低。

同心圆面阵相比于方阵的优点是可以在阵元数基本不变的条件下，同时提高角分辨率及视场范围。

对比文献[2]中的方阵，如果采用阵元数基本相同的同心圆面阵，则在波长布阵时，其角分辨率为0.93°，视场范围为θ∈[0°，90°]，φ∈[0°，360°].

同心圆面阵上，所有的阵元均分布在同心圆环上，相邻圆环的半径差为布阵间距d；每个圆环上，阵元均匀分布，阵元间距同样为d. 一个具有2 362阵元，布阵间距等于波长的同心圆面阵阵型如图2所示。

图2 同心圆面阵阵型Fig.2 Concentric circular array

远场条件下，基于图2的三维成像声纳对位于θ=30°、φ=180°的点目标成像的波束方向图在sinφ=0切片上与方阵的结果对比如图3所示。

图3 同心圆面阵与方阵波束方向图对比Fig.3 Comparison of beam pattern of concentric circular array and square array

方阵的波束方向图在θ=30°、φ=0°出现栅瓣，而同心圆面阵的结果没有出现栅瓣，位于θ=30°、φ=0°的波束高度为-15.21 dB. 因此同心圆面阵的视场范围可以达到θ∈[0°，90°]，φ∈[0°，360°]. 表1为3种阵型的比较。

表1 不同阵型性能比较Tab.1 Comparison of performances of different arrays

图4为同心圆面阵在视场范围内对该点目标成像的波束方向图。

图4 同心圆面阵对点目标成像的波束方向图Fig.4 Beam pattern of point target imaged by concentric circular array

方阵[2]要同时达到角分辨率为1°，视场范围为θ∈[0°,90°],φ∈[0°,360°]，则基阵阵元数需要增加至9 216阵元。这不仅增加了系统复杂度，同时由于布阵间距减小，基阵阵子的直径也同时减小，其内阻增加[12]，会导致数据采集系统输出信号的信噪比降低。

基于同心圆面阵的三维成像声纳则可以在阵元数一定的条件下，同时提高角分辨率及视场范围。

但同心圆面阵的波束形成无法采用FFT运算降低运算量，其波束形成实现一般只能采用编程复杂、灵活性较差的GPU或FPGA，不适合使用DSP. 如果能对同心圆面阵进行改进，使得改进型圆面阵可以使用FFT算法降低运算量，则可以在提高三维成像声纳性能的同时降低算法的实现难度。

2 圆面阵阵型设计

方阵、六角形阵等阵型之所以存在频域FFT波束形成[9]，是因为这两种阵型均可以通过线阵拓扑得到。如果设计一种圆面阵，其阵型同样可以通过线阵拓扑得到，则该圆面阵也可以通过FFT运算降低运算量。

文献[13]中给出一种可以通过线阵拓扑得到的圆面阵，组成圆面阵的每一条线阵为均匀线阵且均不是稀疏的。如果按照这种方式构造圆面阵，则在基阵孔径一定的情况下，基阵只能集成很少的阵元，阵增益较低，不适合三维成像声纳。对这种设计方式进行改进，使设计的圆面阵满足三维成像声纳的应用需求。

三维成像声纳设计时，一般会根据性能指标要求、加工工艺水平、电路规模等限制条件确定声纳基阵能够集成的阵元数及阵元间距。因此改进型圆面阵的设计目标是在给定阵元间距d以及给定基阵的期望阵元数Ke条件下，设计一种圆面阵满足：

1) 圆面阵上任意相邻两个阵元间距不小于d；

2) 圆面阵的阵元数尽可能接近Ke；

3) 圆面阵由同心圆环阵构成，同时可由线阵拓扑而成；

4) 尽可能降低基阵孔径，即阵元间距尽可能接近d.

基于如上要求设计的圆面阵，可以在阵元数一定的条件下尽可能减小基阵面积，从而降低声纳的体积。

为达到设计要求，按照如下方式构造面阵：

1) 假设圆面阵的层数M=3(中心阵元不算做一层)，即圆面阵由M个同心圆环阵组成，第m层圆环阵半径为md；

kmax(M)=3602arcsin(12M)

2) 计算第M层圆环阵阵元间距为d时能容纳的最多阵元数

，其中⎣」表示向下取整。初始化搜索过程中的最外层阵元数变化量Θ=0；

kM(Θ)=2(kmax(M)2-Θ)

αm=2arcsin(12m)α

4) 令m=m-1，计算第m层圆环阵相邻两个阵元方向向量夹角

Km(Θ)=360αm

α，「⎤表示向上取整，即第m层圆环阵相邻两个阵元方向向量夹角必然为α的整数倍。第m层圆阵的阵元数为

；

Θ=kmax(M)22

6) 令Θ=Θ+1,回到第3步，统计并记录K(M,Θ)，直至；

7) 求取K(M,Θ)的最大值K(M,Θmax)，则K(M,Θmax)即为当圆面阵层数为M时，能够容纳的最多阵元数；

8) 如果K(M,Θmax)≥Ke，则M即为圆面阵层数，Θmax对应的kM(Θ)即为圆面阵最外层阵元数，该圆面阵已满足设计要求，搜索过程结束。如果K(M,Θmax)

按照以上步骤构造的圆面阵，是阵元数量大于Ke且最接近Ke的圆面阵，且在该数量下做到了孔径最小，图5为搜索圆面阵层数流程图。

图5 圆面阵层数搜索流程图Fig.5 Flow chart of layer number searching of improved circular array

图6为计算层数为M的圆面阵能容纳阵元数量最大值的算法流程图。

图6 层数为M的圆面阵最大阵元数搜索流程图Fig.6 Flow chart of maximum sensor number searching of M-layer improved circular array

构造过程实际是两个嵌套的穷举过程，首先从圆面阵可能的最小层数开始，计算圆面阵能够包含的最多阵元数。如果阵元数大于期望阵元数Ke，则可以确定使得圆面阵阵元数大于Ke的最小层数。

在计算层数为M的圆面阵能够容纳阵元数的最大值时，首先根据阵元间距计算第M层最多阵元数，穷举第M层阵元数为大于最多阵元数一半的所有偶数。计算对应的圆面阵能够包含的阵元数，阵元数最大值即为层数为M的圆面阵至多能包含的阵元数。这样构造的圆面阵，由层数及最外层阵元数唯一确定。

改进型圆面阵上，位于同一圆环阵上的相邻两个阵元方向向量夹角均为α的整数倍。假设圆面阵中心阵元位于圆面阵坐标系原点，则圆面阵可以看做由L=KM(Θmax)/2条穿过坐标原点的线阵组成。因此该圆面阵可以分解为多个圆环阵，相邻圆阵的半径差为d，也可以分解为多个稀疏线阵，线阵相邻两个阵元间距必然为d的整数倍。

图7 改进型圆面阵实例Fig.7 Example of improved circular array

图7为Ke=48的圆面阵。该圆面阵包含53个阵元，共4层，可以分解为4个等间距圆环阵，同时可以分解为12条线阵，每条线阵相邻阵元间距为d的整数倍。

3 波束方向图波束分布设计

图7所示的圆面阵，是由12条线阵等角度旋转拓扑而成。因此该圆面阵在波束形成时，等效为12条线阵波束形成结果的线性叠加。球坐标系下，层数为M，最外层阵元数为2L的改进型圆面阵在发射单频矩形信号时[2]的远场波束形成公式为

B(u)=

(1)

式中：ωln为第l条线阵上、第n个阵元的权重，如果该位置没有阵元，则权重为0；Sln为该阵元收到的信号经解调后得到的复信号；u=(ux,uy,uz)=(sinθcosφ, sinθsinφ, cosθ)为波束的单位方向矢量。

由三角函数公式，(1)式可化为

(2)

令

(3)

(3) 式表示第l条线阵的波束形成过程。因此改进型圆面阵的波束形成过程，等效为L条线阵波束形成结果的线性叠加。

由于组成改进型圆面阵的每一条线阵均是对一个阵元间距为d的均匀线阵的稀疏抽取，因此只要合理地选取波束方向[9]，就可以使用一维FFT运算来计算每一条线阵的波束方向图，从而降低计算量。

由于组成圆面阵的线阵不是平行的，而是等角度旋转而形成的，如果波束方向选取不合理，会造成不同线阵在波束形成时，其FFT的长度不同，FFT长度越长，计算量越大。为获得最低计算量，需要设计一种波束方向图的波束分布方式，满足如下两个条件：

1)可以使用一次一维FFT运算来计算每一条线阵的波束形成过程；

2)组成圆面阵的每一条线阵波束形成时，所对应的一维FFT运算长度相等。

对于(3)式，假如存在整数k、Nl，使得

(4)

则(3)式可化为

因此当波束分布满足(4)式的约束时，就可以采用一次一维FFT运算优化每一条线阵的波束形成过程，从而满足条件1.

条件2则要求：N0=N1=…=NL-1=N，即对于所有线阵，其FFT运算的长度均为N. 由(4)式可知，Nl的取值与cos(φ-lα)的值域集合相关。

当l取不同值时，如果cos(φ-lα)的值域集合完全相同，则必然有Nl=N. 因此为使波束分布满足条件2，必须对φ的选取进行优化。

若波束方向图上每一个波束，其对应的φ满足：

此时有Nl=N，l=0，2，…，L-1，则(2)式化为

(6)

p、H均为整数，因此当波束分布满足(5)式的约束时，对于任意l，一维FFT运算的长度相等。

除满足(4)式和(5)式外，波束分布同样需要满足视场范围及测向精度的要求。假设需要计算的波束方向图视场角为θ∈[0，θmax]，角分辨率为γ，则可以按照如下步骤，设计波束分布方式：

I=sinθmaxsinγ

1)所有波束分布在I个等间距同心圆环上，即sinθ=0∶sinγ∶sinθmax，则

2)在第i个圆环上，波束为均匀分布，数量为Ji，i=0，1，2，…，I-1，i=0表示波束位于波束方向图中心，对应θ=0°；

3)计算θ=θmax圆环上夹角为γ的相邻波束方位角φ的差值βγ，图8为示意图。

图8 方位角差值计算示意图Fig.8 Sketch of azimuth difference calculation

单位向量OA、OB为两个俯仰角均为θ=θmax的波束方向向量，两向量的夹角为角分辨率γ，向量OC、OD为OA、OB在Oxy平面内的投影向量，则向量OC、OD的夹角即为向量OA、OB所对应方位角的差值。

4)计算β：

β=ααβγ=αkα

，

当βγ<α时，α必为β的整数倍，当βγ≥α时，β=α,因此kα必为正整数；

βi，βi=2arcsin(sin(γ/2)sinθi)ββ=kiβ

Ji=360βi

5)计算位于第i个圆环上相邻波束方位角φ的差值,即每一个波束对应的方位角φ必然为β的整数倍。该圆上包含的波束数量为，则位于第i个圆环上，第j个波束的方向矢量为u(i,j)=(sinθi·cos(jkiβ),sinθisin(jkiβ),cosθi),i=0,1,2,…,I-1,j=0,1,2,…,Ji-1.

至此，确定波束方向图上所有需要计算的波束的方向矢量。

图9为对应图7圆面阵波束方向图的波束分布，视场范围θ∈[0°,8°]，φ∈[0°,360°],角分辨率γ=2°.

图9 优化的波束分布图Fig.9 Optimized beam distribution

4 圆面阵远场频域FFT波束形成

当波束分布满足(4)式和(5)式时，(6)式成立，此时波束形成过程等效为L次N点一维FFT运算结果的线性叠加。

由波束分布设计过程可以看到(5)式是可以严格满足的，而(4)式中的等号很难严格成立，但必然存在一个正整数N,使得在某一精度要求下有

(7)

则

(8)

假设波束方向图共包含Nb个波束，则运算量为L次N点FFT，(L-1)Nb次复数加法，LN次复数乘法。

N=2log2λdsinγsinβ

假设圆面阵阵元间距d=λ，期望阵元数为2 304,则根据图5所示流程搜索到的圆面阵具有阵元T=2 314，由M=30个同心圆环阵组成。同时该圆面阵可分解为L=75条线阵，如图10所示。

图10 2 314阵元圆面阵阵元分布图Fig.10 Distribution map of array elements of improved circular array

假定计算的波束方向图视场范围θ∈[0°,90°]，φ∈[0°,360°],角分辨率γ=1°，则该波束方向图包含Nb=14 222个波束，FFT长度N=8 192.

由于每条线阵原始数据长度为2M+1=61≪8 192，因此可以利用Partial FFT简化运算，将一次8 192点Partial FFT分解为128次64点FFT.

相比于直接计算，采用频域FFT波束形成，减少了92.5%的复数乘法以及85.6%的复数加法。

表2为图10所示圆面阵使用(1)式直接波束形成及FFT频域波束形成运算量比较。

表2 运算量比较Tab.2 Comparison of operation quantity

图11为图10所示圆面阵，在视场范围θ∈[0°,90°]，φ∈[0°,360°]内，对俯仰角为0°方位点目标成像时的波束方向图，主瓣宽度0.963 8°，主副瓣比19.06 dB.

图11 对θ=0°方位点目标成像时的波束方向图Fig.11 Beam pattern of θ=0° point target imaged by improved circular array

由于(8)式中，B(k)与B(u(i,j))并非严格相等，因此在对数波束方向图上，FFT波束形成算法计算的BP(k)与通过(1)式直接计算的BP(u(i,j))之间存在一定的误差。该误差与FFT长度N的取值有关，N越大，误差越小，计算量越大。在声纳设计时，需要根据实际情况取得误差与计算量之间的折衷，以使声纳获得最佳的综合性能。

考察图11中，所有波束高度高于-35 dB的波束与通过(1)式直接计算结果之间的误差σ(i,j)：σ(i,j)=|BP(u(i,j))-BP(k)|，BP(u(i,j))>-35 dB,BP(k)>-35 dB.

图12为σ(i,j)的统计直方图,平均误差为0.016 dB，最高误差为0.058 dB.

图12 波束高度大于-35 dB的波束误差统计Fig.12 Error statistics of beam higher than -35 dB

表3为N取不同值时，误差与计算量的比较。

5 改进型圆面阵栅瓣分析

改进型圆面阵能够在布阵间距等于波长的条件下，对空间点目标成像不出现假目标[4-5]，是因为组成圆面阵的线阵是旋转的，而不是平行的。

表3 N取不同值时误差与计算量的比较Tab.3 Comparison of operation quantity and errors for different N

组成圆面阵的每一条线阵的布阵间距为波长，因此对于每一条线阵而言，其对点目标成像必然会出现栅瓣。但由于线阵是通过旋转组成圆面阵的，因此所有线阵形成的栅瓣能量经叠加后均匀分布在圆周上，而不是位于一点，因此不会出现假目标。

图13为对位于θ=30°、φ=180°的点目标成像的波束方向图。

图13 θ=30°、φ=180°方位点目标成像的波束方向图Fig.13 Beam pattern of θ=30°， φ=180° point target

图13中点目标所形成的栅瓣能量近似均匀分布在圆周上，在平面直角坐标系下，该圆周的方程为

(x-x0)2+(y-y0)2=R2，

图14 改进型圆面阵与方阵波束方向图对比Fig.14 Comparison of beam patterns of improved circular array and square array

方阵波束方向图上，在θ=30°、φ=0°方位出现栅瓣，而改进型圆面阵在该方位附近波束的最大高度为-21.9 dB.

6 结论

为降低采用圆面阵的三维成像声纳的信号处理系统复杂度，提出改进型圆面阵阵型设计及优化波束方位的方法，使得改进型圆面阵的波束形成可等效为多次一维FFT运算结果的线性叠加。

仿真结果表明：

1)相比于传统同心圆面阵的波束形成算法，改进型圆面阵的波束形成计算量大幅度降低。

2)采用改进型圆面阵的三维成像声纳的信号处理系统复杂度较低，易于实现，使得声纳可以同时兼顾视场范围、测向精度、运算量、阵元数4个方面，具有良好的综合性能。

References)

[1] 陈鹏. 相控阵三维成像声纳系统的稀疏阵及波束形成算法研究[D]. 杭州:浙江大学, 2009. CHEN Peng. Research on sparse array and beamforming algorithm for phase array three-dimension imaging sonar system[D]. Hangzhou:Zhejiang University, 2009.(in Chinese)

[2] 王朋. 基于稀疏布阵的三维成像声纳信号处理算法研究[D]. 北京:中国科学院大学, 2015. WANG Peng. Research on signal processing algorithm of three-dimensional acoustical imaging sonar based on sparse planar array[D]. Beijing: University of Chinese Academy of Sciences, 2015. (in Chinese)

[3] 袁龙涛. 相控阵三维摄像声纳系统信号处理关键技术研究[D]. 杭州:浙江大学, 2013. YUAN Long-tao. Research on key technologies of signal processing for phased array three-dimensional imaging sonar system[D]. Hangzhou:Zhejiang University, 2013.(in Chinese)

[4] Tekes C, Karaman M, Degertekin F L. Optimizing circular ring arrays for forward-looking IVUS imaging[J]. IEEE Transactions on Ultrasonics Ferroelectrics & Frequency Control, 2011, 58(12): 2596-2607.

[5] Chen K S, Chen H, Wang L,et al. Modified real GA for synthesis of sparse planar circular arrays[J]. IEEE Antennas and Wireless Propagation Letters, 2016, 15: 274-277.

[6] 杨洪亮. 基于均匀圆阵的测向算法研究[D]. 西安:西安电子科技大学, 2014. YANG Hong-liang. Research on direction finding algorithm based on uniform circular array[D]. Xi’an:Xidian University, 2014.(in Chinese)

[7] Palmese M,Trucco A. Three-dimensional acoustic imaging by chirp zeta transform digital beamforming[J]. IEEE Transaction on Instrumentation and Measurement, 2009, 58(7): 2080-2086.

[8] 王朋, 张扬帆, 黄勇, 等. 基于稀疏布阵的实时三维成像声纳系统[J]. 仪器仪表学报, 2016, 37(4): 843-851. WANG Peng, ZHANG Yang-fan, HUANG Yong, et al. Real-time 3D acoustical imaging sonar system based on sparse planar array[J]. Chinese Journal of Scientific Instrument, 2016, 37(4): 843-851.(in Chinese)

[9] 陈丰. Fourier变换在平面、柱、球阵波束形成和DOA估计中的应用[D]. 杭州:浙江大学, 2006. CHEN Feng. The application of Fourier transform to planar,cylindrical and spherical array beamforming and DOA estimation[D]. Hangzhou:Zhejiang University, 2006.(in Chinese)

[10] 印明明, 刘平香. GPU实现水下三维成像研究[J]. 声学技术, 2014, 33(5): 53-57. YIN Ming-ming, LIU Ping-xiang. Study on 3D image sonar based on GPU[J]. Technical Acoustics, 2014, 33(5): 53-57.(in Chinese)

[11] 胡将. 基于Kintex-7的三维声学成像主信号处理系统硬件设计[D]. 杭州:浙江大学, 2016: 17-47. HU Jiang. Hardware design of three-dimensional acoustic imaging main signal processing system based on Kintex-7 [D]. Hangzhou:Zhejiang University, 2016: 17-47.(in Chinese)

[12] 张沛霖, 张仲渊. 压电测量[M]. 北京: 国防工业出版社, 1983: 63-67. ZHANG Pei-lin, ZHANG Zhong-yuan. Piezoelectric measurement[M]. Beijing:National Defense Industry Press, 1983: 63-67.(in Chinese)

[13] 蒋毅. 天线阵列阵元位置优化方法研究[D]. 哈尔滨:哈尔滨工程大学, 2013. JIANG Yi. Research on the optimum method for element position of antenna arrays[D]. Harbin: Harbin Engineering University, 2013.(in Chinese)

[14] Cooklev T. Generalized and partial FFT[J]. IEICE Transactions on Fundamentals of Electronics, Communications and Computer Sciences, 1995(9):1466-1473.

Design of Fast Beam Forming of Improved Circular Planar Array

YU Di-fei， ZHANG Chun-hua， HUANG Yong

(University of Chinese Academy of Sciences, Beijing 100190, China)

In order to reduce the computational complexity of beam forming of 3D imaging sonar based on circular planar array, a method of improving the array design and optimizing the beam azimuth is proposed. The improved circular array consists of a plurality of linear arrays passing through the center, and the beam forming is equivalent to the superposition of the beam forming of all the linear arrays. By optimizing the beam azimuth, the beam forming process of each linear array is equivalent to a one-dimensional fast Fourier transform operation of equal length. The simulated results show that the beam forming computation of the improved circular array is significantly lower than that of the direct calculation. The use of improved circular array can significantly reduce the complexity of the signal processing system and improve the comprehensive performance of 3D imaging sonar.

acoustics; 3D underwater sonar; circular arrays; beam pattern; fast Fourier transform

2016-08-11

国家自然科学基金项目(11304343)

于涤非(1982—)， 男， 博士研究生。 E-mail: nacassoce@163.com

张春华(1962—)， 男， 研究员， 博士生导师。 E-mail: zch@mail.ioa.ac.cn

TB565+.2; TN911.72

A

1000-1093(2017)05-0959-09

10.3969/j.issn.1000-1093.2017.05.016