# Cutoff phenomenon and entropic uncertainty for random quantum circuits

Sangchul Oh, Sabre Kais† Department of Chemistry, Department of Physics and
Astronomy, and Purdue Quantum Science and Engineering Institute, Purdue
University, West Lafayette, IN, USA  $^†$ Corresponding author:
<EMAIL_ADDRESS>
###### Abstract

How fast a state of a system converges to a stationary state is one of the
fundamental questions in science. Some Markov chains and random walks on
finite groups are known to exhibit the non-asymptotic convergence to a
stationary distribution, called the cutoff phenomenon. Here, we examine how
quickly a random quantum circuit could transform a quantum state to a Haar-
measure random quantum state. We find that random quantum states, as
stationary states of random walks on a unitary group, are invariant under the
quantum Fourier transform. Thus the entropic uncertainty of random quantum
states has balanced Shannon entropies for the computational bases and the
quantum Fourier transform bases. By calculating the Shannon entropy for random
quantum states and the Wasserstein distances for the eigenvalues of random
quantum circuits, we show that the cutoff phenomenon occurs for the random
quantum circuit. It is also demonstrated that the Dyson-Brownian motion for
the eigenvalues of a random unitary matrix as a continuous random walk
exhibits the cutoff phenomenon. The results here imply that random quantum
states could be generated with shallow random circuits.

††: Electronic Structures: Focus Issue on “Quantum Chemistry in the Ages …”

  *

Keywords: Random circuits, quantum computing, cutoff phenomenon, random walks

##  1 Introduction

How many shuffles are enough to ensure that a deck of 52 cards is mixed
randomly? The answer is $\tfrac{3}{2}\log_{2}52$, approximately 8.55, for the
riffle shuffling case, as shown by Aldous and Diaconis [1] and by Bayer and
Diaconis [2]. Fewer than this number are not enough to mix the deck of cards
and more do not significantly improve the randomness. This non-asymptotic
convergence to a steady or equilibrium state is called the cutoff phenomenon
[3] and has been discovered in various fields. Some finite Markov chains
exhibit the cutoff phenomenon. These include the Ehrenfest urn model as a
simplified diffusion model [4], random walks on a hypercube [5], some
Metropolis algorithms [3], Glauber dynamics of Ising models [6], and certain
quantum Markov chains [7].

The card shuffling problem is considered a random walk on the symmetric group
$S_{52}$. The cutoff phenomenon occurs also for random walks on a compact Lie
group. Rosenthal [8] considered that random walk on the special orthogonal
group ${\rm SO}(N)$ given by repeated rotations by a fixed angle through
randomly chosen planes. Rosenthal showed that the random walk on ${\rm
SO}(N)$, after $1/(2(1-\cos\theta))N\log N+cN$ rotations with a fixed angle
$\theta$ and a constant $c$, converges rapidly to the Haar measure in total
variation distance. Porod [9, 10] showed that random walks on ${\rm O}(N)$,
${\rm U}(N)$, and ${\rm Sp}(N)$ given by random reflections converge to Haar
measure in the total variation distances and the cutoff phenomena occurs at
$\frac{1}{2}N\log N$ steps.

Random circuit sampling is a task to sample bit-strings from a probability
distribution defined by random quantum circuits, and considered a good
candidate to demonstrate quantum advantage with noisy-intermediate scale
quantum computers. In 2019, it was implemented on the Sycamore processor with
53 qubits [11] and recently on the Zuchongzhi processor with 56 and 60 qubits
[12, 13]. In both quantum computations, the random circuit was implemented by
applying repeatedly, up to 20 and 24 cycles, randomly-chosen single qubit
gates and the two-qubit gate acting on the nearest neighbor qubits. The number
of cycles plays the same role as the number of times a deck of cards is
shuffled or the number of random rotations on Lie groups. So it is natural to
ask the same questions which have been answered for the card shuffling
problem. How many depths of random quantum circuits, i.e., the number of
cycles, are needed to obtain a Haar-measure random unitary operator or to
transform an initial quantum state into a pure random quantum state? Does the
cutoff phenomenon occur? What is going on in a system after a steady state has
been reached? Is there any tool to measure the closeness to a steady state in
addition to the total variance distance? In this paper, we present partial
answers to these questions using the random circuit implemented on the
Sycamore processor and a time-dependent random Hamiltonian model. The former
is considered a discrete random walk on a unitary group and the latter a
continuous random walk called the Dyson-Brownian motion for eigenvalues of a
unitary operator. A pure random quantum state at the steady state will be
analyzed with the Shannon entropy. Instead of the total variation distance
$||\mu_{k}-\mu_{\rm Haar}||_{\rm TV}$ between the probability distribution
$\mu_{k}$ at the $k$-th time step and the Haar measure distribution $\mu_{\rm
Haar}$, we employ the Wasserstein distance for the distribution of eigenvalues
of a unitary operator.

The paper is organized as follows. In Sec. 2, we will discuss the cutoff
phenomenon for random quantum circuits by calculating the Shannon entropy and
the Wasserstein distance. We show that the Shannon entropy of a random quantum
state is invariant under the quantum Fourier transform. We discuss the
entropic uncertainty relation of random quantum states for the computational
bases and the quantum Fourier transform bases. In Sec. 3, we investigate the
cutoff phenomenon of the Dyson-Brownian motion generated by a time-dependent
random Hamiltonian. Finally, in Sec. 4, we will summarize the result and
present the discussion.

##  2 Cutoff Phenomenon for Random Quantum Circuits

Let us begin with the introduction to the random quantum sampling implemented
on the Sycamore quantum processor and on the Zuchongzhi quantum processor [11,
12]. The task is to sample bit-strings $x=a_{1}a_{2}\cdots
a_{n}\in\\{0,1\\}^{n}$ from the probability
$p(x)=|\left\langle{x}\right|U\left|{0}\right\rangle^{\otimes n}|^{2}$ given
by a random quantum circuit $U$ acting on $n$ qubits where
$\left|{0}\right\rangle^{\otimes n}=\left|{0_{1}\cdots 0_{n}}\right\rangle$ is
the initial state and $\left|{x}\right\rangle\equiv\left|{a_{1}a_{2}\cdots
a_{n}}\right\rangle$ is a computational basis. The random quantum circuit $U$
implemented on both the Sycamore and Zuchongzhi processors is given by
repeatedly applying random unitary operators $U_{k}$ and finally the single-
qubit gates $S$ before the measurement

$U=SU_{m}\cdots U_{2}U_{1}\,,$ (1)

where each random quantum circuit $U_{k}$ is composed of single-qubit gates
$S$ chosen randomly from the set $\\{\sqrt{X},\sqrt{Y},\sqrt{W}\\}$ on all
qubits and two-qubit gates on the pair of qubits selected in the sequence of
the coupler activation patterns of a 2-dimensional array of qubits. Millions
of bit-strings were sampled from the Sycamore and Zuchongzhi processors. The
distributions of bit strings obtained from these noisy quantum processors are
deviated from the ideal distribution. Oh and Kais [14, 15, 16] investigated
this deviation using the random matrix theory and the Wasserstein distances.

Eq. (1) is considered random walks or random rotations on a unitary group. If
the number of cycles $k$ is large, $U^{(k)}\equiv U_{k}\cdots U_{1}$ would
approach a random unitary operator $U_{\rm Haar}$ sampled from the Haar
probability measure on a unitary group ${\rm U}(2^{n})$. Typically, the
convergence to a stationary state is measured by the total variation distance
$||v^{*(k)}-\mu_{\rm Haar}||_{\rm TV}$ where $v^{*(k)}$ is the distribution of
the random walk at step $k$ and $\mu_{\rm Haar}$ is the Haar measure [8, 9,
10]. For the Sycamore processor, the sub-linear convergence, the depth
proportional to $\sqrt{n}$, was claimed by calculating the average entropy of
random quantum states [17]. Emerson et al. [18] studied a pseudo-random
circuit given by repeated applications of $n$ single-qubit random gates
sampled from the Haar measure on ${\rm U}(2)$ and simultaneous two-qubit
interactions. Moreover, Emerson et al. employed the measure of entanglement
for a multipartite system as the measure of the convergence. It was shown that
this random quantum circuit converges to the Haar measure if the circuit depth
$m$ is larger than $m_{c}\equiv{\cal O}(n^{3}N^{2})$ with $N=2^{n}$. This is
larger than the cutoff step $k_{c}\equiv\frac{1}{2}N\ln N={\cal O}(nN)$ for
random rotations or random reflections on Lie groups shown by Rosenthal [8]
and Porod [9, 10].

Different measures have been used to quantify the convergence to the Haar
measure distribution, and give rise to the different cutoff steps. Here, we
employ the Shannon entropy for quantum states and the Wasserstein distance
between the eigenvalue distribution of a random unitary operator sampled from
the Haar measure and those of random quantum circuits. Random unitary matrices
drawn from the Haar measure on a unitary group are called the circular unitary
ensemble. The properties of random quantum states and the eigenvalue
distribution of the circular unitary ensemble are well known from the random
matrix theory [19].

Let us first investigate how close a quantum state at the $k$-th step,
$\left|{\psi^{(k)}}\right\rangle=R^{(k)}\left|{0}\right\rangle=U_{k}\cdots
U_{1}\left|{0}\right\rangle$, is to a random quantum state, a stationary state
of random walks on a unitary group. An immediate question is what a pure
random quantum state is and how to generate it. A random pure state can be
generated in several ways, that is, there are several ways of drawing a random
unitary operator from the Haar measure [20, 21, 22]. A random unitary operator
could be sampled from the Haar measure through the Euler angle method or the
QR decomposition of a complex Gaussian random matrix. Basically, a random
quantum state
$\left|{\psi}\right\rangle=\sum_{i=0}^{N-1}a_{i}\left|{i}\right\rangle$ can be
viewed as a random vector on the $(2N-1)$ sphere and expansion coefficients
$a_{i}=x_{i}+iy_{i}$ are drawn from the normal distributions, i.e.,
$x_{i},y_{i}\sim{\cal N}(0,1)$. The distribution
$P=\\{p_{1},p_{2},\dots,p_{N-1}\\}$ of probabilities $p_{i}=|a_{i}|^{2}$ of a
random quantum state obeys the $\chi^{2}$ distribution with 2 degree of
freedom [19]

${\rm Pr}(p)=(N-1)(1-p)^{N-2}\,.$ (2)

The $\chi^{2}$ distribution for $P$ of $p_{i}$ of random quantum states makes
it possible to calculate the average Shannon entropy. The Shannon entropy for
the probability distribution $P=\\{p_{1},p_{2},\dots,p_{N-1}\\}$ is

$\displaystyle H(P)=H(\left|{\psi}\right\rangle)=-\sum_{i=0}^{N-1}p_{i}\ln
p_{i}\,,$ (3)

where $\sum_{i}p_{i}=1$ and $p_{i}=|a_{i}|^{2}$. The Shannon entropy
$H(\left|{\psi}\right\rangle)$ for a quantum state measures the amount of
uncertainty or the concentration of $P=(p_{0},\dots,p_{N-1})$. The average
Shannon entropy over random quantum states can be calculated with Eqs. (2) and
(3) and is given by

$\langle H(\left|{\psi}\right\rangle)\rangle=\ln N-1+\gamma\,,$ (4)

where $\gamma\approx 0.5772$ is the Euler constant and $\langle\cdots\rangle$
is the average over random quantum states. So the Shannon entropy of a quantum
state could be used as a measure of convergence of random walks on a unitary
group.

