A distributed scatterer InSAR phase enhancement method based on Laplacian operator

Abstract

The application discloses a distributed scatterer InSAR phase enhancement method based on a Laplacian operator, and comprises the following steps: constructing a SAR amplitude information set, identifying a homogeneous pixel set based on the SAR amplitude information set; identifying an approximate deformation range based on the Laplacian operator, and adaptively adjusting the SHP window size of a deformation area 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; combining the sample complex coherence matrix model, performing a split-pile compression phase enhancement calibration processing on a time sequence SAR image set, and obtaining an enhanced phase after calibration. The application can effectively improve the problem that noise suppression and fringe protection are difficult to be considered in the deformation interferogram after the DSInSAR phase enhancement under a complex environment, and improve the quality of the enhanced phase.

Description

Technical Field

[0001] This invention belongs to the field of synthetic aperture radar interferometric data processing technology, and particularly relates to a distributed scatterer InSAR phase enhancement method based on the Laplacian operator. Background Technology

[0002] Interferometric Synthetic Aperture Radar (InSAR) technology has strong application potential in the field of Earth deformation observation. It obtains high-precision surface deformation products by identifying potential high-coherence target points across a wide area of ​​the Earth's surface, performing statistical analysis, and iteratively separating atmospheric delay information. Distributed Scatterer (DS) InSAR technology is widely used due to its ability to effectively improve the spatial coverage of monitoring points and provide richer and more reliable deformation detail information in non-artificial surface areas. However, the random scattering signals within DS cells make their internal phase information extremely unstable. Therefore, phase enhancement methods aimed at improving the phase signal-to-noise ratio and recovering the optimal phase of the DS cells are particularly important. The PTA algorithm was first proposed within the SqueeSAR framework, outputting results through nonlinear calculation of the maximum likelihood estimation model. Subsequently, many extended phase enhancement methods emerged, such as least-squares phase enhancement estimation and eigenvalue decomposition output. The later-developed EMI algorithm can simultaneously maintain the accuracy of maximum likelihood estimation and the efficiency of eigenvalue decomposition. Other methods focus on frequency estimation and principal phase spatial filtering, providing solutions for phase enhancement in dense fringes.

[0003] In evaluating the quality of phase enhancement, the degree of interferogram noise suppression and the integrity of the interferometric fringe signal are two key criteria that require attention. While the aforementioned phase enhancement methods have been widely used in InSAR data processing, the trade-off between phase noise suppression and fringe preservation is rarely considered. In fact, most traditional methods are based on a complex coherence matrix constructed using SHPs within a fixed window as the solution basis. The accurate estimation of this matrix largely determines the quality of the phase output. However, fixed-window solutions are prone to two problems: in non-deformed regions and some deformed regions, the insufficient number of SHPs within the window limits the improvement of the interferogram's signal-to-noise ratio, and the reliability of the complex coherence matrix is ​​questioned; in deformed regions, an excessively large window and a large number of integrated SHPs lead to severe multi-view blurring, resulting in fringe loss. Therefore, how to flexibly and effectively balance the relationship between the degree of interferogram phase noise suppression and fringe signal preservation is particularly important, especially in obtaining high-precision phase enhancement results in complex surface deformation regions, which has become an important research task. Summary of the Invention

[0004] To address the aforementioned technical problems, this invention proposes a distributed scatterer InSAR phase enhancement method based on the Laplacian operator, which can effectively improve the problem of simultaneously achieving noise suppression and fringe protection in the deformed interferogram after DSInSAR phase enhancement under complex environments, thereby improving the quality of the enhanced phase.

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

[0006] Construct a SAR amplitude information set, and identify homogeneous pixel sets based on the SAR amplitude information set;

[0007] The approximate deformation range is identified based on the Laplacian operator, and the size of the SHP window of the deformation region is adaptively adjusted based on the identified approximate deformation range.

[0008] Non-local similarity fusion is performed on the SHP within the adaptive window to construct a sample complex coherence matrix model;

[0009] Based on the aforementioned sample complex coherence matrix model, the time-series SAR image set is subjected to stacked compression phase enhancement calibration processing to obtain the calibrated enhanced phase.

