Quantification of the Effect of Array Element Pitch on Imaging Performance

This paper investigates how the pitch of elements in periodic ultrasonic arrays is related to their imaging performance, with particular emphasis on imaging artifacts (grating lobes) arising from discrete spatial sampling. Although the classical Nyquist rules for array element pitch are well known, they only provide the limiting condition needed to eliminate grating lobes from an array with an infinitely large aperture at a single frequency. Physical arrays have finite-sized apertures and most applications employ broadband pulses. For these reasons, grating lobe artifacts are always present at some level, and practical array design is, therefore, based on suppressing grating lobe artifacts to a level appropriate to a given application. In this paper, a theoretical framework is developed that enables the point spread function of a periodic imaging array to be decomposed into the sum of contributions from a main lobe and different orders of grating lobes, thus allowing grating lobe artifacts to be unambiguously quantified. Numerical simulations are used to analyze the performance of 1-D linear arrays in both far-field (steering only) and near-field (focusing only) scenarios, and design guidelines are deduced. It is shown that in general, the classical Nyquist rules are overly conservative and that the pitch of an array can be increased without significantly compromising image quality, provided that certain constraints on ray angles are implemented in the imaging algorithm. Experimental examples are shown that illustrate the practical application to arrays in two configurations.


I. INTRODUCTION
U LTRASONIC arrays are widely used in medical [1]- [3] and nondestructive testing (NDE) [4]- [9] applications. They can be used simply as a tool for translating an ultrasonic aperture over a target or to perform dynamic steering and/or focusing of an ultrasonic beam. However, for periodic arrays, including the ubiquitous 1-D linear array, there is a lack of quantitative understanding of the effect of array element width and pitch. The concept that the element pitch should satisfy some rules to avoid the appearance of imaging artifacts (grating lobes) is well known but the interpretation of such rules is unclear and frequently violated without significantly deleterious effects. For example, the "half-wavelength rule" [10] for the upper limit on element pitch is the logical extension of Nyquist sampling theory to the spatial domain. However, this rule is often not adhered to even at the ultrasonic wavelength at the center frequency of an array, let alone at the upper limit of the array bandwidth. The objective of this paper is to quantitatively relate the array element pitch to the quality of image obtained and hence determine appropriate rules for specifying array element pitch in practical applications.
There are two drawbacks of using overly conservative array element pitch. First, there is the physical challenge of manufacturing array elements of the necessary size. Second, an increased number of elements is required to populate a given size of spatial aperture and this brings attendant problems of physical connectivity and data throughput. Both problems are especially acute for 2-D arrays.
In this paper, the imaging performance of periodic arrays is quantitatively analyzed. In Section II, a general model is developed that enables the point spread function (PSF) of an array to be written as a sum over contributions due to the main lobe and different orders of grating lobe. Expressions are derived for both far-field and near-field operations. In the special case of far-field, single-frequency operation, it is shown how the analysis relates to classical rules that predict whether grating lobe peaks exist and their position if they do exist. However, even in this case, effects such as the finite size of array aperture, the angular extent of the imaging region, and the nonuniform directivity of array elements mean that actual amplitude of grating lobe artifacts observed is considerably more complicated. In Section III, the models developed in the previous section are used to perform numerical simulations to quantify the imaging performance of arrays in both far-field and near-field operations. In the case of the latter, an imaging algorithm is considered where focusing is performed in transmission and reception over a specified aperture angle range, as conclusions from this case yield important insights that can be readily extended to other cases. The simulations allow quantitative design guidelines to be elucidated that relate element pitch to the amplitude of grating lobe artifacts and (in the case of near-field operation) imaging resolution. The application of the design guidelines to real experimental data in two practical scenarios is presented in Section IV.

II. THEORETICAL DEVELOPMENT
The purpose of the theoretical development presented here is to obtain an expression for the PSF of a periodic array in a form that enables artifacts in the PSF due to spatial undersampling to be explicitly identified.
In general, to obtain a PSF, it is necessary to model the process of elastic wave excitation, detection, and scattering This work is licensed under a Creative Commons Attribution 3.0 License. For more information, see ht. tp://creativecommons.org/licenses/by/3.0/ in an unbounded, homogenous isotropic medium containing a single point target. The necessary modeling could be achieved using any suitable direct numerical simulation technique (e.g., finite elements or finite difference methods). However, even packages such as k-wave [11], POGO [12], and PZ-Flex [13], [14] that are optimized for wave propagation simulations are extremely computationally expensive, and for the simple problem at hand, they are unnecessary. More suitable are programs such as field II [15] and FOCUS [16], [17] which are designed to predict ultrasonic fields and scattering from targets in otherwise homogeneous media. Field II uses the concept of one-way spatial impulse response developed by Tupholme [18] and Stepanishen [19] in conjunction with reciprocity and convolution to determine the overall response of a transducer to a point target when it is excited with a particular input signal. FOCUS uses either the fast nearfield method or the angular spectrum method to calculate the pressure field from single transducers and phased arrays. However, while these are essential tools for the design of arrays and imaging algorithms for specific applications, they do not allow the separation of grating lobe artifacts from other imaging artifacts and this motivates the method described in this section.
As far as possible, the method is developed using a general vector notation that is applicable to either a 2-D wave scattering model describing the operation of a periodic 1-D linear array, or a 3-D wave scattering model describing the operation of a periodic 2-D array. For clarity, the latter case will be described, with any specialization necessary for the 2-D case noted where appropriate.

