Forward and Inverse Studies on Scattering of Rayleigh Wave at Surface Flaws

The Rayleigh wave has been frequently applied in geological seismic inspection and ultrasonic non-destructive testing, due to its low attenuation and dispersion. A thorough and effective utilization of Rayleigh wave requires better understanding of its scattering phenomenon. The paper analyzes the scattering of Rayleigh wave at the canyon-shaped flaws on the surface, both in forward and inverse aspects. Firstly, we suggest a modified boundary element method (BEM) incorporating the far-field displacement patterns into the traditional BEM equation set. Results show that the modified BEM is an efficient and accurate approach for calculating far-field reflection coefficients. Secondly, we propose an inverse reconstruction procedure for the flaw shape using reflection coefficients of Rayleigh wave. By theoretical deduction, it can be proved that the objective function of flaw depth d(x1) is approximately expressed as an inverse Fourier transform of reflection coefficients in wavenumber domain. Numerical examples are given by substituting the reflection coefficients obtained from the forward analysis into the inversion algorithm, and good agreements are shown between the reconstructed flaw images and the geometric characteristics of the actual flaws.


Introduction
The non-destructive testing and structural health monitoring are of great importance in civil, mechanical, and aerospace engineering industries. Corrosions, pittings, and cracks may occur during service years of structures, and cause threats to public life and possession safety. Ultrasonic inspection is one of the fundamental techniques for locating surface and hidden flaws. Compared with traditionally-applied bulk waves, the Rayleigh surface wave has low attenuation and high energy concentration, which makes it suitable for inspection of defects or abnormalities on or close to surface. On the other hand, Rayleigh waves are not dispersive, as compared with other guided waves such as Lamb, SH, and Love waves, thus easier for signal processing [1,2].
Due to the great potential of Rayleigh waves in NDT (Non-Destrctive Testing) applications, it is necessary to examine their scattering and diffraction behaviors, both in the forward and inverse aspects. Researchers have long been interested in forward calculation of scattering phenomena for Rayleigh waves. One of pioneering work can be attributed to Ayster and Auld's paper of calculating the scattered field of Rayleigh waves by surface crack [3], which is followed by Achenbach, who discussed the emission of surface waves from a crack by cyclic loading in detail [4]. Furthermore, Lee and Liu deduced the analytical solution of a scattered Rayleigh wave field by semi-circular canyons [5]. However, theoretical calculations are only applicable for special cases. For more general situations, numerical approaches are required. As a widely applied numerical method, finite element boundary integral equation over the flaw area, it can be proved that the objective function of flaw depth d(x 1 ) is approximately expressed as an inverse spatial Fourier transform of reflection coefficients in wavenumber domain. This inversion procedure is especially suitable for reconstruction of weak and shallow notches, and can also be used as an initial-step image for deeper and steeper canyons.
The paper will be organized as follows: In Section 2, we will discuss about the improved forward problem using modified BEM in Section 2.1, and give some numerical examples and compare with existing results in Section 2.2. In Section 3, we will elaborate on the mathematical formula for inverse reconstruction method for surface-breaking canyon-shaped flaws in Section 3.1, and parametric influences will be discussed in detail in Section 3.2. We will give conclusions in Section 4.

