A Statistical Approach to Surface Metrology for 3D-Printed Stainless Steel

Abstract

Surface metrology is the area of engineering concerned with the study of geometric variation in surfaces. This article explores the potential for modern techniques from spatial statistics to act as generative models for geometric variation in 3D-printed stainless steel. The complex macro-scale geometries of 3D-printed components pose a challenge that is not present in traditional surface metrology, as the training data and test data need not be defined on the same manifold. Strikingly, a covariance function defined in terms of geodesic distance on one manifold can fail to satisfy positive-definiteness and thus fail to be a valid covariance function in the context of a different manifold; this hinders the use of standard techniques that aim to learn a covariance function from a training dataset. On the other hand, the associated covariance differential operators are locally defined. This article proposes to perform inference for such differential operators, facilitating generalization from the manifold of a training dataset to the manifold of a test dataset. The approach is assessed in the context of model selection and explored in detail in the context of a finite element model for 3D-printed stainless steel.

Keywords:

• Finite element
• Gaussian Markov random field
• Natural gradient
• Stochastic partial differential equation

1 Introduction

The rate at which new manufacturing technologies are being developed and the range of components that can now be produced pose a challenge to formal standardization, which has traditionally formed the basis for engineering design and assessment. An important motivation for new technologies is to trade the quality of the component with the speed and cost at which it can be manufactured. This tradeoff results in the introduction of imperfections or defects, in a manner that is intrinsically unpredictable and could therefore be described as statistical. As such, there is a need to develop principled and general statistical methods that can be readily adapted to new technologies, in order to obtain reliable models as the basis for formal standards in the engineering context.

The recent emergence of additive manufacture technologies for stainless steel promises to disrupt traditional approaches based on standardized components and to enable design of steel components of almost arbitrary size and complexity (Buchanan and Gardner 2019). This directed energy deposition technology arises from the combination of two separate technologies; traditional welding modalities and industrial robotic equipment (ASTM Committee F42 2021). In brief, material is produced in an additive manner by applying a layer of steel onto the surface of the material in the form of a continuous weld. To ensure that this process is precisely controlled, the manufacture is performed by a robot on which the weld head is mounted. The complexity of this process is such that the geometric (e.g., thickness and volume) and mechanical (e.g., stiffness and strength) properties of additively manufactured (henceforth “3D-printed”) steel are not yet well-understood. Indeed, conventional models for welding focus on the situation of a localized weld that is typically intended to bind two standard components together. The material properties of elementary welds, as a function of the construction protocol and ambient environment, have been modeled in considerable detail (Kou 2003). In contrast, 3D-printed steel is a single global weld whose construction protocol is highly nonstandard. Furthermore, the limited precision of the equipment results in imperfections in the material geometry, which can range in severity from esthetic roughness of the surface of the material through to macroscale defects, such as holes in the printed component. These imperfections occur in a manner that is best considered as random, since they cannot be easily predicted or controlled.

The variation in the manufactured components precludes the efficient use of raw material, as large factors of safety are necessarily employed. An important task is therefore to gain a refined understanding of the statistical nature of 3D-printed steel, enabling appropriate factors of safety to be provided for engineering design and assessment (Buchanan et al. 2017). It is anticipated that geometric variation, as opposed to mechanical variation, is the dominant source of variation in the behavior of a 3D-printed thin-walled steel component under loading. The focus of this work is therefore to take a statistical approach to surface metrology, with the aim to develop a generative model for geometric variation in 3D-printed steel. Specifically, the following problem is considered:

Problem Statement. Given the notional and actual geometries of a small number of 3D-printed steel components, construct a generative statistical model for the actual geometry of a 3D-printed steel component whose notional geometry (only) is provided.

For example, in Section 4.1 of this article a statistical model is used to generate realistic simulations of geometric variation that might be encountered on a 3D-printed steel cylinder, based only on a training dataset of notionally flat panels of 3D-printed steel; see Figure 1.

Fig. 1 Illustration of the problem statement. Here, the actual geometries of a small number of 3D-printed steel components, each notionally a flat panel, are provided and the aim is to construct a generative statistical model for the actual geometry of a 3D-printed steel component whose notional geometry is cylindrical.

1.1 A Brief Review of Surface Metrology

Surface metrology concerns small-scale features on surfaces and their measurement (Jiang et al. 2007a, 2007b). The principal aim is to accurately characterize the nature of the surface of a given type of material, in order that such material can be recreated in simulations or otherwise analyzed. Typical applications of surface metrology are at the nanometer-to-millimeter scale, but the principles and methods apply equally to the millimeter-to-centimeter scale variation present in 3D-printed stainless steel.

The first quantitative studies of surfaces can be traced back to the pioneering work of Abbot and Firestone (1933), and concern measurement techniques to produce one-dimensional surface cross-sections (or “profiles”), denoted $z:[0,L]\to \mathbb{R}$, and the subsequent description of profiles in terms of a small number of summary statistics, such as amplitude and texture statistics (Thomas 1999). Typical examples of traditional summary statistics are presented in Table 1. The use of summary statistics for surface profiles continues in modern surface metrology (e.g., as codified in standard ISO/TC 213 2015). These include linear filters, for example based on splines and wavelets (Krystek 1996), and segmentation filters that aim to partition the profile into qualitatively distinct sections (Scott 2004). See Jiang et al. (2007a).

Table 1 Summary statistics traditionally used in surface metrology; here $z:[0,L]\to \mathbb{R}$ is the profile through the surface of the material and μ is its average or intercept as per (1).

Parameter	Parameter
name	type	Description	Mathematical definition
Rq	Amplitude	Root mean square roughness	Rq $:=\sqrt{\frac{1}{L}{\int }_0^L{[z(x)-\mu ]}^2\mathrm{d}x}$
Ra	Amplitude	Average roughness	Ra $:=\frac{1}{L}{\int }_0^L|z(x)-\mu |\mathrm{d}x$
HSC	Texture	Number of local peaks per unit length	(NA—this depends on definition of a peak)
${\mathrm{R}}_L$	Texture	Ratio of developed length of a profile compared to nominal length	${\mathrm{R}}_L:=\frac{1}{L}{\int }_0^L\sqrt{1+z\prime {(x)}^2}\mathrm{d}x$

NOTE: Adapted from (Thomas 1999, chap. 7 and 8) .

Although useful for the purposes of classifying materials according to their surface characteristics, these descriptive approaches do not constitute a generative model for the geometric variation in the material. Subsequent literature in surface metrology recognized that surfaces contain inherently random features that can be characterized using a stochastic model (Patir 1978; DeVor and Wu 1972). In particular, it has been noted that many types of rough surface are well-modeled using a nonstationary stochastic model (see, e.g., Majumdar and Tien 1990). Traditional stochastic models for surface profiles fall into two main categories: The first seeks to model the surface profiles as time-series, typified by the autoregressive moving average (ARMA) model introduced in this context in DeVor and Wu (1972) and extended to surfaces with a non-Gaussian height distribution by Watson and Spedding (1982). An ARMA(p, q) model combines a pth order autoregressive model and a qth-order moving average(1) $$z_i=\mu +{ϵ}_i+\sum _{k=1}^p{\varphi }_k{z}_{i-k}+\sum _{l=1}^q{\theta }_l{ϵ}_{i-l},$$(1) where $z_i=z(t_i)$ represents a discretization of the surface profile and the ϵ_i are residuals to be modeled. These techniques have been applied to isotropic surfaces. Patir (1978) had noted that the theory is less well-understood regarding anisotropy and pointed out that observations of surface profiles in differing directions may be used to establish anisotropic models. This represents a significant weakness where simulation of rough surfaces are concerned, as it becomes unwieldly and complex to simulate materials with a layered surface such as 3D-printed steel. Patir (1978) introduced a second class of methods, wherein a two-dimensional surface is represented as an N × M matrix of heights $(z_{i,j})$, assumed to have arisen as a linear transform of an $(N+n)\times (M+m)$ random matrix $({\eta }_{i,j})$ according to$$z_{i,j}=\mu +{ϵ}_i+\sum _{k=1}^n\sum _{l=1}^m{\alpha }_{k,l}{\eta }_{i+k,j+l},$$ where the coefficients $({\alpha }_{k,l})$ are selected in order to approximate the measured autocorrelation function for the $z_{i,j}$. The surface heights are often assumed to be Gaussian as the theory in this setting is well-understood. A number of authors have developed this method: A fractal characterization was explored in Majumdar and Tien (1990), recognizing that the randomness of some rough surfaces is multiscale and exhibits geometric self-similarity. Hu and Tonder (1992) proposed the use of the fast Fourier transform on the autocorrelation function to facilitate rapid surface simulation; see also Newland (1980) and Wu (2000). A possible drawback with this second class of approaches, compared to the former, is that it can be more difficult to account for nonstationarity compared to for example, the ARMA model. The use of wavelets to model nonstationary and multi-scale features has received more recent attention, see Jiang et al. (2007b).

