Multi-baseline INSAR phase unwrapping method, system, device and medium

Abstract

The invention discloses a multi-baseline INSAR phase unwrapping method, system and device and a medium, and the method comprises the steps: constructing a to-be-clustered data set: obtaining the intercept information corresponding to each pixel according to the interferograms corresponding to different vertical baselines, then taking the position information of each pixel as the multi-dimensional clustering feature of the pixel, and obtaining the to-be-clustered data set; taking the pixel as a to-be-clustered target to obtain a to-be-clustered data set; wavelet clustering processing: performing clustering processing on the data set to be clustered based on wavelet clustering; clustering result correction: aiming at pixels in each noise cluster, correcting cluster labels of the pixels; phase unwrapping is carried out cluster by cluster: on the basis of the corrected cluster distribution, a fuzzy vector of each cluster is calculated by using a closed solution formula or a sparse-TSPA method, and the absolute phase of each pixel in the corresponding cluster is calculated on the basis of the fuzzy vector of each cluster. When the method is used for processing a large-size interference pattern, the method has more advantages in the aspects of efficiency, precision and base line ratio adaptability compared with an existing multi-base line phase unwrapping method.

Description

Technical Field

[0001] The present invention relates to the field of remote sensing technology, and in particular to a multi-baseline INSAR phase unwrapping method, system, equipment and medium. Background Art

[0002] Interferometric Synthetic Aperture Radar (InSAR) is an advanced remote sensing technology used to obtain digital elevation models (DEMs) and detect surface deformation. Phase unwrapping is a key step in InSAR data processing. Multi-baseline phase unwrapping overcomes the phase continuity assumption of traditional single-baseline phase unwrapping methods by increasing the ambiguity interval of the interferometric phase, thereby achieving higher accuracy and robustness in complex terrain and noisy environments.

[0003] Although multi-baseline phase unwrapping technology has made significant progress in recent years, existing methods still face many challenges when processing large-scale interferograms, mainly in the following aspects:

[0004] 1. The balance between accuracy and efficiency.

[0005] Existing multi-baseline phase unwrapping methods struggle to balance accuracy and efficiency when processing large-scale interferograms. For example, methods based on parameterization, Kalman filtering, and two-stage spatial planning (TSPA), while demonstrating good noise robustness in certain scenarios, suffer from low efficiency when processing large-scale data, making them difficult to meet the demands of practical applications.

[0006] 2. Limitation of baseline ratio.

[0007] Current state-of-the-art TSPA methods perform poorly when the baseline ratio is less than 2. According to the baseline design criterion for nonlinear mixed integer programming models (NIP criterion), the optimal baseline ratio should meet specific conditions. However, when the baseline ratio is less than 2, a search window that meets the conditions cannot be designed, resulting in the accumulation of phase unwrapping errors, which limits its scope in practical applications.

[0008] 3. Shortcomings of existing clustering algorithms.

[0009] While multi-baseline phase unwrapping methods based on cluster analysis offer unique advantages in efficiency and accuracy, existing clustering algorithms still suffer from unstable performance. For example, the DBSCAN (spatial clustering) algorithm employed by CANOPUS (an improved multi-baseline phase unwrapping algorithm based on cluster analysis) is highly sensitive to user-defined parameters and can easily generate a large number of small clusters, leading to inaccurate clustering results. Furthermore, CANOPUS does not address noise clusters, which can be a significant source of phase unwrapping errors.

[0010] 4. Dimension mismatch problem.

[0011] In multi-baseline phase unwrapping methods based on cluster analysis, there is a dimensionality mismatch between the intercept feature and the position features (row and column). The intercept feature's value range is related to the baseline ratio and measured phase, while the position feature's value range is related to the specific coordinate position of the pixel. This dimensionality mismatch makes it difficult for clustering algorithms to accurately distinguish different clusters, which in turn affects the accuracy of phase unwrapping.

[0012] 5. Insufficient noise robustness.

[0013] Existing multi-baseline phase unwrapping methods exhibit limitations when dealing with noise. For example, parameterized methods are sensitive to noise, and even small amounts of noise can cause noticeable glitches in the unwrapped results. While deep learning-based methods exhibit promising performance in some situations, their generalization and interpretability are limited, and research in this area remains limited.

[0014] In summary, existing technologies for processing large-scale interferograms face challenges such as a difficult balance between accuracy and efficiency, baseline ratio limitations, unstable clustering algorithm performance, dimensionality mismatch, and insufficient noise robustness, limiting their application in complex terrain and large-scale data processing. Therefore, a new approach is needed to overcome these limitations and improve the accuracy, efficiency, and adaptability of multi-baseline phase unwrapping to baseline ratios. Summary of the Invention

[0015] In response to the deficiencies in the prior art, the present invention provides a multi-baseline INSAR phase unwrapping method, system, device and medium, introduces wavelet clustering to achieve more accurate clustering results, and can efficiently process large-scale interferograms, thereby improving the phase unwrapping accuracy.

[0016] In a first aspect, a multi-baseline INSAR phase unwrapping method is provided, comprising the following steps:

[0017] S1: Constructing the dataset to be clustered: According to the interference patterns corresponding to different vertical baselines, the intercept information corresponding to each pixel is obtained. Then, the position information of each pixel is combined as the multidimensional clustering feature of the pixel, and the pixel is used as the target to be clustered to obtain the dataset to be clustered;

[0018] S2: Wavelet clustering processing: clustering the data set based on wavelet clustering;

[0019] S3: Clustering result correction: for each pixel in the noise cluster, correct its cluster label;

[0020] S4: Cluster-by-cluster phase unwrapping: Based on the corrected cluster distribution, the fuzzy vector of each cluster is calculated using a closed-form solution formula or the sparse-TSPA method. Based on the fuzzy vector of each cluster, the absolute phase of each pixel in the corresponding cluster is calculated.

