Multi-scale adaptive time sequence InSAR (Interferometric Synthetic Aperture Radar) earth surface deformation extraction method

Abstract

The invention discloses a multi-scale adaptive time sequence InSAR (Interferometric Synthetic Aperture Radar) earth surface deformation extraction method, and relates to the field of remote sensing signal processing. The method comprises the following steps of: constructing an unwrapped time sequence InSAR phase simulation data set, performing pixel-by-pixel time sequence empirical mode decomposition on a time domain, obtaining a time sequence intrinsic mode function at the same pixel position, and constructing a FastICA wavelet enhancement denoising algorithm on a space domain: performing two-dimensional multi-stage wavelet decomposition on the intrinsic mode function of each dimension by taking a power spectrum entropy as a cut-off condition, the method comprises the following steps of: acquiring a wavelet approximation coefficient of each level and a detail coefficient, performing noise reduction and dimension reduction processing and reconstruction on the wavelet approximation coefficient and the detail coefficient of each level, and finally performing wavelet reconstruction on the updated wavelet coefficient to obtain a processed time sequence InSAR phase simulation data set. The method provided by the invention is based on data self-driving, self-adaptive decomposition is carried out on phase signals from two scales of time and space, earth surface deformation related signals can be effectively captured, other phase signals can be separated and explained, and a scientific theoretical basis is provided for high-precision earth surface deformation monitoring.

Description

Technical Field

[0001] The invention belongs to the field of remote sensing signal processing, and specifically relates to a multi-scale adaptive time-series InSAR surface deformation extraction method. Technical Background

[0002] Surface deformation refers to the deformation and destruction of the ground morphology due to internal and external dynamic geological effects and human activities, which is mainly manifested as geological disasters such as ground collapse and ground fissures, which are mainly characterized by changes in ground elevation. At present, the mainstream means of surface deformation monitoring include point measurement methods marked by high-precision measuring instruments such as total stations, levels, and GPS. Such methods are flexible in operation and high in accuracy, but have low spatial resolution and unstable data acquisition cycles. Synthetic aperture radar interferometry (InSAR) technology has the ability to work all day and all weather, and has the advantages of high-resolution imaging and wide coverage in surface deformation measurement, which can complement point measurement. Conventional radar differential interferometry (D-InSAR) technology is easily affected by factors such as temporal and spatial decorrelation and atmospheric delay disturbance, and the measurement accuracy is not high. The gradually developed time-series InSAR technology can overcome the above defects to a certain extent and obtain deformation information with millimeter-level accuracy.

[0003] When microwave signals propagate from satellites to ground targets through the atmosphere, they are affected by atmospheric refraction, resulting in time delay and path bending, which will cause signal propagation delay. When acquiring SAR images at different times, different atmospheric conditions (temperature, humidity, atmospheric pressure, etc.) lead to inconsistent propagation delays between images, resulting in an atmospheric phase screen in the interferogram. Satellite orbit parameter errors will cause baseline errors, which will appear in the interferogram in the form of residual fringes. Although precise orbit data is used in the data processing process, completely removing residual fringes requires an absolute orbit accuracy of less than 1 mm, which cannot be met at present. In the processed time-series InSAR surface deformation information, the residual atmospheric phase delay, orbit error and other noise phases jointly restrict InSAR high-precision surface deformation monitoring.

[0004] The Global Navigation Satellite System (GNSS) can provide high-precision ground reference points. Adding GNSS observations before and after SAR data collection is effective in correcting atmospheric delays and orbit errors. In addition, space radiometers, ground meteorological data, and data from medium-resolution imaging spectrometers [Moderate-resolution Imaging Spectroradiometer, MODIS, Medium Resolution Imaging Spectrometer, MERIS] have also been used to weaken atmospheric signals and have made some progress. However, methods such as the above-mentioned external data assistance have limited most of their usage scenarios due to their inconsistency with the spatial resolution of SAR data and the asynchrony of measurement time.

[0005] Correction methods that utilize the inherent characteristics of SAR data, such as atmosphere and elevation modeling, common scene stacking (CSS), permanent scatterer method, and independent component analysis (ICA), are widely used. Among them, ICA, as a statistical method, can achieve good signal separation without establishing a function equation between the observed value and the unknown. However, ICA assumes that the source signals are independent of each other and non-Gaussian distributed. In practical applications, there is often a certain correlation between the source signals, such as the aliasing of the InSAR deformation phase and the atmospheric delay phase. When the frequency components of the source signal are similar, resulting in a low signal-to-noise ratio, the separation effect of ICA is poor, which is also a common defect of ICA when facing complex aliased signals. Therefore, first decomposing the source signal at the time and space levels can better analyze and understand the characteristics of the signal and reduce the complexity, providing more favorable conditions for subsequent ICA processing. Summary of the invention

[0006] In view of the above technical problems, the present invention provides a multi-scale adaptive time-series InSAR surface deformation extraction method. By applying empirical mode decomposition in the time dimension and two-dimensional wavelet decomposition in the space dimension to the time-series InSAR simulation data set, the source signal is input into the FastICA algorithm for analysis according to different time frequencies and spatial frequencies, thereby solving the defects of ICA in processing aliased signals.

[0007] To achieve the above technical objectives, a multi-scale adaptive time-series InSAR surface deformation extraction method of the present invention comprises the following steps:

[0008] On each scene image in the real digital elevation model DEM, the unwrapped time series InSAR phase simulation data set is constructed by simulating the deformation phase, atmospheric delay phase, orbit error phase and random noise phase in the InSAR differential interferometry phase; the deformation phase includes linear deformation phase and periodic deformation phase; the atmospheric delay phase includes layered atmospheric phase and turbulent atmospheric phase;

[0009] Using the empirical mode decomposition method, the image data in the time series InSAR phase simulation data set are decomposed pixel by pixel by time series empirical mode EMD to obtain a series of intrinsic mode functions IMF and trend term components Res components arranged from high to low in terms of time frequency;

[0010] The IMF components and Res components obtained by EMD decomposition of the time series InSAR phase simulation data set are converted into a two-dimensional matrix, and multi-level two-dimensional Sym wavelet decomposition is used to obtain the effective wavelet coefficients of each IMF component and Res component arranged from high to low according to spatial frequency at each decomposition level through power spectrum entropy cutoff threshold determination. The effective wavelet coefficients include approximate coefficients and detail coefficients.

