An InSAR deformation data processing method

Abstract

The application discloses an InSAR deformation data processing method, and relates to interferometric synthetic aperture radar. Isoline data is extracted from DEM data of an InSAR monitoring area, and the InSAR monitoring area is divided into a plurality of sub-areas; water body data is obtained by extracting the number of water bodies in the monitoring area according to SAR data intensity maps or optical images; pre-processing is performed on InSAR deformation data, such as abnormal value removal and neighborhood analysis, to obtain pre-processed deformation data; according to the water body data, deformation values distributed in water areas in the InSAR data are removed to obtain new InSAR deformation data; data correction processing is performed on the InSAR deformation data and sub-area distribution data to obtain corrected InSAR deformation data; and the corrected InSAR deformation data is rendered to generate a deformation map of the monitoring area. The method is simple, data is easy to obtain, overall error is reduced, and the precision of InSAR deformation data in a small range and local area is improved.

Description

Technical Field

[0001] This invention relates to interferometric synthetic aperture radar, and more particularly to an InSAR deformation data processing method. Background Technology

[0002] Interferometric Synthetic Aperture Radar (InSAR) is a radar system that uses electromagnetic waves to measure changes in the phase of a monitored area, thereby calculating minute surface deformations in that area.

[0003] Common InSAR deformation monitoring techniques include Differential Interferometric Synthetic Aperture Radar (DInSAR), Persistent Scatterer Interferometric Synthetic Aperture Radar (PSInSAR), and Small Baseline Subset Interferometric Synthetic Aperture Radar (SbasInSAR). These techniques offer millimeter-level accuracy and are widely used in urban subsidence monitoring, coal mine monitoring, earthquake monitoring, crustal plate movement research, and landslide monitoring. However, InSAR deformation data is affected by various factors such as SAR payload positioning bias, spatiotemporal incoherence, topographic effects, atmospheric interference, and crustal movement, which severely interfere with SAR observation accuracy. Correction processing is often necessary to improve the accuracy of InSAR data. Existing methods for correcting InSAR deformation mainly involve using GPS (Global Positioning System) data, GNSS (Global Navigation Satellite System) data, precise coordinates of the SAR payload, decoherence techniques, DEM data, and improved phase unwrapping algorithms to correct InSAR deformation data.

[0004] Inner Mongolia University of Technology proposed a genetic-kriging interpolation model that integrates GPS-InSAR data (Cai Qiuhuan. Analysis and Research on Surface Deformation in Mining Areas Based on GPS and InSAR Data [D]. Inner Mongolia University of Technology, 2023, 20-32) to correct the InSAR deformation data of the mining area surface, and used GPS data that was not involved in the interpolation for verification analysis. Huaneng Lancang River Hydropower Co., Ltd. disclosed a high-precision three-dimensional deformation inversion method suitable for reservoir bank slopes (Li Li. A high-precision three-dimensional deformation inversion method suitable for reservoir bank slopes [P]. Chinese Patent). This method can obtain more comprehensive and accurate surface deformation information by integrating GNSS and InSAR data. Tymofyeyeva et al. (Tymofyeyeva E, Fialko Y. Mitigation of atmospheric phase delays in InSAR data, with application to the eastern California shear zone [J]. Journal of Geophysical Research: Solid Earth, 2015, 120(8): 5952-63) proposed a common scene stacking (CSS) method to correct for temporally uncorrelated atmospheric propagation delays and orbital errors in InSAR data. Central South University proposed a classification error correction process for inter-seismic InSAR time series analysis, which integrates spatiotemporal baseline construction of interferograms, atmospheric correction of interferograms using GACOS (Global Atmospheric Correction for Orbital and Atmospheric Errors in Small Baseline Subset Interferometry) and CSS, long-wave correction of plate models, automatic screening of interferograms based on triangular closed-loop criteria and phase standard deviation, and time series optimization methods (Liu Hongzhi. Study on wide-area inter-seismic deformation of Haiyuan Fault Zone based on time series InSAR error correction [D]. Central South University, 2023. DOI:10.27661 / d.cnki.gzhnu.2023.003159) to obtain high-precision surface deformation data.

[0005] Currently, existing methods primarily correct InSAR data using complex algorithms or high-precision measured deformation data. However, for long-term monitoring areas such as bridges, buildings, and coal mines, deformation data is significantly affected by crustal movement, resulting in inflated deformation data that severely impacts monitoring personnel's judgment. This invention proposes an InSAR deformation data processing method. This method utilizes DEM data (or contour line data) and, based on the contour line distribution characteristics of the terrain, meticulously divides the vast InSAR monitoring area into multiple sub-regions. Then, within these clearly defined sub-regions, the average deformation is systematically estimated and extracted for each region. This average deformation is considered a correction factor and is specifically applied to the InSAR deformation data correction process for each sub-region. The aim is to reduce errors caused by crustal plate movement, accurately extract pure surface deformation information from complex terrain interference, and effectively enhance the accuracy and terrain adaptability of deformation monitoring.

