Distributed scatterer InSAR phase enhancement method based on Laplacian operator

Abstract

The invention discloses a distributed scatterer InSAR phase enhancement method based on a Laplacian operator, and the method comprises the steps: constructing an SAR amplitude information set, and recognizing a homogeneous pixel set based on the SAR amplitude information set; identifying an approximate deformation range based on a Laplacian operator, and adaptively adjusting the window size of a deformation area SHP based on the identified approximate deformation range; performing non-local similarity fusion on the SHP in the adaptive window, and constructing a sample complex coherence matrix model; and in combination with the sample complex coherence matrix model, performing stack compression phase enhancement calibration processing on the time sequence SAR image set to obtain a calibrated enhanced phase. According to the method, the problem that noise suppression and fringe protection are difficult to take into account in a deformation interferogram after DSInSAR phase enhancement in a complex environment can be effectively improved, and the phase enhancement quality is improved.

Description

Technical Field

[0001] The invention belongs to the technical field of synthetic aperture radar interference data processing, and in particular relates to a distributed scatterer InSAR phase enhancement method based on a Laplacian operator. Background Art

[0002] Time series synthetic aperture radar interferometry (InSAR) technology has great application potential in the field of ground deformation observation. It identifies potential high-coherence target points in a wide area of ​​the surface for statistical analysis and iteratively separates atmospheric delay information to obtain high-precision surface deformation products. Among them, distributed scatterer InSAR technology is widely used because it can effectively improve the spatial coverage of monitoring points in non-artificial surface areas and provide more abundant and reliable deformation detail information. However, the random scattering signal in the DS unit makes its internal phase information extremely unstable. Therefore, the phase enhancement method aimed at improving the phase signal-to-noise ratio and restoring the optimal phase of DS is particularly important. The PTA algorithm was first proposed under the SqueeSAR framework. The output result of the maximum likelihood estimation model was obtained by nonlinear solution. Subsequently, many extended phase enhancement methods emerged, such as least squares phase enhancement estimation, eigenvalue decomposition output, etc. The EMI algorithm developed later can simultaneously maintain the accuracy of maximum likelihood estimation and the efficiency of eigenvalue decomposition. There are also methods that focus on frequency estimation and main phase space filtering, which also provides a solution for phase enhancement at dense stripes.

[0003] In terms of phase enhancement quality evaluation, the degree of interference pattern noise suppression and the integrity of interference fringe signals are two standards that need to be focused on. The above phase enhancement methods have been widely used in InSAR data processing, but the contradiction between phase noise suppression and phase fringe protection is rarely considered. In fact, most traditional methods are based on the complex coherence matrix constructed by SHP in a fixed window as the solution basis. The accurate estimation of the matrix largely determines the quality of the phase output. The solution based on the fixed window is prone to the following two problems: the interference pattern signal-to-noise ratio is limited in the non-deformed area and some deformed areas due to the small number of SHP in the window, and the reliability of the complex coherence matrix is ​​questioned; the deformation area has a serious multi-view blur phenomenon due to the large window and the large number of integrated SHPs, and the interference fringes are lost. Therefore, how to flexibly and effectively balance the relationship between the degree of interference pattern phase noise suppression and the protection of interference fringe signals is particularly important, especially in complex surface deformation areas. Obtaining high-precision phase enhancement results has become an important research work. Summary of the invention

[0004] In order to solve the above technical problems, the present invention proposes a distributed scatterer InSAR phase enhancement method based on the Laplacian operator, which can effectively improve the problem of difficulty in balancing noise suppression and fringe protection in the deformation interferogram after DSInSAR phase enhancement in complex environments, and improve the quality of the enhanced phase.

[0005] To achieve the above object, the present invention provides a distributed scatterer InSAR phase enhancement method based on Laplacian operator, comprising:

[0006] Constructing a SAR amplitude information set, and identifying a set of homogeneous pixels based on the SAR amplitude information set;

[0007] Based on the Laplacian operator, the approximate deformation range is identified, and the SHP window size of the deformation area is adaptively adjusted based on the identified approximate deformation range;