Recent research in surface metrology concerns the simulation of rough surfaces for applications in other fields, including tribology (Wang et al. 2015), electronics (Shen et al. 2016), and wave scattering (Choi et al. 2018). Typically, the approaches above are used directly, but alternative methods are sometimes advanced. For instance, Wang et al. (2015) proposed an approach using random switching to generate Gaussian surfaces with prescribed statistical moments and autocorrelations. Clearly, it is necessary to specialise any statistical model to the specific material under study. For example, in the case of 3D-printed steel an approach ought to be taken which is able to natively account for anisotropy due to layering of the steel. More crucially in the present context, existing methods pre-suppose that surface profiles $(z_i)$ or height matrices $(z_{i,j})$ are pre-defined with respect to some regular Cartesian grid on a Euclidean manifold. In applications to 3D-printed steel, where components have possibly complex macroscale geometries, the training and test data can belong to different manifolds and it is far from clear how an existing model, once fitted to training data from one manifold, could be used to simulate realistic instances of geometric variation on a different manifold.

An important application of surface metrology is surface quality monitoring; a classification task where one seeks to determine, typically from a limited number of measurements and in an online environment, whether the quality of a manufactured component exceeds a predefined minimum level (Zang and Qiu 2018a,b; Jin, Pierson, and Liao 2019; Zhao and Del Castillo 2021). Existing research into quality monitoring has exploited Gaussian process models to reconstruct component geometries from a small number of measurements (Xia, Ding, and Wang 2008; Colosimo et al. 2014; Del Castillo, Colosimo, and Tajbakhsh 2015), possibly at different fidelities (Colosimo, Pacella, and Senin 2015; Ding et al. 2020). Going further, performance metrics arising from surface metrology have been used as criteria against which the parameters of a manufacturing protocol can be optimized (Yan, Paynabar, and Pacella 2019).

1.2 Outline of the Contribution

The aim of this article is to explore whether modern techniques from spatial statistics can be used to provide a useful characterization of geometric variation in 3D-printed stainless steel.

To provide a training dataset, a high-resolution laser was used to scan 3D-printed stainless steel panels, on whose surfaces the geometric variation of interest is exhibited. A typical panel is displayed in Figure 2. To characterize geometric variation in the panel dataset, in Section 2 a library of candidate models is constructed, rooted in Gaussian stochastic processes, and in Section 3 formal parameter inference and model selection are performed to identify a most suitable candidate model. Although quite natural from a statistical perspective, this systematic approach to parameter inference and model selection appears to be novel in the surface metrology context. The high-dimensional nature of the dataset introduces computational challenges that are familiar in spatial statistics and the stochastic partial differential equation (SPDE) techniques described in Lindgren, Rue, and Lindström (2011) are leveraged to make the computation practical.

Fig. 2 Laser scan of 3D-printed steel sheet. (a) Photograph of the handheld scanning equipment. (b) Orthographic projection of (a portion of) the scanned sheet, which will be called a panel. The notional thickness of the panel is 3.5mm. [Observe the residual stress induces a slight curvature in the notionally flat panel.]

The aim is to characterize geometric variation not just on panels but on components whose notional geometry, represented abstractly by a two-dimensional manifold $\mathcal{M}$ embedded in ${\mathbb{R}}^3$, may be complicated. (The case study in Section 4.1 considers the case where $\mathcal{M}$ is cylindrical which, whilst not complicated, differs from a flat panel in a fundamental topological way.) This variation, present in the actual geometry of a 3D-printed component, is described by a two-dimensional vector field $u:\mathcal{M}\to {\mathbb{R}}^2$. As illustrated in Figure 3, the value of the first coordinate $u^{(1)}(x)$ of the vector field represents the Euclidean distance of the upper surface of the component from the manifold and similarly the value of $u^{(2)}(x)$ represents the distance of the lower surface of the component from the manifold. The absence of geometric variation corresponds to the functions $u^{(1)}(x)$ and $u^{(2)}(x)$ being constant, but actual 3D-printed steel components exhibit geometric variation that can be statistically modeled. The need to consider diverse notional geometries, represented by different manifolds, precludes conventional parametric (Rasmussen and Williams 2006), compositional (Duvenaud 2014) and nonparametric (Băzăvan, Li, and Sminchisescu 2012; Oliva et al. 2016; Moeller et al. 2016) approaches that aim to infer an appropriate covariance function defined only on the manifold of the training dataset. Indeed, it is well known that a radial covariance function $\phi (d_{{\mathcal{M}}_1}(x,y))$ defined on one manifold $x,y\in {\mathcal{M}}_1$ with geodesic distance $d_{{\mathcal{M}}_1}$ can fail to be a valid covariance function when applied as $\phi (d_{{\mathcal{M}}_2}(x,y))$ in the context of a different manifold $x,y\in {\mathcal{M}}_2$ with geodesic distance $d_{{\mathcal{M}}_2}$. This occurs because positive-definiteness can fail; see Section 4.5 of Gneiting (2013). The key insight used in this article to generalise across manifolds is that the Laplacian ${\Delta }_{{\mathcal{M}}_1}$ on ${\mathcal{M}}_1$ can be associated with the Laplacian ${\Delta }_{{\mathcal{M}}_2}$ on ${\mathcal{M}}_2$. In doing so a suitable differential operator on ${\mathcal{M}}_1$, describing the statistical character of the vector field u via an SPDE, can be inferred and then used to instantiate an equivalent differential operator on ${\mathcal{M}}_2$, defining a Gaussian random field on ${\mathcal{M}}_2$ via an SPDE. This allows, for instance, the simulation of high-resolution random surfaces on a cylinder whose local characteristics are related to the panel training dataset. (Note that this article does not consider the question of when the covariance differential operator/SPDE approach defines a true random field rather than a generalized random field. The practical consequences of this are minimal, however, since the random field will be discretized by projection onto a finite element basis whose coefficients are well-defined linear functionals of the random field.)

Fig. 3 Data are conceptualized as a two-dimensional vector field u on an oriented manifold $\mathcal{M}$. The value $u^{(1)}(x)$ represents the height of the upper surface of the material above $x\in \mathcal{M}$, while the value $u^{(2)}(x)$ represents the height of the lower surface below $x\in \mathcal{M}$.

In Sections 4.1 and 4.2, it is demonstrated, using a held out test dataset consisting of a 3D-printed steel cylinder, that random surfaces generated in this way can capture some salient aspects of geometric variation in the material. However, the conclusions of this article are tempered by the substantial computational challenges associated with the fitting of model parameters in this context. The article concludes with a brief discussion in Section 5.

This research forms part of a wider effort to develop realistic computer modeling tools for applications involving 3D-printed steel; these must also account also for other factors such as mechanical variation and aspects of the printing protocol (Dodwell et al. 2021). Specifically, other controllable factors that affect the quality of the manufactured component include the angle of the weld head, the speed of printing and the order in which the steel layers are printed. It is possible that an improved understanding of the nature of geometric variation can inform the improvement of manufacturing techniques and their protocols (Barclift and Williams 2012; Franco et al. 2017), but these directions are beyond the scope of the present article, which aims only to characterise the geometric variation in the material.