[0011] The FastICA wavelet enhanced denoising algorithm is used to decompose the final effective wavelet coefficients of each IMF component and Res component that are less than the preset cutoff threshold value, with the preset power spectrum entropy value as the cutoff threshold. The final effective wavelet coefficients include the final wavelet approximation coefficients and detail coefficients, and the independent components related to the deformation information, the independent components of the noise term and the corresponding mixing matrix A are obtained, thereby simplifying the complexity of the InSAR phase simulation data set.

[0012] The final effective wavelet coefficients are subjected to denoising and dimensionality reduction processing and reconstruction;

[0013] The FastICA inverse operation method is used to reconstruct the independent components related to the deformation information to obtain the updated final-level effective wavelet coefficients. The final-level effective wavelet coefficients are reconstructed using a two-dimensional wavelet inverse operation to obtain the updated IMF components and Res components. The updated IMF components and Res components are reconstructed to obtain the denoised time-series InSAR surface deformation information.

[0014] Furthermore, the method of constructing the time series InSAR phase simulation data set is as follows: select real DEM data, which covers mountains and plains to meet the needs of atmospheric simulation and terrain changes; simulate linear and periodic sedimentation funnels within the area, and the positions, radii, initial depths and maximum deformation depths of the linear and periodic sedimentation funnels are different;

[0015] The simulated linear deformation function is expressed as:

[0016] def_l=-4*π / λ*d*exp{-[(x-x0). ∧2 +(y-y0). ∧2 ] / r ∧2}*t / T

[0017] Where λ is the radar wavelength; x0 and y0 are the initial coordinates of the deformation center; d is the maximum depth of deformation; r is the deformation radius; t is the current time step; T is the total time step;

[0018] The simulated periodic deformation function is expressed as:

[0019] def_p=-4*π / λ*d*sin(2*π*t / T)*exp{-[(x-x0). ∧2 +(y-y0). ∧2 ] / r ∧2}

[0020] Simulating the stratified atmospheric phase: Assuming that the change in atmospheric density is determined by two factors, time and altitude, a sine function is used to simulate the change in atmospheric density over time. The function is expressed as:

[0021]

[0022] Where t is the current time step; T is the total time step; the function will change periodically between 1 and 1.1 over time;

[0023] Assuming that the atmospheric density decreases exponentially with altitude, the function can be expressed as:

[0024] exp(-1*10 -4 cropped_dem*t)

[0025] Where cropped_dem is the cropped altitude DEM and t is the current time step. As time goes by, the effect of altitude on atmospheric density will become greater.

[0026] Taking into account the two factors of time and altitude, the simulated layered atmosphere function is:

[0027]

[0028] Among them, rh0 is the reference density of the atmosphere; dem_min is the minimum value of the DEM data; H is the atmospheric scale height;

[0029] Simulating turbulent atmospheric phase: Using the self-similarity of fractal functions and the iterative rules of fractal geometry to simulate complex turbulent atmospheric delays, the function is expressed as:

[0030] Hs=(D-1) / 2

[0031] Gs=σ1*randn(N,N)

[0032] [I,J]=meshgrid(1:N,1:N)

[0033] K=|(IJ). ∧ Hs|

[0034] Z = real{ifft2[fft2(Z). * fft2(K)]}

[0035] atm_rn=imgaussfilt(Z,σ2)

[0036] Among them, Hs is the Horst parameter, which is used to describe the autocorrelation and long-term dependence of the time series; D is the fractal dimension. The higher the fractal dimension, the more complex the geometric shape, the rougher the surface, and the stronger the space filling ability; N is the image size; σ is the standard deviation; Gs is the generated Gaussian random variable; randn, meshgrid, ifft2, fft2, and imgaussfilt are the built-in random number generation function, two-dimensional grid generation function, two-dimensional discrete inverse Fourier transform function, two-dimensional discrete Fourier transform function, and two-dimensional Gaussian filter function in MATLAB respectively. The purpose of using imgaussfilt is to smooth the fractal function;

[0037] The simulated orbit error function of the simulated orbit error phase is:

[0038] angle=π / T*sin(t / T*2*π)

[0039] coeff_x = [-0.5*sin(angle)]

[0040] coeff_y=[0.2*cos(angle)]

[0041] orbit=coeff_x+coeff_y

[0042] Where angle is the deviation angle of the orbit error over time; coeff_x is the first-order polynomial defined in the x direction; coeff_y is the first-order polynomial defined in the y direction;

[0043] The simulated random noise function is:

[0044] rn=imnoise('noise',nd)

[0045] Among them, imnoise is the noise adding function in MATLAB; noise is the noise type; nd is the noise density;

[0046] The units of the simulated linear deformation phase, periodic deformation phase, layered atmosphere phase, turbulent atmosphere phase, orbit error phase and random noise phase components are unified into radians and superimposed to obtain the InSAR interferometric phase map after single scene unwrapping:

[0047] Analog_Phase=def_l+def_p+atm_v+atm_rn+orbit+rn

[0048] Each phase component is a function of the current time step t and the total time step T, where T is the number of time series images, thereby obtaining a time series InSAR phase simulation data set.

[0049] Furthermore, the specific process of using empirical mode decomposition to decompose the time series InSAR phase simulation data set pixel by pixel is as follows:

[0050] Convert the time series InSAR phase simulation data set into a two-dimensional matrix AP p*T ={ap1,ap2,...,ap i ,...,ap T}, where ap i The total pixel values ​​corresponding to the i-th scene image are p; T is the number of images, corresponding to the total time step of the time series InSAR phase simulation data set simulated in claim 2;

[0051] For the two-dimensional matrix AP p*T Each row, that is, the same pixel point on each scene image is decomposed by EMD, and it is performed p times in total, so as to achieve pixel-by-pixel decomposition. Each pixel point is decomposed into multiple IMF components and one Res component. The formula is as follows:

[0052]

[0053] Where T is the number of images; m is the current decomposition level; M is the last decomposition level at which iteration stops; at this time, each pixel is arranged according to the time frequency and is decomposed into 1 to M-1 IMF components and 1 Res component; the two-dimensional matrix AP p*T Then there are M-1 IMF components IMF p*T and 1 Res component Res p*T If the pixel point has not been iterated to M-1 times, the subsequent IMF values ​​are filled with zero values; the first IMF p*T Represents the component with the largest frequency change in the time domain in the time series InSAR phase simulation data, to the M-1th IMF p*T The frequency decreases successively; Res p*T Represents the trend of change, with the lowest frequency in the time domain;

