Acoustic Source Tracking Based on Probabilistic Data Association and Distributed Cubature Kalman Filtering in Acoustic Sensor Networks

Abstract

A probabilistic data association-based distributed cubature Kalman filter (PDA-DCKF) method is proposed in this paper, whose performance on tracking single moving sound sources in the distributed acoustic sensor network was verified. In this method, the PDA algorithm is first used to sift the observations from neighboring nodes. Then, the sifted observations are fused to update the state vectors in the CKF. Since nodes in a sensor network have different reliabilities, the final tracking result integrates the estimations from the local nodes, which are weighted with the parameters depending on the mean square error of the estimation and the energy of the received signal. The experimental results illustrated that the proposed PDA-DCKF method is superior to the other DCKF methods in tracking sound sources even under severe noise and reverberant conditions.

1. Introduction

The problem of acoustic source localization and tracking has always been one of the research hotspots in the field of speech processing. It has been widely used in many aspects, such as audio and video conferencing systems, human-computer interaction and speech enhancement, etc. [1,2,3,4]. Traditional acoustic localization and tracking methods usually require the microphone array to have a regular geometric structure, and generally use a centralized data processing method [5]. With the continuous advancement of technology, some traditional microphone arrays gradually show some deficiencies. The distributed microphone network has attracted more and more research work because it has no strict restrictions on the arrangement of microphones, and is a network composed of multiple nodes arbitrarily distributed in space, usually each node contains a set of microphones [6,7,8,9,10].

So far, there have been many studies on acoustic source localization using distributed microphone networks [11]. But they only locate the acoustic source based on the current observations of multiple microphones, which can locate the acoustic source when the background noise and reverberation are small. In noisy and reverberant environments, spurious observations may even mask observations from real acoustic sources, degrading localization performance. To avoid this problem, a Bayesian filter [12] combined the current observation with a series of past observations for current position estimation, which is more effective for dealing with the adverse effects of noise or reverberation. Theoretically, Bayesian filters describe the tracking problem with a state-space model that includes a dynamic model that describes the motion of the target and an observation model that describes the relationship between the observations and the state of the acoustic source. When the state space model is linear and Gaussian, Kalman filter can replace Bayesian filter. However, in acoustic source tracking scenarios, the observation function is usually nonlinear, and some conditions and properties applicable to linear systems no longer hold, and the performance of the Kalman filter may be severely degraded.

Extended Kalman filter (EKF) was first proposed in [13] and is the simplest and most widely used nonlinear filtering method. EKF only intercepts the first-order term in the linearization process to approximate the system function, but its error is relatively large. In order to solve this problem, the literature [14] proposed iterative extended Kalman filter (IEKF), which can improve the accuracy of EKF through several iterations. However, when the nonlinearity of the system model is very strong, whether it is EKF or IEKF, their effect is not good, and there are disadvantages such as poor stability and easy divergence. The particle filter (PF) is a Monte Carlo implementation of a Bayesian filter that approximates the state by a series of weighted particles extracted from the proposal function [15]. PF can handle nonlinear and non-Gaussian situations well, and many PF-based sound source tracking methods have been developed. Vermaak and Blake [16] first introduced PF to sound source tracking. The particle filter method breaks through the limitations of linearity and Gaussian, but its computational load is too large. In practical applications, particle filtering will only be considered when the approximate filter fails. The Sigama point Kalman filter method [17] was similar to the idea of particle filter. It does not use the method of approximating nonlinear systems, but directly uses the real system and observation model by selecting a set of effective deterministic sampling sets, namely Sigama points, which can achieve second-order accuracy. According to the different selection methods of Sigama points, it can be divided into UKF, CKF, QKF and so on. In general, several Gaussian filtering methods introduced above are centralized, that is, the data of all nodes is collected and transmitted to the central processing unit to perform the task of acoustic source tracking. This method is generally unreliable, as any failure of the central processor renders the entire network untraceable.

In order to solve the unreliable problem of centralized methods, many distributed methods have been developed for sound source tracking. No central processor is needed in the distributed method, and all nodes realize the estimation of the global state only by exchanging data with their neighbors. In reference [18], a distributed extended Kalman particle filter (DEKPF) for speaker tracking was developed, which combined the current TDOA observations into EKF to propose particle filter. In reference [19], a distributed particle filter (DPF) was proposed, which applied the improved iterative covariance intersection (MICI) algorithm and interactive multiple model (IMM) to speaker tracking in distributed microphone networks. In reference [20], a distributed iterative EKF was proposed to estimate the time-varying speaker position in the microphone array. In reference [21], a Distributed Unscented Kalman Filter (DUKF) is proposed to overcome the nonlinearity of the measurement model in speaker tracking. The time difference of arrival (TDOA) was used as the observation and then the distributed IMM-UKF was used to track the location of the sound source.

In the actual environment, the existence of noise or reverberation usually produces unreliable observations with false peaks, which may lead to serious performance degradation. Usually, the current observations contained in these methods are only extracted from the largest peak value of a certain observation function. In some bad cases, the peak value related to the real acoustic source may be masked by the stray acoustic source. Therefore, it is more reasonable to extract multiple observations from the observation function, rather than one observation, and then incorporate it into the above tracking scheme. Probabilistic Data Association (PDA) [22] is an effective method to combine multiple observations into Kalman filter state update, which has been proved to be suitable for target tracking in clutter environment. In reference [23], an improved distributed unscented Kalman particle filter (DUKPF) was proposed to track a single moving acoustic source using a distributed microphone network in noise and reverberation environments. This method proposed to extract multiple observations from the observation function of each node and combined them into the status update of UKF through probabilistic data association (PDA) technology, so as to generate PDA-UKF, and then brought in particle filter. In reference [24], a microphone array network distributed multi speaker tracking method based on tasteless particle filter and data association was proposed. The available observations were associated with each speaker at each node using data association technology to track the speaker. Reference [25] proposed a volume information filter based on joint probabilistic data association (JPDA) for multi acoustic source tracking based on distributed acoustic vector sensor (AVS) array, in which JPDA was used to deal with the correlation between observations and targets. Issues related to multi-source tracking are beyond the scope of this article. However, most of particle filter-based methods require excessive computational costs, which limits them in real-time applications. Besides, in existing speaker tracking methods, the PDA algorithm is applied to sift the observations without considering the information from neighboring nodes.

Probabilistic data association with cubature Kalman filtering are combined in this paper, and they are applied to the problem of single-acoustic source tracking in noisy and reverberant environments with distributed acoustic sensor networks. The contributions of this paper are as follows:

• Combining the cubature Kalman filter (CKF) with PDA, the probabilistic data association-cubature Kalman filter (PDA-CKF) was developed. In PDA-CKF, multiple possible observations were merged into the state update of CKF by the PDA technique.
• In this paper, PDA-CKF was applied to the distributed acoustic sensor network, and the probabilistic data association-distributed cubature Kalman filter (PDA-DCKF) was developed by combining the observation information of each node’s neighbor nodes in the network.
• Considering the reliability of the local state, it was proposed to combine the mean square error (MSE) of the position estimation of each node and the received signal energy to adjust the weighting coefficient of distributed acoustic sensor data fusion. In this way, the local state of high-quality nodes is enhanced, and each node can achieve global consistency and good speaker tracking performance.

