一种基于互谱优化的互谱广义逆波束形成方法与流程

本发明涉及一种声源识别算法，主要涉及运用一种基于传声器阵列的互谱优化广义逆波束形成算法，在消除背景噪声的基础上来精确识别定位声源。

背景技术：

现在基于传声器阵列的声源识别方法主要有两种，波束形成与声全息。两者各有优劣势，其中波束形成的中远距离声源识别性能相对于声全息较好，而低频性能相对于声全息较差。

波束形成是近年来快速发展的一种声源识别方法，它基于规则或者非规则传声器阵列(如十字轴和轮型阵列)进行信号采集，然后利用如简单延时求和算法进行声场重构。对比传统方法如主观评价法、选择隔离法等，基于多通道信号处理技术的波束形成具有测量方便、重构声场直观可视化等优点。相比声强测量法，波束形成具有测量效率高和应用环境宽等优点。相对于近年来的发展的声全息方法，波束形成在中高频具有较高的空间分辨率，而且能够实现中远距离的声源重构，并且其能够适应不同的运行工况。由于基于波束形成的噪声源识别方法具有灵活的动态显示范围、干扰能力强、中高频高分辨率高等优点，因此近年受到各国研究人员和学者的高度关注。使得其应用范围也不断地拓展，从航空、声呐等军事领域，逐渐发展到车辆噪声测试等民用领域。与此同时衍生出了各种不同的波束形成算法，如传统的互谱波束形成、反卷积算法以及广义逆波束形成算法，各种不断地推出的各色的波束形成算法各具特色又相互联系，不断推进该方法的应用便利性与实用性。

现有技术中，存在一种基于传声器阵列的声源识别方法，在传统的等效源算法的基础上公开了一种宽带声全息算法。其在靠近声源面较近的等效源平面构造一系列等效源，让这些等效源来替代实际声源分布。通过声源识别理论，建立一个等效源到传声器阵列的声学传递方程。并通过传统的梯度下降算法求解该传递方程，从而迭代求解等效源强度。然后通过该等效源强度，基于声传播过程重构在等效源面与阵列面之间的声源分布。而且该方法在等效源求解过程中，引入一个迭代滤波过程，从而保证迭代过程的收敛效率与精准度。该算法由于引入梯度下降算法和滤波过程，能够提高传统等效源方法在中低频的分辨率。同时该算法也能够实现与其他声全息算法之间的结合或者切换，从而实现在中低频声源的重构。但该算法基于传统的等效源算法重构，其在中低频声源以及较近的测量距离中具有较高的空间分辨率，但是对于高频声源或者较远距离测量，其声源识别的空间分辨率将较低。

技术实现要素：

本发明目的在于提供一种声源重构精度高、鲁棒性较好的声源识别算法，其建构在一般广义逆和互谱波束形成的基础上，具有较高的重建精度，并充分结合声压互谱矩阵的优化，抗噪性能较强。

为实现本发明目的而采用的技术方案是这样的，一种基于互谱优化的互谱广义逆波束形成方法，其特征在于，

搭建一个声源重构系统：

该系统包括位于声源附近的摄像头、传声器阵列、多通道的信号采集器以及后处理计算机。所述摄像头用于采集整个声场的空间影像；所述传声器阵列用于测量声场信号；所述多通道的信号采集器用于测量时域模拟信号，并将其转化成数字信号传输给计算机进行后处理；所述计算机的主要功能主要是存储以及后处理采集到的声场信号，并能高效显示重构声场分布情况；

波束形成过程包括以下步骤：

1)基于广义逆波束形成来计算转向矩阵：

声源面布置n个单位声源，则在传声器阵列中的传声器ri点处的测量得到的声压信号如下：

上式中，为声场中ri处的单位激励在r0处的复数响应，q(ri)为假设的单极子声源强度。式(1)表达的是由声源到传声器阵列之间的传递关系，为后续分析必要，可以将其转化成向量—矩阵形式：

p＝aq(2)

其中a为n×m声场传递矩阵，其元素成分由式(1)决定。p为m维测量声压向量，其元素为对应传感器在单一频率f对应的频域响应。q为n维声源强度列向量，其元素成分代表声源平面中对应点处的强度分布。

为了求解重构声场分布，可以将等式(2)转化为该最小化问题：

