Country-scale InSAR monitoring for settlement and uplift damage calculation in architectural heritage structures

Abstract

The article proposes a methodology for assessing the development of damage in building structures, subjected to differential settlement and uplift, using the analysis of Interferometry Synthetic Aperture Radar (InSAR) data. The proposed methodology is targeted towards general applicability, capable of providing assessment results for measurements over wide geographic areas and for varying structural typologies. The methodology is not limited to ground movement measurements linked to tunnelling, as is the common case. Instead it extends to the monitoring of arbitrary movement in buildings, for example, due to ground consolidation, water table changes or excavation. The methodology is designed for use alongside patrimonial building databases, from which data on individual building geometry and typology are extracted on a region or country scale. Ground movement monitoring data are used for the calculation of the building deformation, expressed in different types of deformation parameters. The combined use of this data with analytical models for settlement damage classification in building structures enables the assessment in patrimonial building structures, at a country scale. The methodology is elaborated and applied on the patrimonial inventory of Belgium for the evaluation of potential settlement and uplift damage on buildings over a period of nearly three decades. The analysis results are compared to on-site observations.

Introduction

State of the art

Ground settlement and uplift can impose severe demands on buildings in terms of applied deformation. It can cause cracking of non-structural and structural elements alike, leading to the loss of structural integrity, loss of aesthetic or cultural value and a reduction of functionality. Both building and infrastructure are subject to these effects. Buildings and installations located in cities are particularly prone to ground movement due to urban densification, development and intensification of land use. Masonry structures, in particular, can be very vulnerable to differential ground movement due to their weight, stiffness and brittleness.^1,2 Urban development in the coming decades is expected to extract a heavy toll on existing buildings, with historic masonry construction, both monumental and vernacular, being at high risk. Assessing the potential damage to buildings affected by tunnelling in urban areas can be extremely costly and time-consuming. ³

Ground settlement and uplift are caused by a wide variety of human activities and natural phenomena. These include soil consolidation, heavy surface construction, tunnelling, excavation, mining, water pumping (or its cessation due to deindustrialization) and changes in land use.

The changes in the ground surface due to settlement or uplift can be measured and monitored using both terrestrial means, such as geodetic and levelling surveys, ⁴ and air- or space-borne methods, such as persistent scatterer interferometric synthetic aperture radar (InSAR for short).^5–7 The former can be time-consuming and involve manual effort but are able to provide accurate measurements over relatively large areas. The latter can be acquired automatically but may require significant processing and interpretation effort. Thus its application has been mostly limited to single buildings or areas of limited expanse.^7–9

InSAR measurements can be used for the reconstruction of the settlement profile over any desired geographical extent, centred around a single building or over an expansive area. In the absence of surface construction or other impediments influencing its shape, this profile is called the greenfield profile. Following the acquisition of monitoring data, the greenfield profile can be approximated through curve fitting or interpolation.

Certain causes of settlement, such as tunnelling ⁷ or mining ¹⁰ with a known trajectory, depth and cross-section, result in a greenfield curve that can be empirically estimated with good accuracy. In such a case, a relatively small number of measurement points can be sufficient for a reasonable determination of the greenfield profile. However, in the case of other causes, such as water table changes, soil consolidation and surface loads, unpredictable greenfield profiles may be produced. A small number of monitoring points may not be sufficient for the accurate determination of the greenfield profile in such a case.

The determination of the greenfield profile is critical for the calculation of the imposed loads on building structures at the surface due to ground movement. Analytical models for damage classification in buildings subjected to ground settlement rely on the calculation of the settlement function along lines representing the façades of the building or an area representing its plan. An analytical expression for the greenfield curve allows a quick and straightforward application of these damage classification models.

In the case of sparse InSAR measurement points, the reliability of the calculated greenfield curve is reduced due to lack of this data and interpolation errors. Simple polynomial curves may exaggerate or underestimate the actual profile, as well as miss the measured peaks in the vertical ground movement. Radial basis function networks capture the peaks of the profile, which are of high interest, but large distances between points severely reduce the accuracy of the interpolation.

The assessment of strain and damage in structures due to ground movement can be calculated analytically or numerically. Empirical models impose limits to the building strain, tilt, deflection ratio or angular distortion for a set of limit states and are regularly employed in analytical modelling. ¹¹ The superstructure is modelled simply, often as a beam or plate, ¹² and its deformation generally follows the greenfield profile, thus disregarding uplift and slip. Other models are used to classify existing damage according to the level of difficulty in their repair. ¹³ The so-called relative stiffness methods take soil–structure interaction into account, but generally rely on design charts assembled for settlement due to tunnelling.^14,15 Simple limits have been adopted for allowable deformations in buildings by design codes. ¹⁶ Numerical solutions are complicated, expensive and uncommon in the literature.^17,18 Finally analytical approaches, calibrated against a numerical/experimental benchmark, have been proposed ¹⁹ and recently expanded for application in different structural typologies. ²⁰

Objectives

This article aims at the development of a method for the determination of damage to buildings due to ground settlement and uplift at the country scale. The greenfield curve around each building of interest is calculated from InSAR monitoring data over a given time period. Patrimonial information databases are used for the determination of the building type. Analytical damage models are used for the evaluation of the ground movement damage.

The method is, in principle, suitable for the analysis of buildings situated in areas of any measurement point density. It is mostly intended, however, for the analysis of areas where only sparse InSAR data are available. The proposed framework is suitable for the assessment of individual structures over geographic extents up to the country scale. Therefore, it aims at being as general and computationally efficient as possible.

First, the data acquisition, processing of model input and modelling approach itself are presented. Second, the results of the developed model are compared to examples from the literature. Third, the models are employed in a sensitivity study using generated structural and loading data. Finally, an outline of the analysis results for Belgium is given, involving thousands of architectural heritage buildings.

In summary, the calculation process consists of the following steps: (1) acquisition and processing of settlement and uplift data, (2) processing of patrimonial building data, (3) calculation of building deformation due to settlement and uplift and (4) calculation and classification of potential damage in buildings.

InSAR data analysis

Data acquisition and processing

InSAR encompasses a set of methodologies that combine at least the phase information of two Synthetic Aperture Radar (SAR) images for the calculation of the distance change between the radar antennas and the scattering objects on the ground. This change of distance can be attributed to ground displacement due to subsidence, uplift, earthquake or other causes. However, as the scattering objects are likely to be human constructions due to the high amplitude of the backscattered signal needed, the detected movement can either be attributed to a movement of the construction with the ground or the construction itself.

Thanks to the systematic acquisition of SAR data since the launch of the European remote-sensing satellite (ERS) programme in 1991 and the constant increase in computational power, new algorithms for processing the time series of SAR data have been developed. These techniques are grouped under the name multi-temporal InSAR (MT-InSAR). Ferretti et al.^21,22 proposed a method for extracting the phase information of a selection of pixels in a stack of SAR images. Those pixels are first selected using their high amplitude of reflection and further filtered relying either on a model of deformation in time ²² or on their correlation in space. ²³ The pixels that correspond to these criteria are called persistent scatterer (PS) and the technique persistent scatterer interferometry (PSI). After processing, a velocity is calculated for each PS from the time series data. Multiplying the velocity with the measurement period yields the vertical displacement of the PS. An additional quality factor called coherence, varying from 0 to 1, is also calculated.

