Wave propagation in random media, parameter estimation and damage detection via stochastic Fourier integral operators

This paper presents a new approach to modelling wave propagation in random, linearly elastic materials, namely by means of Fourier integral operators (FIOs). The FIO representation of the solution to the equations of motion can be used to identify the elastic parameters of the underlying media, as well as their statistical hyperparameters in the randomly perturbed case. A stochastic version of the FIO representation can be used for damage detection. Hypothesis tests are proposed and validated, which are capable of distinguishing between an undamaged and a damaged material, even in the presence of random material parameters. The paper presents both the theoretical fundamentals as well as a numerical experiment, in which the applicability of the proposed method is demonstrated.


Introduction
Modelling wave propagation in elastic solids goes back to the first half of the 19th century. According to the short historical survey in [1], the work of Raleigh, Lamb and Love in the first half of the 20th century sparked a huge development of the subject first in seismology [2] and then in the material sciences [3] and geotechnics [4]. In the first case, elastic waves, often sound waves, are used for detecting the internal structure of the earth crust; in the second case, ultrasonic measurements are used for damage detection, as well as mechanically excited waves. Under the heading of vibration based condition monitoring [5], a wide range of methods has been developed, from modal and spectral analysis to finite element based dynamical models and more. The reader is referred e.g. to the monographs [6,7,8,9] and the survey papers [5,10,11,12,13,14,15,16].
In the past two decades, additional attention has been given to the random structure of soils and materials and the modelling of wave propagation in randomly layered or stochastically perturbed media; a sample of references from seismology, geotechnics, and material sciences is [17,18,19,20,21]. This paper proposes a new approach to modelling wave propagation in random elastic materials, based on Fourier integral operators (FIOs). This type of operators has been introduced by [22] in the 1970s and is now widely used for representing solutions to hyperbolic partial differential equations. The operators are of the form where Φ is the phase function, which encodes the geometry of wave propagation, and a is the amplitude function, representing the frequency content. Typically, f (x) is the source and u(x) the response of the medium; time dependence is not shown explicitly in formula (1). A further novelty of the presented approach is that randomness of the material will be modelled through the stochastic phase functions and amplitudes of the representing Fourier integral operators. This is in contrast to the common approach in which randomness in the material is incorporated by taking the coefficients of the equations of motion as random fields [23,24]. The proposed alternative approach has the advantage that the stochastic system response, described through (1), can be analyzed and simulated in a rather explicit form. Fourier integral operators will be used for (a) propagating acoustic waves through the material and computing the dynamic response to localized excitations; (b) identifying the elastic parameters of the material (Young's modulus, Poisson's ratio) as well as their statistical hyperparameters (variances and correlation lengths); (c) designing statistical tests detecting significant changes in the parameters (due to damage or fatigue).
In short, the approach proceeds as follows. The general assumption is that the material is linearly elastic with inherent randomness. The spatial randomness of the material parameters is caused by small random perturbations of otherwise constant nominal values. The first step is to set up the Fourier integral operators solving the equations of motion with constant material parameters. These parameters appear explicitly in the phase functions and amplitudes of the operators and -in the second step -are replaced by spatial random fields there. The material parameters and their statistical hyperparameters are then estimated by comparing the FIO solution to measured data in a finite number of locations (sensors). This results in a stochastic Fourier integral operator model of the propagation of waves in the given elastic material, interpreted as the undamaged state. In a third step, this model can be used to generate a large number of system responses of the undamaged state by Monte Carlo simulation. This sample can be interpreted as a realization of the null hypothesis of undamaged material, against which later measurements of the possibly damaged material are tested.
Statistical methods for damage detection have been proposed in the literature before [25,26], also in combination with Monte Carlo sampling [27]. Compared to other numerical methods (finite differences, finite elements, modal analysis) the advantage of the FIO based approach lies in the following. Usually, in parameter identification and damage detection one has only a few sensors at certain locations of the material. These sensors measure the dynamic system response, which is a time-dependent signal. Thus one does not have complete spatial information at hand. This setup motivates Fourier integral operators, since Fourier integral operators can simulate the response to an excitation of an elastic material at single points without computing the solution on the whole domain. This makes parameter estimation very efficient. Only setting up and calibrating the baseline model of undamaged material requires an optimization procedure and a larger amount of model evaluations. Generating Monte Carlo samples of the baseline response and comparing it with (features of) measured data is computationally inexpensive.
This paper serves to set up the theoretical concepts and to demonstrate the applicability in the case of isotropic, linearly elastic materials. In this case, the Helmholtz decomposition allows one to decouple the equations of motion into four scalar wave equations for the wave potentials, which in turn are conveniently solved by Fourier integral operators. However, the approach works for other materials (e.g. orthotropic materials) as well. The hypothesis tests will be validated by artificially generated measured data, computed by a finite element model.
The computational advantages of the FIO method come with certain limitations on the range of applicability. In fact, FIOs of the form (1) can model forward propagation of stochastically perturbed waves, but the representation does not take into account internal reflections due to local random changes of the material properties. Thus local scattering or mode conversions are not modelled. Consequently, the presented perturbative approach is valid for random media in the weakly heterogeneous scaling regime only. (It is quite common in the theory of wave propagation in random media that different regimes are analyzed by means of different models, see e.g. [19,21,28].) The structure of the paper is as follows. The second section serves to present the theoretical background, the Helmholtz decomposition and the solution of scalar wave equations by means of Fourier integral operators in three space dimensions. The concept of stochastic phase and amplitude functions will be briefly introduced. Detailed formulas and explicit expressions for the representing Fourier integral operators will be given in the case of a plane strain problem.
In the third section, it is assumed that time-dependent measured data at certain locations are available, and it is shown how to estimate the material parameters and their statistical hyperparameters and how to calibrate the baseline stochastic FIO-model.
The fourth section briefly describes how statistical hypothesis tests can be designed with the help of the baseline stochastic FIO-model.
The last section sets out to validate the proposed procedure in a numerical example. It is assumed that a block of aluminum is excited by a harmonic line force. The laboratory measurements are represented by a finite element simulation with stochastically perturbed Lamé parameters. It is shown that the FIOsolution is able to reproduce the results of the FE-simulation for unperturbed constant Lamé parameters accurately. Furthermore, if the material parameters in the FE-simulation are randomly perturbed, it is demonstrated how the nominal values, standard deviations and correlation lengths of the Lamé parameters can be estimated. This results in the calibrated stochastic FIO-model of the undamaged material. The stochastic FIO-model of the undamaged state is used to generate a Monte Carlo sample of possible system responses of the undamaged material. Four scenarios of material states are considered: undamaged material, stiff material, soft material, or presence of a crack. The latter three are considered as degradation, fatigue, or damage. Three hypothesis tests are applied, and it is shown that they can distinguish between all four scenarios. Finally, the performance of the tests is evaluated with respect to false classification rates. Repeating the tests with different samples of the underlying random fields, it is shown that the tests behave as designed.
A rough quantification of the range of applicability, following the classification of scaling regimes from [19], concludes this section.
The paper ends with the conclusion, followed by an appendix containing additional graphical and tabular illustrations of the results referred to in the text.
The paper is based on the work [29] of the second author. The mathematical background of stochastic Fourier integral operators has been published in [30]. The paper is a continuation of [31], where the propagation of ultrasonic waves in one-dimensional media was addressed. The authors would like to point out that another aspect of Fourier integral operators combined with stochastic processes, even for nonlinear equations, has been worked out in a series of papers [32,33,34]. However, in these papers the phase functions are taken deterministic, the stochastic processes appear in the driving forces.

