CN105301563A - 一种基于一致聚焦变换最小二乘法的双声源定位方法 - Google Patents

Abstract

本发明公开了一种基于一致聚焦变换最小二乘法的双声源定位方法，在该方法中，用先设计的六元圆形麦克风阵列采集声源信号，并获得采集信号的协方差矩阵；利用频率范围的中心频率点定义聚焦变换矩阵，并由最小二乘方法求解聚焦变换矩阵；利用不同带宽内的中心频率点、一致聚焦矩阵及多信号分类方法，获得每个中心频率点所对应的信号空间谱，再利用频率点均值和时间快拍估计方法（MUSIC）求得信号空间谱平均估计值，从而获求声源方位角估计值。本发明方法声源定位估计准确性高，有效克服了方位模糊问题。

Description

一种基于一致聚焦变换最小二乘法的双声源定位方法

技术领域

本发明涉及声源定位领域，特别是一种基于一致聚焦变换最小二乘法的双声源定位方法。

背景技术

在阵列信号处理中，利用波达方向(DirectionofArrival，DOA)估计声源位置是一个新的研究方向，在声呐探测(见文献：王燕,邹男,梁国龙.强多途环境下水听器阵列位置近场有源校正方法[J].物理学报,2015,64(2):0243041-10)，语音识别与追踪(见文献：戚聿波,周士弘,张仁和,任云.一种基于β-warping变换算子的被动声源距离估计方法[J].物理学报,2015,64(7):0743011-6)，在未知环境中机器人的移动(见文献：居太亮.基于麦克风阵列声源定位算法研究[D].博士学位论文(成都：电子科技大学),2006)，以及下一代助听器等(见文献：苏林,马力,宋文化,郭圣明,鹿力成.声速剖面对不同深度声源定位的影响[J].物理学报,2015,64(2):0243021-8)领域中都成为研究的热点，精确的DOA估计在许多应用中都是其中的关键要素之一(见文献：时洁,杨德森,时胜国.基于矢量阵的运动声源柱面聚焦定位方法试验研究[J].物理学报,2012,61(12):1243021-15)。其中，最常见的是在电话会议中，通过知晓说话者扬声器的位置信息来引导相机，或者利用波束形成方法来增强捕获到的源信号，以避免翻领式麦克风带来的不便(见文献：DespoinaPavlidi,AnthonyGriffin,MatthieuPuigt,AthanasiosMouchtaris.Real-timemultiplesoundsourcelocalizationandcountingusingacircularmicrophonearray[J].IEEETransactionsonAudio,Speech,andLanguageProcessing,2013,21(10):2193-2206.)。

在声源定位的早期研究中，主要针对单个声源进行定位的(见文献：AliPourmohammad,SeyedMohammadAhadi.Realtimehighaccuracy3-DPHAT-basedsoundsourcelocalizationusingasimple4-microphonearrangement[J].IEEESystemsJournal,2012,6(3):455-468)。很多方法都是基于不同麦克风对间的到达时差(TimeDifferenceofArrival，TDOA)进行单声源定位的，其中广义互相关变换(GeneralizedCross-CorrelationPHAseTransform，GCC-PHAT)的方法是应用最为广泛的(见文献：QilianLiang,BaojuZhang,ChenglinZhao,YimingPi.TDOAforpassivelocalizationunderwaterversusterrestrialenvironment[J].IEEETransactionsParallelandDistributedSystems,2013,24(10):2100-2108)。文献(JingdongChen,JacobBenesty,YitengHuang.Timedelayestimationinroomacousticenvironments:anoverview[J].EURASIPJournalonAppliedSignalProcessing,2006,26503:1-19)给出了TDOA方位估计技术的综述；在单声源中，由于存在噪声、混响等的干扰，会影响麦克风对间采集的信号，使定位估计结果精度下降，文献(FrancescoNesta,MaurizioOmologo.GeneralizedstatecoherencetransformformultidimensionalTDOAestimationofmultiplesources[J].IEEETransactionsonAudio,Speech,andLanguageProcessing,2012,20(1):246-260.)给出了改进的TDOA估计方法。进一步研究发现，同时定位多个声源是一个更困难的问题。事实上，即使是一个短暂的声音都会使采集到的语音信号存在微小的重叠，可能会导致多个源信号的定位不精确(见文献：JacekP.Dmochowski,JacobBenesty,SofieneAffes.BroadbandMUSIC：Opportunitiesandchallengesformultiplesourcelocalization[C].IEEEWorkshoponApplicationsofSignalProcessingtoAudioandAcoustics,2007:18-21)。对于多声源定位问题，一种是传统的多信号分类(MultipleSignalClassification，MUSIC)方法，属于经典子空间方法，该方法针对窄带信号，根据采集到信号的子空间协方差矩阵，进行特征分解，进而估计源信号方位(见文献：DumiduS.Talagala,WenZhang.BroadbandDOAEstimationUsingSensorArraysonComplex-ShapedRigidBodies[J].IEEETransactionsonAudio,Speech,andLanguageProcessing,2013,21(8):1573-1585.)，由于语音信号的频率带宽一般为[300Hz，3000Hz]，经典子空间方法只是针对窄带信号的，会导致声源定位结果不准确；另一种是基于独立分量分析的源定位方法，该方法首先运用盲源分离(BlindSourceSeparation，BSS)方法获得分离后单个信号，然后再运用对单个源信号进行定位的方法进行定位估计(见文献：AnthonyLombard,YuanhangZheng,HerbertBuchner,WalterKellermann.TDOAEstimationforMultipleSoundSourcesinNoisyandReverberantEnvironmentsUsingBroadbandIndependentComponentAnalysis[J].IEEETransactionsonAudio,Speech,andLanguageProcessing,2011,19(6):1490-1503.)，但此类源定位方法，由于将语音信号先分离再定位，会破坏源信号的完整性，会使定位结果不稳定、定位精确度不高。另外，随着麦克风个数的增加，波达方向(DOA)估计的定位精度也会随之提高，但在实际情况下，通常使用线性麦克风阵列模型，会产生定位模糊问题(无法准确区分声源位于阵列的前方或者后方)。

