Spiral-shaped piezoelectric sensors for Lamb waves direction of arrival (DoA) estimation

Here we focus our attention on methods to estimate the wave direction of arrival (DoA) which are based on the measurement of the difference in time of arrival (DToA) among couples of piezo patches forming a sensor. Assuming a planar wavefront, the DToA (Δt(θ)) for a couple of ideal point-like piezo patches both working as sensors, namely P1 and P2 (see figure 1(a)), is related to the wave DoA denoted by the angle θ by this formula:

$$\begin{eqnarray}&&{\rm{\Delta }}t(\theta )=\displaystyle \frac{d\cdot \cos \theta }{v}\end{eqnarray} \tag{ 1 }$$

where d is the distance among the patches and v is wave velocity³ . Consequently, we can evaluate the DoA as:

$$\begin{eqnarray}&&{\theta }_{{est}}=\arccos \displaystyle \frac{v\cdot {\rm{\Delta }}{t}_{{est}}}{d}\end{eqnarray} \tag{ 2 }$$

where Δt_est is the estimation of the actual DToA computed by processing the measured signals at the two patches. It is interesting to investigate how the errors in the measurements of the DToA affect the estimation of the DoA. This can be done by exploiting the Propagation of Uncertainty framework proposed in [22, 23]. In figure 1(b), the upper and lower bounds of the worst-case error in estimating the actual DoA (θ) are represented. In particular, the Matlab tool by Ridder [23] was exploited for the computation of such quantity for two values of d, considering an uncertainty in Δt_est equal to just 2 μs and an arbitrary wave velocity equal to v = 2000 m s⁻¹.

Zoom In Zoom Out Reset image size

It can be seen in figure 1(b) that the error in the estimation can be very large when θ ≤ 40 deg. The error can be reduced by enlarging the spacing between the two transducers d. However, a larger d implies a quadratically larger Fraunhofer distance, i.e. the source location distance beyond which the plane wave approximation can be considered valid.

We can remark here that the configuration in figure 1(a) does not allow to discriminate between DoA with positive or negative incoming angle θ as DToA(θ) = DToA(−θ), whereas opposite directions of propagation can be discriminated by looking at the sign of the DToA, which can be easily detected by the sign of the signals cross-correlation peak, as DToA(θ) = −DToA(θ + 180 deg).

Let us suppose now to reshape the piezo patches so that there is a linear dependence between the DToA and the DoA:

$$\begin{eqnarray}&&{\rm{\Delta }}t(\theta )=\displaystyle \frac{\alpha \cdot \theta }{v}+{\rm{\Delta }}t(0)\end{eqnarray} \tag{ 3 }$$

where α is an arbitrary constant. Such result can be achieved with a very good approximation if one of the piezo patches (P2 in this case) is shaped like a segment of an Archimedean Spiral, as schematically depicted in figure 2(a).

Zoom In Zoom Out Reset image size

The spiral segment can be expressed in polar coordinates (ρ_s, θ_s) as ρ_s = aθ_s + b, where θ_s ∈ [θ_c, θ_d] and the parameters a and b are set to define the maximum and minimum distance between the spiral segment and the cartesian axes.

By applying again the Propagation of Uncertainty theory with the same input conditions used for the generation of the results in figure 1(b), it is possible to compute the worst-case boundaries in the DoA estimation when the patch P2 has a spiral shape. As can be seen in figure 2(a), in this case the worst-case error is constant over the considered range, and considerably lower than the error for low values of θ that can be achieved for the conventional transducers mimicked by two point-like patches⁴ .

Capturing the DToA of the wavefront impinging on P1 and P2 is thus the mean to infer the DoA. However, accurately tracking the onset of the signal acquired by the two piezoelectric patches is very difficult in practical situations where multiple and dispersive guided modes can be present and acquisitions are affected by measurement noise. In such contexts, more robust DToA estimation procedures are necessary, such as those based on the tracking of the peak of the cross-correlation envelope. Unfortunately, even such procedures fail to produce reliable estimations when the interactions between the piezo patches (having finite dimensions) and the modes wavelengths are not taken into account. Such interactions are usually referred to as wavelength tuning effects [24] and are strongly influenced by the transducer shape.