A. General Framework
Consider a planar periodic array with a total of N identical elements that is operating into a homogeneous isotropic medium as shown in Fig. 1(a). The medium may be either solid or fluid; in both cases, only a single-wave mode in the medium is considered. The array elements are assumed to be reciprocal, in that they have the same angular sensitivity in both transmission and reception. The normal to the plane of the array is denoted by the unit vectorn. Unit vectorsp 1 and p 2 are the lattice vector directions, and p 1 and p 2 are the element pitches in the lattice vector directions.
In its most primitive form, the final electrical output from an array can be represented as a matrix of the time, t, domain responses from each possible transmit-receive element combination f i j (t) or its frequency, ω, domain equivalent F i j (ω). Here, the subscripts i = 1, . . . , N and j = 1, . . . , N are the indices of the transmitter and receiver elements in the array, respectively. In the NDE community, it is increasingly common to experimentally acquire f i j (t) directly, a process referred to as full matrix capture (FMC), and perform all imaging in postprocessing. In medical imaging, differential movement within the target over the timescale of acquisition precludes FMC in many cases. Whether FMC is performed or not, f i j (t) provides the fundamental data from which all other data can be derived (assuming linearity of the ultrasonic field in the medium). It will be assumed here that f i j (t) is analytic (i.e., it comprises both the real physically recorded data and an imaginary component in quadrature as obtained through the Hilbert transform). With this in mind, the output of most linear imaging algorithms can be written in the following time-domain form where |I (r)| is the image amplitude at position r, and τ i j (r) and a i j (r) are the time delay and amplitude weighting that describe the imaging algorithm. Note that because f i j (t) are analytic functions, the underlying image quantity I (r) is also complex valued, but its magnitude |I (r)| is a positive, real-valued envelope suitable for display and interpretation. Imaging algorithms that cannot be described in this form are those in which multiple values at different points in time from the same f i j (t) contribute to the same image point (as is the case when multiple plane-wave excitations are used [9]) and algorithms where one or both of τ i j (r) and a i j (r) are frequency dependent. Both of these additional possibilities for linear imaging can be accommodated if (1) is instead written in the frequency domain where the imaging algorithm is now encapsulated in the complex coefficient A i j (ω; r). For the subset of linear imaging algorithms that can be expressed in the time domain by (1), the associated frequency-domain coefficient Consideration is now confined to the most common type of imaging algorithms where the same transmit and receive focal laws are used. This implies that A i j (r;ω) → A i (r;ω)A j (r;ω). For single scattering that can be described by a scalar function, the PSF can be used to assess the performance of the array and imaging algorithm. The PSF, P(r, q), is the image, I (r), that results from an isolated, omni-directional point target (i.e., a target that scatters equally in all directions for all directions of incident wave) at position q. The scattering coefficient of the target is unity. Although the scattering of elastic waves (even in the single-scattering regime) cannot strictly be described by a scalar function, the PSF remains a useful tool for assessing imaging performance.
To predict the PSF, an expression is needed for the FMC data, F i j (q;ω), resulting from an isolated point target at q. The incident ultrasonic displacement field at q due to point excitation at a position u on the surface of the array is denoted by G(u − q;ω), which can be loosely described as Green's function of the system. This has the form where D 0 (x; ω) describes the directivity of a point source on the surface of the array (throughout this paper, a circumflex accent over a vector denotes a unit vector in the same direction, e.g.,x = x/|x|), and k = ω/c is wavenumber where c is the speed of sound in the medium. For a fluid, D 0 (x; ω) is independent of direction. For an elastic half space (the most relevant for direct-contact NDE applications), exact expressions for D 0 (x; ω) for both longitudinal and shear waves exist [20]; however, a reasonable approximation for the longitudinal wave radiation pattern in an elastic solid is a simple cosine dependence on angle, i.e., D 0 (x; ω) =x ×n. β describes the rate of reduction in amplitude due to beam spreading and is equal to 0.5 in 2-D and 1 in 3-D. The point scatterer at q may be regarded as a secondary source; the response to which at position v on the surface of the array is also given (by reciprocity) as G(v − q; ω). The incident ultrasonic displacement field at q due to excitation by the i th array element is obtained by integrating G(u − q; ω) over u ∈ i , where i is the area of the i th element. Similarly, the response of the j th receiving element is obtained by integrating G(v − q; ω) over v ∈ j , enabling the FMC data for a point target at q to be written as where F 0 (ω) is a function that encapsulates the combined effects of the frequency spectrum of the time-domain electrical signal sent to a transmitting array element, together with the frequency response of the transmitting and receiving elements. Substitution in (2) enables the PSF to be written as Let u i be the position of the center of the i th array element as shown in Fig. 1 and let the elements in the array have identical shapes so i = +u i . The substitutions u i = u −u i and q i = q − u i are made in the integral over the element (so du = du i and G(u − q; ω) = G(u i − q i ; ω) and the area of integration becomes ). In practice, the imaging region of practical interest for an array, while not in the far field of the whole array, is almost invariably in the far field of individual elements. Hence, |q i | |u i | and G(u i − q i ; ω) may be replaced by the far-field approximation [21] The integration over the area of an element yields the farfield element response where D(q i ; ω) and B(q i ; ω) describe, respectively, the overall directivity of an element and loss in amplitude due to beam spreading The result of the integral over the area of an element has analytical solutions in certain cases [22], [23]: for elements in a 1-D array of width a 1 , the result is a 1 sinc((1/2)ka 1 sinq ·p 1 ); for rectangular elements of dimension a 1 × a 2 in a 2-D array, the result is a 1 a 2 sinc((1/2)ka 1 sinq ·p 1 )sinc((1/2)ka 2 sinq ·p 2 ).
Consideration is now restricted to imaging algorithms in which the imaging function conjugates the phase associated with the assumed propagation delay to and from each imaging point, i.e., A i (r; ω) = W i (r) exp(ik|r i |), where r i = u i −r, and W i (r) is a real-valued function describing the weighting of the contribution from each element, sometimes referred to as apodization or aperture shading [24]. In the case of farfield operation, this type of imaging algorithm is simply beam steering. In the case of near-field operation, it corresponds to focusing at the image point in both transmission and reception, as is the case in, e.g., the total focusing method (TFM) [25], [26], inverse wavefield extrapolation method [27], and wavenumber method [28].
Substitution of the far-field element response (7) and the weighting A i (r; ω) = W i (r) exp(ik|r i |) into the general expression for PSF (5) yields the PSF for this type of imaging algorithm For the purposes of transforming the summation over elements to an integration, two new continuous functions are now defined: the sampling function E(u) and the apodization function W(u, r). The sampling function describes the position of element centers in an infinitely large periodic array where u 0 is the position of a reference element. The apodization function W(u, r) has the property that W(u i , r) = W i (r) within the aperture of the array and W(u, r) = 0 outside the aperture; hence, the nonzero values of the product E(u)W(u, r) describes the position of the element centers in a finite-sized array.
Therefore, the general form of the PSF for a periodic array can finally be written as where r = r − u and q = q − u.
To summarize the terms in (11): F 0 (ω) is the combined frequency response of input signal, transmitter array element, and receiver array element; E(u)W(u, r) describes the spatial sampling and apodization of the array; D(q ; ω) describes the physical directivity of array elements; B(q ; ω) describes the beam spread associated with wave propagation; and the exponential term describes the phase shifts associated with wave propagation and the imaging algorithm.
In Sections II-B-II-E, the above expressions are manipulated in a manner that allows explicit separation of grating lobe artifacts due to spatial under-sampling.

