Impact of Urbanization on Seismic Risk: A Study Based on Remote Sensing Data

: The management of seismic risk is an important aspect of social development. However, urbanization has led to an increase in disaster-bearing bodies, making it more difﬁcult to reduce seismic risk. To understand the changes in seismic risk associated with urbanization and then adjust the risk management strategy, remote-sensing technology is necessary. By identifying the types of earthquake-bearing bodies, it is possible to estimate the seismic risk and then determine the changes. For this purpose, this study proposes a set of algorithms that combine deep-learning models with object-oriented image classiﬁcation and extract building information using multisource remote sensing data. Following this, the area of the building is estimated, the vulnerability is determined, and, lastly, the economic and social impacts of an earthquake are determined based on the corresponding ground motion level and fragility function. Our study contributes to the understanding of changes in seismic risk caused by urbanization processes and offers a practical reference for updating seismic risk management, as well as a methodological framework to evaluate the effectiveness of seismic policies. Experimental results indicate that the proposed model is capable of effectively capturing buildings’ information. Through veriﬁcation, the overall accuracy of the classiﬁcation of vulnerability types reaches 86.77%. Furthermore, this study calculates social and economic losses of the core area of Tianjin Baodi District in 2011, 2012, 2014, 2016, 2018, 2020, and 2021, obtaining changes in seismic risk in the study area. The result shows that for rare earthquakes at night, although the death rate decreased from 2.29% to 0.66%, the possible death toll seems unchanged, due to the increase in population.


Introduction
Urbanization refers to the transformation of rural populations into urban populations, the migration of rural people into cities and people no longer working in agriculture [1]. Towns and cities are formed and increase in size with this process. In recent decades, industrialization and modernization have accelerated the process of urbanization, and, as a result, the proportion of the urban population has been increasing worldwide, notably in China. China's urbanization rate has increased steadily over the last decades. In 1950, 13% of people in China lived in cities. By 2010, the urban share of the population had grown to 45% [2]. The Seventh Population Census of China, conducted in 2020, showed that about 63.9 % of the total population lived in cities in 2020.
However, as urbanization continues, the accumulation of the urban population and wealth will directly increase the risk of disaster and pose challenges for disaster mitigation [3][4][5].
is commonly used in remote sensing classification [47][48][49][50].Overall, although the research on urbanization process analysis [51][52][53][54][55][56][57][58] and seismic risk assessment [17][18][19][20][21][22][23][24][25][26][27][28][29] using remote sensing data has achieved many results, the specific research does not address the changes in seismic risk caused by urbanization and the connection between the two. Accordingly, there is no established methodology within the field of remote sensing applications for understanding the seismic risk associated with urbanization. Therefore, in this paper, a practical method of observing the changes in seismic risk under urbanization based on remote sensing data is proposed. The objectives of this study are as follows: (i) To improve the remote sensing data analysis method for earthquake bearing-body detection by integrating deep learning semantic segmentation and ensemble learning classification. (ii) To propose a comprehensive workflow for identifying seismic risk change under urbanization processes using remote sensing data. The remainder of this paper is organized as follows: Section 2 provides an overview of the study area and its materials. A method for identifying structural vulnerabilities that integrates object-oriented classification and deep-learning-based segmentation is described as well. Section 3 presents the experimental results in the study area. Section 4 discusses improvements and future directions. A comprehensive summary is given in Section 5.  Tianjin is in the central part of the North China Plain. The city is a port city and the only megacity in China to have experienced a major earthquake in recent decades. Tianjin's District Baodi is a short drive from Beijing, Tangshan City, and the Tianjin downtown area [59]. The agriculture, industry, and tourism here have flourished over the past few decades, and the district has grown from a village to a county, then to a district [60]. However, the district is at risk of earthquakes. Throughout its history, District Baodi has been affected by many major earthquakes, including those of magnitude 7.8 in 1976 Tangshan and 8 in 1679 Sanhe-Pinggu [8]. These earthquakes caused extensive damage to Baodi, and the earthquake intensity was as high as VIII CSIS (China seismic intensity scale) in most areas [61]. An east-west fault in the region, the Baodi fault, shows evidence of activity in the Quaternary [62]. Over the last ten years, District Baodi was affected by three earthquakes that measured magnitude 3.7, magnitude 4.0, and magnitude 3. 3, on  14 January 2005, 18 June 2012, and 26 August 2012, respectively [8].

