Ground subsidence and disaster risk assessment method, system and equipment based on three-dimensional deformation

Abstract

The application discloses a ground subsidence and disaster risk assessment method, system and equipment based on three-dimensional deformation, relates to the field of ground subsidence disaster prevention, and comprises the following steps: acquiring a synthetic aperture radar image set, underground water level dynamic characteristics and global navigation satellite system reference station data of a research area; determining an unwrapping graph queue by adopting a small baseline subset differential interferometric synthetic aperture radar method based on coherence on the basis of the synthetic aperture radar image set; estimating vertical and horizontal rates according to the unwrapping graph queue and the global navigation satellite system reference station data; selecting high-coherence points and decomposing two-dimensional deformation on direction data, further determining surface displacement of strata in the research area, and describing ground subsidence intensity of the research area by adopting subsidence difference and angular deformation; and according to the surface displacement of strata in the research area and the underground water level dynamic characteristics, the disaster risk of ground subsidence in the research area is evaluated. The application improves the accuracy of ground subsidence and disaster risk assessment.

Description

Technical Field

[0001] The present application relates to the field of ground subsidence disaster prevention, and in particular to a ground subsidence and disaster risk assessment method, system and equipment based on three-dimensional deformation. Background Art

[0002] With socioeconomic and technological development, large areas of flat alluvial plains composed of loose sediments in arid regions have been transformed from desert grasslands to cultivated land. This expansion has intensified the demand for groundwater resources and led to their overexploitation. Long-term overexploitation of groundwater can cause aquifer levels to drop, reducing pore water pressure in the aquifer itself and in the relative aquiclude, leading to soil compaction and ultimately severe land subsidence. This can harm surrounding communities and infrastructure, and cause permanent compaction of the aquifer, reducing its storage capacity. Therefore, measuring surface deformation and controlling and mitigating the hazards associated with land subsidence are crucial to supporting groundwater management and maintaining sustainable groundwater extraction activities.

[0003] Multi-temporal Interferometric Synthetic Aperture Radar (MT-InSAR) technology can measure surface deformation with millimeter-per-year accuracy. However, the spatial distribution and intensity of ground deformation are influenced by complex factors, including groundwater level drawdown, clay layer thickness, tectonic activity, and building loads, all of which contribute to land subsidence. The dual hydrogeological structure of inland arid oases, where single-structured phreatic aquifers transition to phreatic and confined aquifers, further exacerbates the uncertainty in the spatiotemporal evolution of land subsidence and groundwater overexploitation. Currently, the response mechanism of surface deformation in arid oases to ongoing aquifer compaction is unclear, and a systematic assessment of the hazards induced by surface deformation is lacking. Furthermore, Synthetic Aperture Radar (SAR) images of oases in arid regions suffer from severe decoherence, low coherent point density, and difficulty separating noise, resulting in large errors in terrain deformation inversion in such areas.

[0004] Therefore, it is urgent to solve the problems of severe coherence in deconfusion areas and quantification of surface deformation in arid areas, reveal the characteristics of stratum deformation during the development of irrigated agriculture, and clarify the response mechanism of the surface to the continuous compaction of aquifers, so as to achieve high-precision assessment of ground subsidence disaster risks. Summary of the Invention

[0005] The purpose of this application is to provide a three-dimensional deformation-based ground subsidence and disaster risk assessment method, system and equipment, which can achieve high-precision assessment of ground subsidence and disaster risks.

[0006] To achieve the above objectives, this application provides the following solutions:

[0007] In a first aspect, the present application provides a method for land subsidence and disaster risk assessment based on three-dimensional deformation, comprising:

[0008] Acquire a synthetic aperture radar image set of the study area, groundwater level dynamic characteristics, and global navigation satellite system reference station data; the images in the synthetic aperture radar image set include ascending and descending orbit data;

[0009] Based on the synthetic aperture radar image set, a disentanglement map queue is determined using a coherence-based small baseline subset differential interferometry synthetic aperture radar method;

[0010] estimating vertical velocity and horizontal velocity according to the unwrapping graph queue and the global navigation satellite system reference station data to obtain direction data;

[0011] Performing high-coherence point selection and two-dimensional deformation decomposition on the directional data to obtain the radar line-of-sight direction displacement;

[0012] Determining the surface displacement of the strata in the study area according to the radar line of sight displacement;

[0013] According to the surface displacement of the strata in the study area, the ground settlement intensity in the study area is described using settlement difference and angular deformation;

[0014] Based on the surface displacement of the strata in the study area and the dynamic characteristics of the groundwater level, the disaster risk of ground subsidence in the study area is assessed.

[0015] In a second aspect, the present application provides a three-dimensional deformation-based land subsidence and disaster risk assessment system, comprising:

