A Bayesian Method for Material Identification of Composite Plates via Dispersion Curves

Ultrasonic guided waves offer a convenient and practical approach to structural health monitoring and non-destructive evaluation. A key property of guided waves is the fully defined relationship between central frequency and propagation characteristics (phase velocity, group velocity and wavenumber)—which is described using dispersion curves. For many guided wave-based strategies, accurate dispersion curve information is invaluable, such as group velocity for localisation. From experimental observations of dispersion curves, a system identification procedure can be used to determine the governing material properties. As well as returning an estimated value, it is useful to determine the distribution of these properties based on measured data. A method of simulating samples from these distributions is to use the iterative Markov-Chain Monte Carlo (MCMC) procedure, which allows for freedom in the shape of the posterior. In this work, a scanning-laser Doppler vibrometer is used to record the propagation of Lamb waves in a unidirectional-glass-fibre composite plate, and dispersion curve data for various propagation angles are extracted. Using these measured dispersion curve data, the MCMC sampling procedure is performed to provide a Bayesian approach to determining the dispersion curve information for an arbitrary plate. The distribution of the material properties at each angle is discussed, including the inferred confidence in the predicted parameters. The percentage errors of the estimated values for the parameters were 10–15 points larger when using the most likely estimates, as opposed to calculating from the posterior distributions, highlighting the advantages of using a probabilistic approach.


Introduction
This paper focusses on determining dispersive characteristics of Lamb waves in arbitrary plates, with an emphasis on their use in non-destructive evaluation (NDE) and structural health monitoring (SHM). The use of ultrasonic guided waves (UGWs) for SHM strategies [1] can offer a number of distinct advantages, such as range and sizing potential, greater sensitivity and cost effectiveness. There are, commonly, three types of high-frequency stress waves which fall under the category of guided waves: Rayleigh waves, Lamb waves and shear horizontal waves. The former of these types propagate on a surface, whereas the latter propagate in 'thin' plates. Full descriptions and derivations of Rayleigh and Lamb waves are in [2][3][4], although a short introduction to some key concepts will be given here. A particular distinction of Lamb waves is their separation into symmetric As well as an estimation of the most likely values, it is also useful to determine the posterior distribution of these parameters. Some advantages of estimating these distributions include accounting for environmental conditions and for uncertainty propagation. A laudable example of this has been shown, using a genetic algorithm and extending the list of parameters to include a noise term [26]; this generated feasible elastic constants and a distribution based on an assumed Gaussian posterior. However, this assumption of the posterior shape is a shortcoming of the approach, as well as the absence of any possible inference on the cross-correlation between material properties. In addition, the genetic algorithm has a high computational cost [27].
An alternative Bayesian approach to this problem is to simulate samples from the posterior distribution using a Markov-Chain Monte Carlo (MCMC) procedure. This approach allows for an observation of the true shape of the posterior, as well as to simulate the multivariate distributional behaviour of the parameters. Previously, using MCMC would have been impractical, as each iteration requires solving the dispersion curves given the current estimate of the material properties, which has a high computational expense. However, with the faster computation of dispersion-curve solutions (as discussed above), it is now much more feasible to apply such a procedure. In this work, the Legendre polynomial expansion approach [18] has been used, though many other options are available. A key argument for this choice is because, using some numerical manipulations, it was possible to increase the speed of this calculation. Adding further, the Legendre polynomial approach directly use the elastic constants, and so inference on these is more direct, and this allows for more efficient sampling in the MCMC process. The primary increase in computational efficiency is a result of the solution form being that of an eigenvalue problem in which the eigenvalues are the negative of the phase velocity squared; which then allows calculation of only the first modes of interest by using the power iteration method.
In this work, propagation of Lamb waves in a glass-fibre-reinforced-polymer (GFRP) plate are measured and dispersion curve data are returned at various propagation angles. These dispersion curve data are then fed into the MCMC procedure and the posterior distributions of the material properties are analysed. In this initial work, the dispersioncurve solution method used does not include material damping. As the experimental observations of the dispersion-curves are normalised, material damping is not likely to affect simulation over the non-complex elastic constants which are used in the model. The limitations of this is that modelling of attenuation affects is restricted to not including viscoelastic effects. For brevity, the work in this paper considers only the real wavenumber modelling and observations. An important challenge, which this work aimed to address, was to develop a method of material parameter distribution estimation which had freedom in the posterior shape, and was able to estimate the cross-covariance parameters. Overcoming this challenge represents the key novelty of the paper, and was done so with the use of a MCMC sampling procedure. This method allows one to quantify the confidence in the estimates, as well as provide information on the elastic-constant combinations and how these influence each other, with respect to elastic waves. Another challenge with this work was the requirement for an efficient methodology setup, and to ensure it is readily expandable to allow use of more complex material layups and types. This issue was overcome by the embedding of the Legendre polynomial-expansion approach to dispersion-curve solutions.
The next section of this paper begins with explanation of how to obtain dispersion curve solutions, including: numerical solutions, experimental observations, and the Legendre polynomial expansion (LPE) approach for orthotropic materials. Using the LPE approach, a brief sensitivity investigation is included to discuss the effects of each material parameter on the dispersion curve solutions. Section 3 details the experimental method for returning observations of the dispersion curves for the plate. The remainder of Section 3 then details how to estimate the elastic constants given observations, and how this is extended to use the MCMC approach to simulate sampling from the posterior. The paper then finishes by presenting and discussing the results of the procedure, along with a discussion of suggested future work that the authors intend to pursue.