[0006] InSAR is widely used for surface deformation monitoring due to its advantages such as all-day, all-weather, high precision, and wide coverage. However, in urban areas and coal mining areas, when using millimeter-precision InSAR (Synthetic Aperture Radar Interferometry) technology for deformation monitoring, the InSAR deformation data in these areas is easily interfered with by the movement of tectonic plates, which seriously affects the monitoring results. Summary of the Invention

[0007] The purpose of this invention is to propose an InSAR deformation data processing method. By utilizing water body data, dispersion analysis, and DEM data for InSAR data preprocessing and correction, outliers are cleaned up, and deformation data errors caused by overall crustal movement are reduced, thereby improving the accuracy of InSAR deformation data. This invention also corrects for interference factors such as topographic effects and crustal movement by dividing the region and calculating the average deformation, providing more reliable data support for surface deformation monitoring.

[0008] This invention includes the following steps:

[0009] 1) Divide the area into sub-regions based on terrain: Use the elevation information in the DEM data of the InSAR monitoring area to extract contour data, and divide the InSAR monitoring area into multiple sub-regions based on the contour data;

[0010] 2) Water body extraction in the monitoring area: Extract the number of water bodies in the monitoring area based on the SAR data intensity map or optical image to remove deformation data in the water area and obtain water body data in binary map format or shape format;

[0011] 3) InSAR deformation data preprocessing: Based on the distribution information of InSAR deformation data, outliers are removed; based on the location of each deformation data, neighborhood analysis is performed on each deformation data, and the number of deformation values ​​within a certain distance is calculated. If the number of deformation values ​​is small, the deformation value at that point is deleted or the deformation value is changed to a null value, thus obtaining the preprocessed deformation data.

[0012] 4) Remove water body deformation data: Based on the water body data obtained in step 2), remove the deformation values ​​distributed in the water body from the deformation data preprocessed in step 3) to obtain a new InSAR deformation data;

[0013] 5) InSAR deformation data correction: Perform data correction processing on the InSAR deformation data obtained in step 4) and the sub-region distribution data obtained in step 1) to correct the InSAR deformation data and obtain the corrected InSAR deformation data.

[0014] 6) InSAR deformation data display: The corrected InSAR deformation data obtained in step 5) is rendered and the processed deformation data is displayed in graphical form to generate a deformation map of the monitoring area.

[0015] In step 1), the contour data can be obtained by using tools such as ArcGIS and QGis to generate contour data of the monitoring area, or by programming, such as using the contour function in MATLAB; the division into multiple sub-regions involves dividing the monitoring area into areas with higher and lower elevations to obtain several sub-regions.

[0016] In step 3), the InSAR deformation data preprocessing includes setting the discrete distance, setting the number of groups, calculating the number of deformation values ​​around each deformation value, and removing discrete points. The specific steps are as follows:

[0017] 3.1) Set a discrete distance R to calculate the number of deformation values ​​within the discrete distance R of each deformation value. R is a positive number. According to the InSAR deformation data type, if the InSAR deformation data is uniformly distributed, set a square frame with a side length of 2R+1; if the InSAR deformation data is discretely distributed, set an actual distance R.

[0018] 3.2) Set the number of groups; create a matrix N to store the InSAR deformation data analysis results, and set a threshold N0 to determine whether each element in matrix N is greater than N0, where N0 is a positive integer; if it is greater than, the filtering condition is met; if it is not met, the deformation data of that point is removed.

[0019] 3.3) Calculate the number of deformation values ​​surrounding each deformation value;

[0020] 3.4) Remove loosely distributed points;

[0021] In step 3.2), the specific steps for setting the number of groups are as follows;

[0022] 3.2.1) Establish a neighborhood analysis matrix; set up a matrix N, where each element is a positive integer, to record the number of each deformation value; if the InSAR deformation data is uniformly distributed, then matrix N is a two-dimensional matrix, and the size of matrix N is the same as the size of the InSAR deformation data D0 matrix, n ij This represents the number of deformation values ​​surrounding the deformation value in the i-th row and j-th column of matrix N. The initial value of matrix N is 0. If it is discrete InSAR deformation data, then matrix N is an m×1 matrix (column vector), where m represents the number of deformation points (equivalent to the number of deformation points in InSAR deformation data D0, simply put, it is the number of rows in matrix D0). The initial value of matrix N is 0.

[0023] 3.2.2) Set the domain analysis threshold; set a threshold N0 to determine whether each element in matrix N is greater than or equal to N0. N0 is a positive integer. If it is greater than or equal to N0, the filtering condition is met. If it is not met, the deformation data of that point is removed.

[0024] In step 3.3), the specific steps for calculating the number of deformation values ​​around each deformation value are as follows:

[0025] 3.3.1) Establish initial values ​​for the loop; establish initial values ​​for the loop, and calculate the number of deformation values ​​around each deformation value in a loop. If the InSAR deformation data is uniformly distributed, establish two initial values ​​i and j, with an initial value of 1, representing the row and column loops respectively, and their initial values ​​are both 1; if the InSAR deformation data is discretely distributed, establish a loop value i with an initial value of 1.

