反卷积波束形成技术在水声阵列中的应用

2020-07-27孙大军马超梅继丹杨宛珊魏秋雨

哈尔滨工程大学学报 2020年6期

孙大军，马超，梅继丹，杨宛珊，魏秋雨

(1.哈尔滨工程大学水声技术重点实验室，黑龙江哈尔滨 150001；2.海洋信息获取与安全工信部重点实验室(哈尔滨工程大学) 工业和信息化部，黑龙江哈尔滨 150001；3. 哈尔滨工程大学水声工程学院，黑龙江哈尔滨 150001)

阵列信号处理是现代信号处理领域的重要分支，被广泛的应用在雷达、声呐中，高分辨波束形成技术作为提高声呐目标探测性能的重要手段，几十年来一直是国内外学者的研究热点。但是由于声呐阵列规模的逐渐增大和本艇的机动等特殊性的应用环境，使得理想的阵列流型和长时间信号平稳性难以保证，导致了许多高分辨阵列信号处理方法的实际应用效果并不理想。近年来具有高稳健性和实时性(小快拍、短样本)特征的高分辨阵列处理方法已成为了阵列信号处理中新的研究热点。

由于经典的常规波束形成算法受到空域傅里叶变换极限分辨能力的限制，如何突破瑞利限曾一度成为广大学者的重要研究方向。自Capon波束形成技术被提出以来，产生了一批高分辨空间谱估计方法，比如谐波分析法[1]、最大熵法[2]、最小方差无畸变响应法(minimum variance distortionless response, MVDR)[3]、基于ARMA (autoregressive moving average model)的线性预测法，MUSIC (multiple signal classification)[4]算法，ESPRIT (estimating signal parameter via rotational invariance techniques)[5]算法等。由于MVDR和MUSIC以及ESPRIT这类经典的高分辨处理方法无法分辨相干源，之后又出现了一系列的解相干处理方法。空间平滑技术其基本思想是利用均匀线阵的平移不变性，将阵列划分为多个重叠子阵，通过对所有子阵协方差矩阵的平均来解相干。但在水声信号处理领域，由于经典类高分辨处理算法对阵列流形精度要求高、低信噪比性能下降严重等问题，实用性一直受到质疑。学者们又针对具有良好干扰抑制效果的MVDR算法开展了稳健性算法研究，典型的方法有线性约束最小方差算法[6]、对角加载方法[7]，基于特征空间波束形成等方法[8]。但是根据经典的阵列信号处理理论，稳健性、高分辨与处理增益往往不可兼得，高分辨方法稳健性差，而为了提高稳健性对权值进行对角加载这类方法往往又会导致实际处理增益下降，现有方法大多是在这三者之间寻找平衡点，所以每种方法都会有应用局限性。

到了20世纪，人们开始利用其他学科的成果来解决阵处理中的高分辨问题，凸优化理论也被应用到阵列信号处理中，Eldar等[9]运用凸优化理论，将波束形成问题转化为二阶锥优化问题，突破了传统波束形成方法中存在的未能求取全局最优解的局限。随后，压缩感知理论被推广于阵列信号处理领域[10]，相比于前几种高分辨算法，该方法的优点是可应用于小快拍数据，无需预先估计声源个数，可应用于阵元数少于声源数的情形等，但受噪声功率估计误差影响较大，且计算量大，使工程应用受限。

Yang[11-12]将反卷积算法与波束形成相结合，提出一种高增益稳健反卷积波束形成技术，并应用到声压阵和圆阵的处理中，对反卷积波束形成的稳健性、分辨率和增益方面进行了验证，改善效果显著。但文献中所提出的基于Richardson-Lucy(RL)的反卷积波束形成方法只适用于PSF具有移不变特性的阵列中。针对上述问题，文献[13]针对矢量阵PSF的移变性，将非负最小二乘(nonnegative least-squares, NNLS)和DAMAS (deconvolution approach for the mapping of acoustic sources)算法应用到矢量阵中，并提出了一种改进的R-L算法[13-14]，该方法在低信噪比的情况下与NNLS和DAMAS 相比，具有更窄的主瓣宽度，更低的旁瓣和更高的稳健性。多目标情况下的常规波束形成空间谱输出可以看成每一个角度的指向性图和该角度的源强度乘积之和，在数学上可以用一个叠加积分来表示。如果阵列的PSF具有移不变性，此时叠加积分可以被表示为卷积过程。因此，对常规波束输出空间谱和自然指向性函数反卷积可以得到目标函数。