Lamb Wave Propagation in Plates
In order to better apply guided waves for SHM and NDE strategies, prior knowledge of their behaviour is essential. This section aims to introduce the physics of guided waves, the concept of dispersion curves, and important characteristics which are prevalent in this work.

Physics of Lamb Waves
Elastic waves in orthotropic, inhomogeneous media are described by the elastodynamic equation [28], where S is the four-index stiffness tensor, ρ is the material density, u is the displacement field for which the double dot represents double differentiation with respect to time. In bounded media, these waves will show as Lamb waves, which in isotropic elastic media will exhibit two distinct modes: symmetric and antisymmetric. For anisotropic or composite media, there also exists shear-horizontal modes as a solution to (1). When modelling guided waves in isotropic materials, the solutions to the two fundamental equations, derived from (1), given the relationship between frequency ω and wavenumber k. The known frequency and wavenumber can then be used to determine the phase and group velocity, c p and c g respectively, using, As the velocity of the wave is a function of the frequency, the waves are dispersive and plots of the relationship between frequency and wavenumber/velocity are called dispersion curves. For anisotropic materials, the calculation of the group velocity must also account for angular variations of the phase velocity, where k = {k 1 , k 2 , k 3 } is the 3-dimensional wavenumber vector [29], in Cartesian coordinates, the tensor form of the calculation is, where c gx and c gy are the group velocities in the x and y directions respectively. For the rest of this paper, as wave propagation is considered in one direction for an orthotropic material, the partial derivative with respect to the angle vanishes, therefore, the group velocity can estimated using (2).

