基于低秩及几何约束的嵌入式平面麦克风阵列声源测向方法与流程

本发明属于嵌入式麦克风阵列声源测向领域，特别是一种基于低秩及几何约束的嵌入式平面麦克风阵列声源测向方法。

背景技术：

麦克风阵列声源测向方法很多，包括高分辨谱估计、波束形成、波达时间差(tdoa)等。其中，麦克风间tdoa测量结合最小二乘(ls)估计的两级声源测向方法因实现复杂度低、环境适应性好等优点在各种实时声探测任务中得到广泛应用。因受外形设计、空间尺寸、使用运输等多种条件限制，图3所示的便携式枪声探测系统所用的麦克风阵列通常采用嵌入式安装方式，可能导致由于外壳遮挡而造成目标与某个或某些麦克风之间处于非通视(nlos)状态。考虑到经过远距离传播后的枪声膛口波信号能量主要集中在数百hz的低频范围，非通视麦克风受壳体衍射传播效应影响，相关麦克风间的实际tdoa不再与麦克风间距和目标位置保持理想全通视条件下的理论对应关系，进而造成基于ls测向方法的膛口波方向估计结果出现明显偏差。

有关非通视情况下的声源目标测向问题，许多文献主要关注声探测系统整体因障碍阻隔而无法直视声源目标时的解决方案，一般不涉及壳体遮挡对嵌入式阵列目标方向估计的影响及其处理方法。针对仅有部分非通视麦克风的情形，已有研究表明，采用球麦克风阵列声场模型能有效抑制衍射传播干扰、获得精确的声源目标方位估计。但基于模型的衍射传播效应抑制方法一般较为复杂，难以实时实现；同时受成本、结构、使用方式等制约，便携式枪声探测系统通常采用麦克风数很少的平面阵结构，很难严格满足均匀球形阵列甚至圆形阵列的模型要求。

技术实现要素：

本发明所要解决的技术问题在于针对嵌入式平面麦克风阵列的声源测向误差，提供一种实时性高、易于实现、且测向较为精确的嵌入式平面麦克风阵列声源测向方法。

实现本发明目的的技术解决方案为：基于低秩及几何约束的嵌入式平面麦克风阵列声源测向方法，包括以下步骤：

步骤1、通过嵌入式平面麦克风阵列获取实测数据，并对实测数据进行波达时间差估计；

步骤2、利用步骤1中波达时间差估计值，根据波达时间的大小关系选择通视麦克风作为参考麦克风；

步骤3、基于步骤2中选择的参考麦克风构建波达时间差矩阵，进行施加低秩及几何约束条件下的凸优化校正；

步骤4、根据步骤3中校正后的波达时间差结果构建波达时间差向量，进而进行最小二乘方位角估计，完成嵌入式平面麦克风阵列声源测向。

本发明与现有技术相比，其显著优点为：1)本发明中通过选择距离声源最近的麦克风、施加低秩以及几何约束，降低了测向误差的均值及标准差，提高了测向精度及稳健性；2)本发明中选择距离声源最近的麦克风，尽可能避免了壳体遮挡对测向性能的影响；3)本发明采用反对称tdoa矩阵的低秩约束条件，抑制了实测tdoa的测量随机误差；4)本发明采用阵列形状的几何约束条件，抑制了壳体遮挡引起的衍射传播误差；5)本发明中将tdoa的校正转化为带约束条件的凸优化问题，利用拉格朗日乘子法求解凸优化问题的解析解，使得该方法的计算量较小，可满足系统的实时性需要，也降低了成本；6)本发明中通过将tdoa校正转化为带约束条件的凸优化问题，仅通过改变各矩阵的维数和几何约束的表现形式，就可针对性的推广到麦克风数量或阵列形状不同的其它麦克风阵列系统。

下面结合附图对本发明作进一步详细的描述。

附图说明

图1为本发明基于低秩及几何约束的嵌入式平面麦克风阵列声源测向方法流程图。

图2为本发明实施例中实验布设示意图。

图3为本发明实施例所用的嵌入式平面麦克风阵列示意图。

图4为本发明实施例中参考麦克风选择的结构框图。

图5为本发明方法和常规方法的方位角估计误差分布对比示意图。

图6为本发明方法和常规方法的方位角估计误差分布直方图，其中(a)为常规方法的方位角估计误差分布直方图，(b)为本发明方法的方位角估计误差分布直方图。

图7为本实施例中均值补偿后两种方法的绝对方位角估计误差累积概率分布示意图。

具体实施方式

结合图1，本发明基于低秩及几何约束的嵌入式平面麦克风阵列声源测向方法，包括以下步骤：

步骤1、通过嵌入式平面麦克风阵列获取实测数据，并对实测数据进行波达时间差估计。

步骤2、利用步骤1中波达时间差估计值，根据波达时间的大小关系选择通视麦克风作为参考麦克风。

进一步地，步骤2具体为：

步骤2-1、根据波达时间差估计值构建波达时间差估计矩阵，具体为：