The structure of this paper is as follows. Section 2 presents the problem formulation, background knowledge, and some prior knowledge of acoustic source tracking. Section 3 first introduces the single-node PDA-CKF and then details the distributed PDA-DCKF. Section 4 presents the experimental results and discussion. Section 5 summarizes some conclusions.

2. Background Knowledge

2.1. Problem Formulation

Consider a distributed sensor network with $N$ nodes deployed as shown in Figure 1. The positions of the nodes can be obtained in advance by calibration [26]. Each node in the DMA consists of two microphones at distance $L$. All nodes of the network are modeled as vertices of the graph $\mathcal{G}1=(\epsilon ,\upsilon )$, where $\upsilon =\{1,2,\dots ,N\}$ is the vertex set, $\epsilon \subset \{(p,q)|p,q\in \upsilon \}$ is the edge set, and $(p,q)\in \epsilon$ represents the network’s communication constraints, i.e., node $p$ can send information to node $q$, and vice versa. Let ${\mathcal{N}}_{p,k}=\{q\in \upsilon |(p,q)\in \epsilon \}\cup \{p\}$ denote the set of neighbors of node $p$ at time $k$, where a node is a neighbor of itself certainly.

2.2. Signal Model and TDOA Estimation

In acoustic sensor networks, the discrete-time signal acquired by the $l-th$ microphone ($l=1,2$) of node $p$ can be modeled as [23]

$$y_{p,l}(t)=h_{p,l}(t)\ast y(t)+e_{p,l}(t),\forall p\in \upsilon$$

(1)

where $t$ is the discrete-time index, $h_{p,l}(t)$ is the room impulse response (RIR) between the microphone and the acoustic source, $\ast$ denotes the convolution operator, $y(t)$ is the source signal, and $e_{p,l}(t)$ is the additive noise.

Traditionally, the generalized cross-correlation function (GCC) [27] is used for TDOA estimation. Assuming that $Y_1(k)$ and $Y_2(k)$ are the acoustic signal received by a microphone pair at time $k$ and $Y_l(f)=FFT\{Y_l(k)\},l=1,2$ is the frequency domain representation of the corresponding acoustic signal in a time frame, the generalized cross-correlation function of the acoustic signal received by the microphone pair is

$$R_{12}(\tau )={\displaystyle {\int }_{-\infty }^{+\infty }\frac{Y_1(f)Y_2^{*}(f)}{|Y_1(f)Y_2^{*}(f)|}{e}^{j2\pi f\tau }df}$$

(2)

where $Y_1(f)$ and $Y_2(f)$ represent the frequency-domain microphone signals at the node, and $\ast$ represents the complex conjugation operation. Therefore, the delay estimation is [27]

$$\widehat{\tau }={\displaystyle {\int }_{-\infty }^{+\infty }\underset{\tau \in [-\tau \max ,\tau \max ]}{\arg \max }{R}_{12}(\tau )}$$

(3)

$\tau \max$ is the largest time delay estimation.

However, in the real indoor environment, reverberation and noise will bring false maxima of $R_{12}(\tau )$ and obtain invalid TDOA estimation. In order to solve this problem, the local largest of the first $Q$ largest peaks of $R_{12}(\tau )$ are taken as the candidate measurement value of multiple TDOA of node $p$ at time $k$. In this paper, multiple TDOA observations were extracted through a two-step selection process, taking node $p$ as an example [23].

(1)

Select $Q$ delays according to the peak amplitude of the GCC, i.e.,

$$z_{p,k}=[{\tilde{\tau }}_{p,k}^{(1)},{\tilde{\tau }}_{p,k}^{(2)},\dots ,{\tilde{\tau }}_{p,k}^{(Q)}]$$

(4)

where ${\tilde{\tau }}_{p,k}^{(i)}$ is the delay of node $p$ related to the $i-\text{th}$ largest peak of $R_{12}(\tau )$ at time $k$.

(2)

Further, select $m_{p,k}$ observations from (4) as local observations, and the selection rules are shown in Section 3.

$${\stackrel{⌣}{z}}_{p,k}=[{\tilde{\tau }}_{p,k}^{(1)},{\tilde{\tau }}_{p,k}^{(2)},\dots ,{\tilde{\tau }}_{p,k}^{(m_{p,k})}]$$

(5)

where each delay ${\widehat{\tau }}_{p,k}^{(j)},j=1,2,\dots ,m_{p,k}$ is deemed as a TDOA candidate.

2.3. Dynamic Model of Acoustic Source

Without loss of generality, the two-dimensional tracking is considered herein, since the height of a moving acoustic source would usually not change significantly. Speakers move in a room with a distributed acoustic sensor network, and Langevin model [24] can accurately and simply describe the time-varying position of speakers. At time $k$, the state of the speaker is defined as $x_k={[x_k,y_k,{\dot{x}}_k,{\dot{y}}_k]}^T$, where ${(x_k,y_k)}^T$ and ${({\dot{x}}_k,{\dot{y}}_k)}^T$ represent the position and moving speed of the speaker, respectively. In this model, the speaker’s motion in the Cartesian coordinate system is considered to be independent and modeled as [23]

$$x_k=\left[\begin{array}{cc}{I}_2& a\Delta T\otimes I_2\\ 0& a\otimes I_2\end{array}\right]x_{k-1}+\left[\begin{array}{cc}b\Delta T\otimes I_2& 0\\ 0& b\otimes I_2\end{array}\right]u_{k-1}$$

(6)

where $a=e^{-\beta \Delta T}$, and $b=\overline{\upsilon }\sqrt{1-a^2}$; $\beta$ and $\overline{\upsilon }$ are the rate constant and the steady velocity parameter, respectively. $I_s$ denotes the s-order identical matrix, $\otimes$ stands for the Kronecker product, $\Delta T$ is the sampling period for position estimation, and $u_{k-1}$ is the zero-mean white Gaussian noise with identity covariance matrix, which describes the uncertainty of the acoustic source motion.

2.4. Bayesian Framework for Speaker Tracking

Bayesian filtering is the basis of Kalman filtering. This section briefly reviews the basic principles of the Bayesian filtering algorithm.

Assuming that the state variable at time $k$ is $x_k\in {\mathbb{ℝ}}^p$ and its observation value is $y_k\in {\mathbb{ℝ}}^q$, where ${\mathbb{ℝ}}^n$ represents the n-dimensional real vector space, the state equation and observation equation are expressed as [21]:

$$x_{k+1}=f_k(x_k)+{\Gamma }_k{w}_k$$

(7)

$$y_k=h_k(x_k)+v_k$$

(8)

where $f_k(\cdot )$ is the nonlinear state transfer function, $h_k(\cdot )$ is the nonlinear observation function, ${\Gamma }_k$ is the noise transfer matrix, $w_k$ is the process noise, and $v_k$ is the observation noise, which meets [21]