[0016] A data acquisition module is used to obtain a synthetic aperture radar image set of the study area, groundwater level dynamic characteristics, and global navigation satellite system reference station data; the images in the synthetic aperture radar image set include ascending orbit data and descending orbit data;

[0017] a phase unwrapping module for determining an unwrapping map queue based on the synthetic aperture radar image set by using a coherence-based small baseline subset differential interferometry synthetic aperture radar method;

[0018] a rate estimation module, configured to estimate vertical rate and horizontal rate based on the unwrapping graph queue and the GNSS reference station data to obtain direction data;

[0019] A displacement decomposition module is used to select high-coherence points and perform two-dimensional deformation decomposition on the directional data to obtain the radar line of sight direction displacement;

[0020] a displacement determination module, configured to determine the surface displacement of the strata in the study area according to the displacement in the radar sight line direction;

[0021] an intensity determination module for describing the ground subsidence intensity of the study area using settlement differences and angular deformations according to the surface displacement of the strata in the study area;

[0022] The risk assessment module is used to assess the disaster risk of ground subsidence in the study area based on the surface displacement of the strata in the study area and the dynamic characteristics of the groundwater level.

[0023] In a third aspect, the present application provides a computer device comprising: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the above-mentioned three-dimensional deformation-based ground subsidence and disaster risk assessment method.

[0024] According to the specific embodiments provided in this application, this application has the following technical effects:

[0025] The present application provides a method, system and equipment for ground subsidence and disaster risk assessment based on three-dimensional deformation. Based on a synthetic aperture radar image set, a small baseline subset differential interferometry synthetic aperture radar method based on coherence is adopted to determine an unwrapped map queue, which can effectively extract three-dimensional deformation information of the ground at different time points, and then construct a high-precision three-dimensional deformation field. The directional data is further subjected to high-coherence point selection and two-dimensional deformation decomposition to filter out noise sources, and the radar line of sight direction displacement is obtained and the surface displacement of the stratum in the study area is determined. The ground subsidence intensity in the study area is described by using settlement difference and angular deformation. The disaster risk of ground subsidence in the study area is assessed according to the surface displacement of the stratum in the study area and the dynamic characteristics of the groundwater level, thereby improving the assessment accuracy of ground subsidence and disaster risk. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0027] Figure 1 This is an application environment diagram of a ground subsidence and disaster risk assessment method based on three-dimensional deformation in one embodiment of the present application;

[0028] Figure 2 A flowchart of a method for land subsidence and disaster risk assessment based on three-dimensional deformation provided in one embodiment of the present application;

[0029] Figure 3A process flowchart of a method for land subsidence and disaster risk assessment based on three-dimensional deformation provided in one embodiment of the present application;

[0030] Figure 4 This is a schematic diagram of the relationship matrix of hazard levels in one embodiment of the present application;

[0031] Figure 5 A schematic diagram of the functional modules of a three-dimensional deformation-based land subsidence and disaster risk assessment system provided in one embodiment of the present application;

[0032] Figure 6 A schematic diagram of the structure of a computer device provided in one embodiment of the present application. DETAILED DESCRIPTION

[0033] The following will be combined with the drawings in the embodiments of this application to clearly and completely describe the technical solutions in the embodiments of this application. Obviously, the embodiments described are only part of the embodiments of this application, not all of the embodiments. Based on the embodiments in this application, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of this application.

[0034] In order to make the above-mentioned purposes, features and advantages of the present application more obvious and easy to understand, the present application is further described in detail below with reference to the accompanying drawings and specific implementation methods.

[0035] The land subsidence and disaster risk assessment method based on three-dimensional deformation provided in the embodiment of the present application can be applied to Figure 1 In the application environment shown. Among them, the terminal 102 communicates with the server 104 through the network. The data storage system can store the data that the server 104 needs to process. The data storage system can be set up separately, integrated on the server 104, or placed on the cloud or other servers. The terminal 102 can send the synthetic aperture radar image set, groundwater level dynamic characteristics and global navigation satellite system base station data of the study area to the server 104. After receiving the synthetic aperture radar image set, groundwater level dynamic characteristics and global navigation satellite system base station data of the study area, the server 104 uses settlement difference and angle deformation to describe the ground settlement intensity of the study area, and evaluates the disaster risk of ground settlement in the study area. The server 104 can feedback the ground settlement intensity and disaster risk assessment results of the study area to the terminal 102. In addition, in some embodiments, the ground settlement and disaster risk assessment method based on three-dimensional deformation can also be implemented separately by the server 104 or the terminal 102.