其中β为吉鸿诺夫正则化参数，l为n×n正则化对角矩阵。通过迭代算法求解该最小化问题(3)，最终可以得到第(n)步下的转向矩阵的表达式：

w⁽ⁿ⁾＝l⁽ⁿ⁾b⁽ⁿ⁾(4)

式(4)中b⁽ⁿ⁾为n×m矩阵，其行向量可以表达成如下形式：

b(i,:)⁽ⁿ⁺¹⁾＝g(:,i)^half/(g(:,i)^halfg(:,i))(5)

式(5)其主要功能是对传递矩阵进行归一化，从而保证输出的准确性。上式中上标h表示共轭转置，g(:,i)为归一化后的传递矩阵g的第i列向量，其中g的表达式如下：

g＝al⁽ⁿ⁾(6)

式(5)中alf为可逆的m×m维矩阵，其表达式有如下形式：

alf＝(gg^h+βi)^-1(7)

上式(7)中正则化参数β，可以通过l曲线或者gvc法得到。

l⁽ⁿ⁾＝diag(|sp|/||sp||∞)(8)

上式(8)中diag为将向量转化成对角矩阵形式，|*|为取向量的绝对值，||*||∞为求向量的无穷范数。sp为归一化后的广义逆波束形成输出：

sp＝w⁽ⁿ⁾p(9)

在求解波束输出sp的求解过程中也加入如下滤波过程：

式(10)中，(|sp|)max取向量|sp|中的最大元素，其中dk随迭代次数(n)的增加而增加，从而使得滤波门限值随迭代次数增加而逐渐减小：

dk＝0.2+0.5n(11)

求解得到基于广义逆波束形成的一个转向矩阵w⁽ⁿ⁾；

2)对声压互谱矩阵进行优化：

基于传感器阵列测量得到的m维声压向量p，求得一个m×m维声压互谱矩阵，其表达式如下：

c＝pp^h(12)

优化后的声压互谱矩阵形式表达如下：

c优＝c-diag(x)(13)

其中，c优为消除掉对角元素的矩阵，diag(x)为一对角阵，x为噪声分量；

x为的元素x(i,1)被限定在了[0,c(i,i)]之间，这就保持了去对角化元素中的一些优势的同时又能够降低通道自噪声。为得到噪声分量x，可以通过迭代优化算法来优化声压互谱矩阵。假设一系列归一化单位向量uk＝ek(k＝1,2,...,m),则有如下不等式成立：

uk^h(c-diag(x))uk≥0(14)

由上式可以推导出如下线性不等式：

上式中d(k,i)＝|u(i,k)|²(其中u(i,k)为向量uk的第i号元素)，d(k,1)＝uk^hcuk。通过上式可得到一个关于噪声分量x的一个约束条件dx≤d，由于声压互谱矩阵的优化的核心思想就是最小化噪声分量-x，所以可得到一个目标函数minf＝-x。将以上最小化噪声分量过程总结如下，即可以一个线性规划问题：

求解最小化噪声分量-x，并输出优化后的声压互谱矩阵c优。

3)依照所述转向矩阵和优化后的声压互谱矩阵重新构造波束输出：波束输出形式如下：

s优＝w⁽ⁿ⁾c优(w⁽ⁿ⁾)^h(18)

上式中优化后的声压互谱矩阵c优可以替代成如下形式：

s优＝w⁽ⁿ⁾(c-diag(x))(w⁽ⁿ⁾)^h(19)

将上式(19)取对角元素并进行开方，最终得到基于优化互谱矩阵的互谱广义逆波束形成的输出如下：

s输出＝(diag(s优))^1/2(20)。

本发明的技术效果是毋庸置疑的，首先基于传统广义逆波束形成构造转向矩阵。其次是结合线性优化方法对声压互谱矩阵进行优化。最后通过求解得到的转向矩阵和优化的声压互谱矩阵，结合传统的互谱矩阵波束形成思想，构造新的基于优化互谱矩阵的互谱广义逆波束形成算法。该算法在运用广义逆波束形成高分辨率优势的同时，能够充分利用基于互谱优化的互谱波束形成的抗噪性能。

附图说明

图1声源重构示意图；

图2声压互谱矩阵优化流程示意图。

具体实施方式

