InSAR deformation data processing method

Abstract

The invention discloses an InSAR deformation data processing method, and relates to an interferometric synthetic aperture radar. The method comprises the following steps: extracting contour line data by using DEM data of an InSAR monitoring area, and dividing the InSAR monitoring area into a plurality of sub-areas; extracting the number of water bodies in the monitoring area according to the SAR data intensity graph or the optical image to obtain water body data; performing preprocessing such as abnormal value removal and neighborhood analysis on the InSAR deformation data to obtain preprocessed deformation data; according to the water body data, deformation values distributed in a water area in the InSAR data are removed, and new InSAR deformation data are obtained; performing data correction processing on the InSAR deformation data and the sub-region distribution data to obtain corrected InSAR deformation data; and rendering the corrected InSAR deformation data to generate a monitoring area deformation graph. The method is simple, the data is easy to obtain, the integral error is reduced, and the InSAR deformation data precision of a small-range local region is improved.

Description

Technical Field

[0001] The invention relates to an interferometric synthetic aperture radar, and in particular to an InSAR deformation data processing method. Background Art

[0002] Interferometric Synthetic Aperture Radar (InSAR) is a radar that uses electromagnetic waves to measure phase changes in the monitored area to calculate tiny surface deformations in the monitored area.

[0003] Common InSAR deformation monitoring technologies include Differential Interferometric Synthetic Aperture Radar (DInSAR), Persistent Scatterer Interferometric Synthetic Aperture Radar (PSInSAR) and Small Baseline Subset Interferometric Synthetic Aperture Radar (SbasInSAR). The monitoring accuracy of these technologies can reach millimeter level and they are widely used in urban settlement monitoring, coal mine monitoring, earthquake monitoring, crustal plate movement research, landslide monitoring and other fields. However, InSAR deformation data is affected by multiple factors such as SAR (Synthetic Aperture Radar, SAR) payload positioning deviation, spatiotemporal incoherence, terrain effects, atmospheric interference and crustal movement, which seriously interfere with SAR observation accuracy. In order to improve the accuracy of InSAR data, correction processing is often required. The existing correction of InSAR deformation mainly uses GPS (Global Positioning System) data, GNSS (Global Navigation Satellite System) data, SAR payload precise coordinates, decorrelation technology, DEM data, phase unwrapping algorithm improvement, etc. 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 surface InSAR deformation data of mining areas and use GPS data that is not involved in interpolation for verification and analysis. Huaneng Lancangjiang 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 master image iterative estimation method CSS (Common Scene Stacking, CSS) to correct the atmospheric propagation delays and orbit errors in InSAR data that are uncorrelated in time. Central South University proposed a classification error correction process for interseismic InSAR time series analysis, which integrates the spatiotemporal baseline construction of the interferogram, atmospheric correction of the interferogram using GACOS (Global Atmospheric Correction for Orbital and Atmospheric Errors in Small Baseline Subset Interferometry) and CSS, long-wave correction of the plate model, automatic screening of interferograms based on the triangular closed-loop criterion and phase standard deviation, and time series optimization methods (Liu Hongzhi. Study on wide-area interseismic deformation of the 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] At present, the existing methods mainly correct InSAR data through complex algorithms or high-precision measured deformation data. However, for long-term monitoring areas such as bridges, buildings, and coal mines, their deformation data are greatly affected by crustal movement, which causes the deformation data in these areas to be larger, seriously affecting the judgment of monitoring personnel. The present invention proposes an InSAR deformation data processing method, which uses DEM data (or contour data) to divide the vast InSAR monitoring area into multiple sub-areas according to the contour distribution characteristics of the terrain. Then, in these clearly defined sub-areas, the average deformation of each area is systematically estimated and extracted. This average deformation is regarded as a correction factor and is specifically applied to the correction process of InSAR deformation data in each sub-area, aiming to reduce the error caused by crustal plate movement, accurately extract pure surface deformation information from complex terrain interference, and effectively enhance the accuracy of deformation monitoring and terrain adaptability.

