多通道群稀疏线性预测时延估计方法

摘要

本发明公开了一种多通道群稀疏线性预测时延估计方法，采用各个可能的时延值计算多通道互相关系数，将多通道互相关系数的平方最大值所对应的时延作为时延估计值；多通道互相关系数的计算方法为：对每只传声器拾取的声信号分别截取一个长度为L的信号帧，对这些信号帧进行时移后叠放成信号向量以及信号矩阵，根据F/l1,2范数优化准则建立时延估计模型，求解群稀疏预测系数矩阵，然后计算预测误差相关矩阵，根据预测误差相关矩阵计算得到多通道互相关系数。本发明利用预测系数矩阵列向量的群稀疏特性构建一个F/l1,2范数优化准则，统一了多通道互相关系数方法和多通道空时预测方法，联合应用空间和时间线性预测的白化能力，以提高时延估计性能。

说明书

技术领域

本发明属于声源定位技术领域，更为具体地讲，涉及一种多通道群稀疏线 性预测时延估计方法。

背景技术

时延估计是利用传声器阵列拾取的声信号估计声源到各阵元间的到达时间 差，在雷达、声呐以及免提语音通信系统的声源定位和跟踪中扮演着极其重要 的角色。常见的时延估计方法有广义互相关方法、盲通道辨识方法、信息论方 法、基于语音信号的某些特征进行时延估计等。然而在室内声环境下，由于噪 声和混响的不利影响，基于传声器阵列的时延估计是一项具有挑战性的工作。

为了提高时延估计对噪声和混响的鲁棒性，一种多通道互相关系数方法被 提出，其具体方法可参考文献“J.Chen,J.Benesty,andY.Huang,“Robusttime delayestimationexploitingredundancyamongmultiplemicrophones,”IEEETrans. SpeechAudioProcess.,vol.11,no.6,pp.549-557,Nov.2003.”和“J.Benesty,J. Chen,andY.Huang,“Time-delayestimationvialinearinterpolationand cross-correlation,”IEEETrans.SpeechAudioProcess.,vol.12,no.5,pp.509-519, Sep.2004.”。相对于传统的双通道方法，多通道互相关系数方法利用多只传声器 间的空间冗余信息抑制噪声和混响的影响，极大地提高了时延估计对噪声的鲁 棒性。然而，这种多通道方法对混响仍然敏感。多通道空时预测方法对多通道 互相关系数方法进行了推广，其具体方法可参考文献“H.He,L.Wu,J.Lu,X.Qiu, andJ.Chen,“Timedifferenceofarrivalestimationexploitingmultichannel spatio-temporalprediction,”IEEETrans.AudioSpeechLang.Process.,vol.21,pp. 463-475,Mar.2013.”。这种方法以目前最优的方式利用空间和时间信息预白化传 声器信号，由此获得时延估计对混响的鲁棒性。然而，这种方法在低信噪比条 件下性能降低。

线性预测广泛应用于语音处理。传统的线性预测器由长时预测器和短时预 测器级联组成。对于语音信号，这种结构的预测系数向量具有高度的稀疏性。 该特性目前已成功应用于语音编码领域。然而，当语音信号受到噪声污染时， 这种稀疏性降低甚至消失，线性预测器的性能急剧下降。在最近的研究中“H. He,T.Yang,andJ.Chen,“Ontimedelayestimationfromasparselinearprediction perspective,”J.Acoust.Soc.Amer.,vol.137,no.2,pp.1044-1047,Feb.2015.”一种 稀疏线性预测时延估计方法被提出。这种方法在最小二乘基础上引入一个稀疏 正则化项，构成一种l₂/l₁范数优化方法，提高了时延估计对噪声的鲁棒性。但该 方法是在双通道情形下对两个通道的数据分别进行预白化，也即仅仅对数据进 行时间预白化，没有综合利用多通道空间和时间预测的白化能力，从而降低该 方法的估计性能。

发明内容

本发明的目的在于克服现有技术的不足，提供一种多通道群稀疏线性预测 时延估计方法，利用预测系数矩阵列向量的群稀疏特性构建一个F/l_1,2范数优化 准则，统一了多通道互相关系数方法和多通道空时预测方法，联合应用空间和 时间线性预测的白化能力，以提高时延估计性能。

为实现上述发明目的，本发明多通道群稀疏线性预测时延估计方法包括以 下步骤：