[0008] Perform non-local similarity fusion on the SHP within the adaptive window and construct a sample complex coherence matrix model;

[0009] Combined with the sample complex coherence matrix model, the time series SAR image set is subjected to stacking compression phase enhancement calibration processing to obtain the calibrated enhanced phase.

[0010] Technical effect of the present invention: The present invention discloses a distributed scatterer InSAR phase enhancement method based on Laplacian operator, which identifies homogeneous pixels through a confidence interval hypothesis test method; delineates the approximate deformation range of the study area based on the Laplacian operator, and automatically adjusts the window size of the SHP according to the deformation magnitude; at the same time, the present invention makes full use of the integration of similar SHP data in space by using non-local theory. The flexible SHP window allocation mode of the present invention can establish a more accurate complex coherence matrix according to the actual environment monitored in the study area. Similar SHPs in the integrated space can effectively reduce the problem of poor phase reliability caused by too few samples in the low coherence area, and further improve the interference signal-to-noise ratio; and the adaptive window mode enables the deformation area to output a more reasonable window size according to its actual deformation magnitude, avoiding excessive loss of phase fringes and reducing the difficulty of deformation phase interpretation. In addition, the present invention also performs phase enhancement output based on the strategy of stacking and compressing SAR data sets, and uses each stack of data to estimate separately and then unify the benchmark, which can further improve the accuracy of phase estimation and maximize the estimation efficiency. The method proposed in the present invention can effectively protect the detailed information of the interference fringes while suppressing noise with high efficiency and high precision. BRIEF DESCRIPTION OF THE DRAWINGS

[0011] The drawings constituting a part of the present application are used to provide a further understanding of the present application. The illustrative embodiments and descriptions of the present application are used to explain the present application and do not constitute an improper limitation on the present application. In the drawings:

[0012] Figure 1 A schematic diagram of a process of a distributed scatterer InSAR phase enhancement method based on a Laplacian operator according to an embodiment of the present invention;

[0013] Figure 2 Schematic diagram of simulated surface deformation interference diagram and corresponding adaptive window size distribution according to an embodiment of the present invention, wherein (a) is a simulated surface deformation interference phase diagram, and (b) is an adaptive window size diagram identified based on the simulated deformation data;

[0014] Figure 3 Schematic diagram of some real surface deformation phase data EMI, Sequential, SP-Sequential (only SHP is fused) used in the embodiments of the present invention, the phase enhancement result of the method of the present invention, and the local phase enhancement results of SP-Sequential and the present invention, (a) is the phase enhancement result of EMI, (b) is the phase enhancement result of Sequential, (c) is the phase enhancement result of SP-Sequential, (d) is the phase enhancement result of the method of the present invention, (e) is the local enlarged phase enhancement result of SP-Sequential, and (f) is the local enlarged phase enhancement result of the method of the present invention;

[0015] Figure 4 Schematic diagram of the best fit distribution of the real surface deformation phase data EMI, Sequential, SP-Sequential and the phase enhancement method of the present invention and the corresponding DS candidate point statistics, wherein (a) is the best fit and DS candidate point statistics of EMI, (b) is the best fit and DS candidate point statistics of Sequential, (c) is the best fit and DS candidate point statistics of SP-Sequential, and (d) is the best fit and DS candidate point statistics of the method of the present invention. DETAILED DESCRIPTION

[0016] It should be noted that, in the absence of conflict, the embodiments and features in the embodiments of the present application can be combined with each other. The present application will be described in detail below with reference to the accompanying drawings and in combination with the embodiments.

[0017] It should be noted that the steps shown in the flowcharts of the accompanying drawings can be executed in a computer system such as a set of computer executable instructions, and that, although a logical order is shown in the flowcharts, in some cases, the steps shown or described can be executed in an order different from that shown here.

[0018] like Figure 1 As shown, this embodiment provides a distributed scatterer InSAR phase enhancement method based on the Laplacian operator, including:

[0019] Constructing a SAR amplitude information set, and identifying a set of homogeneous pixels based on the SAR amplitude information set;

[0020] Based on the Laplacian operator, the approximate deformation range is identified, and the SHP window size of the deformation area is adaptively adjusted based on the identified approximate deformation range;

[0021] Perform non-local similarity fusion on the SHP within the adaptive window and construct a sample complex coherence matrix model;

[0022] Combined with the sample complex coherence matrix model, the time series SAR image set is subjected to stacking compression phase enhancement calibration processing to obtain the calibrated enhanced phase.

[0023] Furthermore, constructing a SAR amplitude information set, and identifying homogeneous pixels based on the SAR amplitude information set includes:

[0024] Acquire a time series SAR image set, preprocess the time series SAR image set, and construct the SAR amplitude information set;

[0025] By assuming a known sample probability density function, a pixel sample confidence interval is constructed, the distribution difference between the data in the SAR amplitude information set is determined, and a homogeneous pixel set is identified.

[0026] Furthermore, constructing the pixel sample confidence interval includes:

[0027]

[0028] Among them, N SLC is the number of SAR images, α is the interval confidence level, is the temporal mean of the reference pixel, is the time mean of the pixel to be estimated, is the alpha quantile of the standard gamma distribution, and p represents the probability.

[0029] Furthermore, the approximate deformation range identified based on the Laplacian operator includes:

[0030] According to the interference information of long period traversed by Laplacian operator, the original interference pattern with the largest deformation coverage is selected, the deformation magnitude threshold is set, and the approximate deformation range is identified.

[0031] Specifically, the Laplacian operator in the field of image processing is introduced to delineate the deformation range of the interference phase. By integrating the long-period interference information, the original interference pattern with the largest deformation coverage is selected, and the approximate range of the deformation area is located by setting the deformation level threshold. The automatically adjusted SHP window processing method is used to replace the original fixed 15×15 window with an adaptive window to distinguish it from the non-deformed area, effectively weakening the multi-view blur of the phase fringes in the deformation area and flexibly protecting the integrity of the interference fringes.

[0032] According to the Laplacian operator, the approximate deformation range is identified, the deformation threshold is set, and the SHP window size of the deformation area is redistributed; the spatial second-order Laplacian operator can be expressed as:

[0033]

[0034] in, represents the gradient of the scalar field u, represents the second-order derivative of u in the x direction; on this basis, the Laplacian gradient in the homogeneous pixel window is:

[0035]

[0036] in, is the gradient of change in the window, represents the gradient in the x direction, Represents the gradient in the y direction; the Laplacian gradient is discretized and expanded as follows:

[0037]

[0038] Among them, θ is the phase value in the window, r is the window radius (the maximum is 7 and the minimum is 3). By setting the deformation magnitude judgment threshold, the approximate deformation range can be delineated based on the time-series interference pattern of the study area.

[0039] Further, adaptively adjusting the SHP window size based on the identified approximate deformation range includes:

[0040] Determine the initial window size, calculate the Laplacian gradient within the window, set the gradient threshold, determine the size of the current gradient and the threshold, and narrow the window if it exceeds the threshold. Repeat the process until the gradient is less than the threshold or reaches the minimum window, then output the window and complete the adaptive adjustment of the SHP window size in the deformation area.

[0041] Specifically, a 15×15 window is used as the starting value, the Laplacian gradient in the window is calculated, Π / 2 is set as the gradient threshold, the current gradient and the threshold are judged, and the window is narrowed when it exceeds the threshold. The cycle is iterated until the gradient is less than the threshold or reaches the minimum window (7×7), and the window is output.

[0042] The specific processing flow of the automatic adjustment of the deformation area SHP window is as follows:

[0043]

[0044] Among them, step1 represents the data input step, r is the initial window radius, is the Laplacian gradient of the initial window, π / 2 is the discrimination threshold; step2 represents the iterative calculation step of the SHP window radius. First, it is determined whether the Laplacian gradient of the initial window exceeds the threshold. If it does not exceed, it is judged as a non-deformation area and the initial window size is retained. If it exceeds, the window radius is narrowed and the new window gradient is recalculated to continue to compare with the threshold until the window gradient is less than the threshold or reaches the minimum window radius 3. The window is output, and break represents jumping out of the loop process; step3 represents the window radius output step.