Dispersion Curve Solutions for Orthotropic Media
Solutions of the dispersion curves for anisotropic media are more demanding. Here, the work of Lefebvre [18] is followed to formulate a computationally efficient method of solving for these curves. This method has been validated in the works of Othmani [30,31], as well as its improved computational efficiency demonstrated. This method uses a Legendre polynomial expansion to form an eigenvalue problem, utilising the orthonormal basis set for expansion of the field quantities. For orthotropic materials, the generalised Hooke's law can be rewritten as, where σ ij is the stress, ε ij is the strain, and C ij is the elastic constants. The elastic constant tensor C is the inverse of the stiffness matrix S, and the elements are defined using Voigt notation. The relationship between the strain and displacement can be expressed as, where x k is coordinates in direction k = 1, 2, 3. The boundary conditions of zero stresses on the surface can be applied by introducing a rectangular window function π h (x 3 ), the above-mentioned boundary conditions are automatically incorporated in the constitutive relations, and by substituting in the relationship for the strain (6), and transforming the spatial coordinates into dimensionless form q α , the constitutive relations are then, For a wave propagating in the x 1 direction, the displacement components are assumed to be of the form where U i (q 3 ) represent the magnitudes of the fields in the x i direction, and the non-italic i is the imaginary unit. Substituting (9) and (10) into (1) gives, where the superscript (·) refers to the partial derivative with respect to q 3 . It is clear to see that (11b) is independent of the other two equations; in fact, (11b) represents the vibration of the SH waves in orthotropic plates and (11a) and (11c) control propagation of Lamb wave modes. In order to solve the decoupled wave equations, the Legendre polynomial method expands U i (x 3 ) into an orthonormal polynomial basis [18,30,31], where p i m is the expansion coefficient and, where P m (x) is the Legendre polynomial of order m. Theoretically, m runs from 0 to ∞; however, in practice, the summation over polynomials in (12) can be halted at some finite value of m = M, when higher-order terms become negligible.
To retrieve the final equations for solution, one substitutes (12) and (13) into (11), multiplies by Q * j (q 3 ) and integrates over q 3 from 0 to kh, giving, with j and m running from 0 to M, and (·) * indicates the complex conjugate. The definitions of the matrix elements are shown in Appendix A. By separating out (14) into only the coupled Lamb wave modes and the decoupled SH wave mode, the final solution can be arranged as an eigenvalue problem, with eigenvalues −c 2 p and corresponding eigenvectors {p 1 m p 3 m } . Here, 3(M + 1) eigenmodes are generated at order M of the expansion. The only solutions to be accepted are those eigenmodes for which convergence is obtained as M is increased [18,30,31].
These equations are not individually inferrable for respective wave modes, instead the matrix elements must be determined to form the eigenvalue problem in (15), the solutions of which provide dispersion information. Solutions to all available modes must be determined simultaneously using the eigenvalue solution, where the number of modes available is 2(M + 1) for the Lamb wave modes, and M + 1 for the SH wave modes.
Previously, this method has been implemented using a symbolic-programming approach; however, as the expansion forms a series of polynomials, a programmatic approach has been developed here to reduce computation time. This strategy was in fact the first of many numerical tactics employed to reduce computation cost in order to make the method applicable to a probabilistic sampling procedure such as MCMC. More details on all the manipulations developed can be found in Appendix B.
A particularly noteworthy aspect of the numerical manipulations implemented is the use of the power iteration method to determine the eigenvalue solutions. The power method allows one to determine, in sequence, the eigenvalues from largest to smallest, at a lower computational cost [32,33]. Therefore, in order to maintain a computational efficiency, only the first two wave modes (A 0 and S 0 ) are included here.

Prior Exploration of Model
Before progressing to the procedure to identify the material properties governing (14), it is useful to explore the effects on the dispersion curves of changes in the material properties. Figure 1 shows how the curves for the fundamental Lamb wave modes are affected by changes in the elastic constants and density of the orthotropic model. The initial values of each were chosen based on those used in [30].
Changes in all elastic constants appear to be consistently stronger in the solutions for the S 0 mode. The constants C 13 and C 33 appear to have very little effect on the dispersion curve for the A 0 mode, but have significant effects on the S 0 mode. An interesting observation is on the effect of C 33 and C 55 on the S 0 curve; no significant changes appear until after the 'elbow' in the curve-the sharpness of which is unique to more complex models and does not appear in isotropic dispersion curves. For C 55 and ρ, changes appear to be stronger at higher frequencies, whereas the changes appear more consistent across the frequency range for C 11 , C 13 and C 33 .

