A fast phase estimation method for distributed scatterers based on improved PM

Abstract

The present invention discloses a distributed scatterer phase fast estimation method based on improved PM, which belongs to the field of observation data processing technology, including identifying distributed scatterer DS and marking homogeneous pixel SHP of each DS according to the amplitude information of the SAR image after registration; calculating the interference phase mean according to the SHP of DS; constructing a sample covariance matrix C; using the matrix-vector multiplication of PM to iteratively estimate the main eigenvalues ​​and corresponding eigenvectors of C; constructing an adaptive termination strategy for PM, adaptively terminating iteration according to the quality of DS and obtaining the optimal estimation of DS phase. The present invention adopts the above-mentioned distributed scatterer phase fast estimation method based on improved PM, reduces the number of iterations and improves the convergence speed, and by setting an adaptive termination strategy, ensures that the algorithm quickly converges to the optimal solution for high-quality DS, and terminates the iteration in time for low-quality DS, improves the calculation efficiency, effectively solves the problem of excessively time-consuming calculation, and is simple in calculation, and is particularly suitable for parallel calculation.

Description

Technical Field

[0001] The present invention relates to the technical field of observation data processing, and in particular to a distributed scatterer phase fast estimation method based on improved PM. Background Art

[0002] Time-series InSAR technology has become one of the commonly used geodetic technologies due to its advantages such as high precision, wide coverage, and all-weather. The objects of time-series InSAR analysis are a series of special coherent targets, mainly including permanent scatterers (PS) and distributed scatterers (DS). Typical PS are mainly artificial buildings, so they are suitable for urban areas with dense artificial buildings. In non-urban areas where there are no artificial buildings, PS targets are very rare. DS is a natural scatterer widely distributed in nature, such as bare soil, deserts, and sparse low vegetation areas. Therefore, the combination of PS and DS can greatly increase the density of observation points in non-urban areas. However, the backscattering ability of each basic scatterer in the DS pixel is relatively balanced, and there is no dominant strong scattering unit. DS is greatly affected by spatiotemporal incoherence, and the signal coherence is poor and the signal-to-noise ratio is low. To this end, it is necessary to re-estimate the phase of DS according to some methods. A common method is to construct the covariance matrix of DS, and then use the eigendecomposition method to obtain the eigenvector corresponding to the maximum eigenvalue and use it as the DS estimated phase.

[0003] With the launch of new satellites, SAR earth observation has entered an era of high temporal and spatial resolution, multi-polarization, and wide bandwidth, and has also generated unprecedented amounts of data. Taking the Sentinel-1 SAR satellite as an example, a standard image has nearly 1 billion pixels, and the conventional feature decomposition method is too time-consuming to meet the actual production application. PM has the characteristics of simple calculation and small memory usage, but the calculation efficiency and convergence are restricted by the iteration initial value and iteration termination strategy. Summary of the invention

[0004] The purpose of the present invention is to provide a distributed scatterer phase fast estimation method based on improved PM, taking the interference phase mean on homogeneous pixels SHP as the iterative initial value of PM, avoiding PM divergence, reducing the number of iterations and improving the convergence speed, and by constructing an adaptive termination strategy, ensuring that the PM algorithm can quickly obtain solutions for different distributed scatterers DS, greatly improving the computational efficiency and solving the problem of time-consuming computation.

[0005] To achieve the above object, the present invention provides a distributed scatterer phase fast estimation method based on improved PM, comprising the following steps:

[0006] S1. Using the two-sample hypothesis test method, according to the amplitude information of the registered SAR image, the distributed scatterer DS is identified and the homogeneous pixel SHP of each distributed scatterer DS is marked;

[0007] S2. Calculate the mean interference phase of the distributed scatterer DS using homogeneous pixel SHP; the mean interference phase on the homogeneous pixel SHP is set as the initial estimate of the phase of the distributed scatterer DS , and used as the initial value of PM iteration to avoid PM divergence and reduce the number of iterations, the initial value of distributed scatterer DS phase iteration The calculation method is:

[0008] = (1);

[0009] in, Indicates The interference phase on the interference pair, M Indicates the number of homogeneous pixels SHP;

[0010] S3, using homogeneous pixel SHP to construct the sample covariance matrix C; specifically: for each distributed scatterer DS, use its homogeneous pixel SHP to construct the covariance matrix C, that is:

[0011] (2);

[0012] Among them, the superscript represents the conjugate transpose of the matrix, y is the complex reflectance of each homogeneous pixel SHP in the time series, and Ω represents the set of homogeneous pixel SHPs;

