CN106231553B - 基于无线声传感器网络的多节点信息融合声源定位方法 - Google Patents

Abstract

本发明公开了一种基于无线声传感器网络的多节点信息融合声源定位方法，用于解决现有声源定位方法精确度低的技术问题。技术方案是通过卡尔曼滤波器对声源的近场定位结果和远场DOA估计值进行融合，迭代过程中，依据误差协方差矩阵不断校正预测结果，至获取最终声源定位结果，迭代优化过程可以提高定位精度；每次迭代融合远场节点DOA估计值时，是从节点集合中动态选取定位节点，因此可降低单节点估计误差对定位的影响，提高抗噪性。由于采用卡尔曼滤波器方法，加入声源的近场定位结果并与声源的远场DOA估计值进行融合，迭代优化过程有利于提高声源定位精度，使定位精度提高1％～3％。

Description

基于无线声传感器网络的多节点信息融合声源定位方法

技术领域

本发明涉及一种声源定位方法，特别涉及一种基于无线声传感器网络的多节点信息融合声源定位方法。

背景技术

基于无线声传感器网络(Wireless Acoustic Sensor Network,WASN)的声源定位方法，在通信带宽、能量等资源受限条件下，具有资源占用少，简单易实现等显著优势。

文献“Localizing multiple audio sources in a wireless acoustic sensornetwork[J].Signal Processing,2015,107:54-67”公开了一种基于波达方向角(Direction-of-Arrival，DOA)估计值融合的迭代格网声源定位方法。此方法将网络覆盖范围划分成以麦克风阵列节点为顶点的正方形区域，针对声源所在区域迭代地进行网格划分，通过在网格内搜索，从中寻找到各节点DOA估计值误差最小的网格点位置坐标，该坐标作为声源定位结果。文献所述方法是基于声源所在正方形区域的四个节点DOA估计值融合进行声源定位，未考虑当节点距离声源较近时，只能获得声源的近场定位结果而非DOA估计值的实际情况，参与定位信息不全面，定位结果精确度低；此外，若某一节点DOA估计值存在较大偏差，会使定位结果严重偏离真实位置，甚至导致定位失败，因此单节点DOA估计值对定位结果影响大，抗噪性能差。

发明内容

为了克服现有声源定位方法精确度低的不足，本发明提供一种基于无线声传感器网络的多节点信息融合声源定位方法。该方法通过卡尔曼滤波器(Kalman Filter，KF)对声源的近场定位结果和远场DOA估计值进行融合，迭代过程中，依据误差协方差矩阵不断校正预测结果，至获取最终声源定位结果，迭代优化过程可以提高定位精度；每次迭代融合远场节点DOA估计值时，是从节点集合中动态选取定位节点，因此可降低单节点估计误差对定位的影响，提高抗噪性。由于采用卡尔曼滤波器方法，加入声源的近场定位结果并与声源的远场DOA估计值进行融合，迭代优化过程有利于提高声源定位精度，使定位精度提高1％～3％；利用动态选择节点代替固定节点参与DOA估计值融合的远场初步定位结果计算，可降低单节点估计误差影响，抗噪性能增强。

本发明解决其技术问题所采用的技术方案：一种基于无线声传感器网络的多节点信息融合声源定位方法，其特点是包括以下步骤：

步骤一、由m个麦克风阵列节点组成观测区域，融合中心接收目标声源的各节点传送数据，来自第k个节点的数据包括：声源的远场DOA估计值θ_k或近场初步定位结果S_N(k)＝(θ_k，r_k)、节点位置信息P_k＝(x_k，y_k)及声源频谱的特征信息，其中k＝1,2,…m。

步骤二、根据接收数据中距离信息r_k将其分为两类：远场节点数据集合P和近场节点数据集合Q。P,Q集合大小分别为p,q，则有p个节点获得声源的远场DOA估计值(θ₁，θ₂，…θ_p)，q个节点获得声源的近场定位结果((θ₁，r₁)，(θ₂，r₂)，…(θ_q，r_q))，满足p<m,q<m且p+q＝m。