Measuring Observations of Dispersion Curves
The first stage of the process here is to determine a set of measured values on the dispersion curve {ω,k} of the plate in question. Dispersion curves can also be determined from arbitrary plates by the use of a two-dimensional Fourier transform (2DFT); this is done by recording the surface displacement of a Lamb wave, spatially sampled along its propagation path, to generate the time-distance [t-x] space. Passing this through a 2DFT then provides a transformation to the frequency-wavenumber [ω-k] space [34]. The surface displacement of a wave at regularly spaced intervals is measured to form time-distance [t-x] data. The signals at each spatial location are then normalised and the matrix passed through a 2DFT to retrieve the frequency-wavenumber [ω-k] data. Details of the experiment setup are given in Table 1, shown in Figure 2, and detailed in the paragraphs below. Lamb waves were initiated in a glass-fibre reinforced-polymer (GFRP) plate by excitation of a 20 mm diameter piezo-electric transducer (PZT) stack actuator (Physik Instrumente P-016.20P), on the surface of the plate, as shown in Figure 2. The PZT was actuated with a chirp signal of length 1ms and upper frequency of 500 kHz, allowing for broadband excitation. A Polytec PSV-400 scanning-laser vibrometer was used to measure the outof-plane surface displacement of the induced wave-packets along a single propagation path, where the recording state was synchronised with the start of the excitation signal. Retro-reflective tape, 0.5 mm thick, was placed along these propagation paths to improve the signal-to-noise ratio.
Information on the material properties for the plate are shown in Table 2; however, it is important to note that these are for model validation purposes and are not fed into the methodology. The parameters E 11 , E 22 , G 12 , ν 12 , ν 21 , and ρ are provided by the manufacturer of the plate, and have been determined by destructive tests performed on the material. From these parameters, the elastic constants were calculated by taking the inverse of the compliance matrix S, which is formulated from Hooke's law for transversely isotropic materials under plane stress [35]. Unit GPa GPa GPa --kgm −3 GPa GPa GPa GPa A PZT stack actuator was used, as opposed to a disc, in order to improve the signalto-noise ratio; as the plate used is relatively thick, and the epoxy matrix causes rapid attenuation of the waves. The location of the PZT stack actuator was chosen in order to maximise propagation distance before reflection due to boundaries. The stack actuator was placed on the same surface on which the displacement was being measured, this being a result of the physical setup required to mount the plate and to use a stack actuator. It was placed 1/3 of the distance in the 90°direction, which would result in reflections from the upper and lower boundaries imposing on the propagation wavefield at the same time. It would initially be preferable to do the same thing for the horizontal direction also, however, the propagation velocity in this direction will be larger as it is in the direction of the fibres, and so may impose on the vertical propagation wavefield before the wave reflection at the boundary.
Data were recorded along two propagation path directions; 0°and 90°, these path directions are also shown in Figure 2. Specific details of the experimental setup are shown in Table 1, including plate dimensions and acquisition parameters. The [t-x] data were then passed through a 2DFT to form [ω-k] images at each angle; the results for 0°are shown in Figure 3. The [ω-k] image data can be interpreted as a representation of the energy content over the [ω-k] space, where darker areas represent more energy. These image data can then provide an estimate of the dispersion curve in terms of ω and k, where ridges will follow these curves [34]. As the natural frequency of the stack actuator is within the frequency range of the dispersion curves of interest, there will be disparities in energy content along this frequency range. Therefore, to improve contrast of the image, resulting in a more equal distribution of observations along the frequency axis, the image data were normalised with respect to the energy content. This was done by dividing the elements of each frequency bin vector were divided by the sum of the frequency vector content. Where U is the image data,Ũ is the normalised dispersion-curve image data, n k is the length of the frequency vector, k represents the wavenumber index and f the frequency index, The dispersion-curve image data for the GFRP plate in Figure 3 shows strongly both the A 0 and S 0 modes, as well as some information present on the S 1 and S 2 modes. Using the ridge-picking algorithm for this data, as described in [36], only the A 0 and S 0 modes are considered. From the experimental setup, the upper limit of the frequency-thickness bandwidth is 8.192 MHz-mm. However, Section 2.3 showed that the dispersion curve solutions are more sensitive to changes in the material properties at higher frequencies, therefore, the frequency-thickness bandwidth of 4.098 MHz-mm was chosen to include all information available on the A 0 and S 0 modes.