Data Sources
In this study, very high resolution (VHR) images, as well as medium-resolution images were collected. They are mainly from Gaofen-2 and Sentinel-2 sensors. ALOS-1 and WorldView-1 satellite data and Google Earth imagery were used as supplements in the years 2011, 2012, and 2014, as Gaofen-2 and Sentinel-2 satellites were not available over these years. At the same time, the census data and statistical yearbook data of the corresponding years were also collected, and the WorldPop [63] data were used as a reference for the spatial distribution of the population. Table 1 shows the information from the data source. In addition, four separate datasets were built for the four main tasks in this study: (i) footprint segmentation of single buildings, (ii) shadow segmentation of single buildings, (iii) rural building groups, and (iv) vulnerability classification of single buildings. Three of the datasets were instance segmentation datasets, and one was a multi-feature classification dataset.
Based on VHR satellite imagery and ground surveys located approximately 50 km from the study area, we produced data in shapefile format for 62,185 buildings. As shown in Figure 2, the footprint of a building is highly detailed. Building property information is given in the form of vulnerability type, usage, and floor numbers. Using ArcGIS Pro, we turned these data into a dataset that can be used as a basis for training building instance segmentation models, as well as for building vulnerability classification models. Additionally, 720 samples of rural building groups were derived from Sentinel 2 data, and shadow contours of 3250 buildings from GF2. Figure 3 illustrates an example of the instance segmentation dataset. To verify the final results, this study conducted a field survey in the study area during 2019-2020, collected a sample of 823 buildings, as shown in Figure 4, and created a sample of building structures by utilizing Baidu Maps API libraries.

Overall Workflow
In this section, the proposed workflow is described in detail. The seismic risk assessment process involves the quantification of three major input components, namely, seismic hazard intensity, exposure, and vulnerability [15].
Our understanding of the change in seismic risk resulting from urbanization relies on the assumption that the seismic hazard remains relatively stable; thus, the change in seismic risk is due to changes in the hazard-bearing body. When a disaster occurs involving an earthquake, social and economic losses are determined by the destruction of buildings. This form of structural change is the focus of our study. The primary role of remote sensing data in this study is to extract building information from images taken over several years. Objectoriented classification and deep-learning-based instance segmentation are integrated into the pipeline to efficiently accomplish this task.
As shown in Figure 5, the overall workflow includes four main parts. Part 1: Building object segmentation and feature extraction. The footprints of single buildings and rural building groups are extracted from high-and medium-resolution imagery, respectively. The image feature extraction is carried out with the building's footprint as the object unit. Part 2: Calculating the proxies in each object unit, according to the extracted object features, and then conducting vulnerability classification to obtain the disaster-bearing body dataset. Part 3: Calculating structural losses and the resulting economic and population losses at three ground motion intensity levels. Finally, repeat the above work for different years to obtain the results regarding changes in seismic risk during the urbanization process. In addition, to extract the building footprint and features and the vulnerability classification of the structure, multiple machine learning classifiers need to be pre-trained. This part can be regarded as Part 0 of the whole process. Figure 5. The overall workflow consists of four major sections: part 0 to pre-train the segmentation and classification models, part 1 to extract image features, part 2 to classify structural vulnerability, and part 3, a seismic risk assessment.

Building Feature Extraction
The acquisition of building parameters is probably the most time-consuming, tedious, and expensive part of each seismic risk assessment [64]. In this study, three BMask R-CNN classifiers were trained to extract footprints from single buildings and rural building groups, respectively. The footprints of single buildings are taken as object units for multi-feature extraction. The footprint of rural building groups can be directly applied to the classification of vulnerability and the calculation of inventory, according to the survey data.
The third BMask R-CNN classifier extracts the shadows of single buildings to calculate the height and number of floors. Other information on single buildings is extracted by eCognition v9.3, Trimble, CA, USA. A pretrained random forest classifier is used for structure vulnerability classification.