[0054] The process of obtaining the intrinsic mode function (IMF) includes: (1) marking the local extreme points of the time-series InSAR phase simulation data; (2) fitting all local maximum points and local minimum points through cubic spline to obtain the upper envelope and lower envelope; (3) using the input signal to subtract the average value of the upper and lower envelopes to obtain the first component h1:

[0055] If component h1 satisfies the following conditions at the same time: (1) in the entire time range, the number of local extreme points and zero crossing points is equal or differs by at most one; (2) at any time point, the average values ​​of the upper envelope and the lower envelope are zero; at this time, h1 is the first IMF component obtained and has the largest frequency; h1 is repeated according to the above steps to obtain the remaining IMF components with gradually decreasing frequencies; when the standard deviation SD k Less than a given threshold θ:

[0056]

[0057] Or when the maximum number of local extreme values ​​in the residual signal is less than the set maximum number, the iteration is stopped and the remaining trend term component Res is obtained; where h k 、h k-1 are the kth and k-1th IMF components respectively; t is the time step, corresponding to each image in the time series InSAR phase simulation data.

[0058] Furthermore, the specific process of applying multi-level two-dimensional wavelet decomposition to each IMF and Res component is as follows:

[0059] IMF p*T Represents the IMF component, with Res p*T Res component, IMF p*T With Res p*T The steps of multi-level two-dimensional wavelet decomposition are the same as those for Res p*T In the multi-level two-dimensional wavelet decomposition of components: first convert the matrix column by column into a two-dimensional matrix for wavelet decomposition, that is, perform wavelet decomposition on the lowest frequency component in the time domain of T images in the spatial domain: convert the i-th column component Res corresponding to the i-th image p*i Convert to a two-dimensional matrix V J , the two-dimensional matrix V J After a wavelet decomposition, we get the J-th level approximation coefficient V J-1 , horizontal detail coefficient Vertical detail factor Diagonal detail factor And the corresponding index matrix I; similarly, for Res p*T When the rest of Ti column vectors are subjected to multi-level two-dimensional wavelet decomposition;

[0060] The power spectrum entropy is used to determine the cutoff level of decomposition and obtain the effective wavelet coefficients. In this process, the length of the effective wavelet coefficients of each level of decomposition varies in Res p*T The weight is the same, with C n / 2j represents the jth level effective wavelet coefficient, where n / 2 j is the effective wavelet coefficient length after J-level wavelet decomposition; the expression for the decomposition cutoff level determined by the power spectrum entropy is:

[0061]

[0062]

[0063] Among them, I original (n) is the input time domain signal, IMF p*T or Res p*T Choose one;e -j2πft It is a complex exponential function, which represents the amplitude and phase information of the signal. The negative sign indicates the counterclockwise rotation of the phase, which is used to capture the phase change. -j2πft The Fourier transform can convert the time domain signal I original (n) is converted into a frequency domain signal X(f); |X(f)| is the amplitude of the frequency domain signal; is the total energy of the power spectrum density; p(f) is the normalized power spectrum density, that is, the proportion of the energy of the frequency domain signal at frequency f to the total energy and its probability distribution; logp(f) represents the logarithm of the energy proportion of the frequency domain signal at frequency f, reflecting the amount of information of the frequency component; p(f)logp(f) represents the contribution of the energy proportion of the frequency domain signal at frequency f to the total entropy, that is, the weighted amount of information; H is the power spectrum entropy. The larger the entropy value, the more complex and random the frequency domain signal is; the smaller the entropy value, the more ordered and stable the signal is.

[0064] Furthermore, Res p*T The final effective wavelet coefficients obtained by component decomposition are used as the data containing the most surface deformation information and FastICA decomposition is applied to them. The independent components obtained after decomposition include surface deformation information and residual noise components. The independent component matrix related to deformation information is selected by calculating the power spectrum entropy of each independent component, comparing the two-dimensional distribution trend of independent component space, and statistically analyzing the time series curve of the mixing matrix A. Finally, the remaining components in other IMFs are decomposed. p*T The deformation information of the effective wavelet coefficients in the middle is processed according to the same idea, and the independent component matrix related to the deformation information is selected; the details are as follows:

[0065] The specific process of decomposing the final effective wavelet coefficients obtained by wavelet decomposition using the FastICA wavelet enhanced denoising algorithm is as follows:

[0066] Extract independent components related to surface deformation information from the final effective wavelet coefficients; before FastICA, p*T and Res p*T The final effective wavelet coefficients obtained by wavelet decomposition are preprocessed to convert the columns of the A matrix into mutually orthogonal ones, while reducing the noise level in the data set; the principal component analysis method (PCA) is used to perform centering and whitening steps on the final effective wavelet coefficients, and the eigenvalue decomposition method (SVD) is used to retain the principal components of the final effective wavelet coefficients. By giving a cumulative threshold of the eigenvalue proportion, all principal components accumulated before this value are retained;

[0067] The formula for FastICA decomposition of the final effective wavelet coefficients is as follows:

[0068] C n / 2j =A*S

[0069] Among them, C n / 2j is the final effective wavelet coefficient of the jth level decomposition, A is the mixing matrix, each column of which represents the time variation of each independent component; S is the source matrix, each row of which represents the spatial characteristics of each independent component.

[0070] Furthermore, the final effective wavelet coefficients (wavelet approximation coefficients and wavelet detail coefficients at all levels) are updated to obtain the updated IMF and Res components. The updated IMF and Res components are summed and reconstructed to obtain the denoised time series InSAR surface deformation information. The specific process is as follows:

[0071] The selected independent component matrix that meets the deformation phase characteristics is reconstructed according to the following formula, that is, the inverse operation of FsatICA is performed to update the wavelet approximation coefficients or detail coefficients:

[0072] C update (f) = A T*n *icasig n*p

[0073] Among them, C update (f) is the updated final effective wavelet coefficient; A T*n is the mixing matrix; icasig n*p is the screened independent component, T is the number of images, corresponding to the time step; n is the number of screened independent components; p corresponds to the number of pixels;

[0074] Reconstruct the updated final wavelet coefficients: reconstruct the final wavelet approximation coefficients and the wavelet detail coefficients at each level:

[0075]

[0076] Among them, [IMFm (n) / Res M (n)] update are the reconstructed IMF components and Res components; j is the number of wavelet decomposition layers; k represents the temporal translation of the wavelet basis function; C update are the updated wavelet approximation coefficients and detail coefficients; ψ j,k,l represents the two-dimensional wavelet basis function; n represents the number of images;