Formulation of Modified BEM
The problem to be solved is shown in Figure 1. An artificial canyon-shaped flaw has been cut on the surface of an isotropic elastic half-space. The incident Rayleigh wave propagates from left to the right side, and the reflected wave is observed at the far end of the left side.
For accurately computing the scattering wave field, we divide the forward analysis into two major steps, as shown in Figure 2.
(i) A time-harmonic incident Rayleigh wave field in the intact half-space is considered, and the correspondent stress σ inc ij is calculated along the curve Γ − 1 , which is actually the flaw surface. (ii) Since the actual traction on Γ − 1 is zero, the stress σ inc ij is released by exerting a compensation traction t sca The scattering problem is transformed into the question of solving the wave field in the flawed half-space subjected to load t sca i .  Then, the displacement in frequency domain of scattering wave field at an arbitrary point can be expressed in the form of boundary integral equation (BIE), as shown by Equation (1) [34].
where u sca i and t sca i are displacement and boundary traction, respectively. G ip (x, X, Ω) and T ip (x, X, Ω) are fundamental solutions representing the i-th component of the displacement and stress at point x due to a unit time-harmonic concentrated load of angular frequency Ω exerted in the p-th direction at point X in the elastic full-space. The coefficient c pi (X) = δ pi for X ∈ V, c pi (X) = δ pi /2 for point X ∈ Γ which is smooth locally, c pi (X) = 0 for X / ∈ V ∪ Γ. It should be noted that here the boundary Γ includes both the flaw area and the intact infinite surface, as the full-space fundamental solutions have been chosen. Now the whole surface is divided into six parts, Γ − ∞ and Γ + ∞ tending to infinity on the left and right side, respectively; Γ − 0 and Γ + 0 representing the near-field boundaries; Γ − 1 and Γ + 1 the left and right half of the flaw surface, respectively. Consequently, the scattered wave field is rewritten from Equation (1): where the traction-free boundary condition on are meshed, thus, the first item on the right hand side of Equation (2) is omitted. This treatment can cause spurious reflected waves at the truncation points which do not attenuate in the 2D half-space, which are unwanted.
It is noticed that, if the truncation point is far enough from the flaw area, the displacement fields on the infinite boundaries will follow a Rayleigh wave pattern travelling out of the flaw region [18], as expressed by where u ± i . is the displacement field of the unitary Rayleigh wave travelling in positive (+) or negative (−) directions, and R ± the coefficient to be solved. Equation (3) makes it possible to take the far-field integrals in Equation (2) into consideration. Then, for any point on the boundary which is locally smooth, we have where the correction terms A ± p having the definition . They are constants. In paper [18], the coefficients R ± are not explicit in the modified BEM formulation, consequently, the reflection coefficients obtained from x 1 and x 2 components of displacement might not coincide with each other. As an improvement, we will try to solve R ± together with displacement fields.
The correction terms A ± p involve the integral over an infinite boundary, however, by applying the reciprocity theorem, the integral can be transferred to a bounded region [18,35].
where Γ 2 represents an auxiliary vertical boundary of 2-3 Rayleigh wavelengths, see Figure 2. t + i = − σ + i1 n 1 , and t − i = σ − i1 n 1 are defined on Γ 2 with n 1 the unit normal vector pointing in positive x1-direction. σ ± i1 is the stress field of the unitary Rayleigh wave propagating in positive or negative x1-directions, respectively.
By substitution of Equations (5) and (6) into Equation (4), we can obtain an equation set with the displacement and coefficients R ± as the undetermined variables. Here, we adopt the constant elements (see Figure 3) in the BEM discretization, in which one element is represented by the only node in the middle, and obtain the matrix representation in the form of Here, we assume that for the 1st or the N-th node, the displacements are not independent, and can be expressed as , respectively, where x 1 and x N are coordinates of the 1st and the N-th node, respectively. Thus, the submatrices concerning the 1st and the N-th elements in [T] are combined with the correction terms A ± p , and placed as the first and last columns. Consequently, the vector {u sca } contains the displacements of nodes 2~(N − 1), together with coefficients R ± , as shown by where u Pi represents the displacement in ith direction at the Pth node. The matrix [T] is expressed in detail as The elements in the first and last columns are defined as (5) and (6), and the elements where N 0 is the number of nodes on which the tractions are nonzero. The elements in [G] can be written as G PiJk = Γ P T ik x, x J , Ω dΓ, and the elements in {t sca } is t sca Pi x, x J , Ω . Since the matrix [T] is not square-sized, the equation is solved in the least-square sense to achieve error minimization. In the final solution, we can obtain reflection coefficient R − at the first and transmission coefficient R + . at the last element simultaneously, as well as the whole near-field displacements.

Verification of Reflection Coefficients
For verification of the reflection coefficients obtained from Equation (7), we adopt the scattering equation in half-space, which uses the half-space Green's function to calculate the far-field wave pattern, as shown in [36]: where C ijkl is the stiffness tensor, n j is the unit normal vector pointing outward the flaw area, and G H ln (x, X) represents the Green's functions in the half-space. In the far-field, the half-space Green's function can be expressed in the Rayleigh wave form [37]: where in which k L = Ω/c L , k T = Ω/c T , ξ 0 are wavenumbers of S, P, and Rayleigh waves, respectively, and The reflection coefficients derived from Equation (11) should coincide with that obtained directly from Equation (7). A numerical example will be given to show the validity of reflection coefficients. Moreover, the scattering Equation (11) is more useful in the inverse reconstruction procedure of the flaw shape.

