Refined fast large-scale phase unwrapping method compatible with multiple optimization models

Abstract

The invention provides a refined fast large-scale phase unwrapping method compatible with multiple optimization models. Before clustering of the residual points, a mask is generated by using a pseudo-coherence coefficient, and the residual points covered by the mask are used as initial residual point classes to improve the clustering speed; in the clustering process, according to the envelope sparse theorem, clustering analysis is carried out on the residual points and initial residual point classes formed through a mask, it is guaranteed or approximately guaranteed that a locally optimal solution and a globally optimal solution are consistent, and errors of a low-quality area are prevented from being diffused to a high-quality area; an efficient convex hull is adopted to mark a clustering region, and rapid integral solution is carried out on a region outside the convex hull; a proper optimization model can be flexibly selected to locally solve the region in the convex hull according to user requirements; the high-quality area and the low-quality area in the convex hull are divided more finely, solving of the low-quality area in the convex hull is abandoned, and meanwhile it is guaranteed that the high-quality area can still be effectively solved.

Description

Technical Field

[0001] The present invention belongs to the technical field of radar imaging, and in particular relates to an Elaborate and Fast Large-Scale Phase Unwrapping Method with Compatible Multiple Optimization Models (EFLS). Background Art

[0002] Interferometric synthetic aperture radar (InSAR) is a powerful and mature remote sensing technology that can generate digital elevation models and measure surface deformation with high spatial resolution. In practical applications, the measured data obtained by the InSAR system is a wrapped phase map, in which each pixel is represented by a wrapped phase in the range of (-π,π]. The wrapped phase is the main value of the absolute phase that contains elevation and / or deformation information. Phase unwrapping is the step in the InSAR processing chain to obtain the absolute phase from the wrapped phase. So far, two main types of InSAR technology have been developed: single baseline and multi-baseline InSAR.

[0003] For the above two InSAR technologies, the principles used in phase unwrapping are different. Single baseline phase unwrapping relies on a single wrapped phase image to recover the absolute phase. It uses the phase continuity assumption, assuming that the absolute phase gradient of any two adjacent pixels is less than π. Therefore, single baseline phase unwrapping is suitable for terrain with good spatial continuity. Multi-baseline phase unwrapping combines multiple wrapped phase images and uses the relationship between multiple baselines to solve the absolute phase. Since the diversity of baselines can significantly increase the ambiguity interval of the wrapped phase in phase unwrapping, multi-baseline phase unwrapping is more suitable for discontinuous terrain, such as valleys and cliffs.

[0004] With the rapid development of InSAR technology, the amount of data that can be obtained by InSAR systems has increased dramatically, and the demand for large-scale data processing is increasing. In some applications, very large and continuous spatial coverage is required. For phase unwrapping processing, it may be necessary to process long strips of wrapped phase images that are several times larger than the standard wrapped phase image. For traditional single-baseline and multi-baseline phase unwrapping methods, when the amount of input data exceeds a certain limit, two problems may arise: (1) insufficient memory; (2) significantly extended calculation time; this is because the computational complexity of traditional phase unwrapping methods is high, and the computing resources used in large-scale phase unwrapping increase nonlinearly. For example, the minimum cost flow algorithm is a widely used phase unwrapping method, but its computational complexity is highly nonlinear. At the same time, since multi-baseline phase unwrapping requires the joint processing of multiple wrapped phase images, these two problems are more serious in multi-baseline phase unwrapping. Under existing technical conditions, even if the performance of the hardware is improved, such as increasing the number and speed of CPUs and the amount of memory, it cannot meet the rapidly growing computational requirements of large-scale phase unwrapping.

[0005] The basic idea of ​​reducing the computational resources of large-scale data phase unwrapping is to divide and conquer, that is, to divide the large-scale wrapped phase image into sub-wrapped phase images (hereinafter referred to as sub-images) and solve them separately, and then splice the phase unwrapping results of each sub-image together to obtain the overall phase unwrapping result. There are two ways to divide the sub-images: regular division and irregular division. The current more effective method is the irregular sub-image division method based on residual point clustering analysis, combined with the large-scale phase unwrapping method of irregular sub-block solution. In order to improve the solution efficiency, the convex hull technology can be further used to generate a mask (the area inside the mask is not solved). Although the strategy of abandoning the solution of the convex hull area can significantly improve the computing speed and reduce memory consumption, since the area inside the convex hull is not an information-free area, this abandonment will lose effective information. Therefore, in order to avoid the loss of effective information in the area inside the convex hull, it is necessary to develop a more sophisticated solution method. Summary of the invention

[0006] In order to solve the above problems existing in the prior art, the present invention provides a refined fast large-scale phase unwrapping method (EFLS) that is compatible with multiple optimization models. The technical problem to be solved by the present invention is achieved by the following technical solutions:

[0007] A refined, fast, large-scale phase unwrapping method compatible with multiple optimization models includes:

[0008] S100, calculating a pseudo coherence coefficient for a winding phase map to be solved, and generating a mask for a connected domain whose pseudo coherence coefficient is lower than a specified threshold; the winding phase map includes a single-baseline winding phase map and a multi-baseline winding phase map;

[0009] S200, according to different situations of the winding phase diagram, using corresponding strategies to perform gradient estimation, and then calculating residual points to obtain a residual point diagram;

