Strain Monitoring Strategy of Deformed Membrane Cover Using Unmanned Aerial Vehicle-Assisted 3D Photogrammetry

Large structures and high-value assets require inspection and integrity assessment methodologies that ensure maximum availability and operational capabilities. Large membranes are used as floating covers at the anaerobic wastewater lagoons of Melbourne Water’s Western Treatment Plant (WTP). A critical function of this high-value asset pertains to the harnessing of the biogas gas generated at these lagoons as well as protecting the environment from the release of odours and greenhouse gases. Therefore, a proactive inspection and efficient management strategy are required to ensure these expensive covers’ integrity and continued operation. Not only is identifying the state of stress on the floating cover crucial for its structural integrity assessment, but the development of rapid and non-contact inspections will significantly assist in determining the “real-life” performance of the cover for superior maintenance management. This study investigates a strain determination method for WTP floating cover which integrates unmanned aerial vehicle (UAV)-assisted photogrammetry with finite element analyses to determine the structural integrity of these covers. Collective aerial images were compiled to form 3D digital models of the deformed cover specimens, which were then employed in computational and statistical analyses to assess and predict the strain of the cover. The findings complement the future implementation of UAV-assisted aerial photogrammetry for structural health assessment of the large floating covers.


Introduction
Membrane covers are commonly used as floating covers for clean water reservoirs to prevent evaporation and pollution; for landfills to trap hazardous chemicals and unpleasant odour; and for tailing impoundment [1][2][3]. Floating cover materials such as high-density polyethylene (HDPE) geomembranes are very durable, resistant to many different solvents, and have a high strength-to-density ratio, and if well-designed, long service life (in decades) even in harsh environments [2,3]. The Western Treatment Plant (WTP) at Werribee, Victoria, Australia, installed two floating covers to assist with the anaerobic treatment of the raw sewage beneath them, leading to the production of methane-rich biogas, see Figure 1. The WTP floating covers are made from a 2 mm thick HDPE spanning an area of 450 m by 170 m. The WTP cover is held around its perimeter by spanning an area of 450 m by 170 m. The WTP cover is held around its perimeter by clamping strips and mechanical fasteners to provide an airtight seal. In addition to controlling odour, these covers are capable of collecting up to 65,000 m 3 of biogas per day that can generate up to 7 MW of electricity, thereby increasing the treatment plant's supply of renewable electricity from this sustainable source. Without this cover, the methane-rich biogas would be released into the environment as a damaging greenhouse gas. All sewage inflow is unscreened and passes first through an anaerobic reactor where biogas is produced as the raw sewage undergoes anaerobic digestion, which is trapped below the floating cover and harvested for electricity generation. Solidified sewage substance can amass on the reactor surface to form a scum-berg which presses against and elevates the covers to a height of approximately 1 metre. The scum-berg can cause lateral displacement, thus causing changes to the stress on the cover, which is partially due to hydrodynamic effects resulting from the sewage inflow. The movement of sewages and the presence of methane-rich biogas require that the proposed inspection methodology is intrinsically safe and, therefore, underpins the requirement of a proactive and effective inspection approach that ensures the integrity and continued operation of this asset. However, for a structure of this size, the definition of its state of stress is currently beyond the capability of current measurement techniques and, hence, requires a more robust approach. Over the past decades, the acquisition of full-field strain measurement has been an active field for structural health monitoring (SHM). The well-established non-contact strain measurements, such as digital image correlation (DIC), have been extensively investigated and have achieved a high degree of accuracy, primarily for in-plane strain measurement for various applications. However, in using the DIC method at large deformation, the uncertainty in strain measurements increases significantly as the speckle pattern degrades, and decorrelation occurs at such high strain levels [4][5][6][7][8]. A more precise full-field strain determination could be achieved by a 3D DIC method for nonplanar specimens with out-of-plane displacement and rotation using two synchronised cameras, and recently, one single camera [9][10][11][12]. However, this method is often limited by the fixed and small depth of field and highly controlled testing environment [9,13]. Nevertheless, it is known that the ability to accurately predict the strain field from 3D deformation is significantly challenging, especially for large structures in its operating environment. Recently, photogrammetry has become one of the commonly available non-contact and full-field measurement methods for SHM [14,15]. This method, with the aid of photogrammetric software, constructs a digital elevation model (DEM) from sequential 2D digital photographs. Baqersad et al. [16][17][18][19] integrated photogrammetry to predict the 3D dynamic strain experienced by the blade of a wind-turbine using finite element analysis. Pappa et al. [20] used optical techniques with reflective markers and projected dots to define the dynamic motion of membranes. Recently, Luo et al. [21] have investigated the concept of using Over the past decades, the acquisition of full-field strain measurement has been an active field for structural health monitoring (SHM). The well-established non-contact strain measurements, such as digital image correlation (DIC), have been extensively investigated and have achieved a high degree of accuracy, primarily for in-plane strain measurement for various applications. However, in using the DIC method at large deformation, the uncertainty in strain measurements increases significantly as the speckle pattern degrades, and decorrelation occurs at such high strain levels [4][5][6][7][8]. A more precise full-field strain determination could be achieved by a 3D DIC method for non-planar specimens with out-of-plane displacement and rotation using two synchronised cameras, and recently, one single camera [9][10][11][12]. However, this method is often limited by the fixed and small depth of field and highly controlled testing environment [9,13]. Nevertheless, it is known that the ability to accurately predict the strain field from 3D deformation is significantly challenging, especially for large structures in its operating environment. Recently, photogrammetry has become one of the commonly available non-contact and full-field measurement methods for SHM [14,15]. This method, with the aid of photogrammetric software, constructs a digital elevation model (DEM) from sequential 2D digital photographs. Baqersad et al. [16][17][18][19] integrated photogrammetry to predict the 3D dynamic strain experienced by the blade of a wind-turbine using finite element analysis. Pappa et al. [20] used optical techniques with reflective markers and projected dots to define the dynamic motion of membranes. Recently, Luo et al. [21] have investigated the concept of using photogrammetry to measure strain fields of deformed large inflatable structures by combining Delaunay triangulation and finite element methods. However, quantitative evaluations, such as strain and stress measurements, based on photogrammetry are still in their infancy, where many previous studies relied on the monitoring of surface fiducial markers, primarily focused on shape measurement of structural surfaces [21][22][23][24][25][26][27].
Unmanned aerial vehicle (UAV) is a generic uninhabited aerial vehicle that is remotely controlled either manually or autonomously. The ability to install non-contact measurement capabilities on unmanned aerial vehicles (UAV) is creating more rapid and safer opportunities in SHM of large structures such as aircraft wings, dams, vertical tanks, and elevated highways. Recently, there has been significant interest in using UAV systems coupled with photogrammetry to rapidly acquire high temporal and spatial resolution image information for 3D modelling and various disciplines [28][29][30][31][32][33]. This UAV photogrammetry system provides a low-cost commercially available or customisable practical alternative in the close-range aerial domain. This system can be deployed in high risk and inaccessible areas, where it is considered too difficult or impossible for traditional methods and without endangering human life or safety [29]. Furthermore, this would lead to a more time-efficient inspection process for large areas which include monitoring archaeological sites, environmental and agricultural areas, and traffic [29,31]. Despite the advantages in UAV photogrammetry, it is limited by its relatively short flight time and long computational processing due to a significantly large dataset.
In this collaborative project, a single camera mounted on a UAV-based photogrammetry system is deployed as a wide coverage non-contact measurement diagnostic tool to monitor the state of deformation of the WTP membrane covers. The research project proposed to build upon the practical diagnostic tool towards an innovative smarter structure and maintenance, reflecting the digital twin paradigm [34]. One of the research project objectives is to devise a two-stage global-local monitoring strategy where UAV photogrammetry monitoring is first deployed to rapidly evaluate the global strain response of the cover, thereby identifying critical areas where, subsequentially, a localised-detailed inspection is conducted for further quantitative assessment. This two-stage strategy aims to significantly reduce time and cost for inspection of the entire cover and mitigates workers' exposure to a high-risk environment. With the current digital deformation and fiducial features on the WTP membrane cover, it is highly advantageous to acquire the strain field of the membrane for a more meaningful measurement of the asset. However, the floating covers DEMs are highly subjected to artefacts from debris, lighting, and weather. Nevertheless, strain fields are extremely difficult to obtain from small measurements difference and as it is derived from displacement information, it requires highly accurate displacement readings [26]. Thus, the raw DEM requires significant preprocessing prior to further analysis. Furthermore, the WTP regulations, such as flight restriction (minimum and maximum height and duration) and the fact that coating or physical attachment on the cover are not feasible options, limit the ability to enhance the accuracy of acquired measurements through other means, making this problem-specific application significantly challenging.
Our previous works have shown the application of UAV photogrammetry systems into finite element (FE) analysis for strain determination and, furthermore, on-the-field work on the optimisation of flight parameters, which includes overlapping of photos, flight path, and altitude [28,[35][36][37]. In this paper, experimental studies are conducted where UAV photogrammetry is deployed to construct the DEMs of two deformed membrane specimens in a controlled laboratory setting, as a continuation of our investigations on reliable SHM of the WTP floating covers. The DEMs are then preprocessed using a smoothing interpolation with different smoothing parameters and applied as displacement loads for FE analysis to predict its displacement and strain field (see schematic in Figure 2). As uncertainty and erroneous measurements inevitably arise from all stages of the procedure (data collection and preprocessing of data), a probabilistic approach is considered to provide a statistical prediction of the strain field and its degree of certainty. Hence, the resulting smoothed strain fields are then utilised as training data for a variational inference method approach for heteroscedastic modelling, based on Gaussian process (GP) regression. The findings give insights into the development of a strain into the development of a strain monitoring strategy tailored for WTP floating covers and in aiding the decision support for WTP floating cover management.

Figure 2.
Schematic depiction of the procedure used for strain determination in using a UAV-assisted photogrammetry DEM for deformed membrane covers.

UAV Photogrammetry Setup and DEM
An unmanned aerial vehicle DJI Spark with an integrated camera was deployed for all the laboratory studies. The 3D DEMs of the membrane are developed by sequential 2D images from the mounted digital camera with pre-determined flight paths. The specification provided by the manufacturer is listed in Table 1.
Metashape Professional by Agisoft [38] was used to perform image alignment, construction of dense clouds, mesh, and DEM Agisoft adopted the computer vision algorithms (i.e. Structure from Motion 'SfM' and bundle adjustment), described in Barazzetti [39], Westoby [40], and Triggs [41], to develop an autonomous adjustment, realignment of multi-images, and reconstruction of a 3D model with high redundancy. Dominic [42] also reviewed the basic principles of the modern photogrammetric and procedures of conducting a photogrammetry approach. All the captured images were loaded with their corresponding metadata (e.g., GPS location and camera setting) to Metashape Professional for generating the DEM. All parameters (including alignment, dense cloud, and mesh generation) were set to the highest available setting to construct accurate DEMs.

UAV Photogrammetry Setup and DEM
An unmanned aerial vehicle DJI Spark with an integrated camera was deployed for all the laboratory studies. The 3D DEMs of the membrane are developed by sequential 2D images from the mounted digital camera with pre-determined flight paths. The specification provided by the manufacturer is listed in Table 1.
Metashape Professional by Agisoft [38] was used to perform image alignment, construction of dense clouds, mesh, and DEM Agisoft adopted the computer vision algorithms (i.e. Structure from Motion 'SfM' and bundle adjustment), described in Barazzetti [39], Westoby [40], and Triggs [41], to develop an autonomous adjustment, realignment of multi-images, and reconstruction of a 3D model with high redundancy. Dominic [42] also reviewed the basic principles of the modern photogrammetric and procedures of conducting a photogrammetry approach. All the captured images were loaded with their corresponding metadata (e.g., GPS location and camera setting) to Metashape Professional for generating the DEM. All parameters (including alignment, dense cloud, and mesh generation) were set to the highest available setting to construct accurate DEMs.

Noise Filtering of DEM
All DEMs contain random errors and artefacts and their accuracy depends on the topography, DEM generation method, and resolution [44][45][46][47]. Maps and models derived from raw noisy DEM are usually readable after some filtering postprocessing of data to filter noise and improve the data quality. It is important to obtain a high-quality DEM with minimal noise in order to accurately obtain the strain field. Common techniques such as attenuating high-frequency noise in DEM can be performed by spatial filtering based (low-pass filter) on a 2D Fourier transform [44]. Other techniques such as 2D discrete wavelet transform, cutting methods, and smoothing can be found in previous literature [46,48]. The DEM smoothing technique, which includes moving average smoothing, median smoothing, etc., is one of the most popular approaches for DEM denoising and widely used in DEM-based geology studies [45,49]. One example is the moving average method, which smooths data, where the size of a moving window and the number of smoothing iterations are manually selected. However, this approach has been used with caution as it may produce undesirable high-frequency artefacts. A median smoothing filter, which performs through the signal entry by entry, replacing each entry with the median of the neighbour entry, is also a well-known nonlinear noise filtering process to filter impulsive-like (speckle) noise. However, an additional scheme is required when processing entries at or near the boundary. Cubic smoothing splines interpolation is an alternative approach that embodies a surface fitting technique to create a smooth model of a complex profile while removing noisy data and avoiding Runge's phenomenon. In this study, the cubic smoothing splines technique is used to preprocess the raw photogrammetry DEM.
In general, the fitted curve spanning each data interval is presented by a cubic polynomial with the endpoints of the adjacent polynomials matching in location and the first and second derivative [50][51][52]. In its basic form, the objective is to minimise the square error and the curvature; refer to Equation (1) for the 2-dimensional cubic smoothing splines form.
where x i , y i , and z i are the data coordinates and f (x, y) are cubic polynomials, and the weighting parameter p balances the two countervailing constraints and results in a smoothing cubic spline. Succinctly, if p = 1, the smoothing spline passes through all data points resulting in an interpolation spline, and if p = 0, only the curvature of the spline is minimised and hence results in a linear least-squares fit. It should be noted that there are other modified versions of the cubic spline smoothing which include piecewise constant weight function in the curvature measurement and error measure weights. However, for simplicity, only the basic form where p is constant for all directions is considered in this current work.

Material and Experimental Setups
In this study, HDPE materials from Melbourne Water's supply were used. A preliminary investigation of a static pull test using static INSTRON 33R 424 was conducted on six HDPE specimens to validate the material properties, Young's modulus of 125 MPa, Poisson's ratio of 0.42 and stress and strain curve. Two HDPE materials of 2 mm thickness with different geometry sizes were used for two experimental investigations to represent two common deformations of the membrane cover: out-of-plane (OOP) deformation and in-plane deformation (wrinkle formation). The surface of each test specimen is marked with an internal grid with the discretisation of 20 mm to represent the target features. The test specimens are mounted on a test frame rig, refer to Figure 3a, and reference points are marked on the frame to align the DEMs to a 3D coordinate system. Rectangular wood blocks sandwich the edge of the specimens and are then secured using multiple screws and clamps to provide the uniform constraints on the edges for the two tests. For each photogrammetry analysis, 18 aerial In the OOP deformation test, a 1200 mm squared HDPE with an internal 1000 mm squared grid marked on the surface is used, as shown in Figure 4. The membrane cover is mounted on a partially fixed frame, allowing in-plane deformation, and dominantly displaced in the OOP direction by a hydraulic car jack located approximately at the centre of the membrane. The membrane is loaded six times, at first incrementally to a height of 23 mm, 48 mm, and 72 mm, then released to 63 mm and 51 mm, and a final loading to a height of 86 mm, which formed the six OOP deformation tests denoted as Step 1 to 6, respectively, for this study. Optical fibres were installed and aligned with the vertical lines of 1000 mm at 800 mm and 1000 mm horizontally on the HDPE specimen for benchmarking the strains; refer to Figure 4. In the in-plane deformation test, a 1200 mm by 500 mm rectangular HDPE with an internal 1000 mm by 300 mm grid is used, as shown in Figure 5. The membrane cover is fixed on one end while the other end is translated in the in-plane direction towards the fixed edge by 112 mm, 161 mm, and 267 mm. These imposed displacements caused the membrane to form a wrinkle/fold denoted as Step 1 to 3, respectively. The maximum fold heights are 190 mm, 223 mm, and 280 mm for each subsequent In the OOP deformation test, a 1200 mm squared HDPE with an internal 1000 mm squared grid marked on the surface is used, as shown in Figure 4. The membrane cover is mounted on a partially fixed frame, allowing in-plane deformation, and dominantly displaced in the OOP direction by a hydraulic car jack located approximately at the centre of the membrane. The membrane is loaded six times, at first incrementally to a height of 23 mm, 48 mm, and 72 mm, then released to 63 mm and 51 mm, and a final loading to a height of 86 mm, which formed the six OOP deformation tests denoted as Step 1 to 6, respectively, for this study. Optical fibres were installed and aligned with the vertical lines of 1000 mm at 800 mm and 1000 mm horizontally on the HDPE specimen for benchmarking the strains; refer to In the OOP deformation test, a 1200 mm squared HDPE with an internal 1000 mm squared grid marked on the surface is used, as shown in Figure 4. The membrane cover is mounted on a partially fixed frame, allowing in-plane deformation, and dominantly displaced in the OOP direction by a hydraulic car jack located approximately at the centre of the membrane. The membrane is loaded six times, at first incrementally to a height of 23 mm, 48 mm, and 72 mm, then released to 63 mm and 51 mm, and a final loading to a height of 86 mm, which formed the six OOP deformation tests denoted as Step 1 to 6, respectively, for this study. Optical fibres were installed and aligned with the vertical lines of 1000 mm at 800 mm and 1000 mm horizontally on the HDPE specimen for benchmarking the strains; refer to  In the in-plane deformation test, a 1200 mm by 500 mm rectangular HDPE with an internal 1000 mm by 300 mm grid is used, as shown in Figure 5. The membrane cover is fixed on one end while the other end is translated in the in-plane direction towards the fixed edge by 112 mm, 161 mm, and 267 mm. These imposed displacements caused the membrane to form a wrinkle/fold denoted as Step 1 to 3, respectively. The maximum fold heights are 190 mm, 223 mm, and 280 mm for each subsequent In the in-plane deformation test, a 1200 mm by 500 mm rectangular HDPE with an internal 1000 mm by 300 mm grid is used, as shown in Figure 5. The membrane cover is fixed on one end while the other end is translated in the in-plane direction towards the fixed edge by 112 mm, 161 mm, Remote Sens. 2020, 12, 2738 7 of 21 and 267 mm. These imposed displacements caused the membrane to form a wrinkle/fold denoted as Step 1 to 3, respectively. The maximum fold heights are 190 mm, 223 mm, and 280 mm for each subsequent step test. Optical fibres were installed and aligned with the horizontal line of 1000 mm at 180 mm and 320 mm (80 mm and 220 mm horizontal grid lines) vertically on the HDPE specimen for benchmarking; refer to Figure 5.  The preprocessing of experimental data is performed using MATLAB. The photogrammetry DEM is transformed into a readable format, and then the fiducial grid is extracted from each test. A contrast threshold technique is used to identify and extract the grid, and a stationary reference outside of the membrane is used to align the DEM to a global coordinate system. The coordinates and elevation of the grid lines from its corresponding DEM and for each loading case are then extracted in-aid of MATLAB's inbuilt computer vision tool detectSURFFeatures.
The cubic smoothing spline is then employed to smooth the model and filter noise due to artefacts, such as light reflection and dirt on the membrane, as a function of weighting parameter, , from 10 to 10 . Afterwards, the preprocessed displacement fields are used as the displacement grid to load the FE membrane model. In this work, the DEM OOP deformation field is fully applied as the grid load and the in-plane deformation fields are only applied on the top and bottom boundaries of the grid for the OOP deformation test and applied on the left boundary of the grid for the in-plane deformation test while the remaining are excluded and set to freely deform in the inplane directions and to ensure convergence of the solution.
ANSYS 19.2 is used as the FE computational analysis tool to simulate the deformation of the HDPE membrane using the processed DEM displacement field as displacement load conditions. A shell membrane element was used to model the membrane cover of mesh discretisation of 5 mm. Given the quasi-static nature of the membrane, ANSYS Static Structural analysis is considered. ANSYS Parametric Design Language scripts were developed to assign and displace each node on the grid. For nonlinear FE analysis, Large Deflection is set ON and in the experimental investigations, ANSYS Solver Controls Weak Springs is set to Program Controlled to restrain rigid body motion. The simulated smoothed and raw models' strain fields are compared to the strain reading from optical fibre for each deformation test and step.

Statistical Process Approach: Gaussian Process
Significant smoothing of the data or model inevitably results in undesirable amplitude reduction. Photogrammetry data and images are subjected to noise and processed data after smoothing interferes with the precise measurement of the elevation/height causing erroneous strain prediction, especially in areas with high displacement gradients. Nevertheless, unavoidable uncertainties exist throughout the stages from acquisition of data to output of results, therefore, the deterministic prediction of strain values is difficult, and hence, this motivates a probabilistic approach to tackle these inevitable variations. Given the practical application in hand, the only information obtained are the raw data and smoothed models with different weighting parameters constructed from the photogrammetry DEM, and the true strain measurements remain unknown. A statistical The preprocessing of experimental data is performed using MATLAB. The photogrammetry DEM is transformed into a readable format, and then the fiducial grid is extracted from each test. A contrast threshold technique is used to identify and extract the grid, and a stationary reference outside of the membrane is used to align the DEM to a global coordinate system. The coordinates and elevation of the grid lines from its corresponding DEM and for each loading case are then extracted in-aid of MATLAB's inbuilt computer vision tool detectSURFFeatures.
The cubic smoothing spline is then employed to smooth the model and filter noise due to artefacts, such as light reflection and dirt on the membrane, as a function of weighting parameter, p, from 10 −1 to 10 −6 . Afterwards, the preprocessed displacement fields are used as the displacement grid to load the FE membrane model. In this work, the DEM OOP deformation field is fully applied as the grid load and the in-plane deformation fields are only applied on the top and bottom boundaries of the grid for the OOP deformation test and applied on the left boundary of the grid for the in-plane deformation test while the remaining are excluded and set to freely deform in the in-plane directions and to ensure convergence of the solution.
ANSYS 19.2 is used as the FE computational analysis tool to simulate the deformation of the HDPE membrane using the processed DEM displacement field as displacement load conditions. A shell membrane element was used to model the membrane cover of mesh discretisation of 5 mm. Given the quasi-static nature of the membrane, ANSYS Static Structural analysis is considered. ANSYS Parametric Design Language scripts were developed to assign and displace each node on the grid. For nonlinear FE analysis, Large Deflection is set ON and in the experimental investigations, ANSYS Solver Controls Weak Springs is set to Program Controlled to restrain rigid body motion. The simulated smoothed and raw models' strain fields are compared to the strain reading from optical fibre for each deformation test and step.

Statistical Process Approach: Gaussian Process
Significant smoothing of the data or model inevitably results in undesirable amplitude reduction. Photogrammetry data and images are subjected to noise and processed data after smoothing interferes with the precise measurement of the elevation/height causing erroneous strain prediction, especially in areas with high displacement gradients. Nevertheless, unavoidable uncertainties exist throughout the stages from acquisition of data to output of results, therefore, the deterministic prediction of strain values is difficult, and hence, this motivates a probabilistic approach to tackle these inevitable variations. Given the practical application in hand, the only information obtained are the raw data and smoothed models with different weighting parameters constructed from the photogrammetry DEM, and the true strain measurements remain unknown. A statistical process approach using sample data of the smoothed and raw strain field is employed to provide a probabilistic prediction in evaluating the likelihood of strain values of the structure.

Gaussian Process Regression
In this study, GP regression is demonstrated as a probabilistic prediction approach for strain determination. GP, which assumes that the joint probability distribution of model outputs is Gaussian, is a Bayesian approach used in applications for model approximation, multivariate regression, forecasting, and experiment design and as an algorithm for regression and classification for machine learning [53][54][55]. The fundamentals of GP regression are briefly introduced, and a more comprehensive report can be found in previous literature [53][54][55][56].
GP is a stochastic process in which any finite subset through its domain follows a multivariate normal distribution [56]. Consider the following nonlinear regression model with homoscedastic noise (constant variance) of observed data D = (x 1 , y 1 ), . . . , (x n , y n ) , where ε i ∼ N 0, σ 2 ε represents the measurement error, which is independent and identically distributed normal random noises with zero mean and variance σ 2 ε , and f (x) is an unknown function. The GP regression model can then be denoted as by the GP method, f (x) is related as a random function and assumed to have a GP prior with a mean, µ(x), and covariance k(x, , where θ denotes the set of hyperparameters. Hence, the joint distribution of the y 1 , . . . , y n is multivariate normal. where µ has entries µ i = µ(x i ) and Ψ is a n by n matrix whose (i, j) th element is given by where δ ij is Kronecker delta. For a new input x * and y * , the corresponding response value is calculated as where Therefore, the lower and upper limits are, respectively, Remote Sens. 2020, 12, 2738

of 21
A simple assumption is that the correlation between two points decays with the distance between the points according to an exponential function. One prevalent choice of a kernel fulfilling this assumption is the squared-exponential kernel function.
where |x − x * | is the Euclidean distance between two input points and θ denotes the collection of hyperparameters associated with the kernel function, in this caseσ and l.σ 2 is the signal deviation and l is the characteristic length scale, which defines how far apart the input values of x can be for the response values to become uncorrelated. The hyperparameters are usually not assumed to be known but are trained by maximising the marginal log-likelihood from the dataset [54,55],

Heteroscedastic Gaussian Process Regression
It is highly desirable to fit a regression model with non-constant variance than the standard GP, where a constant variance is assumed, which is often unrealistic in many applications. There have been developments in novel GP approaches to regression with input-dependent noise rates and demonstrated higher accurate modelling of real-world datasets compared to other regression methods [57][58][59][60][61]. In this study, variational heteroscedastic GP (VHGP) regression is considered and a brief explanation of the theory will be shown in this section. More detail can be found in the work by Gredilla and Titsias [58].
Similar to the Equation (2), instead of the GP prior on f (x) and the error term, we have where r(x) is an unknown function and the function is defined as r(x) = e g(x) to ensure positivity and without losing generality, and Once the kernel functions are defined, heteroscedastic GP is fully specified and depends only on its hyperparameters (θ, θ g and µ 0 ). Unfortunately, the exact inference in the heteroscedastic GPs is no longer analytically tractable and, thus, a marginalised variational approximation bound [58] for the likelihood function is given by where R is a diagonal matrix with elements R ii = e µ i −Σ ii /2 , and KL is Kullback-Leibler divergence. Furthermore, the stationary equations must be satisfied at any local or global maximum and the two equations are obtained as follows, for some positive semidefinite diagonal matrix Λ. µ and Σ −1 are expressed as a function of Λ. So, the marginalised variational bound can be a function of Λ,

F(µ(Λ), Σ(Λ)) = F(Λ)
which needs to be maximised with respect to the n variational parameters in Σ. By implementing the marginal log likelihood for model selection, we can maximise F with respect to the model hyperparameters for a new test point. This study employs the VHGP regression method for data modelling to deal with the different sources of uncertainty that originated from photogrammetry, preprocessing, and postprocessing analyses. The objective of using an inference statistical approach is to provide descriptive insight and prediction to the strain field. In this work, automatic relevance determination squared-exponential kernel function is defined as and employed for f (x) and g(x). All strain results for all p-values were used as the training data set and the predicted VHGP models are compared to the experimental optical fibre strain data. Figure 6 shows OOP deformation test Step 6 strain fields and their relative residual displacement compared to the raw data (no smoothing). The preprocessing smoothing technique on the DEM significantly improved and denoised its strain profile and uncovered its relevant features which improve the accuracy in identifying the location of the maximum strain, thus providing qualitative information on the evolution of the strain field. However, it is shown that overly smoothed profiles result in flattening of peaks, loss of localised and fine-detailed strain and displacement information.

Out-of-Plane Deformation Test
Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 22 Figure 6 shows OOP deformation test Step 6 strain fields and their relative residual displacement compared to the raw data (no smoothing). The preprocessing smoothing technique on the DEM significantly improved and denoised its strain profile and uncovered its relevant features which improve the accuracy in identifying the location of the maximum strain, thus providing qualitative information on the evolution of the strain field. However, it is shown that overly smoothed profiles result in flattening of peaks, loss of localised and fine-detailed strain and displacement information.  There are significant improvements in the strain profile when is smaller than 10 in comparison with the measured strain and the location of predicted maximum strain is within 8.8 mm of the location of the measured maximum strain value. The strain profile gradually improves as decreases, however, the maximum strain value decreases for all tests. An average maximum strain reduction of 51.3%, 57.0%, 63.0%, and 79.2% was observed for values 10 to 10 , respectively, relative to the measured strain.  There are significant improvements in the strain profile when p is smaller than 10 −3 in comparison with the measured strain and the location of predicted maximum strain is within 8.8 mm of the location of the measured maximum strain value. The strain profile gradually improves as p decreases, however, the maximum strain value decreases for all tests. An average maximum strain reduction of 51.3%, 57.0%, 63.0%, and 79.2% was observed for p values 10 −3 to 10 −6 , respectively, relative to the measured strain.

Out-of-Plane Deformation Test
smoothing spline parameters and its residual out-of-plane displacement relative to raw data (no filter). Figures 7-12 show the strain distributions along the vertical line at 800 mm horizontally for the tests with and without the filtering process and those recorded by the optical fibre (measured strain). There are significant improvements in the strain profile when is smaller than 10 in comparison with the measured strain and the location of predicted maximum strain is within 8.8 mm of the location of the measured maximum strain value. The strain profile gradually improves as decreases, however, the maximum strain value decreases for all tests. An average maximum strain reduction of 51.3%, 57.0%, 63.0%, and 79.2% was observed for values 10 to 10 , respectively, relative to the measured strain.          In view of a practical setting, the initial stress condition of installed membrane cover is unknown and therefore it is also of particular interest to compare the change in strain, ∆ , from different test loading. Figure 13 shows the ∆ of Step 6 relative to the preceding tests. The ∆ profiles are similar compared to the measured strain recorded when the model is smoother. However, a reduction of ∆ value, especially in the vicinity of the maximum/minimum ∆ value, is also observed. In Figure 13, the average maximum ∆ reductions for values 10 to 10 are 29.2%, 40.23%, 49.2%, and 74.9%, respectively, relative to the experimental maximum ∆ along the vertical line at 800mm. In Figure 14, the average maximum ∆ reductions for values 10 to 10 are 59.5%, 66.1%, 63.5%, and 61.8%, respectively, relative to the experimental maximum ∆ along the vertical line at 1000 mm. The findings indicate that the change in strain profiles and the location of the maximum change in strain can also be accurately predicted by smoothing the model. Furthermore, the reduction of maximum ∆ value is also observed.  In view of a practical setting, the initial stress condition of installed membrane cover is unknown and therefore it is also of particular interest to compare the change in strain, ∆ , from different test loading. Figure 13 shows the ∆ of Step 6 relative to the preceding tests. The ∆ profiles are similar compared to the measured strain recorded when the model is smoother. However, a reduction of ∆ value, especially in the vicinity of the maximum/minimum ∆ value, is also observed. In Figure 13, the average maximum ∆ reductions for values 10 to 10 are 29.2%, 40.23%, 49.2%, and 74.9%, respectively, relative to the experimental maximum ∆ along the vertical line at 800mm. In Figure 14, the average maximum ∆ reductions for values 10 to 10 are 59.5%, 66.1%, 63.5%, and 61.8%, respectively, relative to the experimental maximum ∆ along the vertical line at 1000 mm. The findings indicate that the change in strain profiles and the location of the maximum change in strain can also be accurately predicted by smoothing the model. Furthermore, the reduction of maximum ∆ value is also observed. In view of a practical setting, the initial stress condition of installed membrane cover is unknown and therefore it is also of particular interest to compare the change in strain, ∆ε, from different test loading. Figure 13 shows the ∆ε of Step 6 relative to the preceding tests. The ∆ε profiles are similar compared to the measured strain recorded when the model is smoother. However, a reduction of ∆ε value, especially in the vicinity of the maximum/minimum ∆ε value, is also observed. In Figure 13, the average maximum ∆ε reductions for p values 10 −3 to 10 −6 are 29.2%, 40.23%, 49.2%, and 74.9%, respectively, relative to the experimental maximum ∆ε along the vertical line at 800 mm. In Figure 14, the average maximum ∆ε reductions for p values 10 −3 to 10 −6 are 59.5%, 66.1%, 63.5%, and 61.8%, Remote Sens. 2020, 12, 2738 13 of 21 respectively, relative to the experimental maximum ∆ε along the vertical line at 1000 mm. The findings indicate that the change in strain profiles and the location of the maximum change in strain can also be accurately predicted by smoothing the model. Furthermore, the reduction of maximum ∆ε value is also observed.   Figure 15 shows the strain in the horizontal direction of the first in-plane displacement test (Step 1) for different , and similarly to the previous findings, the quality of the strain profile significantly improves as the model is smoothed with the loss of detail strain information. Figures 16-18 show the strain distributions along the horizontal line at 180 mm for the three tests with and without the filtering process and recorded by the optical fibre. There are significant improvements in the strain profile when is smaller than 10 in comparison with the measured strain, and the location of predicted maximum strain is within 29.9 mm of the location of the measured maximum strain value. The strain profile gradually improves as decreases and the maximum strain value gradually decreases for all tests, similar to the previous investigation. There is an average maximum strain reduction of 46.9%, 55.1%, 60.7%, and 65.6% for values 10 to 10 , respectively, relative to the optical fibre-measured experimental strain. It should be noted that the FE analysis of Step 3 raw strain distribution did not converge due to noise of the actual displacement measurement.   Figure 15 shows the strain in the horizontal direction of the first in-plane displacement test (Step 1) for different , and similarly to the previous findings, the quality of the strain profile significantly improves as the model is smoothed with the loss of detail strain information. Figures 16-18 show the strain distributions along the horizontal line at 180 mm for the three tests with and without the filtering process and recorded by the optical fibre. There are significant improvements in the strain profile when is smaller than 10 in comparison with the measured strain, and the location of predicted maximum strain is within 29.9 mm of the location of the measured maximum strain value. The strain profile gradually improves as decreases and the maximum strain value gradually decreases for all tests, similar to the previous investigation. There is an average maximum strain reduction of 46.9%, 55.1%, 60.7%, and 65.6% for values 10 to 10 , respectively, relative to the optical fibre-measured experimental strain. It should be noted that the FE analysis of Step 3 raw strain distribution did not converge due to noise of the actual displacement measurement.  Figure 15 shows the strain in the horizontal direction of the first in-plane displacement test (Step 1) for different p, and similarly to the previous findings, the quality of the strain profile significantly improves as the model is smoothed with the loss of detail strain information. Figures 16-18 show the strain distributions along the horizontal line at 180 mm for the three tests with and without the filtering process and recorded by the optical fibre. There are significant improvements in the strain profile when p is smaller than 10 −3 in comparison with the measured strain, and the location of predicted maximum strain is within 29.9 mm of the location of the measured maximum strain value. The strain profile gradually improves as p decreases and the maximum strain value gradually decreases for all tests, similar to the previous investigation. There is an average maximum strain reduction of 46.9%, 55.1%, 60.7%, and 65.6% for p values 10 −3 to 10 −6 , respectively, relative to the optical fibre-measured experimental strain. It should be noted that the FE analysis of Step 3 raw strain distribution did not converge due to noise of the actual displacement measurement.

In-Plane Displacement Test (Wrinkle Formation)
Remote Sens. 2020, 12, x FOR PEER REVIEW 14 of 22 Figure 15. In-plane deformation test Step 1: HDPE membrane strain in the horizontal direction of (a) raw data and different cubic smoothing spline weighting parameter; (b-g) 10 , 10 , 10 , 10 , 10 , and 10 , respectively. Figures 16-18 show improvement of strain in the vicinity of 400 mm and 800 mm as the model smoothed; however, significant changes to strain as decreases suggest a high level of noise in those particular locations. It is observed that the strain slope is different compared to the optical fibre measurement at 400 mm and 800 mm. Furthermore, the development of these abrupt high compressive strains is earlier seen in Figure 16 in the same vicinities. The reason that it is causing such inaccuracy is due to the reconstructed DEM model at the steep regions (high slope gradient) along the membrane.  however, significant changes to strain as decreases suggest a high level of noise in those particular locations. It is observed that the strain slope is different compared to the optical fibre measurement at 400 mm and 800 mm. Furthermore, the development of these abrupt high compressive strains is earlier seen in Figure 16 in the same vicinities. The reason that it is causing such inaccuracy is due to the reconstructed DEM model at the steep regions (high slope gradient) along the membrane.     Step 3 test relative to the preceding tests. In Figure 19, the average maximum ∆ reduction for values 10 to 10 are 25.7%, 32.7%, 46.0%, and 58.5%, respectively, relative to the experimental maximum ∆ along the horizontal line at 180 mm. In Figure  20, the average maximum ∆ reduction for values 10 to 10 are 38.4%, 60.8%, 64.6%, and 68.6%, respectively, relative to the experimental maximum ∆ along the horizontal line at 320 mm. As expected, the changes in strain profiles are different, especially near the 400 mm and 800 mm vicinity compared to the experimental strain recorded by the optical fibre when the model is smoother. Nevertheless, the development of the strain, as well as the change in strain profiles as a function of , provided qualitative guidance of the strain information, and to further progress in strain evaluation, a statistical approach is considered in the next section.  Figures 16-18 show improvement of strain in the vicinity of 400 mm and 800 mm as the model smoothed; however, significant changes to strain as p decreases suggest a high level of noise in those particular locations. It is observed that the strain slope is different compared to the optical fibre measurement at 400 mm and 800 mm. Furthermore, the development of these abrupt high compressive strains is earlier seen in Figure 16 in the same vicinities. The reason that it is causing such inaccuracy is due to the reconstructed DEM model at the steep regions (high slope gradient) along the membrane. Figures 19 and 20 show the ∆ε of Step 3 test relative to the preceding tests. In Figure 19, the average maximum ∆ε reduction for p values 10 −3 to 10 −6 are 25.7%, 32.7%, 46.0%, and 58.5%, respectively, relative to the experimental maximum ∆ε along the horizontal line at 180 mm. In Figure 20, the average maximum ∆ε reduction for p values 10 −3 to 10 −6 are 38.4%, 60.8%, 64.6%, and 68.6%, respectively, relative to the experimental maximum ∆ε along the horizontal line at 320 mm. As expected, the changes in strain profiles are different, especially near the 400 mm and 800 mm vicinity compared to the experimental strain recorded by the optical fibre when the model is smoother. Nevertheless, the development of the strain, as well as the change in strain profiles as a function of p, provided qualitative guidance of the strain information, and to further progress in strain evaluation, a statistical approach is considered in the next section.  Step 2, at 180 mm of optical fibre and different cubic smoothing spline weighting parameter from 10 −3 to 10 −6 . Figure 19. In-plane deformation test: Change in strain of Step 3 relative to the preceding strain; (a) Step 1 and (b) Step 2, at 180 mm of optical fibre and different cubic smoothing spline weighting parameter from 10 to 10 .

Figure 20.
In-plane deformation test: Change in strain of Step 3 relative to the preceding strain; (a) Step 1 and (b) Step 2, at 320 mm of optical fibre and different cubic smoothing spline weighting parameter from 10 to 10 .

Strain Determination in Using VHGP
The strain profile results from OOP deformation and in-plane deformation tests are used as sample data for VHGP to predict the strain field. Figures 21 and 22 show the measured strain, VHGP maximum predicted mean strain with 95% of certainty (2 * confidence intervals) of the OOP and inplane deformation tests, respectively. The measured maximum strains are larger than the predicted mean strain and are within the 2 * confidence interval, specifically in the upper limit half, of the VHGP model for almost all tests, except for OOP test Step 6 which is approximately +2.15 * from VHGP maximum strain prediction.
There are significantly larger variances in the VHGP models of the in-plane deformation tests which are due to the significant noise (in the vicinity of 400 mm and 800 mm) compared to the OOP tests. It is shown that the VHGP model provided acceptable estimations and uncertainty measurements of the maximum strain of OOP and in-plane deformation tests by incorporating the data of raw and smoothed strain fields.  1 and (b) Step 2, at 320 mm of optical fibre and different cubic smoothing spline weighting parameter from 10 −3 to 10 −6 .

Strain Determination in Using VHGP
The strain profile results from OOP deformation and in-plane deformation tests are used as sample data for VHGP to predict the strain field. Figures 21 and 22 show the measured strain, VHGP maximum predicted mean strain with 95% of certainty (2σ * confidence intervals) of the OOP and in-plane deformation tests, respectively. The measured maximum strains are larger than the predicted mean strain and are within the 2σ * confidence interval, specifically in the upper limit half, of the VHGP model for almost all tests, except for OOP test Step 6 which is approximately +2.15σ * from VHGP maximum strain prediction.
There are significantly larger variances in the VHGP models of the in-plane deformation tests which are due to the significant noise (in the vicinity of 400 mm and 800 mm) compared to the OOP tests. It is shown that the VHGP model provided acceptable estimations and uncertainty measurements of the maximum strain of OOP and in-plane deformation tests by incorporating the data of raw and smoothed strain fields.

Discussion
The UAV photogrammetry inspection method is shown to be effective in monitoring the state of deformation of the floating covers at Melbourne Water WTP [37]. In comparison to non-contact methods on measuring strain fields, although DIC methods may yield higher accuracy, DIC methods require the membrane to be coated with random patterns and multiple fixed cameras continuously

Discussion
The UAV photogrammetry inspection method is shown to be effective in monitoring the state of deformation of the floating covers at Melbourne Water WTP [37]. In comparison to non-contact methods on measuring strain fields, although DIC methods may yield higher accuracy, DIC methods require the membrane to be coated with random patterns and multiple fixed cameras continuously

Discussion
The UAV photogrammetry inspection method is shown to be effective in monitoring the state of deformation of the floating covers at Melbourne Water WTP [37]. In comparison to non-contact methods on measuring strain fields, although DIC methods may yield higher accuracy, DIC methods require the membrane to be coated with random patterns and multiple fixed cameras continuously monitoring to capture the whole asset-these requirements are not practicable for WTP. Furthermore, DIC methods are not effective for monitoring wrinkle and fold formations, which makes photogrammetry methods more suitable for monitoring membrane cover [27]. To further facilitate the development of smart SHM assessment on this highly valuable asset, the digital displacement information of the covers can be further processed into strain measurement to allow a more elaborative quantification of its integrity. Overall, this paper introduces a strain determination method proposed for WTP floating cover and the UAV photogrammetry deployment and its general application for large structures are briefly discussed.
The proposed strain determination method is suitable for global inspection in the two-stage monitoring strategy, as it demonstrates the capability to produce a higher quality strain distribution profile of the deformed membrane for detecting critical areas. Furthermore, the VHGP offers a probabilistic prediction approach that supplements the decision support to justify further localised inspection in concerning areas. The results have shown that the use of cubic smoothing interpolation technique can reliably denoise the DEM, improving the quality of the strain field of deformed membranes to locate high strain areas in both OOP and in-plane deformation tests. However, the smoothing procedure adversely degrades fine detail information and flattens peaks, which is seen for increasing deformation test steps. As a result, there is a difficulty in capturing the numeric strain values at high strain concentrated areas. It should be worth noting that the optimisation of the grid resolution should be carefully considered in both the capability of the photogrammetry, camera, and image resolution and in accurately capturing critical deformations or details of the specimen or structure when implementing this approach for strain determination via UAV photogrammetry [21]. In our application, WTP regulations restrict additional markings of the floating cover and only existing features can be used as targeted fiducial markers, which limits the ability to enhance the resolution for an accurate reading. Furthermore, due to the disturbance of wildlife, fire hazard, and WTP regulations, UAV photogrammetry cannot be performed too close to the floating cover for more detailed measurements. To facilitate the strain evaluation, a probabilistic approach is implemented to predict the strain values with statistical measurement by using the smoothed and raw model strain fields as data samples. A VHGP method is demonstrated to provide uncertainty measurements to assist in evaluating the strain field within the 95% confidence level by using only the available smoothed and raw models. As VHGP searches a distribution over the possible functions that are consistent with the data, large variances occur in the vicinity of strain concentrated areas and noise due to its relatively high changes in strain value as smoothing weighting parameter changes, thereby corresponding to higher uncertainty variance in those areas.
A factor contributing to inaccurate strain, as discussed earlier, is the reconstruction of DEM using UAV photogrammetry in the region with a steep slope. As reported in Wong et al. [36], UAV photogrammetry has difficulties in reconstructing the DEM model with sharp corners and high slope gradients. The photogrammetry process may consider the steep region as a discontinuity in the reconstructed model and automatically apply a smoothing function to construct a smooth DEM. For this laboratory experiment, the aerial images were only taken directly above the region of interest. More viewing angles can improve the accuracy of close-range photogrammetry analysis.
The lack of in-plane direction information on the model may have contributed to the inaccuracy of the strain field. However, it was found that by including the experimental full-field displacement as the applied load led to either, mostly, un-converged solution or highly inaccurate strain results with no resemblance to the measured strain profile and values. It was revealed that there were inconsistent feature detections and extractions of points along the grid line. Some data points were extracted within the 1.5 mm width of the line rather than consistently from the centre. These are artificially induced erroneous in-plane deformations, resulting in extremely noisy and unreliable data in the in-plane direction. Therefore, only the in-plane displacements on the boundaries of the grid were manually identified, selected, and loaded to provide sufficient constraint for FE analysis. Further improvement in the object-identification algorithm can accurately identify and extract the markers. Nevertheless, ongoing future work will investigate improving the quality of DEM, to explore feature identifying and optimize selection on the fiducial markers on the WTP floating covers, as well as other strain prediction algorithms.

Conclusions
In the development of an effective two-stage monitoring strategy for Melbourne Water's WTP floating cover, this paper presented a strain determination technique which incorporates UAV photogrammetry for rapid monitoring in locating critical areas, whilst providing statistical insight for further localised evaluation. The effectiveness of UAV photogrammetry is known to capture larger structures and membrane deformation, such as folds and wrinkles, making it suitable for the large floating covers in WTP. This study has shown a method which utilises the raw DEM of deformed membrane covers from UAV photogrammetry to predict the strain field using FE analysis and a statistical approach. Cubic smoothing interpolation was performed to denoise the DEM, which improves the DEM quality for FE analysis to produce a higher quality strain field that has similar results to the associated experimental strain profile recorded by the optical fibres. However, as the smoothing parameter changes, it is shown that smoothing results in loss of amplitude and detailed information which is most apparent in the vicinity of the maximum strain locations. A VHGP regression method was applied, by using the smoothed model strain fields as training data, in order to further assist in the evaluation of those strain values. It is shown that the predicted strain values are within 95% certainty and critical areas have high variance due to significant noise and high changes of strain values as the smoothing parameter changes. The statistical approach provides the likelihood of strain value estimations and, thereby, vital information for decision-making processes at the operational and management levels. Our future work will continue investigating the strain determination of large membrane covers, which includes improving the accuracy and automation to be applied on WTP floating covers, towards smart structure development. More ongoing on-the-field work on the UAV photogrammetry method is currently underway.