CN104793177A - 基于最小二乘法的麦克风阵列测向方法 - Google Patents

Abstract

基于最小二乘法的麦克风阵列测向方法，包括以下步骤：由N个阵元组成的麦克风阵列接收空间远场声源信号，获取麦克风阵列的采样数据；获取声源信号到达麦克风阵列各阵元的到达时延差；构造麦克风阵列的距离差矢量；构造麦克风阵列的位置矩阵；判断麦克风阵列的阵元是否位于同一平面，如果是则采用最小二乘方法计算声源信号的二维方向矢量，从而得到声源信号的方位角估计值和俯仰角估计值；否则采用约束最小二乘方法计算声源信号的单位方向矢量估计值，从而得到声源信号的方位角估计值和俯仰角估计值。本发明将声源信号的到达角求解转化为求解相应的方向矢量，解决了宽带声源信号的角度估计问题，降低了角度估计的难度，并提高了角度估计的精确性。

Description

基于最小二乘法的麦克风阵列测向方法

技术领域

本发明属于信号处理技术领域，尤其涉及一种基于麦克风阵列实现语音信号的方位角和俯仰角测量的方法，用于解决宽带语音信号的实时测向问题，可应用于通信或雷达等基于阵列天线实现辐射源测向等领域。

背景技术

麦克风阵列是由按照特定位置摆放的一组麦克风组成的声音信号采集装置，根据实际应用，麦克风阵列信号处理主要解决的是从阵列的输出结果中提取出有用信号或用于进行参数估计。麦克风阵列测向作为语音信号处理的重要研究内容得到了广泛的关注。

基于到达时间差进行声源定位的基本原理是利用了声源到达不同阵元的距离差与声源相对于阵列的方向角有关，因此通常也称为距离差测向方法。由于麦克风阵列一般是由多个阵元构成，因此如何利用所有阵元的测量信息以提高声源的角度估计性能成为当前麦克风阵列研究的重点。现有技术中，利用干涉仪测量不同阵元接收信号之间的相位差，再通过解模糊和相位差解算得到相应信号的入射角度。由于语音信号为宽带信号，因此很难直接利用传统的窄带相位干涉仪方法进行语音信号的角度估计，如果利用空间谱估计方法，必须变换到频域，并进行子带处理，不仅实现复杂、运算量大，而且无法保证全频段处理的一致性。

发明内容

本发明的目的是提供一种可有效提高实时性和精确性的麦克风阵列的测向方法。

为了实现上述目的，本发明采取如下的技术解决方案：

基于最小二乘法的麦克风阵列测向方法，包括以下步骤：由N个阵元组成的麦克风阵列接收空间远场声源信号，N≥3，

步骤一、获取麦克风阵列的采样数据；

步骤二、获取声源信号到达麦克风阵列各阵元的到达时延差；

步骤三、构造麦克风阵列的距离差矢量r；

麦克风阵列的距离差矢量r＝[Δr₂,…,Δr_i,…,Δr_N]^T，其中，Δr_i＝cτ_i为声源信号到达麦克风阵列第i个阵元的距离差，τ_i为声源信号到达第i个阵元的到达时延差，c为声音传播速度，i＝2,…,N；

步骤四、构造麦克风阵列的位置矩阵A；

$A=\left[\begin{array}{ccc}{x}_2-x_1& y_2-y_1& z_2-z_1\\ .& .& .\\ .& .& .\\ .& .& .\\ x_n-x_1& y_n-y_1& z_n-z_1\\ .& .& .\\ .& .& .\\ .& .& .\\ x_N-x_1& y_N-y_1& z_N-z_1\end{array}\right],$

其中，x_n、y_n、z_n为第n个阵元的位置坐标，n＝1,…,N；

声源信号的单位方向矢量声源信号的单位方向矢量与距离差矢量之间的关系为：AΘ＝r，其中，θ为声源信号的方位角，为声源信号的俯仰角，[·]^T表示转置操作；

步骤五、判断麦克风阵列的阵元是否位于同一平面，如果是则执行步骤六，否则执行步骤七；

步骤六、采用最小二乘法计算声源信号的二维方向矢量，从而得到声源信号的方位角估计值和俯仰角估计值；

麦克风阵列的阵元位于同一平面时，位置矩阵A和声源信号的单位方向矢量Θ分别退化为二维位置矩阵A₂和二维方向矢量Θ₂，根据计算声源信号的方位角估计值和俯仰角估计值

$\widehat{\mathrm{\θ}}=a\tan \left(\frac{{\mathrm{\μ}}_2}{{\mathrm{\μ}}_1}\right),$