$$E\left\{\left[\begin{array}{c}{w}_k\\ v_k\end{array}\right]{\left[\begin{array}{c}{w}_l\\ v_l\end{array}\right]}^T\right\}=\left[\begin{array}{cc}{Q}_k{\delta }_{k,l}& 0\\ 0& R_k{\delta }_{k,l}\end{array}\right]$$

(9)

where the superscript $T$ represents the transpose of the matrix, $E\{\cdot \}$ represents the expected operator, and ${\delta }_{k,l}$ represents the Kronecker delta function. $Q_k$ and $R_k$ are the covariance matrices of noise $w_k$ and $v_k$, respectively, and it is assumed that they are both positive definite.

The Bayesian filtering problem is to infer the estimated value of the state variable $x_k$ at time $k$ given the observation information $y_{1:k}=\{y_1,\dots ,y_k\}$ at time $k$, i.e., to estimate the posterior probability density $p(x_k|y_{1:k})$. Assuming that the initial probability density function $p(x_0)$ of the state variable is known as prior knowledge, the posterior probability density $p(x_k|y_{1:k})$ can be obtained recursively by the following equations [20]:

$$p(x_k|y_{1:k-1})={\displaystyle \int p(x_k|x_{k-1})p(x_{k-1}|y_{1:k-1})d{x}_{k-1}}$$

(10)

$$p(x_k|y_{1:k})=\frac{p(y_k|x_k)p(x_k|y_{1:k-1})}{p(y_k|y_{1:k-1})}$$

(11)

In Equations (10) and (11), the state transition probability density function $p(x_k|x_{k-1})$ is defined by the state equation; the observation likelihood probability density function $p(y_k|x_k)$ is defined by the observation equation.

3. Improved Distributed Cubature Kalman Filter

In the CKF, the observation corresponding to the largest peak of the observation function is used for the state update. This approach works well under moderate acoustic environments, while its performance degrades in severe noise and reverberation conditions because the spurious peaks from noise or reverberation may cover up the peaks from real acoustic sources. To alleviate this problem, multiple observations are selected from the multiple local maxima of the observation function. A general framework for state updates that integrates multiple possible observations is provided by the probabilistic data association (PDA). Inspired by this idea, the probabilistic data association-cubature Kalman filter (PDA-CKF) was derived in this paper. Next, PDA-CKF was used for acoustic source tracking in distributed acoustic sensor networks, and an improved PDA-DCKF algorithm was developed. The observations of multiple nodes in the neighborhood are filtered by PDA and then merged into the state update of CKF to integrate the information of multiple nodes to realize distributed tracking.

Before introducing PDA-CKF, the preliminary knowledge of cubature Kalman—cubature point set $\{{\xi }_i,{\omega }_i\}$ [28]—should be introduce first.

The standard Gaussian weighted integral is calculated using the spherical-radial cubature rule, i.e., [28]

$$E[x|z]={\displaystyle \underset{R}{\int }f(x)\mathcal{N}(x;0,P)dx\approx {\displaystyle \sum _{i=1}^{2n}{\omega }_{i}f({\xi }_i)}}$$

(12)

In Equation (12), $f(\cdot )$ is the nonlinear state transfer function or observation function, $n$ is the dimension of the state variable, $\mathcal{N}(x;0,P)$ is a Gaussian distribution function with a mean of zero and a variance of $P$, and ${\xi }_i$ is the cubature points.

$${\xi }_i=\sqrt{n}{\left[1\right]}_i,i=1,2,\dots ,2n$$

(13)

$${\omega }_i=\frac{1}{2n}$$

(14)

${\left[1\right]}_i$ represents the point set of $n$ (n-dimensional state) dimensional space, i.e.,

$${\left[1\right]}_i=\left\{\left(\begin{array}{c}1\\ 0\\ ⋮\\ 0\end{array}\right),\left(\begin{array}{c}0\\ 1\\ ⋮\\ 0\end{array}\right),\dots ,\left(\begin{array}{c}0\\ 0\\ ⋮\\ 1\end{array}\right),\left(\begin{array}{c}-1\\ 0\\ ⋮\\ 0\end{array}\right),\left(\begin{array}{c}0\\ -1\\ ⋮\\ 0\end{array}\right),\dots ,\left(\begin{array}{c}0\\ 0\\ ⋮\\ -1\end{array}\right)\right\}$$

(15)

3.1. PDA-CKF Algorithm

(a)

Initialization

When $k=0$, assuming $x_0~\mathcal{N}({\overline{x}}_0,P_0)$, the initial value of the process noise and observation noise matrix are set to Q₀ and R₀, respectively. Then, the optimal initialization of the filter is

$$\begin{array}{l}{\widehat{x}}_{p,0|0}={\overline{x}}_0\\ {\widehat{P}}_{p,0|0}=P_0\end{array}$$

(16)

(b)

State Prediction

For each node $p$, the state estimate and covariance matrix ${\widehat{x}}_{p,k-1|k-1},{\widehat{P}}_{p,k-1|k-1}$ at time $k-1$ are given, and the positive definite noise matrix $Q_{p,k-1},R_{p,k-1}$ are given. Using Equations (13) and (14), the state predicted cubature points ${\chi }_{p,k-1|k-1}^i$ is calculated as:

$$S_{p,k-1|k-1}=\sqrt{{\widehat{P}}_{p,k-1|k-1}}$$

(17)

$${\chi }_{p,k-1|k-1}^i={\widehat{x}}_{p,k-1|k-1}+S_{p,k-1|k-1}{\xi }_i,i=1,2,\dots ,2n$$

(18)

According to the state transition model, the cubature points are propagated nonlinearly, i.e.,

$${\chi }_{p,k|k-1}^i=f({\chi }_{p,k-1|k-1}^i),i=1,2,\dots ,2n,p=1,2,\dots ,N$$

(19)

where $n$ represents the dimension of the state variable, and $N$ represents the number of nodes in the distributed acoustic sensor network. At this time, the state prediction ${\widehat{x}}_{p,k|k-1}$ and its error matrix ${\widehat{P}}_{p,k|k-1}$ are calculated as:

$${\widehat{x}}_{p,k|k-1}=\frac{1}{2n}{\displaystyle \sum _{i=1}^{2n}{\chi }_{p,k|k-1}^i},i=1,2,\dots ,2n,p=1,2,\dots ,N$$

(20)

$${\widehat{P}}_{p,k|k-1}=\frac{1}{2n}{\displaystyle \sum _{i=1}^{2n}({\chi }_{p,k|k-1}^i-{\widehat{x}}_{p,k|k-1}){({\chi }_{p,k|k-1}^i-{\widehat{x}}_{p,k|k-1})}^T}+Q_{p,k},p=1,2,\dots ,N$$

(21)

(c)

Status Update

From the estimated ${\widehat{x}}_{p,k|k-1}$ and variance ${\widehat{P}}_{p,k|k-1}$ at time $k$, the state update cubature points ${\chi }_{p,k|k-1}^i$ is calculated as:

