Spatio-temporal analysis of georeferenced time-series applied to structural monitoring

Abstract

Spatio-temporal (S-T) analysis is not typical in structural monitoring applications of buildings and infrastructure. However, monitoring always includes the temporal component, and observations are often captured in specific locations. In other words, a monitoring dataset could also be considered a spatio-temporal archive, notwithstanding that not all monitoring applications can benefit from S-T processing methods. The paper discusses spatio-temporal analysis using the structural monitoring dataset of the Cathedral of Milan, which has an archive of vertical settlements collected from more than 50 years of measurements. The proposed methods can be adapted and extended for other structural monitoring applications, including single buildings, infrastructure, and the environmental level. The cases of pure temporal (T) and spatial (S) analyses are also discussed, comparing the different approaches, illustrating the pros and cons, and describing the opportunities of the S-T combined workflow. The paper specifically focuses on different typologies of S-T processing: data visualization and exploration techniques, clustering, change detection, prediction, and forecasting. The proposed algorithms were all implemented within the R open-source programming language. They can be replicated (and adapted) for other structural monitoring datasets featuring spatio-temporal correlation.

1 Introduction

Structural monitoring implicitly includes the temporal domain. An example is the acquisition of observations at specific time intervals and their comparison [1]. Nowadays, observations can be acquired with several methods and sensors, and the commercial market offers different solutions depending on the phenomenon under investigation. However, data collection is just an aspect of monitoring applications. According to [2], data processing and feature extraction are necessary to convert raw data into damage information, because sensors do not directly measure damage.

Historical heritage plays a fundamental role in various structural monitoring applications for the specific characteristics of the structures [3]. A multidisciplinary approach involving different specialists is essential to face the different aspects and preserve heritage buildings and sites. A review of relevant monitoring examples is illustrated and discussed in [4].

According to [5], in-depth knowledge is always essential for conserving a monument. The word “knowledge” covers different topics, such as the history of the monument and its complex geometry, the different stratifications and transformations that occurred over time, materials and construction technologies, decay patterns, structural problems, previous restoration activities, and many other aspects. This means that information is not necessarily limited to numerical observations acquired with sensors or other specific measuring tools. Information indicates heterogeneous data and the interpretative work of specialized staff must relate various sources, among which quantitative numeric data and other materials in different formats. Then, monitoring supported by historical research and experimental data is fundamental for the conservation of a monument and for developing intervention strategies [6, 7].

In this respect, the availability of time-series, the use of processing algorithms to generate new information from raw data, and the graphic representations of the results can be used by the specialists involved in the monitoring project.

The manuscript focuses on those monitoring applications that include the (geo)spatial dimension. The case of periodic monitoring using geomatics methods provides input datasets for data processing. Generally speaking, the observed phenomena always have a spatial dimension. They occur in specific points, directions, areas, and structural elements. Monitoring observations can be captured with sensors installed in specific positions, obtaining time-series at given spatial locations. The approach proposed here is based on spatio-temporal (S-T) analysis, assuming that the behavior at spatial locations at time t will almost certainly affect the behavior at nearby spatial locations at time $t+1$ [8].

Some applications can benefit from processing methods capable of combining the time domain and the spatial domain using geostatistical tools, especially those based on long-term monitoring datasets, in which the temporal and (geospatial) domains can be simultaneously combined to generate additional information with automated or semi-automated processing algorithms. According to the author, the combined use of S-T observations could offer novel alternative solutions for data processing and visualization in various applications characterized by spatial and temporal correlation.

The proposed methods can be applied to data captured at the level of a single building (like in the case study used in the manuscript) or extended toward infrastructures and environmental applications. The transition to a wider area probably makes the methods more intuitive than acquiring specific data at the building level. Examples of monitoring applications at the environmental level are those based on synthetic-aperture radar (SAR) [9,10,11,12,13], global navigation satellite system (GNSS) [14,15,16,17], or other monitoring data acquired with geomatics instruments [18,19,20], such as drones, total stations, and levels. Examples with datasets acquired at a more local scale are described in [21,22,23,24,25,26] and focus not only on historical buildings but also bridges, dams, archeological sites, etc.

GIS software could appear as an ideal processing and visualization environment. Tobler’s first law of geography states: “everything is related to everything else, but near things are more related than distant things” [27]. However, most tools implemented in GIS software do not fully exploit S-T information. Observations acquired at different times are often analyzed independently, reducing the S-T dataset to a sequence of independent spatial processes. In this case, data processing is based on spatial analysis methods, i.e., neglecting temporal correlations [28].

Another approach could consider the different time-series as spatially independent, i.e., ignoring spatial correlation. In this case, processing techniques based on univariate or multivariate time-series analysis can be exploited. [29, 30].

A complete spatio-temporal analysis, on the other hand, must consider correlation in space and time without splitting the dataset into “snapshots” at a specific time [13, 31]. In other words, S-T analysis extends the previous concepts considering proximity in both space and time.

This work aims to analyze monitoring observations acquired at specific periods (static and quasi-static monitoring), producing graphic representations, and supporting the interpretation using S-T techniques. Methods for spatial (S), temporal (T), or combined S-T analysis will be compared to understand the pros and cons of a complete S-T approach.

Data processing was carried out with the R open-source programming language, which has advanced tools for geostatistics [32] and allows the implementation and validation of S-T processing methods here extended to structural monitoring datasets. R runs on different UNIX platforms, Windows and MacOS, and it can be downloaded from https://www.r-project.org/. R was used together with the RStudio integrated development environment (IDE) https://posit.co/products/open-source/rstudio/.

The paper is structured as follows: Sect. 2 provides a methodological overview of the proposed methods and describes the case study used in the manuscript to illustrate the implemented S-T algorithms with numerical examples. Section 3 shortly describes the creation of the digital archive for the Cathedral of Milan (the case study presented in the paper), including the origin of the monitoring system and the acquisition of new observations. Section 4 illustrates visualization methods for S-T data. Section 5 describes time-series clustering methods, showing that columns can be grouped depending on their relative displacements. Section 6 extends the temporal domain and introduces techniques for detecting anomalous values in new monitoring campaigns. Finally, Sect. 7 considers spatio-temporal prediction and advanced methods for visualization in space and time.

2 Overview of the proposed method

The manuscript describes spatio-temporal (S-T) analysis and its application to georeferenced time-series for structural monitoring. The temporal (T) and (S) spatial analysis cases are also considered and discussed. Results are compared to clarify when and how an S-T can provide more benefits than other traditional analysis methods.

The structural monitoring application illustrated and discussed in this manuscript closely relates to spatio-temporal data processing, information extraction, and visualization. The paper discusses the static monitoring system that tracks vertical displacements due to subsidence in the Cathedral of Milan (Duomo di Milano). The author does not intend to provide a complete and detailed description of the events that led to the restoration project of the columns of the Cathedral, for which the reader is referred to specific textbooks [33]. The manuscript aims to show the pros and cons of S-T analysis applied to a vertical settlement dataset, including the creation and validation of the dataset and the extraction of additional information.

A description of the different S-T processing algorithms applied to a specific case study can better clarify the outcomes of the different sections, which mainly describe S-T analysis and more traditional temporal and spatial methods. Thus, a (brief) overview of the historical events is useful to present the implemented processing methods with numerical examples from a real case study, providing insight not limited to theoretical aspects.

2.1 Methodological workflow

The structure of the paper follows the workflow shown in Fig. 1, which also corresponds to the different sections of the manuscript. The input is the spatio-temporal dataset of vertical settlements measured in spatial locations $\left\{ \textbf{s}_1\textbf{s}_2,...,\textbf{s}_C\right\}$ at (regular) time intervals $\left\{ t_1,t_2,...,t_T\right\}$. The algorithms can be adapted and extended to other datasets of spatially (geo)referenced time-series. The workflow has an internal loop in which the input dataset can be integrated and updated. Indeed, the data collected can be affected by specific problems, such as anomalous values (outliers) or missing data in specific locations at specific epochs. In the case of the example illustrated in this work (Sect. 3), some discontinuities were also found for the change of instruments and reference points. One of the aims of the proposed work is the identification of anomalous values and the recovery of missing information. Regenerating an incomplete archive is not a trivial operation, and future users must be aware of the technical operations carried out on the dataset.