[0010] S300, performing cluster analysis on the residual points in the residual point graph using the envelope sparsity theorem and the mask to obtain a clustering result of the residual points; wherein the residual points covered by each unconnected mask corresponding to the residual point graph are used as initial residual point classes;

[0011] S400, forming a convex hull of the residual point class according to the clustering result, and marking a convex hull area according to the convex hull points in the residual point map;

[0012] S500, performing phase unwrapping on the area outside the convex hull in the wrapped phase image by using a method for solving the path integral; then, according to the accuracy and efficiency requirements of the solution, using a predetermined quality map and a specified threshold, dividing the area inside the convex hull into high-quality and low-quality areas, and selecting a suitable phase unwrapping algorithm to perform phase unwrapping on the high-quality area inside the convex hull;

[0013] S600, splicing the phase unwrapping results of the area inside the convex hull and the area outside the convex hull to obtain a final phase unwrapping result.

[0014] Beneficial effects:

[0015] (1) After the present invention generates a mask using the pseudo-coherence coefficient, the residual points in each connected domain covered by the mask are used as the initial residual point class, which can save the clustering time required for the residual points to form the initial residual point class and improve the clustering speed; then, according to the envelope sparsity theorem, the residual points in the residual point graph and the initial residual point class formed by the mask are clustered and analyzed to prevent the errors in the low-quality area in the convex hull from spreading to the high-quality area, thereby ensuring the solution accuracy.

[0016] (2) The EFLS method proposed in the present invention first solves the high-quality area outside the convex hull with extremely high efficiency, so that the user can obtain reliable phase unwrapping information of the high-quality area in a very short time; then, the area inside the convex hull is solved. Therefore, this is a phase unwrapping algorithm that has both front-end fast-viewing and back-end fine processing capabilities, and the algorithm is more practical.

[0017] (3) The present invention divides the area in the convex hull into high-quality and low-quality areas more finely. Due to the randomness of the residual point distribution, the convex hull contains not only low-quality areas but also some high-quality areas; therefore, the present invention can finely divide the convex hull area according to the pseudo-coherence coefficient, abandon the solution of the low-quality area in the convex hull to reduce the amount of calculation, and at the same time, the high-quality area in the convex hull can still be effectively solved.

[0018] (4) The present invention can flexibly select appropriate optimization models for solving the phase unwrapping of single baselines and multiple baselines in the convex hull area, for example, ODLS, L 0 Norm, L 1 Norm, L ∞ The invention can be used to perform image stitching correctly, and ensure that the optimal solution of the region in the convex hull is consistent with the global optimal solution. The invention will be further described in detail below with reference to the accompanying drawings and embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0019] Figure 1 A flow chart of a refined, fast, large-scale phase unwrapping method that is compatible with multiple optimization models provided by the present invention;

[0020] Figure 2 It is a schematic diagram of a residual point auxiliary example used in the present invention to illustrate the residual point theory;

[0021] Figure 3 Google Earth optical images and multi-baseline entangled phase diagrams of measured data;

[0022] Figure 4 The comparison results between EFLS clustering method and traditional clustering method are shown;

[0023] Figure 5 The result diagram is a comparison between the solution of the high-quality area in the convex hull of EFLS and the convex hull solution of EFLS;

[0024] Figure 6 Phase unwrapping results and difference maps of the convex hull region solved using a variety of optimization models.

[0025] The present invention is further described in detail below with reference to specific embodiments, but the embodiments of the present invention are not limited thereto.

[0026] Before introducing the present invention, the technical concept of the present invention is first introduced.

[0027] Existing large-scale phase unwrapping methods mainly adopt the "divide and conquer" strategy, that is, dividing the large-scale wrapped phase map into sub-maps for separate solutions and then splicing. The traditional mask generation method based on the quality map is difficult to prevent the solution errors that may occur in the mask area from spreading to the high-quality area; the block method based on clustering analysis and convex hull technology chooses to abandon the strategy of solving the area within the convex hull in order to improve the solution speed, thereby having to lose the utilization rate of the information within the convex hull. In addition, the existing technology uses a single model in a single processing, and cannot process the local different landforms and image features in large-scale data separately. Therefore, the present invention needs to design a block solution method that meets the following conditions:

[0028] (1) After the winding phase map is divided into sub-maps and solved separately, it can be ensured that the solution results of the sub-maps can be correctly spliced, and the solution errors of low-quality areas in the sub-maps will not spread to high-quality areas;

[0029] (2) Based on the existing block solution method based on cluster analysis and convex hull technology, in order to improve information utilization, the high-quality area within the convex hull should also be effectively solved;

[0030] (3) When solving the high-quality region in the convex hull, no additional resource consumption will be added due to the fine division of high-quality and low-quality regions for solution. In other words, the ratio of the computing resources consumed by solving the high-quality region in the convex hull to the computing resources consumed by solving the convex hull region as a whole should at least be less than the ratio of the area of ​​the solved region in the convex hull to the area of ​​the convex hull;

[0031] (4) For the area inside the convex hull, a variety of optimization models can be flexibly selected, and it is ensured that the solution results of the area inside the convex hull can be seamlessly spliced ​​with the surrounding areas.

[0032] Therefore, based on the residual point and envelope sparsity theorems, the present invention adopts a research approach that combines theoretical methods with simulations, and proposes a refined fast large-scale phase unwrapping (EFLS) method that is compatible with multiple optimization models to achieve efficient phase unwrapping solutions for high-quality areas in large-scale interferometric data.