式中的μ₁为二维方向矢量Θ₂的第一项元素，μ₂为二维方向矢量Θ₂的第二项元素；

步骤七、采用约束最小二乘法计算声源信号的单位方向矢量估计值从而得到声源信号的方位角估计值和俯仰角估计值；

根据计算声源信号的方位角估计值和俯仰角估计值

式中的μ₁为单位方向矢量估计值的第一项元素，μ₂为单位方向矢量估计值的第二项元素，μ₃为单位方向矢量估计值的第三项元素，I为单位矩阵，λ为估计参数。

本发明方法更进一步的方案为，所述步骤5中通过检查位置矩阵A的每一列是否为零来判断阵元是否位于同一平面，如有一列为零，则说明麦克风阵列的阵元分布在同一平面上，否则说明麦克风阵列的阵元不是分布在同一平面上。

本发明方法更进一步的方案为，所述步骤5中通过检查位置矩阵A是否列满秩来判断阵元是否位于同一平面，如位置矩阵A不满足列满秩，则说明麦克风阵列的阵元分布在同一平面上，否则说明麦克风阵列的阵元不是分布在同一平面上。

本发明方法更进一步的方案为，所述估计参数λ通过方程r^TA(A^TA+λI)^-2A^Tr＝1求解，计算步骤如下：

对矩阵A^TA进行特征分解，即其中，Λ＝diag(γ₁,γ₂,γ₃)为特征值矩阵，U＝[u₁,u₂,u₃]为之对应的特征矢量矩阵，N_dim为位置矩阵A的列数，N_dim＝3，u_k为特征值γ_k对应的特征矢量，k＝1,2,3，(·)^H表示转置复共轭操作；

根据 $r^{T}A\·U\stackrel{\mathrm{\Δ}}{=}{p}^T=[p_1,p_2,p_3]$ 及 $U^H{A}^{T}r\stackrel{\mathrm{\Δ}}{=}q={[q_1,q_2,q_3]}^T,$ 得到a_k＝p_kq_k；

估计参数λ通过多项式方程c₆λ⁶+c₅λ⁵+c₄λ⁴+c₃λ³+c₂λ²+c₁λ+c₀＝0进行求解，多项式方程中的系数计算公式如下：

c₆＝1，

c₅＝2γ₁+2γ₂+2γ₃，

$\begin{array}{c}{c}_4=({\mathrm{\γ}}_1^2+{\mathrm{\γ}}_2^2+{\mathrm{\γ}}_3^2+4{\mathrm{\γ}}_1{\mathrm{\γ}}_2+4{\mathrm{\γ}}_1{\mathrm{\γ}}_3+4{\mathrm{\γ}}_2{\mathrm{\γ}}_3)\\ -(a_1+a_2+a_3)\end{array},$

$\begin{array}{c}{c}_3=(2{\mathrm{\γ}}_1{\mathrm{\γ}}_2^2+2{\mathrm{\γ}}_1{\mathrm{\γ}}_3^2+2{\mathrm{\γ}}_1^2{\mathrm{\γ}}_2+2{\mathrm{\γ}}_1^2{\mathrm{\γ}}_3+{2\mathrm{\γ}}_2{\mathrm{\γ}}_3^2+{2\mathrm{\γ}}_2^2{\mathrm{\γ}}_3+8{\mathrm{\γ}}_1{\mathrm{\γ}}_2{\mathrm{\γ}}_3)\\ -[a_1({2\mathrm{\γ}}_2+{2\mathrm{\γ}}_3)+a_2({2\mathrm{\γ}}_1+{2\mathrm{\γ}}_3)+a_3({2\mathrm{\γ}}_1+{2\mathrm{\γ}}_2)]\end{array},$

$\begin{array}{c}{c}_2=({\mathrm{\γ}}_1^2{\mathrm{\γ}}_2^2+{\mathrm{\γ}}_1^2{\mathrm{\γ}}_3^2+{\mathrm{\γ}}_2^2{\mathrm{\γ}}_3^2+{4\mathrm{\γ}}_1{\mathrm{\γ}}_2{\mathrm{\γ}}_3^2+{4\mathrm{\γ}}_1{\mathrm{\γ}}_2^2{\mathrm{\γ}}_3+{4\mathrm{\γ}}_1^2{\mathrm{\γ}}_2{\mathrm{\γ}}_3)\\ -[a_1({\mathrm{\γ}}_2^2+{\mathrm{\γ}}_3^2+{4\mathrm{\γ}}_2{\mathrm{\γ}}_3)+a_2({\mathrm{\γ}}_1^2+{\mathrm{\γ}}_3^2+{4\mathrm{\γ}}_1{\mathrm{\γ}}_3)+a_3({\mathrm{\γ}}_1^2+{\mathrm{\γ}}_2^2+{4\mathrm{\γ}}_1{\mathrm{\γ}}_2)]\end{array},$