Fig. 1

Schematic workflow of the different algorithms used in the manuscript

Another important topic is related to the visualization of spatio-temporal data (Sect. 4). Simultaneous visualization and exploitation could be challenging in the case of multiple (geo)spatially referenced time-series containing monitoring measurements. Video animations (i.e., short clips) are very effective but require additional devices such as computers and monitors. Similar considerations apply to interactive graphic representations, in which the user can dynamically navigate a virtual scene, changing parameters, and settings.

Different techniques also exist for generating efficient “static” visualizations of S-T information. They were mainly developed for applications at the environmental level, combining cartographic locations, and the variations of the considered variable. The paper exploits mainly 2D graphic representations, which can be included in technical reports and could activate cognitive processes facilitating the interpretative work of the reader.

Effective visualization is fundamental for the preliminary inspection of the input dataset as well as the production of graphics for additional information (output) produced with further data processing. The implemented geostatistical processing methods are described in the remaining different sections. Given a spatio-temporal process $\left\{ z(\textbf{s},t)\right\}$ with data sampled in spatial locations $\textbf{s}$ and time intervals t, the processing methods discussed in this work can be classified into three main categories:

• Temporal data processing (T), in which spatial correlation is ignored resulting in a pure collection of time-series $z(\textbf{s}_1,t),...,z(\textbf{s}_C,t)$. This does not mean that multiple time-series cannot be simultaneously used in the analysis. Both univariate and multivariate time-series processing methods can be implemented. However, proximity in space is not considered in the analysis;

• Spatial data processing (S), in which proximity in space is taken into consideration and the dataset is split into a set $z(\textbf{s},t_1),...,z(\textbf{s},t_T)$;

• Spatio-temporal processing (S-T), which simultaneously considers correlations in space and time, and the whole dataset $\left\{ z(\textbf{s},t)\right\}$ is used. Of course, this is also the most demanding type of processing.

The manuscript provides different examples in which different strategies (T, S, or S-T) are used to solve the same problem. Although S-T methods could be interpreted as the most general solution, the results achieved by simplified S or T approaches can be sufficient in several practical applications. Fine-tuning the parameters of a pure S-T approach is also more complicated for the intrinsic difficulty in defining the joint spatio-temporal dependence.

Another advantage is the opportunity to process the data with different methods, cross-validating the output. Uncertainty quantification cannot be neglected in the case of displacement monitoring. Uncertainties are associated with input measurements, adjusted quantities using least squares (elevation values), and other quantities derived from additional processing with the proposed methods for information extraction. Quantification of uncertainty in structural monitoring has a direct relationship with the characteristics of the phenomenon under investigation. The estimated value must have sufficiently small uncertainty to highlight changes.

The processing phase described in the paper can be divided into two operations: (1) time-series clustering and (2) time-series prediction/forecasting. Clustering aims at creating groups of time-series based on some similarity measures. As the different time-series correspond to vertical displacements of the different columns caused by subsidence, discovering meaningful groups of columns can identify specific areas inside the Cathedral with similar behavior. The number of clusters is unknown and several methods are available to compute the similarity measures. The examples reported in the next sections will illustrate how to deal with relative (differential) displacements, which is often the case of structural monitoring datasets in which conventional choices are unavailable to materialize a reference system and to compare measurements acquired at different times.

The second main topic is time-series prediction/forecasting, which is related to spatio-temporal statistical modeling and has three principal goals, as described in [31]:

• 1. predicting a plausible value at some location in space within the time span of the observations and reporting the uncertainty of that prediction;

• 2. performing scientific inference about the importance of covariates in the presence of spatio-temporal dependence; and

• 3. forecasting the future value at some location along with the uncertainty of that forecast.

Prediction is here considered as an interpolation to recover missing values in specific sampling locations at a specific time, for instance when some points cannot be accessed during an on-site monitoring campaign. More in general, other possible applications are the estimation of discontinuities for the loss of a monitoring point, which is then replaced with a new one and there is a jump between the measurements. Time-series prediction allows the calculation of values also at unsampled locations in space, which is useful for creating visual representations of vertical settlements. Prediction can recover missing values from incomplete and noisy data in space and time. Processing in the paper considered deterministic methods (without providing estimates of the prediction uncertainty) and kriging, which instead considers the covariability between any two space–time locations.

Forecasting instead focuses on future values and is mainly used as an anomaly detection tool in additional monitoring campaigns. Only the short term is considered (e.g., 1- or 2-step ahead forecasting), so that new values acquired on-site can be compared to the forecasted values using a threshold that depends on the uncertainty of forecasting and precision of adjusted observations.

Several other possible applications can be intended as subtasks for the processing methods shown in the workflow, notwithstanding the distinction between the different levels is not rigid. For instance, anomaly detection could be intended as a subtask of T, S, and S-T analysis. However, anomaly detection is a wide topic and can also be intended as an independent task. The author’s opinion is not to introduce rigid boundaries for the definition of the different operations with spatio-temporal structural monitoring datasets. Other practical applications with different S-T monitoring datasets could require operations depending on the specific research problems and the data under investigation.

2.2 Description of the case study

Figure 2 shows the nave and aisles of the Duomo di Milano and some of the monitoring points installed on the columns. Starting from the second half of the twentieth century, groundwater extraction in Milan’s center became more intense. Between 1956 and 1970, the water table was lowered by 20 m [34], i.e., an average velocity of 1.5 m/year. The lowering of the groundwater level resulted in movements of the ground for the entire center of Milan, breaking a secular equilibrium. The overall movement is not constant throughout the historic center and features differential movements at the local level [35].

The rapid lowering of the groundwater level caused differential movements between the columns of the Cathedral and redistributed internal forces, damaging some architectural elements such as the columns and vaults [36]. The static condition became critical in the 1960s, as confirmed by the rapid development of the crack pattern of the columns, with also detachments of material. The differential settlements of the columns brought the Cathedral very close to collapse, making it necessary to build a jacket of reinforced concrete (in 1969) around the four columns of the tiburium, with a thickness of 0.35 m. The reinforcement avoided the collapse of the columns and allowed a period of analysis and development of the restoration strategy [37].

Fig. 2

The columns of the nave and aisles of the Duomo di Milano. The image was acquired from the pulpit in June 2022. The location of a few monitoring points is also shown

The columns of the Duomo are composed of heterogeneous materials. The external ring is made up of Candoglia marble and has a variable thickness between 0.2 m and 0.9 m, with an average value of about 0.5 m. Lime mortar binds the internal core composed of irregular blocks of serizzo, granite, pieces of marble, and bricks. There is a significant difference between the material properties of the external layer and the internal core of the column.

After a period of study and analysis, the restoration was carried out starting from the top of the column, progressively removing the jacketing, the damaged marble ashlars, and the heterogeneous materials. Voids were filled with new marble ashlars and cement mortar. New marble blocks were cut with millimeter-level accuracy depending on the size of the voids, making new ashlars available on-site after a few hours. The reader is referred to [33] for more details and a complete description of the interventions, which is here briefly introduced.

The four columns of the tiburium (n. 74, 75, 84, 85) were restored between 1981 and 1984, after restoring 21 sary columns of the transept, apse, and central nave in a period of about 10 years. Thanks to a prefectural order of 1972, water withdrawal was greatly reduced with the closure of all wells in the Cerchia dei Navigli, the water ring enclosing (in the past) the medieval historic center of Milan. At the end of the 1980s, the level of the aquifer increased by about 10 m. Since the 1990s, the groundwater level has shown a continuous increase mainly caused by the reduction of withdrawals for industrial uses [38]. In [39], an additional rise of 4–5 m in the central zone of Milan is described, whereas other levels are reached in the northern (8–10 m) and southern areas (2 m) of the city, causing problems at existing underground structures and infrastructures, such as the metro lines.

Among the different actions to mitigate risk and understand the effect of subsidence, the deterioration of the static conditions required the installation of a monitoring system to quantify differential settlements. In 1961, leveling points were installed on the columns of the transept, including the four columns of the tiburium. Subsequently, other monitoring points were installed on several other columns inside the church. In 1969, the system assumed a configuration with 59 monitored columns, which is still available today. Monitoring measurements are still acquired twice a year, producing a technical report with differential displacements of columns as well as monitoring information acquired with several other methods and techniques [40].