如何解决现有技术的不足已成为声源定位领域亟待解决的一大难题。

发明内容

本发明所要解决的技术问题是为了克服现有技术的不足，而提供一种基于一致聚焦变换最小二乘法的双声源定位方法，本发明方法实时、有效地解决多声源定位问题，在一定程度上降低了噪声、混响等对语音信号的干扰，提高了定位的精度与稳定性。

本发明为解决上述技术问题采用以下技术方案：

根据本发明提出的一种基于一致聚焦变换最小二乘法的双声源定位方法，包括如下步骤：

步骤一：采用M个相同麦克风等间距的排列成半径为R的圆形麦克风阵列；

步骤二：根据麦克风阵列采集到的语音信号，得到信号的协方差矩阵，并在给定频率范围内定义中心频率点k_b；

步骤三：根据预设的一定测量范围内的任意角度，存在一个不随角度变化的一致聚焦变换，定义一致聚焦变换矩阵，并通过最小二乘方法求得聚焦变换矩阵；

步骤四：根据步骤二中定义的中心频率点，结合最小二乘法求得的一致聚焦变换矩阵，采用MUSIC方法求得每个中心频率点所对应的信号空间谱，进而求得信号空间谱的均值函数h_array(θ)；

步骤五：采用频率点均值和时间快拍估计的方法求得信号空间谱平均估计值，进而求得声源估计角度。

作为本发明所述的一种基于一致聚焦变换最小二乘法的双声源定位方法进一步优化方案，所述步骤二具体如下：

(201)、对第d个声源发出的语音信号s_d(t)进行采样，得到第m个麦克风采集到的离散信号x_dm(t)，其中，t为时间，d＝1,2，…,D，D表示声源的个数，m＝1,2，…,M；

(202)、根据(201)中采集到的离散信号x_dm(t)，经傅里叶变换后得到输入信号矢量X(k)为X(k)＝(X₁(k),…,X_M(k))^T，X_m(k)＝(X_1m(k),…,X_Dm(k))^T，信号的协方差矩阵为C_X＝E[X(k)X(k)^H]；

其中，k为波数且f为频率，c为声波传播的速度，上标T为转置运算符，X_dm(k)表示x_dm(t)的傅里叶变换，E[*]为期望的运算符，上标H为厄密共轭运算符；

(203)、在给定频率范围内划分频率间隔相等的B个子带，第b个子带的中心频率为k_b，其中，b＝1,…,B。

作为本发明所述的一种基于一致聚焦变换最小二乘法的双声源定位方法进一步优化方案，所述步骤三具体如下：

(301)根据带宽内定义的任意中心频率点k_b以及给定的聚焦频率点k₀，定义一致聚焦变换矩阵T(r,k_b)，得到

V(r,θ,k₀)＝T(r,k_b)V(r,θ,k_b)

其中，V(r,θ,k₀)表示声源距离阵列中心的距离为r、方位角度为θ、聚焦频率点为k₀的导向矢量矩阵，V(r,θ,k_b)表示声源距离阵列中心的距离为r、方位角度为θ、中心频率点为k_b的导向矢量矩阵；

(302)采用最小二乘方法，得到一致聚焦变换矩阵为

T(r,k_b)＝R(r,k_b)L(r,k_b)

其中，T(r,k_b)表示声源距离阵列中心的距离为r，中心频率点为k_b的一致聚焦变换矩阵，R(r,k_b)为矩阵V^H(r,θ,k₀)V(r,θ,k_b)的左奇异矢量矩阵，L(r,k_b)为矩阵V^H(r,θ,k₀)V(r,θ,k_b)的右奇异矢量矩阵。

作为本发明所述的一种基于一致聚焦变换最小二乘法的双声源定位方法进一步优化方案，所述(202)中的c＝343m.s^-1。

作为本发明所述的一种基于一致聚焦变换最小二乘法的双声源定位方法进一步优化方案，所述(203)中的B＝180。

作为本发明所述的一种基于一致聚焦变换最小二乘法的双声源定位方法进一步优化方案，所述步骤四中的h_array(θ)经下面公式得到，