$\begin{array}{c}{c}_1=({2\mathrm{\γ}}_1^2{\mathrm{\γ}}_2^2{\mathrm{\γ}}_3+{2\mathrm{\γ}}_1^2{\mathrm{\γ}}_2{\mathrm{\γ}}_3^2+{2\mathrm{\γ}}_1{\mathrm{\γ}}_2^2{\mathrm{\γ}}_3^2)\\ -[a_1({2\mathrm{\γ}}_2{\mathrm{\γ}}_3^2+{2\mathrm{\γ}}_2^2{\mathrm{\γ}}_3)+a_2({2\mathrm{\γ}}_1{\mathrm{\γ}}_3^2+{2\mathrm{\γ}}_1^2{\mathrm{\γ}}_3)+a_3({2\mathrm{\γ}}_1{\mathrm{\γ}}_2^2+{2\mathrm{\γ}}_1^2{\mathrm{\γ}}_2)]\end{array},$

$c_0={\mathrm{\γ}}_1^2{\mathrm{\γ}}_2^2{\mathrm{\γ}}_3^2-[a_1\left({\mathrm{\γ}}_2^2{\mathrm{\γ}}_3^2\right)+a_2\left({\mathrm{\γ}}_1^2{\mathrm{\γ}}_3^2\right)+a_3\left({\mathrm{\γ}}_1^2{\mathrm{\γ}}_2^2\right)];$

利用多项式求根算法，根据以上系数值求得方程的根，将其中模值最小的根作为所求的估计参数λ。

本发明方法将声源信号的到达角求解转化为求解相应的方向矢量，进而得到了声源方向矢量与声源到达各阵元距离差之间的线性方程，有效解决了宽带声源信号的角度估计问题；通过将所求方向角参数转化为方向矢量，实现了将非线性方程转化为线性方程，降低了角度估计的难度，同时提高了角度估计的精确性，为麦克风数量较多时的冗余观测条件下的角度估计提供了有效的处理方法；通过利用声源方向矢量的恒模特性对最小二乘求解方法进行约束，提高了方位角和俯仰角的估计精度。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中需要使用的附图做简单介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明方法的流程图；

图2为仿真实验中采用不同方法计算估计参数的分布图；

图3为不同多项式根下的角度估计误差对比图；

图4a至图4d分别为最大距离误差对到达角估计结果影响曲线对比图；

图5a至图5d分别为阵元圆环半径对角度估计结果影响曲线对比图；

图6a至图6d分别为阵列Z轴坐标的差异性对角度估计结果影响曲线对比图；

图7a至图7d分别为采用最小二乘法估计到达角的估计误差分布图；

图8a至图8d分别为采用约束最小二乘法估计到达角的估计误差分布图。

具体实施方式

为了让本发明的上述和其它目的、特征及优点能更明显，下文特举本发明实施例，并配合所附图示，做详细说明如下。

参照图1，图1为本发明方法的流程图。本发明方法的步骤如下：由N个阵元组成的麦克风阵列接收空间远场声源信号，N为阵元数，N≥3，

步骤一、获取麦克风阵列的采样数据；

步骤二、获取声源信号到达麦克风阵列的各阵元的到达时延差；

利用步骤一获得的麦克风阵列采样数据中的每一通道数据估计声源信号到达各阵元之间的时延差，根据麦克风阵列结构选择参考阵元，采用广义互相关等时延估计法获取声源信号到达其它阵元与参考阵元之间的到达时延差，根据阵列结构选取参考阵元，参考阵元作为第1个阵元，声源信号到达参考阵元(第1个阵元)的到达时延差τ₁＝0；

步骤三、构造麦克风阵列的距离差矢量r；

麦克风阵列的距离差矢量r＝[Δr₂,…,Δr_i,…,Δr_N]^T，其中，Δr_i＝cτ_i为声源信号到达麦克风阵列中第i个阵元的距离差，τ_i为声源信号到达第i个阵元的到达时延差，c为声音传播速度，i＝2,…,N；

步骤四、构造麦克风阵列的位置矩阵A；

$A=\left[\begin{array}{ccc}{x}_2-x_1& y_2-y_1& z_2-z_1\\ .& .& .\\ .& .& .\\ .& .& .\\ x_n-x_1& y_n-y_1& z_n-z_1\\ .& .& .\\ .& .& .\\ .& .& .\\ x_N-x_1& y_N-y_1& z_N-z_1\end{array}\right],$

