Building seismic capacity analysis method and system based on unmanned aerial vehicle image

Abstract

The invention provides a building seismic capacity analysis method and system based on an unmanned aerial vehicle image, and provides a method for extracting contour vectors and seismic attributes of a single building by combining technologies such as a region growing algorithm, buffer area analysis, a building land index, gray texture, a seismic damage factor method and a machine learning algorithm. And analyzing the earthquake resistance of the building by using two methods of an earthquake damage factor method and a machine learning classifier aiming at the obtained various earthquake resistance attributes. According to the invention, full-automatic building contour segmentation, anti-seismic attribute extraction and anti-seismic capacity analysis are realized, good effects can be obtained, and the requirement of rapid analysis of large-scale building anti-seismic capacity in China is met.

Description

Technical Field

[0001] The present invention relates to the field of remote sensing technology and building seismic analysis, and in particular to a method and system for analyzing the seismic capacity of buildings based on unmanned aerial vehicle images. Background Art

[0002] Conducting a rapid analysis of the seismic resistance of buildings before an earthquake and fully understanding the spatial distribution of buildings in the region and their seismic resistance are of great practical significance for improving the region's earthquake emergency preparedness, enhancing earthquake loss prediction and post-earthquake rapid response capabilities.

[0003] The types of data involved in building earthquake damage investigation and analysis are complex, but my country's current building basic information is not perfect, and many regions have not yet formed a complete building earthquake damage basic database. The traditional household survey method conducts on-site investigations of housing conditions in accordance with a series of standards and processes such as the "Survey Standards for Seismic Fortification of Existing Buildings" and the "Building Seismic Design Code". It not only requires a large amount of human resources and is inefficient; the data obtained is characterized by individualization and fragmentation, which is far from meeting the needs of my country for large-scale building earthquake resistance analysis.

[0004] The drone technology that has emerged in recent years has the advantages of strong adaptability, low acquisition cost, fast response speed, and high image spatial resolution. It has extremely high value in the field of regional information refinement. With the help of camera parameters, positioning attitude system (POS) data and homonymous image points, drone images taken along the planned route can generate extremely high-resolution orthophotos (Digital Orthophoto Map, DOM); in addition, since the route has a certain degree of repetition, there will be drone images from multiple angles in the same area. Using the principle of stereo image pairs, a digital surface model (Digital Surface Model, DSM) can be generated; using the above data, a variety of seismic factors of buildings can be extracted, so as to conduct seismic capacity analysis. However, the current application of drone images in the seismic capacity analysis of buildings is still limited to building extraction and visual interpretation, which still requires high labor costs. Summary of the invention

[0005] The present invention provides a method for quickly analyzing the seismic resistance of buildings based on drone images. The method utilizes high-resolution images acquired by drones and, through deep learning semantic segmentation technology, automatically extracts building outlines and related seismic properties, and then evaluates the seismic resistance of buildings.

[0006] The present invention also provides a method for analyzing the seismic capacity of a building, comprising the following steps:

[0007] Step 1: Extract the outline and seismic attributes of single-family buildings from UAV images to obtain the length, width, area, height, construction age, and building structure of the buildings;

[0008] Step 2: For different situations, select different building seismic capacity analysis methods for processing. For the case where a building basic database already exists, use machine learning methods to analyze the building seismic capacity. For the case where there is a lack of a building basic database, analyze the building seismic capacity based on earthquake damage factors;

[0009] Step 2.1: Establish a building seismic capacity sample library through on-site investigation and visual interpretation, train a machine learning classifier, and establish the relationship between attributes and seismic capacity levels;

[0010] The input of the machine learning classifier is the height, area, width-to-height ratio, length-to-width ratio, building structure, and construction age of the building. The building seismic capacity levels are divided into three levels: estimated seismic capacity up to standard, suspected seismic capacity insufficient, and suspected seismic capacity seriously insufficient;

[0011] Step 2.2: Use the improved earthquake damage factor method to establish the relationship between the damage degree of the building and the earthquake damage factors. The earthquake damage factor method believes that there is no linear relationship between the damage degree of the building and the earthquake damage factors, and the damage degree of the building can be represented by the product of the weights of each earthquake damage factor, and the corresponding weights are calculated by the least squares method;

[0012]

[0013] Where D is the damage index of the building; w is the conversion coefficient of ground motion intensity, d 0 is the correction coefficient; N is the number of types of earthquake damage factors, including building structure, number of building floors, construction age, area, and width-to-height ratio; T is the classification number of the current earthquake damage factor; d ij is the weight value of the jth classification of the ith earthquake damage factor; m ij is the power exponent, which takes 1 when the actual situation of the building conforms to the jth classification of the ith earthquake damage factor, and 0 in other cases;

[0014] After calculating the damage index D of the building through the above formula, classify the seismic capacity of the building according to the following criteria: D < a1 indicates that the estimated seismic capacity meets the standard, a1 ≤ D < a2 indicates that the suspected seismic capacity is insufficient, and D ≥ a2 indicates that the suspected seismic capacity is seriously insufficient.

[0015] Furthermore, the specific implementation method of Step 1 is as follows;

[0016] Step 11: Perform semantic segmentation of ground object types based on DOM and DSM images to obtain the semantic segmentation results;

[0017] Step 12, convert the DOM image in RGB color space to HSV color space, and then normalize the Hue band, DSM and semantic segmentation results to 0-255 to synthesize a pseudo-color image that can simultaneously reflect spectral information, height information and prediction information;

[0018] Step 13, processing the synthesized pseudo-color image by using the Mean Shift clustering algorithm, so that the texture of the pseudo-color image becomes smoother and the edges between regions become sharper;

[0019] Step 14, using a region growing algorithm to extract a single building based on the synthetic pseudo-color image obtained above;