[0021] According to the first aspect, in some possible implementations, when there are M interference patterns corresponding to different vertical baselines, one of them is selected as the main image and is respectively compared with the other interference patterns. Figure 1 A combination of calculations obtains M-1 groups of intercept information, which are combined with the row and column of each pixel to form M+1-dimensional clustering features; where M≥2.

[0022] According to the first aspect, in some possible implementations, step S2 specifically includes:

[0023] Initialize the grid size and neighborhood radius;

[0024] Based on the initial grid size, each dimension in the multidimensional feature space is divided into multiple intervals, thereby dividing the multidimensional feature space into several grids; then the number of pixel objects contained in each grid is counted;

[0025] Applying wavelet transform to the quantized multidimensional feature space to obtain a new transform space, where each grid in the new transform space has a corresponding wavelet coefficient;

[0026] Filter out the significant grids whose wavelet coefficients are greater than or equal to the threshold τ in the new transformed space, and then complete the clustering by finding k-connected significant grids;

[0027] The cluster labels corresponding to each grid in the new transformed space are mapped to the grid in the original multidimensional feature space, and then each grid cluster label in the multidimensional feature space is assigned to all pixel targets in the grid.

[0028] According to the first aspect, in some possible implementations, when initializing the grid size and neighborhood radius, the following conditions must be met:

[0029] Assume that the grid size of the qth intercept dimension is Δb q express;

[0030] but:

[0031]

[0032] Where, represents rounding down, ε represents the neighborhood radius, represents the minimum distance between adjacent cluster centerlines in the qth intercept dimension; therefore, Δb q The maximum value should not exceed At the same time, when the winding phase When the intercept b∈(-1,B q / B1), B1, B q are two vertical baselines respectively; thus:

[0033]

[0034] Where, Indicates the number of intervals into which the qth intercept dimension is divided; Indicates rounding up;

[0035] Assuming that the grid sizes in the row and column dimensions are represented by Δr and Δl respectively, when it is a dual-baseline INSAR system, the number of intervals m2 and m3 divided in the row and column dimensions meet the following conditions:

[0036]

[0037] Where n 2l and n 2s are the number of fringes in the row dimension of the long baseline and short baseline interferograms, n 3l and n 3s are the number of fringes in the column dimension of the long baseline and short baseline interferograms, respectively; N2 and N3 are the number of rows and columns of the interferogram, respectively;

[0038] When it is a multi-baseline INSAR system, the number of intervals m2 and m3 in the row and column dimensions must meet the following conditions:

[0039]

[0040]

[0041] Where n2 max and n3 max are the maximum number of fringes in the row dimension and the maximum number of fringes in the column dimension among all baseline interferograms, respectively.

[0042] According to the first aspect, in some possible implementations, the initial value of the threshold τ for screening significant grids is determined by the following method:

[0043] Screen out all grids with non-zero wavelet coefficients in the new transformed space;

[0044] Determine the proportion of non-significant grids in all grids with non-zero wavelet coefficients based on the estimated noise grid proportion;

[0045] The initial value of the threshold τ is obtained based on the non-salient grid ratio.

[0046] According to the first aspect, in some possible implementations, when applying wavelet transform to the quantized multidimensional feature space, a stationary wavelet transform is adopted.

[0047] According to the first aspect, in some possible implementations, step S3 specifically includes:

[0048] Noise clusters and pseudo clusters with the number of pixel targets within the cluster below the threshold are uniformly regarded as noise clusters;

[0049] For each pixel in the noise cluster, the cluster label with the highest frequency within the preset window centered on the noise cluster is selected as the corrected cluster label of the pixel.

[0050] Secondly, a multi-baseline INSAR phase unwrapping system is provided, comprising:

[0051] The module for constructing the dataset to be clustered is used to obtain the intercept information corresponding to each pixel based on the interference patterns corresponding to different vertical baselines, and then combine the position information of each pixel as the multidimensional clustering feature of the pixel, and use the pixel as the target to be clustered to obtain the dataset to be clustered;

[0052] Wavelet clustering processing module, used to perform clustering processing on the clustering data set based on wavelet clustering;

[0053] The clustering result correction module is used to correct the cluster label of each pixel in the noise cluster;

[0054] The cluster-by-cluster phase unwrapping module is used to calculate the blur vector of each cluster based on the corrected cluster distribution using a closed-form solution formula or the sparse-TSPA method, and calculate the absolute phase of each pixel in the corresponding cluster based on the blur vector of each cluster.

[0055] According to a third aspect, an electronic device is provided, including:

[0056] Memory on which computer programs or instructions are stored;

[0057] A processor is configured to load and execute the computer program or instructions to implement the multi-baseline INSAR phase unwrapping method as described above.

[0058] In a fourth aspect, a readable storage medium is provided, on which a computer program or instruction is stored. When the computer program or instruction is executed by a processor, the multi-baseline INSAR phase unwrapping method as described above is implemented.

[0059] The present invention proposes a multi-baseline INSAR phase unwrapping method, system, device and medium, which uses a wavelet clustering algorithm to cluster pixels with multidimensional clustering features. Through multi-scale analysis of wavelet transform, the problem of incorrect clustering caused by dimensionality mismatch in existing clustering methods (such as DBSCAN) is avoided, thereby improving clustering efficiency and accuracy; according to the unique information of the InSAR data set, appropriate initial grid and neighborhood parameters are selected to solve the dimensionality mismatch problem between intercept features and position features. By optimizing these parameters, the number of incorrect clusters caused by dimensionality mismatch is reduced; cluster correction is introduced during post-processing, which can effectively identify and correct the noise clusters identified by wavelet clustering, further improving the accuracy of phase unwrapping. Therefore, when processing large-scale interferograms, the present invention has more advantages than existing multi-baseline phase unwrapping methods in terms of efficiency, accuracy and adaptability of baseline ratio. BRIEF DESCRIPTION OF THE DRAWINGS