2 Statistical Models for Geometric Variation in 3D-Printed Steel

This section proceeds as follows: Section 2.1 contains the main conceptual contribution, that identification of Laplacians can be used to generalize from the training manifold to a test manifold. Section 2.2 presents the mathematical setting of a two-dimensional Gaussian random field characterized by a differential operator and Section 2.3 presents the resulting collection of candidate models.

2.1 Gaussian Models and Generalization Via the Laplacian

In this article, geometric variation for a 3D-printed component with notional geometry $\mathcal{M}$ is modeled using a Gaussian stochastic process. Letting $\mathbb{E}$ and $\mathbb{C}$ denote expectation and covariance operators, set$$\begin{array}{cc}\mathbb{E}[u^{(r)}(x)]\hfill & =0\hfill \\ \mathbb{C}(u^{(r)}(x),u^{(s)}(y))\hfill & =k_{\mathcal{M}}((x,r),(y,s)),\hfill \end{array}$$where the covariance function $k_{\mathcal{M}}((x,r),(y,s))$ describes the statistical relationship between the upper and lower surfaces ($r,s\in \{1,2\}$) at different spatial locations ($x,y\in \mathcal{M}$). To describe geometric variation in the panel datasets, whose notional flat geometry is denoted ${\mathcal{M}}_1$, any existing method for estimating a suitable covariance function $k_{{\mathcal{M}}_1}$ from the dataset could be employed. Such methods are well-developed and examples include maximum likelihood estimation over a parameterized class of covariance functions (Rasmussen and Williams 2006), inference based on cross-validation with arbitrarily complicated compositions of simple kernels (Duvenaud 2014) and even fully nonparametric inference of a covariance function (Băzăvan, Li, and Sminchisescu 2012; Oliva et al. 2016; Moeller et al. 2016). However, in doing so one would be learning a covariance function that is defined on the manifold ${\mathcal{M}}_1$ of the training dataset, being the Green’s function of a stochastic differential equation defined on ${\mathcal{M}}_1$ (see, e.g., Fasshauer and Ye 2013). As a consequence, one cannot expect a covariance function learned from the training dataset to carry the same meaning when applied on the manifold ${\mathcal{M}}_2$ of the test dataset. Indeed, it is a classic result that the Matérn kernel $\phi$ on ${\mathbb{R}}^3$ with $d_{{\mathbb{R}}^3}(x,y)=⏧x-y⏧$ can fail to be positive definite when applied as $\phi (d_{{\mathbb{S}}^2}(x,y))$ on the sphere $x,y\in {\mathbb{S}}^2\subset {\mathbb{R}}^3$ with $d_{{\mathbb{S}}^2}$ the geodesic distance on ${\mathbb{S}}^2$ (see, e.g., Gneiting 2013, sec. 4.5). Thus, inferential approaches based on the covariance function are unsuitable in this context, since a description is required of geometric variation on a range of manifolds that includes potentially complicated manufacturing geometries.

The key insight in this article is that, whilst covariance functions do not generalize across manifolds, their associated differential operators are only locally defined. In particular, a wide range of models for stochastic vector fields u can be cast as the solutions of a coupled system of SPDEs based on a Laplacian-type differential operator. The specific form of the differential operator, together with suitable boundary conditions, determine the statistical properties of the associated random field. The Laplacian ${\Delta }_{{\mathcal{M}}_1}$ on the manifold ${\mathcal{M}}_1$ of the training dataset can be associated with the Laplacian ${\Delta }_{{\mathcal{M}}_2}$ on the manifold ${\mathcal{M}}_2$ where generative models are to be tested. This provides a natural way in which to decouple the statistical properties of the random vector field u from the manifold on which it is defined. To sound a cautionary note, this assumption is likely to break down in situations where one of the manifolds is highly curved. For example, it is conceivable that the protocols used to control weld head movement may be less effective when the required trajectories are highly curved, leading to higher levels of geometric variation in corresponding regions of the component. As a more extreme example, if two distinct regions of the manifold come within millimetres due to the manifold curving back on itself when embedded in 3D, then through stochastic fluctuation there is a potential for the two corresponding regions of the printed component to become physically fused. The first of these phenomena could be handled by considering Gaussian models that allow covariances to be both position- and curvature-dependent, while an entirely different approach would be needed to model accidental fusing of distinct regions of a component. Neither of these extensions are pursued in the present work.

The SPDE formulation of Gaussian random fields has received considerable attention in the spatial statistics literature following the seminal work of Lindgren, Rue, and Lindström (2011), generating an efficient computational toolkit which has been widely used (e.g., Coveney et al. 2019). Interestingly, the present account appears to describe the first instance of transfer learning using the SPDE framework, where observations made on one manifold are used to provide information on the equivalent stochastic process defined on a different manifold. As noted earlier, no attempt is made to address the technical issue of whether the covariance operator corresponds to a true random field rather than a generalized random field.

2.2 Mathematical Set-Up

The purpose of this section is to rigorously set up a coupled SPDE model for a random vector field u on a manifold $\mathcal{M}$. For completeness, the required regularity assumptions are stated in Section 2.2.1; these are essentially identical to the presentation in Lindgren, Rue, and Lindström (2011) and may be safely skipped if desired. The coupled SPDE model is introduced in Section 2.2.2, following a similar presentation to that of Hu et al. (2013a,b).

2.2.1 Regularity Assumptions

The domain of the printing process is ${\mathbb{R}}^3$ and throughout the article the convention is used that the printing occurs in the direction of the z-axis of ${\mathbb{R}}^3$. That is, each layer of steel is notionally perpendicular to the z-axis according to the printing protocol. The notional geometry of a component is described by a manifold $\mathcal{M}$ embedded in ${\mathbb{R}}^3$. Below are listed some strong regularity conditions on $\mathcal{M}$ that ensure the coupled SPDE model in Section 2.2.2 is well-defined.

Let $\mathcal{M}$ be a 2-dimensional manifold embedded in ${\mathbb{R}}^3$, so that the Hausdorff measure on $\mathcal{M}$, denoted ${\mu }_{\mathcal{M}}$, can be defined, and thus an inner product$${〈\varphi ,\psi 〉}_{\mathcal{M}}:={\int }_{\mathcal{M}}\varphi (x)\psi (x)\mathrm{d}{\mu }_{\mathcal{M}}(x)$$and an associated space $L^2(\mathcal{M})$ of functions $\varphi :\mathcal{M}\to \mathbb{R}$ for which ${〈\varphi ,\varphi 〉}_{\mathcal{M}}<\infty$. It is assumed that $\mathcal{M}$ is endowed with a Riemannian structure, induced by the embedding in ${\mathbb{R}}^3$, which implies that one can define differential operators including the gradient ${\nabla }_{\mathcal{M}}$, the Laplacian ${\Delta }_{\mathcal{M}}$ and the normal derivative ${\partial }_{\mathcal{M},\mathrm{n}}$ on $\mathcal{M}$. It is assumed that $\mathcal{M}$ is compact (with respect to the topology induced by ${\mathbb{R}}^3$), which implies the existence of a countable subset $\{E_k:k=0,1,\dots \}$ of eigenfunctions of the negative Laplacian, $-{\Delta }_{\mathcal{M}}{E}_k={\lambda }_k{E}_k$, which can be taken as a basis for $L^2(\mathcal{M})$. From compactness it follows that $\mathcal{M}$ is bounded and it is further assumed that $\mathcal{M}$ has a boundary $\partial \mathcal{M}$ that is a piecewise smooth 1-dimensional manifold. (From a practical perspective, the assumption that the manifold has a boundary is not restrictive, since manifolds without boundary, such as the sphere, are likely to be challenging to realise as 3D-printed components.) The induced measure on $\partial \mathcal{M}$ is denoted ${\mu }_{\partial \mathcal{M}}$. In the sequel the manifold $\mathcal{M}$ is fixed and the shorthand $\Delta ={\Delta }_{\mathcal{M}}$ is used, leaving the dependence of the Laplacian on $\mathcal{M}$ implicit.