[0020] Step 15, perform morphological closing operation on the grown area to eliminate internal holes, smooth the contour, screen through morphological feature parameters such as area, aspect ratio, and rectangularity, and finally use boundary tracking algorithm to obtain the building outline vector;

[0021] Step 16, extracting the length, width and area of ​​the building according to the building vector, and obtaining the building height through a buffer analysis algorithm with the assistance of the semantic segmentation result;

[0022] Step 17, calculating the enhanced building land index through Landsat time series image data to obtain the age of the building area, and then infer the age of the building;

[0023] Step 18, classify the building roof based on its grayscale texture, and then infer the building structure based on the acquired building attributes.

[0024] Furthermore, in step 13, in the traditional Mean Shift algorithm, each pixel is represented by a five-dimensional vector [X, Y, R, G, B], the first two elements are the coordinates of the pixel in the image, and the last three elements are the color components of the pixel, which indicates that the weights of the three components R, G, and B are consistent. However, the complexity of the three bands in the pseudo-color image is different. Therefore, the color component weight is added on the basis of the traditional Mean Shift algorithm, and a weighted five-dimensional vector [X, Y, Pr×R, Pg×G, Pb×B] is used to represent the pixel, so as to obtain a better clustering effect, where Pr, Pg, and Pb represent the weights of the three components R, G, and B, respectively.

[0025] Furthermore, the calculation process of the building length and width, building area and building height in step 16 is as follows;

[0026] 1) Building length and width: The length and width of the building's minimum circumscribed rectangle are taken as the length and width of the building;

[0027] 2) Building area: Fill the building outline vector with a line width of -1 to obtain the building area without any holes. Count the number of pixels in the building area and combine it with the resolution of the image to obtain the actual area of ​​the building.

[0028] 3) Building height: Generate a circular buffer around the building, find ground points around the building through buffer analysis, and use the average value of the DSM of the ground points as the ground height.

[0029] Furthermore, the specific implementation of step 17 is as follows:

[0030] Use ρ to represent the reflectivity after atmospheric correction and radiation correction, ρ blue represents the blue band, ρ green represents the green band, ρ red represents the red band, ρ nir represents the near-infrared band, ρ mir1 represents the short-wave infrared 1 band, ρ mir2 Indicates shortwave infrared 2 band;

[0031] The modified normalized water index MNDWI is used to extract water body information and generate a mask. The building land extracted by the enhanced exponential building land index EIBI is masked to eliminate the erroneous water body information. The formula is as follows:

[0032]

[0033] Among them, L is the soil adjustment factor, which adjusts the soil brightness difference; V NDBBI is the value of the normalized difference bare land and built-up land index NDBBI. The other Vs have the same meaning and represent the value of the subscript index.

[0034] NDBBI is an index (normalized difference bareness and built-up index, NDBBI), which is used to replace NDBI. Bare land has the highest reflectivity in the shortwave infrared 1 band, so its value in NDBI is also the largest. The building land extracted by NDBI is bound to include bare land. In addition, NDBI will be affected by water shortage and sparse vegetation, resulting in vegetation noise. EBSI is an enhanced bare soil index (EBSI), which is constructed using the property of bare land with the largest difference between the bare soil index (BSI) and MNDWI.

[0035]

[0036] Furthermore, in step 18, the gray-level co-occurrence matrix is ​​first used to calculate and analyze the texture features of the roof, and the category of the entire building roof is judged based on the texture features of the roof, specifically: the roof is sliced, the type of each cropped small area is judged, and finally the type mode of the entire roof area is taken as the entire roof type; finally, the building is classified according to the combination of roof type, building height and construction age, and the building structure is divided into frame structure, brick-concrete structure and brick-wood structure, specifically: buildings with a building height of more than 45 meters are directly judged as frame structures; buildings less than 10 meters, built before 1986 and with old tiles on the roof are judged as brick-wood structures; the remaining buildings are judged as brick-concrete structures.

[0037] Furthermore, the calculation formula of the earthquake intensity conversion coefficient is as follows:

[0038]

[0039] Where a is the peak acceleration of the building, g is the gravitational acceleration, and the peak acceleration a is calculated by the earthquake attenuation model, which is obtained by the following method:

[0040] First, based on the earthquake magnitude, focal depth, epicenter location, and distance between the building and the focal point, the ground peak acceleration at the building's location is estimated using the seismic attenuation relationship, and then the peak acceleration is adjusted using the correction coefficient and the building's dynamic characteristics.

[0041] Furthermore, when using the least squares method to solve the problem, first take the derivatives on both sides of the equation (1) to convert the multiplication nonlinear multiple regression problem into a linear multiple regression problem. The correction coefficient and the earthquake damage factor weight can be solved by the least squares method.

[0042] Furthermore, the values ​​of all earthquake damage factors and weights are as follows:

[0043]

[0044]

[0045] The present invention also provides a building seismic capacity analysis system based on drone images, comprising:

[0046] one or more processors;

[0047] A storage device is used to store one or more programs. When the one or more programs are executed by the one or more processors, the one or more processors implement a method for analyzing the seismic resistance of buildings based on drone images as described in the above technical solution.

[0048] The method for analyzing the seismic capacity of buildings based on drone images provided by the present invention combines regional growing algorithm, buffer zone analysis, building land index, grayscale texture, earthquake damage factor method and machine learning algorithm and other technologies, and proposes a method for extracting the outline vector and seismic properties of a single-family building. For the obtained multiple seismic properties, the seismic capacity of the building is analyzed using the earthquake damage factor method and machine learning classifier. The present invention realizes fully automatic building outline segmentation, seismic property extraction and seismic capacity analysis, which can achieve good results and meet the needs of my country for rapid analysis of large-scale building seismic capacity. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] In order to more clearly illustrate the technical solutions in the present invention or the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0050] Figure 1 The figure is a flow chart of an embodiment of the present invention.

[0051] Figure 2 It is a schematic diagram of a statistical graph of building shape characteristic indexes provided by the present invention;