[0013] S4, using PM's matrix-vector multiplication to iteratively estimate the principal eigenvalues ​​and corresponding eigenvectors of the sample covariance matrix C;

[0014] S5, construct an adaptive termination strategy for PM;

[0015] S6, obtain the DS phase of the distributed scatterer, through Iterations, the result obtained at the end of the iteration That is the optimal estimation of the DS phase of the distributed scatterer ,Right now: .

[0016] Preferably, in S1, the two-sample hypothesis testing method is one of Kolmogorov-Smirnov and Anderson-Darling.

[0017] Preferably, S4 specifically includes the following steps:

[0018] S41, the interference phase mean calculated by formula (1) As an initial estimate of the main eigenvector;

[0019] S42, according to the matrix-vector multiplication of PM, iteratively update the main eigenvalue and eigenvector, after k After iterations:

[0020] (3);

[0021] (4);

[0022] in, Indicates The feature vector obtained by the iteration is Indicates k The eigenvalues ​​obtained by iteration;

[0023] S43, Normalized feature vector :

[0024] (5);

[0025] because , when the number of iterations As it approaches infinity, we have:

[0026] (6);

[0027] (7);

[0028] in, is the theoretical principal eigenvalue of the sample covariance matrix C The corresponding feature vector.

[0029] Preferably, in S5, the number of iteration steps is It is impossible to take infinity, and the termination condition must be set for the iteration. The constructed adaptive iteration termination strategy is:

[0030] (8);

[0031] in, and Respectively represent Second and The main eigenvalue obtained by the iteration; Indicates the rate of change of the main eigenvalue. When the rate of change of the main eigenvalue is less than 1e-6, it means that the main eigenvalue and eigenvector have converged. Indicates that the maximum number of iterations is 50.

[0032] Therefore, the present invention adopts the above-mentioned distributed scatterer phase fast estimation method based on improved PM, and takes the interference phase mean on the homogeneous pixel SHP as the iterative initial value of PM, so as to avoid PM divergence, reduce the number of iterations and improve the convergence speed. By setting an adaptive termination strategy, it is ensured that the algorithm can quickly obtain the optimal solution for different distributed scatterers DS. Compared with the existing methods, the computational efficiency of the method of the present invention is greatly improved, which effectively solves the problem of excessively time-consuming calculations. Moreover, because the calculation is simple, it is particularly suitable for parallel calculations, which can further reduce the calculation time.

[0033] The technical solution of the present invention is further described in detail below through the accompanying drawings and embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] Figure 1 It is a flow chart of an embodiment of a distributed scatterer phase fast estimation method based on improved PM of the present invention;

[0035] Figure 2 is a graph showing the relationship between the number of iterations and the DS quality of an embodiment of a distributed scatterer phase fast estimation method based on improved PM of the present invention;

[0036] Figure 3 It is a comparison of the DS phase accuracy estimated by a distributed scatterer fast phase estimation method based on improved PM of the present invention and a conventional eigendecomposition method;

[0037] Figure 4 It is a time-consuming comparison between a distributed scatterer phase fast estimation method based on improved PM of the present invention and a conventional eigendecomposition method. DETAILED DESCRIPTION

[0038] The technical solution of the present invention is further described below through the accompanying drawings and embodiments.

[0039] Unless otherwise defined, the technical terms or scientific terms used in the present invention should be understood by people with ordinary skills in the field to which the present invention belongs. The words "first", "second" and similar words used in the present invention do not indicate any order, quantity or importance, but are only used to distinguish different components. "Include" or "comprise" and similar words mean that the elements or objects appearing before the word include the elements or objects listed after the word and their equivalents, without excluding other elements or objects. "Connect" or "connected" and similar words are not limited to physical or mechanical connections, but may include electrical connections, whether direct or indirect. "Up", "down", "left", "right" and the like are only used to indicate relative positional relationships. When the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0040] Embodiment 1

[0041] The present invention provides a distributed scatterer phase fast estimation method based on improved PM, the process of which is as follows: Figure 1 As shown, the following steps are included:

[0042] S1. Using two-sample hypothesis test methods such as Kolmogorov-Smirnov and Anderson-Darling, the distributed scatterers DS are identified according to the amplitude information of the registered SAR image and the homogeneous pixels SHP of each distributed scatterer DS are marked.