[0026] 3.3.2) Calculate the distance; calculate the distance between the new deformation point and the surrounding deformation points; based on the InSAR deformation data type, if the InSAR deformation data is uniformly distributed, determine D0. ij If an element in the i-th row and j-th column of matrix D0 is deformed data (i.e., contains data), then extract a neighborhood matrix MR, where iR:i+R represents extracting elements from the iR-th to the i+R-th row of matrix D0, and jR:j+R represents extracting all elements from the jR-th to the j+R-th column of matrix D0. The row and column indices must be greater than 0 and the index value cannot be greater than the number of rows and columns; otherwise, they are discarded.

[0027] (1)

[0028] Among them, MR ijThis represents the neighborhood element (which is a matrix) corresponding to the element in the i-th row and j-th column of matrix D0.

[0029] If D0 ij If the element is not deformable data (i.e., a null value), then directly set MR. ij It is a null matrix;

[0030] If the InSAR deformation data is discretely distributed, then each deformation point is calculated as n i The distance between deformation points is calculated using the following formula:

[0031] (2)

[0032] Among them, R i It is a column vector, where (*,*) represents selecting rows and columns, and a single ":" represents all rows. Alternatively, it can be understood as an m×1 matrix, where m represents the number of rows in matrix D0; n i This represents the number of deformation values ​​surrounding the i-th deformation value;

[0033] 3.3.3) Count the number of adjacent deformation points; for each InSAR deformation value, perform neighborhood analysis and count the number of deformation values ​​within discrete distances. If the InSAR deformation data is uniformly distributed, count the number of deformation points surrounding the deformation point in the i-th row and j-th column of the DW, i.e., calculate the MR. ij The formula for the number of non-empty elements in a matrix is ​​as follows:

[0034] (3)

[0035] The functions `sum()` calculate the summation, `not()` negates the matrix, and `isnan()` checks if any element is null (1 if null, 0 otherwise). ij Let represent the element in the i-th row and j-th column of matrix N, which is a non-negative integer;

[0036] If the InSAR deformation data are discretely distributed, then the R obtained in step 3.3.2) is statistically analyzed. i The number of elements greater than or equal to R in a given set is determined by the following formula:

[0037] (4)

[0038] The sum() function represents summation, and R i ≥R represents the judgment vector R i If an element in the set is greater than or equal to R, the value is 1; otherwise, it is 0.

[0039] 3.3.4) Row loop check; determine whether processing is complete, based on InSAR data type:

[0040] If the InSAR deformation data are uniformly distributed, determine if i is greater than the number of rows in the D0 matrix. If so, set i=1 and proceed to step 3.3.5; otherwise, set i=i+1 and return to step 3.3).

[0041] If the InSAR deformation data is discretely distributed, determine if i is greater than the number of rows in the D0 matrix. If it is, set i=1 and proceed to step 3.4; otherwise, set i=i+1 and return to step 3.3.

[0042] 3.3.5) Column loop check: Column loop check is performed only if the InSAR data type is uniformly distributed deformed data; check if i is greater than the number of rows in the D0 matrix. If it is, set j=1 and proceed to step 3.4); otherwise, set j=j+1 and return to step 3.3).

[0043] In step 3.4), the removal of loosely distributed points, i.e., in the comparative neighborhood analysis, n... ij Similar to setting the population size N0 in step 3.2), if the InSAR deformation data is uniformly distributed, find the elements in matrix N that are less than N0 and assign them the value nan (null value). The formula is as follows:

[0044] (5)

[0045] The find() function is a search function that finds the row and column positions in matrix N where the element is less than N0, assigns the corresponding element to nan, and assigns the new data D0 to D1.

[0046] For discrete InSAR deformation data, find and delete the rows in matrix N where the elements are less than N0; the formula is as follows:

[0047] (6)

[0048] The `find()` function is a search function that finds the row in matrix N where the element is less than N0. The colon ":" indicates that all columns are selected. In other words, it gets the index of the row in matrix N where the element is less than N0. The row index is used to extract all elements in the corresponding row of matrix D0 and assign them the value D1.

[0049] In step 4), the specific steps for removing water body deformation data can be as follows: based on the water body data Shp obtained in step 2). water After removing the deformation values ​​distributed in the water area from the InSAR data, a new InSAR deformation data D2 is obtained; the specific details are as follows, depending on the InSAR data type:

[0050] If the InSAR deformation data is uniformly distributed, the data in binary map format (Shp) will be used. waterPerform a dot product (element-wise multiplication) with the InSAR preprocessed deformation data D1 obtained in step 3.4), where D2 ij This represents the deformation value in the i-th row and j-th column of the new InSAR deformation data D2; the formula is as follows:

[0051] (7)

[0052] If the InSAR deformation data is discretely distributed, the coordinates of each element in the InSAR deformation data are substituted into the shapefile (SHP) format deformation data. water If the point is located in water, then the deformation data for that point (i.e., the row of data containing that point) is deleted, and a new InSAR deformation data D2 is obtained. i This represents the i-th deformation value in the new InSAR deformation data D2, i.e., the i-th row of deformation data;