[0077] Remove [IMF m (n) / Res M (n)] update The useless IMF components in the IMF components, that is, in M-1 IMF p*T Among the components, the components whose wavelet coefficient power spectrum entropy exceeds the threshold, the remaining useful IMF components are summed with the Res components to obtain the denoised time series InSAR surface deformation information:

[0078] I update (n) = [IMF m (n)] update,i +[Res M (n)] update

[0079] Among them, I update (n) is the denoised time series InSAR phase simulation data set; [IMF m (n)] update,i To select the IMF components for reconstruction, n is the number of images, corresponding to the time step.

[0080] Beneficial Effects

[0081] The present invention provides a multi-scale adaptive time-series InSAR surface deformation extraction method, which is easy to implement and has significant effects, mainly reflected in the following four points:

[0082] First, by constructing a simulated data set, compared with real InSAR data, on the one hand, a large amount of data can be generated quickly, which is convenient for rapid iteration and optimization of the algorithm; on the other hand, various influencing factors can be accurately controlled, including linear / periodic deformation, stratified / turbulent atmosphere, orbit error and random noise, so as to better understand the impact of these factors on the phase signal.

[0083] Second, the present invention uses empirical mode decomposition to decompose the data set pixel by pixel, and constructs a time series with the same pixel spatial coordinates and different time steps. This method does not require a preset basis function, and adaptively decomposes the time series of each pixel into several IMF and Res components representing different frequency components, thereby achieving the separation of the original mixed phase in time and frequency.

[0084] Third, the present invention uses two-dimensional wavelet decomposition to perform two-dimensional spatial decomposition of IMF and Res components, and obtains the corresponding spatial high-frequency to low-frequency wavelet coefficients step by step in the order of high frequency to low frequency in time. Each level of decomposition improves the frequency resolution, obtains lower frequency information, corresponding to the approximate coefficient, and higher frequency information, corresponding to the detail coefficient, and further realizes the separation of the original mixed phase in spatial frequency.

[0085] Fourth, the present invention uses the FastICA algorithm to analyze the extracted wavelet coefficients of different time frequencies and different spatial frequencies. Compared with directly using the FastICA algorithm on the original aliased signal, its input signal has better separability in the frequency domain, which effectively solves the defect of ICA in dealing with poor separation effect of non-Gaussian and non-independent signals. BRIEF DESCRIPTION OF THE DRAWINGS

[0086] Figure 1 The present invention is a flowchart of the multi-scale adaptive time-series InSAR surface deformation extraction method.

[0087] Figure 2 are the phase components and the time-series InSAR phase simulation data sets simulated in the embodiments of the present invention, wherein a is the deformation phase time-series InSAR phase simulation data set, b is the layered atmosphere phase time-series InSAR phase simulation data set, c is the turbulent atmosphere phase time-series InSAR phase simulation data set, d is the orbit error phase time-series InSAR phase simulation data set, and e is the time-series InSAR phase simulation data set after the deformation phase, layered atmosphere phase, turbulent atmosphere phase, and orbit error phase are superimposed.

[0088] Figure 3 It is the denoised time-series InSAR surface deformation information obtained after applying empirical wavelet decomposition and ICA to the time-series InSAR simulation data set in an embodiment of the present invention. DETAILED DESCRIPTION

[0089] The following is a further description in conjunction with the embodiments of the present invention.

[0090] if Figure 1 As shown, the present invention discloses a multi-scale adaptive time-series InSAR surface deformation extraction method, the specific steps are as follows:

[0091] Step 1. Build a time series InSAR phase simulation dataset:

[0092] (1) Simulate linear / periodic deformation phase, such as Figure 2 As shown in (a):

[0093] The DEM data covering Reunion Island (20.87°S-21.39°S, 55.21°E-55.84°E) was selected. Because the island has a tropical rainforest climate and the Fournaise volcano is located in the east, the high water vapor content and mountainous conditions throughout the year make the area have obvious atmospheric delay signals, which meet the phase simulation conditions. By intercepting 1000×1000 pixel data, three linear changes and one periodic change sedimentation funnel are simulated, and the simulated linear deformation function is expressed as:

[0094] def_l=-4*π / λ*d*exp{-[(x-x0). ∧2 +(y-y0). ∧2 ] / r ∧2}*t / T

[0095] The initial coordinates x0 and y0 of the linear sedimentation funnel are (720, 530) (500, 300) (250, 750), respectively; the radius r is 100 m, 50 m, and 150 m, respectively; and the maximum sedimentation d is 0.1 m, 0.05 m, and 0.2 m, respectively.

[0096] The simulated periodic deformation function is expressed as:

[0097] def_p=-4*π / λ*d*sin(2*π*t / T)*exp{-[(x-x0). ∧2 +(y-y0). ∧2 ] / r ∧2}

[0098] Among them, the initial coordinates x0 and y0 of the periodic sedimentation funnel are (750, 750), the radius r is 50m, and the maximum sedimentation d is 0.1m.

[0099] The present invention sets the radar wavelength d to 0.056m, that is, simulates the C-band radar signal; sets each step length t to represent 12 days, and the total time step length T is 40, that is, the linear deformation reaches the maximum deformation within 480 days, and the periodic deformation changes with time.

[0100] (2) Simulating the stratified atmospheric phase, such as Figure 2 As shown in (b):

[0101] The present invention assumes that the change of atmospheric density is determined by two factors, time and altitude. Taking these two factors into account, the simulated layered atmosphere function is:

[0102]

[0103] Among them, the atmospheric reference density rh0 is set to 1.225; cropped_dem is the elevation of each pixel after DEM is cropped, dem_min is the minimum value of DEM data, which is automatically obtained by the code; H is the atmospheric scale height, which is set to 8000.

[0104] (3) Simulate turbulent atmospheric phases, such as Figure 2 As shown in (c):

[0105] Based on the self-similarity of fractal functions and the iterative rules of fractal geometry, it is used to simulate complex turbulent atmospheric delays. The function is expressed as:

[0106] Hs=(D-1) / 2

[0107] Gs=σ1*randn(N,N)

[0108] [I,J]=meshgrid(1:N,1:N)

[0109] K=|(IJ). ∧ Hs|

[0110] Z = real{ifft2[fft2(Z). * fft2(K)]}

[0111] atm_rn=imgaussfilt(Z,σ2)

[0112] Among them, the fractal dimension D is set to 3, and the Horst parameter Hs is obtained to be 1; N is the image size, that is, 1000; the standard deviation σ1 is set to 1, and Gaussian random variables are automatically generated by the code; randn, meshgrid, ifft2, fft2, and imgaussfilt are the built-in random number generation function, two-dimensional grid generation function, two-dimensional discrete inverse Fourier transform function, two-dimensional discrete Fourier transform function, and two-dimensional Gaussian filter function in MATLAB, respectively, and the filter function standard deviation σ2 is set to 1.8.

