Refined fast large-scale phase-unwrapping method compatible with multiple optimization models
By generating a mask using pseudocoherence coefficients and performing residual point clustering using the envelope sparsity theorem, combined with path integrals and optimization models, the computational resource requirements for large-scale InSAR data are addressed, achieving efficient and precise phase unwrapping, which is suitable for interferometric synthetic aperture radar data processing.
Patent Information
- Application Number
- CN202411317557.8
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2024-09-20
- Publication Date
- 2025-10-24
- Estimated Expiration
- 2044-09-20
AI Technical Summary
Existing single-baseline and multi-baseline InSAR phase unwrapping methods require significantly more computational resources when processing large-scale data, leading to insufficient memory and extended computation time. Furthermore, traditional methods lose valuable information when discarding the region within the convex hull during the solution process.
A mask is generated using pseudocoherence coefficients. Residual point clustering analysis is performed using the envelope sparsity theorem to divide high-quality and low-quality regions. A suitable optimization model is selected to refine the solution for the high-quality region inside the convex hull. The path integral method is combined to quickly untangle the region outside the convex hull. Finally, the results are spliced together to obtain the final result.
It improves clustering speed, reduces computational resource consumption, ensures high-quality region solution accuracy, and allows flexible selection of optimization models to handle different terrain features, achieving efficient and precise phase unwrapping.
Smart Images

