1 Introduction

Maintenance management of infrastructure networks has been a challenging task for civil asset owners. While the number of deteriorating structures in a network increases exponentially every year, the structural assessment of a single asset, for instance a railway or a motorway bridge, still requires significant effort, time and cost. The structural assessment in everyday engineering practise is a complex process that may include frequent site inspections, sensing and surveying monitoring, non-destructive testing, material characterisation testing and structural modelling.

Optimising this process to evaluate local material deterioration and global structural performance appears particularly challenging for ageing infrastructure. Masonry railway bridges for instance, which comprise around 50–60% of the UK and European rail stock [1, 2], provide an example of civil assets where structural assessment based on common practices is characterised by a high level of uncertainty. One of the main reasons is the unilateral, nonlinear behaviour of masonry due to its inherently discontinuous nature. The difficulty of modelling masonry has led to the development of a wide range of simplified and advanced modelling approaches. These approaches may treat masonry as continuous or discontinuous medium, capturing each time different aspects of its response in the micro or macro scale [3, 4]. Another difficulty is the assessment of the influence of the backfill material, which significantly influences the stability of the arch barrel. As before, a wide range of simplified and advanced modelling approaches exist to account for the backfill effect [5,6,7,8,9]. Moreover, most of these structures have been built more than a century ago, before the enforcements of building codes and design details are unknown. Destructive testing is often avoided in case of structures with cultural value, which limits further engineers to choose reliably material parameters for modelling. Furthermore, the long deformation and loading history of the structure over the decades is unknown due to the absence of monitoring data.

Due to these modelling challenges, measurement-based assessment and maintenance procedures appear attractive. For example, the use of accelerometers for vibration-based damage detection has been at the forefront of research for the SHM of bridges [10,11,12] over the last decades. In the case of masonry bridges though, identifying damage by tracking changes in the modal or frequency domain is more challenging. Masonry arch bridges are relatively stiff, are highly indeterminate, and often contain a significant number of cracks, which poses challenges for modal identification-based damage methods.

Recently, different techniques for distant monitoring have been proposed for these structures: long-term deformation and distortions in geometry can be identified through LiDAR data [1, 13,14,15]; image-based crack detection algorithms can facilitate inspection of heritage masonry structures [16]; field digital image correlation (DIC) can be used to monitor dynamic displacements [17]. These surveying technologies can offer rapid structural assessment solutions and a good first understanding of the overall response of the structure. For continuous long-term deterioration monitoring though, and considering the difficulties to evaluate the structural performance discussed before, higher accuracy is required.

Following this reasoning, a fibre optic monitoring system was trialled to monitor a masonry railway bridge in Leeds, UK [18]. Subsequently, this trial was improved to an autonomous, remote, multi-sensing system to monitor long-term degradation [19,20,21]. The main component of the system is a network of fibre Bragg grating (FBG) sensors that monitors distributed dynamic strain along three arch spans of the bridge. The system is complemented by acoustic emission (AE) sensors and high-sensitivity accelerometers. A preliminary statistical analysis of FBG strain and temperature data over a period of two years indicated changes of dynamic strain due to cyclic temperature variations and material degradation [20, 21]. A more recent FBG deployment on a skewed arch railway bridge in UK showed that with FBG networks, monitoring of multiple aspects of the structural behaviour of masonry railway bridges can be achieved, including visualisation of the flow of forces under train loading [22].

Following these deployments and recent research efforts on self-sensing composite bridges with embedded FBGs combined with statistical modelling [23], this paper examines the use of statistical shape analysis (SSA) [24] to enhance damage detection and structural assessment of ageing railway infrastructure assets monitored with fibre optics. The method is appealing for railway bridge monitoring, where repeated loading from specific train type corresponds to a signature response of the bridge, which results in FBG signals of similar ‘shape’, as will be discussed in the following sections.

In this paper, SSA is used to track changes in the shape of FBG signals at multiple locations in a masonry railway bridge in Leeds, UK. Section 2 discusses the bridge pathology, the sensing system and FBG data, together with previous findings from a train classification algorithm combined with statistical analysis. Section 3 introduces SSA, and more specifically, the ordinary procrustes analysis (OPA), and shows how this method can be adapted for railway asset monitoring, to account separately for changes in the shape of signals due to variations in train speed and dynamic strain. Section 4 presents the results of SSA, complemented by a trained support vector machine (SVM) model. Furthermore, the section investigates how identification of strain changes related to mechanical damage can be distinguished from environmental effects. Based on these findings, the paper discusses that distributed sensing approaches, such as high-sensitivity fibre optic networks, combined with data analytics, may enable a much finer level of deterioration monitoring. Finally, Sect. 5 summarizes the conclusions of the study.

2 Case study: the Marsh Lane Bridge

2.1 Bridge pathology and the fibre optic sensing system

The Marsh Lane Bridge is a typical Victorian viaduct in the UK, located next to the Eastern entrance of Leeds Railway Station, near the city centre. The bridge was constructed between 1865 and 1869 [25]. Today it carries two electrified tracks with a traffic load that exceeds 200 trains per day, ranging from typical passenger trains to multi-wagon freight trains. Figure 1-left shows the Northern side of the investigated section of the bridge, which comprises Arches 37, 38 and 39.

Fig. 1
figure 1

Left—Northern view of the Marsh Lane Viaduct in Leeds, UK. Right—Closer view of the pier between Arches 38 and 39, showing the 2015 retrofit