[0060] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0061] Figure 1 This is a flow chart of a multi-baseline INSAR phase unwrapping method provided by an embodiment of the present invention;

[0062] Figure 2 : is a comparison chart of the experimental results of the first experiment provided by the present invention; wherein, (a) is the reference DEM, (b) is the short baseline simulation interferogram, (c) is the long baseline simulation interferogram, (d) is the intercept obtained by formula (4), (e) is the cluster distribution obtained by CANOPUS, (f) is the cluster distribution obtained by WCRFPU, (g) is the elevation inversion result of (c) by CANOPUS, (h) is the elevation inversion result of (c) by WCRFPU, (i) is the difference between (g) and (a), and (j) is the difference between (h) and (a);

[0063] Figure 3: is a comparison diagram of the experimental results of the second experiment provided by the present invention; wherein, (a) is the reference DEM, (b) is the short baseline simulation interferogram, (c) is the long baseline simulation interferogram, (d) is the intercept diagram obtained by formula (4), (e) is the cluster distribution obtained by the CANOPUS algorithm, (f) is the cluster distribution obtained by the WCRFPU algorithm, (g) is the result of the elevation inversion of (c) by the TSPA algorithm, (h) is the result of the elevation inversion of (c) by the CANet-PU algorithm, (i) is the result of the elevation inversion of (c) by the CANOPUS algorithm, (j) is the result of the elevation inversion of (c) by the WCRFPU algorithm, (k) is the difference between (g) and (a), (l) is the difference between (h) and (a), (m) is the difference between (i) and (a), and (n) is the difference between (j) and (a);

[0064] Figure 4 : is a comparison diagram of the experimental results of the third experiment provided by the present invention; wherein, (a) is a Google Earth image of the study area, (b) is the reference phase unwrapping result of (d) obtained by SRTM, (c) is the short baseline interferogram of TanDEM-X, (d) is the long baseline interferogram of TanDEM-X, (e) is the phase unwrapping result of (d) obtained by CANOPUS, (f) is the phase unwrapping result of (d) obtained by CANet-PU, (g) is the phase unwrapping result of (d) obtained by TSPA, (h) is the phase unwrapping result of (d) obtained by WCRFPU, (i) is the difference between (e) and (b), (j) is the difference between (f) and (b), (k) is the difference between (g) and (b), and (l) is the difference between (h) and (b). DETAILED DESCRIPTION

[0065] To make the objectives, technical solutions, and advantages of the present invention more apparent, the technical solutions of the present invention will be described in detail below. Obviously, the embodiments described are only some of the embodiments of the present invention, not all of them. Based on the embodiments of the present invention, all other implementations obtained by those of ordinary skill in the art without inventive effort are within the scope of protection of the present invention.

[0066] like Figure 1 As shown, the embodiment of the present invention provides a multi-baseline INSAR phase unwrapping method, including the following steps:

[0067] S1: Construct the data set to be clustered: According to the interference patterns corresponding to different vertical baselines, obtain the intercept information corresponding to each pixel and amplify it. Then, combine the position information of each pixel as the multidimensional clustering feature of the pixel, and use the pixel as the target to be clustered to obtain the data set to be clustered.

[0068] In some embodiments, when there are M interference patterns corresponding to different vertical baselines, one of them is selected as the main image and compared with the other interference patterns. Figure 1 A combination of calculations obtains M-1 groups of intercept information, which are combined with the row and column of each pixel to form M+1-dimensional clustering features; where M≥2.

[0069] For ease of understanding, the following is an example of a dual-baseline InSAR system. The viewing angle, wavelength, and vertical baseline of the system are represented by the parameters θ, λ, and B, respectively. i (i=1,2). Therefore, the terrain height h(s) of the s-th pixel can be expressed by this equation:

[0070]

[0071] Among them, r(s) represents the slant distance of the s-th pixel, ψ i (s) is the absolute phase of the s-th pixel in the i-th interferogram. Since the InSAR system can only measure the wrapping phase, the wrapping phase corresponding to the s-th pixel in the i-th interferogram is It is expressed as follows:

[0072]

[0073] Among them, k i (s) is the fuzzy number corresponding to the sth pixel in the i-th interference pattern. According to equations (1) and (2), the elevations corresponding to pixels at the same position in the interference patterns measured by different baselines are the same, so:

[0074]

[0075] Taking the two unknowns k1(s) and k2(s) in the above formula as the independent variable and the dependent variable respectively, formula (3) can be expressed as follows:

[0076]

[0077] According to formula (4), the fuzzy numbers k1(s) and k2(s) corresponding to the same pixel in different interference patterns are linearly related, and the two together constitute the fuzzy vector, This is the intercept of the sth pixel. This intercept and the pixel position information (including row and column) can form a three-dimensional cluster feature. In addition, pixels with the same intercept correspond to the same blur vector, otherwise Equation (4) will be contradictory. Therefore, by grouping pixels with the same intercept into a cluster, each pixel in the cluster has the same blur vector. Multi-baseline phase unwrapping is the process of obtaining the blur vector corresponding to each pixel.

[0078] The above analysis is based on a dual-baseline INSAR system. If it is a multi-baseline INSAR system, the above formula (3) can be rewritten as:

[0079]

[0080] Select the interference pattern corresponding to the vertical baseline B1 as the main image, and the interference pattern corresponding to other vertical baselines Figure 1 One combination can calculate M-1 groups of intercept information, which can be combined with the row and column of each pixel to form M+1 dimensional clustering features.

[0081] S2: Wavelet clustering processing: clustering the data set based on wavelet clustering.

[0082] First, some basic concepts of wavelet clustering are introduced:

[0083] (1) Empty grid: In the quantitative feature space, a grid with a count value of 0 is called an empty grid;

[0084] (2) Non-empty grid: In the quantitative feature space, the grid with non-zero count value is called a non-empty grid;

[0085] (3) Significant grid: In the transformed space, the grid with wavelet coefficients greater than a certain threshold τ is called a significant grid;