Figure 1: The distributions of $p_{i}$ for (a) a uniform superposition of all
basis states $N=20$, (c) a random quantum state generated by a random unitary
matrix ${\rm U}(N)$ with $N=20$, and (e) a random quantum state of $n=14$
qubits generated by a random quantum circuit implemented on the Sycamore
processor [11, 23]. Plots (b), (d), and (f) are the distributions of $p_{k}$
after the quantum Fourier transform of (a), (b), and (c), respectively. The
red lines indicate $p=1/N$ and the Shannon entropy $H$ for each state is
labeled by $H$.

Note that the Shannon entropy $H(\left|{\psi}\right\rangle)$ is defined with
respect to a specific basis set $\\{\left|{i}\right\rangle\\}$. If the same
quantum state is expanded in terms of another basis set
$\\{\left|{k}\right\rangle\\}$,
$\left|{\psi}\right\rangle=\sum_{k=0}^{N-1}b_{k}\left|{k}\right\rangle$, its
Shannon entropy $H(\left|{\psi}\right\rangle)$ will change. For two
orthonormal basis sets, $\\{\left|{i}\right\rangle\\}$ and
$\\{\left|{k}\right\rangle\\}$, the entropic uncertainty relation [24, 25, 26,
27] is written as

$H(P)+H(Q)\geq-2\ln c$ (5)

where $c=\max|\langle i|k\rangle|$, $P=(p_{0},\dots,p_{N-1})$ with
$p_{i}=|a_{i}|^{2}$, and $Q=(q_{0},\dots,q_{N-1})$ with $q_{k}=|b_{k}|^{2}$.
We consider the computational basis and the quantum Fourier transformed basis,
which are mutually unbiased. The quantum Fourier transform (QFT) is the
discrete Fourier transform of the amplitude vector
$(x_{0},x_{1},\dots,x_{N-1})$ to another amplitude vector
$(y_{0},y_{1},\dots,y_{N-1})$, defined by

$y_{k}=\frac{1}{\sqrt{N}}\sum_{j=0}^{N-1}e^{i2\pi\frac{jk}{N}}\,x_{j}\,.$ (6)

The QFT acting on $\left|{\psi_{P}}\right\rangle$ is written as

$\left|{\psi_{P}}\right\rangle=\sum_{i=0}^{N-1}x_{i}\left|{i}\right\rangle\quad\longrightarrow\quad\left|{\psi_{Q}}\right\rangle=\sum_{k=0}^{N-1}y_{k}\left|{k}\right\rangle\,.$
(7)

Fig. 1 illustrates the distributions of probabilities $p_{i}$ for two kinds of
quantum states in the computational basis and probabilities $q_{k}$ in the QFT
basis. First, consider a uniform superposition of all possible basis states in
the computational basis that is given by
$\left|{\psi}\right\rangle=\frac{1}{\sqrt{N}}(\left|{0}\right\rangle+\left|{1}\right\rangle+\cdots+\left|{N-1}\right\rangle)$.
As shown in Fig. 1 (a), its $p_{i}$ is distributed uniformly, $p_{i}={1}/{N}$
with $i=0,1,\dots,N-1$ and its Shannon entropy has the maximum value,
$H(P)=\ln(N)$. The equally-likely distribution among $N$ possible states
implies the largest randomness, and the greatest entropy [28]. However,
$\left|{\psi}\right\rangle$ in the QFT basis set is localized at one site, so
its entropy $H(Q)$ is zero, as depicted in Fig. 1 (d). The entropic
uncertainty is given by $H(P)+H(Q)=\ln(N)$. Next, consider a random quantum
state $\left|{\phi}\right\rangle=U_{\rm Haar}\left|{0}\right\rangle$. Fig. 1
(b) plots the distribution of $p_{i}$ of a random quantum state where $U_{\rm
Haar}$ is generated by the QR decomposition of an $N\times N$ complex Gaussian
random matrix. The entropy of a random quantum state is approximately given by
$H(P)=\ln N-1+\gamma\approx 2.302$ for $N=10$. Fig. 1 (e) plots the
distribution of $q_{k}$ in the QFT basis. Interestingly, the entropy $H(P)$ of
a random quantum state in the computational basis is almost equal to its
entropy $H(Q)$ in its QFT basis. We observe that the average entropy for all
random quantum states in the computational basis, generated by a random
circuit $U$, is equal to that in the QFT basis. A random quantum state may
have the balanced entropic uncertainty, $H(P)=H(Q)=\ln N-1+\gamma$. It is
analogous to a coherent state in the sense that the latter has balanced or
symmetric minimum uncertainty: $\Delta x=\Delta p$. It may be interesting to
understand why a random quantum state is invariant under the QFT. Fig. 1 (c)
depicts the distribution of $p_{i}$ for a random quantum state generated by a
random quantum circuit implemented on the Sycamore processor for $n=12$ and
the cycles $m=14$ [23]. Fig. 1 (f) plots the distribution of $q_{k}$ in the
QFT basis. One can see that the Shannon entropy of a random quantum state in
the computational basis is almost same as that of its QFT state. Note that the
random circuit here is implemented on a classical computer without any noise.
The Shannon entropy in Fig. 1 e is close to the theoretical value, $H=\ln
N-1+\gamma\approx 7.8949$. The Shannon entropy calculated from the Sycamore
data for $n=12$ and $m=14$ is $H\approx 8.217$ and close to $\ln(4096)\approx
8.31$ [29].

Figure 2: For a random quantum circuit for $n=14$ qubits and up to $m=14$
cycles [11], (a) Shannon entropy of the average of quantum states and (b)
Wasserstein distance between the eigenvalues of random circuits and the Haar
random unitary operator are plotted as a function of the number of gates
applied, i.e., the depth of a quantum circuit.

Since a random quantum state is characterized by its Shannon entropy,
$H(P)\approx\ln N-(1-\gamma)$, the Shannon entropy
$H(\left|{\psi^{(k)}}\right\rangle)$ of a quantum state
$\left|{\psi^{(k)}}\right\rangle=U_{k}\dots U_{1}\left|{0}\right\rangle$ at
the $k$-th step of random walks could be used as a convergence measure to see
whether the cutoff phenomenon happens. The sub-linear convergence proportional
to $\sqrt{n}$ was claimed by calculating the average entropy of quantum states
[17]. As shown in Fig. 2 (a), the Shannon entropy of a quantum state converges
to $\ln N-(1-\gamma)$ as the number of gates increases and remains there. An
interesting question is what happens to a quantum state after the Shannon
entropy converges.

The eigenvalue distribution for random unitary operators drawn from the Haar
measure is well known [21], so the distance between the eigenvalue
distributions could be used to measure the convergence of a random walk. We
consider the Wasserstein distance defined by

$W_{p}(P,Q)=\left(\inf_{J\in J(P,Q)}\int||x-y||^{p}\,dJ(x,y)\right)^{1/p}\,,$
(8)

where ${\cal J}(P,Q)$ denotes all joint distributions $J(P,Q)$. If $P$ and $Q$
are the empirical distributions of a data set $\\{x_{1},\dots,x_{n}\\}$ and
$\\{y_{1},\dots,y_{n}\\}$ respectively, then the Wasserstein distance is given
by the distance between order statistics

$W_{p}(P,Q)=\left(\sum_{i=1}^{n}||x_{(i)}-y_{(i)}||^{p}\right)^{1/p}\,,$ (9)

where $x_{(i)}$ is the $i$-th order statistic of a samples, i.e., its $i$-th
smallest value. Fig. 2 (b) plots the Wasserstein distance of order 1 between
the eigenvalue distribution of the circular unitary ensemble and that of a
random quantum circuit as a function of the number of quantum gates applied.
The calculation of the Wasserstein distance supports the cutoff phenomenon for
random quantum circuits as the Shannon entropy does.

##  3 Dyson-Brownian Motions on a Unitary Group

In Sec. 2, the random walk on the unitary group is implemented by applying the
sequence of random quantum circuits, $\\{U_{1},U_{2},\dots\\}$. A random
quantum state after the $k$-th step is $\left|{\psi}\right\rangle=U_{k}\dots
U_{1}\left|{0}\right\rangle$. This may be considered as a discrete process. In
this section, we consider a continuous random walk given by a time-dependent
random Hamiltonian to see how quickly a quantum state converges to a random
quantum state. The time evolution operator at $t+\delta t$ is given by

$\displaystyle U(t+\delta t)$ $\displaystyle=e^{i\frac{\delta
t}{\hbar}H(t)}\,U(t)$ (10) $\displaystyle\approx\Bigl{[}1+i\frac{\delta
t}{\hbar}H(t)-\frac{1}{2}\left(\frac{\delta
t}{\hbar}\right)^{2}H^{2}(t)\Bigr{]}U(t)\,,$ (11)

where $U(0)=I$. The time dependent random Hamiltonian $H(t)$ at time $t$ is
obtained as follows. We draw a complex random matrix $A$ whose real and
imaginary parts of a matrix element $A_{ij}$ are sampled independently from
the normal distribution ${\cal N}(0,\sigma^{2}/2N)$ with the variance
$\sigma^{2}/2N$. The Hermitian property of the Hamiltonian $H$ is fulfilled by
summing $A$ and its conjugate transpose $A^{\dagger}$,
$H=\frac{1}{2}(A+A^{\dagger})$. $H$ is the elements of the Gaussian orthogonal
ensemble.

Figure 3: For the Dyson-Brownian random walk, (a) trajectories of eigenvalues
of a random unitary operator and (b) the Shannon entropy
$H(\left|{\psi}\right\rangle)$ of a quantum state are plotted. Here, we take
$20\times 20$ random Hamiltonian matrices, i.e., $N=20$ and the time step
$\delta t=0.01$. In (b) the red dotted horizontal line represents the Shannon
entropy of a random quantum state, $\ln N-1+\gamma\approx 2.5572$ for $N=20$.

The trajectories of eigenvalues of a random unitary operator $U$ are known as
the Dyson-Brownian motions and do not overlap with others. Eigenvalues repel
each other. We simulate the time evolution with a time-dependent random
Hamiltonian for $N=20$ and $\sigma=1$ so $\sigma^{2}/2N=1/2N$. For simplicity,
we take $\hbar=1$. We set the time step $\delta t=0.01$. The eigenvalues of
$U(t)$ from Eq. (11) are obtained by diagonalizing it and their trajectories
as a function of time are plotted in Fig. 3 (a). The time-evolution of the
Shannon entropy of a quantum state
$\left|{\psi(t)}\right\rangle=U(t)\left|{0}\right\rangle$ is shown in Fig. 3
(b). We observe the non-asymptotic convergence to $\ln N-1+\gamma\approx
2.5572$ for $N=20$. It is interesting to see the eigenvalues and the Shannon
entropy fluctuate after arriving at the steady state. Particularly, the
fluctuation in the Shannon entropy in Fig. 3 (b) is in contrast to no
fluctuation in the Shannon entropy for a finite random walk of a random
quantum circuit, shown in Fig. 1 (a).

##  4 Summary

In this paper, we examined some properties of random quantum states generated
by discrete and continuous random walks on a unitary group. It is found that
the Shannon entropy of a random quantum state generated by random quantum
circuits is invariant under the QFT in the sense that the Shannon entropy does
not change before and after applying the QFT. The entropic uncertainty
relation of a random quantum state for the computational bases and the QFT
bases is balanced, i.e., $H(P)=H(Q)$. This may remind us of the coherent state
with the balanced minimum uncertainty relation, $\delta x=\delta p$. We showed
that the cutoff phenomenon for a random quantum circuit occurs by calculating
the Shannon entropy and the Wasserstein distance for the eigenvalue
distributions. It is a open question whether the cutoff of random walks on a
unitary group scales with the number of qubits $n$ as sub-linear ${\cal
O}(\sqrt{n})$ or ${\cal O}({n}2^{n})$ while the numerical calculations here
seem to support the sub-linear scaling of the cutoff.