S1：M只传声器分别对声源信号进行持续采集，第m只传声器采集的时域 信号记为X_m＝[x_m(0),x_m(1),…,x_m(L-1)]，其中m＝1,2,…,M，x_m(n)表示第m只传声 器在时刻n的采集样本，n＝0,1,…,L-1，L表示每只传声器采集的样本数量；

S2：令时移序号d＝1，声源信号到达第1只和第2只传声器间的初始时移 p₁＝-p_max，p_max表示时延的最大可能值；

S3：分别将第m只传声器采集的信号X_m按照时延f_m(p_d)进行时移，f_m(p_d)表 示声源信号到达第1只和第m只传声器间的相对时延，该时延是关于时移p_d的 函数；时移后的信号中每个样本记为x_m(n,p_d)，将M只传声器在时刻n的时移样 本x_m(n,p_d)进行叠放，得到信号向量x(n,p_d)；

S4：求解以下公式，得到预测系数矩阵A(p_d)：

$\underset{A\left(p_d\right)}{min}\left\{\frac{1}{2}\right||X(0,p_d)-Y(-1,p_d)A\left(p_d\right)|{|}_F^2+\lambda \left|\right|A\left(p_d\right)\left|{|}_{l_{1,2}}\right\}$

其中，||·||_F表示矩阵的F范数，表示矩阵的l₂范数，λ是正则化参数， 取值范围为λ＞0；

X(0,p_d)＝[x(0,p_d)x(1,p_d)…x(K+L-1,p_d)]^T

Y(-1,p_d)＝[y(-1,p_d)y(0,p_d)…y(K+L-2,p_d)]^T

其中：

x(w,p_d)＝[x₁(w,p_d),x₂(w,p_d),…x_M(w,p_d)]^T

$y(w-1,p_d)={\left(\begin{array}{cccc}{x}_1^T(w-1,p_d)& x_2^T(w-1,p_d)& ...& x_M^T(w-1,p_d)\end{array}\right)}^T$

x_m(w-1,p_d)＝[x_m(w-1,p_d),x_m(w-2,p_d),…x_m(w-K,p_d)]^T

w＝0,1,…,K+L-1；K表示预测器阶数，其取值范围为K＜L；对于x_m(q,p_d)， 如果q＜0或q＞L-1，则令x_m(q,p_d)＝0；

S5：计算预测误差矩阵E(0,p_d)：

E(0,p_d)＝X(0,p_d)-Y(-1,p_d)A(p_d)

S6：计算预测误差相关矩阵R(p_d)：

$R\left(p_d\right)=\frac{1}{K+L}{E}^T(0,p_d)E(0,p_d)$

S7：计算多通道互相关系数ρ(p_d)：

${\rho }^2\left(p_d\right)=1-\frac{\det \left(R\right(p_d\left)\right)}{{\Pi }_{m=1}^M{r}_{m,m}\left(p_d\right)}$

其中det(·)表示方阵的行列式，r_m,m(p_d)是矩阵R(p_d)的第m个对角元素；

S8：如果p_d＜p_max，令p_d＝p_d+1，返回步骤S3，否则根据以下公式求得时 延估计值

$\hat{\tau }=\arg \underset{p_d}{max}{\rho }^2\left(p_d\right).$

本发明多通道群稀疏线性预测时延估计方法，采用各个可能的时延值计算 多通道互相关系数，将多通道互相关系数的平方最大值所对应的时延作为时延 估计值；多通道互相关系数的计算方法为：对每只传声器拾取的声信号分别截 取一个长度为L的信号帧，对这些信号帧进行时移后叠放成信号向量以及信号矩 阵，根据F/l_1,2范数优化准则建立时延估计模型，求解群稀疏预测系数矩阵，然 后计算预测误差相关矩阵，根据预测误差相关矩阵计算得到多通道互相关系数。

本发明将基于稀疏线性预测的l₂/l₁范数优化方法推广到多通道，提出了一种 多通道群稀疏线性预测时延估计方法。本发明利用预测系数矩阵列向量的群稀 疏特性构建一个F/l_1,2范数优化准则，从时延估计性能角度统一了多通道互相关 系数方法和多通道空时预测方法。本发明还可以通过调节正则化参数来构建一 组在传声器信号预白化和非预白化间可作折中处理的时延估计器，以适应不同 需求。

附图说明

图1是基于最小二乘的多通道空时预测器对纯净语音预测时预测系数矩阵 的一个列向量；

图2是基于最小二乘的多通道空时线性预测器对含噪语音预测时预测系数 矩阵的一个列向量；