[0086] (4) ε-neighborhood: Assume there are two salient grids (in the transformed space) or non-empty grids (in the quantized feature space) c1 and c2. If the Euclidean distance between c1 and c2 satisfies d(c1,c2)≤ε, then the two grids are each other's ε-neighborhoods, where ε is the radius of the sphere used to define the neighborhood;

[0087] (5) k-ε-neighborhood: If c1 and c2 are both salient grids (in the transformed space) or both non-empty grids (in the quantized feature space), and c1 is one of the k preset ε-neighborhoods of c2, then the grid c1 is called the k-ε-neighborhood of the grid c2;

[0088] (6) k-connected: If there exists a set of meshes p1, p2, ..., p j , so that p1=c1 and p j =c2, and p i+1 It is p i If there is a k-ε-neighborhood (1≤i≤j-1) of c1 and c2, then the grids c1 and c2 are said to be k-connected;

[0089] (7) Cluster: If the significant grid set {c1,c2,…,c m} is k-connected in the transformed space, then the set is considered to be a cluster.

[0090] From the above definitions, we can see that the k-connectivity between salient meshes in the transformed space is defined by the k-ε-neighborhood. If there are two salient meshes c1 and c2 that are k-connected, then we can imagine that there is a connected path between them that allows them to reach each other in the transformed space. These definitions can be used to describe the relationship between meshes in the feature space, thus helping us understand their interactions.

[0091] The following is a detailed description of the wavelet clustering process, which includes:

[0092] S2.1: Initialize the grid size and neighborhood radius.

[0093] Based on the unique information of the InSAR dataset, appropriate initial grid size and neighborhood parameters should be selected. The selection of the initial grid size should consider the feature ranges in the intercept, row, and column dimensions to address the dimensionality mismatch between intercept and position features. Neighborhood parameters (such as the neighborhood radius ε) should be selected to ensure reasonable grid connectivity in the quantized space and avoid clustering errors caused by overly large or small neighborhoods.

[0094] The optimal grid size needs to be found in the grid partitioning scheme. Initial values can be set based on the InSAR dataset information, and then grid parameters can be adjusted based on the clustering results. Wavelet clustering does not require the precise number of clusters to be known in advance, but an estimate of the expected number of clusters helps select an appropriate grid size. A grid that is too small or too large can negatively impact clustering results. When the grid is too small, some pixels that should belong to the same cluster may be assigned to different clusters, resulting in over-quantization. In this case, the number of bins selected for each dimension can be gradually reduced. As this process continues, the grid size gradually increases. When the grid is too large, some pixels from different clusters may be assigned to the same cluster, resulting in under-quantization. In this case, the number of bins selected for each dimension can be gradually increased, resulting in a smaller grid size. Using a heuristic approach, we can visually observe the cluster distribution and whether the distribution and boundaries of each cluster correspond to the intercept distribution. This allows us to determine whether the boundaries between clusters are reasonable. This process is repeated until an appropriate cluster distribution is achieved. Thanks to the multi-resolution nature of the wavelet transform, clustering results at different scales can be observed with the same grid size, reducing the need for grid parameter adjustments. Based on the above analysis, the choice of initial grid size is very important. Based on the proprietary information of the InSAR dataset, the grid size can be initialized in three dimensions: intercept, row, and column.

[0095] The wavelet clustering algorithm is a grid-based clustering algorithm, which requires the feature space to be gridded, and then the grid is used as the basic clustering unit to judge the connectivity of the grid in the quantitative feature space or the transformation space. Taking the dual-baseline InSAR system as an example, it is assumed that the grid sizes in the intercept, row, and column dimensions are represented by Δb, Δr, and Δl, respectively. The connectivity of the grid is determined in the quantitative feature space, and a sphere with a radius of ε is defined to define the neighborhood. It must be ensured that all pixels in all grids in the neighborhood belong to the same cluster. This requires that all pixels in any grid in the neighborhood must belong to the same cluster, and the difference in intercept between pixels in all adjacent grids and pixels in the central grid must not exceed half of the minimum distance between the center lines of adjacent clusters in the intercept dimension, that is, d min / 2. Therefore, the following conditions must be met:

[0096]

[0097] Where, Indicates rounding down. Therefore, the division of intercept dimension should not be too rough during quantization. q The maximum value should not exceed At the same time, it can be seen from formula (4) that when the winding phase When , the intercept b∈(-1,B2 / B1), B1 and B2 are two vertical baselines respectively; therefore:

[0098]

[0099] Where m1 represents the number of intervals into which the intercept dimension is divided; Indicates rounding up.

[0100] On the other hand, the grid sizes Δr and Δl in the row and column dimensions are related to the number of fringes in the long-baseline and short-baseline interferograms. The basic requirement is that the number of intervals m2 and m3 divided in each dimension should not be less than the sum of the maximum number of fringes in the two interferograms in that dimension. Assuming the size of the interferogram is N2 rows and N3 columns, then:

[0101]

[0102] Where n 2l and n 2s are the number of fringes in the row dimension of the long baseline and short baseline interferograms, n 3l and n 3s are the number of fringes in the column dimension of the long baseline and short baseline interferograms, respectively.

[0103] The above is a description of grid parameter initialization for a dual-baseline INSAR system. For a multi-baseline INSAR system, the following conditions must be met when initializing the grid size and neighborhood radius:

[0104] Assume that the grid size of the qth intercept dimension is Δb q express;

[0105] but:

[0106]

[0107] Where, represents rounding down, ε represents the neighborhood radius, represents the minimum distance between adjacent cluster centerlines in the qth intercept dimension; therefore, Δb q The maximum value should not exceed At the same time, when the winding phase When the intercept b∈(-1,B q / B1), B1, B q There are two vertical baselines, B q ∈B, where B is the set of all vertical baselines except the vertical baseline B1 used as the reference; therefore:

[0108]

[0109] Where, Indicates the number of intervals into which the qth intercept dimension is divided; Indicates rounding up;

[0110] Assume that the grid sizes in the row and column dimensions are represented by Δr and Δl respectively, and the number of intervals m2 and m3 divided in the row and column dimensions meet the following conditions:

[0111]

[0112] Where n2 max and n3 max are the maximum number of fringes in the row dimension and the maximum number of fringes in the column dimension after counting all baseline interferograms.

[0113] S2.2: Based on the initial grid size, each dimension in the multidimensional feature space is divided into multiple intervals, thereby dividing the multidimensional feature space into several grids; then counting the number of pixel objects contained in each grid.

[0114] S2.3: Apply wavelet transform to the quantized multidimensional feature space to obtain a new transform space, where each grid in the new transform space has a corresponding wavelet coefficient.

[0115] In this embodiment, a stationary wavelet transform (SWT) is applied to the quantized multidimensional feature space to obtain a new transform space. It should be noted that the method used here is not the discrete wavelet transform (DWT) commonly used in traditional wavelet transforms. This is because although DWT is more efficient, it is not suitable for the high-precision requirements of the PU field. Therefore, this embodiment uses SWT at the expense of slightly reducing efficiency in exchange for a significant improvement in accuracy. Thereafter, each grid in the transform space has a corresponding value, called a wavelet coefficient, so that a grid with wavelet coefficient information can be obtained.

[0116] S2.4: Screen out the significant grids whose wavelet coefficients are greater than or equal to the threshold τ in the new transformed space, and then complete the clustering by finding k-connected significant grids.

[0117] Significant grids with wavelet coefficients greater than or equal to a threshold τ in the new transformed space are selected. The clustering process is to find a set of k-connected significant grids. The threshold τ for significant grids is crucial, as it determines the proportion of noise grids among all non-empty grids. If the threshold is set too low, some noise grids may be mistakenly merged into a cluster; conversely, if the threshold is set too high, some normal grids may be mistakenly identified as noise grids. An initial threshold can be estimated by setting an initial non-significant grid ratio for all grids with non-zero wavelet coefficients based on the estimated proportion of noise grids, based on the distribution of wavelet coefficients. A heuristic method is then used to find the optimal threshold. By gradually adjusting the threshold and observing the cluster boundaries, this helps determine the optimal threshold. Once the threshold τ is determined, grids in the transformed space with wavelet coefficients greater than or equal to the threshold τ are marked as significant, while all other grids are marked as noise grids.

[0118] S2.5: Map the cluster labels corresponding to each grid in the new transformed space to the grid in the original multidimensional feature space, and then assign each grid cluster label in the multidimensional feature space to all pixel targets in the grid.

[0119] In the new transformed space, each grid is assigned a label for the corresponding cluster based on the identified cluster distribution. Since clusters are identified in the transformed space, a lookup table should be created to map them to grids in the original multidimensional feature space. Afterwards, the cluster label of each grid in the multidimensional feature space is assigned to all pixel objects in that grid. Ultimately, we can obtain a 2D cluster distribution of multidimensional clustered feature pixel objects at multiple resolutions. If it is found that this 2D cluster distribution does not look appropriate (judged by whether the distribution of each cluster, the cluster boundaries, etc. correspond to the intercept distribution), the grid size and neighborhood parameters can be adjusted until the user is satisfied with the result.

[0120] S3: Clustering result correction: For each pixel in the noise cluster, correct its cluster label.

[0121] The accuracy of the phase unwrapping results is affected by the noise clusters identified by the wavelet clustering algorithm, which require correction. Furthermore, some small clusters (those with a target number of pixels below a threshold) are artifacts caused by phase noise and can therefore also be considered noise clusters. The basic principle behind this correction is to select, for each pixel in a noise cluster, the cluster label with the highest frequency within a preset window centered around it as the corrected cluster label for that pixel.

[0122] S4: Cluster-by-cluster phase unwrapping: Based on the corrected cluster distribution, the fuzzy vector of each cluster is calculated using a closed-form solution formula or the sparse-TSPA method, which effectively improves efficiency compared to the search-based method in CANOPUS. Based on the fuzzy vector of each cluster, the absolute phase of each pixel in the corresponding cluster can be calculated according to Equation (2).

[0123] In order to ensure that the solution of the present invention is applicable to the processing of large-scale interference patterns, the clustering process should have a high execution efficiency. Assume that there are N (N is very large) pixels (targets) in the interference pattern and the clustering feature has d dimensions (d is very small). The wavelet clustering algorithm first identifies all the data and then divides them into corresponding grids; therefore, the time complexity of this step is O(N). For convenience, assume that each dimension is divided into m equal intervals and the original feature space is divided into K = m d grids. The second step is to apply SWT (Stationary Wavelet Transform) to these grids; its time complexity is O(ldK) = O(dK) = O(K), where l is the filter length used by SWT, which is usually small. If we apply SWT to decompose the number of layers T, and each layer downsamples the space by half, the time required will be:

[0124]

[0125] The above formula shows that the cost of applying the wavelet transform is at most O(4K / 3) = O(K). This shows that the wavelet clustering algorithm can achieve multi-resolution representation very efficiently. The time required to find connected components is O(cK) = O(K), where c is a small constant relative to K. The time required to build the lookup table is O(K). The final step requires assigning corresponding cluster labels to all targets, so the time consumption is O(N). Based on the above analysis, if the time consumption of reading the dataset and assigning cluster labels to targets is not taken into account, the time complexity of wavelet clustering is O(K). Therefore, the intermediate process is only related to the number of grids, not the size of the interference pattern. When the size of the interference pattern is very large, K is usually much smaller than N, and the efficiency advantage of wavelet clustering will be more obvious.

