Pre-Collapse Space Geodetic Observations of Critical Infrastructure: The Morandi Bridge, Genoa, Italy

We present a methodology for the assessment of possible pre-failure bridge deformations, based on Synthetic Aperture Radar (SAR) observations. We apply this methodology to obtain a detailed 15-year survey of the Morandi bridge (Polcevera Viaduct) in the form of relative displacements across the structure prior to its collapse on August 14th 2018. We generated a displacement map for the structure from space-based SAR measurements acquired by the Italian constellation COSMO-SkyMed and the European constellation Sentinel-1A/B over the period 2009–2018. Historical satellite datasets include Envisat data spanning 2003–2011. The map reveals that the bridge was undergoing an increased magnitude of deformations over time prior to its collapse. This technique shows that the deck next to the collapsed pier was characterized since 2015 by increasing relative displacements. The COSMO-SkyMed dataset reveals the increased deformation magnitude over time of several points located near the strands of this deck between 12th March 2017 and August 2018.


Introduction
Over the past 20 years Interferometric Synthetic Aperture Radar (InSAR) and Persistent Scatterer Interferometry (PSI) have established themselves as a methodology capable of monitoring surface deformations in urban environments, thanks to the presence and density of coherent radar targets. The imaging capability of the current generation of SAR constellations (e.g., COSMO-SkyMed/SAOCOM, TerraSAR-X/PAZ and Senintel-1A/B) has been significantly improved, to the extent that SAR images are now provided with meter to sub-meter spatial resolution, and are acquired with a revisit time of a few days, in most cases according to a regular acquisition schedule which over the years can result in hundreds of images per individual site [1].
There is a wealth of scientific literature demonstrating the benefits of InSAR to monitor a wide spectrum of urban processes, such as subsidence [2][3][4][5], foundation settlement [1, 6,7], deformation of new buildings [8], surface effects due to tunneling [9][10][11], ground motions associated to faults and local tectonics [12][13][14]. The implementation on real case studies proves the increasing use of these techniques on a routine basis and how they can influence the management of infrastructure assets and preservation of the built environment [1, 15,16].

The Morandi Bridge
The Polcevera viaduct (also known as Morandi bridge) is a prestressed concrete cable-stayed bridge composed of 11 spans (Figure 1) for a total length of about 1300 m [20]. It connected the motorways Milano-Genova and Genova-Savona, in northern Italy, crossing the Polcevera river on an area composed of marine and alluvial sand. Starting from the southern side, the bridge is supported by two single piers (1 to 2), six V-shaped piers (3 to 8) and three A-shaped frames (9 to 11) ( Figure 1). Piers 3 to 8 consist of two inclined piers connected by a double cantilever deck of variable length [21], while piers 9 to 11 were designed as independent "balanced systems" made of a three-span continuous deck resting on four supports. Each of the three-span decks is supported externally by the anchorage of two prestressed stay-cables passing over a 90 m high mast. The anchorages are connected by transverses located at a distance of 10 m from the cantilever tips. The two internal supports are the top end of two pairs of diverging inclined legs. Each balanced system is supported by a concrete ribbed foundation raft on concrete piles. On August 14th 2018, the system number 9, including an approximately 240 m long section of the adjacent deck suddenly collapsed, causing the death of 43 people.
The bridge was built between 1960 and 1966. In 1993, a detailed inspection of the balanced system number 11 reported oxidation of the metallic membrane protecting the strands together with humidity, absence of grout, slacked or sheared strands, and signs of tensile rupture due to steel cross-section reduction [22]. The 1993 report informed the structure management authorities of the ongoing material degradation and micro-cracking phenomena occurring within the prestressed concrete [23]. The balanced system 11 was subject to retrofitting and many other maintenance activities were carried over during the years [24]. The strands replacement for the balanced system 9 was planned to start two months after the collapse [23]. From West to East, the bridge includes two single piers (1 to 2), six V-shaped piers (3 to 8) and three independent balanced systems (9 to 11). Dimensions in meters. (B) Three-dimensional view of the bridge.

Materials and Methods
The methodology adopted in this paper followed two main steps including SAR multi-temporal interferometric analysis (MT-InSAR) and Markov Chain Monte Carlo three-dimensional inversion.

Multi-Temporal InSAR
Multi-temporal InSAR combines several acquisitions highlighting the temporal behavior of the relative displacement between each acquisition. In this paper surface deformation measurements have been generated using the SARPROZ software extending the standard linear PS technique [25] [26] to solve for non-linear motion with no a priori information. SRTM has been used as a supporting DEM for the PS processing. COSMO-SkyMed (CSK) and Sentinel-1 (SNT) data have been specifically processed for this study while the Envisat time-series are part of a nationwide monitoring program [27]. Table 1 shows information related to the dataset processed and used in this study. We consider a set of N+1 coregistered single look complex (SLC) SAR images and generates a set of N single look differential interferograms (for both SNT and CSK) with respect to a master imaged chosen maximizing the following coherence formula:

Materials and Methods
The methodology adopted in this paper followed two main steps including SAR multi-temporal interferometric analysis (MT-InSAR) and Markov Chain Monte Carlo three-dimensional inversion.

Multi-Temporal InSAR
Multi-temporal InSAR combines several acquisitions highlighting the temporal behavior of the relative displacement between each acquisition. In this paper surface deformation measurements have been generated using the SARPROZ software extending the standard linear PS technique [25,26] to solve for non-linear motion with no a priori information. SRTM has been used as a supporting DEM for the PS processing. COSMO-SkyMed (CSK) and Sentinel-1 (SNT) data have been specifically processed for this study while the Envisat time-series are part of a nationwide monitoring program [27]. Table 1 shows information related to the dataset processed and used in this study. We consider a set of N+1 coregistered single look complex (SLC) SAR images and generates a set of N single look differential interferograms (for both SNT and CSK) with respect to a master imaged chosen maximizing the following coherence formula: where, B ⊥ , T, F D and B ⊥C , T C , F DC are the perpendicular, temporal and doppler baselines and their critical values, respectively. The selected master images for CSK ascending, descending and SNT descending were acquired on 1st January 2015, 11th June 2015 and 20th July 2017, respectively. After selecting a master image, we applied the proposed algorithm only on pixels showing an inverse amplitude stability (i.e., ratio between the mean calibrated image amplitudes and its standard deviation) higher than 0.65, choosing 4 reference points as explained in the results section of this manuscript. The advantage of this approach is to consider only persistent scatterers (PS) that are unaffected by geometrical decorrelation, i.e., PS in interferograms characterized by perpendicular baseline larger than the critical baseline still contain useful phase information [25,26]. Specifically, our output consisted of four different time-series. The interferometric phase difference between two neighbor pixels x r and x s in the interferogram i can be written as the sum of five terms: where ∆C topo and ∆C V are proportional to the DEM errors and mean velocity difference between pixels x r and x s , respectively. ∆w r,s contains the sum of three contributions related to non-linear motion and atmospheric phase screen (APS) noise. Assuming ∆w r,s < π for neighbor pixels it is possible to estimate ∆C topo and ∆C V by maximizing the absolute value of the temporal coherence |ξ PS | for neighbor pixels: The following steps consisted in unwrapping every interferogram and estimating the phase term ∆w i r,s . Once the data were unwrapped we estimated the low pass component using a three sample temporal base-line weighted moving average and assumed this residual phase term as an estimation of the non-linear motion contribution. These steps have been performed considering only neighbor pixels. The same processing was then applied considering only differences between the reference point and all persistent scatterers including the estimated APS in Equation (4).
The accuracy of the estimated velocities and DEM errors can be derived by a regression analysis considering the irregularly sampled data (perpendicular baseline and temporal baseline in the case of the CSK constellation. As in [28] we define the standard deviation of scatterer height (σ 2 ∆h ) and velocity (σ 2 ∆v ): where R is the target-sensor distance, θ is the incidence angle, λ is the sensor wavelength, σ 2 φ is the phase noise variance supposed independent from the acquisition and σ 2 are the variance of the Remote Sens. 2019, 11, 1403 5 of 14 perpendicular and temporal baselines respectively. For further details on the technique we refer the reader to [25,26,28].

Markov Chain Monte Carlo
We used a Markov chain Monte Carlo (MCMC) approach [1] (Algorithm 1) to solve for three-dimensional deformation and their relative errors. We estimated the East North and Vertical velocity components (in mm/yr, Figure 2) for each pillar on the northern and southern side of the bridge, at the deck level and at each pillar ground level. Using Equation (7), we estimated for each incidence angle θ and heading angle ϕ (Table 2), the probability density functions (PDF) from histograms of retained solutions ( Figure 2). We ran the MCMC algorithm iteratively until reaching three million kept solutions. The most probable parameters were estimated from the PDF together with 95% confidence intervals. Figure 2 shows the inverted parameters for the CSK and Sentinel-1A data. When too many parameters are inverted at the same time, the PDFs for several of the parameters do not form smooth, Gaussian-like distributions, but instead are multi-modal. We found a stable solution only when three independent line of sight were available (i.e., CSK ascending, descending and SNT descending). The input LOS parameters are listed in Table 1 while the obtained results are shown in Figures 2 and 3.   Figure 4A. West/East are referred to points located at the deck level (located on the west/east side) near the strands while average refers to the average trend of the deck. "/" indicates no persistent scatterer (PS) coverage in the area. Systems 9-11 refer to the balanced system numbers as listed in Figures 1 and 4. North, South and Ground indicates where the pixels are located on the northern side, southern side, or at the base of the bridge, respectively. 2018-2010 trends for CSK data are presented in Table S2.    Table 2. MCMC outputs are listed in Table 3 and Figure 2.

Results
We processed 130 and 148 images acquired in an ascending and descending geometry, respectively, from the Italian COSMO-SkyMed (CSK) constellation to infer ground deformation   Table 2. MCMC outputs are listed in Table 3 and Figure 2.  Table S2.

Results
We processed 130 and 148 images acquired in an ascending and descending geometry, respectively, from the Italian COSMO-SkyMed (CSK) constellation to infer ground deformation during 2009-2018. . While both the SNT and CSK descending datasets captured the temporal behavior of points located on the bridge deck, the SNT ascending dataset only recorded points located on the ground and for this reason has been discarded from our analysis (Figure 4 and Figure S2). Because of the differential nature of the MT-InSAR measurements, the time series of cumulative deformation were calculated relative to a reference area. Different reference areas were used, in order to highlight different phenomena characterizing the bridge deck and its base. In this analysis we used four reference areas (Figure 4, Panel A, Panel F): the first reference area was located on the west side of the bridge, close to pier 5 and 6, at the deck level (REF1), while the other three were positioned on the north east side of each balanced system 9-11 (REF2-4) ( Figure 1A). REF1 was used as a reference area for analyzing the behavior of the entire bridge and the interaction of each balanced system with the ground, while REF2-4 were used for characterizing the three-dimensional displacement field of the deck, independently for each system. Time series of deformation using REF1 show that average deck deformation recorded from ascending and descending data were characterized by linear trends during the entire satellite acquisition campaign (grey pattern in Movies 1-5, see supplementary material). Movie 1 and Figure 5 show the time-series deformation for the descending line of sight (LOS) of the CSK satellites. The deformation refers to the northern side of the bridge for the three balanced systems at the deck level. In Movie 1 all the PS were located on the north side of the bridge (Figure 4, Figure 5). We found one PS (red dot in Movie 1 and Point CSK DSC 1 in Figure 1), located at the junction between the northwest deck level and the northwest strand, which was characterized by an abrupt acceleration starting March 12th 2017. The LOS deformation rate of this point increased by seven times from 10 ± 2.8 mm/yr to 70 ± 2.8 mm/yr away from the satellite, as shown in Movie 1 and Figure 1. The descending SNT dataset (Movie 2) recorded similar average trends (9.4 ± 2.8 mm/yr) for all the scatterers located on the northern side of the bridge, although there are no persistent scatterers in the SNT dataset covering the area of the bridge characterized by an anomalous speed increase. However, SNT descending scatterers located on the southern side of the balanced system 9 highlighted an area close to the junction between the southern west deck level and the southern west strand that was characterized by a motion toward the satellite at a rate of about 10.3 mm/yr and two deformation transient steps of 1 and 2 cm occurring on January 21st 2017 and April 28th 2018, respectively (Movie 3). These steps were incompatible with phase unwrapping errors only occurring at multiples of half wavelength (2.8 cm for C-band 5.6 cm Sentinel-1A/B wavelength). The balanced system 9 moved away from the satellite at a constant rate of 6.7 ± 2.8 mm/yr. PS recorded by the CSK ascending were located on the southern side of the balanced systems (Movie 3) and moved by about 6.6 ± 5.8 mm/yr away from the satellite, confirming the SNT average deformation trend of the deck southern side of the balanced system 9. As for CSK descending, we found two PS (red and blue dots in Movie 4) located at the junction between the south-west (Point CSK ASC 1 in Figure 1) and east deck level (Point CSK ASC 2 in Figure 1) and the south-west and east strand, respectively, which were characterized by an abrupt acceleration starting March 12th 2017. The LOS deformation rate of these points increases by 10 times from -3.3 ± 2.8 mm/yr to 33.3 ± 2.8 mm/yr away from the satellite for point CSK ASC1, and from -1.6 ± 2.8 mm/yr to 13.3 ± 2.8 mm/yr for point CSK ASC2 (Figure 1). Movies 6-8 show time-series of cumulative deformation at the ground level (using REF1) for CSK descending, SNT descending and CSK ascending, respectively. Given that neither GPS nor other ground truths were available, it is not possible to establish if REF1 was affected by displacement. This made the single balanced system analysis difficult to interpret when using REF1. We interpret this behavior later on when inferring the three-dimensional deformation field. We found that for the deck 9, PS located on the north and southern side moved by 1.2 ± 2.8 mm/yr and 2.8 ± 2.8 mm/yr towards the sensor, respectively, while PS located at the ground level moved by 1.5 ± 2.8 mm/yr away from the satellite. The points located on the north side of both decks 10 and 11 were characterized by a 0.6 ± 2.8 mm/yr motion toward the sensor along the LOS. An insufficient number of PS was found on the southern side of the deck 9 and 10, making the trend analysis not feasible in this area. The ground under pier 10 was moving at 0.5 ± 2.8 mm/yr away from the sensor, similarly to the ground under pier 9.
A comprehensive interpretation of LOS displacements for each dataset remains a difficult task. For this reason, it is fundamental to infer three-dimensional displacement field and reproject our data into East, North and Vertical deformation trends. We analyzed LOS measurements from three independent datasets, using the reference areas REF2-4 to characterize the three-dimensional deformation of each deck prior to collapse and filter out deformation potentially characterizing REF 1. We modeled the multi-geometry inferred deformation using a Markov chain Monte Carlo (MCMC) approach to solve for three-dimensional deformation of the bridge deck. Each pier was designed as an independent balanced system made of a 171 m continuous deck connected to two simply-supported Gerber beam systems of 36 m. It is therefore reasonable to assume that each deck and beam can be approximated as an independent rigid block. Under these assumptions we infer three-dimensional deformation trends of each side (south and north) of the deck and calculate its deformation along three main independent axes. We use an MCMC approach to infer a three-dimensional displacement field at each side of the deck using available LOS directions and Equation (7) presented in the methodology section. The MCMC approach indicated that the northwest side of deck 9 moved (starting February 2015) at a rate of about 3.5 ± 2.1 mm/yr west, 14.1 ± 7.5 mm/yr south, and 4.6 ± 1.1 mm/yr down compared to the northeast side (i.e., REF2). The southwest side of deck 9 moved at a rate of about 7.1 ± 2 mm/yr west, 9.8 ± 6.1 mm/yr north, and 0.9 ± 1.1 mm/yr down compared to the northeast side, while the southeast side of deck 9 moved at a rate of about 7.5 ± 3.5 mm/yr west, 6.15 ± 6 mm/yr north, and 1.4 ± 1.1 mm/yr down compared to REF2. Figure 3 shows a three-dimensional characterization of the bridge deck and Table 3 indicates the numeric values inferred using the MCMC approach (Figure 2), suggesting that the deck 9 was characterized since 2015 by increased deformation accumulation on each side of the deck compared to deck 10 and deck 11. The trend inversion of the three-dimensional displacement also showed that on deck 9 the deformation accumulated mostly in the north-south direction, and decks 9-11 were relatively stable compared to the ground underneath their corresponding pier.

Discussion
This study highlights how the availability of new constellations of SAR could be applied to bridge deformation monitoring. While these techniques cannot decisively distinguish between stress accumulation or material degradation processes, they are useful to detect structural distress signs. A dedicated study of the available literature also confirms how these methods proved to be crucial in characterizing three-dimensional displacements with millimetric accuracy and capable of highlighting anomalous behaviors in near-real time [17,18]. The comprehensive multi-sensor dataset available for the Morandi bridge collapse allows us to make several unique considerations regarding the overall quality of the presented data, the future use of such data for a near-real time monitoring system and bridge construction guidelines.
Our InSAR single line of sight multi temporal analysis coupled with the three-dimensional MCMC inversion scheme show an increased accumulation of relative displacements of deck 9 compared to decks 10 and 11, starting from 2015. The analysis does not show macroscopic

Discussion
This study highlights how the availability of new constellations of SAR could be applied to bridge deformation monitoring. While these techniques cannot decisively distinguish between stress accumulation or material degradation processes, they are useful to detect structural distress signs. A dedicated study of the available literature also confirms how these methods proved to be crucial in characterizing three-dimensional displacements with millimetric accuracy and capable of highlighting anomalous behaviors in near-real time [17,18]. The comprehensive multi-sensor dataset available for the Morandi bridge collapse allows us to make several unique considerations regarding the overall quality of the presented data, the future use of such data for a near-real time monitoring system and bridge construction guidelines.
Our InSAR single line of sight multi temporal analysis coupled with the three-dimensional MCMC inversion scheme show an increased accumulation of relative displacements of deck 9 compared to decks 10 and 11, starting from 2015. The analysis does not show macroscopic deformations characterizing the base of the piers 9-11. The absence of such deformations is also confirmed by the IG18 report by the Commission of the Ministry of Infrastructures and Transportation [29]. Furthermore, there are several PS recorded from different dataset and localized on different areas of deck 9 showing an increased acceleration. These PS are located close to the strands on the north and south west side of deck 9 respectively (Figure 1, Figure 5, Movie 1 and Movie 3) and show an anomalous behavior consisting in an abrupt acceleration starting 12th March 2017 at the junction between the deck and the strand located on the northwest, southwest and east sides of deck 9 respectively.
At the time of writing, the causes of the Morandi bridge collapse have not been officially identified [29]. No consensus has been found yet on the failure mechanism and the conditions which led to that mechanism. Different potential failure mechanisms have been evaluated [30].
Possible collapse modes included: (a) The loss of capacity of one of the stays connecting the deck with the antenna of balanced system 9; (b) the shear failure of the cantilever beam of system 9 or the connecting simply supporting beam; (c) a local failure at the connection between the stay and the deck. These mechanisms have been considered in relation to their ability of triggering the complete failure of one of the stays and the subsequent failure of the whole system 9 [30]. With reference to the pre-existing conditions which could have led to one of these mechanisms, a general material degradation has been documented for both the stays and the deck beams [23,29]. Our results show deformation patterns compatible with both the accelerating degradation of a stay strand, leading to its elongation, and a deterioration of the deck-stay connection [31].
The IG18 analysis does not highlight any of these deformation trends, indicating that previous analyses of CSK satellite-based datasets have been reported to give inconclusive results on pre-collapse variations of the displacement field in time. The discrepancy between the two different results can be ascribable to several reasons such as different temporal coverage (i.e., the IG18 only covers the period 2016-2018 whereas our CSK analysis covers 8 years 2010-2018, different coherence threshold used for the data selection (our analysis only considers pixels with a temporal coherence higher than 0.6 whereas the coherence threshold for IG18 is not known), absence of line of sight diversity and the lack of SNT data, together with the three-dimensional inversion scheme used in this study.

Conclusions
The synergistic use of CSK, SNT and Envisat played a key role in our study. The SNT data have been acquired at a rate nearly three times faster than CSK (6 days repeat vs 16 days repeat for SNT and CSK respectively), while the large CSK orbital tube and better CSK spatial resolution (3 x 3 meters vs. 4 x 22 meters range and azimuth resolution for CSK and SNT respectively) allowed us to better geolocate the analyzed scatterers in the three-dimensional space by a factor of three, as confirmed by independent studies [32,33]. The SNT tight orbital tube introduces a further level of difficulty when geolocating the observed scatterers located on the northern/southern side of the bridge. In our study we used a reference ground control point for georeferencing our data and only considered scatterers along the bridge aligning with optical images. The points not aligning with the optical images are not located at the deck level and hence discarded from our analysis. However, the severe limitations introduced in urban area by a small orbital tube represent a straightforward advantage for geophysical applications in mountainous area characterized by severe slopes introducing geometrical decorrelations. While the Envisat historical data allowed disclosing the presence of any pre-existing long-term trend characterizing the bridge pillars, no data were available for the deck level, principally due to the spatial resolution and the repeat time of the Envisat satellite (25 x 5 meters ground range/azimuth resolution and an average of one acquisition every three months). The SNT dataset revealed to be fundamental for the three-dimensional deformation decomposition when merged with the CSK data. The absence of scatterers located at the deck level for CSK descending (on the southern side of the bridge) and SNT ascending leads to the suggestion of introducing artificial scatterers such as corner-reflectors or transponders embedded within the bridge design. These devices could bridge the gaps between the SNT and CSK constellations, allowing a precise characterization of infrastructures deformation and the capability to monitor specific areas of the bridge according to the system design.
Overall this study shows that either CSK and SNT InSAR measurements are now accurate enough to measure millimeter scale deformation processes, providing information, that can help mitigate deformation-related hazards to large-scale infrastructure.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/11/12/1403/s1, Table S1. LOS trends from different satellites and line of sight recorded during the period 2015-2009. Table S2. Dataset characteristics used in the MCMC inversion. Figure S1. Envisat multi-temporal line of sight (LOS) analysis. Figure S2. Sentinel-1A/B (SNT) ascending multi-temporal Interferometric SAR analysis of the Morandi Bridge. Movie1-8. Time-series deformation for the ascending/descending line of sight (LOS) of the CSK/Sentinel satellites recorded on the north/south side of the bridge, respectively.