[0045] Furthermore, performing non-local similarity fusion on the SHP within the adaptive window and constructing a sample complex coherence matrix model includes: performing similarity fusion on the SHP within the adaptive window according to the non-local theory.

[0046] Furthermore, combined with the central limit theorem, the complex Gaussian distribution PDF of any pixel in a homogeneous pixel set can be expressed as:

[0047]

[0048] Where x is a single pixel observation, Λ represents the original complex matrix built based on the observed pixels, ν represents the mean of the Gaussian distribution, and T represents the transposition operation.

[0049] Furthermore, the original complex coherence matrix of each pixel neighborhood in the homogeneous pixel set is constructed respectively, similarity judgment is performed based on the PDF of the homogeneous pixel set and the original pixel, a weight coefficient is assigned to the output similarity result, and non-local similarity fusion is performed on the SHP within the adaptive window. The construction of the sample complex coherence matrix model includes:

[0050]

[0051] in, is the newly built sample complex coherence matrix model, wt is the weight coefficient of homogeneous pixel similarity, x refis the original pixel, l is the number of homogeneous pixel samples, i is the homogeneous pixel sequence number, x i is the i-th pixel in the homogeneous sample set, is the original complex coherence matrix constructed for the ith pixel neighborhood in the homogeneous sample set, The raw complex coherence matrix constructed for each pixel neighborhood.

[0052] Furthermore, obtaining the complete enhanced phase after calibration includes:

[0053] The time series SAR image set is stacked, the stacked SAR image subsets are input into the sample complex coherence matrix model for phase enhancement, and the phase-enhanced SAR subsets are compressed to obtain compressed image sets of each sub-stack;

[0054] The phase of the compressed image set is estimated again to obtain an estimation result, and the estimation result is added to the enhanced phase of each sub-stack to obtain a calibrated enhanced phase.

[0055] Furthermore, inputting the stacked SAR image subsets into the sample complex coherence matrix model for phase enhancement includes:

[0056]

[0057] in, is the output sub-stack enhanced phase vector, ξ is the real SLC phase vector, is the coherence estimator, is the Hadamard product operation, and H is the transpose operation. The original SAR image set is split, and the phase is enhanced based on the output of a small number of adjacent SAR images to weaken the influence of the decoherence error and reduce the phase enhancement dimension.

[0058] Furthermore, obtaining a compressed image set after compression of each sub-heap includes:

[0059]

[0060] in, is the compressed image, ω n is the compressed variation basis, B is the number of SAR images in the sub-pile, and Z(x,s) is the time-series SLC complex data vector of the s-th pixel.

[0061] Further, obtaining the calibrated enhanced phase includes:

[0062]

[0063] in, is the enhanced phase vector of the kth sub-stack, is the calibration phase of the kth sub-stack, is the enhanced phase after calibration.