[0053] In step 5), the specific steps for InSAR deformation data correction are as follows:

[0054] 5.1) Calculate the sub-region correction factor; calculate the mean of the InSAR deformation data within each sub-region (sub) to obtain the sub-region (sub) correction factor. i The average deformation value Ave within i ;

[0055] 5.2) Correcting InSAR deformation data; subtracting the sub-region data obtained in step 051 from the InSAR deformation data D2 after removing water area deformation data. i The average deformation value Ave within i To obtain the final InSAR deformation data D3 that requires correction.

[0056] In step 6), the InSAR deformation data display can render the corrected InSAR deformation data D3 obtained in step 5). The corrected InSAR deformation data in step 5) can be read in geographic information processing software such as ArcGIS and QGIS, and the deformation data can be rendered to produce a deformation map of the monitoring area.

[0057] InSAR technology processes two types of deformation data. The first type is uniformly distributed InSAR deformation data, such as InSAR deformation data obtained after DInSAR processing. The second type is discretely distributed InSAR deformation data, such as InSAR deformation data obtained after PSInSAR processing.

[0058] The principle of this invention is to utilize DEM data (or contour data) and, based on the contour line distribution characteristics of the terrain, meticulously divide the vast InSAR monitoring area into multiple sub-regions. Then, within these clearly defined sub-regions, the average deformation of each region is systematically estimated and extracted. This average deformation is considered a correction factor and is specifically applied to the InSAR deformation data correction process for each sub-region, aiming to reduce errors caused by tectonic plate movement, accurately extract surface deformation information from complex terrain interference, and improve the accuracy and terrain adaptability of deformation monitoring.

[0059] Compared with the prior art, the advantages of the present invention are as follows:

[0060] Existing technical solutions can be divided into two categories. The first category improves the accuracy of InSAR deformation data by modifying InSAR processing algorithms, thereby reducing errors such as incoherence, baseline, registration, atmospheric delay, topographic, or phase unwrapping. The second category utilizes multi-source data to correct InSAR deformation data, commonly including GPS, GNSS, and inverted image data. This invention uses DEM data for correction, offering the following advantages compared to existing technologies:

[0061] 1. The correction algorithm is simple and does not require modification of complex InSAR algorithms to improve the accuracy of InSAR deformation data;

[0062] 2. Data is easy to obtain; compared to GPS data, GNSS data, and angular reflection data, DEM data is easier to obtain.

[0063] 3. Reduce overall errors caused by crustal movement and improve the accuracy of InSAR deformation data in small local areas, such as coal mines and bridges. Attached Figure Description

[0064] Figure 1 is a flowchart.

[0065] Figure 2 shows the technical principle route.

[0066] Figure 3 shows the sub-regional division of the monitoring area.

[0067] Figure 4 shows the neighborhood analysis of InSAR deformation data. Detailed Implementation

[0068] To make the objectives, technical solutions, and advantages of this invention clearer, the following embodiments will be used in conjunction with the accompanying drawings to further illustrate the invention. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.

[0069] like Figure 1 and 2The specific process of this invention is as follows:

[0070] Step 01: Divide the area into sub-regions based on the terrain; obtain contour data of the monitoring area based on the DEM data of the monitoring area, and then divide it into several sub-regions; the specific process is as follows:

[0071] Step 011: Extract contour lines; Use the elevation information in the DEM data of the monitoring area to extract the contour line data of the monitoring area; You can use tools such as ArcGIS and QGis to generate the contour data of the monitoring area, or you can obtain it through programming, such as using the contour function in MATLAB;

[0072] Step 012: Divide the area into sub-regions; based on the contour data obtained in Step 011, divide the area into sub-regions; according to the contour data of the monitoring area, divide the area into higher and lower elevation regions, obtaining several sub-regions. i ,like Figure 3 As shown in the figure, Figure (a) represents the sub-region division of the monitoring area, Figure (b) represents the sub-region division of uniformly distributed deformation data, and Figure (c) represents the sub-region division of discretely distributed deformation data.

[0073] Step 02: Water body extraction in the monitoring area; extract the number of water bodies in the monitoring area based on the SAR data intensity map or optical image; this step aims to remove deformation data in the water area (outlier handling is done in Step 04); extract water body data in the monitoring area based on the SAR data intensity map or optical image; obtain water body data in binary map format or shape format (Shp). water ;

[0074] Step 03: InSAR deformation data preprocessing; based on the distribution information of InSAR deformation data D0, remove some discrete data (outliers); based on the location of each deformation data point... For each deformation data point, a neighborhood analysis is performed. The number of deformation values ​​*n* within a certain distance is calculated. If the number of deformation values ​​*n* is small, the corresponding deformation value is deleted or replaced with a null value *nan*. Figure 4 As shown, the preprocessed deformation data D1 is obtained, and the specific steps are as follows:

[0075] Step 031: Set the discrete distance; Set a discrete distance R (R is a positive number) to calculate the number of deformation values ​​within a distance R around each deformation value; According to the InSAR deformation data type D0, there are two cases:

[0076] The first type: uniformly distributed InSAR deformation data, then define a rectangular frame (side length 2R+1, e.g., R=10), such as... Figure 4 As shown in Figure (a);

[0077] The second type: Discretely distributed InSAR deformation data, then a set actual distance R (e.g., R=100m) is used. Figure 4 As shown in Figure (b);

[0078] Step 032: Set the number of groups; create an empty matrix N to store the InSAR deformation data analysis results, and then set a threshold N0 (a positive integer) to determine whether each element in matrix N is greater than N0. If it is greater, the filtering condition is met; if it is not met, the deformation data of that point is removed; the details are as follows;

[0079] Step 0321: Establish a neighborhood analysis matrix; set up a matrix N, where each element is a positive integer, to record the number of each deformation value; specifically as follows: First type: uniformly distributed InSAR deformation data, then matrix N is a two-dimensional matrix, and the size of matrix N is the same as the size of the InSAR deformation data D0 matrix, n ij This represents the number of deformation values ​​surrounding the deformation value in the i-th row and j-th column of matrix N. The initial value of matrix N is 0.

[0080] The second type: Discretely distributed InSAR deformation data, then matrix N is an m×1 matrix (column vector), where m represents the number of deformation points (equivalent to the number of deformation points in InSAR deformation data D0, simply put, the number of rows in matrix D0), and n i This represents the number of deformation values ​​surrounding the i-th deformation value, and the initial value of matrix N is 0;

[0081] Step 0322: Set the domain analysis threshold; set the threshold N0 (a positive integer) to determine whether each element in matrix N is greater than or equal to N0. If it is greater than or equal to N0, the filtering condition is met; otherwise, the deformed data of that point is removed.

[0082] Step 033: Calculate the number of deformation values ​​surrounding each deformation value; the specific steps are as follows:

[0083] Step 0331: Establish initial values ​​for the loop; establish initial values ​​for the loop, and repeatedly calculate the number of deformation values ​​around each deformation value, as follows:

[0084] The first type: uniformly distributed InSAR deformation data, then establish two loop initial values ​​i and j, with an initial value of 1, to represent row and column loops respectively, and their initial values ​​are both 1;

[0085] The second type: for discrete InSAR deformation data, a loop value i is established with an initial value of 1;

[0086] Step 0332: Calculate the distance; calculate the distance between the new deformation point and surrounding deformation points; based on the InSAR deformation data type, there are two cases:

[0087] The first method: using uniformly distributed InSAR deformation data to determine D0. ij If an element in the i-th row and j-th column of matrix D0 is deformed data (i.e., contains data), then extract a neighborhood matrix MR, where iR:i+R represents extracting elements from the iR-th to the i+R-th row of matrix D0, and jR:j+R represents extracting all elements from the jR-th to the j+R-th column of matrix D0. The row and column indices must be greater than 0 and the index value cannot be greater than the number of rows and columns; otherwise, they are discarded.

[0088] (1)

[0089] Among them MR ij This represents the neighborhood element (which is a matrix) corresponding to the element in the i-th row and j-th column of matrix D0.

[0090] If D0 ij If the element is not deformable data (i.e., a null value), then directly set MR. ij It is a null matrix;

[0091] The second type: for discrete InSAR deformation data, the deformation point is calculated as n. i The distance between deformation points is calculated using the following formula:

[0092] (2)

[0093] Among them, R i It is a column vector, where (*,*) represents the selection of rows and columns, and a single ":" represents all rows, or it can be understood as an m×1 matrix, where m represents the number of rows in matrix D0;

[0094] Step 0333: Count the number of adjacent deformation points; for each InSAR deformation value neighborhood analysis, count the number of deformation values ​​within discrete distances, which can be divided into two cases:

[0095] The first method involves uniformly distributed InSAR deformation data. This is done by counting the number of deformation points surrounding the deformation point in the i-th row and j-th column of the DW, i.e., calculating the MR. ij The formula for the number of non-empty elements in a matrix is ​​as follows:

[0096] (3)

[0097] The `sum()` function performs a summation, the `not()` function negates the matrix, and the `isnan()` function checks if an element is null (returning 1 if null, 0 otherwise). ij Let represent the element in the i-th row and j-th column of matrix N, which is a non-negative integer;

[0098] The second type: Discretely distributed InSAR deformation data, then the R obtained in step 0332 is statistically analyzed. i The number of elements greater than or equal to R in a given set is determined by the following formula:

[0099] (4)

[0100] The sum() function represents summation, R i ≥R represents the judgment vector R i If an element in the set is greater than or equal to R, the value is 1; otherwise, it is 0.

[0101] Step 0334: Row loop check; determine if processing is complete; based on InSAR data type, it can be divided into two categories, as follows:

[0102] The first type: uniformly distributed InSAR deformation data. Then, determine whether i is greater than the number of rows of the D0 matrix. If it is, set i=1 and proceed to step 0334; otherwise, set i=i+1 and return to step 033.