[0126] The above embodiment provides a multi-baseline INSAR phase unwrapping method, which uses a wavelet clustering algorithm to cluster pixels with multidimensional clustering features. Through multi-scale analysis of wavelet transform, it avoids the problem of erroneous clustering caused by dimensionality mismatch in existing clustering methods (such as DBSCAN), thereby improving clustering efficiency and accuracy. According to the unique information of the InSAR data set, the appropriate initial grid and neighborhood parameters are selected to solve the dimensionality mismatch problem between the intercept feature and the position feature. By optimizing these parameters, the number of erroneous clusters caused by dimensionality mismatch is reduced. Cluster correction is introduced during post-processing, which can effectively identify and correct the noise clusters identified by wavelet clustering, further improving the accuracy of phase unwrapping. Therefore, when processing large-scale interferograms, this method has advantages over existing multi-baseline phase unwrapping methods in terms of efficiency, accuracy, and adaptability of baseline ratio.

[0127] An embodiment of the present invention further provides a multi-baseline INSAR phase unwrapping system, comprising:

[0128] The module for constructing the dataset to be clustered is used to obtain the intercept information corresponding to each pixel based on the interference patterns corresponding to different vertical baselines, and then combine the position information of each pixel as the multidimensional clustering feature of the pixel, and use the pixel as the target to be clustered to obtain the dataset to be clustered;

[0129] Wavelet clustering processing module, used to perform clustering processing on the clustering data set based on wavelet clustering;

[0130] The clustering result correction module is used to correct the cluster label of each pixel in the noise cluster;

[0131] The cluster-by-cluster phase unwrapping module is used to calculate the blur vector of each cluster based on the corrected cluster distribution using a closed-form solution formula or the sparse-TSPA method, and calculate the absolute phase of each pixel in the corresponding cluster based on the blur vector of each cluster.

[0132] It should be understood that the functional unit modules in various embodiments of the present invention can be concentrated in one processing unit, or each unit module can exist physically separately, or two or more unit modules can be integrated into one unit module, and can be implemented in the form of hardware or software.

[0133] An embodiment of the present invention further provides an electronic device, including:

[0134] Memory on which computer programs or instructions are stored;

[0135] A processor is configured to load and execute the computer program or instructions to implement the multi-baseline INSAR phase unwrapping method as described above.

[0136] An embodiment of the present invention further provides a readable storage medium having a computer program or instructions stored thereon. When the computer program or instructions are executed by a processor, the multi-baseline INSAR phase unwrapping method as described above is implemented.

[0137] It can be understood that the same or similar parts of the above embodiments can be referenced to each other, and the contents not described in detail in some embodiments can refer to the same or similar contents in other embodiments.

[0138] The effects of the technical solution of the present invention will be further explained below in conjunction with specific experiments.

[0139] In this experiment, three datasets were used to validate the performance of WCRFPU (the solution of the present invention). The first two were simulated datasets, and the last was a real InSAR dataset. The first experiment verified the efficiency of WCRFPU on large-scale simulated interferograms. The second simulation experiment verified the accuracy of WCRFPU. The third experiment tested the practicality of WCRFPU using a real InSAR dataset. To better demonstrate the performance of WCRFPU, all experiments were compared in detail with other methods.

[0140] In the first experiment, a test was conducted in a simulated scene. The reason for choosing this simulated dataset is that the scene is hilly and undulating, which makes phase unwrapping challenging. In addition, due to the large size of the dataset (3000×3000 pixels), it can effectively verify the efficiency of the algorithm. The height range of the scene is set to 0 meters to 162 meters, and the blur height is set to 60 meters and 36 meters respectively. In order to test the noise robustness of WCRFPU, phase noise with coherence coefficients of 0.8 and 0.7 is added to the short baseline simulated interferogram and the long baseline simulated interferogram, respectively. Reference DEM is as follows Figure 2 As shown in (a), the simulated interference patterns corresponding to different baselines are shown in Figure 2 In (b) and (c). The intercept graph obtained according to formula (4) is as follows Figure 2 As shown in (d).

[0141] For CANOPUS, the radius is set to 4 and the minimum number of points in the neighborhood is set to 29. The final cluster distribution is as follows Figure 2 As shown in (e). For WCRFPU, we divide the three-dimensional feature space into 90×90×90 grids, and the percentage of non-significant grids of all grids with non-zero wavelet coefficients is set to a certain value. The final cluster distribution is shown in Figure 2 By comparing the cluster distributions obtained by the two algorithms, we can observe that the number of clusters identified by WCRFPU is significantly less than that of CANOPUS, and the cluster distribution of WCRFPU is closer to the actual distribution.

[0142] In order to present the final results more fairly and objectively, all algorithms use closed-formula solutions to obtain cluster ambiguity vectors. The phase unwrapping results and corresponding error distributions of the two algorithms are marked with the same color scale range. In order to comprehensively evaluate the phase unwrapping results of different algorithms, the root mean square error (RMSE) and phase unwrapping success rate (PUSR) are calculated. PUSR is defined as the ratio of the number of pixels that successfully recover the ambiguity to the total number of pixels in the long baseline interferogram. The elevation inversion results obtained using CANOPUS and WCRFPU are shown in Figure 2 The error distributions of CANOPUS and WCRFPU are shown in (g) and (h). Figure 2 Table 1 lists the time consumption of the two algorithms, as well as the RMSE and PUSR of the corresponding phase unwrapping results.

[0143] Table 1 Comparison of the two methods

[0144]

[0145] Since the number of grids divided in this experiment is much lower than the target number, and the number of identified clusters is also much lower than CANOPUS, WCRFPU has a significant advantage in efficiency.