Helmholtz decomposition
This subsection serves to motivate the occurrence of FIOs in the solution of the equations of motion of linear elasticity. The displacement u of a linear elastic, homogeneous, isotropic body is governed by the equations of motion where f is the body force, λ and µ the Lamé constants and ρ the density; ∆ denotes the Laplace operator. The classical Helmholtz decomposition allows one to decompose this system of equations into four decoupled wave equations of the form where c is either the lateral wave speed (c 2 = (λ + 2µ)/ρ) or transversal wave speed (c 2 = µ/ρ), and u denotes the pressure potential or one of the components of the shear potential. Details about the decomposition can be found, e.g., in [1]; the case of two space dimensions will be elaborated in Subsection 2.3. For simplicity, assume that F can be written as F (x, t) = g(x)h(t). Taking the Fourier transform with respect to x, yields the ordinary differential equation Here x, ξ denotes the Euclidean scalar product and ξ the Euclidean norm.
Applying the inverse Fourier transform yields the solution operator, represented by means of oscillatory integrals These integrals can be seen as time-integrated Fourier integral operators. For example, the first integral contains 2ic ξ dy dξ (5) which is of the form (1) with Φ(x, y, ξ) = x − y, ξ + c ξ (t − s) and with a(x, y, ξ) = (2ic ξ ) −1 . In general, this integral does not converge as it stands, but has to be interpreted as an oscillatory integral, as e.g. described in [35].
(Actually, an additional regularization at ξ = 0 is required to fully conform with the theory of Fourier integral operators.) However, if g is sufficiently regular so that g is absolutely integrable then (5) can be directly evaluated by iterated integration.

Stochastic perturbations
In the derivation of the solution operator in the previous subsection, the material parameters were assumed to be constant. Generally, even in an undamaged material, the material parameters show spatially random deviations. (See e.g. the ultrasonic measurements of composite plates in [31].) As outlined in the introduction, the random structure of the material will be accounted for by adding random perturbations of the phase function and the amplitude function in the solution operator. These random perturbations take the form of spatial random fields, to be calibrated by measurements. This has the advantage that evaluating the solution by means of stochastically perturbed Fourier integral operators is much simpler than modelling the coefficients in system (2) and then computing the stochastic solution by standard numerical methods.
This representation has a physical interpretation, namely r 1 describes the perturbation geometry of the propagation, while r 2 is the perturbation amplitude of the wave number.
If the input g is not regular enough, for example, a rectangular pulse or even a Dirac delta function, one needs to ensure that the perturbed integral still defines an oscillatory integral. Thus, one has to make sure that r 1 and r 2 are of the right class, i.e., Φ + r 1 must lie in an appropriate space of phase functions (see [30]) and a + r 2 must belong to the Hörmander symbol class S m ,δ with 0 < ≤ 1, 0 ≤ δ < 1 (see e.g. [35]). This difficulty will be avoided in the sequel by assuming that g(y) in (6) is sufficiently smooth and has a sufficiently rapid decay rate at infinity, so that the Fourier integral operators are given by converging iterated integrals, which can be directly numerically evaluated.

