基于概率数据关联的移动声源到达时延估计方法与流程

本发明涉及阵列信号领域，特别是涉及麦克风阵列信号处理和到达时延估计。

背景技术：

本发明的背景基于实际需要产生。近年来，在无人机迅速成为研究热点时，也带来一系列问题，比如无人机黑飞，严重影响区域安全。因此无人机防御正在成为各国政府和军方重点关注的新领域。无人机的声音具有明显的特征，可以有效侦别出空中飞行的物体。通过架设多个麦克风阵列，能有效地收集到无人机的信号，但是实际环境存在的噪声、信号在电路设备传输过程中的噪声以及一些其他的混响、多径干扰等，收集到的声音信号非常嘈杂，导致对无人机进行定位时，计算得到的信号到达麦克风的时延(tdoa)会出现比较大的误差，为了提高tdoa估计的准确性，需要对tdoa估计算法进行改进。

在目前的研究内容中，常用的改进方式是结合前一时刻的tdoa的信息，将本时刻的tdoa估计函数的取值范围限定在一个区间内，从而减小tdoa估计的波动，但这种方法对初值要求较高，且对无人机的飞行速度有要求。因此，急需一种新的方法来真正有效的估计出tdoa，同时能够保证较快的运行速度和系统较低的成本。

技术实现要素：

为了实现对无人机移动时声音信号时延的跟踪，本发明采用麦克风阵列传感器对空中无人机声音信号进行处理，可以有效估计移动中的无人机声音信号到达麦克风时延。

本发明解决其技术问题所采用的技术方案是：一种基于概率数据关联的移动声源到达时延估计方法，包括如下步骤：

(1)根据广义互相关函数计算麦克风阵列t时刻采集的两组时域声信号x1(t)、x2(t)互功率谱密度分布函数

其中为x1(t)和x2(t)的傅立叶变换结果的乘积，为x1(t)和x2(t)的频域滤波器。

(2)提取出函数中每一时刻最大的m个峰值对应的m个时延。

(3)将每一个时刻的m个时延看做到达麦克风时延tdoa的m个量测，并引入tdoa的动态方程：

xk+1＝axk+vk,

zk+1,i＝cxk+wk,i,

其中xk是tdoa在k时刻的状态，vk是均值为0，协方差矩阵为q的高斯白噪声。zk+1,i是在k+1时刻的第i个时延的量测，wk,i是均值为0，协方差矩阵为r的高斯白噪声，且

c＝[10],

(4)利用卡尔曼滤波的预测部分对tdoa状态进行预测。

(5)利用概率数据关联算法计算各个时延是tdoa的概率，然后引入卡尔曼滤波的更新部分对tdoa状态进行更新。

进一步地，步骤(1)中，采用相位变换(phat)加权函数

进一步地，步骤(4)中，引入卡尔曼滤波的预测部分对tdoa状态进行预测：

xk+1|k＝axk|k,

pk+1|k＝apk|ka^t+q,

其中xk|k代表tdoa的状态估计，pk|k代表tdoa的状态估计误差的协方差矩阵，xk+1|k代表一步预测tdoa状态估计，pk+1|k代表一步预测tdoa状态估计误差的协方差矩阵，a^t代表a的转置。

进一步地，步骤(5)中，利用如下公式计算各个时延是tdoa的概率：

其中βⁱk代表在k时刻第i个时延是tdoa的概率，β⁰k代表在k时刻所有时延都不是tdoa的概率；

bk＝λ(2π)^1/2|cpk|k-1c^t+r^-1|^1/2,

其中λ是一个常数，表示第i个时延量测的残差，

进一步地，步骤(5)中，引入卡尔曼滤波的更新部分对tdoa状态进行更新：

kk＝pk|k-1c^t(cpk|k-1c^t+r^-1),

pk|k＝β⁰kpk|k-1+(1-β⁰k)p^ck|k+p^dk,

其中

p^ck|k＝[i-kkc]pk|k-1,

此时tdoa估计值为cxk|k。

本发明提出的基于概率数据关联的移动声源到达时延估计方法，可以估计移动中的无人机声音信号到达麦克风时延，具有成本低廉，计算迅速等特点。本发明具有以下优势：

(1)在不增加成本的条件下，直接利用同一时刻互功率谱密度分布函数提取的多个时延进行融合，提高了tdoa估计的精度。

(2)能工作在低信噪比受多径效应影响严重的环境下，拓宽了声音阵列的工作条件；

(3)本发明方法简便易行，实时性好，易实现。

附图说明

图1为到达麦克风时延图解；

图2为多时刻的最大互功率谱密度分布函数的峰值对应时延的连线；

图3为多时刻的第二大互功率谱密度分布函数的峰值对应时延的连线；

图4为将多个时延融合后的tdoa连续估计结果；

图5为声源或者无人机定位所采用的四面体阵列模型；

图6为无人机飞行定位结果。

具体实施方式

以下结合附图和具体实施例对本发明作进一步详细说明。

本发明提出的一种基于概率数据关联的移动声源到达时延估计方法，包括如下步骤：