在实际阵列信号处理中，由于应用场景和需求的不同，阵列的形式多种多样，除了等间距直线阵、圆阵、矢量阵等典型阵列形式外，不等间距阵、十字阵、面阵等阵列形式也是常采用的阵列形式。本文主要分析了水声常用阵列结构PSF的特点，将多种典型阵列按照卷积模型类型进行了分类综述，将其划分为波束图移变阵列和波束图移不变阵列。然后分别总结了移不变和移变模型阵列的反卷积波束形成典型求解方法，以期为反卷积技术在水声阵列信号处理中的应用提供一定的指导。

1 移不变点扩散函数的反卷积波束形成问题

1.1 移不变的点扩散函数定义

在空间域处理中，有一大类阵列的自然指向性函数直接或者通过简单的变换就能够具有PSF移不变的特点。在这里移不变的定义就是指向某个角度的波束输出是和角度无关的，其可以看做自然指向性函数的移位，即R(θ|ϑ)=R(θ-ϑ)。或者在某一个变换域具有波束输出具有移不变性，不同位置的波束输出可以视指向性函数的移位。这种移不变性可能是在一维，也可能是在二维或高维空间。

1.2 典型移不变模型阵列

1.2.1 圆阵

对于某些阵列，PSF是移不变条件是可以直接满足的，例如，一个等间距圆阵，如果仅考虑二维平面的情况下其自然指向性函数为:

(1)

式中：N为阵元数；r为圆阵半径；λ为波长。则其在角度ϑ的指向性函数可表示为：

(2)

当阵元数N满足N>4πr/λ+2时，指向性函数可以近似为：

(3)

其中，J0(x)为零阶贝塞尔函数：

(4)

此时，圆阵ϑ方向的指向性函数R(θ|ϑ)是θ-ϑ的函数，即R(θ|ϑ)=R(θ-ϑ)，近似具有关于方位ϑ的移不变性质。

常规波束形成功率谱P(θ)可以看成每一个角度的指向性函数R(θ|ϑ)和该角度的信源强度S(ϑ)乘积之和，在数学上可以用叠加积分表示：

(5)

由于圆阵的指向性函数具有移不变特性，波束输出叠加积分的过程可以等效为卷积过程：

(6)

窄带信号频率3 kHz，圆阵的阵元数20，半径0.5 m，改变信号的方向，源1的方位为0°，源2为30°，源3为60°。指向不同角度的圆阵指向性函数如图1所示，可以看出圆阵的指向性函数是自然指向性函数的圆周移位，因此是移不变的。

图1 指向不同角度的圆阵的指向性函数Fig.1 The PSF of a circular array steering to different angles

.2.2 均匀声压线列阵

另外还有一类阵列的指向性函数是需要进行一定的映射变换才能够满足移不变性质的，例如均匀等间距声压线阵、非均匀线阵。以均匀声压线阵为例，均匀线阵几何模型图如图2所示。其自然指向性函数为：

图2 均匀线阵模型Fig.2 Array geometry for uniform line array

(7)

其在角度ϑ方向的指向性函数为：

(8)

这里取平方是为了和CBF的功率谱对应，显然R(θ|ϑ)≠R(θ-ϑ)。因此，叠加过程不能直接用卷积过程代替。但是经过简单的变换，单位线阵的自然指向性函数具有关于cosθ移不变的性质，即有R(cosθ|cos ϑ)=R(cosθ-cos ϑ)。

因此均匀线阵波束输出可以等效为卷积过程：

R(cosθ)*S(cosθ)

(9)

窄带信号，信号频率500 Hz，直线阵的阵元数20，阵元间距半波间距，信号来波方位的余弦值分别为0、0.25和0.5。指向不同角度的声压线阵的指向性函数如图3所示，可以看出指向不同角度的指向性函数和角度无关，具有移不变性，可以看做自然指向性函数的移位。

图3 指向不同角度的声压线阵的指向性函数Fig.3 The PSF of ULA steering to different angles

1.2.3 非均匀声压线列阵

假设N个声压传感器非均匀的分布在同一条直线上，信号从阵列的远场入射，如图4所示。