A plane strain model problem
For the sake of simplicity in presentation and computation, the methods of the previous subsections will be elaborated in detail in the case of two space dimensions. Consider a linearly elastic, homogeneous, isotropic body occupying the whole three-dimensional space. Assume that the force term f in (2) is constant along the z-axis, i.e., depends only on x, y and t. Since the body is assumed to be at rest at t = 0, the solution u to (2) is constant along the z-axis as well, and it is of the form u(x, y, z, t) = [u 1 (x, y, t), u 2 (x, y, t), 0] T .
Without loss of generality, one may assume that the force is applied only in y-direction. More formally, it is assumed that the force is of form where g is the support of the force, and h is the time-dependent force term. The problem is reduced to two dimensions, where the displacements u 1 (x, y, t), u 2 (x, y, t) satisfy the equations In this case, the Helmholtz decomposition is performed as follows, see e.g. [1]. There is a pressure potential V and a shear potential W such that The potentials satisfy the wave equations with c 2 l = (λ + 2µ)/ρ and c 2 s = µ/ρ, where the functions v and w are to be determined from Making the ansatz v(x, y) = ∂ y ϕ(x, y) and w(x, y) = −∂ x ϕ(x, y) for some function ϕ(x, y), one has that Thus, by solving the two-dimensional Poisson equation with the growth condition lim (x,y) →∞ ϕ(x, y) = 0, one uniquely obtains ϕ and thus v and w. In more detail, if g(x, y) is the Dirac delta function δ(x, y), then ϕ is the fundamental solution of the two-dimensional Laplace operator, For general compactly supported g(x, y), the solution is obtained by convolution with the fundamental solution.
To continue solving (7), the spatial Fourier transform is applied in equation (9), similarly to Subsection 2.1. This leads to the ordinary differential equations Solving the ordinary differential equations and applying the inverse Fourier transform one ends up with the following solution operators: and Instead of inverting the Laplace operator as in (10), one can simply assert that Together with the relation from (8), the final solution is obtained as where and These integrals are convergent, since there are no poles left anymore.
Before describing the shape of the stochastic perturbations of the phase functions, and for later reference, a switch of notation from the Lamé parameters λ and µ to the more commonly used elastic material parameters, Young's modulus E and Poisson's ratio ν, is convenient. The parameters are related by respectively. The pressure and shear wave speeds are expressed as respectively. The stochastic perturbations of the phase functions are obtained by making c l and c s space dependent and setting that is, the constant parameters E and ν from formula (16) are replaced by random fields E(x, y) and ν(x, y) of the form Here R E and R ν are homogeneous, mean zero random fields, whose statistical parameters will later have to be calibrated by measurements. Note that the density ρ enters in c l and c s only through the quotient E/ρ. Thus, for the purpose of this analysis, it is no restriction of generality to take ρ constant, putting all randomness formally into E and ν. Strictly speaking, the insertion of non-constant parameters in (17) violates the derivation of the underlying equations in (7); otherwise occurring spatial derivatives of the parameters are neglected. However, this approximation is in line with the general limitation of the method to small random perturbations, as detailed in Subsection 5.5.
Remark. In reality, elastic bodies occupy a finite region. Then system (2) will have to be complemented by boundary conditions. However, in the intended applications, the force is localized in a small spatial region. Thus, for short time, wave propagation can be modelled by the full-space system, as long as the waves do not reach the boundary. (For a discussion of distances and travelling times in the numerical example of Section 5 see Subsection 5.5.)

Stochastic parameter estimation
The random fields (18) associated with the modulus of elasticity and Poisson's ratio are determined by their type of distribution and their statistical parameters. This section presents the proposed procedure how to identify these parameters using the Fourier integral operator representation of the solution.
It will be assumed that the random fields R E and R ν are independent, homogeneous, centered and Gaussian with autocovariance functions where Here σ E , σ ν and L E , L ν denote the respective standard deviations and correlation lengths. However, these assumptions are not obligatory. Different types of autocovariance functions can easily be implemented, and the same methods are also applicable for non-Gaussian random fields. The randomness of the material is thus described by the six parameters E mean , σ E , L E , ν mean , σ ν , L ν . The model scenario is that the time-dependent system output is recorded at n sensor locations. In practice, this will be a record of measured data. For the present study, the record is artificially generated by computer simulation. In any case, it is assumed that the displacements u 1 and u 2 are known at the locations (x j , y j ), j = 1, . . . , n for all times t in the time interval under consideration. The task is to estimate the six parameters by comparing the data with the solution obtained by the Fourier integral operator representation (13).