B. Far-Field Analysis
Consider the case when |r| = |q|, and both are much larger than the extent of the array aperture, W(u, r) = 0, as shown in Fig. 1(b).
If the coordinate origin is taken to be the center of the array aperture, then for all points within the aperture |r| |u|. In this case, the following far-field approximations (see [21]) with respect to the overall array aperture (rather than individual elements) can be made With these approximations, the following expression for the far-field PSF as a sum over lobes can be obtained using the method described in Appendix A: where P FF mn (r, q) In (14), s = q t − r t , where q t and r t are the components of, respectively, the target directionq, and image directionr that lie in the plane of the array [as shown in Fig. 1(b)], and the term contains a scale factor and the phase perturbation associated with the grating lobe order. Each lobe in the response is associated with an integral over the aperture of the array, (· · · )du. For the main lobe, 00 (u) = (4π 2 / p 1 p 2 ) is constant and the integral results in constructive interference when the image direction is close to the target direction, i.e., r ≈q and, hence exp(iks · u) ≈ 1 in (14). For grating lobes, mn (u) adds a phase perturbation across the aperture causing constructive interference to potentially occur in some other image direction, i.e.,r =q. In the far-field case, the integral over the aperture has the form, (· · · ) exp(iks · u)du, of a spatial Fourier transform from the u to the ks domain. Furthermore, the form of the function mn (u) means that in the transform domain, the grating lobes are shifted copies of the main lobe. This means that the result in (14) can also be written as where where W(ks, r) is the spatial Fourier transform of the array aperture function W(u, r) in the ks domain. The function W(ks, r) tends to 0 away from the origin at ks = 0; hence, for m = 0 and/or n = 0, W mn (ks, r) describes an offset version of the same function centered at some other point in the ks domain. In terms of the resulting PSF, W 00 (ks, r) (c) +1 order lobe. (d) Superposition of −1, 0, and +1 order lobes using (18). (e) Direct prediction of PSF using (9). (f) Direct prediction of PSF using the FMC array data set obtained from field II [15].
corresponds to the main lobe, while W mn (ks, r) for m = 0 and/or n = 0 corresponds to a grating lobe. This means that in the far-field case, it is straightforward to predict the position and peak amplitude of grating lobes.

C. Near-Field Analysis
A similar mathematical treatment can be used to express the near-field PSF as a sum over different lobes. This is described in detail in Appendix B and leads to the following: where P mn (r, q) and mn (u) again contains a scale factor and the phase perturbation associated with the grating lobe order. As in the far-field case, there is an integral (· · · )du over the aperture associated with each lobe in the response. Again, for the main lobe, 00 (u) = 4π 2 / p 1 p 2 is constant and the integral results in constructive interference occur when the image point and target are close, i.e., r ≈ q hence exp[ik(|r | − |q |)] ≈ 1.
All nonzero values of m or n result in a phase perturbation in the integral over the aperture causing imaging artifacts due to the spatial under-sampling, which are the near-field equivalent of grating lobes. In contrast to the far-field case, the integral over the aperture cannot be interpreted as a spatial Fourier transform. Hence, numerical integration is required to determine the location and intensity of grating lobe artifacts.

D. Demonstration
An example of the decomposition into different lobe orders is shown in Fig. 2. Note that in (16) and (19), the combined frequency response F 0 (ω) is arbitrary. Here, as an example, the demonstration of the effect of signal bandwidth on image grating lobe performance is performed on an input signal which is a five-cycle, Gaussian-windowed toneburst, with an ultrasonic wavelength at the center frequency of λ c . A 12-element linear array with pitch equal to λ c is used to image a target at a depth of 8λ c . Fig. 2(a)-(c) shows the contributions to the PSF from grating lobe orders −1, 0, and +1, respectively. These are computed using (19), and Fig. 2(d) shows their superposition to form the overall PSF according to (18). For comparison, Fig. 2(e) shows the exact PSF obtained directly from (9), which is in close agreement with the result in Fig. 2(d). In addition, an independent validation was performed to increase the confidence of model derivation. In the validation, an FMC array data set was first simulated using field II [15]. Field II is a 3-D model, hence to obtain results suitable for comparison with a 2-D model, an array with elements that were 15λ c long perpendicular to the imaging plane was modeled. The analytic form of the FMC data set output from field II was obtained by the Hilbert transform. This was used as f i j (t) in (1) to generate the corresponding PSF, with a i j (r) = 1 and