2.2.2 Two-Dimensional Gaussian Random Fields

Hu et al. (2013a,b) considered using coupled systems of SPDEs to model random vector fields, building on the constructions in Lindgren, Rue, and Lindström (2011). In this work, a similar coupled SPDE construction is used, defined next.

Let Ω be an underlying probability space, on which all random variables in this article are defined, and let $\omega \mapsto Z^{(r)}(x;\omega ),r\in \{1,2\}$ be independent, centred, real-valued (true or generalized) random fields, whose distribution is to be specified. Following standard convention, the argument $\omega \in \Omega$ is left implicit throughout. Let ${\mathcal{L}}^{(rs)},r,s\in \{1,2\}$, be differential operators on $\mathcal{M}$, to be specified, and consider then the coupled system of SPDEs(2) $$\left[\begin{array}{cc}{\mathcal{L}}^{(11)}& {\mathcal{L}}^{(12)}\\ {\mathcal{L}}^{(21)}& {\mathcal{L}}^{(22)}\end{array}\right]\left[\begin{array}{c}{u}^{(1)}(x)\\ u^{(2)}(x)\end{array}\right]=\left[\begin{array}{c}{Z}^{(1)}(x)\\ Z^{(2)}(x)\end{array}\right],\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}x\in \mathcal{M},$$(2) whose specification is completed with Neumann boundary conditions on $\partial \mathcal{M}$. The system in (2) is compactly represented in matrix form as $\mathcal{L}\mathbf{u}=\mathbf{Z}$. The distribution of u, the vector field solution of (2), is therefore fully determined by the specification of the matrix differential operator $\mathcal{L}$ and the distribution of the random field Z. In particular, the components $u^{(r)},r\in \{1,2\}$ are independent if ${\mathcal{L}}^{(12)}={\mathcal{L}}^{(21)}=0$ and if $Z^{(1)}$ and $Z^{(2)}$ are independent. Suitable choices for $\mathcal{L}$ and Z are discussed next.

2.3 Candidate Models

A complete statistical description for the random vector field u follows from specification of both the differential operator $\mathcal{L}$ and the noise process Z. Sections 2.3.1–2.3.4 present four candidate parametric models for $\mathcal{L}$, then Sections 2.3.5–2.3.7 present three candidate parametric models for Z. Thus, a total of 12 candidate parametric models are considered for the random vector field u of interest.

It is important to acknowledge that the SPDE framework does not currently enjoy the same level of flexibility for modeling of complex structural features as the covariance function framework. This is for two reasons; first, the indirect relationship between an SPDE and its associated random field presents a barrier to the identification of a suitable differential operator to represent a particular feature of interest and, second, each SPDE model requires the identification or development of a suitable numerical (e.g., finite element) method in order to facilitate computations based on the model. A discussion of the numerical methods used to discretise the SPDEs in this work is reserved for Appendix A.1, with model-specific details contained in Appendix A.2. For these two reasons, the aim in the sequel is not to arrive at a generative model whose realizations are indistinguishable from the real samples of the material. Rather, the aim is to arrive at a generative model whose realizations are somewhat realistic in terms of the salient statistical properties that determine performance in the engineering context. The models were therefore selected to exhibit a selection of physically relevant features, including anisotropy, oscillatory behavior, degrees of smoothness and nonstationarity of the random field. It is hoped that this article will help to catalyze collaborations with the surface metrology community that will lead in the future to richer and more flexible SPDE models for geometric variation in material.

The candidate models presented next are similar to those considered in Lindgren, Rue, and Lindström (2011), but generalized to coupled systems of SPDEs in a similar manner to Hu et al. (2013a,b). Note that formal model selection was not considered in Lindgren, Rue, and Lindström (2011) or Hu et al. (2013a,b).

2.3.1 Isotropic Stationary

For the differential operators ${\mathcal{L}}^{(rs)},r,s\in \{1,2\}$, the simplest candidate model that considered is(3) $${\mathcal{L}}^{(rs)}u(x):=({\eta }^{(rs)}-\Delta )({\tau }^{(rs)}u(x)),\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}x\in \mathcal{M},$$(3) where ${\eta }^{(rs)},{\tau }^{(rs)}>0$ are constants to be specified. This model is stationary, meaning that the differential operator ${\mathcal{L}}^{(rs)}$ does not depend on x, and isotropic, meaning that it is invariant to rotation of the coordinate system on ${\mathbb{R}}^3$. The differential order of the differential operator in (3) is fixed, which can be contrasted with Lindgren, Rue, and Lindström (2011) where powers ${(\eta -\Delta )}^{\alpha /2}$, $\alpha \in \mathbb{N}$, of the pseudo-differential operator were considered. The decision to fix α = 2 was taken because (a) the roughest setting of α = 1 was deemed not to be realistic for 3D-printed steel, since the solutions u to the SPDE are then not well-defined, and (b) over-parameterization will be avoided by controlling smoothness of the vector field u through the smoothness of the driving noise process Z, rather than through the differential order of $\mathcal{L}$. The model parameters in this case are therefore $\theta ={\{{\tau }^{(rs)},{\eta }^{(rs)}\}}_{r,s\in \{1,2\}}$. To reduce the degrees of freedom, it was assumed that ${\eta }^{(12)}={\eta }^{(21)}$, ${\eta }^{(11)}={\eta }^{(22)},{\tau }^{(12)}={\tau }^{(21)}$, and ${\tau }^{(11)}={\tau }^{(22)}$. This encodes an exchangeability assumption on the statistical properties of the upper and lower surfaces of the material. Such an assumption seems reasonable since the distinction between upper and lower surface was introduced only to improve the pedagogy; in reality such a distinction is arbitrary.

2.3.2 Anisotropic Stationary

The visible banding structure in Figure 2(b), formed by the sequential deposition of steel as the material is printed, suggests that an isotropic model is not well-suited to describe the dataset. To allow for the possibility that profiles parallel and perpendicular to the weld bands admit different statistical descriptions, the following anisotropic stationary model was considered:(4) $$\begin{array}{c}{\mathcal{L}}^{(rs)}u(x):=({\eta }^{(rs)}-\nabla ·{\mathbf{H}}^{(rs)}\nabla )({\tau }^{(rs)}u(x)),\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}x\in \mathcal{M}.\end{array}$$(4)

Here (4) differs from (3) through the inclusion of a diffusion tensor ${\mathbf{H}}^{(rs)}$, $r,s\in \{1,2\}$, a positive-definite 3 × 3 matrix whose elements determine the nature and extent of the anisotropy being modeled. The model parameters in this case are $\theta ={\{{\eta }^{(rs)},{\tau }^{(rs)},{\mathbf{H}}^{(rs)}\}}_{r,s\in \{1,2\}}$. To reduce the number of parameters, it was assumed that ${\mathbf{H}}^{(rs)}\equiv \mathbf{H}$ where H is a $3\times 3$ diagonal matrix, so that the nature of the anisotropy is the same for each of the four differential operators being modeled. Moreover, anisotropy was entertained only in the direction of printing (the z-axis of ${\mathbb{R}}^3$) relative to the nonprinting directions, meaning that the model is in fact orthotropic with $\mathbf{H}=\text{diag}(h_1,h_1,h_2)$ for as yet unspecified $h_1,h_2>0$. As before, an exchangeability assumption on the upper and lower surfaces of the material was enforced, so that ${\eta }^{(12)}={\eta }^{(21)},{\eta }^{(11)}={\eta }^{(22)},{\tau }^{(12)}={\tau }^{(21)}$, and ${\tau }^{(11)}={\tau }^{(22)}$.

2.3.3 Isotropic Nonstationary