令远场目标信号传播至麦克风mi的波达时间为ti，i＝1～n；定义嵌入式平面麦克风阵列的波达时间向量为x＝[t1,t2,t3,…,tn]^t，远场目标信号传播至麦克风mi和mj之间的波达时间差τij＝ti-tj，结合维度为n的全1向量1＝[1,1,1,…,1]^t可得相应的反对称波达时间差矩阵m：

m＝x·1^t-1·x^t

结合上述波达时间差矩阵m及实际测量的波达时间向量x获得实测波达时间差矩阵

步骤2-2、记录步骤2-1构建的实测波达时间差矩阵中每列中最小值所在的行序号；

步骤2-3、以记录最多的行序号即众数作为距离目标最近的通视麦克风的序号，并选取该麦克风作为参考麦克风。

进一步地，n的取值为4。

步骤3、基于步骤2中选择的参考麦克风构建波达时间差矩阵，进行施加低秩及几何约束条件下的凸优化校正。

进一步地，步骤3具体为：

步骤3-1、将嵌入式平面麦克风阵列中的参考麦克风记为m1，其余麦克风以逆时针为序依次命名，重新构建实测波达时间差矩阵

步骤3-2、对实测波达时间差矩阵进行凸优化校正，获得校正后的理想全通视波达时间差矩阵估计的解析解；具体为：

施加波达时间差矩阵代数特征的低秩约束，以及对应阵列形状特征的几何约束，将对麦克风全通视情况下的理想波达时间差矩阵m的估计任务转化为有约束凸优化问题：

式中，||·||f为f-范数，约束矩阵d＝[1,g]，包括低秩约束1^tx＝0和几何约束g^tx＝0，前者与理想全通视tdoa矩阵m应为秩2的低秩要求相关，用于抑制高斯随机测量噪声，后者则体现了平行四边形阵列所对应的形状特征几何约束，其中g称为约束向量；

对上式应用拉格朗日乘数法构造函数：

式中，待定系数λ1、λ2为拉格朗日乘子；

考虑：

分别对x、λ1、λ2的一阶导数为0，可得

结合的反对称性质：以及可得

由上获得校正后麦克风全通视情况下的理想波达时间差矩阵m估计的解析解

步骤4、根据步骤3中校正后的波达时间差结果构建波达时间差向量，从而完成最小二乘方位角估计。

进一步地，步骤4具体为：

步骤4-1、构建目标方向对应的麦克风间波达时间差向量b为：

式中，a为与嵌入式平面麦克风阵列尺寸及形状有关的阵列系数矩阵，为远场目标方向矢量，θ、分别为远场目标方位角和俯仰角，τi1为第i个麦克风的坐标向量与第1个麦克风的坐标向量的差值，mi为第i个麦克风的坐标向量，c为声速；

步骤4-2、求取上式向量b的线性最小二乘解为：

式中，为阵列系数矩阵a的摩尔-彭若斯广义逆矩阵，分别为的估计值；

步骤4-3、进行最小二乘方位角估计为：

根据和各自的正负号，最终在0°～360°范围内唯一确定目标方位角θ的估计结果。

下面结合实施例对本发明作进一步详细的描述。

实施例

结合图1，本实施例基于低秩及几何约束的嵌入式平面麦克风阵列声源测向方法，包括以下内容：

步骤1、本实施例的实验在较空旷的室外环境进行，实验布设示意图如图2所示。实施例中采用的嵌入式平面麦克风阵列示意图如图3所示，其包括4个麦克风，嵌入式平面麦克风阵列系统被水平固定在三脚架上，以阵列中心为坐标原点，通过水平顺时针旋转阵列系统等效从小到大改变目标方位角。实施例中作为模拟目标的扬声器与嵌入式平面麦克风阵列高度基本相同，放置在距离阵列中心15m处。本实施例中，目标方位角从0°开始以30°步进递增变化一周，在每个目标方位角位置进行26次独立实验，12个方位角位置共计进行312次独立实验。每次实验中扬声器播放一段含有外场实测的枪声膛口波信号录音。对312次独立实验的实测数据各自进行波达时间差估计。

步骤2、对于一组波达时间差估计值，显然通视麦克风离声源最近，因此根据波达时间的大小关系选择通视麦克风，并作为参考麦克风。结合图4，具体如下：

步骤2-1、令嵌入式平面麦克风阵列系统中的其中一个麦克风为m1，其余麦克风以俯视逆时针为序依次命名m2、m3、m4，远场目标信号传播至麦克风mi的波达时间为ti，i＝1～4；定义阵列的波达时间向量为x＝[t1,t2,t3,t4]^t，远场目标信号传播至麦克风mi和mj之间的波达时间差τij＝ti-tj，i、j＝1～4，结合全1向量1＝[1,1,1,1]^t可得相应的反对称波达时间差矩阵m：

结合上述波达时间差矩阵m及实际测量的波达时间向量x获得实测波达时间差矩阵

步骤2-2、记录步骤2-1构建的实测波达时间差矩阵中每列最小值所在的行序号tminj，j＝1～4，定义行序号向量tmin＝[tmin1,tmin2,tmin3,tmin4]；