[0052] Figure 3 is a schematic diagram of the tracking direction of the boundary tracking algorithm provided by the present invention;

[0053] Figure 4 is a schematic diagram of common building roof categories provided by the present invention;

[0054] Figure 5 It is a statistical diagram of the color tone and texture characteristics of roof samples provided by the present invention;

[0055] Figure 6 It is a schematic diagram of a building structure discrimination decision tree model provided by the present invention;

[0056] Figure 7 It is a schematic diagram of the classification of earthquake resistance levels provided by the present invention. DETAILED DESCRIPTION

[0057] In order to make the purpose, technical solution and advantages of the present invention clearer, the technical solution of the present invention will be clearly and completely described below in conjunction with the drawings of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0058] Step 1: Extract the outline and seismic properties of a single-family building to obtain the length, width, area, height, construction age, and building structure of the building;

[0059] Step 11, perform semantic segmentation of ground object types based on DOM and DSM images to obtain semantic segmentation results;

[0060] Step 12, convert the DOM image in RGB color space to HSV color space, and then normalize the Hue band, DSM and semantic segmentation results to 0-255, so as to synthesize a pseudo-color image that can simultaneously reflect spectral information, height information and network prediction information.

[0061] In the HSV color space, H (Hue) represents hue, which represents the basic hue of the color; S (Saturation) represents saturation, which represents the degree to which the color is close to the spectral color; V (Value) represents brightness, which represents the brightness of the color. After normalizing RGB to 0-1, the specific conversion method is as follows:

[0062] V=max(R,G,B)

[0063]

[0064] The pseudo-color image obtained by observation shows that the spectral information is reflected in the red band. The red components of different roofs are different, which makes the buildings show obvious differences. The height information is reflected in the green band. The higher vegetation appears green in the image, and the green component gradually increases with the height of the building area, and the color of the building gradually becomes lighter. The prediction information of the network is reflected in the blue band. The building area is mostly blue or purple, while the non-building area is mostly ground points, with hue tending to 0, height tending to 0, and no blue component, so it is generally black, which makes the building area and non-building area also have obvious differences. This shows that the band synthesis method proposed in this method can effectively segment the outline of a single building.

[0065] Step 13: Process the synthesized pseudo-color image by using the Mean Shift clustering algorithm, so that the texture of the pseudo-color image becomes smoother and the edges between regions become sharper.

[0066] For three-band pseudo-color images, the bandwidth parameters of the Mean Shift clustering algorithm are the coordinate space bandwidth (i.e., the calculation range of the mean shift vector) and the color space bandwidth (i.e., the degree of merging of colors within the calculation range).

[0067] In the traditional Mean Shift algorithm, each pixel is represented by a five-dimensional vector [X, Y, R, G, B]. The first two elements are the coordinates of the pixel in the image, and the last three elements are the color components of the pixel. This shows that the weights of the three components R, G, and B are consistent. However, the three bands in the pseudo-color image have different complexities and need to be smoothed to different degrees. The impact on clustering should also be different, that is, the weights of the three bands should be different. Therefore, this method adds color component weights to the traditional MeanShift algorithm and uses a weighted five-dimensional vector [X, Y, P r ×R,P g ×G,P b ×B] to represent the pixel points, so as to obtain a better clustering effect, where P r , P g , P b Represent the weights of the three components R, G, and B respectively.

[0068] Step 14: extract the single-family building based on the synthetic pseudo-color image obtained above using a region growing algorithm. The basic idea of ​​the region growing algorithm is to group together pixels with similar features, thereby obtaining a connected region with the same features.

[0069] The region growing algorithm starts from a seed point, compares the pixels in the neighborhood around the seed point with the current region, merges the points that meet the requirements into the current region, and continues to grow outward using these points as seed points until there are no more pixels that can be merged, thus obtaining a fully grown region. The effectiveness of the region growing algorithm depends on the following three aspects:

[0070] 1) Selection of seed points.

[0071] 2) Determination of adjacent pixels and their merging criteria during the growth process.

[0072] 3) Conditions under which growth stops.

[0073] The seed point in this method is randomly selected. A point is randomly selected from the building area in the network prediction result as the seed point to start the region growing algorithm. The non-building area and the area that has been grown are not within the randomly selected range.

[0074] When the seed point grows, the four connected areas of top, bottom, left and right are selected as the growth direction. Compared with the eight-neighborhood, the four-neighborhood will not be stuck, and the grown area will be more compact. The area that has been grown is not within the growth range.

[0075] When judging whether a certain neighborhood pixel can be merged, its pixel value is compared with the mean pixel value of the current area. If it is within a certain threshold range, it is allowed to be merged. The pseudo-color image obtained in the previous section has three bands, which is not convenient for calculating the mean, so it is grayed out and the pixel value is calculated.

[0076] When there are no pixels that can be merged in the current region, the region growing algorithm of the seed point stops and the region belonging to the seed point is obtained. Then the seed point is randomly selected again and the growth process is repeated until there are no new seed points to be selected and the entire region growing algorithm ends.

[0077] Step 15, perform morphological closing operation (dilation first, then erosion) on the grown area to eliminate internal holes, smooth the contour, screen it through morphological feature parameters such as area, aspect ratio, rectangularity, etc., and finally use the boundary tracking algorithm to obtain the building outline vector.

[0078] The area obtained by growth is not necessarily a building area, but may be a small noise point or the edge of a building. Therefore, the following shape feature index is used to eliminate these situations to improve the accuracy of the extraction results. The statistical results of the shape feature index of some buildings in the sample survey are as follows: Figure 2 shown.

[0079] 1) Area S. This is expressed by counting the number of pixels in the area, mainly used to remove noise areas. Buildings usually have a large area. According to statistical results, the area of ​​common buildings is generally more than 20 square meters. Combined with the resolution of TIFF images, the square meter threshold is converted into a pixel number threshold.

