Urban land subsidence risk assessment method and device, storage medium and computer program product

Abstract

The invention discloses an urban land subsidence risk assessment method and device, a storage medium and a computer program product, and relates to the field of wide area land subsidence risk grade assessment. According to the scheme, a danger factor set and a vulnerable factor set are constructed by obtaining the height of a permanent scatterer point in a monitoring area, the deformation rate of the permanent scatterer point and the elevation data of a digital elevation model, and a land subsidence risk membership degree vector is calculated according to the danger factor set and the vulnerable factor set; and evaluating the land subsidence risk grade according to the land subsidence risk membership degree vector. According to the method, rapid assessment and visual expression of the urban wide-area land subsidence risk level can be realized, and reliable foundation support is provided for subsequent rapid mastering of urban wide-area land subsidence risk evolution and early warning of land subsidence key areas.

Description

Technical Field

[0001] The present invention belongs to the field of wide-area land subsidence risk level assessment, and in particular relates to an urban land subsidence risk assessment method, device, storage medium and computer program product. Background Art

[0002] With the rapid development and expansion of urban economy and infrastructure, the ground and underground spaces in cities are being continuously used for development and construction. With the introduction and implementation of new policies such as urban renewal planning, more high-rise and super-high-rise buildings will appear on the ground, and underground hub transportation facilities are also in a high-intensity development stage. With the increasing number of large-scale engineering construction projects concentrated in cities, road subsidence and slope instability are becoming more frequent. There are major safety hazards in the geological conditions such as groundwater level fluctuations and thick soft soil layers in cities. Therefore, regular and rapid urban wide-area ground subsidence monitoring and risk level assessment are of great significance to the development of urban underground projects and the safety of people's production and life.

[0003] At present, the technology used in actual projects is mainly based on manual point observation, and most of them can only passively conduct key observations based on expert opinions after early accidents occur. At this stage, the traditional methods of urban wide-area observation have application limitations such as insufficient manual efficiency, labor-intensive and costly maintenance, and large lags.

[0004] By constructing a systematic and comprehensive assessment method, the risk of wide-area ground subsidence in the city is scientifically assessed and risk levels are divided, providing an effective technical means for active screening and investigation of subsidence risk areas in the city. Summary of the invention

[0005] The main purpose of the present invention is to propose a method for urban land subsidence risk assessment, aiming to solve the application limitations of the current urban wide-area traditional observation methods, such as insufficient labor efficiency, laborious and costly maintenance, and large lag.

[0006] To achieve the above object, the present invention proposes a method for assessing urban land subsidence risk, which comprises:

[0007] Obtain the permanent scatterer point height, permanent scatterer point deformation rate, land development data and digital elevation model elevation data within the monitoring area;

[0008] Constructing a set of hazard factors and a set of vulnerability factors according to the permanent scatterer point heights, permanent scatterer point deformation rates, land development data, and digital elevation model elevation data;

[0009] Calculating a ground subsidence risk membership vector according to the hazard factor set and the vulnerability factor set;

[0010] The ground subsidence risk level is evaluated according to the ground subsidence risk membership vector.

[0011] Optionally, the step of constructing a set of hazard factors and a set of vulnerability factors according to the deformation rate of near-ground permanent scatterer points specifically includes:

[0012] Calculate the relative elevation of each permanent scatterer point based on the digital elevation model elevation data of the monitoring area and the height of the permanent scatterer point;

[0013] The permanent scatterer points are classified into near-ground permanent scatterer points and non-near-ground permanent scatterer points according to their relative elevations;

[0014] Removing non-near-ground permanent scatterer points in the permanent scatterer point classification to obtain near-ground permanent scatterer points;

[0015] The sedimentation rate index is calculated based on the deformation rate of the near-ground permanent scatterer points;

[0016] Calculate the heterogeneity index index and the gradient index index according to the sedimentation rate index;

[0017] Acquire a land development data set, and calculate a land development index indicator according to the land development data set;

[0018] Calculating a digital elevation model elevation index indicator according to the digital elevation model elevation data;

[0019] Constructing a risk factor set according to the sedimentation rate index, the uneven index index and the gradient index index;

[0020] A vulnerability factor set is constructed according to the land development index indicator and the digital elevation model elevation index indicator.

[0021] Optionally, the step of calculating the sedimentation rate index indicator according to the deformation rate of the near-ground permanent scatterer point specifically includes:

[0022] Obtain the actual study area and generate a rectangular grid according to the actual study area;

[0023] A subsidence rate index is generated based on the grid division of the rectangular grid and the deformation rate of the near-surface permanent scatterer points.

[0024] Optionally, the step of calculating the uneven index index and the gradient index index according to the sedimentation rate index specifically includes:

[0025] Generate a power spectrum of the sedimentation data according to the sedimentation rate index, and calculate the radial average data of the power spectrum of the sedimentation data;

[0026] Performing spectrum clipping according to radial average data of the power spectrum of the sedimentation data;

[0027] Inputting the high-frequency reconstructed component after the spectrum is trimmed into the local sedimentation variation index calculation model, and convolutionally calculating the uneven index index;

[0028] The low-frequency reconstruction component after the spectrum is clipped is input into the variation gradient index calculation model, and the gradient index index is calculated by convolution.

[0029] Optionally, the step of obtaining a land development dataset and calculating a land development index indicator according to the land development dataset specifically includes:

[0030] Obtain the land development index calculation model;

[0031] The land development data set is input into a land development index calculation model, and the land development index indicator is calculated according to the grid division.

[0032] Optionally, the step of calculating a digital elevation model elevation index indicator according to the digital elevation model elevation data specifically includes:

[0033] The digital elevation model data of the study area are resampled according to the grid division, and the digital elevation model elevation index index is calculated.

[0034] Optionally, the step of calculating the ground subsidence risk membership vector according to the hazard factor set and the vulnerability factor set specifically includes:

[0035] Calculate the weight of each indicator in the risk factor set and the vulnerability factor set according to the risk factor set and the vulnerability factor set;

[0036] Calculate the membership value of each indicator in the risk factor set and the vulnerability factor set according to the risk factor set and the vulnerability factor set;

[0037] The ground subsidence risk membership vector is calculated based on the weight of each indicator and the membership value of each indicator.

[0038] In addition, to achieve the above-mentioned purpose, the present application also proposes an urban land subsidence risk assessment device, which can implement the steps of the urban land subsidence risk assessment method described above.

[0039] In addition, to achieve the above-mentioned purpose, the present application also proposes a storage medium, which is a computer-readable storage medium, and a computer program is stored on the storage medium. When the computer program is executed by a processor, the steps of the urban land subsidence risk assessment method described above are implemented.

[0040] In addition, to achieve the above-mentioned purpose, the present application also proposes a computer program product, which includes a computer program, and when the computer program is executed by a processor, it implements the steps of the urban land subsidence risk assessment method as described above.

[0041] The present application discloses a method, device, storage medium and computer program product for assessing urban land subsidence risk, which relates to the field of wide-area land subsidence risk level assessment. The scheme obtains the height of permanent scatterer points, deformation rate of permanent scatterer points, land development data and digital elevation model elevation data in the monitoring area; constructs a set of hazard factors and a set of vulnerability factors according to the height of permanent scatterer points, deformation rate of permanent scatterer points, land development data and digital elevation model elevation data; calculates the land subsidence risk membership vector according to the set of hazard factors and the set of vulnerability factors; and assesses the land subsidence risk level according to the land subsidence risk membership vector. The method can realize the rapid assessment and visual expression of the urban wide-area land subsidence risk level, and provide a reliable basic support for the subsequent rapid grasp of the evolution of urban wide-area land subsidence risk and early warning of key land subsidence areas. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in 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 only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on the structures shown in these drawings without paying creative work.

[0043] Figure 1 This is a flow chart of the first embodiment of the urban land subsidence risk assessment method of the present application;

[0044] Figure 2 This is a flow chart of the second embodiment of the urban land subsidence risk assessment method of the present application;

[0045] Figure 3 This is a flow chart of the third embodiment of the urban land subsidence risk assessment method of the present application;

[0046] Figure 4 This is a flow chart of a fourth embodiment of the urban land subsidence risk assessment method of the present application;

[0047] Figure 5 This is a flow chart of the fifth embodiment of the urban land subsidence risk assessment method of the present application;

[0048] Figure 6 This is a schematic diagram of the distribution of settlement rate index for the urban land subsidence risk assessment method applied in this application;

[0049] Figure 7 This is a schematic diagram of the uneven index distribution of the urban land subsidence risk assessment method for this application;

[0050] Figure 8 This is a schematic diagram of the gradient index distribution of the urban land subsidence risk assessment method for this application;

[0051] Fig. 9 This is a schematic diagram of land development index distribution for the urban land subsidence risk assessment method applied in this application;

[0052] Fig.10 This is a schematic diagram of the distribution of elevation index of the digital elevation model for the urban land subsidence risk assessment method in this application;

[0053] Fig.11 This is a schematic diagram of the ground subsidence risk classification assessment results of the urban ground subsidence risk assessment method applied in this application.

[0054] The realization of the purpose, functional features and advantages of the present invention will be further explained in conjunction with embodiments and with reference to the accompanying drawings. DETAILED DESCRIPTION

[0055] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only 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.

[0056] It should be noted that if the embodiments of the present invention involve directional indications (such as up, down, left, right, front, back, etc.), the directional indications are only used to explain the relative position relationship, movement status, etc. between the components in a certain specific posture. If the specific posture changes, the directional indication will also change accordingly.

[0057] In addition, if there are descriptions involving "first", "second", etc. in the embodiments of the present invention, the descriptions of "first", "second", etc. are only used for descriptive purposes and cannot be understood as indicating or implying their relative importance or implicitly indicating the number of technical features indicated. Therefore, the features limited to "first" and "second" may explicitly or implicitly include at least one of the features. In addition, if "and / or" or "and / or" appears in the full text, its meaning includes three parallel solutions. Taking "A and / or B" as an example, it includes solution A, solution B, or solutions that satisfy both A and B. In addition, the technical solutions between the various embodiments can be combined with each other, but it must be based on the ability of ordinary technicians in this field to implement. When the combination of technical solutions is contradictory or cannot be implemented, it should be deemed that such combination of technical solutions does not exist and is not within the scope of protection required by the present invention.

[0058] The present application embodiment provides a method for assessing urban land subsidence risk, referring to Figure 1 , Figure 1 This is a flow chart of the first embodiment of the urban land subsidence risk assessment method of the present application.

[0059] In this embodiment, the urban land subsidence risk assessment method includes steps S10 to S40;

[0060] S10: Obtain the height of permanent scatterer points, deformation rate of permanent scatterer points, land development data and digital elevation model elevation data within the monitoring area;

[0061] S20: constructing a set of hazard factors and a set of vulnerability factors according to the permanent scatterer point height, the permanent scatterer point deformation rate, the land development data and the digital elevation model elevation data;

[0062] S30: Calculating a ground subsidence risk membership vector according to the hazard factor set and the vulnerability factor set.

[0063] S40: Evaluate the land subsidence risk level according to the land subsidence risk membership vector.

[0064] It is understood that a monitoring area refers to a geographical or logical area selected for specific monitoring activities.

[0065] It should be noted that permanent scatterer points are a special term in remote sensing, referring to various ground objects that have strong backscattering of radar waves and are relatively stable in time sequence. These targets dominate the radar pixel unit, and their echo signals have significant stability in the pixel unit. Specifically, permanent scatterer points may include artificial buildings, bridges, exposed rocks, artificially placed corner reflectors, etc. These ground objects maintain relatively unchanged scattering characteristics in SAR image sequences of multiple periods, so they are called permanent scatterers. In radar interferometry (InSAR) technology, especially permanent scatterer differential interferometry (PSInSAR), permanent scatterer points are used to construct mathematical models to effectively remove the influence of atmospheric factors, etc., so as to achieve high-precision monitoring of surface deformation. This technology is widely used to monitor the movement and changes of the earth's surface, such as deformation monitoring of urban buildings, bridges and dams.

