# Theoretical and numerical modeling of Rayleigh wave scattering
by an elastic inclusion

Shan Lia,b a School of Traffic and Transportation Engineering, Central South
University, Changsha, 410075, China
b Department of Mechanical Engineering, Imperial College London, London SW7
2AZ, United Kingdom  Ming Huangb a School of Traffic and Transportation
Engineering, Central South University, Changsha, 410075, China
b Department of Mechanical Engineering, Imperial College London, London SW7
2AZ, United Kingdom  Yongfeng Songa a School of Traffic and Transportation
Engineering, Central South University, Changsha, 410075, China
b Department of Mechanical Engineering, Imperial College London, London SW7
2AZ, United Kingdom  Bo Lanb∗ a School of Traffic and Transportation
Engineering, Central South University, Changsha, 410075, China
b Department of Mechanical Engineering, Imperial College London, London SW7
2AZ, United Kingdom  Xiongbing Lia111Corresponding authors. E-mail addresses:
<EMAIL_ADDRESS>(S. Li<EMAIL_ADDRESS>(M. Huang),
<EMAIL_ADDRESS>(Y. Song<EMAIL_ADDRESS>(B. Lan),
<EMAIL_ADDRESS>(X. Li). a School of Traffic and Transportation
Engineering, Central South University, Changsha, 410075, China
b Department of Mechanical Engineering, Imperial College London, London SW7
2AZ, United Kingdom

###### Abstract

This work presents theoretical and numerical models for the backscattering of
two-dimensional Rayleigh waves by an elastic inclusion, with the host material
being isotropic and the inclusion having arbitrary shape and crystallographic
symmetry. The theoretical model is developed based on the reciprocity theorem
using the far-field Green’s function and the Born approximation, assuming a
small acoustic impedance difference between the host and inclusion materials.
The numerical finite element (FE) model is established to deliver relatively
accurate simulation of the scattering problem and to evaluate the
approximations of the theoretical model. Quantitative agreement is observed
between the theoretical model and the FE results for arbitrarily-shaped
surface/subsurface inclusions with isotropic/anisotropic properties. The
agreement is excellent when the wavelength of the Rayleigh wave is larger
than, or comparable to, the size of the inclusion, but it deteriorates as the
wavelength gets smaller. Also, the agreement decreases with the anisotropy
index for inclusions of anisotropic symmetry. The results lay the foundation
for using Rayleigh waves for quantitative characterization of
surface/subsurface inclusions, while also demonstrating its limitations.

Keywords: Rayleigh wave; Scattering; Born approximation; Finite element.

##  1  Introduction

Rayleigh waves propagating on a solid surface can be scattered by a flaw, such
as a void, crack or inclusion, due to the abrupt change of acoustic
properties. The scattered waves usually carry information about the geometric
and elastic properties of the flaw. Therefore, a thorough understanding of
flaw-induced Rayleigh wave scattering is critical for the nondestructive
evaluation (NDE) of flaw location and the characterization of flaw size,
shape, and other physical characteristics.

Theoretical modeling has been a major method for understanding Rayleigh wave
scattering, and most theoretical studies are based on the Kirchhoff and Born
approximations. The Kirchhoff approximation is mainly used for scatterers of
void or crack type because it treats a scatterer as body sources 1, 2, 3, 4,
5, 6, 7, while the Born approximation uses the incident field to replace the
field inside an inclusion 8, 9 and is therefore suitable for elastic
inclusions. The Born approximation has been extensively used for bulk wave
scattering 10, 11 and demonstrated to be effective and accurate when the
wavelength is larger, or comparable to, the scatterer size. The Born
approximation has also seen applications in studying Rayleigh wave scattering,
and example works include those of Auld 12 and Snieder 13. Apart from the
Kirchhoff and Born approximations, there have been other theoretical efforts
on investigating Rayleigh wave scattering. These include the first-order
perturbation theory for weak subsurface inclusions14, the reciprocity theorem
for linear 15 and nonlinear 16 Rayleigh waves, and the second-order
approximation for multiple shallow cavities 17. However, most theoretical
studies so far have focused on a specific type of inclusion/flaw, while some
studies involve difficult-to-solve equations 14, 15.

The subject of Rayleigh wave scattering has also received considerable
numerical studies. Various numerical schemes have been used, and most common
ones are the boundary element method and the finite element (FE) method. The
boundary element method was mainly used for scatterers that are surface voids
(cavities) 18, 19, 17. By comparison, the FE method is capable of dealing with
all sorts of scatterers and related examples include its uses in simulating
wave propagation 20 and scattering 21 in complex polycrystalline media. For
this reason, the FE method has been used to analyze the interaction of
Rayleigh waves with surface cracks 22, 23, 24, 16, and most recently, it was
successfully applied to predicting the attenuation and velocity dispersion of
Rayleigh waves in polycrystalline materials 25, 26. These studies have
demonstrated the power of the numerical methods (particularly, the FE method)
in realistically simulating Rayleigh wave scattering.

In comparison to the existing studies, this work sets out to study a more
general case of Rayleigh wave scattering, with the scatterer being a
surface/subsurface elastic inclusion with an arbitrary shape and
isotropic/anisotropic property. To achieve this aim, this work contributes to
two aspects. First, this work develops a theoretical model based on the
reciprocity theorem, utilising the Green’s function and the Born
approximation, to calculate the backscattering of Rayleigh waves from an
arbitrarily-shaped surface/subsurface elastic inclusion with an arbitrary
symmetry. Our model is valid for a general elastic inclusion and is not
limited to buried scatterers as reported by a similar prior work 13. Second,
this work also makes use of the proven capability of the FE method to
realistically simulate the same scattering problem. This allows for relatively
accurate results to be obtained, enabling the validation of our general
theoretical model.

The work is organised as follows. Secs. 2 and 3 describe respectively the
theoretical and FE models for the backscattering of Rayleigh waves by an
elastic inclusion. Sec. 4 compares the results of the theoretical and FE
models for a variety of surface and subsurface inclusions of different shapes
and elastic anisotropies. Sec. 5 concludes this work.

##  2  Theoretical model

We consider an isotropic solid with density $\rho_{0}$ and elastic tensor
$c_{pjkl}^{0}$ in the two-dimensional (2D) half-space defined by the $x-z$
coordinates. As shown in Fig. 1, an arbitrarily-shaped inclusion is present on
the surface or subsurface of the host material. The inclusion is defined in
the region $V$ surrounded by the boundary $S$. The inclusion has density
$\rho_{1}\left(\mathbf{x}_{\mathrm{s}}\right)$ and isotropic/anisotropic
property described by the elastic tensor
$c_{pjkl}^{1}\left(\mathbf{x}_{\mathrm{s}}\right)$. A Rayleigh wave
propagating in the host material will be scattered as it encounters the
inclusion. Two types of scattering arises, one from the incident Rayleigh wave
into the same mode and another from the Rayleigh wave into bulk waves. Here we
only consider the Rayleigh-to-Rayleigh scattering as it is more prominent than
the other type, as proved by prior work 3 and supported by our simulation
results in Sec. 3.

Figure 1:  Schematic of (a) theoretical and (b) finite element models. The
yellow dots in (b) are the source line for generating the Rayleigh wave. The
originally generated signal is monitored at the transmitting node T and the
backscattered signal from the surface/subsurface inclusion A is monitored at
the receiving node R. Structured mesh with uniform square elements is used.
$Z_{0}$ is the distance between the upper surface and inclusion. The absorbing
boundary is applied on the left, right and bottom of the surface.

Now we develop a theoretical model to describe the Rayleigh-to-Rayleigh
scattering, and we start from the reciprocity theorem that gives the scattered
Rayleigh wave by15, 9

$\begin{split}{u}^{\text{sc}}_{n}\left(\mathbf{x},\omega\right)&=\int_{S}[c_{pjkl}^{0}n_{l}\left(\mathbf{x}_{\mathrm{s}}\right){u}_{k}^{+}\left(\mathbf{x}_{\mathrm{s}},\omega\right)G_{nj,p}\left(\mathbf{x},\mathbf{x}_{\mathrm{s}},\omega\right)\\\
&\quad\quad-
c_{pjkl}^{0}n_{p}\left(\mathbf{x}_{\mathrm{s}}\right)G_{nj}\left(\mathbf{x},\mathbf{x}_{\mathrm{s}},\omega\right){u}_{k,l}^{+}\left(\mathbf{x}_{\mathrm{s}},\omega\right)]\mathrm{\leavevmode\nobreak\
d}S\textrm{,}\end{split}$ (1)

where the Einstein summation convention over repeated indices ($p$, $j$, $k$,
$l$) from 1 to 3 (or $x$ to $z$) is assumed. $\mathbf{x}$ is the point of
interest at which we evaluate the scattered wave, while
$\mathbf{x}_{\mathrm{s}}$ is a point in the inclusion $V$.
${u}_{l}(\mathbf{x}_{\mathrm{s}})$ is the $l$-th displacement component of the
Rayleigh wave. $G_{nj}(\mathbf{x},\mathbf{x}_{\mathrm{s}})$ is the dyadic
representation of the Green’s function. $n_{k}(\mathbf{x}_{\mathrm{s}})$ is
the $k$-th component of the outward unit normal to the inclusion surface. The
comma derivative notation is used and the derivative is over the scattered
coordinate throughout this work; e.g. ${u}_{k,l}=\partial{u}_{k}/\partial
x_{\mathrm{s}l}$ with $x_{\mathrm{s}l}$ being the $l$-th component of
$\mathbf{x}_{\mathrm{s}}$. The plus and minus superscripts are used to
indicate quantities that are evaluated on the host or inclusion side of the
surface $S$, respectively. From the continuity of displacement and traction
across $S$, we have

${u}_{k}^{+}\left(\mathbf{x}_{\mathrm{s}},\omega\right)={u}_{k}^{-}\left(\mathbf{x}_{\mathrm{s}},\omega\right)\textrm{,}\quad
n_{p}c_{pjkl}^{0}{u}_{k,l}^{+}\left(\mathbf{x}_{\mathrm{s}},\omega\right)=n_{p}c_{pjkl}^{1}{u}_{k,l}^{-}\left(\mathbf{x}_{\mathrm{s}},\omega\right)\textrm{,}$
(2)

so Eq. 1 becomes

$\begin{split}{u}^{\text{sc}}_{n}\left(\mathbf{x},\omega\right)&=\int_{S}[c_{pjkl}^{0}n_{l}\left(\mathbf{x}_{\mathrm{s}}\right){u}_{k}^{-}\left(\mathbf{x}_{\mathrm{s}},\omega\right)\partial
G_{nj,p}\left(\mathbf{x},\mathbf{x}_{\mathrm{s}},\omega\right)\\\ &\quad\quad-
c_{pjkl}^{1}\left(\mathbf{x}_{\mathrm{s}}\right)n_{p}\left(\mathbf{x}_{\mathrm{s}}\right)G_{nj}\left(\mathbf{x},\mathbf{x}_{\mathrm{s}},\omega\right){u}_{k,l}^{-}\left(\mathbf{x}_{\mathrm{s}},\omega\right)]\mathrm{\leavevmode\nobreak\
d}S\textrm{.}\end{split}$ (3)

Applying the divergence theorem $\int_{V}v_{i,k}dV=\int_{S}v_{i}n_{k}dS$ to
the right side of Eq. 3 leads to

$\begin{split}{u}^{\text{sc}}_{n}\left(\mathbf{x},\omega\right)&=\int_{V}\\{\left[c_{pjkl}^{0}{u}_{k}\left(\mathbf{x}_{\mathrm{s}},\omega\right)G_{nj,p}\left(\mathbf{x},\mathbf{x}_{\mathrm{s}},\omega\right)\right]\textrm{,}_{l}\\\
&\quad\quad-\left[c_{pjkl}^{1}G_{nj}\left(\mathbf{x},\mathbf{x}_{\mathrm{s}},\omega\right){u}_{k,l}\left(\mathbf{x}_{\mathrm{s}},\omega\right)\right]\textrm{,}_{p}\\}\mathrm{\leavevmode\nobreak\
d}V\textrm{,}\end{split}$ (4)

where we have dropped the minus superscript since the integration points are
all within the inclusion $V$. Now substituting the equation of motion9

$\begin{split}\left[c_{pjkl}^{0}G_{nj,p}\left(\mathbf{x},\mathbf{x}_{\mathrm{s}},\omega\right)\right]\textrm{,}_{k}+\rho_{0}\omega^{2}G_{nl}\left(\mathbf{x},\mathbf{x}_{\mathrm{s}},\omega\right)&=-\delta_{nl}\delta\left(\mathbf{x}_{\mathrm{s}}-\mathbf{x}\right)\textrm{,}\\\
\left[c_{pjkl}^{1}\left(\mathbf{x}_{\mathrm{s}},\omega\right){u}_{k,l}\left(\mathbf{x}_{\mathrm{s}},\omega\right)\right]\textrm{,}_{p}&=-\rho_{1}\omega^{2}{u}_{j}\left(\mathbf{x}_{\mathrm{s}},\omega\right)\textrm{,}\end{split}$
(5)

into Eq. 4 and using the sampling properties of the delta function, we have

$\begin{split}{u}^{\text{sc}}_{n}\left(\mathbf{x},\omega\right)=\int_{V}&\left[\omega^{2}\Delta\rho
G_{nl}\left(\mathbf{x},\mathbf{x}_{\mathrm{s}},\omega\right){u}_{l}\left(\mathbf{x}_{\mathrm{s}},\omega\right)\right.\\\
&\left.-\Delta
c_{pjkl}(\mathbf{x}_{\mathrm{s}})G_{nj,p}\left(\mathbf{x},\mathbf{x}_{\mathrm{s}},\omega\right){u}_{k,l}\left(\mathbf{x}_{\mathrm{s}},\omega\right)\right]\mathrm{\leavevmode\nobreak\
d}V\textrm{,}\end{split}$ (6)

with

$\begin{gathered}\begin{aligned}
\Delta\rho&=\rho_{1}(\mathbf{x}_{\mathrm{s}})-\rho_{0}\quad\textrm{and}\quad\Delta
c_{pjkl}(\mathbf{x}_{\mathrm{s}})&=c_{pjkl}^{1}(\mathbf{x}_{\mathrm{s}})-c_{pjkl}^{0}\textrm{.}\end{aligned}\end{gathered}$
(7)

Equation 6 provides an exact solution for the scattering amplitude, but the
solution is intractable. We address this difficulty by invoking the Born
approximation with assuming a small property difference between the inclusion
and host materials. Here we do not attempt to define how small the property
difference needs to be for obtaining a reasonably accurate solution, but as we
shall see in Sec. 4, it depends on various factors such as frequency,
inclusion size and material anisotropy. As a result of the Born approximation,
the incident wave is only slightly perturbed by the inclusion, and therefore
the displacement and its derivative in Eq. 6 can be obtained from their values
due to the incident wave only 9, i.e.

$\begin{gathered}{u}_{l}(\mathbf{x}_{\mathrm{s}})\approx{u}_{l}^{\text{in}}(\mathbf{x}_{\mathrm{s}})\quad\textrm{and}\quad{u}_{k,l}(\mathbf{x}_{\mathrm{s}})\approx{u}_{k,l}^{\text{in}}(\mathbf{x}_{\mathrm{s}})\textrm{.}\end{gathered}$
(8)

The formulation so far is valid for both two-dimensional or three-dimensional
cases. Since the 3D dyadic Green’s function is not known yet, here we only
address the 2D case. Now, we assume that the incident Rayleigh wave is a time-
harmonic plane wave propagating along the surface in the $x$-direction, as
shown in Fig. 1(a). In this case, the unit displacement components and their
derivatives can be given by 8, 27, 28

$\begin{split}{u}_{k}^{\text{in}}(\mathbf{x}_{\mathrm{s}})={d}_{k}^{\text{in}}(z_{\mathrm{s}})\exp{\left(\text{i}k_{R}\mathbf{e}^{\text{in}}\cdot\mathbf{x}_{\mathrm{s}}\right)}\textrm{,}&\quad{u}_{k,l}^{\text{in}}(\mathbf{x}_{\mathrm{s}})=\left[{d}_{k,l}^{\text{in}}(z_{\mathrm{s}})+\text{i}{d}_{k}^{\text{in}}(z_{\mathrm{s}})k_{R}{e}^{\text{in}}_{l}\right]\exp{\left(\text{i}k_{R}\mathbf{e}^{\text{in}}\cdot\mathbf{x}_{\mathrm{s}}\right)}\textrm{,}\\\
\end{split}$ (9)

with

$\begin{split}&{d}^{\text{in}}_{k}=\left[U_{R}(z_{\mathrm{s}}),0,\text{i}W_{R}(z_{\mathrm{s}})\right]\textrm{,}\\\
U_{R}(z_{\mathrm{s}})&=\frac{k_{R}}{p}\frac{(2c_{T}^{2}-c_{R}^{2})}{2c_{T}^{2}}\exp(-pz_{\mathrm{s}})-\frac{q}{k_{R}}\exp(-qz_{\mathrm{s}})\textrm{,}\\\
W_{R}(z_{\mathrm{s}})&=\frac{(2c_{T}^{2}-c_{R}^{2})}{2c_{T}^{2}}\exp(-pz_{\mathrm{s}})-\exp(-qz_{\mathrm{s}})\textrm{,}\end{split}$
(10)

where $p=k_{R}\sqrt{1-c_{R}^{2}/c_{L}^{2}}$ and
$q=k_{R}\sqrt{1-c_{R}^{2}/c_{T}^{2}}$. $\mathbf{e}^{\text{in}}=[1,0,0]$ is the
propagation direction of the incident Rayleigh wave. $k_{R}$ and $c_{R}$ are
the wave number and phase velocity of the incident Rayleigh wave. The phase
velocity $c_{R}$ can be calculated by

$(2-c_{R}^{2}/c_{T}^{2})^{2}-4(1-c_{R}^{2}/c_{L}^{2})^{1/2}(1-c_{R}^{2}/c_{T}^{2})^{1/2}=0\textrm{,}$
(11)

where $c_{L}$ and $c_{T}$ are the velocities of the longitudinal and shear
waves in the host material.

For the 2D case considered here, the dyadic Green’s function is given by 13,
29

$\begin{split}G_{nl}\left(\mathbf{x},\mathbf{x}_{\mathrm{s}},\omega\right)=A_{0}\frac{\exp(\text{i}k_{R}r_{2})}{\sqrt[]{2\pi
r_{2}/k_{R}}}[{d}_{l}^{\text{sc}}\left(z_{\mathrm{s}}\right)]^{*}{p}^{\text{sc}}_{n}\left(z\right)\textrm{,}\end{split}$
(12)

with

$\begin{split}&\quad\quad\quad
A_{0}=\frac{1}{4P_{R}c_{R}}\frac{\exp(\text{i}\pi/4)}{k_{R}},\quad\\\
{p}^{\text{sc}}_{n}\left(z\right)=&\left[U_{R}\left(z\right),0,\text{i}W_{R}\left(z\right)\right]\textrm{,}\quad{d}_{l}^{\text{sc}}\left(z_{\mathrm{s}}\right)=\left[-U_{R}(z_{\mathrm{s}}),0,\text{i}W_{R}(z_{\mathrm{s}})\right]\textrm{,}\\\
&P_{R}=\frac{1}{2}\rho_{0}{c}_{g}\int_{o}^{\infty}\left[U_{R}\left(z\right)^{2}+W_{R}\left(z\right)^{2}\right]\mathrm{\leavevmode\nobreak\
d}z\textrm{,}\end{split}$ (13)

where $P_{R}$ represents a normalized power per unit width in the travelling
wave mode. ${c}_{g}$ is the group velocity of the Rayleigh wave. The asterisk
superscript denotes the complex conjugate. In the far field, the distance
$r_{2}$ between the evaluation point $\mathbf{x}$ and a point
$\mathbf{x}_{\mathrm{s}}$ on the inclusion can be approximated as 15, 9

$r_{2}=\mathbf{x}-\mathbf{x_{\mathrm{s}}}\approx
r-\mathbf{e}^{\text{sc}}\cdot\mathbf{x}_{\mathrm{s}}$ (14)

where a fixed point on the inclusion is taken as the origin of the
$\left(x,z\right)$ coordinates. $\mathbf{e}^{\text{sc}}$ is the 2D unit vector
in the $\left(x,z\right)$ plane from the origin $O$ to point $\mathbf{x}$,
representing the propagation direction of the scattered Rayleigh wave. $r$ is
the distance between $O$ and $\mathbf{x}$, i.e., $r=\sqrt{x^{2}+z^{2}}$.

The Green’s function in Eq. 12 applies to any scattering direction. For
simplicity, we only consider the scattering in the backward direction of the
incident wave, namely $\mathbf{e}^{\text{sc}}=[-1,0,0]$. Then, the Green’s
function and its derivative can be written as

$\begin{split}G_{nl}\left(\mathbf{x},\mathbf{x}_{\mathrm{s}},\omega\right)&=A_{0}\frac{\exp(\text{i}k_{R}r)}{\sqrt[]{r}}[{d}_{l}^{\text{sc}}\left(z_{\mathrm{s}}\right)]^{*}{p}^{\text{sc}}_{n}\left(z\right)\exp(-\text{i}k_{R}\mathbf{e}^{\text{sc}}\cdot\mathbf{x}_{\mathrm{s}})\textrm{,}\\\
G_{nj,p}\left(\mathbf{x},\mathbf{x}_{\mathrm{s}},\omega\right)=&A_{0}\frac{\exp(\text{i}k_{R}r)}{\sqrt[]{r}}\left\\{[{d}_{j,p}^{\text{sc}}\left(z_{\mathrm{s}}\right)]^{*}-\text{i}k_{R}{e}^{\text{sc}}_{p}[{d}_{j}^{\text{sc}}\left(z_{\mathrm{s}}\right)]^{*}\right\\}\\\
&\quad\times{p}^{\text{sc}}_{n}\left(z\right)\exp{\left(-\text{i}k_{R}\mathbf{e}^{\text{sc}}\cdot\mathbf{x}_{\mathrm{s}}\right)}\textrm{.}\end{split}$
(15)

Substituting Eq. 15 into Eq. 6 and rearranging the result, we have

${u}^{\text{sc}}_{n}\left(\mathbf{x},\omega\right)=A^{\text{sc}}\left(\omega\right)\frac{\exp(\text{i}k_{R}r)}{\sqrt[]{2\pi
r/k_{R}}}{p}^{\text{sc}}_{n}\left(z\right)\textrm{,}$ (16)

where $A^{\text{sc}}\left(\omega\right)$ is the far-field amplitude of the
backscattered Rayleigh wave, given by

$\begin{split}A^{\text{sc}}\left(\omega\right)=A_{0}\int_{V}&\left\\{\Delta\rho\omega^{2}d^{\text{in}}_{l}\left(\mathbf{x}_{\mathrm{s}}\right)[{d}^{\text{sc}}_{l}\left(\mathbf{x}_{\mathrm{s}}\right)]^{*}-\Delta
c_{pjkl}(\mathbf{x}_{\mathrm{s}})N(\mathbf{x}_{\mathrm{s}})M(\mathbf{x}_{\mathrm{s}})\right\\}\\\
&\times\exp{\left[\text{i}k_{R}\left(\mathbf{e}^{\text{in}}-\mathbf{e}^{\text{sc}}\right)\cdot\mathbf{x}_{\mathrm{s}}\right]}\mathrm{\leavevmode\nobreak\
d}V\textrm{,}\end{split}$ (17)

with

$N(\mathbf{x}_{\mathrm{s}})={d}_{k,l}^{\text{in}}\left(z_{\mathrm{s}}\right)+\text{i}k_{R}{e}^{\text{in}}_{l}{d}_{k}^{\text{in}}\left(z_{\mathrm{s}}\right)\quad\textrm{and}\quad
M(\mathbf{x}_{\mathrm{s}})=[{d}_{j,p}^{\text{sc}}\left(z_{\mathrm{s}}\right)]^{*}-\text{i}k_{R}{e}^{\text{sc}}_{p}[{d}_{j}^{\text{sc}}\left(z_{\mathrm{s}}\right)]^{*}\textrm{.}$
(18)

Eq. 17 can be further evaluated and the resulting final expression is provided
in Appendix A. Note that Einstein summation convention over repeated indices
($p$, $j$, $k$, $l$) is only for 1 and 3 (or $x$ to $z$) in the 2D case.

Although Eq. 17 is derived for isotropic inclusions, it can be modified to
accommodate anisotropic inclusions if the root-mean-square (RMS) of the
backscattering response from a lot of randomly-oriented inclusions is
concerned. Let us consider a case that the host material is the Voigt average
of the anisotropic inclusion, with the elastic constants given by

$\begin{split}c_{11}^{0}=\left<c_{11}\right>&=\frac{3(c_{11}+c_{22}+c_{33})+2(c_{23}+c_{13}+c_{12})+4(c_{44}+c_{55}+c_{66})}{15}\\\
c_{44}^{0}=\left<c_{44}\right>&=\frac{(c_{11}+c_{22}+c_{33})-(c_{23}+c_{13}+c_{12})+3(c_{44}+c_{55}+c_{66})}{15}\textrm{,}\end{split}$
(19)

The fourth-rank elastic tensor $c_{pjkl}$ is written as $c_{ij}$ using the
Voigt index notation where the pairs of indices are contracted to the
following single values: $11\rightarrow 1$, $22\rightarrow 2$, $33\rightarrow
3$, $23\,\mathrm{or}\,32\rightarrow 4$, $13\,\mathrm{or}\,31\rightarrow 5$ and
$12\,\mathrm{or}\,21\rightarrow 6$. The host material has the same density as
the inclusion, namely $\rho_{0}=\rho_{1}$. In this case, the RMS of the
backscattering amplitude from an infinite number of inclusions with random
crystallographic orientations can be calculated by

$\begin{split}A^{\text{sc}}_{rms}\left(\omega\right)=\sqrt{\left<A_{pjkl}^{\text{sc}}\left(\omega\right)A^{sc}_{\alpha\beta\gamma\delta}\left(\omega\right)\right>}=\sqrt{\left<\sum_{pjkl}^{1,3}A_{pjkl}^{\text{sc}}\left(\omega\right)\sum_{\alpha\beta\gamma\delta}^{1,3}A^{sc}_{\alpha\beta\gamma\delta}\left(\omega\right)\right>}\textrm{.}\end{split}$
(20)

The developed theoretical model calculates the far-field amplitude of the
backscattered Rayleigh wave from an arbitrarily-shaped inclusion in 2D. The
size, shape and depth of the inclusion, and its property contrast to the host
material are all incorporated in the integral in Eqs. 17 and 20. For the
inclusions considered in this work, the integral is evaluated by numerical
integration for both regularly and irregularly shaped inclusions. A major
advantage of the model is that it is capable of dealing with arbitrarily-
shaped surface/subsurface inclusions of isotropic or general anisotropic
properties. The use of the Born approximation, however, means that the
property contrast between the inclusion and the host material needs to be
necessarily small; thus, the traditional Kirchhoff approximation is preferable
when the scatterer is a void or crack. We will use numerical simulations to
evaluate how the Born approximation affects the accuracy of the obtained
solution.

##  3 Finite element model

The FE method has been demonstrated recently to be powerful and accurate for
simulating the propagation and scattering of bulk 30, 31, 32, 33, 34, 21, 35
and surface 36, 26 waves in complex solids. Here we report our use of the
method for simulating the same physical problem as addressed by the above
theoretical model, which is the backscattering of a Rayleigh wave by a
surface/subsurface inclusion. Note that we are concerned with the wave
phenomena in the middle of a wide wavefront of a propagating Rayleigh wave,
and the symmetry means that this can be simplified and well simulated by a 2D
FE model.

As schematically shown in Fig. 1, the 2D FE model is based in the $x-z$ plane.
The dimension of the model depends on the modelling frequency and inclusion
type, and detailed parameters are given in Table 1 for the modeled cases. The
model space is discretized with uniform linear square elements, with an edge
size $h$ of about one sixtieth of the center-frequency wavelength of the
Rayleigh wave in the host material to minimize numerical error to $\sim 0.1\%$
30, 31, 35. The inclusion is modelled by an aggregate of elements that are
assigned with the density $\rho_{1}$ and elastic tensor $c_{pjkl}^{1}$ of the
inclusion material. The remaining elements in the middle of the model are
defined as the host material with density $\rho_{0}$ and elastic tensor
$c_{pjkl}^{0}$. The inclusions considered in this work are summarized in Table
2. For isotropic inclusions, the host material is defined as aluminum with
Young’s modulus $E_{0}$ = $70$ GPa, Poisson’s ratio $\nu_{0}$ = $0.35$ and
density $\rho_{0}$ = 2700 kg/m3, and the inclusion has 4% impedance contrast
to the host material caused by density or/and Young’s modulus differences. For
anisotropic inclusions, the host material has the same density as the
inclusion (i.e., $\rho_{0}=\rho_{1}$), and its isotropic elastic tensor
$c_{pjkl}^{0}$ is the Voigt average of $c_{pjkl}^{1}$ calculated by Eq. 19.
The anisotropic inclusion materials considered in this work are listed in
Table 3.

Table 1:  Models used in the simulation. Center frequency of FE modelling $f_{c}$ (MHz), dimensions $d_{x}\times d_{z}$ (mm $\times$ mm), mesh size $h$ (mm), degree of freedom (d.o.f); copper sulfate pentahydrate is abbreviated as CSP. Material | $f_{c}$ | $d_{x}\times d_{z}$ | $h$ | d.o.f |
---|---|---|---|---|---
Aluminum | $0.5$ |  560 $\times$ 56 |  96 $\times$ 10-3 |  6.8 $\times$ 106 |
$1$ |  280 $\times$ 28 |  48 $\times$ 10-3 |  6.8 $\times$ 106 |
Aluminum | $2$ |  140 $\times$ 14 |  24 $\times$ 10-3 |  6.8 $\times$ 106 |
Inconel, Lithium | $2$ |  160 $\times$ 16 |  24 $\times$ 10-3 |  8.9 $\times$ 106 |
CSP | $2$ |  120 $\times$ 12 |  18 $\times$ 10-3 |  8.9 $\times$ 106 |
Aluminum | $4$ |  70 $\times$ 7 |  12 $\times$ 10-3 |  6.8 $\times$ 106 |
Inconel, Lithium | $4$ |  80 $\times$ 8 |  12 $\times$ 10-3 |  8.9 $\times$ 106 |
CSP | $4$ |  60 $\times$ 6 |  9 $\times$ 10-3 |  8.9 $\times$ 106 |
Aluminum | $8$ |  35 $\times$ 3.5 |  6 $\times$ 10-3 |  6.8 $\times$ 106 |
Inconel, Lithium | $8$ |  40 $\times$ 4 |  6 $\times$ 10-3 |  8.9 $\times$ 106 |
CSP | $8$ |  30 $\times$ 3 |  4.5 $\times$ 10-3 |  8.9 $\times$ 106 |
Aluminum | $16$ |  17.5 $\times$ 1.75 |  3 $\times$ 10-3 |  6.8 $\times$ 106 |
$20$ |  14 $\times$ 1.4 |  2.5 $\times$ 10-3 |  6.8 $\times$ 106 |
$32$ |  8.75 $\times$ 0.875 |  1.5 $\times$ 10-3 |  6.8 $\times$ 106 |
Table 2:  Inclusion shape and size used in the simulations. A regularly-shaped inclusion has depth $D$ (mm) and width $L$ (mm), and an irregular inclusion has an equivalent radius, $d=\sqrt{s/\pi}$, with $s$ being the area of the inclusion. | Inclusion | Illustration | $D$ | $L$ | $d$ | $Z_{0}/\lambda$ |
---|---|---|---|---|---|---|---
Regular | Half-circle |  | $0.19$ | $D$ | - | 0 |
half-ellipse |  |  4$D$ | - |
Square |  | $D$ | - |
Rectangle |  |  2$D$ | - |
Circle |  | 0.38 | $D$ | - |  0 $\sim$ 2 |
Irregular |  Inclusion $M$ |  | - | - | 0.35 | 0 |
Table 3:  Properties of anisotropic inclusion materials. Density $\rho$ (kg/m3), equivalent anisotropy index37, 38 $A^{eq}$, elastic constants $c_{ij}$ (GPa), and Voigt averages $c_{11}^{0}$ and $c_{44}^{0}$ (GPa). | Material | $\rho$ | $A^{eq}$ | $c_{11}$ | $c_{12}$ | $c_{44}$ | $c_{11}^{0}$ | $c_{44}^{0}$ |  |  |  |  |  |  |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
Cubic | Aluminum | 2700 | 1.24 | 106.7 | 60.4 | 28.3 | 110.8 | 26.2 |  |  |  |  |  |  |
Inconel | 8260 | 2.83 | 234.6 | 145.4 | 126.2 | 299.9 | 96.6 |  |  |  |  |  |  |
Lithium | 534 | 9.14 | 13.4 | 11.3 | 9.6 | 20.4 | 6.18 |  |  |  |  |  |  |
| Material | $\rho$ | $A^{eq}$ | $c_{11}$ | $c_{12}$ | $c_{13}$ | $c_{14}$ | $c_{15}$ | $c_{16}$ | $c_{22}$ | $c_{23}$ | $c_{24}$ | $c_{25}$ | $c_{26}$ |
Triclinic | CSP | 2286 | 2.37 | 56.5 | 26.5 | 32.1 | -3.3 | -0.8 | -3.9 | 43.3 | 34.7 | -0.7 | -2.1 | 2.0 |
$c_{33}$ | $c_{34}$ | $c_{35}$ | $c_{36}$ | $c_{44}$ | $c_{45}$ | $c_{46}$ | $c_{55}$ | $c_{56}$ | $c_{66}$ | $c_{11}^{0}$ | $c_{44}^{0}$ |  |
56.9 | -4.4 | -2.1 | -1.6 | 17.3 | 0.9 | 0.3 | 12.2 | -2.6 | 10.0 | 54.3 | 12.2 |  |

The elements on the bottom, left, and right sides of the model are used to
define absorbing boundary conditions. The thickness of each absorbing boundary
region in the boundary normal direction is chosen to be at least three times
the wavelength of the Rayleigh wave in the host material. The absorbing
boundary elements in the vicinity of the host material have material
properties close to those of the host material but their material damping
increases gradually towards the model edge. This gradual increase of damping
helps absorb the propagating wave so as to minimise unwanted reflections from
the boundaries 39.

The desired Rayleigh wave is generated by applying two sinusoidal time-domain
signal of 90∘ phase shift to multiple source nodes located on the top surface
of the model (yellow points in Fig. 1). The size of the source is set to be
equal to three center-frequency wavelengths of the simulated Rayleigh wave,
and each source node is assigned a unique amplitude following Eq. (17) in
Sarris et al.36. The simulation is solved using the GPU-accelerated Pogo
program 40 with an explicit time-stepping scheme. A relatively large time step
of $\Delta t=0.9h/c_{L}$, satisfying the Courant-Friedrichs-Lewy condition 41,
is used to minimise numerical error 35.

Over the course of the FE solution, the $z$ \- displacement of the generated
incident wave is monitored at a transmitting node (point T in Fig. 1), while
that of the backscattered wave is recorded at a receiving node (point R). We
emphasize that the transmitting and receiving nodes are placed respectively
far away from the source nodes and the inclusion, in order for the former to
monitor the well-formed incident wave and for the latter to record solely the
scattered Rayleigh wave in the far field. In addition, a reference signal is
obtained at the receiving point using an identical but inclusion-free FE
model, and the reference signal is subtracted from the relatively small raw
signal to minimise the influence of numerical error. The signal $U_{T}(t)$ at
the transmitting node and the corrected signal $U_{R}(t)$ at the receiving
node are Fourier transformed into the frequency domain to obtain the spectra
$U_{T}(f)$ and $U_{R}(f)$. The frequency-dependent amplitude of the
backscattered Rayleigh wave is then calculated by

$U_{R}(f)=U_{T}(f)A^{\text{sc}}(f)\textrm{,}$ (21)

which will be used to evaluate the theoretical model result,
$A^{\text{sc}}(\omega)$, in Sec. 4. It should be noted that this equation is
only applicable when the attenuation of the host material and the diffraction
losses of the Rayleigh wave are not considered15.

Now we present an example to illustrate the simulated wave field and Rayleigh
wave signals. The host material is the aforementioned isotropic aluminum, and
the half-circle inclusion (Table 2) has a 4% impedance contrast to the host
material. The wave field in the model is shown in Fig. 2(a1) shortly after
exciting the source nodes with a signal of 2 MHz center frequency. The wave
field shows the coexistence of multiple wave modes in the model, led by faster
skimming longitudinal, bulk longitudinal and head waves, and followed by
slower bulk shear and Rayleigh waves 24. However, as the waves propagate
further to the transmitting point, the relatively slow Rayleigh wave gets
separated from other wave modes and a rather pure Rayleigh wave is formed, as
can be seen in Fig. 2(a2). Similarly, the inclusion causes multiple
backscattered wave modes, but as illustrated in Fig. 2(a3), a relatively pure
Rayleigh wave is obtained as the waves reach the receiving point. In this
specific case, the scattered Rayleigh wave is about 100 times stronger than
the scattered longitudinal wave, confirming our choice of focusing on the
scattered Rayleigh wave only in this work. The Rayleigh wave signals recorded
at the transmitting and receiving nodes are plotted in Fig. 2(b), and the
respective frequency-domain amplitude spectra are displayed in Fig. 2(c). The
calculated backscattering amplitude will be reported in Sec. 4.

Figure 2:  Example FE modelling of Rayleigh wave scattering by a surface half-
circle inclusion. (a) The FE model setup and the simulated wave fields at (a1)
the point P at an early time of $t=10\,\mu s$, (a2) the point close to the
inclusion before the incident wave is scattered by the inclusion at $t=35\,\mu
s$ , and (a3) the point far away from the inclusion after the wave scattered
by the inclusion at $t=69\,\mu s$. (b) The $z-$ displacements of the incident
Rayleigh wave and the backscattered Rayleigh wave in the time domain. (c)
shows the respective amplitude spectra in the frequency domain.

##  4 Results and discussions

###  4.1 Simple inclusion with different material properties

It is demonstrated by the theoretical model in Eq. 17 that both the density
difference $\Delta\rho$ and elastic tensor contrast $\Delta c_{pjkl}$ between
the inclusion and host materials contribute to Rayleigh wave backscattering.
As an essential first step, here we investigate their individual
contributions. For this purpose, we consider a simple case of a half-circle
inclusion (see Table 2) on the surface of an aluminum host material. The host
material has Young’s modulus $E_{0}$ = $70$ GPa, Poisson’s ratio $\nu_{0}$ =
$0.35$ and density $\rho_{0}$ = 2700 kg/m3. The inclusion is defined in three
distinctive cases that differ from the host material in (1) density with
$\rho_{1}$ = 2500 kg/m3, (2) Young’s modulus with $E_{1}$ = $65$ GPa, and (3)
both density and Young’s modulus with $\rho_{1}$ = 2600 kg/m3 and $E_{1}$ =
$67$ GPa. These three cases have the same acoustic impedance mismatch of 4% to
the host material.

For these three inclusion cases, the amplitudes of the backscattered Rayleigh
waves are plotted in Fig. 3 against the normalized frequency $k_{R}D$ ($D$ is
the radius of the inclusion). The theoretical curves in the figure are
calculated from Eq. 17, while the FE points are simulated using the models in
Table 1. We emphasize that we have achieved a high degree of accuracy for the
FE results. A prominent evidence is the overlapping of the points between the
neighboring models having different model parameters (thus different numerical
errors)35. Therefore, the FE results are well suited to evaluate the
approximations of the theoretical model.

Figure 3:  Rayleigh wave backscattering by a simple half-circle inclusion that
differs from the host aluminum in (a) density, (b) Young’s modulus and (c)
both density and Young’s modulus. The resulting backscattering amplitude $A$
is plotted versus normalised frequency $k_{R}D$ ($D$ is the radius of the
inclusion). The theoretical curves calculated from Eq. 17 are compared with
the FE points obtained from the models in Table 1 with center frequencies of
0.5, 1, 2, 4, 8 and 16 MHz. The inclusion and host materials have a 4%
impedance mismatch for the solid theoretical lines and FE points, and 2%
mismatch for the dash-dotted theoretical lines. Note that the y-axis range in
(b) is half of those in (a) and (c).

The backscattering results for the three inclusion cases in Fig. 3 demonstrate
an oscillating, cyclic behavior. For each case, the peaks of individual cycles
increase gradually with $k_{R}D$, but their cycles are seemingly constant
across different cycles. The overall backscattering amplitude for the first
case of $\Delta\rho\neq 0$ is about twice for the second case of $\Delta
c_{pjkl}\neq 0$, while the third case resides in the middle. Similarly, the
average cycle of the first case is slightly larger than that of the second,
again with the third case lying in between. Considering the differences
between the three cases, we can infer two main results from the theoretical
model in Eq. 17.

First, $\Delta\rho$ and $\Delta c_{pjkl}$ are scaling factors affecting only
the magnitude of backscattering. This is further corroborated by the same
cyclic behavior but different amplitudes of the two theoretical curves in each
of Fig. 3(a) and (b); the extra theoretical curve in each plot is obtained
using a smaller impedance mismatch of 2% by varying $\Delta\rho$ or $\Delta
c_{pjkl}$.

Second, the remaining integral terms in the equation affect not only the
magnitude but also the cyclic period of the backscattering. We can observe
from the two cases in Fig. 3(a) and (b) that, for a given host material and
inclusion, the two integral terms associated with $\Delta\rho$ and $\Delta
c_{pjkl}$ exhibit different cycles. This is also the reason for the small
shift of around 2% between the two theoretical curves in Fig. 3(c). In
addition, we shall see in the subsection below that the cyclic behavior of the
term associated with $\Delta c_{pjkl}$ (presumably for the term associated
with $\Delta\rho$ as well) is much more significantly affected by the geometry
of the inclusion.

For the first two inclusion cases in Fig. 3(a) and (b), the theoretical model
predictions exhibit very good agreement with the FE results at small $k_{R}D$.
As $k_{R}D$ increases, the theoretical curve in Fig. 3(a) tends to have a
noticeably shorter cycle (as if compressed) than the FE results, while this
trend is reversed for the case in Fig. 3(b) with the theoretical curve being
seemingly expanded. As a result, the agreement between the theoretical and FE
results deteriorates as $k_{R}D$ increases. This is particularly evident as we
observe their $k_{R}D$ differences at individual peaks: the difference
increases from 0.1% at the first peak to 3.7% at the fifth peak for the first
case, and it grows from 0.1% to 3.5% at the same peaks for the second case.
Interestingly, for the third case with both $\Delta\rho\neq 0$ and $\Delta
c_{pjkl}\neq 0$, the theoretical model agrees well with the FE results even at
a very large $k_{R}D$, with their difference observable only after $k_{R}D$
$\approx 11$.

Such differences can be understood by investigating the phase change caused by
each inclusion in comparison to the case where the inclusion is absent. The
Born approximation tends to have a larger deviation from the true value when
the inclusion-induced phase change gets bigger 42, 43. In our cases, the
analytically estimated phase changes are $-0.22\pi$, $0.20\pi$ and $0.02\pi$
for Fig. 3(a), (b) and (c) at $k_{R}D$ = 9, which are clear evidence
supporting the observed large (and similar) theoretical-FE differences in the
first two cases and the good agreement in the third one. Apparently, the
inclusion-induced phase change increases with frequency, which leads to the
increased theoretical-FE difference as observed in each figure panel.

We should emphasize that the three inclusions considered have a smaller
density and/or Young’s modulus than the host material. We also study the three
opposite cases (results not shown) with the density and/or Young’s modulus of
the inclusion being larger than those of the host. In comparison to Fig. 3,
their FE results are different but their theoretical predictions remain the
same (obvious from Eq. 17); specifically, the theoretical curves are scaled
(compressed/expanded over $k_{R}D$) with respect to the FE results in an
opposite way to those in Fig. 3.

For simplicity, we shall consider only the second case of $\Delta c_{pjkl}\neq
0$ in the next two subsections. This is because this particular case involves
a considerable approximation in the theoretical solution and is thus
beneficial for us to thoroughly evaluate the theoretical model. In this case,
the inclusion and host materials will have the same density but different
elastic constants.

###  4.2 Isotropic surface inclusion with different shapes

Now we compare the backscattering amplitudes from the four regularly-shaped
inclusions in Table 2. The host material is aluminum (isotropic and with the
above-mentioned material properties) and the inclusion has a 4% impedance
difference to the host caused by $\Delta c_{pjkl}\neq 0$. The FE results
(points) and theoretical predictions (lines) are displayed in Figs.
4(a1)-(a4), plotted against $k_{R}D$. Note that Figure 4(a1) is the same as
Fig. 3(b). All four cases show a good agreement between the theoretical and
numerical results when $k_{R}D$ is small, which is excellent proof of the
validity of the theoretical model.

When comparing Fig. 4(a1) and (a2) (or (a3) and (a4)), we observe that the
results exhibit different cyclic behaviors depending on the lateral dimension
of the inclusion. This is further demonstrated in Figures 4(b1) and (b2),
which reveals that the lateral-to-depth dimension ratio $L/D$ affects the
cycle of the backscattering amplitude curve and a larger $L/D$ ratio
corresponds to an apparently smaller average cycle for both half-ellipse and
rectangle inclusions. This can be explained by that the $z-$ displacement of
Rayleigh waves is non-uniform, and the energy of Rayleigh wave is becoming
smaller with the increase of depth. Meanwhile, by comparing Fig. 3 with Fig.
4, it can be seen that the geometry of inclusion has a more obvious effect on
the cyclic behaviour.

Figure 4:  Rayleigh wave backscattering by a (a1) half-circle , (a2) half-
ellipse , (a3) square and (a4) rectangle inclusion. The theoretical results
(curves) are compared with the FE results (points) with center frequencies of
0.5, 1, 2, 4, 8, 16 MHz. The inclusion is defined to be slightly different
from the host material aluminum only in Young’s modulus ($E_{1}$ = $65$). (b)
represents the relationship between the size of inclusion $L/D$ and the cycle
of the backscattering amplitude curve for a (b1) half-ellipse and (b2)
rectangle inclusion.

Then, we investigate the backscattering amplitude from an irregular inclusion
(depicted in Fig. 5(a) and listed in Table 2). The theoretical prediction and
the FE results are compared in Fig. 5(b), which show very good agreement
between each other, demonstrating that the theoretical solution is accurate to
describe the backscattered wave of arbitrarily shaped inclusions. As expected
from the applicability of the Born approximation, the agreement is becoming
worse with the increase of $k_{R}D$.

Figure 5: (a) Irregularly shaped inclusions modeled with uniform square
elements; (b) comparison of backscattering amplitude by the irregular
inclusion between the FE solution and the theoretical prediction. The host
material and inclusion material properties are the same as those of the above
regularly shaped inclusion simulations.

###  4.3  Isotropic subsurface inclusion

Furthermore, we conduct research to study the backscattering amplitude of a
subsurface inclusion. As shown in Fig. 6(a), we utilize a circular inclusion
(parameters in Table 2) with a depth $Z_{0}=$ 1.45 mm to the surface of the
host material. We set it with $\Delta\rho$ = 0, Young’s modulus $E_{1}$ = 65
GPa and keep other material properties unchanged with respect to the host
material.

The analytical result obtained by Eq. 17 is compared with the FE solution in
Fig. 6(b). The results are plotted against the normalized frequency $k_{R}D$.
In the studied frequency range, the ratio of the depth $Z_{0}$ to the
wavelength $\lambda_{r}$ covers a range of 0 $\sim$ 2 $\lambda_{r}$, which is
displayed in Figure 6(c). Combining Fig. 6(b) and (c), we can see that the
analytical results match well with the FE results when $Z_{0}/\lambda_{r}$
$<$1, demonstrating a good accuracy of the theoretical model in describing
shallow subsurface inclusions. The analytical solution begins to divert from
the FE result at around $Z_{0}$ $\approx$ $\lambda_{r}$. Their difference
increases with $Z_{0}/\lambda_{r}$, with the theoretical result being 10 times
smaller than the FE results at $Z_{0}/\lambda_{r}$ $\approx$ 2, revealing the
limitation of the Born approximation for deeper inclusions. Given that the
energy of Rayleigh surface waves is generally concentrated in the near-surface
region within a depth of about one wavelength 44, the limitation can be
understandable.

Figure 6:  Rayleigh wave scattering by a subsurface inclusion: (a)
illustration of the inclusion; (b) plot of theoretical and FE results versus
normalized frequency $k_{R}D$ on a log scale; (c) the relationship between
$Z_{0}/\lambda_{r}$ and frequency. The host material and inclusion material
properties are the same as those of the above regularly shaped inclusion
simulations.

###  4.4 Anisotropic inclusion

Here we also evaluate the applicability of our theoretical model to an
anisotropic inclusion. This is to identify if the model can be developed
further to describe Rayleigh wave scattering in a polycrystalline material in
the future. Therefore, for the evaluation, we define the inclusion (half-
circle in this case, Table 2) as a single crystal. Similarly to the
polycrystalline material case, we are interested in the average backscattering
response from the inclusion when its crystallographic axis is differently
oriented, which is exactly what Eq. 20 predicts. In the FE modeling, we
randomly rotate the inclusion for 100 times (realizations) and perform FE
simulation for each rotated inclusion, and then take the RMS of the 100
backscattering results as the final result. The host material has the same
density as the inclusion, and its elastic properties are the Voigt averages of
the single inclusion. The single-crystal elastic constants and their Voigt
averages are provided in Table 3 for the materials considered.

We begin by the case of the inclusion being aluminum, which has a cubic
single-crystal symmetry and an anisotropy index of 1.24. The theoretically
predicted RMS backscattering amplitudes $A_{rms}$ are plotted as the grey
curve in Fig. 7(a). The respective FE results are plotted as points in the
figure, which are the RMS over 100 realizations with the inclusion being
randomly rotated in each realization. The error bars show the 95% confidence
interval 45 for the FE points, demonstrating the variation across the
realizations with different crystallographic orientations. It is important to
note that the use of 100 realizations is sufficient for obtaining
statistically converged FE results. This is evidenced by the two convergence
curves at $k_{R}D$ $\approx$ 0.5 and 5.6 (marked as stars in Fig. 7(a)) in
Figs. 7(a1) and (a2) that show gradually stabilizing RMS value as the number
of realizations increases. A general finding from Fig. 7(a) is that the
theoretical and FE results display a good agreement.

Figure 7:  Rayleigh wave scattering by a half-circle inclusion that is cubic
aluminum with an anisotropy index $A^{eq}$ of 1.24. (a) displays the
theoretically predicted RMS backscattering amplitude (line) and the FE
calculated RMS backscattering amplitude (points) of 100 realizations. The
error bar show the 95% confidence interval45. (a1) and (a2) show the
convergences of the FE RMS value with the number of realizations at $k_{R}D$
$\approx$ 0.5 and 5.6.

This good agreement is not surprising for an inclusion material (aluminum) of
small anisotropy. To observe how the agreement changes with anisotropy, we
compare the results of three cubic inclusion materials, namely aluminum,
Inconel and lithium, that have anisotropy indices of 1.24, 2.83 and 9.14.
Their theoretical and FE results are depicted in Figs. 8(a1)-(a3). It is
clearly demonstrated that the agreement between the theoretical and FE results
decreases as anisotropy index increases. This is particularly evident as we
observe the ‘starting point’ (exhibited as yellow dots in Fig. 8) where the
theory starts to deviate from the FE results, which is 5.62 for aluminum, 1.32
for Inconel and 0.98 for lithium. Such results are reasonable because the Born
approximation is expected to gradually fail with the increase of scattering
intensity.

Figure 8:  Rayleigh wave scattering by a simple half-circle inclusion defined
as cubic materials of (a1) aluminum ($A^{eq}$ = 1.24), (a2) Inconel ($A^{eq}$
= 2.83) and (a3) lithium ($A^{eq}$ = 9.14), and triclinic material of CSP
($A^{eq}$ = 2.37). The yellow points are the ‘starting points’ where the
theoretical and FE results start to deviate. (b) Comparison of the theoretical
curves of the four materials.

The inclusion materials addressed above are of the highest cubic symmetry.
Here we also consider the lowest triclinic symmetry to investigate the
influence of single-crystal symmetry. The triclinic material is copper sulfate
pentahydrate (CSP) and it has an anisotropy index of 2.37, which is close to
that of Inconel. Comparing its results in Fig. 8 (a4) with those of Inconel in
Fig. 8 (a3), we can see that the two materials have roughly the same level of
theoretical-FE agreement, and their ‘starting points’ are similar. It
illustrates that the anisotropy factor has a larger effect on the
backscattering amplitude compared with the symmetry of the material.

An extra observation from Figs. 8(a1)-(a4) is that the cycle of the
backscattering amplitude curve only varies subtly among the four materials, as
further illustrated in Fig. 8(b). This supports our earlier finding of the
cyclic behavior being mainly determined by the size and shape (rather than the
material properties) of the inclusion.

All above results show that the theoretical model can be used to evaluate the
backscattering of Rayleigh waves and to establish a direct relation between
Rayleigh backscattering and the material and geometrical properties of the
inclusion.

##  5 Conclusion

In this work, we developed a 2D theoretical model for Rayleigh-to-Rayleigh
backscattering by an inclusion. The model is formulated in the frequency
domain based on the reciprocity theorem using the far-field Green’s function,
and the Born approximation is invoked to derive the final result. The model is
widely valid for a surface or subsurface inclusion with a regular or irregular
geometry, and prominently, the inclusion can be isotropic or anisotropic. A FE
model is established to provide relatively accurate reference data for
evaluating the approximations of the theoretical model. The comparison of the
theoretical and FE results across a range of scattering problems led to
various conclusions, mainly including:

1\. The theory exhibits very good agreement with the FE results for isotropic
inclusions with differently-defined material properties ($\Delta\rho$ or/and
$\Delta c_{pjkl}$) and different shapes (half-circle, ellipse, square and
rectangle, and even irregular) at small $k_{R}D$. As a result of the use of
the Born approximation, the theory starts to break down gradually with the
increase of $k_{R}D$ due to an increasing phase shift caused by the inclusion.
In addition, the sign of the phase shift determines how the theoretical curve
scales (compresses or expands over $k_{R}D$) with respect to the FE results.

2\. The theoretical model agrees well with the FE results for the
backscattering by a subsurface inclusion when the depth of the inclusion is
smaller than the wavelength, and it gradually loses its accuracy as the depth
exceeds a wavelength.

3\. The theory can predict the RMS backscattering amplitude of Rayleigh waves
by an anisotropic inclusion (studied cases include three cubic and one
triclinic materials). The results revealed the good applicability of the
theory to inclusions of weak anisotropy, and also uncovered the larger effect
of the anisotropy, rather than the symmetry, of the inclusion on
backscattering.

4\. The backscattering curve demonstrates an oscillating, cyclic behavior. The
cycle period is mainly determined by the geometry of the inclusion, and it
becomes smaller with the increase of the lateral-to-depth dimension ratio
$L/D$ of the inclusion. The material properties of the inclusion and their
differences to the host material properties only have a very small influence
on the cycle period.

Generally speaking, we have demonstrated the applicability of our theoretical
model to a wide spectrum of inclusion types, especially at low frequencies.
The model would have great potential for further applications in the area of
NDE as well as, for example, for developing models to describe Rayleigh wave
scattering in polycrystals.

## Acknowledgments

This work is supported by the China Scholarship Council and National Natural
Science Foundation of China (Grant No. 92060111). BL gratefully acknowledges
the Imperial College Research Fellowship and MH the generous funding from the
NDE group at Imperial.

## Appendices

##  Appendix A Detailed backscattering equations

The expression of $A^{\text{sc}}\left(\omega\right)$ in Eq. 17 can be further
evaluated, leading to

$\begin{split}A^{\text{sc}}\left(\omega\right)&=\sum_{p,j,k,l}^{1,3}A_{pjkl}^{\text{sc}}\left(\omega\right)=\sum_{l}^{1,3}Aa_{l}^{\text{sc}}-\sum_{p,j,k,l}^{1,3}Ab_{pjkl}^{\text{sc}}\left(\omega\right)\textrm{,}\end{split}$
(A.1)

with

$\displaystyle\begin{split}Aa^{\text{sc}}_{1}\left(\omega\right)=-A_{0}\int_{V}\Delta\rho\omega^{2}U_{R}^{2}\left(\mathbf{x}_{\mathrm{s}}\right)\exp(2\text{i}k_{R}x_{\mathrm{s}})\mathrm{\leavevmode\nobreak\
d}V\textrm{,}\end{split}$ (A.2a)
$\displaystyle\begin{split}Aa^{\text{sc}}_{3}\left(\omega\right)=\leavevmode\nobreak\
A_{0}\int_{V}\Delta\rho\omega^{2}W_{R}^{2}\left(\mathbf{x}_{\mathrm{s}}\right)\exp(2\text{i}k_{R}x_{\mathrm{s}})\mathrm{\leavevmode\nobreak\
d}V\textrm{,}\end{split}$ (A.2b)
$\displaystyle\begin{split}Ab^{\text{sc}}_{1111}\left(\omega\right)=\leavevmode\nobreak\
A_{0}\Delta
c_{1111}\left(\mathbf{x}_{\mathrm{s}}\right)\int_{V}\left[k_{R}^{2}U_{R}(z_{\mathrm{s}})^{2}\right]\exp\left(2\text{i}k_{R}x_{\mathrm{s}}\right)\mathrm{\leavevmode\nobreak\
d}V\textrm{,}\end{split}$ (A.2c)
$\displaystyle\begin{split}Ab^{\text{sc}}_{1113}\left(\omega\right)=\leavevmode\nobreak\
A_{0}\Delta
c_{1113}\left(\mathbf{x}_{\mathrm{s}}\right)\int_{V}\left[\text{i}k_{R}^{2}U_{R}(z_{\mathrm{s}})W_{R}(z_{\mathrm{s}}))\right]\exp\left(2\text{i}k_{R}x_{\mathrm{s}}\right)\mathrm{\leavevmode\nobreak\
d}V\textrm{,}\end{split}$ (A.2d)
$\displaystyle\begin{split}Ab^{\text{sc}}_{1131}\left(\omega\right)=\leavevmode\nobreak\
A_{0}\Delta
c_{1131}\left(\mathbf{x}_{\mathrm{s}}\right)\int_{V}\left[-\text{i}k_{R}U_{R}(z_{\mathrm{s}})\frac{\partial
U_{R}\left(z_{\mathrm{s}}\right)}{\partial
z_{\mathrm{s}}}\right]\exp\left(2\text{i}k_{R}x_{\mathrm{s}}\right)\mathrm{\leavevmode\nobreak\
d}V\textrm{,}\end{split}$ (A.2e)
$\displaystyle\begin{split}Ab^{\text{sc}}_{1133}\left(\omega\right)=\leavevmode\nobreak\
A_{0}\Delta
c_{1133}\left(\mathbf{x}_{\mathrm{s}}\right)\int_{V}\left[k_{R}U_{R}(z_{\mathrm{s}})\frac{\partial
W_{R}\left(z_{\mathrm{s}}\right)}{\partial
z_{\mathrm{s}}}\right]\exp\left(2\text{i}k_{R}x_{\mathrm{s}}\right)\mathrm{\leavevmode\nobreak\
d}V\textrm{,}\end{split}$ (A.2f)
$\displaystyle\begin{split}Ab^{\text{sc}}_{1313}\left(\omega\right)=\leavevmode\nobreak\
A_{0}\Delta
c_{1313}\left(\mathbf{x}_{\mathrm{s}}\right)\int_{V}\left[-k_{R}^{2}W_{R}(z_{\mathrm{s}})W_{R}(z_{\mathrm{s}})\right]\exp\left(2\text{i}k_{R}x_{\mathrm{s}}\right)\mathrm{\leavevmode\nobreak\
d}V\textrm{,}\end{split}$ (A.2g)
$\displaystyle\begin{split}Ab^{\text{sc}}_{1331}\left(\omega\right)=\leavevmode\nobreak\
A_{0}\Delta
c_{1331}\left(\mathbf{x}_{\mathrm{s}}\right)\int_{V}\left[k_{R}W_{R}(z_{\mathrm{s}})\frac{\partial
U_{R}\left(z_{\mathrm{s}}\right)}{\partial
z_{\mathrm{s}}}\right]\exp\left(2\text{i}k_{R}x_{\mathrm{s}}\right)\mathrm{\leavevmode\nobreak\
d}V\textrm{,}\end{split}$ (A.2h)
$\displaystyle\begin{split}Ab^{\text{sc}}_{1333}\left(\omega\right)=\leavevmode\nobreak\
A_{0}\Delta
c_{1333}\left(\mathbf{x}_{\mathrm{s}}\right)\int_{V}\left[\text{i}k_{R}W_{R}(z_{\mathrm{s}})\frac{\partial
W_{R}\left(z_{\mathrm{s}}\right)}{\partial
z_{\mathrm{s}}}\right]\exp\left(2\text{i}k_{R}x_{\mathrm{s}}\right)\mathrm{\leavevmode\nobreak\
d}V\textrm{,}\end{split}$ (A.2i)
$\displaystyle\begin{split}Ab^{\text{sc}}_{3131}\left(\omega\right)=\leavevmode\nobreak\
A_{0}\Delta
c_{1331}\left(\mathbf{x}_{\mathrm{s}}\right)\int_{V}\left[-\frac{\partial
U_{R}\left(z_{\mathrm{s}}\right)}{\partial z_{\mathrm{s}}}\frac{\partial
U_{R}\left(z_{\mathrm{s}}\right)}{\partial
z_{\mathrm{s}}}\right]\exp\left(2\text{i}k_{R}x_{\mathrm{s}}\right)\mathrm{\leavevmode\nobreak\
d}V\textrm{,}\end{split}$ (A.2j)
$\displaystyle\begin{split}Ab^{\text{sc}}_{3133}\left(\omega\right)=\leavevmode\nobreak\
A_{0}\Delta
c_{1333}\left(\mathbf{x}_{\mathrm{s}}\right)\int_{V}\left[-\text{i}\frac{\partial
U_{R}\left(z_{\mathrm{s}}\right)}{\partial z_{\mathrm{s}}}\frac{\partial
W_{R}\left(z_{\mathrm{s}}\right)}{\partial
z_{\mathrm{s}}}\right]\exp\left(2\text{i}k_{R}x_{\mathrm{s}}\right)\mathrm{\leavevmode\nobreak\
d}V\textrm{,}\end{split}$ (A.2k)
$\displaystyle\begin{split}Ab^{\text{sc}}_{3333}\left(\omega\right)=\leavevmode\nobreak\
A_{0}\Delta
c_{3333}\left(\mathbf{x}_{\mathrm{s}}\right)\int_{V}\left[\frac{\partial
W_{R}\left(z_{\mathrm{s}}\right)}{\partial z_{\mathrm{s}}}\frac{\partial
W_{R}\left(z_{\mathrm{s}}\right)}{\partial
z_{\mathrm{s}}}\right]\exp\left(2\text{i}k_{R}x_{\mathrm{s}}\right)\mathrm{\leavevmode\nobreak\
d}V\textrm{,}\end{split}$ (A.2l)
$\displaystyle\begin{split}Ab^{\text{sc}}_{1311}\left(\omega\right)=Ab^{\text{sc}}_{1113}\left(\omega\right),Ab^{\text{sc}}_{3111}\left(\omega\right)=Ab^{\text{sc}}_{1131}\left(\omega\right),Ab^{\text{sc}}_{3113}\left(\omega\right)=Ab^{\text{sc}}_{1331}\left(\omega\right),\\\
Ab^{\text{sc}}_{3311}\left(\omega\right)=Ab^{\text{sc}}_{1133}\left(\omega\right),Ab^{\text{sc}}_{3313}\left(\omega\right)=Ab^{\text{sc}}_{1333}\left(\omega\right),Ab^{\text{sc}}_{3331}\left(\omega\right)=Ab^{\text{sc}}_{3133}\left(\omega\right),\end{split}$
(A.2m)

The detailed expression of the RMS of the backscattering amplitude
$A^{\text{sc}}_{rms}\left(\omega\right)$ in Eq. 20 can be written as

$\begin{split}A^{\text{sc}}_{rms}\left(\omega\right)&=\sqrt{\left<A_{pjkl}^{\text{sc}}\left(\omega\right)A^{sc}_{\alpha\beta\gamma\delta}\left(\omega\right)\right>}=\sqrt{\left<\sum_{pjkl}^{1,3}Ab_{pjkl}^{\text{sc}}\left(\omega\right)\sum_{\alpha\beta\gamma\delta}^{1,3}Ab^{sc}_{\alpha\beta\gamma\delta}\left(\omega\right)\right>}\\\
&=\sqrt{\left<Ab_{1111}^{\text{sc}}Ab_{1111}^{\text{sc}}\right>+\left<Ab_{1111}^{\text{sc}}Ab_{1113}^{\text{sc}}\right>+\left<Ab_{1111}^{\text{sc}}Ab_{1131}^{\text{sc}}\right>+...+\left<Ab_{3333}^{\text{sc}}Ab_{3333}^{\text{sc}}\right>}\textrm{,}\end{split}$
(A.3)

where the individual terms in the square root are obtained by multiplying and
averaging the respective $A$ terms in Eq. A.2. Here, we give an example
expression for one of the terms by

$\allowdisplaybreaks\begin{split}\left<Ab_{1111}^{\text{sc}}Ab_{1111}^{\text{sc}}\right>=A_{0}^{2}&\left<\Delta
c_{1111}\Delta c_{1111}\right>\\\
&\times\int_{V}\int_{V}k_{R}^{4}U_{R}(z)^{2}U_{R}(z_{\mathrm{s}})^{2}\exp\left(2\text{i}k_{R}(x+x_{\mathrm{s}})\right)\mathrm{\leavevmode\nobreak\
d}V\mathrm{\leavevmode\nobreak\ d}V\textrm{.}\end{split}$ (A.4)

The elastic covariance $\left<\Delta c_{pjkl}\Delta
c_{\alpha\beta\gamma\delta}\right>$ involved in Eq. A.3 is given by

$\begin{split}&\left<\Delta c_{pjkl}\Delta
c_{\alpha\beta\gamma\delta}\right>=\left<c_{pjkl}c_{\alpha\beta\gamma\delta}\right>-\left<c_{pjkl}\right>\left<c_{\alpha\beta\gamma\delta}\right>\\\
&=\left(\left<a_{ja}a_{pb}a_{kc}a_{ld}a_{\alpha m}a_{\beta n}a_{\gamma
o}a_{\delta q}\right>-\left<a_{ja}a_{pb}a_{kc}a_{ld}\right>\left<a_{\alpha
m}a_{\beta n}a_{\gamma o}a_{\delta
q}\right>\right)c_{abcd}c_{mnoq}\end{split}$ (A.5)

for materials of any symmetry. $\textbf{a}(\theta,\phi,\xi)$ is the rotation
matrix defined using the Euler angles 46. The final expressions are provided
in the supplementary material.

## References

  * Ayter and Auld 1979 S. Ayter and B. Auld.  Characterization of surface wave scattering by surface breaking cracks.  _Proceedings of the ARPA/AFML Review of Progress in Quantitative NDE; Thompson, DO, Ed_ , pages 498–504, 1979.
  * Achenbach 2000 J. D. Achenbach.  Calculation of wave fields using elastodynamic reciprocity.  _International journal of solids and structures_ , 37(46-47):7043–7053, 2000.
  * Wang et al. 2019 C. Wang, O. Balogun, and J. D. Achenbach.  Scattering of a rayleigh wave by a near surface crack which is normal to the free surface.  _International Journal of Engineering Science_ , 145:103162, 2019.
  * Phan et al. 2013a H. Phan, Y. Cho, and J. D. Achenbach.  A theoretical study on scattering of surface waves by a cavity using the reciprocity theorem.  In _Nondestructive Testing of Materials and Structures_ , pages 739–744. Springer, 2013a.
  * Phan et al. 2013b H. Phan, Y. Cho, and J. D. Achenbach.  Validity of the reciprocity approach for determination of surface wave motion.  _Ultrasonics_ , 53(3):665–671, 2013b.
  * Yang and Achenbach 2017 C. Yang and J. D. Achenbach.  Time domain scattering of elastic waves by a cavity, represented by radiation from equivalent body forces.  _International Journal of Engineering Science_ , 115:43–50, June 2017.
  * wang et al. 2020 C. wang, O. Balogun, and J. D. Achenbach.  Application of the reciprocity theorem to scattering of surface waves by an inclined subsurface crack.  _International Journal of Solids and Structures_ , 207:82–88, 2020.
  * Auld 1990 B. Auld.  _Acoustic Fields and Waves in Solids. Vol. I.(Book)_.  Robert E. Krieger Publishing Company, 435, 1990.
  * Schmerr 2016 L. W. Schmerr.  _Fundamentals of ultrasonic nondestructive evaluation_.  Springer, 2016.
  * Gubernatis et al. 1977 J. Gubernatis, E. Domany, J. Krumhansl, and M. Huberman.  The born approximation in the theory of the scattering of elastic waves by flaws.  _Journal of Applied Physics_ , 48(7):2812–2819, 1977.
  * Kino 1978 G. Kino.  The application of reciprocity theory to scattering of acoustic waves by flaws.  _Journal of Applied Physics_ , 49(6):3190–3199, 1978.
  * Auld 1979 B. Auld.  General electromechanical reciprocity relations applied to the calculation of elastic wave scattering coefficients.  _Wave motion_ , 1(1):3–10, 1979.
  * Snieder 1986 R. Snieder.  3-d linearized scattering of surface waves and a formalism for surface wave holography.  _Geophysical Journal International_ , 84(3):581–605, 1986.
  * Razin 2010 A. Razin.  Scattering of a rayleigh surface acoustic wave by a small-size inhomogeneity in a solid half-space.  _Radiophysics and Quantum Electronics_ , 53(7):417–431, 2010.
  * Schmerr Jr and Sedov 2011 L. W. Schmerr Jr and A. Sedov.  Ultrasonic measurement models for contact testing with bulk waves, surface waves, and plate waves.  _Research in Nondestructive Evaluation_ , 22(3):129–146, 2011.
  * Xu et al. 2022 L. Xu, K. Wang, Y. Su, Y. He, J. Yang, S. Yuan, and Z. Su.  Surface/sub-surface crack-scattered nonlinear rayleigh waves: A full analytical solution based on elastodynamic reciprocity theorem.  _Ultrasonics_ , 118:106578, 2022.
  * Phan et al. 2018 H. Phan, Y. Cho, and W. Li.  A theoretical approach to multiple scattering of surface waves by shallow cavities in a half-space.  _Ultrasonics_ , 88:16–25, 2018.
  * Arias and Achenbach 2004 I. Arias and J. D. Achenbach.  Rayleigh wave correction for the bem analysis of two-dimensional elastodynamic problems in a half-space.  _International journal for numerical methods in engineering_ , 60(13):2131–2146, 2004.
  * Liu et al. 2011 W. Liu, Y. Cho, H. Phan, and J. D. Achenbach.  Study on the scattering of 2-d rayleigh waves by a cavity based on bem simulation.  _Journal of Mechanical Science and Technology_ , 25(3):797–802, 2011.
  * Huang et al. 2021 M. Huang, G. Sha, P. Huthwaite, S. Rokhlin, and M. Lowe.  Longitudinal wave attenuation in polycrystals with elongated grains: 3d numerical and analytical modeling.  _The Journal of the Acoustical Society of America_ , 149(4):2377–2394, 2021.
  * Liu et al. 2019 Y. Liu, A. Van Pamel, P. B. Nagy, and P. Cawley.  Investigation of ultrasonic backscatter using three-dimensional finite element simulations.  _The Journal of the Acoustical Society of America_ , 145(3):1584–1595, 2019.
  * Vu and Kinra 1985 B. Q. Vu and V. K. Kinra.  Diffraction of rayleigh waves in a half-space. i. normal edge crack.  _The Journal of the Acoustical Society of America_ , 77(4):1425–1430, 1985.
  * Jian et al. 2006 X. Jian, Y. Fan, R. Edwards, and S. Dixon.  Surface-breaking crack gauging with the use of laser-generated rayleigh waves.  _Journal of applied physics_ , 100(6):064907, 2006\.
  * Hassan and Veronesi 2003 W. Hassan and W. Veronesi.  Finite element analysis of rayleigh wave interaction with finite-size, surface-breaking cracks.  _Ultrasonics_ , 41(1):41–52, 2003.
  * Grabec et al. 2017 T. Grabec, M. Ryzy, and I. A. Veres.  Numerical modeling of surface elastic wave scattering in polycrystalline materials.  In _2017 IEEE International Ultrasonics Symposium (IUS)_ , pages 1–4. IEEE, 2017.
  * Grabec et al. 2022 T. Grabec, I. A. Veres, and M. Ryzy.  Surface acoustic wave attenuation in polycrystals: Numerical modeling using a statistical digital twin of an actual sample.  _Ultrasonics_ , 119:106585, 2022.
  * Rose 2000 J. L. Rose.  _Ultrasonic waves in solid media_.  Acoustical Society of America, 2000.
  * Phan et al. 2013c H. Phan, Y. Cho, and J. D. Achenbach.  Application of the reciprocity theorem to scattering of surface waves by a cavity.  _International Journal of Solids and Structures_ , 50(24):4080–4088, 2013c.
  * Aki and Richards 2002 K. Aki and P. G. Richards.  _Quantitative seismology_.  University Science, 2002.
  * Van Pamel et al. 2015 A. Van Pamel, C. R. Brett, P. Huthwaite, and M. J. Lowe.  Finite element modelling of elastic wave scattering within a polycrystalline material in two and three dimensions.  _The Journal of the Acoustical Society of America_ , 138(4):2326–2336, 2015.
  * Van Pamel et al. 2017 A. Van Pamel, G. Sha, S. I. Rokhlin, and M. J. Lowe.  Finite-element modelling of elastic wave propagation and scattering within heterogeneous media.  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 473(2197):20160738, 2017.
  * Van Pamel et al. 2018 A. Van Pamel, G. Sha, M. J. Lowe, and S. I. Rokhlin.  Numerical and analytic modelling of elastodynamic scattering within polycrystalline materials.  _The Journal of the Acoustical Society of America_ , 143(4):2394–2408, 2018.
  * Ryzy et al. 2018 M. Ryzy, T. Grabec, P. Sedlák, and I. A. Veres.  Influence of grain morphology on ultrasonic wave attenuation in polycrystalline media.  _The Journal of the Acoustical Society of America_ , 143(1):219–229, 2018.
  * Bai et al. 2018 X. Bai, B. Tie, J.-H. Schmitt, and D. Aubry.  Finite element modeling of grain size effects on the ultrasonic microstructural noise backscattering in polycrystalline materials.  _Ultrasonics_ , 87:182–202, 2018.
  * Huang et al. 2020 M. Huang, G. Sha, P. Huthwaite, S. Rokhlin, and M. Lowe.  Maximizing the accuracy of finite element simulation of elastic wave propagation in polycrystals.  _The Journal of the Acoustical Society of America_ , 148(4):1890–1910, 2020.
  * Sarris et al. 2021 G. Sarris, S. G. Haslinger, P. Huthwaite, P. B. Nagy, and M. J. Lowe.  Attenuation of rayleigh waves due to surface roughness.  _The Journal of the Acoustical Society of America_ , 149(6):4298–4308, 2021.
  * Sha and Rokhlin 2018 G. Sha and S. Rokhlin.  Universal scaling of transverse wave attenuation in polycrystals.  _Ultrasonics_ , 88:84–96, 2018.
  * Huang et al. 2022 M. Huang, S. I. Rokhlin, and M. J. S. Lowe.  Appraising scattering theories for polycrystals of any symmetry using finite elements.  _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 380(2231):20210382, sep 2022.
  * Rajagopal et al. 2012 P. Rajagopal, M. Drozdz, E. A. Skelton, M. J. Lowe, and R. V. Craster.  On the use of absorbing layers to simulate the propagation of elastic waves in unbounded isotropic media using commercially available finite element packages.  _NDT & E INT_, 51:30–40, 2012.
  * Huthwaite 2014 P. Huthwaite.  Accelerated finite element elastodynamic simulations using the gpu.  _Journal of Computational Physics_ , 257:687–707, 2014\.
  * Gnedin et al. 2018 N. Y. Gnedin, V. A. Semenov, and A. V. Kravtsov.  Enforcing the Courant–Friedrichs–Lewy condition in explicitly conservative local time stepping schemes.  _Journal of Computational Physics_ , 359:93–105, 2018.
  * Kak and Slaney 2001 A. C. Kak and M. Slaney.  _Principles of computerized tomographic imaging_.  SIAM, 2001.
  * Huthwaite 2012 P. E. Huthwaite.  _Quantitative imaging with mechanical waves_.  PhD thesis, Imperial College London, 2012.
  * Viktrov 1967 I. Viktrov.  _Rayleigh and Lamb waves: physical theory and applications_.  Plenumn, 1967.
  * Spiegel et al. 2012 M. R. Spiegel, J. J. Schiller, and R. A. Srinivasan.  _Probability and statistics_.  Schaum’s outline series. McGraw-Hill, New York, fourth edi edition, 2012\.  ISBN 9780071795579.
  * Weaver 1990 R. L. Weaver.  Diffusivity of ultrasound in polycrystals.  _Journal of the Mechanics and Physics of Solids_ , 38(1):55–86, 1990.


# On the Sensitivity of Deep Load Disaggregation to Adversarial Attacks

Hafsa Bousbiat  University of Klagenfurt
Klagenfurt, Austria
<EMAIL_ADDRESS>Yassine Himeur  University of Dubai
Dubai, UAE
<EMAIL_ADDRESS>Abbes Amira  University of Sharjah
Sharjah, UAE
<EMAIL_ADDRESS>Wathiq Mansoor  University of Dubai
Dubai, UAE
<EMAIL_ADDRESS>
###### Abstract

Non-intrusive Load Monitoring (NILM) algorithms, commonly referred to as load
disaggregation algorithms, are fundamental approaches for effective energy
management. Despite the success of deep models in load disaggregation, they
face various challenges, particularly those pertaining to privacy and
security. This paper investigates the sensitivity of prominent deep NILM
baselines to adversarial attacks, which have proven to be a significant threat
in domains such as computer vision and speech recognition. Adversarial attacks
entail the introduction of imperceptible noise into the input data with the
aim of misleading the neural network into generating erroneous outputs. We
investigate the Fast Gradient Sign Method (FGSM), a well-known adversarial
attack, to perturb the input sequences fed into two commonly employed CNN-
based NILM baselines: the Sequence-to-Sequence (S2S) and Sequence-to-Point
(S2P) models. Our findings provide compelling evidence for the vulnerability
of these models, particularly the S2P model which exhibits an average decline
of 20% in the F1-score even with small amounts of noise. Such weakness has the
potential to generate profound implications for energy management systems in
residential and industrial sectors reliant on NILM models.

###### Index Terms:

Non-intrusive Load Monitoring, Load Disaggregation, Adversarial Attacks, Smart
Home Energy Management

##  I INTRODUCTION

Providing appliance-level feedback has been widely acknowledged as a promising
approach for achieving significant energy savings [1]. The previous strategy
can lead up to 12% of energy savings in the residential sector [2, 3]. It
allows individuals to make informed decisions about their daily interaction
with in-home electrical devices. In this regard, Non-intrusive Load Monitoring
(NILM) has emerged as a promising technology that utilizes advanced algorithms
to infer the power consumption of individual appliances based on aggregate
household measurements [4, 5]. With minimal hardware requirements, NILM
approaches are an attractive solution for energy management and sustainability
efforts. More precisely, NILM algorithms rely only on a single metering point
reporting the total power consumption of the household to identify operating
appliances [5] allowing thus to avoid complex hardware deployment and
maintance.

Exceptional performance enhancement in NILM scholarship can be observed in
recent years. Particularly, the adoption of deep models for load
disaggregation by Kelly et al. [6] in 2015 has been a turning point [4, 2].
Compared to statistical classical models, deep neural networks have shown
tremendous potential in enhancing the accuracy of energy consumption estimates
[7]. A main advantage of these models is the automatic feature extraction that
provides an end-to-end solution. Nonetheless, the robustness of these models
against adversarial attacks has become a growing concern in recent years [8].
Adversarial attacks can pose a significant threat to NILM’s accuracy and
reliability by injecting small amounts of noise into the input data to deceive
the neural network and generate erroneous outputs [9]. In the case of NILM,
adversarial attacks can cause wrong energy consumption estimates [10], privacy
breaches, instability in home energy management systems [11, 12], and even
safety risks [13] when used in critical infrastructures. This paper aims to
address this gap by evaluating the sensitivity of two popular NILM baselines
to adversarial attacks.

The remainder of this paper is organized as follows: Section II discusses
related work and recent literature providing insights into the state-of-the-
art NILM models and adversarial attacks. Section III describes the different
steps followed to perform the attack, including the generation of adversarial
examples and their evaluation on the target models. Section IV describes the
experimental setup followed by a presentation the results of the experimental
setup in Section V. Finally, Section VI concludes the paper by summarizing the
main contributions of this research and discussing potential future directions
in this area.

##  II RELATED WORK

Figure 1: Flow chart of the performed attack

Non-intrusive load monitoring (NILM) gained significant attention in the past
five years, with neural networks being widely acknowledged as state-of-the-art
models [14, 4]. Popular deep baselines for NILM, such as Sequence-to-Sequence
and Sequence-to-point, are based on convolutional layers [6]. However, other
contributions based on LSTM [15] have also demonstrated the robustness of
these layers in identifying appliance usage. Advanced models, such as the BERT
model [16] and UNET model [17], have also been suggested in the literature.

Despite the success of these models in estimating the power consumption of
household appliances under different conditions [7], they suffer from privacy
and security issues [18]. They can reveal sensitive information about
individuals’ daily routines, as illustrated in [19, 18, 20], leading to
concerns from the consumer’s side. The security of NILM systems is thus a main
concern when deep models are adopted. Yet, it received little attention. A
compact review of adversarial attacks in smart grid scenarios was suggested in
[21], considering deep neural networks. The review study revealed that despite
the potential for extreme harm, these attacks have received very little
attention from the research community. Two types of attacks can be performed:
white-box and black-box attacks.

White-box attacks assume full knowledge of the attacked model [22], and are
gaining interest in related work. They were particularly considered in [9],
targeting IoT devices where the case of residential smart meters was
considered as a case study. The obtained results demonstrated that adversarial
samples are often indistinguishable from real samples, leading to high success
rates of the performed attack. This finding illustrates that adversarial
attacks represent a major threat for energy management systems (EMS) both in
residential and industrial setups, notably in the case of sensitive sectors,
as they can easily be used to fool demand response programs. However, this
attack model makes the assumption of complete knowledge of the model, which is
hard to obtain in real scenarios. In contrast, black-box attacks assume no
knowledge of the model or the training data [22]. The main advantage of
adopting this threat model, compared to the previous one, is its practicality.
The significant influence of this attack on energy analytics and load
disaggregation was explored in [23, 24]. The findings revealed that black-box
attacks can be performed even in cases where the adversary has limited
knowledge of the target system.

##  III METHODOLOGY

The threat model aims to fool neural networks into producing inaccurate
estimates of power consumption, thereby disrupting the energy management
systems of smart buildings. As the attack involves adding perturbations to the
input data, the attacker’s primary challenge would be to determine the optimal
amount of noise to introduce, maximizing the success rate while remaining
undetected by data poisoning detection methods. The current paper’s
methodological design aims to evaluate this aspect while also assessing
whether the attack’s effect is similar for two different disaggregation
models, given appliances with heterogeneous power consumption magnitudes.

The steps involved in the methodology are depicted in Figure 1. Firstly, the
models are trained on the original data and evaluated on a testing data set,
without attack. Secondly, the gradient of the pre-trained model is used to
generate perturbed inputs, which are then fed back into the network to produce
the corresponding predictions.

###  III-A Models

The Seq2Point and Seq2Seq [6] models for NILM have gained significant
attention due to their competitive performance in estimating the power
consumption of household appliances. These models are based on convolutional
layers and have similar structures, differing only in the shape of their
output. The choice of the Seq2Point and Seq2Seq models as representative
models for load disaggregation is a solid choice, given their competitive
performance, wide acknowledgement, and adoption in the research community. To
implement these models, we leveraged the code available from NILMtk [25], a
popular open-source toolkit for NILM research. The training process was
conducted for 150 epochs, using an Adam optimizer with a learning rate of
$10^{-4}$.

###  III-B Attack

We propose to assess the robustness of the baseline models against a widely
recognized white-box adversarial attack, namely the Fast Sign Gradient Sign
Method (FGSM). This attack leverages the gradient of the model to generate an
adversarial sample of the input. For time series data, this technique involves
computing the gradient of the loss function $\mathcal{L}$ with respect to the
input sequence, in order to create a new sequence that increases the loss.
This newly generated sequence is referred to as an adversarial sample
$\tilde{x}$, where $\epsilon$ is a multiplier to ensure that the added
perturbations remain small.

$\tilde{x}=x+\epsilon*sign(\nabla_{x}\mathcal{L}(\theta,x,y))$ (1)

The FGSM has an important feature, namely that the gradient calculation is
performed with respect to each input sequence, enabling the attacker to
identify how each value of the input sequence contributes to the loss
function, with the ultimate goal of misleading the model. It is worth noting
that the FGSM aims to mislead an already pre-trained model, without affecting
its parameters. Hence, it represents a powerful tool to assess the model’s
robustness against adversarial attacks, as it provides insight into the
model’s vulnerabilities without affecting its training process.

##  IV Evaluation

###  IV-A Data

The present study utilizes data from UKDALE, a dataset recorded in the UK,
comprising five buildings, with the first building having a recording period
longer than three years. The pre-trained models are constructed using data
from the first building, where a four-month recording period with a sampling
rate of 30 seconds is considered. The study further considers four appliances
with varying consumption patterns: the washing machine, the kettle, the
microwave, and the fridge. For testing purposes, a period of one month is
considered for each appliance.

###  IV-B Metrics

For a fair evaluation of the models as well as the effect of the FGSM attack,
three metrics are used, the Mean Average Error (MAE), the F1-score and the
Normalised Disaggregation Error (NDE). The MAE for an appliance $i$ is
described as follows:

$\text{MAE}^{(i)}=\frac{1}{N}\cdot\sum_{t=0}^{N-1}|\hat{y}_{t}^{(i)}-y_{t}^{(i)}|$
(2)

where, $y_{t}$ is the ground truth values of the power consumption,
$\hat{y}_{t}$ is the predicted power consumption, and $N$ represents the
number of samples. In addition to the MAE, the NDE is considered as
recommended in recent literature [26]:

$\text{NDE}^{(i)}=\sqrt{\frac{\sum_{t=0}^{N-1}{(\hat{y}_{t}^{(i)}-y_{t}^{(i)}})^{2}}{\sum_{t=0}^{N-1}{(y_{t}^{(i)}})^{2}}}$
(3)

Furthermore, the F1-score is used to assess the performance of the model on
state estimation of different appliances, defined as follows:

$F1-score=\frac{2\cdot Precision\cdot Recall}{Precision+Recall}$ (4)

Where the Precision = TP/(TP+FP), Recall = TP/(TP+FN). A threshold of 10
_watts_ was used to derive states of the fridge and a threshold of 500 _watts_
was used in the case of the remaining appliances.

Figure 2: Seq2Point model for the washing machine Figure 3: Seq2Seq model for
the washing machine

##  V RESULTS & DISCUSSION

The results obtained for both baseline models are presented in Table I. We
report on four scenarios: the performance of the model without attack and the
performance of the model under attack with three different values of
$\epsilon$ (0.01, 0.10, and 0.25). For each of the considered scenarios, we
report on the three metrics for the four appliances considered.

TABLE I: Obtained results for Seq2Point models |  | Without Attack |  | FGSM ($\epsilon=0.01$) |  | FGSM ($\epsilon=0.10$) |  | FGSM ($\epsilon=0.25$)
---|---|---|---|---|---|---|---|---
|  | MAE | F1-score | NDE |  | MAE | F1-score | NDE |  | MAE | F1-score | NDE |  | MAE | F1-score | NDE
Seq2Point
Washing Machine |  | 12.4 | 0.67 | 0.73 |  | 14.20 | 0.63 | 0.75 |  | 126.63 | 0.39 | 1.17 |  | 153.94 | 0.13 | 1.40
Kettle |  | 05.4 | 0.86 | 0.38 |  | 5.5 | 0.86 | 0.39 |  | 9.48 | 0.77 | 0.52 |  | 15.38 | 0.66 | 0.71
Microwave |  | 06.0 | 0.63 | 0.71 |  | 7.57 | 0.60 | 0.79 |  | 11.68 | 0.34 | 1.02 |  | 12.41 | 0.22 | 1.04
Fridge |  | 15.0 | 0.83 | 0.42 |  | 26.4 | 0.74 | 0.61 |  | 39.3 | 0.61 | 0.79 |  | 43.0 | 0.45 | 0.93
Seq2Seq
Washing Machine |  | 12.9 | 0.66 | 0.76 |  | 13.24 | 0.66 | 0.75 |  | 51.55 | 0.66 | 0.75 |  | 71.39 | 0.59 | 0.82
Kettle |  | 06.3 | 0.89 | 0.35 |  | 6.43 | 0.89 | 0.35 |  | 7.24 | 0.88 | 0.36 |  | 9.17 | 0.83 | 0.40
Microwave |  | 07.4 | 0.60 | 0.70 |  | 7.83 | 0.63 | 0.69 |  | 9.98 | 0.65 | 0.75 |  | 12.48 | 0.63 | 0.85
Fridge |  | 17.6 | 0.79 | 0.43 |  | 17.33 | 0.77 | 0.43 |  | 29.59 | 0.62 | 0.55 |  | 36.22 | 0.60 | 0.63

When no attack is performed, both Seq2Point and Seq2Seq models yield
equivalent performance considering the three metrics with the Seq2Point model
providing slightly better results considering all appliances and all metrics.
The latter observation is reflected by a minimal f1-score of 67% and a maximum
of 12.4 _watts_ for the MAE. The previous observation aligns well with
existing literature [7] suggesting the superiority of Seq2Point model.

The results show that as the epsilon value increases, the performance of all
appliances deteriorates with the Kettle performing the best and the Fridge
performing the worst, for both models. At an epsilon value of 0.01, all
appliances show a slight decrease in performance. Nonetheless, starting from
an epsilon value of 0.10, the MAE of all appliances increases significantly,
with the Washing Machine appliance showing the highest increase. Meanwhile,
F1-score decreases for all appliances, with the Microwave appliance showing
the highest decrease. At an epsilon value of 0.25, the MAE, F1-score, and NDE
values continue to deteriorate for all appliances, with the Washing Machine
appliance showing the highest MAE and the Fridge appliance showing the highest
NDE.

When comparing the effect of the added noise on the Seq2Seq and Seq2Point
models, the Seq2Seq shows more robustness against the attack. On one side,
Seq2Point demonstrates higher sensitivity to the added noise with a decrease
in all metrics by a factor of three for frequently used appliances (i.e., the
fridge, the microwave, the kettle) and a decrease by factor of ten (10) in the
case of the washing machine for $\epsilon=0.10$, considered as one of the
major loads. On the other side, the Seq2Seq demonstrates moderate sensitivity
where a decrease with only a factor of 1.5 is recorded for frequently used
appliances and a factor of 5 in the case of the washing machine for the same
values of $\epsilon$. It can be seen that even if high values of $\epsilon$
lead to slight deterioration in the MAE, the seq2seq model still yields
approximately the same F1-score values unlike the Seq2Point model where the
deterioration is recorded in all metrics.

The Figure 2 illustrates one generated activation for the washing machine for
the Seq2Point model. It can be observed that the noise added first reflects on
the OFF states through small fluctuations in the order of hundreds of _watts_
starting from values of $\epsilon$=0.10. Furthermore, the figure shows that
higher values of $\epsilon$ lead to an attenuated signal or undetected states
(e.g., second peak in the figure). Figure 3 illustrates the same activation
generated with Seq2Seq model. It is clear that the model shows more robustness
where it is capable of generating stable signal even with higher values of
added noise with slight deterioration in the estimated values.

The adversarial attack illustrated in this study on NILM models can easily be
leveraged in smart homes to mislead the neural networks. This is particularly
relevant in the case of home automation functions [27]. For example, an attack
can be performed on light control systems that are based on NILM and leads to
considerable losses if unnoticed. It can also lead to the incorrect estimation
of load flexibility in the system. If a device is identified as ”off” although
it is actually ”on,” it will not be considered for automated control, which
will affect the ability to shave demand peaks when a large number of flexible
loads are attacked. For instance, if the system estimates that 75 kW of demand
during peak hours in a day stem from the heating system in winter time and the
NILM algorithm misclassifies the state of this appliance, then the peak
reduction potential is lost. Moreover, the attack can be leveraged at specific
times to mislead occupancy detection systems and easily lead to erroneous
predictions that can be used by malicious parties.

##  VI CONCLUSION

This paper suggested to assess the robustness of the Seq2Seq and Seq2Point
models against adversarial attacks considering the FGSM attack. Even with
smaller values of $\epsilon$, the obtained results reveal high sensitivity of
the Seq2Point model, while demonstrating a moderate sensitivity of the Seq2Seq
model. This study confirms findings from related about the superiority of
Seq2Point in normal scenarios but highlights that Seq2Seq reveal more strength
to injected noise. Adversarial attacks can have destructive effect on EMS
systems and different energy services. In future work, we aim reinforce the
robustness of these models to such attacks using different strategies
including augmenting the set of training data with adversarial samples.

## References

  * [1] Y. Himeur, A. Alsalemi, F. Bensaali, A. Amira, I. Varlamis, G. Bravos, C. Sardianos, and G. Dimitrakopoulos, “Techno-economic assessment of building energy efficiency systems using behavioral change: A case study of an edge-based micro-moments solution,” _Journal of Cleaner Production_ , vol. 331, p. 129786, 2022.
  * [2] Y. Himeur, A. Alsalemi, F. Bensaali, A. Amira, and A. Al-Kababji, “Recent trends of smart nonintrusive load monitoring in buildings: A review, open challenges, and future directions,” _International Journal of Intelligent Systems_ , vol. 37, no. 10, pp. 7124–7179, 2022.
  * [3] R. Gopinath, M. Kumar, C. P. C. Joshua, and K. Srinivas, “Energy management using non-intrusive load monitoring techniques–state-of-the-art and future research directions,” _Sustainable Cities and Society_ , vol. 62, p. 102411, 2020.
  * [4] H. Bousbiat, Y. Himeur, I. Varlamis, F. Bensaali, and A. Amira, “Neural load disaggregation: Meta-analysis, federated learning and beyond,” _Energies_ , vol. 16, no. 2, p. 991, 2023.
  * [5] H. Bousbiat, C. Klemenjak, Y. Himeur, W. Elmenreich, A. Amira, W. Mansoor, and S. Atalla, “Neural nilm learning paradigms: From centralised to decentralised learning,” in _2022 5th International Conference on Signal Processing and Information Security (ICSPIS)_. IEEE, 2022, pp. 138–142.
  * [6] J. Kelly and W. Knottenbelt, “Neural nilm: Deep neural networks applied to energy disaggregation,” in _Proceedings of the 2nd ACM international conference on embedded systems for energy-efficient built environments_ , 2015, pp. 55–64.
  * [7] N. Batra, R. Kukunuri, A. Pandey, R. Malakar, R. Kumar, O. Krystalakos, M. Zhong, P. Meira, and O. Parson, “Towards reproducible state-of-the-art energy disaggregation,” in _Proceedings of the 6th ACM international conference on systems for energy-efficient buildings, cities, and transportation_ , 2019, pp. 193–202.
  * [8] H. Tan, L. Wang, H. Zhang, J. Zhang, M. Shafiq, and Z. Gu, “Adversarial attack and defense strategies of speaker recognition systems: A survey,” _Electronics_ , vol. 11, no. 14, p. 2183, 2022.
  * [9] A. Singh and B. Sikdar, “Adversarial attack and defence strategies for deep-learning-based iot device classification techniques,” _IEEE Internet of Things Journal_ , vol. 9, no. 4, pp. 2602–2613, 2021.
  * [10] L. Ou, Z. Qin, S. Liao, T. Li, and D. Zhang, “Singular spectrum analysis for local differential privacy of classifications in the smart grid,” _IEEE Internet of Things Journal_ , vol. 7, no. 6, pp. 5246–5255, 2020.
  * [11] J. Kotak and Y. Elovici, “Adversarial attacks against iot identification systems,” _IEEE Internet of Things Journal_ , 2022.
  * [12] M. Kaselimi, E. Protopapadakis, A. Voulodimos, N. Doulamis, and A. Doulamis, “Towards trustworthy energy disaggregation: A review of challenges, methods, and perspectives for non-intrusive load monitoring,” _Sensors_ , vol. 22, no. 15, p. 5872, 2022.
  * [13] R. Liaqat and I. A. Sajjad, “An event matching energy disaggregation algorithm using smart meter data,” _Electronics_ , vol. 11, no. 21, p. 3596, 2022.
  * [14] P. Huber, A. Calatroni, A. Rumsch, and A. Paice, “Review on deep neural networks applied to low-frequency nilm,” _Energies_ , vol. 14, no. 9, p. 2390, 2021.
  * [15] M. Kaselimi, N. Doulamis, A. Voulodimos, E. Protopapadakis, and A. Doulamis, “Context aware energy disaggregation using adaptive bidirectional lstm models,” _IEEE Transactions on Smart Grid_ , vol. 11, no. 4, pp. 3054–3067, 2020.
  * [16] Z. Yue, C. R. Witzig, D. Jorde, and H.-A. Jacobsen, “Bert4nilm: A bidirectional transformer model for non-intrusive load monitoring,” in _Proceedings of the 5th International Workshop on Non-Intrusive Load Monitoring_ , 2020, pp. 89–93.
  * [17] A. Faustine, L. Pereira, H. Bousbiat, and S. Kulkarni, “Unet-nilm: A deep neural network for multi-tasks appliances state detection and power estimation in nilm,” in _Proceedings of the 5th International Workshop on Non-Intrusive Load Monitoring_ , 2020, pp. 84–88.
  * [18] H. Bousbiat, G. Leitner, and W. Elmenreich, “Ageing safely in the digital era: A new unobtrusive activity monitoring framework leveraging on daily interactions with hand-operated appliances,” _Sensors_ , vol. 22, no. 4, p. 1322, 2022.
  * [19] P. A. Schirmer, I. Mporas, and A. Sheikh-Akbari, “Identification of tv channel watching from smart meter data using energy disaggregation,” _Energies_ , vol. 14, no. 9, p. 2485, 2021.
  * [20] H. Bousbiat, C. Klemenjak, G. Leitner, and W. Elmenreich, “Augmenting an assisted living lab with non-intrusive load monitoring,” in _2020 IEEE international instrumentation and measurement technology conference (I2MTC)_. IEEE, 2020, pp. 1–5.
  * [21] J. Hao and Y. Tao, “Adversarial attacks on deep learning models in smart grids,” _Energy Reports_ , vol. 8, pp. 123–129, 2022.
  * [22] N. Akhtar and A. Mian, “Threat of adversarial attacks on deep learning in computer vision: A survey,” _Ieee Access_ , vol. 6, pp. 14 410–14 430, 2018.
  * [23] J. Wang and P. Srikantha, “Stealthy black-box attacks on deep learning non-intrusive load monitoring models,” _IEEE Transactions on Smart Grid_ , vol. 12, no. 4, pp. 3479–3492, 2021.
  * [24] J. Wang, “Deep learning on smart meter data: Non-intrusive load monitoring and stealthy black-box attacks,” Ph.D. dissertation, The University of Western Ontario (Canada), 2020.
  * [25] H. Bousbiat, A. Faustine, C. Klemenjak, L. Pereira, and W. Elmenreich, “Unlocking the full potential of neural nilm: On automation, hyperparameters & modular pipelines,” _IEEE Transactions on Industrial Informatics_ , 2022\.
  * [26] C. Klemenjak, S. Makonin, and W. Elmenreich, “Towards comparability in non-intrusive load monitoring: On data and performance evaluation,” in _2020 IEEE power & energy society innovative smart grid technologies conference (ISGT)_. IEEE, 2020, pp. 1–5.
  * [27] H. Bousbiat, C. Klemenjak, G. Leitner, and W. Elmenreich, “Augmenting an assisted living lab with non-intrusive load monitoring,” in _2020 IEEE International Instrumentation and Measurement Technology Conference (I2MTC)_ , 2020, pp. 1–5.