[0080] 2) Aspect ratio R. The ratio of the length and width of the minimum circumscribed rectangle (rotatable) of the growing area is calculated, which is mainly used to remove the elongated areas at the edge of the building. Such areas often have a certain area, and the area threshold cannot handle this situation, but their aspect ratio is generally high, while the aspect ratio of the building is generally less than 7, which is significantly different from the other two.

[0081] 3) Rectangularity F. It is calculated by dividing the area of ​​the region by the area of ​​its minimum circumscribed rectangle. It is mainly used to eliminate the edge area passing through the corner of the building. This type of area is usually the connecting edge of the long side and the short side of a building. The area and aspect ratio are relatively normal, but its rectangularity is low. The rectangularity of a building is generally higher than 0.4.

[0082] After the closing operation and morphological feature index screening, the building area can be vectorized. Since the building is obtained through the region growing algorithm, the entire building area must be a single connected closed figure, so the boundary tracking algorithm can be used to obtain the outer contour of the area. According to the geodetic coordinates corresponding to the upper left corner of the TIFF image and the resolution in the X and Y directions, the pixel coordinates of the contour can be converted into geodetic coordinates, and then the contour vector of the single building can be obtained.

[0083] The boundary tracking algorithm used in this method sets the lower left area point as the starting point, and the initial tracking direction is set to 45 degrees from the upper left corner. The tracking direction sequence is as follows: Figure 3 As shown. Starting from the starting point, detect the region point along the initial tracking direction. If the region point is not in front of the current tracking direction, rotate 45 degrees clockwise and re-detect until the region point is found. Save the current starting point, use the detected point as the new starting point, rotate the tracking direction 90 degrees counterclockwise as the new tracking direction, and repeat the above detection and recording process until the detection point returns to the original starting point, and the boundary tracking algorithm ends.

[0084] The resulting vectors were normalized using the Regularize Building Footprints tool in ArcPY, which uses a polyline compression algorithm to remove artifacts and distortion from the building vectors.

[0085] Step 16, extract the length, width and area of ​​the building according to the building vector, and obtain the building height through the buffer analysis algorithm with the assistance of the semantic segmentation result.

[0086] 1) Building length and width. Since most buildings are relatively regular rectangles, the length and width of the building’s smallest circumscribed rectangle can be roughly considered as the building’s length and width.

[0087] 2) Building area. Morphological closing operation can be used after region growing to eliminate some smaller internal holes. However, due to factors such as incorrect semantic segmentation network prediction results and elevator shafts, large holes may appear inside the building obtained by region growing, which cannot be eliminated by morphological closing operation. Therefore, the area obtained by region growing cannot be directly used as the building area. The building outline vector is filled with a line width of -1 to obtain a building area without any holes. The number of pixels in the area is counted and combined with the resolution of the image to obtain the actual area of ​​the building.

[0088] 3) Building height. The traditional building height extraction method is to calculate the difference between the DSM and DEM average values ​​of the building area, that is, the difference between the absolute building height and the absolute ground height. This method requires the DEM of the corresponding area, but the existing DEM resolution is low and it is difficult to use with the high-resolution DSM; the method of automatically generating DEM using DSM filtering has poor results. It is necessary to manually select the building area and vegetation area and manually adjust the number of filters to obtain a better result, which is contrary to the idea of ​​fully automatic extraction and analysis of this method. Therefore, this method does not use the average DEM value of the building area as the ground height, but searches for ground points around the building through buffer analysis, and uses the average value of the DSM of the ground point as the ground height. The final building height is = (DSM of the building area) - (average value of the DSM of the ground point in the buffer zone of the building area). When generating the buffer zone around the building, since the extraction result is not very accurate at the edge of the building, it is chosen to establish a ring buffer zone to exclude the possible unextracted building edges to avoid being misjudged as ground points.

[0089] Step 17, calculate the enhanced building land index through Landsat time series image data to obtain the age of the building area, and then infer the building age.

[0090] We selected the Landsat satellites with the most historical data so far, and used Google Earth Engine (GEE) to obtain the CollectionLevel2 surface reflectance products of Landsat5 and Landsat8 after geometric correction, radiometric calibration, and atmospheric correction. We used CLOUD_COVER data to screen out images with less than 20% cloud cover, used the QA_PIXEL band to generate cloud masks, and cropped each image. Finally, we averaged and synthesized the images throughout the year to represent the situation in the study area in that year.

[0091] For the convenience of the following description, ρ will be used to represent the reflectivity after atmospheric correction and radiation correction. blue represents the blue band (TM-1, OLI-2), ρ green represents the green band (TM-2, OLI-3), ρ red represents the red band (TM-3, OLI-4), ρ nir represents the near infrared band (TM-4, OLI-5), ρ mir1 Indicates shortwave infrared 1 band (TM-5, OLI-6), ρ mir2 Indicates shortwave infrared 2 band (TM-7, OLI-7).

[0092] The modified normalized difference water index (MNDWI) is used to extract water body information and generate a mask. The building land extracted by the enhanced index-based built-up index (EIBI) is masked to eliminate the erroneous water body information. The formula is as follows:

[0093]

[0094] Among them, L is the soil adjustment factor, which can adjust the soil brightness difference and is generally set to 0.5; V NDBBI is the value of the index NDBBI. The other Vs have the same meaning and represent the value of the subscript index.

[0095] NDBBI is the normalized difference bareness and built-up index (NDBBI), which is used to replace the normalized building index NDBI. The reflectivity of bare land is the highest in the shortwave infrared 1 band, so its value in NDBI is also the largest. The built-up land extracted by NDBI is bound to include bare land. In addition, NDBI is also affected by water shortage and sparse vegetation, resulting in vegetation noise. EBSI is the enhanced bare soil index (EBSI), which is constructed using the property of bare land that has the greatest difference between the bare soil index (BSI) and the MNDWI.