In this study, a processing chain called the Stanford Method for Persistent Scatterers (StaMPS) ²⁴ was used to process the SAR data. In StaMPS, the PS is processed using their correlation in space which has the benefit of extracting more PS in rural areas than the time model. Three datasets of SAR images (ERS 1/2, ENVISAT and Sentinel 1), each covering the entirety of Belgium (geological map shown in Figure 1), were acquired and processed, spanning a total period of 26 years (1992–2018). These datasets are summarized in Table 1. The three radar satellites are operated in C-band (wavelength of 5.6 cm) and scan the earth with a spatial resolution of 8–14 m. Figure 2(a) represents a map of the annual average velocity for the ENVISAT measurement period (2003–2010) over the area of Belgium.

Figure 1.

Geological map of Belgium.

Table 1.

Characteristics of the processed satellite image datasets.

Satellite	Tracks	Pass	No. of scenes	Acquisition period	No. of PS
ERS 1/2	423, 466, 380	Descending	78, 62, 67	1992–2006	617,028
ENVISAT	423, 466, 380	Descending	73, 51, 62	2003–2010	1,047,106
Sentinel 1	37	Ascending	71	2014–2018	1,048,575

Figure 2.

Ground movement monitoring data for Belgium: average annual velocities of PS from (a) ERS data (1992–2006), (b) ENVISAT data (2003–2010), (c) Sentinel 1 data (2014–2018) and (d) total accumulated soil vertical displacement.

Processing results

Based on a minimum coherence of 0.7, the ERS dataset yielded 617,000 PS, the ENVISAT dataset yielded 1,047,106 PS and the Sentinel 1 dataset yielded 1,048,575 PS. The maximum density of PS is observed in cities, Brussels, for example, reaching 800 PS/km² for the Sentinel 1 dataset. This density means that on average a 1250 m² urban area (a square with a side length of 35.5 m) will contain a single PS. The density drops to 35 PS/km² on average over rural areas. The difference of urbanization between the northern part and the southern part of Belgium, especially in the Ardennes region in the area between the cities of Namur, Liege and Arlon, is clearly visible with the reduction of potential targets reflecting the radar signal. For those rural areas, where PS is absent, an inverse distance weighted (IDW) interpolator was used to fill the gaps between the measured values. Velocity grids of 10 m covering the entire country were then calculated by IDW for each measurement period (Figure 2(a) to (c)). The 10-m grid size was selected for maintaining as high as possible a PS resolution in areas with high density data (cities), and, conversely, in areas lacking PS, allowing the clear visualization of the influence of neighbouring PS.

The total soil displacements over the investigated period can be calculated by summing the contributions measured by each sensor. In the cases of temporal overlap, namely between the ERS and ENVISAT data during 2003–2006, the displacements measured by the two sensors are averaged. For the measurement gap period, namely between the ENVISAT and Sentinel 1 data during 2010–2014, the two time series are extended linearly to cover the gap. Through this process, it is possible to construct the ground movement profile in the entire country over the time period 1992–2018 (Figure 2(d)).

Several trends are clearly visible from the PS velocity maps in Belgium (Figure 2). Among them are (a) the settlement registered in the harbour area of Antwerp, ²⁵ (b) the uplift in the de-industrialized northern part of Brussels, ²⁶ (c) the shift from settlement to uplift in the area of Limburg north and east of the city of Hasselt, caused by the shutdown of the collieries in the region^27,28 and (d) the uplift zone between the cities of Mons and Charleroi, also linked to the closure of former collieries.^25,29

Patrimonial data

Building data

Typically in patrimonial inventory databases, every object is accompanied by data on its age, function, architectural style and typology. These data serve to classify the building in terms of its structural type, as well as to define whether settlement damage models are applicable.

In the case of Belgium, more than 700 typologies were registered in the patrimonial databases of the three federal regions, Brussels, Flanders and Wallonia. The databases include numerous types of non-building patrimonial objects, such as natural artefacts and statues. They further include constructions for which ordinary settlement damage models are not applicable, such as bridges, piers and underground structures.

Having filtered the databases according to the applicability of the models to the set of objects, it is necessary to assign the appropriate model. The age and the architectural style of the building, parameters that are often closely associated, are helpful indicators in this regard. In the cases where more than one date entry is encountered, indicating different construction and intervention phases, the oldest date is considered as the representative of the building. While typically databases will not include readily usable data on the specifics of the structural system and structural materials (such as masonry, structural steel and reinforced concrete), the architectural style helps to distinguish between masonry and infilled frame buildings. Conversely, bare frames are related to a limited number of typologies, such as certain industrial buildings and steel frame gasometers.

Building polygon and offset

The plan of each building is delineated by a simple polygon. For building structures, the edges of the polygon generally correspond to façades or frames in the perimeter of the building. Considering the external bays or façades for the calculations is consistent with the majority of actual construction. The external walls in masonry buildings are normally the stiffest and well-founded structural elements, thus more vulnerable to soil movement. Similarly, the external bays of frame structures bear the heaviest infill walls and are thus the most vulnerable part of the structure as regards to differential settlement. In addition, regardless of structural typology, damage in the façade is highly significant from a preservation perspective. The façade is the most visible part of the structure, making even aesthetic damage on its surface a matter worthy of investigation.

It is known that building structures are influenced by ground movement not only within their perimeter but also within a distance from their perimeter. This is considered in the analysis by expanding the area of influence for each building by an offset extending outwards from the building polygon. PS within the polygon and the offset are equally taken into account in the calculation of damage for the building. The offset serves the additional purpose of increasing the chances of a sufficient number of PS being found per building for analysis. An illustration of the building polygon and its offset are presented in Figure 3. The use of a 10-m offset guarantees that in the absence of PS within the area defined by the offset, at least two interpolated grid values will be captured, thus making the damage calculation always possible.

Figure 3.

Simplified building polygon and defined offset.

The basic data on the building polygons were extracted from the digital cadastre of each federal region. The height of the buildings was determined through the subtraction of the digital terrain model (DTM) from the LiDAR-derived digital surface model (DSM) of Belgium as acquired in 2014. These measurements afford a spatial accuracy of 0.50 m and an elevation accuracy of 0.25 m.

A simplification of the polygons is performed by calculating a rotated minimum bounding rectangle. This shape is defined as the smallest area rectangle that contains all the vertices of the polygon. Subsequently, this rectangle is offset outwards by 10 m and the edges are connected with quarter circles. The rectangular shape of the simplified polygons reduces the computational cost, which can be prohibitive when handling tens or hundreds of thousands of cases, and precludes the generation of self-intersecting polygons by the offset, thus avoiding spurious load computations for the edges. However, for the analysis of single case or a limited number of buildings, the polygon simplification step can be skipped, provided the offset is manually checked.

The rectangular shape used for analysis is not intended to represent the actual building plan with every detail but is adopted as a necessary simplification. While the majority of building structures, particularly those made of masonry, present a simple rectangular plan, some building typologies may feature architectural characteristics that deviate from this simple layout. This includes, for example, large church structures with crossings and side chapels. Nevertheless, the rectangular simplification is able to accurately represent the actual building plan for the majority of cases, particularly for masonry buildings. In addition, regardless of actual plan, the building polygons registered in the patrimonial database often feature dozens of redundant edges, which increase the number of calculations without contributing to calculation accuracy and in fact being detrimental to the evaluation of the results. These cases also benefit from the rectangular simplification. For the limited number of cases where geometrical plan complexity deviates from the simple rectangular shape, this approach is still capable of pin-pointing the orientation of the most, and least, loaded external walls.