The static monitoring of relative movements between columns has both spatial and temporal components. The monitoring period covers the years 1970–2021, notwithstanding additional data available for a subset of columns even before this period. The long observation period is another characteristic of the Duomo dataset: more than 50 years of observations. According to [41], this monitoring period could be seen as “a window over historical time allowing possible insight on active processes long-term damaging processes”. Such a period is very long compared to modern dynamic monitoring methods able to measure a huge number of digital data in a relatively short time [42]. On the contrary, this period appears as a fairly short window when compared with the beginning of the construction (in 1386). More in general, monitoring plays an active role in the conservation of historical buildings [43], and the integration of different monitoring systems (such as static and dynamic methods [44,45,46,47]) is essential to prevent damage and plan conservation activities, as reported in numerous works in the technical literature (see, for example, [48,49,50,51,52,53]).

In the case of the Duomo di Milano, the monitoring system includes various sensors and data collection methods with different characteristics: optical-mechanical and digital sensors, static and dynamic measurements, in real time and with a frequency limited to seasonal and annual cycles. The monitoring system has undergone numerous changes and additions in the last 60 years following technological development. More details on the different systems installed in the Duomo are illustrated and discussed in [54,55,56,57]. A detailed description of the evolution of the monitoring systems is out of the scope of this paper, which focuses more on spatio-temporal methods for analyzing the historical archive of differential settlements. Nowadays, data stored in the Duomo’s spatio-temporal digital archive correspond to differential movements $\left\{ z(\textbf{s}_c,t_i)\right\}$ for the spatial locations $\textbf{s}_c=(x_c,y_c)$ with $c=1,...,C=59$ for the 59 monitored columns, and $t_i$ with $i=1,...,T=104$ in the considered monitoring period 1970–2021, with observations always acquired in May and November.

3 Description of the monitoring dataset of vertical movements

3.1 Initial and actual configurations of the monitoring system

The dataset analyzed in the manuscript is a collection of time series of differential vertical movements for the columns of the Cathedral of Milan. A leveling network was installed at the beginning of the ’60 to measure the displacements induced by subsidence in the city center of Milan. The initial configuration of the monitoring system (1961) was limited to some columns along the transept and a benchmark close to the facade. In 1969, the system was integrated to monitor several columns distributed in the Duomo, obtaining a configuration with 59 columns (Fig. 3). Such a scheme is still available today, and monitoring measurements are still acquired twice a year, always in May and November. The red color in the figure shows the columns with a leveling benchmark. The other columns (in gray) complete the plan of the Cathedral. As can be seen, the monitoring configuration is relatively homogeneous except for the Southern part of the nave. The dashed lines represent the connections of leveling measurements. Procedures for acquisition and adjustment are discussed in the next section.

Fig. 3

Map of the columns in the Duomo di Milano and the scheme of the measuring network

Recently (2020–2021), the archive of measurements was reconstructed and digitized to provide a continuous monitoring dataset in the period between 1970 and 2021. Some of the methods presented in this paper were also used to integrate and validate the regenerated time-series dataset. As mentioned in the previous section, the digital archive has a semestral frequency, notwithstanding initial measurements were collected with a high frequency (usually, four monitoring campaigns per year, in different seasons).

Benchmark heights $H_{t,c}$ are calculated using geometric leveling measurements (i.e., height differences between consecutive benchmarks), which are adjusted via least squares, with $t=1$ (May 1970) being the first observation, and $t=T$ (November 2021) the last. In this paper, indexes t and c refer to time and column number, respectively.

The creation of the digital archive required prediction techniques to recover some missing values, which were limited only to about 3% of the entire dataset and mainly refer to two specific epochs, May 1973 and May 2015. The first method used to recover missing value $H_{t,c}$ at a specific epoch is a linear interpolation based on the previous and next values $H_{t,c}=\frac{1}{2}(H_{t-2,c}+H_{t+2,c})$ always acquired in the same season. Results were also cross-checked with geostatistical prediction methods illustrated in this manuscript.

In the case of three columns (43, 52, and 85), it was necessary to reconnect the time-series (in November 1971 for column 43, and November 1984 for both 52 and 85) that had discontinuities. In general, the Duomo di Milano is an active construction site, and Veneranda Fabrica is in charge of preservation with activities from ordinary maintenance to major restoration interventions. In some cases, temporary scaffolding prevented the acquisition of measurements for short periods, requiring interpolation to determine missing values.

Column 39 was selected as a reference at the beginning of monitoring, i.e., it was chosen to calculate vertical movements assuming $H_{t,39}=$ const. The convention is used for historical consistency. The computed elevations are relative values, and the calculated displacements are also relative values (i.e., differential movements). Starting from the adjusted height values at different epochs, we can represent the detected movements using the following notation:

• total height variations $\Delta H_\mathrm{{TOT}}=H_{t,c}-H_{1,c}$;

• annual height variations $\Delta H_{A}=H_{t+2,c}-H_{t,c}$;

• semestral height variations $\Delta H_{S}=H_{t+1,c}-H_{t,c}$.

The collection of time-series considered in this paper is denoted as $z_{t,c}$ and corresponds to total height variations after setting the first value of all time-series to zero $z_{1970,c}=0$. The entire dataset is made up of 59 time-series, and every time-series has 104 values. The time-series for column 39 has only zero values $z_{t,39}=0$, because it is the historical reference for calculating differential movements. The cartographic location of the columns $(x_c,y_c)$ was measured in the reference system UTM-WGS84 (zone 32, North), so that the collection of time-series is also georeferenced. Each individual column is, therefore, denoted using the notation $(x_c,y_c,z_{t,c})$, obtaining a spatio-temporal dataset.

3.2 On-site observation collection and least-squares adjustment

On-site measurements are acquired using high-precision geometric leveling, which allows measuring the differences in elevation $\Delta H_{t,ij}=H_{t,j}-H_{t,i}$ between benchmarks installed on columns i and j. Special (mini) 5-mm leveling rods are connected to the metal benchmarks using a magnetic head. The historical instrument is still a high-precision optical level Zeiss Ni1 with a 5-mm range plane parallel plate micrometer (precision of $\pm 0.2$ mm/km—double run) making it more precise than most modern digital levels today available on the commercial market. The combination of level and rods forms the historical equipment for the monitoring system of the subsidence in the Duomo, and it has been used for several decades. The continuous maintenance of the optical level, benchmarks, and rods allowed the preservation of the monitoring system for more than 50 years. Nowadays, the entire system is still operational and one of the fundamental tasks of the team involved in the monitoring phase is to maintain the system active for future campaigns (over a span of several decades). It is the author’s opinion that the use of optical (mechanical) tools has a fundamental role in guaranteeing the posterity of the system. Probably, more contemporary digital instruments do not offer the same reliability over a long time span.

Nowadays, measurements are collected every 6 months. May and November are the historical months for collecting all measurements, so that the environmental conditions are similar after 1 year. The temperature is also measured to correlate measurements and environmental conditions. The scheme of the geometric leveling network is shown in Fig. 3. The network features several closed loops with multiple intersections, allowing real-time verification of metric accuracy during the measurement phase, because the overall difference in elevation in a loop must be zero. Such value (called loop closure) can be calculated as the sum of consecutive differences in elevation.

As there is no perfect measurement, the closure of a loop is compared to a specific tolerance, which depends on the number of differences in elevation M in a closed loop with M benchmarks. Tolerance is computed as $\pm 0.1\sqrt{M}$ mm. Moreover, a single difference in elevation is the average of multiple measurements, making the average more precise than single measurements. A scheme with multiple closed loops allows the staff involved in the measurement phase to evaluate metric accuracy during the measurement phase and to detect possible outliers in the observations.

After completing the on-site measurement phase (usually an entire working day), the differences in elevation are adjusted with ordinary least squares, obtaining the adjusted elevation values. Benchmark n. 39 is assumed as a reference to remove rank deficiency following the historical conventions of first measurements in the ’60s.