[0033] The idea of ​​achieving the purpose of the present invention is: for single-baseline winding phase diagram and multi-baseline winding phase diagram, the phase continuity hypothesis and the Chinese remainder theorem are used to solve the corresponding residual point diagrams; before performing residual point clustering, a mask is first generated using a quality map (for example, a pseudo-coherence coefficient) and a specified threshold, and then the residual points covered by the mask are used as the initial residual point class for clustering to improve the clustering speed; in the clustering process, according to the envelope sparsity theorem, cluster analysis is performed on the residual points and the initial residual point class formed by the mask to ensure (or approximately ensure) that the local optimal solution is consistent with the global optimal solution, and prevent errors in low-quality areas from spreading to high-quality areas; an efficient convex hull is used to mark the clustering area, and the area outside the convex hull is quickly integrated and solved; according to user needs, a suitable optimization model is flexibly selected to solve the area inside the convex hull; in order to improve information utilization, the high- and low-quality areas in the convex hull are divided more finely, and the solution of the low-quality area in the convex hull is abandoned, while ensuring that the high-quality area can still be effectively solved.

[0034] The technical solution details of the present invention are described in detail below.

[0035] refer to Figure 1 The present invention provides a refined, fast, large-scale phase unwrapping method that is compatible with multiple optimization models, including:

[0036] S100, calculating a pseudo coherence coefficient for a winding phase map to be solved, and generating a mask for a connected domain whose pseudo coherence coefficient is lower than a specified threshold; the winding phase map includes a single-baseline winding phase map and a multi-baseline winding phase map;

[0037] S200, according to different situations of the winding phase diagram, using corresponding strategies to perform gradient estimation, and then calculating residual points to obtain a residual point diagram;

[0038] S300, performing cluster analysis on the residual points in the residual point graph using the envelope sparsity theorem and the mask to obtain a clustering result of the residual points; wherein the residual points covered by each unconnected mask corresponding to the residual point graph are used as initial residual point classes;

[0039] S400, forming a convex hull of the residual point class according to the clustering result, and marking a convex hull area according to the convex hull points in the residual point map;

[0040] S500, performing phase unwrapping on the area outside the convex hull in the wrapped phase image by using a method for solving the path integral, and then dividing the area inside the convex hull into high-quality and low-quality areas more finely by using a predetermined quality map and a specified threshold according to the accuracy and efficiency requirements of the solution, and selecting a suitable phase unwrapping algorithm to perform phase unwrapping on the high-quality area inside the convex hull;

[0041] Therein, the predetermined quality map is formed using pseudo-coherence coefficients.

[0042] S600, splicing the phase unwrapping results of the area inside the convex hull and the area outside the convex hull to obtain a final phase unwrapping result.

[0043] In an optional implementation manner of the present invention, S200 includes:

[0044] S210, for the single baseline phase unwrapping problem, estimating the fuzzy number gradient of the single baseline wrapped phase image using the phase continuity assumption;

[0045] The process of estimating the fuzzy number gradient for the single baseline winding phase image in the present invention is as follows:

[0046] S211, determining, based on a phase continuity assumption, that an absolute phase gradient between adjacent pixels of the single baseline wrapping phase image is less than π;

[0047] S212, based on the requirement that the absolute phase gradient is less than π, the fuzzy number gradient of the single baseline winding phase image is estimated using the phase continuity assumption, expressed as:

[0048]

[0049] in, is the winding phase of any pixel s, is the wrapping phase of the adjacent pixel t to the right of pixel s, is the wrapping phase of the adjacent pixel u below pixel s; is the lateral fuzzy number gradient between adjacent columns of pixel s and pixel t, is the longitudinal fuzzy number gradient between adjacent rows of pixel s and pixel u; r represents the serial number of the interferogram and baseline, r = 1, 2, …, M, M is the number of baselines in the multi-baseline InSAR system, for single baseline phase unwrapping M = 1, so r = 1.

[0050] S220, for the multi-baseline phase unwrapping problem, the Chinese remainder theorem is used to estimate the fuzzy number gradient of the multi-baseline wrapped phase image;

[0051] The process of estimating the fuzzy number gradient for the multi-baseline winding phase image of the present invention is as follows:

[0052] S221, determining the winding phase map of the multi-baseline InSAR system (taking a dual baseline as an example, i.e., M=2), the relationship between pixels at corresponding positions of the winding phase maps of different baselines is expressed as

[0053]

[0054] Among them, B1 and B2 are the winding phases obtained by the dual baseline InSAR system. Figure 1 and winding phase Figure 2 The lengths of the two corresponding vertical baselines (hereinafter referred to as baseline lengths), The winding phase Figure 1 The winding phase at pixel s in the figure, k1(s) is the winding phase Figure 1 The fuzzy number at pixel s in , The winding phase Figure 2 The winding phase at pixel s in the figure, k2(s) is the winding phase Figure 2 The blur number at pixel s in ;

[0055] S222, determine the relationship between the pixel t that is horizontally adjacent to the pixel s and the pixel u that is vertically adjacent to the pixel s, assuming that the arrangement of the adjacent 4 pixels is as follows Figure 2 As shown, the same is true for the horizontally adjacent pixel t and the vertically adjacent pixel u of the pixel s, which are respectively expressed as

