CN105068048A - 基于空间稀疏性的分布式麦克风阵列声源定位方法 - Google Patents

Abstract

本发明提供一种基于空间稀疏性的分布式麦克风阵列声源定位方法，通过两步离散余弦变换(Discrete？Cosine？Transform,DCT)方式进行信号特征提取，然后用该特征构建稀疏定位模型，以便综合利用语音信号的短时和长时特性，同时能降低模型维数；并利用字典学习方式对模型进行修正；最后采用近似l₀范数重构算法对稀疏信号进行重构，得出声源位置。本发明对阵列结构无特殊要求、计算量少，有效地提高了抗噪声和抗混响能力以及定位精度，更适合用于室内声源定位，可广泛应用于视频会议系统、语音识别系统以及智能机器人等各个领域。

Description

基于空间稀疏性的分布式麦克风阵列声源定位方法

技术领域

本发明属于声源定位技术领域，具体涉及一种基于空间稀疏性的分布式麦克风阵列声源定位方法。

背景技术

基于麦克风阵列的声源定位作为定位领域的一个重要的研究方向，具有广泛的应用前景。在民用领域中它被广泛用于电话会议系统、视频会议系统和可视电话系统中，通过对发声者的定位，实现会议现场声源及外接摄像机的实时跟踪。在助听器中的麦克风阵列被用作语音增强的预处理系统，通过检测声源的位置信息，将阵列指向感兴趣的声源方向，从而抑制掉其它方向的干扰和噪声。在军事领域，可利用麦克风阵列声源定位技术实现狙击手定位，地面或低空声目标的定位与追踪等。分布式麦克风阵列由于比传统的单阵列具有更大的空间孔径，可以获得更好的声源定位性能，因此基于分布式麦克风阵列的声源定位方法成为当前麦克风阵列定位领域研究的热点之一。

与声纳和雷达系统等传统的阵列定位技术所处理的对象不同，语音信号为宽带信号，具有短时平稳特性，且所处环境一般还具有高混响，噪声大等特点。目前，常用的声源定位算法一般可划分为三类：一是基于波束成型的方法；二是基于高分辨率谱估计的方法；三是基于时延估计的方法。其中基于波束成形方法通过对麦克风阵列接收信号进行滤波、加权求和，直接控制麦克风阵列指向使波束具有最大输出功率的方向，可实现多声源定位。但该类方法需要知道声源和噪声的先验知识，并存在对初值敏感，计算量大，不利于实时处理等问题。基于高分辨率谱估计的方法在理论上可以对声源的方向进行有效估计，但该类算法是针对窄带信号的，并且该类算法无法处理高度相关的信号，因此混响会给算法的定位精度带来较大影响。基于时延估计的方法是利用广义互相关等时延估计算法求出信号到不同麦克风的相对时延，并利用时延信息与麦克风阵列的空间位置关系估计声源位置。由于这一类方法计算量较小，容易实现实时定位，因此研究得较多，但其同样易受混响因素的影响，特别是用于多声源定位时，性能严重下降。

近年来，压缩感知(CompressiveSensing,CS)理论的成熟和发展给传统的定位方法带来了一场革命性的变革，在声源定位的研究上也揭开了崭新一页。在麦克风阵列声源定位系统中，声源一般可认为是一个点源，同时声源的数目有限，因此声源定位问题本质上是稀疏的。根据定位问题的天然稀疏性，Cevher等人根据声音信号的传播模型将声源定位问题建模为稀疏逼近问题(CevherV,BaraniukR.Compressivesensingforsensorcalibration.5thIEEESensorArrayandMultichannelSignalProcessingWorkshop.2008.175-178.)，展示了较好的定位性能。文献(JiangH,MathewsB,WilfordP.Soundlocalizationusingcompressivesensing.Proceedingsofthe1stInternationalConferenceonSensorNetworks.2012.159-166.)则利用压缩感知方法实现了声源信号达到时间差(timedifferenceofarrival,TDOA)的精确估计，可以获得比广义互相关方式更高的估计精度。Simard等人利用基于格林函数的声源稀疏定位模型进行目标识别和定位，获得了比传统波束成型方法更好的性能(SimardP,AntoniJ.Acousticsourceidentification:Experimentingthel₁minimizationapproach.AppliedAcoustics,2013,74(7):974-986.)。这些工作都显示了压缩感知理论与技术强大的生命力及其在声源定位研究中广阔的应用前景。

然而，与传统麦克风阵列定位方法一样，基于CS理论的声源定位方法依然面临着复杂声学环境的影响，导致其定位精度难有实质性提高。一是噪声因素的影响，由于目前稀疏重构普遍采用l₁范数代替l₀范数进行稀疏约束，而在l₁范数约束模型中，目标函数会对大系数施加更多的约束以保证整个代价函数的收敛性。在定位应用中，大系数往往对应于声源目标，而小系数可能对应于噪声，这样常用的贪婪算法和凸优化算法等可能在重构过程中削弱大系数的贡献，而对小系数没有施加过多约束，从而导致重构准确性下降，在低信噪比时甚至会出现将噪声误判为定位目标。另一方面室内环境中普遍存在混响，而且墙壁、天花板等对声音的反射、散射作用具有时变特性和不可预知性，使得实际信号往往与稀疏定位模型之间存在较大误差，最终影响到定位性能。