E. Classical Far-Field Rules
It is instructive to relate the previous analysis to the classical rules for array design. These are based on the contents of the integral in (16) at a single frequency. For brevity, the factors that are not a function of frequency are dropped to leave the single-frequency far-field lobe contributions aš The range of possible values of |ks| is bounded: because s = q t − r t , |q t | < 1 and |r t | < 1, the largest possible value of |ks| must be less than |ks| max = 2k. Grating lobe artifacts in the PSF occur if there are contributions in the region |ks| < 2k from anyP FF mn (r, q;ω) other thaň P FF 00 (r, q;ω). This can be visualized in the ks plane, examples of which are shown in Fig. 3.
The general condition that must be satisfied to avoid grating lobes is, therefore, for all pairs of m and n except (m, n) = (0, 0). For the cases (m, n) = (±1, 0) and (m, n) = (0, ±1), this condition leads, respectively, to p 1 < λ/2 and p 2 < λ/2; these are the classic half-wavelength rules [10], [29]. For a periodic 2-D array, these must be satisfied by both lattice vectors, which implicitly means that the case (m, n) = (±1, ±1) also satisfies the grating lobe condition. However, for periodic 2-D arrays, an extra condition is required for the case (m, n) = (±1, ∓1) and this is An example where the criteria p 1 < λ/2 and p 2 < λ/2 are satisfied but (23) is not satisfied is shown in Fig. 3(a). In the limiting case of p 1 = p 2 = λ/2, (23) leads top 1p2 = 1/2, which means an angle of 60°between the lattice vectors. This is so-called hexagonal sampling, the most efficient 2-D sampling scheme [10], [29] and is shown in Fig. 3 The direction of potential grating lobes can be readily determined from the ks plane. Consider a target directionq. The main lobe is in image directionr =q (i.e., at ks = 0 in Fig. 3). Assuming the only sampling criteria not satisfied is p 1 < λ/2, the nearest grating lobe is at ks = 2πp 1 / p 1 , which means that the in-plane component of the grating lobe direction r (g) t is related to the in-plane component of the target direction q t by From this, the grating lobe direction itselfr (g) can be obtained. It is instructive to consider the simple case of a spatially under-sampled 1-D array with pitch p 1 > λ/2, lattice vectorp 1 = [1, 0, 0] T , and a target at an angle θ ≥ 0°relative to the array normal. In this case, q t = [sin θ, 0, 0] T and r (g) t = [sin θ g , 0, 0] T , where θ g is the grating lobe direction, which from (24) is then given by θ g = sin −1 (sin θ − λ/ p 1 ).
A simple relaxation to the half-wavelength pitch rule can be achieved by limiting the maximum angle between the image directionr and the array normaln. If the maximum image angle is θ max relative to the array normal direction, then |r t | ≤ sin θ max (note that the target directionq cannot be controlled in a similar manner). Therefore, |ks| max = k(1 + sin θ max ) rather than 2k and the requirement to avoid grating lobes given by (22) and (23) become instead and Here, this is termed the modified Nyquist rule.

A. Quantification of Far-Field Performance
The classical rules for far-field operation summarized in Section II-E are binary conditions: they are either satisfied or they are not. In practice, grating lobe artifacts are always present [since the Fourier transform of the aperture of a finite-sized array W(ks) extends at some level throughout the ks domain]. Hence, practical rules require a subjective decision on the amplitude of grating lobe artifacts that are acceptable and, in general, numerical analysis of the governing equations.
Here, only 1-D arrays are considered, so P mn → P m and p 1 → p. The maximum grating lobe artifact to main lobe amplitude is defined as a metric .
It should be stressed that M FF G (r) is based on specifying the image position r and then searching the space of all possible target locations q to find the one which gives the largest possible grating lobe amplitude at r. This is subtly different to specifying a target location and then searching for the largest grating lobe artifact in the image, which is a more common way of visualizing grating lobes via a PSF. The reason for the choice used here is that the imaging region (i.e., the choice of r) can be controlled whereas the location of possible targets cannot be.

B. Quantification of Near-Field Imaging Performance
Equation (18) enables the PSF of an array to be computed and the contributions from spatial under-sampling identified separately. However, to draw conclusions regarding the optimal pitch for an array, it is necessary to first define metrics that enable the imaging performance to be assessed. Again, only 1-D arrays are considered, so P mn → P m and p 1 → p. Two separate metrics are used to independently assess resolution and image artifacts due to grating lobes. The resolution metric M R (r) uses a measure of the area of the main lobe of the PSF, obtained by dividing the main lobe volume by its peak value Note that M R (q) is a function of main lobe only and is independent of the element pitch; however, M R (q) is a function of element width (due to the element directivity function) and in most real arrays, the element width is intentionally set equal to (or almost equal to) the element pitch.
As in the far-field case, the image artifact metric M G (r) again uses the ratio of the maximum amplitude of the contribution from first-order grating lobes to the amplitude of the main lobe peak

C. Reduction of Parameter Space for Near-Field Imaging
The above metrics allow the quality of a PSF to be represented by two numbers, but the parameter space for near-field imaging remains large: the PSF is dependent on the target position, the number, size, and pitch of elements in the array, the input signal spectrum F 0 (ω), and any apodization applied in the imaging algorithm itself. To draw general conclusions, it is necessary to significantly reduce the parameter space. Here, the following steps are taken to reduce the number of parameters.
1) All distances are normalized to the wavelength λ c at the center frequency of the array. 2) The array element width is assumed equal to the pitch p.
3) The array can have an unlimited number of elements. 4) The input signal spectrum is described by a Gaussian function with a fixed bandwidth. 5) The imaging algorithm is the TFM [25] but with an apodization rule that limits the contributions from array elements to those where the angle to the image point relative to the array normal is less than a specified value α referred to henceforth as the aperture angle limit: Such a rule was introduced previously to suppress backscatter when imaging planar composites [30]. This rule means that the total angle subtended by that active aperture at the imaging point is constant and equal to 2α. This is equivalent to a constant f -number of 1/(2 tan α).
It is worth noting that although the array can have an unlimited number of elements, the number of elements used to form the image at any point is finite due to the aperture angle limit. The steps above reduce the number of independent parameters to three: the array element pitch to wavelength ratio p/λ c , the aperture half angle α, and the image point depth to wavelength ratio z/λ c .