[0006] InSAR is widely used in surface deformation monitoring due to its advantages of all-day, all-weather, high precision and large range. However, for urban areas and coal mining areas, when using millimeter-level precision InSAR (synthetic aperture radar interferometry) technology for deformation monitoring, the InSAR deformation data in these areas are easily affected by the movement of crustal plates, which seriously affects the monitoring results. Summary of the invention

[0007] The purpose of the present invention is to propose an InSAR deformation data processing method, which performs InSAR data preprocessing and correction processing by using water body data, discrete degree analysis and DEM data, cleans outliers, reduces deformation data errors caused by the overall movement of the crust, and improves the accuracy of InSAR deformation data; the present invention corrects interference factors such as terrain effects and crustal movement through regional division and calculation of average deformation, and provides more reliable data support for surface deformation monitoring.

[0008] The present invention comprises the following steps:

[0009] 1) Divide sub-areas according to 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-areas according to the contour data;

[0010] 2) Extraction of water bodies 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 the deformation data in the water area and obtain water body data in binary image format or shape format;

[0011] 3) InSAR deformation data preprocessing: According to the distribution information of InSAR deformation data, outliers are removed; according to the position of each deformation data, each deformation data is calculated for neighborhood analysis, and the number of deformation values ​​within a certain distance around it is calculated. If the corresponding deformation value is small, the deformation value of the point is deleted or the deformation value is changed to a null value to obtain the preprocessed deformation data;

[0012] 4) Removing water body deformation data: According to the water body data obtained in step 2), the deformation data preprocessed in step 3) is removed from the deformation values ​​distributed in the water area in the InSAR data to obtain a new InSAR deformation data;

[0013] 5) InSAR deformation data correction: the InSAR deformation data obtained in step 4) and the sub-region distribution data obtained in step 1) are subjected to data correction processing, the InSAR deformation data are corrected, and corrected InSAR deformation data are obtained;

[0014] 6) InSAR deformation data display: The corrected InSAR deformation data obtained in step 5) are rendered, and the processed deformation data are displayed in a 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, QGis, etc. to generate the contour data Contour of the monitoring area, or by programming, such as using the contour function in Matlab; the division into multiple sub-areas is to divide the areas with higher terrain and lower terrain in the monitoring area to obtain several sub-areas.

[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, which is used to calculate the number of deformation values ​​within the discrete distance R around each deformation value, where R is a positive number; according to the type of InSAR deformation data, if the InSAR deformation data is uniformly distributed, set a square box 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; establish a matrix N to store the InSAR deformation data analysis results, set a threshold N0 to determine whether each element in the matrix N is greater than N0, N0 is a positive integer; if greater than N0, the screening condition is met, if not, the deformation data of the point is removed;

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

[0020] 3.4) Remove points with loose distribution;

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

[0022] 3.2.1) Establish a neighborhood analysis matrix; set a matrix N, each element of which is a positive integer, to record the number of each deformation value; if the InSAR deformation data is uniformly distributed, the matrix N is a two-dimensional matrix, and the matrix size of N is the same as the matrix size of the InSAR deformation data D0, n ij It represents the number of deformation values ​​around the deformation value in the i-th row and j-th column of the matrix N. The initial value of the matrix N is 0. If it is discretely distributed InSAR deformation data, the 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 the InSAR deformation data D0, or simply the number of rows in the D0 matrix). The initial value of the matrix N is 0.

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

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

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

[0026] 3.3.2) Calculate the distance; calculate the distance between the new deformation point and the surrounding deformation points; according to the InSAR deformation data type, if the InSAR deformation data is evenly distributed, determine D0 ij Whether the element (the element in the i-th row and j-th column of the D0 matrix) is deformed data (that is, there is data). If so, extract a neighborhood matrix MR, where iR:i+R means intercepting the elements from the iRth row to the i+Rth row of the matrix D0, and jR:j+R means intercepting all elements from the jRth column to the j+Rth column of the matrix D0; the row and column indexes must be greater than 0 and the index value cannot be greater than the number of rows and columns, otherwise, it will be discarded;