[0036] Terminal 102 may include, but is not limited to, various desktop computers, laptops, smartphones, tablet computers, IoT devices, and portable wearable devices. IoT devices may include smart speakers, smart TVs, smart air conditioners, and smart car devices. Portable wearable devices may include smart watches, smart bracelets, and head-mounted devices. Server 104 may be implemented as a standalone server or a server cluster consisting of multiple servers, or may be a cloud server.

[0037] In an exemplary embodiment, Figure 2 and Figure 3 As shown, a method for land subsidence and disaster risk assessment based on three-dimensional deformation is provided. The method is executed by a computer device, specifically a computer device such as a terminal or a server, or a terminal and a server. In the embodiment of the present application, the method is applied to Figure 1 The server 104 in the example is used for explanation, including the following steps 201 to 207.

[0038] Step 201: Acquire a synthetic aperture radar image set, groundwater level dynamic characteristics, and Global Navigation Satellite System (GNSS) base station data of the study area. The images in the synthetic aperture radar image set include ascending and descending orbit data.

[0039] In one exemplary embodiment, a synthetic aperture radar image set, groundwater level observation data, historical annual runoff, groundwater extraction, groundwater level drawdown data, and global navigation satellite system reference station data are obtained for the study area. Groundwater level dynamic characteristics are determined based on the synthetic aperture radar image set, the groundwater level observation data, the historical annual runoff, the groundwater extraction, and the groundwater level drawdown data.

[0040] This application obtains time-series SAR images of the study area using ascending and descending Sentinel-1A and Sentinel-1S. SAR images of the study area are used to obtain the spatiotemporal characteristics of surface deformation in the study area. The dynamic characteristics of the groundwater level are analyzed by combining groundwater level observation wells and historical annual runoff, groundwater extraction volume, and groundwater level drop. GNSS base station data are also collected as verification data.

[0041] Among them, Sentinel-1's ascending and descending SAR images are obtained in the Interferometric Wideswath (IW) scanning mode, with a total scan length of 250 kilometers. It consists of three sub-scan modes: IW1, IW2 and IW3, each of which consists of a series of pulse trains.

[0042] Step 202 : Based on the synthetic aperture radar image set, a disentanglement map queue is determined by using a coherence-based small baseline subset (SBAS) differential interferometry synthetic aperture radar method.

[0043] In an exemplary embodiment, a coherence-based SBAS differential InSAR method is used to obtain an SBAS dataset, which is a time series of surface deformation, i.e., a disentangled map sequence. Step 202 includes the following steps (21) to (25):

[0044] (21) Selecting a reference image from the synthetic aperture radar image set.

[0045] (22) Using the reference image to geometrically register other images in the synthetic aperture radar image set, and resampling the geometrically registered images to the spatial reference system of the reference image to obtain a resampled image set.

[0046] Geometric registration satisfies the following conditions:

[0047] ①For the VV polarization mode level 1 synthetic space radar scenario.

[0048] ② Phase noise was mitigated using the azimuth and range multifocal tools in the Environment for Visualizing Images (ENVI) remote sensing image processing platform, and the terrain phase component was subtracted using the 30-meter resolution digital surface model of the Shuttle Radar Topography Mission (SRTM).

[0049] ③ The burst interferogram is corrected using the enhanced spectral diversity method, especially the enhanced spectral diversity (ESD) estimation by range and azimuth within the overlapping region of adjacent and consecutive bursts.

[0050] ④ Perform interferometric phase compensation by removing the phase slope. For a dual-channel scene, the differential phase should be the incident angle λ / 4π, and the line-of-sight displacement spacing is 75m.

[0051] (23) Based on the resampled image set, a multi-view interference map is generated according to a preset time baseline threshold and a preset space baseline threshold.

[0052] (24) The multi-view interferogram is filtered using a Goldstein filter to obtain a filtered interferogram. The interferogram is filtered using a Goldstein filter to enhance the coherence of the InSAR time series interferogram.

[0053] (25) The filtered interference graph is phase unwrapped using a minimum cost flow to generate an unwrapped graph queue with a smaller spatial baseline and a shorter time interval.

[0054] Step 203 : estimating vertical velocity and horizontal velocity according to the unwrapping graph queue and the GNSS reference station data to obtain direction data.

[0055] In an exemplary embodiment, since the dewrapped graph queue has phase noise and decoherence problems, the present application performs post-processing and estimates the vertical rate and horizontal rate. Step 203 includes the following steps (31) to (34):

[0056] (31) Low spatial frequency signals are removed from the disentangled map queue, and the radar line of sight (LOS) velocity of each coherent target is converted into a temporary vertical direction.

[0057] Specifically, the low spatial frequency signals due to uncompensated orbit errors and atmospheric effects are removed from the unwrapped map queue (i.e., SBAS dataset), and the LOS velocity of each coherent target is converted to a temporary vertical direction using the following formula;

[0058]