[0064] In order to verify the technical effect of the present invention, the deformation interference data of the mining area was selected and phase enhancement processing analysis was carried out under the same conditions based on EMI, Sequential, SP-Sequential (only SHP was fused) and the method of the present invention. Among them, the schematic diagram of the simulated deformation interference phase and the window size after window redistribution is shown in the figure. Figure 2 As shown; Figure 2 (a) is the simulated surface deformation interference phase map. Figure 2 (b) is the adaptive window size identified based on the simulated deformation data; after adjusting the window in combination with the actual deformation situation, Figure 2 The window size corresponding to the approximate deformation area begins to narrow, while the non-deformation area remains basically unchanged, and the window size and SHP output are more reasonable. The results of the real deformation interferogram after phase enhancement by EMI, Sequential, SP-Sequential and the method of the present invention are as follows: Figure 3 As shown; Figure 3 (a) is the phase enhancement result of EMI. Figure 3 (b) is the phase enhancement result of Sequential. Figure 3 (c) is the phase enhancement result of SP-Sequential. Figure 3 (d) is the phase enhancement result of the method of the present invention, Figure 3 (e) is the local amplification phase enhancement result of SP-Sequential. Figure 3 (f) is the local amplification phase enhancement result of the method of the present invention; for EMI, the number of residual points (×10 5 ) is 1.22, PSD(×10 6 ) is 4.8; for Sequential, the number of residual points is 0.43, and the PSD is 3.11; for SP-Sequential, the number of residual points is 0.23, and the PSD is 2.41; the number of residual points of the method of the present invention is 0.29, and the PSD is 2.47; it can be found that there are still many noise points remaining in the enhanced interference diagrams of EMI and Sequential, and the stripe signals in the approximate deformation area are still relatively fuzzy; SP-Sequential has the strongest signal-to-noise ratio improvement ability, but the deformation stripes have more serious aliasing loss phenomenon; in comparison, the noise filtering ability of the method of the present invention is far better than EMI and Sequential, second only to SP-Sequential, but the overall phase of the edge of the deformation area is smoother, and the continuity of the stripe edge is stronger. The posterior optimal fit distribution of each phase enhancement method in the study area and the corresponding number of DS candidate points are shown in the figure. Figure 4 As shown; Figure 4 (a) is the best fit of EMI and the statistical results of DS candidate points. Figure 4 (b) is the best fit of Sequential and the statistical results of DS candidate points. Figure 4 (c) is the best fit of SP-Sequential and the statistical results of DS candidate points. Figure 4 (d) is the best fit and DS candidate point statistics of the method of the present invention. Among them, Sequential has the worst best fit result, and its number of candidate DS points (183505) is less than EMI (224285), but the fitting results of SP-Sequential and the method of the present invention are relatively enhanced, and the DS candidate points reach 273860 and 258668 respectively. In particular, the method of the present invention, while enhancing the high quality of the results in the non-deformed area, also improves the overall performance of the deformed area (shown in the black box in the figure below). In summary, it is shown in many aspects that the method of the present invention can adapt to the interferometric phase enhancement of surface deformation in complex environments, and has good noise suppression and fringe protection capabilities.

[0065] In summary, the present invention proposes a distributed scatterer InSAR phase enhancement method and system based on the Laplacian operator. Compared with commonly used phase optimization methods such as EMI and Sequential, the optimization accuracy of the method of the present invention is better. Compared with SP-Sequential that only performs SHP fusion, the surface deformation interference fringe information protection capability of the present invention is stronger. First, the SHP neighborhood set of pixel samples is identified based on the parameter hypothesis test method; secondly, the approximate deformation range is identified according to the Laplacian operator, the deformation threshold is set, and the SHP window size of the deformation area is reallocated; thirdly, the SHP in the above-mentioned adaptive window is similarly fused according to the non-local theory; finally, the SAR data set is stacked and compressed for phase enhancement output. Through simulation data and real data experiments, the results show that the distributed scatterer InSAR phase enhancement method and system based on the Laplacian operator proposed in the present invention can effectively improve the problem that it is difficult to simultaneously take into account high noise removal and phase fringe protection during DSInSAR phase enhancement.

[0066] The above are only preferred specific implementations of the present application, but the protection scope of the present application is not limited thereto. Any changes or substitutions that can be easily thought of by a person skilled in the art within the technical scope disclosed in the present application should be included in the protection scope of the present application. Therefore, the protection scope of the present application should be based on the protection scope of the claims.

Claims

1. A distributed scatterer InSAR phase enhancement method based on Laplacian operator, characterized in that: include: Constructing a SAR amplitude information set, and identifying a set of homogeneous pixels based on the SAR amplitude information set; Based on the Laplacian operator, the approximate deformation range is identified, and the SHP window size of the deformation area is adaptively adjusted based on the identified approximate deformation range; Perform non-local similarity fusion on the SHP within the adaptive window and construct a sample complex coherence matrix model; Combined with the sample complex coherence matrix model, the time series SAR image set is subjected to stacking compression phase enhancement calibration processing to obtain the calibrated enhanced phase.