Estimating Elastic Constants
Consider the concept of an elementary system-identification procedure to estimate a set of parameters given a vector y of n observations, y = {y 1 , y 2 , . . . , y n } When using a probabilistic approach to determining elastic constants, one must form a definition of the likelihood of a set of constants, given some observed data. Here, the likelihood is defined based on the choice of noise model, i.e., Gaussian. This likelihood could be used to retrieve a maximum likelihood estimate [37], which is a popular and asymptotically optimal statistical approach to fitting model parameters using data [38]. Assuming the model is of the form, where r i is the mean at point i, and ε i is a white Gaussian noise process. The observations are distributed as y ∼ N (r, σ 2 ). The likelihood is then defined as, For determining the likelihood of some model, the mean can be replaced with a function of the input dimension x of the observations, and some parameters Θ, so the likelihood becomes, In Section 3.1, it was noted that there is a much lower relative resolution in the wavenumber dimension of the dispersion image and in the resulting selected points on the curve. This implies that the Gaussian white-noise distribution is mostly in ω; thus, if one were to estimate the likelihood based on a model of k(ω, Θ), the function would be of the form, Therefore, the problem is formulated as based of a model of ω(k, Θ). The observations are taken as the points on the dispersion curve, whereω i andk i are the measured values of frequency and wavenumber respectively, at point i. The set of observations is grouped into m modes, individually represented by ψ; y = {y ψ 1 , y ψ 2 , . . . , y ψ m }, where y ψ = {ω ψ ,k ψ }. For example, in the case where only the fundamental modes are considered, ψ 1 and ψ 2 represent the A 0 and S 0 modes respectively. The likelihood is then defined as, In this case, ω(k i , Θ) is determined using the methods outlined in Section 2.2. For the work presented here, only the first antisymmetric mode A 0 is considered, and so only solutions for that curve are returned. The parameters are defined as the elastic constants which enter into the equations for the Lamb wave modes in (11a) and (11c), as C 31 = C 13 . Maximising L(y|Θ) provides an estimate of the most likely elastic constants; however, it is also possible to retrieve information on their distribution.

Estimating the Posterior Distributions
The objective at this stage is to determine the distribution of the parameters which could plausibly define the dispersion curve. As the likelihood includes a noise variance term σ, the parameter vector is extended to include this, so that, The distribution of these parameters can be determined by identifying the posterior probability given a set of measured data, p(θ|y). However, this is not directly inferable, so a manipulation is done using Bayes rule, where p(y|θ) is calculated using (24), and p(θ) is the prior, which can be defined using initial knowledge of the parameters. For d parameters, assuming each parameter is independent, the prior is calculated as, Now, the problem is transferred, in that the normalisation term in the denominator is intractable. Instead, a procedure can be done to sample from the posterior with enough repetition that an estimate of the distribution over the parameters can be inferred. One such procedure is the Markov-Chain Monte Carlo (MCMC) method, where subsequent samples depend on assessing their probability with respect to the previous one. An outline of the derivation and procedure for MCMC is given in [39][40][41]. In practice, for computational stability, the probabilities are calculated in the log space, so the marginal likelihood becomes, log(p(θ|y)) = log(p(y|θ)) + log(p(θ)) (29) where, Now, consider how to define this problem for the application to dispersion curve material identification. The first step is to define the likelihood, which is done using (24), which in the log space is, During sampling using MCMC, the size of the random step taken for each parameter is important as too large a step will cause stall, and too small a step will require a large number of iterations. An improvement is made on the standard MCMC procedure, which incorporates Hamiltonian mechanics, to adapt the step size for an optimal simulation, and is so called Hamiltonian Monte Carlo (HMC) [42,43]. For this work, the probabilistic programming language Stan [44], was used to perform the actual sampling procedure.
Next, consider the definition of the priors, which can be done using reasonable knowledge of the material of application. As the prior is a combination of the individual probabilities of each parameter, prior belief on the distribution of these parameters can be used to define each p(θ i ). In this case, the density of the plate is supplied, but no other material properties were provided. Therefore, a tight prior can be given on ρ and priors on the elastic constants are defined to capture reasonable values for the material. Here, a gamma distribution was chosen for the elastic constants as this provides a broad definition of the prior, which can be interpreted as embedding belief on the magnitude of the value, and enforces a positive-only value. As the prior for the density can be defined relatively tightly, this was defined using a normal distribution. The type and definitions of the priors used here are shown in Table 3. Table 3. Definitions of priors for parameters in θ.

Parameter
Distribution Definition