[0056]

[0057]

[0058] Subtracting formula (4) and (5) from formula (3) respectively, we get the following derivation formula:

[0059]

[0060]

[0061] in, s, t and u all represent pixels. The horizontal adjacent pixel of pixel s is t, and the vertical adjacent pixel of pixel s is u; k r (s), k r (t) and k r (u) is the ambiguity number at these pixel positions; r = 1, 2, ..., M, for the dual baseline InSAR system, M = 2, so r = 1, 2;

[0062] S223, using the Chinese Residue Theorem to estimate the fuzzy number gradient of the multi-baseline winding phase image; the optimization model used by the Chinese Residue Theorem to solve the fuzzy number gradient of the multi-baseline winding phase image is:

[0063]

[0064] in,

[0065]

[0066] and are decision variables, which are integers; for some special combinations of (B1, B2), within a certain range, a unique solution can be found.

[0067] S230, calculating the residual points of the winding phase graph using the fuzzy number gradient, and forming a residual point graph;

[0068] This step uses the fuzzy number gradient to calculate the residual points of the winding phase image, assuming that the arrangement of the adjacent 4 pixels is as follows Figure 2 As shown, the calculation formula of the residual point is expressed as:

[0069]

[0070] Where N represents the polarity of the residual point. For single-baseline residual points, the value of N is -1, 0, or 1. For multi-baseline residual points, the value of N is an integer. and They all represent the fuzzy number gradient of the corresponding pixel position.

[0071] In an optional implementation manner of the present invention, S300 includes:

[0072] The residual points in each connected domain covered by the mask in S100 are taken as the initial residual point class, and the residual points in the residual point graph and the initial residual point class formed by the mask are clustered according to the envelope sparsity theorem to obtain the clustering results of the residual points, including:

[0073] 3a), set the initial residual point class: for each residual point res i , set its cluster radius β i The initial cluster radius β i0 Among them, the residual points covered by a connected mask area in the residual point graph are regarded as an initial residual point class; those residual points not covered by the connected mask area, each residual point is regarded as an initial residual point class;

[0074] 3b), calculate the distance between the residual points: in any two different residual point clusters n and cluster m Any residual point res i and res j , calculate res i and res j The distance between ij , where res i ∈cluster n ,res j ∈cluster m ;

[0075] 3c), determine whether the class is merged: determine res i and res j The distance between ij Is it satisfied? If satisfied, res i and res j The two residual point clusters n and cluster m merged into a residual point class; where κ is the norm of the clustering algorithm, if L is used 0 norm, then κ=0. If L 1 norm, then κ=1; α is the approximation factor, which is used to control the accuracy of the envelope sparsity theorem;

[0076] 3d), calculation of the total polarity of the residual point class: for each residual point class cluster n If cluster n If the polarity of the internal residual point is 0, the cluster n The clustering radius of the inner residual points is set to β i =r n , where r n For cluster n The envelope distance refers to the maximum distance between two residual points in the residual point class; if cluster n If the polarity of the internal residual point is not 0, add clustern The clustering radius of the inner residual points is β i , β i Should be greater than r n The value of

[0077] 3e), determine whether the clustering iteration is finished: if there is a residual point class with unbalanced polarity, go directly to 3b) and re-cluster; otherwise, determine again whether any two residual points belonging to different residual point classes are equal. i and res j The distance between ij : If exists If , go to 3b) and re-cluster; if for all the residual points res belonging to different residual point classes i and res j The distance between ij ,That The clustering is completed and the clustering result is obtained.

[0078] In an optional implementation manner of the present invention, S500 includes:

[0079] 5a) For the area outside the convex hull, phase unwrapping is performed using the path integral method;

[0080] 5b) For the area in the convex hull, mark the high-quality area in the convex hull according to the pseudo-coherence coefficient value and the given threshold, and discard the solution of the low-quality area in the convex hull to reduce the amount of calculation;

[0081] 5c) For the high-quality area in the convex hull, according to the accuracy and efficiency requirements of the solution, select a suitable phase unwrapping algorithm to solve and obtain the fuzzy number k of each pixel in the convex hull. r (s);

[0082] 5d) Using the fuzzy number k of each pixel in the convex hull r (s) and winding phase Solve the absolute phase ψ of the pixels in the convex hull r (s); the formula is expressed as:

[0083]

[0084] The phase unwrapping algorithm includes L P Algorithm, assuming that the arrangement of 4 adjacent pixels is as follows Figure 2 As shown, L P The algorithm can be expressed as:

[0085]

[0086] in, and is the weight at pixel s, which can be determined by the pseudo-coherence coefficient value of the input winding phase image; set W is the pixel set of each convex hull area to be solved, s, t and u are all pixel points in set W, where the horizontal adjacent pixel of pixel s is t, and the vertical adjacent pixel of pixel s is u; k r (s), k r (t) and k r (u) is the fuzzy number of pixels s, t, and u, representing the decision variable; is the fuzzy number gradient of the pixel, which is an integer; r = 1, 2, ..., M, represents the serial number of the winding phase map to be solved; when p takes different integer values, it represents different L P Norm algorithms, such as the methods that can be used are L 0 Norm, L 1 Norm, L ∞ Norms, etc.

[0087] In an optional embodiment of the present invention, 5c) includes:

[0088] 5c1) Before solving the area inside the convex hull, expand the area inside the convex hull outward by one pixel;