[0146] The second experiment simulated the MB InSAR interferogram using a DEM (458×157 pixels) of a mountainous area. The reference DEM has an altitude range of 0–136.7 meters, as shown in Figure 2. Figure 3 As shown in (a). The fuzzy heights corresponding to the two different baselines are 73 meters and 43.8 meters respectively, and the baseline ratio is B2 / B1=5 / 3. In order to test the noise robustness of WCRFPU, phase noise with coherence coefficients of 0.85 and 0.8 was added to the short baseline interferogram and long baseline interferogram respectively. The simulated short baseline and long baseline interferograms are shown in Figure 1. Figure 3 As shown in (b) and (c). The intercept corresponding to each pixel can be calculated by formula (4), as Figure 3 As shown in (d). It can be observed that the intercept value is much smaller than the number of rows and columns. Next, we will compare TSPA, CANet-PU (a cluster analysis phase unwrapping method based on the introduction of deep learning) and CANOPUS to verify the performance of WCRFPU. For TSPA, we use a residual point-based algorithm to solve the minimum cost flow (MCF) problem. For CANOPUS, the radius is set to 2, and the minimum number of points in its neighborhood is set to 7. The resulting cluster distribution is shown in Figure 3 As shown in (e). For WCRFPU, the three-dimensional feature space is divided into 40×40×40 grids. For all grids with non-zero wavelet coefficients, the percentage of non-significant grids is set to 60%, and the final cluster distribution is as follows Figure 3The elevation inversion results using TSPA, CANet-PU, CANOPUS and WCRFPU are shown in (f). Figure 3 The elevation inversion errors of TSPA, CANet-PU, CANOPUS and WCRFPU are shown in (g)–(j). Figure 3 As shown in (k)–(n).

[0147] Table 2 lists the RMSE and PUSR corresponding to the phase unwrapping results of the four algorithms. Due to the small baseline ratio (less than 2), the performance of TSPA and CANet-PU is severely limited, resulting in the gradual accumulation of errors in the phase unwrapping process in space. CANOPUS and WCRFPU are able to overcome the baseline ratio limitation and can process interferograms with small baseline ratios. However, CANOPUS identifies too many clusters and tends to calculate incorrect blur vectors at cluster boundaries. Due to dimensionality mismatch, some false clusters are also generated, resulting in many errors in the phase unwrapping results. The WCRFPU algorithm significantly reduces errors by selecting appropriate grid and neighborhood parameters, matching the intercept dimension with the position dimension, and performing additional processing on noisy clusters. In terms of efficiency, due to the small interferogram size, all algorithms can quickly obtain the corresponding phase unwrapping results. However, the WCRFPU algorithm still has a slight advantage.

[0148] Table 2 Performance comparison of four algorithms

[0149]

[0150] The third experiment used a real dual-baseline InSAR dataset acquired from the TanDEM-X satellite in Weinan, Shaanxi Province. This dataset was chosen because the terrain is steep and does not satisfy the phase continuity assumption. Therefore, the SBPU algorithm is not suitable for this type of terrain, which poses a significant challenge to multi-baseline phase unwrapping techniques. Figure 4 (a) shows the long baseline reference phase unwrapping result obtained by SRTM DEM, from which a large number of phase jumps can be observed. Figure 4 (b) and (c) show the interferograms of short baseline and long baseline, respectively. Figure 4 Panels (e)-(h) show the long-baseline phase unwrapping results obtained by CANOPUS, CANet-PU, TSPA, and WCRFPU. For CANOPUS, the radius is set to 4, and the minimum number of points in its neighborhood is set to 35. For TSPA, the MCF problem is solved based on the estimation of the fuzzy number gradient. For WCRFPU, the 3D feature space is divided into a 70×70×70 grid. The percentage of non-significant grids among all grids with non-zero wavelet coefficients is set to 40%. Figure 4(i)-(l) show the error distributions of CANOPUS, CANet-PU, TSPA, and WCRFPU.

[0151] The main parameters of the TanDEM-X satellite system are listed in Table 3. Table 4 lists the time consumption of the four algorithms and the corresponding RMSE and PUSR of the phase unwrapping results. Since the Tandem-X system has a baseline ratio of approximately 2.8, it meets the baseline ratio requirements for TSPA and CANet. Therefore, TSPA and CANet did not experience the error accumulation phenomenon seen in the second experiment. Although TSPA achieved the best phase unwrapping results, it was the most time-consuming and unsuitable for processing large interferograms. CANet requires a large amount of data to train its network model. Insufficient or low-quality data may prevent the model from fully learning the inherent structure of the dataset, resulting in poor clustering performance. Furthermore, CANet requires longer training time and more computing resources. WCRFPU took the shortest time in this experiment. However, the high number of clusters (457 clusters) not only affects the efficiency of cluster-by-cluster PU but also inevitably leads to large errors at cluster boundaries. In contrast, WCRFPU reduces the adverse effects of dimensionality mismatch between cluster features on the clustering results. The 107 clusters identified by WaveCluster are closer to the true spatial distribution of the data. Therefore, WCRFPU exhibits excellent efficiency and accuracy when processing large-scale interferograms.

[0152] Table 3 Main interference parameters of the real data set

[0153]

[0154] Table 4 Performance comparison of four algorithms

[0155]

[0156] Those skilled in the art will appreciate that the embodiments of the present application can be provided as methods, systems, or computer program products. Therefore, the present application can adopt the form of a complete hardware embodiment, a complete software embodiment, or an embodiment in combination with software and hardware. Moreover, the present application can adopt the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to magnetic disk storage, CD-ROM, optical storage, etc.) that contain computer-usable program code.

[0157] The present application is described with reference to the flowcharts and / or block diagrams of the methods, devices (systems), and computer program products according to the embodiments of the present application. It should be understood that each process and / or box in the flowchart and / or block diagram, as well as the combination of the processes and / or boxes in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the steps in the process. Figure 1 a process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0158] These computer program instructions may also be stored in a computer readable memory that can direct a computer or other programmable data processing device to work in a specific manner, so that the instructions stored in the computer readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 a process or multiple processes and / or boxes Figure 1 The function specified in one or more boxes.

[0159] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operational steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing the instructions executed on the computer or other programmable device for implementing the process. Figure 1 a process or multiple processes and / or boxes Figure 1 A step that specifies a function in one or more boxes.