Mask R-CNN Framework
Mask R-CNN is a flexible object instance segmentation framework that efficiently detects objects in an image while simultaneously generating a high-quality segmentation mask for each instance [46]. Based on Faster R-CNN [65] and Fast R-CNN [66], Mask R-CNN adds a branch to predict an object mask while preserving the branch for bounding box recognition, thereby achieving pixel-level instance segmentation. Since Mask R-CNN is easy to generalize to other tasks, it has been widely adopted in remote sensing object classification. As shown in Figure 6, the original network structure of MASK R-CNN includes several components. The multiscale feature maps are extracted from the input image through the backbone part based on ResNet and the feature pyramid network (FPN). These features are shared by the RPN part and the RoIAlign layer. The feature map fed to the RPN is further extracted to generate candidate ROIs. After filtering, the obtained feature maps are used as proposals. The feature maps from the backbone part and the RPN part are properly aligned with the input based on bilinear interpolation [46] through the RoIAlign layer. Finally, the aligned feature map enters two branches in the head part: one is a fully convolutional mask prediction branch, and the other branch is divided into two sub-branches for class prediction and bounding box regression [67].

BMask R-CNN Framework
When performing pixel-level instance segmentation based on Mask R-CNN, predictions are made based on the local information. Although large receptive fields are obtained through the deep framework, which helps to extract features and improve the accuracy of classification, the information details, such as the shape information of the object, are still elusive.
To solve the problem of coarseness and indistinctness in the prediction output, Cheng, T.H. et al. [68] proposed a boundary-preserving Mask R-CNN to exploit boundary information and guide more precise mask prediction.
By adopting boundary features and boundary prediction, BMask R-CNN optimizes the mask head in Mask R-CNN, as illustrated in Figure 7. The new mask head is called the boundary-preserving mask head. Boundary-preserving mask heads synchronously learn object boundaries and masks. First, features from the mask sub-network can provide high-level semantic information for learning boundaries. Then, after obtaining the boundaries, the shape information and abundant location information in boundary features can help to achieve more precise mask predictions [68].
Since a boundary learning head branch is added to the multiple Mask R-CNN tasks, the loss function of the model also needs to accordingly increase a component. Here, ref. [68] proposes a combination of dice loss [69] and binary cross-entropy to optimize the boundary learning is binary cross-entropy loss, with λ as a hyperparameter to adjust the weight. p b and y b representing the predicted boundary for a particular category and the corresponding boundary ground truth, respectively. The Dice coefficient is used to measure the spatial overlap or similarities between the two sets. Here, the consistency between the predicted boundary and the corresponding boundary ground truth is compared. Since Dice loss is insensitive to the number of foreground/background pixels, it alleviates the class-imbalance problem. The calculation formula of Dice loss is as follows: where H and W are the height and width of the predicted boundary map, respectively; i denotes the i-th pixel, and is a smooth term to avoid zero division. Finally, after adding a boundary-preserving branch to Mask R-CNN, the combined multi-task learning loss functions are as follows: where L cls , L box , and L mask represent the loss of classification, localization, and segmentation mask, respectively, which are identical to those in [46].

Post-Processing of Building Footprint
After obtaining the footprints of single buildings and rural building groups from high-resolution and medium-resolution images, respectively, the post-processing of both results must be performed. This mainly includes the following processes: (i) Eliminating non-structural misclassification by setting an area threshold; (ii) Intersecting single buildings and building groups, keeping single buildings in the rural building group, and eliminating redundancy in the two output footprints; (iii) Calculating the actual building area of the rural building group.
In a typical rural residential setting in the study area, a class of simple structures with intact roofs often causes buildings to be misclassified. A good example of this is shown in Figure 8, where the roofs indicated in the red frame generally correspond to actual buildings, whereas other roofs may be simple structures or serve only as shelter from the sun and rain. To estimate the total amount of rural buildings, we apply the empirical formula based on the ratio of the total land area of the rural buildings group to the total area of the rural buildings. For our study area, the relationship between the group area of rural buildings and the number of households is as follows: group area = 652.77 m 2 × number of households + 1.4366 m 2 (4) Further, according to the average floor area of each household, the building stock to be counted can be obtained.