[0103] The second type: discrete InSAR deformation data, then determine whether i is greater than the number of rows of the D0 matrix. If it is, let i=1 and proceed to step 034; otherwise, i=i+1 and return to step 033.

[0104] Step 0334: Column loop judgment; This step is executed only when the InSAR data type is the first type (i.e., the InSAR data is uniformly distributed deformed data); Determine whether i is greater than the number of rows of the D0 matrix. If it is, set j=1 and proceed to step 034; otherwise, set j=j+1 and return to step 033.

[0105] Step 034: Remove discrete points; remove loosely distributed points, i.e., in the comparative neighborhood analysis, n ij In step 032, the population size N0 is set; specifically as follows:

[0106] The first method uses uniformly distributed InSAR deformation data. It identifies elements in matrix N that are less than N0 and assigns them the value nan (null). The formula is as follows:

[0107] (5)

[0108] The find() function is a search function that finds the row and column positions in matrix N where the element is less than N0, then assigns the corresponding element to nan, and assigns the new data D0 to D1.

[0109] The second method involves finding and deleting rows in matrix N containing elements smaller than N0, based on discrete InSAR deformation data. The formula is as follows:

[0110] (6)

[0111] The find() function is a search function that finds the row in matrix N where the element is less than N0. The colon ":" indicates that all columns are selected. In other words, it gets the index of the row in matrix N where the element is less than N0. The row index is used to extract all elements in the corresponding row of matrix D0 and assign them the value D1.

[0112] Step 04: Remove InSAR deformation data from the water body; remove InSAR deformation data from the water body; based on the water body data Shp obtained in Step 02... water After removing the deformation values ​​distributed in the water area from the InSAR data, a new InSAR deformation data D2 is obtained. Based on the InSAR data type, it can be divided into two categories, as follows:

[0113] Category 1: Uniformly distributed InSAR deformation data, converted into binary image format (Shp). water Perform a dot product (element-wise multiplication) with the InSAR preprocessed deformation data D1 obtained in Step 034, where D2 ij This represents the deformation value in the i-th row and j-th column of the new InSAR deformation data D2; the formula is as follows:

[0114] (7)

[0115] The second category: Discretely distributed InSAR deformation data, where the coordinates of each element in the InSAR deformation data are input into the shapefile (SHP) format deformation data. water If the point is located in water, then the deformation data for that point (i.e., the row of data containing that point) is deleted, and a new InSAR deformation data D2 is obtained. i This represents the i-th deformation value in the new InSAR deformation data D2 (i.e., the i-th row of deformation data).

[0116] Step 05: InSAR Deformation Data Correction; Correct the InSAR deformation data D2 obtained in Step 04 and the sub-region distribution data sub obtained in Step 012. i Data correction processing is performed to correct the InSAR deformation data. The specific steps are as follows:

[0117] Step 051: Calculate the sub-region correction factor; calculate the mean of the InSAR deformation data within each sub-region (sub), and obtain the sub-region (sub) correction factor. i The average deformation value Ave within i ;

[0118] Step 052: Correct InSAR deformation data; subtract the sub-region data obtained in Step 051 from the InSAR deformation data D2 after removing water area deformation data. i The average deformation value Ave within i To obtain the final InSAR deformation data D3 that requires correction;

[0119] Step 06: Generate a deformation map of the monitoring area; Render the corrected InSAR deformation data D3 obtained in Step 052 to display and obtain the deformation map of the monitoring area. The specific operation steps are as follows: You can read the InSAR deformation data D3 obtained in Step 052 in geographic information processing software such as ArcGIS and QGIS, and render the deformation data to produce a deformation map of the monitoring area.

[0120] The principle of this invention is to utilize DEM data (or contour data) and, based on the contour line distribution characteristics of the terrain, meticulously divide the vast InSAR monitoring area into multiple sub-regions, such as... Figure 3 As shown. Furthermore, within these clearly defined sub-regions, the average deformation of each region is systematically estimated and extracted. This average deformation is considered a correction factor and is applied specifically to the InSAR deformation data correction process for each sub-region. The aim is to reduce errors caused by tectonic plate movement, accurately extract surface deformation information from complex terrain interference, and improve the accuracy and terrain adaptability of deformation monitoring.

[0121] InSAR technology processes two types of deformation data. The first type is uniformly distributed InSAR deformation data, such as InSAR deformation data obtained after DInSAR processing, as illustrated in the diagram. Figure 3 As shown in Figure (b), gray squares represent InSAR data and white squares represent null values ​​(NaN). The second type is discretely distributed InSAR deformation data, such as InSAR deformation data obtained after PSInSAR processing, as illustrated in the diagram. Figure 3 As shown in Figure (c), each gray area represents an InSAR deformation data point, and each point contains the true coordinates (such as latitude and longitude) and deformation data.

[0122] The above embodiments are merely preferred embodiments of the present invention and should not be considered as limiting the scope of the present invention. All equivalent variations and improvements made within the scope of the present invention should still fall within the patent coverage of the present invention.

