An InSAR monitoring method for continuous surface deformation under vegetation cover

Abstract

The present invention discloses an InSAR monitoring method for continuous deformation of the surface under vegetation coverage conditions, comprising the following steps: S1: identifying pixel units with deep monitoring potential; S2: determining and cutting high potential and areas of interest within the range; S3: aggregating and dividing subsets within a time series; S4: obtaining the cumulative deformation amount and average deformation rate results within the time range of a single time subset based on a matrix singular value decomposition method; S5: for different scenarios, a discrete subset fusion method in the time dimension is used to achieve the acquisition of the average deformation rate of the vegetation-covered surface in the entire time series and high-time frequency surface deformation results. This method makes full use of limited high-quality phase information, solves the problem of limited coverage of results obtained by traditional InSAR methods, and reduces the amount of InSAR calculations, expanding the application scenarios of InSAR technology.

Description

Technical Field

[0001] The invention relates to the technical field of surface exploration, and in particular to an InSAR monitoring method for continuous surface deformation under vegetation coverage conditions. Background Art

[0002] Synthetic Aperture Radar (ISAR) has been widely used in high-precision surface deformation monitoring tasks with large scales and long time series, which is of great help in improving the understanding of surface movement. With different application scenarios and accuracy requirements, a series of improved InSAR technologies have been widely developed, among which the PS-InSAR (Persistent Scatterer InSAR) method for permanent scatterers and the SBAS-InSAR (Small Baseline Subsets InSAR) method for distributed scatterers are the main development trends. However, although the SBAS-InSAR method for natural surfaces can obtain InSAR deformation results with higher coverage by optimizing the interference combination to a certain extent, it is still impossible to obtain effective monitoring results in some low-coherence areas caused by long-term vegetation coverage, which has become an important factor restricting the development of InSAR technology.

[0003] Although it is still impossible to obtain InSAR monitoring results for these long-term stable vegetation coverage areas in principle, long-term observation experiments in these areas show that there are some areas with periodic high-coherence features due to vegetation cover changes or human activities in long time series. In order to effectively mine and utilize these limited high-quality information in long time series, a small number of improved algorithms have been studied recently, including the InSAR method using natural years as the segmentation standard and the ISBAS-InSAR (Intermittent SBAS-InSAR) method. However, these methods always face limitations such as limited coverage of results, large InSAR calculation volume and a large amount of computing power wasted, and many accuracy limiting factors. Summary of the invention

[0004] The present invention aims to solve at least one of the technical problems existing in the prior art. To this end, one purpose of the present invention is to propose an InSAR monitoring method for continuous surface deformation under vegetation cover conditions, which makes full use of limited high-quality phase information to solve the problem of limited coverage of results obtained by traditional InSAR methods, while reducing the amount of InSAR calculations and expanding the application scenarios of InSAR technology.

[0005] In order to solve the above problems, the present invention provides an InSAR monitoring method for continuous surface deformation under vegetation coverage conditions, comprising the following steps:

[0006] S1: Extract the coherence coefficient of the interference pixel from the long time series differential interference image generated by complex conjugate multiplication, combine the maximum value of the pixel time series coherence coefficient, the duration of each high coherence feature time subset of the pixel in the time series, and the complete duration feature of the pixel high coherence feature, to identify the pixel unit with deep monitoring potential;

[0007] S2: Using the pixel monitoring potential results of S1 as prior knowledge, determine and crop the high potential and interesting regions within the range;

[0008] S3: Based on the functional relationship between InSAR monitoring accuracy and the number of interferograms and the high-coherence subset construction strategy, the subsets within the entire time series are aggregated and divided;

[0009] S4: The matrix-based singular value decomposition method is used to obtain the cumulative deformation and average deformation rate results within a single time subset time range;

[0010] S5: Based on S4, for different scenarios, the discrete subset fusion method in the time dimension is adopted to achieve the average deformation rate of the vegetation-covered surface in the entire time series and the acquisition of high-temporal-frequency surface deformation results.

[0011] Preferably, S1 comprises the following sub-steps:

[0012] S11: Perform preliminary inspection on long-time series SAR image data with the same track number and map number, and remove SAR images with obvious offset; use Enhanced Spectral Diversity registration method to achieve accurate sub-pixel registration of sequence images;

[0013] S12: Set appropriate time dimension baseline and space dimension baseline to generate an interference combination with high quality expectations; call the complex data conjugate multiplication function to obtain the interference image result, the formula is expressed as:

[0014] ; (1)

[0015] Where, U1 and U2 are the amplitudes of the main and auxiliary SAR images respectively; Φ1 and Φ2 are the phases of the main and auxiliary SAR images respectively; * represents the complex conjugate operation of complex numbers; j is a complex unit; after removing the phase component generated by the terrain phase, the differential interferometric image is obtained;

[0016] S13: The differential interferometric images obtained above are screened according to the revisit period of the SAR satellite, and are sorted according to the time sequence to obtain an interferometric image stack; in the stack, the data statistics of the coherence coefficient Coh are performed pixel by pixel, and the maximum value Coh of the pixel temporal coherence coefficient is mainly focused on. higest , the duration TH of each high coherence feature time subset of pixels in the time series i and the complete duration TH of the pixel's high coherence feature whole , the three element features are expressed as:

[0017] ; (2)

[0018] In the formula, Coh higest is the maximum value of the pixel unit coherence coefficient in the entire time series, which is used to preliminarily identify the monitoring potential of the pixel; Threshold1 is the coherence threshold used to determine whether the pixel unit has high coherence characteristics; Time is the revisit period of the SAR satellite; Sum is the number statistics operation; SET j To keep the subset of high coherence features, j∈(1,…,M) is the identification number of the high coherence feature set; M is the total number of high coherence time subsets;

[0019] S14: Combined with the judgment conditions, identify the target pixels with further monitoring potential on the vegetation-covered surface within the map sheet. The judgment conditions are expressed as:

[0020] ; (3)

[0021] In the formula, Time range min Under certain interference combination conditions, the product of the number of images required and Time in the high coherence feature time subset to achieve the corresponding InSAR monitoring accuracy; → indicates that on the basis of satisfying the previous condition, the latter condition needs to be further satisfied;

[0022] Combining the relationship between InSAR monitoring accuracy and the number of differential interferograms, the minimum number of SAR images required in a single high-coherence time subset can be fed back; the formula is expressed as:

[0023] ; (4)

[0024] In the formula, ε vel is the standard error of the InSAR monitoring results, which is modified according to the actual application scenario; C is a constant; N min The minimum number of interferometric combinations required to achieve the expected monitoring accuracy can be inferred from the minimum number of SAR images required in the subset.

[0025] Preferably, S2 comprises the following sub-steps:

[0026] Based on the high potential pixel results fed back in step S1, the research area of ​​interest is divided and cropped to reduce the workload of subsequent entanglement phase unwrapping, singular value decomposition least squares deformation rate solution, and average deformation rate and cumulative deformation solution in the time dimension, thereby reducing additional calculations and speeding up operation efficiency.

[0027] Preferably, S3 includes the following sub-steps:

[0028] S31: For potential areas of interest, the proposed high-coherence time subset construction strategy is used to achieve the division of subsets within the entire time series; the construction strategy mainly includes:

[0029] a. Subset start time T start : The first SAR image is used as the initial start time of the subset, and the subset interval will be added to the start time of subsequent subsets;

[0030] b. The shortest time span of the subset T range : Combined with N in S1 min Parameters and time thresholds to obtain the number of SAR images required within the shortest time span of the high coherence feature subset;

[0031] c. Stop time T of the subset stop :Terminate the expansion of the subset according to the decrease in the number of highly coherent pixels in the subset, so as to maximize the utilization of the interference information of the current subset;

[0032] d. Spacing between subsets T step : In order to reduce the number of subsets generated, the time threshold of the interference combination is used as the interval between the start times of adjacent time subsets;

[0033] S32: Use a text document to record the start time and end time of each subset; the specific production process is: loop steps ac in S31 to achieve the construction of a single high-coherence time subset, and combine with step d to achieve the construction of the time subset in the entire time series.

[0034] Preferably, S4 comprises the following sub-steps:

[0035] S41: calling the unwrapping method of the minimum cost flow to obtain the optimal solution of the whole-image ambiguity of the entanglement phase, that is, obtaining the real differential interferometry phase;

[0036] S42: Using the least square method of singular value decomposition, the average deformation rate and cumulative deformation amount of a single high coherence time subset are calculated, and the formula is:

[0037] ; (5)

[0038] In the formula, M j , N j represents the number of interferograms and SAR images in the jth high coherence time subset, Indicates that the pixel is in M j The interference phase in the interference image; N j The differential interferometric phase of the scene SAR image relative to the initial image; Used to represent interference combinations;

[0039] Pair Matrix Perform singular value decomposition, that is:

[0040] ; (6)

[0041] in the formula is the estimated value of the phase; W is the right singular matrix of matrix A, with dimension M j ×M j ; U T is the transpose of the left singular matrix of matrix A, with dimension N j ×N j ; S -1 are the singular values ​​of matrix A;

[0042] The deformation rate V of each time period in a single time subset is calculated, that is:

[0043] ; (7)

[0044] Where V T represents the transpose of the deformation velocity V, and t represents the imaging date of the SAR image;

[0045] The cumulative deformation within a single time subset j is:

[0046] ; (8)

[0047] In the formula, Def j represents the cumulative deformation of a single time subset j; N j is the number of SAR images in the high coherence time subset;

[0048] The average deformation rate within a single time subset, that is:

[0049] ; (9)

[0050] in the formula represents the average deformation rate in the jth time subset; Time (SET j ) is the time span in the jth time subset.

[0051] Preferably, S5 includes the following sub-steps:

[0052] S51: Based on the weighted method of time dimension, the average deformation rate in the entire time series is obtained in discrete time subsets; the formula is expressed as:

[0053] ; (10)