Estimating Floor Numbers
The height of a building is an effective way of assessing its seismic capacity and is essential to the calculation of its area [16]. This can be extracted by applying the shadow length of the building structure in the high-resolution optical image [18][19][20]27], SAR image imaging geometric characteristics [18], LiDAR data [27], or DSM data [25][26][27]. Since highresolution LiDAR data and DSM data were not obtained in this study, the number of floors in the building was inferred from the building shadows in the GF-2 data.
As illustrated in Figure 9, taking a regular building model as an example, according to the angle of solar irradiation and the angle of satellite observation, the geometric relationship between the building and its own shadow mainly presents two situations [70]. When the sun and the satellite are located on the same side of the building, the building itself partially occludes its own shadow. When the satellite and the sun are located on both sides of the building, the shadow of the building can be fully exposed to the viewing direction. According to the basic trigonometric function principle, the formula to calculate the height of the building in two different cases can be obtained [70]: here, H represents the height of the building, L 1 represents the length of the unobstructed shadow, and ω is the sun elevation angle. Considering the situation where the shadow is occluded when the sun and the satellite are on the same side, it is necessary to infer H according to L 2 and L 3 , as follows: where θ is the altitude angle of the satellite, α is the azimuth of the satellite, β is the azimuth of the sun, and ϕ is the azimuth of the building. The computational difficulty is significantly simpler when the satellites and the sun are on different sides of the building.

Occupancy and Population Disaggregation
We assume that the change in the POI of functional facilities is stable and that the increase or decrease in POI in each year is purely dependent on the existence of buildings. Based on the POI obtained in 2018 and 2020, the POI information has been assigned to the building through spatial analysis. By 2020, there were 110 medical institutions of various types, 63 educational institutions at all levels, 24 shopping malls, 39 enterprises exceeding their designated sizes, and 238 manufacturing companies found in the research area.
Samples from the on-site investigation were used to count the density of people in various buildings. These were divided into daytime and nighttime counts. Table 2 represents the Indoor Population Density of different occupancy (people per square meter).

Vulnerability Classification
Based on footprints of single buildings, object feature extraction was carried out in eCognition [71]. Then, a pretrained Random Forest (RF) [72] classifier was adopted to classify the vulnerability type of buildings. The selected features are shown in Table 3.

Loss Assessment
In this study, approaches to loss assessment are based on structural damage. Moreover, this study aims to determine the number of deaths and direct losses due to structural damage caused by ground motions. Explicitly addressing the damage and loss caused by secondary disasters such as surface fault rupture, landslides, soil liquefaction, fire, etc., as well as damage to infrastructures such as bridges and roads, is outside the scope of this paper.

Structural Damage
Structural damage assessment based on vulnerability analysis is the basis for quantifying economic loss and casualty. The vulnerability of structures that are exposed to earthquake loading expresses the likelihood of the occurrence of certain damage levels caused by seismic action [15].
Furthermore, the fragility model can be assumed to be a reliable measurement of damage to a respective set of buildings with similar structural taxonomy of dynamic behavior.
Building damage states are divided into five levels, which are intact, slightly damaged, moderately damaged, severely damaged, or collapsed, according to the damage to structural members or the entire structure. The division points of five kinds of failures correspond to the four limit state divisions of the structure in turn, and those from LS1 to LS4 gradually become more serious. The fragility model relates building response to seismic demand inputs and follows the lognormal assumption [73], as given by Equation (8), which refers to the conditional probability of various limit states of the structure under different earthquakes.
F (x, µ, σ) = P(LS|x) = Φ ((ln (x/µ))/σ) (8) here, Φ is a standard normal cumulative distribution function and x denotes seismic motion intensity, which is taken as peak ground acceleration (PGA). Parameters µ and σ are the fragility model median and standard deviation of ln (x). Table 4 presents the parameters of the fragility model used for different structural typologies in this study. According to the parameters of LSs (s = 1, 2, 3, 4), Equation (8) calculates the probability that the structure will reach LSs. Then, the probability of each DSi (i = 0, 1, 2, 3, 4) of the seismic intensity can be calculated: Based on the hazard level of the study area, we calculated the loss at three different levels of ground motion intensity. The three ground motion levels correspond to rare earthquakes, moderate earthquakes, and frequent earthquakes, respectively. The probability of exceedance during 50 years is 2-3%, 10%, and 63%, respectively.