Figure CN119335531B_ABST
Abstract
Description
Technical Field
[0001] The present invention belongs to the field of radar imaging technology, 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 an InSAR system is a wrapped phase image, in which each pixel is represented by a wrapped phase in the range (-π, π]. The wrapped phase is the principal value of the absolute phase that contains elevation and / or deformation information. Phase unwrapping is the step in the InSAR processing chain that derives the absolute phase from the wrapped phase. To date, two main types of InSAR techniques have been developed: single-baseline and multi-baseline InSAR.
[0003] The two aforementioned InSAR techniques employ different principles for phase unwrapping. Single-baseline phase unwrapping relies on a single wrapped phase image to recover the absolute phase. It employs the phase continuity assumption, assuming that the absolute phase gradient between 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 exploits the relationships between multiple baselines to determine the absolute phase. Since the diversity of baselines can significantly increase the ambiguity range of the wrapped phase during 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 acquired by InSAR system increases dramatically, and the demand for large-scale data processing is increasing. In some applications, very large and continuous spatial coverage is required, which may require processing several times larger than the standard wrapped phase map for phase unwrapping processing. For the traditional single baseline and multi baseline phase unwrapping method, when the input data exceeds a certain limit, two problems may occur: (1) insufficient memory; (2) significantly prolonged calculation time; this is because the calculation complexity of the traditional phase unwrapping method is high, and the calculation resources used in large-scale phase unwrapping are nonlinearly increased. For example, the minimum cost flow algorithm is a widely used phase unwrapping method, but its calculation complexity is highly nonlinear. At the same time, since multi baseline phase unwrapping needs to process multiple wrapped phase maps jointly, these two problems are more serious in multi baseline phase unwrapping. Under the existing technical conditions, even if the performance of the hardware is improved, such as increasing the number and speed of CPUs and increasing the amount of memory, it cannot meet the rapidly growing computing needs of large-scale phase unwrapping.
[0005] The basic idea of reducing the computing resources of large-scale data phase unwrapping is divide and conquer, that is, dividing the large-size wrapped phase map into sub-wrapped phase maps (hereinafter referred to as sub-maps) for solving, and then splicing the phase unwrapping results of each sub-map to obtain the overall phase unwrapping result. The sub-map division method is divided into regular division and irregular division, and the current more effective is the irregular sub-map division method based on residual point clustering analysis, combined with the large-scale phase unwrapping method of irregular sub-block solving. In order to improve the solving efficiency, the convex hull technology can be further used to generate a mask (the area inside the mask is not solved). Although this strategy of discarding the solving of the convex hull area can significantly improve the operation speed and reduce the memory consumption, since the area inside the convex hull is not an information-free area, this discarding will lose effective information. Therefore, in order to avoid the loss of effective information in the area inside the convex hull, a more refined solving method needs to be developed. SUMMARY
[0006] In order to solve the above problems existing in the prior art, the application provides a refined fast large-scale phase unwrapping method compatible with multiple optimization models (EFLS). The technical problems to be solved by the application are solved by the following technical scheme:
[0007] A refined fast large-scale phase unwrapping method compatible with multiple optimization models comprises:
[0008] S100, pseudo-coherence coefficients are calculated for a wrapped phase map to be solved, and a mask is generated for a connected domain with a pseudo-coherence coefficient lower than a specified threshold; the wrapped phase map includes a single baseline wrapped phase map and a multi baseline wrapped phase map;
[0009] S200, according to different cases of the wrapping phase map, gradient estimation is performed using a corresponding strategy, and then residual points are calculated to obtain a residual point map;
[0010] S300, residual points in the residual point map are analyzed by clustering using the envelope sparsity theorem and the mask, to obtain a clustering result of the residual points; wherein each unconnected mask corresponds to residual points covered in the residual point map as an initial residual point class;
[0011] S400, a convex hull of the residual point class is formed according to the clustering result, and a convex hull region is marked in the residual point map according to the convex hull point;
[0012] S500, a method for solving path integration is used to perform phase unwrapping on a region outside the convex hull in the wrapping phase map; according to the accuracy and efficiency requirements, the region inside the convex hull is divided into high-quality and low-quality regions by using a predetermined quality map and a specified threshold, and a suitable phase unwrapping algorithm is selected to perform phase unwrapping on the high-quality region inside the convex hull;
[0013] S600, the phase unwrapping results of the region inside the convex hull and the region outside the convex hull are spliced to obtain the final phase unwrapping result.
[0014] Advantages:
[0015] (1) After the mask is generated by using the pseudo-coherence coefficient, the residual points in each connected domain covered by the mask are used as the initial residual point class, so that the clustering time required for forming the initial residual point class from the residual points can be saved, and the clustering speed is improved; and then the residual points in the residual point map and the initial residual point class formed by the mask are analyzed by clustering according to the envelope sparsity theorem, so that the error diffusion of the low-quality region in the convex hull to the high-quality region is prevented, and the solving accuracy is ensured.
[0016] (2) The EFLS method proposed in the present application firstly solves the high-quality region outside the convex hull with extremely high efficiency, so that the user can obtain reliable phase unwrapping information of the high-quality region in a very short time; and then, the region inside the convex hull is solved. Therefore, this is a phase unwrapping algorithm with both front-end fast viewing and back-end fine processing capabilities, and the algorithm has more practical value.
[0017] (3) The region inside the convex hull is more finely divided into high-quality and low-quality regions in the present application. Due to the randomness of the residual point distribution, the region inside the convex hull not only contains the low-quality region, but also contains part of the high-quality region; therefore, the convex hull region can be finely divided according to the pseudo-coherence coefficient, and the solving of the low-quality region inside the convex hull is abandoned to reduce the amount of calculation, while the high-quality region inside the convex hull can still be effectively solved.
[0018] (4) The present application can flexibly select a suitable optimization model to solve the single baseline and multi-baseline phase unwrapping of the area within the convex hull, for example, ODLS, L 0 norm, L 1 norm, L ∞ norm and MIN, and can correctly perform image stitching to ensure that the optimal solution and the global optimal solution of the area within the convex hull are consistent. The present application will be further described in detail below in combination with the drawings and embodiments. BRIEF DESCRIPTION OF DRAWINGS
[0019] Figure 1 A flowchart of the fine and fast large-scale phase unwrapping method compatible with multiple optimization models provided by the present application;
[0020] Figure 2 An auxiliary residual point example diagram used when explaining the residual point theory of the present application;
[0021] Figure 3 Google Earth optical images and multi-baseline wrapped phase images of actual measurement data;
[0022] Figure 4 Comparison results of the clustering method of EFLS and the traditional clustering method;
[0023] Figure 5 Result diagram of the comparison of the solution of the high-quality area within the convex hull of EFLS and the solution of the convex hull of EFLS;
[0024] Figure 6 Phase unwrapping results and difference diagrams of the solution of the convex hull area by using multiple optimization models.
[0025] The present application will be further described in detail below in combination with specific embodiments, but the embodiments of the present application are not limited thereto.
[0026] Before introducing the present application, the technical concept of the present application is first introduced.
[0027] The existing large-scale phase unwrapping method mainly adopts the strategy of "divide and conquer", that is, the large-size wrapped phase image is divided into sub-images for separate solution, and then spliced. The traditional mask generation method based on the quality map cannot avoid the diffusion of the solution errors in the mask area to the high-quality area; the block method based on clustering analysis and convex hull technology selects the strategy of discarding the solution of the area within the convex hull to improve the solution speed, thereby losing the utilization rate of the information within the convex hull. In addition, the existing technology adopts a single model in a single processing, and cannot separately process different local topographies and image features in large-scale data. Therefore, the present application needs to design a block solution method that meets the following conditions:
[0028] (1) The wrapped phase map is divided into sub-maps for solving respectively, so that the solving results of the sub-maps can be correctly spliced, and the phenomenon that the solving error of a low-quality area in a sub-map diffuses to a high-quality area is avoided;
[0029] (2) On the basis of the existing block solving method based on clustering analysis and convex hull technology, in order to improve information utilization, the high-quality area in the convex hull is also effectively solved;
[0030] (3) When the high-quality area in the convex hull is solved, the additional resource consumption caused by the fine division of the high-quality area and the low-quality area is avoided, in other words, the ratio of the calculation resources consumed by the solving of the high-quality area in the convex hull to the calculation resources consumed by the overall solving of the convex hull is at least less than the ratio of the area of the solved area in the convex hull to the area of the convex hull;
[0031] (4) For the area in the convex hull, a plurality of optimization models can be flexibly selected, and the solving result of the area in the convex hull can be seamlessly spliced with the surrounding areas.
[0032] Therefore, the present application takes the residual point, envelope sparsity theorem and the like as the theoretical basis, adopts the research idea of combining the theoretical method with simulation, and proposes an EFLS method which is compatible with multiple optimization models and is fine and fast, so that the high-quality area in large-scale interference data is efficiently phase unwrapped.
[0033] The idea for realizing the present application is as follows: for the single baseline wrapped phase map and the multi-baseline wrapped phase map, the residual point maps corresponding to the single baseline wrapped phase map and the multi-baseline wrapped phase map are respectively solved by using the phase continuity assumption and the Chinese remainder theorem; before clustering the residual points, a mask is generated by using a quality map (for example, a pseudo-coherence coefficient) and a specified threshold, and the residual points covered by the mask are taken as the initial residual point classes for clustering, so as to improve the clustering speed; in the clustering process, the residual points and the initial residual point classes formed by the mask are analyzed according to the envelope sparsity theorem, so that the local optimal solution and the global optimal solution are consistent (or approximately consistent), and the error diffusion from the low-quality area to the high-quality area is prevented; the clustering area is marked by using an efficient convex hull, and the area outside the convex hull is quickly integrated and solved; according to the user demand, a suitable optimization model is flexibly selected to solve the area in the convex hull; in order to improve the information utilization, the high-quality area and the low-quality area in the convex hull are more finely divided, and the solving of the low-quality area in the convex hull is abandoned, while the high-quality area can still be effectively solved.
[0034] The technical scheme details of the present application are described in detail below.
[0035] Reference Figure 1 The present application provides a fine and fast large-scale phase unwrapping method compatible with multiple optimization models, which comprises the following steps:
[0036] S100, calculating pseudo-coherence coefficients for the wrapped phase map to be solved, and generating a mask for connected domains with pseudo-coherence coefficients lower than a specified threshold; the wrapped phase map includes a single baseline wrapped phase map and a multi-baseline wrapped phase map;
[0037] S200, performing gradient estimation using a corresponding strategy according to different cases of the wrapped phase map, and then calculating residual points to obtain a residual point map;
[0038] S300, performing cluster analysis on residual points in the residual point map by using an envelope sparsity theorem and the mask to obtain a cluster result of the residual points; wherein residual points covered by each unconnected mask in the residual point map are taken as initial residual point classes;
[0039] S400, forming a convex hull of a residual point class according to the cluster result, and marking a convex hull region in the residual point map according to convex hull points;
[0040] S500, performing phase unwrapping on a region outside the convex hull in the wrapped phase map by using a method for solving path integrals, and then performing phase unwrapping on a high-quality region in the convex hull by using a predetermined quality map and a specified threshold to divide the region into high-quality and low-quality regions according to accuracy and efficiency requirements, and selecting a suitable phase unwrapping algorithm to perform phase unwrapping on the high-quality region in the convex hull;
[0041] The predetermined quality map is formed by using pseudo-coherence coefficients.
[0042] S600, splicing phase unwrapping results of the region inside the convex hull and the region outside the convex hull to obtain a final phase unwrapping result.
[0043] In an optional embodiment of the present application, S200 includes:
[0044] S210, for a single baseline phase unwrapping problem, estimating a blur number gradient of the single baseline wrapped phase map by using a phase continuity assumption;
[0045] The present application estimates a blur number gradient of a single baseline wrapped phase map as follows:
[0046] S211, determining, according to a phase continuity assumption, that an absolute phase gradient between adjacent pixels of the single baseline wrapped phase map is less than π;
[0047] S212, estimating a blur number gradient of the single baseline wrapped phase map by using the phase continuity assumption based on the requirement that the absolute phase gradient is less than π, and the estimation is expressed as:
[0048]
[0049] wherein, is a wrapped phase of an arbitrary 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 horizontal 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 sequence number of the interferogram and baseline, r = 1, 2, …, M, M is the number of baselines in the multi-baseline InSAR system, and 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 in the present invention is as follows:
[0052] S221, determining the winding phase map of a multi-baseline InSAR system (taking a dual baseline as an example, i.e., M=2), the relationship between pixels at corresponding positions in the winding phase map 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 middle, 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, 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 pixel s, which are expressed as
[0056]
[0057]
[0058] Subtracting formula (4) and (5) from formula (3) respectively, we can get the following derivation formula:
[0059]
[0060]
[0061] wherein, s, t and u all represent pixels, t is a horizontally adjacent pixel of pixel s, and u is a vertically adjacent pixel of pixel s; k r (s), k r (t) and k r (u) are ambiguity numbers at these pixel positions; r = 1, 2, …, M, for a dual-baseline InSAR system, M = 2, so r = 1, 2;
[0062] S223, estimating ambiguity number gradients of the multi-baseline wrapped phase map by using the Chinese remainder theorem; the optimization model for solving the ambiguity number gradients of the multi-baseline wrapped phase map adopted by the Chinese remainder theorem is as follows:
[0063]
[0064] wherein,
[0065]
[0066] and are decision variables, which are integers; for some special combinations of (B1, B2), within a certain range, a unique
[0067] S230, calculating residual points of the wrapped phase map by using the ambiguity number gradients, and forming a residual point map;
[0068] This step calculates residual points of the wrapped phase map by using the ambiguity number gradients, assuming that the arrangement of four adjacent pixels is as shown in Figure 2 The calculation formula of the residual points is as follows:
[0069]
[0070] wherein, N represents the polarity of the residual point, for a single-baseline residual point, the value of N is -1, 0, 1; for a multi-baseline residual point, the value of N is an integer; and all represent ambiguity number gradients of corresponding pixel positions.
[0071] In an alternative embodiment of the present application, S300 comprises:
[0072] taking the residual points in each connected domain covered by the mask in S100 as initial residual point classes, and performing clustering analysis on the residual points in the residual point map and the initial residual point classes formed by the mask according to the envelope sparsity theorem to obtain clustering results of the residual points, comprising:
[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 in the residual point graph that are covered by a connected mask area are regarded as an initial residual point class; those residual points that are 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 d 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 d ij Is it satisfied If satisfied, res i and res j Belongs to the two residual point clusters n and cluster m Merge into a residual point class; where κ is the norm of the clustering algorithm, if L 0 norm, then κ=0, if L is used 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 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, then add clustern Cluster radius β of internal residual points 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 d ij : If exists If , go to 3b) and re-cluster; if for all residual points res belonging to different residual point classes i and res j The distance d ij ,That The clustering is completed and the clustering result is obtained.
[0078] In an optional embodiment of the present invention, S500 includes:
[0079] 5a) For the region outside the convex hull, phase unwrapping is performed using the path integral method.
[0080] 5b) For the area within the convex hull, mark the high-quality area within the convex hull according to the pseudo-coherence coefficient value and a given threshold, and discard the solution of the low-quality area within the convex hull to reduce the amount of computation;
[0081] 5c) For the high-quality area in the convex hull, according to the accuracy and efficiency requirements of the solution, select the appropriate 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 The weight at the pixel s can be determined by the pseudo-coherence coefficient value of the input wrapped phase map; the set W is a pixel set of each to-be-solved convex hull region, s, t and u are all pixels in the set W, wherein the horizontally adjacent pixel of the pixel s is t, and the vertically adjacent pixel of the pixel s is u; k r (s), k r (t) and k r (u) are the blur numbers of the pixels s, t and u, and represent decision variables; The blur number gradient of the pixel is an integer; r = 1, 2,..., M, represents the serial number of the to-be-solved wrapped phase map; when p takes different integer values, different L P norms can be adopted, for example, L 0 norm, L 1 norm, L ∞ norm, etc.
[0087] In an alternative embodiment of the present application, 5c) comprises:
[0088] 5c1) Before solving the in-convex-hull region, the in-convex-hull region is expanded by one pixel outward as a whole;
[0089] 5c2) When solving the in-convex-hull region, the pixels on the convex hull boundary are taken as the equality constraint condition of the algorithm, and then the phase unwrapping algorithm is solved. Since the out-of-convex-hull region is preferentially phase unwrapped, the absolute phase value of the pixel on the convex hull boundary is known; for the case of L P norm (p > 2), in order to make the values on the boundary between the in-convex-hull region and the out-of-convex-hull region equal, a weighted method can be used for the constraint condition on the boundary.
[0090] The solving result of the in-convex-hull region of the present application can be directly used for image stitching with the out-of-convex-hull region, without additional calculation. Since the out-of-convex-hull region is preferentially phase unwrapped in the present application, in order to ensure that the pixels in and out of the convex hull can be correctly stitched, when solving the in-convex-hull high-quality region, the pixels on the convex hull boundary are taken as the equality constraint condition of the algorithm, and then the solving is performed.
[0091] The device implementation part of the present application is consistent with the details 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 dataset to be solved is as follows: Figure 3 As shown in Figures (b), (c), and (d), they are multi-baseline twisted phase images with baseline lengths of -406.2m, -957.72m, and -1862.3m, respectively. Comparing the twisted phase image with the optical image, it can be seen that the area covered by the lake surface forms a clear 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 this 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. Verification of clustering speed improvement:
[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 traditional solution methods that do not divide the high / low quality regions, the computational resources required by the EFLS method are significantly reduced.
[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 is verified from the above three aspects in turn.
[0104] Experiment 1:
[0105] In this experiment, the clustering time of the traditional clustering method and the EFLS method of the application is compared to verify that the residual point clustering method of the EFLS method improves the clustering speed.
[0106] As Figure 4 Fig. (a) in the middle shows Figure 3 Fig. (b) in the middle is a pseudo-coherence coefficient map, and the threshold is set to 0.7. The area with a pseudo-coherence coefficient less than 0.7 is used as a mask, and the generated mask map is shown in Figure 4 Fig. (b) in the middle, 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 Fig. (c) in the middle is a residual point map obtained by jointly solving three wrapped phase maps (as shown in Figure 3 Figs. (b)-(d) in the middle). By comparing the mask map and the residual point map, it can be seen that the area with dense residual points is basically overlapped with the mask area. Therefore, in this experiment, the residual points covered by the mask in Figure 4 Fig. (b) in the middle form an initial residual point class. In Figure 4 Fig. (c) in the middle, the convex hull of the initial residual point class is marked with a black solid line (for the sake of clear image display, 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 takes each residual point as an initial residual point class, while the EFLS method of the application forms an initial residual point class by using the pseudo-correlation coefficient 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 number of clustering iterations of the traditional clustering method is 1735, and the time consumption is 1846.79 seconds. The number of clustering iterations of the EFLS method is 1593, and the time consumption is 1197.05 seconds, and the clustering time consumption is reduced to 64.8% of the traditional method. Of course, the improvement of the clustering speed is at the expense of the clustering precision. Figure 4 Fig. (d) in the middle is a convex hull map of the clustering result of the EFLS method, Figure 4The convex hulls in Figs. (d) and (e) are marked in the figures. The clustering results of the two methods are basically the same, so the EFLS method effectively improves the clustering efficiency without affecting the clustering results, and the superiority of the method will be more obvious when dealing with real data with a large range of noise. Figure 4 The convex hulls in Figs. (d) and (e) are marked in the figures. The clustering results of the two methods are basically the same, so the EFLS method effectively improves the clustering efficiency without affecting the clustering results, and the superiority of the method will be more obvious when dealing with real data with a large range of noise.
[0111] Experiment Two:
[0112] Based on the convex hull generation results of Experiment One, Experiment Two will verify that the EFLS method has fine solving capability and reasonable consumption of computing resources.
[0113] For easy comparison and testing, this experiment only tests the lake and its surrounding area, as shown in Fig. (a). Figure 5 As shown in Fig. (a), the largest convex hull and the convex hulls connected to it are marked in the pseudo-correlation coefficient figure with a white line frame. The area surrounded by these convex hulls is the area to be tested in this experiment, which is referred to as the lake area in the following text. Obviously, in the lake area, there is also a high-quality image area, which is the blue area within the white line frame 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 experiments compare this area.
[0114] Figure 5 Fig. (b) is Figure 3 The reference absolute phase of Fig. (b) is used to compare the solving accuracy. Figure 5 Fig. (c) is the L 1 The unwrapping result of the lake area solved by the L Figure 5 Fig. (d) is the difference between Figure 5 Fig. (b) and Fig. (c). Figure 5 Fig. (e) is the L 1The norm solution result, the low quality area is not solved and the value is 100 to observe. Figure 5 Fig. (f) is a diagram of Figure 5 Fig. (b) is a difference diagram of Fig. (e), wherein the difference of the unsolved area is 0. Figure 5 The color bar of Fig. (d) and Fig. (f) is 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, and after calculation, Figure 5 The mean square error of the high quality area in the lake of Fig. (d) and Fig. (f) is 0.55 and 0.52 respectively, which shows that the solution accuracy of the method proposed in the application is ensured. Table 3 compares the calculation resources consumed by the refined solution.
[0115] Table 3 Comparison table of calculation resource consumption of EFLS method refined solution
[0116]
[0117] As can be seen from the data in Table 3, after the improvement of the refined solution method, the solution time is reduced to 27.62% of the convex hull full solution, and the peak memory is also reduced to 23.19% of the original. The decrease of the calculation resource consumption is because the area of the solved area is reduced, but it should be noted that the implementation of the refined solution needs to consume additional calculation resources, and the refined solution method proposed in the application does not increase too much calculation resources, because the number of pixels to be solved is 38.73% of the number of pixels in the lake area, that is, the decrease of the resource consumption is greater than the decrease of the solved area. Therefore, this experiment verifies that the refined solution is effective in improving the solution accuracy and speed.
[0118] Experiment three:
[0119] The observed area used in this experiment has various terrains, in addition to the lake area concerned in the above experiment, such as Figure 3 As shown in Fig. (a), the upper left corner area is a complex mountain terrain, and the results obtained by using different phase unwrapping optimization models will have differences, so in the actual data processing, different optimization models may be used in different areas to obtain higher processing accuracy. Therefore, in this experiment, the L 1 norm and MIN two optimization models are used for solution, to verify the flexibility of the L P norm model selection in the convex hull area.
[0120] As shown in Fig. (a), the upper left corner area is a complex mountain terrain, and the results obtained by using different phase unwrapping optimization models will have differences, so in the actual data processing, different optimization models may be used in different areas to obtain higher processing accuracy. Therefore, in this experiment, the L 1 norm and MIN two optimization models are used for solution, to verify the flexibility of the L P norm model selection in the convex hull area. Figure 6The mask image (a) is shown in the figure. 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 more serious and more 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 rest of the convex hulls are solved by 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 plot, the mean squared errors in the mountainous region are 24.97 and 24.46, respectively. The two methods take 3697.87 seconds and 6897.24 seconds to solve, respectively, with peak memory usage of 263.11Mb and 408.62Mb. Different optimization models yield different accuracy and consume different computational resources. Therefore, the EFLS phase unwrapping algorithm proposed in this paper allows users to flexibly select the appropriate optimization model for local solution in different regions based on application requirements, without affecting the correct image stitching of different convex hull regions.
[0121] From the above three experiments, it can be seen that when performing residual point clustering, the EFLS method proposed in the present invention first uses the pseudo-coherence coefficient to generate a mask, 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 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 within 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] Furthermore, the terms "first", "second", etc. are used only for descriptive purposes and do not connote or imply relative importance. Thus, a feature defined with "first", "second", etc. can implicitly or explicitly include one or more of the features.
[0124] Although the present application has been described in connection with various embodiments thereof, it will be understood that other variations and modifications of the described embodiments can be made by those skilled in the art upon reading the description of the application set forth above. In the claims, the word "comprising" does not exclude other components or steps not mentioned in the claims. The word "a" or "an" preceding an element does not exclude the presence of a plurality of such elements.
[0125] The above description is further to make further detailed description of the present application in connection with specific preferred embodiments, and cannot be deemed to limit the specific implementation of the present application to these descriptions. For those skilled in the art, without departing from the concept of the present application, a number of simple deductions or substitutions can be made, which should be regarded as falling within the protection scope of the present application.
Claims
1. A multi-optimization model compatible refined fast large-scale phase unwrapping method, characterized in that, The method comprises the following steps: S100, pseudo-coherence coefficients are calculated for the wrapped phase map to be solved, and a mask is generated for a connected domain with a pseudo-coherence coefficient lower than a specified threshold value; The wrapped phase map comprises a single-baseline wrapped phase map and a multi-baseline wrapped phase map; S200, gradient estimation is performed using a corresponding strategy according to different cases of the wrapped phase map, and then residual points are calculated to obtain a residual point map; S300, clustering analysis is performed on the residual points in the residual point map by using an envelope sparsity theorem and the mask to obtain a clustering result of the residual points; wherein each unconnected mask corresponds to residual points covered in the residual point map as an initial residual point class; S400, a convex hull of the residual point class is formed according to the clustering result, and a convex hull region is marked in the residual point map according to the convex hull points; S500, a phase unwrapping method for solving path integration is used to perform phase unwrapping on a region outside the convex hull in the wrapped phase map; and according to accuracy and efficiency requirements of the solving, a predetermined quality map and a specified threshold value are used to divide a region inside the convex hull into a high-quality region and a low-quality region, and a phase unwrapping algorithm is selected to perform phase unwrapping on the high-quality region inside the convex hull; S600, the phase unwrapping results of the region inside the convex hull and the region outside the convex hull are spliced to obtain a final phase unwrapping result. S200 comprises: S210, for a single-baseline phase unwrapping problem, a blur number gradient of the single-baseline wrapped phase map is estimated by using a phase continuity assumption; S220, for a multi-baseline phase unwrapping problem, a blur number gradient of the multi-baseline wrapped phase map is estimated by using a Chinese remainder theorem; S230, residual points of the wrapped phase map are calculated by using the blur number gradient, and a residual point map is composed.
2. The fine-scale fast large-scale phase unwrapping method compatible with multi-optimization models according to claim 1, wherein, S210 comprises: S211, according to the phase continuity assumption, it is determined that an absolute phase gradient between adjacent pixels of the single-baseline wrapped phase map is less than π; S212, based on the requirement that the absolute phase gradient is less than π, the blur number gradient of the single-baseline wrapped phase map is estimated by using the phase continuity assumption, and is expressed as: where, is the wrapping phase of the pixel s, is the wrapping phase of the neighboring pixel t right to the pixel s, is the wrapping phase of the neighboring pixel u below the pixel s; is the horizontal blur number gradient between the pixel s and the neighboring column of the pixel t, is the vertical blur number gradient between the pixel s and the neighboring row of the pixel u; r represents the serial number of the wrapping phase map and the baseline, r = 1, 2, …, M, M is the number of baselines of the multi-baseline InSAR system, for single baseline phase unwrapping M = 1, therefore r = 1.
3. The fine-scale fast large-scale phase unwrapping method compatible with multi-optimization models of claim 2, wherein, S220 comprises: S221, a relationship between pixels at corresponding positions of wrapped phase maps of different baselines in a wrapped phase map of a multi-baseline InSAR system is determined, and is expressed as wherein B1 and B2 are the lengths of the two perpendicular baselines corresponding to wrapped phase map 1 and wrapped phase map 2 respectively obtained by the dual-baseline InSAR system, is the wrapped phase at pixel s in wrapped phase map 1, and k1(s) is the ambiguity number at pixel s in wrapped phase map 1, is the wrapped phase at pixel s in wrapped phase map 2, and k2(s) is the ambiguity number at pixel s in wrapped phase map 2. S222, relationships between pixels at pixels s transversely adjacent to the pixel s and pixels u longitudinally adjacent to the pixel s are determined, and are respectively expressed as The deduction formula is obtained by subtracting formula (4) and formula (5) from formula (3) respectively as follows: wherein, s, t and u all represent pixels, t is a horizontally adjacent pixel of pixel s, and u is a vertically adjacent pixel of pixel s; k r (s), k r (t) and k r (u) are the blur numbers at these pixel positions; r = 1, 2, …, M, for a dual-baseline InSAR system, M = 2, thus r = 1, 2; S223, an optimization model of the Chinese remainder theorem is used to solve the blur number gradient of the multi-baseline wrapped phase map, and the optimization model is: wherein, and are decision variables, which are integers; for certain special combinations of (B1, B2), within a certain range, a unique solution can be found 4. The fine-scale fast large-scale phase unwrapping method compatible with multi-optimization models of claim 3, wherein, S230 comprises: The residual points of the wrapped phase map are calculated by using the blur number gradient to compose a residual point map, and a calculation formula of the residual points is expressed as: Wherein, N represents the polarity of the residual point, for single baseline residual point, N takes the value of -1, 0, 1; for multi-baseline residual point, N takes the value of integer; and Both represent the blur number gradient of the corresponding pixel position.
5. The fine-scale fast large-scale phase unwrapping method compatible with multi-optimization models according to claim 4, wherein, S300 comprises: The residual points in each connected domain covered by the mask in S100 are taken as initial residual point classes, and clustering analysis is performed on the residual points in the residual point map and the initial residual point classes formed by the mask according to an envelope sparsity theorem to obtain a clustering result of the residual points.
6. The fine-scale fast large-scale phase unwrapping method compatible with multi-optimization models according to claim 5, wherein, S300 comprises: 3a), setting initial residual point classes: for each residual point res i , set its clustering radius β i as the initial clustering radius β i0 ; wherein, the residual points covered by one connected mask region in the residual point map are regarded as one initial residual point class; those residual points not covered by the connected mask region, each residual point is regarded as an initial residual point class; 3b), calculate the distance between residual points: take any two different residual points res i and res j from cluster n and cluster m respectively, and calculate the distance d ij between res i and res j , where res i ∈ cluster n , res j ∈ cluster m ; 3c), determine whether to merge the residual point classes: determine whether the distance d i between res j and res ij satisfies If yes, merge the two residual point classes cluster i and cluster j to which res n and res m belong; where κ is the norm selection of the clustering algorithm, κ = 0 if L 0 norm is used, and κ = 1 if L 1 norm is used; and α is an approximation factor to control the accuracy of the envelope sparsity theorem. 3d), total polarity calculation of the cluster of residual points: for each cluster of residual points n if the polarity sum of the residual points in cluster n is 0, set the cluster radius of the residual points in cluster n to β i = r n , where r n is the envelope distance of cluster n , the envelope distance refers to the maximum distance between two residual points in the cluster of residual points; if the polarity sum of the residual points in cluster n is not 0, increase the cluster radius β n of the residual points in cluster i , β i should take a value greater than r n . 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 d ij : If exists If , go to 3b) and re-cluster; if for all residual points res belonging to different residual point classes i and res j The distance d ij ,That The clustering is completed and the clustering result is obtained.
7. The fine-scale fast large-scale phase unwrapping method compatible with multi-optimization models according to claim 6, wherein, S500 comprises: 5a) for a region outside the convex hull, a phase unwrapping method for solving path integration is used to perform phase unwrapping thereon; 5b) For the area inside the convex hull, mark the high quality area inside the convex hull according to the pseudo-coherence coefficient value and a given threshold value, and discard the solution for the low quality area inside the convex hull to reduce the amount of calculation; 5c) For high quality area within the convex hull, select a phase unwrapping algorithm according to the accuracy and efficiency requirement of the solution, get the ambiguity number k of each pixel within the convex hull r (s); 5d) Using the blur number k for each pixel inside the convex hull r (s) and the wrapped phase Solving for the absolute phase ψ of the pixels inside the convex hull r (s) is expressed as:
8. The multi-optimization model compatible refined fast large-scale phase unwrapping method according to claim 7, wherein, The phase unwrapping algorithm comprises an algorithm based on L P norm, then L P The algorithm can be expressed as wherein, and is the weight at pixel s, which can be determined by the pseudo-coherence coefficient value of the input wrapped phase map; W is the pixel set of each to-be-solved convex hull region, s, t and u are all pixel points in the set W, wherein the horizontally adjacent pixel of pixel s is t, and the vertically adjacent pixel of pixel s is u; k r s), k r t) and k r u) are the blur numbers of pixels s, t and u, representing the decision variables; are the blur number gradients of the pixels, which are integers; r = 1, 2, …, M, representing the serial number of the to-be-solved wrapped phase map; when p takes different integer values, it represents different L P norm algorithms.
9. The multi-optimization model compatible refined fast large-scale phase unwrapping method of claim 8, wherein, 5c) comprising: 5c1) before solving the area inside the convex hull, expand the area inside the convex hull as a whole by one pixel outward; 5c2) when solving the area inside the convex hull, the pixels of the convex hull boundary are taken as the equation constraint condition of the algorithm, and then the phase unwrapping algorithm is solved.
Citation Information
Patent Citations
Regional division based multi-strategy InSAR (Interferometric synthetic aperture) radar phase unwrapping method
CN104316922A
InSAR (Interferometric Synthetic Aperture Radar) interferometric phase two-step unwrapping method combining quality map and minimum cost flow
CN113311433A