[0089] 5c2) When solving the area inside the convex hull, the pixels on the convex hull boundary are used as the equality constraint of the algorithm, and then the phase unwrapping algorithm is used to solve it. Since the phase unwrapping is performed on the area outside the convex hull first, the absolute phase values ​​of the pixels on the convex hull boundary are known; for L P In the case of the norm (p ≥ 2), in order to make the values ​​on the boundary between the convex hull area and the area outside the convex hull equal, the constraint conditions on the boundary can be implemented by using a weighted method.

[0090] The solution result of the inner region of the convex hull of the present invention can be directly used for image stitching with the outer region of the convex hull without the need for additional calculations. Because the present invention preferentially performs phase unwrapping on the outer region of the convex hull, in order to ensure that the pixels inside and outside the convex hull can be correctly stitched, when solving the high-quality region in the convex hull, the pixels on the boundary of the convex hull are used as the equality constraint of the algorithm before solving.

[0091] The implementation details of the device implementation part of the present invention are consistent with those of the method part, and will not be repeated here.

[0092] experiment:

[0093] This experiment uses a real InSAR dataset, which is acquired by the Phased Array type L-band Synthetic Aperture Radar (PALSAR) InSAR system of the Advanced Land Observing Satellite (ALOS). The main radar parameters are shown in Table 1.

[0094] Table 1 Main parameters of the measured data set

[0095]

[0096] The observation area is Crowley Lake and its surrounding areas in the United States. Figure 3 Figure (a) shows the Google Earth optical image of the area. The interferometric phase data set to be solved is as follows: Figure 3 As shown in Figures (b), (c), and (d), they are multi-baseline entangled phase images with baseline lengths of -406.2m, -957.72m, and -1862.3m, respectively. Comparing the entangled phase image with the optical image, it can be seen that the area covered by the lake surface forms an obvious noise area in the interference phase image. Such a large area of ​​noise area will consume a lot of computing resources in the traditional solution method. The mountains around the lake also pose a challenge to the solution of the data. Next, this experiment will verify the superiority of the EFLS method proposed in this invention over the traditional method from the following three aspects:

[0097] 1. Clustering speed improvement verification:

[0098] Comparing the clustering time of the traditional clustering method and the EFLS method, it is verified that the residual point clustering method proposed in the present invention improves the clustering speed.

[0099] 2. Refined solution and verification of computing resource consumption:

[0100] The EFLS method can accurately divide the high / low quality regions in the winding phase diagram. Compared with the traditional solution method that does not divide the high / low quality regions, the EFLS method requires significantly reduced computing resources.

[0101] 3. Verification of local area processing flexibility:

[0102] For complex terrain, different data features may be suitable for different algorithms. The EFLS method allows the use of different L P Norm or other optimization model methods independently process different convex hull areas in an interference pattern, which can meet the user's needs for flexible data processing.

[0103] The following will verify the above three aspects in turn.

[0104] Experiment 1:

[0105] This experiment compares the clustering time of the traditional clustering method and the EFLS method of the present invention, and verifies that the residual point clustering method of the EFLS method improves the clustering speed.

[0106] like Figure 4 Figure (a) shows Figure 3 The pseudo coherence coefficient map in (b) is shown in Figure 1. In this experiment, the threshold is set to 0.7, and the area with a pseudo coherence coefficient less than 0.7 is used as a mask. The generated mask map is shown in Figure 1. Figure 4 As shown in the middle (b), the mask area is marked with 1, that is, the value of the low-quality area in the image is 1, and the value of the high-quality area is 0.