图4 非均匀线阵模型Fig.4 Array geometry for non-uniformly linear array

假设入射信号是单频信号，第1个阵元H1接收到的信号是Acos(2πft)，把1号阵元当做参考阵元，那么第i个传感器接收到的信号为：

si(t)=Acos(2πft+φi)，i=1,2,…，N

(10)

其中φi为声程差HiPi导致的相位差:

(11)

将阵元接收到的信号相加，整个阵列的输出为：

(12)

将s(t)取平方并归一化得到的指向性函数为：

(13)

相似地，其在角度ϑ方向的指向性函数为：

(14)

显然，非均匀线阵的指向性函数是sinθ-sin ϑ的函数，因此非均匀线阵的PSF是关于角度的正弦值的移不变函数，R(θ|ϑ)=R(sinθ-sin ϑ)。非均匀线阵常规波束形成功率可以写做如下卷积过程:

R(sinθ)*S(sinθ)

(15)

为了说明非均匀线阵移不变性，以20元非均匀阵为例，各阵元的位置分别选取半波间距阵列的[1,2,3,4,5,6,8,10,12,14,16,18,20,22,24,26,27,30,31,35]号阵元构成非等间距阵列。信号是500 Hz单频信号，其中λ对应信号频率的波长。

图5 指向不同角度的非均匀声压线阵的指向性函数Fig.5 The PSF of NLA steering to different angles

1.2.3 远场平面阵

一个任意分布的N元平面阵，如图6所示放置于xoy平面，假设信号从(θ0,φ0)方位以远场平面波入射，N个阵元位置坐标分别为[xn,yn]，n=1,2,…,N。则第n个阵元与参考阵元之间的声程差为：

图6 任意面阵模型Fig.6 Array geometry for planar array

βn=k(xnsinφcosθ+ynsinφsinθ)

(16)

其中，k为波数，k=ω/c=2π/λ，λ为信号的波长，定义ux=sinφcosθ，uy=sinφsinθ，那么基阵的阵列流形向量为a(k)=(e(-jk(uxx1+uyy1)),…,e(-jk(uxxN+uyyN)))T，常规波束形成器的加权向量为wc=a(θ0,φ0)/M。因此常规波束形成的指向性函数为：

Rp(ux,ux0,uy,uy0)=|wHa(θ,φ)|2=

Rp(ux-ux0,uy-uy0)

(17)

故面阵的波束可以看做是一个关于ux=sinφ·cosθ和uy=sinφsinθ的二维移不变函数，且ux和uy都在[-1,1]范围内，因此可以将矩形面阵归为PSF移不变阵的二维反卷积求解问题。以矩形阵为例说明平面阵的二维反卷积问题，矩形阵仿真条件：X轴阵元数8，Y轴阵元数8，信号形式为频率500 Hz单频信号，阵元间距为1.5 m。

图7是矩形面阵在不同角度所对应的单目标指向性函数。由图可见，2个指向性函数形状相同，只是进行了平移，因此面阵的指向性函数可以看成是移不变的。

图7 矩形阵指向不同角度的三维指向性函数Fig.7 PSF of rectangle array steering to different source angles

.3 PSF移不变模型反卷积求解方法

对于移不变的点扩散函数，常用的卷积方法有CLEAN[15]、DAMAS[16]、NNLS[17]、R-L[18-19]算法等。在其基础上又衍生出适用于相干信号的CLEAN-C算法，基于傅里叶变换的快速DAMAS2[20]和DAMAS3[20]算法，基于稀疏信号重构的SC-DAMAS[21]算法，基于小波变换计算网格的快速DAMAS算法[22]，基于傅里叶变换的FFT-NNLS算法等。Klaus在文献[17]中比较了DAMAS，NNLS,基于傅里叶变化的DAMAS2, FFT-NNLS和R-L算法，结果表明DAMAS和NNLS不要求PSF是移不变的，而DAMAS2,FFT-NNLS和RL算法适用于PSF移不变的反卷积问题中，且RL算法和DAMAS,NNLS相比具有更低的计算量。文献[23]将CLEAN算法，DAMAS,NNLS和R-L算法应用到三维空气声源定位中，结果表明几种方法都能显著提高声源分辨能力，但在处理相干源时R-L的效果最好。文献[24]分析比较了DAMAS,NNLS和ISCA 3种反卷积算法，结果表明DAMAS的定位误差精度最差， ISCA有最高的定位精度和最高的稳健性。在这里简要介绍几种常用的反卷积算法。