A system of linear equations can be written using the matrix notation $A_tH_t=\delta _t$, where $A_t$ is the design matrix with only zeros and ones, $H_t=(H_{t,1},H_{t,3},H_{t,5}...,)^\mathrm{{T}}$ the vector of unknowns, and $\delta _t$ the observation vector. Column 39 is not included in the unknown vector. The least-squares solution is calculated as $\hat{H_t}=(A_t^T A_t)^{-1}A_t^\mathrm{{T}}\delta _t$. The system of equations does not use a weight matrix, assuming that observations have equal weights for the relatively short and similar distances between the columns. Residuals $\hat{v_t}$ are computed as $\hat{v_t}=A\hat{H_t}-\delta _t$, obtaining the standard deviation for the set of observations $\sigma _{0t}=\sqrt{\hat{v_t}^\mathrm{{T}}\hat{v_t}/r}$ with r redundant measurements. The variance–covariance matrix is given by $C_{\hat{H_t}\hat{H_t}}=\hat{\sigma }_{0t}^2(A_t^\mathrm{{T}} A_t)^{-1}$ and contains the variances of elevation values along the diagonal, and the covariances (off-diagonal elements) between all pairs of elevations.

Specific in-house software was implemented in MATLAB to adjust the Duomo leveling network. The software has a similar implementation if compared to the software used in the past, which was implemented in Fortran. The new software offers more flexibility and tools while providing the same numerical results with a small numerical discrepancy($<0.0001$ mm), which is smaller than the expected precision achievable with leveling measurements. Least-squares adjustment provides elevation values with an average precision of about $\pm 0.1$ mm.

After completing a measurement campaign, a complete report showing the computed elevation values and the achieved uncertainty is delivered. Total, annual, and relative variations are also computed and illustrated with graphs and schemes.

The network geometry shown in Fig. 3 is the reference scheme for any measurement campaign. However, Duomo di Milano is a dynamic site with continuous activities. In some cases, temporary scaffolding or restoration activities prevented the installation of the rods on a few benchmarks, which were impossible to include during some monitoring periods, resulting in a few missing data. The network scheme could require local variations to include as many points as possible, notwithstanding the presented scheme remains the reference configuration. Predicting the value in an unsampled location at a specific epoch is another task considered in this manuscript.

4 Visualization of monitoring data archived using georeferenced time-series

4.1 Direct spatio-temporal visualization methods

Geospatial time-series datasets can be stored using advanced solutions, such as space-time cubes [58] or S-T databases [59]. However, S-T data are essentially numbers in one or more tables, which can be related using common fields to connect temporal and spatial information. Additional metadata are also necessary to establish the geospatial reference system, using geographic (latitude–longitude) or cartographic (East, North) coordinates. After completing a monitoring campaign and adjusting measurements via least squares, the archive is updated with new data.

In the case of multiple time-series datasets, data visualization starts with multiple time plots, i.e., graphs showing the profiles of the time-series over time. A collection of time-series can be visualized on a single graph using different line properties (colors, thickness, type of line, etc.) to separate the various profiles. However, the Duomo dataset has 59 time-series, and a single plot does not provide sufficient clarity.

Data visualization is a fundamental aspect of structural monitoring applications. Graphs, schemes, and other (usually) 2D products facilitate data interpretation, which is usually performed by expert (human) operators with the support of analytical tools and additional knowledge related to the structure. Indeed, the problem is not only the identification of anomalies but also the interpretation of possible causes and solutions. In this phase of the work, the user utilizes information from the monitoring dataset (existing data and new observations) and a variety of additional products, such as existing reports, technical drawings, photographs, 3D models, and various schemes, to name a few. Moreover, a deep understanding of the structure cannot be neglected [43], including the historical aspects, materials, construction technologies, damages, structural problems, deterioration patterns, and previous interventions, among the other main factors. This list is only indicative and not exhaustive, and multidisciplinary knowledge involving various specialists is required.

One of the aims of the work proposed in this paper is the use of both traditional and innovative processing and visualization methods for spatio-temporal datasets. As mentioned, a single graph showing the time-series is not sufficiently clear, and an alternative plot was created using trellis graphics, which is specifically designed to visualize multivariate data [60]. Figure 4 is a trellis plot, i.e., a rectangular array of plots with different time-series. Although the figure does not provide any information about the spatial location of the columns (the user must use the profiles together with the scheme of the columns reported in Fig. 3), a visual inspection of the graphs indicates the presence of differential positive and negative movements.

The time-series also shows variable patterns, including trend, stationarity, volatility, and fluctuation. In several cases, the slope is greater at the beginning of the monitoring period (1970–1980) than in more recent years (2010–2020), confirming the subsidence that affected the Duomo (and the city center of Milan) starting from the second half of the twentieth century. The closure of the wells in the city center (in the Cerchia dei Navigli) that occurred in the ’70s was beneficial, resulting in a progressive reduction of the slope of the curves.

Fig. 4

Trellis plot with different charts for all the columns of Duomo di Milano from 1970 to 2021

An alternative visualization that combines spatial and temporal dimensions can be obtained with a lattice representation of the columns (as points or other filled symbols) and a color palette associated with the relative displacement. Such layers can be plotted in the R environment using packages that can handle geospatial datasets. The input dataset is a point shapefile with East and North coordinates for the columns. The relative displacements measured at different epochs are stored in various fields.

A subset of the field values was performed to limit visualization every 5 years, starting from 1971 and completing the entire observation period. The final trellis representation is shown in Fig. 5. The author believes that such a visualization method is still insufficient to describe the subsidence’s spatio-temporal evolution in 50 years. However, a progressive change from the initial configuration (1970) can be visually detected, especially along the nave of the Duomo, with maximum variations in the North-West corner of the facade.

Fig. 5

Trellis plot combined with a spatial representation of all the columns of Duomo di Milano

It is challenging to create a single picture representing such a complex spatio-temporal phenomenon with sufficient clarity. The visualization methods proposed in this section can be rapidly generated with good knowledge of the R language. They can provide some valuable indications, but they must be coupled with other information, including alternative visualization methods using the same data and a variety of knowledge that cannot be directly combined with georeferenced time-series.

As mentioned in the introduction, this paper aims to find alternative processing and visualization methods, facilitating interpretation and considering the main events and major restorations in the latest 50 years. The aim is not to replace the methods currently used to create the semestral report but to complement the work with other solutions that can provide more information. In the following sections, particular attention will be paid to this aspect, introducing additional visualization methods that require further processing.

4.2 Trend analysis in the case of differential displacements with a conventionally fixed benchmark

In the case of differential movements measured over time, particular attention must be paid to interpreting the trend visible when the time-series are plotted. Trend analysis allows the visualization of monotonic upward or downward differential movements throughout the observation period. Most time-series used in this manuscript exhibit nonlinear trends, i.e., profiles have a non-constant positive or negative slope.

Instead of representing trends in traditional 2D charts (like in Fig. 4) using smoothing techniques (such as moving average or exponential smoothing) or decomposition methods (e.g., classical decomposition, STL, X11, etc.), the proposed work aims at establishing relative trends between the columns as an indicator of differential movements inside the Cathedral. Therefore, the goal is a spatial representation of the trend, from which additional considerations can be derived using the relative positions of the columns.

Figure 6 (top-left corner) shows the result of a trend analysis conducted using the Mann–Kendall test, which assesses whether a time-series has a monotonic upward or downward trend. The test uses the difference in signs between each observation and all subsequent observations. The figure shows four columns without a trend (39, 43, 75, 88—green color). Blue columns have a negative trend, and rose columns have a positive trend. Such results must always be interpreted as relative trends, in which column 39 has no variation over time. Therefore, the figure shows the trend of the different columns in the hypothesis that 39 remains perfectly stable over time: $z_{t,39}=0$. Alternatively, the figure can be interpreted as an indicator for columns with positive, negative, or similar trends compared to 39. Figure 6 also indicates the spatial distribution of relative trends in the Duomo. Nave and choir show mainly positive trends, especially toward the facade and apse, respectively, suggesting large subsidence values in the transept area.

Fig. 6