[0010] Technical advantages of this invention: This invention discloses a distributed scatterer InSAR phase enhancement method based on the Laplacian operator. It identifies homogeneous pixels using a confidence interval hypothesis testing method; delineates the approximate deformation range of the study area based on the Laplacian operator; and automatically adjusts the window size of the SHP (Synchronous Hierarchical Positive Spectroradiometer) based on the deformation level. Simultaneously, this invention utilizes nonlocal theory to fully fuse and utilize similar SHP data in space. The flexible SHP window allocation mode of this invention can establish a more accurate complex coherence matrix based on the actual environment monitored in the study area. Integrating similar SHP in space can effectively mitigate the problem of poor phase reliability caused by insufficient samples in low-coherence regions, further improving the interferometric signal-to-noise ratio. Furthermore, the adaptive window mode allows the deformation region to output a more reasonable window size based on its actual deformation level, avoiding excessive loss of phase fringes and reducing the difficulty of deformation phase interpretation. In addition, this invention also uses a strategy of stacked compressed SAR datasets for phase enhancement output. By estimating each stack of data separately and then unifying the benchmark, the accuracy of phase estimation can be further improved, and the estimation efficiency can be maximized. The method proposed in this invention can effectively protect the detailed information of interference fringes while suppressing noise with high efficiency and high precision. Attached Figure Description

[0011] The accompanying drawings, which form part of this application, are used to provide a further understanding of this application. The illustrative embodiments and descriptions of this application are used to explain this application and do not constitute an undue limitation of this application. In the drawings:

[0012] Figure 1 This is a schematic flowchart of a distributed scatterer InSAR phase enhancement method based on the Laplacian operator according to an embodiment of the present invention.

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

[0014] Figure 3 The diagram shows the partial real surface deformation phase data used in the embodiments of the present invention, including EMI, Sequential, SP-Sequential (only SHP is fused), phase enhancement results of the method of the present invention, and 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 locally magnified phase enhancement result of SP-Sequential, and (f) is the locally magnified phase enhancement result of the method of the present invention.

[0015] Figure 4 This diagram illustrates the distribution of the EMI, Sequential, and SP-Sequential phase data of real surface deformation used in the embodiments of the present invention, as well as the distribution of the posterior best fit of the phase enhancement method of the present invention and the corresponding number of DS candidate points. In this diagram, (a) shows the best fit of EMI and the statistical results of DS candidate points, (b) shows the best fit of Sequential and the statistical results of DS candidate points, (c) shows the best fit of SP-Sequential and the statistical results of DS candidate points, and (d) shows the best fit of the method of the present invention and the statistical results of DS candidate points. Detailed Implementation

[0016] It should be noted that, unless otherwise specified, the embodiments and features described in this application can be combined with each other. This application will now be described in detail with reference to the accompanying drawings and embodiments.

[0017] It should be noted that the steps shown in the flowchart in the accompanying drawings can be executed in a computer system such as a set of computer-executable instructions, and although a logical order is shown in the flowchart, in some cases the steps shown or described may be executed in a different order than 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] Construct a SAR amplitude information set, and identify homogeneous pixel sets based on the SAR amplitude information set;

[0020] The approximate deformation range is identified based on the Laplacian operator, and the size of the SHP window of the deformation region is adaptively adjusted based on the identified approximate deformation range.

[0021] Non-local similarity fusion is performed on the SHP within the adaptive window to construct a sample complex coherence matrix model;

[0022] Based on the aforementioned sample complex coherence matrix model, the time-series SAR image set is subjected to stacked 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 confidence interval for pixel samples is constructed to determine the distribution differences among data in the SAR amplitude information set and to identify homogeneous pixel sets.

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

[0027]

[0028] Where, N SLC Let α be the number of SAR images and α be the interval confidence level. For reference pixel time average, The time average of the pixels to be estimated. Let α be the quantile of the standard gamma distribution, and p represent the probability.

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

[0030] Based on the Laplacian operator traversing long-period interference information, the original interferogram with the largest deformation coverage is selected, and a deformation level threshold is set to identify the approximate deformation range.

