CN102969000B - 一种多通道语音增强方法 - Google Patents

Abstract

本发明公开了一种F范数表征信号子空间维度的多通道语音增强方法，该方法包括：步骤1：通过N个麦克风组成的麦克风阵列采集带噪声的多路语音信号y(t)，并计算该多路语音信号y(t)的带噪语音互相关矩阵R_yy，其中，t表示离散时间点；步骤2：利用噪声估计算法估计加性噪声互相关矩阵 步骤3：利用所述带噪语音互相关矩阵R_yy和加性噪声互相关矩阵估计纯净语音互相关矩阵 步骤4：利用所述纯净语音互相关矩阵估计信号子空间的维度 步骤5：对进行广义特征值分解，并结合信号子空间的维度和拉格朗日乘子μ获得时域约束线性信号估计器；步骤6：利用时域约束线性信号估计器对多路语音信号y(t)进行滤波获得增强后的语音

Description

一种多通道语音增强方法

技术领域

本发明涉及语音增强和噪声消除的语音信号处理领域，特别涉及一种多通道语音增强方法。

背景技术

信号子空间方法(Signal Subspace Approach，SSA)主要基于矩阵正交分解理论，利用特征值分解(Eigenvalue Decomposition，EVD或奇异值分解(Singular Value Decomposition，SVD)，将带噪语音信号的向量空间分解为彼此正交的两部分：信号子空间(信号+噪声)和噪声子空间(只有噪声)。其中，EVD利用信号的自相关矩阵，实现由时域到特征域的变换，被视为依赖于信号的变换，其经常被称作Karhunen-Loeve变换(KLT)。本质上，基于SSA的语音增强是将噪声子空间置零，同时去除信号子空间中的噪声成分。

SSA首先被引入到单通道语音增强领域(参考文献1：Y.Ephraim andH.L.Van Trees，“A signal subspace approach for speech enhancement，”IEEETrans.Speech Audio Process.，vol.3，no.4，pp.251-266，Jul.1995)，在Karhunen-Loeve域上设计了多种有约束的线性滤波器，很好的处理了白噪声污染的语音信号。Hu等人(参考文献2：Y.Hu and P.C.Loizou，“Ageneralized subspace approach for enhancing speech corrupted by colorednoise，”IEEE Trans.Speech Audio Process.，vol.11，no.4，pp.334-341，Jul.2003)利用广义特征值分解(Generalized Eigenvalue Decomposition，GEVD)的方法，提高了SSA处理有色噪声的能力。程宁等人(参考文献3：N.Cheng，W.Liu and L.Wang，“Masking Property Based Microphone Array Post-filterDesign”，Interspeech，2010)将人耳的听觉掩蔽效应应用到了SSA中，精确了约束条件，使算法性能进一步得到提升。

Doclo等人(参考文献4：S.Doclo and M.Moonen，“GSVD-basedoptimal filtering for single and multimicrophone speech enhancement，”IEEETrans.Signal Process.，vol.50，no.9，pp.2230-2244，Sep.2002)将SSA应用到多通道语音增强算法中，它以广义奇异值分解(Generalized SingularValue Decomposition，GSVD)为工具，推导了在多通道环境下SSA的最优线性滤波器。Kim等人为提高语音识别系统的性能，先后提出了子频带独立插值方法(参考文献5：G.Kim and N.I.Cho，Principal subspacemodification for multi-channel wiener filter in multi-microphone noisereduction，Proc.IEEE Int.Conference on Acoustics，Speech and SignalProcessing，2008)和将多干扰源分离到不同的子空间中的子空间建模方法(参考文献6：G.Kim and P.C.Loizou，Improving Speech Intelligibility inNoise Using Environment-Optimized Algorithms，IEEE Trans.Audio，Speech，And Language Process.，vol.18，no.8，pp.2080-2090，2010)。

多通道SSA也存在着诸多缺陷，主要包括：在实际环境中应用KLT时，信号子空间维度的难以精确估计；处理有色噪声时的广义特征向量矩阵并不是正交阵，因此在特征值域上设计的线性信号估计器会引入语音畸变；多通道条件下GEVD的计算复杂度极高，不利于SSA在实时系统中应用，而现有的快速算法并不能保证算法的性能。

从先验知识的角度来审视SSA，其基于这样一个假设：一个长为p的随机序列可以由长为Q＜p的随机序列的线性变换得到，称之为低秩模型。由自回归(AutoRegressive，AR)模型模型和线性预测系数(Linear PredictionCoefficients，LPC)等知识，我们知道这个假设对于语音信号是成立的。但将这个假设应用到语音增强算方法中，却必须非常小心，对信号子空间(低秩模型)的维度必须精确估计。一方面，对信号子空间维度的过估计会导致信号子空间中噪声成分增加，从而在增强后的语音中保留了过多残余噪声；另一方面，对信号子空间的欠估计又会导致语音成分的丢失，从而增大语音畸变，降低语音感知质量。因此，精确的信号子空间维度估计是SSA优异性能的保证，并直接影响到算法后续过程的有效性。但在实际应用，由于房间回声和非平稳噪声等因素的影响，信号子空间的维度往往难以精确估计，因此低秩模型的先验知识难以准确运用。

信号子空间维度估计的传统方法的根本缺陷在于过分依赖噪声估计，而忽略了信号子空间维度自身的变化规律。当基于SSA的语音增强算法工作在真实环境中时，由于噪声的干扰，SVD或EVD都会产生大量接近于零的奇异值或特征值，而且其值变化具有连续性，这也信号子空间维度的确定增加了难度。

发明内容

(一)要解决的技术问题

针对传统的信号子空间维度估计方法在真实环境中难以有效工作，其原因是过分依赖噪声的估计，本发明充分利用语音离散傅里叶变换(Discrete Fourier Transform，DFT)数幅值的超高斯分布特性和麦克风之间的相关性，提出了一种新的信号子空间维度估计方法，目的在于精确估计信号子空间的维度。

(二)技术方案

本发明为解决上述问题而提出的一种F范数表征信号子空间维度的多通道语音增强方法，其包括下列步骤：

步骤1：通过N个麦克风组成的麦克风阵列采集带噪声的多路语音信号y(t)，并计算该多路语音信号y(t)的带噪语音互相关矩阵R_yy，其中，t表示离散时间点；

步骤2：利用噪声估计算法估计加性噪声互相关矩阵

步骤3：利用所述带噪语音互相关矩阵R_yy和加性噪声互相关矩阵估计纯净语音互相关矩阵

步骤4：利用所述纯净语音互相关矩阵估计信号子空间的维度

步骤5：对进行广义特征值分解，并结合信号子空间的维度和拉格朗日乘子μ获得时域约束线性信号估计器；

步骤6：利用时域约束线性信号估计器对多路语音信号y(t)进行滤波获得增强后的语音

(三)有益效果

本发明充分利用了语音DFT系数幅值的超高斯分布特性和麦克风之间的相关性。首先利用超高斯分布模型来更好地拟合各个频点上语音幅值谱的统计柱状图。其次，利用麦克风之间的相关性和正交变换不改变矩阵F范数的特点，获得了语音互相关矩阵F范数的高斯分布模型。然后，利用一种基于最大化原则的估计策略，在接受原假设的前提下最大化信号子空间的维度。通过以上步骤，本发明能够克服真实环境中噪声能量波动的不利影响，更精确地估计信号子空间的维度，使SSA可以在语音畸变和噪声消除之间建立更合理的折中。

附图说明

图1是本发明中F范数表征信号子空间维度的多通道语音增强方法流程图；

图2是本发明中基于F范数的估计信号子空间维度的方法流程图；

图3是本发明中TDC线性信号估计器的设计流程图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明作进一步的详细说明。

本发明提供一种基于F范数的信号子空间维度估计方法，该方法的具体步骤如下：

步骤S1：通过N个麦克风组成的麦克风阵列采集多路带噪语音信号y(t)，具体包括：

步骤S11：采集信号所用的麦克风阵列，其由N个等间距的麦克风组成。假设第n个麦克风在时刻t采集到的带噪语音信号y_n(t)为x_n(t)与v_n(k)之和，即

y_n(t)＝g_n*s(t)+v_n(t)＝x_n(t)+v_n(t)    (1)

其中，g_n是声源s(t)到第n个麦克风的冲击响应，x_n(t)是第n个麦克风采集到的纯净语音，v_n(t)是第n个麦克风采集到的加性噪声，t表示离散时间点。

步骤S12：利用长为L的汉宁窗对信号进行截短。信号模型可表示为向量形式：

y_n＝x_n+v_n，n＝1，2，…，N    (2)

其中，y_n＝[y_n(t)y_n(t-1)...y_n(t-L+1)]^T，其包含了带噪语音信号y_n(t)最近的L个采样，v_n和x_n类比定义。那么，所有麦克风采集到的带噪语音信号的全局向量形式为：

y＝x+v    (3)

其中，

$y={\left[\begin{array}{cccc}{y}_1^T& y_2^T& ...& y_N^T\end{array}\right]}^T$

$x={\left[\begin{array}{cccc}{x}_1^T& x_2^T& ...& x_N^T\end{array}\right]}^T---\left(4\right)$

$v={\left[\begin{array}{cccc}{v}_1^T& v_2^T& ...& v_N^T\end{array}\right]}^T$

同时，利用长度为75ms的带噪语音信号通过短时平均计算互相关矩阵R_yy，

R_yy＝E{yy^T}        (5)

其中，E{·}是求期望算子，[·]^T是转置算子，R_yy的阶数为NL×NL。

步骤S2：计算当前帧的后验信噪比(Signal-to-Noise Ratio，SNR)γ：

$\mathrm{\γ}=10\·\log \frac{R_{\mathrm{yy}}\left(\mathrm{1,1}\right)}{{\hat{R}}_{\mathrm{vv},k-1}\left(\mathrm{1,1}\right)}---\left(6\right)$

其中，log{·}是求自然对数算子，是前一帧中加性噪声互相关矩阵中第1行第1列元素的估计。

步骤S3：利用噪声估计算法，例如语音激活检测(Voice ActivityDetection，VAD)方法，更新加性噪声互相关矩阵的估计具体包括：

根据有声段的SNR的下限阈值γ_th，更新加性噪声互相关矩阵的估计值如下：

${\hat{R}}_{\mathrm{vv},k}=\left\{\begin{array}{cc}{\hat{R}}_{\mathrm{vv},k-1}& \mathrm{if}{\mathrm{\γ}}_k>{\mathrm{\γ}}_{\mathrm{th}}\\ 0.98\·{\hat{R}}_{\mathrm{vv},k-1}+R_{\mathrm{yy}}& \mathrm{else}\end{array}\right.---\left(7\right)$

其中，γ_th为SNR的下限阈值，γ_k为第k帧的SNR；所述是前一帧的加性噪声互相关矩阵的估计值，R_yy是当前帧的带噪语音的互相关矩阵，k是帧号。

步骤S4：利用麦克风间纯净语音互相关矩阵F范数的概率密度分布模型，估计信号子空间的维度具体包括：

步骤S41：初步估计纯净语音的互相关矩阵并进行特征值分解(EVD)，获得特征值和特征向量其中，q＝1，2，...，NL；具体包括：

步骤S411：利用步骤S12中带噪语音互相关矩阵R_yy和步骤S3中加性噪声的互相关矩阵估计纯净语音的互相关矩阵如下：

${\hat{R}}_{\mathrm{xx}}=R_{\mathrm{yy}}-{\hat{R}}_{\mathrm{vv},k}---\left(8\right)$

步骤S412：对进行特征值分解(EVD)，得到特征值集合 ${\widehat{\mathrm{\λ}}}_{x,q}(q=\mathrm{1,2},...,\mathrm{NL})$ 和特征向量集合 ${\hat{b}}_{x,q}(q=\mathrm{1,2},...,\mathrm{NL}).$

步骤S42：计算当前帧的纯净语音互相关矩阵的拒绝水平δ，并计算相对应的F范数阈值。具体步骤如下：

步骤S421：利用广义超高斯分布模型来拟合纯净语音离散傅里叶变换DFT系数幅值的统计柱状图，以得到纯净语音离散傅里叶变换DFT系数幅值的超高斯分布模型。

其中，所述纯净语音离散傅里叶变换DFT系数幅值的超高斯分布模型定义为：

$p\left({\mathrm{\α}}_x\right)=\frac{{\mathrm{\η}}^{\mathrm{\θ}+1}}{\mathrm{\Γ}(\mathrm{\θ}+1)}\frac{{\mathrm{\α}}_x^{\mathrm{\θ}}}{{\mathrm{\σ}}_x^{\mathrm{\θ}+1}}\exp \{-\mathrm{\η}\frac{{\mathrm{\α}}_x}{{\mathrm{\σ}}_x}\}---\left(9\right)$

其中，θ和η是可调参数，α_x是纯净语音离散傅里叶变换DFT系数的幅值，σ_x是纯净语音的标准差，Γ是Gamma函数，exp{·}是指数算子。

本发明的一个实施例中通过1小时的纯净语音数据训练得到纯净语音DFT系数幅值的统计柱状图p_o。并通过将超高斯分布模型和纯净语音DFT系数幅值的统计柱状图之间的Kullback-Leibler距离(Kullback-LeiblerDivergence，KLD)最小化，得到最佳参数组合(θ，η)。其中，KLD J(α_x：o)可如下计算：

$J({\mathrm{\α}}_x:o)=\underset{0}{\overset{M}{\∫}}(p_o\left(m\right)-p_{{\mathrm{\α}}_x}\left(m\right))\log \left(\frac{p_o\left(m\right)}{p_{{\mathrm{\α}}_x}\left(m\right)}\right)\mathrm{dm}---\left(10\right)$

其中，p_o(m)是纯净语音DFT系数幅值的统计柱状图，是所述纯净语音离散傅里叶变换DFT系数幅值的超高斯分布模型，log{·}是求自然对数算子，m是纯净语音DFT系数幅值的统计柱状图中纯净语音DFT系数幅值的可能值，其取值范围从0到M。

通过在纯净语音DFT的每个频点上最小化(10)式，我们可以得到每个频点上纯净语音DFT系数幅值的超高斯分布参数θ和η，也就是说，得到公式(9)表示的纯净语音离散傅里叶变换DFT系数幅值的超高斯分布模型的精确表达。

步骤S422：推导麦克风间互相关矩阵(也称作全局互相关矩阵)F范数与单通道纯净语音DFT系数幅值的函数关系。

由于不同麦克风的自功率谱矩阵具有相等的F范数，同时，不同组合的麦克风间的互相关矩阵的F范数相等，但小于自相关矩阵的F范数。那么可以得到每一个自/互相关矩阵的F范数为：

${\left|\right|R_{\mathrm{xx},\mathrm{ij}}\left|\right|}_f=\left\{\begin{array}{cc}{\left|\right|R_{\mathrm{xx},11}\left|\right|}_f,& \mathrm{ifi}=j\\ \sqrt{{\mathrm{\β}}_1}{\left|\right|R_{\mathrm{xx},11}\left|\right|}_f,& \mathrm{else}\end{array}\right.---\left(11\right)$

其中，R_xx，11是第一个麦克风的自相关矩阵，||·||_f是矩阵F范数算子，β₁是麦克风间相关系数。

那么全局纯净语音互相关矩阵的F范数可如下计算：

${\left|\right|R_{\mathrm{xx}}\left|\right|}_f=\sqrt{\underset{i,j=1}{\overset{N}{\mathrm{\Σ}}}{\left|\right|R_{\mathrm{xx},\mathrm{ij}}\left|\right|}_f^2}$

$=\sqrt{N(N-1){\mathrm{\β}}_1+N}{\left|\right|R_{\mathrm{xx},11}\left|\right|}_f$

$=\sqrt{[N(N-1){\mathrm{\β}}_1+N]L}\·\sqrt{\underset{\mathrm{\ω}=0}{\overset{L-1}{\mathrm{\Σ}}}{R}_{\mathrm{xx},11}^2\left(\mathrm{\ω}\right)}---\left(12\right)$

$={\mathrm{\β}}_2\sqrt{\underset{\mathrm{\ω}=0}{\overset{L-1}{\mathrm{\Σ}}}{\mathrm{\α}}_x^4\left(\mathrm{\ω}\right)}$

其中，∑{·}是求和算子，α_x(ω)是在长为L的纯净语音DFT系数中频点ω处的幅值， ${\mathrm{\β}}_2=\sqrt{[N(N-1){\mathrm{\β}}_1+N]L}.$

步骤S423：利用步骤S421中所述纯净语音DFT系数幅值的超高斯分布模型和步骤S423中推导得出的麦克风间互相关矩阵F范数表示的麦克风间的相关性，获得纯净语音互相关矩阵F范数的概率密度分布模型，并用高斯分布拟合。具体过程如下：

首先，为简化表达，定义f_x＝||R_xx||_f。那么，f_x的概率密度分布函数p(f_x)可以由(8)式中的p(α_x(ω))结合(11)式得到。但此计算的变量太多，难以得到清晰的解析解。马尔柯夫链蒙特卡洛(Markov chain Monte Carlo，MCMC)方法恰到好处地用来计算p(f_x)的统计解。通过MCMC算法产生1,000,000个服从p(α_x(ω))分布的观测值，从而得到相同数量的f_x计算结果。然后，为清晰的描述此分布，本发明将f_x对应到对数域内，即定义f_log＝log₁₀f_x。从而获得了对数域上的统计柱状图。f_log的统计柱状图具有明显的高斯特性，因此本发明利用高斯分布来拟合f_log的统计柱状图，得到纯净语音互相关矩阵F范数的高斯分布模型p_k(f_log)：

$p_k\left(f_{\log }\right)=\frac{1}{\sqrt{\mathrm{\π}}{\mathrm{\σ}}_{\log }}\exp \{-\frac{{(f_{\log }-{\widehat{\mathrm{\μ}}}_{\log }\left(k\right))}^2}{{\mathrm{\σ}}_{\log }^2}\}---\left(13\right)$

其中，f_log是多通道纯净语音互相关矩阵R_xx的F范数的对数，σ_log是f_log高斯分布模型的方差，π是圆周率，k表示当前帧号，是数学期望，如下计算：

${\widehat{\mathrm{\μ}}}_{\log }\left(k\right)=\mathrm{mean}\{\log {\left|\right|{\hat{R}}_{\mathrm{xx}}\left(k\right)\left|\right|}_f,...,\log {\left|\right|{\hat{R}}_{\mathrm{xx}}(k-D+1)\left|\right|}_f\}---\left(14\right)$

其中，log{·}是取对数算子，||·||_f是取F范数算子，mean{·}是取均值算法，D是取均值的帧数，是对纯净语音互相关矩阵R_xx的初步估计。步骤S424：计算当前帧语音互相关矩阵的拒绝水平δ，其利用了步骤S2估计的SNRγ。方法如下：

$\mathrm{\&delta;}=\left\{\begin{array}{ccc}{\mathrm{\&delta;}}_{\max }& \mathrm{for}& \mathrm{\&gamma;}<{\mathrm{\&gamma;}}_{\min }\\ \frac{{\mathrm{\&delta;}}_{\max }-{\mathrm{\&delta;}}_{\min }}{{10}^{\frac{{\mathrm{\&gamma;}}_{\max }}{10}}-{10}^{\frac{{\mathrm{\&gamma;}}_{\min }}{10}}}& \mathrm{for}& {\mathrm{\&gamma;}}_{\min }\&le;\mathrm{\&gamma;}\&le;{\mathrm{\&gamma;}}_{\max }\\ {\mathrm{\&delta;}}_{\min }& \mathrm{for}& \mathrm{\&gamma;}>{\mathrm{\&gamma;}}_{\max }\end{array}\right.---\left(15\right)$

其中，γ_max＝10dB，γ_min＝-10dB，δ_max＝0.5，δ_min＝0.05。

步骤S425：利用(13)式获得当前帧的纯净语音互相关矩阵F范数的高斯分布，并利用(15)给定的阈值，解下面的方程，获得F范数的上限阈值fδ：

${\&Integral;}_{f_{\mathrm{\&delta;}}}^{+\&infin;}{p}_k\left(f_{\log }\right){\mathrm{df}}_{\log }=\mathrm{\&delta;}---\left(16\right)$

其中，δ是拒绝水平，p_k(f_log)是F范数的高斯分布模型。

步骤S43：利用最大化估计器来估计信号子空间维度。具体为：

基于信号子空间应尽可能最大限度地包含所有语音信息的估计策略，利用最大化信号子空间维度估计器从过估计的方向逼近真实维度，并基于如下考虑：由于信号子空间维度的过估计所带来的残留噪声可以被滤波器部分消除，而低估计引起的语音畸变却难以修复。因此，信号子空间的维度可如下估计：

$\hat{Q}=\arg \underset{1\&le;q\&le;\mathrm{NL}}{\max }\{\frac{1}{2}\log \left(\underset{i=1}{\overset{q}{\mathrm{\&Sigma;}}}{\widehat{\mathrm{\&lambda;}}}_{x,i}^2\right)\&le;f_{\mathrm{\&delta;}}\}---\left(17\right)$

其中，表示对进行特征值分解(EVD)后得到的第i个特征值，argmax{·}是取最大值坐标的算法，∑[·]是有限范围求和算子，q为在1到NL之间取值的整数。

步骤S5：对进行广义特征值分解(GEVD)，并结合信号子空间的维度和拉格朗日乘子μ设计时域约束(Time Domain Constraint，TDC)线性信号估计器，具体包括：

步骤S51：设计TDC线性信号估计器。具体包括：

步骤S51：对进行广义特征值分解GEVD：

利用步骤S3中加性噪声互相关矩阵和步骤S411中纯净语音互相关矩阵进行矩阵乘法运算，得到

${\hat{R}}_{\mathrm{vx}}={\hat{R}}_{\mathrm{vv}}^{-1}{\hat{R}}_{\mathrm{xx}}---\left(18\right)$

对进行广义特征值分解GEVD，得到特征向量矩阵W和特征值矩阵Λ。

步骤S52：利用步骤S2中得到的SNR计算拉格朗日乘子μ，如下：

$\mathrm{\&mu;}=\left\{\begin{array}{ccc}5& \mathrm{for}& \mathrm{\&gamma;}\&le;-5\mathrm{dB}\\ 4.2-\frac{\mathrm{\&gamma;}}{6.25}& \mathrm{for}& -5\mathrm{dB}<\mathrm{\&gamma;}<20\mathrm{dB}\\ 1& \mathrm{for}& \mathrm{\&gamma;}\&GreaterEqual;20\mathrm{dB}\end{array}\right.---\left(19\right)$

步骤S53：设计TDC线性信号估计器：

计算阶增益矩阵G_TDC，其每个元素可如下获得：

$g_{\mathrm{TDC},i}=\frac{{\mathrm{\&lambda;}}_{\mathrm{xv},i}}{{\mathrm{\&lambda;}}_{\mathrm{xv},i}+\mathrm{\&mu;}},\mathrm{for\; i}=\mathrm{1,2},...,\hat{Q}---\left(20\right)$

其中，λ_xv，i是对进行广义特征值分解GEVD后得到的第i个特征值，μ是拉格朗日乘子。

按照TDC规则生成线性信号估计器，如下：

$H={\mathrm{UW}}_1{G}_{\mathrm{TDC}}{W}_1^T---\left(21\right)$

其中，矩阵W₁是特征向量矩阵W中的前个特征向量组成的列满秩矩阵，U＝[I_L×L0_L×L…0_L×L]是L×NL的矩阵，I_L×L是L阶单位阵。

步骤S6：利用TDC线性信号估计器对多路语音信号y(t)进行滤波，如下式所示：

$\hat{x}\left(t\right)=\mathrm{Hy}\left(t\right)---\left(22\right)$

步骤S7：对滤波后语音信号进行帧间重叠获得增强后语音，其中重叠率为75％。

本发明的实施例使用卡内基美隆大学(Carnegie Mellon University，CMU)语音数据库中的8通道的阵列信号(10个说话人，每人13句话)中进行了测试。SNR大约在-6dB左右，采样率f_s为8KHz。采用oSNR，LLR和PESQ作为语音增强算法的客观评价性能指标。对比试验表明所提算法可以更精确地估计信号子空间的维度，在噪声消除和语音畸变之间进行更好的折中，从而得到了最高的PESQ得分。

表1在CMU阵列数据库(inputSNR≈-6dB)下的算法性能

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种F范数表征信号子空间维度的多通道语音增强方法，其特征在于，包括下列步骤：

步骤1：通过N个麦克风组成的麦克风阵列采集带噪声的多路语音信号y(t)，并计算该多路语音信号y(t)的带噪语音互相关矩阵R_yy，其中，t表示离散时间点；

步骤2：利用噪声估计算法估计加性噪声互相关矩阵

步骤3：利用所述带噪语音互相关矩阵R_yy和加性噪声互相关矩阵估计纯净语音互相关矩阵

步骤4：利用所述纯净语音互相关矩阵估计信号子空间的维度

步骤5：对进行广义特征值分解，并结合信号子空间的维度和拉格朗日乘子μ获得时域约束线性信号估计器；

步骤6：利用时域约束线性信号估计器对多路语音信号y(t)进行滤波获得增强后的语音

其中，步骤4中所述估计信号子空间的维度具体包括：

步骤41：计算纯净语音离散傅里叶变换系数幅值的超高斯分布模型：

$p\left({\mathrm{\&alpha;}}_x\right)=\frac{{\mathrm{\&eta;}}^{\mathrm{\&theta;}+1}}{\mathrm{\&Gamma;}(\mathrm{\&theta;}+1)}\frac{{\mathrm{\&alpha;}}_x^{\mathrm{\&theta;}}}{{\mathrm{\&sigma;}}_x^{\mathrm{\&theta;}+1}}\exp \{-\mathrm{\&eta;}\frac{{\mathrm{\&alpha;}}_x}{{\mathrm{\&sigma;}}_x}\}$

其中，θ和η是可调参数，α_x是纯净语音离散傅里叶变换系数的幅值，σ_x是纯净语音的标准差，Γ是Gamma函数，exp{·}是指数算子；

步骤42：利用所述纯净语音离散傅里叶变换系数幅值的超高斯分布模型获得纯净语音互相关矩阵的F范数：

${\left|\right|R_{\mathrm{xx}}\left|\right|}_f={\mathrm{\&beta;}}_2\sqrt{\underset{\mathrm{\&omega;}=0}{\overset{L-1}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_x^4\left(\mathrm{\&omega;}\right)}$

其中，∑{·}是求和算子，α_x(ω)是纯净语音在频点ω处的离散傅里叶变换系数的幅值， ${\mathrm{\&beta;}}_2=\sqrt{[N(N-1){\mathrm{\&beta;}}_1+N]L}$

步骤43：获取纯净语音互相关矩阵F范数的高斯分布模型：

$p_k\left(f_{\log }\right)=\frac{1}{\sqrt{\mathrm{\&pi;}}{\mathrm{\&sigma;}}_{\log }}\exp \{-\frac{{(f_{\log }-{\widehat{\mathrm{\&mu;}}}_{\log }\left(k\right))}^2}{{\mathrm{\&sigma;}}_{\log }^2}\}$

其中，f_log是纯净语音互相关矩阵R_xx的F范数的对数，即f_x＝||R_xx||_f，f_log＝log₁₀f_x，σ_log是f_log高斯分布模型的方差，π是圆周率，k表示当前帧号，是数学期望，如下计算：

${\widehat{\mathrm{\&mu;}}}_{\log }\left(k\right)=\mathrm{mean}\{\log {\left|\right|{\hat{R}}_{\mathrm{xx}}\left(k\right)\left|\right|}_f,\&CenterDot;\&CenterDot;\&CenterDot;,\log {\left|\right|{\hat{R}}_{\mathrm{xx}}(k-D+1)\left|\right|}_f\}$

其中，log{·}是取对数算子，||·||_f是取F范数算子，mean{·}是取均值算法，D是取均值的帧数；

步骤44：估计信号子空间维度

$\hat{Q}=\arg \underset{1\&le;q\&le;\mathrm{NL}}{\max }\{\frac{1}{2}\log \left(\underset{i=1}{\overset{q}{\mathrm{\&Sigma;}}}{\widehat{\mathrm{\&lambda;}}}_{x,i}^2\right)\&le;f_{\mathrm{\&delta;}}\}$

其中，表示对纯净语音互相关矩阵的估计进行特征值分解(EVD)得到的第i个特征值，q为在1到NL之间取值的整数，argmax{·}是取最大值坐标的算法，∑[·]是有限范围求和算子，f_δ是当前高斯分布模型下f_log的上限阈值，如下计算：

${\&Integral;}_{f_{\mathrm{\&delta;}}}^{+\&infin;}{p}_k\left(f_{\log }\right){\mathrm{df}}_{\log }=\mathrm{\&delta;}$

其中，δ是纯净语音互相关矩阵的估计的拒绝水平，p_k(f_log)是所述纯净语音互相关矩阵F范数的高斯分布模型；

所述步骤5具体包括以下内容：

步骤51、对进行广义特征值分解(GEVD)，得到特征向量矩阵W和特征值矩阵Λ：

步骤52、利用信噪比SNRγ计算拉格朗日乘子μ：

$\mathrm{\&mu;}=\left\{\begin{array}{ccc}5& \mathrm{for}& \mathrm{\&gamma;}\&le;-5\mathrm{dB}\\ 4.2-\frac{\mathrm{\&gamma;}}{6.25}& \mathrm{for}& -5\mathrm{dB}<\mathrm{\&gamma;}<20\mathrm{dB}\\ 1& \mathrm{for}& \mathrm{\&gamma;}\&GreaterEqual;20\mathrm{dB}\end{array}\right.$

步骤S53、生成TDC线性信号估计器：

$H={\mathrm{UW}}_1{G}_{\mathrm{TDC}}{W}_1^T$

其中，矩阵W₁是前个特征向量组成的列满秩矩阵，U＝[I_L×L 0_L×L … 0_L×L]是L×NL的矩阵，I_L×L是L阶单位阵，G_TDC为阶增益矩阵G_TDC，其每个元素如下计算：

$g_{\mathrm{TDC},i}=\frac{{\mathrm{\&lambda;}}_{\mathrm{xv},i}}{{\mathrm{\&lambda;}}_{\mathrm{xv},i}+\mathrm{\&mu;}},\mathrm{fori}=\mathrm{1,2},...,\hat{Q}$

其中，λ_xv，i是对进行广义特征值分解(GEVD)后得到的第i个特征值，μ是拉格朗日乘子。

2.如权利要求1所述的方法，其特征在于，所述麦克风阵列由N个等间距的麦克风组成，利用长为L的汉宁窗对语音信号进行截短，所述带噪语音互相关矩阵R_yy的阶数为NL×NL。

3.如权利要求1所述的方法，其特征在于，所述步骤2中加性噪声互相关矩阵的估计通过下式得到：

${\hat{R}}_{\mathrm{vv},k}=\left\{\begin{array}{cc}{\hat{R}}_{\mathrm{vv},k-1}& {\mathrm{if\&gamma;}}_k>{\mathrm{\&gamma;}}_{\mathrm{th}}\\ 0.98\&CenterDot;{\hat{R}}_{\mathrm{vv},k-1}+R_{\mathrm{yy}}& \mathrm{else}\end{array}\right.$

其中，γ_th为SNR的下限阈值，γ_k为第k帧的SNR；所述是前一帧加性噪声互相关矩阵的估计，R_yy是当前帧带噪语音互相关矩阵，k是帧号。

4.如权利要求1所述的方法，其特征在于，所述步骤41还包括：

利用纯净语音数据训练得到纯净语音离散傅里叶变换系数幅值的统计柱状图p_o；

通过将所述纯净语音离散傅里叶变换系数幅值的超高斯分布模型和纯净语音离散傅里叶变换系数幅值的统计柱状图p_o之间的KLD距离最小化，得到最佳参数组合(θ，η)，其中，KLD距离如下计算：

$J({\mathrm{\&alpha;}}_x:o)=\underset{0}{\overset{M}{\&Integral;}}(p_o\left(m\right)-p_{{\mathrm{\&alpha;}}_x}\left(m\right))\log \left(\frac{p_o\left(m\right)}{p_{{\mathrm{\&alpha;}}_x}\left(m\right)}\right)\mathrm{dm}$

其中，是所述纯净语音离散傅里叶变换DFT系数幅值的超高斯分布模型，p_o是纯净语音离散傅里叶变换系数幅值的统计柱状图，log{·}是求自然对数算子，m是纯净语音DFT系数幅值的统计柱状图中纯净语音DFT系数幅值的可能值，其取值范围从0到M。

5.如权利要求1所述的方法，其特征在于，所述当前帧语音互相关矩阵的拒绝水平δ如下计算：

$\mathrm{\&delta;}=\left\{\begin{array}{ccc}{\mathrm{\&delta;}}_{\max }& \mathrm{for}& \mathrm{\&gamma;}<{\mathrm{\&gamma;}}_{\min }\\ \frac{{\mathrm{\&delta;}}_{\max }-{\mathrm{\&delta;}}_{\min }}{{10}^{\frac{{\mathrm{\&gamma;}}_{\max }}{10}}-{10}^{\frac{{\mathrm{\&gamma;}}_{\min }}{10}}}\&CenterDot;({10}^{\frac{{\mathrm{\&gamma;}}_{\max }}{10}}-{10}^{\frac{\mathrm{\&gamma;}}{10}})& \mathrm{for}& {\mathrm{\&gamma;}}_{\min }\&le;\mathrm{\&gamma;}\&le;{\mathrm{\&gamma;}}_{\max }\\ {\mathrm{\&delta;}}_{\min }& \mathrm{for}& \mathrm{\&gamma;}>{\mathrm{\&gamma;}}_{\max }\end{array}\right.$

其中，γ_max＝10dB，γ_min＝-10dB，δ_max＝0.5，δ_min＝0.05，γ为当前帧的后验信噪比SNR。

6.如权利要求3-5中任一项所述的方法，其特征在于，所述当前帧的后验信噪比SNRγ如下计算：

$\mathrm{\&gamma;}=10\&CenterDot;\log \frac{R_{\mathrm{yy}}\left(\mathrm{1,1}\right)}{{\hat{R}}_{\mathrm{vv},k-1}\left(\mathrm{1,1}\right)}$

其中，log{·}是求自然对数算子，是前一帧加性噪声互相关矩阵中第1行第1列元素的估计。

7.如权利要求1所述的方法，其特征在于，所述步骤6具体为：利用所述TDC线性信号估计器H对所述多路语音信号y(t)进行滤波，并对得到的滤波后的语音信号进行帧间重叠获得增强后语音，其中滤波如下公式所示：

$\hat{x}\left(t\right)=\mathrm{Hy}\left(t\right).$

8.如权利要求1所述的方法，其特征在于，所述步骤3中纯净语音互相关矩阵如下估计：

${\hat{R}}_{\mathrm{xx}}=R_{\mathrm{yy}}-{\hat{R}}_{\mathrm{vv},k}.$