图3是采用本发明对含噪语音预测时预测系数矩阵的一个列向量；

图4是本发明多通道群稀疏线性预测时延估计方法的流程图；

图5是基于增广拉格朗日乘子交替方向法求解预测系数矩阵的流程图；

图6是噪声和混响环境下各方法时延估计性能随传声器个数变化的曲线图；

图7是噪声环境下各方法时延估计性能随混响时间变化的曲线图；

图8是轻度混响环境下各方法时延估计性能随信噪比SNR变化的曲线图；

图9是中度混响环境下各方法时延估计性能随信噪比SNR变化的曲线图；

图10是轻度混响环境下本发明时延估计性能随参数δ变化的曲线图；

图11是中度混响环境下本发明时延估计性能随参数δ变化的曲线图。

具体实施方式

下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员 更好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和 设计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。

实施例

为了更好地说明本发明的技术内容，首先对本发明所使用的信号模型及方 法推导进行说明。

假设一宽带声源辐射平面波信号，M只传声器构成的阵列用于接收声信号。 如果选择第1只传声器作为参考点，则第m只传声器在时刻n拾取的信号样本 x_m(n)可表示为：

x_m(n)＝α_ms[n-t-f_m(τ)]+w_m(n)(1)

其中，m＝1,2,…,M，α_m是声信号传输效应导致的衰减因子，s(n)是未知的 零均值因果宽带声源信号，t是从声源到第1只传声器间的传输时间，w_m(n)是 第m只传声器上的加性噪声，假设该噪声与声源信号以及其它传声器的噪声不 相关，τ是声源信号到达第1只和第2只传声器间的时延，f_m(τ)是声源信号到 达第1只和第m只传声器间的相对时延，它是时延τ的函数。本实施例为了易于 推导与实验，假设使用等间距的线性阵列，则在远场条件下f_m(τ)＝(m-1)τ。

在上述信号模型基础上，时延估计的目的是给定M只传声器采集的样本信 号来估计时延τ。对于一个假定的时延p，利用时移信号x_m(n+f_m(p))对齐传声 器信号。当p＝τ时，M只传声器信号对齐，此时M只传声器拾取的信号间相似 性最大。为简化符号，将x_m(n+f_m(p))记为x_m(n,p)。

将M只传声器在时刻n拾取的信号样本表示成一个向量：

x(n,p)＝[x₁(n,p),x₂(n,p),…x_M(n,p)]^T(2)

其中，x_m(n,p)表示第m只传声器在时刻n的信号样本的时移信号，(·)^T表示 向量或矩阵的转置。

为了便于推导，对第m只传声器的信号在时刻n-1定义另一个向量：

x_m(n-1,p)＝[x_m(n-1,p),x_m(n-2,p),…x_m(n-K,p)]^T(3)

其中K表示预测器阶数，其取值范围为K＜L。

现在利用M个通道过去的信号向量x₁(n-1,p)、x₂(n-1,p)、…、x_M(n-1,p)来 预测x(n,p)，即预测值表示为：

$\hat{x}(n,p)=A_1\left(p\right)x_1(n-1,p)+A_2\left(p\right)x_2(n-1,p)+...+A_M\left(p\right)x_M(n-1,p)---\left(4\right)$

其中是多通道前向预测器的系数矩阵。

于是，可得预测误差向量e(n,p)的表达式为：

$e(n,p)=x(n,p)-\hat{x}(n,p)=x(n,p)-A^T\left(p\right)y(n-1,p)---\left(5\right)$

e(n,p)可以记为：

e(n,p)＝[e₁(n,p),e₂(n,p),…e_M(n,p)]^T(6)

A(p)是KM×M的多通道前向预测系数矩阵，可以表示为：

A(p)＝[A₁(p)A₂(p)…A_M(p)]^T(7)

y(n-1,p)是M只传声器在时刻n-1拾取的信号向量，可以表示为：

$y(n-1,p)={\left(\begin{array}{cccc}{x}_1^T(n-1,p)& x_2^T(n-1,p)& ...& x_M^T(n-1,p)\end{array}\right)}^T---\left(8\right)$

根据以上推导过程可以看出，公式(5)是采用M个通道时刻n-K到时刻 n-1的样本预测M个通道时刻n的样本，以期获得对应的预测系数矩阵。在实际 应用中，为了利用预测稀疏矩阵的群稀疏特性以提高时延估计的准确度，需要 将(5)写成矩阵方程。因此利用M只传声器采集的长度分别为L的数据帧，根 据公式(5)扩展得到：