[0059] Among them, V Uj is the vertical velocity of the coherent target j, V LOSj is the coherent target j in the vertical direction of LOS, U LOSj is the cosine of the coherent target j in the horizontal direction of LOS.

[0060] (32) Adjusting the position of each coherent target according to the global navigation satellite system base station data to obtain a small baseline subset point data set.

[0061] Specifically, combined with the GNSS base station data, the following formula is used to recalculate the t of each coherent target: n The original position of the time series at the initial time t0 is adjusted for displacement:

[0062] d Uj (t n )=d Uj' (t n )-V U' *{t n -t0};

[0063] Among them, t n is the current time, t0 is the initial time, d Uj (t n) is the current position of the coherent target j, d Uj' (t n ) is the original position of the coherent target j, V U' Vertical velocity offset resulting from reference point homogenization to remove low spatial frequency trends.

[0064] (33) Resampling and interpolating the small baseline subset point data set to obtain a plurality of grid elements.

[0065] Specifically, the SBAS point datasets were initially resampled using a grid with a spacing of 100 meters to describe the inherent differences in the distribution and location of targets in the two geometries. Nearest neighbor interpolation was then used to transfer the attributes of each dataset to the grid. In particular, no interpolation was performed in areas without targets.

[0066] (34) Performing cosine vector inversion on each grid element to obtain an east-west velocity component, a north-south velocity component, and a vertical velocity component. The directional data includes the east-west velocity component, the north-south velocity component, and the vertical velocity component.

[0067] Specifically, for each grid element i, the direction cosines of the LOS of Sentinel-1 IW in the east-west, north-south, and vertical directions are considered in each acquisition mode, and the following formula is used to invert the two-variable linear equations containing three unknown variables to obtain the east-west velocity component, the north-south velocity component, and the vertical velocity component:

[0068] V Ai =E Ai *V Ei +N Ai *V Ni +U Ai *V Ui ;

[0069] V Di =E Di *V Ei +N Di *V Ni +U Di *V Ui ;

[0070] Among them, V Ai is the annual velocity of the IW band of the ascending orbit at grid element i, E Ai is the known value of the east-west cosine of the IW band of the ascending orbit at grid element i, N Ai is the known value of the north-south direction cosine of the IW band of the ascending orbit at grid element i, U Ai is the known value of the vertical direction cosine of the IW band of the ascending orbit at grid element i, V Diis the annual velocity of the descending orbit IW band at each grid element i, E Di is the known value of the east-west direction cosine of the descending orbit IW band at each grid element i, N Di is the known value of the north-south direction cosine of the descending orbit IW band at each grid element i, U Di is the known value of the vertical direction cosine of the descending orbit IW band at each grid element i, V Ei is the east-west velocity component, V Ni is the north-south velocity component, V Ui is the vertical velocity component.

[0071] Step 204 : performing high coherence point selection and two-dimensional deformation decomposition on the directional data to obtain LOS displacement.

[0072] In an exemplary embodiment, the directional data is imported into the Mintpy open source software, and noise is filtered out using high coherence point selection and two-dimensional deformation decomposition. After eliminating various noise sources, linear regression is applied to the resolved displacement time series to obtain the LOS displacement. Step 204 includes the following steps (41) to (44):

[0073] (41) Construct pixel stability index and coherence index:

[0074]

[0075]

[0076] Where X is the pixel stability index, Υ is the coherence index, N is the number of available interferograms, and φ n is the phase value of the pixel of interest, is the mean of the phase values ​​at multiple moments at the pixel of interest, J() is the Bessel function, z1 and z2 represent the magnitudes of the two negative coherences, and the superscript * represents the complex conjugate.

[0077] (42) The directional data is imported into the Mintpy open source software, and high coherence points are selected based on the pixel stability index and the coherence index to obtain an interferometric synthetic aperture radar image.

[0078] Specifically, pixels whose frequency distributions of pixel stability index and coherence index are both in the top 80% are regarded as high coherence points.

[0079] (43) The three-dimensional surface deformation of the interferometric synthetic aperture radar image in the LOS direction is projected onto a two-dimensional plane to obtain a two-dimensional surface deformation time series.

[0080] Specifically, by constructing a multi-dimensional spatial projection transformation relationship, the three-dimensional surface deformation in the LOS direction of the InSAR image is projected onto a two-dimensional plane.

[0081] (44) Taking the spatial coverage of the ascending orbit data as a reference and combining it with the descending orbit data of the corresponding date, the LOS displacement of the two-dimensional surface deformation time series is decomposed into the vertical direction and the east-west direction to obtain the LOS displacement in the vertical direction and the LOS displacement in the east-west direction.

[0082] Specifically, the following formulas are used to determine the LOS displacement in the vertical direction and the LOS displacement in the east-west direction:

[0083] d U cosθ-d E cosαsinθ+d N sinαsinθ=d LOS +δ LOS ;

[0084]

[0085] Among them, d U is the LOS displacement in the vertical direction, d E is the LOS displacement in the east-west direction, d N is the LOS displacement in the north direction, α is the heading angle of the satellite, θ is the incidence angle of the satellite, and d LOS is the observed LOS direction position, δ LOS is the observed LOS direction error. The superscript A represents the parameter obtained from the ascending data, and the superscript D represents the parameter obtained from the descending data.

[0086] Step 205: Determine the surface displacement of the stratum in the study area based on the LOS displacement.

[0087] In an exemplary embodiment, step 205 includes the following steps (51) to (52):

[0088] (51) The strata in the study area are considered as homogeneous, isotropic, and flat cylinders in an elastic half-space. The total displacement of a point due to the contraction of the flat cylinder is obtained by superimposing the displacements of countless infinitesimally small strain units.

[0089] (52) Based on the flat cylinder and the LOS displacement, the surface displacement of the strata in the study area is determined using the following formula:

[0090]

[0091] Where L is the surface displacement of the stratum in the study area, represents the surface displacement in the vertical component, v is the Poisson's ratio, Δh is the changing thickness of the flat cylinder, R is the radius of the flat cylinder, r is the center distance of the flat cylinder, D is the depth of the flat cylinder, J1() is the zero-order Bessel function, J2() is the first-order Bessel function, and α' is the variable used for the integration operator, which represents the LOS displacement.

[0092] Step 206 : Based on the surface displacement of the strata in the study area, the ground settlement intensity in the study area is described using settlement difference and angular deformation.

[0093] In an exemplary embodiment, the relationship between the settlement difference and the angular deformation is:

[0094]

[0095] Where β is the angular deformation, Δd ui is based on the settlement difference between two points, i.e. vertical displacement, l res is the image resampling resolution of the LOS displacement, which is determined according to the surface displacement of the strata in the study area.

[0096] Step 207 : Evaluate the disaster risk of land subsidence in the study area based on the surface displacement of the strata in the study area and the dynamic characteristics of the groundwater level.

[0097] In an exemplary embodiment, based on the sedimentation gradient, combined with the exposure and vulnerability information of the risk factors, the hazard level is divided using the optimal threshold of the empirical statistical model. The risk factors include population size, transportation infrastructure, housing density, and groundwater level dynamic characteristics. Step 207 includes the following steps (71) to (75):

[0098] (71) Determine multiple subsidence gradients based on the surface displacement of the strata in the study area. Specifically, the different subsidence gradients are divided according to the surface displacement l of the strata in the study area, and the total vertical displacement Δd of each 100-meter grid cell compared to its adjacent cells is calculated as ui The maximum percentage of

[0099] (72) Obtain the population, transportation infrastructure, and housing density of the study area.

[0100] (73) The dynamic characteristics of groundwater level, population, transportation infrastructure and housing density in the study area were weighted to obtain a comprehensive risk index.

[0101] (74) The support vector machine (SVM) classification model was used, with the sedimentation gradient and comprehensive risk index as input variables and the hazard level as the output variable, to find the optimal threshold for dividing the hazard level.

[0102] (75) Based on the multiple subsidence gradients and the comprehensive risk index, multiple hazard levels are divided according to the optimal threshold to assess the disaster risk of ground subsidence in the study area. Figure 4 As shown in the figure, the hazard levels are specifically divided into R0, R1, R2, R3 and R4.

[0103] This application acquires a surface deformation dataset through time-series InSAR analysis. High-coherence point selection and two-dimensional deformation decomposition are used to filter out noise. Linear regression is applied to the analyzed displacement time series to obtain LOS displacements. High-coherence point extraction and LOS displacement decomposition are used to filter out noise sources. The Geertsma analytical method is used to model surface displacements in elastic media, and the intensity of ground subsidence is described by settlement differences and angular deformation. Based on the settlement gradient and combined with the exposure and vulnerability information of risk factors, the optimal threshold of the SVM classification model is used to classify the hazard into five categories: R0, R1, R2, R3, and R4, thereby achieving a risk assessment of the hazard of ground subsidence. The three-dimensional deformation field obtained through SAR image analysis is used to model ground subsidence in the study area and conduct a disaster risk assessment, which has certain ecological and social value.

[0104] In summary, this application has the following advantages:

[0105] (1) The SBAS differential InSAR method is used to process multi-time synthetic aperture radar data, which can effectively extract the three-dimensional deformation information of the ground at different time points, including vertical settlement and horizontal displacement, and then construct a high-precision three-dimensional deformation field.