[0027] (1)

[0028] Among them, MR ijRepresents the neighborhood element corresponding to the element in the i-th row and j-th column of the matrix D0 (it is a matrix);

[0029] If D0 ij If the element is not a deformable data (that is, a nan value), then MR ij is a null value matrix;

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

[0031] (2)

[0032] Among them, R i is a column vector, (*,*) indicates the selection of rows and columns, and a single ":" indicates all rows, or it can be understood as an m×1 matrix, where m indicates the number of rows in the D0 matrix; n i Represents the number of deformation values ​​around the i-th deformation value;

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

[0034] (3)

[0035] The sum() function represents the sum function, the not() function represents the matrix inversion function, and the isnan() function represents the function that determines whether an element is null. If it is, it is 1, otherwise it is 0. ij Represents the element in the i-th row and j-th column of the matrix N, which is a non-negative integer;

[0036] If the InSAR deformation data is discretely distributed, the R obtained in step 3.3.2) is statistically i The number of elements in is greater than or equal to R, and the formula is as follows:

[0037] (4)

[0038] Among them, the sum() function represents the sum, R i ≥R indicates the judgment vector R i Whether the element in is greater than or equal to R, if it is satisfied, it is 1, otherwise it is 0;

[0039] 3.3.4) Loop judgment: judge whether the processing is completed according to the InSAR data type:

[0040] If the InSAR deformation data is uniformly distributed, determine whether 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 whether i is greater than the number of rows in the D0 matrix. If so, 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 judgment: Only when the InSAR data type is uniformly distributed deformation data, perform column loop judgment; judge whether i is greater than the number of rows of the D0 matrix. If so, set j=1 and go 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, that is, the comparison of n in the neighborhood analysis ij The size of the population number N0 is set in step 3.2). If the InSAR deformation data is uniformly distributed, find the elements in the matrix N that are less than N0 and assign them to nan (null value). The formula is as follows:

[0044] (5)

[0045] The find() function represents a search function, which finds the row and column positions where the elements in the matrix N are less than N0, assigns the corresponding elements to nan, and assigns the obtained new D0 data to D1;

[0046] If the InSAR deformation data is discretely distributed, find the row where the element in the matrix N is less than N0 and delete it; the formula is as follows:

[0047] (6)

[0048] Among them, the find() function represents a search function, which searches for the row where the element in the matrix N is less than N0. The “:” indicates selecting all columns, that is, obtaining the row index where the element in the matrix N is less than N0. The row index extracts all elements in the corresponding row in the D0 matrix and assigns them to D1.

[0049] In step 4), the specific steps of removing the water body deformation data may be: according to the water body data Shp obtained in step 2) water , remove the deformation values ​​distributed in the water area in the InSAR data, and obtain a new InSAR deformation data D2; according to the InSAR data type, the details are as follows:

[0050] If the uniformly distributed InSAR deformation data is converted into binary image format data Shp waterPerform dot multiplication (multiply corresponding elements) with the InSAR preprocessed deformation data D1 obtained in step 3.4), where D2 ij Represents the deformation value data of the i-th row and j-th column in 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 converted into the deformation data Shp format. water If the point is in the water, the deformation data of the point (i.e., the row of data where the point is located) is deleted to obtain a new InSAR deformation data D2, where D2 i represents the i-th deformation value data in the new InSAR deformation data D2, that is, the i-th row of deformation data;

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

[0054] 5.1) Calculate the sub-region correction factor; calculate the mean value of the InSAR deformation data in each sub-region sub, and obtain the sub-region sub i The average deformation value Ave i ;

[0055] 5.2) Correct the InSAR deformation data; subtract the sub-area sub obtained in step 051 from the InSAR deformation data D2 after removing the water deformation data. i The average deformation value Ave i , obtain the final InSAR deformation data D3 that needs to be corrected.

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