1.3.1 DAMAS

将系统的输入输出描述成线性方程组的形式：

P=RS+N

(18)

式中：P代表常规波束形成的输出；R代表基阵的指向性函数PSF；S代表源的分布函数；N代表噪声。DAMAS是在忽略噪声的条件下用高斯赛德尔方法求解线性方程组。对于线性移不变的PSF，矩阵R是Toeplitz 矩阵且是一个半正定对称阵。

高斯赛德尔方法通过引入正约束来解决式(18)，迭代过程：

(19)

(20)

1.3.2 NNLS算法

最小二乘问题，其基本原理可以推广到阵列反卷积处理，具体方法是在常规波束输出，矢量阵点扩散函数字典，声源目标函数之间建立差函数方程组，通过最小化差函数的原则来实现对目标函数S的求解，获取声源的分布，矩阵为：

φ=‖RS-P‖2

(21)

式中：P是矢量阵常规波束形成输出空间谱；R是所有角度矢量阵的指向性函数形成的指向性函数字典矩阵；S是反映声源方位和强度信息的目标函数。对于上述方程组求解可以采用梯度投影法，该方法的核心是负梯度方向指向标量场下降最快的方向，通过在φ关于S的负梯度方向上按特定步长反复迭代搜索来求得S的结果。计算的具体步骤可以参考文献[13]。

NNLS算法并不局限于移不变反卷积模型求解，只需要获得阵列所有方位的指向性函数矩阵即可。需要注意的是，对于移不变模型采用NNLS算法时可采用FFT变换域计算的方式实现快速计算，但这并不适用于移变模型阵列。该方法并不局限于矢量阵。对结构固定不变的阵列该方法都适用。

1.3.3 R-L算法

R-L算法是一种基于贝叶斯理论的迭代方法，该方法被广泛用于图像复原中。R-L算法中要求PSF是移不变的，在阵列信号处理中，就是指基阵的波束图是和角度无关的(指向某一个角度的波束图都可以看做自然指向性函数的移位)，R-L算法的迭代方程为：

R(sinθ)*Si(sinθ)

(22)

(23)

dsinθ≤ε

(24)

1.3.4 仿真对比

以均匀分布声压线阵为例，分析阵列的PSF是移不变的反卷积问题求解。由式(9)可知，ULA的波束输出为：

P(cosθ)=R(cosθ)*S(cosθ)

(25)

基于阵列的参数可以得到阵列的自然指向性函数，对接收到的数据做常规波束形成可以得到波束输出的方位谱，对常规波束输出空间谱P(θ)和自然指向性函数R(θ)反卷积可以得到目标函数S(θ)，其理想情况下是δ函数，必然是对目标方位的高分辨估计结果。关于PSF是移不变的反卷积问题求解引用文献[14]中所采用的R-L方法。

为了说明反卷积波束形成算法在PSF是移不变的阵列的应用效果，以圆阵为例进行仿真分析。

仿真条件同图1，为理想无噪声条件下，双目标的仿真结果图。图8说明常规波束输出的空间谱是阵列的PSF和信号的目标强度函数的卷积。与CBF相比，反卷积波束形成方法有更窄的主瓣和更低的旁瓣。

图8 双目标圆阵仿真结果Fig.8 Simulation results of two target for circular array

2 移变点扩散函数的反卷积波束形成问题

2.1 移变的点扩散函数

在空间域处理中，有一大类阵列的自然指向性函数具有PSF移变的特点，移变的定义就是指向某个角度的波束输出是和角度有关的，即R(θ|ϑ)≠R(θ-ϑ)，不可以看做自然指向性函数的简单移位。但与一般的时域系统相比，空间域上的移变问题有一个典型特点：在水声中一般是将具有固定拓扑结构指的阵列去采集声信号，对于任意形状的阵列，当阵列拓扑结构固定后，PSF无论是移不变还是移变都是可以预先测定得到的，如果采用测量获得的PSF函数由于其包含了对阵型误差等因素等影响的校正会有更好的反卷积求解效果。

2.1.1 矢量阵

以均匀矢量线阵为例来说明这一问题，角度定义同图1。与声压阵相似，直线声矢量阵的波束输出可以写为[16]：

