Ambient vibration‐based finite element model updating of an earthquake‐damaged masonry tower

This paper presents a vibration‐based model updating procedure for historical masonry structures which have suffered severe damage due to seismic events. This allows gathering in‐depth insights on the current condition of damaged buildings, which can be beneficial for the knowledge of their actual structural behaviour and, consequently, for the design of repairing and strengthening interventions. The methodology, based on the experimentally identified modal parameters, is tested on the San Felice sul Panaro medieval fortress, which was heavily damaged by the 2012 Emilia earthquake. The finite element mesh of the structure in its post‐quake condition is generated by means of a nonstandard semi‐automatic mesh generation procedure based on a laser scanner points cloud. Ambient vibration testing is performed on the main tower of the fortress. Mechanical properties of the tower and the level of connections with the rest of the fortress in its current damaged state are investigated. To fully characterize the actual behaviour of the tower in operational condition, mesh elements corresponding to the damaged masonry are identified and different material properties are assigned to them. This allows to account for the effect of damage and cracks, which appeared essential in the calibration process. The updating procedure is carried out by means of an advanced surrogate‐assisted evolutionary algorithm designed for reducing the computational effort.


INTRODUCTION
an evolutionary strategy with a quadratic response surface to evaluate the structural parameters that potentially minimize the difference between numerical and experimental results.
The paper is organized as follows. First, the San Felice sul Panaro fortress is described in Section 2. Ambient vibration testing and modal parameter identification are presented in Section 3, whereas Section 4 describes the FE model of the fortress generated by means of the CLOUD2FEM procedure. Finally, the updating of the FE model using the DE-S algorithm is presented and discussed in Section 5.

DESCRIPTION OF THE CASE-STUDY
The case under study is the main tower of the San Felice sul Panaro medieval fortress (Figure 1), located near the city of Modena. Such a tower is called Mastio because of its dominant dimensions compared to the rest of the building (see Figure 1a). The tower, 32 m high with an almost square plan with sides of slightly more than 10 m, is composed of seven levels: three cross vaults and three timber decks as well as a timber trussed roof. The thickness of the tower trunk walls ranges from 2.5 m at the bottom up to 1.25 m at the top. Several openings of different sizes are irregularly placed along the tower. Furthermore, it is characterized by the presence of a crowning in its upper part, which presents a larger plan, achieved through masonry corbels, which allows to enlarge the perimeter overhanging from the wall below. Along its west side, the tower is adjacent to the south entrance building, whereas in its north side, it is linked to the curtain wall.
In 2012, the investigated structure was hit by the Emilia earthquake with two main shocks of magnitude M W = 5.86 (May 20) and M W = 5.66 (May 29). [40] The epicentres of the first and the second main shocks were located at approximately 10 and 5 km far from San Felice sul Panaro, respectively. After the seismic sequence, the collapse of the four minor tower roofs was observed and cracks of different relevance appeared on all the fortress structural elements extensively, as shown in Figure 1b. In Cattari et al., [4] an accurate description of the Emilian medieval fortresses damage mechanisms is reported. Particularly, the Mastio most relevant damage consists of diagonal cracks, clearly visible in the lower half of the south and north front (Figure 2a,b). In order to prevent further collapses of the structure after the seismic sequence, the San Felice sul Panaro municipality commissioned first-aid safety interventions. In particular, the cracks of the main tower (and of the remaining structures) have been partially filled with lime and polyurethane (Figure 2a,b). Thereby, 15-mm-diameter steel strands have been diffusely inserted into the main tower sections (highlighted with green circles in Figure 2c), inside 60 mm large holes drilled into the walls and then grappled through lime mortar and pozzolan.