Economic Loss
In this paper, the method used to calculate direct economic loss according to the damage state, loss ratio, and replacement price refers to the provisions of China Code GB/T 18208.4-2011 (Seismic Field Work Part IV: Disaster Direct Loss Assessment) [77]. According to this code, the formula proposed in this study for the calculation of direct earthquake economic loss for H types of structures is as follows: (1) The direct economic loss L A of H building structure types and D damage levels in a certain area is calculated as follows: A d denotes the total area of the h-type structure with damage stated; R d is the loss ratio when the h-type structure has a damage state of d; and P h is the replacement price of the h-type structure.
(2) The direct economic loss L B of the indoor property is calculated as: where T d is the ratio of the indoor property loss when the damage level of the h-type structure is d; P h is the replacement unit price of the h-type structure; and µ 1 is the ratio of the indoor property value of the building structure to the structural replacement price.
(3) The direct economic loss L C of decoration damage is calculated as follows: S d denotes the decoration loss ratio when the damage level of the h-type structure is d; Q h denotes the decoration price of the h-type structure; γ 1 represents the correction factor, considering the difference in economic development levels; γ 2 represents the correction factor considering building occupancy; and γ 3 represents the proportion of high-level decoration. γ 1 and γ 2 values are specified in [77].

Social Loss
In this study, death and injury calculations were carried out according to the relationship between the damage state of the house and the casualty rate of people, without distinguishing between structure types; the calculation formula is as follows: N d and N I denote the numbers of dead people and injured people, RD d and RI d represent the death and injury rate of people under different damage levels, A d is the area of the building structure under the damage level, and ρ is the density of people in the room. The fatality rate values used in this study are listed in Table 5.

Evaluation Indicators
The performance of the proposed algorithms was assessed using six metrics, namely, IoU, precision, recall, F1-score for segmentation tasks, and Overall Accuracy and Kappa coefficients for vulnerability classification.
Overall Accuracy = TP + TN TP + TN + FP + FN (23) where TP is the value of the true positives, FP is the value of the false positives, TN is the value of the true negatives, and FN is the value of the false negatives. P 0 is the relative observed agreement among raters. P e is the hypothetical probability of chance agreement [78].

Results
In this section, the results obtained by applying the workflow and method mentioned in Section 3 to the study area in Section 2 are presented. The results of the extraction of disaster-bearing body information and the estimation of seismic risk will be described, as well as the evaluation of the outputs.

Building Information Extraction Result
The building segmentation of the study area data is performed by the BMask-RCNN model presented in Section 3. Table 6 shows the extraction accuracy of single buildings and rural building groups. In terms of IoU, precision, recall, and F1 score, rural group buildings exhibit a lower extraction accuracy than single buildings. Extraction examples of single buildings and rural group buildings are presented in Figures 10 and 11. The results show a small number of errors and missing records.  The accuracy of height and floor estimation is shown in Figure 12a,b. An excellent linear fitting relationship exists between the height estimate and the real value, and the error of the height estimate is within one meter. The difference between the estimated number of floors and the actual number of floors may be as large as three stories, and the error for from 5-to 10-story buildings is larger.
According to field sampling data, the Overall Accuracy for the building vulnerability classification reached 86.77%, and the Kappa coefficient was 0.6538. Figure 13 shows the distribution of building structures over different years. In Figure 13, the brick-wood structures represented by the color red are gradually disappearing in the study area. In the extraction results for 2011, 2012, and 2014, brick-wood structures occupied the largest portion of the study area. According to the extraction results, only a few brick-wood structures survived in 2016.  Nevertheless, the overall pattern at the center of the study area has not changed dramatically, and the road between building blocks remains generally unchanged. This area of buildings remains unchanged and includes many functional facilities, including schools, hospitals, malls, offices, etc. Furthermore, based on the above results, the construction area of various structures in different years can be estimated. In 2020, the construction area of shear wall structure buildings reached 15 times that of the area in 2011, from 0.842 km 2 to 12.643 square km 2 . The proportion changed from 9.22% to 64.66%. Figure 14 illustrates how the area of each structure type changes over the years of this study. From 2011 to 2014, the construction area of various types of buildings remained stable, while the area for shear wall structures steadily increased. The construction area of brick wood buildings dramatically declined in 2016. However, the shear wall structure continued to grow into 2018 and saw a surge in 2020.
Grids are used to estimate the density of the building area.    Table 7 shows the estimated total population for each year.