[0054] In the formula, represents the average deformation rate in the entire time series; j is the sequence number of the high coherence time subset; is the average deformation rate of the jth time subset; t j is the time span of the jth time subset;

[0055] The obtained results are output as GeoTIFF and GRD files of standard scale;

[0056] S52: Based on discrete time subsets, the continuous cumulative deformation in the entire time series is obtained to achieve the goal of high-frequency continuous surface deformation monitoring;

[0057] a. When multiple adjacent high-coherence time subsets overlap in time range, first normalize the deformation of the adjacent subsequent subsets at the starting time, and then use the least squares adjustment method to correct the cumulative deformation in the previous and next high-coherence subsets for the overlapping stage; repeat the above steps to achieve the fusion of the cumulative deformation in the entire time series; the formula is expressed as:

[0058] ; (11)

[0059] In the formula, j represents the serial number of the overlapping subset interval, and i represents the cumulative deformation at time i in the jth subset; Def j,i is the cumulative deformation at the i-th moment in the j-th high coherence time subset; is the estimated value of the cumulative deformation after correction;

[0060] b. When the two adjacent high-coherence time subsets do not overlap in time range, assuming that the deformation rate of the previous stage remains unchanged, the fusion of the cumulative deformation variables of the subsequent subsets is achieved; the formula is expressed as:

[0061] ; (12)

[0062] Def T1 ,Def T2 are the accumulated deformation at time T1 and T2 respectively; is the average deformation rate of the previous subset in the discontinuous time subset;

[0063] The obtained results are output as GeoTIFF and GRD files of standard scale.

[0064] A device used in an InSAR monitoring method for continuous surface deformation under vegetation cover conditions, comprising:

[0065] Data storage unit: used to store the original data of SAR image data, digital elevation model data, precise orbit data and atmospheric delay correction data, as well as the intermediate data in the processing process and the processed InSAR result data;

[0066] Interference processing unit: Based on the interference combination, it realizes the calculation and screening of differential interference images and builds the interference image stack; at the same time, it realizes the calculation of pixel monitoring potential and divides the area of ​​interest according to the spatial clustering method;

[0067] Subset division and solution unit: Combined with the relationship between InSAR monitoring accuracy and the number of differential interferograms and the division strategy of high coherence time subsets, the division of subsets within the entire time series is realized; at the same time, the singular value decomposition method based on the matrix is ​​used to obtain the cumulative deformation amount and average deformation rate results within a single time subset;

[0068] Deformation rate and cumulative deformation fusion unit: use the weighted method of time dimension to calculate the average deformation rate in the entire time series; according to the actual situation, use the least squares adjustment method or conditional assumption method to realize the fusion of cumulative deformation;

[0069] Results output unit: outputs the text document of the above-mentioned high-coherence time subset division results and the average deformation rate and cumulative deformation of the subset, and outputs the average deformation rate and cumulative deformation in the entire time series.

[0070] An electronic device used in an InSAR monitoring method for continuous surface deformation under vegetation cover conditions, the electronic device comprising a data storage device, a data processor and a graphic interaction device;

[0071] The data storage device serves as a data storage unit, realizing corresponding data storage and file transmission and exchange functions;

[0072] The data processor is used to store the computer program and is capable of executing the above-mentioned computer program;

[0073] The graphic interaction device is used to identify high-potential pixels and demarcate regions of interest, smoothly promoting the progress of InSAR processing, and to intuitively observe the InSAR deformation results obtained by the proposed method.

[0074] The advantages of the present invention compared with the prior art are:

[0075] The present invention constructs an InSAR monitoring method for continuous surface deformation under vegetation cover conditions, which can realize fast and high-time-frequency InSAR monitoring of continuous surface deformation in a long time series under a vegetation-covered surface. The method identifies pixel units with deep monitoring potential by constructing a differential interference image stack, combining the maximum value of the pixel time series coherence coefficient, the duration of each high-coherence feature time subset of the pixel in the time series, and the complete duration characteristics of the pixel high-coherence feature; adopts a spatial aggregation method to achieve the division of high-potential and interesting areas; on this basis, combines the relationship between the accuracy of InSAR monitoring results and the number of interference patterns and the proposed high-coherence time subset construction strategy to divide the SAR images in the entire time series; based on the singular value decomposition method, the accumulated deformation amount and the average deformation rate in a single time subset are solved; finally, according to different situations, a discrete subset fusion method in the time dimension is adopted to obtain the average deformation rate of the vegetation-covered surface in the entire time series and the surface deformation results with high time frequency.

[0076] The technical effects of the present invention can be mainly manifested in the following aspects:

[0077] First, the identification of high-potential pixels in the study area and the division of the area of ​​interest reduce the amount of task calculation, avoid unnecessary waste of computing power, and improve the efficiency of obtaining InSAR deformation results.

[0078] Second, compared with the standard SBAS-InSAR method, the proposed method can effectively obtain InSAR monitoring results with a larger coverage area, improve the monitoring potential of the InSAR algorithm, and expand the application scenarios of the InSAR algorithm.