Results of Forward Problem
Example 1. To verify the modified BEM method, we now try to solve the classical Rayleigh wave scattering problem presented by Kawase [10], for which the solutions have already been obtained from discrete wavenumber method applying half-space Green functions.
As shown by Figure 4a, an artificial half-circular canyon-shaped flaw is set on the surface of an isotropic elastic half-space, with Poisson's ratio of 1/3. The incident wave is a unitary Rayleigh wave travelling from left to right, whose time-domain waveform is expressed as in which Ω = Ωr/c T is the normalized circular frequency, where r is the radius of canyon, and c T is the shear wave velocity. T = tc T /r represents the normalized time array. The time-domain signal is shown in Figure 5, when Ω = 2π and T 0 = 2π.  We perform a Fourier transform to the incident wave signal in time-domain to obtain its frequency spectrum, and calculate the displacement field and reflection coefficients at each frequency point by the modified BEM. As an example, the whole displacement field along the surface at frequency Ω = 2π is plotted in Figure 6, and is compared with the results obtained by Kawase [10]. It is seen that the BEM displacements show good agreement with the previous results. Figure 6. Comparisons of the surface displacements obtained by the modified BEM and those in Kawase [10]. Displacement fields at other frequency points have also been calculated in the same manner. An inverse Fourier transform is then applied to the frequency-domain displacements to attain time-domain wave signals for different points along the surface, as shown in Figure 7a,b, from which we can see that the spurious reflection waves are not generated. The circular dots in Figure 8 shows the reflection coefficients obtained directly from the modified BEM program, and the rectangular dots represents those derived from half-space scattering Equation (11). The good agreement illustrates the validity of the forward analysis. This example will be used in the following inverse reconstruction.

Formulation of Inverse Reconstruction
We now try to reconstruct the position and geometric properties of the canyon-shaped flaw. In order to quantitatively depict the flaw's geometric characteristics, an objective function d(x 1 ) is defined in x 1 domain, which means the flaw depth at coordinate x 1 , whereas in the intact area, its value is zero.
The scattering Equation (11) in half-space is adopted for inversion. In this paper, a shallow flaw is considered, which can be thought of as a weak scatterer. In this case, the Born approximation [36] can be used for linearization, which replaces the total wave field by the incident wave field on flaw boundaries. e.g., u tot 1 ≈ u inc 1 . Then Equation (11) becomes It is noticed that, on the original half-space surface Γ ′ 1 (see the dotted line in Figure 1), the traction induced by half-space Green functions is zero. Thus, the surface for integration can be extended to Since surface Γ − 1 ∪ Γ + 1 ∪ Γ ′ 1 is closed, we can perform Gauss's theorem to the flawed region and obtain It is assumed that a pure incident Rayleigh wave mode is sent from left to right, and the reflected wave is observed at the far field on the left side. Thus, the incident wave with wavenumber ξ 0 can be expressed as where p i (x 2 ) (i = 1, 2) is defined in Equation (13), and A inc represents the complex amplitude of the incident wave. We use the far-field form of half-space Green's functions (12), and substitute them into Equation (17), and obtain the reflected wave of the same wave number ξ 0 , as shown by Equation (19): If one takes a thorough observation on Equation (19), he can find that p n (X 2 ) exp(−iξ 0 ·X 1 ) represents a unitary reflected Rayleigh wave travelling from right to left, and the integral over V is taken with respect to coordinate x of field points, and can be viewed as a complex amplitude of that unitary reflected Rayleigh wave. Also, the integrand is separated into the part in the bracket concerning x 2 coordinate, and an exponential item exp(+iξ 0 ·2x 1 ).
This inspires us to rewrite the expression of reflected wave into the following form where the Fun(x 2 , ξ 0 ) is the expression in the bracket in Equation (19), and the complex amplitude for the reflected wave is termed as Also, the volume integral is further rewritten as a dual integral: Fun(x 2 , ξ 0 )dx 2 (22) where d(x 1 ) is the flaw depth at coordinate x 1 . The integral range has been extended from minus infinity to infinity, since d(x 1 ) = 0 in the intact areas. After organization, the reflection coefficient is expressed as where c TT , c TL , and c LL are functions of wavenumber ξ 0 , and are expressed in detail as If we perform an integral over x 2 first, we can rewrite Equation (23) as Since the flaw depth is considered small, it is natural to apply Taylor expansion for the exponential functions in Equation (25) at d(x 1 ) = 0, i.e., Then, it is found that the objective function d(x 1 ) can be reconstructed using a spatial inverse Fourier transform, as illustrated by It should be noted that the inversion formula adopts a series of assumptions for linearization. Thus, parametric analysis should be performed to show the validity range of this method.