A more detailed inspection of the material in Figure 2(b) reveals that (at least on the visible surface) the vertical centre of the panel exhibits greater variability in terms of surface height compared to the other portions of the panel. It can therefore be conjectured that the printing process is nonstationary with respect to the direction of the z-axis, in which case this nonstationarity should be encoded into the stochastic model. To entertain this possibility, a differential operator of the form$${\mathcal{L}}^{(rs)}u(x):=({\eta }^{(rs)}(x_3)-\Delta )({\tau }^{(rs)}(x_3)u(x)),\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}x\in \mathcal{M},$$was considered, which generalizes (3) by allowing the coefficients ${\eta }^{(rs)}$ and ${\tau }^{(rs)}$ to depend on the z-coordinate x₃ of the current location x on the manifold. Here, it was assumed that ${\tau }^{(rs)}(x_3)=c_{\tau }^{(rs)}\tau (x_3)$ and ${\eta }^{(rs)}(x_3)=c_{\eta }^{(rs)}\eta (x_3)$ for some unspecified $c_{\tau }^{(rs)},c_{\eta }^{(rs)}$ and some functions τ, η, to be specified, whose scale is fixed. The nature of the nonstationarity, as characterized by τ and η, will differ from printed component to printed component and this higher order variation could be captured by a hierarchical model that can be considered to have given rise to each specific instance of the $\tau ,\eta$.

Recent work due to Roininen et al. (2019) considered the construction of nonstationary random fields that amounted to using a nonparametric Gaussian random field for both τ and η. However, this is associated with a high computational overhead that prevents one from deploying such an approach at scale in the model selection context. Therefore, the functions τ and η were instead parameterized using a low-order Fourier basis, for τ taking(5) $$\tau (x)=1+\sum _{i=1}^4\sum _{j=0}^1{\gamma }_{\tau }^{(ij)}\hspace{0.17em}\text{sin}\hspace{0.17em}\left(\frac{2\pi i}{100}{x}_3+\frac{j\pi }{4}\right),$$(5) where the ${\gamma }_{\tau }^{(·)}$ are parameters to be specified. The units of x₃ are millimeters, so (5) represents low frequency nonstationarity in the differential operator. For η an identical construction was used, with coefficients denoted ${\gamma }_{\eta }^{(·)}$. The model parameters in this case are $\theta =\{c_{\tau }^{(rs)},c_{\eta }^{(rs)},{\gamma }_{\tau }^{(ij)},{\gamma }_{\eta }^{(ij)}\}$. To reduce the number of parameters an exchangeability assumption was again used, setting $c_{\tau }^{(11)}=c_{\tau }^{(22)},c_{\tau }^{(12)}=c_{\tau }^{(21)}$, $c_{\eta }^{(11)}=c_{\eta }^{(22)}$, and $c_{\eta }^{(12)}=c_{\eta }^{(21)}$.

To circumvent the prohibitive computational cost of fitting a hierarchical model for the Fourier coefficients ${\gamma }_{\tau }^{(·)}$ and ${\gamma }_{\eta }^{(·)}$, joint learning of these parameters across the panel datasets was not attempted. This can be justified by the fact that we have only 6 independent samples, drastically limiting the information available for any attempt to learn a hierarchical model from the dataset.

2.3.4 Anisotropic Nonstationary

The final model considered for the differential operator is the natural combination of the anisotropic model of Section 2.3.2 and the nonstationary model of Section 2.3.3. Specifically, the differential operator$$\begin{array}{cc}{\mathcal{L}}^{(rs)}u(x)\hfill & :=({\eta }^{(rs)}(x_3)-\nabla ·{\mathbf{H}}^{(rs)}\nabla )({\tau }^{(rs)}(x_3)u(x)),\hfill \\ \hfill & \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}x\in \mathcal{M}\hfill \end{array}$$was considered, where the 3 × 3 matrices ${\mathbf{H}}^{(rs)}$ and the functions ${\tau }^{(rs)},{\eta }^{(rs)}$ are parameterized in the same manner, with the same exchangeability assumptions, as earlier described. Note that the matrices ${\mathbf{H}}^{(rs)}$ are not considered to be spatially dependent in this model.

2.3.5 White Noise

Next, attention turns to the specification of a generative model for the noise process Z. In what follows three candidate models are considered, each arising in some way from the standard white-noise model. The simplest of these models is just the white-noise model itself, described first.

Let $Z^{(j)}:=W^{(r)}$ be a spatial random white-noise process on $\mathcal{M}$, defined as an $L^2(\mathcal{M})$-bounded generalized Gaussian field such that, for any test functions $\{{\phi }_i\in L^2(\mathcal{M}):i=1,\dots ,n\}$, the integrals ${〈{\phi }_i,W^{(j)}〉}_{\mathcal{M}},i=1,\dots ,n$, are jointly Gaussian with(6) $$\begin{array}{c}\mathbb{E}{〈{\phi }_i,W^{(j)}〉}_{\mathcal{M}}=0\\ \mathbb{C}({〈{\phi }_i,W^{(j)}〉}_{\mathcal{M}},{〈{\phi }_j,W^{(j)}〉}_{\mathcal{M}})={〈{\phi }_i,{\phi }_j〉}_{\mathcal{M}}.\end{array}$$(6)

It is not necessary to introduce a scale parameter for the noise process since the ${\tau }^{(rs)}$ parameters perform this role in the SPDE.

For this choice of noise model, the random field u of interest is typically continuous but not mean-square differentiable, which may or (more likely) may not be an appropriate mathematical description of the material.

2.3.6 Smoother Noise

It is difficult to determine an appropriate level of smoothness for the random field u based on visual inspection, except to rule out models that cannot be physically realized (see Section 2.3.1). For this reason, a smoother model for the noise process Z was also considered. In this case, the noise process is itself the solution of an SPDE(7) $$\begin{array}{c}({\eta }^{(r)}-\nabla ·\mathbf{H}\nabla )Z^{(r)}(x)=W^{(r)}(x)x\in \mathcal{M}\\ {\partial }_{\mathrm{n}}{Z}^{(r)}(x)=0x\in \partial \mathcal{M}\end{array}$$(7) where $W^{(r)}$ is the white-noise process defined in Section 2.3.5. To limit scope, it was assumed that the matrix H is the same as the matrix used (if indeed one is used) in the construction of the differential operator in Sections 2.3.2 and 2.3.4. The only free parameters are therefore ${\eta }^{(1)},{\eta }^{(2)}>0$ and the further exchangeability assumption was made that these two parameters are equal.

For this choice of noise model, the random field u of interest is typically twice differentable in the mean square sense, which may or may not be more realistic than the white-noise model.

2.3.7 Smooth Oscillatory Noise

The final noise model allows for the possibility of oscillations in the realizations of Z. As explained in Appendix C4 of Lindgren, Rue, and Lindström (2011), such noise can be formulated as arising from the complex SPDE$$\begin{array}{cc}(\eta e^{i\pi \theta }-\nabla ·\mathbf{H}\nabla )[Z^{(1)}(x)+i{Z}^{(2)}(x)]\hfill & =W^{(1)}(x)+i{W}^{(2)}(x)\hfill \\ \hfill & \hspace{1em}x\in \mathcal{M}\hfill \\ {\partial }_n[Z^{(1)}(x)+i{Z}^{(2)}(x)]\hfill & =0\hfill \\ \hfill & \hspace{1em}x\in \partial \mathcal{M}.\hfill \end{array}$$

Here, $\eta >0$ and $0\le \theta <1$ are to be specified, with θ = 0 corresponding to the smoother noise model of Section 2.3.6. The matrix H is again assumed to be equal to the corresponding matrix used (if indeed one is used) in the construction of the differential operator in Sections 2.3.2 and 2.3.4, and $W^{(1)}$ and $W^{(2)}$ are independent white-noise processes as defined in Section 2.3.5.