$h_{array}\left(\mathrm{\θ}\right)=\frac{1}{B}\underset{b=1}{\overset{B}{\Σ}}\frac{1}{{\left(R\right(r,k_b\left)L\right(r,k_b\left)V\right(r,\mathrm{\θ},k_b\left)\right)}^H{U}_N\left(r\right)U_N^H\left(r\right)\left(R\right(r,k_b\left)L\right(r,k_b\left)V\right(r,\mathrm{\θ},k_b\left)\right)}$

其中，U_N(r)表示声源位置矢量的噪声子空间矩阵，下标N表示噪声，下标array表示阵列。

作为本发明所述的一种基于一致聚焦变换最小二乘法的双声源定位方法进一步优化方案，所述步骤五具体如下：

(501)、采用时间快拍估计方法得到协方差矩阵C_X的近似值

其中，t′＝T₀,2T₀…表示时间系数，T₀表示时间间隔，W表示窗口的长度，表示输入信号在时间间隔为T₀时、第l个傅里叶级数的近似值，为取整运算符；

(502)、根据求得信号空间谱平均估计值即

${\hat{h}}_{array}\left(\mathrm{\θ}\right)=\frac{1}{B}\underset{b=1}{\overset{B}{\Σ}}\frac{1}{{\left(R\right(r,k_b\left)L\right(r,k_b\left)V\right(r,\mathrm{\θ},k_b\left)\right)}^H{\hat{U}}_N\left(r\right){\hat{U}}_N^H\left(r\right)\left(R\right(r,k_b\left)L\right(r,k_b\left)V\right(r,\mathrm{\θ},k_b\left)\right)}$

其中，表示声源位置矢量的噪声子空间矩阵U_N(r)的近似值；

信号空间谱平均估计值的峰值所对应的角度即为声源估计角度

作为本发明所述的一种基于一致聚焦变换最小二乘法的双声源定位方法进一步优化方案，所述步骤三中一定测量范围内的任意角度为-90°～90°。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：

(1)本发明构建了六元圆形麦克风阵列，运用基于一致聚焦变换最小二乘法的宽带信号MUSIC方法，对室内近场双声源进行定位；为了确保声源定位的稳定性与精确性，求出麦克风对采集到的语音信号的协方差矩阵，并在频率范围内定义中心频率点，通过一致聚焦变换和最小二乘法，在不同的中心频率点上，运用多信号分类的方法求得每个中心频率点所对应的信号空间谱，再利用频率点均值和时间快拍估计的方法求得信号空间谱平均估计值，进而估计求得声源方位；

(2)本发明方法降低了噪声、混响等对语音信号的干扰，提高了定位的精度与稳定性；

(3)本发明在仿真与实测实验中都能准确地确定多声源方位，定位精度高与稳定性强，在语音信号处理领域，具有较强的实用性。

附图说明

图1是本发明的流程图。

图2是本发明设计的六元麦克风房间仿真模型。

图3是本发明在相同混响时间(T₆₀)，不同信噪比(SNR)条件下声源定位结果：(a)为T₆₀＝250ms，SNR＝20dB条件下声源定位结果；(b)为T₆₀＝250ms，SNR＝15dB条件下声源定位结果；(c)为T₆₀＝250ms，SNR＝10dB条件下声源定位结果；(d)为T₆₀＝250ms，SNR＝5dB条件下声源定位结果；(e)为T₆₀＝250ms，SNR＝0dB条件下声源定位结果。

图4是本发明在不同混响时间、相同信噪比条件下声源定位结果。

图5是本发明与传统的MUSIC、BSS-TDOA方法的声源定位结果的比较。

图6是本发明在全消声实验室声源定位结果：(a)为三维图；(b)为侧视图。

图7是本发明在非消声实验室声源定位结果：(a)为三维图；(b)为侧视图。

图8是本发明仿真实验与实测实验进行声源定位结果的对比图。

具体实施方式

下面结合附图对本发明的技术方案做进一步的详细说明：

本发明方法的流程参见图1，本发明是一种基于一致聚焦变换最小二乘法的麦克风阵列双声源定位方法，利用六元麦克风阵列，结合语音信号特性进行声源定位，其具体实施步骤如下：

步骤一：建立圆形麦克风阵列模型；

建立圆形麦克风阵列模型，如图2所示，由M个相同的麦克风等间距的排列组成，阵列的半径为R，M个阵元围绕旋转成一个圆阵，这里约定：当且仅当i≡j(modM)时(mod表示数学中的求余数运算)，第i个麦克风与第j个麦克风为同一个麦克风。声源S(t)位于近场条件下，满足

$r\≤\frac{2L^2}{\mathrm{\λ}}---\left(1\right)$

式中，r为声源距离阵列中心的距离，L为阵列的最大尺寸(这里L＝2R)，λ为语音信号的波长。

步骤二：求麦克风阵列采集到的语音信号的协方差矩阵，并在频率范围内定义中心频率点；