E(n,p)＝X(n,p)-Y(n-1,p)A(p)(9)

其中：

E(n,p)＝[e(n,p)e(n+1,p)…e(n+K+L-1,p)]^T(10)

X(n,p)＝[x(n,p)x(n+1,p)…x(n+K+L-1,p)]^T(11)

Y(n-1,p)＝[y(n-1,p)y(n,p)…y(n+K+L-2,p)]^T(12)

其中，X(n,p_d)是(K+L)×M的矩阵，Y(n-1,p_d)是(K+L)×KM的矩阵。本 发明中，将传声器采集的长度为L的样本中第一个样本序号作为时刻n来构建 X(n,p_d)和Y(n-1,p_d)，即n＝0。由于数据帧的长度为L，那么在构建X(n,p)和 Y(n-1,p)时，对于时刻q的信号样本x_m(q,p)，如果q＜0或q＞L-1，则令 x_m(q,p)＝0。

传统线性预测器由长时预测器和短时预测器级联组成，由此预测器系数向 量具有高度的稀疏性。因此在多通道线性预测中，预测器系数矩阵也具有稀疏 性。以预测器阶数为80、传声器采集的语音帧长度为1024个采样点、4只传声 器用于时延估计为例，对多通道预测器的系数矩阵进行分析。图1是基于最小 二乘的多通道空时预测器对纯净语音预测时预测系数矩阵的一个列向量。从图1 可以看出，基于最小二乘的多通道空时线性预测系数矩阵的列向量是稀疏的， 而且预测系数矩阵不同列向量的大能量样本基本上处于相同的位置，这说明所 有列向量具有类似的群稀疏特性。图2是基于最小二乘的多通道空时线性预测 器对含噪语音预测时预测系数矩阵的一个列向量。图2中语音信号的信噪比 SNR＝5dB。从图2可以看出，当存在噪声时，预测器系数矩阵的稀疏特性降低。 既然对于纯净语音信号，预测器系数矩阵是群稀疏的，则可利用这一特性来提 高噪声环境下线性预测器的鲁棒性。为此，本发明在最小二乘基础上引入一个 预测系数矩阵的群稀疏正则化约束，强制预测系数矩阵的列向量具有相同的稀 疏结构。于是，提出如下F/l_1,2范数优化准则来预处理传声器信号：

$\underset{A\left(p_d\right)}{min}\left\{\frac{1}{2}\right||X(n,p_d)-Y(n-1,p_d)A\left(p_d\right)|{|}_F^2+\lambda \left|\right|A\left(p_d\right)\left|{|}_{l_{1,2}}\right\}---\left(13\right)$

其中，||·||_F表示矩阵的F范数，λ是正则化参数，取值范围为λ＞0。可表示为：

$\left|\right|A\left(p\right)|{|}_{l_{1,2}}=\underset{i=1}{\overset{MK}{\Sigma }}|\left|A^i\left(p\right)\right|{|}_{l_2}---\left(14\right)$

其中，表示矩阵或向量的l₂范数，Aⁱ(p)表示预测系数矩阵A(p)的第i行。

本发明中，令n＝0，求解公式(13)即可得到预测系数矩阵A(p)，再根据 预测系数矩阵A(p)即可计算出时延。

图3是采用本发明对含噪语音预测时预测系数矩阵的一个列向量。图3中 语音信号的信噪比SNR＝5dB，从图3可以看出，在噪声环境下通过本发明所提 出的如下F/l_1,2范数优化准则可以获得稀疏的预测系数矩阵。

基于以上推导，可以得到本发明多通道群稀疏线性预测时延估计方法。图4 是本发明多通道群稀疏线性预测时延估计方法的流程图。如图4所示，本发明 多通道群稀疏线性预测时延估计方法的具体步骤包括：

S401：信号采集：

M只传声器分别对声信号进行持续采集，第m只传声器采集的时域信号记 为X_m＝[x_m(0),x_m(1),...,x_m(L-1)]，其中m＝1,2,...,M，x_m(n)表示第m只传声器在时 刻n的采集样本，n＝0,1,...,L-1，L表示每只传声器采集的样本数量。

S402：初始化时移序号：

令时移序号d＝1，声源信号到达第1只和第2只传声器间的初始时移 p₁＝-p_max，p_max表示时延的最大可能值。

S403：信号时移：