[0031] Specifically, the Laplacian operator from the field of image processing is introduced to delineate the deformation range of the interference phase. By synthesizing long-period interference information, the original interferogram with the largest deformation coverage is selected. By setting a deformation level threshold, the approximate range of the deformation region is located. Furthermore, an automatically adjusted SHP window processing method is used, replacing the original fixed 15×15 window with an adaptive window to distinguish it from the non-deformed region, effectively reducing the multi-view blurring of the phase fringes in the deformation region and flexibly protecting the integrity of the interference fringes.

[0032] Based on the Laplacian operator, an approximate deformation range is identified, a deformation threshold is set, and the SHP window size of the deformation region is reallocated. The spatial second-order Laplacian operator can be expressed as:

[0033]

[0034] in, The gradient represents the scalar field u. Let represent the second derivative of u in the x-direction; based on this, the Laplacian gradient within the homogeneous pixel window is:

[0035]

[0036] in, For the gradient of change in the window, Represents the gradient in the x-direction. The Laplacian gradient represents the gradient along the y-direction; its discretized expansion is shown in the following equation:

[0037]

[0038] Where θ is the phase value within the window and r is the window radius (maximum 7, minimum 3). By setting the deformation level judgment threshold, the approximate deformation range can be delineated based on the temporal interferogram of the study area.

[0039] Furthermore, the adaptive adjustment of 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 a gradient threshold, compare the current gradient with the threshold, and if the gradient exceeds the threshold, narrow the window. Iterate until the gradient is less than the threshold or the minimum window size is reached, then output the window to complete the adaptive adjustment of the SHP window size in the deformation region.

[0041] Specifically, starting with a 15×15 window, the Laplacian gradient within the window is calculated. Π / 2 is set as the gradient threshold. The current gradient is compared with the threshold. If the gradient exceeds the threshold, the window is narrowed. This process is repeated until the gradient is less than the threshold or the minimum window size (7×7) is reached, at which point the window is output.

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

[0043]

[0044] Where step1 represents the data input step, and r is the initial window radius. π / 2 is the initial window's Laplacian gradient, and π / 2 is the discrimination threshold. Step 2 represents the SHP window radius iterative calculation step. First, it determines whether the initial window's Laplacian gradient exceeds the threshold. If it does not exceed the threshold, it is determined to be a non-deformation region, and the initial window size is retained. If it exceeds the threshold, the window radius is narrowed, and the new window gradient is recalculated and compared with the threshold again until the window gradient is less than the threshold or reaches the minimum window radius of 3. Then, the window is output. Break means to exit the loop process. Step 3 represents the window radius output step.

[0045] Furthermore, nonlocal similarity fusion is performed on the SHP within the adaptive window to construct the sample complex coherence matrix model, which includes: performing similarity fusion on the SHP within the aforementioned adaptive window based on nonlocal theory.

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

[0047]

[0048] In the formula, 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 transpose operation.

[0049] Furthermore, the original complex coherence matrix of the neighborhood of each pixel in the homogeneous pixel set is constructed respectively. Similarity judgment is performed between the PDF of the homogeneous pixel set and the original pixels. Weight coefficients are assigned to the output similarity results. Non-local similarity fusion is performed on the SHP within the adaptive window. The sample complex coherence matrix model is constructed as follows:

[0050]

[0051] in, For the newly created sample complex coherence matrix model, wt is the similarity weight coefficient of homogeneous pixels, and x refHere, x represents the original pixel, l represents the number of homogeneous pixel samples, i represents the homogeneous pixel index, and x represents the original pixel. i For the i-th pixel in a homogeneous sample set, Let be the original complex coherence matrix constructed for the neighborhood of the i-th pixel in a homogeneous sample set. The original complex coherence matrix is ​​constructed for the neighborhood of each pixel.

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

[0053] The time-series SAR image set is split into clusters, and the clustered SAR image subsets are input into the sample complex coherence matrix model for phase enhancement. The phase-enhanced SAR subsets are then compressed to obtain compressed image sets of each cluster.

[0054] Phase estimation is performed again on the compressed image set to obtain estimation results. The estimation results are then added to the phase enhancement of each sub-stack to obtain the calibrated phase enhancement.