Estimation of Seismic Risk Changes
Based on the methods outlined in Section 3 and the results presented in Section 3.1, estimates of losses under different ground motion levels, including direct losses in economic terms and deaths due to seismic activity during night and daytime, were further obtained for the study area. Table 8 provides a summary of the losses.  Figure 16 illustrates the trend in earthquake losses over many years based on different levels of ground motion intensity. The estimated results of losses include death tolls when earthquakes occur during the daytime as well as at night, as well as a death toll per 10,000 people when earthquakes occur at night. A direct economic loss is also incurred. In general, for frequent earthquakes, there was little change in losses from 2011 to 2021, except for an increase in direct economic losses after 2018. In terms of loss estimates from moderate earthquakes, as well as those from rare earthquakes, the trends were almost identical in recent years.
Firstly, it should be noted that the number of deaths caused by earthquakes occurs during the day and night, with one obvious trend being that the results in 2016 decreased compared with those before 2014 and began increasing thereafter. Additionally, after 2016, the death toll when the earthquake occurred at night was much lower than the death toll when it occurred during the day. A rare earthquake occurring at night will cause 1300 fewer deaths in 2021 than it did in 2011, which represents a reduction of one-third. Second, the death rate has declined from 2011 to 2021. For example, the death rate per 10,000 people resulting from a rare earthquake declined from 229 in 2011 to 66 in 2021. It can be said that the death rate decreased from 2.29% to 0.66%. Similarly, the death rate for every 10,000 people affected by the moderate earthquake declined from 93 to 27.
Third, the economic losses showed a slight upward trend between 2011 and 2014 and began to decline in 2016. There was an inflection point in 2018, and economic losses began to rise after that year. Figure 17 demonstrates the death toll of a rare earthquake that occurred at night. According to the results of the death toll estimation, the spatialized results for 2011 to 2014 indicate that this stage is mainly characterized by lower values uniformly distributed within space. In the 2016 and 2018 results, some low-value areas disappeared, while the original high-value areas remained. Additionally, the results for 2020 and 2021 indicate several increasing potential deaths.