Results
In this section, samples from the posterior distributions of the parameters are shown in both univariate and bivariate distributions, and a kernel density estimate is used to estimate the probability density function of the bivariate distributions. Also shown are samples of the dispersion curves drawn from the samples of the posterior distributions, overlaid onto the observed dispersion-curve image data taken from the two-dimensional Fourier transform of the measured surface displacement. This section is split into three subsections; firstly, the simulated posterior distributions of the parameters are shown, followed by a simulated distribution of dispersion-curve data based on these parameters. Lastly, the first two statistical moments of the univariate samples are calculated; which are then used to display the estimated mean and variance of each parameter for each plate. At this stage, one is not dealing with time-signals, and is instead utilising only observations from frequency-domain data, therefore, the use of the term 'samples' is with reference to samples (or iterations) drawn from the MCMC procedure.

Posterior Distribution of the Material Parameters
The results of 20,000 accepted samples of the sampling procedure for propagation angles of 0°and 90°are shown in Figures 4 and 5 respectively. The first observation that can be made is of the evidence of correlation between all material parameters, whereas the distribution of the noise parameter appears to converge to an independent distribution. This result is anticipated, as the elastic properties which form the stiffness matrix are described by a series of inseparable equations.   Figure 4 indicates the univariate and bivariate distributions for a propagation angle of 0°. There is an apparent 'edge' on the scatter correlation plots between certain parameters, in particular between C 13 and all other elastic constants. As a condition of the solution to the dispersion curve equations is that λ < 0, any solutions where this is not the case are rejected. The edge may indicate a region of forbidden parameter combinations which cannot exist, given a real elastic material.
There is evidence of a particularly strong correlation between C 55 and ρ, which appears to be a linear relationship. This could be explained by comparison to the isotropic case; for an isotropic material their relationship can be defined as C 55 = ρc 2 T . For an orthotropic material, the transverse-wave velocity would remain the same when rotating around the axis in the direction of wave propagation. This property could be used to reduce the number of parameters, increasing performance of the simulation.
One observation made in the results of the sample plots for the 90°propagation, in comparison to 0°, is the less apparent hard edge caused by the rejection parameter. This may indicate that the posterior space of valid elastic constants is less discontinuous when modelling Lamb-wave propagation through fibres. Another difference between the two sets of results is that there is a less strong correlation in the parameter pairs C 11 -C 13 and C 13 -C 33 . This could be a result of the fibres no longer acting as a secondary wave guide, and instead shear forces through the fibres have more of an influence on wave propagation.

Distribution of Dispersion Curve Models
Using the parameters at each sample point, a distribution of dispersion curves was generated, and is shown in Figure 6, along with observations taken from the [ω-k] image data. For the propagation angle of 0°, the darker areas of the image data, as well as the observation points, lie within the distribution well for both fundamental wave modes. This result shows that the method works well for obtaining dispersion characteristics of Lamb waves. Although the coupon used here has unidirectional fibre, an orthotropic model was still used for the data at a propagation angle of 90°to test its applicability to all directions. Figure 6b indicates that for determining dispersion curves of the A 0 mode, this model still provides a useful solution. However, the curves generated for the S 0 mode become mismatched from the image data, as well as the points taken from the ridge-picking algorithm, at a frequency-thickness greater than 1.2 MHz-mm. This indicates that the model used is insufficient for modelling the S 0 mode, however, it still generates a reasonable model for the dispersion curves for the A 0 mode.
As stated at the beginning of this paper, a key advantage of the method shown here is the freedom in the posterior distribution, as no assumption is made as to its shape. In an engineering context, this allows freedom in the material type to be modelled, so long as the model of the dispersion curve solutions is accurate. For both propagation angles here, the univariate distributions of the parameters do not all appear to be of the same shape. In fact, the elastic constants and density appear to converge to a Gamma distribution of varying skewness, and the noise parameter appears to converge to a normal distribution. In Figures 4 and 5, all elastic constants appear to converge to Gamma distributions; this indicates that the true posterior of the elastic constants should converge to a skewed distribution.