发明内容

本发明的目的是提供一种基于空间稀疏性的分布式麦克风阵列声源定位方法，该方法能够有效提高定位精度，抗混响能力强，自适应能力强。

本发明的目的是针对现有技术中存在的不足，从构造语音信号特征入手，提出一种鲁棒的基于CS理论的声源定位方法。该方法首先通过两步离散余弦变换(DiscreteCosineTransform,DCT)方式进行信号特征提取，然后用该特征构建稀疏定位模型，以便综合利用语音信号的短时和长时特性，并能降低模型维数。为了克服l₁范数代替l₀范数进行稀疏约束时可能造成的约束错误问题，本发明提出一种近似l₀范数稀疏重构算法，通过构造合理的l₀范数近似函数，直接在l₀范数意义下进行稀疏重构，可以有效提高噪声条件下的重构精度。同时利用字典学习技术，根据接收信号不断调整冗余字典，使之能够克服实际信号与稀疏模型之间的失配问题，使得稀疏重构结果能够精确反映出位置信息。

本发明提供了如下的技术方案：

一种基于空间稀疏性的分布式麦克风阵列声源定位方法，包括如下步骤：

S1、定位系统建立：分布式麦克风定位系统由M个已知自身位置的麦克风和K个待定位声源组成，采用分步格点划分方法将整个定位区域均匀划分为若干个格点；每个麦克风分别接收声源发出的信号，并传送给定位中心；

S2、语音信号特征提取：定位中心将任一麦克风接收到的语音信号首先经过加窗处理，把原来长度为P×1的语音信号向量r分解为J个长度为Q×1的短帧，即：

$r={\left[\begin{array}{c}{r}_1\\ r_2\\ .\\ .\\ .\\ r_P\end{array}\right]}_{P\×1}----->Z={\left[\begin{array}{c}{z}_1\\ z_2\\ .\\ .\\ .\\ z_J\end{array}\right]}_{J\×Q}---\left(1\right)$

其中，r_i(i＝1,…,P)表示输入信号向量r的每一个分量，Z是一个维数为J×Q的矩阵，每一行代表经过加窗处理后的1帧数据；

接下来对这J帧信号分别进行一次DCT变换：

式中D(·)表示DCT变换，表示z_i(i＝1,…,J)经过DCT变换后的结果；对变化后的每一帧数据通过除以该帧的最大值实现归一化处理，然后对每一帧数据进行求平均计算：

接下来考虑连续多帧信号的长时特性，为此对求平均后的向量s再进行一次DCT变换，得到

向量的长度仍为J×1，也即基于两步DCT变换的特征提取方法将计算复杂度从样本长度级降低到帧长数量级；

S3、建立稀疏定位模型：当定位中心接收到各麦克风采集的信号后，按S2分别进行特征提取，构成新的测量向量其中[·]^T表示转置计算，从而稀疏定位模型可表示为：

y＝Φx+v(5)

其中，x为N×1的稀疏向量，Φ为一个M×N的矩阵，表示冗余字典，v为M×1的向量，表示噪声干扰；Φ中的任一项1≤i≤M,1≤j≤N,表示第i个麦克风收到声源在第j个格点处发出声波信号按S2进行提取后获得的特征量；

一旦将定位区域内划分为N个格点，则声源在空间上的位置可精确地用一个N×1的稀疏索引向量x表示，声源位置所对应格点处索引值为1，而其他格点对应索引值为0，即

x＝[0,1,…,0,1,0…0]^T(6)

如此一来，定位问题就转变为依据接收信号判断稀疏向量x中非零值所在位置的问题；

S4、模型失配修正：修正式为：

y＝(Φ+Γ)x+v＝Hx+v(7)

其中v为M×1的向量，表示噪声干扰；H＝(Φ+Γ)表示真实的冗余字典，其中Γ是预先未知的；

采用字典学习：

$\begin{array}{cc}\arg \underset{H,X}{min}\left|\right|Y-HX|{|}_F^2& s.t.\left|\right|x_i|{|}_0\≤K,i=1,...,L\end{array}---\left(8\right)$

其中‖·‖_F表示Frobenius范数，Y＝[y₁y₂…y_L]表示训练样本，X＝[x₁x₂…x_L]为稀疏矩阵，其分量x_i(i＝1,…,L)表示对应训练向量y_i的稀疏向量；

采用正则化方式进行字典更新，计算公式如下：

H＝YX^T(XX^T+βI)^-1(9)

其中I表示单位矩阵，β为正则化系数；