Numerical Examples
In this sub-section, numerical examples are given to show the validity of the inversion procedure. Also, parametric studies are performed to illustrate which kind of flaw is the most suitable for this reconstruction method.
Here, we mainly focus on the arc-shaped canyons as representatives of surface flaws. It is believed that other shapes can lead to similar conclusions.

Choice of Parameters
The radius r, maximum depth d max in Figure 4a, as well as the distance l in Figure Figure 4b for the double-flaw case, are the key geometric parameters for arc-shaped canyons; the maximum frequency Ω max is an important parameter of input data.
For generality, we perform dimensionless analysis in this paper. The normalized radius R = r/r is kept constant as 1; the depth is normalized as D = d/r; and frequency Ω = Ωr/c T , where c T is the shear wave velocity in the half-space. Meanwhile, the normalized coordinates are written as The readers can readily verify that, in dimensionless analysis, the model (i) with radius r, maximum depth d max , and frequency range 0 ∼ Ω max , is equivalent to model (ii) with radius 2r, maximum depth 2d max , and frequency range 0 ∼ 0.5Ω max . Thus, the absolute value of radius is of little importance. Instead, the depth-radius ratio d max /r (or D max /R) is a more essential parameter showing the "steepness" of the canyon. Another important factor is the normalized frequency range, which reflects the incident Rayleigh wavelength in comparison with the geometric flaw scale.

Reconstruction Results for Models of Different Depth-Radius Ratio
We established three canyon models with depth-radius ratios d max /r equaling 0.1, 0.2, and 0.5, respectively. For each model, we calculated the reflection coefficients of Rayleigh wave within the frequency range from 0 to Ω max = 5.0, with a uniform interval of dΩ = 0.2, as can be illustrated in Figure 8. These coefficients are substituted into the inversion Equation (27) to obtain the flaw depth function d(x 1 ).
It should be noted, Equation (27) takes an integral from minus to positive infinity. However, the reflection coefficients are only available at positive wavenumber points. Since d(x 1 ) is a real function, The functions c TT (ξ 0 ), c TL (ξ 0 ), and c LL (ξ 0 ) are also defined similarly. The whole integrand function is defined as 0 at ξ 0 = 0. Moreover, the integration range is reduced to 0 ∼ Ω max .
After discretization, the integral of Equation (27) becomes where ∆x 1 = π/N∆ξ 0 is the increment in x 1 coordinate. Integers p and k are indices in spatial and wavenumber series, respectively. The discrete inverse Fourier transform is easily performed by FFT (Fast Fourier Transform) algorithm.
The results are plotted in Figure 9. It is observed that when ratio d max /r is as small as 0.1 (see Figure 9a), the geometric properties (such as location, depth, and width) and shape of the flaw can be reconstructed with good precision. When d max /r value increases to 0.5 (see Figure 9c), the location of flaw can still be precisely detected, and the width and depth are slightly larger than the reality. However, the reconstructed shape is not as accurate as its shallower counterparts. The image underestimates the bottom region, and looks as if there were two local maximum depth points. The reason why the shallow and flat flaws are reconstructed better can be attributed to the application range of Born approximation, which is especially suitable for weak scatterers. Even the reconstructed image is not accurate for deeper and steeper canyons, the geometric characteristics (e.g., location, depth, and width) can be obtained with acceptable accuracy, so the image can be used as the "initial solution" in the following iterative reconstruction algorithm.