D. Far-Field Results
First, some example far-field, single-frequency results of a 1-D linear array with aperture size 20λ c obtained using (20) are presented in Fig. 4. The "perfect" far-field PSF is shown in Fig. 4(a): the array sees the target in the correct direction with no artifacts over the complete range of angles from −90°to 90°. The effect of modest under-sampling is shown in Fig. 4(b), in which the effect of array element directivity is ignored (i.e., elements are assumed to be omni directional). In this case, grating lobes become visible at −90°when the image angle exceeds a certain value. This value can be found to be 19.5°by rearranging the classical far-field expression (24) and setting θ g to −90°. Image angle limits of ±19.5°a re shown by the dotted lines in Fig. 4(b)-(d). However, it can be seen in Fig. 4(b) that imaging operation within these limits only ensures that the peak of potential grating lobes is excluded. The finite size of the array aperture means that the grating lobes have finite width. In this case, the limit on maximum imaging angle to maintain grating lobe artifacts redicted by (24) and the dashed lines indicate the actual angle limits required to keep the grating lobe level to below 40 dB relative to the main lobe.
below, for example, −40 dB, is the considerably smaller range of ±11.5°and is indicated by the dashed lines Fig. 4(b). Fig. 4(c) shows the results for the same array but this time including the effect of element directivity resulting from the elements having finite width (equal to the element pitch) and the inherent directivity associated with out-of-plane excitation of a solid half space. This results in a natural suppression of signals from targets at large angles, which, in contrast to the effect of finite-size aperture, acts to increase the imaging angle range. In this case, for a main lobe to grating lobe ratio of 40 dB, the maximum imaging angle range is ±21.5°i ndicated by dashed lines.
Finally, the effect of broadband pulsed operation can be included. Throughout this paper, a standard broadband pulse is used that is a Gaussian-windowed five-cycle toneburst, with the number of cycles defined by the −40-dB points of the window. In this case, it is necessary to form the PSF as a function of both angle and radial distance using (16). This is because the maximum grating lobe artifacts do not necessarily occur at the same radial distance as the target. In order to present the resulting PSF in a consistent format with the singlefrequency cases, the radial dimension of the PSF is collapsed and replaced with the peak amplitude along that dimension. The result is shown in Fig. 4(d). Because the grating lobe position is a function of frequency, the effect of using a broadband pulse is to blur the grating lobe artifacts out and reduce their peak amplitude relative to the main lobe. This enables the maximum imaging angle to be further extended to ±24.5°, indicated by dashed lines in Fig. 4(d).  5 shows a graph of maximum allowable steering angle versus element pitch in wavelengths for various cases. The classical half-wavelength rule (22) is indicated by the vertical gray line and the modified Nyquist rule (25) by the gray dash-dotted line. The solid black line shows the numerically calculated far-field steering angle limit required to maintain the grating lobe level 40 dB below the main lobe as defined by (27). This calculation is for an array with a finite-size aperture of length 20λ and takes account of the element directivity. Two competing effects cause this curve to depart from the curve of the modified Nyquist rule (25), in which it intersects at ∼0.7λ and a steering angle of ∼25°. For element pitches greater than ∼0.7λ, the physical directivity of the elements provides a natural suppression of grating lobe artifacts and hence enables a greater steering angle to be obtained than that predicted by (25). However, for pitches below ∼0.7λ, the finite width of grating lobes cause artifacts that exceed the 40-dB criterion to occur before the center of the grating lobe comes into view. This results in a reduction in the maximum steering angle compared to that predicted by (25). It is interesting to note that this curve actually passes the classical Nyquist criterion at ∼68°; in other words, a pitch of less than half a wavelength is required to steer beyond this angle and suppress grating lobe artifacts below 40 dB. Finally, the black dotted line shows the equivalent, numerically calculated result when the array is under broadband pulsed operation (a five-cycle, Gaussian-windowed toneburst) rather than single-frequency operation. Again, this curve intersects curve for the modified Nyquist rule (25) at a pitch of ∼0.7λ c and a steering angle of ∼25°, with a similar but larger departure to the single-frequency result on either side of this point. For pitches above ∼0.7λ c , the finite bandwidth causes a blurring and consequent suppression of the peak grating lobes, which is increased compared to the single-frequency case, as evident by comparison of Fig. 4(c) and (d). At pitches below ∼0.7λ c , performance becomes limited by the finite width of grating lobes.