[0043] S2. Use homogeneous pixels SHP to calculate the mean interference phase of the distributed scatterer DS. The essence of the distributed scatterer DS phase estimator is a filter that aims to improve the coherence quality of the distributed scatterer DS phase. Spatial mean filtering is a simple and effective method to improve the quality of interference phase. Its biggest drawback is that the filtering error is large when the scattering types of each pixel in the filter window are different. This defect can be suppressed if mean filtering is performed on homogeneous pixels SHP. Therefore, the mean interference phase on homogeneous pixels SHP is set as the initial estimate of the distributed scatterer DS phase. , and used as the initial value of PM iteration to avoid PM divergence and reduce the number of iterations, the initial value of distributed scatterer DS phase iteration The calculation method is:

[0044] = (1);

[0045] in, Indicates The interference phase on the interference pair, M Indicates the number of homogeneous pixels SHP.

[0046] S3. Use homogeneous pixel SHP to construct the sample covariance matrix C. For each distributed scatterer DS, use its homogeneous pixel SHP to construct the covariance matrix C, that is:

[0047] (2);

[0048] Among them, the superscript represents the conjugate transpose of the matrix, y is the complex reflectance of each pixel on each homogeneous pixel SHP in the time series, and Ω represents the set of homogeneous pixel SHP.

[0049] S4, using PM's matrix-vector multiplication, iteratively estimate the main eigenvalues ​​and corresponding eigenvectors of the sample covariance matrix C. The specific process is:

[0050] S41, the interference phase mean calculated by formula (1) As an initial estimate of the main eigenvector;

[0051] S42, according to the matrix-vector multiplication of PM, iteratively update the main eigenvalue and eigenvector, after After iterations:

[0052] (3);

[0053] (4);

[0054] in, Indicates The feature vector obtained by the iteration is Indicates The eigenvalues ​​obtained by iteration;

[0055] S43, Normalized feature vector :

[0056] (5);

[0057] because , when the number of iterations As it approaches infinity, we have:

[0058] (6);

[0059] (7);

[0060] in, is the theoretical principal eigenvalue of the sample covariance matrix C The corresponding feature vector.

[0061] S5. Set the adaptive termination strategy. Iteration steps It is impossible to take infinity, so the termination condition must be set for the iteration. The constructed adaptive iteration termination strategy is:

[0062] (8);

[0063] in, and Respectively represent Second and The main eigenvalue obtained by the iteration; Indicates the rate of change of the main eigenvalue. When the rate of change of the main eigenvalue is less than 1e-6, it means that the main eigenvalue and eigenvector have converged. It shows that the maximum number of iterations is 50. For high-quality distributed scatterers, DS can quickly converge to the exact solution, while for low-quality distributed scatterers, DS can terminate the iteration in time.

[0064] set up The basis is as follows:

[0065] In practical applications, after estimating the phase of the distributed scatterer DS, the quality of the estimation needs to be evaluated:

[0066] (9);

[0067] in, and denote the observed phase and estimated phase of the distributed scatterer DS, respectively. N represents the number of SAR images, Represents the real part operation. indicates the quality of the estimate, the larger the value, the better the estimate; When it is less than a certain value (such as ), indicating that the estimation quality is poor, and the distributed scatterer DS should be discarded. Actual SAR data show that, if Figure 2 Shown: For those For the distributed scatterer DS, only a few iterations (about 35 times) are needed to converge; The distributed scatterer DS is a low-quality target that should be discarded. Even if it does not converge, it will not affect the subsequent results.

[0068] S6, obtain the DS phase of the distributed scatterer, obtained when the iteration ends That is the optimal estimation of the DS phase of the distributed scatterer ,Right now: .

[0069] In order to verify the effectiveness and computational efficiency of the method described in this embodiment, the method described in this embodiment and the conventional eigendecomposition method are used to estimate the DS phase of the distributed scatterer, and the phase estimation accuracy and computational efficiency of the two methods are compared.

[0070] 1. Accuracy assessment

[0071] The simulation data contains 100 scenes of SAR data. According to the imaging parameters of Sentinel-1A, the time interval is set to 12 days, the wavelength is set to 55.6 mm, the deformation rate is set to 1 mm / a, and the SHP is set to 300. The method described in this embodiment and the conventional feature decomposition method are used to estimate the DS phase of the distributed scatterer, and each is repeated 100,000 times, and then the root mean square error (RMSE) is calculated to evaluate the accuracy. The experimental results are shown in Figure 2. Figure 3As shown, it can be found that the accuracy of the method described in this embodiment is completely consistent with the conventional method, which shows that the DS phase of the distributed scatterer estimated by the method described in the present invention is reliable.