[0096]

[0097] In addition, due to the low spatial resolution of Landast images, building land may be missed in areas such as farmland and woodland. In the above identification and extraction of building outlines based on drone images, small houses in farmland or woodland may be detected. In this case, the buildings will lack the corresponding building area age, and the corresponding building age cannot be inferred. Therefore, unknown numbers are used instead and they are not involved in subsequent analysis.

[0098] Step 18, classify the building roof based on its grayscale texture, and then infer the building structure based on the acquired building attributes.

[0099] Due to the difficulty of using drone images to identify building structures, it is difficult to make detailed divisions. Therefore, the building structures are simply divided into frame structures, brick-concrete structures, and brick-wood structures with reference to the current standards. The frame structure uses reinforced concrete as the main load-bearing structure, which has the toughness of steel bars and the compressive resistance of concrete. It has high structural rigidity and excellent bearing capacity, can effectively prevent structural deformation and damage, and has good earthquake resistance; the brick-concrete structure mainly uses brick walls to bear vertical loads, with poor toughness. When the earthquake intensity is 6-7 degrees, the wall will crack, and it will collapse when the intensity is above 9 degrees; the existing brick-wood structure is generally an older house, using wooden beams and wooden columns to bear the vertical load of the building, supporting the stability of the brick wall, and the earthquake resistance is extremely poor. The wood and brick wall cannot work together effectively, and the wood is prone to large bending deformation, which will collapse when the earthquake intensity is 6-7 degrees.

[0100] There is often a close relationship between building roofs and building structures. Buildings of the same period and structure tend to use the same building roofs. Common building roofs are Figure 4 As shown, the colors and textures of the roofs of these buildings are different and can be classified by hue and texture features.

[0101] Texture refers to the grayscale distribution pattern of pixels within a certain range in space. In order to describe the texture characteristics of the roof, the Gray-Level Co-occurrence Matrix (GLCM) is used for calculation and analysis. Each element in the grayscale co-occurrence matrix represents the co-occurrence probability between pairs of pixels at fixed distances and directions. The probability value of the elements at the diagonal of the image with slower texture changes is larger, while the probability value of the elements on both sides of the diagonal of the image with faster texture changes is larger. Therefore, the entire grayscale co-occurrence matrix contains the spatial relationship information between various grayscale levels, reflecting the distribution and position characteristics of grayscale, from which the local pattern and arrangement rules of the image can be analyzed. Due to the large amount of data in the grayscale co-occurrence matrix, it is not convenient to directly compare and classify. The following derived statistics are usually used as texture classification features:

[0102] 1) Energy, also known as angular second moment (ASM). Energy is the square root of ASM and is used to describe the thickness and uniformity of grayscale texture. When the elements in the grayscale co-occurrence matrix are discretely distributed or the element values ​​are similar, the energy value is small, indicating that the texture is fine or irregular.

[0103] ASM=∑ i ∑ j p(i,j) 2

[0104] Where p(i,j) represents the element value of the gray-level co-occurrence matrix at (i,j).

[0105] 2) Contrast: It indicates the difference in grayscale levels between adjacent pixels and is used to describe the clarity of image texture. The greater the texture contrast, the clearer the effect and the greater the contrast.

[0106] Contrast=∑ i ∑ j p(i,j)×(ij) 2

[0107] 3) Entropy: It is used to describe the complexity of an image. The more complex the image, the greater the entropy value.

[0108] Entropy=∑ i ∑ j p(i,j)×ln(p(i,j))

[0109] 4) Homogeneity: It is used to describe the homogeneity of image texture. The greater the variation between image texture regions, the smaller the value of inverse variance.

[0110]

[0111] 5) Correlation: It is used to describe the similarity of the horizontal or vertical distribution of the grayscale of an image. The greater the correlation between local grayscales, the greater its value.

[0112]

[0113] Mean=∑ i ∑ j p(i, j)×i

[0114] Variance=∑ i ∑ j p(i, j)*(i-Mean) 2

[0115] In actual calculation, since the grayscale of the image 0-255 is too large, the calculated grayscale co-occurrence matrix will be too large and the calculation efficiency will be low. Therefore, it is necessary to compress the grayscale of the image, compress the image grayscale to 0-25, a total of 26 levels, and then calculate the grayscale co-occurrence matrix in four directions (0°, 45°, 90° and 135°), statistically analyze the texture classification features, and take the extreme values ​​in the four directions as the final result.

[0116] Figure 4 The normalized results of the average values ​​of the hue and texture features of the nine roof samples are shown in Figure 5As shown in the figure, the color and texture classification features of different roofs are obviously different, which can be classified by machine learning. Since these attributes are continuous real-valued features, they are more in line with the applicable scope of Gaussian Bayesian classifier, so the naive Gaussian Bayesian classifier is selected for classification. The classification accuracy of the final test set is 88.24%, the precision is 90.63%, the recall rate is 89.38%, and the F1 score is 90.00. When judging the category of the entire building roof, since the entire roof situation is relatively complex, there may be multiple roof types, resulting in inaccurate classification. Therefore, the roof is sliced, and the type judgment is performed on each small area cut out. Finally, the type mode of the entire roof area is taken as the entire roof type.

[0117] Through the summary of roofs of different building structures, it is found that brick-wood buildings all use older tile roofs. The roof texture is also related to the building structure. Figure 6 In the building structure identification decision tree shown, buildings with a height of more than 45 meters are directly judged as frame structures; buildings less than 10 meters, built before 1986 and with old-style tiles on the roof are judged as brick-wood structures; the rest of the buildings are judged as brick-concrete structures.

[0118] Step 2: Select different building seismic capacity analysis methods for different situations. For situations where there is already a basic building database, use machine learning methods to analyze the building seismic capacity. For situations where there is no basic building database, use earthquake damage factors to analyze the building seismic capacity.