分别将第m只传声器采集的信号X_m按照时延f_m(p_d)进行时移，f_m(p_d)表示 声源信号到达第1只和第m只传声器间的相对时延，该时延是关于时移p_d的函 数。时移后的信号样本即为x_m(t+f_m(p_d))，为简化表达，记第m只传声器在时刻 n的时移信号为x_m(n,p_d)。将M只传声器在时刻n的时移信号x_m(n,p_d)进行叠放， 得到信号向量x(n,p_d)。

S404：求解预测系数矩阵：

求解以下公式，得到预测系数矩阵A(p_d)：

$\underset{A\left(p_d\right)}{min}\left\{\frac{1}{2}\right||X(0,p_d)-Y(-1,p_d)A\left(p_d\right)|{|}_F^2+\lambda \left|\right|A\left(p_d\right)\left|{|}_{l_{1,2}}\right\}---\left(15\right)$

其中，A(p_d)是KM×M的矩阵，||·||_F表示矩阵的F范数，表示矩阵或向 量的l₂范数，λ是正则化参数，取值范围为λ＞0，

X(0,p_d)＝[x(0,p_d)x(1,p_d)...x(K+L-1,p_d)]^T(16)

Y(-1,p_d)＝[y(-1,p_d)y(0,p_d)...y(K+L-2,p_d)]^T(17)

其中：

x(w,p_d)＝[x₁(w,p_d),x₂(w,p_d),…x_M(w,p_d)]^T(18)

$y(w-1,p_d)={\left(\begin{array}{cccc}{x}_1^T(w-1,p_d)& x_2^T(w-1,p_d)& ...& x_M^T(w-1,p_d)\end{array}\right)}^T---\left(19\right)$

x_m(w-1,p_d)＝[x_m(w-1,p_d),x_m(w-2,p_d),…x_m(w-K,p_d)]^T(20)

w＝0,1,…,K+L-1；K表示预测器阶数，其取值范围为K＜L。对于构建矩 阵X(0,p_d)和Y(-1,p_d)的时刻q的信号样本x_m(q,p_d)，如果q＜0或q＞L-1，则令 x_m(q,p_d)＝0。显然根据公式(16)至公式(20)可知，q的取值范围为[-K,K+L-1]。

可表示为：

$\left|\right|A\left(p_d\right)|{|}_{l_{1,2}}=\underset{i=1}{\overset{MK}{\Sigma }}|\left|A^i\left(p_d\right)\right|{|}_{l_2}---\left(21\right)$

其中，表示矩阵或向量的l₂范数，Aⁱ(p_d)表示预测系数矩阵A(p_d)的第i 行。

S405：计算预测误差矩阵：

将步骤S404得到的预测系数矩阵，代入(9)式，即可得到预测误差矩阵 E(0,p_d)，即：

E(0,p_d)＝X(0,p_d)-Y(-1,p_d)A(p_d)(22)

S406：计算预测误差相关矩阵：

根据预测误差矩阵E(0,p_d)计算预测误差相关矩阵R(p_d)：

$R\left(p_d\right)=\frac{1}{K+L}{E}^T(0,p_d)E(0,p_d)---\left(23\right)$

S407：计算多通道互相关系数：

利用预测误差相关矩阵R(p_d)计算多通道互相关系数ρ(p_d)：

${\rho }^2\left(p_d\right)=1-\frac{\det \left(R\right(p_d\left)\right)}{{\Pi }_{m=1}^M{r}_{m,m}\left(p_d\right)}---\left(24\right)$

其中det(·)表示方阵的行列式，r_m,m(p_d)是矩阵R(p_d)的第m个对角元素。

S408：判断是否p_d＜p_max，如果是，进入步骤S409，否则进入步骤S410。

S409：令p_d＝p_d+1，返回步骤S403。

S410：计算得到时延估计值：

根据以下公式求得时延估计值

$\hat{\tau }=\arg \underset{p_d}{max}{\rho }^2\left(p_d\right)---\left(25\right)$

公式(25)的含义是，在[-p_max,p_max]范围内的各个时移p_d中，当其对应的 ρ²(p_d)最大，即将该p_d作为时延估计值。

显然，在步骤S404中，(15)式是一个凸优化问题，可用多种方法求解， 例如线性规划、内点法、原始-对偶内点法等。增广拉格朗日乘子交替方向法能 有效利用多个变量可分离(可解耦)的特性，因此本实施例中采用该方法求解 这一问题。引入一个辅助矩阵Z(p_d)，(13)式可等价地表示为：