This completes the specification of candidate models for $\mathcal{L}$ and Z, constituting 12 different models for the random vector field u in total. These will be denoted M_i, $i=1,\dots ,12$ in the sequel. Section 3 seeks to identify which of these models represents the best description of the panel training dataset.

3 Data Preprocessing and Model Selection

This section proceeds as follows: In Section 3.1 the training dataset is described and Section 3.2 explains how it was pre-processed. Section 3.3 describes how formal statistical model selection was performed to identify the “best” model from the collection of 12 models just described. The fitting of model parameters raises substantial computational challenges, which this article does not fully solve but discusses in detail.

3.1 Experimental Set-Up

The data were obtained from large sheets of notional 3.5 mm thickness manufactured by the Dutch company MX3D (https://mx3d.com) using their proprietary printing protocol. The dataset consisted of six panels of approximate dimensions 300 mm × 300 mm, cut from a large sheet. These panels were digitized using a laser scanner to form a data point cloud; the experimental set-up and a typical dataset are displayed in Figure 2. The measuring process itself can be considered to be noiseless, due to the relatively high (0.1mm) resolution of the scanning equipment. For each panel a dense data point cloud was obtained, capturing the geometry of both sides of the panel.

3.2 Data Preprocessing

Before data were analyzed they were preprocessed and the steps taken are now briefly described.

3.2.1 Cropping of the Boundary

First, as is clear from Figure 2(b), part of the boundary of the panel demonstrates a high level of variation. This variation is not of particular interest, since in applications the edges of a component can easily be smoothed. Each panel in the dataset was therefore digitally cropped using a regular square window to ensure that boundary variation was removed. Note that no reference coordinate system was available and therefore the orientation of the square window needed to be determined (i.e., a registration task needed to be solved). To identify the direction perpendicular to the plane of the panel a principal component analysis of the data point cloud was used to determine the direction of least variation; this was taken to be the direction perpendicular to the plane of the panel. To orient the square window in the plane a coordinate system was chosen to maximize alignment with the direction of the weld bands. The weld bands are clearly visible (see Figure 2(b)) and the average intensity of a Gabor filter was used to determine when maximal agreement had been reached (see, e.g., Weldon, Higgins, and Dunn 1996).

3.2.2 Global De-Trending

As a weld cools, a contraction in volume of the material introduces residual stress in the large sheet (Radaj 2012). As a consequence, when panels are cut from the large sheet these stresses are no longer contained, resulting in a global curvature in each panel as internal forces are re-equilibrated. This phenomenon manifests as a slight curvature in the orthographic projections of the panel shown in Figure 2(b). The limited size of the dataset, in terms of the number of panels cut from the large sheet, precludes a detailed study of geometric variation due to residual stress in this work. Therefore residual stresses were removed from each panel dataset by subtracting a quadratic trend, fitted using a least-square objective in the direction perpendicular to the plane of the panel. This narrows focus to the variation that is associated with imprecision in the motion of the print head as steel is deposited. The data at this stage constitute a point cloud in ${\mathbb{R}}^3$ that is centred around a two-dimensional, approximately square manifold $\mathcal{M}$, aligned to the x and y axes in ${\mathbb{R}}^3$.

3.2.3 High-Frequency Filter

Closer inspection of Figure 2(b) reveals small bumps on the surface of the panel. These are referred to as “splatter” and result from small amounts of molten steel being sprayed from the weld head as subsequent layers of the sheet are printed. Splatter is not expected to contribute much to the performance of the material and therefore is of little interest. To circumvent the need to build a statistical model for splatter, a preprocessing step was used in order to remove such instances from the dataset. For this purpose a two-dimensional Fourier transform was applied independently to each surface of each panel and the high frequencies corresponding to the splatter were removed. Since data take the form of a point cloud, the data were first projected onto a regular Cartesian grid, with grid points denoted $x_i\in \mathcal{M}$, in order that a fast Fourier transform could be performed. The frequencies for removal were manually selected so as to compromise between removal of splatter and avoidance of smoothing of the relevant variation in the dataset. The resolution of the grid was taken to be 300 × 300 (i.e., millimeter scale), which was fine enough to capture relevant variation whilst also enabling the computational analyses described in the sequel.

3.2.4 Preprocessed Dataset

The result of preprocessing a panel dataset is a collection of triples$$(x_i,u_i^{(1)},u_i^{(2)})$$where $x_i\in {\mathcal{M}}_1$ represents a location in the manifold representation ${\mathcal{M}}_1$ of the panel, $u_i^{(1)}=u^{(1)}(x_i)\in \mathbb{R}$ represents the measured Euclidean distance of the upper surface of the panel from $x_i\in {\mathcal{M}}_1$ and $u_i^{(2)}=u^{(2)}(x_i)\in \mathbb{R}$ represents the measured Euclidean distance of the lower surface of the panel from $x_i\in {\mathcal{M}}_1$. The subscript on ${\mathcal{M}}_1$ indicates that the training data are associated to this manifold; a test manifold ${\mathcal{M}}_2$ will be introduced in Section 4.1. The latent two-dimensional vector field u describes the continuous geometry of the specific panel and it is the statistical nature of these vector fields which will be modeled. In the sequel a collection of candidate models for u are considered, each based on a Gaussian stochastic process.

The preprocessed data for a panel are collectively denoted ${\mathcal{D}}^{(j)}$, where $j\in \{1,\dots ,6\}$ is the index of the panel. The collection of all 6 of these panel datasets is denoted ${\mathcal{D}}_1$, the subscript indicating that these constitute the training dataset.

3.3 Model Selection

To perform formal model selection among the candidate models an information criterion was used, specifically the Akaike information criterion (AIC)(8) $$\text{AIC}(M)=2\text{dim}({\widehat{\theta }}_M)-2\hspace{0.17em}\text{log}\hspace{0.17em}p({\mathcal{D}}_1|{\widehat{\theta }}_M,M),$$(8) where $p({\mathcal{D}}_1|{\theta }_M,M)$ denotes the likelihood associated with the training dataset ${\mathcal{D}}_1$ under the SPDE model M, whose parameters are θ_M, and ${\widehat{\theta }}_M$ denotes a value of the parameter θ_M for which the likelihood is maximized. Small values of AIC are preferred. The parameter vector θ_M can be decomposed into parameters that are identical for each panel dataset ${\mathcal{D}}^{(j)}$, denoted ${\theta }_M^{\text{id}}$, and those that are specific to a panel dataset ${\mathcal{D}}^{(j)}$, denoted ${\theta }_M^{(j)}$. The latter are essentially random effects and consist of the Fourier coefficients ${\gamma }_{\tau }^{(·)}$ and ${\gamma }_{\eta }^{(·)}$ in the case where M is a nonstationary model (Sections 2.3.3 and 2.3.4), otherwise ${\theta }_M^{(j)}$ is empty and no random effects are present. The panel datasets ${\mathcal{D}}^{(j)}$ that together comprise the training dataset ${\mathcal{D}}_1$ are considered to be independent given θ_M and M, so that(9) $$p({\mathcal{D}}_1|{\theta }_M,M)=\prod _{j=1}^{6}p({\mathcal{D}}^{(j)}|{\theta }_M^{\text{id}},{\theta }_M^{(j)},M).$$(9)

The maximization of (9) to obtain ${\widehat{\theta }}_M$ requires a global optimization routine and this will in general preclude the exact evaluation of the AIC for the statistical models considered. Following standard practice, ${\widehat{\theta }}_M$ was therefore set to be the output of a numerical optimization procedure applied to maximization of (9). A drawback of the SPDE approach is that implementation of an effective numerical optimization procedure is difficult due to the considerable computational cost involved in differentiating the likelihood. To obtain numerical results, this work arrived, by trial-and-error, at a procedure that contains the following ingredients:

Gradient ascent: In the case where M is a stationary model, so the ${\theta }_M^{(j)}$ are empty, the maximum likelihood parameters ${\widehat{\theta }}_M^{\text{id}}$ were approximated using the natural gradient ascent method. This extracts first order gradient information from (9) and uses the Fisher information matrix as a surrogate for second-order gradient information, running a quasi-Newton method that returns a local maxima of the likelihood. Full details are provided in Appendix A.4.

Surrogate likelihood: The use of gradient information constitutes a computational bottleneck, due to the large dense matrices involved, which motivates the use of a surrogate likelihood within the gradient ascent procedure whose Fisher information is more easily computed (Tajbakhsh, Aybat, and Del Castillo 2018). For this article a surrogate likelihood was constructed based on a subset of the dataset, reduced from the 300 × 300 grid to a 50 × 50 grid in the central portion of the panel. Natural gradient ascent was feasible for the surrogate likelihood. It is hoped that minimization of this surrogate likelihood leads to a suitable value ${\widehat{\theta }}_M$ with respect to maximizing the exact likelihood, and indeed the values that reported in the main text are evaluations of the exact likelihood at ${\widehat{\theta }}_M$.

Two-stage procedure for random effects: The case where M is nonstationary is challenging for the natural gradient ascent method, even with the surrogate likelihood, due to the much larger dimension of the parameter vector when the random effects ${\theta }_M^{(j)}$ are included. In this case a two-stage optimization approach was used, wherein the common parameters ${\theta }_M^{\text{id}}$ were initially fixed equal to their values under the corresponding stationary model and then optimized over the random effects ${\theta }_M^{(j)}$ using natural gradient ascent, which can be executed in parallel.

These approximation techniques combined to produce the most robust and reproducible parameter estimates that the authors were able to achieve, but it is acknowledged that they do not guarantee finding a local maximum of the likelihood and this analytic/numerical distinction can lead to pathologies that are discussed next. The computations that we report were performed in Matlab R2019b and required approximately one month of CPU time on a Windows 10 laptop with an Intel® Core^TMi7-6600U CPU and 16GB RAM. The number of iterations of natural gradient ascent required to achieve convergence ranged from 10 to 80 and was model-dependent. Models with $\ge 100$ parameters required several days of computation, while the simpler models with < 10 parameters typically required less than one day of computation in total.

The twelve candidate models are presented in the first three columns of Table 2, together with the dimension $\text{dim}({\theta }_M)$ of their parameter vector, the (approximately) maximized value of the surrogate and exact log-likelihood and the AIC. The main conclusions to be drawn are as follows:

Table 2 Model selection based on the training dataset ${\mathcal{D}}_1$.

Isotropic	Stationary	Noise	Log-likelihood
			# Parameters	Surrogate	Exact	AIC
✓	✓	White	4	21,187	1,041,766	–2,083,524
✓	✓	Smoother	5	21,775	1,154,950	–2,309,890
✓	✓	sm. oscil.	6	22,752	1,208,810	–2,417,608
✓	✗	White	100	24,404	774,539	–1,548,878
✓	✗	Smoother	101	25,900	1,023,561	–2,046,920
✓	✗	sm. oscil.	102	29,118	1,042,737	–2,085,270
✗	✓	White	6	22,005	1,100,358	–2,200,704
✗	✓	Smoother	7	23,166	1,164,952	–2,329,890
✗	✓	sm. oscil.	8	25,151	1,233,936	–2,467,856
✗	✗	White	102	25,385	914,913	–1,829,622
✗	✗	Smoother	103	26,336	1,056,735	–2,113,264
✗	✗	sm. oscil.	104	28,150	1,111,449	–2,222,690

NOTE: The first three columns indicate whether the differential operator $\mathcal{L}$ was isotropic and stationary and the driving noise process Z that was used. The remaining columns indicate the number of free parameters in each model, the approximately maximized values of the surrogate and exact log-likelihood and the Akaike information criterion (AIC). [The log-likelihood and AIC values are rounded to the nearest integer.]

Isotropy: In all cases, the anisotropic models for $\mathcal{L}$ out-performed the corresponding isotropic models according to both the (approximately) maximized exact log-likelihood and the AIC. (The same conclusion held for 5 of 6 comparisons based on the surrogate likelihood.) Since the anisotropic model includes the isotropic model as a special case, these higher values of the likelihood are to be expected, but the AIC conclusion is nontrivial. (Inspection of the log-likelihood does, however, provide a useful validation of the optimization procedure that was used.)

Stationarity: In all cases, the nonstationary models for $\mathcal{L}$ out-performed the corresponding stationary models according to the (approximately) maximized surrogate log-likelihood, which is to be expected since these models are nested. However, the conclusion was reversed when the parameters ${\widehat{\theta }}_M$, estimated using the surrogate likelihood, were used to evaluate the exact likelihood. These results suggest that the two-stage numerical optimization procedure for estimation of random effects was robust but that the surrogate likelihood was lossy in respect of information in the dataset that would be needed to properly constrain parameters corresponding to nonstationarity in the model.

Noise model: In all cases, the smooth oscillatory noise models for Z led to substantially lower values of the surrogate and the exact log-likelihood compared to the white-noise and smooth noise models. This is to be expected since these models are nested, but the fact that this ordering was recovered in both the values of the surrogate and exact log-likelihood indicates that the numerical procedure worked well in respect of fitting parameters pertinent to the noise model.

Best model: The best-performing model, denoted $M^{*}$ in the sequel, was defined as the model that provided the smallest value of AIC, computed using the (approximately) maximized exact likelihood. This was the anisotropic stationary model $\mathcal{L}$ with smooth oscillatory noise Z, which had 8 parameters and achieved a substantially lower value of AIC compared to each of the 11 candidate models considered. Although the use of AIC in this context is somewhat arbitrary and other information criteria could also have been used, the substantially larger value of the log-likelihood for the best-performing model suggests that use of an alternative information criterion is unlikely to affect the definition of the best-performing model.

If instead the surrogate likelihood had been used to produce values for AIC, it would have declared that the best-performing model also includes nonstationary behavior. This suggests that the aforementioned “best” model is more accurately described as the “best model that could be fitted” and that, were the computational challenges associated with parameter estimation in the SPDE framework surmounted, a model that involves nonstationarity may have been selected. This suggests that additional models of further complexity may provide better descriptions of the data, at least according to AIC, but that fitting of model parameters presents a practical barrier to the set of candidate models that can be considered. The remainder investigates whether predictions made using the simple model $M^{*}$ may still be useful.

4 Transfer Learning

This section aims to assess the capacity of our selected statistical model $M^{*}$ to perform transfer learning, and proceeds as follows: Section 4.1 assesses the performance of the statistical model in the context of a 3D-printed cylinder, which constitutes a held-out test dataset. Then, Section 4.2 implements finite element simulations of compressive testing based on the fitted model and compares these against the results of a real compressive experiment.

4.1 Predictive Performance Assessment

To investigate the generalization performance of the model $M^{*}$ identified in Section 3.3, a test dataset, denoted ${\mathcal{D}}_2$ and displayed in Figure 4(a), was obtained as a laser scan of a 3D-printed cylinder of diameter $\approx$170 mm and length $\approx$581 mm. The notional thickness of the cylinder was 3.5 mm, in agreement with the training dataset. The manifold representation of the panel is denoted ${\mathcal{M}}_1$ and the manifold representation of the cylinder is denoted ${\mathcal{M}}_2$. The fitted model $M^{*}$ induces a bivariate Gaussian field on any suitably regular manifold, in particular ${\mathcal{M}}_2$. A sample from this model is displayed in Figure 4(b). At a visual level, the “extent” of the variability is similar between Figure 4(a) and Figure 4(b), but the experimentally produced cylinder contains more prominent banding structure compared to the cylinder simulated from $M^{*}$.

Fig. 4 (a) Laser scan of a section from a real 3D-printed cylinder. (b) Sample from the fitted model $M^{*}$ generated based on a cylindrical manifold. (c) The distribution over the thickness of the wall of the real (blue, shaded) and the simulated (red, solid) 3D-printed cylinder.

The thickness of the cylinder wall is relevant to its mechanical performance when under load. Recall that $M^{*}$ was not explicitly trained as a predictive model for thickness; rather, thickness is a derived quantity of interest. It is therefore interesting to investigate whether the distribution of material thickness, accouring to $M^{*}$ instantiated on ${\mathcal{M}}_2$, agrees with the experimentally obtained test dataset. Results in Figure 4(c) indicate that the fitted model agrees with the true distribution of material thickness in as far as the modal value and the right tail are well-approximated. However, the true distribution is positively skewed and the model $M^{*}$, under which thickness is necessarily distributed as a Gaussian random field, is not capable of simultaneously approximating also the left tail. From an engineering perspective, the model $M^{*}$ is conservative in the sense that it tends to under-predict wall thickness, in a similar manner to how factors of safety might be employed.

The predictive likelihood of a held-out test dataset is often used to assess the predictive performance of a statistical model. However, there are two reasons why such assessment is unsuitable in the present context: First, the preprocessing of the training dataset, described in Section 3.2, was adapted to the manifold ${\mathcal{M}}_1$ of the training dataset and cannot be directly applied to the manifold ${\mathcal{M}}_2$ of the test dataset. This precludes a fair comparison, since the fitted statistical model is not able to explain the high-frequency detail (such as weld splatter, visible in Figure 4(a)) that are present in the test dataset ${\mathcal{D}}_2$ but were not included in the training dataset. Second, recall that the purpose of the statistical model is limited to capturing aspects of geometric variation that are consequential in engineering. The predictive likelihood is not well-equipped to determine if the statistical model is suitable for use in the engineering context.

For a more detailed assessment of the statistical model and its predictive performance in the engineering context, finite element simulations of cylinders under load were performed, described next.

4.2 Application to Compressive Testing

As discussed in Section 2.3, the aim of this work is to model those aspects of variability that are salient to the performance of components manufactured from 3D-printed steel. To this end, predictions were generated for the outcome of a compressive test of a cylinder. This test, which was also experimentally performed, involved a cylinder of the same dimensions described in Section 4.1, which was machined to ensure that its ends are parallel to fit the testing rig. The test rig consisted of a pair of end plates which compress the column ends axially. The cylinder was placed under load until a set amount of displacement has been reached. Both strain gauge and digital image correlation measurements of displacement under a given load were taken; further details may be found in Buchanan, Wan, and Gardner (2018).

To produce a predictive distribution over possible outcomes of the test, three geometries were simulated for the cylinder using $M^{*}$ and finite element simulations of the compressive test were performed, using these geometries in the construction of the finite element model. To facilitate these simulations a simple homogeneous model for the mechanical properties of the material was used, whose values were informed by the coupon testing undertaken in Buchanan, Wan, and Gardner (2018). This is not intended as a realistic model for the mechanical properties of 3D-printed steel. Indeed, as discussed in Section 1, it is expected that 3D-printed steel exhibits local variation in mechanical properties such as stiffness; the focus of this article is surface metrology (only) and joint modeling of geometric and mechanical variation for 3D-printed steel will be addressed in future work.

Three typical end points from the finite element simulation are displayed in Figure 5, panels (a)-(c). Notice that local variation in geometry leads to nontrivial global variation in the location of the buckling point. Panel (d) in Figure 5 displays the buckled states of two real 3D-printed steel cylinders, experimentally obtained. At a visual level there is some qualitative agreement in the appearence of these failure modes. The approach taken in this article can be contrasted with approaches in the engineering literature that do not generalize across manifolds and do not report full probability distributions (Wagner, Hühne, and Elishakoff 2020). However, an important caveat is that the precise quantitative relationship between load and displacement in these simulations is not yet meaningful due to the simplistic model used for material properties. Further research will be required to incorporate these aspects into a complete statistical model for 3D-printed steel.

Fig. 5 Predicting the outcome of a compressive test on a 3D-printed cylinder, using the fitted generative model $M^{*}$ for the geometry of the cylinder and a finite element simulation of the experiment. (a)–(c) Samples from the generative model for the geometry of the cylinder after the simulated compressive test. (d) and (e) Buckled cylinders, experimentally produced. (Reproduced with permission from Buchanan, Wan, and Gardner 2018).

5 Discussion

This article identified surface metrology as a promising new application area for recent advances in spatial statistics. Though attention was focussed on 3D-printed stainless steel, the techniques and methods have the potential for application to a wide variety of problems in this engineering subfield.

Our approach sought to decouple the local statistical properties of the geometric variation from the global shape of the component, which offered the potential for transfer learning; that is, conclusions drawn from one component to be carried across to a possibly different component. This assumption is intuitively reasonable for components that do not involve high levels of curvature, but its appropriateness for general components remains to be tested. To achieve this we exploited the SPDE framework. However, these was a practical limitation on the complexity of SPDE that could be reliably fit, since gradients of the likelihood required large dense matrices to be computed. As a result, it is possible that models that are too simple but whose parameters can be successfully optimized may be preferred to more parsimonious models whose parameters cannot easily be optimized. In other applications of surface metrology where generalization across manifolds is not required, model fitting and model selection could be more easily achieved using the traditional covariance function framework (e.g., Rasmussen and Williams 2006; Duvenaud 2014; Băzăvan, Li, and Sminchisescu 2012; Oliva et al. 2016; Moeller et al. 2016).

The principal mathematical assumption that made transfer learning possible in this work was the identification of Laplacians across different manifolds, corresponding to the different components that can in principle be manufactured. This assumption (together with knowledge of the printing direction) enabled the meaning of a statistical model to be unambiguously interpreted across different manifolds, despite the fact that in general there are an infinitude of different geometries with which a given manifold can be equipped. Gaussian process models arise naturally in this context, since in specific cases they can be conveniently expressed in terms of the Laplacian via the SPDE framework. Alternative statistical models could also be considered: For example, the eigenfunctions $\{E_k:k=0,1,\dots \}$ of the negative Laplacian $-{\Delta }_{\mathcal{M}}$ form a basis for $L^2(\mathcal{M})$, and one can consider more general generative models for the random fields $u_1,u_2\in L^2(\mathcal{M})$ based on linear combinations of the E_k with non-Gaussian coefficients (see, e.g., Hosseini and Nigam 2017; Sullivan 2017). However, computational methods for non-Gaussian models remain relatively under-developed at present.

The scope of this investigation was limited to describing the geometric variation of 3D-printed steel and there are many interesting questions still to be addressed. These include statistical modeling of the mechanical properties of the material, optimization of the printing protocol to maximize the performance of (or reduce the variability in performance of) a given component, and adaptive monitoring and control of the manufacturing process to improve the quality of the printed component. These important problems will need to be collaboratively addressed; see Dodwell et al. (2021). The same printing protocol that we studied was recently used to construct an entire pedestrian bridge using 3D-printed steel (Zastrow 2020). To ensure the safety of the bridge for pedestrians, extensive engineering assessment and testing were required. A long-term goal of this research is to enable such assessment to be partly automated in silico before the component is printed. This would in turn reduce the design and testing costs associated with components manufactured with 3D-printed steel.

Acknowledgments

The authors thank the editor, associate editor, and two reviewers for their valuable comments on the article.

Supplementary Materials

The electronic supplement contains a detailed description of how computation was performed, together with the Matlab R2019b code that was used.

Engineering and Physical Sciences Research Council;

The authors are grateful to Craig Buchanan, Pinelopi Kyvelou and Leroy Gardner for providing the laser-scan datasets and to MX3D, from whom the material samples were provided. In addition, the authors are grateful for discussions with Alessandro Barp, Tim Dodwell, Mark Girolami, Peter Gosling and Rob Scheichl.

Additional information

Funding

This research was supported by the Lloyd’s Register Foundation programme on data-centric engineering and the Alan Turing Institute (under the EPSRC grant EP/N510129/1). WSK gratefully acknowledges support from the UK EPSRC (grant EP/K031066/1).