E. Near-Field Results
The following near-field results are all obtained by simulations, based on (18) and the numerical integration of (19). As described above, the number of independent parameters governing near-field performance are reduced to three: p/λ c , α, and z/λ c . Again, the input signal is assumed to be a Gaussian-windowed, five-cycle toneburst. The range of p/λ c considered is 0.4 to 3 in 0.1 increments and the range of α considered is 5°to 45°in 1°increments.
There is an important, counter-intuitive consequence of the aperture angle limit that warrants some explanation before presenting the results. For a given aperture angle limit α, the aperture size D is proportional to the depth of the image point according to D = 2z tan α, while the near-field depth z NF of an aperture is proportional [21] to D 2 and given by z NF = D 2 /(4λ c ) = (z 2 tan 2 α)/λ c . However, near-field operation requires z < z NF . Hence, the image depth at which near-field operation is possible is obtained by setting Fig. 5. Graph of maximum far-field steering angle and near-field aperture angle limit as functions of element pitch to wavelength ratio for various rules. The finite-sized aperture array considered has an aperture size of 20λ (20λ c in the case of broadband pulsed operation), and the pulse considered is a five-cycle, Gaussian-windowed toneburst.
Hence, when an aperture angle limit is applied, near-field operation is only possible at depths z > z min (α). This in contrast to imaging with a fixed aperture size D, where nearfield operation requires z < z max (D) where z max (D) = D 2 /(4λ c ).
Since this paper is concerned with near-field operation, the main consideration must be image depths where z > z min (α). In the first instance, imaging performance is considered at a depth that is a prescribed multiple of z min (α), which reduces the number of independent parameters to two: p/λ c and α. The smallest z min (α) occurs at the largest aperture angle α = 45°, where z min (α)/λ c = 1. However, the PSF of a target at this depth would not be fully separated from the plane of the array (as the depth extent of the PSF is approximately equal to the number of cycles in the time-domain signal 5 multiplied by half the wavelength). To allow the PSF to be separated from the array at α = 45°, the target depth is set to be z = 4z min (α). The resulting behavior of M R and M G as functions of p/λ c and α is shown in Fig. 6(a) and (b). The result in Fig. 6(b) shows an increased level of image artifacts for increased p/λ c or increased α. From the point of view of determining a practically useful design rule, a relationship between the aperture angle α and the pitch p max needed to achieve a given value of M G is desirable. Plane waves of wavelength λ incident on an array at angle α have a spatial period on the plane of the array equal to λ/ sin α. For this reason, the proposed form of an empirical expression relating pitch to aperture angle limit is where c 1 is a constant that depends on the desired level of grating suppression. For M G = −40 dB, it has been empirically determined that c 1 = 0.5. The corresponding line is superposed on the results in Fig. 6(a) and (b) and can be seen to provide a good fit to the −40-dB threshold in Fig. 6(b). This line is also shown as the black dashed line in Fig. 5 to compare the predication from the other models. It is shown that the aperture angle limit in near field can be relaxed based on (32). Having established relationships for M R and M G at specific depths related to the near-field length, it is now appropriate to see how these results are affected by the actual depth of the image point. To this end, α and z/λ c are now used as independent parameters, the latter over the range 5-50 in increments of 1, and the pitch p/λ c is set to p max (α)/λ c according to (32). The results are shown in Fig. 6(c) and (d), which also show the minimum image depth for near-field operation z min (α) and the actual image depth 4z min (α) used to generate the results in Fig. 6(a) and (b). Considering Fig. 6(c) first, it is clear that the resolution metric M R is almost completely independent of depth for z > z min . Conversely, the artifact metric M G shown in Fig. 6(d) is a function of depth; however, it decreases monotonically with depth for z > z min . At z = 4z min , M G has dropped to −40 dB, which is expected as these were the image depth and value of M G used to compute the maximum pitch p max (α). Fig. 6(e) shows the resolution M R plotted as a function of aperture angle for all the points in Fig. 6(a) and (c) that satisfy p ≤ p max (α) and z ≥ z min . For these points, M R has negligible dependence on p/λ c or z and is essentially a function of α only. With this in mind, additional near-field results have been generated for the extended aperture angle range, 45°< α ≤ 60°, using p = p max (α) and z = 12z min , and these points are also included in Fig. 6(e). From the point of view of developing a design rule, it is desirable to obtain a simple empirical relationship between M R and α. Noting that the resolution of an imaging system is inversely proportional to the numerical aperture sin α, the proposed form of such an expression is where r 1 = 0.3 has been determined by fitting to the data points in Fig. 6(e). The resulting curve is plotted as a solid line. This provides a good fit to the data points up to α ∼ = 55°( equivalent to f -numbers greater than 0.35). Beyond this point, the resolution actually worsens slightly. Up to this point, the lateral size of the PSF is controlled mainly by interference effects and axial size of the PSF is controlled mainly by the spatial length of the broadband pulse. However, at larger aperture angles, the lateral size of the PSF also increasingly influenced and ultimately compromised by the spatial length of the pulse. It is worth noting that such large aperture angles can rarely be achieved in real applications as the maximum depth that could be imaged with such a large aperture angle from a given length of array is low. Of much more practical importance is the resolution performance at angles up to ∼ = 25°( f -numbers greater than unity).

IV. EXPERIMENTAL EXAMPLES
The previous analysis has quantified the effect of grating lobes due to spatial under-sampling in both the far-and nearfields of an array. For the latter, it was chosen to impose an aperture angle limit on the imaging algorithm to reduce the parameter space. In this section, this technique is applied to experimental data obtained in two configurations.

A. Equipment and Sample
An aluminum sample of depth 50 mm is used for all tests with the geometry shown in Fig. 7(a). The longitudinal wave speed in the sample is c = 6300 m · s −1 . The salient feature of the sample is a row of 1-mm-diameter side-drilled throughhole (SDH) targets with a pitch of 10 mm at a depth of 20 mm. Perpendicular to the imaging plane, the thickness of the sample is 20 mm. A 5-MHz, 64-element array manufactured by Imasonic, Besançon, France, with a nominal center frequency f c = 5 MHz is used for all experiments. FMC is performed using a commercial array controller device (MicroPulse FMC manufactured by Peak NDT Ltd., Derby, U.K.) which is then processed off-line in MATLAB (Mathworks, Natick, MA, USA). The physical pitch of the array is 0.63 mm, which is equal to half of the wavelength, 0.5λ c , of longitudinal waves in aluminum at the nominal center frequency of the array. However, by collecting FMC data and summing the data from groups of contiguous elements, it is possible to synthesize FMC data from virtual arrays with 32 elements at 2 × 0.63 = 1.26 mm (1λ c ) pitch and 21 elements at 3 × 0.63 = 1.89 mm (1.5λ c ) pitch.

B. Normal-Incidence Direct-Contact Imaging
In the first example, the array is used in direct contact with the sample as shown in Fig. 7(b) with coupling provided by a thin layer of ultrasonic coupling gel. In this and subsequent diagrams and images, the dotted and dashed lines are drawn from the extremes of the array aperture inclined at the aperture angle limit. If the aperture angle limit is imposed in the imaging algorithm, the dashed lines indicate the extent of the area that can be imaged; for points outside this region, there are no ray paths to any element in the array at angles that satisfy the aperture angle limit. The triangular area between the dotted lines in the sample is the area where points are imaged with the maximum possible aperture permitted by the aperture angle limit. The resolution within this region should be constant and determined by (33). Between the dotted and dashed lines are regions where points are imaged with less than the maximum possible aperture so the resolution deteriorates.  Fig. 8(a)-(c) shows the results of processing the data obtained from the 0.5λ c , 1λ c , and 1.5λ c pitch arrays with the conventional TFM algorithm (i.e., using the full array aperture at each imaging point). The image in Fig. 8(a) shows the best resolution of the point targets including those outside the footprint of the array. It can be seen from the images in Fig. 8(b) and (c) that the performance of the conventional TFM algorithm deteriorates rapidly as the array element pitch is increased first to 1λ c and then to 1.5λ c , with substantial artifacts appearing. In Fig. 8(b), most of the artifacts are grating lobe effects associated with the reflection from the back wall of the sample, while in Fig. 8(c), artifacts from the SDHs themselves are also visible at depths of z < 20 mm directly below the array. Such images are unusable for most practical purposes.
The effect of applying aperture angle limits is now considered. Equation (32) is used to calculate the appropriate angle limits shown in Table I, and the resulting images are shown in Fig. 8(d)-(f). Note that in the case of the array with 0.5λ c pitch, the aperture angle limit is 90°, i.e., no aperture angle  limit is required. Since this would lead to an identical result to that shown in Fig. 8(a), an image based on an aperture angle limit of 30°is instead shown in Fig. 8(d). The imposition of this aperture angle limit reduces the observable region, but the central three targets in Fig. 8(d) are in or very close to the region where they are visible over the complete specified aperture angle range. As expected from (33), the resolution of these targets with a fully sampled array and an aperture angle limit of 30°is very similar to that obtained without a limit, i.e., Fig. 8(a). For the arrays with pitches of 1λ c and 1.5λ c , comparison of Fig. 8(e) and (f) with Fig. 8(b) and (c) shows that the grating lobe artifacts have been eliminated throughout the observable region. In Fig. 8(e), the resolution of the visible targets is indistinguishable from that in Fig. 8(d), despite the array having twice the pitch and half the number of elements. The result in Fig. 8(f) does show a noticeable loss of resolution as the required aperture angle to suppress grating lobe artifacts at this element pitch is now only 19.5°.