[0106] (2) Aiming at the problem of severe deterioration of SAR image interferometric phase quality (decoherence) in oasis areas with intense agricultural activities, the pixel stability index and coherence index were combined to extract high coherent points. Based on the coherent point density and GNSS records, the threshold for coherent point selection was adjusted, and a coherent point selection method suitable for oases in arid areas was constructed. This method solved the problems of severe decoherence, low coherent point density, and difficulty in separating noise in oasis areas, and enabled effective surface deformation monitoring of oases in arid areas.

[0107] (3) Using formation deformation information as a constraint, we modeled the surface displacement of elastic media caused by changes in underground reservoir pore pressure and inverted hydrogeological parameters. We innovatively integrated the groundwater level drop and clay layer thickness to assess the accuracy of the inversion parameters and constructed an accurate surface displacement model based on settlement differences and angular deformation.

[0108] Based on the same inventive concept, the embodiments of the present application also provide a three-dimensional deformation-based ground subsidence and disaster risk assessment system for implementing the three-dimensional deformation-based ground subsidence and disaster risk assessment method mentioned above. The implementation solution provided by this system is similar to the implementation solution described in the above-mentioned method. Therefore, the specific limitations in one or more embodiments of the three-dimensional deformation-based ground subsidence and disaster risk assessment system provided below can be found in the limitations of the three-dimensional deformation-based ground subsidence and disaster risk assessment method above, and will not be repeated here.

[0109] In an exemplary embodiment, Figure 5 As shown, a ground subsidence and disaster risk assessment system based on three-dimensional deformation is provided, including: a data acquisition module 501, a phase unwrapping module 502, a rate estimation module 503, a displacement decomposition module 504, a displacement determination module 505, an intensity determination module 506 and a risk assessment module 507.

[0110] The data acquisition module 501 is used to acquire a synthetic aperture radar image set of the study area, groundwater level dynamic characteristics, and global navigation satellite system reference station data. The images in the synthetic aperture radar image set include ascending orbit data and descending orbit data.

[0111] The phase unwrapping module 502 is configured to determine an unwrapping map sequence based on the synthetic aperture radar image set by adopting a coherence-based small baseline subset differential interferometry synthetic aperture radar method.

[0112] The velocity estimation module 503 is configured to estimate the vertical velocity and the horizontal velocity according to the unwrapping graph queue and the GNSS reference station data to obtain direction data.

[0113] The displacement decomposition module 504 is used to perform high coherence point selection and two-dimensional deformation decomposition on the direction data to obtain LOS displacement.

[0114] The displacement determination module 505 is used to determine the surface displacement of the stratum in the study area according to the LOS displacement.

[0115] The intensity determination module 506 is used to describe the ground subsidence intensity of the study area using settlement difference and angular deformation according to the surface displacement of the strata in the study area.

[0116] The risk assessment module 507 is used to assess the disaster risk of land subsidence in the study area based on the surface displacement of the strata in the study area and the dynamic characteristics of the groundwater level.

[0117] In an exemplary embodiment, a computer device is provided. The computer device may be a server or a terminal. The internal structure diagram thereof may be as follows: Figure 6As shown. The computer device includes a processor, a memory, an input / output interface (Input / Output, I / O) and a communication interface. The processor, the memory and the input / output interface are connected via a system bus, and the communication interface is connected to the system bus via the input / output interface. The processor of the computer device is used to provide computing and control capabilities. The memory of the computer device includes a non-volatile storage medium and an internal memory. The non-volatile storage medium stores an operating system, a computer program and a database. The internal memory provides an environment for the operation of the operating system and the computer program in the non-volatile storage medium. The database of the computer device is used to store a set of synthetic aperture radar images of the study area, dynamic characteristics of groundwater levels and global navigation satellite system base station data. The input / output interface of the computer device is used to exchange information between the processor and an external device. The communication interface of the computer device is used to communicate with an external terminal via a network connection. When the computer program is executed by the processor, a ground subsidence and disaster risk assessment method based on three-dimensional deformation is implemented.

[0118] Those skilled in the art will understand that Figure 6 The structure shown in the figure is only a block diagram of a part of the structure related to the solution of the present application, and does not constitute a limitation on the computer device to which the solution of the present application is applied. The specific computer device may include more or fewer components than shown in the figure, or combine certain components, or have a different component arrangement.

[0119] In an exemplary embodiment, a computer-readable storage medium is provided, storing a computer program. When the computer program is executed by a processor, the steps in the above-mentioned method embodiments are implemented.

[0120] In an exemplary embodiment, a computer program product is provided, including a computer program. When the computer program is executed by a processor, the steps in the above method embodiments are implemented.

