一种基于混合时频域的旋转声源识别方法与流程

本发明属于声场识别技术领域，具体涉及一种基于混合时频域的旋转声源识别方法。

背景技术：

对以一定角速度旋转的声源，时域追踪das与反卷积方法组合形成混合时频域方法，可以高分辨率成像结果。该方法的第一步是：利用时域追踪识别声源；第二步是：基于频域反卷积方法清晰化声源识别结果。但是，混合时频域方法需要高采样频率实现过采样，以实现准确的时间延迟。当采样频率较低时，时间延迟计算不精确，则时域das波束形成结果不精确，后续反卷积声源识别结果的空间收敛性欠佳。

2015年，lylloff等提出了一种新颖的提高效率的声源识别反卷积算法，即基于傅里叶变换的快速迭代收缩阈值算法(fastiterativeshrinkage-thresholdingalgorithm,fista)，该方法相比于基于傅里叶变换的nnls空间收敛性更好、计算效率更高。但是，基于傅里叶变换的快速迭代收缩阈值算法仍需要较高的采样频率。

技术实现要素：

针对现有技术存在的技术问题，本发明所要解决的技术问题就是提供一种基于混合时频域的旋转声源识别方法，它能降低时域采样频率，且在较低采样频率下，本发明的空间收敛、旁瓣抑制及计算效率的综合性能比先有技术更优。

本发明所要解决的技术问题是通过这样的技术方案实现的，它包括：

步骤1、使用das识别声源

假设t时刻旋转声源在第m个麦克风位置rm处，m＝1,2,3,…,m为麦克风的序号，产生的时域声压为p(rm,t)；

声源扫描面上第j个扫描网格点时域追踪das波束形成输出bj(t)的时间信号为：

式中，j＝1,2,3,…,j，j是扫描网格点的总数，rj(t)是t时刻第j个扫描网格点的位置矢量，rmj(t)表示t时刻第m个麦克风到第j个扫描网格点的距离，c是声速；tmj(rm,rj(t))表示t时刻第m个麦克风和第j个扫描网格点之间的传递函数；

步骤2、基于频域反卷积方法清晰化声源识别结果

波束形成输出和点传播函数之间的关系式为：

b＝aq

subjecttoqj≥0

式中，b是时间间隔上所有扫描网格点输出的自谱向量，q为未知的声源强度列向量，qj是向量q的第j个元素；a是阵列点传播函数矩阵；

将波束形成输出和点传播函数之间的关系式转化为无约束最优化式：

式中，||·||2表示l2范数，||·||1表示l1范数，λ为促进声源稀疏的非负的正则化参数，表示声源强度

选择正则化参数λ＝min{0.001,0.001||b||∞}，其中||·||∞表示无穷范数；

求声源强度q的迭代流程如下：

初始迭代设置有，声源强度q0＝0，辅助向量y1＝q0，步长τ1＝1，

第一步，使用梯度信息并将结果在非负象限进行投影，第l步迭代的zl如下：

zl＝p+(yl-a^t(ayl-b)/l)

式中，zl表示未经稀疏促进的源强向量，p+表示在非负象限上的欧几里得投影，yl表示经过第l次迭代的辅助向量，l表示lipschitz常数，l等于a^ta的最大特征值；

第二步，使用收缩算子更新声源强度ql如下：

ql＝sign(zl)max{|zl|-λ,0}

第三步，更新中间迭代步长τl+1如下：

第四步，更新辅助向量yl+1如下：

以上四步为一个循环，将yl+1代入第一步求得zl进入下一次循环，如此反复迭代直到迭代次数或者终止条件满足时，得到声源强度q。

本发明的技术效果是：

在fista算法基础上，将声源稀疏正则化引入声源识别，在采用较低的时域采样频率就能达到现有混合时域方法在高时域采样频率下的声源识别性能。

附图说明

本发明的附图说明如下：

图1为仿真布置示意图；

图2为参考采样频率下不同反卷积方法仿真结果；

图3为不同采样频率下、不同反卷积方法仿真结果；

图4为仿真下不同采样频率下、不同反卷积方法声源主瓣点数、积分及峰值列表；

图5为试验布置图；

图6为两个采样频率下、不同反卷积方法试验结果；

图7为试验下两个采样频率下、不同反卷积方法声源主瓣点数、积分及峰值列表。

具体实施方式

下面结合附图和实施例对本发明作进一步说明：

本发明包括以下步骤：