(201)在室内环境中有D个指向性声源，同时也存在着无指向性的环境噪声，由M＞D个全指向性麦克风采集声场中的语音信号。第d(d＝1，2，…,D)个声源的位置矢量在极坐系中表示为r_d＝(r_d,θ_d)^T，r_d表示第d个声源距离阵列中心的距离，θ_d表示第d个声源的方位角度，上标T表示转置运算符。设声波传播的速度c＝343m.s^-1。

(202)第d个声源的语音信号为s_d(t)，则其频域值为

S_d(k)＝∫s_d(t)e^-2jπftdt＝∫s_d(t)e^-jkctdt(2)

式中，j为虚数单位，e表示自然指数，f表示频率，表示波数，则源信号矢量为S(k)＝(S₁(k),…,S_D(k))^T，S_d(k)表示第d个声源的频域信号。

圆形麦克风阵列中，第m个麦克风采集到的第d个声源的语音信号为x_dm(t)(d＝1，2，…,D；m＝1,…M)，则输入信号矢量为X(k)＝(X₁(k),…,X_M(k))^T，X_m(k)＝(X_1m(k),…,X_Dm(k))^T，X_dm(k)表示x_dm(t)的傅里叶变换，且

X(k)＝V(r₁,…,r_D,k)S(k)+B(k)(3)

式中

V(r₁,…,r_D,k)＝(V(r₁,k),…,V(r_D,k))(4)

是由与第d个声源相关的导向矢量矩阵V(r_d,k)构成的M×D阶矩阵，r_d表示第d个声源的位置矢量。B(k)＝(B₁(k),…,B_M(k))^T表示每个麦克风上的加性噪声，B_m(k)第m个麦克风上的加性噪声，假设噪声为零均值、稳定的白噪声，每个麦克风上的噪声能量是相等的，并且噪声信号与声源信号间是相互独立的，则

$E\[B\left(k\right)B{\left(k\right)}^H\]={\mathrm{\σ}}_N^2{I}_M---\left(5\right)$

E[B(k)(V(r_d,k)S(k))^H]＝0(6)

式中，E[*]表示期望的运算符，Ι_M表示M×M的单位矩阵，H表示厄密共轭运算符，V(r_d,k)表示导向矢量矩阵，表示噪声方差。

在近场条件下，需要考虑每个声源与麦克风阵列间的距离，则导向矢量矩阵V(r_d,k)＝V(r_d,θ_d,k)，其中，第d个分量为

$V_d(r_d,{\mathrm{\θ}}_d,k)=r_d{e}^{{\mathrm{jkr}}_d}\frac{e^{-jk\sqrt{r_d^2-2r_d{\mathrm{cos\θ}}_d}}}{\sqrt{r_d^2-2r_d{\mathrm{cos\θ}}_d}}---\left(7\right)$

根据输入信号X(k)，可求得信号的协方差矩阵C_X，为

C_X＝E[X(k)X(k)^H](8)

根据式(3)、(5)、(6)，C_X可进一步表示为

$C_X=V(r_1,...,r_D,k)C_S{V}^H(r_1,...,r_D,k)+{\mathrm{\σ}}_N^2{I}_M---\left(9\right)$

式中，C_S为源信号D×D阶协方差矩阵，C_B为噪声信号的协方差矩阵，

C_S＝E[S(k)S(k)^H](10)

$C_B=E\[B\left(k\right)B{\left(k\right)}^H\]={\mathrm{\σ}}_N^2{I}_M---\left(11\right)$

M×M阶矩阵C_Y＝V(r₁,…,r_D,k)C_SV^H(r₁,…,r_D,k)，矩阵C_Y满足埃尔米特对称、半正定，因此可得M个实的、非负的第m个特征值λ_m，以及相关的第m个正交特征向量U_m(m＝1,…M)。这里假定声源信号间是相互独立的，V(r₁,…,r_D,k)为满秩矩阵，C_Y的秩为D那么其特征值满足λ₁≥λ₂≥…≥λ_D＞λ_D+1＝…＝λ_M＝0。由上述推导可以注意到向量U₁,…,U_D与V(r₁,…,r_D,k)生成空间的范围是一致的，因此可根据导向矢量所形成的D维子空间S来估计声源位置，且D维子空间S被称作信号子空间。

根据

$C_X=C_Y+{\mathrm{\σ}}_N^2{I}_M---\left(12\right)$

可得

式中，U_S＝(U₁…U_D)∈R^M×D为信号子空间S的矩阵，是由上述的特征向量所构成的D阶矩阵，且与特征值相关；U_N＝(U_D+1…U_M)∈R^M×(M-D)为噪声子空间N的矩阵，是由余下的特征向量构成的M-D阶矩阵，且其特征值为

信号子空间矩阵U_S与噪声子空间矩阵U_N有如下关系，

(U_S|U_N)^H(U_S|U_N)＝I_M(14)

(203)由于语音信号的频率带宽一般为[300Hz，3000Hz]，为了在处理时确保语音信号的完整性，在给定频率范围内划分频率间隔相等的B个子带，第b个子带的中心频率为k_b，其中，b＝1,…,B，这里取B＝180。