[0121] It should be noted that the user information (including but not limited to user device information, user personal information, etc.) and data (including but not limited to data used for analysis, stored data, displayed data, etc.) involved in this application are all information and data authorized by the user or fully authorized by all parties, and the collection, use and processing of relevant data must comply with relevant regulations.

[0122] In this application, all actions to obtain signals, information or data are carried out in compliance with the relevant data protection laws and policies of the country where they are located and with the authorization given by the owner of the corresponding device.

[0123] Those skilled in the art will understand that all or part of the processes in the above-mentioned embodiment methods can be implemented by instructing the relevant hardware through a computer program, and the computer program can be stored in a non-volatile computer-readable storage medium. When the computer program is executed, it can include the processes of the embodiments of the above-mentioned methods. Among them, any reference to memory, database or other media used in the embodiments provided in this application may include at least one of non-volatile and volatile memory. Non-volatile memory may include read-only memory (ROM), magnetic tape, floppy disk, flash memory, optical memory, high-density embedded non-volatile memory, resistive random access memory (ReRAM), magnetic random access memory (MRAM), ferroelectric random access memory (FRAM), phase change memory (PCM), graphene memory, etc. Volatile memory may include random access memory (RAM) or external cache memory, etc. By way of illustration and not limitation, RAM may be in various forms, such as static random access memory (SRAM) or dynamic random access memory (DRAM).

[0124] The databases involved in the various embodiments provided herein may include at least one of a relational database and a non-relational database. Non-relational databases may include, but are not limited to, distributed databases based on blockchains. The processors involved in the various embodiments provided herein may include, but are not limited to, general-purpose processors, central processing units, graphics processing units, digital signal processors, programmable logic units, data processing logic units based on quantum computing, and the like.

[0125] The technical features of the above embodiments can be combined arbitrarily. To make the description concise, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0126] This document uses specific examples to illustrate the principles and implementation methods of this application. The description of the above examples is only intended to help understand the method and core concept of this application. At the same time, for those skilled in the art, based on the concept of this application, there may be changes in the specific implementation methods and application scope. In summary, the content of this specification should not be understood as limiting this application.

Claims

1. A land subsidence and disaster risk assessment method based on three-dimensional deformation, characterized in that: The three-dimensional deformation-based land subsidence and disaster risk assessment method includes: Acquire a synthetic aperture radar image set of the study area, groundwater level dynamic characteristics, and global navigation satellite system reference station data; the images in the synthetic aperture radar image set include ascending and descending orbit data; Based on the synthetic aperture radar image set, a disentanglement map queue is determined using a coherence-based small baseline subset differential interferometry synthetic aperture radar method; estimating vertical velocity and horizontal velocity according to the unwrapping graph queue and the global navigation satellite system reference station data to obtain direction data; Performing high-coherence point selection and two-dimensional deformation decomposition on the directional data to obtain the radar line-of-sight direction displacement; Determine the surface displacement of the strata in the study area based on the radar line of sight displacement, specifically including: The strata in the study area are considered as flat cylinders in a homogeneous, isotropic, and elastic half-space; Based on the displacement of the flat cylinder and the radar line of sight, the surface displacement of the stratum in the study area is determined using the following formula: Where l is the surface displacement of the stratum in the study area, v is the Poisson's ratio, Δh is the thickness change of the flat cylinder, R is the radius of the flat cylinder, r is the center distance of the flat cylinder, D is the depth of the flat cylinder, J1( ) is the zero-order Bessel function, J2( ) is the first-order Bessel function, and α' is the variable used for the integration operator, which represents the displacement in the radar line of sight direction; According to the surface displacement of the strata in the study area, the ground settlement intensity in the study area is described using settlement difference and angular deformation; Based on the surface displacement of the strata in the study area and the dynamic characteristics of the groundwater level, the disaster risk of ground subsidence in the study area is assessed.

2. The method for land subsidence and disaster risk assessment based on three-dimensional deformation according to claim 1, characterized in that: Acquire synthetic aperture radar imagery, groundwater level dynamics, and GNSS base station data for the study area, including: Acquire synthetic aperture radar imagery, groundwater level observation data, historical annual runoff, groundwater extraction, groundwater level drawdown data, and global navigation satellite system reference station data for the study area; The dynamic characteristics of the groundwater level are determined based on the synthetic aperture radar image set, the groundwater level observation data, the historical annual runoff, the groundwater extraction volume and the groundwater level drop data.

3. The method for land subsidence and disaster risk assessment based on three-dimensional deformation according to claim 1, characterized in that: Based on the synthetic aperture radar image set, a disentanglement map queue is determined using a coherence-based small baseline subset differential interferometry synthetic aperture radar method, specifically comprising: selecting a reference image from the synthetic aperture radar image set; performing geometric registration on other images in the synthetic aperture radar image set using the reference image, and resampling the geometrically registered images to a spatial reference system of the reference image to obtain a resampled image set; Based on the resampled image set, generating a multi-view interferogram according to a preset time baseline threshold and a preset space baseline threshold; The multi-view interferogram is filtered using a Goldstein filter to obtain a filtered interferogram; The filtered interference graph is phase unwrapped using a minimum cost flow to generate an unwrapped graph queue.