S5、稀疏重构：在稀疏恢复阶段，字典H保持不变，通过重构算法计算稀疏向量；根据CS理论，准确描述稀疏约束的是l₀范数，即

min‖x_i‖₀s.t.y_i＝Hx_i,i＝1,…,L(10)

由于l₀范数的求解是NP难问题，因此常用l₁范数来代替l₀范数进行稀疏约束，即

min‖x_i‖₁s.t.y_i＝Hx_i,i＝1,…,L(11)

用如下组合函数来作为l₀范数的近似函数，组合函数表达式如下：

$f_{\mathrm{\σ}}\left(x\right)=\mathrm{\λ}(1-\frac{{\mathrm{\σ}}^2}{x^2+{\mathrm{\σ}}^2})+(1-\mathrm{\λ})(1-e^{-x^2/{\mathrm{\λ\σ}}^2})---\left(12\right)$

其中0<λ<1；显然，f_σ(x)是连续函数，并且具有如下性质

$\underset{\mathrm{\&sigma;}\&RightArrow;0}{\lim }{f}_{\mathrm{\&sigma;}}\left(x\right)=\left\{\begin{array}{cc}1& x\&NotEqual;0\\ 0& x=0\end{array}\right.---\left(13\right)$

令 $F_{\mathrm{\&sigma;}}\left(x_i\right)=\underset{j=1}{\overset{N}{\&Sigma;}}{f}_{\mathrm{\&sigma;}}\&lsqb;x_i\left(j\right)\&rsqb;,$ 则当σ很小时有

F_σ(x_i)≈‖x_i‖₀(14)

其中x_i(j)表示x_i的第j个分量；因此(10)式所示问题可以变为

$\arg \underset{x_i}{min}{F}_{\mathrm{\&sigma;}}\left(x_i\right),s.t.y_i={\mathrm{Hx}}_i---\left(15\right)$

为了进一步逼近l₀范数的约束效果，引入加权约束思想，通过加权，使重构信号中的大系数和小系数获得同等约束，此时(15)式改为

$\arg \underset{x_i}{min}{\mathrm{WF}}_{\mathrm{\&sigma;}}\left(x_i\right),s.t.y_i={\mathrm{Hx}}_i---\left(16\right)$

其中W为N维对角加权矩阵，该矩阵仅对角线上元素不为零，其它元素绝为零，其对角线上元素取为w_j＝1/(|x_i(j)|+η)，η是一个小量，防止w_j出现奇异值；

最小化式(16)的问题可以通过拉格朗日乘子法转化为无约束最优化问题进行求解，即

$\underset{x_i}{min}L\left(x_i\right)={\mathrm{WF}}_{\mathrm{\&sigma;}}\left(x_i\right)+\mathrm{\&gamma;}\left|\right|y_i-{\mathrm{Hx}}_i|{|}_2^2---\left(17\right)$

其中γ为约束参数；

当参数σ较小时问题(17)的解容易收敛于局部最优解，为了使问题(17)能够尽量收敛于全局最优解，本发明取参数σ为下列一组下降序列

σ＝[σ₁,…,σ_T](18)

其中σ₁为较大正常数，σ_T为接近于零的正常数；迭代求解过程中将参数σ＝σ_m-1(m＝2,…,T)时得到的式(17)的最优解作为参数σ＝σ_m时求解最优化问题(17)的初始解，从而使算法逐步收敛于参数为σ＝σ_T时的全局最优解。

进一步的，S1中采用分步格点划分技术，具体步骤如下：

用较大的间隔将定位空间均匀划分为N₁个格点，并实现粗定位，得到定位目标对应的格点位置；然后将粗定位结果所对应格点空间用较小的间隔再次均匀划分为N₂个格点，并在这N₂个格点上再次实现稀疏重构，此结果作为最后的输出结果。

优选的，S2中所加的窗为汉明窗。

本发明的有益效果是：本发明采用基于两步DCT变换的特征提取方法将计算复杂度从样本长度级降低到帧长数量级，这样可以大大节省计算量，并且利用分帧获得的短时平稳特性能够提高声源特征的抗噪性；利用模型修正技术动态调整字典，克服稀疏模型与实际信号之间的失配问题，增强稀疏模型的鲁棒性，提高在混响条件下的适用性。近似l₀范数稀疏重构算法，通过构造合理的l₀范数逼近函数，直接在l₀范数意义下进行稀疏重构，可以有效提高噪声条件下的重构精度并具有较快的重构速度，不仅适用于单目标定位，而且适用于多目标定位。较传统的声源定位方法而言，本发明对阵列结构无特殊要求、计算量少，有效地提高了抗噪声和抗混响能力以及定位精度，更适合用于室内声源定位，可广泛应用于视频会议系统、语音识别系统以及智能机器人等各个领域。

附图说明

附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。在附图中：

图1是本发明声源定位方法的结构框图；

图2是本发明实施例中逼近函数的比较图(σ＝0.1,λ＝1/16)；

图3是本发明实施例中定位环境示意图；

图4是本发明实施例中无混响条件下三种算法在不同信噪比条件下的RMSE对比图；

图5是本发明实施例中混响时间对定位精度的影响。

具体实施方式

如图1所示，一种基于空间稀疏性的分布式麦克风阵列声源定位方法，包括如下步骤：

S1、定位系统建立：分布式麦克风定位系统由M个已知自身位置的麦克风和K个待定位声源组成，并将整个定位区域均匀划分为若干个格点。每个麦克风分别接收声源发出的信号，并传送给定位中心；

本发明采用分步格点划分技术，具体步骤如下：首先用较大的间隔将定位空间均匀划分为N₁个格点，并实现粗定位，得到定位目标对应的格点位置；然后将粗定位结果所对应格点空间用较小的间隔再次均匀划分为N₂个格点，并在这N₂个格点上再次实现稀疏重构，此结果作为最后的输出结果。这样不仅可以有效减少计算量，而且可以避免因格点划分过多而造成冗余字典各列之间相关性变大，从而导致稀疏恢复性能下降。

S2、语音信号特征提取：

定位中心将任一麦克风接收到的语音信号首先经过加窗(采用汉明窗)处理，把原来长度为P×1的语音信号向量r分解为J个长度为Q×1的短帧，即：

其中r_i(i＝1,…,P)表示输入信号向量r的每一个分量，Z是一个维数为J×Q的矩阵，其中每一行代表经过加窗处理后的1帧数据。

接下来对这J帧信号分别进行一次DCT变换：

式中D(·)表示DCT变换，表示z_i(i＝1,…,J)经过DCT变换后的结果。对变化后的每一帧数据通过除以该帧的最大值实现归一化处理，然后对每一帧数据进行求平均计算：

尽管语音信号一般来说在长时间内是非平稳的，但经过上述分帧处理后，时间片段相对较短，可以认为每一帧数据是短时平稳的，因此每一帧数据的统计特征可近似认为具有短时不变特性，而噪声特别是高斯噪声在求平均后，均值一般是减小的。因此，这样做不仅可以降低模型维数，而且具有一定的抑制噪声作用。接下来考虑连续多帧信号的长时特性，为此对求平均后的向量s再进行一次DCT变换，得到

注意此时向量的长度仍为J×1，也即基于两步DCT变换的特征提取方法将计算复杂度从样本长度级降低到帧长数量级，这样可以大大节省稀疏重构的计算量，并且经过再次DCT变换后，中大系数相对集中，有利于稀疏重构的实现。

S3、建立稀疏定位模型：

当定位中心接收到各麦克风采集的信号后，按S2分别进行特征提取，构成新的测量向量其中[·]^T表示转置计算，从而稀疏定位模型可表示为：

y＝Φx+v(5)

其中x为N×1的稀疏向量，Φ为一个M×N的矩阵，表示冗余字典，v为M×1的向量，表示噪声干扰；Φ中的任一项1≤i≤M,1≤j≤N,表示第i个麦克风收到声源在第j个格点处发出声波信号按S2进行提取后获得的特征量。由于格点一旦划定，格点位置为已知，再加上麦克风的位置是事先已知的，所以冗余字典Φ中的每个元素所对应的语音波形可根据文献(CevherV,BaraniukR.Compressivesensingforsensorcalibration.5thIEEESensorArrayandMultichannelSignalProcessingWorkshop.2008.175-178.)中声音信号传播模型计算得到，再经过S2就可以得到特征量。

由于声源仅在定位区域中的一个或几个格点处(或其附近)，也即在某一特定时间声源所在的位置在空间域上是稀疏的。于是，一旦定位区域内的格点位置划定，声源在空间上的位置可精确地用一个N×1的稀疏索引向量x表示，声源位置所对应格点处索引值为1，而其他格点对应索引值为0，即

x＝[0,1,…,0,1,0…0]^T(6)

如此一来，定位问题就转变为依据接收信号判断稀疏向量x中非零值所在位置的问题。

S4、模型失配修正：

以上稀疏定位模型是根据声音信号在自由空间的传播模型得到的，忽略了实际定位环境中干扰因素的影响。事实上，实际测量中不仅含有噪声，而且还存在由于墙壁、天花板等障碍物对声音的反射、散射所引起的混响效应。为了表示理想稀疏模型与实际情况之间的差别，本文用矩阵Γ表示模型误差，因此实际情况下(5)式修正为：

y＝(Φ+Γ)x+v＝Hx+v(7)

其中v为M×1的向量，表示噪声干扰；H＝(Φ+Γ)表示真实的冗余字典，其中Γ是预先未知的。显然，实际模型不仅含有噪声的影响，而且真实字典H和理想字典Φ之间也存在着明显差别。为了解决这一问题，本发明采用字典学习技术，根据接收信号不断调整冗余字典，使之能够自适应实际信号的动态改变，克服稀疏模型与实际信号之间的失配问题，使得稀疏重构结果能够精确得出位置信息。字典学习一般可以归结为如下优化问题：

$\begin{array}{cc}\arg \underset{H,X}{min}\left|\right|Y-HX|{|}_F^2& s.t.\left|\right|x_i|{|}_0\&le;K,i=1,...,L\end{array}---\left(8\right)$

其中‖·‖_F表示Frobenius范数，Y＝[y₁y₂…y_L]表示训练样本，X＝[x₁x₂…x_L]为稀疏矩阵，其分量x_i(i＝1,…,L)表示对应训练向量y_i的稀疏向量。

字典学习一般包括稀疏恢复和字典更新两个部分，并且两个部分交替进行，即在字典更新时固定稀疏向量不变，而在稀疏恢复时，采用上一步已更新过的字典。因此在模型失配修正阶段，本发明固定稀疏向量不变，采用正则化方式进行字典更新，计算公式如下：

H＝YX^T(XX^T+βI)^-1(9)

其中I表示单位矩阵，β为正则化系数。

S5、稀疏重构：

在稀疏恢复阶段，字典H保持不变，通过重构算法计算稀疏向量。根据CS理论，准确描述稀疏约束的是l₀范数，即

min‖x_i‖₀s.t.y_i＝Hx_i,i＝1,…,L(10)

由于l₀范数的求解是NP难问题，因此常用l₁范数来代替l₀范数进行稀疏约束，即

min‖x_i‖₁s.t.y_i＝Hx_i,i＝1,…,L(11)

虽然l₁范数优化问题是一个严格凸问题，可以求得其全局唯一解，但该问题只有在一定条件下才与l₀范数问题等价，也即此时所求得的解才是原问题的解。在很多实际应用中，特别是在低信噪比情况下，采用l₁范数逼近l₀范数仍存在一定的误差。正如式(10)所示的l₀范数约束模型，对于任何一个向量，其l₀范数仅存在0、1之分，因此大系数和小系数对目标函数的贡献是相同的。然而在l₁范数约束模型中，求解对象为模值最小情况下对应的解，此时向量x_i中大系数和小系数对目标函数的贡献是不一样的，即大系数模值大，相应贡献大，反之亦然。因此在求解过程中，目标函数会对大系数施加更多约束以保证整个代价函数的收敛性。然而在定位应用中，大系数往往对应于定位目标，而小系数可能对应于噪声。由于l₁范数约束模型对待大、小系数的不公平性，则可能在优化过程削弱大系数的贡献，而对小系数没有施加过多约束，从而导致重构准确性下降，在低信噪比时甚至会将噪声误判为定位目标。

针对此问题，用如下组合函数来作为l₀范数的近似函数，既克服直接用l₀范数进行约束难以计算的问题，又可以避免使用l₁范数造成的对大、小系数约束不公平性问题，从而得到近似l₀(approximatel₀norm,AL0)范数逼近算法。

组合函数表达式如下：

$f_{\mathrm{\&sigma;}}\left(x\right)=\mathrm{\&lambda;}(1-\frac{{\mathrm{\&sigma;}}^2}{x^2+{\mathrm{\&sigma;}}^2})+(1-\mathrm{\&lambda;})(1-e^{-x^2/{\mathrm{\&lambda;\&sigma;}}^2})---\left(12\right)$

其中0<λ<1。显然，f_σ(x)是连续函数，并且具有如下性质

$\underset{\mathrm{\&sigma;}\&RightArrow;0}{\lim }{f}_{\mathrm{\&sigma;}}\left(x\right)=\left\{\begin{array}{cc}1& x\&NotEqual;0\\ 0& x=0\end{array}\right.---\left(13\right)$

因此，该组合函数不仅符合文献(MohimaniH,BabaieZM,JuttenC.AFastApproachforOvercompleteSparseDecompositionBasedonSmoothedl₀Norm.IEEETransactionsonSignalProcessing,2009,57(1):289-301.)中指出的平滑函数所需满足的条件，而且如图2所示，该函数与标准高斯函数相比在x＝0附近具有更大的陡峭性，因此其对于l₀范数的逼近更精确。

令 $F_{\mathrm{\&sigma;}}\left(x_i\right)=\underset{j=1}{\overset{N}{\&Sigma;}}{f}_{\mathrm{\&sigma;}}\&lsqb;x_i\left(j\right)\&rsqb;,$ 则当σ很小时有

F_σ(x_i)≈‖x_i‖₀(14)

其中x_i(j)表示x_i的第j个分量。因此，(10)式所示问题可以变为

$\arg \underset{x_i}{min}{F}_{\mathrm{\&sigma;}}\left(x_i\right),s.t.y_i={\mathrm{Hx}}_i---\left(15\right)$

为了进一步逼近l₀范数的约束效果，引入加权约束思想，通过加权，使重构信号中的大系数和小系数获得同等约束，此时(15)式改为

$\arg \underset{x_i}{min}{\mathrm{WF}}_{\mathrm{\&sigma;}}\left(x_i\right),s.t.y_i={\mathrm{Hx}}_i---\left(16\right)$

其中W为N维对角加权矩阵，该矩阵仅对角线上元素不为零，其它元素绝为零，其对角线上元素取为w_j＝1/(|x_i(j)|+η)，η是一个小量，防止w_j出现奇异值。显然，当重构系数x_i(j)较小时，所对应的权值w_j较大，反之亦然。如此一来，刚好可以拉近大、小系数在目标函数贡献性上的差距，因此加权过程相当于在优化过程中实现了对大、小系数的平衡“惩罚”。

最小化式(16)的问题可以通过拉格朗日乘子法转化为无约束最优化问题进行求解，即

$\underset{x_i}{min}L\left(x_i\right)={\mathrm{WF}}_{\mathrm{\&sigma;}}\left(x_i\right)+\mathrm{\&gamma;}\left|\right|y_i-{\mathrm{Hx}}_i|{|}_2^2---\left(17\right)$

其中γ为约束参数。综合考虑计算速度和计算复杂度，本发明选用最优化理论中的BFGS算法进行求解，该算法可以获得比最陡梯度法更快的收敛速度，又没有牛顿法那样大的计算复杂度。值得注意的是，当参数σ较小时问题(17)的解容易收敛于局部最优解。为了使问题(17)能够尽量收敛于全局最优解，本发明取参数σ为下列一组下降序列

σ＝[σ₁,…,σ_T](18)

其中σ₁为较大正常数，σ_T为接近于零的正常数。迭代求解过程中将参数σ＝σ_m-1(m＝2,…,T)时得到的式(17)的最优解作为参数σ＝σ_m时求解最优化问题(17)的初始解，从而使算法逐步收敛于参数为σ＝σ_T时的全局最优解。

为了更好地理解本发明的技术方案，以下将结合附图及具体实施例对本发明的工作流程及有益效果进行详细说明。

在本实施例中，声源信号选择NTT中文数据库中的语音信号，截取长度L＝65536的语音片段，分成J＝128帧，每帧长度为512，采样频率F_s＝16kHz。实验中采用Image模型产生不同混响时间下的房间冲激响应，进而通过声源信号与房间冲激响应的卷积得到麦克风的接收信号，噪声采用Noisex92噪声库的白噪声。选择文献(CevherV,BaraniukR.Compressivesensingforsensorcalibration.5thIEEESensorArrayandMultichannelSignalProcessingWorkshop.2008.175-178.)中基于l₁范数约束的定位方法(称为对比方法1)和文献(JiangH,MathewsB,WilfordP.Soundlocalizationusingcompressivesensing.Proceedingsofthe1stInternationalConferenceonSensorNetworks.2012.159-166.)中基于正交匹配追踪(OrthogonalMatchingPursuit，OMP)的定位方法(称为对比方法2)进行性能对比。每种算法均重复进行50次实验，每次实验在定位区域内随机选择若干个目标位置。定位计算由笔记本电脑(配置为因特尔的i5CPU，3.1GHz主频，4GB内存)实现，负责处理收集数据和运行定位算法。定位性能用均方根误差(rootmeansquareerror,RMSE)来衡量，其中RMSE的计算公式为

$RMSE=\frac{1}{QK}\underset{q=1}{\overset{Q}{\&Sigma;}}\underset{i=1}{\overset{K}{\&Sigma;}}\sqrt{\left|\right|{\hat{P}}_q\left(i\right)-P_q\left(i\right)|{|}_2^2}---\left(19\right)$

其中Q为测试次数，K为声源个数，为第i个目标第q次测量时的估计位置，P_q(i)为第i个目标第q次测量时的真实位置。RMSE可以衡量不同定位算法的平均定位精度，但在多目标定位时有时会出现多数目标的定位精度很高，而因1、2个目标定位误差较大拉低RMSE的情况。为了公平起见，引入定位失效率指标，并且规定当目标估计位置与真实位置直接距离相差大于0.3m时，定位失效。定位失效率用下式来表示

考虑三维定位效果，设房间大小为6m×6m×2.5m，8个麦克风高度均为1m，在xy平面上等间隔放在定位区域的边缘，即相邻麦克风之间距离均为3m(如图3所示)。为了实现压缩感知定位，将定位空间均匀划分为若干格点。首先进行粗格点方式划分，在xyz三个方向均按0.5m的间隔进行划分，也即对应N₁＝720；在粗定位完成后，在目标对应的0.5m×0.5m×0.5m大小格点范围内再按xyz三个方向均为0.1m精度进行划分，此时N₂＝125。参数σ的下降序列取为[1,0.5,0.2,0.1,0.05,0.02,0.01，0.001]，λ＝1/16。设置目标数为3，目标位置随机选取。首先考虑无混响条件下的定位性能。表1给出了无混响时三种方法的定位结果对比，信噪比为15dB。

表1无混响条件下三种方法定位性能比较

从表上可以看出，在无混响且信噪比较高的情况下，三种定位方法的定位误差均在可接受范围内。但本发明方法的定位误差要小于其它两种方法，RMSE为0.12m，接近格点划分的最小精度，失效率也低于5％。对比方法1的失效率与本发明方法接近，这说明基于l₁范数约束的方法能够保证大多数目标的定位精度在小于0.3m的有效范围内，但对比方法1的均方根误差要高于本发明方法，说明其失效目标的偏离误差是比较大的。对比方法2采用贪婪算法进行稀疏重构，由于OMP算法的抗噪声能力相对较弱，因此对比方法2的定位精度最低。

图4给出了无混响条件下三种方法在不同信噪比条件下的RMSE对比图。如图4所示，由于经过两步DCT处理得到的信号特征对噪声不敏感，再加上重构时改进的平滑l₀范数重构算法采用连续函数来逼近l₀范数约束，克服了l₁范数约束模型对待大、小系数的不公平性，从而能较好地减轻噪声的影响，确保重构的准确性，所以本发明方法的定位性能要优于其它两种算法。在低信噪比时定位性能的改善尤为明显，能将定位精度提高约0.2m。当信噪比大于20dB时，本发明方法的定位误差基本不再变化，且定位精度都在0.1m左右。

接下来分析混响条件下的定位性能，混响时间设为250ms，信噪比为15dB，其它条件不变。定位结果如表2所示，对比表2可以看出在混响条件下三种方法的定位性能都有不同程度的下降，尤其对比方法1和对比方法2在混响条件下的RMSE分别要比无混响时增大0.12m和0.16m。这是因为根据声音信号的传播模型构建的字典与混响条件下的实际接收信号之间存在失配情况，从而导致稀疏重构精度大大下降。而本发明方法中由于采用在线字典学习方式实时调整字典使之与实际接收信号相适应，所以定位性能下降最少，失效率仍能控制在5％之内，这意味着即使在混响条件下采用本发明方法绝大多数时候也能满足有效定位的需求。

表2混响条件下三种方法定位性能比较(混响时间250ms)

进一步，图5给出了不同混响条件下三种方法定位性能的对比情况，信噪比仍为15dB。从图上可以看出，对比方法1和对比方法2的均方根误差随着混响时间的增大而增加，而本发明方法的均方根误差的变化则很轻微。这是因为随着混响时间增大，对比方法1和对比方法2中字典的失配情况变得越发严重，而本发明方法通过字典学习调整字典与实际情况相匹配，从而可以减轻混响的影响。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (3)

1.一种基于空间稀疏性的分布式麦克风阵列声源定位方法，其特征在于，包括如下步骤：

S1、定位系统建立：分布式麦克风定位系统由M个已知自身位置的麦克风和K个待定位声源组成，采用分步格点划分方法将整个定位区域均匀划分为若干个格点；每个麦克风分别接收声源发出的信号，并传送给定位中心；

S2、语音信号特征提取：定位中心将任一麦克风接收到的语音信号首先经过加窗处理，把原来长度为P×1的语音信号向量r分解为J个长度为Q×1的短帧，即：

其中，r_i(i＝1,…,P)表示输入信号向量r的每一个分量，Z是一个维数为J×Q的矩阵，每一行代表经过加窗处理后的1帧数据；

接下来对这J帧信号分别进行一次DCT变换：

式中D(·)表示DCT变换，(i＝1,…,J)表示z_i(i＝1,…,J)经过DCT变换后的结果；对变化后的每一帧数据通过除以该帧的最大值实现归一化处理，然后对每一帧数据进行求平均计算：

接下来考虑连续多帧信号的长时特性，为此对求平均后的向量s再进行一次DCT变换，得到

向量的长度仍为J×1，也即基于两步DCT变换的特征提取方法将计算复杂度从样本长度级降低到帧长数量级；

S3、建立稀疏定位模型：当定位中心接收到各麦克风采集的信号后，按S2分别进行特征提取，构成新的测量向量其中[·]^T表示转置计算，从而稀疏定位模型可表示为：

y＝Φx+v(5)

其中，x为N×1的稀疏向量，Φ为一个M×N的矩阵，表示冗余字典，v为M×1的向量，表示噪声干扰；Φ中的任一项1≤i≤M,1≤j≤N,表示第i个麦克风收到声源在第j个格点处发出声波信号按S2进行提取后获得的特征量；

一旦将定位区域内划分为N个格点，则声源在空间上的位置可精确地用一个N×1的稀疏索引向量x表示，声源位置所对应格点处索引值为1，而其他格点对应索引值为0，即

x＝[0,1,…,0,1,0…0]^T(6)

如此一来，定位问题就转变为依据接收信号判断稀疏向量x中非零值所在位置的问题；

S4、模型失配修正：修正式为：

y＝(Φ+Γ)x+v＝Hx+v(7)

其中v为M×1的向量，表示噪声干扰；H＝(Φ+Γ)表示真实的冗余字典，其中Γ是预先未知的；

采用字典学习：

$\begin{array}{cc}\arg \underset{H,X}{min}\left|\right|Y-HX|{|}_F^2& s.t.\left|\right|x_i|{|}_0\&le;K,i=1,...,L\end{array}---\left(8\right)$

其中||·||_F表示Frobenius范数，Y＝[y₁y₂…y_L]表示训练样本，X＝[x₁x₂…x_L]为稀疏矩阵，其分量x_i(i＝1,…,L)表示对应训练向量y_i的稀疏向量；

采用正则化方式进行字典更新，计算公式如下：

H＝YX^T(XX^T+βI)^-1(9)

其中I表示单位矩阵，β为正则化系数；

S5、稀疏重构：在稀疏恢复阶段，字典H保持不变，通过重构算法计算稀疏向量；根据CS理论，准确描述稀疏约束的是l₀范数，即

min||x_i||₀s.t.y_i＝Hx_i,i＝1,…,L(10)

由于l₀范数的求解是NP难问题，因此常用l₁范数来代替l₀范数进行稀疏约束，即

min||x_i||₁s.t.y_i＝Hx_i,i＝1,…,L(11)

用如下组合函数来作为l₀范数的近似函数，组合函数表达式如下：

$f_{\mathrm{\&sigma;}}\left(x\right)=\mathrm{\&lambda;}(1-\frac{{\mathrm{\&sigma;}}^2}{x^2+{\mathrm{\&sigma;}}^2})+(1-\mathrm{\&lambda;})(1-e^{-x^2/{\mathrm{\&lambda;\&sigma;}}^2})---\left(12\right)$

其中0<λ<1；显然，f_σ(x)是连续函数，并且具有如下性质

$\underset{\mathrm{\&sigma;}\&RightArrow;0}{\lim }{f}_{\mathrm{\&sigma;}}\left(x\right)=\left\{\begin{array}{cc}1& x\&NotEqual;0\\ 0& x=0\end{array}\right.---\left(13\right)$

令 $F_{\mathrm{\&sigma;}}\left(x_i\right)=\underset{j=1}{\overset{N}{\&Sigma;}}{f}_{\mathrm{\&sigma;}}\&lsqb;x_i\left(j\right)\&rsqb;,$ 则当σ很小时有

F_σ(x_i)≈||x_i||₀(14)

其中x_i(j)表示x_i的第j个分量；因此(10)式所示问题可以变为

$\arg \underset{x_i}{min}{F}_{\mathrm{\&sigma;}}\left(x_i\right),s.t.y_i={\mathrm{Hx}}_i---\left(15\right)$

为了进一步逼近l₀范数的约束效果，引入加权约束思想，通过加权，使重构信号中的大系数和小系数获得同等约束，此时(15)式改为

$\arg \underset{x_i}{min}{\mathrm{WF}}_{\mathrm{\&sigma;}}\left(x_i\right),s.t.y_i={\mathrm{Hx}}_i---\left(16\right)$

其中W为N维对角加权矩阵，该矩阵仅对角线上元素不为零，其它元素绝为零，其对角线上元素取为w_j＝1/(|x_i(j)|+η)，η是一个小量，防止w_j出现奇异值；

最小化式(16)的问题可以通过拉格朗日乘子法转化为无约束最优化问题进行求解，即

$\underset{x_i}{min}L\left(x_i\right)={\mathrm{WF}}_{\mathrm{\&sigma;}}\left(x_i\right)+\mathrm{\&gamma;}\left|\right|y_i-{\mathrm{Hx}}_i|{|}_2^2---\left(17\right)$

其中γ为约束参数；

当参数σ较小时问题(17)的解容易收敛于局部最优解，为了使问题(17)能够尽量收敛于全局最优解，本发明取参数σ为下列一组下降序列

σ＝[σ₁,…,σ_T](18)

其中σ₁为较大正常数，σ_T为接近于零的正常数；迭代求解过程中将参数σ＝σ_m-1(m＝2,…,T)时得到的式(17)的最优解作为参数σ＝σ_m时求解最优化问题(17)的初始解，从而使算法逐步收敛于参数为σ＝σ_T时的全局最优解。

2.根据权利要求1所述的基于空间稀疏性的分布式麦克风阵列声源定位方法，其特征在于，S1中采用分步格点划分技术，具体步骤如下：

用较大的间隔将定位空间均匀划分为N₁个格点，并实现粗定位，得到定位目标对应的格点位置；然后将粗定位结果所对应格点空间用较小的间隔再次均匀划分为N₂个格点，并在这N₂个格点上再次实现稀疏重构，此结果作为最后的输出结果。

3.根据权利要求1所述的基于空间稀疏性的分布式麦克风阵列声源定位方法，其特征在于，S2中所加的窗为汉明窗。