似P范数特征值分解高分辨率声源定位识别方法研究

刘月婵 ，何元安，商德江, 孙 超

(1.哈尔滨工程大学 水声技术重点实验室，哈尔滨 150001；2.船舶系统工程部，北京 100036)

为了突破常规波束形成(conventional beamforming, CBF)在多声源探测中受阵列物理孔径的限制所导致的“瑞利限”，众多高分辨率谱估计法应运而生[1-2]，其中以子空间多重信号分类(MUSIC)算法为代表。该算法通过对基阵接收数据协方差矩阵进行特征分解，从而得到与信号分量相对应的信号子空间与噪声子空间，利用其正交特性构造出“针状”空间谱峰，相较于以往的各种方法大大的提高了定位精度及空间分辨率。然而该方法仅能得到声源方位的估计结果，无法真实反映声源贡献的相对大小，不能达到声源识别的目的。近场声全息技术[3-6](near field acoustic holography, NAH)作为另一种重要的噪声源定位识别方法虽然能获得声源的真实贡献量，还可对声场进行重建进而获得更丰富的声场信息，然而该方法在水下辐射噪声源定位方面的应用较少。究其原因，水下(航行器辐射噪声通常被视为分布式声源，且多集中在中低频段；噪声源辐射形式未知；噪声源间相关性未知。这些问题导致以往方法获得的定位结果缺乏可靠性，并且很难获得噪声源辐射强度的精确值。

鉴于特征分解类算法的高分辨率特性，本文提出了一种基于似P范数特征分解的噪声源定位识别方法。通过特征分解获得各阶特征子空间及子空间响应函数，对各阶特征子空间采用预设声源的方法建立子空间声源向量重构模型，利用似P范数稀疏性约束[7-9]求解该重构模型，迫使子空间声源向量能更快，更有效的向真实稀疏源位置收敛，获得能量更加集中的稀疏解。仿真及水池试验结果显示该方法可实现对声源位置及幅度的准确估计，空间分辨率高。具有很高的工程应用价值。

1 基本理论

(1)

(2)

由于特征子空间与特征矢量具有相同的正交特性，式(2)可写成：

(3)

式(3)说明阵列对空间声源的响应可表示成各阶特征子空间响应函数的叠加。若Si为第i阶特征子空间对应的i阶声源向量，H为声源向量与接收阵元间的传递函数矩阵，那么每阶特征子空间均可看成是一个网络系统，其输出为子空间响应函数ui：

HSi=ui

(4)

与常规波束形成声源定位识别方法相同，对于任意阶特征子空间，首先将扫描域离散化为N个格点区域，节点处为波束扫描点，各扫描点均存在一预设声源，该预设声源可以是点声源，多极子声源等，或不同类型声源的组合形式，文中选定点声源为预设声源，其它情况将在后续研究中讨论。除目标声源外，其余为假想声源，所有预设声源幅值构成了子空间内N×1维声源向量Si，认为该子空间响应函数ui是由这些预设源辐射的声场叠加产生，不同预设声源类型联合构成了扫描点与阵元间的M×N维传递函数矩阵H。通过求解式(4)可计算出声源向量Si,即该子空间内波束扫描输出。考虑式(4)为声学逆问题，通常采用正则化技术[10-11]获得问题的近似解：

(5)

其中α称为正则化参数[12]。其取值与输入的信噪比有关,通常为矩阵HH-1(或H-1H)最大特征值的0.1%～10%。

由于假想声源并不真实存在，其幅值在向量Si中表现为零或近似为零，而目标声源相较与假想源通常具有较大的幅值，这说明Si具有信息局部化的特点，可将其看成是一组稀疏信号。正则化技术虽然解决了逆问题的不适定性，但并不能获得最优的稀疏解。为了得到能量集中的稀疏解，本文采用基于似P范数的迭代加权算法，将式(4)转化为如下线性约束优化问题：

(6)

利用Lagrange乘子的方法求解上式,最小化的代价函数为：

(7)

解出关于声源向量Si和Lagrange乘子λ的形式解如下：

(8)

(9)

综上得到求解声源向量Si的迭代公式如下(k≥1)：

(10)

2 仿真分析

为了验证本文提出的基于似P范数特征值分解的噪声源定位算法(P-EVD)的有效性,将其声源定位识别结果与经典的常规波束形成技术(CBF),子空间MUSIC算法及基于幅度补偿的MVDR(ac-MVDR)算法[14]进行比较分析。采用由声压水听器构成的7×7均匀网格阵，阵元间距0.3 m，目标声源位于x-y平面上，扫描域 8 m×8 m，扫描间隔0.25 m，基阵面与扫描面相距z=5 m，设水中声速为1 500 m/s，背景噪声为高斯白噪声，除特别说明，信噪比snr=20 dB，Monte-Carlo实验数为50，范数p的取值将在后续研究中重点分析，这里不做分析，仿真中取p=1，正则化参数α取最大特征值的1%。

算例1: 两个非相干单极子声源A和B分别位于(-1,0,0)和(1,0,0)，发射频率为1.5 kHz幅值分别为aA=2，aB=1的单频信号。

图1 非相干单极子声源定位

图1所示是四种算法对非相干单极子声源的定位结果，由于常规波束形成算法的空间分辨率低，并且旁瓣起伏较大(图1(a))，从图中仅能判断出在(-1,0,0)处存在幅度为1.944的声源，不能分辨出两个声源的空间位置,定位失效；与CBF方法相类似，ac-MVDR算法的空间分辨率虽然有所提高，然而对于文中采用的小孔径阵列而言，其性能相较于(b)图的MUSIC算法并无改善；图1(b)中的MUSIC算法是一种高分辨率空间谱估计算法，其空间分辨率明显高于CBF，具有尖锐的聚焦峰和低旁瓣级，图中两声源空间位置清晰可见，该方法可实现对非相干声源的准确定位。但对声源幅值的估计不尽理想(图中两声源幅值为aA=0.055,aB=0.056)，这是由于MUSIC空间谱估计算法是利用信号子空间与噪声子空间的正交性估计目标声源方位，并不能获得目标声源的幅值。由上节的阐述可知本文提出的P-EVD算法是以P范数为约束获得信号能量高度集中的稀疏解，因此在完成对声源精确定位的基础上还可实现对声源幅值的估计，使这一问题得以解决。图1(d)显示的是采用P-EVD算法输出的空间谱，图中无旁瓣，具有更尖锐的聚焦峰尺度，声源位置清晰可见，可获取高精度高分辨率声源定位结果。图中显示两声源幅值为aA=1.991，aB=0.991,估计相对误差为0.67%，估计精度良好，该方法可获得对声源方位及幅值满意的估计结果。

算例2: 考虑两相干单极子声源A和B，幅值aA=aB=2，其他声源参数同算例1。

图2 相干单极子声源定位

图3 分布式声源定位

从图2显示的结果可以看出：① 与算例1的仿真结果相似，CBF的低分辨率使其对两相干声源的扫描结果看起来更像是一个声源，依然无法分辨出声源的个数及空间位置,定位失效。② 虽然MUSIC算法在非相干声源定位中展现出高分辨率的优越性，但对于相干声源，由于无法准确估计声源个数，因此很难实现对相干声源的检测(图3(b))，且幅度估计结果与实际声源相差较大，其定位结果同样失效。③ ac-MVDR具有比CBF和MUSIC算法更低平稳的旁瓣级，然而其空间分辨率不足以区分出两相距较近的相干声源。④ P-EVD算法依然具有优异的空间分辨率和精确的幅度估计能力，aA=aB=1.997，其幅度估计相对误差为0.3%，克服了其他方法不能处理相干声源的不足。算例1同算例2说明P-EVD方法可实现对声源位置及幅度的精确估计而无需考虑其相干性。

算例3: 目标声源为一条有限长线源A和两个非相干单极子声源B和C，线源A与x轴夹角为45°，长度1.4 m，其几何中心位于(-1,1.5,0)，仿真中将其看成是有多个单极子声源构成的分布式声源。声源B与声源C分别位于 (2,0,0)和(-1,1,0)处，声源A，B和C均发射频率为1.5 kHz的单频信号，幅度为aA=aB=aC=1。声源布放情况及扫描结果如图3所示，图中＊标注处为声源，□标注处为测量阵列。

算例3仿真结果显示：① CBF，MUSIC及ac-MVDR算法对分布式声源的扫描输出看起来更像是集中声源，图中不能准确分辨出声源数，声源定位及幅度估计值均存在偏差。MUSIC 与ac-MVDR算法仅仅是消除了CBF中部分旁瓣干扰，依然不能获得精准的定位结果及幅度估计。② P-EVD的扫描结果与图中标定的分布式声源吻合的很好，可同时实现集中式声源及分布式声源的定位识别，输出声源幅度为aA=0.986,aB=0.979，aC=0.996,估计误差1.3%。通过该仿真可知，该方法可实现对分布式声源形态，方向，位置的有效识别及其幅度的准确估计。

算例4：本算例将讨论分析频率及测距阵列孔径比对本文算法(P-EVD)目标方位及声源幅度估计的影响，分别给出各自的均方根误差曲线。频率变化范围100～15 000 Hz，测距阵列孔径比值0.1～30，每个分析频率及测距孔径比的Monte-Carlo实验数为50。

图4 P-EVD性能分析

图4给出了P-EVD声源定位性能受分析频率及测距孔径影响的曲线。通过对比分析可以得到如下结论：① P-EVD算法克服了常规波束形成低频受限的缺点，可实现在低，中，高全频段内声源定位，并不受测量距离及阵列孔径的约束，定位准确有效；② 相对于P-EVD算法对声源方位估计的宽容性，其声源幅度估计仅在中高频段内有效，且测距孔径比需小于15方可获得满意的幅度估计结果。

由于文中采用的是迭代求解计算，且扫描域较大，因此计算速率受到一定程度的制约，表1给出了上述算例中各类算法的计算时间。多次迭代导致P-EVD相对于其他两种方法计算耗时有所增长，但并不冗长，综合其在目标方位和幅度估计方面优良性能，该方法依然具有显著的优势。

表1 计算时间(单位:秒)

3 水池试验验证

为了验证本文提出算法的可行性，于消声水池中开展了水下噪声源定位识别验证试验，如图5所示。由上节对P-EVD算法的性能分析结果，确定试验中各项参数的选取。试验中声源，水听器阵列及扫描面布放情况如图6所示。其中扫描面定义为吊放声源所在的平面，与池底垂直，扫描面大小设定为-4 m≤x≤4 m，-4 m≤y≤4 m，扫描间隔0.25 m。测量面为水听器阵列所在的平面，与扫描面平行，测量间距z=3 m，考虑到布放条件及水听器个数的限制，测量系统由10元垂直线阵构成,阵元间距0.2 m，该阵列沿x方向每间隔0.2 m进行一次采样构成10×10个测点的测量平面。吊放声源采用球型换能器作为单极子声源，声源频率f=2 kHz，采样率fs=32.768 kHz，p=1。

图5 消声水池及扫描系统

图6 声源定位试验示意图

图7 两非相干单极子声源试验结果

工况1：两非相干单极子声源A和B入水深度为4 m，x方向间距1 m,其在扫描面上的坐标分布为A：(-0.5,0,0)和B：(0.5,0,0)，声源B输出信号幅值是声源A的2倍，试验结果以声源A为基准做归一化处理，结果如图7所示。

图8 两相干单极子声源试验结果

图7为非相干声源试验处理结果。试验结果显示P-EVD具有比常规波束形成算法更高的空间分辨率，并且在声源相距较近时其空间分辨能力要优于子空间MUSIC算法及ac-MVDR算法。图中所示的声源幅度估计的相对归一化结果验证了本文算法可同时实现目标方位及声源幅度满意的识别结果。

工况2:两相干单极子声源A和B入水深度为4 m，x方向间距1.6 m,其在扫描面上的坐标分布为A: (-0.8,0,0)和B: (0.8,0,0)，两声源输出信号幅值相同，试验结果以声源A为基准做归一化处理,结果如图8所示。

图8中利用似P范数特征值分解(P-EVD)算法对两相干单极子声源的试验结果显示，与其他方法相比,定位旁瓣得到了明显的抑制,分辨率高,两声源幅值估计结果近似相等. 试验结果表明P-EVD克服了其他方法不能识别相干声源的不足，计算过程中可无需考虑声源间相干性。

4 结 论

提出了一种基于似P范数特征值分解的声源定位识别方法。该方法通过在扫描点预设不同的声源类型，认为阵列接受信号由所有预设声源辐射声场叠加产生，建立声源重构模型。利用似P范数获得最优稀疏解，进而得到高分辨率，背景更为平滑的声源定位效果。计算机仿真结果及水池试验结果显示该方法有效消除了传统算法分辨率不足，存在栅瓣效应以及相干声源处理能力不佳的问题，具有良好的背景噪声抑制能力，可有效实现对单极子声源及分布式声源的准确定位及幅度估计而不受声源间相干性的影响。

[1]Harry L, Van T著.汤俊译. 最优阵列处理技术[M]. 北京：清华大学出版社，2008.

[2]陈志菲，孙进才，侯宏. 宽带DOA估计的类MUSIC波束形成算法[J]. 电子学报，2011，39(6)：1257-1260.

CHENG Zhi-fei，SUN Jin-cai，HOU Yong. MUSIC-Like beamforming method for wideband DOA estimation[J]. Acta Electronica Sinica，2011，39(6)：1257-1260.

[3]杨德森，陈欢，时胜国. 基于最大似然估计的辐射噪声源近场定位方法性能分析[J]. 振动与冲击，2010，29(12)：67-70.

YANG De-sen，CHEN Huan，SHI Sheng-guo. Orientation measure of radiant noise source at short distance based on maximum likelihood estimation[J]. Journal of Vibration and Shock，2010，29(12)：67-70.

[4]徐亮，毕传兴，陈剑，等. 基于波叠加法的patch近场声全息及其实验研究[J]. 物理学报，2007，56(5)：2776-2783.

XU Liang，BI Chuan-xing，CHEN Jian, et al. Algorithm and experimental investigation of patch nearfield acoustic holography based on wave superposition approach[J]. Acta Physica Sinica，2007，56(5)：2776-2783.

[5]张小正，毕传兴，徐亮，等. 基于波叠加法的近场声全息空间分辨率增强方法[J]. 物理学报，2010，59(8)：5564-5571.

ZHANG Xiao-zheng，BI Chuan-xing，XU Liang, et al. Resolution enhancement of nearfield acoustic holography by the wave superposition approach[J].Acta Physica Sinica, 2010，59(8)：5564-5571.

[6]Christophe L，Manuel M，Alexandre G. Boundary element method for the acoustic characterization of a machine in bounded noisy environment[J]. Journal of the Acoustial Socity of America，2007，121(5)：2750-2757.

[7]Christophe L，Manuel M，Alexandre G. Measurement of confined acoustic sources using near-field acoustic holography[J]. Journal of the Acoustial Socity of America，2009，126(3)：1250-1256.

[8]Figueiredo M A T，Nowak R D. A bound optimization approach to wavelet-based image disconsolation[J]. IEEE Trans Image Processing，2005，2(2)：782-787.

[9]Figueiredo M A T，Bioucas-Dias J M，Nowak R D. Majorization-Minimization algorithms for wavelet-based image restoration[J]. IEEE Trans Image Processing，2007，16(12)：2980-2991.

[10]Chartrand R，Wotao Y. Iteratively reweighted algorithms for compressive sensing[J]. IEEE Trans Acoust，Speech，Signal Processing，2008，34(6)：3869-3872.

[11]Nelson P A，Yooh S H. Estimation of acoustic source strength by inverse methods: Part1 conditioning of the inverse problem[J]. Journal of Sound and Vibration，2000，233(4)：639-664.

[12]Williams E G. Regularization method for near-field acoustic holography[J]. Journal of the Acoustical Society of America. 2001，110(4)：1976-1988.

[13]Candes E J，Wakin M B，Boyd S P. Enhancing sparsity by reweighted l1 minimization[J]. Journal of Fourier Analysis and Application. 2008，14(5-6)：877-905.

[14]陈欢，何良，杨德森，等.基于幅度补偿的MVDR水下噪声源近场定位识别方法研究[J]. 振动与冲击，2012，31(2)：51-54.

CHEN Huan，HE Liang，YANG De-sen，et al. Underwater noise sources identification in near-field location based on MVDR method with amplitude compensation [J]. Journal of Vibration and Shock，2012，31(2)：51-54.