4. The method for land subsidence and disaster risk assessment based on three-dimensional deformation according to claim 1, characterized in that: Estimating vertical velocity and horizontal velocity according to the unwrapping graph queue and the global navigation satellite system reference station data to obtain direction data, specifically including: Eliminating low spatial frequency signals from the disentangled image queue and converting the radar line-of-sight velocity of each coherent target into a temporary vertical direction; Adjusting the position of each coherent target according to the global navigation satellite system reference station data to obtain a small baseline subset point data set; resampling and interpolating the small baseline subset point data set to obtain a plurality of grid elements; Cosine vector inversion is performed on each grid element to obtain an east-west velocity component, a north-south velocity component, and a vertical velocity component; the direction data includes the east-west velocity component, the north-south velocity component, and the vertical velocity component.

5. The method for land subsidence and disaster risk assessment based on three-dimensional deformation according to claim 1, characterized in that: The directional data is subjected to high-coherence point selection and two-dimensional deformation decomposition to obtain the radar line-of-sight direction displacement, specifically including: Construct pixel stability index and coherence index; Importing the directional data into Mintpy open source software, and selecting high coherence points based on the pixel stability index and the coherence index to obtain an interferometric synthetic aperture radar image; Projecting the three-dimensional surface deformation of the interferometric synthetic aperture radar image in the radar line of sight onto a two-dimensional plane to obtain a two-dimensional surface deformation time series; Taking the spatial coverage of the ascending orbit data as a reference and combining it with the descending orbit data of the corresponding date, the radar line of sight displacement of the two-dimensional surface deformation time series is decomposed into the vertical direction and the east-west direction to obtain the radar line of sight displacement.

6. The method for land subsidence and disaster risk assessment based on three-dimensional deformation according to claim 1, characterized in that: The relationship between settlement difference and angular deformation is: Where β is the angular deformation, Δd ui is based on the difference in settlement between two points, l res is the image resampling resolution of the radar line of sight displacement, which is determined according to the surface displacement of the strata in the study area.

7. The method for land subsidence and disaster risk assessment based on three-dimensional deformation according to claim 1, characterized in that: Based on the surface displacement of the strata in the study area and the dynamic characteristics of the groundwater level, the disaster risk of land subsidence in the study area is assessed, including: determining a plurality of subsidence gradients based on surface displacements of strata in the study area; Obtain the population, transportation infrastructure, and housing density of the study area; The groundwater level dynamic characteristics, population, transportation infrastructure and housing density of the study area were weighted to obtain a comprehensive risk index; A support vector machine classification model was used, with the settlement gradient and comprehensive risk index as input variables and the hazard level as output variable, to find the optimal threshold for hazard level classification. According to the multiple settlement gradients and the comprehensive risk index, multiple hazard levels are divided according to the optimal threshold to evaluate the disaster risk of ground subsidence in the study area.

8. A land subsidence and disaster risk assessment system based on three-dimensional deformation, applied to the land subsidence and disaster risk assessment method based on three-dimensional deformation according to any one of claims 1 to 7, characterized in that: The three-dimensional deformation-based land subsidence and disaster risk assessment system includes: A data acquisition module is used to obtain a synthetic aperture radar image set of the study area, groundwater level dynamic characteristics, and global navigation satellite system reference station data; the images in the synthetic aperture radar image set include ascending orbit data and descending orbit data; a phase unwrapping module for determining an unwrapping map queue based on the synthetic aperture radar image set by using a coherence-based small baseline subset differential interferometry synthetic aperture radar method; a rate estimation module, configured to estimate vertical rate and horizontal rate based on the unwrapping graph queue and the GNSS reference station data to obtain direction data; A displacement decomposition module is used to select high-coherence points and perform two-dimensional deformation decomposition on the directional data to obtain the radar line of sight direction displacement; a displacement determination module, configured to determine the surface displacement of the strata in the study area according to the displacement in the radar sight line direction; an intensity determination module for describing the ground subsidence intensity of the study area using settlement differences and angular deformations according to the surface displacement of the strata in the study area; The risk assessment module is used to assess the disaster risk of ground subsidence in the study area based on the surface displacement of the strata in the study area and the dynamic characteristics of the groundwater level.

9. A computer device comprising: A memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the three-dimensional deformation-based ground subsidence and disaster risk assessment method according to any one of claims 1 to 7.

Images