2. The distributed scatterer InSAR phase enhancement method based on Laplacian operator according to claim 1, characterized in that: Constructing a SAR amplitude information set, and identifying homogeneous pixels based on the SAR amplitude information set includes: Acquire a time series SAR image set, preprocess the time series SAR image set, and construct the SAR amplitude information set; By assuming a known sample probability density function, a pixel sample confidence interval is constructed, the distribution difference between the data in the SAR amplitude information set is determined, and a homogeneous pixel set is identified.

3. The distributed scatterer InSAR phase enhancement method based on Laplacian operator according to claim 1, characterized in that: Constructing a pixel sample confidence interval involves: Among them, N SLC is the number of SAR images, α is the interval confidence level, is the temporal mean of the reference pixel, is the time mean of the pixel to be estimated, g α / 2；NSLC is the alpha quantile of the standard gamma distribution, and p represents the probability.

4. The distributed scatterer InSAR phase enhancement method based on Laplacian operator according to claim 1, characterized in that: The approximate deformation range identified based on the Laplacian operator includes: According to the interference information of long period traversed by Laplacian operator, the original interference pattern with the largest deformation coverage is selected, the deformation magnitude threshold is set, and the approximate deformation range is identified.

5. The distributed scatterer InSAR phase enhancement method based on Laplacian operator according to claim 1, characterized in that: Adaptive adjustment of the SHP window size based on the identified approximate deformation range includes: Determine the initial window size, calculate the Laplacian gradient within the window, set the gradient threshold, determine the size of the current gradient and the threshold, and narrow the window if it exceeds the threshold. Repeat the process until the gradient is less than the threshold or reaches the minimum window, then output the window and complete the adaptive adjustment of the SHP window size in the deformation area.

6. The distributed scatterer InSAR phase enhancement method based on Laplacian operator according to claim 1, characterized in that: The non-local similarity fusion of SHP within the adaptive window and the construction of the sample complex coherence matrix model include: in, is the newly built sample complex coherence matrix model, wt is the weight coefficient of homogeneous pixel similarity, x ref is the original pixel, l is the number of homogeneous pixel samples, i is the homogeneous pixel sequence number, x i is the i-th pixel in the homogeneous sample set, is the original complex coherence matrix constructed for the ith pixel neighborhood in the homogeneous sample set, The raw complex coherence matrix constructed for each pixel neighborhood.

7. The distributed scatterer InSAR phase enhancement method based on Laplacian operator according to claim 1, characterized in that: Obtaining the complete enhanced phase after calibration includes: The time series SAR image set is stacked, the stacked SAR image subsets are input into the sample complex coherence matrix model for phase enhancement, and the phase-enhanced SAR subsets are compressed to obtain compressed image sets of each sub-stack; The phase of the compressed image set is estimated again to obtain an estimation result, and the estimation result is added to the enhanced phase of each sub-stack to obtain a calibrated enhanced phase.

8. The distributed scatterer InSAR phase enhancement method based on Laplacian operator as claimed in claim 7, characterized in that: Inputting the stacked SAR image subsets into the sample complex coherence matrix model for phase enhancement includes: in, is the output sub-stack enhanced phase vector, ξ is the real SLC phase vector, is the coherence estimator, is the Hadamard product operation, and H is the transpose operation.

9. The distributed scatterer InSAR phase enhancement method based on Laplacian operator according to claim 7, characterized in that: The compressed image set obtained after compression of each sub-heap includes: in, is the compressed image, ω n is the compressed variation basis, B is the number of SAR images in the sub-pile, and Z(x,s) is the time-series SLC complex data vector of the s-th pixel.

10. The distributed scatterer InSAR phase enhancement method based on Laplacian operator according to claim 7, characterized in that: Obtaining the calibrated enhanced phase includes: in, is the enhanced phase vector of the kth sub-stack, is the calibration phase of the kth sub-stack, is the enhanced phase after calibration.

Images