The bridge appears significantly damaged due to the increase of the operational train loads over the last decades and due to environmental effects. For this reason, it was repaired in 2015 by filling the relieving arches of the piers with concrete and installing steel ties to arrest the transverse movements of the piers and spandrel walls, as shown in Fig. 1-right. Figure 2 is the plan view of the bridge showing with red lines the main cracks observed under the arches. The most severe damages are concentrated over the relieving arches at the centre of the piers due to a spreading mechanism that forces the relieving arch keystone to descend and the pier walls to bow outwards [18]. The pier between Arches 37 and 38 is surrounded by transverse cracks at the level where the backing meets the backfill material (see Fig. 1-right). These cracks are due to the intense out-of-plane rocking response of the pier under train loading before the 2015 retrofit [18]. In the longitudinal direction, separation cracks between the spandrel walls and the arch barrels are observed in Arches 37 and 38 as well as a longitudinal crack in Arch 37 develops below the North railway track.

Fig. 2
figure 2

Plan view of the investigated section of the bridge, showing the main damages and location of the fibre Bragg grating network

In July 2016, a few months after the bridge repair, a network of FBG strain sensors was installed underneath Arches 37 and 38, allowing for the detailed study of their dynamic deformation [18]. In November, 2017, this proof-of-concept, temporary monitoring installation was replaced by a permanent FBG monitoring installation [19,20,21] underneath Arches 37, 38 and 39, shown in Fig. 2, to assess the long-term effectiveness of the 2015 retrofitting intervention and to monitor bridge deterioration. In Fig. 2, the blue lines along the longitudinal and transverse direction of the bridge correspond to the FBG network of strain sensors and the black triangles to FBG temperature sensors.

A close view of the FBG strain and temperature sensors is offered in Fig. 3. In the case of strain sensors, the fibre optic cable is pre-tensioned and the Bragg gratings are located between aluminium clamps. In that way, every FBG sensor measures strain variation between two fixed points, which are about 1 m apart.

Fig. 3
figure 3

Fibre Bragg grating strain and temperature sensor

There are various ways of measuring temperature with FBGs. A common method is to have a loose bare fibre in an enclosure (tube), which can expand thermally without any influence of external constraints. This is appropriate when sensors are embedded into concrete structures (see for instance [26] and references reported therein). Another method is to have the FBG bonded with a uniform solid, usually metal, which expands thermally and transfers the resulting strain to the fibre (such as the commercially available Temperature Compensation FBG sensors that are adjusted to steel ‘carriers’). These are appropriate for surface mounting applications.

The FBG temperature sensor developed for this study (see Fig. 3-bottom left) consists of two machined rectangle aluminium plates with a tight groove fitted around the cable, screwed together. The clamps are long enough to transfer strain due to thermal expansion of the aluminium, while the cable that exits the clamp is loose to omit any influence of structural deformation.

The temperature clamps were calibrated in a TX 150 Grant water bath at the laboratory of the Centre for Smart Infrastructure and Construction (CSIC), University of Cambridge, UK. The relationship between the normalised wavelength change and temperature was linear and the slope obtained was 24.25 10–6/oC, which is a typical value of the thermal expansion coefficient of aluminium, as expected.

The FBG network shown in Fig. 2 comprises: 47 strain sensors in the longitudinal direction; 17 strain sensors in the transverse direction; 4 sensors attached to steel wires that connect the springings of Arch 37 and Arch 39, which are located underneath the longitudinal strain sensors arrays, to measure variations in arch span (see also Fig. 4-top); and lastly, 5 temperature sensors.

Fig. 4
figure 4

FBG strain signals at four bridge location for the most common

The dynamic strain variations in the bridge are calculated based on the following equation:

$$ \Delta \varepsilon = \frac{{\Delta} {\lambda}}{{\lambda_{{\text{o}}} \left( {1 - \rho} \right)}}, $$
(1)

where Δε is the relative strain change, λο is the original transmitted wavelength, which is constant for each FBG strain sensor, Δλ is the instantaneous wavelength shift, and ρ is the photo-elastic constant.

A typical theoretical value of the photo-elastic constant of the FBG sensors at 1550 nm wavelength is ρ = 0.22 [27]. In practise, some deviations may exist. To derive the actual value of the fibre optic sensors installed in the bridge, calibration tests were performed at the laboratory of CSIC, University of Cambridge, UK. In particular, the sensors were batch calibrated at constant temperature of 20 °C on a strain rig, which consists of a 1.6-m rigid aluminium I-beam mounted with a Pi M-414 1PD linear stage, a PIC-863 DC motor controller and clamps specifically design for the 2 mm cable used here. The linear stage has an accuracy of 1 μm giving an overall estimated measurement accuracy of 10 μm. Strain was applied with the linear stage and the corresponding wavelength shift recorded to derive the strain coefficient, p. As expected, the relationship between the normalised wavelength change and strain was linear with a slope p = 0.77, and the derived photo-elastic constant used in Eq. (1) was estimated to be ρ = 1 − p = 0.23.

In this study, only transient dynamic strain variation during a train loading is investigated. For the short period of one train loading, the temperature effect on wavelength variation is negligible [18].

The FBG sensors used in this study are draw tower gratings (DTGs), with 6 μm core diameter, 125 μm cladding, 250 μm Omocer coating, 20% mean reflectivity, FWHM 0.1 to 0.2 nm, and manufactured at wavelength range 1512–1588 nm. A 4-channel sm130 Optical Sensing Interrogator of Micron Optics, Inc. was permanently installed in the bridge, offering 1 kHz sampling rate per sensor, and the achieved signal resolution is ± 2 micro-strain (με). A fibre optic cable of 20 FBG sensors was connected to each one of the four channels, offering 80 FBG sensors in total.

The monitoring system is autonomous and powered by a solar panel. The system data-logger is triggered by train loading and the data for each train passage are transmitted through 4G connection to the main server of the Department of Engineering at the University of Cambridge (see also [19]).

2.2 Monitoring data

Figure 4 offers an overview of the bridge response and the various shapes of FBG signals under the most common train loadings. Figure 4-top presents schematically a typical strain sensor array on the North side of the bridge in the longitudinal direction. Four sensors are numbered, which correspond to the four columns of the graphs in Fig. 4. Locations 1, 2 and 3 correspond to the strain sensor between clamps A1 and A2 (near the arch springing), clamps A3 and A4 (near the quarter span), and A5 and A6 (near the half span), respectively. Location 4 corresponds to the sensor that measures strain of the steel wire which connects the arch springings, measuring the dynamic variation of the arch span.