Trend computed considering the historical conventions $z_{t,39}=0$ (top-left), and the case with reference benchmark moved to point 18 (top-right), point 74 (bottom-left), and point 30 (bottom-right)

A fixed benchmark for adjusting the leveling network affects trend analysis. If we change the historical convention from point 39 to another column, trend analysis will only show the relative trend of the remaining column compared to a different reference point. Changing the historical convention can be carried out by subtracting the values of the new column chosen as a reference from all the remaining time-series. For example, assuming column 18 as the new reference, the new set of time-series can be calculated as $z_{t,c}^*=z_{t,c}-z_{t,18}$. The trend analysis results for the new set of time-series are illustrated in Fig. 6 (top-right corner) and show another relative trend with a new reference. If the reference is placed on column 74, trend analysis indicates that all columns exhibit a positive trend. Indeed, column 74 reaches maximum negative displacements as shown in Fig. 6, bottom-left corner.

The last case in Fig. 6 (bottom-right) represents results fixing column 30, which is the point reaching the largest positive displacements compared to column 39. One may expect the opposite result compared to the previous analysis (fixing point 74), i.e., overall negative trends. However, the figure shows some positive trends, mainly toward the facade, and negative trends in other areas. Such a result can be explained by considering that the Mann–Kendall test considers sign variations, not the total displacement’s magnitude. For this reason, some columns with smaller displacements could show a positive trend and vice versa.

To summarize the previous examples, the calculation of relative vertical displacement is the only possible solution in monitoring projects without a stable reference system. As no stable benchmark can be installed in the Duomo di Milano or the surrounding area, trend analysis outcomes are affected by the convention established at the beginning of monitoring activities. In general, several monitoring applications are based on initial conventions to define (and materialize) a reference system choosing specific points. If the aim is the generation of results not affected by conventions in the choice of the reference system, other methods more robust against the use of relative displacements are required. Trend analysis becomes more relevant for testing a specific column against the others instead of providing an integrated indicator for all columns. Section 5 describes some alternative processing methods invariant against initial conventions in the choice of the reference system.

5 Spatio-temporal clustering with relative displacement observations

The time-series dataset exhibits different behavior regarding short- and long-term trends, autocorrelation, seasonality, stationarity, and volatility. This section aims at (1) clustering relative movements, (2) interpreting the result of cluster analysis together with the spatial distribution of the measured points, and (3) exploiting clustering solutions robust against the use of relative measurements. Similarity measures quantify the distance between different time-series and create clusters with variable numbers of columns.

Spatio-temporal clustering considers the spatial distribution of the clusters, introducing both spatial and temporal constraints. Spatial relations between the columns are encapsulated into a spatial weight matrix W, a $59\times 59$ squared matrix reflecting spatial relationships between columns [61].

The spatial weight matrix can be generated using different methods. The input is a vector point layer representing the spatial location columns and relative displacement stored in multiple fields as attributes. W was constructed with two different methods:

• A: Delaunay triangulation, which tends to maximize the angles of triangles;

• B: k-Nearest Neighbor (k-NN) using $k=8$ points.

The results with both methods are shown in Fig. 7. As can be seen, the k-NN method provides more connections than the Delaunay triangulation, because it forces each column to reach a minimum number of connections with the others.

Fig. 7

Visualization of spatial weight matrices created with Delaunay triangulation (A) and 8-Nearest Neighbor (B)

After constructing the weight matrices, spatio-temporal clustering was carried out with the approach proposed by [62], which is known under the name Skater: Spatial “k” luster analysis by tree edge removal. Skater includes the spatial weight matrix in the clustering process as additional spatial constraints. The temporal dissimilarity is evaluated using the Euclidean distance, using the values from November 1970 to November 2021 (values in May 1970 are all zeros and are excluded from processing). A predefined number of four classes was set as input, making the comparison of results more direct.

Clustering results are summarized in Fig. 8. The colors used in the representation are indicative, and there is no direct correspondence between the figures. In other words, the reader must look only at the overall distribution achieved without connecting the use of a particular color in the figures. As can be seen, the different methods tend to identify specific parts of the Duomo, notwithstanding differences in the different clusters. These parts can be identified as:

• the nave and aisles (the Western part);

• the transept, which is split into the southern/central part (the area with maximum negative displacements) and the northern transept;

• the apse and the choir (East), and the central part of the transept.

Fig. 8

Comparative visualization of trend analysis results with different methods. The colors used in the different figures do not correspond to the same clusters (colour figure online)

The general advice in the case of cluster analysis with relative displacements consists of using different methods and interpreting results with a visual representation showing the different clusters on 2D maps. As explained in the previous section, processing methods could be very sensitive to the conventional choices to set a relative reference system and possible distortions in time. Different factors could highly affect results, leading to wrong interpretations if just a single clustering algorithm is chosen. The intention is not to finely tune the parameters used as input in the analysis, which is a complicated task in cluster analysis.

The results show that the subdivision previously obtained is rather consistent. The final clusters aggregate the four columns of the tiburium (n. 74, 75, 84, 85) with the southern parts of the transept. A visual inspection of the profiles of these columns shows a relatively constant negative trend, which assumes minimum values for column 74. This area is historically the zone with the most significant depression due to subsidence. A second area can be identified around the previous one, surrounded by a third group of columns, including those of the northern transept and the apse. Finally, a final cluster includes a large part of the nave.

6 Time-series forecasting and its application to monitoring

The monitoring dataset can be considered a set of time-series with a common temporal domain between May 1970 ($t=1$) and November 2021 ($t=T)$. Time-series forecasting allows the estimation of future values at times $T+1, T+2,..., T+F$, where F is the forecast horizon. In the case of the Cathedral of Milan, time-series forecasting is used with a specific goal: the identification of anomalous values during a new monitoring campaign. An automatic comparison can be made by evaluating the difference between the forecasted value $\hat{z}_{T+1,c}$ (one step ahead) and the measured value $z_{T+1,c}$ after acquiring and adjusting new on-site observations.

Usually, such comparison is carried out using the last observed value, applying the so-called naive forecasting method $\hat{z}_{T+1,c}=\hat{z}_{T,c}$. However, the method does not include drift and gives a horizontal straight line for multiple forecasted values. There exist different advanced solutions for univariate time-series forecasting and some methods were implemented and tested using a single script. The aim was to understand which method (or which combination of methods) performs better in forecasting future values. In the case of reliable results, a new strategy for the comparison of newly acquired data can be employed in future monitoring campaigns.

Section 6.1 briefly presents the implemented forecasting strategies, whereas Sect. 6.2 describes the achieved results and quantifies metric accuracy and uncertainty.

6.1 The proposed approach for univariate time-series forecasting

Different methods for univariate time-series forecasting were tested. The implementation is based on the forecast package [63], which contains different forecasting functions. An overview of the methods is illustrated and discussed using the data $\left\{ z_{t,49} \right\}$ of column $c=49$. Such time-series is addressed as $z_t$ in the rest of this section.

The graph is shown in Fig. 9. The upward trend reveals a positive movement for column 49 compared to the reference column 39. The autocorrelation function (ACF) shows very slow decay. The partial autocorrelation function (PACF) shows that only the first lag value is statistically significant.

Fig. 9

Graphs of differential displacements for column 49 (top-left) and second-order difference (bottom-left), the autocorrelation, and partial autocorrelation plots (right)

The time-series is not stationary and the first difference (the change between consecutive observations in the original series) was computed to make the time-series stationary using the relationship $\nabla z_t = z_t - z_{t-1}$. However, the differenced data are still not stationary, and it was necessary to calculate the second-order differences $\nabla ^2 z_{t}$, which are the first differences of first differences. The KPSS test (Kwiatkowski–Phillips–Schmidt–Shin [64]) indicates that the time-series of second differences $\nabla ^2 z_{t}$ is stationary.

Most time-series of the Duomo collection are not stationary and show a slowly decreasing ACF. Stationarity is a fundamental requirement for using ARMA (autoregressive moving average) models [65], whereas other forecasting methods do not require stationary data.

The next subsection discusses four methods for time-series forecasting, which were integrated into a single R script to obtain automatic forecasts for the different time-series of the Duomo dataset. Column 49 is still used as an illustrative example. Results with all columns are described in Sect. 6.2.