下面结合实施例对本发明作进一步说明，但不应该理解为本发明上述主题范围仅限于下述实施例。在不脱离本发明上述技术思想的情况下，根据本领域普通技术知识和惯用手段，做出各种替换和变更，均应包括在本发明的保护范围内。

由于本实施例涉及的是一个基于声源识别的波束形成算法，其建立在声源识别系统的基础上，所以首先将简要介绍该声源识别系统。其中图1描述了一个声源重构系统示意图，该系统包括位于声源附近的摄像头、传声器阵列、多通道的信号采集器以及后处理计算机。

所述摄像头用于采集整个声场的空间影像；，通过波束形成算法计算得到的波束输出与空间影响相重合即可以得到整个声场的声源分布情况。摄像头采用一般的数字摄像头即可，可通过一般的usb接口传输到计算机。

所述传声器阵列用于测量声场信号；，其结构包括了声压或者声速传感器、传声器支架以及相关线束。在传声器阵列中，传感器的分布可以是规则化的或者是非规则化一、二、三维空间分布，这个根据实际测量需求而选择。传感器数目的选择和传声器间距的选择则根据需要测量的频率范围而决定。

所述多通道的信号采集器用于测量时域模拟信号，并将其转化成数字信号传输给计算机进行后处理；其应该包括多通道的前置放大器以及多通道的具有高精度的模拟/数字转换器，甚至应该包括相关的滤波器，从而保证输出到计算机的数字信号足够精确。

所述计算机的主要功能主要是存储以及后处理采集到的声场信号，并能高效显示重构声场分布情况；。其中计算机要求的内存512mb以上，建议是2g，cpu主频2.0g以上比较好，这样能够充分保证计算的高效性。

波束形成过程包括以下步骤：

假设以上声源重构系统中包括m个声压传声器，通过多通道信号采集器采集到的时域信号传输并储存到计算机，并将采集得到的时域信号通过快速傅里叶变换转化到频域中。实施例中，波束形成算法主要涉及中高频声源的重构问题，所以我们将得到的阵列频域信号中的某个单一频率f成分所对应的复数信号提取出来，并将其转化为后面算法所要处理的测量声压信号p。

由于该算法是基于广义逆波束形成，互谱矩阵优化以及互谱波束形成所开发。

1)基于广义逆波束形成来计算转向矩阵：

，其中的转向矩阵为归一化之后的转向矩阵，其基本作用于传统互谱波束形成中的转向向量类似。为求解该归一化的转向矩阵，声源面布置n个单位声源，则在传声器阵列中的传声器ri点处的测量得到的声压信号(频域单频成分)如下：

上式中，为声场中ri处的单位激励在r0处的复数响应，q(ri)为假设的单极子声源强度。式(1)表达的是由声源到传声器阵列之间的传递关系，为后续分析必要，可以将其转化成向量—矩阵形式：

p＝aq(2)

其中a为n×m声场传递矩阵，其元素成分由式(1)决定。p为m维测量声压向量，其元素为对应传感器在单一频率f对应的频域响应。q为n维声源强度列向量，其元素成分代表声源平面中对应点处的强度分布。

为了求解重构声场分布，可以将等式(2)转化为该最小化问题：

其中β为吉鸿诺夫正则化参数，l为n×n正则化对角矩阵。通过迭代算法求解该最小化问题(3)，最终可以得到第(n)步下的转向矩阵的表达式：

w⁽ⁿ⁾＝l⁽ⁿ⁾b⁽ⁿ⁾(4)

式(4)中b⁽ⁿ⁾为n×m矩阵，其行向量可以表达成如下形式：

b(i,:)⁽ⁿ⁺¹⁾＝g(:,i)^half/(g(:,i)^halfg(:,i))(5)

式(5)其主要功能是对传递矩阵进行归一化，从而保证输出的准确性。上式中上标h表示共轭转置，g(:,i)为归一化后的传递矩阵g的第i列向量，其中g的表达式如下：

g＝al⁽ⁿ⁾(6)

式(5)中alf为可逆的m×m维矩阵，其表达式有如下形式：

alf＝(gg^h+βi)^-1(7)

上式(7)中正则化参数β，可以通过l曲线或者gvc法得到。本专利中推荐使用的正则化参数β为gg^h最大特征值的倍(其中f为聚焦声源频率)。