(1)根据广义互相关函数计算麦克风阵列t时刻采集的两组时域声信号x1(t)、x2(t)互功率谱密度分布函数

其中为x1(t)和x2(t)的傅立叶变换结果的乘积，为x1(t)和x2(t)的频域滤波器。

可以采用相位变换(phat)加权函数

(2)提取出函数中每一时刻最大的m个峰值对应的m个时延。

(3)将每一个时刻的m个时延看做到达麦克风时延tdoa的m个量测，并引入tdoa的动态方程：

xk+1＝axk+vk,

zk+1,i＝cxk+wk,i,

其中xk是tdoa在k时刻的状态，vk是均值为0，协方差矩阵为q的高斯白噪声。zk+1,i是在k+1时刻的第i个时延的量测，wk,i是均值为0，协方差矩阵为r的高斯白噪声，且

c＝[10],

(4)利用卡尔曼滤波的预测部分对tdoa状态进行预测，具体如下：

xk+1|k＝axk|k,

pk+1|k＝apk|ka^t+q,

其中xk|k代表tdoa的状态估计，pk|k代表tdoa的状态估计误差的协方差矩阵，xk+1|k代表一步预测tdoa状态估计，pk+1|k代表一步预测tdoa状态估计误差的协方差矩阵，a^t代表a的转置。

(5)利用概率数据关联算法计算各个时延是tdoa的概率，然后引入卡尔曼滤波的更新部分对tdoa状态进行更新。

计算各个时延是tdoa的概率的公式如下：

其中βⁱk代表在k时刻第i个时延是tdoa的概率，β⁰k代表在k时刻所有时延都不是tdoa的概率；

bk＝λ(2π)^1/2|cpk|k-1c^t+r^-1|^1/2,

其中λ是一个常数，表示第i个时延量测的残差，

引入卡尔曼滤波的更新部分对tdoa状态进行更新的公式如下：

kk＝pk|k-1c^t(cpk|k-1c^t+r^-1),

pk|k＝β⁰kpk|k-1+(1-β⁰k)p^ck|k+p^dk,

其中

p^ck|k＝[i-kkc]pk|k-1,

此时tdoa估计值为cxk|k。

附图1是一个麦克风到达时延的一个图解，在实际环境中，麦克风与无人机或者移动声源的距离相对于麦克风阵列阵元之间的间距要大很多，所以声源信号传播到麦克风阵列不同阵元可以看做平行入射。从图中可以明显的看到声源信号到达第一个麦克风比到达第二个麦克风要多传播的距离l’，利用这个距离除以声速就可以得到到达时延(tdoa)。

附图2是连续采集几分钟数据后，将每秒的数据利用广义互相关算法求出互功率谱密度分布函数，然后直接将互功率谱密度分布函数上最大的峰值对应的时延当作声音信号到达麦克风时延(tdoa)，再将连续估计的tdoa画成连线图。从图中可以看出，如果直接取最大的峰值对应的时延当作tdoa，得到tdoa估计曲线会有比较大的波动点，这些波动点是由噪声或者多径效应产生的，

附图3是连续采集几分钟数据后，将每秒的数据利用广义互相关算法求出互功率谱密度分布函数，然后直接将互功率谱密度分布函数上第二大的峰值对应的时延当作声音信号到达麦克风时延(tdoa)，再将连续估计的tdoa画成连线图。从图中可以看到，直接利用互功率谱密度分布函数第二大峰值对应时延画出来的连线的波动点相较于附图2是不同的，在附图2中是波动点的地方，附图3中却可能是平滑曲线。这说明由于噪声等干扰，tdoa不一定对应互功率谱密度分布函数的最大峰值，可能第二大，第三大，甚至第四大峰值对应的时延才是真实的tdoa。因此，我们可以引入概率数据关联(pda)算法来计算前几个最大峰值对应时延点是真实tdoa的概率，再利用上述概率将时延加权起来，这样可以提高tdoa估计的准确性。

附图4展示了用本发明提出的时延估计算法进行tdoa估计后得到的tdoa估计连线，可以看到，相较于附图2与附图3，附图4的曲线没有任何的波动，曲线非常平滑，消除了前面由噪声和多径效应产生的干扰点。

附图5是声源或者无人机定位所采用的四面体阵列模型，图中l是阵列原点到麦克风的距离，r1是1号麦克风到达无人机的距离，s1,s2,s3,s4分别是1，2，3，4号麦克风，令(x0,y0,z0)代表无人机的空间坐标，d1i代表声源信号到达1号麦克风与到达i号麦克风的路程差，然后引入下列方程：

x0²+y0²+(z0-l)²＝r1²,

x0²+(y0-l)²+z0²＝(r1+b12)²,

通过上述方程，四个方程四个未知量(x0,y0,z0,r1)，可以定位无人机。

附图6展示了我们利用上述阵列进行定位后的结果，并与无人机飞行的gps轨迹进行了对比，可以看到，在大部分区域都能较好的跟踪无人机。

上述实施例用来解释说明本发明，而不是对本发明进行限制，在本发明的精神和权利要求的保护范围内，对本发明作出的任何修改和改变，都落入本发明的保护范围。