$\underset{A\left(p_d\right),Z\left(p_d\right)}{\min }\left\{\frac{1}{2}\left|\right|X\left(0,p_d\right)-Y\left(-1,p_d\right)A\left(p_d\right)|{|}_F^2+\lambda |\left|Z\left(p_d\right)\right|{|}_{l_{1,2}}:A\left(p_d\right)-Z\left(p_d\right)=0\right\}---\left(26\right)$

其增广拉格朗日子问题可表示为

$\begin{array}{c}\underset{A\left(p_d\right),Z\left(p_d\right)}{min}\left\{\frac{1}{2}\right||X(0,p_d)-Y(-1,p_d)A\left(p_d\right)|{|}_F^2+\lambda \left|\right|Z\left(p_d\right)|{|}_{l_{1,2}}\\ +<\theta \left(p_d\right),A\left(p_d\right)-Z\left(p_d\right)>+\frac{\beta }{2}\left|\right|A\left(p_d\right)-Z\left(p_d\right)\left|{|}_F^2\right\}\end{array}---\left(27\right)$

式中是线性约束拉格朗日乘子矩阵；β是惩罚参数，其取值范 围为β＞0；〈θ(p_d),A(p_d)-Z(p_d)〉是矩阵内积，其表达式为：

〈θ(p_d),A(p_d)-Z(p_d)〉＝tr{θ^T(p_d)[A(p_d)-Z(p_d)]}(28)

其中，tr(·)表示矩阵的迹；

(27)式中大括号内的第四项、即增广项，用于确保目标函数是严格凸的。 (27)式的迭代求解过程中，给定Z_k(p_d)，θ_k(p_d)，可利用(27)式交替保持一 个未知矩阵固定对另一个矩阵最小化求得A_k+1(p_d)，Z_k+1(p_d)，θ_k+1(p_d)。

首先，当Z(p_d)＝Z_k(p_d)，θ(p_d)＝θ_k(p_d)固定，(27)式对A(p_d)最小化等价 于：

$\underset{A\left(p_d\right)}{min}\left\{\frac{1}{2}\right||X(0,p_d)-Y(-1,p_d)A\left(p_d\right)|{|}_F^2+\frac{\beta }{2}\left|\right|A\left(p_d\right)-Z_k\left(p_d\right)+{\theta }_k\left(p_d\right)/\beta \left|{|}_F^2\right\}---\left(29\right)$

其解为：

$\begin{array}{c}{A}_{k+1}\left(p_d\right)={[Y^T(-1,p_d)Y(-1,p_d)+\beta I]}^{-1}\\ \times \left[Y^T\left(-1,p_d\right)X\left(0,p_d\right)+{\mathrm{\beta Z}}_k\left(p_d\right)-{\theta }_k\left(p_d\right)\right]\end{array}---\left(30\right)$

式中I表示KM×KM的单位矩阵。

当A(p_d)＝A_k+1(p_d)，θ(p_d)＝θ_k(p_d)固定，(27)式对Z(p_d)最小化等价于：

$\underset{Z\left(p_d\right)}{min}\left\{\lambda \right|\left|Z\left(p_d\right)\right|{|}_{l_{1,2}}+\frac{\beta }{2}\left|\right|A_{k+1}\left(p_d\right)-Z\left(p_d\right)+{\theta }_k\left(p_d\right)/\beta \left|{|}_F^2\right\}---\left(31\right)$

Z(p_d)的第i行Zⁱ(p_d)的解为：

$Z_{k+1}^i\left(p_d\right)=soft\left(A_{k+1}^i\left(p_d\right)+{\theta }_k^i\left(p_d\right)/\beta ,\lambda /\beta \right)---\left(32\right)$

其中表示辅助矩阵Z_k+1(p_d)的第i行，i＝1,2,…,MK；

soft函数定义为：

最后，线性约束拉格朗日乘子矩阵θ(p_d)更新如下：

θ_k+1(p_d)＝θ_k(p_d)+β(A_k+1(p_d)-Z_k+1(p_d))(34)

因此通过迭代计算(30)、(32)和(34)式获得(13)式的解。图5是基 于增广拉格朗日乘子交替方向法求解预测系数矩阵的流程图。如图5所示，本 实施例中求解预测系数矩阵A(p_d)的具体过程为：

S501：初始化参数：

令迭代次数k＝1，初始化KM×M的辅助矩阵Z₁(p_d)和拉格朗日乘子矩阵 θ₁(p_d)。

S502：计算预测系数矩阵：