[0113] (4) Simulate orbit error phase, such as Figure 2 As shown in (d):

[0114] The present invention uses a first-order polynomial to simulate the orbit error and make its angle change slightly over time. The simulated orbit error function is:

[0115] angle=π / T*sin(t / T*2*π)

[0116] coeff_x = [-0.5*sin(angle)]

[0117] coeff_y=[0.2*cos(angle)]

[0118] orbit=coeff_x+coeff_y

[0119] The orbit error deviation angle angle is determined by t and T, and the first-order polynomial coeff_x defined in the x direction and the first-order polynomial coeff_y defined in the y direction are further obtained.

[0120] (5) Simulating random noise

[0121] In order to reflect the high random characteristics of random noise in space and time, the present invention uses MATLAB built-in salt and pepper noise to simulate random noise:

[0122] rn=imnoise('salt&pepper',nd)

[0123] Among them, imnoise is the noise adding function in MATLAB; the noise density nd is set to 0.05.

[0124] The units of the simulated 4 types and 6 phase components are unified into radians and superimposed to obtain the InSAR interferometric phase map after single scene unwrapping, namely:

[0125] Analog_Phase=def_l+def_p+atm_v+atm_rn+orbit+rn

[0126] Each phase component is a function of "current time step t and total time step T". The first scene image is taken as the main image, and the maximum simulated time baseline is 480 days (12*T, T=40). A single main time series interferogram is obtained, that is, the time series InSAR phase simulation data set, such as Figure 2 As shown in (e):

[0127] Step 2. Use empirical mode decomposition to decompose the time series InSAR phase simulation data set pixel by pixel:

[0128] First, convert the data set into a two-dimensional matrix AP 1000000*40 ={ap1,ap2,...,ap i ,...,ap 40}, where ap i All pixel values ​​corresponding to the i-th image, the number is 1*10 6 .

[0129] In this example, for the two-dimensional matrix AP 1000000*40 Apply empirical mode decomposition to limit the maximum number of extreme values ​​of the residual signal to 1, according to the formula:

[0130]

[0131] From high frequency to low frequency in the time dimension, different numbers of IMF components and 1 Res component are obtained, among which the maximum M value of the decomposed IMF components is 8 and the minimum is 3. The 8 IMF components are screened, and the 5th to 8th groups of IMF components show most of the random noise points and no useful signals, so they are eliminated (this step is based on the fact that the power spectrum entropy of the wavelet coefficients in step 3 exceeds the set threshold).

[0132] Step 3. Apply two-dimensional multi-level wavelet decomposition to each IMF and Res component:

[0133] Set AP of 1-4 IMF components 1000000*40 Each column vector of the matrix is ​​converted into 40 two-dimensional matrices V 1000*1000 , decomposed using Sym6 wavelet (Sym6 wavelet is a type of Symlets wavelet, has approximate symmetry characteristics, can effectively avoid phase distortion; has 6th-order vanishing moment, can better represent complex polynomial signals; has compact support, helps to reduce boundary effects; has good regularity, reduces the discontinuity of the reconstructed image).

[0134] After applying Sym6 wavelet decomposition to each layer, the two-dimensional matrix is ​​converted into a one-dimensional vector and the power spectrum entropy is calculated based on the index of the I matrix to the dimension of the wavelet approximation coefficient matrix of the current layer; for the detail coefficients in the three directions, they are converted into one-dimensional vectors, added and combined, and the power spectrum entropy is calculated. In this example, the normalized power spectrum entropy threshold is set to 0.75. After 5-level decomposition, the average values ​​of the power spectrum entropy of the wavelet coefficients of the first group of IMF components are 0.27, 0.36, 0.45, 0.46, and 0.57, respectively; after 6-level decomposition, the average values ​​of the power spectrum entropy of the wavelet coefficients of the second group of IMF components are 0.14, 0.21, 0.22, 0.29, 0.58, and 0.73, respectively; after 4-level decomposition, the average values ​​of the power spectrum entropy of the wavelet coefficients of the third group of IMF components are 0.13, 0.36, 0.45, and 0.67, respectively; after 4-level decomposition, the average values ​​of the power spectrum entropy of the wavelet coefficients of the fourth group of IMF components are 0.24, 0.42, 0.41, and 0.55, respectively; after 5-level decomposition, the average values ​​of the power spectrum entropy of the wavelet coefficients of the Res component are 0.11, 0.21, 0.34, 0.52, and 0.68, respectively. Furthermore, the approximate coefficient matrices of the four groups of IMF components and Res components are obtained, which are: Similarly, the horizontal, vertical, diagonal and merged detail coefficient matrices of the corresponding levels are obtained, and the dimensions are the same as the approximate coefficients, which will not be repeated here.

[0135] Step 4. Use the FastICA algorithm to decompose the nearly effective wavelet coefficients obtained by wavelet decomposition:

[0136] The cumulative proportion of eigenvalues ​​is set to 90%, and the power spectrum entropy threshold is set to 0.7. First, the final approximate coefficient matrix of the Res component is processed. After PCA, the proportion of the eigenvalue of the first principal component is 96.34%, so only the first principal component is retained for ICA to obtain the corresponding independent component. By analyzing its spatial pattern, it is almost completely consistent with the three linear sedimentation funnels simulated in the time series InSAR simulation data set. The column vector corresponding to its A matrix (actually 40*1 columns due to only one independent component) shows a linear change trend over time. In addition, the power spectrum entropy mean of this component is 0.28, which is far below the threshold, and there are obvious peaks at sampling points 3, 14, and 20, indicating that there are three "deterministic" signals in the entire signal. Therefore, the independent component is determined to be a linear surface deformation related component.

[0137] For the final approximate coefficient of the IMF1 component, after PCA, the eigenvalue of the first principal component accounts for 41.59%, the second principal component accounts for 15.88%, and the 12th principal component accounts for 90.35% of the eigenvalue. ICA is performed on these 12 principal components, and the number of independent components is set to 12 (for the subsequent number of independent components, it is uniformly set to the same number as the number of principal components), and the decomposition results are obtained. By analyzing the spatial pattern, the sixth independent component is almost completely consistent with the simulated periodic sedimentation funnel, and the column vector corresponding to its A matrix conforms to the change trend of the descending part of the sine function, and the mean value of the power spectrum entropy is 0.25, and there is an obvious peak at sampling point 22. Therefore, it is determined that the independent component is a component related to periodic surface deformation. Among the remaining 11 independent components, they all show the characteristics of random distribution in the spatial pattern, which is consistent with the simulated turbulent atmosphere. The column vector corresponding to the A matrix has no correlation with time, and the mean value of the power spectrum entropy is generally between 0.75-0.9. Therefore, it is preliminarily determined that these independent components are delayed components of the turbulent atmosphere.