Discussion
In this paper, we describe a new task regarding the perception of changes in seismic risk due to urbanization, based on remote sensing data. By applying multi-source remote sensing image data for different years and combining auxiliary information, we were able to monitor the earthquake-bearing body changes within the study area as the urbanization process continued and subsequently estimate seismic risk changes. Our work is based on the assumption that the seismic hazard remains relatively stationary over ten years.
Building structure information was extracted using the improved Mask R-CNN instance segmentation model, and the random forest classification algorithm was applied after instance segmentation to obtain a classification result. Furthermore, a method of earthquake disaster loss assessment was used to calculate the societal and economic losses over different years based on the three earthquake intensity levels.
First, building object types were divided into single buildings and rural building groups. Based on the BMask R-CNN model, we obtained relatively reliable results for both building object types.
Several factors affect the extraction accuracy of single buildings: The shadows cast by tall buildings and tall plants block the view of the buildings; these are the major reasons for misclassifications. Additionally, some low buildings that blend into the background environment are not identified. A third reason is that we use annotated data with very fine edges and an edge-preserving model to construct a finer outline of the building geometry. However, since the sample data have such a large number of edges, it is still challenging to achieve the same level of building segmentation.
Factors that influence building segmentation accuracy in rural areas include: In some cases, the plastic sheds surrounding rural buildings are similar in color, tone, and geometric size to the buildings, which confuses their classification. As in single buildings, some backgrounds, such as barren land, cause misclassifications and omissions. Moreover, some independent buildings are also incorrectly classified as rural group buildings.
A redundant processing operation is carried out at the intersection between rural group buildings and single buildings, based on the results of the single building extraction.
Although the results for estimations of the height of the building are relatively accurate, there is a certain amount of variance in the estimation of the number of floors of the building. This can be attributed to several factors: (i) There is a wide range of story heights in factory buildings with a single floor, ranging from 3 m to 10 m. (ii) Additionally, classifications of building use types do not always reflect the attributes of each specific building use. An example would be a university campus, which contains a wide variety of buildings, all of which are classified as one occupation, or a gymnasium that may be in a high-rise, one-story structure. Furthermore, the height of the teaching building in schools is often higher than that of the office. (iii) The third point is that the top floor of some buildings has a decorative roof, and the number of floors is rounded off incorrectly.
Confusion between structures primarily exists among confined masonry, reinforced concrete, and shear wall structures. This is particularly true for low-rise dwellings, which usually have a variety of structural types and a close geometrical arrangement. Since high-rise buildings are commonly shear wall structures, the number of stories plays an important role in improving their classification accuracy.
As seen from the result, there has been a significant change in the building stock within the study area from 2011 to 2021. This change could be attributed to several factors. First, the gross floor area of buildings decelerated between 2014 and 2016 but subsequently grew rapidly. A second characteristic of the change in building types is the decreased number of brick wood structures, followed by a significant rise in the number of shear wall structures, while other types of buildings continue to steadily change. Moreover, most of the buildings that were reduced or added were residential. Third, the renewal of buildings occurs in a variety of areas throughout the city, from the center to the edge.
The above analysis results are consistent with the period of demolition and reconstruction in this area in the past. According to the statistical yearbook of the Baodi District, the district launched a relocation and reconstruction project in October 2015, involving 25,000 households and 66,000 residents, with an investment of 43.3 billion RMB yuan (approximately 6551 million EURO). A total of 35 urban villages and 21 dormitories were demolished in 2016, and, in 2020, 35 replacement communities were built, which are currently in operation.
In this study, the calculation to estimate building area change was based on the extraction of the building's footprint, the estimation of the number of floors, and the classification of building types. Additionally, a change in the building use is implied. This is a major difference from previous methods of monitoring urbanization expansion by identifying impervious layers. We also consider this to be a fundamental issue in our research, and it forms the basis for the determination of variations in specific earthquake losses.
There are still some limitations in our study that need to be acknowledged. During remote sensing analysis, first, we adopted a large private dataset to train the image segmentation model and classification model, which is a time-consuming procedure that needs to be performed carefully. This localized dataset has proven to be indispensable in our research. Second, even though high-rise buildings and one-story rural dwellings are classified with high accuracy, building structure type confusion exists among confined masonry and reinforced concrete.
In calculations of population and economic losses, no earthquake simulation would enable the results to be verified, as in the case of [79]. This means that the results of the calculation presented in this paper are only theoretically reliable. Second, since the study area is relatively small, we perform calculations by setting a consistent PGA value rather than using the conventional seismic hazard model [80,81].
In future work, we will study transfer learning techniques to reduce the dependence on sample size and examine the attention module. We will also study the subdivision of the vulnerability function, as well as the spatialization of population data based on remote sensing.

Conclusions
Identifying and mastering the changes in earthquake risk due to urbanization is crucial to effectively adjust countermeasures to reduce earthquake losses and execute earthquake emergency preparedness in a targeted manner. To achieve this goal, we propose an integrated workflow incorporating deep learning and ensemble learning methods for remote sensing image analysis. As an example, the urbanization process of the Baodi core area in Tianjin was studied using remote sensing images taken in recent years. The type of building and changes in building area were analyzed from the extracted building information, and the economic and social losses resulting from these changes were calculated at three ground-motion intensity levels. Yet, the seismic loss calculation was made based on certain parameters; additional factors such as geology, soil type, distance from the seismogenic fault, seismic wave propagation, and secondary hazards were not taken into consideration.
According to our research, the study area has changed from a county seat to a dormitory town of a big city as a result of urbanization. In the study area, the development of commercial real estate eliminated low-quality housing and nearly doubled the number of people that could be accommodated, but the growth rate of functional infrastructure was relatively low.
New residential buildings have resulted in a reduction in the rate of night-time earthquake deaths in this area. The spatial distribution of the probability of a fatality also changes accordingly. However, there is a possibility that the number of deaths caused by rare earthquakes may not significantly decrease, given the dramatic increase in population. This may be overlooked in the context of urbanization.
Earth observation data complement ground-collected data and play a pivotal role in risk assessment and reduction [82]. Moreover, our study demonstrates that remote sensing data can be a valuable resource to observe changes in seismic risk as a result of urbanization. Although there are still some areas for improvement, the method based on remote sensing data can be used as a tool for updating seismic risk management plans and urban planning in other parts of China in the future.