[0055] Furthermore, inputting the subset of SAR images after stacking into the sample complex coherence matrix model for phase enhancement includes:

[0056]

[0057] in, Let ξ be the output sub-heap enhancement phase vector, and ξ be the true SLC phase vector. As a coherence estimator, This is the Hadamard product operation, where H is the transpose operation. The original SAR image set is split, and phase enhancement is performed on the output portions of a smaller number of neighboring SAR images to reduce the impact of decoherence errors and decrease the phase enhancement dimensionality.

[0058] Furthermore, the compressed image sets obtained after compressing each sub-heap include:

[0059]

[0060] in, To compress the image, ω n Let B be the number of SAR images in the sub-pile, and Z(x,s) be the temporal SLC complex data vector of the s-th pixel.

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

[0062]

[0063] in, This is the enhancement phase vector for the k-th sub-heap. For the calibration phase of the k-th sub-heap, This is the enhanced phase after calibration.

[0064] To verify the technical effect of this invention, deformation interferometry data from a mining area were selected, and phase enhancement processing analysis was performed under the same conditions based on EMI, Sequential, SP-Sequential (SHP fusion only), and the method of this invention. The schematic diagrams of the simulated deformation interferometry phase and the window size after window redistribution are shown below. Figure 2 As shown; Figure 2 (a) is the simulated surface deformation interferometric phase diagram. Figure 2 (b) represents the adaptive window size identified based on the simulated deformation data; after adjusting the window according to the actual deformation, Figure 2 The window size corresponding to the approximately deformed region begins to narrow, while the non-deformed region remains basically unchanged, making the window sizes more reasonable in relation to the SHP output. The results of the real deformed interferogram after phase enhancement by EMI, Sequential, SP-Sequential, and the method of this invention are as follows: Figure 3 As shown; Figure 3 (a) shows the phase enhancement result of EMI. Figure 3 (b) shows the phase enhancement results for Sequential. Figure 3 (c) shows the phase enhancement results of SP-Sequential. Figure 3 (d) represents the phase enhancement result of the method of the present invention. Figure 3 (e) represents the locally magnified phase enhancement result of SP-Sequential. Figure 3 (f) represents the local amplified phase enhancement result of the method of the present invention; for EMI, the number of residual points (×10) 5 The value is 1.22, and the PSD (×10) 6 The signal-to-noise ratio (SNR) of the method in this invention 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 in this invention is 0.29 and the PSD is 2.47. It can be found that the enhanced interferograms of EMI and Sequential still have a lot of noise points, and the fringe signal in the near-deformation region is still relatively blurred; SP-Sequential has the strongest SNR improvement capability, but the deformation fringes show a more serious aliasing loss phenomenon; in contrast, the noise filtering capability of the method in this invention is far better than EMI and Sequential, second only to SP-Sequential, but the overall phase of the deformation region edge is smoother and the fringe edge continuity is stronger. The posterior best fit distribution of each phase enhancement method in the study area and the corresponding number of DS candidate points are as follows. Figure 4 As shown; Figure 4 (a) shows the best fit of EMI and the statistical results of DS candidate points. Figure 4 (b) shows the best fit of Sequential and the statistical results of DS candidate points. Figure 4 (c) represents the best fit of SP-Sequential and the statistical results of DS candidate points. Figure 4 (d) represents the optimal fit and statistical results of DS candidate points for the method of this invention. Among them, Sequential showed the worst optimal fit, with fewer candidate DS points (183505) than EMI (224285). However, the fitting results of SP-Sequential and the method of this invention were relatively enhanced, with DS candidate points reaching 273860 and 258668 respectively. In particular, the method of this invention showed high-quality enhancement in non-deformation regions while also improving the overall performance in deformation regions (shown by the black box in the figure below). In summary, the method of this invention can adapt to surface deformation interferometric phase enhancement under complex environments and has good noise suppression and fringe protection capabilities.