[0079] Third, the proposed fusion method of discrete subsets in the time dimension for different situations enriches the data information that can be provided by the above-mentioned InSAR technology and expands the application scenarios of the InSAR algorithm. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0081] Figure 1 A flow chart of a method for implementing the present invention;

[0082] Figure 2 InSAR monitoring potential map of the study area and pixels;

[0083] Figure 3 Optical images of cropped high potential and research areas of interest;

[0084] Figure 4 Pixel monitoring potential maps for cropped high potential, study areas of interest;

[0085] Figure 5 A deformation result map of the research area obtained according to the inventive method of the present invention;

[0086] Figure 6 This is the deformation result map obtained by standard SBAS-InSAR in the study area;

[0087] Figure 7 To show the instance points of the deformation rate time dimension fusion results, the coherence changes in the entire time series;

[0088] Figure 8 The cumulative deformation amount (coincidence condition) obtained according to the method proposed in the present invention;

[0089] Fig. 9 It is the accumulated deformation amount (non-overlapping condition) obtained according to the method proposed in the present invention. DETAILED DESCRIPTION

[0090] The embodiments of the present application are described in detail below, and examples of the embodiments are shown in the accompanying drawings, wherein the same or similar reference numerals throughout represent the same or similar elements or elements having the same or similar functions. The embodiments described below with reference to the accompanying drawings are exemplary and are only used to explain the present application, and cannot be understood as limiting the present application.

[0091] In the description of this application, it should be noted that, unless otherwise clearly specified and limited, the terms "installed", "connected", and "connected" should be understood in a broad sense, for example, it can be a fixed connection, a detachable connection, or an integral connection; it can be a mechanical connection or an electrical connection; it can be a direct connection or an indirect connection through an intermediate medium, it can be the internal connection of two elements or the interaction relationship between two elements. For ordinary technicians in this field, the specific meanings of the above terms in this application can be understood according to specific circumstances.

[0092] The present invention is further described in detail below in conjunction with the accompanying drawings.

[0093] The present invention provides an InSAR monitoring method for continuous surface deformation under vegetation cover conditions, and the implementation process is as follows: Figure 1 As shown,

[0094] The specific operation process will be given below, and the technical method of the present invention will be further explained and supplemented.

[0095] In this embodiment, the following steps are mainly included:

[0096] Step S1:

[0097] S11: In this embodiment, 141 Sentinel-1SAR image data with track number: 149 and map number: 588 between January 1, 2018 and December 31, 2022 are processed, and the preparation of digital elevation model data, precise orbit data, and atmospheric delay correction data is not described in detail. After removing the SAR images with obvious offset, the Enhanced Spectral Diversity registration method is used to achieve accurate sub-pixel registration of the sequence images. The Enhanced Spectral Diversity registration method is an existing algorithm and can be called directly.

[0098] S12: Set appropriate time dimension baseline and space dimension baseline to generate an interference combination with high quality expectations; call the complex data conjugate multiplication function to obtain the interference image result, the formula is expressed as:

[0099] ; (1)

[0100] Where, U1 and U2 are the amplitudes of the main and auxiliary SAR images respectively; Φ1 and Φ2 are the phases of the main and auxiliary SAR images respectively; * represents the complex conjugate operation of complex numbers; j is a complex unit; after removing the phase component generated by the terrain phase, the differential interferometric image is obtained;

[0101] In this embodiment, the actual application scenario is taken into consideration, and the baseline of the time dimension is 24 days, and the baseline of the space dimension is 500 days.

[0102] S13: The differential interferometric images obtained above are screened according to the revisit period of the SAR satellite, and are sorted according to the time sequence to obtain an interferometric image stack; in the stack, the data statistics of the coherence coefficient Coh are performed pixel by pixel, and the maximum value Coh of the pixel temporal coherence coefficient is mainly focused on. higest , the duration TH of each high coherence feature time subset of pixels in the time series i and the complete duration TH of the pixel's high coherence feature whole , the three element features are expressed as:

[0103] ; (2)

[0104] In the formula, Coh higestis the maximum value of the pixel unit coherence coefficient in the entire time series, which is used to preliminarily identify the monitoring potential of the pixel; Threshold1 is the coherence threshold used to determine whether the pixel unit has high coherence characteristics; Time is the revisit period of the SAR satellite; SET j To keep the subset of high coherent features, j∈(1,…,M) is the identification number of the high coherent feature set; M is the total number of high coherent time subsets.

[0105] In this embodiment, considering the actual experimental results, the value of Threshold1 is 0.1, and Time is the playback period of the Sentinel-1 radar satellite, 12 days.

[0106] S14: Combined with the judgment conditions, target pixels with further monitoring potential on the vegetation-covered surface within the map sheet can be identified. The judgment conditions can be expressed as:

[0107] S14: Combined with the judgment conditions, identify the target pixels with further monitoring potential on the vegetation-covered surface within the map sheet. The judgment conditions are expressed as:

[0108] ; (3)