$$S_{p,k|k-1}=\sqrt{{\widehat{P}}_{p,k|k-1}}$$

(22)

$${\chi }_{p,k|k-1}^i={\widehat{x}}_{p,k|k-1}+S_{p,k|k-1}{\xi }_i$$

(23)

${\chi }_{p,k|k-1}^i$ is propagated through the observation equation,

$${\widehat{Z}}_{p,k|k-1}^i=h({\chi }_{p,k|k-1}^i),i=1,2,\dots ,2n,p=1,2,\dots ,N$$

(24)

Further, the observation prediction ${\widehat{z}}_{p,k|k-1}$ and the observation prediction error variance $P_{p,k|k-1}^{zz}$ are, respectively, obtained by

$${\widehat{z}}_{p,k|k-1}=\frac{1}{2n}{\displaystyle \sum _{i=1}^{2n}{\widehat{Z}}_{p,k|k-1}^i}$$

(25)

$$P_{p,k|k-1}^{zz}=\frac{1}{2n}{\displaystyle \sum _{i=1}^{2n}({\widehat{Z}}_{p,k|k-1}^i-{\widehat{z}}_{p,k|k-1}){({\widehat{Z}}_{p,k|k-1}^i-{\widehat{z}}_{p,k|k-1})}^T+R_{p,k}}$$

(26)

Then, according to the probabilistic data association, the verification area of node $p$ can be constructed by [29]:

$$\{z_{p,k}:{(z_{p,k}-{\widehat{z}}_{p,k|k-1})}^T{(P_{p,k|k-1}^{zz})}^{-1}(z_{p,k}-{\widehat{z}}_{p,k|k-1})\le \gamma \}$$

(27)

where $\gamma$ is the gate threshold. Suppose $m_{p,k}(m_{p,k}\ge 0)$ observations fall into the validated region (27) at time $k$. Define validate observations ${\stackrel{⌣}{z}}_{p,k}$, i.e.,

$${\stackrel{⌣}{z}}_{p,k}=z_{p,k}^{(j)},j=1,2,\dots ,m_{p,k}$$

(28)

Actually, only one of the above observations is related to the real source; the others are due to noise or reverberation, or none of them are related to the real source. Correspondingly, for $m_{p,k}$ validated observations, there maybe be $m_{p,k}+1$ possible hypothesis, i.e.,