The sensors are labelled according to (i) the arch number, (ii) the location of the sensor array, and (iii) the name of the adjacent clamps. For instance, the sensor in the longitudinal direction at the South side of Arch 38, which is between clamps A2 and A3, is labelled as 38SA2A3; the sensor in the transverse direction at the West side of Arch 37, which is between clamps T4 and T5, is labelled as 37WT4T5, and so forth. The sensors monitoring the span opening are indicated with the letters SP. For instance, the sensor that measures the opening of the North span of Arch 39 is labelled as 39NSP.

The signals presented in Fig. 4 are for the sensor array on the North side of Arch 37, which is the most severely damaged arch. In particular, the first column corresponds to sensor 37NA1A2, the second to 37NA3A4, the third to 37NA5A6 and the fourth to 37NSP. In all plots, the trains are heading towards the East, using the North track. The first row of graphs corresponds to the most common type of passenger trains with two cars that cross the bridge, which are the Class 144 (red line) and Class 155/158 (black line). The second row corresponds to passenger trains with 3 cars, which are the Class 155/158 (green line), Class 170 (red line) and Class 185 (black line). The third row corresponds to the most common passenger train with 4 cars, which is the Class 222. There are more types of passenger trains crossing the bridge, but these signals represent the vast majority of the data. The most common train type is the 3-carriage Class 185, which represents approximately the 50% of the records.

The graphs in the three first columns and three first rows have the same scale to allow direct comparison. All signals represent trains with approximately the same velocity. Negative values of the relative strain correspond to contraction and positive to expansion. Moving from location 1 to location 3, the response gradually changes from fully compressive at the springings to mostly tensile at the keystone. The passenger trains with 3 cars, and especially the most frequent train Class 185, cause the largest dynamic deformation. In the last column, the 37NSP sensor measures strain variations of the steel wire connecting the springings, and the response tends to be symmetric, with the negative and positive values corresponding to span closing and opening. The last row offers an example of a heavy multi-wagon freight train.

Figure 4 shows that each train has a signature response, which results in signals of characteristic shape in each location. The number of peaks are associated with the number of cars. For instance, for all passenger train events in column 3, the first positive peak corresponds to the moment when the first axle (of the first car) passes above the keystone, the last positive peak corresponds to the last axle (of the last car), and the intermediate double peaks correspond to the intermediate pairs of axles.

2.3 Train classification algorithm and previous results

Preliminary statistical analysis results of the dynamic strain variation per sensor location over a period of two years, between July 2016 and June 2018, are presented in [20]. The analysis is combined with a train classification algorithm that categorizes the signals based on the number of peaks and amplitude. In particular, the paper presents mean value and standard deviation of the maximum and minimum signal peaks that correspond to the 3-carriage Class 185 train heading at both West and East directions. The results are summarised in Fig. 5.

Fig. 5
figure 5

Mean value and standard deviation of the maximum and minimum peaks of FBG signals for 42 sensors underneath Arches 37 and 38, for the Class-185 3-carriage train that is heading East (top graph) and West (bottom graph)

The analysis identifies five locations where the dynamic strain has notably increased over the two year monitoring period. The locations are indicated with dashed boxes, both in Figs. 2 and 5, which correspond to sensors 37NA6A7-37NA7A8, 38NA2A3, 37SA4A5, 37SA7A8, and 38SA2A3. In other words, the dynamic deformation has increased in all four symmetric quarter span locations around the Arch 37–38 pier, where the transverse cracks are observed, together with the South keystone area of the most severely damaged Arch 37.

2.4 The effect of ambient temperature

Apart from these five cases highlighted in Figs. 2 and 5, the response remained consistent in all other sensor locations (see Fig. 5). In particular, a uniform seasonal pattern was observed, which was further studied in a subsequent work by the authors [21]. The work shows that temperature and dynamic strain are linearly related. When temperature increases, the bridge dynamic deformation decreases uniformly in the longitudinal direction. In the transverse direction, where the thermal expansion is not constrained, the decrease is smaller. However, the increase in dynamic strain at the five identified locations is beyond what can be explained by temperature variation.

In general, the measured temperature is used to compensate for the photo-thermal effect affecting the fibre core density and resulting in wavelength shifts. The resulting total strain following this compensation is the strain due to mechanical action, as well as thermal expansion/contraction of the structure. In theory, these two effects could be easily decoupled if the structure had a uniform material with well-known thermal expansion coefficient. However, due to the heterogeneity of material (e.g. historic masonry), complexity in geometry (e.g. barrel vaults, spandrel walls, backfill), and variations in exposure (e.g. due to sun, rain), the structure experiences uneven and localised thermal expansion which cannot be fully taken into account by considering air temperature alone. Furthermore, depending on the environmental exposure of the fibre optic network installed externally on an existing structure, the photo-thermal effect could also unevenly affect the FBG sensors.

Acknowledging these difficulties, this study proposes an alternative damage detection approach, which is based on the study of transient dynamic strain variations at the instant of train loading using Eq. (1). During the short duration of the train loading, temperature variations are negligible, and the ‘shape’ of FBG signals depends exclusively on the dynamic deformation of the bridge, inherently omitting any photo-thermal effect in the fibre. Consequently, differences in signal shapes (e.g. variations in signal amplitude) can be associated with either daily or seasonal effects (e.g. change of dynamic behaviour of the bridge due to thermal expansion/contraction of the bridge) or due to permanent changes that are associated with material degradation.

Decoupling mechanical damage from normal seasonal effect is of critical importance for the development of reliable early-warning structural alert systems. This paper presents how statistical shape analysis can be used to improve train classification, enable automated damage detection considering environmental effects, and study multiple aspects of the structural behaviour.