In the following section, we will introduce a design strategy based on the Radon transform which allows to define the geometry of the piezo patches such that the effect of the wavelength tuning can be controlled and exploited for the DoA estimation.

The ground basis of the piezo patches' shaping strategy described in this paper is the formulation of the equations that govern the sensing of Lamb waves presented in [25]. According to such formulation, the voltage V_P(ω) generated by a surface mounted piezo patch of arbitrary shape Ω_P(x, y) in presence of a plane wave propagating at angle θ can be expressed as:

$$\begin{eqnarray}&&{V}_{P}(\omega )={jU}(\omega ){k}_{0}(\omega ){H}_{p}(\theta ){D}_{p}(\omega ,\theta )\end{eqnarray} \tag{ 4 }$$

where U(ω) denotes the amplitude and the polarization of the wave component relevant to the piezo properties of the patch at the considered frequency ω, k₀(ω) is the wave vector which characterizes the propagation, H_p(θ) is a quantity related to the material properties of the piezo-structure system, and finally D_p(ω, θ) is the sensor Directivity function which can be computed by the following integral:

$$\begin{eqnarray}&&{D}_{p}(\omega ,\theta )={\int }_{{{\rm{\Omega }}}_{P}}{e}^{{{jk}}_{0}(\omega )(x\cos \theta +y\sin \theta )}{\phi }_{p}(x,y){dxdy}\end{eqnarray} \tag{ 5 }$$

where ϕ_p(x, y) is referred to as shape function. Without lack of generality, here we consider the case of a single piezoelectric patch, with a single polarization, constant piezoelectric properties and ϕ_p(x, y) as a step function equal to 1 if the point of coordinates (x,y) belongs to the area of the piezoelectric patch Ω_p(x, y), and 0 elsewhere.

From equation (5) follows that D_p(ω, θ) can be computed from the coefficients in the direction θ of a bidimensional (spatial) Fourier transform (FT) of ϕ_p(x, y). It is worth noting that the Projection-slice theorem [26] states that the bidimensional FT of the initial function along a line at inclination angle θ is equal to the mono-dimensional FT of the Radon transform (acquired at angle θ) of that function. The Radon transform at angle θ of a generic shape function ϕ(x, y) defined in domain Ω can be computed by the following formula:

$$\begin{eqnarray}&&{{ \mathcal R }}_{\theta }(\varrho )[\phi (x,y)]={\int }_{{\rm{\Omega }}}\phi (x,y)\delta (\varrho -x\cos \theta -y\sin \theta ){dxdy}\end{eqnarray} \tag{ 6 }$$

and it consists of multiple line-integrals, as schematically depicted in figure 3.

Zoom In Zoom Out Reset image size

From equation (5) and the Projection-slice theorem we can conclude that if two piezoelectric patches ϕ₁ and ϕ₂ have the same Radon transform along a given direction θ apart from a spatial shift ϱ₀, i.e. ${{ \mathcal R }}_{\theta }(\varrho )[{\phi }_{1}]={{ \mathcal R }}_{\theta }(\varrho +{\varrho }_{0})[{\phi }_{2}]$, their directivity functions D₁ and D₂ (when the patches are solicited by the same plane wave) differs only by a phase shift which is directly related to the spatial shift, i.e. ${D}_{1}(\omega ,\theta )={D}_{2}(\omega ,\theta ){e}^{-{{jk}}_{0}(\omega ){\varrho }_{0}}$.

Similarly, assuming an undamped propagation the amplitude of the frequency response of the two patches is equal, a part from a scaling factor due to the possible variations associated to the term H(θ) which does not depends on the frequency.