C. Normal-Incidence Immersion Imaging
A common inspection configuration in industrial inspections for NDE is to use a layer of liquid couplant between the array and the component under inspection (the sample) to allow the array to be scanned. An equivalent algorithm to the TFM for immersion imaging, the immersion TFM, can be defined that is focused in both transmission and reception at every image point. The necessary phase delays for such an algorithm are obtained using the technique described in [31] to efficiently calculate the ray path between each array element and each image point, taking account of refraction at the couplantsample interface. In most cases, the ultrasonic waves in the coupling liquid have a considerably lower ultrasonic velocity than those in the sample. According to Snell's law, this means that at the interface between couplant and sample, there is a significant amplification of the angle of a ray relative to the surface normal. The following propositions are made for extending the previously described techniques to immersion inspection.
1) The presence of grating lobe artifacts in the image depends on the angle of rays in the couplant relative the array normal and to the ultrasonic wavelength in the couplant according to (32). 2) The image resolution is determined by the angular range of rays reaching the image point in the sample and the ultrasonic wavelength in the sample according to (33). 3) Beyond the relevant critical angle, there are no physical ray paths between an array element and a point in the sample; this sets an implicit limit on the angles of rays that are considered in reconstruction, even without the imposition of an explicit aperture angle limit. These propositions could be applied to determine either the required pitch of an array to achieve a given resolution, or the appropriate aperture angle limit to apply to an array with a given pitch. Here, the latter is considered, and water is used as the couplant, in the normal-incidence configuration shown in Fig. 7(c). At the center frequency of the array, f c = 5 MHz, the wavelength in water is λ W = c W / f c = 0.3 mm, where c W is the speed of sound in water and is taken to be 1480 m · s −1 . The appropriate aperture angle limits are shown in Table II. The first critical angle in water at a water-aluminum interface For the particular case of immersion imaging at normal incidence, there is an interesting phenomenon that arises from the form of the empirical equations (32) and (33) which govern the aperture angle limit to avoid grating lobes in the couplant and the resolution obtained for a given aperture angle limit in the sample. At normal incidence, both of these are functions of the sine of the angle of rays relative to the surface normal; however, the ratio of the sines of these angles is in turn governed by Snell's law. Assuming a couplant with a lower ultrasonic, a smaller aperture angle limit must be imposed on a given array when it is used in immersion rather than in direct contact. However, Snell's law amplifies the angle in the sample and the aperture angle achieved in the sample remains exactly the same as would be achieved in direct contact.
The resulting images are shown in Fig. 9, using the same format as in Fig. 8 with the top row of images corresponding to direct application of the immersion TFM algorithm and the lower row corresponding to the immersion TFM algorithm with appropriate aperture angle limits. Although the array with 2.1λ W pitch is significantly spatially under-sampled in water, the image without aperture angle limits in Fig. 9(a) exhibits negligible grating lobe artifacts. This is because the maximum angle of rays relative to the array normal required for the normal-incidence immersion TFM algorithm cannot exceed the first critical angle in water of 13.6°associated with the water-aluminum interface. In other words, the immersion TFM algorithm provides an implicit aperture angle limit anyway which substantially suppresses grating lobe artifacts in the sample compared to what might be expected for the pitch of array. When an explicit aperture angle limit of 13.6°i s applied, the result is unchanged as shown in Fig. 9(d). For the arrays with pitches of 4.2λ W and 6.3λ W , the results are similar to those for direct contact. Without an explicit aperture angle limit, significant grating lobe artifacts are present and these are completely suppressed when the appropriate aperture angle limit is imposed. Furthermore, the resolution achieved is almost identical to that achieved by the same array in direct contact for the reasons explained above. For the 4.2λ W pitch array, the loss in resolution compared to a 2.1λ W is minimal but for a 6.3λ W pitch array, the loss in resolution is noticeable.

V. CONCLUSION
An analytical model has been developed that enables ultrasonic array images to be decomposed into a sum over different grating lobe orders. The model has enabled explicit separation of artifacts due to spatial under-sampling in both far-and nearfield cases under both single frequency and broadband (pulsed) operation. Numerical simulations have enabled the following conclusions to be drawn.