3 Statistical shape analysis (SSA)

In many applications the size and shape of objects is of interest. For instance, in biology it may be of interest to investigate if there is any difference between the sexes by looking at the size and shape of monkey skulls. Statistical shape analysis (SSA) provides a framework to reason about the size and shape of objects. Other applications of SSA, part from the biology example [28], are in chemistry [29] and image analysis [30]. A thorough introduction to SSA can be found in [24]. In the following, a brief introduction of SSA for use for the Marsh Lane Viaduct data is presented.

First, “shape” has the following definition, as given in [31],

Shape is all the geometrical information that remains when location, scale and rotational effects are removed from an object

that is, a shape is invariant to the Euclidean similarity transformations of translation, scaling and rotation.

To describe a shape, a finite number of points on a shape will be used. These points are referred to as landmarks. A landmark is a location of correspondence on each object that matches between and within populations. For instance, if the objects were outlines of human hands, then landmarks could correspond to the ends of each digit. Landmarks may be selected in a variety of ways corresponding to scientifically meaningful locations, anatomical locations or mathematical locations.

In the following statistical analysis we consider the FBG data from sensors 37NA6A7 and 38NA3A4. Sensor 37NA6A7 is located at the 1/3 of the Arch 37 span at the North side of the bridge (see Fig. 2). This is one of the five locations where an amplification of the dynamic strain has been identified in [20], as discussed in the previous Sect. 38NA3A4 is its symmetric location in Arch 38, at the 1/3 of the arch span (see Fig. 2), where no significant variations in the dynamic strain was observed [20].

For this analysis, the dataset consists of M = 1151 3-carriage train passage events heading East on the viaduct that occurred between July 2016 and March 2019. These events were extracted using the train classification algorithm discussed in [20]. In particular, the number of cars is identified based on the number of peaks of smoothed signals, whereas the direction of the train by cross-correlating signals from different locations. For instance, trains heading East will first trigger the sensors in Arch 37 and consequently in Arch 39. The opposite happens with trains heading towards West.

The M = 1151 dataset used in this analysis comprises a one-day dataset from the preliminary, proof-of-concept FBG installation in July 2016. No temperature sensors were installed in this case. Then, there is absence of monitoring data for about 16 months, followed by datasets in November 2017 and March 2018 from the permanent FBG installation, which includes temperature monitoring. After July 2018, the bridge was continuously monitored, as remote connection and power (solar panel) were established.

Initially, the train classes were “manually” labelled as (i) Class 185 and (ii) Class 170 or 155/158 (non-Class 185 trains). Class 185 trains are heavier and more frequent, and easier to be distinguished. Note that these labels are subject to human error, but are the closest representation of truth, based on the available information from site spotting to distinguish the subclasses of trains.

Consider the single train passage event data recorded at sensor 37NA6A7, presented in Fig. 6. The landmarks that are used correspond to the red triangle points in Fig. 6. The first and last landmarks correspond to the start and end of the train passage event. The procedure to landmark all train passage events was automated for all sensing locations and different train loadings—it involved applying a bandpass filter in the frequency domain of the data, then finding the turning points of the smoothed version. For the specific case of sensors 37NA6A7 and 38NA3A4 and 3-carriage passenger train loading, the ‘shapes’ for the data from both sensing locations consist of 15 landmarks.

Fig. 6
figure 6

FBG sensor data (centred) for a single train passage event at 37NA6A7. Triangle points represent landmarks

We shall mathematically denote the landmarks as a matrix \(X \in {\mathbb{R}}^{15 \times 2}\), where Xj1 and Xj2 correspond to the time (x) and microstrain (y) coordinate of the jth landmark. The matrix X is called a configuration matrix.

Let's denote the train passage event configurations for sensor location 37NA6A7 as \(X_{1}^{{\left( {37} \right)}}\),…, \( X_{M}^{{\left( {37} \right)}}\), where the events are chronologically ordered. Similarly, denote the event configurations for sensor location 38NA3A4 as \(X_{1}^{{\left( {38} \right)}}\),…, \( X_{M}^{{\left( {38} \right)}}\), again ordered.

Before proceeding, we centre the configuration matrices as follows. More precisely, we work with CX as opposed to X, where \(C = I_{15} - \frac{1}{15}1_{15} 1_{15}^{T}\), where \(1_{15}\) is an \(15 \times 1\) vector with all entries one and \(I_{15}\) is the \(15 \times 15\) identity matrix. This transformation centres a shape onto the origin.

In general, one way to compare two shapes is to match their configurations, X onto Y say, as closely as possible up to similarity transformations. This transformation can be done using Ordinary Procrustes Analysis (OPA), which is mathematically achieved by minimising

$$ D_{{{\text{OPA}}}}^{2} \left( {X,Y} \right) = \left\| {Y - \beta_{{{\text{OPA}}}} X{\Gamma } - 1_{15} \gamma^{T} } \right\|^{2} , $$
(2)

where \( \left\| X \right\| = \sqrt {{\text{trace}}\left( {X^{T} X} \right)}\) is the Euclidean norm, Γ is a \(2 \times 2\) rotation matrix, \(\beta_{{{\text{OPA}}}} > 0\) is a scale parameter, and γ is a \(2 \times 1\) location vector. The minimisation is over all possible Γ, β and γ.

The OPA procedure is not appropriate for our application to the Marsh Lane Viaduct data as we now explain. First, the β term scales both time and space (strain) at the same magnitude e.g. a value of \(\beta = 2\) would scale both time and space by a factor of 2. This is not appropriate since the timing of the train passage event and the amount of measured strain should be allowed to vary separately. This suggests using separate scaling terms for time and space. Second, rotation of shapes will again combine time and space—this suggests not using rotations at all. For our application, we will match two configurations, X onto Y by minimising