The aim of this section is to provide a method to synthesize the shape of a sensor such that its directivity along a given direction is equal to a prescribed one apart from a phase-shift linearly dependent with the direction itself.

By doing so, two goals can be achieved:

• the uncertainty in the DoA estimation is minimized, because having the phase-shift linearly dependent with the DoA is equivalent to have a linear dependence between DToA and DoA proposed in the previous section;
• the detrimental effect of wavenumber tuning is avoided with the directivity function design, and, consequently, the similarity between the frequency responses of two patches can be used to estimate the difference in time of arrival DToA with procedures based on dispersion compensation and cross-correlation as those presented in [27].

The shape synthesis procedure exploits the fact that the Radon transform can be inverted with suitable algorithms [28], and consists of the following steps:

• (i)  
  define the geometry of the piezo patch P1, by designing its shape function ϕ₁(x, y) (a regular geometry such as that of a disc can be selected);
• (ii)  
  compute the Radon transform ${ \mathcal R }{1}_{\theta }(\varrho )$ of the shape ϕ₁(x, y) associated to the first piezo patch;
• (iii)  
  design the desired Radon transform of the shape ϕ₂(x, y) associated to the second piezo patch so that:

  $$\begin{eqnarray}&&{ \mathcal R }{2}_{\theta }(\varrho )={ \mathcal R }{1}_{\theta }(\varrho -{\varrho }_{0}-\alpha \theta )\end{eqnarray} \tag{ 7 }$$

  in a predetermined interval [θ₁, θ₂], it is worth noting that selection of the parameter α directly influences⁵ the sensitivity of the sensor to the variation of the DoA θ;
• (iv)  
  compute the inverse Radon transform (IRT) of ${ \mathcal R }{2}_{\theta }(\varrho )$ to obtain ${\hat{\phi }}_{2}(x,y);$
• (v)  
  apply a binary quantization procedure⁶ which transforms ${\hat{\phi }}_{2}(x,y)$ in a step function ϕ₂(x, y).

It is worth noting that the sensor design strategy just presented is not based on the knowledge of the waveguide geometrical and mechanical properties (dispersion curves), i.e. it is material agnostic and can be applied both to isotropic and anisotropic waveguides. An example of the results which can be achieved with the shape synthesis procedure is represented in figure 4(c). In particular, P1 was assumed to be a piezoelectric circular disc (diameter equal to 10 mm), then the shape of the second patch P2 was designed following steps 2 to 5, and using the following parameters: ρ₀ = 3.5 mm, α = 0.56 mm/degrees, θ₁ = 0° and θ₂ = 90°. Such parameters generate a transducer whose geometrical dimensions are similar those of the Archimedean spiral depicted in figure 2(a), as it can be seen from figure 4(c) where the Archimedean spiral of figure 2(a) has been over imposed as dashed line. The resulting ${ \mathcal R }{2}_{\theta }(\varrho )$ is depicted in figure 4(a). The IRT of ${ \mathcal R }{2}_{\theta }(\varrho )$ (i.e. ${\hat{\phi }}_{2}(x,y)$) is then depicted in figure 4(b). Finally, the piezopatch shape shown in figure 4(c) is obtained by selecting ϕ₂(x, y) as the region of the domain in which ${\hat{\phi }}_{2}(x,y)$ is greater than the 10% of its maximum value. By comparing subplots (a) and (d), it can be observed that the quantization procedure somehow degrades the desired RT. However, in the next section, we will show that such degradation has a negligible impact in the DoA performances.

Zoom In Zoom Out Reset image size

In order to assess numerically the performances of the proposed methodology to shape transducers, we have simulated a Lamb waves setup which consists of an aluminum plate 3-mm thick (Young's modulus 70 GPa, Poisson's coefficient 0.3 and material density 2700 kg m⁻³). The response of a shaped piezo-patch to an impact was simulated using the Greens function formalism adopted in [29]. In particular, the plate out-of-plane displacement w_i at the i − th location due to a point source characterized by a frequency spectrum g₀(ω) is estimated as:

$$\begin{eqnarray}&&{w}_{i}(x,y,\omega )={g}_{0}(\omega )\cdot {G}_{i}(x,y,\omega )\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{G}_{i}(x,y,\omega )=j{\pi }^{2}{H}_{0}^{(1)}(k(\omega ){d}_{i}(x,y))\end{eqnarray} \tag{ 9 }$$

where G_i(x, y, ω) is the Green function, defining the response to a unit source, ${H}_{0}^{(1)}$ is the Hankel function of the first kind and zero order and d_i(x, y) is the distance of the generic point of the plate at coordinates (x, y) to the impact location. Multimodal dispersive propagation is included by exploiting superposition and taking into account the dispersive property of each Lamb mode by means of the proper wavenumber-frequency k(ω) relation.

The overall signal recorded by a given patch thus becomes:

$$\begin{eqnarray}&&{s}_{t}(\omega )=\displaystyle \sum _{i=1}^{{N}_{r}}{w}_{i}(x,y,\omega )\end{eqnarray} \tag{ 10 }$$

where the w_i(x, y, ω) is the signal at each of the N_r discrete points belonging to the patch geometry.

Signals s_t(ω) were computed for the circular piezo patch P1 and the spiral shaped patch P2 as described in the previous subsection for different possible impact locations obtained by varying the distance (d_i = 100, 200, 400, 800 mm) and the DoA (θ = 0, 5,...,90°). The patches P1 and P2 and the impact locations, schematically illustrated as circles, are represented in figure 5(a). For each impact location the signal s_t(ω) was computed considering as g₀(ω) the spectrum of the pulse in figure 5(b) and a number of integration points N_r = 52 and N_r = 356 points, for patches P1 and P2, respectively.

Zoom In Zoom Out Reset image size

Then, the DToA Δt (θ) between the two piezo patches was computed by applying the cross-correlation procedure detailed in [27] which includes a dispersion compensation step based on the warped frequency transform [30]. Such procedure is compatible with the implementation of low-cost and low power devices [31]. Finally, the DoA was estimated by inverting equation (3).

The estimated DoA (θ_est) is represented in figure 6 as a function of its actual value (i.e. θ). In particular, figures 6(a) and (b) report the results related to the A₀ and S₀ mode, respectively. As can be seen the DoA estimation is very accurate for both modes. For all the considered impact points, the DoA estimated using both the A₀ and S₀ modes is clearly between the upper and lower bounds derived with the Propagation of Uncertainty theory for the ideal patches (as those shown in figures 1(b) and 2(b) which do not suffer the wavelength tuning effect), and the standard deviation of $| {\theta }_{{est}}-\theta |$ is well below 2 degrees in both cases.

Zoom In Zoom Out Reset image size

It is worth noting that: i) the simulated sampling frequency is f_s = 1 MHz, and the uncertainty in ToA was supposed to be 2/f_s; ii) the upper and lower bound errors are different in the two cases owing to the fact that the average group velocity of the two modes is considerably different in the frequency region which corresponds to the spectrum of the actuated pulse⁷ .

We can remark here that the spirally shaped sensor has intrinsically the capability to estimate the DoA also for waves characterized by an incoming angle within the third quadrant, i.e. 180 ≤ θ ≤ 270 degrees. In that cases Δt(θ) is negative rather than positive. A test in which the impact points were distributed on the third quadrant was considered. In particular, figure 7(a) shows the P1 and P2 patches and the location of the impacts points distributed on a circle at d_i = 100 mm from the center of P2. Figure 7(b) shows the estimated DoA with respect to the actual DoA.

Zoom In Zoom Out Reset image size

From this latest result follows that the DoA can be estimated over the 360 degrees by adding a second spirally shaped patch P3 rotated by 90 degree, around the center of P1, with respect to P2.

Finally, by combining the information gained by two of this clusters (P1-P2-P3), it is possible to infer the impact location through triangulation techniques, as those used in [6] and [18].