[0107] Figure 4 Figure (c) is composed of three entangled phase diagrams (such as Figure 3 The residual point map obtained by joint solution (as shown in Figures (b) to (d) in the figure). Comparing the mask map and the residual point map, it can be seen that the area with dense residual points basically overlaps with the mask area. Therefore, this experiment uses Figure 4 The residual points covered by the mask in (b) form the initial residual point class. Figure 4 In the middle (c) figure, the convex hull of the initial residual point class is marked with a black solid line (for the sake of image clarity, only the initial residual point class containing more than 20 residual points is marked, and the subsequent convex hull marking is also the same). The traditional clustering algorithm uses each residual point as an initial residual point class, while the EFLS method of the present invention uses the pseudo correlation coefficient to form the initial residual point class to save the number of clustering iterations. The clustering speed and clustering results of the two are shown in Table 2:

[0108] Table 2 Comparison of clustering speed and clustering results

[0109]

[0110] The traditional clustering method has 1735 clustering iterations and takes 1846.79 seconds, while the EFLS method has 1593 clustering iterations and takes 1197.05 seconds, reducing the clustering time to 64.8% of the traditional method. Of course, the improvement in clustering speed comes at the expense of clustering precision. Figure 4 Figure (d) is the convex hull diagram of the clustering result of the EFLS method. Figure 4Figure (e) is the clustering result of the traditional method. By comparing the convex hulls in the two figures, we can see that the residual point class has increased relatively. As shown in Table 2, the average number of residual points in the traditional residual point class is 15.75, while that of the EFLS method is 20.43. It is not difficult to understand that since all residual points must be contained in the convex hull, the fewer the number of convex hulls, the larger the area covered by the convex hull (the proportion of the number of pixels in the convex hull obtained by the traditional clustering method to the number of pixels in the entire image is 20.65%, while the proportion of the convex hull of the EFLS method is 21.74%). Since the traditional method will abandon the convex hull solution, the increase in the proportion of the convex hull will lead to a decrease in information utilization. However, the increase in the total area covered by the convex hull will hardly affect the information utilization of EFLS, because the total area of ​​the solved area depends on the generation of the mask, not on the size of the convex hull. Comparison Figure 4 The convex hull marks in Figure (d) and Figure (e) show that the clustering results obtained by the two methods are basically the same. Therefore, the EFLS method effectively improves the clustering efficiency without having a significant impact on the clustering results. This improvement will be more obvious when processing measured data with a large range of noise.

[0111] Experiment 2:

[0112] Based on the convex hull generation results of Experiment 1, Experiment 2 will verify that the EFLS method has refined solution capabilities and reasonable computing resource consumption.

[0113] To facilitate comparative testing, this experiment only tested the lake and its surrounding areas, such as Figure 5 As shown in Figure (a), the convex hull with the largest area and the convex hulls connected to it are marked in the pseudo-correlation coefficient diagram with white wireframes. The area surrounded by these convex hulls is the area to be tested in this experiment. In the following text, the marked area is called the lake area. Obviously, there are also high-quality image areas in the lake area, that is, the blue area in the white wireframe in the figure. The EFLS method can solve the high-quality area in the irregular lake area while ensuring (or approximately ensuring) that the local solution is consistent with the global solution. The following is an experimental comparison of this area.

[0114] Figure 5 (b) is Figure 3 The reference absolute phase in Figure (b) is used to compare the solution accuracy. Figure 5 The middle (c) figure is L 1 The unwrapping result of the lake area solved by the norm method. In order to keep the image intact and facilitate comparison, other convex hulls outside the lake area are also solved by this method. Figure 5 The middle (d) figure is Figure 5 The difference between (b) and (c). Figure 5 The middle figure (e) shows the high-quality area in the lake area using L 1Norm solution results. The low-quality area is not solved and its value is set to 100 for observation. Figure 5 Figure (f) shows Figure 5 The difference between (b) and (e), where the difference in the unsolved area is set to 0. Figure 5 The color bars of (d) and (f) are set to the same value range. The accuracy of the high-quality area after the refined solution should not be lower than the accuracy of the convex hull solution. After calculation, Figure 5 The mean square errors of the high-quality areas in the lakes in (d) and (f) are 0.55 and 0.52 respectively, which shows that the solution accuracy of the method proposed in the present invention is guaranteed. Table 3 compares the computing resources consumed by the refined solution.

[0115] Table 3 Comparison of computational resource consumption of EFLS method for refined solution

[0116]

[0117] It can be seen from the data in Table 3 that after the improvement to the refined solution method, the solution time is reduced to 27.62% of the total solution of the convex hull, and the peak memory is also reduced to 23.19% of the original. The decrease in computing resource consumption is due to the reduction in the area of ​​the solved region. However, it should be noted that the realization of refined solution requires the consumption of additional computing resources, and the refined solution method proposed in the present invention does not increase too many computing resources, because the number of pixels solved is 38.73% of the number of pixels in the lake area, that is, the reduction in resource consumption is greater than the reduction in the area solved. Therefore, this experiment verifies that refined solution is effective in both improving solution accuracy and speed.

[0118] Experiment 3:

[0119] The observed area used in this experiment has a variety of terrains, in addition to the lake area focused on in the above experiments, such as Figure 3 As shown in Figure (a), the upper left corner is a complex mountainous terrain. The results obtained using different phase unwrapping optimization models will be different. Therefore, in the actual data processing, different optimization models may be used in different regions to obtain higher processing accuracy. Therefore, this experiment uses L 1 The two optimization models of norm and MIN are solved to verify the EFLS method in the convex hull region for L P Flexibility in norm model selection.

[0120] like Figure 6The mask image is shown in Figure (a). The area surrounded by the convex hull marked by multiple white solid lines in the upper left corner is a mountainous area with complex terrain. It is referred to as the mountainous area in the following text. Considering that the noise in this area is relatively serious and scattered, resulting in the formation of many small-area high-quality independent areas (called islands) between the masks, this experiment adopts the overall solution of the convex hull area to ensure the continuity of the solved area. This experiment uses L 1 The two optimization models of norm and MIN are used to solve the problem, and the remaining convex hulls are all based on L 1 Norm to solve. Figure 6 Figures (b) and (c) are L 1 The results of the norm and MIN methods for mountainous areas, Figure 6 Figures (d) and (e) are Figure 5 In the figure (b) and their difference graph, the mean square errors in the mountainous area are 24.97 and 24.46 respectively. The time taken for the two methods to solve is 3697.87 seconds and 6897.24 seconds respectively, and the peak memory is 263.11Mb and 408.62Mb respectively. Different optimization models are used for solving, and the accuracy obtained is different, and the computing resources consumed are also different. Therefore, the EFLS phase unwrapping algorithm proposed in the present invention allows users to flexibly select appropriate optimization models for local solutions in different regions according to application requirements, and this does not affect the correct image stitching of different convex hull regions.

[0121] It can be concluded from the above three experiments that the EFLS method proposed in the present invention first uses the pseudo-coherence coefficient to generate a mask when clustering residual points, and then uses the residual points covered by the mask as the initial class for clustering analysis, which improves the clustering speed compared with the traditional clustering algorithm; the EFLS method proposed in the present invention performs phase unwrapping on more high-quality areas, thereby improving data utilization; the EFLS method proposed in the present invention divides the convex hull into high / low quality areas and abandons the solution of the low-quality areas in the convex hull, which not only reduces the consumption of computing resources but also ensures the solution accuracy of most high-quality areas; the EFLS method proposed in the present invention can flexibly select different optimization models to perform local solutions to the areas in the convex hull, and different convex hull areas can perform correct image stitching.

[0122] In summary, experiments based on measured data have verified the correctness, effectiveness and reliability of the present invention.

[0123] In addition, the terms "first" and "second" are used for descriptive purposes only and should not be understood as indicating or implying relative importance or implicitly indicating the number of the indicated technical features. Therefore, the features defined as "first" and "second" may explicitly or implicitly include one or more of the features. In the description of the present invention, the meaning of "plurality" is two or more, unless otherwise clearly and specifically defined.

[0124] Although the present application is described herein in conjunction with various embodiments, in the process of implementing the claimed application, those skilled in the art may understand and implement other variations of the disclosed embodiments by reviewing the drawings, the disclosure, and the appended claims. In the claims, the word "comprising" does not exclude other components or steps, and "a" or "an" does not exclude a plurality of components or steps.

[0125] The above contents are further detailed descriptions of the present invention in combination with specific preferred embodiments, and it cannot be determined that the specific implementation of the present invention is limited to these descriptions. For ordinary technicians in the technical field to which the present invention belongs, several simple deductions or substitutions can be made without departing from the concept of the present invention, which should be regarded as falling within the protection scope of the present invention.

Claims

1. A refined, fast, large-scale phase unwrapping method compatible with multiple optimization models, characterized in that: include: S100, calculating a pseudo coherence coefficient for the winding phase image to be solved, and generating a mask for the connected domains whose pseudo coherence coefficients are lower than a specified threshold; The winding phase diagram includes a single-baseline winding phase diagram and a multi-baseline winding phase diagram; S200, according to different situations of the winding phase diagram, using corresponding strategies to perform gradient estimation, and then calculating residual points to obtain a residual point diagram; S300, performing cluster analysis on the residual points in the residual point graph using the envelope sparsity theorem and the mask to obtain a clustering result of the residual points; wherein the residual points covered by each unconnected mask corresponding to the residual point graph are used as initial residual point classes; S400, forming a convex hull of the residual point class according to the clustering result, and marking a convex hull area according to the convex hull points in the residual point map; S500, performing phase unwrapping on the area outside the convex hull in the wrapped phase image by using a method for solving the path integral; then, according to the accuracy and efficiency requirements of the solution, using a predetermined quality map and a specified threshold, dividing the area inside the convex hull into high-quality and low-quality areas, and selecting a phase unwrapping algorithm to perform phase unwrapping on the high-quality area inside the convex hull; S600, splicing the phase unwrapping results of the area inside the convex hull and the area outside the convex hull to obtain a final phase unwrapping result.

2. The refined, fast, large-scale phase unwrapping method compatible with multiple optimization models according to claim 1, characterized in that: S200 includes: S210, for the single baseline phase unwrapping problem, estimating the fuzzy number gradient of the single baseline wrapped phase image using the phase continuity assumption; S220, for the multi-baseline phase unwrapping problem, the Chinese remainder theorem is used to estimate the fuzzy number gradient of the multi-baseline wrapped phase image; S230, using fuzzy number gradient to calculate the residual points of the winding phase graph and form a residual point graph.

3. The refined, fast, large-scale phase unwrapping method compatible with multiple optimization models according to claim 2 is characterized in that: S210 includes: S211, determining, based on a phase continuity assumption, that an absolute phase gradient between adjacent pixels of the single baseline wrapping phase image is less than π; S212, based on the requirement that the absolute phase gradient is less than π, the fuzzy number gradient of the single baseline winding phase image is estimated using the phase continuity assumption, expressed as: in, is the winding phase of any pixel s, is the wrapping phase of the adjacent pixel t to the right of pixel s, is the wrapping phase of the adjacent pixel u below pixel s; is the lateral fuzzy number gradient between adjacent columns of pixel s and pixel t, is the longitudinal fuzzy number gradient between adjacent rows of pixel s and pixel u; r represents the serial number of the wrapped phase map and baseline, r = 1, 2, …, M, M is the number of baselines in the multi-baseline InSAR system, for single baseline phase unwrapping M = 1, so r = 1.

4. The refined, fast, large-scale phase unwrapping method compatible with multiple optimization models according to claim 3 is characterized in that: S220 includes: S221, determining the relationship between pixels at corresponding positions of the winding phase images of different baselines in the winding phase image of the multi-baseline InSAR system, expressed as Among them, B1 and B2 are the lengths of the two vertical baselines corresponding to the winding phase image 1 and the winding phase image 2 obtained by the dual-baseline InSAR system, is the winding phase at pixel s in the winding phase image 1, k1(s) is the fuzzy number at pixel s in the winding phase image 1, is the winding phase at pixel s in the winding phase image 2, k2(s) is the fuzzy number at pixel s in the winding phase image 2; S222, determine the relationship between the pixel t that is horizontally adjacent to the pixel s and the pixel u that is vertically adjacent to the pixel s, which are expressed as Subtracting formula (4) and (5) from formula (3) respectively, we get the following derivation formula: in, s, t and u all represent pixels. The horizontal adjacent pixel of pixel s is t, and the vertical adjacent pixel of pixel s is u; k r (s), k r (t) and k r (u) is the ambiguity number at these pixel positions; r = 1, 2, ..., M, for the dual baseline InSAR system, M = 2, so r = 1, 2; S223, using the optimization model of the Chinese remainder theorem to solve the fuzzy number gradient of the multi-baseline winding phase diagram, the optimization model is: in, and are decision variables, which are integers; for some special combinations of (B1, B2), within a certain range, a unique solution can be found.

5. The refined fast large-scale phase unwrapping method compatible with multiple optimization models according to claim 4 is characterized in that: S230 includes: The residual points of the winding phase graph are calculated using the fuzzy number gradient to form a residual point graph. The calculation formula of the residual points is expressed as: Where N represents the polarity of the residual point. For single-baseline residual points, the value of N is -1, 0, or 1. For multi-baseline residual points, the value of N is an integer. and They all represent the fuzzy number gradient of the corresponding pixel position.

6. The refined fast large-scale phase unwrapping method compatible with multiple optimization models according to claim 5 is characterized in that: S300 includes: The residual points in each connected domain covered by the mask in S100 are taken as initial residual point classes, and cluster analysis is performed on the residual points in the residual point graph and the initial residual point classes formed by the mask according to the envelope sparsity theorem to obtain a clustering result of the residual points.

7. The refined fast large-scale phase unwrapping method compatible with multiple optimization models according to claim 6 is characterized in that: S300 includes: 3a), set the initial residual point class: for each residual point res i , set its cluster radius β i is the initial cluster radius β i0 ; Among them, the residual points covered by a connected mask area in the residual point graph are regarded as an initial residual point class; those residual points not covered by the connected mask area, each residual point is regarded as an initial residual point class; 3b), calculate the distance between the residual points: in any two different residual point clusters n and cluster m Any residual point res i and res j , calculate res i and res j The distance between ij , where res i ∈cluster n ,res j ∈cluster m ; 3c), determine whether the residual point class is merged: determine res i and res j The distance between ij Is it satisfied? If satisfied, res i and res j The two residual point clusters n and cluster m merged into a residual point class; where κ is the norm of the clustering algorithm, if L is used 0 norm, then κ=0. If L 1 norm, then κ=1; α is the approximation factor, which is used to control the accuracy of the envelope sparsity theorem; 3d), calculation of the total polarity of the residual point class: for each residual point class cluster n If cluster n If the polarity of the internal residual point is 0, the cluster n The clustering radius of the inner residual points is set to β i =r n , where r n For cluster n The envelope distance refers to the maximum distance between two residual points in the residual point class; if cluster n If the polarity of the internal residual point is not 0, add cluster n The clustering radius of the inner residual points is β i , β i Should be greater than r n The value of 3e), determine whether the clustering iteration is finished: if there is a residual point class with unbalanced polarity, go directly to 3b) and re-cluster; otherwise, determine again whether any two residual points belonging to different residual point classes are equal. i and res j The distance between ij : If exists If , go to 3b) and re-cluster; if for all the residual points res belonging to different residual point classes i and res j The distance between ij ,That The clustering is completed and the clustering result is obtained.