[0119] Step 2.1, establish a sample library of building seismic resistance through field investigation and visual interpretation, train machine learning classifiers such as multi-layer perceptron, and establish the relationship between attributes and seismic resistance level.

[0120] The input of the machine learning classifier is the building's height, area, aspect ratio, length-to-width ratio, building structure, and construction age. The building's seismic capacity is divided into three levels: estimated seismic capacity meets the standard, suspected seismic capacity is insufficient, and suspected seismic capacity is seriously insufficient.

[0121] This method divides the seismic capacity of buildings into three levels: estimated seismic capacity that meets the standard, suspected seismic capacity that is insufficient, and suspected seismic capacity that is seriously insufficient, according to the seismic capacity classification method of the China Earthquake Disaster Prevention Center. Figure 7 shown.

[0122] Buildings that are estimated to meet the earthquake resistance standards are generally newer communities in cities and towns, large public buildings, special buildings, and earthquake-resistant housing in rural areas. These buildings are mostly frame structures and brick-concrete structures, built in strict accordance with earthquake-resistant design specifications, and adopt a variety of earthquake-resistant fortification measures. Newly built residential communities are mostly regular rectangles in drone images, with consistent spacing between buildings, neat distribution, and basically the same height; large public buildings, special buildings, etc. are mostly single large buildings with neat streets around them.

[0123] Buildings suspected of insufficient earthquake resistance are generally low-rise houses in towns, public buildings, small factories and mid-rise houses in rural areas. These buildings are basically brick-concrete structures. No earthquake-resistant design was carried out before construction, and basically no additional earthquake-resistant reinforcement measures were taken. In the image, these houses are spaced at different distances, distributed chaotically, and have uneven heights. The colors of the roofs are mainly red, blue, gray, and white, and the arrangement is messy. Public buildings and small factories are mostly small single buildings. There is no unified planning of the surrounding streets, and there are a large number of small roads interspersed.

[0124] Buildings suspected of being seriously inadequate in earthquake resistance are generally single-story buildings in towns and low-rise self-built houses in rural areas. Such buildings include old brick-concrete and brick-wood structures, with no earthquake-resistant measures at all. In the images, these buildings are closely connected, low in height and small in area, with old-style gray bricks and tiles on the roofs, mostly distributed in old urban areas or rural areas.

[0125] Step 2.2: Use the improved earthquake damage factor method to establish the relationship between the degree of damage to a building and the earthquake damage factor. The earthquake damage factor method assumes that there is no linear relationship between the degree of damage to a building and the earthquake damage factor. The degree of damage to the building can be represented by the product of the weights of each earthquake damage factor, and the corresponding weight is calculated by the least squares method. The earthquake damage factor method assumes that there is generally no linear relationship between the degree of damage to a building and these earthquake damage factors. The degree of damage to the building can be represented by the product of the weights of each earthquake damage factor.

[0126]

[0127] Where D is the damage index of the building; w is the seismic intensity conversion coefficient, a is the peak acceleration of the building, and g is the gravitational acceleration; d 0 is the correction coefficient; N is the number of types of earthquake damage factors, including building structure, number of floors (calculated using the building height obtained in step 16, the average floor height of the building is 3 meters, and the approximate number of floors of the building can be obtained by simply dividing the building height by the average floor height), building age, area and aspect ratio; T is the number of categories of the current earthquake damage factor; d ij is the weight value of the jth classification of the ith earthquake damage factor; m ijIt is a power index, which takes the value of 1 when the actual situation of the building meets the jth classification of the ith earthquake damage factor, and takes the value of 0 in other cases.

[0128] Among them, the peak acceleration a of the earthquake motion suffered by the building is calculated by the earthquake motion attenuation model. First, according to the magnitude of the earthquake, the focal depth, the epicenter location, and the distance between the building and the focal point, the ground peak acceleration (PGA) at the location of the building is estimated using the earthquake motion attenuation relationship. This value is then adjusted by the site correction factor (the geological conditions of the site such as soil type, topography, etc.), which may amplify or attenuate the earthquake motion. Finally, the dynamic characteristics of the building (such as structural type, vibration frequency) will also affect the actual peak acceleration of the earthquake motion.

[0129] When solving the problem, we can find the derivatives on both sides of the equal sign of the above formula, convert the multiplication nonlinear multiple regression problem into a linear multiple regression problem, and solve the correction coefficient and earthquake damage factor weight by the least squares method. Combined with the above solution method, the weight value of the new factor is calculated. All earthquake damage factors and weight values ​​are shown in the following table.

[0130] Table 1 Earthquake damage factors and their weight values

[0131]

[0132]

[0133] All the building properties involved in the table have been extracted before. After calculating the building's damage index D through the formula, the building's seismic resistance is graded according to the table below.

[0134] Table 2 Comparison table of damage index and seismic resistance level

[0135]

[0136] Step 3: Use the extracted buildings and their attributes to analyze the earthquake resistance of the buildings through the earthquake damage factor method or the trained machine learning classifier, and evaluate the accuracy of the analysis results.

[0137] Based on the previously extracted building attributes, six attributes including height, structure, area, aspect ratio, length-to-width ratio and age were selected for training in accordance with the "Code for Seismic Design of Buildings" and expert recommendations, and a classifier was used to judge the seismic resistance of the building.

[0138] The seismic capacity analysis of the extracted results in the study area is carried out. When evaluating the accuracy, the intersection of the two vectors is taken for statistics to obtain the confusion matrix of the seismic capacity analysis results of the seismic damage factors in the study area.

[0139] Results of building seismic capacity analysis experiments:

[0140] The seismic capacity analysis method of buildings based on seismic damage factors was used to analyze the seismic capacity of the test set in the sample, and the confusion matrix obtained is shown in Table 3.

[0141] Table 3 Test set earthquake damage factor seismic capacity analysis accuracy evaluation table

[0142]