$\begin{array}{c}{A}_{k+1}\left(p_d\right)={\left[Y^T\left(-1,p_d\right)Y\left(-1,p_d\right)+\beta I\right]}^{-1}\\ \times \left[Y^T\left(-1,p_d\right)X\left(0,p_d\right)+{\mathrm{\beta Z}}_k\left(p_d\right)-{\theta }_k\left(p_d\right)\right]\end{array}---\left(35\right)$

S503：更新辅助矩阵：

$Z_{k+1}^i\left(p_d\right)=soft\left(A_{k+1}^i\right(p_d)+{\theta }_k^i(p_d)/\beta ,\lambda /\beta ),i=1,2,...,MK---\left(36\right)$

S504：更新拉格朗日乘子矩阵：

θ_k+1(p_d)＝θ_k(p_d)+β(A_k+1(p_d)-Z_k+1(p_d))(37)

S505：判断是否k＜Q，Q表示最大迭代次数，如果是，进入步骤S506，否 则进入步骤S507。

S506：令k＝k+1，返回步骤S502。

S507：得到预测系数矩阵：

将当前的预测系数矩阵A_k+1(p_d)作为最终的预测系数矩阵A(p_d)，即 A(p_d)＝A_k+1(p_d)。

此外，从(15)式可看出，参数λ在控制预测器系数矩阵的稀疏程度方面扮 演着重要角色。该参数主要受传声器信号影响，也即受矩阵X(0,p)和Y(-1,p)的 影响，因此可以通过如下方式确定参数λ：

$\lambda =\delta \left|\right|X^T(0,p)Y(-1,p)|{|}_{l_{\infty }}---\left(38\right)$

其中表示求取矩阵的无穷范数，δ是一个正常数，根据需要进行设置。

在实际应用中，通常需要持续实时地对时延进行估计，因此采用本发明对 M只传声器分别采集的长度为L的数据帧进行时延估计后，M只传声器需要分 别采集长度为L的新数据帧，再利用本发明进行下一次时延估计。每次进行时延 估计时，都以第1个样本(即时刻0的样本)作为基准来构建矩阵，以求解预 测系数矩阵。

为了更好地说明本发明的技术效果，采用本发明(记为MCSTGSP)对一个 具体的实施例进行实验验证，并以多通道互相关系数(MCCC)方法、具有预白 化的多通道互相关系数(预白化MCCC)方法、多通道空时预测(MCSTP)方 法作为对照算法对实验结果进行对照。

实验场地为一个7m×6m×3m的房间。使用6只等间距全指向传声器构成的 直线阵列拾取传声器信号，阵元间距为0.1m。假设房间地面西南角为坐标原点， 房间内任意位置的坐标表示为(x,y,z)。阵列的第1只和第6只传声器分别放置 于(3.25,3.00,1.40)和(3.75,3.00,1.40)。声源位于(2.49,1.27,1.40)。

声源信号是一段预先录制的采样率为16kHz的纯净语音信号，其长度约为 1分钟。声源到6只传声器间的脉冲响应通过Image模型产生，Image模型可参 考文献“J.B.AllenandD.A.Berkley,“Imagemethodforefficientlysimulating small-roomacoustics,”J.Acoust.Soc.Amer.,vol.65,pp.943-950,Apr.1979.”脉冲 响应长度为2048个采样点。声源信号与对应通道脉冲响应进行卷积获得传声器 信号，在其中加入零均值高斯白噪声以控制信噪比SNR。

在实验中，将传声器信号分割成长度为64ms互不重叠的信号帧。对每一帧 信号用Hamming窗函数加窗后用于时延估计。采用畸变估计概率和非畸变估计 的根均方误差这两种准则评价时延估计方法的性能(有关这两种准则的定义和 分类方法参见文献“H.He,L.Wu,J.Lu,X.Qiu,andJ.Chen,“Timedifferenceof arrivalestimationexploitingmultichannelspatio-temporalprediction,”IEEETrans. AudioSpeechLang.Process.,vol.21,pp.463-475,Mar.2013.”，“J.P.Ianniello, “Timedelayestimationviacross-correlationinthepresenceoflargeestimation errors,”IEEETrans.Acoust.,Speech,SignalProcess.,vol.ASSP-30,pp.998-1003, Dec.1982.”和“B.Champagne,S.Bedard,andA.Stephenne,“Performanceof time-delayestimationinthepresenceofroomreverberation,”IEEETrans.Speech AudioProcess.,vol.4,pp.148-152,Mar.1996.”)。用于性能统计的语音帧总数为 936帧(帧长为1024个采样点)。声源到前两只传声器间的真实时延为2.0个采 样间隔。