[0057] There are two types of deformation data processed by InSAR technology. 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 the present invention is to use DEM data (or contour data) to divide the vast InSAR monitoring area into multiple sub-areas according to the contour distribution characteristics of the terrain. Then, in these clearly defined sub-areas, the average deformation of each area is systematically estimated and extracted. This average deformation is regarded as a correction factor and is specifically applied to the correction process of InSAR deformation data of each sub-area, aiming to reduce the error caused by the movement of crustal plates, accurately extract surface deformation information from complex terrain interference, and improve the accuracy of deformation monitoring and terrain adaptability.

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

[0060] Existing technical solutions can be divided into two categories. The first category is to improve the accuracy of InSAR deformation data by improving the InSAR processing algorithm, and to reduce the InSAR data decoherence error, baseline error, registration error, atmospheric delay effect error, terrain error or phase unwrapping error by improving the algorithm; the second category is to use multi-source data to correct InSAR deformation data, and the common ones are GPS data, GNSS data, and angle inversion data. The present invention uses DEM data for correction, which has the following advantages compared with the existing technology:

[0061] 1. The correction algorithm is simple, and there is no need to improve the complex InSAR algorithm to improve the accuracy of InSAR deformation data;

[0062] 2. Data is easy to obtain. Compared with GPS data, GNSS data, and angular inversion data, DEM data is easy to obtain;

[0063] 3. Reduce the overall error caused by crustal movement and improve the accuracy of InSAR deformation data in small local areas, such as coal mines, bridges, etc. BRIEF DESCRIPTION OF THE DRAWINGS

[0064] Figure 1 is a flow chart.

[0065] Figure 2 shows the technical principle route.

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

[0067] Figure 4 shows the neighborhood analysis of InSAR deformation data. DETAILED DESCRIPTION

[0068] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the following embodiments will further illustrate the present invention in conjunction with the accompanying drawings. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0069] like Figure 1 and 2, the specific process of the present invention is as follows:

[0070] Step 01: Divide sub-areas according to terrain; obtain the contour data of the monitoring area according to the DEM data of the monitoring area, and then divide it into several sub-areas; 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; ArcGIS, QGis and other tools can be used to generate the contour line data of the monitoring area, or it can be obtained by programming, such as using the contour function in Matlab;

[0072] Step 012: Divide the sub-areas; Divide the sub-areas according to the contour data Contour obtained in step Step 011; Divide the areas with higher terrain and lower terrain in the monitoring area according to the contour data Contour of the monitoring area to obtain several sub-areas sub i ,like Figure 3 As shown, 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: Extract water bodies 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 (outliers are processed in step Step 04); extract water body data in the monitoring area based on the SAR data intensity map or optical image; obtain water body data Shp in binary image format or shape format water ;

[0074] Step 03: InSAR deformation data preprocessing: remove some discrete data (outliers) according to the distribution information of InSAR deformation data D0; , calculate each deformation data to perform neighborhood analysis, calculate the number of deformation values ​​n within a certain distance around, and if the number n is smaller, the corresponding deformation value will be deleted or the deformation value will be changed to a null value nan, such as Figure 4 As shown, the preprocessed deformation data D1 is obtained, which is specifically as follows: The specific steps are as follows:

[0075] Step 031: Set discrete distance; Set a discrete distance R (R is a positive number) to calculate the number of deformation values ​​within the surrounding distance R of each deformation value; according to the InSAR deformation data type D0, it can be divided into two cases:

[0076] The first type: uniformly distributed InSAR deformation data, set a square box (side length is 2R+1, such as: R=10), such as Figure 4 As shown in Figure (a);

[0077] The second type: For discretely distributed InSAR deformation data, an actual distance R (such as R = 100m) is set. Figure 4 As shown in Figure (b);

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

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

[0080] The second type: discretely distributed InSAR deformation data, the 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 the InSAR deformation data D0, in simple terms, the number of rows in the D0 matrix), and n i Represents the number of deformation values ​​around the i-th deformation value, and the initial value of the 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 the matrix N is greater than or equal to N0. If it is greater than or equal to N0, the screening condition is met. If not, the deformation data of the point is removed.

[0082] Step033: Calculate the number of deformation values ​​around each deformation value; the specific steps are as follows:

[0083] Step0331: Establish the initial value of the loop; Establish the initial value of the loop, and loop to calculate the number of deformation values ​​around each deformation value, as follows:

[0084] The first type: For uniformly distributed InSAR deformation data, two cycle initial values ​​i and j are established, and their initial values ​​are 1, representing row and column cycles respectively, and their initial values ​​are both 1;

[0085] The second type: For discretely distributed InSAR deformation data, a number cycle 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 the surrounding deformation points; according to the type of InSAR deformation data, it can be divided into two cases:

[0087] The first type: Uniformly distributed InSAR deformation data, determine D0 ij Whether the element (the element in the i-th row and j-th column of the D0 matrix) is deformed data (that is, there is data). If so, extract a neighborhood matrix MR, where iR:i+R means intercepting the elements from the iRth row to the i+Rth row of the matrix D0, and jR:j+R means intercepting all elements from the jRth column to the j+Rth column of the matrix D0; the row and column indexes must be greater than 0 and the index value cannot be greater than the number of rows and columns, otherwise, it will be discarded;

[0088] (1)

[0089] Among them, MR ij Represents the neighborhood element corresponding to the element in the i-th row and j-th column of the matrix D0 (it is a matrix);

[0090] If D0 ij If the element is not a deformable data (that is, a nan value), then MR ij is a null value matrix;

[0091] The second type: discretely distributed InSAR deformation data, each deformation point is calculated as n i The distance between deformation points is as follows:

[0092] (2)

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

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

[0095] The first type: For uniformly distributed InSAR deformation data, the number of deformation points around the deformation point in the i-th row and j-th column in DW is counted, that is, MR is calculated. ij The number of non-empty elements in the matrix is ​​as follows:

[0096] (3)

[0097] The sum() function represents the sum function, the not() function represents the matrix inversion function, and the isnan() function represents the function that determines whether an element is a null value. If it is, it is 1, otherwise it is 0. ij Represents the element in the i-th row and j-th column of the matrix N, which is a non-negative integer;

[0098] The second type: discretely distributed InSAR deformation data, the R obtained in Step 0332 is statistically i The number of elements in is greater than or equal to R, and the formula is as follows:

[0099] (4)

[0100] The sum() function represents the sum, R i ≥R indicates the judgment vector R i Whether the element in is greater than or equal to R, if it is satisfied, it is 1, otherwise it is 0;

[0101] Step 0334: Loop judgment: judge whether the processing is completed; according to the 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 so, set i=1 and go to Step 0334; otherwise, set i=i+1 and return to Step 033.

[0103] The second type: discretely distributed InSAR deformation data, determine whether i is greater than the number of rows of the D0 matrix. If so, set i=1 and go to Step 034; otherwise, set i=i+1 and return to Step 033.

[0104] Step0334: 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 deformation data); determine whether i is greater than the number of rows of the D0 matrix. If so, set j=1 and enter step 034; otherwise, set j=j+1 and return to step 033;

[0105] Step 034: Remove discrete points; remove points with loose distribution, that is, compare the n points in the neighborhood analysis. ij The size of the group number N0 is set in step Step032; the details are as follows:

[0106] The first method: uniformly distributed InSAR deformation data, find the elements in the matrix N that are less than N0 and assign them to nan (null value) as follows:

[0107] (5)

[0108] The find() function represents a search function, which finds the row and column positions where the elements in the matrix N are less than N0, and then assigns the corresponding elements to nan, and assigns the obtained new D0 data to D1;

[0109] The second method: For discretely distributed InSAR deformation data, find the row where the element in matrix N is less than N0 and delete it; the formula is as follows:

[0110] (6)

[0111] The find() function is a search function that finds the row where the element in matrix N is less than N0. The “:” indicates selecting all columns, that is, obtaining the row index where the element in matrix N is less than N0. The row index extracts all elements in the corresponding row in matrix D0 and assigns them to D1.

[0112] Step04: Remove the InSAR deformation data in the water area; remove the InSAR deformation data in the water area; according to the water body data Shp obtained in Step02 water , remove the deformation values ​​distributed in the water area in the InSAR data, and obtain a new InSAR deformation data D2; according to the InSAR data type, it can be divided into two categories, as follows:

[0113] The first category: uniformly distributed InSAR deformation data, the binary image format data Shp water Perform dot multiplication with the InSAR preprocessed deformation data D1 obtained in step Step034 (multiply the corresponding elements), where D2 ij Represents the deformation value data of the i-th row and j-th column in the new InSAR deformation data D2; the formula is as follows:

[0114] (7)

[0115] The second category: discretely distributed InSAR deformation data, which converts the coordinates of each element in the InSAR deformation data into the deformation data Shp format water If the point is in the water, the deformation data of the point (i.e., the row of data where the point is located) is deleted to obtain a new InSAR deformation data D2, where D2 i Represents the i-th deformation value data (i.e., the i-th row of deformation data) in the new InSAR deformation data D2;

[0116] Step 05: InSAR deformation data correction; the InSAR deformation data D2 obtained in step Step 04 and the sub-region distribution data sub i Perform data correction processing to correct the InSAR deformation data. The specific steps are as follows:

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

[0118] Step 052: Correct the InSAR deformation data; subtract the sub-area sub obtained in step Step 051 from the InSAR deformation data D2 after removing the water area deformation data. i The average deformation value Ave i , obtain the final InSAR deformation data D3 that needs to be corrected;

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

[0120] The principle of the present invention is to use DEM data (or contour data) to divide the vast InSAR monitoring area into multiple sub-areas according to the contour distribution characteristics of the terrain, such as Figure 3 As shown. Then, in these clearly defined sub-areas, the average deformation of each area is systematically estimated and extracted. This average deformation is regarded as a correction factor and is specifically applied to the correction process of InSAR deformation data in each sub-area, aiming to reduce the error caused by crustal plate movement, accurately extract surface deformation information from complex terrain interference, and improve the accuracy of deformation monitoring and terrain adaptability.

[0121] There are two types of deformation data processed by InSAR technology. The first type is uniformly distributed InSAR deformation data, such as the InSAR deformation data obtained after DInSAR processing, as shown in the schematic diagram. Figure 3 As shown in Figure (b), the gray grid indicates that InSAR data is stored, and the white grid indicates a null value nan; the second type is discretely distributed InSAR deformation data, such as the InSAR deformation data obtained after PSInSAR processing, as shown in the schematic diagram Figure 3 As shown in Figure (c), each gray point represents an InSAR deformation data, and each point contains real coordinates (such as longitude and latitude) and deformation data.

[0122] The above embodiments are only preferred embodiments of the present invention and cannot be considered to limit the scope of the present invention. All equivalent changes and improvements made within the scope of the present invention should still fall within the scope of the present invention.

Claims

1. An InSAR deformation data processing method, characterized in that The following steps are involved: 1) Divide sub-areas according to 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-areas according to the contour data; 2) Extraction of water bodies 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 the deformation data in the water area and obtain water body data in binary image format or shape format; 3) InSAR deformation data preprocessing: remove outliers based on the distribution information of InSAR deformation data; According to the position of each deformation data, calculate each deformation data for neighborhood analysis, and calculate the number of deformation values ​​within a certain distance around it. If the corresponding deformation value is smaller, delete the deformation value of the point or change the deformation value to a null value to obtain the preprocessed deformation data; 4) Removing water body deformation data: According to the water body data obtained in step 2), the deformation data preprocessed in step 3) is removed from the deformation values ​​distributed in the water area in the InSAR data to obtain a new InSAR deformation data; 5) InSAR deformation data correction: the InSAR deformation data obtained in step 4) and the sub-region distribution data obtained in step 1) are subjected to data correction processing, the InSAR deformation data are corrected, and corrected InSAR deformation data are obtained; 6) InSAR deformation data display: The corrected InSAR deformation data obtained in step 5) are rendered, and the processed deformation data are displayed in a graphical form to generate a deformation map of the monitoring area.