[0065] In summary, this 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 method of this invention achieves better optimization accuracy. Compared with SP-Sequential, which only performs SHP fusion, this invention has stronger protection capability for surface deformation interferometric fringe information. First, the SHP neighborhood set of pixel samples is identified based on the parameter hypothesis testing method; second, the approximate deformation range is identified according to the Laplacian operator, a deformation threshold is set, and the SHP window size of the deformation region is reallocated; third, the SHP within the above adaptive window is fused based on nonlocal theory; finally, the SAR dataset is compressed and the phase enhancement output is performed. Through simulated 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 this invention can effectively improve the problem of simultaneously achieving high noise removal and phase fringe protection in DSInSAR phase enhancement.

[0066] The above are merely preferred embodiments of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A distributed scatterer InSAR phase enhancement method based on the Laplacian operator, characterized in that, include: Construct a SAR amplitude information set, and identify homogeneous pixel sets based on the SAR amplitude information set; The approximate deformation range is identified based on the Laplacian operator, and the size of the SHP window of the deformation region is adaptively adjusted based on the identified approximate deformation range. Non-local similarity fusion is performed on the SHP within the adaptive window to construct a sample complex coherence matrix model; Based on the aforementioned sample complex coherence matrix model, the time-series SAR image set is subjected to stacked compression phase enhancement calibration processing to obtain the calibrated enhanced phase; 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 confidence interval for pixel samples is constructed to determine the distribution differences among data in the SAR amplitude information set and to identify homogeneous pixel sets. Constructing confidence intervals for pixel samples includes: ， in, For SAR image numbers, For interval confidence levels, For reference pixel time average, The time average of the pixels to be estimated. For standard gamma distribution quantiles, Represents probability; The approximate deformation range identified based on the Laplacian operator includes: Based on the Laplacian operator traversing long-period interference information, the original interferogram with the largest deformation coverage is selected, and a deformation level threshold is set to identify the approximate deformation range. 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 a gradient threshold, compare the current gradient with the threshold, and if the gradient exceeds the threshold, narrow the window. Iterate until the gradient is less than the threshold or the minimum window size is reached, then output the window to complete the adaptive adjustment of the SHP window size in the deformation region.

2. The distributed scatterer InSAR phase enhancement method based on the Laplacian operator as described in claim 1, characterized in that, Nonlocal similarity fusion is performed on the SHP within the adaptive window to construct the sample complex coherence matrix model, including: ， in, For the newly created sample complex coherence matrix model, The similarity weighting coefficient for homogeneous pixels. Original pixels, l The number of homogeneous pixel samples. i The serial number of the homogeneous pixel. The first in a homogeneous sample set i 1 pixel, The first in a homogeneous sample set i The original complex coherence matrix constructed from the neighborhood of pixels. The original complex coherence matrix is ​​constructed for the neighborhood of each pixel.

3. The distributed scatterer InSAR phase enhancement method based on the Laplacian operator as described in claim 1, characterized in that, The complete enhanced phase obtained after calibration includes: The time-series SAR image set is split into clusters, and the clustered SAR image subsets are input into the sample complex coherence matrix model for phase enhancement. The phase-enhanced SAR subsets are then compressed to obtain compressed image sets of each cluster. Phase estimation is performed again on the compressed image set to obtain estimation results. The estimation results are then added to the phase enhancement of each sub-stack to obtain the calibrated phase enhancement.

4. The distributed scatterer InSAR phase enhancement method based on the Laplacian operator as described in claim 3, characterized in that, The phase enhancement process involves inputting a subset of the split SAR images into the sample complex coherence matrix model, including: ， in, Enhance the phase vector of the output sub-heap. For the true SLC phase vector, As a coherence estimator, For Hadamard product operation, This is a transpose operation.

5. The distributed scatterer InSAR phase enhancement method based on the Laplacian operator as described in claim 3, characterized in that, The compressed image set obtained after each sub-heap compression includes: ， in, To compress the image, For compressible variation basis, B This represents the number of SAR images in the sub-stack. For the first s A temporal SLC complex data vector of pixels.

6. The distributed scatterer InSAR phase enhancement method based on the Laplacian operator as described in claim 3, characterized in that, The enhanced phase obtained after calibration includes: ， in, For the first k The enhanced phase vector of each sub-heap For the first k The calibration phase of each sub-pillar This is the enhanced phase after calibration.

Images