[0138] For the final approximate coefficient of the IMF2 component, after PCA, the eigenvalue of the first principal component accounts for 83.37%, the second principal component accounts for 8.63%, and the cumulative eigenvalue accounts for 92.00%. ICA is performed on these two principal components. By analyzing the spatial pattern, the first independent component shows the edge of a linear sedimentation funnel, and the corresponding column vector of the A matrix is ​​also roughly correlated with time. The mean power spectrum entropy is 0.54, and there is a significant peak at sampling point 10. It is determined that this independent component is a supplement to the linear surface deformation related component; the second independent component is almost completely consistent with the simulated periodic sedimentation funnel, and the corresponding column vector of the A matrix conforms to the trend of the rising part of the sine function. The mean power spectrum entropy is 0.48, and there is a significant peak at sampling point 22, which is determined to be a periodic surface deformation related component.

[0139] For the final approximate coefficient of the IMF3 component, after PCA, the eigenvalue of the first principal component accounts for 73.70%, and the eigenvalue of the fourth principal component accounts for 92.41%. ICA is performed on these four principal components. By analyzing the spatial pattern, the first independent component is almost consistent with the simulated periodic sedimentation funnel. The column vector corresponding to its A matrix conforms to the periodic change trend of the sine function. The mean value of the power spectrum entropy is 0.52, and there is an obvious peak at sampling point 22. Therefore, it is determined that the independent component is a component related to periodic surface deformation. The remaining three independent components all show characteristics consistent with the turbulent atmosphere.

[0140] For the IMF4 component, from the first to the last approximate coefficients, after PCA dimensionality reduction and ICA, none of the independent components showed a spatial pattern similar to the simulated linear / periodic sedimentation funnel, showing a mixture of a large number of random noise points and an approximate turbulent atmosphere, and was determined to be a noise component.

[0141] It is worth mentioning that the layered atmosphere is reflected in the Res component. Specifically: after PCA decomposition, the second principal component eigenvalue of the final approximate coefficient of the Res component accounts for 3.4%. After ICA, it is found that among the two independent components, one is the aforementioned linear surface deformation component, and the other is almost consistent with the simulated layered atmosphere in terms of spatial pattern. The column vector corresponding to its A matrix shows a linear trend with time and has a strong correlation. The mean power spectrum entropy is 0.65, and there is no obvious peak, indicating that the overall uncertainty of this independent component is strong. The orbit error is in the IMF2 component. Specifically: after PCA decomposition of the third-level approximate coefficient of the IMF2 component, the third and fourth principal components are retained for ICA. The other two independent components of the four independent components are consistent with the simulated orbit error in terms of spatial pattern. The column vector corresponding to its A matrix has no correlation with time. The mean power spectrum entropy is 0.57 and 0.64, respectively, without obvious peaks and flat, indicating that the signal spectrum is evenly distributed.

[0142] For the detail coefficients of the four groups of IMF components and Res components, this example decomposes them according to the level of ICA decomposition of the corresponding approximate coefficients. The detail coefficients correspond to spatial high-frequency information. By analyzing the spatial pattern of each group of component detail coefficients, the four groups of IMF components all contain high-frequency noise independent components, that is, they show dense noise points, and the power spectrum entropy of the wavelet coefficients exceeds the threshold. The relevant components are eliminated, and only the remaining components are retained to participate in the reconstruction of the wavelet detail coefficients.

[0143] Step 5. Update the final effective wavelet coefficients (wavelet approximation coefficients and wavelet detail coefficients at all levels) to obtain the updated IMF and Res components. Sum and reconstruct the updated IMF and Res components to obtain the denoised time series InSAR surface deformation information:

[0144] Multiply the mixing matrix A corresponding to the independent component by the independent component icasig, and reconstruct according to the following formula to obtain the updated wavelet coefficients (taking the independent component that meets the linear deformation phase characteristics selected in step 4 as an example):

[0145] C update (f) = A 40*1 *icasig 1*1681

[0146] It should be noted that the independent components subsequently reconstructed are all components related to linear / periodic surface deformation.

[0147] All updated wavelet approximation coefficients are reconstructed from the last level (i.e., the level at which ICA is run), and updated wavelet detail coefficients are reconstructed from all levels (i.e., all lower levels before the level at which ICA is run). Specifically, the waverec2 function in MATLAB, C++, and update (f) and the index matrix I to obtain the reconstructed [IMF(40)] update,1,2,3 Components and Res(40)] update .

[0148] Sum the reconstructed IMF component and the Res component, that is:

[0149] I update (40) = [IMF (40)] update,1,2,3 +[Res(40)] update

[0150] The denoised time series InSAR surface deformation information is obtained, such as Figure 3 shown.

Claims

1. A multi-scale adaptive time-series InSAR surface deformation extraction method, characterized in that: Here are the steps: On each scene image in the real digital elevation model DEM, the deformation phase, atmospheric delay phase, orbit error phase and random noise phase in the InSAR differential interferometry phase are simulated to construct an unwrapped time-series InSAR phase simulation data set. The deformation phase includes linear deformation phase and periodic deformation phase; The atmospheric delay phase includes the stratified atmospheric phase and the turbulent atmospheric phase; Using the empirical mode decomposition method, the image data in the time series InSAR phase simulation data set are decomposed pixel by pixel by time series empirical mode EMD to obtain a series of intrinsic mode functions IMF and trend term components Res components arranged from high to low in terms of time frequency; The IMF components and Res components obtained by EMD decomposition of the time series InSAR phase simulation data set are converted into a two-dimensional matrix, and multi-level two-dimensional Sym wavelet decomposition is used to obtain the effective wavelet coefficients of each IMF component and Res component arranged from high to low according to spatial frequency at each decomposition level through power spectrum entropy cutoff threshold determination. The effective wavelet coefficients include approximate coefficients and detail coefficients. The FastICA wavelet enhanced denoising algorithm is used to decompose the final effective wavelet coefficients of each IMF component and Res component that are less than the preset cutoff threshold value, with the preset power spectrum entropy value as the cutoff threshold. The final effective wavelet coefficients include the final wavelet approximation coefficients and detail coefficients, and the independent components related to the deformation information, the independent components of the noise term and the corresponding mixing matrix A are obtained, thereby simplifying the complexity of the InSAR phase simulation data set. The final effective wavelet coefficients are subjected to denoising and dimensionality reduction processing and reconstruction; The FastICA inverse operation method is used to reconstruct the independent components related to the deformation information to obtain the updated final-level effective wavelet coefficients. The final-level effective wavelet coefficients are reconstructed using a two-dimensional wavelet inverse operation to obtain the updated IMF components and Res components. The updated IMF components and Res components are reconstructed to obtain the denoised time-series InSAR surface deformation information.