步骤三、在近场节点数据集合Q中，第i个节点获得声源的近场定位结果为S_N(i)＝(θ_i，r_i)，对其进行坐标变换得S_N(i)＝(x_i，y_i)(i＝1,2…q)，计算声源的近场初步定位结果S_N＝(x_N,y_N)，并将其作为KF的初始状态输入值。S_N的计算公式为：

步骤四、针对远场节点数据集合P，计算声源的远场初步定位结果S_F＝(x_F,y_F)。具体步骤为：

(a)从P中选取四个节点组成远场定位节点集合P_select＝{m_j}(j＝1,2,3,4)，节点位置坐标为P_j＝(x_j,y_j)，节点DOA估计值θ_j，并且有其中集合P_select的选取需满足两个条件：

条件1：四个节点围成凸四边形，保证声源的近场初步定位结果在该四边形凸包内。

条件2：四个节点中任意两个节点两两组合，每组节点的最小角度距离函数A(X,Y)满足γ阈值条件。函数A(X,Y)∈(0,π)为X和Y之间的最小角度距离：

其中，X,Y∈(0，2π)分别为每组节点的DOA估计值θ_j1，θ_j2(j1≠j2)。设定平行阈值γ(0<γ<π)，则A(θ_j1，θ_j2)需满足：

γ<A(θ_j1，θ_j2)<π-γ (3)

若上式不成立则任意剔除其中一个节点，重新选择一个新节点加入到集合P_select中，重复上述验证过程，直至满足要求获得P_select。

(b)利用线性最小二乘法融合P_select中节点的DOA估计值，计算公式为：

x_Fsinθ-y_Fcosθ＝x_jsinθ-y_jcosθ(j＝1,2,3,4) (4)

上式的矩阵表示形式记作AS_F ^T＝b，其中 A,b分别为系数矩阵和常量矩阵，S_F ^T为S_F的转置，计算求得声源的远场初步定位结果：

S_F ^T＝(A^TA)^-1A^Tb (5)

(c)输出声源的远场初步定位结果S_F＝(x_F,y_F)，并将其作为KF的第k次观测值赋值给Z(k)＝(x_k,y_k)，则Z(k)＝S_F。

步骤五、利用KF融合声源的远、近场初步定位结果S_F和S_N。在不考虑控制作用的情况下，线性离散***的状态方程和观测方程为：

其中，X是状态变量，Φ是状态转移矩阵，Z是观测变量，H是观测矩阵；W,V为***过程噪声和观测噪声，其协方差矩阵分别为Q(k)和R(k)； u_k(x),u_k(y)分别为x,y方向的观测噪声，两分量相互独立，服从零均值，方差为的高斯白噪声分布；Q(k),W(k)的定义与观测噪声一致。

状态预测方程：

X(k|k-1)＝ΦX(k-1|k-1) (7)

P(k|k-1)＝ΦP(k-1|k-1)Φ^T+Q(k-1) (8)

其中，Φ^T是Φ的转置矩阵，P是X的协方差矩阵。

观测更新方程：

K_g(k)＝P(k|k-1)H^T/(HP(k|k-1)H^T+R) (9)

X(k|k)＝X(k|k-1)+K_g(k)(Z(k)-HX(k|k-1)) (10)

P(k|k)＝(I-K_gH)P(k|k-1) (11)

其中，K_g为卡尔曼增益，I为单位矩阵。

此外，将近场初步定位结果作为参考声源位置S_R＝(x_R,y_R)，以此计算对应的协方差矩阵P，步骤三中声源的近场初步定位结果S_N的协方差矩阵为P_N，步骤四中声源的远场初步定位结果S_F的协方差矩阵为P_F，P的计算公式为：