[0109] In the formula, Time range min Under certain interference combination conditions, in order to achieve the corresponding InSAR monitoring accuracy within the high coherence feature time subset, the product of the number of images required and Time; → indicates that on the basis of satisfying the previous condition, the latter condition must be further satisfied.

[0110] Combining the relationship between InSAR monitoring accuracy and the number of differential interferograms, the minimum number of SAR images required in a single high-coherence time subset can be fed back; the formula is expressed as:

[0111] ; (4)

[0112] In the formula, ε vel is the standard error of the InSAR monitoring results, which is modified according to the actual application scenario; C is a constant; N min The minimum number of interferometric combinations required to achieve the desired monitoring accuracy can be inferred from the minimum number of SAR images required in the subset.

[0113] In this embodiment, ε ve The value is 2mm / year. In addition, the constant C is set to 11, and the minimum number of SAR images in a single high coherence time subset is 16. The obtained full-image high potential pixel results can be seen Figure 2 .

[0114] Step S2:

[0115] Based on the high potential pixel results fed back in step S1, the research area of ​​interest is divided and cropped to reduce the workload of subsequent entanglement phase unwrapping, singular value decomposition least squares deformation rate solution, and average deformation rate and cumulative deformation solution in the time dimension, thereby reducing additional calculations and speeding up operation efficiency.

[0116] The high potential and interesting regions cut out in this embodiment are shown in Figure 3 , the optical image of the area of ​​interest is shown in Figure 4 .

[0117] Step S3:

[0118] S31: For potential areas of interest, the proposed high-coherence time subset construction strategy is used to achieve the division of subsets within the entire time series. The construction strategy mainly includes:

[0119] a. Subset start time T start : The first SAR image is used as the initial start time of the subset, and the subset interval will be added to the start time of subsequent subsets.

[0120] b. The shortest time span of the subset T range :Combined with S1 in N min Parameters and time thresholds are used to obtain the number of SAR images required within the shortest time span of the high coherence feature subset.

[0121] c. Stop time T of the subset stop :Terminate the expansion of the subset according to the decrease in the number of high coherent pixels in the subset to maximize the utilization of the interference information of the current subset.

[0122] d. Spacing between subsets T step : In order to reduce the number of subsets generated, the time threshold of the interference combination is used as the interval between the start times of adjacent time subsets.

[0123] S32: In this embodiment, a text document is used to record the start time and end time of each subset. The specific production process is: looping S31 steps ac to achieve the construction of a single high-coherence time subset, and combined with step d to achieve the construction of the time subset in the entire time series.

[0124] Step S4:

[0125] S41: Call the unwrapping method of the minimum cost flow to obtain the optimal solution of the whole-graph ambiguity of the entanglement phase, that is, to obtain the real differential interferometer phase.

[0126] The phase unwrapping algorithm of the minimum cost flow is an existing algorithm and can be directly called in this embodiment, where the pixel unwrapping threshold is 0.1 and the jump phase threshold is 40.

[0127] S42: Using the least square method of singular value decomposition, the average deformation rate and cumulative deformation amount of a single high coherence time subset j are calculated, and the formula is:

[0128] ; (5)

[0129] In the formula, M j , N j represents the number of interferograms and SAR images in the jth high coherence time subset, Indicates that the pixel is in M j The interference phase in the interference image; N j The differential interferometric phase of the scene SAR image relative to the initial image; Used to represent interference combinations;

[0130] Pair Matrix Perform singular value decomposition, that is:

[0131] ; (6)

[0132] in the formula is the estimated value of the phase; W is the right singular matrix of matrix A, with dimension M j ×M j ; U T is the transpose of the left singular matrix of matrix A, with dimension N j ×N j ; S -1 are the singular values ​​of matrix A;

[0133] The deformation rate V of each time period in a single time subset is calculated, that is:

[0134] ; (7)

[0135] Where V T represents the transpose of the deformation velocity V, and t represents the imaging date of the SAR image;

[0136] The cumulative deformation within a single time subset j is:

[0137] ; (8)

[0138] In the formula, Def j represents the cumulative deformation of a single time subset j; N j is the number of SAR images in the high coherence time subset;

[0139] The average deformation rate within a single time subset, that is:

[0140] ; (9)

[0141] in the formula represents the average deformation rate in the jth time subset; Time (SET j ) is the time span in the jth time subset.

[0142] The least squares algorithm for singular value decomposition is an existing algorithm and can be called directly.

[0143] Step S5:

[0144] S51: Based on the weighted method of the time dimension, the average deformation rate in the entire time series is obtained in discrete time subsets. The formula can be expressed as:

[0145] ; (10)

[0146] In the formula, represents the average deformation rate in the entire time series; j is the sequence number of the high coherence time subset; is the average deformation rate of the jth time subset; t j is the time span of the jth time subset;

[0147] The obtained results are output as GeoTIFF and GRD files of standard scale, which are stored in physical media such as hard disks in this embodiment. The average deformation rate results in the study area obtained in this embodiment are as follows: Figure 5 At the same time, in order to intuitively demonstrate the effectiveness of the method proposed by the present invention, in the same image area and under the same interference conditions, the results obtained based on the standard SBAS-InSAR method are as follows Figure 6 shown.