MODAL IDENTIFICATION FROM AMBIENT RESPONSE
Ambient vibration tests were conducted on the Mastio of the San Felice sul Panaro fortress in July 2016, to measure the dynamic response in operational conditions and identify its modal properties. As the fortress was significantly damaged by the 2012 earthquake, the modal identification refers to the dynamic behaviour of the Mastio in damaged conditions.   Instrumented levels (L1-L4) and layout of the accelerometers (A1-A10) during the ambient vibration testing The dynamic acquisition system was composed of 10 uniaxial piezoelectric accelerometers (seven PCB/393B12 and three PCB/393B31), with a dynamic range of ±0.5 g, a bandwidth ranging from 0.15 to 1000 Hz and a resolution of 8 g (PCB/393B12) and 1 g (PCB/393B31). The accelerometers were connected to a National Instruments acquisition system for data storage and system management. The sampling frequency was set to 200 Hz. The response of the tower was measured simultaneously in seven points belonging to four levels along the height of the main tower. In each measuring point, one or two accelerometers were placed, for a total of 10 measurement channels (A1-A10 in Figure 3). The accelerometers were installed on the inner walls by means of metal plates and screws, as shown in Figure 4. Figure 5a presents a typical acceleration time series recorded at the upper instrumented level (L4). The measured acceleration ranges between ±15 mg (corresponding to ±0.15 m/s 2 ), stating the low level of ambient excitation during the test. The corresponding power spectral density (PSD) function is shown in Figure 5b. The modal identification is performed applying the Enhanced Frequency Domain Decomposition method to the acquired accelerations. [41,42] First, the PSD matrix of the acquired accelerations is calculated and decomposed with the singular value decomposition (SVD) method. Each singular vector of the PSD matrix represents the jth mode shape whereas the corresponding singular value is the amplification factor, that is, the structural response amplification under dynamic loads. The jth natural frequency  is identified from the peak of the PSD graph. Finally, the damping ratio is estimated through the logarithmic decrement. The reader is referred to other works [41,42] for all the details about the method. Table 1 presents the first five estimated natural frequencies and damping ratios whereas the corresponding mode shapes are shown in Figure 6. Two closely spaced modes are identified around 1.75 Hz. These modes are dominant bending and involve flexure in W-E direction ( Figure 6a

FE MESH GENERATION AND NUMERICAL MODELLING
The FE model of the structure under study has been generated by means of a non-standard mesh generation procedure called CLOUD2FEM. [12,38] Such an innovative method semi-automatically transforms 3D points cloud of complex buildings into 3D FE meshes, minimizing the user-time investment. The procedure details are briefly revisited in this section as the FE mesh is used in the model updating process.
Given a points cloud, obtained by means of laser scanner or photogrammetric surveys of the inner and outer surfaces of a building, it can be processed by reducing the points density and by generating the triangular irregular network (TIN)  mesh, which is a standard routine. Such a mesh of the surfaces is broken down by means of a subdivision of the 3D domain into bi-dimensional subdomains by slicing it perpendicularly to an opportune direction (typically the vertical direction) with a certain step. The slicing step is chosen according to the complexity of the building along the slicing direction. By using a concave hull algorithm, [43] the boundary polygon that encloses the points of each slice is computed. Then a filled region for each slice of the building is extracted and idealized as a digital image composed of pixels with a certain resolution. Such slices are stackable as the digitalization is performed on a fixed space region. The subsequent stacking of the slices generates voxels. Finally, each voxel is automatically transformed into an eight-node hexahedral finite element and, therefore, the structure is completely discretized as an unique continuum composed by evenly spaced eight-node hexahedral elements.
As far as the investigated monument is concerned, the municipality of San Felice sul Panaro, after the first-aid safety interventions, commissioned a fine terrestrial laser scanning (TLS) survey of the damaged building in order to acquire a snapshot of the after-quake structural condition (see Figure 7). The surveyed raw points cloud (composed of over 40 millions of points) is shown in Figure 7a, whereas the processed TIN mesh is shown in Figure 7b. On this mesh, a slicing operation has been carried out (see Figure 7c), using a vertical gap of z = 25 cm. One hundred twenty-one slices have been extracted and, after their processing with a concave hull algorithm, digitalized (see some examples in Figure 7d) with a bi-dimensional resolution in the horizontal plane of 25 cm× 25 cm, as suggested in Castellazzi et al. [44] Indeed, in Castellazzi et al., [44] the authors carried out a natural frequencies comparison between different CLOUD2FEM-based mesh sizes of the main tower (as isolated tower), obtaining a very good performance of the 25 cm × 25 cm × 25 cm mesh. Furthermore, the effectiveness of the meshing approach has been also investigated in Castellazzi et al. [38] through a comparison with a very detailed CAD-based model. As a result, the resolution 25 cm × 25 cm × 25 cm was found to be the best compromise between result accuracy and computational effort. In fact, although this mesh dimension does not accurately reproduce each small architectural detail, it guarantees a good accuracy in terms of global dynamic response having at least five hexahedral FEs in the thickness of the main tower walls. [12] The resulting mesh, depicted in Figure 8, is characterized by 409,300 hexahedral finite elements (each one 25 cm × 25 cm × 25 cm) and 1,512,444 dofs.
It has to be pointed out that, although the laser scanner survey has been carried out on the damaged fortress, its capability to collect information about the cracks is substantially limited. Indeed, the chance of picking points in the depth of a crack is extremely challenging and strongly depends on the crack width and on the relative angle between the surveyed surface and the laser direction. Furthermore, most of the cracks was partially filled in the first-aid interventions. In addition, during the idealization of each slice as a digital image composed of 25-cm-large pixels, every geometric detail smaller than this dimension is approximated.
Modelling the floors and vaults of monumental masonry structures has always been a challenging task. [45] In particular, it is commonly accepted that in order to perform 3D seismic analyses of masonry buildings, equivalent diaphragms can be used to model vaults. [46] In the adopted FE model, floors and vaults are automatically meshed by means of a jagged 3D representation of the original geometry and the mesh accuracy can be considered satisfactory aiming at the global structural response. Particularly, an isotropic elastic material has been assumed to roughly model deck timber elements with the values 8,000 MPa, 2,918 MPa, and 415 kg/m 3 for Young's modulus, shear modulus, and density, respectively. [47] The fortress was surrounded by a moat, and, therefore, the outer ground level is approximately 3.5 m lower than the inner one. Accordingly, clamped boundary conditions have been applied in all the nodes located at the moat level, whereas the elements located within the courtyard have been modelled through an elastic continuum to coarsely take into account the presence of terrain. In particular, a linear elastic material with Young's modulus, shear modulus, and density equal to 935 MPa, 316 MPa, and 1,200 kg/m 3 , respectively, has been considered for terrain. [48] For the masonry elements, a linear elastic material with density and Poisson's ratio equal to 1,800 kg/m 3 and 0.2, respectively, has been adopted. According to the Italian Codes, [49] the Young's modulus of rubble masonry elements made of solid bricks and lime mortar is expected to lie in the range [850-1,250 MPa]. However, the actual value of the masonry elastic modulus is identified from the model updating procedure in Section 5.
Finally, roofing wooden structures have been considered as concentrated mass on the Mastio top elements, that is, where the roof is borne.
Although in this study the attention is focused on the Mastio modal identification and model updating, the modelling of the whole fortress is considered essential by the authors as the tower under study is clearly non-isolated and the FIGURE 8 Generated finite element model: The magnified portion shows the mesh discretization interaction between adjacent structures can play a considerable role in the evaluation of natural frequencies and modal shapes. In particular, in previous studies, [50,51] it arises that the stiffness of the adjacent parts and their mutual level of connection with the tower strongly influences its dynamic response. This aspect is further investigated in the following section. Therefore, by modelling the whole monument, the actual stiffness of the adjacent parts of the tower is directly accounted for.
To conclude, it has to be particularly stressed that the the adopted mesh generation approach directly considers the after-quake geometry of the structure, including for instance partial collapses, that is, the configuration of the structure when the ambient vibration testing was carried out.

MODEL UPDATING AND STRUCTURAL PARAMETER IDENTIFICATION
The FE model of the fortress is calibrated with respect to the experimental results so that the modal properties of the main tower agree as close as possible with the experimental ones. To this aim, a set of unknown structural parameters is evaluated from the minimization of an objective function defined as the difference between measured data and numerical predictions. Due to the complexity of the numerical model, the evaluation of numerical modal parameters is highly time-consuming and the success of the optimization problem strongly depends on the efficiency of the optimization algorithm. Hence, an improved surrogate-assisted evolutionary strategy is adopted to reduce the number of objective function evaluations and find a compromise between local and global search. The optimization algorithm is described in Section 5.1 whereas the identification of the structural parameters is presented in Section 5.2.

The DE-S algorithm
Genetic and evolutionary algorithms are widely used to solve global optimization problems. Their architecture is designed for large-scale problems and allows to avoid local minima. The main drawback is that a large number of function evaluations is often required to reach the convergence. Surrogate-assisted evolutionary strategies [52] use efficient computational models, such as response surfaces (RSs), high polynomial functions, or Kriging models, [53] to approximate the objective function. They aim at evaluating those individuals that potentially have a good prediction of the objective function value. The introduction of a second-order surrogate in the differential evolution (DE) algorithm is proposed in Vincenzi and Savoia. [54] In the DE-S algorithm, proposed by Vincenzi and Gambarelli, [39] a proper infill sampling strategy is introduced to further reduce the number of objective function evaluations. The candidate points are chosen trying to find a compromise between local and global search, that is, enhancing both the accuracy in the region of the optimum predicted by the surrogate (local exploitation) and the global exploration. The DE-S algorithm is summarized in the following.
The optimization process is based on the simultaneous adoption of NP vectors x i,G , with i = 1, ..NP. Each vector x i,G has dimension D, being D the number of unknown structural parameters, whereas G indicates the Gth population of vectors. First, the objective function values are evaluated for the initial population of vectors, randomly chosen in the search space. At each subsequent iteration and for each vector x i,G , a trial vector v i,G is generated. To this aim, NP subsets of NS vectors are built, with NS < NP. Each subset contains the vector x i,G and other NS − 1 vectors randomly selected among the remaining vectors of the Gth population. Each subset is used to calibrate a second-order response surfaceĤ as a local approximation of the objective function H: where Q is a D × D coefficient matrix collecting the quadratic terms, L is a D-dimensional vector of linear terms, and 0 is a constant. Applying the least square estimation methods from the NS evaluated points, the parameters calibrating the RS model are evaluated. Depending on the shape of the approximating functionĤ, two possibilities occur. If the response surface is convex, the mutant vector v i,G is defined as the minimizer x * of the approximating function, that is, On the contrary, if the response surface is nonconvex, the mutant vector v i,G is obtained from the Mutation operation, based on a linear combination of vectors [55] : where r 1 , r 2 , r 3 are integer numbers randomly selected in the range [1-NP], with r 1 ≠ r 2 ≠ r 3 , and F is a positive constant controlling the amplitude of the mutation, normally assumed smaller than 2. The Crossover operation is applied to increase the diversity of vectors and escaping local optima. Each mutant vector represents a candidate for a new evaluation. For each candidate v i,G , the score is defined as the weighted sum of two criteria. The first one depends on the objective function value predicted by the surrogate. The second criterion depends on the distances of the candidate from the points where the objective function has been already evaluated.
The two scoring criteria are combined as follows: where V R (v i,G ) and V D (v i,G ) are the scores and w R and w D the weights associated to the two criteria. Candidate points with low vales of W are preferred. The score V R depends on the objective function value predicted by the surrogate as follows: where s f is the minimum of the quadratic approximationĤ (that is the minimum estimated by the surrogate), H best is the best objective function value at the current iteration, whereas is the minimum of the evaluations of the NS points (x s,G ) selected for the calibration of the response surface. The score V D is derived from the distances of the candidate point to the points where the objective function has already been evaluated: where min and max are the minimum and maximum weighted distances of the candidate point to the n points where the objective function is evaluated: The weighted distances D are calculated from the Euclidean distances d and the objective function values as The weights w R and w D are chosen in the range 0 ≤ w ≤ 1, with w R + w D = 1. If w D is close to 1, the global exploration prevails on the local one and candidate points placed in a rather unexplored region of the parameter domain are preferred. On the contrary, if w R is high, candidate points with low objective function values are preferred (local exploitation).
Results of analyses performed in Vincenzi and Gambarelli [39] suggest the adoption of the following weight factors: w R = 2∕3, w D = 1∕3 if the surrogate presents a minimum (Case A in the flowchart of Figure 9) and w R = 0, w D = 1 when the surrogate search fails (Case B). Indeed, in Case A, it is worth investigating the region close to the optimum predicted by the surrogate whereas in the second case, the global exploration of the parameter domain is preferred.
To reduce the number of the objective function evaluations, a subset of NH candidates (with NH < NP) is selected. Candidates with the lowest score are preferred and only for this subset of points the objective function is evaluated. Finally, in the Selection operation, all points (candidates and points of the current population) are ordered depending on their objective function value. The NP points with the lowest value of the objective function are selected to be the members of the next population. The flowchart of the DE-S algorithm is reported in Figure 9. The reader is referred to Vincenzi and Gambarelli [39] for all details about it.

Identification of the structural parameters
The objective function H is a measure of the difference between numerical (f num , num ) and experimental (f sper , sper ) natural frequencies and mode shapes: MAC( num, , s er, ) where (0 ≤ ≤ 1) is a weighting factor that determines the relative importance between frequency and mode shape residuals, x is the vector of unknown structural parameters, and N is the number of modes considered in the calibration procedure. Finally, MAC is the modal assurance criterion, [56] representing the correlation between two modal vectors. First, it is important to identify which modes should be included in the calibration process. Indeed, in ambient vibration testing, higher frequencies are often obtained with less accuracy than the lower order ones. Minimizing the error between experimental and numerical modal properties for higher modes may prevent matching the lower modes of vibration. [57] In this study, the first three natural modes (Modes 1-3 in Table 1) are accounted for in the calibration procedure, whereas Modes 4 and 5 are adopted for the purpose of validation of the updated FE model. This is in-line with the fact that generally, only the first natural modes are relevant for seismic analyses.
In Forghieri et al., [25] the influence of the weighting factor on the identified structural parameters is investigated evaluating the optimal solutions forming the Pareto front. Adopting two different selection criteria, an optimal value of equal to 0.4 is obtained, although results are observed to be similar for in the range [0.1-0.8]. Hence, in this study, the weighting factor is set to 0.4.
Moreover, a proper selection of structural parameters is fundamental for the success of the calibration process. In particular, those parameters whose values are uncertain and that potentially have a considerable effect on the vibration response of the structure are to be selected.
The FE model of the fortress has been generated from the 3D points cloud obtained after a fine TLS survey of the fortress. This means that the physical dimensions of the masonry elements and floors are defined with reasonable accuracy. On the contrary, the material properties of the major structural components are more uncertain. Indeed, due to the serious damage of the fortress and its age of construction, the effective stiffness of the masonry is uncertain.
The first attempt to update the numerical model (Model updating #1) is carried out considering a homogeneous distribution of the masonry elastic properties. In particular, the structural parameter considered in the optimization procedure is the equivalent elastic modulus of the cracked masonry of the Mastio, indicated as E M . Table 2 presents the optimization analysis results, including the range of variation for each updating parameter and the identified value, the average frequency error, and the average MAC value. The average frequency error and MAC value are calculated with reference to Modes 1-3 and 4-5. The first are the natural modes against which the FE model is calibrated whereas the second are used to evaluate the goodness of fit of the calibrated model. The numerical frequency and MAC value for each natural mode are reported in Table 3, together with a comparison to the experimental results. As far as the calibrated modes (1-3) are concerned, a pretty good match between numerical and experimental frequencies is obtained, with errors lower than 6%. Modes 1 and 3 show MAC values higher than 0.90, whereas a slightly lower correlation is observed for Mode 2, characterized by a MAC value of 0.86. Larger differences between numerical and experimental results are obtained for Modes 4 and 5, in relation to both natural frequencies (Mode 4) and mode shapes (Mode 5).  With regard to the identified structural parameters, an equivalent elastic modulus of the masonry of about 825 MPa is obtained. The identified value is lower than the expected one due to the presence of significant cracks in the masonry.
As the updated model does not describe the actual behaviour of the tower with respect to Modes 4 and 5, additional structural parameters are accounted for in the optimization procedure. Due to the complex historical evolution of the fortress, its monolithic behaviour and the efficiency of connections between adjacent walls cannot be assured. Furthermore, some important damage can be observed in correspondence of the interfaces between the Mastio and the rest of the fortress. Therefore, the elastic moduli of the portions of walls connecting the Mastio to the fortress on the west (E CW ) and north (E CN ) side are considered as additional unknown parameters (highlighted in Figure 10).
A sensitivity analysis is performed to evaluate how each structural parameter influences the dynamic behaviour. The analysis is carried out varying the structural parameters one-by-one and setting the others to 825 MPa (i.e., the elastic modulus obtained from model updating #1). Figure 11 presents the value of the objective function evaluated from Equation 11 as well as the partial contributions of natural frequency and mode shape residuals to the objective function. Figure 11a,b is obtained considering Modes 1-3, whereas Figure 11c considering Modes 4 and 5.
As expected, the variation in the masonry elastic modulus E M affects the natural frequencies and not the mode shapes (Figure 11a). On the contrary, the variation in the elastic moduli E CW and E CN of the connection elements causes changes in modes shapes and slightly affects natural frequencies. In particular, the objective function values for different E CW calculated considering Modes 1-3 and 4-5 are reported in Figures 11b,c, respectively. When E CW decreases, values of the objective function calculated from Modes 1-3 increase whereas those calculated from Modes 4-5 decrease. This highlights the non-negligible contribution of the interface element stiffness.
Hence, a second calibration (Model updating #2) is carried out considering as uncertain parameters the elastic modulus of the cracked masonry E M and the elastic moduli of the connection elements E CW and E CN . The updated model shows an equivalent elastic modulus of the masonry equal to 932 MPa whereas the connection elements are character-  ized by low (<50 MPa) elastic moduli ( Table 2). Compared to the model updating #1, the elastic modulus of masonry E M increases to make up for the reduction of the connection stiffness. Moreover, the calibrated modes (Modes 1-3) show a slight improvement in the mode shapes against a slight worsening in terms of natural frequencies. Similar to the model updating #1, a low correlation between experimental and numerical results is observed for Modes 4 and 5.
From the previous calibrations, it follows that the equivalent elastic moduli of the cracked masonry is not the appropriate parameter to accurately characterize the Mastio dynamic behaviour. Indeed, the Mastio presents a severe crack pattern whose contribution cannot be accounted for simply defining an equivalent elastic modulus of the cracked masonry. To fully describe the actual behaviour of the main tower in operational conditions, an improvement of the FE model is proposed in this paper. After an accurate survey of damage and cracks, mesh elements corresponding to the damaged masonry are identified in the FE model (highlighted in Figure 12) and a different elastic modulus E D is assigned to them. This allows to account for the effect of damage and cracks in operational conditions, that is, when the external actions are not such as to involve non-linear behaviour and damage and cracks only imply a local stiffness reduction. As the calibration is performed with reference to the modal properties of the Mastio, only the cracks of the Mastio are introduced in the FE model.
The third calibration process (Model updating #3) is carried out considering as unknown structural parameters the elastic moduli of the damaged E D and undamaged E U masonry. The updated model presents elastic moduli of the damaged and undamaged masonry equal to 700 and 892 MPa, respectively ( Table 2). As expected, the elastic modulus of the damaged masonry is lower than the one of the undamaged of masonry. However, the elastic modulus of the damaged masonry is representative of the cracks that, for the grater part, have been partially filled with lime and polyurethane after the 2012 earthquake. Once the FE model is improved accounting for the actually damaged portions of masonry, a much better consistency between numerical and experimental results is achieved. Indeed, the calibrated FE model presents frequency errors between 3.51% and 7.72% and MAC values in the range [0.75-0.95] for the five natural modes.

CONCLUDING REMARKS
In this paper, a vibration-based model updating procedure for historical masonry structures that have suffered severe damage due to seismic events has been presented.
The methodology developed for the calibration of the numerical model, based on the experimentally identified modal parameters, has been carried out on the earthquake-damaged main tower of San Felice sul Panaro medieval fortress. The FE model of the structure in its post-quake condition has been generated by means of a non-standard semi-automatic mesh generation procedure. Ambient vibration testing was performed on the main tower of the fortress. Mechanical properties of the tower in its current damaged state have been investigated. To fully characterize the actual behaviour of the tower in operational conditions, mesh elements corresponding to the damaged masonry have been identified in the FE model and different material properties have been assigned to them, allowing to account for the effect of damage and cracks. This improvement was found to be fundamental aiming at consistently describing the dynamic behaviour of the tower in operational conditions by means of numerical modelling. The updated FE model of the structure was able to guarantee a good accuracy of modal parameters of the concerned modes, which were in close agreement with the experimental results, still preserving the physical meaning of updated parameters.