The following subsections give a very synthetic description of the different methods. The reader is referred to [66] for more technical details.

6.1.1 ARIMA models

An AutoRegressive Integrated Moving Average model—ARIMA(p, d, q)—has the form

$$\begin{aligned} (1-\phi _1B - \cdots - \phi _p B^p) (1-B)^d y_{t} = c + (1 + \theta _1 B + \cdots + \theta _q B^q)\varepsilon _t, \end{aligned}$$

(1)

where p is the order of the autoregressive part, d is the order of differencing, and q is the order of the moving average part. B is the lag operator defined as $B^j z_{t} = z_{t - j}$.

The auto.arima function in the forecast package automatically fits an ARIMA model to univariate time-series, selecting the values (p, d, q), and returning the best ARIMA model using the corrected Akaike’s Information Criterion (AICc).

In the case of column 49, an ARIMA(1,2,1) was automatically chosen. After computing the model, the bias-adjusted mean squared error (MLE) of the innovations variance was 0.018 mm$^2$. The same function can also handle seasonal ARIMA models, which could occur in the Duomo dataset, because monitoring campaigns are carried out in May and November. The different environmental conditions in different seasons could affect structural behavior.

A seasonal ARIMA can be written as ARIMA$(p,d,q)\times (P,Q,D)_s$. It has a similar structure to non-seasonal ARIMA while introducing seasonal autoregressive, differencing, and moving average components.

6.1.2 Neural network autoregression

Neural network autoregression for time-series forecasting is an extension of feed-forward networks in which p lagged values $z_{t-1},z_{t-2},...,z_{t-p}$ are used as input for forecasting value $z_p$. The nnetar function in the forecast package uses a single hidden layer with k neurons. It can also handle univariate seasonal time-series using the last values from the same season as input parameters. The notation NNAR$(p,P,k)_s$ indicates a network with k neurons in the hidden layer and values $(y_{t-1},y_{t-2},\dots ,y_{t-p},y_{t-s},y_{t-2s},\dots , y_{t-Ps})$ as inputs.

6.1.3 ETS

ETS models form a family of forecasting methods based on a measurement equation that describes the observed data and some state equations that describe how the unobserved components or states (level, trend, and seasonal) change over time [66]. The book [67] provides a complete description of ETS models, also called space-state models. The ets function in the forecast package determines model parameters. After fitting a model, forecasts can be generated. The function can also automatically select model parameters with the AICc or other information criteria.

6.1.4 Natural cubic smoothing splines

The use of natural cubic smoothing splines for time-series forecasting is described in [68]. The method is an extension of cubic splines, often used to interpolate noisy data. Additional constraints are imposed with this method, so the spline can provide better forecasts without compromising the fit [66]. The function splinef of the forecast package was included in the script to fit a natural cubic smoothing spline and calculate forecasted values.

6.2 Accuracy evaluation

The proposed time-series forecasting methods were implemented in a single script to process the entire Duomo dataset simultaneously. In other words, each column was processed with the four forecasting methods and the results were cross-checked. Accuracy evaluation is essential to understand if the forecasted values are representative of future relative displacements. The analysis of uncertainties after fitting a model is not always a reliable solution, because good fitting statistics do not necessarily mean that the forecasts will also be accurate.

Accuracy evaluation was performed by splitting each time-series into training and test data. We indicate the observation period as $(t_1, t_2,..., T)$, in which $t_1$ corresponds to May 1970 and T to November 2021. The model is estimated using only training data $(z_1,z_2,...,z_p)$, then forecasts $(\hat{z}_{p+1},\hat{z}_{p+2},...,\hat{z}_T)$ are generated and compared with the test data $(z_{p+1},z_{p+2},...,z_T)$. The root-mean-square error $RMSE_c$ for a generic column c can be estimated as

$$\begin{aligned} \mathrm{{RMSE}}_c=\sqrt{\frac{1}{T-p}\sum _{t=p+1}^{T} \left( z_{t,c}-\hat{z}_{t,c} \right) ^2}. \end{aligned}$$

(2)

Accuracy evaluation was carried out with multiple training and test datasets by moving index p in time, starting from 2011 (i.e., training dataset 1970–2011, and test dataset 2012–2021), covering 10 years. The procedure is then repeated moving p one epoch ahead (training dataset 1970–2012, test dataset 2013–2021, 9-year period) in an iterative way. The last iteration covers a single year (training dataset 1970–2020, test dataset 2021).

The results are shown in Table 1, which illustrates only the data in even years. Another clarification is required to understand the meaning of the reported values. As the analysis was carried out for all columns (except 39), 58 RMSE$_c$ values can be calculated for each testing period. The mean RMSE value and standard deviation were estimated to synthesize evaluation results.

A graphic visualization for the different columns (in different testing periods) is shown in Fig. 10. Table 1 also shows two rows called Av.1 and Av.2. Such values correspond to the average values estimated with the four methods (Av.1) and the first three methods (Av.2).

Fig. 10

RMSE values of single columns in different testing periods. Results are in mm

Averaging forecasts calculated from different methods could provide overall better results [69]. In the case of Av.2, results with the cubic smoothing spline method were removed. In fact, the analysis of separated values suggested that the smoothing spline method provides the worst forecasting results.

Slightly better results were achieved with the ETS method, followed by Av.1. All the methods provide similar accuracy for short testing periods. However, the computed accuracy is worse than the expected precision of the geometric leveling network, which is about $\pm 0.1$ mm. We must also consider that such precision is determined via least-squares adjustment and refers to elevation value standard deviation. The precision of differential variation can be assumed as $\pm 0.1\sqrt{2}=\pm 0.14$ mm.

If we suppose to use only the results of the ETS method and a short-term forecast (1 year, i.e., two forecasted values), the expected accuracy is 0.18 mm ± 0.09 mm, which is slightly larger than the precision of variations. A 3-year forecast shows an accuracy of 0.24 mm ± 0.11 mm, and the value remains relatively constant after increasing the testing period. The worst forecasting results were achieved mainly for those columns in the Northwestern corner of the Duomo.

As mentioned at the beginning of the section, the main goal of time-series forecasting is the comparison of forecasts with actual relative displacement detected during a new monitoring campaign instead of the direct use of naïve forecasts. The idea is to develop a strategy for anomaly detection based on two-step ahead forecasts, which means an additional year (May–May, November–November). In this case, accuracy evaluation results and the expected metric accuracy after adjusting leveling measurements and computing relative variations can be compared to find discontinuities exceeding a $\pm 0.2$ mm threshold.

Such a comparison could also integrate the traditional semestral monitoring report transmitted to Veneranda Fabrica del Duomo di Milano. It provides a quick numerical check to identify possible anomalies after adjusting the observations and calculating relative differences.

Table 1 Accuracy evaluation using variable validation datasets, showing mean value and standard deviation of RMSE values for all the columns (except 39)

7 Spatio-temporal prediction and area-based visualization techniques

This section aims to explore alternative visualization methods able to describe with more clarity the effect of subsidence in the Duomo di Milano during half a century of measuring activities. The different types of analysis described in the previous sections used localized points related to columns’ relative displacements. However, subsidence affects the whole area. Extending data processing and visualization to the entire surface provides alternative ways to understand differential movements, which are here considered a spatio-temporal phenomenon.

Spatio-temporal processes in R can be managed with the spacetime and the gstat packages. The second one uses space-time classes to estimate the spatio-temporal variogram and perform spatio-temporal prediction [70]. Additional packages used in this section are the raster, terra, and rasterVis packages, which provide enhanced visualization tools for raster datasets.

Data processing was carried out in two ways: (1) pure (geo)spatial analysis, in which the different measurement periods are considered independent datasets (7.1), and (2) complete space-time analysis, in which temporal correlation is also exploited (7.3).

7.1 Spatial interpolation of relative subsidence measurements

Interpolation allows prediction of relative displacements $z_{t,s}$ in (unmeasured) locations $\textbf{s}=(x,y)$ starting from $z_{t,c}$ values for the 59 columns at spatial locations $\textbf{s}_c=(x_c,y_c)$. The set of unmeasured locations corresponds to the whole area enclosed by the external walls of the Cathedral, which is discretized with a 0.5 m $\times$ 0.5 m grid. The considered observation period is still 1970–2021, i.e., time does not consider additional years here, whereas space has a new bounded domain given by the perimeter wall. The extension of the temporal domain (2022, 2023,...) will be discussed in Sect. 7.3.