8. The refined fast large-scale phase unwrapping method compatible with multiple optimization models according to claim 7 is characterized in that: S500 includes: 5a) For the area outside the convex hull, phase unwrapping is performed using the path integral method; 5b) For the area in the convex hull, mark the high-quality area in the convex hull according to the pseudo-coherence coefficient value and the given threshold, and discard the solution of the low-quality area in the convex hull to reduce the amount of calculation; 5c) For the high-quality area in the convex hull, according to the accuracy and efficiency requirements of the solution, the phase unwrapping algorithm is selected to solve the fuzzy number k of each pixel in the convex hull. r (s); 5d) Using the fuzzy number k of each pixel in the convex hull r (s) and winding phase Solve the absolute phase ψ of the pixels in the convex hull r (s), expressed as:

9. The refined fast large-scale phase unwrapping method compatible with multiple optimization models according to claim 8, characterized in that: The phase unwrapping algorithm includes L P norm algorithm, then L P The algorithm can be expressed as in, and is the weight at pixel s, which can be determined by the pseudo-coherence coefficient value of the input winding phase image; set W is the pixel set of each convex hull area to be solved, s, t and u are all pixel points in set W, where the horizontal adjacent pixel of pixel s is t, and the vertical adjacent pixel of pixel s is u; k r (s), k r (t) and k r (u) is the fuzzy number of pixels s, t, and u, representing the decision variable; is the fuzzy number gradient of the pixel, which is an integer; r = 1, 2, ..., M, represents the serial number of the winding phase map to be solved; when p takes different integer values, it represents different L P Norm algorithm.

10. The refined fast large-scale phase unwrapping method compatible with multiple optimization models according to claim 9, characterized in that: 5c) Includes: 5c1) Before solving the area inside the convex hull, expand the area inside the convex hull outward by one pixel; 5c2) When solving the area within the convex hull, the pixels on the boundary of the convex hull are used as the equality constraints of the algorithm, and then the phase unwrapping algorithm is used to solve it.

Images