[0143] As can be seen from the table, the earthquake damage factor method rarely classifies buildings suspected of insufficient or severely insufficient seismic resistance as buildings with estimated seismic resistance that meet the standards. It is easy to misjudge buildings suspected of insufficient seismic resistance. The overall accuracy is 84.23%, and the overall classification effect is good, which shows that the earthquake damage factor method is applicable in different regions.

[0144] The seismic capacity analysis of the extraction results of the three study areas was carried out. Since there are certain differences between the automatically extracted buildings and the manually outlined buildings, the intersection of the two vectors was taken for statistics during the accuracy assessment, and the confusion matrix of the seismic capacity analysis results of the earthquake damage factors in the study area was obtained, as shown in Table 4.

[0145] Table 4. Table of accuracy assessment of earthquake damage factor analysis of seismic capacity in the study area

[0146]

[0147] Compared with the accuracy assessment results of the test set, the overall accuracy and Kappa coefficient of the earthquake damage factor method were reduced when the actual earthquake resistance analysis was carried out in the study area, but the misclassification trend it showed was consistent with the test set. The main reason for the decline in accuracy was that the buildings extracted from the complex area and the artificially outlined samples were quite different, which made the extracted building area and aspect ratio inconsistent with the actual ones. However, its overall accuracy was still 73.40%, which could still reflect the spatial distribution of buildings in the area and their earthquake resistance to a certain extent.

[0148] Then, the machine learning-based building seismic capacity analysis method was tested. The two trained multi-layer perceptrons were used to analyze the seismic capacity of the study area. The same accuracy assessment method was used, and the confusion matrix obtained was shown in Table 5.

[0149] Table 5. Analysis accuracy of multi-layer perceptron seismic capability in the study area

[0150]

[0151] Also affected by the building extraction error, the accuracy of the multi-layer perceptron in determining the seismic capacity of the study area is also reduced, and the overall accuracy is reduced to 79.49%, which is slightly higher than the earthquake damage factor method.

[0152] Those skilled in the art will appreciate that embodiments of the present invention may be provided as methods, systems, or computer program products. Therefore, the present invention may take the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware. Moreover, the present invention may take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0153] The present invention is described with reference to the flowcharts and / or block diagrams of the methods, devices (systems), and computer program products according to the embodiments of the present application. It should be understood that each process and / or block in the flowchart and / or block diagram, as well as the combination of the processes and / or blocks in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to generate a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowchart and / or block diagram. Figure 1 A process or multiple processes and / or boxes Figure 1 These computer program instructions can also be stored in a computer-readable memory that can guide a computer or other programmable data processing device to work in a specific way, so that the instructions stored in the computer-readable memory produce a product including an instruction device, which implements the functions specified in the process. Figure 1 A process or multiple processes and / or boxes Figure 1 A function specified in one or more boxes.

[0154] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operating steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing instructions for implementing the process. Figure 1 A process or multiple processes and / or boxes Figure 1 The steps for the functions specified in one or more boxes.

[0155] Through the description of the above implementation methods, those skilled in the art can clearly understand that each implementation method can be implemented by means of software plus a necessary general hardware platform, and of course, it can also be implemented by hardware. Based on this understanding, the above technical solution is essentially or the part that contributes to the prior art can be embodied in the form of a software product, and the computer software product can be stored in a computer-readable storage medium, such as ROM / RAM, a disk, an optical disk, etc., including a number of instructions for a computer device (which can be a personal computer, a server, or a network device, etc.) to execute the methods described in each embodiment or some parts of the embodiments.

[0156] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some of the technical features therein. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for analyzing the seismic capacity of buildings based on drone images, characterized in that: The steps are as follows: Step 1: Extract the outline and seismic resistance attributes of single-family buildings based on UAV images to obtain the length, width, area, height, construction age, and building structure of the buildings; Step 2: For different situations, select different building seismic capacity analysis methods for processing. For the case where a building foundation database already exists, use machine learning methods to analyze the building seismic capacity. For the case where there is a lack of a building foundation database, analyze the building seismic capacity based on earthquake damage factors; Step 2.1: Establish a building seismic capacity sample library through on-site investigations and visual interpretation, train a machine learning classifier, and establish the relationship between attributes and seismic capacity levels; The input of the machine learning classifier is the height, area, width-to-height ratio, length-to-width ratio, building structure, and construction age of the building. The building seismic capacity levels are divided into three levels: estimated seismic capacity up to standard, suspected seismic capacity insufficient, and suspected severe seismic capacity insufficient; Step 2.2: Use the improved earthquake damage factor method to establish the relationship between the damage degree of the building and the earthquake damage factors. The earthquake damage factor method believes that there is no linear relationship between the damage degree of the building and the earthquake damage factors, and the damage degree of the building can be represented by the product of the weights of each earthquake damage factor, and the corresponding weights are calculated by the least squares method; where D is the damage index of the building; w is the seismic ground motion intensity conversion coefficient, d0 is the correction coefficient; N is the number of types of earthquake damage factors, including building structure, number of building floors, construction age, area, and width-to-height ratio; T is the classification number of the current earthquake damage factor; d ij is the weight value of the jth classification of the ith earthquake damage factor; m ij is a power index, which takes the value of 1 when the actual situation of the building meets the jth classification of the ith earthquake damage factor, and takes the value of 0 in other cases; After calculating the building damage index D through the above formula, classify the seismic capacity of the building according to the following criteria: D < a1 indicates that the estimated seismic capacity meets the standard, a1 ≤ D < a2 indicates that the seismic capacity is suspected to be insufficient, and D ≥ a2 indicates that the seismic capacity is suspected to be severely insufficient.