图6是噪声和混响环境下各方法时延估计性能随传声器个数变化的曲线图。 图6展示了噪声(SNR＝10dB)和混响(混响时长T₆₀＝300ms)环境下，本发明 方法与对比方法的时延估计性能随传声器个数变化的关系。其中对于本发明方 法，δ＝0.001。可以看出，随着传声器个数增加，四种时延估计方法的畸变估计 概率和非畸变估计的根均方误差基本上都降低，说明适当增加传声器个数能有 效改善时延估计的鲁棒性。当使用两只传声器进行时延估计时，所有具有预白 化能力的时延估计方法都未获得明显的优势，而原始的MCCC方法具有一定的 鲁棒性。当使用多只传声器时，尽管MCCC方法获得较小的畸变估计概率，但 对应的根均方误差较大。尽管预白化MCCC和多通道空时预测MCSTP方法获 得适当的根均方误差，然而其畸变估计概率较大。相比之下，本发明提出的 MCSTGSP方法在畸变估计概率和非畸变估计的根均方误差间获得了良好的折 中，这表明在噪声和混响环境下本发明用于时延估计是有效的。

图7是在噪声环境下各方法时延估计性能随混响时间变化的曲线图。图7 中传声器信号的信噪比SNR＝10dB。其中对于本发明方法，δ＝0.001。如图7 所示，当混响时长T₆₀＝0ms时，MCCC获得了最好的性能，表明该方法对噪声 具有鲁棒性。随着T₆₀增加，尽管MCCC方法的畸变估计概率较小，但根均方误 差较大，表明MCCC对混响鲁棒性较差。尽管MCSTP方法对混响具有良好的 鲁棒性，但对噪声不具有鲁棒性。对于预白化MCCC、MCSTP和MCSTGSP这 三种具有预白化能力的时延估计方法而言，MCSTGSP方法获得了最好的性能， 表明提出的多通道群稀疏预测时延估计方法在不同混响环境是有效的。

图8是轻度混响环境下各方法时延估计性能随信噪比SNR变化的曲线图。 图9是中度混响环境下各方法时延估计性能随信噪比SNR变化的曲线图。图8 中混响时长T₆₀＝120ms，图9中混响时长T₆₀＝300ms。其中对于本发明方法， δ＝0.001。根据图8和图9可以看出原始的MCCC算法对噪声具有最好的鲁棒 性，尤其在低信噪比条件下，然而在高信噪比条件下MCCC算法对混响最敏感。 相对于MCCC，预白化MCCC对混响获得了更好的鲁棒性。MCSTP和预白化 MCCC在低信噪比条件下获得相当的性能，然而，由于其最优的预白化能力， MCSTP对混响更具有鲁棒性。尽管这两种具有预白化能力的时延估计方法的性 能在混响条件下获得了较大的提高，但在噪声影响下其性能降低。本发明 MCSTGSP在MCCC和MCSTP间获得了良好的折中，尤其是MCSTGSP增强 了对噪声的免疫力，这证实了利用预测系数矩阵的群稀疏特性可提高多通道空 时预测对噪声的鲁棒性。

图10是轻度混响环境下本发明时延估计性能随参数δ变化的曲线图。图11 是中度混响环境下本发明时延估计性能随参数δ变化的曲线图。图10中混响时 长T₆₀＝120ms，图11中混响时长T₆₀＝300ms。从图10和图11可以看出，一方面 正则化参数δ越小，MCSTGSP越接近MCSTP方法，极限情况是δ＝0，即 MCSTGSP退化成MCSTP。另一方面，随着δ增加，预测系数矩阵越来越稀疏， 因此MCSTGSP对噪声越鲁棒，对混响越敏感；极限情况是传声器信号未进行 预白化处理直接用于计算MCCC，因此对应的MCSTGSP退化成原始的MCCC。 因此调节正则化参数δ的不同取值，MCSTGSP方法可构成一组具有不同程度预 白化能力的时延估计器。

以上实验结果表明，提出的MCSTGSP方法在噪声和混响环境下获得了良 好的鲁棒性，在MCCC和MCSTP两种方法间获得了有效的折中。而且，正则 化参数的不同取值使MCSTGSP方法构成一组具有不同预白化能力的时延估计 器，可以根据用户的实际需要来进行调节。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域 的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对 本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定 的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发 明创造均在保护之列。