其中，噪声方差分别为

KF融合过程中参数对应关系：初始状态输入值X₀＝S_N，P₀＝P_N；迭代过程中，第k次的观测值Z(k)＝S_F，第k次的预测值X(k|k-1)，对应的估计值X(k|k)；迭代终止，最终声源定位结果记为S^*＝(x^*,y^*)。

具体实现流程：

①输入初始状态值X₀＝S_N，P₀＝P_N；

②利用公式(7)、公式(8)预测第k次迭代的声源位置，预测结果为X(k|k-1)＝(x_k,y_k)，对应的协方差矩阵为P(k|k-1)；

③利用步骤四获得声源的远场初步定位结果S_F＝(x_F,y_F)，由公式(12)计算对应的协方差矩阵P_F；

④输入第k次的观测值Z(k)＝S_F，R(k)＝P_F，利用公式(9)、公式(10)校正预测结果，获得第k次的最优估计结果X(k|k)＝(x_k,y_k)，更新公式(11)获得协方差矩阵P(k|k)；

⑤对第k次的估计结果，判断是否满足终止条件。若满足则终止迭代，输出估计结果S^*＝X(k|k)；否则令k＝k+1继续迭代，返回步骤②。

本发明的有益效果是：该方法通过卡尔曼滤波器(Kalman Filter，KF)对声源的近场定位结果和远场DOA估计值进行融合，迭代过程中，依据误差协方差矩阵不断校正预测结果，至获取最终声源定位结果，迭代优化过程可以提高定位精度；每次迭代融合远场节点DOA估计值时，是从节点集合中动态选取定位节点，因此可降低单节点估计误差对定位的影响，提高抗噪性。由于采用卡尔曼滤波器方法，加入声源的近场定位结果并与声源的远场DOA估计值进行融合，迭代优化过程有利于提高声源定位精度，使定位精度提高1％～3％；利用动态选择节点代替固定节点参与DOA估计值融合的远场初步定位结果计算，可降低单节点估计误差影响，抗噪性能增强。

下面结合附图和具体实施方式对本发明作详细说明。

附图说明

图1是本发明基于无线声传感器网络的多节点信息融合声源定位方法的流程图。

图2是图1中选择远场定位节点集合的子流程图。

图3是本发明方法实施例8m×8m监测区域内节点和声源位置分布示意图。

图4是背景技术方法与本发明方法进行声源定位的定位均方根误差与参考信噪比变化关系的仿真图。

具体实施方式

参照图1-4。本发明基于无线声传感器网络的多节点信息融合声源定位方法具体步骤如下：

三角形表示目标声源真实位置S＝(x,y)，本实施例设为S＝(2.6,3.0)。黑色圆点表示各麦克风阵列节点，对应编号为m_k(k＝1,2,…9)并组成观测节点集合M＝{m_k}(k＝1,2…9)。声源的观测区域由四个小正方形区域组成，边长V＝4。节点m_k的位置坐标P_k＝(x_k,y_k)对应正方形各顶点的坐标，如节点m₁和m₅的位置坐标分别为P₁＝(x₁,y₁)＝(0,0)，P₅＝(x₅,y₅)＝(4,4)。本实施例中，以观测区域中心位置m₅节点处的信噪比作为参考信噪比计算其他各节点信噪比，参考信噪比变化范围为0～20dB。节点处DOA观测噪声η_θ服从零均值，协方差为的高斯分布根据节点信噪比SNR与DOA估计误差关系式：δ_θ(SNR)＝1.979e^-0.2875(SNR)+1.884，可产生各节点m_k(k＝1,2,…9)的观测数据。

利用卡尔曼滤波器进行声源定位，具体过程为：

1、接收各节点数据并进行分类。

