Optimal Array Design and Directive Sensors for Guided Waves DoA Estimation

Abstract

The estimation of Direction of Arrival (DoA) of guided ultrasonic waves is an important task in many Structural Health Monitoring (SHM) applications. The aim is to locate sources of elastic waves which can be generated by impacts or defects in the inspected structures. In this paper, the array geometry and the shape of the piezo-sensors are designed to optimize the DoA estimation on a pre-defined angular sector, from acquisitions affected by noise and interference. In the proposed approach, the DoA of a wave generated by a single source is considered as a random variable that is uniformly distributed in a given range. The wave velocity is assumed to be unknown and the DoA estimation is performed by measuring the Differences in Time of Arrival (DToAs) of wavefronts impinging on the sensors. The optimization procedure of sensors positioning is based on the computation of the DoA and wave velocity parameters Cramér-Rao Matrix Bound (CRMB) with a Bayesian approach. An efficient DoA estimator is found based on the DToAs Gauss-Markov estimator for a three sensors array. Moreover, a novel directive sensor for guided waves is introduced to cancel out undesired Acoustic Sources impinging from DoAs out of the given angles range. Numerical results show the capability to filter directional interference of the novel sensor and a considerably improved DoA estimation performance provided by the optimized sensor cluster in the pre-defined angular sector, as compared to conventional approaches.

1. Introduction

Acoustic source (AS) localization from the measurements of passive sensors is a widely-investigated problem in structural health monitoring (SHM) [1,2]. An AS could consist of undesired impacts on the monitored structure. Alternatively, Acoustic Emissions (AE) can be generated by the growth of defects, such as cracks and corrosion. Acoustic sources in plate-like structures generate Lamb waves [3], i.e., guided stress waves (GW). The detection by means of piezoelectric sensors and the subsequent analysis of AS signals allow to infer whether an impact has occurred and to localize it. Although several methods for AS localization were proposed and validated at laboratory scale, these approaches rarely satisfy the following stringent constraints which are posed by field applications:

• Minimal monitoring system weight, including cabling and circuitry;
• Minimal power consumption, to be compatible with wireless systems. This means reducing the computational cost for the signals local processing and/or the amount of data to be wirelessly transferred;

Among all methods of AS localization, three main typologies can be identified:

• Inverse methods [4,5];
• Hyperbolic positioning [6];
• Methods which use Directions of Arrival (DoAs) estimation [7,8,9];

The first typology requires very accurate modeling of ultrasonic propagation. It can be performed by knowing previously the structural parameters or by measuring them. However, the main limitations of inverse methods, are that they require a high computational cost and, consequently, are incompatible with lightweight and low-power embedded systems. The second typology estimates the AS location as the intersection point of geometric hyperbola derived by the time delays between sensors. The hyperbolic positioning methods require a large number of distant and well-synchronized sensors in order to reduce the estimation uncertainty. Unfortunately, accurate synchronization is difficult to implement in a wireless system. For these reasons, the methods of the third typology, which estimate DoAs by using clusters of closely located sensors, are often the only feasible approach to AS localization. Indeed, these methods have a computational cost compatible with embedded systems and do not require synchronization between distant sensors.

Different strategies which use one or more clusters for the AS estimation have been studied and tested [7,8,9]. Among these, the first two [7,8] use just one cluster to locate an AS. Nevertheless, in order to estimate the AS distance, they require to detect two modes (typically the A0 and S0 modes) and compute their time difference on the sensors. Conversely, the technique proposed by Kundu [9] is suitable even when one only wave mode is detected. This is the case of AS generated by an impact, when, typically, the Lamb waves A0 fundamental mode has a much higher amplitude w.r.t. other modes. Such technique uses–at least–two closely spaced sensors clusters placed apart to locate sources both in isotropic and slightly anisotropic plates. The DoA estimation performance of a simple cluster of three circular sensors placed on the vertices of an isosceles right triangle, with unknown wave velocity, was investigated in [10,11]. The DoA estimation is performed by means of the Differences in Time of Arrival (DToAs) estimations, via simple Cross-Correlation (CC) procedures. In this work, this cluster will be referred to as Standard Cluster (SC). In [12,13], it was shown that it is possible to use multiple clusters to estimate the AS localization even in the case of heavily anisotropic structures (i.e., characterized by elliptical or rhomboidal wavefronts).

Although the SC using, for AS location, was well-validated at laboratory scale, it is well-known that in realistic field deployment, several problems arise due to noise and directional interference. Typical noisy physical sources are given by: structural vibration (e.g., due to turbulence on an aircraft), scattered wavefield, noisy acquisition channels, or noisy electronic devices. Whereas the directive interference is due to the waves reverberation within bounded structures. In particular, we can distinguish coherent interference w.r.t the signal to be detected and incoherent signals. The first ones are due to edge-reflections of the impact/defect AE to be detected. The second ones are due to reflections produced by different impacts or acoustic events.

Several studies have been already presented to address the noise and the reverberation issues in DoA estimation and AS localization. In [14], Oktel and Moses proposed a sensors cluster design procedure to increase the DoA estimation performance in presence of noise. It was based on the Bayesian approach (or global) of the Cramér Rao bound (CRB), which depends on the sensors positioning and defines the lower bound of any unbiased estimator. However, in this work, the wave velocity is supposed to be known. Conversely, in many applications, such an assumption is not verified and results in a loss in accuracy.

Regarding the directional interference issue, several estimation methods were proposed to filter or distinguish the undesired corruption from the useful signal. In particular, among the DoA estimators based on the DToAs estimations, the Generalized CC (GCC) procedure with phase transform (GCC-PHAT), has been shown to be a suitable alternative in reverberant scenarios [15,16]. However, its performance can only be considered optimal for a high signal-to-noise ratio (SNR) (as shown in [17]). As a consequence, GCC-PHAT is not suited for the applications considered in this work.

A different DoA estimation method, the Multiple Signal Classification (MUSIC) [18], is able to estimate up to N-1 DoAs due to different sources with N-sensor arrays. Originally designed to estimate the number and DoAs of uncorrelated signals, modified versions [19,20] have been proposed to estimate also the DoAs of coherent signals for the multipath environment. However, with a simple 3-sensor cluster, just 2 coherent signals can be detected. This means that other directional interferences may cause wrong estimations. A more robust MUSIC algorithm for reverberant scenarios is proposed in [21]. MUSIC algorithms are also limited by the assumption of accurate knowledge or estimation of wave velocity. Therefore, an additional iterative wave velocity estimation procedure is needed (as shown in [22]), which increases consistently the computational cost.

This work provides original solutions to tackle the detrimental effects both of noise and of directional interference. The first problem is addressed by means of a novel strategy for sensor cluster design. Unlike [14], we considered the velocity to be unknown. In our approach, the cluster design procedure is based on the computation of the CRB of DoA in case of the unknown velocity of propagation (CRB${}_{u-v}$), and on the usage of its Bayesian average (assuming that the DoA is a random variable with a known probability distribution) as a cost function for the optimal design. This allows us to minimize both the DoA lower bound and the DoA accuracy loss due to the unknown wave velocity.