Building data analysis

Overall, the percentage of object typologies in Belgium, registered in the databases, for which settlement damage models are applicable, was roughly 70%. The percentage varies between federal regions due to differences in database structure and object inclusivity. This brings the total number of analysis cases to 269,194 buildings.

The number of PS per building for the three satellite image datasets used in this study is evaluated over the case studies of the Brussels region. The distribution of these numbers is illustrated in Figure 4. For the Sentinel 1 dataset, which featured the highest number of PS overall, the number of PS per building is generally low, with only 9.9% of cases having two points or more, even in urban areas where the PS density is highest.

Figure 4.

Cumulative frequency graph of the number of PS per building in Brussels.

Definition of ground deformation curve

For the definition of the ground deformation curve $f(x,y)$ in three-dimensional (3D) space, a conservative approach is adopted. Practical limitations are imposed by the number of PS situated in the area of the polygon and the offset. The coordinates of each PS are defined by its location in the $(x,y)$ plane and its vertical displacement $z$ , equal to the PS velocity times and the length of the measurement period. Absence of PS or existence of a single PS within the area of influence results in inability to determine a curve. Two PS can define a single plane whose dip direction coincides with the horizontal direction of the vector connecting the two PS. Three PS uniquely define a single plane in 3D space. For four or more PS, a least-squares fit of a plane or polynomial surface can be easily calculated. This fit, however, may miss by a large margin that actually measured vertical displacements at the PS location.

As shown from the PS density in urban areas and in the joint analysis of InSAR and patrimonial data, the number of PS per case study is generally limited. Therefore, for every building, the two PS defining the maximum tilt plane are used for the determination of the ground deformation profile in its area of influence. In the absence of PS, the 10-m interpolated grid data are used instead, treating the value of the grid as a standard PS. In the following expressions involving linear algebra, matrices are shown in bold uppercase, vectors in bold lowercase, matrix elements in lowercase with a double subscript index and vector elements with a single subscript index.

The function of the maximum tilt plane is

$$f(x,y)=a+\mathrm{bx}+\mathrm{cy}$$

(1)

Tilt $\omega$ between two points is generally defined as

$$\omega =\frac{s}{d}$$

(2)

where $s$ is the height difference of the points and $d$ is the horizontal distance. For a set of $n$ PS, $P_1,\dots ,P_n$ , each with a pair of coordinates $(x_p,y_p),p\in \{1,\dots ,n\}$ in the $R^2$ Euclidean space, the distance matrix $\mathit{D}$ is

$$\mathit{D}=\left[\begin{array}{cccc}0\hfill & d_{1,2}\hfill & \cdots \hfill & d_{1,n}\hfill \\ \multicolumn{1}{c}{d_{2,1}}& 0\hfill & \cdots \hfill & d_{2,n}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{d_{n,1}}& d_{n,2}\hfill & \cdots \hfill & 0\hfill \end{array}\right]$$

(3)

where

$$d_{i,j}=\Vert P_i-P_j{\Vert }_2$$

(4)

with $\Vert .{\Vert }_2$ being the 2-norm of the $R^2$ Euclidean space.

Each PS has a vertical displacement $z_p,p\in \{1,\dots ,n\}$ obtained from InSAR data. The vertical relative displacement matrix between the PS is

$$\mathit{S}=\left[\begin{array}{cccc}0\hfill & s_{1,2}\hfill & \cdots \hfill & s_{1,n}\hfill \\ \multicolumn{1}{c}{s_{2,1}}& 0\hfill & \cdots \hfill & s_{2,n}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{s_{n,1}}& s_{n,2}\hfill & \cdots \hfill & 0\hfill \end{array}\right]$$

(5)

where

$$s_{i,j}=z_i-z_j$$

(6)

The tilt matrix containing the absolute value of the tilt between each pair of PS $i,j$ is defined as the Hadamard, or entry-wise, division

$$\mathbf{\Omega }=|\mathit{S}ø(\mathit{D}+{\mathit{I}}_n)|=\left[\begin{array}{cccc}0\hfill & {\omega }_{1,2}\hfill & \cdots \hfill & {\omega }_{1,n}\hfill \\ \multicolumn{1}{c}{{\omega }_{2,1}}& 0\hfill & \cdots \hfill & {\omega }_{2,n}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{{\omega }_{n,1}}& {\omega }_{n,2}\hfill & \cdots \hfill & 0\hfill \end{array}\right]$$

(7)

where ${\mathit{I}}_n$ is the identity matrix, whose inclusion precludes division with the zero elements of the diagonal of the $\mathit{D}$ matrix. Therefore, the non-diagonal elements of the tilt matrix are

$${\omega }_{i,j}=\left|\frac{s_{i,j}}{d_{i,j}}\right|$$

(8)

and all diagonal elements ${\omega }_{i,i}$ are zero. The maximum tilt is simply defined as

$${\omega }_{\max }=\underset{i,j=\{1,\dots ,n\}}{max}\left({\omega }_{i,j}\right)$$

(9)

The two PS from which the maximum tilt ${\omega }_{\max }$ is defined are designated as $P_{s1}$ and $P_{s2}$ with coordinates $(x_{s1},y_{s1})$ and $(x_{s2},y_{s2})$ and vertical displacements $z_{s1}$ and $z_{s2}$ , respectively. For $z_{s1}<z_{s2}$ , the numerical parameters of the maximum tilt plane function (equation (1)) are equal to

$$\begin{array}{c}b=\frac{(x_{s1}-x_{s2})(z_{s1}-z_{s2})}{{(y_{s1}-y_{s2})}^2+{(x_{s1}-x_{s2})}^2}\hfill \\ c=\frac{(y_{s1}-y_{s2})(z_{s1}-z_{s2})}{{(y_{s1}-y_{s2})}^2+{(x_{s1}-x_{s2})}^2}\hfill \end{array}$$

(10)

The numerical parameter $a$ of the maximum tilt plane equation affects the absolute but not the vertical relative displacement between PS on the plane and is thus disregarded for simplicity. The tilt surface and the PS from which it is derived are graphically illustrated in Figure 5.

Figure 5.

Maximum tilt surface and defining points ${\mathit{P}}_1$ and ${\mathit{P}}_2$ .

The entire area of Belgium cannot be covered by a single InSAR track, necessitating the use of multiple tracks which partially overlap. This overlap can lead to spurious tilt results from PS at a short distance with drastically different velocities which could not be resolved during the automatic processing. In the areas of track overlap, the above calculations are performed for the PS of each track individually. The most conservative result for each measurement period in terms of tilt between tracks is considered for the calculation of the maximum tilt.

Definition of building deformation

Polygon vertex and edge data

Each building polygon consists of a set of $m$ vertices $V_1,\dots ,V_m$ , each with a pair of coordinates $(x_v,y_v),v\in \{1,\dots ,m\}$ in the $R^2$ Euclidean space. In the present investigation, the number of vertices is four for all cases, but the methodology is applicable to any closed polygon. The vertical displacements $z_v$ of the vertices of a building polygon can be calculated by the substitution of its coordinates in the maximum tilt plane function (equation (1)). The vector of relative vertical displacements between consecutive vertices is calculated according to

$$\mathit{s}=(|z_1-z_2|,\dots ,|z_m-z_1|)$$

(11)

The vector of lengths of the $m$ polygon edges are calculated according to