其中，x_n、y_n、z_n为第n个阵元的位置坐标，n＝1,…,N，阵元的位置坐标根据麦克风阵列的设置参数确定；

声源信号的单位方向矢量 为单位方向矢量Θ的第一项元素的关于方位角和俯仰角的函数关系式，为单位方向矢量Θ的第二项元素的关于方位角和俯仰角的函数关系式，为单位方向矢量Θ的第三项元素的关于方位角和俯仰角的函数关系式，声源信号的单位方向矢量与距离差矢量之间的关系为：AΘ＝r；其中，θ为声源信号的方位角，为声源信号的俯仰角，[·]^T表示转置操作，方位角为入射波达方向在XY平面上的投影与X轴正向之间的夹角，其取值范围为(-π,π]，俯仰角为入射波达方向与XY平面的夹角，其取值范围为[0，π/2]；

步骤五、判断麦克风阵列的阵元是否位于同一平面，如果是则执行步骤六，否则执行步骤七；

判断麦克风阵列的阵元是否位于同一平面，可以通过检查位置矩阵A的每一列是否为零进行判断，如有一列为零，则说明阵元分布在同一平面上，否则说明麦克风阵列的阵元不是分布在同一平面上；或者通过检查位置矩阵A是否列满秩，若位置矩阵A的秩r(A)＝2，即不满足列满秩，则说明麦克风阵列的阵元分布在同一平面上，否则，即r(A)＝3时，位置矩阵A列满秩，则说明麦克风阵列的阵元不是分布在同一平面上；根据不同的阵列分布选择不同的角度估计方法；

步骤六、采用最小二乘法计算声源信号的二维方向矢量，从而得到声源信号的方位角估计值和俯仰角估计值；

当麦克风阵列的阵元位于同一平面时，则有一维坐标完全相同，位置矩阵A相应列的元素即为0，此时声源信号的位置矩阵A和单位方向矢量Θ分别退化为二维位置矩阵A₂和二维方向矢量Θ₂，本实施例以Z轴坐标相同为例进行说明，则

$A_2=\left[\begin{array}{cc}{x}_2-x_1& y_2-y_1\\ .& .\\ .& .\\ .& .\\ x_n-x_1& y_n-y_1\\ .& .\\ .& .\\ .& .\\ x_N-x_1& y_N-y_1\end{array}\right],$

相应的，当X轴坐标相同时， $A_2=\left[\begin{array}{cc}{y}_2-y_1& z_2-z_1\\ .& .\\ .& .\\ .& .\\ y_n-y_1& z_n-z_1\\ .& .\\ .& .\\ .& .\\ y_N-y_1& z_N-z_1\end{array}\right],$ 当Y轴坐标相同时， $A_2=\left[\begin{array}{cc}{x}_2-x_1& z_2-z_1\\ .& .\\ .& .\\ .& .\\ x_n-x_1& z_n-z_1\\ .& .\\ .& .\\ .& .\\ x_N-x_1& z_N-z_1\end{array}\right],$

根据计算声源信号的方位角估计值和俯仰角估计值

$\widehat{\mathrm{\θ}}=a\tan \left(\frac{{\mathrm{\μ}}_2}{{\mathrm{\μ}}_1}\right),$

式中的μ₁为二维方向矢量Θ₂的第一项元素，μ₂为二维方向矢量Θ₂的第二项元素，r为麦克风阵列的距离差矢量，A₂为麦克风阵列的二维位置矩阵，本实施例中μ₁和μ₂对应方位角和俯仰角的表达式分别为：

步骤七、采用约束最小二乘法计算声源信号的单位方向矢量估计值从而得到声源信号的方位角估计值和俯仰角估计值；