式(6)中的正则化矩阵l⁽ⁿ⁾为一重要参数，在迭代过程中随广义逆波束形成的输出的变化而变化，从而保证重构过程的精确性。一般可以将其定义为如下形式：

l⁽ⁿ⁾＝diag(|sp|/||sp||∞)(8)

上式(8)中diag为将向量转化成对角矩阵形式，|*|为取向量的绝对值，||*||∞为求向量的无穷范数。sp为归一化后的广义逆波束形成输出：

sp＝w⁽ⁿ⁾p(9)

在传统的广义逆波束形成输出中，由于噪声的影响，将在非声源方向出现干扰旁瓣。为进一步抑制这些干扰旁瓣的影响，从而实现更为精确声源定位，因此在求解波束输出sp的求解过程中也加入如下滤波过程：

该滤波过程的核心思想主要是通过假定一个随迭代变化门限值，在迭代过程中不断地将波束输出向量中小于该阈值的元素强制置为零，进而抑制噪声的影响。式(10)中(|sp|)max取向量|sp|中的最大元素，其中dk随迭代次数(n)的增加而增加，从而使得滤波门限值随迭代次数增加而逐渐减小：

dk＝0.2+0.5n(11)

值得说明的是，步骤1)最终目的不是为了求解波束输出，而是为了求解得到基于广义逆波束形成的一个转向矩阵w⁽ⁿ⁾，即为求解到如式(2)的转向矩阵。

2)对声压互谱矩阵进行优化：

基于传感器阵列测量得到的m维声压向量p，求得一个m×m维声压互谱矩阵，其表达式如下：

c＝pp^h(12)

传统基于声压互谱矩阵的波束形成为了去除通道自噪声的影响，通常通过消除掉互谱矩阵的对角元素，即令c＝c-(c)对角。由于声压互谱矩阵中的对角元素是同通道的功率谱，同通道中的信号的功率谱可能引入通道自噪声，所以通过去去除对角元素，确实能够提高波束输出的抗噪能力。不过该方法也具有较大的副作用，如导致波束输出谱中会出现负功率，甚至可能导致声源定位出现较大的偏差。当然实际测量中我们当然想要在保持去对角化方法抗噪优势的同时，又能够减小定位偏差等负面情况的影响。针对这种情况，我们基于互谱矩阵的广义逆波束形成输出，对声压互谱矩阵进行优化改进。

假设声压互谱矩阵中的对角元素中包括了有功分量和噪声分量，如果在声压互谱矩阵的计算过程中能够在保存有功分量的基础上剔除噪声分量，则此时就能够充分运用优化后的声压互谱矩阵进行声场重构。根据以上假设，我们可以将优化后的声压互谱矩阵形式表达如下：

c优＝c-diag(x)(13)

所以可通过优化方法剔除噪声分量而保持有功分量。上式类似于声压互谱矩阵去掉对角元素，其中，c为消除掉对角元素的矩阵，diag(x)为一对角阵，其代表噪声分量，是需要在对角元素中去除的部分，该对角阵为其对角元素对应向量x中的元素，例如x为的元素x(i,1)被限定在了[0,c(i,i)]之间，c(i,i)指的是矩阵c第i个对角元素，这就保持了去对角化元素中的一些优势的同时又能够降低通道自噪声。为得到噪声分量x，我们可以通过迭代优化算法来优化声压互谱矩阵。假设一系列归一化单位向量

uk＝ek(k＝1,2,...,m),则有如下不等式成立：

uk^h(c-diag(x))uk≥0(14)

由上式可以推导出如下线性不等式：

上式中d(k,i)＝|u(i,k)|²(其中u(i,k)为向量uk的第i号元素)，d(k,1)＝uk^hcuk。通过上式可得到一个关于噪声分量x的一个约束条件dx≤d，(d为一矩阵，其元素值由d(k,i)＝|u(i,k)|²决定；d为一向量，其元素值由d(k,1)＝uk^hcuk决定；)由于声压互谱矩阵的优化的核心思想就是最小化噪声分量-x，所以可得到一个目标函数minf＝-x。将以上最小化噪声分量过程总结如下，即可以一个线性规划问题：

当然求解线性规划问题(17)的方法很多，我们可以根据不同的需要选择不同的优化方法。最后通过选取一般的线性迭代算法求解最小化噪声分量-x，并输出优化后的声压互谱矩阵c优。