Quantifying the Posterior Distributions
For each of the parameters, the expected value and variance are calculated as the first two arithmetic statistical moments. Using a kernel density estimate, the mode of the distributions is also calculated and designated as the most likely estimate. The values were calculated for the samples from the posterior distributions for both angles, the results of which are shown in Table 4. A notable observation here is the relatively large discrepancy between the mean values for density; for a propagation angle of 90°it is predicted to be much lower. This discrepancy may be as a result of the model having to 'counteract' any conflict between the model form being fit and the data values.
Relative to the mean value, the standard deviations of each parameter are similar, which aligns well with the observed posterior distributions that are seen in Figures 4 and 5. This observation can be interpreted as the level of uncertainty being similar for each parameter, meaning that discrepancies are not confined to a single parameter, but are instead in the combination of parameters. It is important to note that, in the prior, each parameter is treated as independent, whereas the posterior shows that there is strong co-dependence between the parameter's dispersion-curve solutions. When comparing the results at 0°to those shown in Table 2, some further discussion can be made on the result and the advantages of the methodology. The results compare reasonably well with the provided values, however, there is still some discrepancy. This may be a result of there not being enough observations provided from the higher-frequency range, which was shown in Section 2.3 to be more sensitive to changes in the elastic constants. It is, however, important to note the much less accurate values obtained if one were to choose the most likely estimate, in comparison to choosing the mean value. This, once more, shows the advantage of estimating the posterior distribution, as opposed to simply determining the most likely elastics constants obtained by maximising (24).
As the aim of the work here is to determine the dispersion characteristics useful for NDE/SHM strategies, a key motivation of which is to find the group velocity of the waves, it is also useful to look at the distribution of curves for this attribute. During the same curve sample-drawing procedure as above, the value of c g was also calculated as the slope of the generated [ω-k] curves. The distributions of the [ω-c g ] curves for propagation angles of 0°and 90°are shown in Figure 7. Much like the curves seen in Figure 6, the distribution of the A 0 mode is much tighter than that of the S 0 mode. For the propagation angle of 0°, the curves have no discontinuities and appear to have a low uncertainty.
In Figure 7b, there are some discontinuities of the curve in the range 1.25 < f h < 1.8 MHz-mm. From (2), the group velocity is taken as the gradient of the [ω-k] curves. By inspecting Figure 6b, one can see the gradient of the curve changes rapidly, which is the model solution with the highest likelihood. From the results shown, the method presented here returns accurate and precise models of the dispersion curves for an arbitrary orthotropic plate. The objective problem of the work here was to determine dispersion-curve information on the fundamental modes, as this is the necessary information required for guided wave-based localisation. For the purposes of determining the dispersion characteristics from the data provided, it performs well.

On the Confidence of the Results
In this paper, the capability of the method has been shown with respect to the initial objective; to determine accurate dispersion curve information for an unknown material, or if the dispersion curve information needs updating. It is important to note, that even with confident estimates of these dispersion curves, the estimated variance in each parameter is still reasonably large. This indicates that using alignment of estimated dispersion curves, using the most likely estimate of the parameters, is not enough to reasonably state confidence in these parameters. In Section 2.3, it can be observed that relatively large changes in the material properties result in small changes in the dispersion curve. The work here shows the importance of obtaining the posterior distribution of the parameters, as opposed to just the most likely estimate.
Furthermore, the results for the dispersion curve observations for the 90°data show that the model used here in unsuitable for this propagation direction. The model used in the solution equation is for propagation in the direction of the fibres, and these results show there is no combination of material properties that can accurately model dispersion curves for propagation through the fibres. Therefore, even though the posterior distribution appear to have settled to a reasonable shape and values, by observing the dispersion curve solutions, one can see the inapplicability of the model. These two remarks show the importance of using both the posterior distributions, estimated values/variances, and the dispersion curve solutions to properly assess the results of this procedure.
Another important point to note regards the information supplied to the sampling procedure. In this paper, a ridge-picking algorithm was used to generate a vector of {ω, k} observations; from Figure 6, it can be observed that the distribution of the generated curves lies within the distribution of these observations. The use of only the A 0 and S 0 modes-although important for maintaining computational efficiency-decreases the confidence in the posterior estimates, as a result of the stronger influence of the parameters on higher-order modes. However, this point also highlights the advantage of this method; where the there is a lack of confidence from the observations, it is important to consider the posterior distribution, as opposed to simply the most likely value.