[0160] Although the embodiments of the present invention have been shown and described above, it will be understood that the above embodiments are illustrative and are not to be construed as limitations on the present invention. A person skilled in the art may change, modify, replace and modify the above embodiments within the scope of the present invention.

Claims

1. A multi-baseline INSAR phase unwrapping method, characterized in that: The steps include: S1: Constructing the dataset to be clustered: According to the interference patterns corresponding to different vertical baselines, the intercept information corresponding to each pixel is obtained. Then, the position information of each pixel is combined as the multidimensional clustering feature of the pixel, and the pixel is used as the target to be clustered to obtain the dataset to be clustered; S2: Wavelet clustering processing: clustering the data set based on wavelet clustering; S3: Clustering result correction: for each pixel in the noise cluster, correct its cluster label; S4: Cluster-by-cluster phase unwrapping: Based on the corrected cluster distribution, the fuzzy vector of each cluster is calculated using a closed-form solution formula or the sparse-TSPA method. Based on the fuzzy vector of each cluster, the absolute phase of each pixel in the corresponding cluster is calculated.

2. The multi-baseline INSAR phase unwrapping method according to claim 1, characterized in that: In step S1, when there are M interference patterns corresponding to different vertical baselines, one of them is selected as the main image, and is combined with other interference patterns one by one to calculate M-1 groups of intercept information, which is combined with the row and column of each pixel to form an M+1-dimensional clustering feature; Where M≥2.

3. The multi-baseline INSAR phase unwrapping method according to claim 1, characterized in that: Step S2 specifically includes: Initialize the grid size and neighborhood radius; Based on the initial grid size, each dimension in the multidimensional feature space is divided into multiple intervals, thereby dividing the multidimensional feature space into several grids; then the number of pixel objects contained in each grid is counted; Applying wavelet transform to the quantized multidimensional feature space to obtain a new transform space, where each grid in the new transform space has a corresponding wavelet coefficient; Filter out the significant grids whose wavelet coefficients are greater than or equal to the threshold τ in the new transformed space, and then complete the clustering by finding k-connected significant grids; The cluster labels corresponding to each grid in the new transformed space are mapped to the grid in the original multidimensional feature space, and then each grid cluster label in the multidimensional feature space is assigned to all pixel targets in the grid.

4. The multi-baseline INSAR phase unwrapping method according to claim 3, characterized in that: When initializing the grid size and neighborhood radius, the following conditions must be met: Assume that the grid size of the qth intercept dimension is Δb q express; but: Where, represents rounding down, ε represents the neighborhood radius, represents the minimum distance between adjacent cluster centerlines in the qth intercept dimension; therefore, Δb q The maximum value should not exceed At the same time, when the winding phase When the intercept b∈(-1,B q / B1), B1, B q are two vertical baselines respectively; thus: Where, Indicates the number of intervals into which the qth intercept dimension is divided; Indicates rounding up; Assuming that the grid sizes in the row and column dimensions are represented by Δr and Δl respectively, when it is a dual-baseline INSAR system, the number of intervals m2 and m3 divided in the row and column dimensions meet the following conditions: Where n 2l and n 2s are the number of fringes in the row dimension of the long baseline and short baseline interferograms, n 3l and n 3s are the number of fringes in the column dimension of the long baseline and short baseline interferograms, respectively; N2 and N3 are the number of rows and columns of the interferogram, respectively; When it is a multi-baseline INSAR system, the number of intervals m2 and m3 in the row and column dimensions must meet the following conditions: Where n2 max and n3 max are the maximum number of fringes in the row dimension and the maximum number of fringes in the column dimension among all baseline interferograms, respectively.

5. The multi-baseline INSAR phase unwrapping method according to claim 3, characterized in that: The initial value of the threshold τ for filtering significant grids is determined by the following method: Screen out all grids with non-zero wavelet coefficients in the new transformed space; Determine the proportion of non-significant grids in all grids with non-zero wavelet coefficients based on the estimated noise grid proportion; The initial value of the threshold τ is obtained based on the non-salient grid ratio.

6. The multi-baseline INSAR phase unwrapping method according to claim 3, characterized in that: When applying wavelet transform to the quantized multidimensional feature space, stationary wavelet transform is adopted.

7. The multi-baseline INSAR phase unwrapping method according to claim 1, characterized in that: Step S3 specifically includes: Noise clusters and pseudo clusters with the number of pixel targets within the cluster below the threshold are uniformly regarded as noise clusters; For each pixel in the noise cluster, the cluster label with the highest frequency within the preset window centered on the noise cluster is selected as the corrected cluster label of the pixel.

8. A multi-baseline INSAR phase unwrapping system, characterized by: include: The module for constructing the dataset to be clustered is used to obtain the intercept information corresponding to each pixel based on the interference patterns corresponding to different vertical baselines, and then combine the position information of each pixel as the multidimensional clustering feature of the pixel, and use the pixel as the target to be clustered to obtain the dataset to be clustered; Wavelet clustering processing module, used to perform clustering processing on the clustering data set based on wavelet clustering; The clustering result correction module is used to correct the cluster label of each pixel in the noise cluster; The cluster-by-cluster phase unwrapping module is used to calculate the blur vector of each cluster based on the corrected cluster distribution using a closed-form solution formula or the sparse-TSPA method, and calculate the absolute phase of each pixel in the corresponding cluster based on the blur vector of each cluster.

9. An electronic device, characterized in that: include: Memory on which computer programs or instructions are stored; A processor, configured to load and execute the computer program or instructions to implement the multi-baseline INSAR phase unwrapping method according to any one of claims 1 to 7.

10. A readable storage medium having a computer program or instruction stored thereon, characterized in that: When the computer program or instruction is executed by a processor, the multi-baseline INSAR phase unwrapping method according to any one of claims 1 to 7 is implemented.

Images