$$\mathit{l}=({\Vert V_1-V_2\Vert }_2,\dots ,{\Vert V_m-V_1\Vert }_2)$$

(12)

The tilt ${\omega }_v$ of the each of the edges of the polygon can be calculated by the entry-wise division of the vertical displacements with the edge length

$$\mathit{\omega }=\mathit{s}ø\mathit{l}=(\frac{s_1}{l_1},\dots ,\frac{s_m}{l_m})$$

(13)

Polygon edge tilt

The calculation of the tilt for each individual edge allows the evaluation of damage at different locations of the structure. Edges parallel to the dip direction of the maximum tilt surface suffer the most damage, whereas edges perpendicular to it do not suffer damage. In the case of piecewise analysis defined by accumulated ground deformation from different measurement periods, this approach allows the evaluation of the damage in different parts of the building over time. Since this is normally the case with InSAR data, the analysis results from different periods constitute the damage history of the building, which can be compared with previous inspection records and patrimonial database information.

The absolute tilt of each edge is accumulated in every subsequent measurement period. Damage from previously registered tilt is considered not reversible. Practically, this means that when opposite sign tilts are registered at different periods, only the sum of their absolute values is considered for the accumulated tilt of the edge. The repeated unloading and reloading of building structures subjected to soil settlement have not been previously investigated experimentally, analytically or numerically. However, empirical damage assessment models are based not only on the magnitude of strain or crack width from monotonic settlement loading but also on parameters, such as crack distribution and extent. ¹³ Since opposite sign tilt may cause damage to different parts of the structures still retaining their stiffness, the conservative approach of the accumulation of absolute tilt is adopted here.

Definition of load parameters

For a number of settlement damage models for masonry, the accumulated tilt is not a sufficient parameter to determine the loading. Since the maximum tilt surface only provides information on the tilt of an edge, further manipulation is required for the indirect determination of other loading parameters.

In brief, the loading parameters are the angular distortion $\beta$ , the deflection ratio $\delta$ , the bending strain $\epsilon$ and the tilt $\omega$ . The tilt $\omega$ is defined as the rigid body rotation of the whole superstructure. The angular distortion $\beta$ is the rotation of the straight line joining two reference points relative to the tilt. The deflection ratio $\delta$ is defined as the quotient of relative deflection and the corresponding length. ¹⁶ Finally, the bending strain $\epsilon$ does not describe ground deformation but is indirectly calculated for the superstructure according to the other load parameters.

As shown above, the vertices of an edge with a horizontal distance $l$ and a relative vertical displacement $s$ define two points on a plane perpendicular to the geographic $(x,y)$ plane. To calculate the deflection ratio $\delta$ , a Gaussian curve is fitted to these two points. This fitting is intended to represent the distorted shape of a shell element resting on a plane characterized by linearly variable settlement. ¹⁵ This deformation profile is presented schematically in Figure 6. The equation of the Gaussian curve along the edge, with $x$ here corresponding to the local axis of the edge for brevity, reads

$$z_g\left(x\right)=s·\exp (-\frac{x^2}{2x_i^2}),x\in [0,\infty )$$

(14)

where $x_i$ is the inflection distance of the curve. A Gaussian curve may be roughly approximated by a cosine expression. One such expression is

$$z_c\left(x\right)=\frac{s}{2}(1+\cos \left(\frac{\pi x}{l}\right)),x\in [0,l]$$

(15)

While this approximation is not very accurate, it helps in defining the point beyond which equation (14), which has no root, is very close to zero. The root of equation (15) is $l$ . Defining the inflection point distance of the Gaussian curve to be at $x_i=l/\pi \cong 0.318·l$ , the difference between equations (14) and (15) at $x=l$ is reduced to less than $0.1\%$ of $s$ , which may be considered negligible. Therefore, the sagging length is equal to $l_s=x_i$ and the hogging length is equal to $l_h=l-x_i$ . The deflection at the inflection point is $z_g\left(x_i\right)\cong 0.6065·s$ and the deflection at $l$ is $z_g\left(l\right)\cong 0$ . The piecewise linear interpolation between the maximum deflection point, the inflection point and the zero deflection point is defined by the expression

