一种适用于球麦克风阵列的超分辨声源定位方法与流程

本发明涉及一种适用于球麦克风阵列的超分辨声源定位方法，属于声源定位技术领域。

背景技术：

与自由空间中的信号源定位问题不同，房间混响环境下存在信号多径传播和系统欠定问题(散射体个数大于麦克风数目)，因此经典环境下的doa估计方法在高混响环境下存在性能下降。为了解决上述问题，rafaely等人针对球麦克风阵列，提出基于直接路径检测的声源定位方法(见参考文献[1]o.nadiriandb.rafaely.localizationofmultiplespeakersunderhighreverberationusingasphericalmicrophonearrayandthedirect-pathdominancetest[j].ieee/acmtransactionsonaudio,speech,andlanguageprocessing,22(10),2014，见文献[2]liormadmoniandboazrafaely.directionofarrivalestimationforreverberantspeechbasedonenhanceddecompositionofthedirectsound[j].ieeejournalofselectedtopicsinsignalprocessing,2018:1-1)，引入频率平滑以缓解信号多径传播的影响，同时利用语音信号的时频稀疏特性，选择直接路径信号占主体的局部时频块，使用现有子空间类(music，espreit等)方法估计声源位置，从而解决系统欠定问题。但是，频率平滑只能缓解信号多径传播导致的自相关矩阵秩损问题，并不能保证其主特征值仅对应直接路径分量(见文献[2])，可能存在信号强相关甚至相干的情况。而传统的子空间类方法在强相关、强噪声、短时信号等条件下性能迅速下降。因此，现有的直接路径检测声源定位方法对恶劣环境(强噪声，强混响)下的多声源定位问题性能有待进一步提升。

技术实现要素：

本发明所解决的技术问题是，克服现有方法所遇到的局限，提出一种适用于球麦克风阵列的超分辨声源定位方法，用于解决现有直接路径检测声源定位方法在恶劣环境下性能差的问题。

为了实现上述目的，本发明采用的技术方案为：

一种适用于球麦克风阵列的超分辨声源定位方法，包括如下步骤：

步骤一、球阵采集高阶声场进行球谐波分解和短时傅里叶变换，计算局部时频块的自相关矩阵，并进行频率平滑；

步骤二、将步骤一获得的局部时频块自相关矩阵，代入半定规划问题求解，获得恢复后的球谐波域协方差矩阵；

步骤三、利用步骤二求解得到的球谐波域协方差矩阵作直接路径检测，即对求解得到的球谐波域协方差矩阵作奇异值分解，若最大奇异值与次大奇异值的比值大于阈值，则认为该时频块通过检测；

步骤四、对通过检测的时频块计算空间谱，即利用步骤三中球谐波域协方差矩阵大奇异值对应奇异值矢量作为信号子空间计算空间谱，并将各时频块空间谱叠加，搜索叠加后空间谱的谱峰，获得声源位置。

具体的，为介绍本发明中所述具体步骤，首先给出如下设定：

考虑一个半径为r由i个全向阵元组成的球阵，球阵的第i个阵元位于ri＝(r,φi)，其中其中θ和分别为仰角和方位角，接收信号表示为x(t)＝[x1(t),…,xi(t)]^t，t表示时间。

本发明中步骤一具体步骤为：对接收信号x(t)作短时傅里叶变换，得到接收信号的时频序列x(v,l)，其中v，l分别为时间和频率点索引。再对x(v,l)作球谐波变换(见参考文献[3]b.rafaely,“analysisanddesignofsphericalmicrophonearrays,”ieeetransactionsonspeechandaudioprocessing,,vol.13,no.1,pp.135-143,2005.)：

式(1)中diag[·]表示对角矩阵，κl＝λ/2π，λl表示频点l对应的波长。

hn和jn分别是二类球汉克尔函数和一类球贝塞尔函数，h′n和j′n分别是hn和jn的导数。w为与球阵结构有关的固定加权系数矩阵(见参考文献[3])，为(n+1)²×i维矩阵，其第i列表示为:

式(3)中[·]^t表示矩阵转置，为球谐波函数

为连带勒让德多项式，n为对球阵采集信号作球谐波分解的最高球谐波阶数。将时频序列p(v,l)划分为若干局部时频块。令第v个时间、第l个频率上的局部时频块的包含qv和ql个时间和频率点，计算该局部时频块上的信号自相关矩阵：

本发明中步骤二具体步骤为：步骤一中式(5)计算得到的代入下列半正定规划问题求解，