[0148] S52: Based on discrete time subsets, the continuous cumulative deformation in the entire time series is obtained to achieve the goal of high-frequency continuous surface deformation monitoring.

[0149] a. When multiple adjacent high-coherence time subsets overlap in time range, first normalize the deformation of the adjacent subsequent subset (i.e., the second subset) at the start time, and then use the least squares adjustment method to correct the cumulative deformation in the overlapping stage (before the start time of the third subset). Repeat the above steps to achieve the fusion of the cumulative deformation in the entire time series. The formula can be expressed as:

[0150] ; (11)

[0151] In the formula, j represents the serial number of the overlapping subset interval, and i represents the cumulative deformation at time i in the jth subset; Def j,i is the cumulative deformation at the i-th moment in the j-th high coherence time subset; is the estimated value of the cumulative deformation after correction;

[0152] c. When the two adjacent high-coherence time subsets do not overlap in time range, assuming that the deformation rate of the previous stage remains unchanged, the fusion of the cumulative deformation variables of the subsequent subsets is achieved; the formula is expressed as:

[0153] ; (12)

[0154] Def T1 ,Def T2 are the accumulated deformation at time T1 and T2 respectively; is the average deformation rate of the previous subset in the discontinuous time subset;

[0155] The obtained results are output as GeoTIFF and GRD files of standard image format, which are stored in physical media such as hard disks in this embodiment. In this embodiment, example points are selected for the above two different scenes. The continuous coherence results of pixels in the two different scenes are as follows: Figure 7 As shown in the figure, the cumulative deformation results obtained according to the above method are as follows Figure 8 , Fig. 9 shown.

[0156] The results of the embodiment show:

[0157] First, the standard SBAS-InSAR method has very limited effect on deformation monitoring in the natural subsurface and is often concentrated in stable high-coherence areas, which limits the application scope of InSAR technology.

[0158] Second, in a wide range of natural surfaces, there are intermittent high coherence features in a certain area due to changes in vegetation coverage or human activities. The method proposed in the present invention can effectively realize the full utilization of limited phase information.

[0159] Third, the deformation rate of the study area during 2018-2022 is between -339 mm / year and 239 mm / year. The coverage obtained by the proposed method is 16.4 times that of the standard SBAS-InSAR method. At the same time, the calculation of the area of ​​interest is carried out in a targeted manner, reducing the workload.

[0160] Fourth, the average deformation rate and time series deformation results obtained by the time dimension fusion method further expand the application prospects of InSAR technology.

[0161] In addition to the above algorithm principle itself, this embodiment also includes a device for implementing the above steps:

[0162] Data storage unit: used to store raw data such as SAR image data, digital elevation model data, precise orbit data, and atmospheric delay correction data, as well as intermediate data during the processing process and processed InSAR result data.

[0163] Interference processing unit: Based on the interference combination, it realizes the calculation and screening of differential interference images and builds the interference image stack. At the same time, it realizes the calculation of pixel monitoring potential and divides the area of ​​interest according to the spatial clustering method.

[0164] Subset division and solution unit: Combined with the relationship between InSAR monitoring accuracy and the number of differential interferograms and the division strategy of high coherence time subsets, the division of subsets in the entire time series is realized. At the same time, the singular value decomposition method based on the matrix is ​​used to obtain the cumulative deformation and average deformation rate results within a single time subset.

[0165] Deformation rate and cumulative deformation fusion unit: The weighted method of the time dimension is used to calculate the average deformation rate in the entire time series. According to the actual situation, the least squares adjustment method or conditional assumption method is used to realize the fusion of cumulative deformation.

[0166] Results output unit: outputs the text document of the above-mentioned high-coherence time subset division results and the average deformation rate and cumulative deformation of the subset, and outputs the average deformation rate and cumulative deformation in the entire time series.

[0167] This embodiment also designs an electronic device: the device includes a data storage device, a data processor and a graphic interaction device.

[0168] The data storage device is used as a data storage unit to realize the corresponding data storage and file transfer and exchange functions.

[0169] The data processor is used to store the computer program, and is characterized by being able to execute the above computer program, that is, any one of the above methods.

[0170] The graphic interaction device is used to identify high-potential pixels and demarcate regions of interest, smoothly promoting the progress of InSAR processing, and to intuitively observe the InSAR deformation results obtained by the proposed method.

[0171] The present invention and its embodiments are described above, and such description is not restrictive. The drawings show only one embodiment of the present invention, and the actual structure is not limited thereto. In short, if ordinary technicians in the field are inspired by it, without departing from the purpose of the invention, they can design a structure and embodiment similar to the technical solution without creativity, which should belong to the protection scope of the present invention.

Claims