Estimating the nominal values
The first task is to estimate the nominal values E mean and ν mean based on the data given in the sensor locations. This will be done by fitting FIO solutions (13) with constant parameters E 0 , ν 0 in (16).
To shorten notation, let g j (t) = u meas (x j , y j , t) be the measured signal in (x j , y j ), recorded over the time interval [0, Goodness of fit is measured in terms of the weighted norm and E 0 , ν 0 are estimated as solutions to the minimization problem The weights w 1 and w 2 are chosen as the maximum of the measured displacements in x-direction and in y-direction, respectively, at each sensor location and over the whole time interval [0, T max ]. This ensures that both directions enter with equal importance in the error norm. If a sample of size N of the measured signals u meas (x j , y j , t) is available, one can increase the robustness of the algorithm by repeating the calibration N times, resulting in the estimates E Numerical experiments showed that E * mean and ν * mean are good estimators of E mean and ν mean , respectively, even with small sample size N (see Subsection 5.2).

Estimating the standard deviations
The next task is to estimate the standard deviations σ E and σ ν . To obtain the required data, the nominal values are first estimated for each sensor location separately: The estimator of the variance of E 0,opt and ν 0,opt will be denoted by σ 2 E0,opt and is computed as follows: invoking the hypothesis of homogeneity so that the expectations E(E 0,opt,j ) = E mean may be assumed to be independent of j.
It turned out that the estimator σ E0,opt has a bias depending on the correlation length of the underlying random field E(x, y) in (18). In fact, σ E0,opt underestimates σ E for short correlation lengths, the underestimation getting more distinct for smaller correlation lengths.
As detailed in Subsection 5.5, the effect can be explained as follows. When the correlation length gets smaller and smaller, one is entering the quasi-homogeneous regime. The medium behaves like a deterministic one with almost constant E [28]. Therefore, the estimator σ E0,opt tends to zero as well. On the other hand, when the correlation length gets larger, the underestimation disappears because then the random field R E (x, y) tends to become a constant random variable R E with standard deviation σ E , for which σ E0,opt is an unbiased statistical estimator. The result is a nonlinear dependence of σ E0,opt on L E , as shown in Figure 3.
Numerical experiments showed that σ E0,opt is a linear function of σ E with a slope depending on L E , that is, It is proposed here that h(L E ) can be obtained as a calibration curve through numerical experiments as follows. First, one fixes a standard deviation σ E and models the dependence of the estimator (23) on the correlation length, that is one tries to find a functional relationship σ E0,opt = f(L E ). For this, one selects a set of correlation lengths (L E ) i , i = 1, . . . , m. Then, for each correlation length one produces a Monte Carlo sample of signals u meas (x j , y j , t) by means of FE-simulations, where the random field has the fixed standard deviation σ E and correlation length (L E ) i . For the i th sample, one estimates (σ 2 E0,opt ) i by the procedure described above, resulting in the values (σ E0, The curve f(L E ) is obtained through regression by means of a power function. The linear dependence of σ E0,opt on σ E allows one to compute the slope by h(L E ) = f(L E )/σ E . An explicit example of the calibration procedure will be given in Section 5.
This calibration curve can be used for any later experiment. Thus, if one has to determine the standard deviation of a new sample of the same material, one can estimate it by This estimation procedure, however, is only possible if one knows L E . The next subsection indicates how the correlation length L E can be estimated. Similarly, the variance of ν 0,opt is estimated as: Numerical experiments showed that σ ν0,opt is a sufficiently accurate estimator of σ ν . In any case, in the later application the coefficient of variation of ν was taken rather small (around 1%), so a possible bias of the estimator for σ ν was not further investigated, and was taken as estimate for the standard deviation of Poisson's ratio.

Estimating the correlation length
Since data at different sensor locations are available, one can also estimate the correlation between the sensors and from there the overall correlation length. It is assumed that the parameters in the sensor locations derive from a stationary random field obeying the autocovariance function is the Euclidean distance between sensors j 1 and j 2 . The sensors are paired into groups with the same Euclidean distance, denoting the l th group of sensor pairs with Euclidean distance r l by P l . Then the covariance at lag r l is empirically estimated by The correlation length L E0 is then estimated by fitting the autocovariance function (27) to the empirical covariances, resulting in the estimate Similar to the observed bias of the estimator σ E0,opt , the estimator L E0,opt overestimates L E when L E gets small. Numerical experiments suggested, though, that the estimator depends on the correlation length through a functional relationship L E0,opt = g(L E ). To determine this relationship, one can use the same simulation data as used to calibrate f(L E ). This time, for every i = 1, . . . , m, one computes (L E0,opt ) i and obtains the relation (L E0,opt ) i = g((L E ) i ). Then, by regression, one can fit the curve g(L E ) and for new samples one may obtain the estimator L * E by inverting g and setting Finally, L * E is inserted in place of L E in (24) to get the estimate σ * E of the standard deviation.

Testing for damage
This section serves to present a proposal for statistical tests for damage, indicating whether the properties of a given material have undergone some changes or not. The starting point is a calibrated FIO model of the undamaged, random material. The response of the undamaged material is modelled through equations (13), (14), (15) with wave speeds (17) given by the random fields (18), (19). The parameters of the random fields are calibrated by the procedure outlined in Section 3. The sensor locations remain fixed.
The strategy is as follows. Using the stochastic Fourier integral operator model of the undamaged material, a sample of signals in the sensor locations is generated. Thereby, percentiles of the distribution of certain features of the undamaged material are obtained. These data are compared with measured data of the possibly damaged material. The procedure can be cast in the form of a hypothesis test, in which the null hypothesis (undamaged material) is accepted or rejected, depending on whether the measured features remain within certain bounds or not. The generation of the Monte Carlo sample can be done in advance of the testing. Generating a large Monte Carlo sample of system responses in finitely many sensor locations through the FIO representation is computationally inexpensive.
Details will be presented in Section 5. Here is a summary of the testing strategy.
Step 1. A Monte Carlo sample of size N MC of the random fields (18), (19) is generated. The corresponding solutions (13) are computed in the sensor locations, using the Fourier integral operator representation.
Step 2. One selects features of each signal at each sensor location. This choice is critical, since this has a very strong influence on the efficiency of the test.
Step 1 delivers a Monte Carlo sample of size N MC of each feature at each sensor location.
Step 3. The null hypothesis (undamaged material) is cast in the form of critical percentile ranges for the features.
Step 4. Data of the specimen to be tested are acquired by measuring the signals of the sensors. Then, for each feature one computes the p-value with respect to the Monte Carlo sample.
Step 4 results in a tuple of p-values for each sensor location. An overall p-value can be aggregated by e.g. computing the overall minimum or mean. In this paper, the mean p-value for each sensor location was computed and the minimum of these values was taken as the overall p-value.
Step 6. The null hypotheses is rejected if the overall p-value is smaller than a given threshold, for example, 1% or 5%.
Remark. The aggregation of p-values is a much debated issue [36,37,38,39], and there are certainly more sophisticated approaches. However, the method chosen in Step 5 worked well. In fact, the performance of the test was validated in Subsection 5.4 where it was found to obey the designed acceptance and rejection rates.

Numerical example
In this section, an application of the concepts to a prototypical example, a plane strain problem, will be presented. The section starts with the description of the model and the material assumptions, followed by a presentation of the proposed parameter calibration method. Next, statistical tests for damage are put to work in four scenarios. The section ends with a validation and performance analysis of the devised tests.

The numerical model
The material structure is an aluminum block with the following nominal material parameters: Young's modulus is E = 70 GPa, Poisson's ratio is ν = 0.35, and the density is ρ = 2.70g/cm 3 . The cross section of the block was assumed to be a square of size 10 × 10 cm (henceforth referred to as the numerical domain).
The geometric configuration can be seen in Figure 1. There are n = 8 (virtual) sensors located around the centered force.  As in Subsection 2.3 it is assumed that the force term f is independent of z and has the special form where δ(x, y) is the Dirac delta function. The displacements u 1 (x, y), u 2 (x, y) at point (x, y) are then given by (13). When evaluating this formula, some numerical issues have to be addressed. Since the Fourier transform of δ(x, y) is δ(ξ, η) ≡ 1, the integrals in (13) are not convergent. To remedy this problem, the Dirac delta function was replaced by a regularized version, namely by where ε is a small regularization parameter (chosen here as ε = 0.1, i.e., about three times the grid length of the numerical domain). The Fourier transform is then explicitly given by δ ε (ξ, η) = exp − ε 2 (ξ 2 + η 2 ) . The numerical domain was discretized with 256×256 points. The considered time interval was 7 µs with ∆t = 0.05 µs. The integrals were computed by the rectangular quadrature formula for all sensor locations. For sufficient accuracy, higher order integration schemes were not required.
At ξ = η = 0 one faces the additional problem that the integrands in (14) and (15) are not defined. This does not matter analytically since the value at one point does not change the integral. However, it is a problem for numerical evaluations. So in order to compute the solution for this point one has to go back to equation (7). If one takes the Fourier transform and sets ξ = η = 0 then using the fact that δ ε (0, 0) = 1. Since the material is at rest at t = 0, it follows that u 1 (0, 0, t) = 0 and u 2 (0, 0, t) =

Numerical validation
In the deterministic case of constant material parameters, the solution procedure via the Fourier integral operator formula (13) was validated by comparison with a finite element simulation of the plane strain model, solved by (7), on the same domain with the same material parameters. The same force was applied, but in non-regularized form as a point force. For the simulation, a standard Abaqus routine was used, with the implicit Euler scheme for time integration.
The (time-dependent) solution in the eight locations (x j , y j ), j = 1, . . . , 8 was computed at all times t. The evaluation of the FE-model took approximately 800 seconds on a workstation with 5 CPU cores used.
The comparison of the FIO-solution and the FE-solution in x-and y-direction can be found in Figures 6 and 7 in the Appendix and showed satisfactory accuracy. The relative error is with respect to the maximum displacement occurring in the FE-solution in the respective directions. One should point out that there is a systematic error in the FIO-solution, since the point force was regularized. Therefore, the force term acts on an area and not at a point. That means that the FIO-signal is slightly blurred. In addition, this also means that the signal arrives at the sensor locations slightly earlier than the FE-signal. Of course, if one refines the mesh, the regularization of the point force can be more sharp and the error vanishes. The finite element model was also used to simulate artificial data, meant to represent laboratory measurements in the next subsections.

Parameter estimation
This subsection serves to demonstrate the calibration procedure of Section 3 in the numerical example. Artificial data were generated by the FE-model of subsection 5.1. Young's modulus and Poisson's ratio were entered as random fields of the form (18) with autocovariance functions (19). Input parameter values were E mean = 70 GPa, ν mean = 0.35, σ E = 3.5 GPa (5% of its nominal value), σ ν = 0.005, and L E = L ν = 3 cm. For the parameter calibration, the output of the FE-calculation was extracted at the eight sensor locations (Figure 1). This output served as measured data for the FIO-based calibration procedure.
Concerning the estimation of the nominal values E mean and ν mean , three scenarios were tested, according to Table 1. (II) Full random field model for the elastic parameters, no repetitions in formula (22), N = 1.
(III) Full random field model for the elastic parameters, means estimated by formula (22) with N = 10.
The minimization of (21) was done by means of the Nelder-Mead algorithm [40,41]. To show that the estimate works well even with a bad initial guess for the algorithm, E ini = 50 GPa and ν ini = 0.3 was chosen. The results of the optimization procedure as well as the computing time can be found in Table 1 In order to estimate the standard deviations, formulas (23) and (25) were used with N = 100 finite element simulations to compute σ E0,opt and σ ν0,opt . As indicated in (26), σ ν0,opt could be taken as estimator σ * ν of the standard deviation of Poisson's ratio. (The present simulation resulted in the estimate σ * ν = 0.0036.) In Subsection 3.2, the bias of σ E0,opt has been pointed out. As noted there, further numerical experiments had shown that it can be described by a relation of the form σ E0,opt = h(L E )σ E . For example, the linear dependence of σ E0,opt on σ E at fixed correlation length L E = 3 cm can be seen in Figure 2.   Figure 3. Finally, one computes h(L E ) = f(L E )/3.5. So far, if the correlation length is assumed to be known, one can obtain the estimate σ * E of the standard deviation of Young's modulus through formula (24).
It remains to estimate the correlation lengths. The procedure has been outlined in Subsection 3.3. To obtain the empirical covariances at different distances, the sensors were paired as shown in Table 2. Fitting the autocovariance   function (27) results in the estimator (28), which was observed to be biased in Subsection 3.3, obeying a functional relation L E0,opt = g(L E ).
Again, a power function of the form g(L E ) = γ 0 L γ1 E showed a satisfactory fit. The data and fitted function can be seen in Figure 4.
Finally, an estimate L * E of the correlation length L E can be achieved through formula (29). The results of the present numerical simulation are summarized in Table 3.  It can be observed that the obtained estimates of the correlation length and the standard deviation are satisfactory, but not as accurate as the estimates  Table 1. Larger sample sizes in the determination of the calibration curves could improve the estimates, but at the price of considerably increased computational cost. The tests for damage in Subsection 5.3 are based on the nominal value E mean of Young's modulus which can be accurately estimated even by a single measurement, with N = 1 in formula (22).

Damage detection
As outlined in Section 4, the testing strategy is based on a stochastic Fourier integral operator model of the undamaged material, which is used to generate a sample of signals in the sensor locations. This sample serves as a realization of the null hypothesis. Measured data of the possibly damaged material are tested against this sample. As before, the parameter values for the undamaged state were assumed to be E mean = 70 GPa, ν mean = 0.35, σ E = 3.5 GPa, σ ν = 0.005, and L E = L ν = 3 cm. These data were used to generate the random fields (18), (19) producing the wave speeds (17), which entered the Fourier integral operator representation (13) of the displacements. The coefficients of variation were so small that no problems with the square roots in (17) could arise, and further cut-offs were not needed.
In practice, estimated values would be entered in (18), (19) according to Section 3. The purpose of this subsection is to demonstrate how the test procedure works, so this step was skipped here. The finite element model was used to generate artificial responses of the material having undergone various kinds of damage, represented through changes in the nominal value E mean , or through geometric deterioration. The standard deviations and correlation lengths were kept fixed throughout.
The null hypothesis was that the material is in the undamaged state, i.e., it has the parameters described above. As initial step, a sample of size N = 10000 of system responses was generated, distributed according to the null hypothesis, by means of the FIO-solution operator with wave speeds (17).
In order to test the material one has to choose characteristic features of the signal. This choice is quite critical, since a bad choice of features can lead to a poor distinction between damaged and undamaged material.
Note that for the tests, a full record of the measured signal (duration 7 µs) is available. Using the information captured in the Fourier spectrum of the signal suggests itself, because the amplitude of a selected frequency relates to attenuation, while the phase is directly related to propagation speed. Experiments with various sets of frequencies of the Fourier spectrum were undertaken. For the present purpose, the amplitude and the phase angle of the two most dominant frequencies of the recorded signal turned out to work best.
As described in Subsection 5.1, there were eight sensor locations with signals in x-and y-direction. The four signals in x-direction at sensors 2, 4, 5, 7 were close to zero. (Due to the special excitation, no pressure waves arrive at sensors 4 and 5, and no shear waves arrive at sensors 2 and 7.) Thus 12 signals remained for the analysis. Of each signal 2 × 2 features were extracted: the phase angle and amplitude at the first two nonzero frequencies in the DFT-spectrum (here ω 1 = 2π/7 ≈ 0.9 MHz and ω 2 = 4π/7 ≈ 1.8 MHz). Thus a total number of 48 features were used for damage detection.
Three tests with the following decision procedures were implemented: (I) A left-sided test for the phase angle of the signal, detecting small phase angles (corresponding to late arrival of the signal). For this test, the phase angles need to be ordered: The phase angle space was determined by the mean phase angle of the FIO sample plus/minus 180 degrees.
For each sensor location, the p-value for each single feature (i.e. the phase angle of one signal for one frequency) was determined, and the average p-value (over all features) at each sensor location was calculated. The overall p-value was taken as the minimum over all sensor locations.
(II) A right-sided test for the phase angle of the signal, detecting large phase angles, using the same principle as in (I).
(III) A two-sided test for the amplitude of the dominant frequencies. The aggregated p-value (over features and sensor locations) was computed as in (I).
These tests were applied to the following scenarios: S1: The material is having the desired properties (undamaged state).
S2: The material is not having the desired properties: E mean is too small (60 MPa).
S3: The material is not having the desired properties: E mean is too large (80 MPa).
S4: The material is having the desired properties, but is suffering a crack. The crack is modelled by a small region in the FE-model having a very small Young's modulus. The modelled crack lies between the applied point force and Sensor 3. Figure 5 shows realizations of the random field E(x, y) in Scenarios S1-S4. A confidence level of 99% was adopted. This meant rejection of the null hypothesis for p-values below α = 1%. (In Subsection 5.4, α = 5% will be considered for comparison as well.) The confidence bounds of the features corresponding to the chosen α-levels were determined through the sample of size 10000 of simulated stochastic FIO-solutions of the undamaged material. For a typical random field realization the tests showed the following results. Test (I) S1: The computed p-value lies well above the chosen α-level, which means that undamaged material will be accepted as such.
S2: The p-value is well below 1%, the phase angles are to the left of the lower confidence bounds. The test detects that E mean is small. In fact, for small E mean , the wave speed is low, the signal arrives later than in the undamaged state, and this corresponds to small phase angles.
S3: The p-value is close to one. The left-sided test does not reject the null hypothesis. In fact, Young's modulus is high, the wave speed is large, the signal arrives early and the phase angles are actually large.
S4: The crack has no significant influence on the phase angles.

Test (II)
This shows the same result as Test (I) with the exception that Scenario S2 has a large p-value and Scenario S3 has a low p-value. This is due to the fact that this is a right-sided test.
Test (III) S1: The p-value of the amplitude lies above the chosen α-level, and therefore the null hypothesis is not rejected.
S2: The amplitude also changes with Young's modulus, but the test is not suitable for detecting this scenario.
S3: The amplitude also changes with the Young's modulus, but the test is not suitable for detecting this scenario.
S4: The crack influences the signal in Sensor 3, since it is in between Sensor 3 and the wave source. As a consequence, the amplitude changes as well.
The null hypothesis is rejected in Sensor 3. In the other sensors, the behavior is varying, but the null hypothesis is usually accepted in Sensor 6, where the signal is not influenced by the crack.
It is interesting to note that the different behavior of Test (III) in the different sensors allows one to draw conclusions about the location of the crack. The three tests are most expressive when considered together. It also should be noted that, when repeating the tests, errors of the first kind (rejection of true null hypothesis) and of the second kind (acceptance of false alternative) do occur. However, as will be seen in the next subsection, the error rate of the first kind is around or below the chosen α-level. Figures 8 to 11 show the comparison of the Monte Carlo sample of size 10000 (representing the null hypothesis) with the FE-computed signal, based on one typical realization of the random field of the corresponding scenario in a certain sensor location. The histograms depict the marginal distributions of the respective features generated by the Monte Carlo sample. Also shown are the two-sided 99% bounds of the amplitude and the one-sided 99% upper and lower bounds of the phase angle for the indicated frequency, as well as the value of the feature obtained from the FE-computation. The top figures always depict the time dependent displacement in the indicated direction at the indicated sensor. The gray band visualizes the collection of the solution curves from the Monte Carlo sample whose features lie in the 99% confidence intervals.
The asymmetries in the histograms of the amplitudes in Figures 9 to 11 are caused by a strong nonlinear correlation of amplitude and phase angle at frequency ω 1 , resulting in an inversely U-shaped scatterplot (not shown here).

Performance evaluation of the tests
In order to check whether the tests obey the designed rates of false negative and false positive classifications, tests were repeatedly applied to FE-simulated data with different realization of the random fields. More precisely, the tests were applied to a sample of 100 FE-simulations of each scenario. The numbers in Table 4 suggest that the tests perform as designed (at α-levels of 5% and 1%, respectively). The performance of the tests, based on the 100 FE-simulations, is visualized in Figure 12 by means of a scatter plot 1 and a violin plot 2 .  Table 4: Number of p-values below the α-level among 100 repetitions.
As a conclusion, one can say that the performance analysis shows that the proposed procedure is capable of detecting damaged material, and furthermore the three tests are capable of distinguishing between the scenarios.
The influence of L E and σ E As was seen from the parameter estimation, the correlation length of E plays a critical role in estimating the parameters. So the question arises if it also influences the p-values of the tests. To check this point, the tests from above were applied to an undamaged material (i.e. E mean = 70 GPa and σ E = 3.5 GPa) with different correlation lengths L E . For each correlation length, a sample of finite element simulations with sample size 100 was tested against the Monte Carlo sample of size 10000 from Subsection 5.3. The results can be seen in Figure 13. One may observe that the correlation length does play a role, but does not influence the results too critically. The number of rejections can be found in Table 5.
In contrast to the correlation length, the standard deviation severely influences the measured signals. The tests from above were again applied to an undamaged material (i.e. E mean = 70 GPa and L E = 3 cm) with different σ E . For each standard deviation value, a sample of finite element simulations with sample size 100 was tested against the Monte Carlo sample from Subsection 5.3. The results can be seen in Figure 14. Up to a coefficient of variation of E of 5% all p-values were above the chosen α-level, at least in Test (I). If the coefficient of variation was larger than 15%, the null hypothesis was falsely rejected for almost the whole sample of size 100. However, since a large standard deviation can be interpreted as damaged material, too, this can be seen as a desirable property of the tests. The number of rejections can be found in Table 6.

Advantages and limitations
Advantages of the method The advantage of the described method is its low computational cost. The evaluation of 10000 FIOs took approximately 11 hours (= 10000 · 4 sec). If one wanted to compute a Monte Carlo sample of responses of the undamaged material of the same size by the finite element method, one could expect a computation time of 93 days (10000 · 800 sec).
The computation of the p-values is so fast that it can be done online. The computation of the features of the Monte Carlo sample of size 10000 representing the null hypothesis can be done in advance. So if one gets a measured signal one just has to compute its features and compare them with the Monte Carlo sample. This can be done in real time.

Limitations and range of applicability
The following limitation of the presented procedure should be pointed out. If the stochastic FIO representation has the form as described in Section 2, it does not incorporate any internal reflections (scattering, mode conversion) of the waves due to local changes of material parameters, as was already noted in the introduction and in Subsection 3.2. The experiments subsumed in Table 6 suggest that the applicability of the method is mostly influenced by the coefficient of variation of the material parameters, which should be sufficiently small.
In order to give a rough quantification, also with respect to frequency, wavelength, correlation length and propagation distance, one may invoke the classification of scaling regimes in [19, Section 5.1] 3 ; see also [28]. Denoting the typical frequency by f 0 , the important parameters for the occurrence or nonoccurrence of scattering are: the angular frequency ω 0 = 2πf 0 ; the wave number k 0 = ω 0 /c 0 , where c 0 is the typical propagation speed, the coefficient of variation σ of the underlying random fields, the scale of heterogeneities (here taken as the correlation length), and the typical propagation distance L 0 . For the present purpose, the following regimes listed in [19] are of interest: The effective medium: k 0 1, k 0 L 0 ∼ 1, σ 1 or σ ∼ 1. The medium is quasi-homogeneous, random scattering is weak, and no backscattering occurs.
The weakly homogeneous regime: k 0 ∼ 1, k 0 L 0 1, σ 1. The coupling between the wave and the medium is weak due to small σ. Significant scattering can only occur at large propagation distances.
Weakly homogeneous regime, subcase σ 2 k 0 L 0 ∼ 1: In this case, mode coupling and backscattering are of order one.
This puts the setup into the weakly homogeneous regime; at the considered propagation distance, at most weak scattering is expected. With regard to the applicability of the hypothesis tests, Table 5 shows the results of experiments with = L E between 0.05 cm and 20 cm. The smaller , the smaller k 0 and the closer one gets to the quasi-homogeneous regime. This confirms the observation that in the considered frequency range, the correlation length has little influence on the hypothesis tests.
On the contrary, the coefficient of variation σ does have an influence. In fact, for σ ≈ 0.1, the subcase σ 2 k 0 L 0 ∼ 1 is reached, where scattering of order one is to be expected. This in turn is confirmed by Table 6, which shows failure of the hypothesis tests beyond this value.
Summarizing, it is expected that the FIO approximation is of acceptable accuracy as long as one remains in the quasi-homogeneous or the weakly homogeneous regime and the coefficient of variation is significantly smaller than one.
This holds unrestrictedly for the hypothesis tests. For parameter estimation, calibration curves as described in Section 3 and Subsection 5.2 are required. Finally, the dimensions of the domain beyond the sensors must be chosen large enough so that the signal does not reach the domain's boundary in the measuring period (see the remark at the end of Subsection 2.3).

Conclusion
This article presented a Fourier integral operators based approach to modelling wave propagation in random linearly elastic materials. This FIO representation of the solution to the equations of motion is numerically fast, if one is only interested in the time-dependent solution in certain locations. For example, this is the case if one needs to compare the solution with sensor data in certain locations. This efficient simulation procedure can be used for identification of the material parameters. Even if the material parameters are randomly perturbed by a Gaussian random field, the estimation of the nominal value is stable and accurate. Furthermore, using a calibration curve one can also estimate the stochastic parameters of the random fields, such as standard deviation or correlation length.
The obtained material parameters can be used for constructing a stochastic FIO, describing wave propagation in the material with its given random properties. Since the evaluation of the FIO is fast, large Monte Carlo samples of the stochastic FIO solution can be generated with relatively small effort. This sample can be used to design hypothesis tests, to determine if the material has the desired properties or not (or if the properties have changed over time). In a numerical example, four scenarios were considered: the material has the assumed properties; the material is stiff; the material is soft; a crack is present. For the decision procedures of the designed tests certain features of the time-dependent solution were considered: the phase angles and the amplitudes of the most dominant frequency of the signals. It was shown that the three implemented tests are capable of distinguishing between the four scenarios. Furthermore, it was also shown by additional simulations that the implemented tests behave as designed with regard to the error rates of false classifications.