式(6)中tr(·)表示矩阵的迹，(·)^h表示矩阵的共轭转置，re(·)表示取实部，对于矩阵m，m≥0表示矩阵m半正定，||·||2表示欧几里得范数，i为单位阵，表示第v个时间、第l个频率上的局部时频块自相关矩阵，令为kronecker积，则θk维第k个对角线为1其他元素为0的toeplitz矩阵，g为固定系数矩阵，η为人工设定拟合误差，n为对球阵采集高阶声场进行球谐波分解的最高球谐波阶数。

通过cvx软件计算式(6)，获得ρ的对偶变量rv,l作为恢复后的球谐波域协方差矩阵。

本发明中步骤三具体步骤为：对第v个时间、第l个频率上的局部时频块，按照步骤二求解得到的球谐波域协方差矩阵rv,l作奇异值分解，得到奇异值若σ1/σ2大于给定阈值α，则认为该时频块通过检测。

本发明中步骤四具体步骤为：将通过步骤三检测的局部时频块的时频点{v,l}集合记为a。对于时频点{v,l}上的球谐波域协方差矩阵rv,l，将步骤三奇异值分解中得到的最大奇异值对应的奇异值向量记为us(v,l)，计算融合空间谱步骤为：

式(7)中其中θ和分别为仰角和方位角。搜索式(7)空间谱的谱峰，即得到该时间段内声源的位置。

本发明所公开的一种超分辨球阵声源定位方法，相对于已有的子空间类声源定位方法对噪声和混响具有更好鲁棒性，对邻近声源分辨力更好。具体为：

(1)本发明步骤二中公开一种新的球谐波域信号协方差矩阵恢复方法，与传统子空间类方法相比，对强相关、强噪声、短时信号等条件具有更好的鲁棒性。与球上超分辨方法(见参考文献[4]t.bendory,s.dekel,anda.feuer,“super-resolutiononthesphereusingconvexoptimization,”ieeetransactionsonsignalprocessing,,vol.63,no.9,pp.2253-2262,2015)相比具有更小的半正定约束维数，因而具有更小的计算量。

(2)基于本发明步骤二求解得到的球谐波域信号协方差矩阵，进行直接路径检测和空间谱计算，使声源定位结果对噪声和混响具有更好鲁棒性，对邻近声源分辨力更好。

附图说明

图1为本发明方法流程图。

图2为本发明方法房间混响环境下邻近声源空间谱对比图。

图3为本发明方法多固定声源定位结果对比图。

图4为本发明方法单运动声源定位结果对比图。

具体实施方式

下面结合附图对本发明作进一步说明。

一种适用于球麦克风阵列的超分辨声源定位方法，其特征在于，包括如下步骤：

步骤一、球阵采集高阶声场进行球谐波分解和短时傅里叶变换，计算局部时频块的自相关矩阵，并进行频率平滑；

步骤二、将步骤一获得的局部时频块的自相关矩阵代入半定规划问题求解，获得恢复后的球谐波域协方差矩阵；

步骤三、利用步骤二求解得到的球谐波域协方差矩阵作直接路径检测，即对求解得到的球谐波域协方差矩阵作奇异值分解，若最大奇异值与次大奇异值的比值大于设定的阈值，则认为该时频块通过检测；

步骤四、对通过检测的时频块计算空间谱，即利用步骤三中球谐波域协方差矩阵最大奇异值对应奇异值矢量作为信号子空间计算空间谱，并将各时频块空间谱叠加，搜索叠加后空间谱的谱峰，获得声源位置。

考虑一个半径为r由i个全向阵元组成的球阵，球阵的第i个阵元位于ri＝(r,φi)，其中其中θ和分别为仰角和方位角，接收信号表示为x(t)＝[x1(t),…,xi(t)]^t，t表示时间；

步骤一具体步骤为：对接收信号x(t)作短时傅里叶变换，得到接收信号的时频序列x(v,l)，其中v，l分别为时间和频率点索引，再对x(v,l)作球谐波变换：

式(8)中diag[·]表示对角矩阵，κl＝λ/2π，λl表示频点l对应的波长；

hn和jn分别是二类球汉克尔函数和一类球贝塞尔函数，h′n和j′n分别是hn和jn的导数；w为与球阵结构有关的固定加权系数矩阵，维矩阵，其第i列表示为:

式(3)中[·]^t表示矩阵转置，为球谐波函数；

为连带勒让德多项式，n为对球阵采集信号作球谐波分解的最高球谐波阶数；将时频序列p(v,l)划分为若干局部时频块；令第v个时间、第l个频率上的局部时频块的包含qv和ql个时间和频率点，计算该局部时频块上的信号自相关矩阵：

所述步骤二中采用如下半定规划问题：