融合中心接收各节点m_k(k＝1,2,…9)传送数据：声源的远场DOA估计值θ_k或近场定位结果S_N(k)＝(θ_k，r_k)以及节点位置坐标P_k＝(x_k，y_k)；并依据距离信息r_k将接收数据分为两类：远场节点数据集合P＝{m_j}(j＝1,2…p)和近场节点数据集合Q＝{m_i}(i＝1,2…q)，p,q分别为远、近场节点数据集合P,Q大小，本实施例中q＝2,p＝7，则有两个近场节点m₄和m₅，其他节点m_j(j＝1,2,3,6,7,8,9)为远场节点。

2、计算声源的近场初步定位结果S_N＝(x_N,y_N)。

对Q中节点的近场定位结果S_N(i)＝(θ_i，r_i)进行坐标变换得S_N(i)＝(x_i，y_i)(i＝1,2…q)，利用公式(1)计算S_N＝(x_N,y_N)，并赋值给KF的初始状态值X₀，有X₀＝S_N。

以近场初步定位结果S_N作为参考声源位置S_R＝(x_R,y_R)。根据S_R分别计算x,y方向的方差并由公式(2)计算S_N对应的协方差矩阵P_N，将其赋值给KF的初始协方差矩阵P₀，有P₀＝P_N。

3、计算声源的远场初步定位结果S_F＝(x_F,y_F)。

从远场节点数据集合P中动态选取远场定位节点集合P_select，具体步骤为：

①从P中随机选取4个节点作为P_select＝{m_j}(j＝1,2,3,4)，节点m_j的DOA估计值为θ_j，位置信息为P_j＝(x_j,y_j)，并记录对应的剩余节点集合其中m_j∈P_select，且满足转至步骤④；

②判断是否为空，若则从当前组合中任意剔除一个节点m_j1或m_j2，转至步骤③；否则返回步骤①；

③从中随机选取一个新节点，添加到P_select并更新转至步骤④，

④对P_select使其能够围成凸四边形，并验证声源的近场初步定位结果S_N是否在该四边形的凸包内，若满足条件，则转至步骤⑤，否则返回步骤①；

⑤P_select中任意2个节点m_j1,m_j2两两组合，对每组节点根据其DOA估计值θ_j1,θ_j2，计算最小角度距离函数A(θ_j1,θ_j2)(j1≠j2)；

⑥根据平行阈值γ判断A(θ_j1,θ_j2)是否满足阈值条件：γ<A(θ_j1,θ_j2)<π-γ，本实施例中γ设为20°。若满足条件，则转至步骤⑦，否则转至步骤②；

⑦节点选择成功，输出选择结果P_select。

之后利用输出的P_select，通过公式(3)计算声源的远场初步定位结果S_F＝(x_F,y_F)。

表示为矩阵形式：AS_F ^T＝b

其中，A,b分别为系数矩阵和常量矩阵，S_F ^T为S_F的转置。公式(4)的计算结果为：

S_F ^T＝(A^TA)^-1A^Tb (4)

根据S_R分别计算x,y方向的观测噪声方差由公式(5)计算S_F对应的协方差矩阵P_F，将其赋值给KF的观测协方差矩阵R，有R＝P_F。

4、利用卡尔曼滤波器进行信息融合，获取声源定位结果，具体过程为：

①利用公式(1)(2)初始化***状态方程(6)，开始迭代。

X(k)＝Φ(k-1)X(k-1)+W(k-1) (6)

②利用公式(7)(8)，预测第k次的状态变量X(k)和协方差矩阵P(k)。

X(k|k-1)＝ΦX(k-1|k-1) (7)

P(k|k-1)＝ΦP(k-1|k-1)Φ^T+Q(k-1) (8)

其中，将声源位置S作为***状态变量X(k)＝(x_k,y_k)，P(k)为X(k)的协方差矩阵，初始化输入值为：X(0)＝S_N，P(0)＝P_N。状态转移矩阵ΦT是Φ的转置矩阵。W为***过程噪声，其协方差矩阵设为恒定值，本实施例中取Q＝diag(10^-4，10^-4)。