$$ D^{2} \left( {X,Y} \right) = \left\| {Y - \beta \circ X - 1_{15} \gamma } \right\|^{2} . $$
(3)

Over all \(\beta = \left( {\beta_{1} ,\beta_{2} } \right)^{T} \in {\mathbb{R}}^{2 \times 1}\) and γ. Note that \(\circ\) denotes the Hadamard product, which is an element-wise operation such that for \( \beta \in {\mathbb{R}}^{2 \times 1}\) and \(X \in {\mathbb{R}}^{15 \times 2}\)

$$ \beta \circ X = \left( {\begin{array}{*{20}c} {\begin{array}{*{20}c} {\beta_{1} X_{11} } & {\beta_{2} X_{12} } \\ {\beta_{1} X_{21} } & {\beta_{2} X_{22} } \\ \end{array} } \\ {\begin{array}{*{20}c} \vdots & \vdots \\ {\beta_{1} X_{15,1} } & {\beta_{2} X_{15,2} } \\ \end{array} } \\ \end{array} } \right). $$
(4)

The minimisation of Eq. (3) is straightforward as the optimal values of β and γ can be computed analytically (closed form). We shall denote the optimal values of \( \beta = \left( {\beta_{1} ,\beta_{2} } \right)^{T}\) as \(\hat{\beta } = \left( {\hat{\beta }_{{1}} ,\hat{\beta }_{{2}} } \right)^{T}\), which are given by

$$ \hat{\beta } = \left( {\frac{{\mathop \sum \nolimits_{j = 1}^{15} X_{j1} Y_{j1} }}{{\mathop \sum \nolimits_{j = 1}^{15} X_{j1}^{2} }},\frac{{\mathop \sum \nolimits_{j = 1}^{15} X_{j2} Y_{j2} }}{{\mathop \sum \nolimits_{j = 1}^{15} X_{j2}^{2} }}} \right)^{T} . $$
(5)

Note that the estimate of the location vector γ is not of interest. To illustrate the difference between these two types of transformation, consider transforming an event, say event B depicted in Fig. 7 top-right, onto event A (Fig. 7 top-left). The transformation of B onto A that minimises \( D_{{{\text{OPA}}}}^{2}\) is depicted in Fig. 7 bottom-left. Notice that the transformation of B does not correspond well to A. The total Euclidean distance between the landmarks is 256.3. The estimated transformation parameters are

$$ {\hat{\Gamma }} = \left[ {\begin{array}{*{20}c} {0.99999474} & { - 0.00324237} \\ {0.00324237} & {0.99999474} \\ \end{array} } \right]\quad {\text{and}}\quad \hat{\beta }_{{{\text{OPA}}}} = 0.9561869 $$
(6)
Fig. 7
figure 7

Illustration of transformation B (top right) onto A (top left). Transformed version of B using \(D_{{{\text{OPA}}}}^{2}\) and \(D^{2}\) are presented bottom left and right respectively

The transformation of B onto A that uses our modified version of OPA is depicted in Fig. 7 bottom-right. The separate scaling factors in D2 allows event B to be mapped onto A with visually good results. The total Euclidean distance between the landmarks for this transformation is 144.2. The estimated scaling parameters are \( \hat{\beta } = \left( {0.3159522, 0.9771222} \right)^{T}\). Note that the estimated train speed of events A and B is 17.45 m/s and 5.11 m/s. These velocities are estimated by finding the time-lag between signals through cross-correlation from corresponding sensing locations of Arch 37 and 39, where their mean distance in 17.1 m. The \(\hat{\beta }_{{1}} = 0.3159522\) term is the estimated ratio of the velocities based on mapping the shapes, which is close to \(5.11/17.45 = 0.29\). In general, the \(\hat{\beta }_{{1}}\) term is the estimated ratio of train velocities, and similarly, the \(\hat{\beta }_{{2}}\) term is the estimated ratio of the landmark strain values.

4 Results

4.1 Data overview and trends identification

Returning to our application, consider mapping each configuration, from a particular sensor, onto the first chronological event by minimising Eq. (3). More precisely, minimise \(D^{2} \left( {X_{1} ,X_{j} } \right)\) for \(j = 1, \ldots ,M\), to obtain the optimal values \( \hat{\beta }\left( 1 \right), \ldots ,\hat{\beta }\left( M \right)\). Note that the interpretation of \( \hat{\beta }_{1} \left( j \right)\) and \(\hat{\beta }_{2} \left( j \right)\) is as follows: \(\hat{\beta }_{1} \left( j \right)\) is the amount of time-scaling required to map the first chronological event onto the jth train passage event. Similarly,\( \hat{\beta }_{2} \left( j \right)\) is the amount of strain-scaling required to map the first chronological event onto the jth train passage event. We choose to transform all events onto the first chronologically, as it represents the state of the bridge at the beginning of the monitoring period, in July 2016. Thus, any observed changes in the transformation may be attributable to the passage of time.

Figure 8 displays the estimates of \(\hat{\beta }_{1}\) and \(\hat{\beta }_{2}\) for all train passage events at sensor location 37NA6A7. The data from the preliminary FBG installation in July 2016 are shown separately from the permanent FBG installation, after November 2017. Data from the preliminary installation do not contain information from temperature sensors, as in the permanent installation. Note that there is a 16 months gap between the two installations. The data from the permanent installation, after November 2017, have been grouped to Class 185 and non-Class 185 (i.e. Class 170 or Class 155/158), as discussed previously. There is not a discernible correlation between \(\hat{\beta }_{1}\) and \(\hat{\beta }_{2}\). In other words, train speeds (over the limited range of speeds measured) do not seem to affect the level of dynamic deformation at 37NA6A7. Similar observations can be made in the other sensor locations.

Fig. 8
figure 8

Estimates of \(\hat{\beta }_{1}\) and \(\hat{\beta }_{2}\) for train passage events at sensor 37NA6A7