1. An InSAR monitoring method for continuous surface deformation under vegetation cover conditions, characterized in that: The following steps are involved: S1: Extract the coherence coefficient of the interference pixel from the long time series differential interference image generated by complex conjugate multiplication, combine the maximum value of the pixel time series coherence coefficient, the duration of each high coherence feature time subset of the pixel in the time series, and the complete duration feature of the pixel high coherence feature, to identify the pixel unit with deep monitoring potential; S2: Using the pixel monitoring potential results of S1 as prior knowledge, determine and crop the high potential and interesting regions within the range; S3: Based on the functional relationship between InSAR monitoring accuracy and the number of interferograms and the high-coherence subset construction strategy, the subsets within the entire time series are aggregated and divided; S4: The matrix-based singular value decomposition method is used to obtain the cumulative deformation and average deformation rate results within a single time subset time range; S5: Based on S4, for different scenarios, the discrete subset fusion method in the time dimension is adopted to achieve the average deformation rate of the vegetation-covered surface in the entire time series and the acquisition of high-temporal-frequency surface deformation results.

2. According to the InSAR monitoring method for continuous surface deformation under vegetation cover conditions as described in claim 1, it is characterized by: The S1 comprises the following sub-steps: S11: Perform preliminary inspection on long-time series SAR image data with the same track number and map number, and remove SAR images with obvious offset; use Enhanced Spectral Diversity registration method to achieve accurate sub-pixel registration of sequence images; S12: Set appropriate time dimension baseline and space dimension baseline to generate an interference combination with high quality expectations; call the complex data conjugate multiplication function to obtain the interference image result, the formula is expressed as: (1) In the formula, U 1 and U 2 are the amplitudes of the main and auxiliary SAR images respectively; Φ 1 and Φ 2 are the phases of the main and auxiliary SAR images respectively; * represents the complex conjugate operation of complex numbers; j is a complex unit; after removing the phase component generated by the terrain phase, the differential interferometric image is obtained; S13: The differential interferometric images obtained above are screened according to the revisit period of the SAR satellite, and are sorted according to the time sequence to obtain an interferometric image stack; in the stack, the coherence coefficient is calculated pixel by pixel. Coh Data statistics, focusing on the maximum value of the pixel temporal coherence coefficient Coh higest , the duration of each high coherence feature time subset of pixels in the time series TH i and the full duration of the pixel's high coherence feature TH whole , the three element features are expressed as: (2) In the formula, Coh higest It is the maximum value of the pixel unit coherence coefficient in the entire time series, which is used to preliminarily identify the monitoring potential of the pixel; Threshold 1 is a coherence threshold used to determine whether a pixel unit has a high coherence feature; Time is the revisit period of the SAR satellite; Sum It is a number counting operation; SET j To keep the subset of high coherence features, j∈(1,…,M) is the identification number of the high coherence feature set; M is the total number of high coherence time subsets; S14: Combined with the judgment conditions, identify the target pixels with further monitoring potential on the vegetation-covered surface within the map sheet. The judgment conditions are expressed as: (3) In the formula, Time range min To achieve the corresponding InSAR monitoring accuracy in the high coherence feature time subset under certain interference combination conditions, the number and Time → indicates that on the basis of satisfying the previous condition, the next condition needs to be satisfied; Combining the relationship between InSAR monitoring accuracy and the number of differential interferograms, the minimum number of SAR images required in a single high-coherence time subset can be fed back; the formula is expressed as: (4) In the formula, ε vel It is the standard error of InSAR monitoring results, which is modified according to the actual application scenario; C is a constant; N min The minimum number of interferometric combinations required to achieve the expected monitoring accuracy can be inferred from the minimum number of SAR images required in the subset.

3. According to the InSAR monitoring method for continuous surface deformation under vegetation coverage conditions as described in claim 1, it is characterized by: The S2 comprises the following sub-steps: Based on the high potential pixel results fed back in step S1, the research area of ​​interest is divided and cropped to reduce the workload of subsequent entanglement phase unwrapping, singular value decomposition least squares deformation rate solution, and average deformation rate and cumulative deformation solution in the time dimension, thereby reducing additional calculations and speeding up operation efficiency.

4. The InSAR monitoring method for continuous surface deformation under vegetation coverage conditions according to claim 1, characterized in that: The S3 comprises the following sub-steps: S31: For potential areas of interest, the proposed high-coherence time subset construction strategy is used to achieve the division of subsets within the entire time series; the construction strategy mainly includes: a. Subset start time T start : The first SAR image is used as the initial start time of the subset, and the subset interval will be added to the start time of subsequent subsets; b. Shortest time span of subset T range :Combined with S1 N min Parameters and time thresholds to obtain the number of SAR images required within the shortest time span of the high coherence feature subset; c. Stop time of subset T stop :Terminate the expansion of the subset according to the decrease in the number of highly coherent pixels in the subset, so as to maximize the utilization of the interference information of the current subset; d. Spacing between subsets T step : In order to reduce the number of subsets generated, the time threshold of the interference combination is used as the interval between the start times of adjacent time subsets; S32: Use a text document to record the start time and end time of each subset; the specific production process is: loop steps ac in S31 to achieve the construction of a single high-coherence time subset, and combine with step d to achieve the construction of the time subset in the entire time series.