③利用公式(4)(5)将声源的远场初步定位结果作为观测值，输入观测方程(9)中。

Z(k)＝H(k)X(k)+V(k) (9)

其中，观测变量Z(k)＝S_F，观测矩阵V为观测噪声，对应的观测误差R由公式(5)计算协方差矩阵P_F获得，有R＝P_F。

④利用公式(10)～(12)，校正预测结果。

K_g(k)＝P(k|k-1)H^T/(HP(k|k-1)H^T+R) (10)

X(k|k)＝X(k|k-1)+K_g(k)(Z(k)-HX(k|k-1)) (11)

P(k|k)＝(I-K_gH)P(k|k-1) (12)

其中，K_g为卡尔曼增益，I为单位矩阵。

⑤设置KF的最大迭代次数N，每次迭代时进行验证，本实施例中设N＝10。若k<N，则重复公式(6)～(12)对应过程；否则终止迭代，输出当前估计结果X(k|k)＝(x_k,y_k)，赋值给S^*＝X(k|k)，即获得最终声源定位结果S^*＝(x^*,y^*)。此外，设置T为蒙特卡洛仿真次数，重复进行多次实验，本实施例中取T＝20。

分别对背景技术的格网定位法GB及本发明的KF方法进行仿真，计算各自仿真定位结果与真实声源位置的定位均方根误差，并求其相对于正方形边长的百分比实现标准化，以此作为定位精度的衡量标准。由图4仿真结果可以看出，通过融合声源的近场定位结果和远场DOA估计值，本发明方法能够降低声源定位误差，在较高信噪比时，定位精度可提高1％，较低信噪比时，提高1.5％～3％；迭代优化的定位结果抗噪性能增强，特别是在较低信噪比的条件下，对于改善定位结果，降低声源定位误差具有显著效果。

Claims (1)

1.一种基于无线声传感器网络的多节点信息融合声源定位方法，其特征在于包括以下步骤：

步骤一、由m个麦克风阵列节点组成观测区域，融合中心接收目标声源的各节点传送数据，来自第k个节点的数据包括：声源的远场DOA估计值θ_k或近场初步定位结果S_N(k)＝(θ_k，r_k)、节点位置信息P_k＝(x_k，y_k)及声源频谱的特征信息，其中k＝1,2,…m；

步骤二、根据接收数据中距离信息r_k将其分为两类：远场节点数据集合P和近场节点数据集合Q；P,Q集合大小分别为p,q，则有p个节点获得声源的远场DOA估计值(θ₁，θ₂，…θ_p)，q个节点获得声源的近场定位结果((θ₁，r₁)，(θ₂，r₂)，…(θ_q，r_q))，满足p<m,q<m且p+q＝m；

步骤三、在近场节点数据集合Q中，第i个节点获得声源的近场定位结果为S_N(i)＝(θ_i，r_i)，对其进行坐标变换得S_N(i)＝(x_i，y_i)(i＝1,2…q)，计算声源的近场初步定位结果S_N＝(x_N,y_N)，并将其作为卡尔曼滤波器KF的初始状态输入值；S_N的计算公式为：

步骤四、针对远场节点数据集合P，计算声源的远场初步定位结果S_F＝(x_F,y_F)；具体步骤为：

(a)从P中选取四个节点组成远场定位节点集合P_select＝{m_j}(j＝1,2,3,4)，节点位置坐标为P_j＝(x_j,y_j)，节点DOA估计值θ_j，并且有其中集合P_select的选取需满足两个条件：

条件1：四个节点围成凸四边形，保证声源的近场初步定位结果在该四边形凸包内；

条件2：四个节点中任意两个节点两两组合，每组节点的最小角度距离函数A(X,Y)满足γ阈值条件；函数A(X,Y)∈(0,π)为X和Y之间的最小角度距离：