当各阵元不是处于同一平面时，此时的位置矩阵A是列满秩，即r(A)＝3，为了进一步提高声源信号的测向精度，利用声源信号方向矢量的恒模特性进行约束，利用公式 $\left\{\begin{array}{c}\min \left|\right|\mathrm{A\Θ}-r\left|\right|\\ s.t.\left|\right|\mathrm{\Θ}\left|\right|=1\end{array}\right.$ 对方向矢量进行优化求解，得到单位方向矢量估计值 $\widehat{\mathrm{\Θ}}={(A^{T}A+\mathrm{\λI})}^{-1}{A}^{T}r;$

根据计算声源信号的方位角估计值和俯仰角估计值

式中的I为单位矩阵，μ₁为单位方向矢量估计值的第一项元素，μ₂为单位方向矢量估计值的第二项元素，μ₃为单位方向矢量估计值的第三项元素，μ₁，μ₂和μ₃对应方位角和俯仰角的表达式分别为： r为麦克风阵列的距离差矢量，A为麦克风阵列的位置矩阵，λ为估计参数。

俯仰角估计值可用公式或公式计算，也可同时采用前述两个公式计算俯仰角，然后取平均值作为最后的俯仰角估计值，以取得更为准确、可靠的估计结果。

本发明的估计参数λ可通过方程r^TA(A^TA+λI)^-2A^Tr＝1求解，计算步骤如下：

对矩阵A^TA进行特征分解，即其中，Λ＝diag(γ₁,γ₂,γ₃)为特征值矩阵，U＝[u₁,u₂,u₃]为之对应的特征矢量矩阵，N_dim为位置矩阵A的列数，由于A是列满秩的，因此N_dim＝3，u_k为特征值γ_k对应的特征矢量，k＝1,2,3，(·)^H表示转置复共轭操作；

根据 $r^{T}A\·U\stackrel{\mathrm{\Δ}}{=}{p}^T=[p_1,p_2,p_3]$ 及 $U^H{A}^{T}r\stackrel{\mathrm{\Δ}}{=}q={[q_1,q_2,q_3]}^T,$ 得到a_k＝p_kq_k，k＝1,2,3，p_k为矢量p的第k项元素，q_k为矢量q的第k项元素；

本发明的估计参数λ通过多项式方程c₆λ⁶+c₅λ⁵+c₄λ⁴+c₃λ³+c₂λ²+c₁λ+c₀＝0进行求解，多项式方程中各系数计算公式如下：

c₆＝1，

c₅＝2γ₁+2γ₂+2γ₃，

$\begin{array}{c}{c}_4=({\mathrm{\γ}}_1^2+{\mathrm{\γ}}_2^2+{\mathrm{\γ}}_3^2+4{\mathrm{\γ}}_1{\mathrm{\γ}}_2+4{\mathrm{\γ}}_1{\mathrm{\γ}}_3+4{\mathrm{\γ}}_2{\mathrm{\γ}}_3)\\ -(a_1+a_2+a_3)\end{array},$

$\begin{array}{c}{c}_3=(2{\mathrm{\γ}}_1{\mathrm{\γ}}_2^2+2{\mathrm{\γ}}_1{\mathrm{\γ}}_3^2+2{\mathrm{\γ}}_1^2{\mathrm{\γ}}_2+2{\mathrm{\γ}}_1^2{\mathrm{\γ}}_3+{2\mathrm{\γ}}_2{\mathrm{\γ}}_3^2+{2\mathrm{\γ}}_2^2{\mathrm{\γ}}_3+8{\mathrm{\γ}}_1{\mathrm{\γ}}_2{\mathrm{\γ}}_3)\\ -[a_1({2\mathrm{\γ}}_2+{2\mathrm{\γ}}_3)+a_2({2\mathrm{\γ}}_1+{2\mathrm{\γ}}_3)+a_3({2\mathrm{\γ}}_1+{2\mathrm{\γ}}_2)]\end{array},$

$\begin{array}{c}{c}_2=({\mathrm{\γ}}_1^2{\mathrm{\γ}}_2^2+{\mathrm{\γ}}_1^2{\mathrm{\γ}}_3^2+{\mathrm{\γ}}_2^2{\mathrm{\γ}}_3^2+{4\mathrm{\γ}}_1{\mathrm{\γ}}_2{\mathrm{\γ}}_3^2+{4\mathrm{\γ}}_1{\mathrm{\γ}}_2^2{\mathrm{\γ}}_3+{4\mathrm{\γ}}_1^2{\mathrm{\γ}}_2{\mathrm{\γ}}_3)\\ -[a_1({\mathrm{\γ}}_2^2+{\mathrm{\γ}}_3^2+{4\mathrm{\γ}}_2{\mathrm{\γ}}_3)+a_2({\mathrm{\γ}}_1^2+{\mathrm{\γ}}_3^2+{4\mathrm{\γ}}_1{\mathrm{\γ}}_3)+a_3({\mathrm{\γ}}_1^2+{\mathrm{\γ}}_2^2+{4\mathrm{\γ}}_1{\mathrm{\γ}}_2)]\end{array},$

$\begin{array}{c}{c}_1=({2\mathrm{\γ}}_1^2{\mathrm{\γ}}_2^2{\mathrm{\γ}}_3+{2\mathrm{\γ}}_1^2{\mathrm{\γ}}_2{\mathrm{\γ}}_3^2+{2\mathrm{\γ}}_1{\mathrm{\γ}}_2^2{\mathrm{\γ}}_3^2)\\ -[a_1({2\mathrm{\γ}}_2{\mathrm{\γ}}_3^2+{2\mathrm{\γ}}_2^2{\mathrm{\γ}}_3)+a_2({2\mathrm{\γ}}_1{\mathrm{\γ}}_3^2+{2\mathrm{\γ}}_1^2{\mathrm{\γ}}_3)+a_3({2\mathrm{\γ}}_1{\mathrm{\γ}}_2^2+{2\mathrm{\γ}}_1^2{\mathrm{\γ}}_2)]\end{array},$

$c_0={\mathrm{\γ}}_1^2{\mathrm{\γ}}_2^2{\mathrm{\γ}}_3^2-[a_1\left({\mathrm{\γ}}_2^2{\mathrm{\γ}}_3^2\right)+a_2\left({\mathrm{\γ}}_1^2{\mathrm{\γ}}_3^2\right)+a_3\left({\mathrm{\γ}}_1^2{\mathrm{\γ}}_2^2\right)];$

上式中，γ_k为特征值矩阵Λ的第k个特征值，k＝1,2,3，利用多项式求根算法，根据以上系数值求得方程的根，由于多项式方程共有六个根λ₁、λ₂、λ₃、λ₄、λ₅、λ₆，将其中模值最小的根作为所求的估计参数λ。

本发明的效果可以通过以下的仿真结果进一步说明：

实验一、

仿真条件如下：麦克风阵列为阵元在XY平面为均匀圆环阵，相邻阵元的Z轴坐标具有一定的差异，仿真中的误差是通过加入到目标与阵元之间的距离差测量值中，其中误差的最大值为R/10，阵元数N＝8，所有阵元的X、Y坐标分布在R＝0.2m的圆环上，Z轴的坐标不同，其中奇数序号阵元Z＝0，而偶数序号阵元Z＝R/5。声源信号的方位角和俯仰角为(70°,30°)。

图2所示为采用不同方法计算估计参数λ的分布图。其中标注为Fun-Opt的+为采用约束条件搜索得到的估计参数λ，标注为Err-Opt的○为采用全局搜索最小误差法得到的估计参数λ，标注为Roots的×为采用本发明方法的多项式求根法获得的估计参数λ。

从图2可以看出，基于多项式求根方法计算最优λ时，R6的模值最小，距离(0,0)最近，结合图3，应该选取模值最小的根作为所求的最优λ值。然而基于全局搜索方位和俯仰角估计误差最小而搜索到的最优λ为复数，且与利用约束条件和多项式求根法得到的λ具有一定的差别，这是由于两者的优化目标不同造成的。而基于约束条件获得的λ与多项式求根法获得的λ之间的差异性则是由搜索过程中的步长或搜索精度引起的，其中基于多项式求根法得到的最优λ最准确，而且速度更快。

图3为不同多项式根下的角度估计误差对比图，图3中CLSMMerr和CLSMfun的曲线重合，从图3可以看出，基于不同多项式根作为λ的角度估计误差，第六个根，也就是模值最小的根对应的误差最小。

实验二、

实验二的麦克风阵列阵元的排列与实验一相同。实验二将最小二乘法与约束最小二乘法进行对比。由于当阵元的Z轴坐标不一致时，俯仰角的估计有两个公式可以利用，因此，将利用反余弦函数计算的方法标注为-F1，将利用反正弦函数计算的方法标注为-F2，将两者计算结果平均的方法标注为-MF12，采用最小二乘法的曲线标注为LSM，采用约束最小二乘法的曲线标注为CLSM。

图4a为最大距离误差对方位角估计结果影响曲线对比图，图4b至图4d分别为最大距离误差对俯仰角估计结果影响对比曲线图。从图4a至图4d可以看出，模值约束对方位角和俯仰角估计结果具有较大的改善，其中基于反正弦函数计算结果改善最明显，从该仿真结果也可看出，约束最小二乘法中的三种俯仰角估计结果相同，而最小二乘法中利用反余弦函数计算的俯仰角精度最高。

图5a至图5d分别为阵元圆环半径对角度估计结果影响曲线对比图。从图5a至图5d可以看出，随着圆环半径的增加，最小二乘和约束最小二乘的方位角和俯仰角估计误差逐渐变小，约束最小二乘法对方位角和俯仰角的估计结果具有一定改善，均优于最小二乘法对方位角和俯仰角的估计结果，其中利用反正弦函数方法进行估计的结果改善最明显。

图6a至图6d分别为阵列Z轴坐标的差异性对角度估计结果影响曲线对比图，从图6a至图6d可以看出，随着Z轴坐标的差异性增加，估计误差逐渐变小，而模值约束对方位角估计结果具有一定改善，对俯仰角估计改善比较大，尤其在阵列Z轴差异性比较小时，利用反正弦方法估计结果改善最明显。

图7a至图7d分别为采用最小二乘法估计到达角的估计误差分布图。图8a至图8d分别为采用约束最小二乘法估计到达角的估计误差分布图。通过对比可以发现，模值约束不仅可以提高方位角的估计精度，对于俯仰角的估计精度提高最明显，相比于最小二乘算法，方位角和俯仰角估计精度都有较大提升，尤其是两种俯仰角估计方法获得的误差分布几乎一样。

通过实验验证，本发明方法不仅测向速度快，而且精度高，尤其当利用约束加权最小二乘方法时，具有较高的方位角和俯仰角估计性能，可应用于语音信号的定位与跟踪场景。

以上所述，仅是本发明的较佳实施例而已，并非对本发明做任何形式上的限制，虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容做出些许更动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。

Claims (4)

1.基于最小二乘法的麦克风阵列测向方法，其特征在于，包括以下步骤：由N个阵元组成的麦克风阵列接收空间远场声源信号，N≥3，

步骤一、获取麦克风阵列的采样数据；

步骤二、获取声源信号到达麦克风阵列各阵元的到达时延差；

步骤三、构造麦克风阵列的距离差矢量r；

麦克风阵列的距离差矢量r＝[Δr₂,…,Δr_i,…,Δr_N]^T，其中，Δr_i＝cτ_i为声源信号到达麦克风阵列第i个阵元的距离差，τ_i为声源信号到达第i个阵元的到达时延差，c为声音传播速度，i＝2,…,N；

步骤四、构造麦克风阵列的位置矩阵A；

其中，x_n、y_n、z_n为第n个阵元的位置坐标，n＝1,…,N；

声源信号的单位方向矢量声源信号的单位方向矢量与距离差矢量之间的关系为：AΘ＝r，其中，θ为声源信号的方位角，为声源信号的俯仰角，[·]^T表示转置操作；

步骤五、判断麦克风阵列的阵元是否位于同一平面，如果是则执行步骤六，否则执行步骤七；

步骤六、采用最小二乘法计算声源信号的二维方向矢量，从而得到声源信号的方位角估计值和俯仰角估计值；

麦克风阵列的阵元位于同一平面时，位置矩阵A和声源信号的单位方向矢量Θ分别退化为二维位置矩阵A₂和二维方向矢量Θ₂，根据计算声源信号的方位角估计值和俯仰角估计值

$\widehat{\mathrm{\θ}}=a\tan \left(\frac{{\mathrm{\μ}}_2}{{\mathrm{\μ}}_1}\right),$

式中的μ₁为二维方向矢量Θ₂的第一项元素，μ₂为二维方向矢量Θ₂的第二项元素；

步骤七、采用约束最小二乘法计算声源信号的单位方向矢量估计值从而得到声源信号的方位角估计值和俯仰角估计值；

根据计算声源信号的方位角估计值和俯仰角估计值

式中的μ₁为单位方向矢量估计值的第一项元素，μ₂为单位方向矢量估计值的第二项元素，μ₃为单位方向矢量估计值的第三项元素，I为单位矩阵，λ为估计参数。

2.根据权利要求1所述的基于最小二乘法的麦克风阵列测向方法，其特征在于：所述步骤5中通过检查位置矩阵A的每一列是否为零来判断阵元是否位于同一平面，如有一列为零，则说明麦克风阵列的阵元分布在同一平面上，否则说明麦克风阵列的阵元不是分布在同一平面上。

3.根据权利要求1所述的基于最小二乘法的麦克风阵列测向方法，其特征在于：所述步骤5中通过检查位置矩阵A是否列满秩来判断阵元是否位于同一平面，如位置矩阵A不满足列满秩，则说明麦克风阵列的阵元分布在同一平面上，否则说明麦克风阵列的阵元不是分布在同一平面上。

4.根据权利要求1所述的基于最小二乘法的麦克风阵列测向方法，其特征在于：所述估计参数λ通过方程r^TA(A^TA+λI)^-2A^Tr＝1求解，计算步骤如下：

对矩阵A^TA进行特征分解，即其中，

Λ＝diag(γ₁,γ₂,γ₃)为特征值矩阵，U＝[u₁,u₂,u₃]为之对应的特征矢量矩阵，N_dim为位置矩阵A的列数，N_dim＝3，u_k为特征值γ_k对应的特征矢量，k＝1,2,3，(·)^H表示转置复共轭操作；

根据 $r^{T}A\·U\stackrel{\mathrm{\Δ}}{=}{p}^T=[p_1,p_2,p_3]$ 及 $U^H\·A^{T}r\stackrel{\mathrm{\Δ}}{=}q=[q_1,q_2,q_3{]}^T,$ 得到a_k＝p_kq_k；

估计参数λ通过多项式方程c₆λ⁶+c₅λ⁵+c₄λ⁴+c₃λ³+c₂λ²+c₁λ+c₀＝0进行求解，多项式方程中的系数计算公式如下：

c₆＝1，

c₅＝2γ₁+2γ₂+2γ₃，

$\begin{array}{c}{c}_4=({\mathrm{\γ}}_1^2+{\mathrm{\γ}}_2^2+{\mathrm{\γ}}_3^2+{4\mathrm{\γ}}_1{\mathrm{\γ}}_2+4{\mathrm{\γ}}_1{\mathrm{\γ}}_3+4{\mathrm{\γ}}_2{\mathrm{\γ}}_3)\\ -(a_1+a_2+a_3)\end{array},$

$\begin{array}{c}{c}_3=(2{\mathrm{\γ}}_1{\mathrm{\γ}}_2^2+2{\mathrm{\γ}}_1{\mathrm{\γ}}_3^2+2{\mathrm{\γ}}_1^2{\mathrm{\γ}}_2+2{\mathrm{\γ}}_1^2{\mathrm{\γ}}_3+2{\mathrm{\γ}}_2{\mathrm{\γ}}_3^2+2{\mathrm{\γ}}_2^2{\mathrm{\γ}}_3+8{\mathrm{\γ}}_1{\mathrm{\γ}}_2{\mathrm{\γ}}_3)\\ -[a_1(2{\mathrm{\γ}}_2+2{\mathrm{\γ}}_3)+a_2(2{\mathrm{\γ}}_1+2{\mathrm{\γ}}_3)+a_3(2{\mathrm{\γ}}_1+2{\mathrm{\γ}}_2)]\end{array},$

$\begin{array}{c}{c}_2=({\mathrm{\γ}}_1^2{\mathrm{\γ}}_2^2+{\mathrm{\γ}}_1^2{\mathrm{\γ}}_3^2+{\mathrm{\γ}}_2^2+{\mathrm{\γ}}_3^2+4{\mathrm{\γ}}_1{\mathrm{\γ}}_2^2{\mathrm{\γ}}_3+4{\mathrm{\γ}}_1^2{\mathrm{\γ}}_2{\mathrm{\γ}}_3)\\ -[a_1({\mathrm{\γ}}_2^2+{\mathrm{\γ}}_3^2+4{\mathrm{\γ}}_2{\mathrm{\γ}}_3)+a_2({\mathrm{\γ}}_1^2+{\mathrm{\γ}}_3^2+4{\mathrm{\γ}}_1{\mathrm{\γ}}_3)+a_3({\mathrm{\γ}}_1^2+{\mathrm{\γ}}_2^2+4{\mathrm{\γ}}_1{\mathrm{\γ}}_2)]\end{array},$

$\begin{array}{c}{c}_1=(2{\mathrm{\γ}}_1^2{\mathrm{\γ}}_2^2{\mathrm{\γ}}_3+2{\mathrm{\γ}}_1^2{\mathrm{\γ}}_2{\mathrm{\γ}}_3^2+2{\mathrm{\γ}}_1{\mathrm{\γ}}_2^2{\mathrm{\γ}}_3^2)\\ -[a_1(2{\mathrm{\γ}}_2{\mathrm{\γ}}_3^2+2{\mathrm{\γ}}_2^2{\mathrm{\γ}}_3)+a_2(2{\mathrm{\γ}}_1{\mathrm{\γ}}_3^2+2{\mathrm{\γ}}_1^2{\mathrm{\γ}}_3)+a_3(2{\mathrm{\γ}}_1{\mathrm{\γ}}_2^2+2{\mathrm{\γ}}_1^2{\mathrm{\γ}}_2)]\end{array},$

$c_0={\mathrm{\γ}}_1^2{\mathrm{\γ}}_2^2{\mathrm{\γ}}_3^2-[a_1\left({\mathrm{\γ}}_2^2{\mathrm{\γ}}_3^2\right)+a_2\left({\mathrm{\γ}}_1^2{\mathrm{\γ}}_3^2\right)+a_3\left({\mathrm{\γ}}_1^2{\mathrm{\γ}}_2^2\right)];$

利用多项式求根算法，根据以上系数值求得方程的根，将其中模值最小的根作为所求的估计参数λ。