Processing is carried out by generating separated prediction maps at different epochs using traditional spatial interpolation [28]. Spatial interpolation assumes that the measured relative displacement at a given time $t=t^*$ is sufficient to predict relative displacement in other nearby locations at the same time $t^*$. The sample semi-variogram $\hat{\gamma }(d)$ describes the relationship between data and distance

$$\begin{aligned} \hat{\gamma }(d) = \frac{1}{2 N(d)} \sum _{N(d)}^{} (z(\textbf{s}_i) - z(\textbf{s}_j))^2, \end{aligned}$$

(3)

in which the variance between two points $\textbf{s}_i$ and $\textbf{s}_i$ depends on the distance d

$$\begin{aligned} d=\sqrt{(x_i-x_j)^2+(y_i-y_j)^2}, \end{aligned}$$

(4)

and N(d) is the set of all pairs of locations separated by d.

Prediction is carried out with ordinary kriging, which is a widely used method for data interpolation in geostatistics [71, 72]. Ordinary kriging predicts values at not sampled point locations (i.e., $\hat{z}(\textbf{s}_0)$) based on the observed data and the estimated semi-variogram, using a linear combination of observed values $z(\textbf{s}_i)$ and weights $\lambda _i$

$$\begin{aligned} \hat{z}(\textbf{s}_0)=\sum _{i=1}^{n} \lambda _i z(\textbf{s}_i). \end{aligned}$$

(5)

Automation in prediction is achieved with the automap package, which automatically fits the semi-variogram model to the data and generates the interpolated map. Such operation was repeated with an automated script for the observation period $t_1,t_2,\dots ,t_T$, considering the different epochs as independent events.

Results with the data from 1971 to 2021 are shown in Fig. 11, showing the interpolated displacement maps with kriging every 3 years. The lattice-based visualization was achieved with rasterVis, which can print single raster files or raster stacks. A stack is a collection of raster images with the same spatial extent and resolution organized using a single object. The raster stack was also cropped using the perimeter of the Duomo.

Fig. 11

Interpolated maps of the predicted differential displacements using spatial prediction based on kriging. The axes show UTM coordinates: East and North

The sequence of interpolated maps shows the temporal evolution of the displacements. According to the historical convention (point 39 is the fixed benchmark to calculate relative displacements), the maximum negative displacement area is the southern transept. Differential displacements between the different columns of the Cathedral are visible, and the most significant relative values are located between the northwest side of the facade and the southern transept.

Geospatial prediction allows the creation of useful representations to illustrate the effect of subsidence over time better. Obviously, an additional work is necessary to show the distribution for the whole surface of the Duomo, especially if compared to the direct visualization of values stored in the archive. However, automation is possible, because operations become repetitive for the different epochs. In this case, the implemented function performed four steps: (1) extraction of the measured displacements at a specific time $t^*$, (2) estimation of the sample semi-variogram, (3) model fitting, and (4) generation of the interpolated map.

It is also important to mention that spatial interpolation is not carried out to predict relative displacement for columns without a benchmark (see the scheme in Fig. 3), notwithstanding this is a natural consequence of data processing. The main goal is not recovering the values for columns without a measuring point, but the creation of visualizations that can clarify the effect of subsidence with enhanced visual representation. The raster stack is the main output of this phase, but it can also be used to create additional products as described in the next section.

7.2 Alternative visualization methods for relative displacement acquired over time

As the interpolated maps are stored using a raster stack, it is possible to create an animated plot, i.e., a sequential plot that can be exported as a video clip (for instance, an animated GIF). Another solution is the creation of a 3D cube with the cubeview package. The cube can be interactively rotated and translated in space, and slices through each dimension can be visualized to quantify displacement. Another way of representing digital elevation models (DEMs) is a sequence of 3D plots after scaling relative displacements. 3D plots may provide a graphic visualization of surface deformations starting from the initial flat surface in 1970.

All the previous methods require a computer for visualization and are useful for interactive exploration of the data, and they are not suitable for printing 2D representations. Indeed, the use of spatial and temporal information makes the production of 2D plots more complicated, requiring multiple maps structured in a grid, such as in the case of Fig. 11.

An interesting 2D graphic representation of spatio-temporal phenomena is the Hovmöller diagram, which is very common for interpreting meteorological data. The Hovmöller diagram is generated from a raster stack using the average values in single columns (North direction) or rows (East direction), storing the values on a new x axis. The y axis instead represents time.

The Hovmöller diagrams for the Duomo di Milano are shown in Fig. 12, which depicts different behaviors in the two directions. The Cathedral is oriented to the East (meaning that the altar is on the East side of the church), longitudinal (along the nave), and transversal (along the transept) axes are aligned to cartographic axes. The first figure shows the increasing differential movement between the nave and transept. Further movement is also visible between the transept and the apse and has a smaller magnitude.

The second Hovmöller diagram instead cuts the Duomo in the longitudinal direction. The red area at the bottom indicates the transept’s South area, where maximum negative displacements are visible. The maximum is reached in the top-right corner, indicating the progressive movement over time and the presence of differential movements. A comparison of the two diagrams shows an overall lack of symmetric behavior, with local positive or negative concentrations. The different magnitude of the legend must also be considered during the visual interpretation of the two diagrams.

Fig. 12

Hovmöller diagrams for the east–west (façade-apse) and north–south (along the transept) axes of the Duomo; coordinates in UTM, east, and north, respectively

Additional visualizations can be created using raster algebra applied to the raster stack with relative displacements. A useful representation is a raster plot of the velocity field and its evolution over time. The average velocity in different periods (1970–75; 1981–85; 1991–95; 2015–21) is shown in Fig. 13. The same color convention was applied to all images, and velocity is expressed in mm/year. As can be seen, the velocity at the beginning of the monitoring period was much larger than the values in recent years. Velocity is also a relative value due to the relative displacements calculated using column 39 as a stable reference.

The relative velocity between 1970 and 1975 reaches values larger than 1.4 mm/year between the columns of the transept and the southwest corner of the Duomo. Relative velocity is much more uniform after a decade, as demonstrated by the figure in the period 1981–1975. The magnitude is also significantly reduced. Such results can also be related to the decision to close all the wells in the Cerchia dei Naviglia, starting from 1972 [33], recovering several meters in terms of aquifer level and reducing differential displacements.

Nowadays, the average velocity is relatively small and more uniform, notwithstanding differential movements are still present. The average is generally lower than 0.1 mm/year for the entire area inside the Duomo, except around column 28. Historically, the column with the most significant positive displacement has always been column 30 (the one in the northwestern corner of the cathedral). However, the (positive) relative variations of column 28 increased in 2015.

Fig. 13

Average velocity in different periods computed from differential movements

7.3 From spatial to spatio-temporal prediction

The previous Sect. 7.1 considered the measured relative displacements $z_{t,c}$ at different observation periods t as independent observations, creating a set of interpolated maps without exploiting temporal aggregation.

Spatio-temporal interpolation uses the spatio-temporal semi-variogram [73, 74]

$$\begin{aligned} \hat{\gamma }(d,u) = \frac{1}{2 N(d,u)} \sum _{N(d,u)}^{} (z(\textbf{s}_i,t_i) - z(\textbf{s}_j,t_j))^2, \end{aligned}$$

(6)

where N(d, u) is the set of all pairs of locations separated by d and $u=t_i-t_j$. The estimation of $\hat{\gamma }(d,u)$ requires preliminary conversion of input data into an STFDF class (a class for spatio-temporal data based on a full space-time grid [75]), which can handle spatio-temporal data when the spatial location of points remains constant over time. Then, variogram modeling is carried out with the procedure described in [76]. Finally, spatio-temporal kriging can be run on a grid with space and time components.

Instead of running spatio-temporal prediction with the whole dataset, two particular cases were considered in this work: (1) an alternative spatio-temporal procedure to recover missing values, and (2) the calculation of future values, i.e., an extension of the temporal domain.