步骤1、使用时域追踪延时求和(delayandsum,das)识别声源

假设t时刻旋转声源在第m个麦克风位置rm处(m＝1,2,3,…,m是麦克风的序号)产生的时域声压为p(rm,t)。

声源扫描面上第j个扫描网格点时域追踪das波束形成输出bj(t)的时间信号为：

式(1)中，j＝1,2,3,…,j，j是扫描网格点的总数，rj(t)是t时刻第j个扫描网格点的位置矢量，rmj(t)表示t时刻第m个麦克风到第j个扫描网格点的距离，c是声速；tmj(rm,rj(t))表示t时刻第m个麦克风和第j个扫描网格点之间的传递函数。

步骤2、基于频域反卷积方法清晰化声源识别结果

为了使用反卷积方法清晰化时域追踪das声源识别结果，选取一段时间间隔内所有扫描网格点输出的自谱向量b(ω)(bj(t)表示在时域下波束形成输出，b(ω)为bj(t)在频域下取自谱得到)作为传统波束形成输出，将该自谱赋值给该时间间隔中间时刻所对应的位置，以该位置计算点传播函数并进行后续反卷积计算。

反卷积方法通过在传统波束形成输出中移除阵列点传播函数的影响获得真实声源强度。传统波束形成输出和点传播函数之间构建如下线性方程组：

式(2)中，b是时间间隔上所有扫描网格点输出的自谱向量，q为未知的声源强度列向量，qj是向量q的第j个元素；a是阵列点传播函数矩阵。

根据参考文献becka,teboullem.afastiterativeshrinkage-thresholdingalgorithmforlinearinverseproblems[j].siamjournalonimagingsciences,2009,2(1):183-202.(becka,teboullem.一种求解线性逆问题的快速迭代收缩阈值算法[j].siamjournalonimagingsciences,2009,2(1):183-202.)：

线性方程(2)可以转化为如下无约束最优化问题：

式(3)中，||·||2表示l2范数，||·||1表示l1范数，λ为促进声源稀疏的非负的正则化参数，表示声源强度；

当λ＝0，未使用稀疏约束，没有收缩算子的fista声源识别反卷积方法不使用收缩算子(为方便区分，标记为fistaλ＝0)，通过快速迭代求解式(3)以获得声源强度。

当λ≠0为完整的fista声源识别反卷积方法，通过选择合适的λ促使声源稀疏以提高声源空间收敛性和抑制旁瓣，本发明选择正则化参数λ＝min{0.001,0.001||b||∞}，其中||·||∞表示无穷范数。

本发明的完整fista声源识别反卷积方法采用迭代求声源强度q，根据参考文献：becka,teboullem.afastiterativeshrinkage-thresholdingalgorithmforlinearinverseproblems[j].siamjournalonimagingsciences,2009,2(1):183-202(becka,teboullem.一种求解线性逆问题的快速迭代收缩阈值算法[j].siamjournalonimagingsciences,2009,2(1):183-202.)，迭代流程如下：

初始迭代设置有，声源强度q0＝0，辅助向量y1＝q0，步长τ1＝1，

第一步，使用梯度信息并将结果在非负象限进行投影，第l步迭代的zl如下：

zl＝p+(yl-a^t(ayl-b)/l)(4)

式(4)中，zl表示未经稀疏促进的源强向量，p+表示在非负象限上的欧几里得投影，yl表示经过第l次迭代的辅助向量，l表示lipschitz常数，l等于a^ta的最大特征值，a是阵列点传播函数矩阵。

第二步，使用收缩算子更新声源强度ql如下：

ql＝sign(zl)max{|zl|-λ,0}(5)

第三步，更新中间迭代步长τl+1如下：

第四步，更新辅助向量yl+1如下：

以上四步为一个循环，将yl+1代入第一步求得zl进入下一次循环，如此反复迭代直到迭代次数或者终止条件满足时，得到声源强度q。

仿真模拟试验

如图1所示，仿真设置如下：双声源绕z轴逆时针旋转，初始坐标分别为(0.5,0)，(-0.5,0)，它们分别被标记为source1、source2。声源频率为2khz，声源的旋转速度为300r/min，两声源强度同为100db。扫描平面为1.5m×1.5m，扫描网格点间距为0.025m×0.025m。采用半径为0.65m的b&k36通道扇形轮阵列，阵列平行于声源面并距离声源面1m。仿真信号中加入20db高斯白噪声，该仿真中均进行1000次迭代，显示动态范围设置为15db。