1) Both the classical (half-wavelength rule) and modified
Nyquist spatial sampling criteria are inadequate for the purposes of practical array design. They both provide binary conditions that can only be fully satisfied by an infinitely long array operating at a single frequency. For an array with a finite-sized aperture under broadband operation, grating lobe artifacts are always present at some level. 2) For arrays operating under far-field conditions where only beam steering is performed, both the classical and modified Nyquist criteria are confounded because: i) finite aperture size increases the angular extent of grating lobes; ii) broadband operation causes angular smearing of grating lobes; and iii) directivity of finitewidth array elements suppresses high-angle grating lobes. Under near-field focused operation, the situation is complicated further as focusing implies a steering angle that varies over the array aperture. 3) Metrics have been defined that quantify the maximum amplitude of grating lobe artifacts in the imaging region and (in the case of near-field focused operation) the imaging resolution. 4) Simulations have shown that in general, larger array element pitches can be employed than predicted by the Nyquist criteria without significantly compromising imaging performance.
Practical guidelines for array design are summarized below. These are based on broadband pulsed operation using a fivecycle Gaussian-windowed toneburst and assume a 40-dB main lobe to grating lobe artifact ratio is desired. In these guidelines, wavelength refers to the wavelength of ultrasonic waves at the center frequency of the toneburst. 1) Under far-field operation, the relationship between element pitch and maximum steering angle for an array with a 20-wavelength aperture is indicated by the black dotted line in Fig. 5. As a rule of thumb, an array with half-wavelength pitch permits steering up to about 45°a nd one with wavelength-pitch permits steering up to 15°.
2) Under near-field focused operation, an array with halfwavelength pitch will not produce grating lobe artifacts anywhere in the imaging region. For larger pitch arrays, an upper aperture angle limit α (equivalent to a lower limit on f -number) should be imposed in the imaging algorithm to ensure grating lobe artifacts are eliminated. However, this also has the effect of limiting the region that can be imaged. (It is analogous to the maximum permissible steering angle in far-field operation.) The relationship between aperture angle and element pitch is given empirically by (32) and is shown by the black dashed line in Fig. 5. 3) The aperture angle limit governs the resolution of the imaging algorithm according to the empirical relation (33), with a larger aperture angle providing better resolution. However, the rate of resolution improvement decreases with increasing aperture angle and above about 30°(equivalent to f -numbers below unity), the improvement is somewhat marginal. The element pitch associated with an aperture angle limit of 30°is, from (32), equal to one wavelength. To put the last number into context, a periodic 1-D array with a pitch of 1 rather than 0.5 wavelengths requires 50% fewer elements to populate an aperture of a certain size and the resulting FMC data set is 75% smaller. For a periodic 2-D array, the equivalent reduction in element count is 75% and the reduction in FMC data set size is 94%, i.e., more than one order of magnitude. A pitch of one wavelength is proposed as a sensible compromise in array design in that it enables high resolution images to be achieved with significantly lower element counts than a fully sampled array.
For an array that is immersion coupled to a sample, the imaging resolution is determined by the angular aperture achieved at the imaging point in the sample. Conversely, the presence of grating lobe artifacts is determined by the angle of rays in the couplant relative to the array normal and their wavelength. In most practical scenarios, the couplant has a lower ultrasonic velocity than the sample. For an array with a given pitch, the shorter ultrasonic wavelength in the couplant means that a smaller aperture angle limit must be imposed in immersion to avoid grating lobes than would be required in a direct-contact configuration. The appropriate angle limit can be calculated according to (32) using the ultrasonic wavelength in the couplant. However, the angle of rays passing from the couplant into the sample is amplified due to refraction at the interface according to Snell's law. This means that the aperture angle achieved at an imaging point in the sample is higher than the aperture angle limit imposed in the couplant. Hence, even if an array is significantly under-sampled relative to the ultrasonic wavelength in the couplant, high-resolution imaging in immersion is still possible. In fact, the ratio of array pitch to ultrasonic wavelength in the sample remains the governing parameter for image resolution in a normalincidence immersion configuration.

A. Far-Field Analysis
The general PSF (11) can be written as P(r, q) = F 0 (ω)P(r, q;ω)dω, whereP(r, q;ω) can be regarded as the PSF associated with continuous-wave operation at single frequency. With the far-field assumptions given in (12)P(r, q;ω) becomeš The target directionq can be written asq = q n +q t , where q n = (q·n)n is the component normal to the plane of the array and q t is in the plane of the array, as shown in Fig. 1(b). In the complex exponent in (A.1),q · u can be written as q · u = (q n + q t ) · u = q t · u. Similarly, the look directionr can be split intonormal and in-plane components, r n and r t , so thatr · u i = r t · u i . Defining s = q t − r t means thať P FF (r, q;ω) The integral with respect to u can be recognized as a spatial Fourier transform, F{·} = (·) exp(iks · u)du, from Strictly there is a leading exponential term, exp(−iks · u 0 ) in eq. (A.5) that results in a phase shift depend on where the reference array element is relative to the coordinate origin. If the origin is at the center of the array, it will correspond to either an array element position, or the midpoint between two adjacent element positions in one or both of the lattice vector directions. This means that the reference element position can always be written as u 0 = a 1 p 1 /2 + a 2 p 2 /2, where a 1 and a 2 are integers. Consequently, at the position of every delta function in the ks plane, exp(−iks · u 0 ) = 1, so this term is dropped.) Substitution of the expressions for E(ks) and W(ks − ξ , r) into the convolution integral enables the shifting property of the delta functions to be exploited (A.14) For comparison with the near-field PSF, it is also instructive to express the far-field result in the same form

B. Near-Field Analysis
As in the far-field case, the general PSF (11) is first written as P(r, q) = F 0 (ω)P(r,q;ω)dω, whereP(r,q;ω) is the PSF associated with continuous-wave operation at single frequency. Next, the array aperture is divided into subapertures of area j . The integral over the array aperture is then written as the sum of integrals I j (r, q,ω) over the subapertures, so thať in which W j (u, r) is the apodization function of the j th subaperture and is equal to the overall apodization function W(u), for u ∈ j and 0 elsewhere. The size of each aperture is sufficiently small that both r and q are in its far field. Let u j be the center of the j th aperture and define r j = r − u j and q j = q − u j . The integration variable over a particular subaperture is changed to u = u − u j and the far-field condition implies |u | |q−u j | and |u | |r−u j |. The usual far-field approximations then enable the integral over a subaperture to be written as As the area j of the subaperture becomes smaller W j (ks j ) → j W(u j ) and hence the subaperture integral becomes I j (r, q;ω) = D(q j ; ω)B(q j ; ω) j W(u j , r) In the limit as i is allowed to tend to 0, the summation over j once again becomes an integral over u to leavě P(r, q;ω) =