其中tr(·)表示矩阵的迹，(·)^h表示矩阵的共轭转置，re(·)表示取实部，对于矩阵m，m≥0表示矩阵m半正定，||·||2表示欧几里得范数，i为单位阵，表示第v个时间、第l个频率上的局部时频块自相关矩阵，令为kronecker积，则θk维第k个对角线为1其他元素为0的toeplitz矩阵，g为固定系数矩阵，η为人工设定拟合误差，n为对球阵采集高阶声场进行球谐波分解的最高球谐波阶数；

通过计算(13)，获得ρ的对偶变量rv,l作为恢复后的球谐波域协方差矩阵。

所述步骤三中采用如下直接路径检测方法：

对给定时频块自相关矩阵即，对第v个时间、第l个频率上的局部时频块，按照步骤二求解得到的球谐波域协方差矩阵rv,l作奇异值分解，得到奇异值若σ1/σ2大于给定阈值α，则认为该时频块通过检测。

将通过步骤三检测的局部时频块的时频点{v,l}集合记为a；对于时频点{v,l}上的球谐波域协方差矩阵rv,l，将步骤三奇异值分解中得到的最大奇异值对应的奇异值向量记为us(v,l)，计算融合空间谱步骤为：

式(14)中其中θ和分别为仰角和方位角；搜索式(14)空间谱的谱峰，即得到该时间段内声源的位置。

本发明方法是基于球麦克风阵列，按照图1所述流程图实现声源定位。

步骤一中，频率平滑和自相关矩阵计算应选择合适频带，频率太低则受噪声影响大，精度不足；频率太高会导致球谐波变换出现空间混叠，一般应使n>κr，且i>(n+1)²。

步骤二中，固定系数矩阵g由如下方式计算：连带勒让德多项式可表示为

式(15)中{βn,k,l}为唯一确定系数。故

式(16)中令g＝[g00,g1(-1),g10,…,gnn]^t，利用式(16)，则gm,n由如下公式确定

式(17)中

步骤二中，人工设定拟合误差η一般设定为(n+1)²σ²，σ为rv,l最小奇异值开根号。实际中，可根据具体球阵结构、球谐波变换残余误差、麦克风幅相误差等情况适当调整。

步骤三中，给定阈值α一般选取2-6，根据语音信号、房间混响等条件适当调整。

步骤四中，一般假设通过直接路径检测的局部时频块中仅包含单一声源，即由rv,l奇异值分解计算信号子空间us(v,l)时，仅选取最大奇异值对应的奇异值矢量。

下面结合一些具体实施实例，对本发明作具体说明：

实施实例1

考虑一个矩形房间，房间尺寸为4×6×3m，麦克风阵列位于[x,y,z]＝[2.6,3.05,0.7]，两个声源位于[x,y,z]＝[1.69,4.05,1.57]和[x,y,z]＝[1.2,3.5,1.7]，相对俯仰方位角分别是和房间混响时间t60＝0.2秒，信噪比30db，对接收语音信号作短时傅里叶变换后选择单频点计算协方差矩阵，每帧快拍数ql＝5，分别计算dpd-music(见参考文献1)与本发明所公开方法的空间谱，如图2所示，其中图中□代表真实值：

可见，在上述房间混响环境下，dpd-music方法只能显示出一个谱峰，误差比较大，而本发明公开方法能够准确定位两个说话人位置，表现出更好的空间分辨率。

实施实例2

采用来自locata竞赛(见参考文献5h.w.c.evers,a.schmidt,h.mellmann,h.barfuss,p.a.naylor,andw.kellermann.thelocatachallengedatacorpusforacousticsourcelocalizationandtracking[c].ieee10thsensorarrayandmultichannelsignalprocessingworkshop(sam):410-414,2018.)固定多声源场景中的eigenmike(见参考文献6themhacoustics.em32eigenmikemicrophonearrayreleasenotes[k].25summitave,summit,usa,2013.)麦克风阵列数据进行实验。

房间尺寸为7.1×9.8×3m，其大概混响时间是t60＝0.55s。这里，将本文提出了球谐波原子范数方法与locata竞赛所提供的的基准方法music(见参考文献[7]christineevers,heinrichw.heinrichmellmann,alexanderschmidt,hendrikbarfuss,patricka.naylorandwalterkellermann.locatachallenge-evaluationtasksandmeasures[c].16thinternationalworkshoponacousticsignalenhancement(iwaenc),2018:565-569.)进行性能比较。本实例中短时傅里叶变换采用帧长为1440点的汉宁窗，采样频率为48khz，fft的点数为1024点，每帧帧长30ms，做频率平滑频点数15个，时间积累为15个快拍，主路径检测阈值α设置为4。图2和图3分别给出了多固定声源和单运动声源场景下的测试结果。

测试结果表明，本发明公开方法在多固定声源场景下优于参考文献[7]方法，在单运动声源场景下优于参考文献[1]方法。