In addition to the demonstration of quantum advantage, random quantum circuits
may be applicable to solving interesting problems, for example, in randomized
linear algebra. The set of random quantum states
$\left|{\psi_{i}}\right\rangle=U_{\rm Haar}\left|{i}\right\rangle$ with
$i=0,1,\dots,2^{n}-1$ form the orthonormal basis set. The trace of a matrix
$A$ could be calculated using ${\rm
Tr}(A)\approx\frac{1}{M}\sum_{i=1}^{M}\left\langle{\psi_{i}}\right|A\left|{\psi_{i}}\right\rangle$
where $\left|{\psi_{i}}\right\rangle$ is a random quantum state [30].

## Acknowledgment

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, National Quantum Information Science Research Centers. We
also acknowledge the National Science Foundation under award number 1955907.

## References

## References

  * [1] Aldous D and Diaconis P 1986 The American Mathematical Monthly 93 333–348 ISSN 00029890, 19300972 URL http://www.jstor.org/stable/2323590
  * [2] Bayer D and Diaconis P 1992 The Annals of Applied Probability 2 294 – 313 URL https://doi.org/10.1214/aoap/1177005705
  * [3] Diaconis P 1996 Proceedings of the National Academy of Sciences 93 1659–1664 URL https://www.pnas.org/doi/pdf/10.1073/pnas.93.4.1659
  * [4] Ehrenfest P and Ehrenfest-Afanassjewa T 1907 Über zwei bekannte Einwände gegen das Boltzmannsche H-Theorem (Hirzel)
  * [5] Diaconis P, Graham R L and Morrison J A 1990 Random Structures & Algorithms 1 51–72 (Preprint https://onlinelibrary.wiley.com/doi/pdf/10.1002/rsa.3240010105) URL https://onlinelibrary.wiley.com/doi/abs/10.1002/rsa.3240010105
  * [6] Lubetzky E and Sly A 2013 Inventiones mathematicae 191 719–755 ISSN 1432-1297 URL https://doi.org/10.1007/s00222-012-0404-5
  * [7] Kastoryano M J, Reeb D and Wolf M M 2012 Journal of Physics A: Mathematical and Theoretical 45 075307 URL https://dx.doi.org/10.1088/1751-8113/45/7/075307
  * [8] Rosenthal J S 1994 The Annals of Probability 22 398 – 423 URL https://doi.org/10.1214/aop/1176988864
  * [9] Porod U 1996 The Annals of Probability 24 74–96 ISSN 00911798 URL http://www.jstor.org/stable/2244833
  * [10] Porod U 1996 Probability Theory and Related Fields 104 181–209 ISSN 1432-2064 URL https://doi.org/10.1007/BF01247837
  * [11] Arute F, Arya K, Babbush R, Bacon D, Bardin J C, Barends R, Biswas R, Boixo S, Brandao F G S L, Buell D A, Burkett B, Chen Y, Chen Z, Chiaro B, Collins R, Courtney W, Dunsworth A, Farhi E, Foxen B, Fowler A, Gidney C, Giustina M, Graff R, Guerin K, Habegger S, Harrigan M P, Hartmann M J, Ho A, Hoffmann M, Huang T, Humble T S, Isakov S V, Jeffrey E, Jiang Z, Kafri D, Kechedzhi K, Kelly J, Klimov P V, Knysh S, Korotkov A, Kostritsa F, Landhuis D, Lindmark M, Lucero E, Lyakh D, Mandrà S, McClean J R, McEwen M, Megrant A, Mi X, Michielsen K, Mohseni M, Mutus J, Naaman O, Neeley M, Neill C, Niu M Y, Ostby E, Petukhov A, Platt J C, Quintana C, Rieffel E G, Roushan P, Rubin N C, Sank D, Satzinger K J, Smelyanskiy V, Sung K J, Trevithick M D, Vainsencher A, Villalonga B, White T, Yao Z J, Yeh P, Zalcman A, Neven H and Martinis J M 2019 Nature 574 505–510 ISSN 0028-0836, 1476-4687 URL http://www.nature.com/articles/s41586-019-1666-5
  * [12] Wu Y, Bao W S, Cao S, Chen F, Chen M C, Chen X, Chung T H, Deng H, Du Y, Fan D, Gong M, Guo C, Guo C, Guo S, Han L, Hong L, Huang H L, Huo Y H, Li L, Li N, Li S, Li Y, Liang F, Lin C, Lin J, Qian H, Qiao D, Rong H, Su H, Sun L, Wang L, Wang S, Wu D, Xu Y, Yan K, Yang W, Yang Y, Ye Y, Yin J, Ying C, Yu J, Zha C, Zhang C, Zhang H, Zhang K, Zhang Y, Zhao H, Zhao Y, Zhou L, Zhu Q, Lu C Y, Peng C Z, Zhu X and Pan J W 2021 Phys. Rev. Lett. 127(18) 180501 URL https://link.aps.org/doi/10.1103/PhysRevLett.127.180501
  * [13] Zhu Q, Cao S, Chen F, Chen M C, Chen X, Chung T H, Deng H, Du Y, Fan D, Gong M, Guo C, Guo C, Guo S, Han L, Hong L, Huang H L, Huo Y H, Li L, Li N, Li S, Li Y, Liang F, Lin C, Lin J, Qian H, Qiao D, Rong H, Su H, Sun L, Wang L, Wang S, Wu D, Wu Y, Xu Y, Yan K, Yang W, Yang Y, Ye Y, Yin J, Ying C, Yu J, Zha C, Zhang C, Zhang H, Zhang K, Zhang Y, Zhao H, Zhao Y, Zhou L, Lu C Y, Peng C Z, Zhu X and Pan J W 2022 Science Bulletin 67 240–245 ISSN 2095-9273 URL https://www.sciencedirect.com/science/article/pii/S2095927321006733
  * [14] Oh S and Kais S 2022 The Journal of Physical Chemistry Letters 13 7469–7475 pMID: 35939529 URL https://doi.org/10.1021/acs.jpclett.2c02045
  * [15] Oh S and Kais S 2022 Phys. Rev. A 106(3) 032433 URL https://link.aps.org/doi/10.1103/PhysRevA.106.032433
  * [16] Oh S and Kais S 2023 Phys. Rev. A 107(2) 022610 URL https://link.aps.org/doi/10.1103/PhysRevA.107.022610
  * [17] Boixo S, Isakov S V, Smelyanskiy V N, Babbush R, Ding N, Jiang Z, Bremner M J, Martinis J M and Neven H 2018 Nature Physics 14 595–600 ISSN 1745-2481 URL https://doi.org/10.1038/s41567-018-0124-x
  * [18] Emerson J, Weinstein Y S, Saraceno M, Lloyd S and Cory D G 2003 Science 302 2098–2100 ISSN 0036-8075 URL https://science.sciencemag.org/content/302/5653/2098
  * [19] Haake F, Gnutzmann S and Kuś M 2010 Quantum Signatures of Chaos (Springer-Verlag Berlin Heidelberg) ISBN 978-3-642-05428-0
  * [20] Mezzadri F 2007 Notices of the American Mathematical Society 54 592 – 604 ISSN 0002-9920
  * [21] Meckes E S 2019 The Random Matrix Theory of the Classical Compact Groups Cambridge Tracts in Mathematics (Cambridge University Press)
  * [22] Ozols M 2009 unpublished essay on http://home. lu. lv/sd20008
  * [23] Martinis et al J M 2022 Quantum supremacy using a programmable superconducting processor, Dryad, Dataset https://doi.org/10.5061/dryad.k6t1rj8
  * [24] Deutsch D 1983 Phys. Rev. Lett. 50(9) 631–633 URL https://link.aps.org/doi/10.1103/PhysRevLett.50.631
  * [25] Kraus K 1987 Phys. Rev. D 35(10) 3070–3075 URL https://link.aps.org/doi/10.1103/PhysRevD.35.3070
  * [26] Maassen H and Uffink J B M 1988 Phys. Rev. Lett. 60(12) 1103–1106 URL https://link.aps.org/doi/10.1103/PhysRevLett.60.1103
  * [27] Coles P J, Berta M, Tomamichel M and Wehner S 2017 Rev. Mod. Phys. 89(1) 015002 URL https://link.aps.org/doi/10.1103/RevModPhys.89.015002
  * [28] Ambegaokar V and Clerk A A 1999 American Journal of Physics 67 1068–1073 ISSN 0002-9505 (Preprint https://pubs.aip.org/aapt/ajp/article-pdf/67/12/1068/10115996/1068_1_online.pdf) URL https://doi.org/10.1119/1.19084
  * [29] Oh S and Kais S (unpublished)
  * [30] Oh S and Kais S (in preparation) Estimating trace of a matrix with random quantum states


# Wave based damage detection in solid structures using artificial neural
networks

Frank Wuttke1,∗, Hao Lyu1, Amir S. Sattari1, Zarghaam H. Rizvi1
1 Geomechanics and Geotechnics Group, Kiel University, Germany
∗ Corresponding Author: Frank Wuttke, Email<EMAIL_ADDRESS>

###### Abstract

The identification of structural damages takes a more and more important role
within the modern economy, where often the monitoring of an infrastructure is
the last approach to keep it under public use. Conventional monitoring methods
require specialized engineers and are mainly time consuming. This research
paper considers the ability of neural networks to recognize the initial or
alteration of structural properties based on the training processes. The
presented work here is based on Convolutional Neural Networks (CNN) for wave
field pattern recognition, or more specifically the wave field change
recognition. The CNN model is used to identify the change within propagating
wave fields after a crack initiation within the structure. The paper describes
the implemented method and the required training procedure to get a successful
crack detection accuracy, where the training data are based on the dynamic
lattice model. Although the training of the model is still time consuming, the
proposed new method has an enormous potential to become a new crack detection
or structural health monitoring approach within the conventional monitoring
methods.

Keywords— Damage Detection, Structural Health Monitoring, Artificial Neural
Networks, Convolutional Neural Networks, Lattice Element Method

##  1 Introduction

For the permanent use of existing structures in urban areas, as well as for
lifelines and for safety structures such as dams, the successful monitoring of
structures is of the highest priority. Usually, conventional methods of
structural dynamics are used to analyse the state of structures and to find
existing or propagating damages, while wave based method are less used in
ordinary structural dynamics. These methods are more usual in the field of
non-destructive testing (NDT) (Kaewunruen and Remennikov, 2006; Farhangdoust
and Mehrabi, 2019) under use of very high excitation frequencies and shorter
wave lengths. Beside the type of analysis method - structural dynamics with
long wave length or ultra-sound methods in NDT with extreme short wave length
- the matter of the structural analysis is based on the active analysis of
excited vibrations or wave fields until the structural damage is detected.
Without any knowledge of preexisting damage zones, the analysis can take lots
of time. The present development shows a new strategy from a numerical case
study to analyse structures under use of artificial neural networks. The
artificial neural networks (ANN) (Sha and Edwards, 2007) are used to learn the
change of structural response patterns of damaged structures in opposition to
non-damaged structures.

With the booming development in large-scale data generation and computation
power, deep learning algorithms, especially deep convolutional networks (known
as CNNs or ConvNets), are developing by leaps and bounds. Deep learning
methods have been applied to various many fields, such as computer vision and
natural language processing, and often outperforms conventional methods (LeCun
et al., 2015). The multi-layer structured deep models can learn patterns from
data at multiple levels of abstraction. Convolution operation performs an
important role in CNN layers. ConvNets are particularly suitable to learn from
array-like data, such as audio, image, and video. These data are uniformly
sampled from signals in spatial and/or temporal domains. In many applications,
ConvNets have been used as the backbone for pattern extraction, e.g.,
recognising objects from images (Girshick, 2015), understanding text (Kim,
2014), and synthesising audio (van den Oord et al., 2016). In recent years,
attempts have been made to apply ConvNets to damage detection as well
(Abdeljaber et al., 2017; Rautela and Gopalakrishnan, 2019).

Deep learning methods have long been expected to be promising for wave-based
damage detection (Avci et al., 2021). Compared to the hand-engineered-feature-
based methods, the deep-learning-based method uses deep neural networks as a
feature extractor to learn representations from wave fields(Guo et al., 2020).
1D CNNs and RNNs are two popular structures to recognize patterns from 1D
signals. Abdeljaber et al. trained multiple 1D CNNs to detect whether damage
exists at specific locations (joints) (Abdeljaber et al., 2017). Their model
uses the acceleration signal at each joint as input and requires extra
workload to segment the signal into frames. Considering damage detection as a
binary prediction problem, i.e., predicting whether a crack exists from input
data, 1D CNN, RNN, and LSTM models all can achieve high accuracy (Rautela and
Gopalakrishnan, 2019). In one following paper, the authors developed a two-
stage damage detection method (Rautela and Gopalakrishnan, 2020). The method
determines whether a sample is damaged or not at the first stage and then
predict the location and length of the damage with another regressor network.
However, the regressor network deals only with the damage that is orthogonal
to the sample’s surface. Khan et al. transformed the structural vibration
responses in the time domain into two-dimensional spectral frame
representation and then applied a 2D CNN to distinguish between the undamaged
and various damaged states (Khan et al., 2019). Besides using CNNs for wave-
based damage detection, there are also methods using different input data or a
combination with other learning schemes. For example, in (Gulgec et al.,
2019), the authors proposed to use 2D CNN to predict the bounding box of a
crack from raw strain field. By adding noise and changing loading patterns to
augment the data, the model could achieve some robustness. Nunes and her
colleagues developed a hybrid strategy to detect undefined structural changes
(Nunes et al., 2020). They proposed to apply a unsupervised k-means clustering
method to the CNN learned features, and thus the features that are extracted
from the samples with undefined changes are expected to fall out of these
clusters. Our work differs from the above mentioned works by training an end-
to-end model to predict both crack shapes and locations from large amounts of
simulated wave fields.

The numerical treatments for the paper are done by use of a meso-scale method
as it is better suited to capture the effects as initial cracking and crack
propagation depending on the material parameter and the initial- & boundary
conditions without pre-definition of damaged patches, and it is also
applicable for 2D and 3D problems. The Lattice Element Method is a class of
discrete models in which the structural solid is represented as a 3D assembly
of one-dimensional elements (Wong et al., 2014; Rizvi et al., 2019; Sattari et
al., 2017). This idea allows one to provide robust models for propagation of
discontinuities, multiple crack interactions, or cracks coalescence even under
dynamic loads and wave fields. Different computational procedures for lattice
element methods for representing linear elastic continuum have been developed.
Beside different mechanical, hydro-mechanical and multi-physical developments,
the extension and basics for a new dynamic Lattice Element Method was
presented in (Rizvi et al., 2018, 2020). This development will be used in the
given paper for the health monitoring of structures.

To perform the damage detection as a suitable numerical software, pattern
indicators and specific designed neural networks are needed. The numerical
simulation is realized under use of Dynamic Lattice-Element Method, where the
advantage of the discontinuum method in opposition to continuum methods
related to the damage detection will be discussed in the methodology section.
The implemented artificial neural networks are also described in this section.
Based on the considered numerical and DNN (deep neural networks) models, a
case study of a 2D plane is performed to show the developments and results of
the new approach.

##  2 Methodology

###  2.1 Dynamic Lattice Approach

The assembly of the heterogeneous and homogeneous material will be generated
by specific meshing algorithms in LEM. The Lattice Element Models with the
lattice nodes can be considered as the centers of the unit cells, which are
connected by beams that can carry normal force, shear force and bending
moment. Because the strain energy stored in each element can be exceeded by a
given threshold, the element is either removed for cracking or assigned a
lower stiffness value. The method is based on minimizing the stored energy of
the system. The size of the localized fracture process zone around the static
or propagating crack plays a key role in failure mechanism, which is observed
in various models of linear elastic fracture mechanics and multi-scale
theories or homogenization techniques. Normally this propagating crack process
needs a regularization, however, an efficient way of dealing with this kind of
numerical problem is by introducing the embedded strong discontinuity into
lattice elements, resulting in mesh-independent computations of failure
response. The generation of the lattice elements are done by Voronoi cells and
Delaunay itself (Sattari et al. (2017); Moukarzel and Herrmann (1992)). With
the performance of this procedure an easy algebraic equation is generated for
the static case. To develop the dynamic LEM for simulation of a propagating
wave field, a more complex extension of the LEM is needed. The following
solution of the dynamic LEM is solved as a transient solution in the time
domain.

####  2.1.1 Equation of motion

To solve the dynamic LEM, the static LEM needs an extension of the equation of
motion. The general equation of motion without the damping term is defined by

$\mathbf{M}\ddot{u}+\mathbf{K}u=F(t)$ (1)

where $\mathbf{M}$ and $\mathbf{K}$ are the mass and the stiffness matrices
terms and $F(t)$ is the applied time-dependent force. Both matrices, the mass
and stiffness matrix, have to be defined in terms of the LEM definition.

####  2.1.2 Mass Matrix generation

The mass matrix or the consistent mass matrix (CMM) is generated either by
lumping the mass at the nodes or by following the variation mass lumping (VMM)
scheme. The VMM scheme is also implemented in the finite element method for
dynamic simulations. The element mass $M^{e}$ is computed using the following
equation

$M_{e}=\int\rho\left\lfloor N_{v}^{e}\right\rfloor^{T}N_{v}d\Omega$ (2)

If the shape functions are identical, that is, $N_{v}^{e}=N^{e}$, the mass
matrix is called the consistent mass matrix (CMM) or $M_{c}^{e}$.

$M_{c}^{e}=\int_{0}^{l}\rho A\left[N^{e}\right]^{T}N^{e}dx=\frac{1}{4}\rho
lA\int_{-1}^{1}\begin{bmatrix}1-\epsilon\\\
1+\epsilon\end{bmatrix}\begin{bmatrix}1-\epsilon&1+\epsilon\end{bmatrix}d\epsilon$
(3)

Where, $\rho$ is the density assigned to the Voronoi cells and $A$ and $l$ are
the area and the length of the lattice elements
The elemental mass matrix is symmetric, physically symmetric, and complies
with the condition of conservation and positively. To obtain the global mass
matrix, a congruent transformation is applied. In contrast to the stiffness
matrix, translational masses never vanish. All the translational masses are
retained in the local mass matrix. The global transformation is achieved
through the following equation.

$\bar{M}_{c}^{e}=\left[T^{e}\right]^{T}\left[M_{c}^{e}\right]\left[T^{e}\right]$
(4)

$M_{c}^{e}=\frac{1}{2}m^{e}\begin{bmatrix}1&0&0&0\\\ 0&1&0&0\\\ 0&0&1&0\\\
0&0&0&1\end{bmatrix}$ (5)

####  2.1.3 Element Stiffness Matrix

The force displacement component of a truss element is given by the spring
relation

$\left\\{F\right\\}=\left[K\right]\left\\{U\right\\}$ (6)

The vectors $\left\\{F\right\\}$ and $\left\\{U\right\\}$ are the member joint
force and member joint displacement, respectively shown in figure 2. The
member stiffness matrix or the local stiffness matrix is [K]. For a truss
element it is given by

$\left[K\right]=\frac{EA}{L}\begin{bmatrix}1&0&-1&0\\\ 0&0&0&0\\\ -1&0&1&0\\\
0&0&0&0\end{bmatrix}$ (7)

After applying the congruent transformation, the member stiffness matrix in
global coordinates are given as

$\left[K^{e}\right]=\left[T^{e}\right]^{T}\left[K\right]\left[T^{e}\right]$
(8) $K^{e}=\frac{E^{e}A^{e}}{L^{e}}\begin{bmatrix}l^{e}&lm&-l^{2}&-lm\\\
lm&-m^{2}&-lm&-m^{2}\\\ -l^{2}&-lm&l^{2}&lm\\\
-lm&-m^{2}&lm&m^{2}\end{bmatrix}$ (9)

Where $l=cos\phi^{e},m=sin\phi^{e}$ and with $\phi^{e}$ as the orientation
angle as shown in Figure 2.

####  2.1.4 Time domain solution of Equation of Motion

The equation of motion for the linear system of equations is solved with the
Newmark beta method due to its unconditional stability. The displacement and
the velocity terms for the next time step are calculated as follows:

$u_{t}=u_{t-\Delta t}+\Delta t\dot{u}_{t-\Delta
t}+\left(\frac{1}{2}-\beta\right)\Delta t^{2}\ddot{u}_{t}$ (10)
$\dot{u}_{t}=\dot{u}_{t-\Delta t}+(1-\gamma)\Delta t\dot{u}_{t-\Delta
t}+\gamma\Delta t\ddot{u}_{t}$ (11)

We follow the average acceleration approach with $\beta=\frac{1}{4}$ and
$\gamma=\frac{1}{2}$
The Newmark beta method solves the algebraic form of the equation of motion
(EOM) of undamped forced vibration at the end time interval $t+\Delta t$

$F_{t+\Delta t}=\mathbf{M}\ddot{u}_{t+\Delta t}+\mathbf{K}u_{t+\Delta t}$ (12)

The stiffness and the mass matrices are computed in the following fashion to
reduce in the form of equation (9)

$\mathbf{\hat{K}}=\mathbf{K}+a_{0}\mathbf{M}$ (13)

Where $\hat{K}$ is the effective stiffness matrix and
$a_{0}=\frac{6}{\gamma\Delta t^{2}}$
Similarly, the effective load vector at time $t+\Delta t$ is calculated as in
(17).

$\hat{F}_{t+\Delta t}=F_{t+\Delta
t}+\mathbf{M}(a_{0}u_{t}+a_{2}\dot{u}_{t}+a_{3}\ddot{u}_{t})$ (14)

Here, $a_{2}=\frac{1}{\gamma\Delta t}$ and $a_{3}=\frac{1}{2\gamma}$.
The above simplification leads to the algebraic form.

$\left\\{\hat{F}_{t+\Delta
t}\right\\}=\left[\mathbf{\hat{K}}\right]\left\\{U_{t+\Delta t}\right\\}$ (15)

From the above equation, displacement of each node is calculated for every
time step.
The natural frequency of the system is calculated as given below.

$\omega^{2}[M]\Phi=[K]\Phi$ (16) $\omega^{2}=eig([M]^{-1}[K])$ (17)

The detailed description of the theory and implementation of the dynamic
Lattice-Element method with validation and verification of the method by
analytical and numerical benchmarks is given in (Rizvi et al., 2018, 2020).

###  2.2 Wave field identification by Convolutional Neural Networks

The basic idea of deep learning system damage detection in the given case of
propagating wave fields is the identification of wave field patterns
respective of the change in wave field patterns during the damage evolution.
To apply the idea, the excitation and receiver points keep constant during the
monitoring process.

The damage evolution process covers the initial static case on the given
plate. After a change of surrounding, static stress condition damages on
different positions in the plate area can be created depending on the stress
condition and the material parameter. Before and after a damage scenario, a
small strain wave field is excited to propagate through the plate. Because of
the damage / crack, within the plate the pattern of the propagating wave field
will be modified.

The interaction of the wave field within the crack is essential for
identifying the correct wave field. Under the assumption of an open crack,
neglecting shear slipping and crack growth under dynamic loads, the crack will
produce a mode conversion and a scattering of the propagating wave field. That
phenomena is studied in (Rizvi et al., 2018, 2020). It becomes obvious that
the transient solution provides that phenomena.

For future application of the DNN-Damage-Detection of real structures, the
virtual tool has to be optimized and needs sufficient training data. The
numerical data should provide the base of the training data of the CNN.

###  2.3 Training of 1D CNN-Based Detector for Damage Detection

####  2.3.1 The Damage Detection Dataset

##### An Overview of Training Data Set - The Wave Propagation Data

The wave field data was produced by numerical simulation on virtual 2D plates
as study examples. After the application of a Dirac input, the elastic wave
propagates from the source point into the plate. Adding source points at a
different location of the plate generates different wave patterns, which are
recorded by a set of receiver points. Particularly when the plate has a crack
inside, the wave patterns are very different from the one in non-crack plates.
Because of the crack, the wave field shadow as well as the reflection becomes
visible and change the pattern of the wave field around the damaged region.

The plate has fixed boundary conditions at the bottom and no constrains at
other sides. The wave source can change its location around the plate
boundary. Receivers that record wave displacement are assumed to be along the
free boundary as well as on the inner surface of the plate.

The sequential measures of time-dependent displacement amplitudes are used as
the data source for the planned damage detection network. The instantaneous
load added at a chosen excitation point causes high displacement amplitudes at
the wave field and decreases rapidly after several time steps to the wave coda
at a smaller strain level. The observable surface wave front, as well as at
interference of back scattered and reflected wave field is the result of the
wave propagation and reflection in the plate. The described phenomenon are
clearly visible in an example in Section 3.1.

In this paper, the whole content of the wave field, the initial wave front,
the coherent part and the diffusive part will be used for the damage detection
in the time domain. There is no selection respective analysis of harmonic wave
modes in the part coherent wave field (Wuttke et al., 2012b) or application of
the interferometric method at the code in the diffusive part (Wuttke et al.,
2012a) of the wave field yet.

##### The Damage Detection Dataset

The dataset is generated by running the numerical simulation repeatedly for
randomly generated plates with or without a crack. In this study, the size of
all plates is set to $0.01(m)\times 0.01(m)$ and the lower-left plate corner
is defined as the Cartesian coordinate system origin. The wave field is
generated by an initial Dirac impulse (in 1 times step span). The simulation
runs for 2000 samples with time steps by $1e-9$ seconds. With the excitation
of the wave field, the recording at all receiver points starts as the base of
the deep detection network. During the method development, possible negative
effects by larger displacement values on deep learning model are analysed. The
plate, receivers and excitation points are shown in Figure 1. The resulted
displacement wave field consists of time histories in the X- and Y-direction
at 81 receiver positions with 2000 time steps.

To validate the damage detection method with randomly generated cracks, the
crack itself is described by 3 parameters, i.e., crack length $l$, orientation
$\alpha$, and start position $(x,y)$, with their values being randomly chosen
from the following description,

$\displaystyle l$ $\displaystyle\in(0,\frac{1}{2}\min(e_{x},\,e_{y})],$ (18)
$\displaystyle\alpha$ $\displaystyle\in[0,360],$ (19) $\displaystyle x$
$\displaystyle\in[s_{x},\,e_{x}-s_{x}],$ (20) $\displaystyle y$
$\displaystyle\in[s_{y},\,e_{y}-s_{y}],$ (21)

where $e_{x}$ and $e_{y}$ are the length of the sample edges along the x-axis
and y-axis, $s_{x}$ and $s_{y}$ are the distance between two receivers in the
X-axis and Y-axis (see Figure 1). If one randomly generated crack stretches
out of the sample plate, the excess part is discarded. The plate particles
that correspond to the crack are marked as removed for the Lattice-Element
model calculation.

To start the identification scenario, i.e., detecting the damage in a given
plate, a binary image of the damage was provided as its label for comparison
between the identified structure and original structure. The binary image
covers the plate’s surface and indicates the location where a crack exists.
The label image can be obtained originally of an 100x100 resolution, where
each pixel covers an area similar to the size of a particle. When the model is
adjusted to refine or enlarge predictions, the resolution of the label image
can be changed accordingly. Figure 2 gives out two label images of different
resolutions for the same plate. The image of 16x16 pixels is resized and
binarized from the image of 100x100 pixels. We use the image of the reduced
resolution (16x16) as the supervision signal for training the detector
network. As the proposed model does predictions for each pixel, the label
image of 16x16 resolution restricts the problem scale while still maintaining
the model’s applicability.

Figure 1: The receiver arrangement on the surface of a plate. Figure 2: Plate
with crack and its binarised labels at different resolutions. (A). The
generated sample with one crack inside in particle view. (B). The 100x100
binary label image for the plate (white pixels indicate the existence of the
crack). (C). The 16x16 binary label image for the plate.

The numerical simulation randomly generates plate particles and cracks within
the plate domain using pre-defined boundary conditions. These randomnesses in
sample generation and crack generation are called Type-N samples.
Additionally, plates without any crack inside are also generated randomly as
reference samples, marked as Type-R. It’s also possible to generate samples of
different plates with similar cracks (Type-S) or the same plates with
different cracks (Type-C). In this work, 3040 samples for training and 320
samples for testing were generated. In the training dataset, nearly half of
the cases are Type-N; Type-R, Type-S, and Type-C cases consist of the
remaining portion of the training data. The testing data consists of 144
Type-N samples, 8 Type-R samples, 160 Type-C samples, and 8 Type-S samples.
The Type-C samples are generated from 10 different samples and 16 different
cracks for each sample (see Figure 20 in the Appendix). Among all test
samples, we intentionally generated 7 random samples, and each one has its
counter case in the training dataset in terms of the same crack (no-crack
cases are excluded in this case). It is worth emphasizing that these samples
are not repeated ones. Because of the randomness in the generation process,
the diversity of the interior particles and their wave field patterns is
ensured.

####  2.3.2 Crack Detection Models with CNN Detector

The proposed damage detection model is trained to detect the exciting crack in
the 2D plate and the damage location on hand of the wave field pattern. The
training covers a series of wave fields in randomly generated plates with or
without cracks. The receivers are placed inside the plate and along the free
boundaries to record the coordinate-depending wave field time histories. Here,
at the receiver points, the displacements in x- and y-directions are recorded.
These typical 1D Euclidean time series are the data base for the 1D
convolution filters within the feature extractors. The proposed deep network
model consists of three components: a set of 1D-CNN layers acting as a wave
pattern (WP) extractor to select WPs from input displacement time histories;
two fully convolutional layers to handle the time dimension and on the
receiver’s dimension to fuse wave patterns; and a predictor module taking the
fused features and making predictions of crack existence (Fig 3).

Figure 3: Conceptual explanation of the proposed crack detection model with
1D-CNN detector. (A). Diagram of a 1D-CNN that transforms a discretized wave
input and produces a discrete feature sequence. (B). Internal structure of a
1D-CNN layer, consisting of a trainable weight $W^{(h)}$, a bias $b^{(h)}$.
Activation function is represented by $\sigma$.(C) Diagram of the complete
structure of the proposed model.

Let $N$ denote the number of receivers that are placed on the plate, they
record displacements for $T$ time steps after a load excitation is added. The
reading of the $i^{\mathrm{th}}$ receiver at time $t$ is denoted as
$s_{k}^{(t)}$, where $k\in{1,2,\cdots,N}$ and $t\in{1,2,\cdots,T}$. The
surface of the plate is decomposed into regular grids. Each cell within the
grid covers a small area of the plate surface. The number of cells indicates
the spatial resolution of the prediction model. The size of the cell was
chosen as 10 times smaller than the wave lengths. Finer resolution requires a
larger number of cells and each cell covers a smaller area on the sample
surface. In contrast, a coarser resolution results in less number of cells and
larger coverage for a single cell. The resulting cells can be denoted by the
column and row index as $c_{i,j}$. The model makes binary classification on
every grid cell to decide whether damage exists in the cell. To summarize, the
model can be written as the following equation, where $f$ presents the
proposed model and $\theta$ are all trainable parameters; $p_{i,j}$ is the
probability of damage existence in $c_{i,j}$.

$f_{\theta}(s_{1}^{(1)},\cdots,s_{k}^{(t)},\cdots,s_{N}^{(T)})=\\{p_{i,j}\\}$
(22)

##### CNN

is a feed-forward network that consists of trainable multistage feature
extractors. The feature extractors are trained in an optimization procedure,
where the gradient of an objective function with respect to the weights of the
feature extractors is calculated and back-propagated (Rumelhart et al., 1986).

CNN is particularly useful to analyse natural signals which can be represented
in arrays of different dimension modalities, sequential signals including
language and audio as 1D arrays; images as 2D; and video and volumetric images
as 3D. The core operation of CNN is to calculate the convolution of input
signals with a set of trainable filters. The transformation produces different
kinds of features from input signals in spatial (e.g., images) or temporal
(e.g., audio) domain (LeCun et al., 1995).

In the implementation, CNN differentiates itself from other ANNs by using the
local connection and weight sharing strategy. In CNNs, one “neuron” connects
locally with only a restricted number of “neurons” in its previous layers, and
the connection weights are shared among all the “neurons”. Local connectivity,
weight sharing, and pooling, are CNN’s three key properties for dealing with
natural signal (LeCun et al., 2015). A typical CNN architecture is composed of
firstly some convolutional layers (ConvLayer), and then some more ConvLayer,
or fully connected layers. ConvLayer detects local features from the previous
layer by transforming the input signal with a set of filters. It produces
different kinds of feature maps with respect to its filters, then an
activation operation is applied to the feature maps. The non-linear activation
functions “squeeze” the values of a feature maps into a certain range, mostly
$[0,1]$, $[-0.5,0.5]$ or $[0,+\infty)$. Sometimes, a pooling layer is used for
down-sampling the feature maps by taking local average or local maximum
values. The pooling layer merges semantically similar features into a higher
level (LeCun et al., 2015). The Pooling layer can be intentionally replaced by
setting a larger stride in the convolution layer (Springenberg et al., 2015).

Figure 4 depicts a schematic drawing of applying n 1d-convolution to N sensory
input of T steps. The kernel size is denoted as m. For each receiver’s input,
the convolution produces n features. The red mark indicates the data patch
used for convolution and the corresponding results in feature maps.

Figure 4: The schematic drawing of 1D convolution on N receiver data in T
steps.

##### The Network Structure

The detailed implementation of the proposed model is described in this part.
The WP-Extractor consists of three 1D-ConvLayer blocks. Each block has two
1D-ConvLayers followed by a BatchNormalization layer ((Ioffe and Szegedy,
2015)) and an activation layer - using the LeakyReLU function. ((He et al.,
2015)). The BatchNormalization layer standardizes the output of the ConvLayer
by re-scale according to the samples in a batch. It prevents dramatic change
of the spread and distribution of inputs during training, and has the effect
of stabilizing and speeding-up the training process. LeakyReLU is a variation
of Rectified Linear Unit activation function (ReLU). The ReLU function results
in zero when the input is less than zero and keeps the input unchanged when
the input is above or equal to zero. The LeakyReLU function “squeezes” the
value when the input is less than zero and thus allows a small, non-zero
gradient when the unit is not active.

The configurations of the two ConvLayers in the same block are identical while
the kernel size and number of filters for ConvLayers vary between blocks. To
reduce the size of the extracted wave patterns, the output of each
1D-ConvBlock is passed through a MaxPooling layer. Because the 1D-ConvLayers
operate only on the time dimension, the WP-Extractor selects only the receiver
patterns. For each case, the input data has shape in $N\times T$, where $N$
represents the number of receivers and $T$ the time steps. The output of the
WP-Extractor has shape $N\times\hat{T}\times C$, where $C$ is the filter
number of the last ConvBlock.

The first fusion layer transforms the receiver data into a 1D vector by
passing through a fully convolutional layer. The convolution kernel of this
layer is $1\times\hat{T}$, so the transformed field of the ConvLayer covers
the whole time domain. The resulted data transformations are given by an
$N\times 1\times C^{\prime}$ array. After the transformed data are reshaped
according to the position of the receivers and has the size of an $9\times
9\times C^{\prime}$ array ($N=9\times 9$ as shown in Figure 1). The second
fusion layer, a 2D fully ConvLayer is used to save the information of all
receivers. The 2D ConvLayer employs a kernel of $9\times 9$ and thus produces
an $1\times 1\times C^{\prime\prime}$ array.

The core module of the predictor is composed of two Transpose-Convolutional
layers (TransConvLayer, sometimes called as deconvolution). The two
TransConvLayers are used to up-sample the saved transformed data.
TransConvLayers are widely used in image generation tasks, such as DCGAN
((Radford et al., 2016)). The saved transformed data are first reshaped to
$4\times 4\times c$, and then passed through the two TransConvLayers. The
ready-for-predict data shape is $16\times 16\times 4$. Finally, it’s passed
through a 2D ConvLayer of $1\times 1$ kernels, and a single output channel.
The sigmoid function is used in this ConvLayer as activation to make sure the
prediction ranges from 0 to 1.

The layer configuration and wave pattern shapes are also shown in Figure 5, K
refers to the kernel size, S is the steps, F indicates the number of filters.
The step of Max-Pooling is set to 4. The 1D convolution can be implemented
using 2D convolution by fixing the kernel size of the receiver dimension to 1.

Figure 5: Detailed design of the damage detection model. (A). The network
architecture. (B) Shapes of input and output of each component.

##### The Loss Function

Training the proposed model is an optimization procedure, and relies on the
objective function (also known as loss function). In this work, the proposed
model makes a prediction on crack existence for each single patch. In other
words, the model makes multiple binary predictions. For single binary
classification problems, cross entropy (CE) loss (see Equation 23) is probably
the most commonly used loss function. However, considering the “has crack”
patches actually consist of only a small portion of total patches, cross
entropy (CE) loss can introduce bias towards “no crack” predictions, i.e.,
simply predicting all patches as “no crack” already result a rather low loss
value. To tackle such extreme class imbalance, we select Focal Loss (FL) as
our loss function. FL was originally proposed to address the extreme class
imbalance in object detection (Lin et al., 2017). FL is a variation of CE loss
by adding a penalty term to reduce the loss value of already correctly
predicted training cases. The penalty term $(1-p_{t})^{\gamma}$ ($\gamma\geq
0$) re-weights between difficult and easy examples. During training, if a
sample is already predicted correctly with a high probability, it is called an
“easy” case. The penalty term reduces its loss and thus focus on “hard” cases,
where correct predictions are made with a much lower probability.

$L_{\mathrm{CE}}=-\frac{1}{N}\sum^{N}_{i=1}log(p_{t}^{(i)}),\quad\mathrm{where}\;p_{t}=\begin{cases}p&\;\mathrm{if}\;y=1,\\\
1-p&\;\mathrm{otherwise}\end{cases}$ (23)

where $y$ is the label (1 and 0).

To adjust the loss values of the two binary classes, a weighting factor
$\alpha\in[0,1]$ can be added. Similar to defining $p_{t}$, $\alpha_{t}$ can
be defined as $\alpha$ for class 1, and $1-\alpha$ for class 0. The focal loss
is written as:

$\mathrm{FL}(p_{t})=-\alpha_{t}(1-p_{t})^{\gamma}log(p_{t}),$ (24)

For the crack detection case, let $\boldsymbol{y}$ denote the binary image and
$\hat{\boldsymbol{y}}={p_{i,j}}$ a predicted image. The average FL on N cases
is calculated by:

$\mathrm{FL}(p_{t})=\frac{1}{N}\sum_{n=1}^{N}\sum_{i=1}^{U}\sum_{j=1}^{V}-\alpha_{t}(1-p_{t\\{n,i,j\\}})^{\gamma}log(p_{t\\{n,i,j\\}}),$
(25)

where $U$, $V$ are the number of columns and rows of a grid.

Two hyper-parameters are introduced by Focal loss, $\alpha$ and $\gamma$. When
$\gamma=0$, FL is equivalent to weighted CE. When $\gamma$ increases, the
modulating factor uses a lower standard to define easy examples and has
greater power to down-weights well-classified examples. In practice, the
optimal $\alpha$ and $\gamma$ are found out by empirical studies.

####  2.3.3 Training Process

##### Data Pre-processing

In this simulation, the recorded data at the edge receivers are discarded to
avoid any possible effects caused by the extremely large values. Thus, in
total the records at 81 receivers are used for both training and testing. Then
the wave displacements are normalized between -1 and 1 according to each
sample’s maximum and minimum value. The resulting input data for the CNN model
is $2000\times 81\times 2$ matrix for each case.

##### Training Configurations

The Adam optimization algorithm (Kingma and Ba, 2015) is chosen as the
optimizer; it’s a commonly used first-order gradient-based optimization
algorithm for training deep networks. When updating model parameters during
each training step, the algorithm adjusts the learning strength according to
the previous gradients. An overview of the optimizers in deep learning can be
found in (Ruder, 2016). The initial learning rate is set to 0.0002. The
training epochs were set to 150 for all experiments to ensure sufficient
training steps for the models to converge. The best model with respect to the
evaluation metric is saved for evaluation.

As mentioned in the previous section, the focal loss has two hyper-parameters,
$\alpha$ and $\gamma$. To determine suitable hyper-parameters, a set of models
were trained with different alpha and gamma. The detailed selection of alpha
and gamma is listed in section 3.3.1.

The model and simulation code is implemented in Python with Tensorflow 2.1
Keras. The simulations are performed on a workstation of Windows 10 platform
with Nvidia GPU.

##### Evaluation Metrics

In this section, we introduce the Dice similarity coefficient (DSC) metric
(defined as $M_{DSC}=2\cdot\frac{recall\cdot precision}{recall+precision}$)
and the IoU-based accuracy for model evaluation. Although the loss values
indicate the quality of the prediction on the patch basis, predicting cracks
is more focused. We expect the model to find out the “has crack” patches cover
as much damaged patches as possible, and as few “no crack” patches as
possible.

In the damage detection study, if the model correctly predicts a cell has (has
no) damage, the result is marked as TP (PN); if the model makes wrong
predictions on damage existence, the result is either FP or FN as shown in
Table 1. Based on this, we can calculate precision ($TP/(TP+FP)$) and the
recall ($TP/(TP+FN)$),

| Right Prediction (T) | Wrong Prediction (F)
---|---|---
Has Damage (P) | TP | FP
No Damage (N) | TN | FN
Table 1: Prediction typology in a binary-classification-based damage
detection.

then the DSC metric can be calculated with Equation 26. Similarly, the IoU
calculates the ratio between the intersection and the union of ground truth
damage area (TDA) and predicted damage area (PDA). As illustrated in Figure 6,
TDA is bound by dashed red lines; PDA is marked by solid red lines. Their
intersection is the TP set, where the model makes correct predictions. The
remaining part of PDA is wrongly identified as damaged area, i.e., the FP set;
the remaining part of TDA that has damage inside but has not been found forms
the FN set. For a single case, the IoU can be calculated as Equation 27 by
counting the number of each set.

Both DSC and IoU metric range between 0 and 1. If there is no overlap between
PDA and TDA, both metrics are equal to 0. When PDA is closer to TDA, the
intersection area becomes larger and the union area becomes smaller, resulting
in a value closer to 1. When PDA covers TDA exactly, DSC and IoU metrics reach
their upper limit of 1. When the sample has no damage and the model makes
correct predictions, both the intersection and the union become 0. In this
special case, both DSC and IoU metrics are assigned the value 1.

$M_{\mathrm{DSC}}=2\cdot\frac{\mathrm{area}(TP)}{2\cdot\mathrm{area}(TP)+\mathrm{area}(FP)+\mathrm{area}(FN)}$
(26)

$\mathrm{IoU}=\frac{\mathrm{area}(TP)}{\mathrm{area}(TP)+\mathrm{area}(FP)+\mathrm{area}(FN)}$
(27)

Figure 6: Conceptual example of calculating IoU.

One underlying assumption of data generation is that each sample has maximum
“one crack” inside. Based on the assumption, we can define the accuracy using
IoU values. For the prediction of a single case, we can consider it as a
“correct” prediction if its IoU is greater than a given threshold. Given the
threshold, the accuracy on the whole dataset is calculated as the ratio of the
number of samples whose IoU value is greater than the threshold, to the total
number of evaluated dataset.

##  3 Results

###  3.1 Simulated Displacement Wave Field Using Dynamic Lattice

The dynamic lattice method described in section 2.1 is used to simulate the
wave fields in a 2D plate. The considered boundary conditions with different
excitation points and crack conditions are shown in Fig.7. Fig.8 shows the
simulated wave fields in lateral direction for a boundary condition of Fig.7a.
The simulated wave fields with a generated crack (Fig.7b) are shown in Fig.9.
For the conditions shown in Fig.7c - excitation point in upper middle boundary
- the wave field is plotted Fig.10. The results clearly show wave shadows
behind the crack as well as the reflection of the wave field from the defined
cracked surface.

Figure 7: Boundary conditions: (a) horizontal excitation without generated
crack, (b) horizontal excitation with generated crack, and (c) vertical
excitation with generated crack Figure 8: The 6 frames (100 time steps
interval, from left to right) of a displacement ($u_{x}$) wave propagation
inside the defined plate in Fig.7a. Figure 9: The 6 frames (100 time steps
interval, from left to right) of a displacement ($u_{x}$) wave propagation
inside the defined plate in Fig.7b. Figure 10: The 6 frames (100 time steps
interval, from left to right) of a displacement ($u_{y}$) wave propagation
inside the defined plate in Fig.7c.

The time histories ($u_{x}$) of the reference points ($R_{1:9}$) inside the
plate (Fig.11) are shown in Fig.12. Two boundary conditions are considered:
one with discontinuity (crack), and one without. The plate dimension is 10x10
$cm$ and the load excitation is at the left middle boundary. The applied
rectangular impulse load with a magnitude of 1 $kN$ is kept for 10 time steps,
where $\Delta t=0.00000001s$. The Young’s modulus of a plate is assigned to 5
$GPa$.

Fig.12 clearly shows the arrival time of the wave fields at each reference
point. The closest reference point ($R_{2}$) has the maximum amplitude and
minimum arrival time. In Fig.12a the arrival time of the wave field to
$R_{1}\approx R_{3}$, $R_{4}\approx R_{6}$ and $R_{7}\approx R_{9}$. Due to
the generated discontinuity (crack) in Fig.12b, the first arrival times of the
wave field to $R_{3}$, $R_{6}$ and $R_{9}$ are delayed. Theses reference
points are located in the shadow field behind the generated discontinuity.
Having a closer look at the $R_{2}$, it is obvious that due to the wave
reflection from the generated discontinuity, the arrival of the second wave
field happens sooner than the first boundary condition, approximately
$1.45\times 10^{-6}s$. The length, location and orientation of the
discontinuities affect the wave field in the domain. The simulated wave fields
at the reference points are used for training and developing the artificial
neural network model.

Figure 11: The boundary conditions and assigned reference points: (a) without
crack, (b) with a generated crack Figure 12: The time histories ($u_{x}$) of
the reference points inside the plate: (a) without crack, (b) with a generated
crack

###  3.2 The Trained Damage Detection Model

####  3.2.1 Identified Optimal Hyper-parameters for Focal Loss

FL introduces two hyper-parameters $\gamma$ and $\alpha$. They are used to
adjust the loss value of a prediction during the training, so the training can
focus on specific types of training cases. Figure 13 shows the FL value with
different $\gamma$ and $\alpha$ values. Figure 13 (A) is a remake of Figure 1
in (Lin et al., 2017). With use of the penalty term, the loss value is reduced
with the probability of making correct predictions increasing. $\gamma$
controls the decay strength, and larger $\gamma$ ensures the loss to decrease
faster. For example, when $\gamma=5$, the predictions where $p_{t}>0.4$ can
hardly contribute to the loss. In contrast, when $\gamma=1$, the predictions
where $p_{t}>0.6$ still contribute to the loss. Meanwhile $\alpha$ can also be
used to re-weight the binary classes (has crack and has no crack) (Figure 13
(B)-(D)). When $\alpha$ is used for one class, the other class is re-weighted
by $(1-\alpha)$. Choosing a small $\alpha$ for a class will obviously decrease
the contribution of the whole class to the loss. For example, if $\alpha=0.1$
is chosen for the “has crack” class and $\gamma=5$, the predictions can be
hardly improved when it is greater than $0.3$. Particularly, when $\gamma=0$
and $\alpha=0.5$, the FL is equivalent to the cross-entropy (CE) loss.

Figure 13: The focal loss values with different $\gamma$s and $\alpha$s. (A).
Focal loss without weight $\alpha$. (B). Focal loss with different $\alpha$
value with fixed $\gamma=0$. (C). Focal loss with different $\alpha$ value
with fixed $\gamma=2$. (D). Focal loss with different $\alpha$ value with
fixed $\gamma=5$

In this paper, $\gamma$ and $\alpha$ control the model’s learning strength on
“no crack” class and “has crack” class. As $\alpha$ increases, the model is
driven to focus on damaged cells, because the false predictions of damaged
cells contributes more to the overall loss. As $\gamma$ increases, the model
is trained to focus on ”hard” cells, where the model can’t make predictions
with high confidence. Because a higher $\gamma$ value forces the model to pay
more attention to the ”hard” cells. The evaluations of models that are trained
with different $\alpha$s are given in Table 2. When assigning larger
weights(larger $\alpha$) to the “has crack” class for CE loss, the trained
model tends to have lower precision and higher recall. This can be interpreted
as the model’s tendency to give more “has crack” predictions. On the contrary,
using small weights(smaller $\alpha$) for the “has crack” class results in
higher precision but lower recall. This means the models tend to give less
“has crack” predictions. When $\alpha=0.9$, the CE loss gives the model of the
highest accuracy, however, the model is also characterised as having low
precision and high recall; having high DSC metric value but not the optimal
one. The results using focal loss are shown in Table 3. By adding the penalty
term, with carefully chosen $\gamma$, the trained models have balanced the
precision and recall, and thus result in an increasing in IoU and DSC metrics.
The accuracy is also improved compared to the models trained with CE loss. We
consider two sets of $\alpha$ and $\gamma$ combinations, $\alpha=0.35$
$\gamma=0.2$ and $\alpha=0.9$ $\gamma=0.4$, have balanced precision and
recall, and achieve high DSC and accuracy at the same time.

$\alpha$ | prec. | recall | IoU | DSC | accu.
---|---|---|---|---|---
0.1 | 0.852 | 0.606 | 0.549 | 0.709 | 0.600
0.2 | 0.917 | 0.666 | 0.629 | 0.772 | 0.697
0.25 | 0.864 | 0.603 | 0.551 | 0.710 | 0.609
0.3 | 0.895 | 0.645 | 0.600 | 0.750 | 0.672
0.35 | 0.906 | 0.705 | 0.656 | 0.793 | 0.706
0.5 | 0.882 | 0.675 | 0.619 | 0.765 | 0.666
0.75 | 0.756 | 0.762 | 0.611 | 0.759 | 0.706
0.9 | 0.839 | 0.734 | 0.643 | 0.783 | 0.738
0.95 | 0.834 | 0.737 | 0.643 | 0.783 | 0.691
Table 2: The evaluation results (including IoU, DSC, and accuracy) of models trained by varying $\alpha$ for CE loss($\gamma$= 0). $\gamma$ | $\alpha$ | prec. | recall | IoU | DSC | accu.
---|---|---|---|---|---|---
0.1 | 0.25 | 0.916 | 0.668 | 0.630 | 0.773 | 0.703
0.2 | 0.35 | 0.880 | 0.730 | 0.664 | 0.798 | 0.756
0.4 | 0.9 | 0.803 | 0.779 | 0.654 | 0.791 | 0.769
1 | 0.5 | 0.884 | 0.711 | 0.650 | 0.788 | 0.731
2 | 0.75 | 0.846 | 0.746 | 0.657 | 0.793 | 0.725
4 | 0.95 | 0.816 | 0.748 | 0.640 | 0.781 | 0.722
Table 3: The evaluation results (including IoU, DSC, and accuracy) of models
trained by varying $\gamma$ for FL loss (optimal $\alpha$).

####  3.2.2 The Selected Thresholds

The accuracy is calculated dependently with two threshold settings: the
threshold of crack existence in a pixel, and the threshold of correct
predictions of a sample. The first threshold defines the probability value,
above which a pixel can be considered to have a crack inside. In this work, it
is also referred to as binarizing threshold ($T_{bin}$). The second threshold
is set as a “tolerance” (marked as $T_{tol}$) to the prediction. The
“tolerance” allows a prediction to be “correct” when the predicted “has crack”
pixels cover a certain area of the crack, i.e., its IoU score is greater than
the threshold. The very strict criteria requires that the predicted “has
crack” patches cover the true crack-existing area, i.e., the $\mathrm{IoU}=1$,
to be a “correct prediction”.

The FL function pushes the predicted probabilities of “has crack” and “no
crack” towards opposite extremes, because a sample with IoU value that is
close to 0.5 will have a large penalty during training. This fact is also
illustrated in Figure 14 and 15, which are resulted from the recommended model
trained with $\alpha=0.9$ $\gamma=0.4$ and $\alpha=0.35$ $\gamma=0.2$. The
sub-figures A of both two figures suggest that most damaged cells are
correctly predicted with a probability above 0.5, while still-minor “hard
cases” get a border prediction around 0.5, with about 45 cases that both
models can not properly handle. In both sub-figures B, the accumulated
histograms show a clearer comparison on the quality of predictions for
different $T_{bin}$ values. They show that different $T_{bin}$ produce similar
accumulative histogram curves. This suggests that most ”no crack” cells and
many ”no crack” cells are predicted with very high confidences. We can choose
$T_{bin}=0.5$ as it also fits the configuration of FL loss. The curves begin
to rise when IoU value reaches 0.5. This suggest us to chose $T_{tol}$ for
evaluation, so that the number of cases with IoU values between 0 to 0.5 are
relatively small and accumulate quickly when
$\mathrm{value}_{\mathrm{IoU}}>0.5$.

Figure 14: The distribution of IoU values among 320 test cases. The IoU values
are calculated by the model trained with $\gamma=0.4$ and $\alpha=0.9$. A: the
6 histograms of the IoU values that are calculated with 6 different
binarization threshold ($T_{bin}$=0.3, 0.4, 0.5, 0.6, 0.7, 0.8); B: the
accumulated histograms of the IoU values that are calculated with 6 different
threshold ($T_{bin}$=0.3, 0.4, 0.5, 0.6, 0.7, 0.8). Figure 15: The
distribution of IoU values among 320 test cases. The IoU values are calculated
by the model trained with $\gamma=0.2$ and $\alpha=0.35$. A: the 6 histograms
of the IoU values that are calculated with 6 different binarization threshold
($T_{bin}$=0.3, 0.4, 0.5, 0.6, 0.7, 0.8); B: the accumulated histograms of the
IoU values that are calculated with 6 different threshold ($T_{bin}$=0.3, 0.4,
0.5, 0.6, 0.7, 0.8).

####  3.2.3 Discussion on Model Performance

The histograms on the distribution of test data (Figure 14 and 15) indicate
that both models are not good at detecting a minor set of damaged cases in
test data. We first examine the distribution of crack size (Crack size is the
pixel count or the percentage of “has crack” pixels in the 100x100 labeling
image.) in training data and test data. As shown in Figure 16, the samples
with small cracks consist of a larger portion in test data then in training
data. To further explore the relation between crack size and model
performance, we plot the histogram of crack size against IoU values for test
data (Figure17). We can easily find out that the IoU values can be very low
for tiny crack samples, whereas the IoU values for larger crack samples are
mostly above 0.5. Then the accuracy, adjusted by excluding samples with small
crack size from test data, is shown in Figure 18. It shows that the proposed
model is particularly good at identifying larger cracks. IoU values are
calculated from the model predictions with $T_{bin}=0.5$. It indicates that
the most low-quality predictions are made for samples with tiny cracks, while
cases of larger crack sizes generally have better predictions. This means the
developed model can easily distinguish between damaged cases and non-damage
cases for large cracks but is not good at detecting tiny cracks. If we only
count the cases with crack size greater than 0.002, the accuracy leaps by
around 0.1. When excluding the cases with small cracks (crack size less than
0.003), the accuracy of the proposed model is already beyond 0.9. If we only
count the cases with larger cracks (crack size greater than 0.004), the
accuracy of both models can reach 0.95.

Figure 16: Histogram of crack size distribution in training data and testing
data. A: crack size distribution in training data; B: crack size distribution
in testing data. Figure 17: Histogram of crack size distribution and IoU
values for test data. A: results from model trained with $\alpha=0.9$
$\gamma=0.4$; B: results from model trained with $\alpha=0.35$ $\gamma=0.2$.
Figure 18: The adjusted accuracy calculated for test data after excluding tiny
crack cases. A: predictions are made by the model with $\alpha=0.9$
$\gamma=0.4$; B: predictions are made by the model with $\alpha=0.35$
$\gamma=0.2$. The line plot presents the accuracy that is re-calculated when
excluding cases with crack size less than 0.001, 0.002, 0.003, and 0.004; the
bar chat shows the number of samples after excluding the samples with tiny
cracks against the total number of samples in test dataset.

##  4 Conclusion

The paper presents a new approach to detect damages by wave pattern
recognition models. The major development is a learning CNN to detect on-hand
the visible wave pattern of the damaged zone within a solid structure. To
generate the cracked structure, a new dynamic Lattice Element method was used.
The major advantage of this method is the application to heterogeneous
structures under mechanical, hydraulically, thermal field influence and local
chemical changes to describe the evolution of damages in solid structures. The
use of new generation deep CNNs to analyse the time dependency within the
changed wave pattern is promising. With the described method, a stable
detection of 90 percent of the generated large cracks was possible. The next
steps will be the reduction of the used number of receivers and increasing the
model’s ability of tiny crack detection .

## Acknowledgement

The presented work is supported by the Kompetenzzentrum Geo-Energie (KGE)
(https://www.kge.uni-kiel.de/, last access: 14 February 2021). We gratefully
acknowledge the funding of KGE by the European Union - European Regional
Development Fund (EFRE). We also thanks the funding of the project ’Seismic
identification of concrete structures’ funded by the Federal Ministry of
Economic Affairs and Industry - BMWI and the German Federation of Industrial
Research Associations - ZIM/ AIF with funding code ZF4016813FF8. Furthermore,
we would like to acknowledge the thoughtful reviews of the reviewers and the
editor, and their constructive comments supporting the manuscript revision.

## References

  * Abdeljaber et al. [2017] Osama Abdeljaber, Onur Avci, Serkan Kiranyaz, Moncef Gabbouj, and Daniel J Inman.  Real-time vibration-based structural damage detection using one-dimensional convolutional neural networks.  _Journal of Sound and Vibration_ , 388:154–170, 2017.
  * Avci et al. [2021] Onur Avci, Osama Abdeljaber, Serkan Kiranyaz, Mohammed Hussein, Moncef Gabbouj, and Daniel J. Inman.  A review of vibration-based damage detection in civil structures: From traditional methods to machine learning and deep learning applications.  _Mechanical Systems and Signal Processing_ , 147:107077, 2021.  ISSN 0888-3270.  doi: https://doi.org/10.1016/j.ymssp.2020.107077.  URL http://www.sciencedirect.com/science/article/pii/S0888327020304635.
  * Farhangdoust and Mehrabi [2019] Saman Farhangdoust and Armin Mehrabi.  Health monitoring of closure joints in accelerated bridge construction: A review of non-destructive testing application.  _Journal of Advanced Concrete Technology_ , 17(7):381–404, 2019.  doi: 10.3151/jact.17.381.
  * Girshick [2015] Ross Girshick.  Fast r-cnn.  In _Proceedings of the IEEE international conference on computer vision_ , pages 1440–1448, 2015.
  * Gulgec et al. [2019] Nur Sila Gulgec, Martin Takáč, and Shamim N. Pakzad.  Convolutional neural network approach for robust structural damage detection and localization.  _Journal of Computing in Civil Engineering_ , 33(3):04019005, 2019.  doi: 10.1061/(ASCE)CP.1943-5487.0000820.
  * Guo et al. [2020] Tian Guo, Lianping Wu, Cunjun Wang, and Zili Xu.  Damage detection in a novel deep-learning framework: a robust method for feature extraction.  _Structural Health Monitoring_ , 19(2):424–442, 2020.  doi: 10.1177/1475921719846051.  URL https://doi.org/10.1177/1475921719846051.
  * He et al. [2015] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Delving deep into rectifiers: Surpassing human-level performance on imagenet classification.  In _Proceedings of the IEEE international conference on computer vision_ , pages 1026–1034, 2015.
  * Ioffe and Szegedy [2015] Sergey Ioffe and Christian Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  In _International Conference on Machine Learning_ , pages 448–456, 2015.
  * Kaewunruen and Remennikov [2006] S. Kaewunruen and A. Remennikov.  Non-destructive testing (ndt): A tool for dynamic health monitoring of railway track structures.  2006\.
  * Khan et al. [2019] Asif Khan, Dae-Kwan Ko, Soo Chul Lim, and Heung Soo Kim.  Structural vibration-based classification and prediction of delamination in smart composite laminates using deep learning neural network.  _Composites Part B: Engineering_ , 161:586 – 594, 2019\.  ISSN 1359-8368.  doi: https://doi.org/10.1016/j.compositesb.2018.12.118.  URL http://www.sciencedirect.com/science/article/pii/S1359836818325411.
  * Kim [2014] Yoon Kim.  Convolutional neural networks for sentence classification.  In _Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , pages 1746–1751, 2014.
  * Kingma and Ba [2015] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  In Yoshua Bengio and Yann LeCun, editors, _3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings_ , 2015.
  * LeCun et al. [1995] Yann LeCun, Yoshua Bengio, et al.  Convolutional networks for images, speech, and time series.  _The handbook of brain theory and neural networks_ , 3361(10):1995, 1995.
  * LeCun et al. [2015] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep learning.  _nature_ , 521(7553):436–444, 2015.
  * Lin et al. [2017] Tsung-Yi Lin, Priya Goyal, Ross Girshick, Kaiming He, and Piotr Dollar.  Focal loss for dense object detection.  In _The IEEE International Conference on Computer Vision (ICCV)_ , Oct 2017.
  * Moukarzel and Herrmann [1992] C. Moukarzel and H. J. Herrmann.  A vectorizable random lattice.  _J. Stat. Phys._ , 68:911–923, 1992.
  * Nunes et al. [2020] Lorena Andrade Nunes, Rafaelle Piazzaroli Finotti Amaral, Flávio de Souza Barbosa, and Alexandre Abrahão Cury.  A hybrid learning strategy for structural damage detection.  _Structural Health Monitoring_ , 2020.  doi: 10.1177/1475921720966943.
  * Radford et al. [2016] Alec Radford, Luke Metz, and Soumith Chintala.  Unsupervised representation learning with deep convolutional generative adversarial networks.  In Yoshua Bengio and Yann LeCun, editors, _4th International Conference on Learning Representations, ICLR 2016, San Juan, Puerto Rico, May 2-4, 2016, Conference Track Proceedings_ , 2016.
  * Rautela and Gopalakrishnan [2019] Mahindra Rautela and S Gopalakrishnan.  Deep learning frameworks for wave propagation-based damage detection in 1d-waveguides.  In _Proceedings of the 11th International Symposium NDT in Aerospace_ , 2019.
  * Rautela and Gopalakrishnan [2020] Mahindra Rautela and S. Gopalakrishnan.  Ultrasonic guided wave based structural damage detection and localization using model assisted convolutional and recurrent neural networks.  _Expert Systems with Applications_ , page 114189, 2020.  ISSN 0957-4174.  doi: https://doi.org/10.1016/j.eswa.2020.114189.  URL http://www.sciencedirect.com/science/article/pii/S0957417420309234.
  * Rizvi et al. [2018] Zarghaam Haider Rizvi, Frank Wuttke, and Amir Shorian Sattari.  Dynamic analysis by lattice element method simulation.  In Wei Wu and Hai-Sui Yu, editors, _Proceedings of China-Europe Conference on Geotechnical Engineering_ , pages 405–409, Cham, 2018. Springer International Publishing.
  * Rizvi et al. [2019] Zarghaam Haider Rizvi, Mijo Nikolić, and Frank Wuttke.  Lattice element method for simulations of failure in bio-cemented sands.  _Granular Matter_ , 21, 2019.  doi: https://doi.org/10.1007/s10035-019-0878-6.
  * Rizvi et al. [2020] Zarghaam Haider Rizvi, Syed Husain Mustafa, Amir Shorian Sattari, Shahbaz Ahmad, Peter Furtner, and Frank Wuttke.  Dynamic lattice element modelling of cemented geomaterials.  In Amit Prashant, Ajanta Sachan, and Chandrakant S. Desai, editors, _Advances in Computer Methods and Geomechanics_ , pages 655–665, Singapore, 2020. Springer Singapore.
  * Ruder [2016] Sebastian Ruder.  An overview of gradient descent optimization algorithms, 2016.
  * Rumelhart et al. [1986] David E Rumelhart, Geoffrey E Hinton, and Ronald J Williams.  Learning representations by back-propagating errors.  _nature_ , 323(6088):533–536, 1986.
  * Sattari et al. [2017] A.S. Sattari, Z.H. Rizvi, H.B. Motra, and F. Wuttke.  Meso-scale modeling of heat transport in a heterogeneous cemented geomaterial by lattice element method.  _Granular Matter_ , 19(66), 2017.
  * Sha and Edwards [2007] W. Sha and K.L. Edwards.  The use of artificial neural networks in materials science based research.  _Materials & Design_, 28(6):1747 – 1752, 2007\.  ISSN 0261-3069.  doi: https://doi.org/10.1016/j.matdes.2007.02.009.  URL http://www.sciencedirect.com/science/article/pii/S0261306907000520.
  * Springenberg et al. [2015] J Springenberg, Alexey Dosovitskiy, Thomas Brox, and M Riedmiller.  Striving for simplicity: The all convolutional net.  In _ICLR (workshop track)_ , 2015.
  * van den Oord et al. [2016] Aaron van den Oord, Sander Dieleman, Heiga Zen, Karen Simonyan, Oriol Vinyals, Alex Graves, Nal Kalchbrenner, Andrew Senior, and Koray Kavukcuoglu.  Wavenet: A generative model for raw audio, 2016.
  * Wong et al. [2014] John Kam-wing Wong, Kenichi Soga, Xiaomin Xu, and Jean-Yves Delenne.  Modelling fracturing process of geomaterial using lattice element method.  01 2014.
  * Wuttke et al. [2012a] F Wuttke, M Asslan, and T Schanz.  Time-Lapse Monitoring of Fabric Changes in Granular Materials by Coda Wave Interferometry.  _Geotechnical Testing Journal_ , 35(2):353–362, feb 2012a.  ISSN 0149-6115.  doi: 10.1520/GTJ103546.  URL https://www.astm.org/DIGITAL_LIBRARY/JOURNALS/GEOTECH/PAGES/GTJ103546.htm.
  * Wuttke et al. [2012b] F Wuttke, K Markwardt, and T Schanz.  Dispersion Analysis in Geotechnical Laboratory Tests: Time-frequency and Time-scale Signal Transforms.  _Geotechnical Testing Journal_ , 35(5):703–714, aug 2012b.  ISSN 0149-6115.  doi: 10.1520/GTJ103727.  URL https://www.astm.org/DIGITAL_LIBRARY/JOURNALS/GEOTECH/PAGES/GTJ103727.htm.

##  Appendix A Example of Predictions

The following figures show true crack occurrence and predictions made by our
best performing model for 320 valid cases in 16 columns and 20 rows. In Figure
19, the 5 different case types (see Section 2.3 “The Damage Detection
Dataset”) are marked with different coloured letters. Figure 20 is the true
occurrence of crack in full resolution (100 x 100) and Figure 21 gives the
true occurrence in a reduced resolution (16 x 16). The models in this work are
trained based on the labels of reduced resolution. The output in Figure 22 \-
25 are made by two models which are trained with focal loss with $\alpha=0.9$
$\gamma=0.4$ and $\alpha=0.35$ $\gamma=0.2$, the two recommended alpha and
gamma values we found in our experiments. Figure 22 and 24 show the predicted
probability of crack existence for each pixel, a brighter pixel indicates a
higher probability of crack inside. While in Figure 23 and 25, binary
predictions with the threshold that the pixel with probability greater than
0.5 is considered as having a crack inside.

Figure 19: The category for 320 test cases. The test cases are categorised
into four types: 1). randomly generated samples with randomly generated cracks
(Type-N), 2). randomly generated samples with no crack (Type-R), 3). randomly
generated with similar cracks (Type-S), and 4). the same sample with different
cracks (Type-C). They are marked by the colored marks. The special cases
(marked as “SC” in yellow) are the 7 cases we intentionally generated with the
same cracks in training data but from different samples. Figure 20: The true
crack occurrence in reduced resolution (100 x 100) for 320 testing cases.
Figure 21: The true crack occurrence in reduced resolution (16 x 16) for 320
testing cases. Figure 22: The predicted probability of crack existence in
pixels for 320 testing cases, the brighter a pixel’s color is indicates the
higher probability of crack existence inside the pixel. It is made by the
model trained with $\alpha=0.9$ and $\gamma=0.4$ Figure 23: The binary
predictions made by the model trained with $\alpha=0.9$ and $\gamma=0.4$,
where the pixel with a probability greater than a certain threshold (0.5) is
considered as having a crack inside. Figure 24: The predicted probability of
crack existence in pixels for 320 testing cases, the brighter a pixel’s color
is indicates the higher probability of crack existence inside the pixel. It is
made by the model trained with $\alpha=0.35$ and $\gamma=0.2$ Figure 25: The
binary predictions made by the model trained with $\alpha=0.35$ and
$\gamma=0.2$, where the pixel with probability greater than a certain
threshold (0.5) is considered as having a crack inside.