$$\{\begin{array}{l}{\mathrm{H}}_{p,0}=\left\{\mathrm{All} \mathrm{observations} \mathrm{are} \mathrm{independent} \mathrm{of} \mathrm{real} \mathrm{sound} \mathrm{sources} \right\},j=0\\ {\mathrm{H}}_{p,j}=\left\{z_{p,k}^{(j)} \mathrm{is} \mathrm{associated} \mathrm{with} \mathrm{the} \mathrm{true} \mathrm{source}\right\},j=1,2,\dots ,m_{p,k}\end{array}$$

(29)

According to Equation (29), the equation for calculating ${\widehat{x}}_{p,k|k}^{(j)},j=0,1,\dots ,m_{p,k}$ is as follows:

$${\widehat{x}}_{p,k|k}={\displaystyle \sum _{j=0}^{m_{p,k}}\mathrm{E}\{x_{p,k}|{\mathrm{H}}_{p,j},z_{p,1:k}\}p({\mathrm{H}}_{p,j}|z_{p,1:k})}={\displaystyle \sum _{j=0}^{m_{p,k}}{\beta }_{p,k}^{(j)}{\widehat{x}}_{p,k|k}^{(j)}}$$

(30)

where ${\beta }_{p,k}^{(j)}\triangleq p({\mathrm{H}}_{p,j}|z_{p,1:k})$ is the prior probability of event ${\mathrm{H}}_{p,j}$, $0\le {\beta }_{p,k}^{(j)}\le 1$, and ${\displaystyle \sum _{j=0}^{m_{p,k}}{\beta }_{p,k}^{(j)}=1}$,${\widehat{x}}_{p,k|k}^{(j)}\triangleq \mathrm{E}\{x_{p,k}|{\mathrm{H}}_{p,j},z_{p,1:k}\}$ is the updated estimate conditioned on the event ${\mathrm{H}}_{p,j},j=0,1,\dots ,m_{p,k}$, and

$${\widehat{x}}_{p,k|k}^{(0)}={\widehat{x}}_{p,k|k-1}$$

(31)

$${\widehat{x}}_{p,k|k}^{(j)}={\widehat{x}}_{p,k|k-1}+K_{p,k}{v}_{p,k}^{(j)},j=1,2,\dots ,m_{p,k}$$

(32)

where $v_{p,k}^{(j)}=z_{p,k}^{(j)}-{\widehat{z}}_{p,k|k-1}$ is the innovation related to the observation $z_{p,k}^{(j)}$, $K_{p,k}$ is the Kalman gain of node $p$, and

$$P_{p,k|k-1}^{xz}=\frac{1}{2n}{\displaystyle \sum _{i=1}^{2n}({\chi }_{p,k|k-1}^i-{\widehat{x}}_{p,k|k-1}){({\widehat{Z}}_{p,k|k-1}^i-{\widehat{z}}_{p,k|k-1})}^T}$$

(33)

$$K_{p,k}=P_{p,k|k-1}^{xz}{(P_{p,k|k-1}^{zz})}^{-1}$$

(34)

where $P_{p,k|k-1}^{xz}$ is the cross covariance between the state and observation $z_{p,k}$ of node $p$.

Given the innovation $v_{p,k}^{(j)}$ and its covariance $P_{p,k|k-1}^{zz}$, the probability ${\beta }_{p,k}^{(j)}$ is generally computed as [30]

$$\begin{array}{l}{\beta }_{p,k}^{(j)}=\{\begin{array}{c}\frac{e_{p,j}}{b_p+{\displaystyle {\sum }_{i=1}^{m_{p,k}}{e}_{p,i}}},j=1,2,\dots ,m_{p,k}\\ \frac{b_p}{b_p+{\displaystyle {\sum }_{i=1}^{m_{p,k}}{e}_{p,i}}},j=0\end{array}\\ e_{p,j}=e^{-\frac{1}{2}{(v_{p,k}^{(j)})}^T{(P_{p,k|k-1}^{zz})}^{-1}{v}_{p,k}^{(j)}}\\ b_p=\lambda |2\pi P_{p,k|k-1}^{zz}{|}^{\frac{1}{2}}\frac{1-P_{p,D}{P}_G}{P_{p,D}}\end{array}$$

(35)

where $\lambda$ is the spatial probability, $P_{p,D}$ is the probability that the acoustic source is detected by sensor $p$, and $P_G$ is the gate probability.

Finally, the state estimate value ${\widehat{x}}_{p,k|k}$ and error covariance ${\widehat{P}}_{p,k|k}$ can be obtained by

$${\widehat{x}}_{p,k|k}={\widehat{x}}_{p,k|k-1}+K_{p,k}{v}_{p,k}$$

(36)

$${\widehat{P}}_{p,k|k}={\beta }_{p,k}^{(0)}{\widehat{P}}_{p,k|k-1}+(1-{\beta }_{p,k}^{(0)}){\dot{P}}_{p,k|k}+{\ddot{P}}_{p,k|k}$$

(37)

where $v_{p,k}={\displaystyle {\sum }_{j=1}^{m_{p,k}}{\beta }_{p,k}^{(j)}{v}_{p,k}^{(j)}}$ is the probability weighted innovation, and the covariances ${\dot{P}}_{p,k|k}$ and ${\ddot{P}}_{p,k|k}$ are respectively given by [29,30]

$${\dot{P}}_{p,k|k}={\widehat{p}}_{p,k|k-1}-K_{p,k}{P}_{p,k|k-1}^{zz}{K}_{p,k}^T$$

(38)

$${\ddot{P}}_{p,k|k}=K_{p,k}\{{\displaystyle \sum _{j=1}^{m_{p,k}}{\beta }_{p,k}^{(j)}(v_{p,k}^{(j)}-v_{p,k}){(v_{p,k}^{(j)}-v_{p,k})}^T\}{K}_{p,k}^T}$$

(39)

To summarize, the pseudo-code of the PDA-CKF method of using the observations from a single node is depicted in Algorithm 1.

Algorithm 1: PDA-CKF Algorithm

Initialization: ${\widehat{x}}_{p,0|0}={\overline{x}}_0,{\widehat{P}}_{p,0|0}=P_0$

Input: ${\widehat{x}}_{p,k-1|k-1},{\widehat{P}}_{p,k-1|k-1},z_{p,k}$

Output: ${\widehat{x}}_{p,k|k},{\widehat{P}}_{p,k|k}$

Iteration: for $k=1,2,\dots$

1: Prediction step:

2: Compute the state predicted cubature points ${\chi }_{p,k-1|k-1}^i$ at time $k-1$ with (18).

3: Compute the predicted estimate ${\widehat{x}}_{p,k|k-1}$ and covariance ${\widehat{P}}_{p,k|k-1}$ with (20) and (21), respectively.

4: Update step: 5: Compute the state update cubature points ${\chi }_{p,k|k-1}^i$ with (23).

6: Compute the predicted observations ${\widehat{z}}_{p,k|k-1}$ with (25).

7: Compute the innovation covariance $P_{p,k|k-1}^{zz}$ with (26).

8: Select the validated observations ${\stackrel{⌣}{z}}_{p,k}$ according to (28).

9: Compute the cross-covariance $P_{p,k|k-1}^{xz}$ with (33). 10: Compute the Kalman gain $K_{p,k}$ with (34).

11: Compute the association probability ${\beta }_{p,k}^{(j)}$ with (35), $j=1,2,\dots ,m_k$

12: Compute the covariances ${\dot{P}}_{p,k|k}$ and ${\ddot{P}}_{p,k|k}$ with (38) and (39), respectively. 13: Compute the updated estimate ${\widehat{x}}_{p,k|k}$ and covariance ${\widehat{P}}_{p,k|k}$ with (36) and (37), respectively.

The PDA-CKF algorithm makes full use of the observation information of the node itself, which improves the tracking accuracy. However, this algorithm will fail when a node is damaged or the environmental noise and reverberation are severe. Therefore, this paper generalized PDA-CKF to a distributed version that can be used in distributed sensor networks. The improved method was named the probabilistic data association-based distributed cubature Kalman filter (PDA-DCKF). The specific process is shown in Section 3.2.

3.2. PDA-DCKF Algorithm

3.2.1. PDA-DCKF

The neighborhood information of nodes are fused in PDA-DCKF to form local node networks. Then, the local state estimations and error covariances for the local node networks are calculated separately. Finally, the local results are fused to obtain the global state estimation.

On the basis of the above steps, the following is defined:

$$\begin{array}{l}{\widehat{Z}}_{Np,k|k-1}^i={[{\widehat{Z}}_{p,k|k-1}^i;{\widehat{Z}}_{q,k|k-1}^i]}_{num({\mathcal{N}}_{p,k})\times 2n},i=1,2,\dots ,2n,p=1,2,\dots ,N,\\ q\in {\mathcal{N}}_{p,k}=\{q\in \upsilon |(p,q)\in \epsilon \}\cup \{p\}\end{array}$$

(40)

where $q$ represents the neighborhood nodes adjacent to node $p$, $\upsilon =\{1,2,\dots ,N\}$ is the vertex set, $\epsilon \subset \{(p,q)|p,q\in \upsilon \}$ is the edge set of the distributed acoustic sensor network, $num({\mathcal{N}}_{p,k})$ indicates the number of nodes in the neighborhood of node $p$. ${\mathcal{N}}_{p,k}=\{q\in \upsilon |(p,q)\in \epsilon \}\cup \{p\}$ denotes the set of neighbors of node $p$ at time $k$, where a node is a neighbor of itself certainly.

Further, the resulting observations are fused into a matrix. Then, the observed prediction and prediction error variance are, respectively, given by

$${\widehat{z}}_{Np,k|k-1}=\frac{1}{2n}{\displaystyle \sum _{i=1}^{2n}{\widehat{Z}}_{Np,k|k-1}^i}$$

(41)

$$P_{Np,k|k-1}^{zz}=\frac{1}{2n}{\displaystyle \sum _{i=1}^{2n}({\widehat{Z}}_{Np,k|k-1}^i-{\widehat{z}}_{Np,k|k-1}){({\widehat{Z}}_{Np,k|k-1}^i-{\widehat{z}}_{Np,k|k-1})}^T+{[R_{P,k};R_{q,k}]}_{num({\mathcal{N}}_{p,k})\times num({\mathcal{N}}_{p,k})}}$$

(42)

For single node $p$, $v_{p,k}^{(j)}=z_{p,k}^{(j)}-{\widehat{z}}_{p,k|k-1}$ is the innovation vector related to observation $z_{p,k}^{(j)}$, and $K_{p,k}$ is the Kalman gain of node $p$. As far as multiple nodes are concerned, the information of node $p$ and surrounding nodes $q$ is fused to obtain

$$P_{Np,k|k-1}^{xz}=\frac{1}{2n}{\displaystyle \sum _{i=1}^{2n}({\chi }_{p,k|k-1}^i-{\widehat{x}}_{p,k|k-1}){({\widehat{Z}}_{Np,k|k-1}^i-{\widehat{z}}_{Np,k|k-1})}^T}$$

(43)

$$K_{Np,k}=P_{Np,k|k-1}^{xz}{(P_{Np,k|k-1}^{zz})}^{-1}$$

(44)

where $P_{Np,k|k-1}^{xz}$ is the cross covariance between the state and the observed value of node $p$ after fusing the information of neighboring nodes, and $K_{Np,k}$ is the Kalman gain of node $p$ at time $k$ after the fusion.

The probability weighted innovation vector of local nodes is defined as

$$v_{NP,k}={[v_{p,k};v_{q,k}]}_{num({\mathcal{N}}_{p,k})\times 1},p=1,2,\dots ,N,q\in {\mathcal{N}}_{p,k}=\{q\in \upsilon |(p,q)\in \epsilon \}\cup \{p\}$$

(45)

The following is defined as

$${\beta }_{Np,k}^{(0)}=({\beta }_{p,k}^{(0)}+{\displaystyle \sum _{q=1}^{num({\mathcal{N}}_{p,k})-1}{\beta }_{q,k}^{(0)}})/num({\mathcal{N}}_{p,k})$$

(46)

where $\{{\displaystyle \sum _{j=1}^{m_{p,k}}{\beta }_{p,k}^{(j)}(v_{p,k}^{(j)}-v_{p,k}){(v_{p,k}^{(j)}-v_{p,k})}^T\}}$ is defined in the covariance ${\ddot{P}}_{p,k|k}$ of node $p$ as $w_p$; when the information of node $p$ and surrounding nodes is fused, the expression of $w_p$ is computed as

$$\begin{array}{l}{w}_{Np}=\mathrm{diag}{(w_p,w_q)}_{num({\mathcal{N}}_{p,k})\times num({\mathcal{N}}_{p,k})},p=1,2,\dots ,N,\\ q\in {\mathcal{N}}_{p,k}=\{q\in \upsilon |(p,q)\in \epsilon \}\cup \{p\}\end{array}$$

(47)

where $w_q=\{{\displaystyle \sum _{j=1}^{m_{q,k}}{\beta }_{q,k}^{(j)}(v_{q,k}^{(j)}-v_{q,k}){(v_{q,k}^{(j)}-v_{q,k})}^T\}}$.

Finally, the state estimate ${\widehat{x}}_{Np,k|k}$ and the error covariance ${\widehat{P}}_{Np,k|k}$ for node $p$ are expressed as

$${\widehat{x}}_{Np,k|k}={\widehat{x}}_{p,k|k-1}+K_{Np,k}{v}_{Np,k}$$

(48)

$${\widehat{P}}_{Np,k|k}={\beta }_{Np,k}^{(0)}{\widehat{P}}_{p,k|k-1}+(1-{\beta }_{Np,k}^{(0)}){\dot{P}}_{Np,k|k}+{\ddot{P}}_{Np,k|k}$$

(49)

$${\dot{P}}_{Np,k|k}={\widehat{p}}_{p,k|k-1}-K_{Np,k}{P}_{Np,k|k-1}^{zz}{K}_{Np,k}^T$$

(50)

$${\ddot{P}}_{Np,k|k}=K_{Np,k}{w}_{Np}{K}_{Np,k}^T$$

(51)

3.2.2. Fusion Strategy

After calculating the estimation of each local node in the distributed acoustic sensor network, these data need to be fused to obtain a global estimate. Since nodes in a sensor network have different reliabilities, the final tracking result integrates the estimations from the local nodes, which are weighted with the parameters depending on the mean square error of the estimation and the energy of the received signal.

(a)

Energy

The energy of the signal received by each node in the acoustic sensor network is calculated [31], and the equation is described as:

$$E_p=\underset{T\to \infty }{\lim }{\displaystyle {\int }_{-T}^T|x_p(t){|}^{2}dt}$$

(52)

where $x_p(t)$ represents the sound signal received by node $p$. In practice, analog signal $x(t)$ is converted into digital signal $x(n)$, and $x(n)$ needs to be framed and windowed. Then, the framed signal is donated by $x(n)\cdot \omega (n)$. In this paper, the Hamming window was selected for the window function $\omega (n)$. Further, the energy of each frame can be obtained by

$$E_{p,n}={\displaystyle \sum _{m=-\infty }^{\infty }[x_p(m)\omega {(n-m)]}^2}={\displaystyle \sum _{m=-\infty }^{\infty }{x_p}^2(m)h(n-m)}={x_p}^2(n)\ast h(n)$$

(53)

where $h(n)={\omega }^2(n)$, and $E_{p,n}$ represents the short-term energy of node $p$ when the window function starts at the $n-th$ point of the signal. The short-term energy can be regarded as the output of the square of the speech signal passing through a linear filter, and the unit impulse response of the linear filter is $h(n)$.

(b)

MSE

In Equation (48), the local estimate ${\widehat{x}}_{Np,k|k}$ of node $p(p=1,2,\dots ,N)$ is calculated, and ${\widehat{r}}_{p,k}=[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\end{array}]{\widehat{x}}_{Np,k|k}$ is expressed as the estimated acoustic source position of node $p$ at time k. The following is defined:

$${\hat{\mathit{r}}}_{N,k}=\frac{1}{N}\sum _{p=1}^N{\hat{\mathit{r}}}_{p,k}$$

(54)

where ${\hat{\mathit{r}}}_{N,k}$ represents the global position estimation result weighted with the average consensus coefficients and calculates the MSE between the position obtained by each local node and ${\widehat{r}}_{N,k}$, defined as

$$M_p=\left|\right|{\widehat{r}}_{p,k}-{\widehat{r}}_{N,k}|{|}^2$$

(55)

After calculating the energy $E_p$ and the mean square error $M_p$ of node $p$ at time $k$, the following is defined:

$$C_p=\frac{E_p}{M_p}$$

(56)

$${\eta }_p=\frac{C_p}{{\displaystyle {\sum }_{p=1}^N{C}_p}}$$

(57)

where ${\eta }_p$ represents the weight of node $p$ during global fusion. A global consistency analysis is performed on the results obtained by each node according to ${\eta }_p,p=1,2,\dots ,N$:

$${\widehat{x}}_{k|k}={\displaystyle \sum _{p=1}^N{\eta }_p{\widehat{x}}_{Np,k|k}}$$

(58)

$${\widehat{P}}_{k|k}={\displaystyle \sum _{p=1}^N{\eta }_p{\widehat{P}}_{Np,k|k}}$$

(59)

To summarize, the PDA-DCKF is depicted in Algorithm 2.

Algorithm 2: PDA-DCKF Algorithm

Initialization: ${\widehat{x}}_{p,0|0}={\overline{x}}_0,{\widehat{P}}_{p,0|0}=P_0$

Input: ${\widehat{x}}_{k-1|k-1},{\widehat{P}}_{k-1|k-1},z_k$

Output: ${\widehat{x}}_{k|k},{\widehat{P}}_{k|k}$

Iteration: for $k=1,2,\dots$ For any node $p(p=1,2,\dots ,N)$ in sensor network

1: Prediction step:

2: Compute the state predicted cubature points ${\chi }_{p,k-1|k-1}^i$ at time $k-1$ with (18).

3: Compute the predicted estimate ${\widehat{x}}_{p,k|k-1}$ and covariance    ${\widehat{P}}_{p,k|k-1}$ with (20) and (21), respectively.

4: Update step: 5: Compute the state update cubature points ${\chi }_{p,k|k-1}^i$ with (23).

6: Compute the observed values of predicted local nodes ${\widehat{z}}_{Np,k|k-1}$ with (41).

7: Compute the innovation covariance of predicted local nodes $P_{Np,k|k-1}^{zz}$ with (42).

8: Select the validated observations ${\stackrel{⌣}{z}}_{p,k}$ according to (28), $p=1,2,\dots ,N$.

9: Compute the cross-covariance of predicted local nodes $P_{Np,k|k-1}^{xz}$ with (43). 10: Compute the Kalman gain $K_{Np,k}$ with (44). 11: Compute the probability weighted innovation vector of local nodes $v_{Np,k}$ with (45) 12: Compute the association probability ${\beta }_{p,k}^{(j)}$ with (35), $j=1,2,\dots ,m_{p,k}$. 13: Compute the association probability ${\beta }_{Np,k}^{(0)}$ with (46). 14: Compute $v{1}_{Np}$ with (47). 15: Compute the updated estimate ${\widehat{x}}_{Np,k|k}$ and covariance ${\widehat{P}}_{Np,k|k}$ of local nodes with (48) and (49), respectively. 16: Compute the weight ${\eta }_p$ of node $p(p=1,2,\dots N)$ with (57). 17: Compute the updated estimate ${\widehat{x}}_{k|k}$ and covariance ${\widehat{P}}_{k|k}$ with (58) and (59), respectively.

The advantages of probabilistic data association and distributed acoustic sensor networks are combined in the PDA-DCKF proposed in this paper. In this method, the PDA algorithm is used to sift the observations from neighboring nodes. Then, the sifted observations are fused to update the state vectors in the CKF. This method not only makes the observation value obtained by each node more accurate, but also makes full use of the information of neighborhood nodes.

Meanwhile, a weighted fusion method based on local node-received signal energy and position estimation mean square error was proposed. This dynamic weighted consistency fusion considers the reliability of the local state of the nodes and provides a good global estimation performance.

4. Experiments and Results Discussion

To verify the performance of the proposed speaker tracking method, the evaluations are performed in a simulated room environment. Under the same conditions, the comparative experiments between PDA-DCKF and current methods are carried out, including centralized method (CCKF), DUKF, DCKF, iteration based DCKF [20] (DICKF) and DEKF. The results obtained by all methods are the average of 100 Monte Carlo runs.

The root mean square error (RMSE) is used here to evaluate the tracking performance. $r_k$ is expressed as the ground truth value of time $k$, and ${\widehat{r}}_{N,k}$ represents the global consistency position calculated by the acoustic sensor network at this time. The RMSE is defined as [32]

$$\mathrm{RMSE}=\sqrt{\frac{1}{K}{\displaystyle \sum _{k=1}^K\left|\right|r_k-{\widehat{r}}_{N,k}|{|}^2}}$$

(60)

where $K$ denotes the number of frames. Generally, the smaller the RMSE, the better the tracking result.

4.1. Simulation Setups

The simulation environment was a typical room of size 6 m × 6 m × 3 m, with an acoustic sensor network of 12 nodes ($N=12$). Each node contained a pair of microphones 0.5 m apart. The communication diagram of the distributed acoustic sensor network is shown in Figure 2, where the communication radius is 2.5 m, and each circle represents a node. The simulated trajectory 1 was a line from (0.5,0.8) to (2.5, 2.8), and trajectory 2 was an arc from (1, 2) to (4.86, 2.1), as shown in Figure 3. In different experiments, the speech sampled at the frequency of $F_s=16\mathrm{KHZ}$ was used as the acoustic source signal; the speech was a female recording, and the waveform and spectrum of the signal are shown in Figure 4a. The sound speed was $c=342 \mathrm{m}/\mathrm{s}$. The microphone signals were simulated with the Image method [33]. Specifically, different RIRS are generated by virtual sound source method to reflect different reverberation time. These RIRSs were convolved with the speech signal and then added to the Gaussian white noise with a determined mean and covariance to produce a received microphone signal with a mixture of reverberation and noise. The different covariance of Gaussian noise determines the different value of the signal-to-noise ratio (SNR), which reflects different environmental noise conditions. The microphone signal was divided into different signal frames along the sound source track, where the frame length of speech signal was $N_f=512$ and each signal frame was used for state estimation. Taking node 1 as an example, Figure 4b shows the waveform and spectrum of the speech signal received by the first microphone of node 1. For the observation TDOA, a total of eight time delays were chosen according to the magnitude of the GCC peak. From these delays, further TDOA observations were selected, where the relevant parameters were set as $\lambda =10$, $\gamma =4$, $P_G=0.93$, and $P_D=0.95$. The standard deviation of TDOA measurement error was $\sigma =50 \mathsf{\mu }\mathrm{s}$. In the acoustic dynamical model, the parameters were $\beta =10 {\mathrm{s}}^{-1}$ and $\overline{\upsilon }=1 {\mathrm{ms}}^{-1}$. In the average consistency calculation of the global state estimation and its error covariance, the Metropolis weight was used, the number of consistency iterations [34] was $N_{con}=10$, and the number of iterations in the iterative CKF was 3.

This paper conducted four experiments to evaluate the tracking performance of PDA-DCKF. In Experiment 1, trajectory 1 was used as the acoustic source trajectory. The initial prior $p(x_0)$ of the acoustic source position was set as a Gaussian distribution with mean $x_0={[0.5,0.8,0.02,0.02]}^T$ and covariance $P_0=diag([0.05,0.05,0.0025,0.0025])$. In experiment 2, the sound source signal and track were the same as experiment 1. Using simple average fusion rules, the influence of fusion rules on PDA-DCKF tracking performance was discussed. Experiment 3 discussed the robustness of the algorithm. The acoustic source and trajectory were the same as the previous two experiments. In Experiment 4, trajectory 2 was used as the acoustic source track to check the tracking results of the acoustic source when the track was nonlinear.

4.2. Simulation Results

4.2.1. Experiment 1

In this experiment, the tracking performance was evaluated under different ambient and reverberant conditions. First, the impact of environmental noise on tracking performance was investigated. Figure 5 depicts the RMSE results as a function of SNR for a reverberation time of T₆₀ = 200 ms. In Figure 5, it is observed that the RMSE of all methods decreases with the increase of SNR, which means that the tracking accuracy increases with the increase of SNR. This is because when the SNR becomes larger, the microphone signal is less affected by ambient noise, resulting in better tracking performance. In addition, under the same SNR, PDA-DCKF performs better than traditional distributed Kalman filtering, such as extended Kalman filtering, unscented Kalman filtering, and cubature Kalman filtering. Since only one time-delayed observation of the GCC largest peak is used in traditional methods, peaks associated with real sources may be masked by spurious peaks caused by noise or reverberation, resulting in erroneous state estimates. In contrast, multiple time-difference observations of multiple largest peaks of GCC are employed in PDA-DCKF, resulting in ideal tracking performance. At the same time, compared with DICKF in this experiment, the results show that the effect of PDA-DCKF is better than that of DICKF. Because DICKF is aimed at the DCKF method, and DCKF has problems such as slow response speed and low tracking accuracy. However, the tracking performance and convergence speed of the algorithm can be improved through several local iterations in DICKF. However, still only one time-delay observation of the GCC largest peak is used in DICKF, which also causes it to be inaccurate, but as can be seen from Figure 5, as the SNR increases, the gap between DICKF and PDA-DCKF becomes smaller because the observations are more reliable when the SNR becomes larger. In addition, Figure 5 shows that PDA-DCKF is not as good as CCKF because the observation information of all nodes is used in CCKF, but PDA-DCKF achieved an effect very close to the CCKF effect, and its computational cost and the burden of the network is less than that of CCKF.

The effect of reverberation on tracking performance was also studied in this paper. Figure 6 depicts the RMSE results as a function of T₆₀ with SNR = 20 dB. From the results, we can observe that the RMSEs of all the methods increased as T₆₀ became larger, which signifies the degradation of the tracking accuracies. This may be because the microphone signal is more affected by reverberation as T₆₀ becomes larger, the time difference observations extracted from only the largest peak or multiple largest peaks are not reliable, and the tracking performance of these methods deteriorates. In addition, it can be found from Figure 6 that the tracking performance of PDA-DCKF is better than DEKF, DUKF, DCKF, and DICKF. In fact, in traditional methods, the time-difference observations included in the scheme are only extracted from the largest peak of the GCC, while the peaks associated with the true hypocenter may be masked by false peaks caused by reverberation. In contrast, PDA-DCKF incorporates TDOA observations of multiple largest peaks of GCC into the scheme, which can alleviate the adverse effects of reverberation to a certain extent. Furthermore, the effect is not as good as CCKF showed in Figure 6, but it also achieves a very close effect.

4.2.2. Experiment 2

The effect of the fusion strategy proposed in this paper on the results is discussed in Experiment 2. When PDA-DCKF adopts a simple average fusion rule, it is called PDA-DCKF-avg. In this section, different SNR and different reverberations are used to test the effectiveness of the fusion strategy. The experimental results are shown in Figure 7 and Figure 8.

As depicted in Figure 7, with the increase of the SNR, the RMSEs for the PDA-DCKF methods with both these two fusion strategies decrease, but the proposed one is more effective. Figure 8 also shows that, with the increase of the reverberation time, the error also increases. In addition, only under 50 ms reverberation, the error of the average fusion strategy is smaller than that proposed in this paper, and the fusion strategy proposed in this paper was better than the average fusion effect under 100–600 ms. Comparing Figure 5, Figure 6, Figure 7 and Figure 8, it can be found that, even if the average fusion strategy is used, the PDA-DCKF in this paper is still smaller than the error obtained by the above comparison test, which further proves the effectiveness of the method in this paper.

4.2.3. Experiment 3

In practical applications, a network may be damaged by nodes, and when a node in a network is damaged, whether the network can still work normally will test the robustness of the system. In this subsection, the node damage in the distributed acoustic sensor network is simulated, and the tracking results of the acoustic source after the damage are compared with those before the damage. When node 1 in the network is damaged, it is called graph $\mathcal{G}2$, as shown in Figure 9a. When node 1 and node 6 in the acoustic sensor network are damaged, it is called graph $\mathcal{G}3$, as shown in Figure 9b. The experimental results are shown in

Table 1. RMSE versus SNR under different graphs with T₆₀ = 200 ms.

SNR (dB)	$\mathcal{G}1$	$\mathcal{G}2$	$\mathcal{G}3$
5	0.3363	0.3521	0.3637
10	0.1735	0.1803	0.2057
15	0.1457	0.1511	0.1786
20	0.1201	0.1284	0.1543
25	0.1161	0.1203	0.1457
30	0.1101	0.1169	0.1376

and

Table 2. RMSE versus T₆₀ under different graphs with SNR = 20 dB.

T60 (ms)	$\mathcal{G}1$	$\mathcal{G}2$	$\mathcal{G}3$
50	0.0992	0.1013	0.1164
100	0.1108	0.1195	0.1351
200	0.1201	0.1284	0.1543
300	0.1393	0.1501	0.1754
400	0.1824	0.1903	0.2158
500	0.2187	0.2305	0.2439
600	0.4703	0.4897	0.5062

.

It can be seen from Table 1 and Table 2 that the acoustic source can still be tracked in the case of node damage. Although the accuracy has decreased, the amplitude of the drop is not large and the acoustic source can still be tracked accurately. This can prove that the method proposed in this paper has good robustness under this network.

4.2.4. Experiment 4

In order to further verify the effectiveness of the algorithm in this paper, the semicircle of trajectory 2 was used as the acoustic source trajectory, and comparative experiments were carried out under different SNR and reverberation. The experimental data are shown in

Table 3. RMSE versus SNR with T₆₀ = 200 ms.

SNR (dB)	RMSE
5	0.3751
10	0.1869
15	0.1423
20	0.1299
25	0.1214
30	0.1167

and

Table 4. RMSE versus T₆₀ with SNR = 20 dB.

T60 (ms)	RMSE
50	0.1174
100	0.1251
200	0.1299
300	0.1322
400	0.1635
500	0.2216
600	0.5027

. Figure 10 shows the tracking results with SNR = 15 dB and T₆₀ = 400 ms.

From the above Table 3 and Table 4 and Figure 10, it can be seen that the algorithm in this paper can still accurately track the sound source in the face of such a strong nonlinear trajectory.

5. Conclusions

An improved PDA-DCKF method was proposed in this paper, which proved to be able to solve the problem of tracking a single mobile acoustic source with distributed acoustic sensor networks in the noise and reverberation environment. First, in order to reduce the adverse effects of noise and reverberation, the prediction value of observation is obtained by using the prediction state and the observation model of distributed nodes. Secondly, the actual observations are screened according to the predicted value. Multiple TDOA observations are extracted at each node and incorporated into the status update of CKF through PDA to generate PDA-CKF. PDA-CKF was applied to distributed acoustic sensor networks, and PDA-DCKF was further developed. In PDA-DCKF, the PDA algorithm is first used to sift the observations from neighboring nodes. Then, the sifted observations are fused to update the state vectors in the CKF. Each node runs PDA-DCKF for local state estimation and TDOA observation. Then, a new fusion strategy is proposed using energy and MSE to merge all single local estimates in a distributed manner for global state estimation. In order to apply the improved PDA-DCKF to the acoustic source tracking problem, the Langevin model was used to model the acoustic source dynamics, and a method to extract the time difference observation was proposed. Finally, a distributed acoustic source tracking framework was obtained. In order to evaluate the effectiveness of PDA-DCKF in acoustic source tracking, comparative experiments were carried out with existing methods (DCKF, DUKF, DEKF, and DICKF) under different ambient noise and reverberation conditions. The results show that the PDA-DCKF has better tracking performance than DCKF, DUKF, DEKF, and DICKF under most noise and reverberation conditions. In addition, the PDA-DCKF achieved the same tracking performance as the centralized CKF. Furthermore, it can even track the acoustic source stably in the case of node damage.