2. The method for analyzing the seismic capacity of a building based on drone images according to claim 1, characterized in that: The specific implementation method of Step 1 is as follows; Step 11: Perform semantic segmentation of ground object types based on DOM and DSM images to obtain the semantic segmentation result; Step 12: Convert the DOM image in the RGB color space to the HSV color space, and then normalize the Hue band, DSM, and semantic segmentation results to 0 - 255 to synthesize a pseudo-color image that can simultaneously reflect spectral information, height information, and prediction information; Step 13: Process the synthesized pseudo-color image through the Mean Shift clustering algorithm to make the texture of the pseudo-color image smoother and the edges between regions sharper; Step 14: Use the region growing algorithm to extract single-family buildings based on the synthesized pseudo-color image obtained above; Step 15: Perform morphological closing operations on the grown regions to eliminate internal holes and smooth the contours, screen through morphological feature parameters such as area, aspect ratio, and rectangularity, and finally use the boundary tracking algorithm to obtain the building outline vector; Step 16: Extract the building length, width, and area according to the building vector, and obtain the building height through the buffer analysis algorithm with the assistance of the semantic segmentation result; Step 17: Calculate the enhanced built-up land index through Landsat time-series image data to obtain the building area age, and then infer the construction age; Step 18, classify the building roof based on its grayscale texture, and then infer the building structure based on the acquired building attributes.

3. The method for analyzing the seismic capacity of a building based on drone images as claimed in claim 1, characterized in that: In step 13, in the traditional Mean Shift algorithm, each pixel is represented by a five-dimensional vector [X, Y, R, G, B]. The first two elements are the coordinates of the pixel in the image, and the last three elements are the color components of the pixel. This indicates that the weights of the three components R, G, and B are consistent. However, the complexity of the three bands in the pseudo-color image is different. Therefore, the color component weight is increased on the basis of the traditional Mean Shift algorithm, and a weighted five-dimensional vector [X, Y, Pr×R, Pg×G, Pb×B] is used to represent the pixel, so as to obtain a better clustering effect, where Pr, Pg, and Pb represent the weights of the three components R, G, and B, respectively.

4. The method for analyzing the seismic capacity of a building based on drone images according to claim 1, characterized in that: The calculation process of the building length and width, building area and building height in step 16 is as follows; 1) Building length and width: The length and width of the building's minimum circumscribed rectangle are taken as the length and width of the building; 2) Building area: Fill the building outline vector with a line width of -1 to obtain the building area without any holes. Count the number of pixels in the building area and combine it with the resolution of the image to obtain the actual area of ​​the building. 3) Building height: Generate a circular buffer around the building, find ground points around the building through buffer analysis, and use the average value of the DSM of the ground points as the ground height.

5. The method for analyzing the seismic capacity of buildings based on drone images according to claim 1, characterized in that: The specific implementation of step 17 is as follows; Use ρ to represent the reflectivity after atmospheric correction and radiation correction, ρ blue represents the blue band, ρ green represents the green band, ρ red represents the red band, ρ nir represents the near-infrared band, ρ mir1 represents the short-wave infrared 1 band, ρ mir2 Indicates shortwave infrared 2 band; The modified normalized water index MNDWI is used to extract water body information and generate a mask. The building land extracted by the enhanced exponential building land index EIBI is masked to eliminate the erroneous water body information. The formula is as follows: Among them, L is the soil adjustment factor, which adjusts the soil brightness difference; V NDBBI is the value of the normalized difference bare land and built-up land index NDBBI. The other Vs have the same meaning and represent the value of the subscript index. NDBBI is an index (normalized difference bareness and built-up index, NDBBI), which is used to replace NDBI. Bare land has the highest reflectivity in the shortwave infrared 1 band, so its value in NDBI is also the largest. The building land extracted by NDBI is bound to include bare land. In addition, NDBI will be affected by water shortage and sparse vegetation, resulting in vegetation noise. EBSI is an enhanced bare soil index (EBSI), which is constructed using the property of bare land with the largest difference between the bare soil index (BSI) and the MNDWI.

6. The method for analyzing the seismic capacity of a building based on drone images according to claim 1, characterized in that: In step 18, the gray-level co-occurrence matrix is ​​first used to calculate and analyze the texture features of the roof, and the category of the entire building roof is judged according to the texture features of the roof, specifically: the roof is sliced, the type of each cropped small area is judged, and finally the type mode of the entire roof area is taken as the entire roof type; finally, the building is classified according to the roof type, building height and construction age, and the building structure is divided into frame structure, brick-concrete structure and brick-wood structure, specifically: buildings with a building height of more than 45 meters are directly judged as frame structures; buildings with a height of less than 10 meters, a construction age before 1986 and old-style tiles on the roof are judged as brick-wood structures; The remaining buildings were determined to be brick-concrete structures.

7. The method for analyzing the seismic capacity of buildings based on drone images as claimed in claim 1, characterized in that: The calculation formula of the earthquake intensity conversion coefficient is as follows: Where a is the peak acceleration of the building, g is the gravitational acceleration, and the peak acceleration a is calculated by the earthquake attenuation model, which is obtained by the following method: First, based on the earthquake magnitude, focal depth, epicenter location, and distance between the building and the focal point, the ground peak acceleration at the building's location is estimated using the seismic attenuation relationship, and then the peak acceleration is adjusted using the correction coefficient and the building's dynamic characteristics.

8. The method for analyzing the seismic capacity of buildings based on drone images according to claim 1, characterized in that: When using the least squares method to solve the problem, first take the derivatives on both sides of the equation (1) to convert the multiplication of nonlinear multiple regression problems into linear multiple regression problems. The correction coefficient and earthquake damage factor weight can be solved by the least squares method.

9. The method for analyzing the seismic capacity of buildings based on drone images as claimed in claim 1, characterized in that: The values ​​of all earthquake damage factors and weights are as follows:

10. A building seismic capacity analysis system based on drone images, characterized by: one or more processors; A storage device for storing one or more programs, which, when executed by the one or more processors, enables the one or more processors to implement a method for analyzing the seismic resistance of buildings based on drone images as described in any one of claims 1 to 9.

Images