Claims

1. An InSAR deformation data processing method, characterized in that The method comprises the following steps: 1) Sub-regional division according to terrain: using InSAR to monitor the elevation information in the DEM data of the region, extracting contour data, and dividing the InSAR monitoring region into multiple sub-regions according to the contour data; 2) Water body extraction in the monitoring region: extracting the water body number in the monitoring region according to the SAR data intensity diagram or optical image to remove the deformation data in the water area, and obtaining the water body data in the binary diagram format or shape format; 3) InSAR deformation data preprocessing: removing outliers according to the distribution information of the InSAR deformation data; According to the position of each deformation data, the neighborhood analysis of each deformation data is calculated, and the number of deformation values within a certain distance is calculated. The smaller the number, the smaller the deformation value, and then the point deformation value is deleted or the deformation value is changed to a null value, to obtain the preprocessed deformation data; 4) Removing deformation data in water area: according to the water body data obtained in step 2), the deformation data in the water area in the preprocessed deformation data in step 3) is removed to obtain a new InSAR deformation data; 5) InSAR deformation data correction: performing data correction processing on the InSAR deformation data obtained in step 4) and the sub-regional distribution data obtained in step 1), correcting the InSAR deformation data, and obtaining corrected InSAR deformation data; The specific steps of the InSAR deformation data correction are as follows: 5.1) Calculate sub-region correction factor; calculate the average of InSAR deformation data within each sub-region sub, to obtain the average deformation value Ave within the sub-region sub i i ;​ 5.2) Correcting the InSAR deformation data; subtracting the average deformation value Ave within the sub-region sub obtained in step 5.1) from the InSAR deformation data D2 processed by removing the water area deformation data i i to obtain the final InSAR deformation data D3 required for correction;​ 6) InSAR deformation data display: rendering the corrected InSAR deformation data obtained in step 5), displaying the processed deformation data in the form of graphics, and generating a deformation map of the monitoring region.

2. The InSAR deformation data processing method of claim 1, wherein In step 1), the contour data is generated by using Arcgis, QGis tools to generate monitoring region contour data Contour, or by using the contour function in matlab; the multiple sub-regions are divided into regions with higher and lower terrain in the monitoring region to obtain a plurality of sub-regions.

3. The InSAR deformation data processing method of claim 1, wherein In step 3), the InSAR deformation data preprocessing comprises setting a discrete distance, setting a group number, calculating the number of deformation values around each deformation value, and removing discrete points, and the specific steps are as follows: 3.1) Set a discrete distance R for calculating the number of deformation values within a distance R around each deformation value, R is a positive number; according to the type of InSAR deformation data, if the InSAR deformation data is uniformly distributed, a square box with a side length of 2R+1 is set; if the InSAR deformation data is discretely distributed, an actual distance R is set; 3.2) Set the group number; A matrix N is established for storing the InSAR deformation data analysis results, and a threshold N0 is set for judging whether each element in the matrix N is greater than N0, N0 is a positive integer; if it is greater than, it meets the screening condition, if it does not meet, the point deformation data is removed; 3.3) Calculate the number of deformation values around each deformation value; 3.4) Remove the points with loose distribution.

4. The InSAR deformation data processing method of claim 3, wherein In step 3.2), the specific steps of setting the group number are as follows: 3.2.1) Establishing the neighborhood analysis matrix; setting a matrix N, each element of which is a positive integer, for recording the number of each deformation value; if it is uniformly distributed InSAR deformation data, the matrix N is a two-dimensional matrix, the matrix size of N is the same as that of the InSAR deformation data D0 matrix, n ij The number of deformation values in the neighborhood of the deformation value in the i-th row and the j-th column of the matrix N is represented, and the initial value of the matrix N is 0. If the InSAR deformation data is discrete distribution, the matrix N is an m x 1 matrix, namely a column vector; wherein, m represents the number of deformation points, namely the number of rows of the D0 matrix; the initial value of the matrix N is 0; 3.2.2) Setting neighborhood analysis threshold; setting threshold N0 for judging whether each element in the matrix N is greater than or equal to N0, N0 is a positive integer, if greater than or equal to N0, it meets the screening condition, if not, the deformation data of the point is removed.