Reconstruction Results Using Different Frequency Range
It is necessary to discuss how the frequency range affects the inverse reconstruction results, which helps us understand the proper choice of incident wavelength in experiments. We have recalculated the inverse reconstruction for the flaws with d max /r equaling 0.2 and 0.5, using frequency range Ω = 0 ∼ 5.0, 0 ∼ 2.0, and 0 ∼ 1.0, respectively.
The reconstruction results are plotted in Figure 10. For the case of d max /r = 0.2 (see Figure 10a), the frequency range of Ω = 0 ∼ 1.0 (see the dash-dot curve) only gives a vague image, yet it still obtains the correct location. When the frequency range increases to Ω = 0 ∼ 2.0, the flaw shape seems to agree well with the reality, however, the image underestimates the flaw depth and overestimates the width. As the full data of reflection coefficients with Ω = 0 ∼ 5.0 are used in the inversion process, the flaw image becomes clearer, and the width and depth are more accurate, although the image for part of the bottom is deviated from reality. For the case of d max /r = 0.5 (see Figure 10b), if only the low frequency parts of reflection data (e.g., Ω = 0 ∼ 1.0 and Ω = 0 ∼ 2.0) are taken into the inversion process, the results imply the existence of the flaw but the images are very blurred. It is only when higher frequency portions are taken into consideration that the flaw image becomes clearer and width and depth are much more accurate.
From these two examples, it can be concluded that the low-frequency or long-wavelength components of the reflection signals (e.g., Ω = 0 ∼ 1.0) control the canyon's basic information, such as existence and location; while the usage of higher frequency components (e.g., Ω = 1.0 ∼ 5.0) increases the clarity of image and gives more detailed geometric information, such as depth and width. It should be noted that the effect of high frequency components of enhancing image resolution is limited by Born approximation. As the Rayleigh wave becomes shorter in wavelength, it will penetrate less deeply into the elastic half-space, thus, the canyon becomes a large scatterer, and Born approximation loses validity. Actually, in the previous cases, if the high frequency components for Ω > 5.0 are taken into consideration, the reconstructed images change little from the ones obtained from data with Ω = 1.0 ∼ 5.0.

Reconstruction Results of Double-Canyon Type Flaws
In order to examine the validity of our inversion approach in multiple flaw cases as seen in Figure 4b, we introduce here two models with double-canyon type flaws. For each model, two canyons with different geometric parameters are cut on the surface of the half-space with a distance l in between. The geometric parameters are set as (i) normalized radius: R 1 = R 2 = 1; (ii) normalized depth: d max,1 = 0.2, d max,2 = 0.1; (iii) normalized edge-to-edge distance between two flaws: 1 L = l/r = 2.0, 2 L = l/r = 0.2.
The inversion results for the two models are plotted in Figure 11, in which input data with frequency ranges Ω = 0 ∼ 1.0, Ω = 0 ∼ 2.0, and Ω = 0 ∼ 5.0 are shown in different curve types. It can be seen that, when the two canyons are distant (see Figure 11a), even the data with low frequency range can give correct locations of both two flaws and rough estimations of their depth and width. As the input frequency range increases to Ω = 0 ∼ 5.0, the geometric parameters of both two flaws are reconstructed with better accuracy.
On the other hand, if the two flaws are close with each other (see Figure 11b), the vague image reconstructed from low frequency range data cannot distinguish the two canyons. As the frequency range increases, the resolution becomes higher, and the images of two canyons are clearly separated. Thus, it is concluded that the frequency range Ω = 0 ∼ 5.0 is necessary for reconstruction of multiple flaws.

Conclusions
In this paper we present the studies of forward and inverse analysis on Rayleigh wave scattering by surface-breaking canyon-shaped flaws.
For the forward problem, we show that the modified boundary element method is an effective numerical technique to calculate the scattered wave field. By introduction of modification items concerning far-field displacement patterns in the BEM matrix, one can effectively suppress artificial reflected waves induced by model truncation. In our improved modified BEM version, the far-field reflection and transmission coefficients are solved simultaneously with the whole displacement field in the least-square sense. The improved version corrects Arias' original formulation, in which reflection coefficients obtained from x 1 and x 2 components might not coincidence with each other.
For the inverse problem, we present an inversion procedure to reconstruct the geometric properties of canyon-shaped flaw by usage of reflection coefficients of Rayleigh wave at a series of frequencies.
From scattering theory for elastic waves, the scattered Rayleigh wave field is expressed as a boundary integral equation over the flawed area. By introduction of far-field form of half-space Green's function and the Born approximation, it can be proven that the objective function of flaw depth is reconstructed by inverse Fourier transform of reflection coefficients. This method gives a direct and simple inversion approach for reconstructing locations and shapes of surface flaws, which does not require baseline data and iteration.
For numerical verifications, the reflection coefficients obtained from forward analysis are used as input data in the inversion procedure. It is shown that the inversion method is effective for both single and multiple flaw cases. The method can reconstruct images especially accurately for shallow and flat notches; while for deeper flaws, the method can still represent their geometric properties (e.g., the location, depth and width), and can produce reliable "initial images" for further usage.