2. The InSAR deformation data processing method according to claim 1, characterized in that In step 1), the contour data is obtained by using ArcGIS, QGis tools to generate the contour data Contour of the monitoring area, or by using the contour function in Matlab; the division into multiple sub-areas is to divide the areas with higher terrain and lower terrain in the monitoring area to obtain several sub-areas.

3. The InSAR deformation data processing method according to claim 1, characterized in that 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: 3.1) Set a discrete distance R, which is used to calculate the number of deformation values ​​within the discrete distance R around each deformation value, where R is a positive number; according to the type of InSAR deformation data, if the InSAR deformation data is uniformly distributed, set a square box with a side length of 2R+1; if the InSAR deformation data is discretely distributed, set an actual distance R; 3.2) Set the number of groups; Establish a matrix N to store the InSAR deformation data analysis results, set a threshold N0, and judge whether each element in the matrix N is greater than N0, where N0 is a positive integer; if it is greater than N0, the screening condition is met, and if it is not met, the deformation data of the point is removed; 3.3) Calculate the number of deformation values ​​around each deformation value; 3.4) Remove points with loose distribution.

4. The InSAR deformation data processing method according to claim 3, characterized in that In step 3.2), the specific steps of setting the number of groups are as follows; 3.2.1) Establish a neighborhood analysis matrix; set a matrix N, where each element is a positive integer, used to record the number of each deformation value; if the InSAR deformation data is uniformly distributed, the matrix N is a two-dimensional matrix, and the size of the matrix N is the same as the size of the InSAR deformation data D0 matrix, n ij It represents the number of deformation values ​​around the deformation value in the i-th row and j-th column in the matrix N. The initial value of the matrix N is 0; If the InSAR deformation data is discretely distributed, the matrix N is an m×1 matrix, that is, a column vector; where m represents the number of deformation points, that is, the number of rows of the D0 matrix; the initial value of the matrix N is 0; 3.2.2) Set the domain analysis threshold; Set the threshold N0 to determine whether each element in the matrix N is greater than or equal to N0, where N0 is a positive integer. If it is greater than or equal to N0, the screening condition is met; if not, the deformation data of the point is removed.

5. The InSAR deformation data processing method according to claim 3, characterized in that In step 3.3), the specific steps of calculating the number of deformation values ​​around each deformation value are as follows: 3.3.1) Establish loop initial values; Establish loop initial values, and loop to calculate the number of deformation values ​​around each deformation value. If the InSAR deformation data is uniformly distributed, establish two loop initial values ​​i and j, whose initial values ​​are 1, representing row and column loops respectively, and their initial values ​​are both 1; if the InSAR deformation data is discretely distributed, establish 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 InSAR deformation data type, if the InSAR deformation data is evenly distributed, determine D0 ij Whether the element is deformation data, if so, extract a neighborhood matrix MR, where iR:i+R means intercepting the elements from the iRth row to the i+Rth row of the matrix D0, and jR:j+R means intercepting all elements from the jRth column to the j+Rth column of the matrix D0; the row and column indexes must be greater than 0 and the index value cannot be greater than the number of rows and columns, otherwise, they are discarded; （1） Among them, MR ij Represents the neighborhood element corresponding to the element in the i-th row and j-th column of matrix D0; D0 ij Represents the element in row i and column j in matrix D0; If D0 ij If the element is not a deformation data, that is, a nan null value, then MR is directly set ij is a null value matrix; If the InSAR deformation data is discretely distributed, then each deformation point is calculated as i The distance between deformation points is as follows: （2） Among them, R i is a column vector. (*,*) indicates the selection of rows and columns. A single ":" indicates all rows, or it can be understood as an m×1 matrix, where m indicates the number of rows in the D0 matrix; n i Represents the number of deformation values ​​around the i-th deformation value; 3.3.3) Count the number of adjacent deformation points; perform neighborhood analysis on each InSAR deformation value and count the number of deformation values ​​within the discrete distance. If the InSAR deformation data is uniformly distributed, count the number of deformation points around the deformation point in the i-th row and j-th column of DW, that is, calculate MR ij The number of non-empty elements in the matrix is ​​as follows: （3） The sum() function represents the sum function, the not() function represents the matrix inversion function, and the isnan() function represents the function that determines whether an element is null. If it is, it is 1, otherwise it is 0. ij Represents the element in the i-th row and j-th column of the matrix N, which is a non-negative integer; If the InSAR deformation data is discretely distributed, the R obtained in step 3.3.2) is statistically i The number of elements in is greater than or equal to R, and the formula is as follows: （4） Among them, the sum() function represents the sum, R i ≥R indicates the judgment vector R i Whether the element in is greater than or equal to R, if it is satisfied, it is 1, otherwise it is 0; 3.3.4) Loop judgment: judge whether the processing is completed according to the InSAR data type: If the InSAR deformation data is uniformly distributed, determine whether 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.1); If the InSAR deformation data is discretely distributed, determine whether i is greater than the number of rows in the D0 matrix. If so, set i=1 and proceed to step 3.4); otherwise, set i=i+1 and return to step 3.3.1); 3.3.5) Column loop judgment: Only when the InSAR data type is uniformly distributed deformation data, perform column loop judgment; judge whether i is greater than the number of rows of the D0 matrix. If it is satisfied, set j=1 and proceed to the next step; otherwise, set j=j+1 and return to step 3.3.1).