$$z_i\left(x\right)=\{\begin{array}{cc}(z_g\left(x_i\right)-s)\frac{x}{x_i}+s\hfill & x\in [0,x_i]\hfill \\ \multicolumn{1}{c}{(z_g\left(l\right)-z_g\left(x_i\right))\frac{x-x_i}{l-x_i}+z_g\left(x_i\right)}& x\in [x_i,l]\hfill \end{array}$$

(16)

Figure 6.

Illustration of building edge deformation profile in blue based on vertex settlements. Single vertex settlement (left) and superposition of multi-edge settlement (right).

The three curves are illustrated in Figure 7. The sagging and hogging deflections are calculated as

$$\begin{array}{c}{\mathrm{\Delta }}_s\left(x\right)=|z_g\left(x\right)-z_i\left(x\right)|,x\in [0,x_i]\hfill \\ \multicolumn{1}{c}{{\mathrm{\Delta }}_h\left(x\right)=|z_g\left(x\right)-z_i\left(x\right)|,x\in [x_i,l]}\end{array}$$

(17)

and their maximum values are designated as $\max \left({\mathrm{\Delta }}_s\left(x\right)\right)={\mathrm{\Delta }}_s$ and $\max \left({\mathrm{\Delta }}_h\left(x\right)\right)={\mathrm{\Delta }}_h$ , respectively. The sagging and the hogging deflection ratios are finally calculated as

$$\begin{array}{c}{\delta }_s=\frac{{\mathrm{\Delta }}_s}{l_s}\cong 0.0809·\frac{s}{x_i}\hfill \\ \multicolumn{1}{c}{{\delta }_h=\frac{{\mathrm{\Delta }}_h}{l_h}\cong 0.1914·\frac{s}{(l-x_i)}}\end{array}$$

(18)

Figure 7.

Edge deformation profile and derivation of the deflection ratios. Cosine approximation ${\mathit{z}}_{\mathit{c}}\left(\mathit{x}\right)$ and piecewise linear interpolation ${\mathit{z}}_{\mathit{i}}\left(\mathit{x}\right)$ are also shown.

The location of maximum sagging is $x_s\cong 0.137·l$ , and the location of maximum hogging is $x_h\cong 0.629·l$ . It is now possible to calculate the deflection ratio. For the assumed deflected shape, the hogging deflection ratio is greater. Therefore, the deflection ratio for the edge is $\delta ={\delta }_h$ .

The angular distortion $\beta$ and the bending strain $\epsilon$ are calculated from the expressions ¹¹

$$\beta =3\delta \frac{1+4R\frac{h^2}{l^2}}{1+6R\frac{h^2}{l^2}}$$

(19)

and ¹³

$$\epsilon =\frac{3\frac{l}{h}\delta }{\frac{1}{4}{\left(\frac{l}{h}\right)}^2+\frac{3}{2}R}$$

(20)

where $h$ is the height of the edge and $R$ is the ratio of the Young’s modulus over the shear modulus

$$R=\frac{E}{G}=2(\nu +1)$$

(21)

with $\nu$ being the Poisson’s ratio of the structural element represented by the edge.

Damage classification models

Damage models

Damage assessment due to settlement is calculated according to models for three different main structural typologies: masonry, infilled frames and bare frames. In the literature, the majority of work on sensitivity of buildings to soil settlement has been carried out for masonry buildings. This is motivated by their high stiffness, low tensile and shear strength and low deformation ductility. A limited amount of work has been carried out for infilled frames, whose sensitivity against soil movement is mostly due to the masonry infills and not the main load-bearing frame. Bare frames have been studied even less due to their relatively low sensitivity but represent a not insignificant part of industrial building heritage.

In this study, six different types of analytical damage models (a, b, c, d, e and f) are considered. Each is provided with a qualitative damage level description and a numerical damage index $I_{\{a,b,c,d,e,f\}}$ . The subscripts $a$ and $b$ refer to the Boscardin and Cording angular distortion and bending strain limits, respectively, ¹¹ $c$ refers to the tilt limits by Fischer, ³⁰ $d$ to the deflection ratio model by Giardina et al. ¹⁹ and $e$ and $f$ to the tilt limits for infilled and bare frames.

A summary of load limits for the a, b and c masonry models considered is provided in Table 2.

Table 2.

Summary of a, b and c damage models for masonry.

Damage level	Boscardin and Cording 11		Fischer 30
	$\beta (\times {10}^{-3})$	$\epsilon (\%)$	$\omega (\times {10}^{-3})$
Negligible	0.0–1.1	0.000–0.050	0–1
Very slight	1.1–1.6	0.050–0.075	1–2
Slight	1.6–3.3	0.075–0.150	2–3
Moderate to severe	3.3–6.7	0.150–0.300	3–5
Severe to very severe	>6.7	>0.300	>5

Model d for masonry walls is based on the deflection ratio $\delta$ and a polynomial expression of a damage index. ¹⁹ This model is based on the allowable crack widths for different damage levels proposed by Burland et al. ³¹ These limits are indicated in Table 3. The polynomial expression of the damage index reads

$$I_d=b_1+b_2\delta +b_3{\delta }^2+b_4{\delta }^3+\sum _{i=1}^6{a}_i{\overline{x}}_i$$

(22)

Table 3.

Settlement-induced crack width limits according to Burland et al. (Model d).

Damage level	Crack width(mm)
Negligible	0.0–0.1
Very slight	0.1–1.0
Slight	1.0–5.0
Moderate	5.0–15.0
Severe	15.0–25.0
Very severe	>25.0

The numerical parameters $a_i$ and $b_i$ of the damage model are explained in detail by Giardina et al. ¹⁹ and the values of ${\overline{x}}_i$ are the functions of the material and geometric properties of the masonry wall. Therefore, this approach allows the analysis of masonry buildings on a case-by-case basis. When these parameters are not known, reference values may be used instead.

In the case of infilled and bare frames, the limits in tilt for each damage level are acquired according to the relative allowable limits for drift given by ASCE 41. ³² These drift limits represent the capacity of the building for deformation based on its typology. The resulting tilt limits are presented in Table 4. Infilled frames are considered to have the same tilt limits as masonry structures, and bare frames have these tilt limits doubled. In the case of the tilt limits for infilled frames, however, it should be noted that the damage level refers primarily to the infill and not the frame. This infill damage should normally be repairable or the damaged element replaceable.

Table 4.

Damage models for infilled (Model e) and bare frames (Model f).

Damage level	Infilled frame $\omega ({10}^{-3})$	Bare frame $\omega ({10}^{-3})$
Negligible	0–1	0–2
Very slight	1–2	2–4
Slight	2–3	4–6
Moderate to severe	3–5	6–10
Severe to very severe	>5	>10

Load limit normalization

To maintain a uniform output, the predicted building damage from all the above models is expressed on a four-tier damage scale: null/negligible (NN) for $I\in [0,1)$ , slight/light (SL) for $I\in [1,2)$ , moderate/severe (MS) for $I\in [2,3)$ and very severe (VS) for $I=3$ . Therefore, the above expressions for the damage indices must be normalized to the desired scale. This is done according to the following linear interpolation and scaling expressions

$$I_a\left(\beta \right)=\{\begin{array}{ccc}0.75(0+\frac{\beta -0.0000}{0.0011-0.0000})\hfill & \mathrm{if}\hfill & 0.0000\le \beta \le 0.0011\hfill \\ \multicolumn{1}{c}{0.75(1+\frac{\beta -0.0011}{0.0016-0.0011})}& \mathrm{if}\hfill & 0.0011\le \beta \le 0.0016\hfill \\ \multicolumn{1}{c}{0.75(2+\frac{{\beta }_j-0.0016}{0.0033-0.0016})}& \mathrm{if}\hfill & 0.0016\le \beta \le 0.0033\hfill \\ \multicolumn{1}{c}{0.75(3+\frac{\beta -0.0033}{0.0067-0.0033})}& \mathrm{if}\hfill & 0.0033\le \beta \le 0.0067\hfill \\ \multicolumn{1}{c}{3}& \mathrm{if}\hfill & 0.0067\le \beta \hfill \end{array}$$

(23)

$$I_b\left(\epsilon \right)=\{\begin{array}{ccc}0.75(0+\frac{\epsilon -0.00000}{0.00050-0.00000})\hfill & \mathrm{if}\hfill & 0.00000\le \epsilon \le 0.00050\hfill \\ \multicolumn{1}{c}{0.75(1+\frac{\epsilon -0.00050}{0.00075-0.00050})}& \mathrm{if}\hfill & 0.00050\le \epsilon \le 0.00075\hfill \\ \multicolumn{1}{c}{0.75(2+\frac{\epsilon -0.00075}{0.00150-0.00075})}& \mathrm{if}\hfill & 0.00075\le \epsilon \le 0.00150\hfill \\ \multicolumn{1}{c}{0.75(3+\frac{\epsilon -0.00150}{0.00300-0.00150})}& \mathrm{if}\hfill & 0.00150\le \epsilon \le 0.00300\hfill \\ \multicolumn{1}{c}{3}& \mathrm{if}\hfill & 0.00300\le \epsilon \hfill \end{array}$$

(24)

$$I_c\left(\omega \right)=\{\begin{array}{ccc}0.75(0+\frac{\omega -0.000}{0.001-0.000})\hfill & \mathrm{if}\hfill & 0.000\le \omega \le 0.001\hfill \\ \multicolumn{1}{c}{0.75(1+\frac{\omega -0.001}{0.002-0.001})}& \mathrm{if}\hfill & 0.001\le \omega \le 0.002\hfill \\ \multicolumn{1}{c}{0.75(2+\frac{\omega -0.002}{0.003-0.002})}& \mathrm{if}\hfill & 0.002\le \omega \le 0.003\hfill \\ \multicolumn{1}{c}{0.75(3+\frac{\omega -0.003}{0.005-0.003})}& \mathrm{if}\hfill & 0.003\le \omega \le 0.005\hfill \\ \multicolumn{1}{c}{3}& \mathrm{if}\hfill & 0.005\le \omega \hfill \end{array}$$

(25)

$$I_d\left(\delta \right)=\frac{\max (0,\min (b_1+b_2\delta +b_3{\delta }^2+b_4{\delta }^3+\sum _{i=1}^6{a}_i{\overline{x}}_i,6))}{2}$$

(26)

$$I_e\left(\omega \right)=\{\begin{array}{ccc}0.75(0+\frac{\omega -0.000}{0.001-0.000})\hfill & \mathrm{if}\hfill & 0.000\le \omega \le 0.001\hfill \\ \multicolumn{1}{c}{0.75(1+\frac{\omega -0.001}{0.002-0.001})}& \mathrm{if}\hfill & 0.001\le \omega \le 0.002\hfill \\ \multicolumn{1}{c}{0.75(2+\frac{\omega -0.002}{0.003-0.002})}& \mathrm{if}\hfill & 0.002\le \omega \le 0.003\hfill \\ \multicolumn{1}{c}{0.75(3+\frac{\omega -0.003}{0.005-0.003})}& \mathrm{if}\hfill & 0.003\le \omega \le 0.005\hfill \\ \multicolumn{1}{c}{3}& \mathrm{if}\hfill & 0.005\le \omega \hfill \end{array}$$

(27)

$$I_f\left(\omega \right)=\{\begin{array}{ccc}0.75(0+\frac{\omega -0.000}{0.002-0.000})\hfill & \mathrm{if}\hfill & 0.000\le \omega \le 0.002\hfill \\ \multicolumn{1}{c}{0.75(1+\frac{\omega -0.002}{0.004-0.002})}& \mathrm{if}\hfill & 0.002\le \omega \le 0.004\hfill \\ \multicolumn{1}{c}{0.75(2+\frac{\omega -0.004}{0.006-0.004})}& \mathrm{if}\hfill & 0.004\le \omega \le 0.006\hfill \\ \multicolumn{1}{c}{0.75(3+\frac{\omega -0.006}{0.010-0.006})}& \mathrm{if}\hfill & 0.006\le \omega \le 0.010\hfill \\ \multicolumn{1}{c}{3}& \mathrm{if}\hfill & 0.010\le \omega \hfill \end{array}$$

(28)

Application examples

Summary of calculation process

The entire process for potential damage calculation according to the developed methodology is illustrated in the flowchart of Figure 8. In summary, the patrimonial building data are processed for the calculation of the building polygons and their 10-m offset. Next, each InSAR data period $i$ is processed in a loop for the calculation of the PS for each period (PSi). The 10-m velocity grids are subsequently calculated through inverse distance weighing for each period (IDW_i). Next, for each building $j$ , a check is made on the number of PS within the area of the polygon plus its offset ${(\mathrm{PS}}_j^i)$ . If two or more PS are found, then the soil maximum tilt plane (during the period $i$ near building $j$ ( $f_i^j(x,y)$ from equation (1)) is calculated based on the PS data, otherwise the calculation resorts to the use of the grid data. The total accumulated tilt of the building $j$ edges at the end of the time period $i$ $({\omega }_i^j)$ is equal to the tilt increment calculated in the current period ( $\mathrm{\Delta }{\omega }_i^j$ from equation (9)) added to the total accumulated tilt from all previous periods $({\omega }_{i-1}^j)$ . Finally, the potential damage $I_i^j$ of building $j$ at the end of period $i$ can be calculated as a function of the accumulated tilt.

Figure 8.

Calculation flowchart for proposed method.

Verification of damage model

The results of the damage model are compared to the related results from the literature on settlement damage to buildings as assembled and summarized by Namazi and Mohamad. ¹² Settlement in these cases was induced by a variety of causes, including tunnelling, self-weight and adjacent excavation. The cited sources are used for the determination of the dimensions and type of the building as well as of the deformation loads in terms of the deflection ratio. The Poisson’s ratio for masonry buildings (Cases 1, 2, 4, 6, 7, 11–15) was assumed to be equal to 0.3, as was for a number of RC frame cases with masonry infills (8–10). Other RC frame buildings (Cases 3, 5, 16) were analyzed using a Poisson’s ratio of 5.25, which suggests the substantial shear deformability of a bare frame. The criteria proposed by Burland et al. ³¹ were used in the cited article for classifying damage to the buildings (Table 3). The results are presented in Table 5. For the RC frame cases, both the infilled and bare frame damage indices are presented.

Table 5.

Building case studies from the literature: comparison of reported with predicted damage.

No.	Reference	Type	$R$	Width(m)	Length(m)	Height(m)	δ(%)	Reported damage	$I_a\left(\beta \right)$	$I_b\left(\epsilon \right)$	$I_c\left(\omega \right)$	$I_d\left(\delta \right)$	$I_e\left(\omega \right)$	$I_f\left(\omega \right)$
1	Netzel 33	Masonry	2.6	8.0	13.0	5.0	0.015	Negligible	0.235	0.314	0.596	1.220	–	–
2	Skempton and MacDonald 34	Masonry	2.6	23.0	32.0	20.0	0.833	Severe	3.000	3.000	3.000	3.000	–	–
3	Skempton and MacDonald 34	RC frame	12.5	55.0	61.0	20.0	0.580	Severe	–	–	–	–	3.000	3.000
4	Breth and Chambosse 35	Masonry	2.6	12.0	22.0	10.0	0.400	Very severe	3.000	3.000	3.000	3.000	–	–
5	Standing 36	RC frame	12.5	18.0	80.0	43.4	0.010	Negligible	–	–	–	–	0.398	0.199
6	Bloodworth et al. 37	Masonry	2.6	4.4	7.6	6.0	0.015	Very slight	0.214	0.199	0.596	1.220	–	–
7	Anketell-Jones and Burland 38	Masonry	2.6	30.0	70.0	20.0	0.120	Slight	2.086	1.960	3.000	2.407	–	–
8	Finno et al. 39	RC frame	2.6	35.0	42.0	15.3	0.040	Slight	–	–	–	–	1.590	0.795
9	Macleod and Paul 40	RC frame	2.6	17.6	74.7	12.6	0.051	Slight	–	–	–	–	2.028	1.014
10	Macleod and Paul 40	RC frame	2.6	8.3	44.0	13.4	0.005	Slight	–	–	–	–	0.199	0.099
11	Macleod and Paul 40	Masonry	2.6	8.0	18.0	7.6	0.275	Moderate	2.896	3.000	3.000	2.981	–	–
12	Macleod and Paul 40	Masonry	2.6	8.0	24.5	7.6	0.081	Moderate	1.652	1.956	2.735	2.084	–	–
13	Macleod and Paul 40	Masonry	2.6	14.6	23.7	4.0	0.042	Slight	0.796	1.158	1.670	1.629	–	–
14	Li and Cameron 41	Masonry	2.6	8.0	16.0	3.0	0.164	Moderate	2.479	1.942	3.000	2.650	–	–
15	Li 42	Masonry	2.6	15.0	17.0	3.0	0.370	Severe/very severe	3.000	3.000	3.000	3.000	–	–
16	Fjeld 43	RC frame	12.5	25.0	34.2	13.5	0.400	Severe/very severe	–	–	–	–	3.000	3.000

Accurate predictions shaded in green, overestimation of damage in blue, underestimation of damage in red.

Overall, the numerical damage indices from the analysis show good agreement with the reported damage for the majority of the cases and for all structural typologies. The masonry models appear to overestimate the damage in Cases 7 and 11, although the magnitude of the deflection ratio in both cases suggests that the reported damage should normally be higher. The tilt model index $\left(I_c\right)$ overestimates the damage for high loads, and the deflection ratio model index $\left(I_d\right)$ overestimates it for low loads. The angular distortion $\left(I_a\right)$ and strain $\left(I_b\right)$ model indices produce similar values for most cases, intermediate in magnitude to the other models.

The cases where the potential damage was underestimated are here investigated in greater detail. The underestimation of the potential damage in Case 8 only occurs when considering a bare frame case $(I_f)$ , while the prediction is accurate when considering an infilled frame $(I_e)$ , as is the actual case. In Case 10, some of the reported damage (vertical cracks beneath windows) is not typical of soil movement and was attributed to shrinkage of the masonry. Finally, only a small underestimation is obtained in Case 12 when using the $I_a$ and $I_b$ models, while the $I_c$ model gives a more accurate, if slightly conservative, prediction. Overall, the $I_c$ model has not underestimated the damage in any of the cases and has overestimated it fewer times than the $I_d$ model.

An additional verification of the proposed model is made through a comparison of the damage curve described by equation (25), based on tilt, with an extensive field study, conducted by Peduto et al. ⁹ in the Netherlands on a large number of structures damaged by differential settlement. In the cited source, numerous buildings with shallow and piled foundations were inspected. For each building, the damage levels were observed, as defined by Burland et al., ³¹ and the rotation angle, which for the small angle values considered is nearly equal to the slope or tilt, was recorded. It was thus possible to establish mean rotation angles for each damage level. The comparison of the resulting damage curve, with damage levels quantified on a numerical scale, with the proposed tilt damage model is shown in Figure 9. No distinction was made between severe and very severe damage levels by Peduto et al., limiting the extent of the curve. Overall, the proposed tilt model is shown to be unconservative compared to the average findings of the field study in the Netherlands. However, the conservative way in which the maximum tilt for each building is defined in the proposed model (equation (9)) can potentially offset this difference.

Figure 9.

Comparison of tilt damage model with damage curve determined by Peduto et al.

Sensitivity analysis

The masonry damage models a, b, c and d are used in a sensitivity study, the main purpose being to quantify the differences in the predicted damage over a variety of cases. Both the loading data, determined by the position and vertical displacement of the PS, and the building dimensions and properties are generated according to specified constraints. The data are used for the calculation of the relevant damage indices. In total, 100,000 cases were generated and analyzed.

The parameters used in the analysis are presented in Table 6. The values for the building parameters are distributed uniformly between the minimum and maximum values provided and are representative of historic masonry buildings. The coordinates of the PS are generated randomly, and the values of the vertical displacement follow a skew normal distribution. The skewness is introduced to reflect the fact that in a given geographic area around a single building, the InSAR points will tend to be mostly uplifting or subsiding as a whole. The shape parameter used for this skewness is equal to 100 times the standard deviation. The values used for the distribution are in part consistent with the data obtained from the actual InSAR data from Belgium but are not intended to be wholly representative of them.

Table 6.

Values for sensitivity analysis.

Building parameter	Minimum value	Reference value	Maximum value
Width/length (m)	15	–	50
Height (m)	3	–	15
Openings percentage (%)	0	30	50
Tensile fracture energy (N/mm)	1	10	50
Tensile strength (N/mm2)	0.05	0.1	0.20
Young’s modulus (N/mm2)	1000	3000	6000
Foundation vertical stiffness (N/m3)	3.5e + 07	7e + 07	14e + 07
Shear behaviour of foundation	0	0	1
Ratio of Young’s and shear moduli	2.6	2.6	12.5
Loading	Mean value	Standard deviation
Vertical displacement of PS (m)	0.00	0.10

The results of this analysis are presented as fragility curves in terms of the cumulative percentage of cases for the damage indices of each model, illustrated in Figure 7. The analyses are performed in two sets: first, by assuming the reference values for the foundation stiffness, foundation shear behaviour and the superstructure $R$ parameter, and second, with a variation of these parameters as shown in Table 6. In the second approach, the openings percentage and the $R$ parameter are linearly linked due to the increase in the shear deformability when the openings percentage is increased. These parameters are linked according to

$$R=2.6+19.8\times \left[\mathrm{openings}\mathrm{percentage}\right]$$

(29)

resulting in a value for $R=2.6$ for no openings and $R=12.5$ for 50% openings. The former value corresponds to a Poisson’s ratio of masonry of $\nu =0.3$ , which is a typically accepted design value for the material. ⁴⁴ The latter value reflects the reduction of shear stiffness due to the presence of openings. ¹³ The variation of these parameters, which strongly affects the soil–structure interaction characteristics of the buildings, does not affect the results for $I_a\left(\beta \right)$ and $I_c\left(\omega \right)$ .

For the first approach (Figure 10(a)), the tilt model (c) is the most conservative of the models, with 12.0% of the cases having a damage index of more than 2. The deflection ratio model (d) is the least conservative of the four, with 1.6% of buildings having a damage index of 2 or more. The angular distortion (a) and bending strain (b) models are nearly equivalent and lie between the other two models, with 7.0% of cases having a damage index of 2 or more.

Figure 10.

Fragility curves for generated cases according to four masonry models. (a) Using reference values. (b) Using varying properties for foundation and superstructure stiffness.

In the second approach (Figure 10(b)), the tilt model (c) remains the most conservative. However, the bending strain model (b) approximates it for the entire range of results. The deflection ratio model (d) is more conservative than in the first approach throughout the range of results, particularly in the range of lower damage indices. The angular distortion model (a) becomes the most conservative for low-damage indices and roughly at the centre of the envelope of the predicted results generated from all the models for indices above 1. Overall, the four models in the second approach provide a much narrower results envelope of predicted damage compared to the first approach.

The differences between the two approaches reflect the importance of considering soil–structure interaction effects in differential settlement analysis for damage prediction in buildings. However, the superstructure characteristics that influence these effects are generally not found in patrimonial building databases and are thus not practical to implement in a country- or region-scale analysis. Therefore, the tilt model (c), which proved to be the most conservative for the reference value approach, is adopted for the presentation of the country-scale results for Belgium.

Application in Belgium

Damage calculation

The results of the damage analysis are summarized in Table 7. These results indicate a total of 3715 structures potentially having sustained differential settlement damage in the period 1992–2018. The majority of potentially damaged structures are located in the Flanders region, yet the percentages for each category are similar for the three regions.

Table 7.

Summary of potential damage calculation for Belgium in 1992–2018: number of potentially damaged buildings per damage level, summarized per federal region.

Damage	Brussels	%	Flanders	%	Wallonia	%	Total	%
Negligible	633	71.69	1658	70.98	370	74.60	2661	71.63
Slight	86	9.74	303	12.97	69	13.91	458	12.33
Light	45	5.10	123	5.27	19	3.83	187	5.03
Moderate	27	3.06	96	4.11	7	1.41	130	3.50
Severe	22	2.49	71	3.04	9	1.81	102	2.75
Very severe	70	7.93	85	3.64	22	4.44	177	4.76
Total	883	100.00	2336	100.00	496	100.00	3715	100.00

Proportion of buildings in each damage level over total number of potentially damaged buildings.

This potential damage, as a function of the ground movement profile gradient in the area of the building, highlights specific structures in need of further investigation.

Buildings on moving soil per time period

In addition to the quantitative damage assessment of the Belgian patrimonial building stock, a further analysis is carried out on all analyzed buildings. As estimated by Peduto et al., ⁴⁵ an area subsiding or uplifting at a rate greater than 2 mm/year may be considered to be moving at a noticeable pace. While this criterion cannot be used for damage assessment, as it does not take differential settlements within the building into account, it provides an indication of possible regions with the risk of settlement damage.

The regions exceeding the velocity threshold of 2 mm/year are identified by the velocity values of the interpolation grid. Building polygons and their offset in contact with interpolation cells having a velocity greater than 2 mm/year are considered to be located on moving soil. The results of this analysis are summarized in Table 8. The period 1992–1997 indicates a total of 21,220 buildings on moving soil, corresponding to 7.9% of the patrimonial building stock. The overall number of heritage buildings on moving soil tends to decrease in later periods, with a small increase noted in Flanders and Wallonia for 2005–2010 and an increase in Brussels during 2016–2018.

Table 8.

Number of buildings on moving soil according to the 2 mm/year threshold, summarized per federal region and time period.

Period	Brussels	%	Flanders	%	Wallonia	%	Total	%
1992–1997	5159	10.46	13,292	7.69	2769	9.28	21,220	8.42
1996–2001	156	0.32	9904	5.73	326	1.09	10,230	4.06
2001–2006	33	0.07	1157	0.67	208	0.70	1398	0.55
2005–2010	47	0.10	2291	1.32	717	2.40	3055	1.21
2016–2018	123	0.25	1667	0.96	235	0.79	2025	0.80

Proportion of buildings on moving soil over total number of buildings.

Results verification and evaluation

The analysis results were compared to site-visit findings from selected case studies. Two methods were employed for the selection of the case studies: (a) selection based on high registered tilt values during the measurement periods and (b) selection according to the analysis result, that is, buildings with high potential damage. Additional attention was paid to factors, such as accessibility. Overall, 18 buildings were investigated, mostly located in the areas around the closed collieries of Limburg and in the city of Leuven. The cases, analysis results and inspection outcomes are summarized in Table 9.

Table 9.

Potential damage calculation in selected cases.

No.	Case type	City/Town	Building type	$I(\omega )$	Inspection	Calculation
1	Residential building block	Hasselt	Masonry	0.19	Negligible	IDW / PS / IDW
2	Large school building	Hasselt	Masonry	0.12	Negligible	IDW / IDW / IDW
3	Courthouse	Hasselt	Masonry	0.32	Slight	IDW / IDW / PS
4	Public building (government)	Hasselt	Masonry	0.65	Slight	PS / PS / IDW
5	Town house	Hasselt	Masonry	0.07	Slight	IDW / IDW / IDW
6	Church and monastery	Hasselt	Masonry	0.38	Slight	IDW / PS / PS
7	Hospital and monastery	Hasselt	Masonry	0.24	Negligible	IDW / PS / PS
8	Industrial heritage (mining concession)	Genk	Masonry + RC	0.15	Moderate	PS / PS / PS
9	Public swimming pool	Genk	RC	0.08	Negligible	PS / IDW / PS
10	Church	Oudsbergen	Masonry	0.21	Light	PS / IDW / IDW
11	Industrial heritage (mining concession)	Beringen	RC	0.17	Negligible	PS / PS / PS
12	Industrial heritage (mining concession)	Beringen	RC	0.31	Negligible	IDW / PS / PS
13	Church with cemetery	Herk-de-Stad	Masonry	0.05	Negligible	IDW / IDW / IDW
14	Town house	Leuven	Masonry	3.00	Severe	IDW / PS / IDW
15	Workers’ housing	Leuven	Masonry	3.00	Severe	PS / IDW / IDW
16	Coupled town houses	Leuven	Masonry	2.48	Moderate	IDW / IDW / PS
17	Farmhouse	Peer	Masonry	0.01	Negligible	IDW / IDW / IDW
18	Town house	Hasselt	Masonry	2.56	Moderate	IDW / PS / IDW

The calculation column indicates the method used for calculating the tilt in the three measurement periods: the interpolated grid (IDW) or the persistent scatterers (PS).

In total, 12 buildings were selected according to Method (a). According to the analysis results, all the buildings suffered negligible damage during the period 1992–2018. For 50% of the cases, on-site investigation did not reveal any differential movement damage and is thus in agreement with the analysis results. The remaining 50% presented damage, ranging from slight to moderate, but with evidence indicating that this damage was induced and repaired before the beginning of the measurement period.

Two buildings were highlighted according to Method (b) in the city of Hasselt, one with light and one with moderate calculated potential damage. Both cases have undergone renovation interventions during the measurement period. It was not possible to determine whether the interventions were related to settlement damage.

Three buildings were identified according to Method (b) in the city of Leuven, one with moderate and two with severe calculated potential damage. The moderately damaged and one of the severely damaged buildings have undergone renovation interventions. The other severely damaged building indeed presents substantial differential movement damage: wall cracks typical of ground settlement, roof parapet collapse and lintel failure.

Overall, distinguishing between movement partially registered due to structural changes, such as alterations in the roof, and movement due to ground settlement and uplift is not straightforward. The effects of the former are smeared out in the 10-m interpolation grid and thus do not affect the velocity threshold criterion. However, they may affect the potential damage calculation in the absence of prominent PS caused by actual ground movement.

It is possible to distinguish between movement in PS on the ground and in PS on buildings simply by evaluation of their coordinates against nearby building polygons. The movement of the former PS can be in most cases assumed to be caused by ground movement, although pavement renovation and topsoil changes cannot be ruled out with certainty. In the latter case, without the use of purpose-built reflectors installed on rigid points in selected buildings, ⁴⁶ it is impossible to absolutely distinguish between movement caused by ground movement and by structural alterations. Wide-scale assessment is not conducive to this targeted strategy, as it would require a prohibitively large number of reflectors. Therefore, the adoption of more advanced interpolation methods, coupled with InSAR processing methods capable of producing a higher density of PS, presents a more realistic strategy. However, this approach would require piecewise application in geographic regions of sufficiently small dimensions to manage computational costs.

Conclusion

A method for the calculation of settlement-induced damage based on InSAR monitoring data on the country scale is presented in this article. This monitoring data are used in the combination with processed patrimonial and cadastral data for the calculation of imposed settlements on buildings, focusing, but not being limited to, historic masonry structures.

The damage calculation method was tested against a number of case studies from the literature, demonstrating sufficient accuracy in terms of damage classification. The verification against these case studies shows the potential of the calculation method to be in principle applicable at any analysis scale. The different damage indexing methods provide a range of results for common structural typologies.

The analysis results for the entirety of Belgium are presented, in terms of potential damage suffered in the period 1992–2018 and in terms of the amount of buildings on subsiding or uplifting soil. The former criterion, derived from the ground movement gradient in the area of the building, can be used for identifying potential building in need of repairs. The latter criterion, based on the ground movement velocity in the area of the building, can be used for prioritizing large-scale inspection efforts in areas of interest.

While the proposed method shows sufficient flexibility in terms of defining the properties of the building and loading parameters, the matter of soil–structure interaction during the application of arbitrary soil movement profiles remains open. The next step in the refinement of the methodology developed here would be the incorporation of this aspect in a generalized way. The marked influence of the soil type and the shear stiffness of the superstructure on the calculated damage index as highlighted in the sensitivity analysis further motivates this investigation.

Practical limitations currently restrict the use of InSAR data to medium-resolution sets. As a subject of future work, a two-step analysis can be designed, in which a wide-area assessment of ground movement is conducted using medium-resolution data, allowing the identification of areas of high interest. These areas can, in turn, be analyzed locally using high-resolution data, such as those provided by the TerraSAR-X or Cosmo-SkyMed satellites.

Finally, a coupling of damage analysis results with wider site-inspection efforts is necessary and needs to be emphasized in future work. This coupling serves the further validation of this and other similar analysis methods for ground movement damage prediction over large areas and can assist in the establishment of new and the enrichment of existing empirical relations between ground movement measured in the country scale and damage to buildings.