5. The InSAR deformation data processing method of claim 3, wherein In step 3.3), the specific steps of calculating the number of deformation values around each deformation value are as follows: 3.3.1) Establishing loop initial value; establishing loop initial value, loop calculating the number of deformation values around each deformation value, if the InSAR deformation data is uniformly distributed, establishing two loop initial values i and j, the initial values of which are 1, respectively representing row and column loops, the initial values of which are both 1; if the InSAR deformation data is discrete distribution, establishing a number loop value i with an initial value of 1; 3.3.2) Calculate the distance; Calculate the distance between the new deformation point and the surrounding deformation points; According to the type of InSAR deformation data, if the InSAR deformation data is uniformly distributed, judge D0 ij whether the element is deformation data, if yes, extract a neighborhood matrix MR, where i-R:i+R represents the elements of the i-Rth row to the i+Rth row of the intercepted matrix D0, j-R:j+R represents all elements of the j-Rth column to the j+Rth column of the intercepted matrix D0; The row and column indexes must be greater than 0 and the index values cannot be greater than the number of rows and columns, otherwise, they are discarded; MR ij = D0(i-R:i+R,j-R:j+R) (1) where MR ij denotes the neighborhood element corresponding to the element in the ith row and jth column of matrix D0; D0 ij denotes the element in the ith row and jth column of matrix D0; If D0 ij If the element is not deformation data, i.e. nan, then MR ij is the null matrix. If the InSAR deformation data are discrete, the distance between each deformation point and n i deformation points is calculated, as follows: where R i is a column vector, (, ) indicates the selection of rows and columns, and " alone indicates all rows or is understood as an m x 1 matrix, where m represents the number of rows of the D0 matrix; n i represents the number of peripheral deformation values of the i-th deformation value. 3.3.3) Count the number of adjacent deformation points; for each InSAR deformation value neighborhood analysis, count the number of deformation values within the discrete distance, if the InSAR deformation data is uniformly distributed, then count the number of deformation points around the i-th row and j-th column deformation point in DW, that is, calculate MR ij The number of non-empty elements in the matrix is as follows: n ij = sum(not(isnan(MR ij ))) (3) Wherein, sum() function represents the sum function, not() function represents the matrix inversion function, isnan() function represents the judgment element is null, is 1, otherwise 0; n ij The element in the i-th row and j-th column of matrix N is a non-negative integer. If the InSAR deformation data are discrete, then the R i element is greater than or equal to the number of R elements, and the formula is as follows: n i = sum(R i ≥R) (4) wherein the sum() function represents summation, R i ≥R represents a judgment whether the element in the vector R i is greater than or equal to R, and if it is satisfied, it is 1, otherwise it is 0. 3.3.4) Row loop judgment; judging whether the processing is completed, according to the type of InSAR data: If the InSAR deformation data is uniformly distributed, judging whether i is greater than the number of rows of the D0 matrix, if yes, setting i = 1, and entering step 3.3.5) at the same time; otherwise, setting i = i + 1, and returning to step 3.3.1); If the InSAR deformation data is discrete distribution, judging whether i is greater than the number of rows of the D0 matrix, if yes, setting i = 1, and entering step 3.4) at the same time; otherwise, setting i = i + 1, and returning to step 3.3.1); 3.3.5) Column loop judgment; only when the type of InSAR data is uniformly distributed deformation data, the column loop judgment is executed; judging whether i is greater than the number of rows of the D0 matrix, if yes, setting j = 1, and entering the next step at the same time; otherwise, setting j = j + 1, and returning to step 3.3.1).

6. The InSAR deformation data processing method of claim 3, wherein In step 3.4), the points with loose distribution are removed, i.e. in the comparison neighborhood analysis, n ij As in step 3.2), the size of the group number N0is set, if the InSAR deformation data is uniformly distributed, the elements in the matrix N smaller than N0are found, and they are assigned as nan, the formula is as follows: Wherein, the find() function represents a search function, finding the row and column positions of the elements in the matrix N that are less than N0, and assigning the corresponding elements as nan; the obtained new D0 data is assigned to D1; If the InSAR deformation data is discrete distribution, finding the row in which the elements in the matrix N are less than N0 and deleting it; the formula is as follows: D1 = D0(find(N ≥ N0), :) (6) Wherein, the find() function represents a search function, searching for the row in which the elements in the matrix N are less than N0, and ": " represents selecting all columns, namely obtaining the row index in which the elements in the matrix N are less than N0, and extracting all elements in the corresponding row in the D0 matrix, and assigning them to D1.

7. The InSAR deformation data processing method of claim 3, wherein In step 4), the specific steps for removing the water area deformation data are as follows: according to the water body data Shp obtained in step 2) water , removing the deformation values distributed in the water area in the InSAR data to obtain a new InSAR deformation data D2; according to the type of InSAR data, the specific steps are as follows: If the InSAR deformation data is uniformly distributed, the binary graph format data Shp water Point multiplication (corresponding elements are multiplied) is performed on the InSAR preprocessed deformation data D1 obtained in step 3.4), and D2 ij represents the deformation value data of the i-th row and the j-th column in the new InSAR deformation data D2; the formula is as follows: D2 = D1 * Shp water (7) If the discrete distributed InSAR deformation data, each element coordinate of the InSAR deformation data is brought into the shape format deformation data Shp water If the point is in the water area, the deformation data of the point is deleted, that is, the data of the row where the point is located is deleted, to obtain a new InSAR deformation data D2, wherein D2 i represents the i-th deformation value data in the new InSAR deformation data D2, that is, the i-th row deformation data.

8. The InSAR deformation data processing method of claim 1, wherein In step 6), the InSAR deformation data display, reading the corrected InSAR deformation data D3 in step 5) in Arcgis, QGis geographic information processing software for rendering, producing deformation map of the monitoring area.

Images