其中将声压互谱矩阵的优化过程总结如图2所示，其主要包括四个过程，初始化、第一循环、第二循环以及第三循环。首先开始进行的初始化过程主要是为了给一些变量赋值。对于优化的声压互谱矩阵c优的初始值我们假定为没有进行任何处理的声压互谱矩阵，并令u为由一系列的归一化单位列向量组成的矩阵。其中b为向量x的幅值因子，x下和x上分别代表向量x在迭代求解过程中的下上限值，主要是为了约束其收敛的范围，保证求解的准确性，其中t1和t2为约束变量，主要是为了提高该优化过程的收敛速度，其中0≤t1<t2≤1，最常规的做法可以将t1令为0且将t2设为1。通过设置t1和t2数值，能够有效地在优化声压互谱矩阵和去对角化方法之间切换。其中迭代步数n1是一个非常重要的参数，其取值随不同的测量环境的不同而不同，理论上其取值越大越好，不过其取值过大也会造成计算效率的降低，所以实际测量中必须要选择可以接受的合理数值。

其中优化过程的第二和第三循环主要是为了计算约束条件中的相关参数，其中第三循环通过u去计算矩阵d。第二循环主要是通过矩阵u和声压互谱矩阵c去计算向量d。其中最重要的属于第一循环，第一循环的开始部分主要是计算第二循环的循环次数，并根据矩阵u的列的维数重新赋予矩阵d以及向量d应有的维度以及初始值。当第二和第三循环完成时，已经计算出新的矩阵d以及向量d。此时第一循环的后半部分将通过常规的线性规划方法，在约束住噪声分量x的上下限的基础上进行优化求解，从而计算出新的噪声分量x。再将噪声分量向量x转化成对角矩阵，然后通过声压互谱矩阵减去该对角矩阵就可以得到优化后的声压互谱矩阵，这相当于在声压互谱矩阵的对角元素中减去噪声分量，从而提高波束输出结果的抗噪性能。当然在第一循环的最后还有一个重要的步奏，就是将上一步得到的矩阵u的列向量与优化后的声压互谱矩阵的特征向量做并集，并将得到并集重新赋值给新矩阵u，这样随着迭代次数的增加，矩阵u的列数将逐渐增加。

完成以上循环过程以后，则可以输出优化后声压互谱矩阵c优。该优化后的声压互谱矩阵可以运用到基于互谱的波束形成，从而保证输出的精准性，以下步骤3)将通过基于广义逆的互谱波束形成进行输出，从而进行声源定位识别。

3)依照所述转向矩阵和优化后的声压互谱矩阵重新构造波束输出：依照传统互谱波束形成理论，其波束输出可以表达成如下形式：

sout＝wcw^h(17)

其中的w为由转向向量组成的转向矩阵，其列向量为转向向量。其中转向向量主要作用是通过将波束聚焦到声源方向，从而自声源方向的波束输出得到聚焦而增加，在非声源方向的波束输出从而受到抑制，这样通过观察聚焦的主瓣方向即可以识别定位声源方向。所以基于以上思想，在得到的基于广义逆波束形成转向矩阵与优化的声压互谱矩阵基础上，可以重新构造一个基于广义逆波束形成的互谱波束形成，其波束输出形式如下：

s优＝w⁽ⁿ⁾c优(w⁽ⁿ⁾)^h(18)

上式中优化后的声压互谱矩阵c优可以替代成如下形式：

s优＝w⁽ⁿ⁾(c-diag(x))(w⁽ⁿ⁾)^h(19)

上式w⁽ⁿ⁾就是式(4)中求得的基于广义逆波束形成求得的转向矩阵，通过该转向矩阵将能够聚焦波束输出方向，从而实现声源识别定位。上式中输出的s优形式上类似于传统互谱波束形成的输出形式，但是由于其转向矩阵经过了广义逆波束形成进行了优化，进一步提高了该算法的重建精度。同时由于引入了优化后的声压互谱矩阵，进而提高了该声源重构方法的抗噪性能。为简化上式的输出，将上式(19)取对角元素并进行开方，最终得到基于优化互谱矩阵的互谱广义逆波束形成的输出如下：

s输出＝(diag(s优))^1/2(20)。