The second problem, i.e., the negative impact of directional interference, is tackled by means of novel directive piezo-sensors, suitable for guided propagation structures. The shaping of piezoelectric transducers has already been used as a powerful means to detect DoA in non-reverberant contexts [23,24]. Here, we show how the shaping can be used to filter out all directional interference, either coherent or incoherent w.r.t. the useful wave signal. It is worth noting that the usage of these transducers is beneficial whenever a limited angular sector has to be monitored, without regard to the number of sensors and the adopted DoA estimator.

The design procedure for the directive transducer proposed in this paper draws inspiration from the work of Senesi and Ruzzene [25] which showed how to relate the transducer shape to its directivity. However, the transducers considered in [25] are characterized by symmetric beam patterns which are unsuited to distinguish sources related to opposite directions. To achieve this capability, a novel complex (i.e., multi-phase) transducer has to be implemented [26]. The proposed Directive Complex sensor (DCS) consists of five piezoelectric patches and allows to suppress the lobes out of a 90° monitoring sector. The experimental validation of the piezo-patches shaping design is beyond the scope of this paper. Nevertheless, future developments are aimed to realize the proposed DCS, by exploiting the available manufacturing techniques. For example, a laser-cut can be used to shape piezo-electrodes on the upper surface of metalized PVDF (polyvinylidene fluoride) sheets, as proposed in [27]. Alternatively, a metallic printing technique can be used on PVDF films to obtain the desired shape patches [28].

The details of the cluster design strategy, the adopted DoA estimator, the novel directional transducer concept and, finally, the numerical validations are thoroughly illustrated in the following sections.

2. System Model and Cramér-Rao Matrix

Let us assume that the sensors array consists of three identical sensors: P1, P2, and P3. The sensors are located at ${\mathbf{r}}_i={[x_i,y_i]}^T$ for $i\in [1,2,3]$. Following [6,29], we adopt a model with a single co-planar far-field source which generates the wavefield impinging the 3 sensors array. The signal at the ith sensor is $s(t-d_i)$, where $s\left(t\right)$ is the signal at a reference point near the array and $d_i$ is the delay at the ith sensor w.r.t. the reference point. Without loss of generality, the reference point is assumed to be coincident with the location of the first element in the array, P1, so that $d_1=0$. We assume also that the sensors are near enough so that the amplitude gradient across the array and the effect of wave dispersion are negligible. The output signal of the ith sensor can be expressed as:

$$x_i\left(t\right)=s(t-d_i)+n_i\left(t\right)$$

(1)

where $n_i\left(t\right)$ is the additive sensor noise at ith sensor. In order to estimate the DoA of wavefront, we first estimate the vector of DToAs, $\mathbf{d}={[d_2,d_3]}^T$. In the discrete Fourier domain, the 3 × 1 measurements vector at kth frequency ${\omega }_k$, is given by:

$$\mathbf{x}\left({\omega }_k\right)={\mathbf{a}}_{\theta }\left({\omega }_k\right)s\left({\omega }_k\right)+\mathbf{n}\left({\omega }_k\right)$$

(2)

where ${\mathbf{a}}_{\theta }\left({\omega }_k\right)$ is the steering vector, defined as:

$${\mathbf{a}}_{\theta }\left({\omega }_k\right)={\left[1,e^{j{\omega }_k{d}_2\left(\theta \right)},e^{j{\omega }_k{d}_3\left(\theta \right)}\right]}^T$$

(3)

where $d_i\left(\theta \right)=({\mathbf{u}}^T\left(\theta \right)·{\mathbf{r}}_i-{\mathbf{u}}^T\left(\theta \right)·{\mathbf{r}}_1)/v$ is the DToA between the ith sensor and the reference sensor P1, v is the wave velocity and $\mathbf{u}\left(\theta \right)={\left[cos\left(\theta \right),sin\left(\theta \right)\right]}^T$ is the unit vector pointing toward the signal source. The followings hypotheses are formulated:

• The noises are stationary Gaussian processes with zero mean. The signal is a Gaussian process with zero means, approximately stationary. The last hypothesis assumes that, for closely spaced sensors, the dispersion effect can be neglected;
• The signal and noises are mutually uncorrelated and uncorrelated between themselves;

Under the previous hypotheses, Hahn and Tretter in [29], derived the Cramér-Rao Matrix Bound (CRMB), $\mathbf{Q}$, for the delays. We suppose that an optimal estimator is used to estimate the DToAs, such as the Maximum Likelihood (ML) estimator proposed in [29], therefore the covariance matrix is equal to the CRMB. Moreover, we suppose that noises have identical covariance matrix. This means that the noises have identical spectrum and in case of white noises, they have the same noise level. When the noises have same spectrum, the covariance matrix $\mathbf{Q}$ assumes the following simple form:

$$\mathbf{Q}={\sigma }_{d_i}^2\left[\begin{array}{cc}1\hfill & \hfill 1/2\\ 1/2\hfill & \hfill 1\end{array}\right]$$

(4)

where the variance of time-delays, ${\sigma }_{d_i}^2$, is given by:

$${\sigma }_{d_i}^2=\frac{2}{3}\frac{2\pi }{T_S{\displaystyle \underset{0}{\overset{B_S}{\int }}}2{\omega }^2\frac{S^2/N^2}{1+3(S/N)}d\omega }$$

(5)

where $S\left(\omega \right)$ and $N\left(\omega \right)$ are the power spectra of signal and noise, $B_S$ is the signal band and $T_S$ is the signal time duration. Thanks to the asymptotically Gaussian property of a ML estimator, the conditional probability density function of $\mathbf{d}$ is:

$$f(\mathbf{d};\theta )=\frac{1}{{\left(2\pi \right)}^2\sqrt{det\mathbf{Q}}}exp\left\{-\frac{1}{2}{\left(\mathbf{d}-\frac{\mathbf{r}}{v}\right)}^T{\mathbf{Q}}^{-1}\left(\mathbf{d}-\frac{\mathbf{r}}{v}\right)\right\}$$

(6)

with $\mathbf{r}={\left[r_2\left(\theta \right),r_3\left(\theta \right)\right]}^T$ the vector of Differences in Distance of Arrival (DDoAs) of the wavefront between the sensors and the reference sensor P1. Note that $r_i\left(\theta \right)$ ($i=2,3$) are so that the “true” values of $d_i\left(\theta \right)$ are $d_i\left(\theta \right)=r_i\left(\theta \right)/v={({\mathbf{u}}^T\left(\theta \right)·{\mathbf{r}}_i-{\mathbf{u}}^T\left(\theta \right)·{\mathbf{r}}_1)}_i/v$. Let’s assume that the wave velocity is unknown. So, the unknown parameters are $\theta$ and v. In [30], Malagò and Pistone provided the Fisher information Matrix (FIM) for a Gaussian distribution when the vector of the means $\mathbf{\mu }$ and the covariance matrix $\mathbf{Q}$ are both functions of a set of parameters $\mathbf{\gamma }={({\gamma }_1,{\gamma }_2,...,{\gamma }_K)}^T$. The expression, specialized in our case (with ${\gamma }_1=\theta$ and ${\gamma }_2=v$), is the following:

$$I_{m,n}\left(\mathbf{\gamma }\right)=\frac{\partial \left[r_2\left(\theta \right)/v,r_3\left(\theta \right)/v\right]}{\partial {\gamma }_m}{\mathbf{Q}}^{-1}\frac{\partial \left[\begin{array}{c}{r}_2\left(\theta \right)/v\hfill \\ r_3\left(\theta \right)/v\hfill \end{array}\right]}{\partial {\gamma }_m}$$

(7)

with ${\mathbf{Q}}^{-1}$, from (4) is given by:

$${\mathbf{Q}}^{-1}=\frac{4}{3{\sigma }_{d_i}^2}\left[\begin{array}{cc}1\hfill & -1/2\\ -1/2\hfill & 1\end{array}\right]$$

(8)

The inverse of the FIM provided in Equation (7) is the sought CRMB for unknown parameters $\theta$ and v, which defines a lower bound of the covariance matrix for any unbiased estimator of two parameters.

3. Cost Function and Array Geometry Design

The D-criterion (see [31]) uses the determinant of CRMB (equal to the inverse of FIM determinant, often called generalized covariance bound) as a cost function to be minimized to obtain the optimal design. However, in the current application domain, only the DoA estimation performance has to be optimized. Therefore, the following function has been considered:

$$\frac{det(\mathrm{FIM})}{I_{22}}=I_{11}-\frac{I_{12}{I}_{21}}{I_{22}}$$

(9)

It consists of two terms: the first one is related only to the Fisher Information (FI) on $\theta$ (when v is known) and the second one is related to information on $\theta$ and v, when they are simultaneously estimated, divided by the FI on v (when $\theta$ is known). Finally, adopting a Bayesian (or global) approach, similarly to [14], the following CRB${}_{v-u}$ cost function is defined:

$$\begin{array}{c}{J}_C\left(\mathbf{r}\right)=\mathrm{E}\left[I_{22}det(\mathrm{CRMB})\right]=\mathrm{E}\left[\frac{I_{22}}{det(\mathrm{FIM})}\right]\\ =\frac{1}{2{\theta }_0}{\displaystyle \underset{-{\theta }_0}{\overset{{\theta }_0}{\int }}}\frac{I_{22}\left(\mathbf{r}\left(\theta \right)\right)}{det\left(\mathrm{FIM}\right(\mathbf{r}\left(\theta \right)\left)\right)}f\left(\theta \right)d\theta \end{array}$$

(10)

where $f\left(\theta \right)$ is the probability density function (pdf) of $\theta$, thought as random variable and $[-{\theta }_0,{\theta }_0]$ is its domain, supposed compact. The column vector $\mathbf{r}\left(\theta \right)$ is function of sensors locations coordinates as:

$$\mathbf{r}\left(\theta \right)=\left[\begin{array}{c}(x_2-x_1)cos\left(\theta \right)+(y_2-y_1)sin\left(\theta \right)\hfill \\ (x_3-x_1)cos\left(\theta \right)+(y_3-y_1)sin\left(\theta \right)\hfill \end{array}\right]$$

(11)

We define the CRB${}_{u-v}$-optimal array ${\mathbf{r}}_C$ as the one whose elements location is given by:

$$\begin{array}{c}{\mathbf{r}}_C=argmin{J}_C\left(\mathbf{r}\right)\hfill \\ \mathrm{with}\left\{\begin{array}{c}\sqrt{x_1^2+y_1^2}\le d\hfill \\ \sqrt{x_2^2+y_2^2}\le d\hfill \\ \sqrt{x_3^2+y_3^2}\le d\hfill \end{array}\right.\hfill \end{array}$$

(12)

where d is the radius of the circular domain where the sensor elements are constrained to lie in.

The general problem statement can be specified to the case of uniform pdf in a 90° sector, i.e., $\left[-\pi /4,\pi /4\right]$. Computing the terms $I_{m,n}$ of (7), the cost function (10) is:

$$J_C\left(\mathbf{r}\right)={\sigma }_{d_i}^2\frac{2v^2}{\pi }{\displaystyle \underset{-\pi /4}{\overset{\pi /4}{\int }}}\frac{r_2^2\left(\theta \right)+r_3^2\left(\theta \right)-r_2\left(\theta \right)r_3\left(\theta \right)}{{\left({r_2}^{\prime }\left(\theta \right)r_3\left(\theta \right)-r_2\left(\theta \right){r_3}^{\prime }\left(\theta \right)\right)}^2}d\theta$$

(13)

Its minimization is obtained with a symmetric configuration with respect to 90°-axis (i.e., Y-axis), and a half-opening angle $\beta$ equal to 23° (see Figure 1b). Conversely, when the wave velocity is known, the optimized array would be symmetric w.r.t the X-axis [14]. This particular configuration is due to the minimization of the DoA accuracy loss due to the unknown wave velocity [32], via Equation (13). The obtained 3-sensor cluster will be referred to as Designed Cluster (DC) while the configuration of Figure 1a, which is related to the cluster proposed by Kundu et al. in [10,11], will be referred to as Standard Cluster (SC) in this paper.

Finally, it is worth noting that the previous design procedure of the optimal sensor positioning is still valid for a generic number M of sensors, with an appropriate covariance matrix (4).

4. A DoA Efficient Estimator

As anticipated, an efficient time delays estimator has to be used in order to match the DToAs covariance matrix with the CRMB (4). The Maximum Likelihood (ML) DToAs estimator, asymptotically efficient, was proposed by Hahn and Tretter in [29]. The technique consists of measuring the DToAs for all possible sensors pairs by Generalized Cross Correlation (GCC) and then calculating the Gauss-Markov (GM) (weighted) estimate of the DToAs with respect to the first sensor. The GCC procedure consists of computing the Cross-Correlation between the acquired signals, filtered first by an appropriate filter. The Optimal Filtering to attain the time delays CRMB is defined by:

$${\left|F_{OPT\left(\omega \right)}\right|}^2=\frac{S\left(\omega \right)/N^2\left(\omega \right)}{1+M\left(S\right(\omega )/N(\omega \left)\right)}$$

(14)

where M is the number of the sensor, whereas S and N are the power spectra of respectively no-noisy signal and noise. In practice, the optimal filter requires knowledge or estimation of the signal and noise spectra. A simple estimation method consists in measuring the noise spectrum and computing $S\left(\omega \right)$ by subtracting the noise spectrum from the noisy signal spectrum. However, due to random variations of noise, spectral subtraction can result in negative estimates of the short-time magnitude or power spectrum. Different methods for reducing and removing the distortions due to the rectification process are proposed in [33].

Under the hypothesis that the noises have the same spectrum for each sensors, the Gauss-Markov estimator coefficients, for the case of three sensors array (regardless of what filter is used for the GCC procedures), are given by:

$$\left[\begin{array}{c}{d}_{2GM}\hfill \\ d_{3GM}\hfill \end{array}\right]=\frac{2}{3}\left[\begin{array}{cc}1 1/2\hfill & -1/2\hfill \\ 1/2 1\hfill & 1/2\hfill \end{array}\right]\left[\begin{array}{c}{d}_{12GCC}\hfill \\ d_{13GCC}\hfill \\ d_{23GCC}\hfill \end{array}\right]$$

(15)

where $d_{ijGCC}$ are the DToAs between sensor i and j estimated by using the GCC procedure, whereas $d_{iGM}$ are the time delays with respect to the reference sensor estimated with the GM estimator. The weights of the Gauss-Markov estimator have a more complex form (expressed by a ratio of time delays variances) only when the noise spectra for each sensor, $N_i$, are not all equal.

In this paper, for testing the DoA estimation performance with the two clusters of Figure 1, we assumed equal white noise spectra and a flat signal spectrum within a Band $B_s$, to emulate the narrow-band impulsive signals due to an impact. In this case, the Optimum Filter (14) is equal to an arbitrary constant within the Band $B_s$, and 0 elsewhere:

$${\left|F_{OP{T}_{Flat}}\left(\omega \right)\right|}^2=\left\{\begin{array}{c}1\hspace{1em}\omega \in B_s\hfill \\ 0\hspace{1em}\mathrm{elsewhere}\hfill \end{array}\right.$$

(16)

The Band $B_s$, is estimated by using the spectral subtraction technique, assuming to know the white noise level. The distortions induced by the rectification of negative values of the estimated power spectrum S are neglected. This assumption is justified when SNR values are sufficiently high, while the non-linear distortions are not negligible when the signal-to-noise ratio decreases.

Finally, an optimal DoA estimation function from the estimated time delays has to be found to attain the CRB${}_{u-v}$ (the inverse of (9)). Given the designed array geometry (DC in Figure 1b), the following result is obtained:

$$d_2=\frac{d}{v}2cos\left(\alpha \right)cos\left(\theta \right),\hspace{1em}{d}_3=\frac{d}{v}\left(cos\left(\alpha \right)cos\left(\theta \right)-(1+sin(\alpha \left)\right)sin\left(\theta \right)\right)$$

(17)

$\theta$ can be computed by using the ratio $\eta \stackrel{.}{=}{d}_2/d_3$ and inverting with respect to $\theta$. The “true” relation between $\theta$ and $\eta$ can be used as estimation function:

$${\widehat{\theta }}_{DC}=\mathrm{atan}\left(\frac{tan\beta (\eta -2)}{\eta }\right),\hspace{1em}\mathrm{with}\beta =23°$$

(18)

Note that the estimator ${\widehat{\theta }}_{DC}$ is a function only the ratio $d_2/d_3$. From the Theory of Uncertainty Propagation [34], the mean square error and variance (when the time delays $d_2,d_3$ are random variables with covariance matrix (4)) can be computed by expanding the estimation function in the first order Taylor series:

$$e_{ms}={\sigma }_{{\widehat{\theta }}_{DC}}^2=\mathrm{E}[\widehat{\theta }]-\theta =\left[{\left(\frac{\partial f}{\partial d_2}\right)}^2+{\left(\frac{\partial f}{\partial d_3}\right)}^2+\left(\frac{\partial f}{\partial d_2}\right)\left(\frac{\partial f}{\partial d_3}\right)\right]{\sigma }_{d_i}^2$$

(19)

this formula is valid around the point of expansion $(\mathrm{E}\left[d_2\right],\mathrm{E}\left[d_3\right])$ of function f in Taylor series. Expanding Equation (19), we have:

$$\begin{array}{c}{\sigma }_{{\widehat{\theta }}_{DC}}^2={\left(\frac{\partial f}{\partial \eta }\right)}^2\left[\frac{d_2^2+d_3^2-d_2{d}_3}{d_3^4}\right]{\sigma }_{d_i}^2\\ =\frac{d_2^2\left(\theta \right)+d_3^2\left(\theta \right)-d_2\left(\theta \right)d_3\left(\theta \right)}{{\left({d_2}^{\prime }\left(\theta \right)d_3\left(\theta \right)-d_2\left(\theta \right){d_3}^{\prime }\left(\theta \right)\right)}^2}{\sigma }_{d_i}^2\end{array}$$

(20)

where the last equivalence is provided by the inverse function theorem. The function in (20) is precisely the integrand of (13), i.e., $I_{22}/det(\mathrm{FIM})$, expressed in terms of DToAs $d_i\left(\theta \right)$, instead of elements $r_i\left(\theta \right)$, namely the Differences in Distance of Arrival (DDoAs).

5. Numerical Performance Comparison between Clusters and DoAs Estimators

In order to validate the design procedure of the array geometry and the performance of the proposed DoA estimator, a numerical analysis was performed. The estimation functions of the two clusters of Figure 1 are respectively:

$${\widehat{\theta }}_{SC}=\mathrm{atan}\left(\frac{d_2}{d_3}\right),\hspace{1em}{\widehat{\theta }}_{DC}=\mathrm{atan}\left(\frac{tan(23°)(d_2/d_3-2)}{d_2/d_3}\right)$$

(21)

In particular, impact waves propagating in aluminium plate 1 [mm] thick (Young’s modulus 70 [GPa], Poisson’s coefficient 0.3 and material density 2700 [kg/m${}^3$]) were simulated with the Greens function formalism adopted in [35]. The impulse response of a bandpass Batterworth filter (10th order) with different bandwidths and center frequencies was used in order to simulate the impact signal. The two DToAs with respect to the first sensor $d_2,d_3$ were computed from the simulated acquired signals by using three different estimation modalities:

• locating the peaks of the cross-correlation ([36]);
• by using the Gauss-Markov estimator (15) consisting in three cross-correlation procedures;
• by combining the Gauss-Markov estimator (15) with three GCC procedures (filtering first the signals with the filter (16);

The last modality is the optimum one in the considered case which involves (additive) white (zero-mean Gaussian) noise (AWGN) and quasi-flat signal spectrum S in a fixed band. The power spectrum $S\left(\omega \right)$ and its band $B_S$ for the filter (16) are estimated by the spectral subtraction technique.

Simulations were performed for multiple impact locations obtained by varying the true DoA with $5^{°}$ steps ($\theta =-45°,-40°,-35°,\dots ,45°$), the distance from the cluster center being 0.8 [m]. The results achieved by the clusters of Figure 1 are given in

Table 1. Comparison of Standard Deviation (SD) and Maximum Error (ME) values (in [Deg]) between the arrays of Figure 1 for noise-affected measurements at different PSNR values. Actuated pulse band: [30–40] [kHz].

	Standard Cluster						Designed Cluster
	CC		GM-CC		GM-GCC		CC		GM-CC		GM-GCC
PSNR	SD	ME	SD	ME	SD	ME	SD	ME	SD	ME	SD	ME
60 dB	0.45	1.05	0.44	1.05	0.43	1.05	0.35	0.67	0.29	0.63	0.29	0.63
40 dB	0.94	3.14	0.80	3.30	0.46	1.86	0.84	2.94	0.59	1.84	0.32	1.16
35 dB	1.53	4.99	1.27	5.05	0.52	1.86	1.46	6.1	0.98	3.6	0.37	1.19
30 dB	2.50	8.13	2.14	7.50	0.67	2.58	2.35	8.8	1.62	5.82	0.51	1.92
28 dB	3.09	10.10	2.55	10.19	0.79	3.48	2.86	11.09	1.97	7.73	0.59	2.25
27 dB	3.50	12.41	2.90	12.46	1.80	8.31	3.24	10.82	2.25	7.79	1.73	6.02
26 dB	3.85	12.30	3.18	11.20	2.74	10.47	3.68	12.78	2.52	9.45	2.09	7.88
24 dB	4.78	15.30	3.98	14.44	3.64	12.85	4.57	20.67	3.1	11.17	2.83	9.87

,

Table 2. Comparison of Standard Deviation (SD) and Maximum Error (ME) values (in degrees) between the arrays of Figure 1 for noise-affected measurements at different PSNR values. Actuated pulse band: [50–60] [kHz].

	Standard Cluster						Designed Cluster
	CC		GM-CC		GM-GCC		CC		GM-CC		GM-GCC
PSNR	SD	ME	SD	ME	SD	ME	SD	ME	SD	ME	SD	ME
60 dB	0.58	1.05	0.52	1.05	0.52	1.05	0.44	0.92	0.28	0.83	0.26	0.83
40 dB	1.22	4.09	1.02	4.03	0.55	1.97	1.1	4.31	0.77	2.49	0.37	1.26
35 dB	1.93	6.99	1.61	6.95	0.63	2.91	1.83	8.7	1.26	4.31	0.44	2.06
30 dB	3.31	11.22	2.72	11.05	0.84	3.03	3.03	13.41	2.11	8.58	0.62	2.31
28 dB	4.14	14.2	3.41	13.05	1	4.88	3.76	14.42	2.58	9.6	0.93	3.51
27 dB	4.52	16.04	3.72	14.96	2.01	8.77	4.26	15.55	2.93	12.2	2.52	9.18
26 dB	5.08	16.14	4.23	15.77	3.92	15.97	4.79	16.89	3.23	12.34	2.86	9.79
24 dB	6.36	21.69	5.28	20.28	4.86	20.9	5.84	20.7	4.01	15.1	3.78	12.78

and

Table 3. Comparison of Standard Deviation (SD) and Maximum Error (ME) values (in degrees) between the arrays of Figure 1 for noise-affected measurements at different PSNR values. Actuated pulse band: [30–60] [kHz].

	Standard Cluster						Designed Cluster
	CC		GM-CC		GM-GCC		CC		GM-CC		GM-GCC
PSNR	SD	ME	SD	ME	SD	ME	SD	ME	SD	ME	SD	ME
60 dB	0.35	0.85	0.34	0.85	0.34	0.85	0.45	0.72	0.4	0.72	0.4	0.72
40 dB	0.48	1.87	0.45	1.87	0.37	1.42	0.46	1.62	0.37	1.08	0.38	0.72
35 dB	0.73	2.28	0.64	2.3	0.41	1.87	0.63	2.31	0.46	1.76	0.37	1.08
30 dB	1.14	4.05	0.96	4.16	0.48	1.87	1.04	4.03	0.72	3.11	0.38	1.08
28 dB	1.35	5.1	1.15	4.83	0.49	1.87	1.25	4.69	0.87	3.15	0.39	1.12
27 dB	1.5	5.07	1.25	5.2	0.54	2.14	1.39	4.92	0.96	4.35	0.51	2.08
26 dB	1.67	5.53	1.4	5.46	1.1	4.53	1.54	6.08	1.07	3.75	0.82	3.75
24 dB	2.05	6.61	1.69	6.61	1.59	6.21	1.88	7.09	1.3	4.67	1.2	4.57

, for different center frequencies and bandwidths of the impact signal and different peak signal to noise ratios (PSNR). To assess the Standard Deviation (SD) of DoA estimations, 200 simulations, on the entire 90° sector, were performed. Furthermore, the maximum error (ME) over all simulations was considered. The simulations were run simulating the propagation of the A0 Lamb mode, and considering: (i) circular piezo sensors with radius equal to 5 [mm], (ii) the radius of the array d equal to 2 [cm], (iii) sampling frequency (Fs) equal to 2 [MHz].

As shown by the Table 1, Table 2 and Table 3, the SD and the ME values obtained with the designed cluster are, as expected, smaller w.r.t. the SC, for PSNR values higher than 24 dB. In particular, the best performances are achieved when the optimal DoA estimator, based on the GM-GCC time delays estimator, is used. In this case, the variances almost equate to the CRB${}_{u-v}$. Due to the wave dispersion of the A0 mode, the higher is the considered center frequency, the higher is the wave (group) velocity v. It is important to consider this fact because the DoA CRB in Equation (13) and the variance of the DoA estimator in Equation (20) increase as $v^2$. Conversely, for the case of quasi-flat signal spectrum in a given band $B_S$, the DToA variance term ${\sigma }_d^2$ of the Covariance Matrix (the CRMB in Equation (4) for an optimal estimator) decreases as $B_S^{3/4}+3B_S^{}{\omega }_c^2$, where ${\omega }_c$ is the central frequency of $B_S$ (see (5)). Therefore, the two terms, wave velocity and central frequency have an opposite influence on the DoA estimation performance (see the SD values of Table 1 and Table 2). Finally, it can be noted that the higher is the bandwidth, the smaller is the DToA and DoA variance (compare Table 1, Table 2 and Table 3 which are characterized by a 10 and 30 [kHz] Bandwidth, respectively).

6. Directive Base Sensor Design

In order to filter ASs with a DoAs out of considered angles range, a novel directive base sensor is investigated in this paper. The sensor beam pattern is ideally equal to 1 in a given range and 0 elsewhere. Without lack of generality, we refer to the [0°,90°] range as the one where the beampattern is equal to 1. The beampattern of a sensor is linked to its shape as described by the model proposed in [37] which estimates the frequency response of a piezo sensor impinged by a Lamb wave mode as:

$$V_P\left(\omega \right)=jU\left(\omega \right)k_0\left(\omega \right)H_P\left(\theta \right)D_P(\omega ,\theta )$$

(22)

in this equation, $U\left(\omega \right)$ denotes the amplitude and the polarization of the wave component relevant to the piezo properties of the patch at the considered frequency, $k_0\left(\omega \right)$ is the wave vector that characterizes the propagation and $H_P\left(\theta \right)$ is a quantity related to the material properties of the piezo-structure system. Without lack of generality, here we consider the case of piezo patches with a single polarization and constant piezoelectric properties. Finally, the only function which depends by the DoA $\theta$ and by the piezo-sensor shape is the $D_P(\omega ,\theta )$ function. It defines the frequency response for all possible angles of arrival $\theta$. Therefore, it is called Directivity function and can be computed by the following integral:

$$D_P(\omega ,\theta )={\displaystyle \underset{{\mathsf{\Omega }}_P}{\int }}{e}^{-j{k}_0\left(\omega \right)(xcos\theta +ysin\theta )}{\varphi }_P(x,y)dxdy$$

(23)

where ${\varphi }_P(x,y)$ is the function that describes the geometry of sensor and is referred to as shape function. Defining ${\mathsf{\Omega }}_P$ the area of the piezoelectric path, the shape function is equal to 1 when $(x,y)\in {\mathsf{\Omega }}_P$ and 0 elsewhere.

Considering as an example a circular piezo-sensor of radius R, the (23) provides:

$$D_P(\omega ,\theta )=2\pi R^2\frac{J_1\left(R{k}_0\left(\omega \right)\right)}{R{k}_0\left(\omega \right)}\approx 2\pi R^2\mathrm{sinc}\left(R{k}_0\left(\omega \right)\right)$$

(24)

where $J_1(·)$ is the first kind Bessel function of first order and $k_0\left(\omega \right)$ is the wave vector of the propagation mode of Lamb waves (e.g., A0 or S0 mode). Observe that the (24) doesn’t depend by $\theta$ but only by frequency, so the directive properties of a piezo-disk are the same for all angles, i.e., the disk is omnidirectional. We define the base sensor beampattern at frequency $\omega$ as;

$$d(\omega ,\theta )=\frac{\left|D_P(\omega ,\theta )\right|}{{max}_{\theta }\left[\left|D_P(\omega ,\theta )\right|\right]}$$

(25)

For a circular sensor then $d(\omega ,\theta )=1$. Thus, the Directivity function $D_P(\omega ,\theta )$ (23) is equal to the bi-dimensional spatial Fourier Transform (2D-FT) at angle $\theta$ of the shape function. Then, the shape function which corresponds to a given desired directivity can be determined with an Inverse Fourier Transform. In our approach, we impose the same Directivity function (and so the same frequency response) of a piezo-disk of given radius in the [0–90]° angles-sector and 0 elsewhere. Therefore, we compute the 2D-FT of a disk, set to 0 all values of 2D-FT in the k-space domain between [90;360]°, and finally get back in the space-domain, via the 2D-Inverse FT (IFT) to obtain the desired shape function (see Figure 2).

Note that a real shape function is achieved if and only if the Directivity function is symmetric with respect to the origin. Instead, our procedure tolerates the generation of a complex shape function, with a real part and an imaginary part, both having positive and negative values (Figure 2c,d). This higher complexity allows us to have a beam pattern that is not symmetrical, i.e., without lobes in the [180:270]° range. Unfortunately, the described procedure produces a continuously modulated shape function that cannot be realized in practice. So, a quantization procedure is applied to the computed shape function. In particular, the phase of the complex shape function is quantized as detailed in the

Table 4. Phase quantization scheme used in the complex shape function implementation.

Phase Interval	Quantized Value	Patch Color in Figure 3a
$[-\pi /4+\pi /4)$	0	Yellow
$[\pi /4+3/4\pi )$	$\pi /2$	Red
$[3/4\pi +5/4\pi )$	$\pi$	Blue
$[5/4\pi +7/4\pi )$	$3/2\pi$	Black

.

Then, the absolute values greater than a certain positive threshold are set to 1 and others to 0 (see Figure 3). Areas associated with the same quantized values define the shape of the electrodes of the piezo patches used as sensors. A distance gap of at least 0.5 [mm] between the patches has been imposed to be compliant with the typical geometrical limitations that are associated with patch manufacturing.

Such a procedure generates the Directive Complex Sensor (DCS). As can be seen in Figure 3a, such a sensor consists of five piezo patches each one depicted with a different color. It is worth noting that the two blue patches related to the quantized phase $\pi$ can be short-circuited. Moreover, the three patches corresponding to the opposite phases 0 and $\pi$ (referred to as the real part) correspond to regions where the computed shape function has almost equal absolute average value. The same applies for the two patches related to phases $\pi /2$ and $3/2\pi$ (the imaginary part). This implies that piezo-patches related to the real part require just one differential acquisition channel and a second differential channel is required by the patches related to the imaginary part. In order to generate the (complex) time-signal, a weighted sum of the two acquired differential signals has to be performed, in which the signal related to the imaginary part is multiplied by the factor $i·W_{Im}$, where i is the imaginary unit and $W_{Im}$ is a suitable weight. Finally, the anti-analytic signal of the complex acquired signal is computed and used to feed the DToA estimator.

Figure 3b shows the actual 2D-FT (absolute value) in the k-space domain, i.e., the 2D-FT computed after the quantization of the shape function. Due to quantization, the values of the 2D-FT, out of [0°,90°], are not 0, but are still smaller than values in the monitored angular sector.

By using Equation (25), the DCS sensor theoretical beampatterns were computed at different values of frequency, considering the wave vector values of A0 mode when propagating in an aluminum plate with 1 [mm] thickness (i.e., for a known dispersion curve $k_0\left(\omega \right)$). As shown in Figure 4, useful directive beampatterns are achieved in the [10–60] [kHz] frequencies band. However, the best directional behavior is achieved in the [30–60] [kHz] frequency range.

The directivity properties at each frequency can be expressed by the Average Directional Attenuation (ADA) parameter. It is computed by the beampatterns values, as the ratio between the average beampattern value in the monitored angular sector (i.e., [0–90]°) and out of that one. In Figure 4, the ADA values for different frequency values are shown. The ADA value is above 8.3 dB in the [10–60] [kHz] band, and 13.0 dB in the [30–60] [kHz] band.

Regarding the $W_{Im}$ parameter, in order to find an optimal value a suitable cost function $J_{W_{\mathrm{Im}}}$ was defined:

$$J_{W_{\mathrm{Im}}}\left(W\right)=\frac{1}{N}{\displaystyle \sum _{i=1}^N}\frac{\varsigma (W,k_i)}{D{A}_{\alpha }^{\beta }(W,k_i)}$$

(26)

where $k_i$ are N wave-numbers in the considered (spatial) spectral region, $D{A}_{\phi }^{\gamma }(W,k_i)$ is the directional attenuation at angle $\phi$ w.r.t the one at angle $\gamma$ and $\varsigma (W,k_i)$ is equal to highest sidelobe level of beampatterns at each frequency $k_i$. By minimizing $J_{W_{\mathrm{Im}}}$ for $k_i=353,409,459,505 [\mathrm{r}\mathrm{a}\mathrm{d}/\mathrm{m}]$ (corresponding to the frequencies $f_i=30,40,50,60 \left[\mathrm{k}\mathrm{H}\mathrm{z}\right]$ of the A0 mode in the considered setup plate), $\phi$ = 120° and $\gamma$ = 90°, the value $W_{\mathrm{Im}-Opt}=6.81$ is obtained. It is worth noting that the DA of the proposed DCS is sufficiently high to mask directional interferences in a given angular range $\left[{\psi }_1,{\psi }_2\right]$, subset of [90°,360°]. For example, the DA is above 11.3 dB within [30–60] [kHz], for all undesired DoAs within [130–320] deg. This is due to DCS beampattern non-idealities, in particular to the non-sharp mainlobe cut-off near 0° and 90°. The DCS DA is also limited by the highest sidelobe level. The previous facts justified the cost function (26).

Such non-idealities can be attributed to the detrimental effect of binary quantization. More specifically, the DCS shows the better directivity properties in the wave vector k values range [353–505] [rad/m] (corresponding to the 30–60 [kHz] beampatterns shown in Figure 4). It is worth noting that the relationship between wave vector function and frequency $k_0\left(\omega \right)$ depends on the monitored structure characteristics (material, thickness, etc.). In other words, the optimal frequency range of the DCS can be found by taking into account the dispersion curves of the monitored medium.

7. Numerical Results of the Directive Complex Sensor

In order to validate the design procedure and the directivity properties of the DCS, in the following subsections, the beampatterns obtained from the Finite Element Method (FEM) are shown and compared with the theoretical ones (see Figure 4). Furthermore, a numerical comparison between the performance of DC (rotated by 45° to work within [0–90]°) of conventional disk-sensors and the same cluster with DCS (see Figure 5) is given, considering undesired AS with DoAs out of [0–90]°.

7.1. Finite Element Simulation Using COMSOL Multiphysics

To validate the DSC performance, the theoretical beampatterns predicted by the model have been compared with the ones resulting from finite element (FE) simulations. To do so, a three-dimensional COMSOL^® [38] based FE model of the proposed DSC was built. In the numerical model, an aluminum plate with dimensions of 500 [mm] × 60 [mm] and thickness 1 [mm] was chosen as the propagation medium. Since it is sufficient to shape just one metalization of the DCS (top or bottom) to achieve the desired directive behavior, the DCS was modeled using the geometry obtained in the design procedure, below which a small disk of piezoelectric material with a radius of 12 [mm] was defined. Then, the DCS was attached to the plate as shown in Figure 6.

The excitation for the A0 mode was simulated using a line load in a way that a plane wave is generated within the plate. It should be noted that the excitation signal is considered as a sine-wave with a combination of four different frequencies of 30, 40, 50 and 60 kHz. Unlike the common procedure to compute the directivity pattern, which includes a number of point sources around a fixed transducer, a different approach was utilized here: at each simulation run, the sensor was rotated of 5 degrees, while the excitation load was fixed, as depicted in Figure 6. The motivation of using such a method was to reduce the computational cost of the numerical model. Furthermore, in order to prevent wave back-reflection at side boundaries, the Low-Reflecting Boundary option was utilized. Two physics, Structural Mechanics and Electrostatics were coupled by the Multiphysics-Piezoelectric Effect to take the solid mechanics of the aluminum plate and the electrical feature of the piezoelectric sensor into consideration. The simulation results including the sensor response and the generated wavefield at different times for $\theta =0$ are given in Figure 7 and Figure 8, respectively. The beam patterns obtained from the FE simulation are compared to that of the theoretical model in Figure 9. A notably good agreement between them is observed, indicating the effectiveness of the proposed complex sensor.

7.2. DoA Estimation with DCS Clusters in Reverberant Environments

As discussed in the Introduction, in realistic reverberant environments, coherent reflections and incoherent reflections may hamper the DoA estimation. This interference can be viewed as waves generated by virtual Image Sources (ISs), due to the mirroring produced by the boundaries of the monitored structure [39]. In the most general case, undesired AE is given by multiple ISs. Let us suppose that the AS signal is corrupted by undesired reflections due to an edge closely spaced w.r.t. the sensor cluster, as shown in Figure 10. In this example, the AS is placed in the position specified by the blue circle (DoA equal to 90°), while Edge 1 generates an IS ($I{S}_1$) which, in turn, generates coherent interference on the signal acquired by the DCS cluster.

Additional incoherent interference may be generated by ISs of previous acoustic events. Both coherent and incoherent directional interference could have a detrimental effect on DoA estimation. Improved DoA estimation performance is achieved when the DoAs of ISs occur in angular ranges $\left[{\psi }_1,{\psi }_2\right]$ filtered by the beampattern of the DCS, which ensures a minimum DA level.

In order to evaluate the DoA estimation performance of DCS clusters in realistic simulation setups, the cases of coherent and incoherent reflection interferences are considered. For both cases, in the [0–90]° angular-sector, the wave to be detected were simulated by changing their orientation with a step of 5°.

At first, coherent interference was simulated, the cluster was placed at $d_c$ = 17 [cm] from the edge, whereas the AS location distance was set equal to 40 [cm]. The directional interference $I{S}_1$ is produced by the mirroring of the AS induced by edge reflections (see the Figure 10: the AS DoA is equal to 90°, whereas the corresponding $I{S}_1$ DoA is equal to 130.36°). Considering a sampling frequency $f_s$ equal to 2 [MHz], a 200 samples Tukey window (i.e., a rectangular window with the first and last 47.5 percent of the samples equal to parts of a cosine) filtered by using a Butterworth filter (10th order) with a bandpass equal to [30–40] [kHz] was used as impact signal and as IS.

Considering the Design Array configuration depicted in Figure 5, the DoAs estimated by the processing of the simulate response (via the GM-GCC time delays estimator) of piezo disk-sensors and DCSs, together the actual AS DoA and the corresponding IS DoA for 19 different simulated angles cases are reported in

Table 5. Comparison of DoA estimation performance between the two arrays depicted in Figure 5 for measurements affected by a coherent edge reflection due to an Image Source (see Figure 10). In the table, the nominal DoA value, the direction of the interferring source, and the estimated values are reported in degrees (ASs band: [30–40] [kHz]; Impact distance d = 40 [cm]; Cluster distance from the edge $d_c$ = 17 [cm]).

		Disks	DCSs			Disks	DCSs
		GM-GCC	GM-GCC			GM-GCC	GM-GCC
AS DoA	IS DoA	Estimated DoA	Estimated DoA	AS DoA	IS DoA	Estimated DoA	Estimated DoA
0	180	2.71	−0.93	50	160.27	60.14	49.85
5	178.24	3.93	3.93	55	157.7	46.95	56.98
10	176.48	11.01	9.19	60	154.91	53.05	58.91
15	174.68	16.14	14.94	65	151.87	66.43	66.18
20	172.86	22.81	19.58	70	148.51	72.26	73.16
25	170.98	27.23	25.52	75	144.77	97.54	74.07
30	169.04	28.96	30.34	80	140.56	91.82	82.38
35	167.02	30.71	36.31	85	135.8	14.23	83.89
40	164.9	37.06	40.15	90	130.36	104.67	93.33
45	162.66	47.31	45.95	-	-	SD 18.48	SD 1.47
						ME 70.76	ME 3.33

. In these conditions, the Standard Deviation and the Maximum Error values are equal to $18.48 deg$ and $70.76 deg$ for the piezo-disk cluster, and $1.47 deg$ and $3.33 deg$ for the DCSs cluster, respectively.

Two examples of acquired time signals, distinguishing the signal component related to the wave to be detected and the reflection, on a piezo Disk and on a DCS for the same DoAs are illustrated in Figure 11 (more specifically, for the DCS, the anti-analytic real part of the complex time signal is plotted). As can be seen, the AS signal and the IS signal are overlapped in time, hence a very unfavorable condition for DoA estimation. However, the DCS clearly shows the capability to strongly attenuate the spurious component.

In order to assess the DoA estimation performance in even more challenging conditions, the case of measurements affected both by directional interference and diffuse noise (AWGN) was considered. 200 simulations of AWGN, on the entire 90° sector, were performed for different PSNR values. The Standard Deviation and the Maximum values are given in

Table 6. DoA estimation performance (SD and ME in degress) for the DCS cluster shown in Figure 5) when the measurementes are affected by a coherent edge reflection (simulation setup of Table 5) and diffuse noise (AWGN).

PSNR	60 dB	30 dB	20 dB	15 dB	10 dB	9 dB	8 dB	7 dB
SD	1.5	1.52	1.65	1.9	2.47	2.63	2.85	3.11
ME	3.33	4.58	5.69	6.69	9.24	9.57	11.91	12.99

.

Then, we considered the case of an additional incoherent interference with random DoA in the angular-sector ([169.69–180]°) (to simulate an undesired incoherent component due to an $I{S}_2$ of a previous impact/defect).

In particular, we have simulated AS and incoherent spurious waves impinging on the sensors simultaneously, because these are the most critical conditions to perform the DoA estimation. The impulse response of a [30–40] [kHz] band Chebyshev Type I filter (10th order with a passband ripple of 4 dB) with a 3 dB amplification (to simulate a high energy impact) was used as actuating signal of the incoherent component.

The DoA estimation results in the presence both of a coherent and incoherent interference by means of the cluster of disk sensors and DCSs (Figure 5), for 19 different simulated cases are reported in

Table 7. Comparison of DoA (in degrees) estimation performance between two designed arrays (Figure 5) for measurements affected by a coherent edge reflection due to an IS (see the $I{S}_1$ in the Figure 10), and a incoherent spurious signal due to a second IS of a previous impact/defect with a random DoA within the range [169–180] (see the $I{S}_2$ in the Figure 10). AS and ISs band: [30–40] [kHz]; (Impact distance d = 40 [cm]; Cluster distance from the edge $d_c$ = 17 [cm]; Structure length $l_s$ = 1.1 [m]).

			Disks	DCSs				Disks	DCSs
			GM-GCC	GM-GCC				GM-GCC	GM-GCC
AS DoA	IS DoA	IS-2 DoA	Estimated DoA	Estimated DoA	AS DoA	IS DoA	IS-2 DOA	Estimated DoA	Estimated DoA
0	180	178	16.66	−2.03	50	160.27	171	69.3	49.85
5	178.24	179	7.62	5.42	55	157.7	179	28.44	58.49
10	176.48	171	17.76	10.07	60	154.91	179	116.24	58.91
15	174.68	179	32.01	15.4	65	151.87	175	20.26	68
20	172.86	176	29.71	20.23	70	148.51	178	83.38	73.16
25	170.98	171	30.2	25.38	75	144.77	171	83.52	73.62
30	169.04	172	23.41	28.9	80	140.56	174	27.38	81.15
35	167.02	175	25.85	34.42	85	135.8	179	17.65	79.78
40	164.9	179	61.1	40.15	90	130.36	178	110.89	92.33
45	162.66	179	92.82	46.87	-	-		SD 31.49 ME 67.35	SD 2.03 ME 5.22

. The Standard Deviation and the Maximum Error values are equal to $31.4 deg$ and $67.3 deg$ for the piezo-disk cluster, and $2.03 deg$ and $5.22 deg$ for the DCSs cluster, respectively. Three examples of acquired time-signals, distinguishing the signal component related to the wave to be detected and the spurious ones, on a piezo Disk and on a DCS for the same DoAs are given in Figure 12. Also, in this case, the spurious components are strongly attenuated by the directional sensor.

Finally, the DoA estimation performance was evaluated when the measurements are affected also by diffuse noise (AWGN). 200 simulations of AWGN were performed for different PSNR values. The Standard Deviation and the Maximum values, shown in

Table 8. DoA estimation performance (SD and ME in degress) by means of a DCS designed cluster (Figure 5)) when the measurementes are affected by a coherent edge reflection (simulation setup of Table 5) and diffuse noise (AWGN).

PSNR	60 dB	30 dB	20 dB	15 dB	10 dB	9 dB	8 dB	7 dB
SD	2.04	2.04	2.13	2.53	3.38	3.54	4.26	4.04
ME	5.22	5.67	7.14	9.28	13.15	16.07	16.04	18.05

, clearly validate the ability of the DCS to cancel out multiple spurious interferences.

It is worth noting that the analyzed configuration is representative of many realistic scenarios. In all the considered cases, there is a clear advantage in the performance achieved by the DCS cluster with respect to the conventional piezo disks. The performance of the DCS is, however, degraded when the ISs and the AS to be detected are generated at closely spaced locations. This is due to the fact that the directional selectivity of the DCS is not perfect, and, in case of interferences whose direction of arrival is slightly larger than 90° or slightly less than 0°, the attenuation is poor.

It is also worth noting that the proposed DCS cluster provides good results in DoA estimation both for coherent and incoherent directional interference even for low PSNR values because the DoA estimator (21), based on the GM-GCC estimator, is the optimal estimator in presence of noise, for any SNR value. Viceversa, the GCC-PHAT signal processing is useful just for high SNR values because is based uniquely on the signal phase information. Furthermore, the DCS sensors allow to work on all signal time-lapse, whereas the commonly-adopted selection of smaller time windows reduces the DoA estimation accuracy, in particular when the non-impulsive and long-lasting signal is to be detected.

8. Conclusions

In this work, an optimal N-sensor array design procedure for DoA estimation is discussed. Specifically, the minimization of the CRB of the DoA with a Bayesian approach is used as an optimality criterion, considering the wave velocity to be unknown.

The general procedure was applied to design a 3-sensor cluster to monitor a 90° sector. Two of such clusters can be used to detect and locate Acoustic Sources such as crack growth emissions or impacts over a rectangular area.

An efficient DoA estimator was found, based on the GCC-Gauss Markov estimator of the DToAs. It was shown that the optimal GCC is a band-pass filter, for the case of white noise and narrow-band quasi-flat signal spectrum whose band can be estimated via the subtraction of noise spectrum. The Designed Cluster and DoA estimator were validated by numerical simulations.

Despite the good performances achieved, it must be considered that, in realistic scenarios, many physical sources can generate directional interference. This is the case of waves due to reverberation, for example. To filter out this interference, a novel directive passive sensor was designed to replace the piezo-disks which are conventionally adopted in this application field. The novel sensor exploits its shape as a means to attenuate spurious waves coming from directions that are not included in the monitored angular sector. It consists of five-piezo patches whose output is collected by two differential channels. The new sensor is able to filter directional interference in a known k wave vector bandwidth ([353–505] [rad/m]). In this band, the average attenuation for spurious waves is 12.89 dB. Future developments are aimed to improve the quantization procedures to achieve even better directivity in a larger bandwidth and to provide experimental results of the proposed DCS concept in different practical scenarios.