Optimal Passive Source Localization for Acoustic Emissions

Abstract

Acoustic emission is a non-destructive testing method where sensors monitor an area of a structure to detect and localize passive sources of elastic waves such as expanding cracks. Passive source localization methods based on times of arrival (TOAs) use TOAs estimated from the noisy signals received by the sensors to estimate the source position. In this work, we derive the probability distribution of TOAs assuming they were obtained by the fixed threshold technique—a popular low-complexity TOA estimation technique—and show that, if the sampling rate is high enough, TOAs can be approximated by a random variable distributed according to a mixture of Gaussian distributions, which reduces to a Gaussian in the low noise regime. The optimal source position estimator is derived assuming the parameters of the mixture are known, in which case its MSE matches the Cramér–Rao lower bound, and an algorithm to estimate the mixture parameters from noisy signals is presented. We also show that the fixed threshold technique produces biased time differences of arrival (TDOAs) and propose a modification of this method to remove the bias. The proposed source position estimator is validated in simulation using biased and unbiased TDOAs, performing better than other TOA-based passive source localization methods in most scenarios.

1. Introduction

Acoustic emission testing is a non-destructive testing method used to detect several kinds of faults in structures such as piping, bridges and aerospace structures. Its main advantages over other non-destructive testing methods are the high sensibility, the wide coverage area and the possibility of monitoring structures in real time [1,2,3,4,5,6]. In the acoustic emission framework, sensors are spread over the surface to be monitored, and elastic waves emitted by passive sources (such as cracks or rivets under stress) propagate through the medium and are detected by these sensors. The segment of signal acquired by a sensor that is detected as a wave is called a hit. In order to reduce memory requirements, the sampled signals are usually not saved, only the estimated times of arrival (TOAs) are recorded and later used to estimate the source position.

Many TOA estimation algorithms are presented in the literature [7,8,9,10], but the fixed threshold method [11,12,13] is a very popular technique due to its simplicity and low complexity. Defining $r_i\left[n\right]$ as the signal sampled with period T by the i-th sensor at instant $t=nT$, this method estimates the TOA $t_i$ by comparing the absolute value of the signal to a fixed threshold K:

$$t_i=n_{i}T,n_i=\underset{n}{min}\{n:|r_i\left[n\right]|\ge K\}.$$

(1)

Thus, the estimated TOA is the first instant the absolute value of the signal crosses the threshold. The value of K should be high enough to be larger than the noise level most of the time, but not too high so that hits will be missed.

There are many approaches in the literature to estimate the source position given the TOAs [14,15,16,17,18,19,20], usually based on the minimization of different choices of cost function that may receive as input the measured TOAs or the time-differences of arrival (TDOAs, defined as the difference of TOAs measured by different sensors). One of the main contributions of this paper is to compare different cost functions and show that they do not usually lead to optimal estimates of source localization in general, except in the particular situation of small noise variance. We concentrate here in the isotropic case, i.e., we assume that the wave velocity is known and independent of the direction, but there are algorithms that are applicable to the anisotropic and unknown velocity case [14,21,22].

Passive source localization algorithms based on TOAs usually assume a Gaussian distribution for the uncertainties in the TOA estimates and a zero-mean Gaussian distribution for TDOA estimates [20,21,22,23,24,25,26]. We prove theoretically that this assumption holds approximately for the threshold method in the case of low noise. In a more general setting, the TOA distribution depends on the method used to estimate it, and we prove in Section 3 that for larger noise variances the fixed threshold method produces TOAs that can be modeled as a mixture of Gaussians (see also [27]). Furthermore, TDOAs estimated by the fixed threshold method may present a large bias since signals reach sensors with different amplitudes due to attenuation (see Section 4 and [27]).

The fixed threshold method only works if the threshold K is sufficiently higher than the noise level so that noise has a low probability of triggering the threshold without an incoming wave. However, choosing a high value for K increases the TOA bias, given by the difference between the expected value of the estimated and the actual time of arrival of the wave. If all sensors received hits with similar amplitudes, the TOA bias would be constant accross sensors, and there would be no TDOA bias. However, because sensors lie at different distances from the acoustic emission source, they usually receive hits with different amplitudes due to attenuation. As such, TOA bias varies accross sensors, generating a bias in TDOA estimates, which may bias the position estimate. An illustration of the fixed threshold method and the bias in TOA and TDOA estimates is presented in Figure 1.

The contributions of this paper are:

• To derive the probability distribution for TOA estimates obtained through the fixed threshold method. We show that the TOAs can be well approximated by a mixture of Gaussian distributions when (as is usually the case) the sensor has a bandpass frequency response, such that the observed signals may be well described by a model of the form

  $$s\left(t\right)=e\left(t\right)sin(2\pi f_{0}t+\varphi ),$$

  (2)

  where $e\left(t\right)$ is a low-pass envelop function and $f_0$ is the peak of the sensor frequency response. The mixture of Gaussians is a reasonable approximation if the sampling rate $f_s$ satisfies $f_s\gg 4f_0$. We also show that the Gaussian distribution usually assumed is a good approximation in the low-noise regime.
• To propose a TOA estimation method that generates unbiased Time Differences of Arrival (TDOAs) for structures where all frequencies propagate in the medium with the same velocity.
• To derive the Cramér–Rao lower bound for the source position assuming TOAs follow a mixture of Gaussian distributions.
• To derive the optimal TOA-based source position estimator on a surface (that is, the estimator whose covariance matrix is the inverse of the Fisher information matrix) assuming TOAs follow a mixture of Gaussian distributions whose parameters are known given the source position. A procedure to estimate these parameters from noisy signals is also presented. The optimal estimator is tested in simulation in several scenarios and compared with other passive source localization methods based on TOAs.
• To show that the optimal estimator reduces to one of the algorithms available in the literature if the TOA estimates are unbiased and the measurement noise is sufficiently small, in which case TOAs are approximately Gaussian-distributed.

Notation: We denote the expected value operator as $\mathbb{E}\{·\}$, the probability of an event $\mathcal{A}$ as $\mathbb{P}\left\{\mathcal{A}\right\}$, the multivariate normal distribution with mean vector $\mu$ and covariance matrix $\mathbf{C}$ as $\mathcal{N}(\mu ,\mathbf{C})$, the covariance matrix of a random vector $\mathbf{x}$ as $\mathrm{cov}\left(\mathbf{x}\right)$, the trace of matrix $\mathbf{A}$ as $\mathrm{tr}\left\{\mathbf{A}\right\}$. $\mathrm{diag}(a_1,\cdots ,a_N)$ is the diagonal matrix $\mathbf{A}$ such that ${\mathbf{A}}_{ii}=a_i$, ${\mathbf{1}}_{\mathcal{A}}\left(x\right)$ is the indicator function, equal to 1 if $x\in \mathcal{A}$ and 0 otherwise, and ${\mathcal{A}}^c$ is the complement of the set $\mathcal{A}$. Vectors and matrices are written in boldface, and the derivative of the k-th element of a vector $\mathbf{v}$ with respect to x is denoted as ${\mathbf{v}}_{k,x}$.

2. Source Localization Methods Based on TOAs

Hits identified at different sensors with similar TOAs are grouped into a single event related to a possible source [12]. Assuming the grouping is correct, the source position can be estimated by solving the system of equations:

$$t_k=t+{\tau }_k\left(\mathbf{x}\right),k=1,2,\cdots ,N,$$

(3)

where $t_k$ is the TOA at sensor k, N is the number of sensors, t is the instant the event occurred (an unknown variable) and ${\tau }_k\left(\mathbf{x}\right)$ is the time the wave on the surface under test takes to propagate from the tentative source position $\mathbf{x}={\left[\begin{array}{cc}x& y\end{array}\right]}^T$ to the position ${\mathbf{x}}_k$ of the k-th sensor assuming the propagation is isotropic with known velocity c:

$${\tau }_k\left(\mathbf{x}\right)=\frac{1}{c}∥\mathbf{x}-{\mathbf{x}}_k∥.$$

(4)

The system of Equations (3) has three unknowns ($x,y$ and t) and N equations, thus it is overdetermined for $N>3$ unless there are co-linear sensors. Since due to noise the system usually has no solution, the problem is recast as an optimization problem using an appropriate cost function to obtain an approximate solution. Three popular cost functions for TOA-based source localization are (5), (8) and (9) below [12,20,28]. The quality of the position estimates depends on the choice of the cost function, as well as on the distribution of the uncertainty in the TOA estimates $t_i$.

The least squares solution of (3) is obtained by minimizing the cost function $J_{\mathrm{TOA}}(\mathbf{x},t)$, defined as

$$J_{\mathrm{TOA}}(\mathbf{x},t)=\sum _{k=1}^N{(t_k-t-{\tau }_k\left(\mathbf{x}\right))}^2.$$

(5)

$J_{\mathrm{TOA}}(\mathbf{x},t)$ is a generalization of the cost function $J_{\mathrm{TOA}}^{\prime }\left(\mathbf{x}\right)=c^2{\sum }_{k=1}^N{(t_k-{\tau }_k\left(\mathbf{x}\right))}^2$ [23] that is used in applications where the instant of emission of the wave $t^{*}$ is known, in which case it is common to use $t^{*}=0$ without loss of generality. In our application, $t^{*}$ is unknown and it must be taken into account to localize the source. Fortunately, it is possible to obtain a cost function equivalent to (5) that only depends on $\mathbf{x}$ by solving for t as a function of $\mathbf{x}$. In fact, setting $\frac{\partial }{\partial t}{J}_{\mathrm{TOA}}(\mathbf{x},t)=-2{\sum }_{k=1}^N\left(t_k-t-{\tau }_k\left(\mathbf{x}\right)\right)$ to zero and isolating t yields $t=\frac{1}{N}{\sum }_{k=1}^N(t_k-{\tau }_k\left(\mathbf{x}\right))$.

Hence, minimizing $J_{\mathrm{TOA}}(\mathbf{x},t)$ is equivalent to minimizing:

$$J_{\mathrm{TOA}}\left(\mathbf{x}\right)=\sum _{k=1}^N{(t_k-{\tau }_k\left(\mathbf{x}\right)-t\left(\mathbf{x}\right))}^2,$$

(6)

where $t\left(\mathbf{x}\right)=\frac{1}{N}{\sum }_{k=1}^N(t_k-{\tau }_k\left(\mathbf{x}\right))$. Another approach to approximately solve (3) is to subtract the first equation from the others to eliminate t, yielding a modified system of $N-1$ equations in two unknowns whose k-th equation is:

$$t_k-t_1={\tau }_k\left(\mathbf{x}\right)-{\tau }_1\left(\mathbf{x}\right),2\le k\le N.$$

(7)

The least squares solution of (7) is obtained by minimizing the cost function $J_{\mathrm{TDOA}}\left(\mathbf{x}\right)$ below, which only depends on the measured time differences of arrival (TDOAs) $t_k-t_1$ [12,29]:

$$J_{\mathrm{TDOA}}\left(\mathbf{x}\right)=\sum _{k=2}^N{(t_k-t_1-({\tau }_k\left(\mathbf{x}\right)-{\tau }_1\left(\mathbf{x}\right)))}^2.$$

(8)

It is also possible obtain an approximate solution to (3) using a constrained least squares approach, as done in [28]. Without loss of generality, consider that sensor 1 is positioned at $(0,0)$. Rewriting each equation from (7) as ${(t_k-t_1+{\tau }_1\left(\mathbf{x}\right))}^2={\tau }_k{\left(\mathbf{x}\right)}^2$, substituting (4), expanding the squares and performing simple algebraic manipulations, we obtain $d_k∥\mathbf{x}∥+{\mathbf{x}}_k^T\mathbf{x}=\frac{{∥x_k∥}^2-d_k^2}{2}$, where $d_k=c(t_k-t_1)$ is the measured range difference between sensors k and 1.

In order to find an approximate solution for this over-determined system of equations, Ref [28] defines the auxiliary unknown variable $\mathbf{z}={\left[\begin{array}{cc}∥\mathbf{x}∥& \mathbf{x}\end{array}\right]}^T$, which is estimated by defining the cost function $J_{CLS}\left(\mathbf{z}\right)$ (where the acronym ’CLS’ stands for ’constrained least squares’) as

$$J_{CLS}\left(\mathbf{z}\right)={∥\mathbf{A}\mathbf{z}-\mathbf{b}∥}^2,$$

(9)

where $\mathbf{A}=\left[\begin{array}{cc}{d}_2& {\mathbf{x}}_2^T\\ ⋮& ⋮\\ d_N& {\mathbf{x}}_N^T\end{array}\right]$ and $\mathbf{b}=\frac{1}{2}{\left[\begin{array}{ccc}{∥{\mathbf{x}}_2∥}_2^2-d_2^2& \cdots & {∥{\mathbf{x}}_N∥}_2^2-d_N^2\end{array}\right]}^T$, and solving the constrained optimization problem:

$$\widehat{\mathbf{z}}=arg\underset{\mathbf{z}}{min}{J}_{CLS}\left(\mathbf{z}\right)\mathrm{s}.\mathrm{t}.\left\{\begin{array}{c}{z}_1^2-z_2^2-z_3^2=0\hfill \\ z_1>0\hfill \end{array}\right.,$$

(10)

in which we denote by $z_i$ the i-th element of vector $\mathbf{z}$. A low-complexity method for obtaining an approximate solution for (10) is proposed in [28], but since the focus of this paper is on the quality of the estimators, in our simulations we solve (10) directly.

3. TOA Probability Distribution

In [27], the probability mass function (pmf) for TOAs obtained by the fixed threshold algorithm was derived considering the measurement noise to be white. In this section, we extend that result, deriving the TOA pmf for non-white noise, and obtaining a simplified expression in the case where the noise is white.

3.1. TOA Pmf for Correlated Noise

Let us model the signal $r\left[n\right]$ sampled by a sensor as the sum of a deterministic component $s\left[n\right]$ and a zero-mean noise $w\left[n\right]$, that is, $r\left[n\right]=s\left[n\right]+w\left[n\right]$. To simplify the notation, the index i representing the sensor is omitted.

The fixed threshold method estimates the TOA as the first instant the signal $\left|r\right[n\left]\right|$ crosses the threshold. Thus, denoting as $p\left[n\right]$ the probability of the estimated TOA being the instant n, we have:

$$p\left[n\right]=\mathbb{P}\left\{\right|r\left[n\right]|\ge K\mathrm{and}|r\left[m\right]|<K\forall m<n\}.$$

(11)

Note that the threshold level should be adjusted so that the noise cannot trigger the threshold by itself [11]. Assume then that the threshold is well adjusted (i.e., $\mathbb{P}\left\{\right|w_k\left[m\right]|\ge K\}\approx 0\forall m$). Then, $\left|r\right[n\left]\right|\ge K$ if and only if:

$$w\left[n\right]sgn\left(s\left[n\right]\right)\ge q_n,q_n\stackrel{\mathrm{def}}{=}K-\left|s\left[n\right]\right|.$$

(12)

Let us prove (12): If $\left|r\right[m\left]\right|\ge K$, then $sgn\left(r\right[n\left]\right)=sgn\left(s\right[n\left]\right)$ because $w\left[n\right]\in [-K,+K]$. Thus,

$$\left|r\right[n\left]\right|=\left(s\right[n]+w[n\left]\right)sgn\left(r\right[n\left]\right)=\left|s\right[n\left]\right|+w\left[n\right]sgn\left(w\right[n\left]\right)\ge K,$$

(13)

leading to:

$$w\left[n\right]sgn\left(s\left[n\right]\right)\ge q_n.$$

(14)

To prove the converse of (12), assume that $w\left[n\right]sgn\left(s\right[n\left]\right)\ge K-\left|s\right[n\left]\right|$. If $s\left[n\right]\ge 0$, this implies $s\left[n\right]+w\left[n\right]\ge K$, thus $w\left[n\right]sgn\left(s\left[n\right]\right)\ge q_n$. On the other hand, if $s\left[n\right]<0$, we have $-w\left[n\right]-s\left[n\right]=w\left[n\right]sgn\left(s\right[n\left]\right)+\left|s\right[n\left]\right|\ge K$, also leading to $w\left[n\right]sgn\left(s\left[n\right]\right)\ge q_n$.

Our goal is to write (11) in terms of the noise joint cumulative distribution. For this, define the intervals ${\mathcal{I}}_m=\left]-\infty ,q_m\right[$ and the modified noise $\tilde{w}\left[m\right]=sgn\left(s\left[m\right]\right)w\left[m\right]$ for $0\le m\le n$. Using these definitions, we conclude from (12) that $\left|r\right[n\left]\right|<K$ if and only if $\tilde{w}\left[n\right]<q_m$, thus $p\left[n\right]$ is the probability of $(\tilde{w}\left[n\right],\tilde{w}[n-1],\cdots ,\tilde{w}\left[0\right])$ belonging to ${\mathcal{I}}_n^c\times {\mathcal{I}}_{n-1}\times \cdots \times {\mathcal{I}}_0$. Define the vector of noise samples $\mathbf{w}\left[n\right]={\left[\begin{array}{cccc}w\left[0\right]& w\left[1\right]& \cdots & w\left[n\right]\end{array}\right]}^T$. We consider that the joint noise pdf $f_{\mathbf{w}\left[n\right]}(w_0,\cdots ,w_n)$ is an even function with respect to each variable, that is, $f_{\mathbf{w}\left[n\right]}(w_0,\cdots ,w_j,\cdots ,w_n)$ = $f_{\mathbf{w}\left[n\right]}(w_0,\cdots ,-w_j,\cdots ,w_n)$ for all j (e.g., a zero-mean Gaussian distribution). In this case, $(\tilde{w}\left[n\right],\cdots ,\tilde{w}\left[0\right])$ has the same distribution as $\mathbf{w}\left[n\right]$. Hence, we have $p\left[n\right]=\mathbb{P}\{\mathbf{w}\left[n\right]\in {\mathcal{I}}_n^c\times {\mathcal{I}}_{n-1}\times {\mathcal{I}}_{n-2}\times \cdots \times {\mathcal{I}}_0\}$, which can be expanded into:

$$\begin{array}{cc}\hfill p\left[n\right]& ={\int }_{{\mathcal{I}}_0\times \cdots \times {\mathcal{I}}_{n-1}\times {\mathcal{I}}_n^c}{f}_{\mathbf{w}\left[n\right]}(w_0,\cdots ,w_n)\mathrm{d}{w}_0\cdots \mathrm{d}{w}_n\hfill \\ \hfill & ={\int }_{-\infty }^{+\infty }{\int }_{{\mathcal{I}}_0\times \cdots \times {\mathcal{I}}_{n-1}}{f}_{\mathbf{w}\left[n\right]}(w_0,\cdots ,w_n)\mathrm{d}{w}_0\cdots \mathrm{d}{w}_n\hfill \\ \hfill & -{\int }_{{\mathcal{I}}_0\times \cdots \times {\mathcal{I}}_{n-1}\times {\mathcal{I}}_n}{f}_{\mathbf{w}\left[n\right]}(w_0,\cdots ,w_n)\mathrm{d}{w}_0\cdots \mathrm{d}{w}_n.\hfill \end{array}$$

(15)

As the second integral in (15) is equal to the joint cumulative distribution of $\mathbf{w}\left[n\right]$ evaluated at $(q_0,q_1,\cdots ,q_n)$ and the first integral is the cumulative distribution of $\mathbf{w}[n-1]$ evaluated at $(q_0,q_1,\cdots ,q_{n-1})$, we obtain:

$$p\left[n\right]=F_{\mathbf{w}[n-1]}(q_0,\cdots ,q_{n-1})-F_{\mathbf{w}\left[n\right]}(q_0,\cdots ,q_n),$$

(16)

where $F_{\mathbf{w}\left[n\right]}(·)$ is the joint cumulative distribution of $\mathbf{w}\left[n\right]$.

3.2. TOA Pmf for Independent and Identically Distributed Noise

If the noise is an independent and identically distributed (i.i.d.) stochastic process, a much simpler expression for (16) can be obtained. In this case, denoting the cumulative density function (cdf) of a single noise sample as $F_W\left(w\right)$, the joint cumulative distribution $F_{\mathbf{w}\left[n\right]}\left(\mathbf{w}\right)$ can be written as the product of the cdfs of each noise sample, that is, $F_{\mathbf{w}\left[n\right]}\left(\mathbf{w}\right)={\prod }_{\ell =0}^n{F}_W\left({\mathbf{w}}_{\ell }\right)$, and thus:

$$p\left[n\right]=(1-F_W\left(q_n\right))\prod _{\ell =0}^{n-1}{F}_W\left(q_{\ell }\right),q_{\ell }=K-\left|s[\ell ]\right|.$$

(17)

Note that $q_n$ and $p\left[n\right]$ depend on the noiseless signal $s\left[n\right]$, which is unknown in practical applications. In Section 5.3, $q_n$ is estimated using the noisy signal $r\left[n\right]$ instead of $s\left[n\right]$, and (17) is used to build a source position estimator.

3.3. Approximating the TOA Pdf as a Gaussian Mixture

Many authors consider that TOAs follow a Gaussian distribution to derive passive source localization methods [23,24,25,26], but (17) is not only discrete, but may also be multimodal. In this section, we investigate in which cases the Gaussian distribution is a good approximation to the TOA pmf (17). We show via an example that, for high sampling rate and small noise variance, TOAs can be well approximated by a Gaussian distribution, and for high sampling rate and high noise variance they can be approximated by a mixture of Gaussian distributions. This approximation is exploited in Section 5.1, where the expression for the optimal TOA-based estimator for TOAs distributed according to a mixture of Gaussian distributions is derived.

Using a sampling rate of $10\mathrm{MHz}$ (sampling period $T=0.1\mathsf{\mu }\mathrm{s}$), the received signal was generated following (2) as $r\left[n\right]={sin}^2\left(\frac{\pi (n-1)}{1000}\right)sin\left(\frac{2\pi (n-1)}{67}\right)+w\left[n\right]$ (a von Hann window modulated with a frequency of $f_0=\frac{1}{67\mathsf{\mu }\mathrm{s}}\approx 150\mathrm{kHz}$, a common resonance frequency for piezoelectric sensors used in acoustic emission), where $w\left[n\right]$ is a Gaussian white noise with standard deviation $\frac{\mathsf{\mu }\mathrm{s}K}{5}=0.06$, and $K=0.3$ is the threshold. We then applied the fixed threshold method to $r\left[n\right]$ to obtain the experimental TOAs and generate the empirical TOA pmf and compare with the theoretical expression (17).

The obtained TOA pmf and cdf for this hit is shown in Figure 2. The pmf $p\left[n\right]$ was fitted into a Gaussian Mixture Model (GMM), leading to the pdf $p_{\mathrm{GMM}}\left(t\right)$. The details of the fitting technique we used are described in Appendix C. The cdf of the fitted distribution was plotted along with the TOA cdf, and in order to allow the comparison between the fitted pdf $p_{\mathrm{GMM}}\left(t\right)$ and $p\left[n\right]$, we plotted $T\times p_{\mathrm{GMM}}\left(t\right)$ along with $p\left[n\right]$.

We conclude from Figure 2 that the TOA distribution is well approximated by a GMM. If the noise level is high enough, the threshold may be triggered in different oscillation periods of the signal, creating secondary lobes spaced by the oscillation period of $\left|r\right[n\left]\right|$ (half the oscillation period of $r\left[n\right]$), explaining the different Gaussian components of $p\left[n\right]$. On the other hand, if the noise were small, $p\left[n\right]$ would present a single lobe, reducing to a Gaussian distribution.

Since we are approximating a discrete probability distribution by a continuous distribution, this approximation is only reasonable if the sampling rate is high enough so that each oscillation period of $\left|r\right[n\left]\right|$ contains multiple samples, leading to multiple samples for each lobe of $p\left[n\right]$. Furthermore, each lobe of $p\left[n\right]$ should correspond to a single oscillation period of $\left|r\right[n\left]\right|$, thus there must be at least one sample below the threshold between oscillation periods of $\left|r\right[n\left]\right|$. Hence, the sampling rate $f_s$ must satisfy $f_s\gg 4f_0$, where $f_0$ is the carrier frequency. Nevertheless, if the sampling rate is small enough so that each lobe of $p\left[n\right]$ contains only one sample, $p\left[n\right]$ may be approximated by a mixture of Gaussian distributions with very small variances. It is worth noting that if both the sampling rate and the noise are small, TOAs may become deterministic, in which case any unbiased estimator produces the same outcome.

4. Estimating Unbiased TDOAs

In general, sensors receive hits generated by the same source with different amplitudes, since the distance of propagation is different for each sensor and waves attenuate as they propagate. Moreover, the attenuation may be frequency-dependent, so that the signal waveforms may differ at each sensor. Because of the difference of attenuation among sensors, the lag between the TOA estimated by the fixed threshold method and the actual time of arrival is not the same for all sensors (that is, $\mathbb{E}\{t_i-t_j\}\ne {\tau }_i-{\tau }_j$ even for noise variance tending to zero). As such, the fixed threshold method produces biased TDOA estimates, which will lead to biased position estimates. The fixed threshold method thus only produces unbiased estimates in scenarios with small attenuation difference between hits received by sensors.

We now derive a low-complexity method to obtain unbiased TDOAs in the noiseless case that consists of using different thresholds for each sensor, extending the fixed threshold method to a scenario where attenuation does not necessarily produce biased TDOA estimates. For this, we assume that waveforms at different sensors have the same shape with different amplitudes. In this model, the analog signals received by sensors $r_i\left(t\right)$ are attenuated and delayed versions of each other:

$$r_i\left(t\right)=a_i\psi (t-{\tau }_i),$$

(18)

$a_i$ is a known function of the distance between the source and the i-th sensor (denoted as $d_{i,s}$) that models the attenuation and $\psi \left(t\right)$ is an auxiliary function that controls the shape of the received signals. In Section 6, we use $a_i=\frac{e^{-\alpha d_{i,s}}}{\sqrt{d_{i,s}}}$, where the term $e^{-\alpha d_{i,s}}$ models the energy loss and the term $\frac{1}{\sqrt{d_{i,s}}}$ is due to the wave dispersion in the 2D medium [30]. Note that (18) assumes the frequency response of all sensors to be similar. Furthermore, if the structure generates reflections of the source signal, then (18) may be valid only for instants t smaller than the instant the first reflection arrives at the sensor.

Define $E_i=\sqrt{\frac{1}{M}{\sum }_n{r}_i^2\left[n\right]}$ as the RMS value of the hit corresponding to the i-th sensor, and ${\beta }_i=\frac{{\sum }_{k=1}^N{E}_k}{E_i}$. We propose to obtain TOAs in real time by comparing $r_i\left[n\right]$ with the modified threshold $K_i=\frac{K}{{\beta }_i}$ instead of using a fixed threshold. This way, the normalized threshold $K_i$ is used to estimate the TOAs given that a hit was detected.

We now show that if signals follow the propagation model (18), this procedure produces unbiased TDOAs in the noiseless case. Since $r_i\left[n\right]=r_i\left(nT\right)=a_i\psi (nT-{\tau }_i)$, we have $E_i=\sqrt{\frac{1}{M}{\sum }_n{a}_i^2{\psi }^2(nT-{\tau }_i)}=a_i{E}_0$, where $E_0=\sqrt{\frac{1}{M}{\sum }_n{\psi }^2(nT-{\tau }_i)}$. Hence, ${\beta }_i=\frac{{\sum }_{k=1}^N{a}_k{E}_0}{a_i{E}_0}=\frac{{\sum }_{k=1}^N{a}_k}{a_i}$.

The estimated TOA is $t_i=T{min}_n\{n:|r_i\left[n\right]|\ge \frac{K}{{\beta }_i}\}=T{min}_n\{n:{\beta }_i{a}_i|\psi (nT-{\tau }_i)|\ge K\}$, and substituting ${\beta }_i=\frac{{\sum }_{k=1}^N{a}_k}{a_i}$, we obtain:

$$\begin{array}{cc}\hfill t_i& =T\underset{n}{min}\{n:|\psi (nT-{\tau }_i)|\sum _{k=1}^N{a}_k\ge K\}\hfill \\ \hfill & =T\underset{n}{min}\{n:|\psi \left(nT\right)|\sum _{k=1}^N{a}_k\ge K\}-{\tau }_i.\hfill \end{array}$$

(19)

The term $\left|\psi \left(nT\right)\right|{\sum }_{k=1}^N{a}_k$ in (19) does not depend on i, thus the TDOA between hits measured by sensors i and j is $t_i-t_j={\tau }_i-{\tau }_j$, which is an unbiased TDOA measurement in the noiseless case.

Note that even though the waveform is required to compute the modified threshold $K_i$, the operations needed to calculate them have very low complexity and can be easily implemented in real time, thus waveforms do not need to be stored to apply this method. There are other TOA estimation methods that generate unbiased TDOAs, as methods based on Akaike information criterion (AIC) [8,31], but these methods present higher computational complexity than the fixed threshold method (which only needs comparisons) and our TOA estimation technique, whose complexity is dominated by the calculation of the energy of the signal.

This method only works in structures where waves propagate according to (18), thus it will not be effective if hits contain wave reflections. In this case, hits must be windowed to avoid adding the energy of reflections into the energy of the hit. The proposed TOA estimation method is employed in Section 5.1 to compare the source position estimators presented in this work. The performance of these estimator using the fixed threshold method and the proposed TOA estimator are also compared.

5. The Optimal Source Position Estimator

In this section, we derive an expression for the optimal cost function (that is, the one that yields the minimum mean-square localization error among all possible TOA-based estimators) assuming TOAs follow a Gaussian mixture model (GMM), which is a good approximation for high sampling rate as discussed in Section 3.3. We prove that $J_{\mathrm{TOA}}$ is the optimal cost function when the noise is small and independent and identically distributed (i.i.d.) in time and across sensors. In order to determine the optimal estimator, the Fisher information matrix must be calculated and compared to the inverse of the covariance matrix of our position estimator. It is worth noting that the Fisher information matrix for passive localization problem using Gaussian TOA or TDOA measurements (that is, a Gaussian mixture with $M=1$ components) was calculated in [32].

We consider that TOAs follow a mixture of M Gaussian distributions with weights $w_1,\cdots ,w_M$, covariance matrices ${\mathbf{C}}_1,\cdots ,{\mathbf{C}}_M$ and means ${\mathbf{v}}_1,\cdots ,{\mathbf{v}}_M$. In practice, all these parameters depend on the source position, but in order to derive the optimal estimator we assume that ${\mathbf{C}}_k$ and $w_k$ are independent of the source position, while ${\mathbf{v}}_k$ is modeled as ${\mathbf{v}}_k=\tau \left(\mathbf{x}\right)+{\mathbf{b}}_k\left(\mathbf{x}\right)+t\mathbf{1}$, where $\mathbf{b}\left(\mathbf{x}\right)$ is a TOA bias model that may depend on $\mathbf{x}$, $\mathbf{1}={\left[\begin{array}{ccc}1& \cdots & 1\end{array}\right]}^T$, $\tau \left(\mathbf{x}\right)={\left[\begin{array}{c}{\tau }_1\left(\mathbf{x}\right),\cdots ,{\tau }_N\left(\mathbf{x}\right)\end{array}\right]}^T$ and t is the instant of emission of the wave. Thus, defining the vector of measured TOAs $\mathbf{u}={\left[\begin{array}{c}{t}_1,t_2,\cdots ,t_N\end{array}\right]}^T$, the TOA pdf is given by:

$$f(\mathbf{u};\mathbf{x},t)=\sum _{k=1}^M\frac{w_k{e}^{-\frac{1}{2}{(\mathbf{u}-{\mathbf{v}}_k)}^T{\mathbf{C}}_k^{-1}(\mathbf{u}-{\mathbf{v}}_k)}}{\sqrt{{\left(2\pi \right)}^{N}det{\mathbf{C}}_k}}.$$

(20)

We also assume that for each $\mathbf{u}$ there is only one dominant term in the sum (20), that is, there is an index $m=m(\mathbf{u};\mathbf{x},t)$ such that for each set $\mathbf{u}$ of possible measured TOAs, the only term in (20) that is not approximately zero is the m-th component, so:

$$f(\mathbf{u};\mathbf{x},t)\approx \frac{w_m{e}^{-\frac{1}{2}{(\mathbf{u}-{\mathbf{v}}_m)}^T{\mathbf{C}}_m^{-1}(\mathbf{u}-{\mathbf{v}}_m)}}{\sqrt{{\left(2\pi \right)}^{N}det{\mathbf{C}}_m}}.$$

(21)

This is a reasonable hypothesis in the case where the sampling rate is such that there will be at least one sample below the threshold in (2) after each maximum, that is, the sampling rate should satisfy $f_s\gg 4f_0$, as in the example given in Section 3.3.

5.1. Optimal Source Position Estimator

Let ${\mathcal{S}}_1,{\mathcal{S}}_2,\cdots ,{\mathcal{S}}_M$ be a partition of ${\mathbb{R}}^N$ such that the k-th mixture component is different than zero for $\mathbf{u}\in {\mathcal{S}}_k$ and approximately zero for $\mathbf{u}\notin {\mathcal{S}}_k$. Each Gaussian mixture component lies in a single set ${\mathcal{S}}_k$, thus ${\mathcal{S}}_k$ can be viewed as the "support" of the k-th Gaussian component. The proposed optimal estimator defined below consists on estimating m such that $\mathbf{u}\in {\mathcal{S}}_m$ and estimating the source position with the conditional maximum likelihood estimator given that $\mathbf{u}\in {\mathcal{S}}_m$ (we show below how to estimate m even without knowing $\mathbf{x}$ and t).

The conditional probability distribution of $\mathbf{u}$ given that $\mathbf{u}\in {\mathcal{S}}_k$ is $f(\mathbf{u};\mathbf{x},t|\mathbf{u}\in {\mathcal{S}}_k)=\frac{1}{\sqrt{{\left(2\pi \right)}^{N}det{\mathbf{C}}_k}}{e}^{-\frac{1}{2}{(\mathbf{u}-{\mathbf{v}}_k)}^T{\mathbf{C}}_k^{-1}(\mathbf{u}-{\mathbf{v}}_k)}$, thus the conditional maximum likelihood estimator given that $\mathbf{u}\in {\mathcal{S}}_k$ is obtained by maximizing $log\left(f(\mathbf{u};\mathbf{x},t|\mathbf{u}\in {\mathcal{S}}_k)\right)$, or equivalently minimizing $J_k(\mathbf{x},t)=\frac{1}{2}{(\mathbf{u}-{\mathbf{v}}_k)}^T{\mathbf{C}}_k^{-1}(\mathbf{u}-{\mathbf{v}}_k)$. $J_k(\mathbf{x},t)$ can be rewritten in terms of only two variables by setting $\frac{\partial }{\partial t}{J}_k(\mathbf{x},t)=0$ and isolating t in terms of $\mathbf{x}$, yielding the equivalent cost function:

$$\begin{array}{cc}\hfill J_k\left(\mathbf{x}\right)& =\frac{1}{2}{(\mathbf{u}-{\mathbf{v}}_k)}^T{\mathbf{C}}_k^{-1}(\mathbf{u}-{\mathbf{v}}_k),\hfill \\ \hfill {\mathbf{v}}_k& =\tau \left(\mathbf{x}\right)+{\mathbf{b}}_k\left(\mathbf{x}\right)+t\left(\mathbf{x}\right)\mathbf{1},\hfill \\ \hfill t\left(\mathbf{x}\right)& =\frac{1}{N}{\mathbf{1}}^T\mathbf{u}-\frac{1}{N}{\mathbf{1}}^T({\mathbf{b}}_k\left(\mathbf{x}\right)+\tau \left(\mathbf{x}\right)).\hfill \end{array}$$

(22)

Assume we could choose $k=m^{*}\stackrel{\mathrm{def}}{=}m({\mathbf{x}}^{*},t^{*})$, that is, the index of the mixture component that most contributes to $f(\mathbf{u};{\mathbf{x}}^{*},t^{*})$, where ${\mathbf{x}}^{*}$, $t^{*}$ are the actual source position and instant of emission. Hence, recalling that $\mathbf{u}$ is the vector of measured TOAs,

$$\begin{array}{cc}\hfill J_{\mathrm{opt}}\left(\mathbf{x}\right)& =J_{m^{*}}\left(\mathbf{x}\right),\hfill \\ \hfill m^{*}& =arg\underset{k}{max}\frac{w_k{e}^{-\frac{1}{2}{(\mathbf{u}-{\mathbf{v}}_k^{*})}^T{\mathbf{C}}_k^{-1}(\mathbf{u}-{\mathbf{v}}_k^{*})}}{\sqrt{{\left(2\pi \right)}^{N}det{\mathbf{C}}_k}},\hfill \end{array}$$

(23)

where ${\mathbf{v}}_k^{*}$ is ${\mathbf{v}}_k$ computed at the actual source position ${\mathbf{x}}^{*}$, i.e., ${\mathbf{v}}_k^{*}=\mathbb{E}\left\{\mathbf{u}\right|\mathbf{u}\in {\mathcal{S}}_k\}=\tau \left({\mathbf{x}}^{*}\right)+{\mathbf{b}}_k+t^{*}\mathbf{1}$. Even though $({\mathbf{x}}^{*},t^{*})$ are unknown variables, in Section 5.3 we show how to obtain ${\mathbf{v}}_k^{*}$ from the TOA pmf, which is estimated using the noisy signals acquired by sensors. $m^{*}$ can be correctly chosen even if an imprecise value for ${\mathbf{v}}_k^{*}$ is used because $m^{*}$ is chosen by picking the maximum of a finite set, thus a perturbation on these values only affects $m^{*}$ if it is large enough to change the index of the maximum value.

We claim that the cost function $J_{\mathrm{opt}}\left(\mathbf{x}\right)$ yields the optimal estimator when minimized if the noise level is not very high and the Hessian of $J_{\mathrm{opt}}\left(\mathbf{x}\right)$ is approximately constant around the actual source position, in which case $J_{\mathrm{opt}}\left(\mathbf{x}\right)$ can be approximated by its second-order Taylor polynomial around ${\mathbf{x}}^{*}$. In order to prove this, we calculate the Fisher information matrix in Appendix A and the covariance matrix of our estimator in Appendix B. The covariance matrix (A11) coincides with the inverse of the submatrix associated with the spatial components of the Fisher information matrix (A2), which proves that the proposed estimator is indeed optimal.

It must be highlighted that we made the following hypotheses to compute the covariance matrix:

• The noise variance is not high, so that the Taylor approximation (A3) works. This is also necessary to guarantee that the estimator is approximately unbiased.
• The parameters of the TOA pdf $w_k,{\mathbf{b}}_k$ and ${\mathbf{C}}_k$ are known.

If these conditions do not hold, (23) may not be optimal. However, we show in Section 5.3 how to obtain an approximation to the optimal estimator when the assumptions are not true, and in Section 6 we show that the resulting estimator usually outperforms other estimators.

5.2. The Maximum Likelihood Estimator

The Maximum likelihood estimator (MLE) is obtained by finding the values of $\mathbf{x}$ that maximize (20). Given the definition of $m=m(\mathbf{x},t)$, this can be obtained by solving:

$$\widehat{\mathbf{x}},\widehat{t},m=arg\underset{\mathbf{x},t,k}{min}\frac{1}{2}{(\mathbf{u}-{\mathbf{v}}_k)}^T{\mathbf{C}}_k^{-1}(\mathbf{u}-{\mathbf{v}}_k)+{\kappa }_k,$$

(24)

where ${\kappa }_k=-log\left(\frac{w_k}{\sqrt{{\left(2\pi \right)}^{N}det{\mathbf{C}}_k}}\right)$. Since ${\kappa }_k$ does not depend on $\mathbf{x}$, the MLE coincides with the arguments that minimize $J_m\left(\mathbf{x}\right)$ for $m=arg{min}_k{min}_{\mathbf{x}}{J}_k\left(\mathbf{x}\right)+{\kappa }_k$.

Therefore, the optimal estimator may coincide with the MLE if they share the same value for m, but they are different estimators in general. However, their estimates always coincide if the noise is small enough so that the number of components in the Gaussian mixture is $M=1$, as in this case both estimators always pick $m=1$. Note that if $M=1$ and TOAs are unbiased and i.i.d, both the MLE and the optimal estimator coincide with $J_{\mathrm{TOA}}$.

5.3. The Gaussian Mixture TOA Estimator

In order to show that (23) is optimal, we assumed that the parameters of the Gaussian mixture are known and that m is such that $\mathbf{u}\in {\mathcal{S}}_m$, but these parameters must be estimated in practical problems. In this section, we describe a procedure to choose m and to estimate the mixture parameters using the noisy signals received by the sensors, leading to a suboptimal estimator we call the Gaussian mixture TOA estimator (GMTOA). To simplify the explanation of our method, we assume that the signals are corrupted by white noise and that the noise at different sensors is independent.

We estimate the TOA pmf for each sensor using (17), but instead of the noiseless waveform we compute $q_n$ using the noisy signal, that is,

$$p_i\left[n\right]=(1-F_W\left({\widehat{q}}_n\right))\prod _{k=0}^{n-1}{F}_W\left({\widehat{q}}_k\right),{\widehat{q}}_n=K-\left|r_i\left[n\right]\right|.$$

(25)

Note that evaluating this pmf does not depend on previous knowledge of the source position, which we use to extract the GMM parameters. Thus, we convert the pmf $p_i\left[n\right]$ to a pdf $p_i\left(t\right)=\frac{1}{T}{p}_i\left[\mathrm{round}(t/T)\right]$ and then we fit $p_i\left(t\right)$ into a Gaussian mixture model, as done in Section 3.3. Then, the means ${\mu }_{i,1},\cdots ,{\mu }_{i,M_i}$, variances ${\sigma }_{i,1}^2,\cdots ,{\sigma }_{i,M_i}^2$ and weights $w_{i,1},\cdots ,w_{i,M_i}$ are extracted for each mixture component $k=1,2,\cdots ,M$ and sensor $i=1,2,\cdots ,N$. Note that, since we are assuming the noise at different sensors to be independent, the TOAs in different sensors will be independent and the full mixture model will be the product of the distributions for each sensor, and have in total $M=M_1{M}_2\cdots M_N$ components.

Now, it is necessary to choose the index $m_i$ in the range $1\le m_i\le M_i$. Since the TOAs are independent, it is possible to choose the active index $m_i$ for each sensor independently. Index $m_i$ is chosen by finding the Gaussian component that most contributes to the fitted pdf:

$$m_i=arg\underset{k}{max}\frac{w_{i,k}}{\sqrt{2\pi {\sigma }_{i,k}^2}}{e}^{-\frac{{(u_i-{\mu }_{i,k})}^2}{2{\sigma }_{i,k}^2}}.$$

(26)

Finally, the TOA bias for the i-th sensor $b_i$ must be extracted from the mean of the chosen GMM component ${\mu }_{m_i}$. Ideally, ${\mu }_{m_i}$ and $b_{m_i}$ would be related by ${\mu }_{m_i}={\tau }_i\left({\mathbf{x}}^{*}\right)+t^{*}+b_{m_i}$, but ${\tau }_i\left({\mathbf{x}}^{*}\right)$ and $t^{*}$ are not known. As such, $b_{m_i}$ is estimated by approximating ${\tau }_i\left({\mathbf{x}}^{*}\right)+t^{*}$ as the first instant the estimated probability distribution $p_i\left[n\right]$ crosses a very small fixed threshold $K_p$, that is, the instant where it assumes a non-negligible value (in our simulations, we use $K_p={10}^{-3}$). Hence, $b_{m_i}$ is estimated as:

$$b_{m_i}={\mu }_{m_i}-\overline{n}T,\overline{n}=\underset{n}{min}\{n:p_i\left[n\right]\ge K_{\mathrm{pdf}}\}.$$

(27)

Therefore, the source position estimated by our proposed estimator is:

$$\widehat{\mathbf{x}}=arg\underset{\mathbf{x}}{min}{(\mathbf{u}-\mathbf{v})}^T{\mathbf{C}}^{-1}(\mathbf{u}-\mathbf{v}),$$

(28)

where $\mathbf{v}=\tau \left(\mathbf{x}\right)+\mathbf{b}+t\left(\mathbf{x}\right)\mathbf{1}$, $\mathbf{C}=\mathrm{diag}({\sigma }_{m_1}^2,{\sigma }_{m_2}^2,\cdots ,{\sigma }_{m_N}^2)$, $t\left(\mathbf{x}\right)=\frac{1}{N}{\mathbf{1}}^T\mathbf{u}-\frac{1}{N}{\mathbf{1}}^T(\mathbf{b}+\tau )$ and $\mathbf{b}={\left[\begin{array}{ccc}{b}_{m_1}& \cdots & b_{m_N}\end{array}\right]}^T$. The proposed algorithm is summarized in Algorithm 1.

Note that our method requires the signals $r_i\left[n\right]$ sampled by each sensor to localize the source, thus it does not use only TOAs as the cost functions $J_{\mathrm{TOA}}$, $J_{\mathrm{TDOA}}$ and $J_{CLS}$. Storing waveforms may be an issue for long acoustic emission tests that last weeks or months, as the waveforms may require an enormous storage capacity. However, it is possible to use our method without the need to store the whole waveforms because it uses $r_i\left[n\right]$ only to calculate $p_i\left[n\right]=(1-F_W\left({\widehat{q}}_n\right)){\prod }_{k=0}^{n-1}{F}_W\left({\widehat{q}}_k\right)$, which can be computed recursively: Defining $Q_i\left[n\right]={\prod }_{k=0}^{n-1}{F}_W\left({\widehat{q}}_k\right)$ and recalling that ${\widehat{q}}_n=K-\left|r_i\left[n\right]\right|$, we have:

$$\begin{array}{cc}\hfill p_i\left[n\right]& =(1-F_W\left({\widehat{q}}_n\right))Q_i\left[n\right],\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill Q_i\left[n\right]& =Q_i[n-1]F_W\left({\widehat{q}}_{n-1}\right).\hfill \end{array}$$

(30)

Since the computational cost of computing $F_W\left({\widehat{q}}_n\right)$ is small, it is possible to compute $p_i\left[n\right]$ in real time and discard the signals $r_1\left[n\right],\cdots ,r_N\left[n\right]$ instead of storing them. Storing $p_i\left[n\right]$ instead of the signals is much more efficient because $p_i\left[n\right]$ assumes a non-negligible value for only a small number of samples, thus it is possible to store only these non-negligible samples. Therefore, our method does not present the problem of storing a huge number of waveforms in long acoustic emission tests as methods that explicitly require the sampled signals as [33,34,35].

Algorithm 1: GMTOA source position estimator

6. Simulations

In this section, we assess the performance of our proposed source position estimator through simulations. Four sensors were positioned at the coordinates $(0,0),(0,1),(1,0)$ and $(1,1)$ in order to estimate the position of sources located inside the convex polygon whose vertices are the positions of the sensors (in this case, this polygon is a square) in several scenarios. In all cases, the cost functions were minimized using the mean position of the sensors $(0.5,0.5)$ as initial condition, except for the MLE which was implemented by finding the global maximum of the exact Gaussian mixture pdf (20) using M different initial conditions, each one obtained by minimizing $J_k\left(\mathbf{x}\right)$ for $k=1,2,\cdots ,M$.

6.1. Unbiased i.i.d. Gaussian TOAs, Different Source Positions

In order to test the optimality of our algorithm in ideal conditions, we generated TOAs as Gaussian random variables with no bias and standard deviation of $1\mathsf{\mu }$s and $10\mathsf{\mu }$s. Denoting the source position as ${\mathbf{x}}^{*}={\left[\begin{array}{cc}x& y\end{array}\right]}^T$, y was fixed at $y=0.4\mathrm{m}$ and x was swept from 0 to 1. For each position, the Mean Square Error (MSE) of the estimators was calculated using ${10}^5$ samples and used to compute the efficiency of the estimator, defined as $\eta =\frac{\mathrm{tr}\left\{\mathcal{I}{\left({\mathbf{x}}^{*}\right)}^{-1}\right\}}{\mathrm{tr}\left\{\mathbb{E}\left\{(\widehat{\mathbf{x}}-{\mathbf{x}}^{*}){(\widehat{\mathbf{x}}-{\mathbf{x}}^{*})}^T\right\}\right\}}$. Note that $\eta =1$ if the estimator is optimal.

The efficiency in terms of the source position is shown in Figure 3. For $\sigma =1\mathsf{\mu }\mathrm{s}$, the efficiency of the cost function $J_{\mathrm{TOA}}$ is very close to $100\%$, while $J_{\mathrm{TDOA}}$ and $J_{CLS}$ show worse performance for all source positions, indicating that $J_{\mathrm{TOA}}$ is indeed the optimal estimator for unbiased i.i.d. Gaussian TOAs.

On the other hand, for $\sigma =10\mathsf{\mu }\mathrm{s}$, the estimates obtained through $J_{\mathrm{TOA}}$ are very close to optimal when the source lies inside the region of interest (the convex polygon whose vertexes are the positions of the sensors). However, outside this polygon, $J_{\mathrm{TOA}}$ is no longer the optimal cost function because the approximations used to prove the estimator is optimal are not precise in this region. This happens because the Taylor approximation (A3) does not hold when ${\mathbf{x}}^{*}$ is in a region where the Hessian of $J_{\mathrm{TOA}}\left(\mathbf{x}\right)$ is not approximately constant, which may happen if the source is close to a sensor.

Even though $J_{\mathrm{CLS}}$ has worse accuracy than $J_{\mathrm{TOA}}$ in most cases, its efficiency does not fall down abruptly for x outside the region of interest, especially in the case of high noise level, in which the performance of the other methods is close to zero outside the region between sensors.

It is worth noting that in practice the sensors are scattered around the monitored region, thus the precision of localization algorithms is more important for sources lying inside the region of interest, in which case $J_{\mathrm{TOA}}$ performs better than the other methods for both noise variances.

6.2. Fixed Source Position and Different SNRs

Next, we performed simulations estimating the TOAs from noisy signals, as happens in practice. Note that:

• For low SNR, the Taylor approximation (A3) (Appendix B) does not hold and the estimator may be biased.
• The MSE of GMTOA depends on the quality of the estimated GMM parameters, being optimal only if these estimates match the actual ones.
• The GMM assumption for TOAs is an approximation, so an estimator may show lower MSE than the derived CRLB when empirical TOAs are used.

We have done four simulations in similar conditions as before, except that the source position is fixed at ${\mathbf{x}}^{*}={\left[\begin{array}{cc}0.3& 0.4\end{array}\right]}^T$. In all simulations, hits were generated according to (18) using $a_i=\frac{e^{-\alpha d_{i,s}}}{\sqrt{d_{i,s}}}$ and $\psi \left(nT\right)={sin}^2\left(\frac{\pi n}{1000}\right)sin\left(\frac{2\pi n}{67}\right)$ (a modulated von Hann window), with white Gaussian noise of different variances. The sampling rate is $10\mathrm{Msps}$, hence the frequency of the carrier is $\frac{{10}^7}{67}\approx 150\mathrm{kHz}$, the resonance frequency for some acoustic emission sensors.

In these simulations, TOAs were calculated using the TOA estimation technique proposed in Section 4 with a fixed threshold $K=0.3$. The GMTOA estimator was implemented using Algorithm 1, and the parameters obtained by this algorithm were used in the MLE.

6.2.1. TOAs Following a GMM and Known Parameters

In this simulation, the noisy hits were used to obtain auxiliary TOAs through the fixed threshold method. Then, the mean and variance of these auxiliary TOAs were computed and used to generate TOAs following a Gaussian mixture distribution with known parameters. Note that ${\mathbf{v}}_k$ is not known because it depends on the tentative source position.

The MSE in terms of noise standard deviation for the presented methods is shown in Figure 4a,b. The MSE of the GMTOA estimator is very close to the CRLB for all noise variances, confirming our derivations. Moreover, the maximum likelihood estimator does not coincide with the optimal one if $\sigma$ is high, but the MLE is optimal for low noise level.

6.2.2. Empirical TOAs and Unknown Parameters

In this simulation, the actual TOAs obtained from the fixed threshold method are used to localize the source. Hence, TOAs are not exactly distributed according to a Gaussian mixture in this case. The Cramér–Rao lower bound was calculated based on the parameters estimated by the GMTOA Estimator from noisy hits using Algorithm 1.

Figure 4c,d shows that if the noise variance is very high, the GMTOA estimator does not perform better than the others, but it does have a significantly better performance for intermediate variances where the TOA pdfs cannot be approximated as single Gaussian distributions.

On the other hand, the proposed estimator has lower MSE than the other estimators for very low variances, but it may perform worse than them if the variance is not very low, but low enough so that the TOA pdf is nearly Gaussian. This is because the MSE of $J_{\mathrm{TOA}}$, $J_{\mathrm{TDOA}}$ and $J_{CLS}$ fall abruptly when the TOA distribution becomes approximately Gaussian (or a mixture with only one component), while the MLE and the GMTOA estimator do not because they depend on the noisy parameters of the estimated Gaussian pdf. For smaller variances, the parameters from the pdf approximate well the actual ones, thus the GMTOA estimator becomes again better than $J_{\mathrm{TOA}}$ and $J_{\mathrm{TDOA}}$.

Note that the MSEs of the estimators fall below the Cramér–Rao lower bound if $\sigma$ is small. This happens because we derived the CRLB for TOAs following a known Gaussian mixture distribution, but the true distribution is discrete, not exactly a mixture of Gaussians.

6.2.3. TOAs Obtained from Filtered Hits and Unknown Parameters

The main issue with the maximum likelihood and GMTOA estimators is that they depend on the noise pdf, which is not accessible in practice and must be estimated. If the noise variance is large, these parameters can be poorly estimated and cause performance loss in the localization methods, as verified in the last simulation. To test the robustness of these estimators in a scenario of non-white noise, a lowpass filter is applied to hits before calculating the threshold in this simulation. The filter used was a Butterworth filter of order 15 and digital cutoff frequency ${\omega }_c=0.1\pi$ (or 1 MHz), but the TOA pmf was calculated assuming white Gaussian noise.

The results are shown in Figure 4e,f. With the filter, the parameters of the GMM are better estimated, and the GMTOA estimator has better performance than all other methods even for high variance. The MLE also presents much lower MSE than $J_{\mathrm{TOA}}$ and $J_{\mathrm{TDOA}}$ for high noise level.

Since in this simulation GMTOA assumes white noise to calculate the TOA pmf (which is used to extract the GMM parameters), it may have worse performance using filtered hits (which are corrupted by non-white noise) instead of the original hits for some values of $\sigma$, but in most cases using filtered hits implied in lower MSE.

6.2.4. Using Biased TDOAs

In all previous simulations, TOAs were obtained using a proposed modification of the fixed threshold method described in Section 4, which generates unbiased TDOAs. Figure 5 shows a simulation with the same setup as in Section 6.2.2 and Section 6.2.3, but the fixed threshold method was used to obtain TOAs instead of the proposed modification. As such, TDOAs are biased in this simulation.

Comparing Figure 4c,e with Figure 5, the proposed TOA estimation method significantly decreases the MSE of TOA-based localization methods. These figures also show that the GMTOA estimator achieves better performance than other methods in most cases for TOAs obtained from original hits and in all cases if hits are filtered before TOA extraction.

It is worth noting that if biased TOAs are used directly for estimation, the MSE may decrease for higher noise variance because bias may be reduced.

7. Conclusions

In this work, we derived the expression for the TOA pmf obtained through the threshold method, and showed that it is not Gaussian-distributed as usually assumed. We show that the TOA pmf can only be approximated as a Gaussian distribution if the sampling rate is high and the noise level is small. The distribution can be modeled by a mixture of Gaussian distributions for higher noise level. We also presented a TOA estimation method based on the fixed threshold algorithm that generates unbiased TDOAs in structures where hits are approximately delayed and attenuated versions of each other.

The optimal source position estimator was derived assuming TOAs follow a mixture of Gaussian distributions of known parameters. This estimator coincides with the maximum likelihood estimator for small noise level, in which case TOAs are Gaussian. Since the mixture parameters are not known in practical applications, we proposed a method to estimate them from noisy waveforms, resulting in the GMTOA estimator, a suboptimal method to localize the source. We verified in simulation that GMTOA yields better results than other methods for both biased and unbiased TDOAs in most cases, even if the TOA pmf is not exactly a mixture of Gaussian distributions, and its performance is enhanced if a low-pass filter is applied to hits before extracting TOAs.

Even though the proposed estimator performed better than other methods in simulation, it is necessary to further validate it with data collected from an acoustic emission test. Further work also includes extending our TOA estimation method to generate unbiased TDOAs for more complex propagation models, in particular for structures where hits contain many reflections of the wave emitted by the source.