[0066] It should be noted that the height of a permanent scatterer point (PS point) is not a directly defined concept, because a PS point itself refers to a ground object with stable backscattering characteristics for radar waves, such as buildings, exposed rocks, etc., and their height is usually determined by their physical characteristics and geographical location. However, in radar interferometry (InSAR) technology, the height information of PS points can be indirectly obtained through observation data from radar satellites.

[0067] Specifically, after the microwave pulses emitted by the radar satellite hit the ground, they will be reflected by the scatterers on the ground (including the PS point) to form echo signals. These echo signals contain information such as the position, shape, and height of the scatterers. By processing and analyzing these echo signals, the elevation information of the ground can be obtained, and then the height of the PS point relative to a certain reference surface can be obtained.

[0068] In general, although the height of permanent scatterer points is not directly defined, its height information can be indirectly obtained through observation data from radar satellites. In PSInSAR technology, by constructing a ground displacement model and optimizing the algorithm, the height change of PS points can be monitored and analyzed.

[0069] It can be understood that the deformation rate of a permanent scatterer point (PS point) refers to the deformation rate of these ground objects with stable scattering characteristics (such as buildings, exposed rocks, etc.) relative to a reference benchmark (such as a point on the earth's surface or the entire surface) within a certain time range. This deformation rate is usually measured in millimeters per year (mm / a).

[0070] In the PS-InSAR (Permanent Scatter Synthetic Aperture Radar Interferometry) technology, the phase information of the PS point can be obtained by interferometrically processing the SAR images of the same area at different time points. Since the PS points are stable in time sequence, their phase changes mainly reflect the surface deformation information. By unwrapping and analyzing these phase information, the deformation rate of the PS point can be calculated.

[0071] It should be noted that the elevation data in the digital elevation model (DEM) refers to the elevation value of the ground point, that is, the altitude of the terrain surface. These elevation data are represented in the form of an ordered numerical array (X, Y, Z), where X and Y represent the plane coordinates of the ground point, and Z represents the elevation value of the point (i.e., the altitude).

[0072] In DEM, elevation data is obtained through different methods, including satellite remote sensing, aerial photogrammetry and ground measurement, etc. These elevation data are used to build a digital model of the terrain surface, thereby realizing the digital expression of the terrain surface morphology.

[0073] It should be noted that the risk factor set mainly reflects the possibility and intensity of geological disasters.

[0074] It should be noted that the vulnerability factor set mainly reflects the losses and impacts that may be caused when geological disasters occur.

[0075] In this embodiment, the risk factor set includes a sedimentation rate index, a non-uniformity index, and a gradient index, and the vulnerability factor set includes a land development index and a digital elevation model elevation index. The following is a definition of each index:

[0076] Subsidence rate index: Subsidence rate refers to the speed of land settlement in underground projects, and is an important indicator for measuring soil deformation and foundation settlement. The subsidence rate index is an index for quantitative evaluation of this rate, usually calculated by continuous monitoring and measurement of land surface settlement. This index can reflect the speed and degree of land settlement, and is an important basis for assessing the risk of geological disasters.

[0077] Inhomogeneity index: The inhomogeneity index is an indicator that reflects the degree of difference or dispersion of the overall unit mark values. In geological hazard assessment, this indicator can be used to describe the inhomogeneity or variability of surface deformation. When the surface deformation of a certain area shows great inhomogeneity or variability, it may mean that the geological structure of the area is unstable and there is a high risk of geological hazards.

[0078] Gradient index indicator: In geological disaster assessment, this indicator can be compared to the measurement of surface deformation gradient or geological structure gradient. When the surface deformation gradient of a certain area is large, it may mean that the geological structure of the area has a large stress concentration or deformation trend, thereby increasing the risk of geological disasters.

[0079] Land development index: Land development index is an indicator used to reflect the degree of land development and utilization. It usually includes specific indicators such as land development area and land development investment, which are used to measure the development and utilization of land resources in a region. In geological disaster assessment, land development index can be used as a basis for assessing vulnerability. When the degree of land development in a certain area is high, it may mean that the building and population density in the area is high, thereby increasing the potential losses when geological disasters occur.

[0080] Digital elevation model elevation index: Digital elevation model (DEM) is a physical ground model that represents ground elevation in the form of a set of ordered numerical arrays (X, Y, Z). The digital elevation model elevation index is an index for quantitatively evaluating the elevation information in this model. This index can reflect the elevation changes and terrain characteristics of the ground, and is one of the important bases for assessing the vulnerability of geological disasters. When the elevation changes in a certain area are large or the terrain is complex, it may mean that the soil and geological structure in the area are relatively fragile and easily affected by geological disasters.

[0081] It should be noted that the land subsidence risk membership vector is a vector composed of multiple memberships, each of which represents the probability or degree of the land subsidence risk belonging to a specific risk level. This vector is usually used to quantitatively assess the size and possibility of land subsidence risk and provide a scientific basis for decision makers.

[0082] It should be noted that the land subsidence risk level is defined based on multiple factors such as the severity of land subsidence, the scope of impact, the development trend, and the possible harm caused. The specific level classification may vary depending on the standards of different regions and industries, and this embodiment does not limit this.

[0083] In this embodiment, the method for obtaining the height of permanent scatterer points and the deformation rate of permanent scatterer points in the monitoring area is to retrieve sufficient heavy-track SAR images of the study area, and use the PS-InSAR algorithm after operations such as spatiotemporal baseline screening and preprocessing and registration to calculate the height of permanent scatterer points and the deformation rate of permanent scatterer points in the monitoring area within the study area.

[0084] In specific implementation, the height of permanent scatterer points, deformation rate of permanent scatterer points and digital elevation model elevation data in the monitoring area are obtained; the relative elevation of each permanent scatterer point is calculated according to the digital elevation model elevation data and the height of the permanent scatterer points in the monitoring area; the permanent scatterer points are classified into near-ground permanent scatterer points and non-near-ground permanent scatterer points according to the relative elevation of the permanent scatterer points; the non-near-ground permanent scatterer points in the permanent scatterer point classification are removed to obtain the near-ground permanent scatterer points; a hazard factor set and a vulnerability factor set are constructed according to the deformation rate of the near-ground permanent scatterer points, the hazard factor set includes a settlement rate index index, a non-uniformity index index and a gradient index index, and the vulnerability factor set includes a land development index index and a digital elevation model elevation index index; the ground subsidence risk membership vector is calculated according to the hazard factor set and the vulnerability factor set; the ground subsidence risk level is evaluated according to the ground subsidence risk membership vector.

[0085] This embodiment obtains the ground subsidence risk membership vector by calculation, which can realize the rapid assessment and visual expression of the urban wide-area ground subsidence risk level, and provide a reliable basic support for the subsequent rapid grasp of the urban wide-area ground subsidence risk evolution and early warning of key ground subsidence areas.

[0086] Based on the first embodiment of the present application, in the second embodiment of the present application, the same or similar contents as those in the first embodiment can refer to the above introduction, and will not be repeated in the following. Figure 2 , Figure 2 This is a flow chart of the second embodiment of the urban land subsidence risk assessment method of the present application.

[0087] In this embodiment, the step S20 further includes steps S21 to S29;

[0088] S21: Calculate the relative elevation of each permanent scatterer point based on the digital elevation model elevation data of the monitoring area and the height of the permanent scatterer point;

[0089] S22: Classify the permanent scatterer points into near-ground permanent scatterer points and non-near-ground permanent scatterer points according to their relative elevations;

[0090] S23: removing non-near-ground permanent scatterer points in the permanent scatterer point classification to obtain near-ground permanent scatterer points;

[0091] S24: Calculate the sedimentation rate index based on the deformation rate of near-surface permanent scatterer points.

[0092] S25: Calculate the uneven index index and the gradient index index according to the sedimentation rate index.

[0093] S26: Acquire a land development data set, and calculate a land development index indicator according to the land development data set.

[0094] S27: Calculating a digital elevation model elevation index indicator according to the digital elevation model elevation data.

[0095] S28: constructing a risk factor set according to the sedimentation rate index, the uneven index index and the gradient index index.

[0096] S29: constructing a vulnerability factor set according to the land development index indicator and the digital elevation model elevation index indicator.

[0097] It should be noted that the relative elevation of a permanent scatterer point (PS point) refers to the vertical height difference of these stable scattering ground objects (such as buildings, rocks, etc.) relative to a selected reference point or reference surface.

[0098] In this embodiment, the calculation formula for the relative elevation of the permanent scatterer point is as follows:

[0099] Relative elevation = permanent scatterer point elevation - ground elevation

[0100] The ground elevation refers to the ground elevation value of the corresponding position of the DEM.

[0101] It can be understood that, since this embodiment only studies the near-ground data, it is necessary to filter the near-ground data and remove the non-near-ground data.

[0102] This embodiment uses the Kriging spatial interpolation method to calculate the ground elevation at the latitude and longitude coordinates of each PS point according to the DEM elevation data, and then calculates the relative elevation of each PS point in turn according to the regional DEM data and classifies the PS points into non-near-ground building PS points and near-ground PS points. The specific judgment formula is as follows:

[0103]

[0104] Among them, H PS ,D PS They represent the estimated height of the PS point and the result of Kriging spatial interpolation of the DEM model at the plane coordinates of the point. Similarly, H GCP ,D GCP It means the height of the origin or ground control point is solved and the DEM elevation interpolation result, and the measuring points located near the ground are selected as the research measuring point set.

[0105] It should be noted that the land development dataset is a dataset specifically used to describe and analyze land development activities. This dataset usually contains various information related to land development, such as land use type, land cover status, land development intensity, geographical location, etc.

[0106] Land use type: Information describing how land is used, such as agricultural land, industrial land, residential land, etc. This information helps to understand the distribution and utilization of land resources.

[0107] Land cover status: reflects the different types of cover on the earth's surface, such as forests, grasslands, farmlands, soil, glaciers, lakes, etc. Land cover data helps to assess the ecological value and environmental status of land resources.

[0108] Land development intensity: describes the degree and scope of land development activities, such as building density, volume ratio, etc. This information helps to understand the development potential and degree of land resources.

[0109] Geographic location: Provides geographical location information of land development activities, such as latitude and longitude, administrative divisions, etc. This information helps determine the spatial distribution and regional characteristics of land development activities.

[0110] The land development dataset has a wide range of data sources, including satellite remote sensing data, geographic information system (GIS) data, government statistics, etc. These data are collected, sorted and analyzed to form data resources that can be used for scientific research. This embodiment does not limit the method of obtaining the land development dataset.

[0111] In a specific implementation, the settlement rate index is calculated according to the deformation rate of the near-ground permanent scatterer point; the uneven index and the gradient index are calculated according to the settlement rate index; a land development data set is obtained, and the land development index is calculated according to the land development data set; the digital elevation model elevation index is calculated according to the digital elevation model elevation data. A risk factor set is constructed according to the settlement rate index, the uneven index and the gradient index; and a vulnerability factor set is constructed according to the land development index and the digital elevation model elevation index.

[0112] Based on the second embodiment of the present application, in the third embodiment of the present application, the same or similar contents as those of the second embodiment can be referred to the above introduction, and will not be repeated in the following. Figure 3 , Figure 3 This is a flow chart of the third embodiment of the urban land subsidence risk assessment method of the present application.

[0113] In this embodiment, the step S24 further includes steps S241 to S242; the step S52 further includes steps S251 to S254;

[0114] S241: Acquire the actual study area, and generate a rectangular grid according to the actual study area.

[0115] S242: Generate a sedimentation rate index according to the grid division of the rectangular grid and the deformation rate of the near-ground permanent scatterer points.

[0116] S251: Generate a power spectrum of the sedimentation data according to the sedimentation rate index, and calculate radial average data of the power spectrum of the sedimentation data.

[0117] S252: Perform spectrum clipping according to radial average data of the power spectrum of the sedimentation data.

[0118] S253: Inputting the high-frequency reconstructed component after the spectrum is trimmed into the local sedimentation variation index calculation model, and convolutionally calculating the uneven index indicator.

[0119] S254: Inputting the low-frequency reconstruction component after the spectrum is clipped into the variation gradient index calculation model, and convolutionally calculating the gradient index index.

[0120] In this embodiment, the actual research area refers to the area of ​​interest, also known as AOI. Its selection is mainly based on subjective demarcation by humans, and the specific method can be positioning by corner coordinates or drawing polygon vector frames using GIS software. After determining the research area, all relevant data need to be clipped to obtain data that matches the size and scope of the research area.

[0121] In this embodiment, the specific method of generating a rectangular grid is as shown in formula (2).

[0122]

[0123] in, Represents the x, y coordinates of point P (upper left corner of the study area) in the UTM coordinate system (unit: meter), Represents the x and y coordinates of the same point P in the latitude and longitude system (unit: degree). proj() represents the coordinate reprojection function. gridnum represents the above coordinates of the lower right corner of the study area. x , gridnum y Indicates the specific number of grids divided in the study area in the x and y directions, resolution x , resolution y Indicates the size of the grid (unit: meter).

[0124] Furthermore, the steps of generating a sedimentation rate index based on the grid division of the rectangular grid and the deformation rate of the near-ground permanent scatterer points are as follows: based on the geographic coordinates and the cell size, the study area is gridded and resampled using the Kriging interpolation method, and the resampled data is used as the sedimentation rate index in the risk factor system.

[0125] For example, the east-west length of the study area selected in this implementation case is about 2.7km, and the north-south length is about 2.6km. Combined with the general engineering practice of urban ground subsidence risk, and to facilitate subsequent data processing, the same number of grids is set in the north-south and east-west directions, which are 101 cells. The setting of this cell scale can not only meet the general requirements of ground subsidence risk assessment, but also ensure that each cell contains a sufficient number of InSAR monitoring points. Using the above-mentioned geographic coordinates and cell size as a benchmark, the study area can be divided into 101*101 grid areas, and the Kriging method is used to interpolate and resample the above-mentioned grids to generate subsidence rate indicators in the risk factor system.

[0126] In a specific implementation, an actual study area is obtained, and a rectangular grid is generated according to the actual study area; a settlement rate index is generated according to the grid division of the rectangular grid and the deformation rate of near-ground permanent scatterer points.

[0127] It should be noted that the power spectrum is the abbreviation of the power spectrum density function, which is defined as the signal power within a unit frequency band. It can be understood that the subsidence rate power spectrum represents the power of signals at different frequencies after the subsidence data is spectrally converted.

[0128] It can be understood that in the power spectrum, the frequency components in each direction are averaged to obtain radial average data, which helps to reduce the interference of directional noise and highlight the main frequency features.

[0129] It can be understood that one or more clipping thresholds are set on the power spectrum, which are used to distinguish between signals and noise, ensuring that only meaningful frequency components are retained.

[0130] As can be understood, the inverse Fourier transform is used to reconstruct the sedimentation signal based on the clipped power spectrum. This process can retain the key frequency components in the original signal while removing unnecessary noise.

[0131] It is understandable that by comparing the sedimentation data before and after clipping, it is possible to verify whether the clipping effect meets expectations. The signal smoothness, noise level, and retention of key frequency components can be checked.

[0132] It is understandable that the trimmed subsidence data is further analyzed to extract useful spectral features, which can be used for subsidence trend prediction, anomaly detection or other related applications.

[0133] Exemplarily, in this implementation case, the sedimentation rate index of 101*101 is regarded as a two-dimensional array, which can be converted into a complex domain spectrum array using the fast Fourier decomposition algorithm. The power spectrum can be further calculated based on the amplitude spectrum obtained by decomposition, thereby calculating its radial average data and setting a reasonable cutoff frequency threshold according to the calculation situation. In this case, the cutoff frequency threshold is set to 0.1Hz based on the actual calculation result, and the frequency can be converted into the subscript number corresponding to the frequency value under discrete coordinates, that is, after the spectrum is centered, the spectrum center is used as the origin, and the spectrum data within a radius of 4 is divided into a low-frequency interval. On the contrary, other spectrum areas under the ideal filter are high-frequency areas, thereby obtaining low-frequency and high-frequency component spectra.

[0134] In the specific implementation, the sedimentation rate index is used as input data, and a two-dimensional fast Fourier transform (fft2) is performed. At the same time, an ideal bandpass high and low frequency filter group is constructed and the high-frequency local component and the low-frequency overall component are reconstructed respectively. The radial average data of the power spectrum of the sedimentation data is calculated by intercepting the two-dimensional Fourier transform and setting the threshold to determine the cutoff frequency and complete the spectrum clipping. The filter group construction is shown as follows:

[0135]

[0136] Among them, highfilter and lowfilter represent the high-frequency part and low-frequency part of the sedimentation spectrum data (frequency domain), respectively, and f shift (x,y) represents the spectrum data after center frequency shift, {(x,y)|x 2 +y 2 >cut threshold}, x and y represent the row and column numbers of the spectrum two-dimensional matrix, respectively. threshold The cutoff frequency threshold is obtained by calculating the radial average data of the power spectrum. The high-frequency component is all the components in the spectrum that are higher than the cutoff frequency, and correspondingly, the low-frequency component is all the components in the spectrum that are lower than or equal to the cutoff frequency.

[0137] Furthermore, the high-frequency reconstruction component is input into the local settlement variation index calculation model, and the convolution calculation is used to obtain the uneven index index in the risk factor system; the low-frequency reconstruction component is input into the change gradient index calculation model, and the convolution calculation is used to obtain the gradient index index in the risk factor system.

[0138] The specific calculation formulas of the uneven index and gradient index are as follows:

[0139]

[0140] Among them, γ(x,y) and η(x,y) represent the sedimentation rate variation unevenness index and sedimentation rate gradient index at the coordinates, respectively. d(x,y) represents the sedimentation rate at the location, and d(x+w,y+w) represents the sedimentation rate in the neighborhood of the point. w can generally be 1 or 2, representing a 9-square grid or a 25-square grid with the cell as the center, obtained by taking the natural logarithm after taking the square sum. grad(x) and grad(y) represent the sedimentation rate gradient values ​​in the x and y directions at the cell, respectively. The gradient index is obtained by taking the root of the square sum.

[0141] Based on the third embodiment of the present application, in the fourth embodiment of the present application, the same or similar contents as those in the third embodiment can refer to the above description, and will not be repeated in the following. Figure 4 , Figure 4 This is a flow chart of the fourth embodiment of the urban land subsidence risk assessment method of the present application.

[0142] In this embodiment, the step S26 further includes steps S261 to S262; the step S27 further includes step S271.

[0143] S261: Obtaining land development index calculation model;

[0144] S262: Input the land development data set into a land development index calculation model, and calculate the land development index indicator according to grid division.

[0145] S271: Resample the digital elevation model data of the study area according to the grid division and calculate the digital elevation model elevation index indicator.

[0146] It is understandable that in order to ensure that the acquired data set is more reliable, the resolution of the land development data set in this embodiment should not be less than 3m. These data may include information such as land cover type, land use status, topography and other information, which is not limited in this embodiment.

[0147] It is understandable that according to the research purpose and data characteristics, a suitable land development index calculation model is selected, which is a method based on multi-criteria decision analysis (MCDA), machine learning or deep learning. This embodiment does not limit the selection of the land development index calculation model.

[0148] In this embodiment, the resampling uses a bilinear interpolation method, which is a commonly used image scaling and resampling method that calculates the value of a new pixel by weighted averaging four surrounding pixels.

[0149] For example, the open source land development dataset used in this implementation case has a resolution of 1m, and the new grid w=25 is determined according to the evaluation grid size, that is, a 51*51 grid in the land development dataset is 1 evaluation cell. The ratio of the number of cells classified as buildings or infrastructure to the number of cells that are not unclassified is used as the land development index of the evaluation grid. From this definition, it can be seen that the index has a normalized feature. Similarly, the DEM elevation data is re-divided according to the new grid size, and the resampled elevation interpolation data is used as the DEM elevation index of the corresponding cell.

[0150] In specific implementation, use external reliable land development or classification related data sets, the resolution of which should not be less than 3m. Input it into the land development index calculation model, and perform convolution calculation according to the divided grid size. The calculation formula is as follows to obtain the land development index index in the vulnerability factor system. Resample the DEM data of the study area according to the divided grid size, and the interpolation method can choose bilinear interpolation to obtain the DEM elevation index index in the vulnerability factor system.

[0151]

[0152] Based on the fourth embodiment of the present application, in the fifth embodiment of the present application, the same or similar contents as those in the fourth embodiment can refer to the above description, and will not be repeated in the following. Figure 5 , Figure 5 This is a flow chart of the fifth embodiment of the urban land subsidence risk assessment method of the present application.

[0153] In this embodiment, step S30 further includes steps S31 to S33.

[0154] S31: Calculating the weight of each indicator in the risk factor set and the vulnerability factor set according to the risk factor set and the vulnerability factor set.

[0155] S32: Calculating the membership value of each indicator in the risk factor set and the vulnerability factor set according to the risk factor set and the vulnerability factor set.

[0156] S33: Calculate the ground subsidence risk membership vector according to the weight of each indicator and the membership value of each indicator.

[0157] It is understandable that weight is a relative concept, which refers to the relative importance or influence of a certain indicator or factor in the overall evaluation or decision-making. In mathematics and statistics, weight is usually quantified as a numerical value that reflects the relative contribution or importance of each variable in the overall.

[0158] In an evaluation system with multiple indicators or factors, each indicator or factor may have a different impact on the final result. In order to quantify this impact, we assign a weight to each indicator or factor. The weight value range is usually 0 to 1 (or in some cases, any positive number, but the sum is 1 or a fixed value for normalization), where 0 means that the indicator or factor has no impact on the final result, and 1 (or the maximum value after normalization) means that the indicator or factor has the greatest impact on the final result.

[0159] It can be understood that the membership value is an important concept in fuzzy set theory, which is used to indicate the degree to which an element belongs to a fuzzy set. In land subsidence risk assessment or other similar fuzzy assessment systems, the membership value is usually quantified as a value between 0 and 1 to reflect the degree to which each indicator or factor belongs to a certain risk level.

[0160] It can be understood that the characteristics of the membership value are as follows:

[0161] Range limitation: The value range of membership is strictly limited to between 0 and 1. Among them, 0 means that the element does not belong to the fuzzy set at all, and 1 means that the element completely belongs to the fuzzy set.

[0162] Continuous variation: The degree of membership can take any real number between 0 and 1, which allows elements to belong to the fuzzy set to a certain degree, rather than simply being divided into either belonging or not belonging.

[0163] Relativity: Membership is relative to a specific fuzzy set. The same element may have different membership values ​​in different fuzzy sets.

[0164] It is understandable that in land subsidence risk assessment, the determination of weights is crucial for calculating the risk membership vector. By reasonably allocating the weights of each indicator in the hazard factor set and the vulnerability factor set, the size and possible level of land subsidence risk can be more accurately assessed.

[0165] It can be understood that the risk membership vector is a mathematical tool used to quantify the degree of subordination of risk factors to different risk levels or evaluation criteria. The risk membership vector means that in the process of risk assessment or evaluation, for a certain risk factor, its degree of subordination to each risk level or evaluation standard is quantified as a set of values. The vector composed of these values ​​is the membership vector of the risk factor. The risk membership vector is usually composed of multiple elements, each element corresponds to a risk level or evaluation standard, and the value range of the element is between 0 and 1. Specifically, if the risk level is divided into n levels, then for a certain risk factor, its membership vector can be expressed as:

[0166] R=(r1,r2,r3,...,rn);

[0167] Among them, ri (i = 1, 2, ..., n) represents the degree to which the risk factor belongs to the i-th risk level and satisfies the following conditions:

[0168] 0≤ri≤1, and Σri=1 (in some cases, Σri may not be 1, but represents the relative size of each membership value).

[0169] It should be noted that the weight of each indicator may be calculated using an entropy weight method or an empirical expert method, which is not limited in this embodiment.

[0170] In this embodiment, if the entropy weight method is used to calculate the weight of each indicator, the specific formula is as follows:

[0171]

[0172] The matrix Q is a standardized matrix obtained after each of the z factors is standardized. The number of data for each factor is n, and the q in the matrix represents the standardized value. j That is, it represents the information entropy of the j-th element, w j Represents the weight of the j-th element after information entropy allocation.

[0173] It can be understood that the membership degree represents the fuzzy relationship between the evaluation criteria and the factors. The common trapezoidal function quantification method is used to determine the membership value and construct a suitable membership matrix. The general form of the membership matrix of the land subsidence risk assessment index of the three hazard factors and the two vulnerability factors in the present invention can be expressed as follows:

[0174]

[0175] where h 11 ~h 33 Indicates the membership value of the three indicators in the risk factor set, v 11 ~v 23 It represents the membership value of the two indicators in the vulnerability factor set. After determining the membership function and setting the interval, the risk assessment membership vector is calculated for the assessment factor set data after each cell is aligned. The specific calculation formula is as follows:

[0176]

[0177] Among them, R1~R3 ​​represent the membership values ​​of the cell as low risk, medium risk and high risk after fuzzy comprehensive evaluation. According to the maximum membership principle, the risk level corresponding to the maximum value is taken as the evaluation result. By evaluating each cell in turn, the distribution of ground subsidence risk in the entire study area can be obtained.

[0178] In order to quickly complete the risk assessment, this application can also use the expert method to determine the weight of each assessment indicator in combination with the actual project.

[0179] For example, the weight set is set to W = {ω1, ω2, ω3, ω4, ω5} = {0.4, 0.175, 0.225, 0.1, 0.1}, and the comment set is V = {low risk, medium risk, high risk}. Therefore, the evaluation membership function is selected as a trapezoidal fuzzy distribution, and the specific expression is as follows, where a, b, c, d are the set key thresholds:

[0180]

[0181] Taking the calculation of the fuzziness of the sedimentation rate index as an example, the above formulas give the membership function expressions of the low risk, medium risk and high risk of the index respectively. In this case, the thresholds a, b, c, and d of the index are set to 0, 1, 2, 4, and 8 respectively. Therefore, for the index results calculated for the cell, the above functions can be used to calculate the risk membership values ​​and normalize them. The thresholds for the normalized variation unevenness index, gradient index, and land development index are set to 0, 0.4, 0.5, 0.7, and 0.95. The corresponding thresholds for the DEM elevation index are 0, 3, 9, 15, and 30. According to the specific data of each evaluation index of the cell, the above-mentioned piecewise functions and thresholds can be referred to for calculation in turn to obtain the membership matrix of each evaluation index of the cell:

[0182]

[0183] After completing the setting of the evaluation factor weight set W and the calculation of the membership matrix, the fuzzy comprehensive evaluation membership calculation can be performed. In order to better display the calculation results, this case selects the actual data of a cell in the study area for display. The calculation process of the remaining cells is exactly the same and will not be shown in detail here. The specific data and formula of this cell are as follows:

[0184]

[0185] This means that after the fuzzy comprehensive evaluation, the cell is classified as low risk, medium risk, and high risk at 0.80, 0.05, and 0.15, respectively. Based on the maximum membership evaluation principle, the low risk with a membership of 0.80 is selected as the ground subsidence risk assessment level of the cell. The subsidence risk of all cells can be assessed according to the above calculation and evaluation ideas.

[0186] In order to more specifically demonstrate the practical effect of the present invention in practical application, Figure 6 This is a schematic diagram of the settlement rate index distribution of the urban land subsidence risk assessment method in this application, where the darker the color in the area, the more obvious the settlement, and the lighter the color, the more obvious the uplift; Figure 7 This is a schematic diagram of the uneven index distribution of the urban land subsidence risk assessment method in this application, where the brighter the color in the area, the more serious the uneven deformation in the area; Figure 8 This is a schematic diagram of the gradient index distribution of the urban land subsidence risk assessment method in this application, where the brighter the color in the area, the greater the gradient of change in the area; Fig. 9 This is a schematic diagram of the land development index distribution of the urban land subsidence risk assessment method for this application, where the brighter the color in the area, the greater the proportion of the area used for construction and municipal facilities development, that is, the greater the load on the ground in the area; Fig.10 This is a schematic diagram of the distribution of the digital elevation model elevation index of the urban land subsidence risk assessment method in this application, where the brighter the color in the area, the higher the digital model elevation value (height); Fig.11 This is a schematic diagram of the ground subsidence risk grading assessment results of the urban ground subsidence risk assessment method of the present application, wherein the ground subsidence risk grading assessment results schematic diagram includes three colors: black, white and gray. Black represents that the area is at a low risk level, gray represents that the area is at a medium risk level, and white represents that the area is at a high risk level.

[0187] It can be seen from the assessment results that medium- and high-risk areas are basically concentrated in potential danger areas such as high settlement rates and large gradients. The comprehensive assessment results can objectively represent the potential settlement danger areas in the study area, which has guiding significance for the actual wide-area settlement assessment and analysis application.

[0188] In specific implementation, the weight of each indicator in the hazard factor set and the vulnerability factor set is calculated based on the hazard factor set and the vulnerability factor set; the membership value of each indicator in the hazard factor set and the vulnerability factor set is calculated based on the hazard factor set and the vulnerability factor set; the ground subsidence risk membership vector is calculated based on the weight of each indicator and the membership value of each indicator.

[0189] In addition, to achieve the above-mentioned purpose, the present application also proposes a ground subsidence risk assessment device, which can implement the steps of the urban ground subsidence risk assessment method described above.

[0190] In addition, to achieve the above-mentioned purpose, the present application also proposes a storage medium, which is a computer-readable storage medium, and a computer program is stored on the storage medium. When the computer program is executed by a processor, the steps of the urban land subsidence risk assessment method described above are implemented.

[0191] The present application also provides a computer program product, including a computer program, which implements the steps of the urban land subsidence risk assessment method as described above when executed by a processor.

[0192] The above description is only for the present invention and does not limit the patent scope of the present invention. All equivalent structural changes made by using the contents of the present invention specification and drawings under the concept of the present invention, or directly / indirectly applied in other related technical fields are included in the patent protection scope of the present invention.

Claims

1. A method for assessing urban land subsidence risk, characterized in that: The urban land subsidence risk assessment method comprises: Obtain the permanent scatterer point height, permanent scatterer point deformation rate, land development data and digital elevation model elevation data within the monitoring area; Constructing a set of hazard factors and a set of vulnerability factors according to the permanent scatterer point heights, permanent scatterer point deformation rates, land development data, and digital elevation model elevation data; Calculating a ground subsidence risk membership vector according to the hazard factor set and the vulnerability factor set; The ground subsidence risk level is evaluated according to the ground subsidence risk membership vector.

2. The urban land subsidence risk assessment method according to claim 1, characterized in that: The step of constructing a set of hazard factors and a set of vulnerability factors according to the point deformation rate of near-ground permanent scatterers specifically includes: Calculate the relative elevation of each permanent scatterer point based on the digital elevation model elevation data of the monitoring area and the height of the permanent scatterer point; The permanent scatterer points are classified into near-ground permanent scatterer points and non-near-ground permanent scatterer points according to their relative elevations; Removing non-near-ground permanent scatterer points in the permanent scatterer point classification to obtain near-ground permanent scatterer points; The sedimentation rate index is calculated based on the deformation rate of the near-ground permanent scatterer points; Calculate the heterogeneity index index and the gradient index index according to the sedimentation rate index; Calculating a land development index indicator based on the land development data set; Calculating a digital elevation model elevation index indicator according to the digital elevation model elevation data; Constructing a risk factor set according to the sedimentation rate index, the uneven index index and the gradient index index; A vulnerability factor set is constructed according to the land development index indicator and the digital elevation model elevation index indicator.

3. The urban land subsidence risk assessment method according to claim 2, characterized in that: The step of calculating the subsidence rate index according to the point deformation rate of the near-ground permanent scatterer specifically includes: Obtain the actual study area and generate a rectangular grid according to the actual study area; A subsidence rate index is generated based on the grid division of the rectangular grid and the deformation rate of the near-surface permanent scatterer points.

4. The urban land subsidence risk assessment method according to claim 3, characterized in that: The step of calculating the uneven index index and the gradient index index according to the sedimentation rate index specifically includes: Generate a power spectrum of the sedimentation data according to the sedimentation rate index, and calculate the radial average data of the power spectrum of the sedimentation data; Performing spectrum clipping according to radial average data of the power spectrum of the sedimentation data; Inputting the high-frequency reconstructed component after the spectrum is trimmed into the local sedimentation variation index calculation model, and convolutionally calculating the uneven index index; The low-frequency reconstruction component after the spectrum is clipped is input into the variation gradient index calculation model, and the gradient index index is calculated by convolution.

5. The urban land subsidence risk assessment method according to claim 4, characterized in that: The step of obtaining a land development data set and calculating a land development index indicator according to the land development data set specifically includes: Obtain the land development index calculation model; The land development data set is input into a land development index calculation model, and the land development index indicator is calculated according to the grid division.

6. The urban land subsidence risk assessment method according to claim 5, characterized in that: The step of calculating the digital elevation model elevation index indicator according to the digital elevation model elevation data specifically includes: The digital elevation model data of the study area are resampled according to the grid division, and the digital elevation model elevation index index is calculated.

7. The urban land subsidence risk assessment method according to claim 6, characterized in that: The step of calculating the ground subsidence risk membership vector according to the hazard factor set and the vulnerability factor set specifically includes: Calculate the weight of each indicator in the risk factor set and the vulnerability factor set according to the risk factor set and the vulnerability factor set; Calculate the membership value of each indicator in the risk factor set and the vulnerability factor set according to the risk factor set and the vulnerability factor set; The ground subsidence risk membership vector is calculated based on the weight of each indicator and the membership value of each indicator.

8. An urban land subsidence risk assessment device, characterized in that: The urban land subsidence risk assessment device can implement the steps of the urban land subsidence risk assessment method described in any one of claims 1-7.

9. A storage medium, characterized in that: The storage medium is a computer-readable storage medium, and a computer program is stored on the storage medium. When the computer program is executed by a processor, the steps of the urban land subsidence risk assessment method according to any one of claims 1 to 7 are implemented.

10. A computer program product, characterized in that The computer program product comprises a computer program, and when the computer program is executed by a processor, the steps of the urban land subsidence risk assessment method according to any one of claims 1 to 7 are implemented.

Images