步骤2-3、以记录最多的行序号即tmin中的众数作为距离目标最近的通视麦克风序号，并选取该麦克风作为参考麦克风。

步骤3、基于步骤2中选择的参考麦克风构建波达时间差矩阵，进行施加低秩及几何约束条件下的凸优化校正，具体如下：

步骤3-1、将嵌入式平面麦克风阵列中的参考麦克风记为m1，其余麦克风以逆时针为序依次命名m2、m3、m4，重新构建实测波达时间差矩阵为：

步骤3-2、对实测波达时间差矩阵进行凸优化校正，获得校正后的理想全通视波达时间差矩阵估计的解析解，具体为：

施加波达时间差矩阵代数特征的低秩约束，以及对应阵列形状特征的几何约束，将对麦克风全通视情况下的理想波达时间差矩阵m的估计任务转化为有约束凸优化问题：

式中，||·||f为f-范数，约束矩阵d＝[1,g]，包括低秩约束1^tx＝0和几何约束g^tx＝0，前者与理想全通视波达时间差矩阵m应为秩2的低秩要求相关，用于抑制高斯随机测量噪声，后者则体现了平行四边形阵列所对应的形状特征几何约束，由于τ12＝τ43，因而式中g＝[1,-1,1,-1]^t。

针对上式应用lagrange乘数法构造函数：

式中，待定系数λ1、λ2为拉格朗日乘子；

考虑：

分别对x、λ1、λ2的一阶导数为0，可得

结合的反对称性质：以及可得

最后可得校正后麦克风全通视情况下的理想波达时间差矩阵m估计的解析解

步骤4、根据步骤3中校正后的波达时间差结果构建波达时间差向量，从而完成最小二乘方位角估计，具体如下：

步骤4-1、对于嵌入式阵列n＝4，构建目标方向对应的麦克风间波达时间差向量b为：

式中，a为与阵列尺寸及形状有关的阵列系数矩阵，为远场目标方向矢量，θ和分别为远场目标方位角和俯仰角，τi1为第i个麦克风的坐标向量与第1个麦克风的坐标向量的差值，mi为第i个麦克风的坐标向量，c为声速；

步骤4-2、求取上式向量b的线性最小二乘解为：

式中，为阵列系数矩阵a的摩尔-彭若斯广义逆矩阵，分别为的估计值；

步骤4-3、进行最小二乘方位角估计为：

arctan(·)表示四象限反正切函数，即需要根据和各自的正负号，最终在0°～360°范围内唯一确定目标方位角的估计结果。具体如下：

当

当

当

当

步骤5、采用理想全通视直线传播模型下基于tdoa的常规ls测向算法(下文简称常规方法)作为性能对比，并定义方位角估计误差为目标方位角估计值与真实值之差。

图5为本发明方法和常规方法在12个测试方位角上各26次独立实验的方位角估计误差分布。图中显示，常规方法的方位角估计误差在一些方位角处散布较大，而经过本发明方法校正后，各方位角的方位角估计误差分布则已趋向集中。

图6比较了两种方法的方位角估计误差总体分布直方图。对比图6(a)、(b)可以看出，常规方法的方位角估计误差主要集中分布在–2.03°～5.14°范围内(上下四分位数误差范围)，分布均匀程度较高，且误差绝对值的最大值高达23°。而本发明方法方位角估计误差则主要集中分布在–1.38°～2.84°范围内(上下四分位数误差范围)，相应频度随着误差绝对值增大而很快下降，误差绝对值的最大值也降为10°。

此外，常规方法和本发明方法的方位角估计误差统计均值±方差分别为1.73±5.87°和0.90±3.55°，统计中值分别为1.88°和0.87°。实际目标方位角通常未知且随机，通常假设实际目标方位角在0°～360°之间均匀分布。为使全方位的总体测角性能相对均衡，同时考虑到两种方法方位角估计误差的统计均值与统计中值基本相等，实际应用时可预先补偿上述统计均值。两种方法在均值补偿后的方位角估计误差绝对值所对应的累积概率分布对比如图7所示。需要强调的是，就枪声探测系统实战环境下的技术指标而言，当前以英国神听猛击(earsswats)ii代、美国飞镖战警-x(boomerangwarrior-x)系统为代表的国外同类便携式产品的枪声膛口波测向性能通常以±7.5°为测向精度考核标准。通过考察图7可知，针对方位角估计误差绝对值在7.5°以内的累积概率，常规方法仅为0.86，而与之形成鲜明对比的是，本发明方法达到0.94，测向性能提升明显。

综上所述，本发明的方法通过在低秩约束的基础上增加阵列形状几何约束来构建凸优化问题并选取距离目标最近的通视麦克风作为参考麦克风，能够同时有效抑制tdoa随机测量噪声与壳体遮挡衍射传播效应对基于tdoa的ls测向方法的不利影响，方位角估计精度得到明显提高。需要说明的是，仅通过改变各矩阵或向量的维数以及几何约束的表现形式，上述测向方法就可直接推广到麦克风数量或阵列形状各不相同的其它嵌入式麦克风阵列系统。