[0072] 2. Computational efficiency test

[0073] The computation time is related to the number of SAR images. Five groups of SAR data sets with 100, 150, 200, 250 and 300 scenes were simulated. Then the DS phase was estimated using the method described in this embodiment and the conventional feature decomposition method. The calculation was repeated 10,000 times and the cumulative time was calculated. The experimental results are shown in Figure 2. Figure 4 As shown, it can be found that as the number of SARs increases, the computational efficiency of the method described in this embodiment is significantly higher than that of the conventional eigendecomposition method, and the computational time is only 1 / 10 of that of the conventional eigendecomposition method.

[0074] Therefore, the present invention adopts the above-mentioned distributed scatterer phase fast estimation method based on improved PM, and takes the interference phase mean on the homogeneous pixel SHP as the iterative initial value of PM, so as to avoid PM divergence, reduce the number of iterations and improve the convergence speed. By setting an adaptive termination strategy, it is ensured that the algorithm can quickly obtain the optimal solution for different distributed scatterers DS. Compared with the existing methods, the computational efficiency of the method described in the present invention is greatly improved, which effectively solves the problem of excessively time-consuming calculations, and is particularly suitable for parallel computing.

[0075] Finally, it should be noted that the above embodiments are only used to illustrate the technical solution of the present invention rather than to limit it. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that they can still modify or replace the technical solution of the present invention with equivalents, and these modifications or equivalent replacements cannot cause the modified technical solution to deviate from the spirit and scope of the technical solution of the present invention.

Claims

1. A distributed scatterer phase fast estimation method based on improved PM, characterized in that: The following steps are involved: S1. Using the two-sample hypothesis test method, according to the amplitude information of the registered SAR image, the distributed scatterer DS is identified and the homogeneous pixel SHP of each distributed scatterer DS is marked; S2. Calculate the mean interference phase of the distributed scatterer DS using homogeneous pixel SHP; the mean interference phase on the homogeneous pixel SHP is set as the initial estimate of the phase of the distributed scatterer DS , and used as the initial value of PM iteration to avoid PM divergence and reduce the number of iterations, the initial value of distributed scatterer DS phase iteration The calculation method is: = (1)； in, Indicates The interference phase on the interference pair, M Indicates the number of homogeneous pixels SHP; S3, using homogeneous pixel SHP to construct the sample covariance matrix C; specifically: for each distributed scatterer DS, use its homogeneous pixel SHP to construct the covariance matrix C, that is: (2)； Among them, the superscript represents the conjugate transpose of the matrix, y is the complex reflectance of each homogeneous pixel SHP in the time series, and Ω represents the set of homogeneous pixel SHPs; S4, using PM's matrix-vector multiplication to iteratively estimate the principal eigenvalues ​​and corresponding eigenvectors of the sample covariance matrix C; S5, construct an adaptive termination strategy for PM; S6, obtain the DS phase of the distributed scatterer, through Iterations, the result obtained at the end of the iteration That is the optimal estimation of the DS phase of the distributed scatterer ,Right now: .

2. The method for fast phase estimation of distributed scatterers based on improved PM according to claim 1, characterized in that: In S1, the two-sample hypothesis testing method is one of Kolmogorov-Smirnov and Anderson-Darling.

3. The method for fast phase estimation of distributed scatterers based on improved PM according to claim 1, characterized in that: S4 specifically includes the following steps: S41, the interference phase mean calculated by formula (1) As an initial estimate of the main eigenvector; S42, according to the matrix-vector multiplication of PM, iteratively update the main eigenvalue and eigenvector, after After iterations: (3)； (4)； in, Indicates The feature vector obtained by the iteration is Indicates The eigenvalues ​​obtained by iteration; S43, Normalized feature vector : (5)； because , when the number of iterations As it approaches infinity, we have: (6)； (7)； in, is the theoretical principal eigenvalue of the sample covariance matrix C The corresponding feature vector.

4. The method for fast phase estimation of distributed scatterers based on improved PM according to claim 3 is characterized in that: In S5, the number of iterations It is impossible to take infinity, and the termination condition must be set for the iteration. The constructed adaptive iteration termination strategy is: (8)； in, and They represent the main eigenvalues ​​obtained at the kth and k-1th iterations respectively; It represents the rate of change of the main eigenvalue. When the rate of change of the main eigenvalue is less than 1e-6, it means that the main eigenvalue and eigenvector have converged. k=50 indicates that the maximum number of iterations is 50.

Images