声源识别过程：1、根据假定旋转目标声源信息计算阵列上各传声器接收的声压信号；2、通过时域追踪das得到初步结果；3、进行频域反卷积计算清晰化结果；4、将输出结果转化成声压级并进行成像。在频域反卷积的方法中，选取了damas、nnls、fistaλ＝0与本方法发明对比。

damas通过高斯赛徳尔迭代求解；nnls利用最小化残差平方和来求解；fistaλ＝0的正则化参数λ＝0，通过快速迭代求解；本方法发明是一种加速梯度投影算法，首先利用梯度下降信息，并在非负象限进行投影，算法更具鲁棒性；其次，本发明使用收缩算子，是典型的稀疏促进的求解方法；最后，使用新的步长和辅助向量，本发明有更高的计算效率。

首先，混合时频域方法采样频率设置为1024khz，并标记为参考。以相同时间长度的时域追踪das波束形成结果变换获得自谱，频率分辨率为16hz。在参考采样频率下混合时频域方法的声源识别结果见图2，从图2看出：damas、nnls、fistaλ＝0和本发明，四种反卷积方法均能够有效识别声源位置，空间分辨率佳，幅值精度damas最高，本发明次之，fistaλ＝0和nnls稍差。

实际中采样频率难以达到参考采样频率1024khz，在采样频率降低分别至16.384khz、65.536khz的声源识别结果见图3，从图3看出：在16.384khz、65.536khz两种较低采样频率下，本发明收敛性最优；nnls和fistaλ＝0收敛性较差，主瓣较宽；damas对低采样频率导致的误差敏感，声源识别结果不精确，在空间上扩散且偏离真实声源位置。

上述四种反卷积方法在有限次的迭代下，主瓣峰值未能收敛到真实声源强度，声源强度扩散到声源的相邻区域，分别获得两个声源主瓣区域覆盖点数、主瓣区域源强积分及峰值，见图4列表中的数值。扩散区域源强积分同真实声源强度非常接近，表明上述方法均能够准确的识别声源总强度。本发明的峰值收敛性在低采样下与damas相当，在高采样频率下较damas稍差，本发明和damas的峰值收敛性明显高于nnls和fistaλ＝0。所有采样频率下，本发明主瓣区域覆盖点数稍优于damas，且明显比nnls及fistaλ＝0少。

综上分析，在低时域采样频率下，本方法发明可达到damas的高时域采样频率的声源识别性能，补偿现有混合时频域方法对高采样频率的要求。在相同较低采样频率下，本方法发明的计算效率、空间收敛、幅值收敛的综合性能较damas、nnls、fistaλ＝0的混合时频域方法更优。

试验验证

试验布局如图5所示，将两个直径为0.025m的同频率不等强度的压电式蜂鸣器安装在旋转臂的两端，强度相差约10db，强源和弱源分别标记为s1和s2，蜂鸣器发声频率为3.15khz，两声源旋转半径为0.15m，旋转臂由可调速电机驱动，该次试验中电机逆时针转动，转速为375r/min，转速传感器置于旋转臂后方获得转速信息。采用半径为0.65m的b&k36通道扇形轮阵列，阵列平行于声源面并距离声源面1m。计算平面被设置为0.3m×0.3m，网格间距0.01m×0.01m。迭代次数1000次，显示动态范围设置为15db。

两个采样频率分别为16.384khz、65.536khz。以相同时间长的数据获得自谱，频率分辨率均为16hz，成像频带范围为δf＝3136-3168hz。两个采样频率下不同反卷积方法声源识别结果如图6所示：两个采样频率下，damas在低采样频率16.384khz下声源识别结果空间收敛性欠佳，存在旁瓣污染，主瓣区域成扩散状，声源识别结果随着采样频率提高空间收敛性提高，旁瓣被部分抑制，但依旧未能完全收敛到真实声源位置。nnls和fistaλ＝0均能有效识别声源，但主瓣较宽且存在少量旁瓣。本发明在两种采样频率下都能够准确的识别声源。

图7列出了所有方法的声源主瓣点数、积分及峰值。在两种采样频率下，本发明与nnls、fistaλ＝0相比，声源识别性能最优，主瓣区域聚焦点数最少，峰值最高，强源主瓣区域积分相当，弱源区域积分稍差。而damas先呈现扩散状，后部分收敛。

试验与仿真结果一致，证明了仿真结果的可靠性。