2. The multi-scale adaptive time-series InSAR surface deformation extraction method according to claim 1, characterized in that: The method of constructing the time series InSAR phase simulation data set is as follows: select real DEM data, which covers mountains and plains to meet the needs of atmospheric simulation and terrain changes; simulate linear and periodic sedimentation funnels within the area, and the positions, radii, initial depths and maximum deformation depths of the linear and periodic sedimentation funnels are different; The simulated linear deformation function is expressed as: def_l=-4*π / λ*d*exp{-[(x-x0). ∧2 +(y-y0). ∧2 ] / r ∧2 }*t / T Where λ is the radar wavelength; x0 and y0 are the initial coordinates of the deformation center; d is the maximum depth of deformation; r is the deformation radius; t is the current time step; T is the total time step; The simulated periodic deformation function is expressed as: def_p=-4*π / λ*d*sin(2*π*t / T)*exp{-[(x-x0). ∧2 +(y-y0). ∧2 ] / r ∧2 } Simulating the stratified atmospheric phase: Assuming that the change in atmospheric density is determined by two factors, time and altitude, a sine function is used to simulate the change in atmospheric density over time. The function is expressed as: Where t is the current time step; T is the total time step; the function will change periodically between 1 and 1.1 over time; Assuming that the atmospheric density decreases exponentially with altitude, the function can be expressed as: exp(-1*10 -4 cropped_dem*t) Where cropped_dem is the cropped altitude DEM and t is the current time step. As time goes by, the effect of altitude on atmospheric density will become greater. Taking into account the two factors of time and altitude, the simulated layered atmosphere function is: Among them, rh0 is the reference density of the atmosphere; dem_min is the minimum value of the DEM data; H is the atmospheric scale height; Simulating turbulent atmospheric phase: Using the self-similarity of fractal functions and the iterative rules of fractal geometry to simulate complex turbulent atmospheric delays, the function is expressed as: Hs=(D-1) / 2 Gs=σ1*randn(N,N) [I,J]=meshgrid(1:N,1:N) K=|(IJ). ∧ Hs| Z=rdal{ifft2[fft2(Z). * fft2(K)]} atm_rn=imgaussfilt(Z,σ2) Among them, Hs is the Horst parameter, which is used to describe the autocorrelation and long-term dependence of the time series; D is the fractal dimension. The higher the fractal dimension, the more complex the geometric shape, the rougher the surface, and the stronger the space filling ability; N is the image size; σ is the standard deviation; Gs is the generated Gaussian random variable; randn, meshgrid, ifft2, fft2, and imgaussfilt are the built-in random number generation function, two-dimensional grid generation function, two-dimensional discrete inverse Fourier transform function, two-dimensional discrete Fourier transform function, and two-dimensional Gaussian filter function in MATLAB respectively. The purpose of using imgaussfilt is to smooth the fractal function; The simulated orbit error function of the simulated orbit error phase is: angle=π / T*sin(t / T*2*π) coeff_x = [-0.5*sin(angle)] coeff_y=[0.2*cos(angle)] orbit=coeff_x+coeff_y Where angle is the deviation angle of the orbit error over time; coeff_x is the first-order polynomial defined in the x direction; coeff_y is the first-order polynomial defined in the y direction; The simulated random noise function is: rn=imnoise('noise',nd) Among them, imnoise is the noise adding function in MATLAB; noise is the noise type; nd is the noise density; The units of the simulated linear deformation phase, periodic deformation phase, layered atmosphere phase, turbulent atmosphere phase, orbit error phase and random noise phase components are unified into radians and superimposed to obtain the InSAR interferometric phase map after single scene unwrapping: Analong_Phase=def_l+def_p+atm_v+atm_rn+orbit+rn Each phase component is a function of the current time step t and the total time step T, where T is the number of time series images, thereby obtaining a time series InSAR phase simulation data set.

3. The multi-scale adaptive time-series InSAR surface deformation extraction method according to claim 1, characterized in that: The specific process of pixel-by-pixel decomposition of the time-series InSAR phase simulation data set using empirical mode decomposition is as follows: Convert the time series InSAR phase simulation data set into a two-dimensional matrix AP p*T ={ap1,ap2,...,ap i ,...,ap T }, where ap i The total pixel values ​​corresponding to the i-th image are p; T is the number of images, corresponding to the total time step of the time series InSAR phase simulation data set simulated in claim 2; For the two-dimensional matrix AP p*T Each row, that is, the same pixel point on each scene image is decomposed by EMD, and it is performed p times in total, so as to achieve pixel-by-pixel decomposition. Each pixel point is decomposed into multiple IMF components and one Res component. The formula is as follows: Where T is the number of images; m is the current decomposition level; M is the last decomposition level at which iteration stops; at this time, each pixel is arranged according to the time frequency and is decomposed into 1 to M-1 IMF components and 1 Res component; the two-dimensional matrix AP p*T Then there are M-1 IMF components IMF p*T and 1 Res component Res p*T If the pixel point has not been iterated to M-1 times, the subsequent IMF values ​​are filled with zero values; the first IMF p*T Represents the component with the largest frequency change in the time domain in the time series InSAR phase simulation data, to the M-1th IMF p*T The frequency decreases successively; Res p*T Represents the trend of change, with the lowest frequency in the time domain; The process of obtaining the intrinsic mode function (IMF) includes: (1) marking the local extreme points of the time-series InSAR phase simulation data; (2) fitting all local maximum points and local minimum points through cubic spline to obtain the upper envelope and lower envelope; (3) using the input signal to subtract the average value of the upper and lower envelopes to obtain the first component h1: If component h1 satisfies the following conditions at the same time: (1) in the entire time range, the number of local extreme points and zero crossing points is equal or differs by at most one; (2) at any time point, the average values ​​of the upper envelope and the lower envelope are zero; at this time, h1 is the first IMF component obtained and has the largest frequency; h1 is repeated according to the above steps to obtain the remaining IMF components with gradually decreasing frequencies; when the standard deviation SD k Less than a given threshold θ: Or when the maximum number of local extreme values ​​in the residual signal is less than the set maximum number, the iteration is stopped and the remaining trend term component Res is obtained; where h k 、h k-1 are the kth and k-1th IMF components respectively; t is the time step, corresponding to each image in the time series InSAR phase simulation data.