7.3.1 Case n.1: recovering missing values with spatio-temporal prediction

The previous sections described the creation of the Duomo dataset and some of the operations used to recover missing data, which are mainly based on univariate time-series analysis. In other words, temporal interpolation for recovering missing values did not consider spatial correlations.

Spatio-temporal interpolation is a relatively flexible approach in this sense. Recovering missing values can be carried out with spatio-temporal kriging after creating a suitable dataset with data before and after the missing period.

As described in Sect. 3, all the values for the monitoring campaign in May 1973 were recovered using the average of measurements in May 1972 and May 1974, which is a basic linear interpolation suitable when velocity is constant. We decided to recalculate such values using spatio-temporal prediction, considering the measurements acquired in May from 1970 to 1976, assuming the values in 1973 as additional unknowns. Such an approach considers more epochs and can deal with variations in the velocity field.

After creating the STFDF class with a script, the empirical variogram was generated. Variogram modeling was a more complicated task requiring several trials with the different models available and choosing the best one with a visual inspection and analysis of mean square errors. Results achieved using the metric variogram model are discussed here.

The variogram has the form

$$\begin{aligned} \gamma _m(d,u)=\gamma _{joint} \left( \sqrt{d^2+(ku)^2} \right) , \end{aligned}$$

(7)

where $\gamma _{joint}$ is any known variogram that may include a nugget effect, and k is an anisotropy parameter between spatial and temporal units.

After fitting the theoretical model, spatio-temporal kriging prediction [77] can be carried as

$$\begin{aligned} \hat{z}(\textbf{s}_0,t_0)=\sum _{i=1}^{n} \lambda _i z(\textbf{s}_i,t_i). \end{aligned}$$

(8)

The prediction was initially carried on a spatial grid limited to the locations of the columns, whereas the temporal domain included May 1973 (Fig. 14). The implementation provides an exact interpolation, meaning that each interpolated value is estimated minimizing the prediction error for that point, and the predicted values in sampled locations are the observed values.

Fig. 14

Results of spatio-temporal prediction evaluated on the different column locations. The red rectangle represents the recovered missing value in 1973 (colour figure online)

Figure 15 shows a comparison between the relative variations computed with independent interpolation of the time-series (the values used in the historical archive), and the alternative solution from local spatio-temporal interpolation. As can be seen, differences are relatively small, with an average value of − 0.05 mm and a standard deviation of ± 0.13 mm. The difference for column 39 is 0.12 mm instead of zero (historical convention), and the most significant differences correspond to columns 51 and 67, with a discrepancy of about 0.4 mm.

Fig. 15

Difference between values in May 1973 recovered by univariate time-series analysis and spatio-temporal prediction

The method can also be used on an extended spatial $0.5m\times 0.5m$ grid covering the entire surface of the Duomo, obtaining the results in Fig. 16. Such a representation is comparable with the pure spatial interpolation shown in the previous section, notwithstanding a different color scale was used to distinguish pure spatial from spatio-temporal prediction.

Fig. 16

Spatio-temporal interpolation extended to a dense grid covering the surface of the Duomo

The opportunity to recover missing values can be considered an alternative to temporal interpolation based on univariate time series analysis. In addition, the method provides the opportunity to validate metric accuracy using a set of measured relative displacements for some input locations excluded from processing.

7.3.2 Case n.2: spatio-temporal forecasting

The temporal domain of the dataset (1970–2021) can be extended to estimate future values with spatio-temporal kriging. We consider the relative measurements taken in May from 2012 to 2021 (10 years) to predict values for May 2022. The values can be organized into another STFDF class, from which another spatio-temporal empirical variogram is generated.

Variogram modeling was again carried out with the metric model, and prediction with kriging was extended to 2022, forecasting future variations. Results are shown in Fig. 17 together with the values in the previous years. The method can also be extended to predict values for the whole surface, notwithstanding such an aspect was not considered in this manuscript.

Finally, results predicted with spatio-temporal kriging can be compared to the forecasts using ETS. Figure 18 shows a comparison of the two methods using the difference in the computed values. The average value of such differences is $-0.09$ mm, and the standard deviation is $\pm 0.12$ mm. It is interesting how the largest differences are located around columns 26, 28, 30, 90, 92, and 94, located in the area with maximum (positive) relative displacements.

Fig. 17

The predicted values in May 2022 using spatio-temporal kriging (red rectangle)

Fig. 18

Difference between values in May 2022 forecasted by univariate time-series analysis and spatio-temporal prediction

8 Conclusions

The paper presented different methods for spatial (S), temporal (T), and spatio-temporal (S-T) analysis applied to georeferenced time-series for structural monitoring. Examples were illustrated and discussed using the Duomo di Milano monitoring archive, which features more than 50 years of differential movements induced by subsidence. However, the proposed methods can be used in several practical applications not limited to historic structures.

Spatio-temporal methods could be extended to other structural monitoring applications where S-T datasets are collected at regular or irregular sampling intervals. Examples could be satellite-based measurements using GNSS, SAR, or georeferenced monitoring information captured with other sensors and tools, resulting in time series for specific points for static and quasi-static monitoring. These methods are normally used in large areas beyond the scale of the building (for example, at the territorial level), but the technical literature reports numerous applications for structures, infrastructures, and at the environmental level [78,79,80]. Other methods able to capture S-T information at the scale of the building could be geodetic solutions, such as total stations and leveling networks. Automation is also feasible with robotic total stations [81] and hydrostatic leveling systems [82].

After creating the archive, one of the aims of the proposed work was to develop solutions for validating, extracting, and visualizing S-T information. Problems, such as clustering, prediction, forecasting, recovery of missing values, and anomaly detection, were illustrated and discussed. Data uncertainty and accuracy were also evaluated, considering the expected precision of adjusted monitoring measurements and the results from S, T, or S-T processing.

Although the paper considered the vertical movements of the columns of the Duomo di Milano, other S-T monitoring datasets are available for the cathedral. For instance, starting from the ’60, leveling measurements are also acquired in the square and on the main buildings (Galleria Vittorio Emanuele II, the Royal Palace, the Church of San Gottardo, etc.) around the Cathedral, forming a large monitoring network covering an area of about 160,000 m$^2$. The external network is still measured once a year, always in May. Another example is the main spire, which is monitored with automatic and manual methods. Differential vertical movements are measured (since 2008) around the octagonal tiburium when the scaffolding was built around the main spire for restoration activities.

Overall, spatio-temporal geostatistics can help solve specific tasks, such as the recovery of missing values using prediction approaches able to consider both spatial and temporal dimensions. Pure temporal (T) or spatial (S) strategies are also helpful for several tasks proposed in this work. Missing values were initially recovered with pure temporal approaches without quantifying uncertainty, and cross-comparison with the alternative approaches presented in the manuscript was considered a solution for data validation. It is the author’s opinion to combine different methods instead of concentrating only on a specific algorithmic solution. Defining optimal input parameters is also a challenging task. S-T analysis is also more computationally expensive, especially space-time kriging with large datasets.

The results achieved with various methods were not the same, offering the opportunity to interpret the phenomenon under investigation from different perspectives, considering the pros and cons of the applied methods and their combined use. The interpretative work can also benefit from additional solutions to graphically represent spatio-temporal changes, which must be related to historical events and deep knowledge of the monument in any monitoring application with heritage buildings.

The paper also considered the use of clustering methods robust against the conventional choices in defining a reference system, which is common in monitoring projects and often provides relative (differential) displacements. Time-series cluster analysis (including the spatial component) can be carried out with dissimilarity measures invariant to the conventional choices used to remove rank deficiency in least-squares adjustment, enabling the identification of time-series groups independently of the solutions implemented to calculate relative displacements.

Finally, extending S-T processing to specific applications at the level of individual structural elements is possible when measurements are captured with techniques capturing surfaces or dense sets of points. Examples could be the photogrammetric monitoring [83] of architectural elements (such as columns, arches, and vaults) starting from images acquired with synchronized cameras. Digital image correlation (DIC) [84], image-based condition monitoring [85], or vision metrology [86] are other applications in which combined spatio-temporal geostatistics can provide valuable information.