Surface deformation risk assessment method based on insar and multi-source data
By combining InSAR with multi-source data and physical mechanism models, the problems of causal analysis and trend prediction in surface deformation risk assessment were solved, enabling dynamic quantification of risk and scientific decision support.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- CHONGQING THREE GORGES UNIV
- Filing Date
- 2026-02-11
- Publication Date
- 2026-06-05
Smart Images
Figure CN122153575A_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of geological disaster monitoring and risk assessment technology, specifically a method for assessing surface deformation risk based on InSAR and multi-source data. Background Technology
[0002] Landslides, ground subsidence, and surface uplift pose serious threats to infrastructure, the ecological environment, and community safety. In recent years, Synthetic Aperture Radar Interferometry (InSAR) technology, especially temporal InSAR techniques (such as Small Baseline Subset SBAS-InSAR), has become a core tool in surface deformation monitoring due to its ability to acquire millimeter-level temporal data of surface deformation over a wide area with high precision and non-contact operation. However, single deformation monitoring data cannot directly reveal the underlying physical mechanisms and causes of deformation, nor can it quantitatively assess future risks. Therefore, how to integrate multi-source data to move from "deformation phenomenon monitoring" to "deformation mechanism analysis and risk assessment" has become a key research challenge.
[0003] To improve the intelligence level of deformation analysis, existing technologies are beginning to combine InSAR technology with artificial intelligence methods. For example, Chinese patent document CN120993412A discloses "A method and system for detecting and classifying surface deformation by integrating CNN and SBAS-InSAR technologies." This method first uses SBAS-InSAR technology to acquire the surface deformation velocity field and combines it with multi-source data such as optical images and topography. Second, it performs preliminary classification of deformation areas (such as stable, subsidence, landslide, etc.) using simple threshold rules based on deformation velocity and slope. Finally, it uses these preliminary classification labels as supervision signals to train a convolutional neural network model to achieve more refined automatic identification of deformation types.
[0004] While this method has made progress in the automated classification of deformation phenomena, it still has the following significant shortcomings: the method can only classify deformation phenomena (such as landslides and subsidence) without deeply analyzing the physical mechanisms that cause deformation (such as groundwater and engineering activities), and cannot answer the questions of "why deformation occurs" and "how it will develop in the future". Furthermore, its output is a discrete category map, which cannot quantitatively classify and rank risks, making it difficult to support accurate risk management decisions. The process stops at classification and map generation, lacking mechanism verification, trend prediction, and result validation based on physical models, thus limiting the reliability and practicality of the results.
[0005] In summary, current technologies have achieved high-precision monitoring and intelligent type identification of surface deformation. However, significant technological gaps remain in the transition to deformation causal mechanism analysis, trend prediction, and quantitative risk assessment. Therefore, an innovative methodology is urgently needed. This methodology should deeply integrate time-series InSAR deformation data with multi-source disaster-causing factor data, systematically combining five key aspects: data monitoring, mechanism analysis, model validation, trend prediction, and quantitative rating. This would construct a complete and coherent technological chain, fundamentally elevating the understanding from "deformation description" to "risk cognition," and providing a scientific basis for precise geological disaster prevention and control.
[0006] The information disclosed in the background section is only intended to enhance the understanding of the background of this disclosure, and therefore may include information that does not constitute prior art known to those skilled in the art. Summary of the Invention
[0007] The purpose of this invention is to provide a method for assessing the risk of land deformation based on InSAR and multi-source data, so as to solve the problems mentioned in the background art.
[0008] To achieve the above objectives, the present invention provides the following technical solution: A method for assessing surface deformation risk based on InSAR and multi-source data, specifically including: A method for assessing surface deformation risk based on InSAR and multi-source data, characterized by the following specific steps: Step 1: Acquire synthetic aperture radar images covering the area to be evaluated at equal time intervals to form a monitoring image set. Perform time-series InSAR processing on the monitoring image set to obtain the time-series deformation data of each pixel location and establish the mapping relationship between pixels and geographic coordinates. Step 2: Divide the area to be evaluated into several grid cells. For each grid cell, calculate the average deformation rate and spatial standard deviation of the deformation rate based on the temporal deformation data of all pixels within it. Step 3: Based on the preset settlement threshold and non-uniformity threshold, mark the grids with average deformation rate exceeding the settlement threshold and spatial standard deviation of deformation rate exceeding the non-uniformity threshold as high-risk areas in the initial screening. Step 4: For each high-risk area in the initial screening, obtain multi-source disaster-causing factor data for its corresponding range, and conduct collaborative analysis with the temporal deformation data of the area to infer the deformation cause based on physical mechanisms and predict its development trend. Step 5: Based on the causes and trends of deformation, and combined with the average deformation rate of the high-risk areas in the initial screening, assess the risk level of each high-risk area in the initial screening, and uniformly assign the lowest risk level representing stability to the remaining grid units that were not marked as high-risk areas in the initial screening, thereby obtaining the risk level of all grid units.
[0009] Furthermore, the synthetic aperture radar images within the monitoring image set are distributed at equal intervals in time, and the various synthetic aperture radar images within the monitoring image set correspond to each other spatially. The steps for obtaining the temporal deformation data of each pixel location include: The earliest synthetic aperture radar image is set as the initial image. Each synthetic aperture radar image in the monitoring image set is paired with the initial image to form several image pairs. The acquisition time of the synthetic aperture radar image is marked as the time tag of the image pair. Differential interferometry is performed on each image pair to obtain a differential interferometric phase map and phase unwrapping is performed to calculate the relative surface deformation map between each image pair. For the relative surface deformation map of each image pair, the phase difference of each pixel at each time tag is obtained. The phase difference of each pixel is sorted according to the time tag order to form the temporal deformation data of each pixel.
[0010] Furthermore, the mapping relationship between pixels and geographic coordinates is established, and the specific steps include: Obtain the geographic coordinates of the top-left pixel in the synthetic aperture radar image. The geographic coordinates of the bottom right pixel And the total number of rows of pixels inside the synthetic aperture radar image. Total number of columns ; For synthetic aperture radar images located at the first line, number The latitude of a column of pixels is calculated using the following linear interpolation formula. With longitude : in, and These are the row index and column index of the cell, respectively. The value range is 0 to , The value range is 0 to ; The row and column indices of each pixel in the synthetic aperture radar image can be obtained using the above calculation formula. Its geographic coordinates The mapping relationship between them.
[0011] Furthermore, the area to be evaluated is divided into regular grid cells. Specific steps include: Determine the scale of the grid cells, and use the smallest inscribed rectangle of the region to be evaluated as the starting range for grid division, setting the lower left corner of the rectangle as the starting origin for grid division; Starting from the origin, the grid is translated at equal intervals along the horizontal and vertical directions according to the grid scale to generate regularly arranged grid boundary lines in sequence. Several regular grid cells are formed by adjacent horizontal and vertical grid boundary lines, from which only grid cells that have effective overlap with the actual spatial range of the area to be evaluated are retained; Traverse all cells with acquired temporal deformation data and determine whether they fall within the spatial boundary of a valid grid cell based on their geographic coordinates. If a pixel falls into a certain grid cell, then the temporal deformation data corresponding to that pixel is associated with that grid cell.
[0012] Furthermore, for each grid cell, based on the temporal deformation data of all pixels within it, the average deformation rate and spatial standard deviation of the deformation rate of that grid cell are calculated, including: For each grid cell, extract the temporal deformation data of all its associated cells; Using the observation time corresponding to the time-series deformation data as the independent variable and the value of the time-series deformation data as the dependent variable, a linear fit is performed. The slope of the fitted line is determined as the deformation rate value of the pixel, and its unit is length per time. Calculate the arithmetic mean of the deformation rates of all pixels within the grid cell, and use it as the average deformation rate of the grid cell. Calculate the sample standard deviation of the deformation rate values of all pixels within the grid cell, and use it as the spatial standard deviation of the deformation rate of the grid cell.
[0013] Furthermore, based on the calculated average deformation rate and spatial standard deviation of deformation rate for each unit grid, and compared with the preset preliminary risk screening threshold, grids with an average deformation rate exceeding the settlement threshold and a spatial standard deviation of deformation rate exceeding the non-uniformity threshold are marked as high-risk areas in the initial screening. Specific steps include: Set the settlement threshold and non-uniformity threshold used for judgment; For each grid cell, perform the following judgment: When the average deformation rate of the grid cell exceeds the settlement threshold, it is determined that it meets the first risk condition; When the spatial standard deviation of the deformation rate of the grid cell exceeds the non-uniformity threshold, it is determined that it meets the second risk condition. For grid cells that simultaneously meet the first and second risk conditions mentioned above, they are marked as high-risk areas in the initial screening.
[0014] Furthermore, for each initially screened high-risk area, multi-source disaster-causing factor data for its corresponding range are obtained and analyzed in conjunction with the temporal deformation data of that area to infer the causes of deformation based on physical mechanisms and predict its development trend, including: For each high-risk area in the initial screening, the monitoring time series of multiple pre-defined disaster-causing factors within the area are obtained. The disaster-causing factors include physical mechanism disaster-causing factors, soil and rock time-effect disaster-causing factors, and groundwater-related disaster-causing factors. The monitoring time series of each category of disaster-causing factor consists of at least two monitoring sequences of specific physical quantities that reflect the disaster-causing mechanism of that category. Extract the temporal deformation data of all pixels in the region, calculate the average value of the deformation of all pixels at the same time point, and connect these average values in time sequence to generate a regional deformation time series that represents the overall deformation of the region. For each type of disaster-causing factor, calculate the Pearson correlation coefficient between each specific monitoring time series and the regional deformation time series to obtain the Pearson correlation coefficient subset for each category. For each type of disaster-causing factor, the one with the largest absolute value is selected from the subset of its corresponding Pearson correlation coefficients and used as the representative correlation coefficient between that type of disaster-causing factor and regional deformation. By comparing the absolute values of the correlation coefficients of all categories, the category with the largest absolute value is identified as the main disaster-causing factor category for the initially screened high-risk area; based on this main disaster-causing factor category, a preliminary hypothesis on the formation cause of the area is proposed.
[0015] Furthermore, for the high-risk areas identified in the initial screening, their development trends are predicted based on the results of the collaborative analysis and deformation data, including: Based on the regional deformation time series, the overall average deformation rate over the entire monitoring period is calculated as the first rate; Based on the data of the regional deformation time series within the most recent preset time period, its recent average deformation rate is calculated as the second rate; Compare the absolute value of the second speed with the absolute value of the first speed: If the absolute value of the second rate is greater than or equal to the preset first threshold, it is determined to be a fast trend; If the absolute value of the second rate is less than the first threshold and greater than or equal to the preset second threshold, it is determined to be in a trend. If the absolute value of the second rate is less than the second threshold, it is determined to be a slow trend; Wherein, the first threshold is greater than the second threshold.
[0016] Furthermore, based on the causes and trends of deformation, and combined with the average deformation rate of the initially screened high-risk areas, the risk level of each initially screened high-risk area is assessed. Specific steps include: Based on the final deformation causation conclusions output, three causal categories—physical mechanism-induced disaster, soil and rock age deformation, and groundwater-related disaster—are assigned causal weight values of 0.5, 0.3, and 0.2, respectively. Based on the deformation development trend prediction results, the development trend is characterized into three types: fast, medium and slow, and assigned trend coefficient values of 0.9, 0.6 and 0.3 respectively; The comprehensive risk index of the region is calculated as the product of the region's average deformation rate, causal weight value, and trend coefficient value. The comprehensive risk index is compared with a preset first comprehensive risk threshold and a second comprehensive risk threshold to determine the risk level, wherein the first comprehensive risk threshold is less than the second comprehensive risk threshold; If the overall risk index of the high-risk area in the initial screening is not greater than the first overall risk threshold, then the risk level of the corresponding grid unit is low risk level. If the comprehensive risk index of a high-risk area in the initial screening is greater than the first comprehensive risk threshold but not greater than the second comprehensive risk threshold, then the risk level of the corresponding grid unit is a certain risk level. If the comprehensive risk index of the initially screened high-risk area is greater than the second comprehensive risk threshold, then the risk level of the corresponding grid unit is a certain level, which is a high-risk level.
[0017] For all other grid cells that were not initially identified as high-risk areas, their risk levels are uniformly assigned the lowest risk level, representing stability.
[0018] Compared with the prior art, the beneficial effects of the present invention are: This invention achieves a comprehensive mechanistic analysis of surface deformation by synergistically analyzing time-series InSAR deformation data with monitoring sequences of multiple disaster-causing factors and establishing a correlation-based quantification model of the impact weights of these factors. This overcomes the limitations of traditional methods, which can only identify deformation types but cannot trace their causes. It enables the assessment process to accurately pinpoint the main physical factors leading to deformation (such as groundwater extraction and engineering activities), providing a direct, physically based scientific basis for subsequent targeted risk management.
[0019] This invention uses an integrated physical prediction model to extrapolate deformation trends and designs a set of quantitative rating rules that integrate deformation causes, trends, and average deformation rates, enabling dynamic quantification and classification of regional risks. This changes the situation where risk information in traditional risk maps is discrete and incomparable, allowing the risk levels of different regions to be measured and ranked through a unified comprehensive index. This provides clear and actionable decision support for prioritizing disaster prevention and mitigation and accurately allocating resources.
[0020] This invention significantly improves the reliability of risk assessment results and the robustness of the system by constructing a closed-loop process that includes physical model verification and logical consistency verification of the final results. This mechanism ensures that every step of reasoning from data to conclusions is verified, and the final output risk thematic map maintains logical consistency with the original data and intermediate conclusions, enhancing the interpretability of the technical results and their credibility and practical value in complex real-world engineering scenarios. Attached Figure Description
[0021] Figure 1 This is a functional model diagram of the present invention; Figure 2 This is a schematic diagram of the overall method flow of the present invention; Figure 3 This is a time-series co-analysis diagram of the average deformation of area A and the groundwater extraction volume in the example; Figure 4 This is a time-series co-analysis diagram of the average deformation and average groundwater level depth in area A of the embodiment. Figure 5 This is a time-series co-analysis diagram of the average deformation and crustal horizontal movement rate in area A of the embodiment; Figure 6 This is a time-series co-analysis diagram of the average deformation of region B and the construction load index in the embodiment. Figure 7 This is a time-series co-analysis diagram of the average deformation and soil volumetric moisture content in area B of the example; Figure 8 This is a time-series co-analysis diagram of the average deformation of area B and the amount of groundwater extraction in the example; Figure 9 This is a time-series co-analysis diagram of the average deformation and crustal horizontal movement rate in region C of the example; Figure 10 This is a time series analysis diagram of the regional average deformation and monthly seismic energy release in area C in the example; Figure 11 This is a time-series co-analysis diagram of the average deformation of area C and the construction load index in the example. Detailed Implementation
[0022] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to specific embodiments.
[0023] It should be noted that, unless otherwise defined, the technical or scientific terms used in this invention should have the ordinary meaning understood by one of ordinary skill in the art to which this invention pertains. The terms "first," "second," and similar terms used in this invention do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "comprising" or "including" mean that the element or object preceding the word encompasses the elements or objects listed following the word and their equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" are used only to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.
[0024] Example: Please see Figures 1 to 11 The present invention provides a technical solution: The method for assessing surface deformation risk based on InSAR and multi-source data includes the following steps: Step 1: Acquire synthetic aperture radar images covering the area to be evaluated at equal time intervals to form a monitoring image set. Perform time-series InSAR processing on the monitoring image set to obtain the time-series deformation data of each pixel location and establish the mapping relationship between pixels and geographic coordinates.
[0025] The setting of the equal time interval depends on the monitoring requirements and the surface deformation rate of the area to be evaluated. In this embodiment, the equal time interval is 14 days. Data is acquired using satellites with synthetic aperture radar imaging capabilities. Synthetic aperture radar images within the specified time range are downloaded through a relevant data service platform, ensuring that the images cover the area to be evaluated.
[0026] From the acquired synthetic aperture radar (SAR) imagery, images that meet the time interval requirements are selected to form a monitoring image set. This monitoring image set should include multiple images within the set time interval to ensure the accuracy of the time series analysis.
[0027] The earliest synthetic aperture radar (SAR) image in the monitoring image set is set as the initial image. Then, each other SAR image in the monitoring image set is paired with this initial image to form several interferometric image pairs. The acquisition time of each SAR image is tagged as the time label of the corresponding image pair.
[0028] For each pair of interferometric images, differential interferometric phase maps are generated using interferometry. These phase maps reflect the relative displacement of the Earth's surface at each observation time point relative to the initial time. Then, the generated differential interferometric phase maps are unwrapped to calculate the relative surface deformation map between each pair of images. For each relative surface deformation map, the phase difference of each pixel at that time label is extracted.
[0029] The phase difference of each pixel under all time tags is arranged in chronological order of the time tags to form the temporal deformation data of that pixel.
[0030] For each pair of synthetic aperture radar images acquired at adjacent times in the monitoring image set, differential interferometric phase maps are generated using interferometry. These phase maps reflect the relative displacement information of the Earth's surface between the two observation time points. Then, the generated differential interferometric phase maps are subjected to phase unwrapping processing to calculate the relative deformation map of the Earth's surface between each pair of adjacent time points.
[0031] Obtain the geographic coordinates of the top left starting point of the synthetic aperture radar image. Geographic coordinates in the lower right corner and the total number of rows of images Total number of columns ; For the image located at the first line, number The latitude of a column of pixels is calculated using the following linear interpolation formula. With longitude : in, and These are the row index and column index of the cell, respectively. The value range is 0 to , The value range is 0 to .
[0032] The geographic coordinates corresponding to each pixel in the image can be obtained by calculating using the above formula. For each pixel in the image, assign it a row number. , column number As a unique identifier, it is related to the calculated geographic coordinates. Associativity as a single record. Organize all such association records for all cells into a mapping table and store it.
[0033] Using the above mapping table, spatial vector points are created for the geographic coordinates of each pixel using a suitable GIS library. The geographic coordinates of each pixel are converted into spatial vector point data. The spatial vector point data is then associated and stored with the temporal deformation data of each pixel location to form a spatiotemporally integrated deformation dataset.
[0034] Step 2: Divide the area to be evaluated into several grid cells. For each grid cell, calculate the average deformation rate and spatial standard deviation of the deformation rate based on the temporal deformation data of all pixels within it.
[0035] Based on the need to balance monitoring accuracy and computational efficiency, the scale of the grid cells is determined. In this embodiment, it is 500 meters. A standard grid size of 500 meters is used. The specific grid generation process is as follows: Based on the spatial extent formed by all pixels in the area to be evaluated, calculate its minimum inscribed rectangle. Set the lower left vertex of this rectangle as the starting point for mesh generation. Starting from the origin, the grid is translated at equal intervals of 500 meters along both the horizontal and vertical directions to generate a series of grid boundary lines parallel to the coordinate axes. Several regular rectangular grid cells are formed by adjacent horizontal and vertical boundary lines. Then, spatial filtering is performed, retaining only grid cells that effectively overlap with the actual spatial extent of the area to be evaluated, and removing invalid cells completely outside the area, resulting in a set of valid grid cells. Iterate through each pixel record in the integrated data table, extracting its geographic coordinates, namely the longitude value Lon and the latitude value Lat. For each grid cell in the effective grid cell set, its spatial boundary is uniquely defined by the lower left corner coordinates (MinLon, MinLat) and the upper right corner coordinates (MaxLon, MaxLat). These coordinate values were calculated and stored during grid division based on the starting origin and grid scale. For the current pixel, check whether its coordinates satisfy the boundary conditions of a certain grid cell. The judgment logic is: if both MinLon and Lat are satisfied... Lon MaxLon and MinLat Lat MaxLat determines whether a cell falls within the spatial boundary of this grid cell. Once a cell is determined to fall within a valid grid cell, its unique identifier (row and column index) and its corresponding complete temporal deformation data in the table are added as a sub-record to the data structure of that grid cell. Each grid cell will ultimately associate all deformation temporal information of zero, one, or more cells, completing the spatial merging and data organization from continuous cell points to discrete grid cells.
[0036] After completing the spatial association of all pixels with grid cells, for each grid cell in the effective grid cell set, based on the temporal deformation data of all its associated pixels, the average deformation rate and spatial standard deviation of the deformation rate of that grid cell are calculated. The specific steps are as follows: For the grid cell to be calculated, obtain a list of all cells belonging to that cell from its associated data structure. Based on the unique identifier, i.e., the row and column index, of each cell in the list, extract its complete temporal deformation data sequence from the comprehensive data table; For each extracted pixel's temporal deformation data, the least squares linear fitting method is used for analysis to obtain the straight line that best represents the trend of deformation over time. The slope of this fitted line is determined as the deformation rate value of that pixel, in millimeters per year. A positive slope indicates an uplift trend, while a negative slope indicates a subsidence trend. Collect the deformation rate values of all pixels in the current grid cell, calculate the arithmetic mean of these values, and use it as the average deformation rate of the grid cell to characterize the overall average trend and intensity of deformation in the region. Based on all the above pixel deformation rate values and their calculated average values, the sample standard deviation of these rate values relative to the average value is calculated as the spatial standard deviation of the deformation rate of the grid cell. This standard deviation is used to quantify the spatial dispersion of the deformation rate within the grid cell; a larger value indicates that the deformation is more non-uniform in space. The calculated average deformation rate and spatial standard deviation of deformation rate for each grid cell are stored as attributes of that grid cell.
[0037] After calculating the average deformation rate of all grid cells, these grid cells are used as basic mapping units to spatially visualize their average deformation rate attribute values. Specifically, the average deformation rate value of each grid cell is rendered onto its corresponding geographic location using a preset color band (in this embodiment, a gradient from blue to red represents the rate change from uplift to subsidence), based on a pre-defined color scheme. This generates a spatial distribution map of the deformation rate field, using grid cells as units, covering the entire area to be evaluated. This map is stored along with the grid cell vector dataset as an important intermediate result for subsequent verification.
[0038] Step 3: Based on the preset settlement threshold and non-uniformity threshold, the grids with average deformation rate exceeding the settlement threshold and spatial standard deviation of deformation rate exceeding the non-uniformity threshold are marked as high-risk areas in the initial screening.
[0039] Based on the geological background, engineering safety standards, and historical deformation data of the area to be evaluated, key settlement and non-uniformity thresholds are set. The settlement threshold is used to identify grids with significant deformation rates. In this embodiment, based on the accuracy requirements and experience of regional settlement monitoring, this threshold is set to -10 mm / year (negative values indicate settlement). Grids with an average deformation rate lower than this value are considered to have a significant settlement trend. The non-uniformity threshold is used to identify grids where deformation is spatially unevenly distributed. In this embodiment, this threshold is set to 5 mm / year. Grids with a spatial standard deviation of deformation rate higher than this value are considered to have significant spatial differences in internal deformation.
[0040] Iterate through all valid mesh cells. For each mesh cell, read its stored average deformation rate and spatial standard deviation of deformation rate attributes, and then perform the following logical checks in sequence: Determine the first risk condition: Compare the average deformation rate of the grid with the preset settlement threshold. If the average deformation rate is less than the settlement threshold, the grid cell is determined to meet the first risk condition; otherwise, it is not met. Determine the second risk condition: Compare the spatial standard deviation of the deformation rate of the grid cell with the preset non-uniformity threshold. If the spatial standard deviation of the deformation rate is greater than the non-uniformity threshold, the grid cell is determined to meet the second risk condition; otherwise, it is not met.
[0041] After determining the above two conditions, a comprehensive assessment of the risk status of each grid cell is performed. If a grid cell simultaneously meets both the first and second risk conditions, it is determined to belong to the initial high-risk area. If a grid cell fails to meet both conditions simultaneously (i.e., meets only one condition or neither condition), it is determined not to belong to the high-risk area defined in this step. For grid cells determined to be in the initial high-risk area, a dedicated risk marker field is added to their attributes.
[0042] The spatial extent and attributes of all grid cells marked as high-risk areas in the initial screening are summarized as the output of this step. The purpose of this process is to automatically and efficiently screen out potential high-risk areas from massive grid data that not only have high overall settlement rates but also extremely uneven internal deformation, greatly narrowing the scope of subsequent refined analysis and improving overall assessment efficiency.
[0043] Step 4: For each high-risk area in the initial screening, obtain multi-source disaster-causing factor data for its corresponding range, and conduct collaborative analysis with the temporal deformation data of the area to infer the causes of deformation based on physical mechanisms and predict its development trend.
[0044] For a specific high-risk area to be analyzed, based on the spatial boundaries of the high-risk area, monitoring time series of three preset categories—physical mechanism disaster-causing factors, soil and rock time-dependent disaster-causing factors, and groundwater-related disaster-causing factors—are extracted from a pre-established disaster-causing factor database. Each category of monitoring time series consists of at least two specific physical quantity sequences. Physical mechanism disaster-causing factor sequence: crustal horizontal movement rate, the data source is the continuous operating reference station network covering the region, the northward and eastward component velocity is extracted and the composite result is used to generate the monthly average movement rate time series; seismic energy release, from the National Earthquake Science Data Center, the earthquake catalog within a certain buffer range of the region is obtained, and the time series is generated according to the monthly cumulative energy release (joules). Time-dependent disaster-causing factor sequence in soil and rock: Soil volumetric moisture content, which is obtained by extracting monthly average soil moisture data at a certain depth below the surface of the area and generating a time series; Construction load index, which is obtained by extracting the total building area of all new and expanded projects in the area, summarizing it by quarter, and converting it into monthly data through linear interpolation, as a proxy time series for changes in engineering load; Groundwater-related disaster-causing factor sequence: groundwater extraction volume, obtaining the monthly total groundwater level of the corresponding hydrogeological unit in the area. Extraction volume (ten thousand cubic meters) time series; average groundwater level depth: from the groundwater monitoring well network in this area, the monitoring wells located within the boundary are selected, and their monthly average water level depth (meters) is calculated to generate a time series representing the overall water level change in this area.
[0045] For the vector boundary of the current high-risk area, the temporal deformation data of all pixels in the area are extracted from the comprehensive data table, the average value of the deformation of all pixels at each time point is calculated, and the data are connected in time sequence to generate a regional deformation time series.
[0046] To ensure that all disaster-causing factor monitoring sequences can be directly and effectively correlated with regional deformation time series, strict time alignment and resampling processing were performed on all sequences. The specific steps are as follows: The time points of the regional deformation time series are used as the reference. This series has clearly defined start and end points, time intervals (e.g., months), and total length. (a set of time points). Define the set of reference time points as... .
[0047] For each original disaster-causing factor monitoring sequence, the following judgments and operations are performed sequentially: For each time series of monitoring of disaster-causing factors (such as a monthly groundwater extraction series), the following judgments and operations are performed sequentially: Case 1: Time points are perfectly aligned. If the original set of time points in the sequence is exactly the same as... If they are exactly the same, no processing is needed; use them directly. Scenario 2: The time resolution is higher than the baseline (more dense). If the original data is daily data and the baseline is monthly data, then the original data is aggregated on a monthly basis: for each month in the baseline... Take the arithmetic mean of all daily observations within the month as the corresponding disaster factor value for that month; Scenario 3: The time resolution is lower than the baseline (sparser). If the original data is quarterly data and the baseline is monthly data, then interpolation is performed on the original data: using the quarterly value as a known node, linear interpolation is used on the time axis to calculate the baseline month for each month. Corresponding estimated values of disaster-causing factors; Scenario 4: Data is missing at a certain baseline time point after the above processing. If no valid value is found, linear interpolation is used to calculate and fill the gap using the values of the nearest valid time points before and after the given point.
[0048] After the above processing, each catastrophic factor sequence is converted into a sequence of length [length missing]. The sequence, its first... Each value corresponds to a reference time point .
[0049] Ultimately, a standardized time-series dataset containing regional deformation sequences and all disaster-causing factor sequences is obtained, where all sequences are strictly one-to-one correspondent in the time dimension.
[0050] Based on the standardized time-series dataset described above, for each type of disaster-causing factor, every specific monitoring time series under it is traversed. The Pearson correlation coefficient between each specific monitoring time series and the regional deformation time series is calculated. This calculation measures the degree of consistency in the linear trends of the two series during the same period. The closer the absolute value of the correlation coefficient is to 1, the stronger the linear correlation between the two over time. The correlation formula is as follows: in, The time series is a regional deformation data, and its length is... (i.e., shared) (at each time point), the set of observations for this sequence is denoted as . , ,in For the sequence at the 1st Each benchmark time point The observed values (i.e., the average deformation of the region at that time point) and the sample mean of all observations in the series are: The number to be calculated The monitoring time series of specific disaster-causing factors are as follows: Its length is also The set of observations for this sequence is denoted as , ,in For the sequence at the 1st Each benchmark time point The observed value represents the physical quantity value of the disaster-causing factor at that time point, and the sample mean of all observed values in this sequence is... . For regional deformation sequence With disaster-causing factor sequence The Pearson correlation coefficient between them.
[0051] By calculating the monitoring sequence of each specific disaster-causing factor according to the above formula, a set of Pearson correlation coefficients can be obtained. each The value ranges from -1 to 1. The absolute value represents the strength of the linear correlation between regional deformation and the specific hazard-causing factor over time. The closer the absolute value is to 1, the stronger the linear correlation; the closer the absolute value is to 0, the weaker the linear correlation. The sign indicates the direction of the relationship: a positive value indicates that the two sequences change in the same direction (when the value of the catastrophic factor increases, the deformation tends to increase); a negative value indicates that they change in opposite directions (when the value of the catastrophic factor increases, the deformation tends to decrease).
[0052] After obtaining the set of Pearson correlation coefficients for monitoring sequences of all specific disaster-causing factors. Subsequently, the following processing was performed to quantify the combined impact of various disaster-causing factors: Determine representative correlation coefficients for each category: For the three major disaster factor categories—physical mechanism, soil and rock aging, and groundwater—process them separately. For all correlation coefficients belonging to the same category... Calculate its absolute value The maximum value is selected from the absolute values of all correlation coefficients for that category. The original correlation coefficient corresponding to this maximum value, retaining its sign, is then used as the representative correlation coefficient for that category. The representative correlation coefficients for the three categories mentioned above are denoted as follows: (Physical Mechanisms) (Soil and rock aging type) (Groundwater related category).
[0053] Compare the absolute values of the representative correlation coefficients for the three categories. The category with the largest absolute value is identified as the potential primary disaster-causing category for the high-risk area. Subsequently, based on the identified potential primary disaster-causing categories, corresponding physical prediction models are selected from a pre-established mapping relationship library of disaster-causing factor categories and physical prediction models for further validation. This model library is based on the principles of geomechanics and hydrogeology, and its core correspondences and models are briefly described below: For groundwater-related models—the seepage-stress coupling model: This model is used to simulate the interaction between groundwater pressure changes and soil deformation. This model typically involves the following equations: The seepage equation is Darcy's law: in, The seepage rate vector, Let be the permeability coefficient. For water head; Effective stress principle: in, For the effective stress tensor, For the total stress tensor, Pore water pressure, The tensor is a unit tensor. The deformation and strength of soil are determined by the effective stress. control; Models for time-dependent deformation of soil and rock – soil and rock creep constitutive models: This model describes the continuous deformation (creep) of soil over time under stress, and typically involves the following equations: Creep constitutive equation: in, For strain rate, For the applied stress, and These are material constants; Stress-strain relationship: in, For elastic modulus, The viscosity coefficient; For models targeting disasters caused by physical mechanisms: Based on the specific physical quantities causing the disaster, the following sub-models can be further selected: Dislocation model: used to simulate elastic surface deformation caused by internal crustal dislocations such as fault slip. in, For displacement, For dislocation strength, Shear modulus The distance between the dislocation source and the observation point; Elastic rebound model: The elastic rebound model is used to describe the recovery process of geological structures after being subjected to stress and then releasing the stress. It is applicable to the rebound behavior of the earth's surface after an earthquake. in, For stress, The elastic modulus of the material, In response to the situation.
[0054] The system uses the potential primary disaster-causing category as the query key to automatically search the aforementioned model library. For example, if the disaster is identified as groundwater-related, the solver of the seepage-stress coupling model is loaded; if the disaster is caused by a physical mechanism and the data indicates tectonic activity, the dislocation model is selected first. Subsequently, the generated standardized time-series data is formatted into the input format required by the selected model, ready for verification calculations.
[0055] The regional deformation time series and related disaster-causing factor driving sequences generated in the above steps are input into the selected physical prediction model. Simultaneously, initial parameters are set for the model based on regional geological data. The model's parameter inversion function is run, automatically adjusting the model parameters to ensure the simulated deformation time series matches the observed regional deformation time series as closely as possible. This process aims to find a set of model parameters that best fits the historical observation data. The fitting effect of the model after parameter inversion is checked. If the model can reproduce the historical deformation process well with reasonable parameter values, then the physical mechanism is considered to reasonably explain the observed deformation. Therefore, the potential primary disaster-causing category identified in the above steps is verified as the actual primary disaster-causing cause in the region, and the final deformation causation conclusion is output. At the same time, a set of calibrated model parameters is obtained for subsequent predictions.
[0056] Based on the standardized time series data of Area A shown in Table 1 and Figure 3 Figure 4 Figure 5 It can be seen that groundwater-related disaster-causing factors and regional deformation exhibit a highly consistent and coordinated trend in their temporal evolution. Specifically, when groundwater extraction (F1) increases, the average groundwater level depth (F2) increases accordingly, and the subsidence amplitude of the regional average deformation (D) increases significantly. Conversely, when extraction decreases, the groundwater level depth becomes shallower, and the subsidence trend slows down. The curves of the three factors are highly consistent in terms of key inflection points and magnitude of change, intuitively demonstrating the close causal relationship between groundwater dynamics and surface subsidence.
[0057] In contrast, Figure 3 and Figure 4In the study, the time-series curves of soil and rock time-related disaster-causing factors and physical mechanism disaster-causing factors showed significant differences from the regional deformation curve in terms of variation patterns, phases, and amplitudes, failing to exhibit similar synchronicity and consistency. This comparison further confirms, from a visualization perspective, the quantitative analysis results based on the Pearson correlation coefficient, that is, groundwater activity is the main disaster-causing factor of deformation in this region.
[0058] Based on the rules set in step 4, the correlation coefficient with the largest absolute value is selected as the representative from each category: Groundwater correlation category represents correlation coefficient ; Representational correlation coefficients for soil and rock time-dependent applications ; Physical mechanism category represents correlation coefficient ; Comparing the absolute values of the three, we can see that: 0.98 0.82 The value is 0.15, therefore, the potential primary hazard category for Area A is identified as groundwater-related. This initially suggests that the continued subsidence in this area is likely mainly caused by the drop in water level due to groundwater over-extraction. Based on this preliminary assessment, a seepage-stress coupled physical model will be used for verification and trend prediction.
[0059] Table 1: Standardized Time Series Data of Groundwater-Related Disaster Categories in Area A According to the standardized time-series data of Area B shown in Table 2 and the corresponding time-series comparison chart, the construction load index (F6) and the regional average deformation (D) exhibit significant step-like coordinated changes in their time-series evolution. Specifically, this is manifested as follows: Figure 6 As shown, when the construction load index experiences a step increase at specific time points (such as April and July 2025), the regional settlement deformation immediately exhibits a synchronous acceleration trend. The changes of the two are highly coupled in time, and their magnitudes correspond. This clearly reveals the immediate triggering and continuous driving effect of engineering loading activities on surface settlement.
[0060] like Figure 7 and Figure 8 As shown, although the soil volumetric water content (F5) curve exhibits seasonal fluctuations, its change pattern is smooth and continuous, which does not match the significant step acceleration characteristic in the deformation curve. The curves of groundwater-related disaster-causing factors (extraction volume F1, water level depth F2) also show relatively smooth changes, failing to explain the dramatic response of deformation at the point of sudden load increase. The lack of consistency between these factors and the deformation curve at key change nodes further corroborates the quantitative analysis conclusions from a visualization perspective, namely, that soil and rock time-dependent deformation (driven primarily by engineering loads) is the main disaster-causing category of deformation in this region.
[0061] The relevant correlation coefficients for each category are selected according to the rules: Groundwater correlation category represents correlation coefficient ; Representational correlation coefficients for soil and rock time-dependent applications ; Physical mechanism category represents correlation coefficient ; Comparing the absolute values of the three, 0.9 0.85 Therefore, the potential primary cause of disaster in Area B is identified as soil and rock time-dependent deformation, with a value of 0.15. This indicates that the settlement development in this area is highly coupled with the increased load from construction projects within the area, and soil compression and creep caused by engineering activities are the main causes of deformation. Based on this preliminary assessment, a soil and rock creep constitutive model will be selected for verification and trend prediction.
[0062] Table 2: Standardized Time Series Data of Time-Induced Disaster Categories in Area B (Soil and Geotechnical Engineering) According to the standardized time-series data of region C shown in Table 3 and the corresponding time-series comparison chart, the crustal horizontal movement rate (F3) and the regional average deformation (D) exhibit a highly consistent linear growth trend in time. Figure 9 As shown, the two curves are almost parallel and maintain a stable synchronous rise over time, which intuitively reflects the decisive control of regional continuous tectonic stress loading on surface deformation.
[0063] like Figure 10 and Figure 11 As shown, the monthly seismic energy release (F4) remains at a low background fluctuation level, and its curve differs significantly from the deformation curve pattern, which exhibits a marked linear trend, showing no clear correlation. The curves of soil and rock aging and groundwater-related disaster-causing factors either show gentle fluctuations or other variation patterns, none of which can explain the stable and continuous linear uplift characteristics exhibited by the deformation. This visual contrast strongly supports the quantitative analysis results, indicating that physical mechanisms of disaster (specifically manifested as tectonic activity) are the primary cause of disasters caused by deformation in this region.
[0064] The relevant correlation coefficients for each category are selected according to the rules: Groundwater correlation category represents correlation coefficient ; Representational correlation coefficients for soil and rock time-dependent applications ; Physical mechanism category represents correlation coefficient ; Comparing the absolute values of the three, 1.00 0.90 The value is 0.55, therefore, the potential primary cause of disaster in Area C is identified as a physical mechanism-induced disaster. This indicates that the deformation in this area is mainly controlled by crustal tectonic movements, manifesting as a continuous uplift synchronized with the accumulation of tectonic stress. Based on this preliminary judgment, subsequent verification and trend prediction will prioritize either a dislocation model or an elastic rebound model.
[0065] Table 3: Standardized Time Series Data of Disaster Categories Based on Physical Mechanisms in Area C After validating the physical model of the deformation causes, the deformation development trend of the initially screened high-risk areas is predicted and qualitatively assessed. This prediction does not rely on the aforementioned physical model, but is based on the obtained deformation time series of the region itself, and is achieved through the following steps: First rate calculation: The deformation amount and time of the deformation time series of the region are linearly fitted, and the slope of the fitted line is taken as the overall average deformation rate of the region during the entire monitoring period. Second rate calculation: Select the most recent preset time period in the deformation time series of the region. In this embodiment, the preset time period is one-third of the total monitoring time. Perform linear fitting and use the obtained slope as the recent average deformation rate of the region. The trend determination threshold is preset based on the overall average deformation rate. In this embodiment, the following is set: First threshold = absolute value of overall average deformation rate × 1.5; The second threshold = the absolute value of the overall average deformation rate × 0.8; Among them, the first threshold is greater than the second threshold; The absolute value of the calculated recent average deformation rate is compared with the two thresholds mentioned above, and the development trend of the region is qualitatively determined based on the comparison results: If the absolute value of the recent average deformation rate is greater than or equal to the first threshold, the deformation development trend of the region is determined to be fast. If the absolute value of the recent average deformation rate is less than the first threshold but greater than or equal to the second threshold, the deformation development trend of the region is determined to be medium. If the absolute value of the recent average deformation rate is less than the second threshold, the deformation development trend of the region is determined to be slow. Thus, for each initially screened high-risk area, the above steps have yielded the final causal conclusions regarding its formation and the qualitative results of its development trend (fast, medium, slow). These two results will serve as input parameters for subsequent comprehensive risk level assessment.
[0066] Step 5: Based on the deformation causal conclusions and development trend qualitative results obtained in Step 4, and combined with the average deformation rate of the initially screened high-risk areas, assess the risk level of each initially screened high-risk area, and uniformly assign the remaining grid cells that were not marked as initially screened high-risk areas the lowest risk level representing stability, thereby obtaining the risk level of all grid cells.
[0067] For each region in the initial high-risk region set, the following quantitative rating process is performed: Determine the causal weights: Query the final causal conclusions obtained in step 4 for the high-risk area in the initial screening. In this embodiment, based on the category to which the conclusion belongs, the corresponding fixed causal weight value is obtained from the table below: Table 4: Causal Weight Values The aforementioned weight values are fixed parameters used for subsequent risk index calculations. Their assignment is based on the fact that deformations caused by different factors exhibit varying degrees of controllability, suddenness, and potential hazard. In this embodiment, deformations caused by physical mechanisms (such as tectonic activity) are assigned the highest weight of 0.5 because they are typically difficult to intervene in and have far-reaching effects; deformations caused by soil and rock age-related factors (such as engineering creep) are assigned the next highest weight of 0.3; and deformations related to groundwater (such as pumping settlement) have relatively clear causes and can be mitigated through management measures, thus receiving a basic weight of 0.2. This weighting system is an important component of the risk quantification model of this invention, aiming to transform the differences in attributes caused by different factors into calculable numerical factors.
[0068] Determine the trend coefficient value: Query the qualitative trend result (fast, medium, slow) obtained in step 4 for the high-risk area in the initial screening. Based on this qualitative result, find the corresponding trend coefficient value in the table below: Table 5: Trend Coefficient Values The trend coefficient is a dynamic factor used to quantify the impact of future deformation trends on the current risk status. In the risk assessment model, the current deformation rate represents the urgency of the current situation, while the causal weights represent the inherent severity of the problem. The trend coefficient introduces a dynamic adjustment over a time dimension: when the predicted development trend is rapid, it indicates that the deformation may intensify in the short term, and the risk dynamically increases, thus a larger amplification coefficient (0.9) is assigned; when the development trend is medium or slow, it indicates that the deformation tends to stabilize or ease, thus a smaller coefficient (0.6 or 0.3) is assigned.
[0069] The average deformation rate (unit: mm / year) of the high-risk area in the initial screening is obtained, and its absolute value is taken. Based on the final deformation causal conclusion of the area, the corresponding unique causal weight value is obtained from Table 4. Next, based on the qualitative result of the development trend of the area, the corresponding unique trend coefficient value is obtained from Table 5. Finally, the absolute value of the average deformation rate, the causal weight value, and the trend coefficient value are multiplied together to obtain the comprehensive risk index of the area.
[0070] The calculated comprehensive risk index is compared with a preset risk level threshold range to determine the risk level of the area. In this embodiment, the threshold range is determined by the following method: Calculate the statistical distribution of the comprehensive risk index for all initially screened high-risk areas within the entire assessment area. Sort the values in ascending order, and take the comprehensive risk index values ranked in the top 20% as the high-risk threshold; take the comprehensive risk index values ranked in the top 50% as the medium-risk threshold. This method ensures that the risk level classification is based on the relative distribution of risks within the area.
[0071] If the comprehensive risk index value of a region is greater than or equal to the high-risk threshold, it is assessed as a high-risk region; if the comprehensive risk index value of a region is less than the high-risk threshold but greater than or equal to the medium-risk threshold, it is assessed as a medium-risk region; if the comprehensive risk index value of a region is less than the medium-risk threshold, it is assessed as a low-risk region.
[0072] After the above steps, all initially screened high-risk areas were assigned three risk levels: high, medium, and low. For all other grid cells not marked as initially screened high-risk areas, their risk levels were uniformly assigned the lowest level, indicating stability. The advantage of this operation is that every spatial grid cell within the area to be evaluated (whether it's an initially screened high-risk area or another area) now has a clearly defined risk level attribute.
[0073] According to Table 6 and Figure 5 As can be seen, this embodiment integrates deformation rate, physical causes, and development trends to quantify and generate a comprehensive risk index. For example, both regions A-2 and B-2 have high average deformation rates, but A-2 has a lower cause weight (groundwater related, 0.2) but a qualitative result of "fast" development trend, while B-2 has a higher cause weight (geotechnical aging, 0.3) but a qualitative result of "medium" development trend. The resulting risk indices are close (4.50 vs 3.24), but are classified as "high risk" and "medium risk" respectively based on a threshold division based on relative distribution within the region. This demonstrates the comprehensiveness and discriminative power of the risk quantification model in this method. Although region C-1 is caused by a physical mechanism (weight 0.5), its deformation rate is low and its development trend is slow, resulting in the lowest risk index, and it is classified as "low risk".
[0074] Finally, combining the initial screening results obtained from the above steps, all remaining grid units that were not marked as high-risk areas in the initial screening were uniformly assigned the "stable" level. This data was then integrated with the results in Table 6 to generate risk level data covering the entire area to be assessed, which was used for the subsequent production of thematic risk maps.
[0075] Table 6: Calculation and Level Assessment Results of Risk Index for Each High-Risk Area in the Initial Screening Based on the obtained risk level data of all grid units in the region, and in a dedicated mapping template according to a preset visualization scheme, in this embodiment, red represents high risk, orange represents medium risk, yellow represents low risk, and green represents stable areas. Different risk levels are rendered onto the corresponding spatial locations of the grid units. Simultaneously, map elements such as legends, scale bars, and north arrows are added to generate an intuitive thematic map of surface deformation risk covering the entire area to be assessed.
[0076] The above formulas are all dimensionless calculations. The formulas are derived from software simulations based on a large amount of collected data to obtain the most recent real-world results. The preset parameters in the formulas are set by those skilled in the art according to the actual situation.
[0077] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented in software, the above embodiments can be implemented, in whole or in part, as a computer program product. Those skilled in the art will recognize that the units and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented by electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution.
[0078] The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of this embodiment, depending on actual needs.
[0079] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.
Claims
1. A method for assessing surface deformation risk based on InSAR and multi-source data, characterized in that, The specific steps include: Step 1: Acquire synthetic aperture radar images covering the area to be evaluated at equal time intervals to form a monitoring image set. Perform time-series InSAR processing on the monitoring image set to obtain the time-series deformation data of each pixel location and establish the mapping relationship between pixels and geographic coordinates. Step 2: Divide the area to be evaluated into several grid cells. For each grid cell, calculate the average deformation rate and spatial standard deviation of the deformation rate based on the temporal deformation data of all pixels within it. Step 3: Based on the preset settlement threshold and non-uniformity threshold, mark the grids with average deformation rate exceeding the settlement threshold and spatial standard deviation of deformation rate exceeding the non-uniformity threshold as high-risk areas in the initial screening. Step 4: For each high-risk area in the initial screening, obtain multi-source disaster-causing factor data for its corresponding range, and conduct collaborative analysis with the temporal deformation data of the area to infer the deformation cause based on physical mechanisms and predict its development trend. Step 5: Based on the causes and trends of deformation, and combined with the average deformation rate of the high-risk areas in the initial screening, assess the risk level of each high-risk area in the initial screening, and uniformly assign the lowest risk level representing stability to the remaining grid units that were not marked as high-risk areas in the initial screening, thereby obtaining the risk level of all grid units.
2. The surface deformation risk assessment method based on InSAR and multi-source data according to claim 1, characterized in that: The synthetic aperture radar (SAR) images within the monitoring image set are distributed at equal intervals in time, and the individual SAR images within the monitoring image set correspond to each other spatially. The steps for obtaining the temporal deformation data for each pixel location include: The earliest synthetic aperture radar image is set as the initial image. Each synthetic aperture radar image in the monitoring image set is paired with the initial image to form several image pairs. The acquisition time of the synthetic aperture radar image is marked as the time tag of the image pair. Differential interferometry is performed on each image pair to obtain a differential interferometric phase map and phase unwrapping is performed to calculate the relative surface deformation map between each image pair. For the relative surface deformation map of each image pair, the phase difference of each pixel at each time tag is obtained. The phase difference of each pixel is sorted according to the time tag order to form the temporal deformation data of each pixel.
3. The surface deformation risk assessment method based on InSAR and multi-source data according to claim 2, characterized in that: Establishing the mapping relationship between pixels and geographic coordinates involves the following steps: Obtain the geographic coordinates of the top-left pixel in the synthetic aperture radar image. The geographic coordinates of the bottom right pixel And the total number of rows of pixels inside the synthetic aperture radar image. Total number of columns ; For synthetic aperture radar images located at the first line, number The latitude of a column of pixels is calculated using the following linear interpolation formula. With longitude : in, and These are the row index and column index of the cell, respectively. The value range is 0 to , The value range is 0 to ; The row and column indices of each pixel in the synthetic aperture radar image can be obtained using the above calculation formula. Its geographic coordinates The mapping relationship between them.
4. The surface deformation risk assessment method based on InSAR and multi-source data according to claim 3, characterized in that: The area to be evaluated is divided into regular grid cells. The specific steps include: Determine the scale of the grid cells, and use the smallest inscribed rectangle of the region to be evaluated as the starting range for grid division, setting the lower left corner of the rectangle as the starting origin for grid division; Starting from the origin, the grid is translated at equal intervals along the horizontal and vertical directions according to the grid scale to generate regularly arranged grid boundary lines in sequence. Several regular grid cells are formed by adjacent horizontal and vertical grid boundary lines, from which only grid cells that have effective overlap with the actual spatial range of the area to be evaluated are retained; Traverse all the pixels with obtained temporal deformation data and determine whether they fall within the spatial boundary of a certain effective grid cell based on their geographic coordinates. If a pixel falls into a certain grid cell, then the temporal deformation data corresponding to that pixel is associated with that grid cell.
5. The surface deformation risk assessment method based on InSAR and multi-source data according to claim 4, characterized in that: For each grid cell, based on the temporal deformation data of all pixels within it, calculate the average deformation rate and spatial standard deviation of the deformation rate for that grid cell, including: For each grid cell, extract the temporal deformation data of all its associated cells; Using the observation time corresponding to the time-series deformation data as the independent variable and the value of the time-series deformation data as the dependent variable, a linear fit is performed. The slope of the fitted line is determined as the deformation rate value of the pixel, and its unit is length per time. Calculate the arithmetic mean of the deformation rates of all pixels within the grid cell, and use it as the average deformation rate of the grid cell. Calculate the sample standard deviation of the deformation rate values of all pixels within the grid cell, and use it as the spatial standard deviation of the deformation rate of the grid cell.
6. The surface deformation risk assessment method based on InSAR and multi-source data according to claim 5, characterized in that: Based on the calculated average deformation rate and spatial standard deviation of deformation rate for each cell grid, and compared with the preset preliminary risk screening threshold, grids with an average deformation rate exceeding the settlement threshold and a spatial standard deviation of deformation rate exceeding the non-uniformity threshold are marked as high-risk areas in the initial screening. Specific steps include: Set the settlement threshold and non-uniformity threshold used for judgment; For each grid cell, perform the following judgment: When the average deformation rate of the grid cell exceeds the settlement threshold, it is determined that it meets the first risk condition; When the spatial standard deviation of the deformation rate of the grid cell exceeds the non-uniformity threshold, it is determined that it meets the second risk condition. For grid cells that simultaneously meet the first and second risk conditions mentioned above, they are marked as high-risk areas in the initial screening.
7. The surface deformation risk assessment method based on InSAR and multi-source data according to claim 6, characterized in that: For each high-risk area identified in the initial screening, data on multi-source disaster-causing factors within its corresponding range are acquired and analyzed in conjunction with the temporal deformation data of that area to infer the causes of deformation based on physical mechanisms, including: For each high-risk area in the initial screening, the monitoring time series of multiple pre-defined disaster-causing factors within the area are obtained. The disaster-causing factors include physical mechanism disaster-causing factors, soil and rock time-effect disaster-causing factors, and groundwater-related disaster-causing factors. The monitoring time series of each category of disaster-causing factor consists of at least two monitoring sequences of specific physical quantities that reflect the disaster-causing mechanism of that category. Extract the temporal deformation data of all pixels in the region, calculate the average value of the deformation of all pixels at the same time point, and connect these average values in time sequence to generate a regional deformation time series that represents the overall deformation of the region. For each type of disaster-causing factor, calculate the Pearson correlation coefficient between each specific monitoring time series and the regional deformation time series to obtain the Pearson correlation coefficient subset for each category. For each type of disaster-causing factor, the one with the largest absolute value is selected from the subset of its corresponding Pearson correlation coefficients and used as the representative correlation coefficient between that type of disaster-causing factor and regional deformation. By comparing the absolute values of the correlation coefficients of all categories, the category with the largest absolute value is identified as the main disaster-causing factor category for the initially screened high-risk area. Based on this main disaster-causing factor category, a preliminary hypothesis on the formation cause of the area is proposed.
8. The surface deformation risk assessment method based on InSAR and multi-source data according to claim 7, characterized in that: For the high-risk areas identified in the initial screening, their development trends are predicted based on the results of the collaborative analysis and deformation data, including: Based on the regional deformation time series, the overall average deformation rate over the entire monitoring period is calculated as the first rate; Based on the data of the regional deformation time series within the most recent preset time period, its recent average deformation rate is calculated as the second rate; Compare the absolute value of the second speed with the absolute value of the first speed: If the absolute value of the second rate is greater than or equal to the preset first threshold, it is determined to be a fast trend; If the absolute value of the second rate is less than the first threshold and greater than or equal to the preset second threshold, it is determined to be in a trend. If the absolute value of the second rate is less than the second threshold, it is determined to be a slow trend; Wherein, the first threshold is greater than the second threshold.
9. The surface deformation risk assessment method based on InSAR and multi-source data according to claim 8, characterized in that: Based on the causes and trends of deformation, and combined with the average deformation rate of the initially screened high-risk areas, the risk level of each initially screened high-risk area is assessed. The specific steps include: Based on the final deformation causation conclusions output, three causal categories—physical mechanism-induced disaster, soil and rock age deformation, and groundwater-related disaster—are assigned causal weight values of 0.5, 0.3, and 0.2, respectively. Based on the deformation development trend prediction results, the development trend is characterized into three types: fast, medium and slow, and assigned trend coefficient values of 0.9, 0.6 and 0.3 respectively; The comprehensive risk index of the region is calculated as the product of the region's average deformation rate, causal weight value, and trend coefficient value. The comprehensive risk index is compared with a preset first comprehensive risk threshold and a second comprehensive risk threshold to determine the risk level, wherein the first comprehensive risk threshold is less than the second comprehensive risk threshold; If the overall risk index of the high-risk area in the initial screening is not greater than the first overall risk threshold, then the risk level of the corresponding grid unit is low risk level. If the comprehensive risk index of a high-risk area in the initial screening is greater than the first comprehensive risk threshold but not greater than the second comprehensive risk threshold, then the risk level of the corresponding grid unit is a certain risk level. If the comprehensive risk index of the high-risk area in the initial screening is greater than the second comprehensive risk threshold, then the risk level of the corresponding grid unit is a certain level, which is a high-risk level. For all other grid cells that were not initially identified as high-risk areas, their risk levels are uniformly assigned the lowest risk level, representing stability.