其中，X,Y∈(0，2π)分别为每组节点的DOA估计值θ_j1，θ_j2(j1≠j2)；设定平行阈值γ(0<γ<π)，则A(θ_j1，θ_j2)需满足：

γ<A(θ_j1，θ_j2)<π-γ (3)

若上式不成立则任意剔除其中一个节点，重新选择一个新节点加入到集合P_select中，重复上述验证过程，直至满足要求获得P_select；

(b)利用线性最小二乘法融合P_select中节点的DOA估计值，计算公式为：

x_Fsinθ-y_Fcosθ＝x_jsinθ-y_jcosθ(j＝1,2,3,4) (4)

上式的矩阵表示形式记作AS_F ^T＝b，其中 A,b分别为系数矩阵和常量矩阵，S_F ^T为S_F的转置，计算求得声源的远场初步定位结果：

S_F ^T＝(A^TA)^-1A^Tb (5)

(c)输出声源的远场初步定位结果S_F＝(x_F,y_F)，并将其作为卡尔曼滤波器KF的第k次观测值

赋值给Z(k)＝(x_k,y_k)，则Z(k)＝S_F；

步骤五、利用卡尔曼滤波器KF融合声源的远、近场初步定位结果S_F和S_N；在不考虑控制作用的情况下，线性离散***的状态方程和观测方程为：

其中，X是状态变量，Φ是状态转移矩阵，Z是观测变量，H是观测矩阵；W,V为***过程噪声和观测噪声，其协方差矩阵分别为Q(k)和R(k)； u_k(x),u_k(y)分别为x,y方向的观测噪声，两分量相互独立，服从零均值，方差为的高斯白噪声分布；Q(k),W(k)的定义与观测噪声一致；

状态预测方程：

X(k|k-1)＝ΦX(k-1|k-1) (7)

P(k|k-1)＝ΦP(k-1|k-1)Φ^T+Q(k-1) (8)

其中，Φ^T是Φ的转置矩阵，P是X的协方差矩阵；

观测更新方程：

K_g(k)＝P(k|k-1)H^T/(HP(k|k-1)H^T+R) (9)

X(k|k)＝X(k|k-1)+K_g(k)(Z(k)-HX(k|k-1)) (10)

P(k|k)＝(I-K_gH)P(k|k-1) (11)

其中，K_g为卡尔曼增益，I为单位矩阵；

此外，将近场初步定位结果作为参考声源位置S_R＝(x_R,y_R)，以此计算对应的协方差矩阵P，步骤三中声源的近场初步定位结果S_N的协方差矩阵为P_N，步骤四中声源的远场初步定位结果S_F的协方差矩阵为P_F，P的计算公式为：

其中，噪声方差分别为

卡尔曼滤波器KF融合过程中参数对应关系：初始状态输入值X₀＝S_N，P₀＝P_N；迭代过程中，第k次的观测值Z(k)＝S_F，第k次的预测值X(k|k-1)，对应的估计值X(k|k)；迭代终止，最终声源定位结果记为S^*＝(x^*,y^*)；

具体实现流程：

①输入初始状态值X₀＝S_N，P₀＝P_N；

②利用公式(7)、公式(8)预测第k次迭代的声源位置，预测结果为X(k|k-1)＝(x_k,y_k)，对应的协方差矩阵为P(k|k-1)；

③利用步骤四获得声源的远场初步定位结果S_F＝(x_F,y_F)，由公式(12)计算对应的协方差矩阵P_F；

④输入第k次的观测值Z(k)＝S_F，R(k)＝P_F，利用公式(9)、公式(10)校正预测结果，获得第k次的最优估计结果X(k|k)＝(x_k,y_k)，更新公式(11)获得协方差矩阵P(k|k)；

⑤对第k次的估计结果，判断是否满足终止条件；若满足则终止迭代，输出估计结果S^*＝X(k|k)；否则令k＝k+1继续迭代，返回步骤②。