Figure 9 presents the \(\hat{\beta }_{1}\) and \(\hat{\beta }_{2}\) values for the train passage events against time at sensor 37NA6A7. There seems to be a trend of the \(\hat{\beta }_{2}\) values over time. On the other hand, there is no pattern in the \(\hat{\beta }_{1}\) estimates over time, which is logical since trains speeds are essentially random.

Fig. 9
figure 9

Estimates of \(\hat{\beta }_{1}\) (top row) and \(\hat{\beta }_{2}\) (bottom row) for chronologically ordered train passage events at sensor 37NA6A7

The \(\hat{\beta }_{2} \) values are plotted against temperature in Fig. 10. Note that the first 31 train passage events did not have recorded temperature readings. We generate temperature data for these events, using information from the National Centre for Atmospheric Science database, available online (https://sci.ncas.ac.uk/leedsweather/Archive/). The sampling rate in the database is 5 min, and the data were interpolated to estimate in-between temperature values. These events with “indicative” temperature values are included in Fig. 10.

Fig. 10
figure 10

Estimates of \(\hat{\beta }_{2}\) against temperature for train passage events at sensor 37NA6A7

It is expected that most of the 31 train passage events (2016 preliminary data) correspond to the most frequent Class 185 train. With this in mind, there is a clear shift in \(\hat{\beta }_{2}\) between the 2016 preliminary events and the remaining values. This will be examined separately in Sect. 4.3. By omitting for now the 2016 preliminary data, in Fig. 10 we see a clear trend, with \(\hat{\beta }_{2}\) decreasing with temperature. Further, there seems to be various levels of the trend – more precisely, there seems to be a three different linear relationships exhibited in the data. The statistical relationship between \(\hat{\beta }_{2}\) and temperature may be modelled by for example a linear model. However, if the different groups are not accounted for in the data will lead to inflated variation in the parameter estimates. Therefore, before fitting a model, we start by investigating the different groups exhibited by the data. The groups seem to correspond to a large extent with the different train classes. Before modelling the data, we first consider constructing a classifier that can be used to predict train classes based on the \(\hat{\beta }_{2}\) value and temperature of an event.

4.2 Classification of train passage events

As suggested by Fig. 10, there may be a way to classify train classes using temperature and the corresponding \(\hat{\beta }_{2}\) value. For this classification problem, a support vector machine (SVM) [30] is used. Given the visible linear division between the train classes, i.e. the data are nearly linearly separable for the train classes, a support vector machine with a linear kernel would be appropriate. Note that, if the data were not linearly separable, then other types of kernels may be appropriate—see Sect. 5.4 in [32]. A support vector machine with a linear kernel for this 2-dimensional (temperature and \(\hat{\beta }_{2}\)) problem would seek to find a line furthest from the data in the sense that it separates train Class 185 points from the train Class 155/158–170 points and maximizes the distance to the closest point. Since the train classes are not linearly separable, then one allows overlap between the groups. This line (and more generally in higher dimensions) is called the maximum margin hyperplane.

A SVM is trained using randomly selected 100 datapoints across the two different classes of train. Note that the events with uniformly generated temperature values were not used, since the generated temperature data are not directly comparable with the data acquired by the sensor network. These training datapoints are indicated by the blue circles in Fig. 11. The resulting partition of the data-space, indicated by the dashed black line, separates the two classes very well. The classification results are summarised in Table 1.

Fig. 11
figure 11

SVM results. Maximum margin hyperplane (dashed black line), training datapoints

Table 1 Classification performance of SVM fitted to \(\hat{\beta }_{2}\) and temperature data

Table 1 states that the fitted SVM predicted 850 Class-185 trains and 266 of other class (these numbers correspond to the rows of Table 1). Of the 850 Class-185 train predictions, 10 did not agree with the labels. Similarly, of the 266 Class 155/158 or 170 trains, only 2 did not agree with the labels. Clearly the fitted SVM is a very good classifier of train types.

Note that the 155/158–170 train class consist of two subclasses of trains, for which we have limited information. More precisely, we only have data from three instances where site spotting distinguished between the 155/158 and 170 trains (see Fig. 12); there are two events corresponding to Class 170 trains and one for a Class 155/158 train. Though the data is limited, the data points do fall in the expected clusters based on train weights.

Fig. 12
figure 12

Estimates of \(\hat{\beta }_{2}\) against temperature for train passage events at sensor 37NA6A7, with train subclasses

To train an SVM for this 3 class problem, training points are artificially labelled corresponding to train subclasses. This labelling is done to (i) demonstrate that SVMs can be easily extended to multiple class problems and (ii) propose a statistical model explaining the relationship between \(\hat{\beta }_{2}\) and temperature that takes into account the different classes of train. Note, however, that we do not subscribe to this artificial labelling of the subclasses as the true train types. Further investigation into the labelling of subclasses is warranted to verify the suggested labels. The training points are indicated in Fig. 13 by the black, green and blue points. The maximum margin hyperplanes are indicated by the dashed black lines.

Fig. 13
figure 13

SVM results using train subclasses. Data are represented by the grey points. Circled points are those datapoints used to train the SVM. The dotted lines represent the hyperplanes that separate the class predictions of the SVM method

Let's assume that these predicted labels from this SVM are true labels (see Fig. 14-top for the predicted labels) and proceed to fit the following linear model with the 3 classes of train:

$$ Y_{i} = \hat{\mu } + \hat{\alpha }_{1} {\text{temp}}_{i} + \hat{\alpha }_{2} V_{1,i} + \hat{\alpha }_{3} V_{2,i} + \varepsilon_{i} , \ldots ,\varepsilon_{i} \sim N\left( {0,\hat{\sigma }^{2} } \right), $$
(7)

where the \(\varepsilon_{i}\) are independent, \({\text{temp}}_{i}\) is the temperature reading from the ith train event, \(Y_{i}\) is the random variable corresponding to \(\hat{\beta }_{2}\) for the ith train event, and \(V_{1,i}\) and \(V_{2,i}\) are dummy variables taking the value 0 or 1 depending on the train class. Specifically, if train i has predicted label 1 corresponding to train class 185, then \(V_{1,i} = 0\) and \(V_{2,i} = 0\); for Class 170, \(V_{1,i} = 1\) and \(V_{2,i} = 0\); for Class 155/158, \(V_{1,i} = 1\) and \(V_{2,i} = 1\). The estimated values of the unknown parameters are given in Table 2 (left columns):

Fig. 14
figure 14

Top: Expected \(\hat{\beta }_{2}\) value (lines) given temperature and train class given by the fitted linear model of Eq. (7) using data from sensor location 37NA6A7. The train classes are the predicted classes given by the SVM. Bottom: Fitted linear model to estimates of \(\hat{\beta }_{2}\), from 37NA6A7, and temperature data with train classes from sensor location 38NA3A4

Table 2 Linear model parameters estimates for Eq. (7) using data from (i) sensor location 37NA6A7 (left columns), and (ii) using data from 38NA3A4 and predicted class types based on information for 37NA6A7 (right columns)

This model gives an R2 value (a measure of goodness-of-fit) of 95%. The fitted model, i.e. the expected value of \(\hat{\beta }_{2}\) given the temperature and train class, is represented in Fig. 14-top.

The interpretation of the estimated parameters of fitted model is as follows: There is an \(\hat{\alpha }_{2} = 0.396\) decrease in relative strain between Class 185 and 170 train events; further there is an \(\hat{\alpha }_{3} = 0.6655\) decrease from Class 185 to Class 155/158 train events.

Now consider the data from Arch 38 processed in a similar manner to the above. Figure 14-bottom presents the \(\hat{\beta }_{2}\) values against temperature for the data from sensor 38NA3A4 where the colours denote the predicted classes based on data from sensor 37NA6A7. Note that the train classes again seem separable. We again fit the model of Eq. (7) to the data presented in Fig. 14-bottom. The estimated parameter values for this model are given in Table 2 (right columns).

4.3 SSA for damage detection

Figure 14 shows that SSA, combined with a trained SVM, offers a sensitive and robust way to classify train events and obtain expected behaviour curves for each sensor location. This expected behaviour can take into account multiple parameters that show a trend on the statistical data. In this study, the temperature was shown to affect the dynamic deformation of the bridge in a consistent linear manner for all predicted train classes.

The expected behaviour curves per sensor location can be now used for damage detection. An example is offered below.

Figure 10 shows a deviation between the 2016 preliminary installation data (indicative temperature) and the events recorded from the permanent installation (measured temperature) for sensor 37NA6A7. Figure 15 compares the deviation between the \(\hat{\beta }_{2}\) estimates of the 2016 preliminary events and the predicted behaviour curve that corresponds to the most frequent Class 185 train for 2017–2019, for all four 1/3-span locations of the North side of Arches 37 and 38. Apart from the 37NA6A7 and 38NA3A4 locations, which were studied in the previous section, two more sensors, 37NA3A4 and 38NA6A7, have been added in this example. The term \({\Delta }\hat{\beta }_{2}^{\rm T}\), shown in the 37NA6A7 plot, expresses the deviation between the mean \(\hat{\beta }_{{2,{\text{prelim}}}}^{\rm T}\) value of the preliminary events and the predicted behaviour curve \(\hat{\beta }_{{2,{\text{predict}}}}^{\rm T}\) for the given temperature range, T = [17–20 °C]. Clearly, the location that shows the larger deviation \({\Delta }\hat{\beta }_{2}^{\rm T}\), is sensor 37NA6A7.

Fig. 15
figure 15

Comparison between the \(\hat{\beta }_{2}\) estimates for the 2016 preliminary events and the predicted Class-185 curve from SSA and SVM in four sensor locations

Table 3 offers the estimated \({\Delta }\hat{\beta }_{2}^{\rm T}\) values for all four locations, together with the peak-to-peak strain amplitude of the first chronological (reference) signal \(\overline{{\varepsilon_{0} }}\), where by definition, \(\hat{\beta }_{2} = 1\). Table 3 also includes values for the term \(\overline{\varepsilon }_{{{\text{prelim}}}}^{T}\), which is the mean peak-to-peak strain amplitude of the 2016 preliminary results, and the term \(\overline{\varepsilon }_{{{\text{predict}}}}^{T}\), which is the mean peak-to-peak strain amplitude of the predicted curve for the given temperature range. The expression that relates the above terms is the following:

$$ {\Delta }\hat{\beta }_{2}^{\rm T} = \hat{\beta }_{{2,{\text{predict}}}}^{\rm T} - \hat{\beta }_{{2,{\text{prelim}}}}^{\rm T} \approx \frac{{\overline{\varepsilon }_{{{\text{predict}}}}^{T} - \overline{\varepsilon }_{{{\text{prelim}}}}^{T} }}{{\overline{{\varepsilon_{0} }} }} $$
(8)
Table 3 Peak-to-peak strain amplitude of the first chronological signal \(\overline{\varepsilon }_{{\text{o}}}\), and deviation of \(\hat{\beta }_{2}\) between the 2016 preliminary events and the predicted Class-185 curve from SSA and SVM in four sensor locations

For instance, in the most severely damaged arch of the bridge, Arch 37, location 37NA3A4 shows a slight deviation between the preliminary and predicted \(\hat{\beta }_{2}\) equal to 0.095, or 9.5%, which corresponds to change in the peak-to-peak amplitude dynamic strain equal to \(9.5{\text{\% }}\overline{{\varepsilon_{0} }} = 6.75\;\mu {\upvarepsilon } \approx 78.296\;\mu {\upvarepsilon } - 72.071\;\mu {\upvarepsilon }\). In other words, the mean peak-to-peak amplitude of the recorded events from the permanent installation, appears amplified by approximately 7 με. Note that this identified change is in an order comparable to the low signal noise of the FBGs (± 2 με). Slight deviations may exist due to incompatibility of the temperature data, or because not all the preliminary events are from Class 185 trains. In another location of Arch 37, 37NA6A7, the biggest \(\hat{\beta }_{2}\) deviation among the four sensors is observed, equal to 0.682, or 68.2%, while the peak-to-peak amplitude appears amplified by \(68.2{\text{\% }}\overline{\varepsilon }_{o} = 31.89\;\mu {\upvarepsilon }\). This observation is in agreement with previous findings [20].

It is noted that some dispersion of data points around the fitted curves in Fig. 14 is expected. This is because not all trains of the same class are identical, and their exact weight may vary depending on the number and distribution of passengers. Other seasonal variations and environmental factors not considered in this study might also have caused some data dispersion. Consequently, any identified changes within this data dispersion range is not a strong indication of damage, which is the case of sensor 37NA3A4. On the other hand, any clear deviation that cannot be explained by usual data dispersion, indicates change in structural behaviour, such as sensor 37NA6A7. Moreover, structural analysis could help to further assess the effect of the observed strain variations on masonry and overall structural performance.

In the same way that the 2016 data where compared with the expected behaviour curves derived from the SSA, future projections are possible, by simply replacing the \(\hat{\beta }_{{2,{\text{prelim}}}}^{\rm T}\) with the new \(\hat{\beta }_{2}^{\rm T}\) for a given temperature range and rewrite Eq. (8) as: \({\Delta }\hat{\beta }_{2}^{\rm T} = \hat{\beta }_{{2,{\text{new}}}}^{\rm T} - \hat{\beta }_{{2,{\text{predict}}}}^{\rm T}\).

These results enable a much finer level of deterioration monitoring. High-sensitivity distributed sensing approaches, such as fibre optic networks, either embedded or externally installed on new or existing structures, may provide ‘good’ quality data, which is indispensable for ‘good’ data analytics. The FBG monitoring system installed in the Marsh Lane Bridge, has low signal noise and proved to be resilient over time. Statistical procedures of low computational cost that are able to cope with streaming monitoring data, may detrend the data from selective parameters (e.g. train speed and train loading variations, environmental effect) to allow a more accurate monitoring of structural deterioration.

To this end, it is clarified that although the proposed system appears sensitive in detecting changes in the dynamic deformation of the structure, damage assessment based on this data is not a straightforward process. This work and previous studies of the authors [19,20,21] have shown that with the exemption of a few locations, the overall structural response of the bridge has not experienced significant permanent changes. This finding indicates that the extensive retrofitting intervention back in 2015 with the installation of steel ties and the filling of the relieving arches with concrete successfully delayed deterioration of the Marsh Lane Bridge during the following 3-year monitoring period. However, the results also clearly define locations where local deformations are increasing and minor damage is occurring. The availability of continuous long-term monitoring data will be a key element for continued tracking of gradual changes at these locations as well as across the whole bridge, to identify when this minor damage or other new damage may become an issue. Multi-sensing information, for instance the combination of acoustic emission sensors for local cracking monitoring with distributed strain monitoring approaches, together with statistics, may also contribute to this effort, enabling a more detailed damage mode identification and potentially automated structural performance assessment of ageing infrastructure [33].

5 Conclusions

In this work, the deterioration monitoring of a masonry arch bridge with a fibre Bragg grating (FBG) sensors network is enhanced with the statistical shape analysis (SSA). The method is appealing for railway bridge monitoring where repeated loading from specific train type corresponds to a signature response of the structure, which results in FBG signals of similar ‘shape’.

In particular, the ordinary procrustes analysis (OPA) was adapted to use two separate scaling parameters \(\hat{\beta }_{{1}}\) and \(\hat{\beta }_{{2}}\) which can be directly linked to train speed and dynamic strain. As a result, the method provides a computationally inexpensive tool to directly evaluate the impact of train speed on the dynamic deformation of the bridge, which is a major concern for railway infrastructure owners. For Marsh Lane Bridge, the method indicates that train speed has no effect on deformation, although the range of train speeds is limited to around 5–18 m/s (20–65 km/h) because the bridge is located near the Leeds central railway station.

Secondly, a clear correlation between \(\hat{\beta }_{{2}}\) (dynamic deformation) and ambient temperature was identified. In particular, the dynamic strain appears to decrease with temperature increase, verifying previous observations [21]. This can be related to thermal expansion of the bridge, which further confines its structural members during warmer periods.

The identified linear pattern between dynamic strain and temperature was found to be consistent regardless of the train loading. Plotting \(\hat{\beta }_{{2}}\) against temperature revealed distinguishable trends from three subclasses of 3-carriege trains, which cause subtle differences in the response of the bridge, and were hardly separable before. An automated train classification process is presented, after training a support vector machine (SVM) to identify the three different classes of trains. Train classification appears very accurate.

The paper shows that SSA, combined with a trained SVM, offers a robust way to obtain expected behaviour curves for each sensor location, by taking into account multiple parameters (in this study, only train type and temperature variations were considered). These behaviour curves can be used for a more accurate deterioration monitoring.

An example is offered, where monitoring data from a preliminary FBG installation in 2016, a few months after the retrofitting of the bridge, are compared with the main set of the data from the permanent FBG installation. Changes in the dynamic deformation are expressed as differences in the \(\hat{\beta }_{{2}}\) values. It is shown that, even tiny dynamic strain variations, comparable to the low noise of FBG signals (± 2 με), result in clearly identifiable \(\hat{\beta }_{{2}}\) shifts.

The introduction of the SSA has significantly improved the clarity of the results and accuracy in deterioration detection. High-sensitivity sensors networks, combined with statistical modelling, may leverage automated railway asset management, even for ageing masonry structures, the structural assessment of which remains a major challenge for civil engineers.