4. The multi-scale adaptive time-series InSAR surface deformation extraction method according to claim 3, characterized in that: The specific process of applying multi-level two-dimensional wavelet decomposition to each IMF and Res component is as follows: IMF p*T Represents the IMF component, with Res p*T Res component, IMF p*T With Res p*T The steps of multi-level two-dimensional wavelet decomposition are the same as those for Res p*T In the multi-level two-dimensional wavelet decomposition of components: first convert the matrix column by column into a two-dimensional matrix for wavelet decomposition, that is, perform wavelet decomposition on the lowest frequency component in the time domain of T images in the spatial domain: convert the i-th column component Res corresponding to the i-th image p*i Convert to a two-dimensional matrix V J , the two-dimensional matrix V J After a wavelet decomposition, we get the J-th level approximation coefficient V J-1 , horizontal detail coefficient Vertical detail factor Diagonal detail factor And the corresponding index matrix I; similarly, for Res p*T When the rest of Ti column vectors are subjected to multi-level two-dimensional wavelet decomposition; The power spectrum entropy is used to determine the cutoff level of decomposition and obtain the effective wavelet coefficients. In this process, the length of the effective wavelet coefficients of each level of decomposition varies in Res p*T The weights are the same, represents the jth level effective wavelet coefficient, where n / 2 j is the effective wavelet coefficient length after J-level wavelet decomposition; the expression for the decomposition cutoff level determined by the power spectrum entropy is: Among them, I original (n) is the input time domain signal, IMF p*T or Res p*T ;e -j2πft It is a complex exponential function, which represents the amplitude and phase information of the signal. The negative sign indicates the counterclockwise rotation of the phase, which is used to capture the phase change. -j2πft The Fourier transform can convert the time domain signal I original (n) is converted into a frequency domain signal X(f); |X(f)| is the amplitude of the frequency domain signal; is the total energy of the power spectrum density; p(f) is the normalized power spectrum density, that is, the proportion of the energy of the frequency domain signal at frequency f to the total energy and its probability distribution; logp(f) represents the logarithm of the energy proportion of the frequency domain signal at frequency f, reflecting the amount of information of this frequency component; p(f)logp(f) represents the contribution of the energy proportion of the frequency domain signal at frequency f to the total entropy, that is, the weighted amount of information; H is the power spectrum entropy. The larger the entropy value, the more complex and random the frequency domain signal is; the smaller the entropy value, the more ordered and stable the signal is.

5. The multi-scale adaptive time-series InSAR surface deformation extraction method according to claim 4, characterized in that: Res p*T The final effective wavelet coefficients obtained by component decomposition are used as the data containing the most surface deformation information and FastICA decomposition is applied to them. The independent components obtained after decomposition include surface deformation information and residual noise components. The independent component matrix related to deformation information is selected by calculating the power spectrum entropy of each independent component, comparing the two-dimensional distribution trend of independent component space, and statistically analyzing the time series curve of the mixing matrix A. Finally, the remaining components in other IMFs are decomposed. p*T The deformation information of the effective wavelet coefficients in the middle is processed according to the same idea, and the independent component matrix related to the deformation information is selected; the details are as follows: The specific process of decomposing the final effective wavelet coefficients obtained by wavelet decomposition using the FastICA wavelet enhanced denoising algorithm is as follows: Extract independent components related to surface deformation information from the final effective wavelet coefficients; before FastICA, p*T and Res p*T The final effective wavelet coefficients obtained by wavelet decomposition are preprocessed to convert the columns of the A matrix into mutually orthogonal ones, while reducing the noise level in the data set; the principal component analysis method (PCA) is used to perform centering and whitening steps on the final effective wavelet coefficients, and the eigenvalue decomposition method (SVD) is used to retain the principal components of the final effective wavelet coefficients. By giving a cumulative threshold of the eigenvalue proportion, all principal components accumulated before this value are retained; The formula for FastICA decomposition of the final effective wavelet coefficients is as follows: in, is the final effective wavelet coefficient of the jth level decomposition, A is the mixing matrix, each column of which represents the time variation of each independent component; S is the source matrix, each row of which represents the spatial characteristics of each independent component.

6. The multi-scale adaptive time-series InSAR surface deformation extraction method according to claim 5, characterized in that: Update the final effective wavelet coefficients (wavelet approximation coefficients and wavelet detail coefficients at all levels) to obtain the updated IMF and Res components. Sum and reconstruct the updated IMF and Res components to obtain the denoised time series InSAR surface deformation information. The specific process is as follows: The selected independent component matrix that meets the deformation phase characteristics is reconstructed according to the following formula, that is, the inverse operation of FsatICA is performed to update the wavelet approximation coefficients or detail coefficients: C update (f)=A T*n *icasig n*p Among them, C update (f) is the updated final effective wavelet coefficient; A T*n is the mixing matrix; icasig n*p is the screened independent component, T is the number of images, corresponding to the time step; n is the number of screened independent components; p corresponds to the number of pixels; Reconstruct the updated final wavelet coefficients: reconstruct the final wavelet approximation coefficients and the wavelet detail coefficients at each level: Among them, [IMF m (n) / Res M (n)] update are the reconstructed IMF components and Res components; j is the number of wavelet decomposition layers; k represents the temporal translation of the wavelet basis function; C update are the updated wavelet approximation coefficients and detail coefficients; ψ j,k,l represents the two-dimensional wavelet basis function; n represents the number of images; Remove [IMF m (n) / Res M (n)] update The useless IMF components in the IMF components, that is, in M-1 IMF p*T Among the components, the components whose wavelet coefficient power spectrum entropy exceeds the threshold, the remaining useful IMF components are summed with the Res components to obtain the denoised time series InSAR surface deformation information: I update (n)=[IMF m (n)] update,i +[Res M (n)] update Among them, I update (n) is the denoised time series InSAR phase simulation data set; [IMF m (n)] update,i To select the IMF components for reconstruction, n is the number of images, corresponding to the time step.

Images