5. The InSAR monitoring method for continuous surface deformation under vegetation coverage conditions according to claim 1, characterized in that: The S4 comprises the following sub-steps: S41: calling the unwrapping method of the minimum cost flow to obtain the optimal solution of the whole-image ambiguity of the entanglement phase, that is, obtaining the real differential interferometry phase; S42: Using the least square method of singular value decomposition, the average deformation rate and cumulative deformation amount of a single high coherence time subset are calculated, and the formula is: (5) in the formula M j , N j represents the number of interferograms and SAR images in the jth high coherence time subset, Indicates that the pixel is M j The interference phase in the interference image; express N j The differential interferometric phase of the scene SAR image relative to the initial image; Used to represent interference combinations; for matrices Perform singular value decomposition, that is: (6) In the formula is the estimated value of the phase; W is the right singular matrix of matrix A, with dimension M j ×M j ; U T is the transpose of the left singular matrix A, with dimensions N j ×N j ; S -1 are the singular values ​​of matrix A; Calculate the deformation rate of each time period in a single time subset V ,Right now: (7) Among them V T represents the transpose of the deformation velocity V, t Indicates the imaging date of the SAR image; Single time subset j The cumulative deformation within is: (8) In the formula, Def j Represents a single time subset j The cumulative deformation of N j is the number of SAR images in the high coherence time subset; The average deformation rate within a single time subset, that is: (9) In the formula Indicates j The average deformation rate within a time subset; Time(SET j ） For the j The time span within a time subset.

6. The InSAR monitoring method for continuous surface deformation under vegetation coverage conditions according to claim 1, characterized in that: The S5 comprises the following sub-steps: S51: Based on the weighted method of time dimension, the average deformation rate in the entire time series is obtained in discrete time subsets; the formula is expressed as: (10) In the formula, Represents the average deformation rate in the entire time series; j All are sequence numbers of high coherence time subsets; For the j The average deformation rate of a time subset; t j For the j The time span of the time subset; The obtained results are output as GeoTIFF and GRD files of standard scale; S52: Based on discrete time subsets, the continuous cumulative deformation in the entire time series is obtained to achieve the goal of high-frequency continuous surface deformation monitoring; a. When multiple adjacent high-coherence time subsets overlap in time range, first normalize the deformation of the adjacent subsequent subsets at the starting time, and then use the least squares adjustment method to correct the cumulative deformation in the previous and next high-coherence subsets for the overlapping stage; repeat the above steps to achieve the fusion of the cumulative deformation in the entire time series; the formula is expressed as: (11) In the formula, j The sequence number representing the overlapping subset interval, i Representative in the j In subset i The cumulative deformation at the moment; Def j,i For the j The first i The cumulative deformation at the moment; is the estimated value of the cumulative deformation after correction; b. When the two adjacent high-coherence time subsets do not overlap in time range, assuming that the deformation rate of the previous stage remains unchanged, the fusion of the cumulative deformation variables of the subsequent subsets is achieved; the formula is expressed as: (12) Def T1 , Def T2 At the time T 1 , T 2 The cumulative deformation of is the average deformation rate of the previous subset in the discontinuous time subset; The obtained results are output as GeoTIFF and GRD files of standard scale.

7. A device used in the InSAR monitoring method for continuous surface deformation under vegetation cover conditions as claimed in claim 1, characterized in that: include: Data storage unit: used to store the original data of SAR image data, digital elevation model data, precise orbit data and atmospheric delay correction data, as well as the intermediate data in the processing process and the processed InSAR result data; Interference processing unit: Based on the interference combination, it realizes the calculation and screening of differential interference images and builds the interference image stack; at the same time, it realizes the calculation of pixel monitoring potential and divides the area of ​​interest according to the spatial clustering method; Subset division and solution unit: Combined with the relationship between InSAR monitoring accuracy and the number of differential interferograms and the division strategy of high coherence time subsets, the division of subsets within the entire time series is realized; at the same time, the singular value decomposition method based on the matrix is ​​used to obtain the cumulative deformation amount and average deformation rate results within a single time subset; Deformation rate and cumulative deformation fusion unit: use the weighted method of time dimension to calculate the average deformation rate in the entire time series; according to the actual situation, use the least squares adjustment method or conditional assumption method to realize the fusion of cumulative deformation; Results output unit: outputs the text document of the above-mentioned high-coherence time subset division results and the average deformation rate and cumulative deformation of the subset, and outputs the average deformation rate and cumulative deformation in the entire time series.

8. An electronic device used in the InSAR monitoring method for continuous surface deformation under vegetation cover conditions as claimed in claim 1, characterized in that: The electronic device includes a data storage device, a data processor and a graphic interaction device; The data storage device is used as a data storage unit to realize corresponding data storage and file transmission and exchange functions; The data processor is used to store a computer program and is capable of executing the computer program, wherein the computer program is used to execute the InSAR monitoring method for continuous surface deformation according to claim 1; The graphic interaction device is used to identify high-potential pixels and demarcate regions of interest, smoothly promoting the progress of InSAR processing, and to intuitively observe the InSAR deformation results obtained by the proposed method.

Images