Conclusions
The aim of this paper was to develop a Bayesian approach to material identification for the purposes of determining dispersion curve models for an orthotropic plate. The Bayesian approach allows for inference on the posterior distribution of the material properties, and for total freedom of the distribution shape, as well as inferring any multivariate correlations between parameters. By determining the multi-variate posterior distributions of the elastic constant space, it was shown that it is also possible to generate distributions of the dispersion-curve solutions. This is another key advantage of the method, as uncertainty bounds on the curve can be propagated through directly to uncertainty in measurements made using these curves-such as wave source localisation. The initial objective for the problem was to determine the dispersion characteristics of the fundamental modes that are important for the purposes of damage localisation in NDE/SHM strategies. The results of the curve distributions indicate that the method works well to achieve this objective. Future work intended by the authors has been discussed following analysis of the results.

Future Work
As discussed, the method presented here returned reliable and robust results for the objective problem. For further application of this method to more exhaustive material identification procedures-such as full elastic-constant identification or SH dispersioncurve information-some additions are necessary. The key aim for these improvements is to increase the fidelity of the information provided to the procedure in order to allow inference from the additional wave modes. An important consideration of this, however, will be the increased computational cost of calculating solutions for more modes. The authors intend to explore a number of approaches to address this objective; one such approach is to include rotation of the stiffness matrix combined with using multiple angles for a single observation set, and run a single parameter identification routine. Adding further, the Legendre polynomial expansion approach can readily include damping characteristics by using a complex-valued stiffness matrix; by obtaining complex observations of the wavenumber, it would be possible to extend the identification procedure to simulate the posterior distributions on the real-imaginary pairs for each elastic constant. Another approach is to develop a multi-dimensional prior definition for the elastic-constant space, which would improve sampling efficiency. The final approach to be explored is to adapt the likelihood to use full 2DFT image data, rather than individual observations taken from the image data using a ridge-picking algorithm.

Conflicts of Interest:
The authors declare no conflict of interest. Section 2.2 shows the method of solving dispersion curves for anisotropic materials using the Legendre polynomial expansion approach. The final system of equations which forms an eigenvalue problem are given in (14), the matrix elements of which are given below, where, and the sifting property of a Dirac delta function states, which leads to, NT 2 (m, j, n) = f n (0) − f n (kh) (A13)

Appendix B.2
During solution of the dispersion-curve equations, the matrix A is produced. This matrix is formed from 4 sub-matrices A 11 , A 31 , A 13 and A 33 , all of size 2(M + 1) × 2(M + 1), where A 11 , A 33 ∈ RM ,M and A 13 , A 31 ∈ CM ,M where Re(A 13 ) = Re(A 31 ) = 0. The use of complex numbers greatly increases computational cost. However, by some manipulation, the requirement for complex variable types could be eradicated, decreasing computational load. The eigendecomposition problem is formulated such that, |A − λI| = 0 (A15) whereM = M + 1, λ 1 = λ{1 :M} and λ 2 = λ{M + 1 : 2M}, as A 13 , A 31 are purely imaginary, the following manipulation is then applied, |A 13 ||A 31 | = |iIm(A 13 )||iIm(A 31 )| = | − Im(A 13 )||Im(A 31 )| (A18) and thus the matrix A can be reformulated as, The returned eigenvectors λ of the matrix A, give solutions for 2M modes where, Each eigenvalue can be related to the phase velocity of each mode by, where j = 2M − i + 1. The only physically possible solutions are c p ∈ R >0 ; therefore, only solutions where λ ∈ R <0 are viable. The minimum-magnitude eigenvalue will represent the smallest-value solution of c p . As here, data points for the A 0 mode are being used, only the smallest value of c p is needed (and thus smallest-magnitude eigenvalue). This leads to the use of the power iteration method [32], which is a method that returns only the dominant eigenvalue of a matrix, at reduced computational expense. Unfortunately, the dominant eigenvalue of a matrix A is the largest magnitude, which would return the largest value of c p . However, the smallest magnitude eigenvalue can be returned by determining the dominant eigenvalue of A −1 .