步骤三：根据一定范围内的任意角度，存在一个不随角度变化的一致聚焦变换，定义聚焦变换矩阵，并通过最小二乘方法求解；

(301)对于一定测量范围内的任意角度(这里设定测量的范围为-90°～90°)，存在一个不随角度变化的一致聚焦变换。根据带宽内定义的任意中心频率点k_b以及给定的聚焦频率点k₀，定义一致聚焦变换矩阵T(r,k_b),b＝1,…,B，任取(r,θ)，有

V(r,θ,k₀)＝T(r,k_b)V(r,θ,k_b)(15)

则变换T(r,k_b)称为一致聚焦变换。

(302)利用最小二乘方法，对式(15)求解

$\left\{\begin{array}{c}min\left|\right|V(r,\mathrm{\θ},k_0)-T(r,k_b)V(r,\mathrm{\θ},k_b)|{|}_F^2\\ T^H(r,k_b)T(r,k_b)=I\end{array}\right.---\left(16\right)$

可得，一致聚焦变换矩阵

T(r,k_b)＝R(r,k_b)L(r,k_b)(17)

式中，R(r,k_b)为矩阵V^H(r,θ,k₀)V(r,θ,k_b)的左奇异矢量矩阵，L(r,k_b)为矩阵V^H(r,θ,k₀)V(r,θ,k_b)的右奇异矢量矩阵；。

步骤四：根据步骤二中定义的中心频率点，结合最小二乘法求得的一致聚焦变换矩阵，利用MUSIC方法求得每个中心频率点所对应的信号空间谱，进而求信号空间谱的均值函数；

在每个中心频率点k_b，计算T(r,k_b)X(k_b)的二阶统计量，再求和，结合式(12)，得输入信号的聚焦协方差矩阵，

${\mathrm{\Γ}}_X\left(r\right)=\underset{b=1}{\overset{B}{\Σ}}T(r,k_b)C_X\left(k_b\right)T^H(r,k_b)---\left(18\right)$

式中，C_X(k_b)表示在中心频率点k_b时，输入信号的协方差矩阵。

根据式(15)、式(17)，式(18)得

$\begin{array}{c}{\mathrm{\Γ}}_X\left(r\right)=V(r,\mathrm{\θ},k_0){\mathrm{\Γ}}_S{V}^H(r,\mathrm{\θ},k_0)+{\mathrm{\σ}}_N^2{\mathrm{\Γ}}_N\left(r\right)\\ =\underset{b=1}{\overset{B}{\Σ}}\left(T\right(r,k_b\left)V\right(r,\mathrm{\θ},k_b\left)\right){\mathrm{\Γ}}_S{\left(T\right(r,k_b\left)V\right(r,\mathrm{\θ},k_b\left)\right)}^H+{\mathrm{\σ}}_N^2{\mathrm{\Γ}}_N\left(r\right)\\ =\underset{b=1}{\overset{B}{\Σ}}\left(R\right(r,k_b\left)L\right(r,k_b\left)V\right(r,\mathrm{\θ},k_b\left)\right){\mathrm{\Γ}}_S{\left(R\right(r,k_b\left)L\right(r,k_b\left)V\right(r,\mathrm{\θ},k_b\left)\right)}^H+{\mathrm{\σ}}_N^2{\mathrm{\Γ}}_N\left(r\right)\end{array}---\left(19\right)$

式中，

${\mathrm{\Γ}}_S=\underset{b=1}{\overset{B}{\Σ}}{C}_S\left(k_b\right)---\left(20\right)$

式中，C_S(k_b)表示在中心频率点k_b时，源信号的协方差矩阵。

根据式(19)，可得噪声信号聚焦协方差矩阵为

$\begin{array}{c}{\mathrm{\Γ}}_N\left(r\right)=\underset{b=1}{\overset{B}{\Σ}}T(r,k_b)C\left(k_b\right)T^H(r,k_b)\\ =\underset{b=1}{\overset{B}{\Σ}}\left(R\right(r,k_b\left)L\right(r,k_b\left)\right)C\left(k_b\right){\left(R\right(r,k_b\left)L\right(r,k_b\left)\right)}^H\end{array}---\left(21\right)$

式中

$E\[B\left(k_b\right)B^H\left(k_b\right)\]={\mathrm{\σ}}_N^{2}C\left(k_b\right)---\left(22\right)$

由输入信号的聚焦协方差矩阵Γ_X(r)与噪声信号的聚焦协方差矩阵Γ_N(r)，可构成矩阵对(Γ_X(r),Γ_N(r))，其第m个特征值为μ_m，且μ_m＞μ_m+1，第m个特征向量为U_m,m＝1,…,M。那么U_S(r)＝(U₁,…,U_D)，U_N(r)＝(U_D+1,…,U_M)，且有

${\mathrm{\μ}}_{D+1}=...={\mathrm{\μ}}_M={\mathrm{\σ}}_N^2---\left(23\right)$

V^H(r,θ,k₀)U_N(r)＝0(24)

基于上述分析，可得信号的空间谱函数为

$\begin{array}{c}{h}_b(\mathrm{\θ},k_b)=\frac{1}{V^H(r,\mathrm{\θ},k_0)U_N\left(r\right)U_N^H\left(r\right)V(r,\mathrm{\θ},k_0)}\\ =\frac{1}{{\left(T\right(r,k_b\left)V\right(r,\mathrm{\θ},k_b\left)\right)}^H{U}_N\left(r\right)U_N^H\left(r\right)\left(T\right(r,k_b\left)V\right(r,\mathrm{\θ},k_b\left)\right)}\\ =\frac{1}{{\left(R\right(r,k_b\left)L\right(r,k_b\left)V\right(r,\mathrm{\θ},k_b\left)\right)}^H{U}_N\left(r\right)U_N^H\left(r\right)\left(R\right(r,k_b\left)L\right(r,k_b\left)V\right(r,\mathrm{\θ},k_b\left)\right)}\end{array}---\left(25\right)$

式中，U_N(r)表示声源位置矢量的噪声子空间矩阵，下标N表示Noise(噪声)。

然后，根据式(25)，可得信号空间谱的均值函数，

$\begin{array}{c}{h}_{array}\left(\mathrm{\θ}\right)=\frac{1}{B}\underset{b=1}{\overset{B}{\Σ}}{h}_b(\mathrm{\θ},k_b)\\ =\frac{1}{B}\underset{b=1}{\overset{B}{\Σ}}\frac{1}{{\left(R\right(r,k_b\left)L\right(r,k_b\left)V\right(r,\mathrm{\θ},k_b\left)\right)}^H{U}_N\left(r\right)U_N^H\left(r\right)\left(R\right(r,k_b\left)L\right(r,k_b\left)V\right(r,\mathrm{\θ},k_b\left)\right)}\end{array}---\left(26\right)$

式中，下标array表示阵列，

因为式(26)正交特性趋向于0，则平均空间谱函数的峰值所对应的角度θ，即为声源位置角度估计值。

步骤五：结合实际情况：仅有麦克风采集到的语音信号可用，运用频率点均值和时间快拍估计的方法求得信号空间谱平均估计值，进而求得声源估计角度；

在实际中，由于C_X是未知的，仅可以利用麦克风采集到的语音信号x(t)，而且矢量信号X(k)的复包络值也不能准确的确定。因此需要计算协方差矩阵C_X、矢量信号X(k)的近似值，这里采用时间快拍估计方法来求近似值。设时间系数t′＝T₀,2T₀…，T₀表示时间间隔，一方面x(t)以(l为整数)速率进行采样，因此在每个t′上，在快拍为时，通过傅里叶级数(FS)估计X(k)的近似值为另一方面，在t′上估计C_X，运用W长度的滑动窗口在T₀空间进行采样再由加权求和方法求得的估计值替代定义的期望值。在W长度的窗口中，期望C_X的近似值是完全基于得到的，且二者是相互独立的，这排除了快拍使用重叠的可能性，即

式中，表示输入信号在时间间隔为T₀时，第l个傅里叶级数的近似值，为取整运算符；

根据上述分析，式(26)可进一步表示为

${\hat{h}}_{array}\left(\mathrm{\θ}\right)=\frac{1}{B}\underset{b=1}{\overset{B}{\Σ}}\frac{1}{{\left(R\right(r,k_b\left)L\right(r,k_b\left)V\right(r,\mathrm{\θ},k_b\left)\right)}^H{\hat{U}}_N\left(r\right){\hat{U}}_N^H\left(r\right)\left(R\right(r,k_b\left)L\right(r,k_b\left)V\right(r,\mathrm{\θ},k_b\left)\right)}---\left(28\right)$

式中，表示声源位置矢量的噪声子空间矩阵U_N(r)的近似值，下标N表示Noise(噪声)。

仿真环境为5.5m×3.3m×2.3m的房间冲激响应模型，运用含有6个麦克风的圆形阵列，相邻麦克风间的夹角为60°，阵列的直径为40cm，声速c＝343m/s，混响时间T₆₀＝250ms。根据近场条件，声源响应在r＝0.6m～1.6m范围内，据此设定声源S₁的角度为θ₁＝60°，距离阵列中心距离为0.7m，声源S₂的角度为θ₂＝-20°，距离阵列中心距离为1.2m；声源与阵列在同一个平面上，且二个声源信号相互独立、能量相等。环境噪声SNR分别取0dB、5dB、10dB、15dB、20dB。对于采集的语音信号，设定帧长为512点，帧移为160点，FFT的长度为1024点，采样率为16000Hz，窗函数选择汉明窗，窗长取150点。

实测环境为全消声实验室、非消声实验室，房间尺寸(5.5m×3.3m×2.3m)、阵列摆放位置与仿真环境相同，声源高度、阵列高度都为1.2m，阵列为6个麦克风的圆形阵列。实验器材：数据采集设备为16通道的PXIE-4496数据采集卡、配套PC机(Intel2GHzCorei7CPU，2GBRAM)；声源为AM012人工嘴、便携式音箱；麦克风为的简易声音传感器模块(全向性，工作电压5V)。由于人工嘴在通电工作时产生的嘶嘶声、房间换气扇转动时产生的呼呼声，实测环境下的信噪比平均为20dB。

图3是本发明方法在相同混响时间(T₆₀)，不同信噪比(SNR)条件下声源定位结果。图4是本发明方法在不同混响时间、相同信噪比条件下声源定位结果。图5是本发明方法与传统的MUSIC、BSS-TDOA方法的声源定位结果的比较。图6是本发明方法在全消声实验室声源定位结果。图7是本发明方法在非消声实验室声源定位结果。图8是本发明方法仿真实验与实测实验进行声源定位结果的比较。图3、图4、图5、图6、图7都是用来说明本发明定位效果。

相同混响时间(T₆₀)、不同信噪比(SNR)条件下声源定位结果：

混响时间为T₆₀＝250ms，信噪比(SNR)不同分别为20dB、15dB、10dB、5dB、0dB条件下，六元麦克风阵列的声源定位结果。

图3表明，在声源与阵列间距离的增加情况下，随着信噪比的降低，声源定位结果的精确度降低，在近场条件下，图3中的(a)、图3中的(b)能准确、稳定地反应声源定位的结果，图3中的(c)、图3中的(d)、图3中的(e)能反应出声源定位的结果，但是随着信噪比降低会出现伪峰、出现局部衰减影响声源定位。图3中的(a)在20dB时，曲线很平滑；图3中的(b)在15dB时，在-40°附近产生幅值较小的伪峰；图3中的(c)在10dB时，在-40°附近产生幅值较小的伪峰，在1.0m附近产生局部衰减；图3中的(d)在5dB时，在40°、-30°、-70°附近产生伪峰，在0.8m、1.0m附近产生局部衰减；图3中的(e)在0dB时，在80°、30°、0°、-60°附近产生伪峰，在0.7m～1.0m附近产生局部衰减；但由图3中的(d)、图3中的(e)可看出，虽然有局部衰减与伪峰的影响，但仍能较为准确的得到声源位置。因此从总体上看，本发明方法能准确、稳定的确定声源位置结果。

不同混响时间、相同信噪比条件下声源定位结果：

图4表明，在相同信噪比下，混响时间对声源定位结果的影响较小。

本发明方法与传统的MUSIC、BSS-TDOA方法的声源定位结果的对比：

声源S₁、S₂分别位于{[10°、-10°]，[20°、-20°]，[30°、-30°]，[40°、-40°]，[50°、-50°]，[60°、-60°]，[70°、-70°]，[80°、-80°]，[90°、-90°]}位置。

图5表明，在信噪比为20dB、混响时间为250ms的条件下，本发明方法能较为准确的确定声源S₁、S₂的位置，传统的MUSIC方法的估计误差基本在8°左右；而BSS-TDOA的方法，由于存在了盲源分离与声源定位二个步骤，对声源定位的精确度造成影响，估计误差基本在10°左右。

在全消声实验室声源定位结果：

图6是本发明方法在全消声实验室声源定位结果：图6中的(a)为三维图；图6中的(b)为侧视图。图6表明，在全消声实验室中，由于仅存在噪声因素的影响，因此本发明方法能准确的测得声源位置，声源S₁为-21.4°、S₂为61.5°。

在非全消声实验室声源定位结果：

图7是本发明方法在非消声实验室声源定位结果：图7中的(a)为三维图；图7中的(b)为侧视图。图7表明，在非消声实验室中，由于存在混响、噪声等因素的影响，本发明方法仍能较为准确的测得声源位置声源S₁为-19.4°、S₂为58.7°。

图8是本发明方法仿真实验与实测实验进行声源定位结果的对比图，仿真实验与实测实验进行声源定位结果的对比：

声源S₁、S₂分别位于{[10°、-10°]，[20°、-20°]，[30°、-30°]，[40°、-40°]，[50°、-50°]，[60°、-60°]}位置。

由于存在实际环境噪声、混响以及采集设备A/D转换等影响因素，实测结果与仿真结果存在一定偏差，在仿真实验中声源定位平均绝对估计误差S₁为0.7°、S₂为1.1°，在实测实验中，全消声实验室情况下平均绝对估计误差S₁为1.3°、S₂为1.5°，非消声实验室情况下平均绝对估计误差S₁为1.9°、S₂为2.3°。

本发明方案所公开的技术手段不仅限于上述实施方式所公开的技术手段，还包括由以上技术特征任意组合所组成的技术方案。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也视为本发明的保护范围。

Claims (8)

1.一种基于一致聚焦变换最小二乘法的双声源定位方法，其特征在于，包括如下步骤：

步骤一：采用M个相同麦克风等间距的排列成半径为R的圆形麦克风阵列；

步骤二：根据麦克风阵列采集到的语音信号，得到信号的协方差矩阵，并在给定频率范围内定义中心频率点k_b；

步骤三：根据预设的一定测量范围内的任意角度，存在一个不随角度变化的一致聚焦变换，定义一致聚焦变换矩阵，并通过最小二乘方法求得聚焦变换矩阵；

步骤四：根据步骤二中定义的中心频率点，结合最小二乘法求得的一致聚焦变换矩阵，采用MUSIC方法求得每个中心频率点所对应的信号空间谱，进而求得信号空间谱的均值函数h_array(θ)；

步骤五：采用频率点均值和时间快拍估计的方法求得信号空间谱平均估计值，进而求得声源估计角度。

2.根据权利要求1所述的一种基于一致聚焦变换最小二乘法的双声源定位方法，其特征在于，所述步骤二具体如下：

(201)、对第d个声源发出的语音信号s_d(t)进行采样，得到第m个麦克风采集到的离散信号x_dm(t)，其中，t为时间，d＝1,2，…,D，D表示声源的个数，m＝1,2，…,M；

(202)、根据(201)中采集到的离散信号x_dm(t)，经傅里叶变换后得到输入信号矢量X(k)为X(k)＝(X₁(k),…,X_M(k))^T，X_m(k)＝(X_1m(k),…,X_Dm(k))^T，信号的协方差矩阵为C_X＝E[X(k)X(k)^H]；

其中，k为波数且f为频率，c为声波传播的速度，上标T为转置运算符，X_dm(k)表示x_dm(t)的傅里叶变换，E[*]为期望的运算符，上标H为厄密共轭运算符；

(203)、在给定频率范围内划分频率间隔相等的B个子带，第b个子带的中心频率为k_b，其中，b＝1,…,B。

3.根据权利要求2所述的一种基于一致聚焦变换最小二乘法的双声源定位方法，其特征在于，所述步骤三具体如下：

(301)根据带宽内定义的任意中心频率点k_b以及给定的聚焦频率点k₀，定义一致聚焦变换矩阵T(r,k_b)，得到

V(r,θ,k₀)＝T(r,k_b)V(r,θ,k_b)

其中，V(r,θ,k₀)表示声源距离阵列中心的距离为r、方位角度为θ、聚焦频率点为k₀的导向矢量矩阵，V(r,θ,k_b)表示声源距离阵列中心的距离为r、方位角度为θ、中心频率点为k_b的导向矢量矩阵；

(302)采用最小二乘方法，得到一致聚焦变换矩阵为

T(r,k_b)＝R(r,k_b)L(r,k_b)

其中，T(r,k_b)表示声源距离阵列中心的距离为r，中心频率点为k_b的一致聚焦变换矩阵，R(r,k_b)为矩阵V^H(r,θ,k₀)V(r,θ,k_b)的左奇异矢量矩阵，L(r,k_b)为矩阵V^H(r,θ,k₀)V(r,θ,k_b)的右奇异矢量矩阵。

4.根据权利要求2所述的一种基于一致聚焦变换最小二乘法的双声源定位方法，其特征在于，所述(202)中的c＝343m.s^-1。

5.根据权利要求2所述的一种基于一致聚焦变换最小二乘法的双声源定位方法，其特征在于，所述(203)中的B＝180。

6.根据权利要求3所述的一种基于一致聚焦变换最小二乘法的双声源定位方法，其特征在于，所述步骤四中的h_array(θ)经下面公式得到，

$h_{array}\left(\mathrm{\θ}\right)=\frac{1}{B}\underset{b=1}{\overset{B}{\Σ}}\frac{1}{{\left(R\right(r,k_b\left)L\right(r,k_b\left)V\right(r,\mathrm{\θ},k_b\left)\right)}^H{U}_N\left(r\right)U_N^H\left(r\right)\left(R\right(r,k_b\left)L\right(r,k_b\left)V\right(r,\mathrm{\θ},k_b\left)\right)}$

其中，U_N(r)表示声源位置矢量的噪声子空间矩阵，下标N表示噪声，下标array表示阵列。

7.根据权利要求6所述的一种基于一致聚焦变换最小二乘法的双声源定位方法，其特征在于，所述步骤五具体如下：

(501)、采用时间快拍估计方法得到协方差矩阵C_X的近似值

其中，t′＝T₀,2T₀…表示时间系数，T₀表示时间间隔，W表示窗口的长度，表示输入信号在时间间隔为T₀时、第l个傅里叶级数的近似值，为取整运算符；

(502)、根据求得信号空间谱平均估计值即

${\hat{h}}_{array}\left(\mathrm{\θ}\right)=\frac{1}{B}\underset{b=1}{\overset{B}{\Σ}}\frac{1}{{\left(R\right(r,k_b\left)L\right(r,k_b\left)V\right(r,\mathrm{\θ},k_b\left)\right)}^H{\hat{U}}_N\left(r\right){\hat{U}}_N^H\left(r\right)\left(R\right(r,k_b\left)L\right(r,k_b\left)V\right(r,\mathrm{\θ},k_b\left)\right)}$

其中，表示声源位置矢量的噪声子空间矩阵U_N(r)的近似值；

信号空间谱平均估计值的峰值所对应的角度即为声源估计角度

8.根据权利要求1所述的一种基于一致聚焦变换最小二乘法的双声源定位方法，其特征在于，所述步骤三中一定测量范围内的任意角度为-90°～90°。