(26)

式中：Rv(θ|ϑ)矢量阵指向ϑ的指向性函数；S(ϑ)是以方位为变量的目标信源强度分布函数；θ和ϑ分别代表扫描角度和目标信号方位。

矢量阵在90°方向指向性函数为：

(27)

式中：Rp(θ)为声压阵的自然指向性函数，可以看出，矢量阵的自然指向性函数Rv(θ)为声压阵自然指向性函数Rp(θ)与一个角度指向因子u(θ)=(1+cosθ)/2的乘积。

指向ϑ方向的矢量阵指向性函数为：

Rv(θ|ϑ)=

(28)

比较Rv(θ)和Rv(θ|ϑ)，不难发现Rv(θ|ϑ)≠Rv(θ-ϑ)，也不可通过简单的正弦，余弦变换将其变为移不变函数。因为，u(θ)=(1+cosθ)/2是关于θ移不变的函数，而Rp(θ)是关于cosθ的移不变函数，二者相乘后无法通过简单的变量代换将其等效为某个统一变量的移不变模型。

矢量阵仿真：20元均匀分布矢量阵，阵元间距3.75 m，为信号频率半波间距，信号频率200 Hz，采样频率20 kHz。改变信号的方向，指向不同角度的矢量阵的指向性函数如图9所示，可以看出不同角度的矢量阵的指向性函数是移变的，不是自然指向性函数的移位。而且尽管矢量有一定的抗左右模糊作用，但是其效果与目标方位有关，在端射方向附近效果仍然较差，会在对称方向形成伪峰。

图9 指向不同角度的矢量阵指向性函数Fig.9 PSF of VSLA steering to different source angles

2.1.2 共形阵

共形阵是传感器附着在非平面载体表面上的阵列，简称共形阵。由于其具有良好的空气动力学特性，可以充分利用导流罩的纵向空间，增加基阵有效孔径，共形阵正被当做一种新型布阵方式应用在新型声呐系统中。

假设各阵元的位置pn是随机分布，如图10所示，第nth个传感器的位置坐标为pn=[xn,yn,zn]。K个平面波从阵列的远场入射，简单起见，假设K=1。θ和φ分别是方位角和俯仰角，定义信号传播方向的单位向量可以表示为：u=-[cosθcosφ,sinθcosφ,sinφ]T，把原点当做参考点时，第nth个传感器相对于参考点的时延差为：

图10 任意分布阵列模型Fig.10 Arbitrary conformal array geometry

(29)

式中c为声速。

若入射信号为Acos(2πft)，那么第nth个传感器接收到的信号为：

sn(t)=Acos[2πf(t+τn)]

(30)

将接收到的信号求和，并将幅度归一化并且取平方，得到阵列的指向性函数为：

相似的，其在角度ϑ方向的指向性函数为：

R(θ,φ)=

(32)

显然，共形阵的指向性函数不具有规律可言，是移变函数。为了对该问题简要说明，这里采用半椭圆阵列来说明共形阵反卷积波束形成问题。对于共形阵声呐(非成像声呐)，探测目标为阵所在平面内的远距离目标，即只考虑目标方位方向上的一维求解问题，此时是一维移变反卷积问题。

将41个传感器如图11所示，放置在椭圆轨迹上，椭圆方程是 4y2+x2=100; 每个传感器在x轴的位置是xn,其中xn=[0:0.5:10]∪[9.5:-0.5:0]。

图11 椭圆形共形阵模型Fig.11 Oval conformal array geometry

信号频率为1 kHz，改变信号的方向，指向不同角度的共形阵的指向性函数如图12所示，可以看出不同角度的共形阵的指向性函数是移变的，不是自然指向性函数的移位。

图12 指向不同角度的共形阵指向性函数Fig.12 PSF of conformal array steering to different source angles

2.2 PSF是移变模型的反卷积求解方法

文献[14]针对矢量阵分析比较了3种反卷积波束形成方法，DAMAS，NNLS和改进的R-L算法。改进的R-L算法在低信噪比的情况下与NNLS和DAMAS 相比，具有更窄的主瓣宽度，更低的旁瓣和更高的稳健性。因此针对空间域反卷积问题点扩散函数移变情况的特点——移变但可预知的特点，应用文献[14]中提出的改进R-L反卷积算法。