6. The InSAR deformation data processing method according to claim 3, characterized in that In step 3.4), the removal of loosely distributed points, that is, the comparison of n in the neighborhood analysis ij The size of the population number N0 is set in step 3.2). If the InSAR deformation data is uniformly distributed, find the elements in the matrix N that are less than N0 and assign them to nan. The formula is as follows: （5） The find() function represents a search function, which finds the row and column positions where the elements in the matrix N are less than N0, assigns the corresponding elements to nan, and assigns the obtained new D0 data to D1; If the InSAR deformation data is discretely distributed, find the row where the element in the matrix N is less than N0 and delete it; the formula is as follows: （6） Among them, the find() function represents the search function, which finds the row where the element in matrix N is less than N0. ":" means selecting all columns, that is, getting the index of the row where the element in matrix N is less than N0. The row index extracts all elements in the corresponding row in matrix D0 and assigns them to D1.

7. The InSAR deformation data processing method according to claim 1, characterized in that In step 4), the specific steps of removing the water body deformation data are as follows: according to the water body data Shp obtained in step 2) water , remove the deformation values ​​distributed in the water area in the InSAR data, and obtain a new InSAR deformation data D2; according to the InSAR data type, the details are as follows: If the uniformly distributed InSAR deformation data is converted into binary image format data Shp water Perform dot multiplication (multiply corresponding elements) with the InSAR preprocessed deformation data D1 obtained in step 3.4), where D2 ij Represents the deformation value data of the i-th row and j-th column in the new InSAR deformation data D2; the formula is as follows: （7） If the InSAR deformation data is discretely distributed, the coordinates of each element in the InSAR deformation data are converted into the deformation data Shp format. water If the point is in the water, the deformation data of the point, that is, the row of data where the point is located, is deleted to obtain a new InSAR deformation data D2, where D2 i It represents the i-th deformation value data in the new InSAR deformation data D2, that is, the i-th row of deformation data.

8. The InSAR deformation data processing method according to claim 1, characterized in that In step 5), the specific steps of the InSAR deformation data correction are: 5.1) Calculate the sub-region correction factor; calculate the mean value of the InSAR deformation data in each sub-region sub, and obtain the sub-region sub i The average deformation value Ave i ; 5.2) Correct the InSAR deformation data; subtract the sub-area sub obtained in step 051 from the InSAR deformation data D2 after removing the water deformation data. i The average deformation value Ave i , obtain the final InSAR deformation data D3 that needs to be corrected.

9. The InSAR deformation data processing method according to claim 1, characterized in that In step 6), the InSAR deformation data is displayed, and the InSAR deformation data D3 corrected in step 5) is read in ArcGIS and QGis geographic information processing software for rendering to produce a deformation map of the monitoring area.

Images