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

Abstract

The application discloses a multi-baseline INSAR phase unwrapping method, system, device and medium, wherein the method comprises the following steps: constructing a to-be-clustered data set: acquiring intercept information corresponding to each pixel according to interferograms corresponding to different vertical baselines, then combining position information of each pixel as multi-dimensional clustering features of the pixel, and taking the pixel as a to-be-clustered target to obtain the to-be-clustered data set; wavelet clustering processing: performing clustering processing on the to-be-clustered data set based on wavelet clustering; cluster result correction: correcting cluster labels of pixels in each noise cluster; cluster-by-cluster phase unwrapping: based on the corrected cluster distribution, calculating a blur vector of each cluster by using a closed solution formula or a sparse-TSPA method, and calculating absolute phases of each pixel in a corresponding cluster based on the blur vector of each cluster. When large-size interferograms are processed, the application has more advantages than existing multi-baseline phase unwrapping methods in terms of efficiency, accuracy and adaptability of a baseline ratio.

Description

Technical Field

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

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

[0003] Despite significant progress in multi-baseline phase unwrapping techniques in recent years, existing methods still face numerous challenges when processing large-scale interferograms, mainly in the following aspects:

[0004] 1. The issue of balancing accuracy and efficiency.

[0005] Existing multi-baseline phase unwrapping methods struggle to effectively balance accuracy and efficiency when processing large-scale interferograms. For instance, while parametric methods, Kalman filtering, and two-stage programming (TSPA) exhibit good noise robustness in certain scenarios, they are inefficient when processing large-scale data, making them unsuitable for practical applications.

[0006] 2. Limitations on baseline ratio.

[0007] The current state-of-the-art TSPA method performs poorly when the baseline ratio is less than 2. According to the baseline design criterion (NIP criterion) of the nonlinear mixed integer programming model, the optimal baseline ratio should meet specific conditions. However, when the baseline ratio is less than 2, it is impossible to design a search window that meets the conditions, leading to the accumulation of phase unwrapping errors and limiting 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 terms of efficiency and accuracy, existing clustering algorithms still suffer from performance instability. For example, the DBSCAN (spatial clustering) algorithm used by CANOPUS (an improved multi-baseline phase unwrapping algorithm based on cluster analysis) is highly sensitive to user-defined parameters, easily generating a large number of small clusters and leading to inaccurate clustering results. Furthermore, CANOPUS does not handle noisy clusters, which can be a significant source of phase unwrapping errors.

[0010] 4. Dimension mismatch issue.

[0011] In multi-baseline phase unwrapping methods based on clustering analysis, a dimensionality mismatch exists between intercept features and positional features (rows and columns). The range of intercept values ​​is related to the baseline ratio and the measured phase, while the range of positional features is related to the specific coordinates of the pixels. This dimensionality mismatch makes it difficult for clustering algorithms to accurately distinguish different clusters, thus affecting the accuracy of phase unwrapping.

[0012] 5. Insufficient noise robustness.

[0013] Existing multi-baseline phase unwrapping methods exhibit certain limitations when handling noise. For example, parameter-based methods are sensitive to noise, and even small amounts of noise can cause noticeable "glitch" in the unwrapping results. While deep learning-based methods show better performance in some cases, their generalization ability and interpretability are relatively weak, and research in this area remains limited.

[0014] In summary, existing technologies for processing large-size interferograms suffer from several drawbacks, including difficulty in balancing accuracy and efficiency, limitations in baseline ratio, unstable clustering algorithm performance, dimensionality mismatch, and insufficient robustness to noise. These limitations restrict their application in complex terrains and large-scale data processing. Therefore, a new method is needed to overcome these shortcomings and improve the accuracy, efficiency, and adaptability to baseline ratios of multi-baseline phase unwrapping. Summary of the Invention

[0015] To address the shortcomings of existing technologies, this invention provides a multi-baseline INSAR phase unwrapping method, system, device, and medium. It introduces wavelet clustering to achieve more accurate clustering results and can efficiently process large-size interferograms, thereby improving the phase unwrapping accuracy.

[0016] Firstly, a multi-baseline INSAR phase unwrapping method is provided, including the following steps:

[0017] S1: Construct the dataset to be clustered: Based on the interferograms corresponding to different vertical baselines, obtain the intercept information corresponding to each pixel, 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.

[0018] S2: Wavelet clustering processing: Clustering the dataset to be clustered based on wavelet clustering;

[0019] S3: Clustering result correction: Correct the cluster label for each pixel in the noisy cluster;

[0020] S4: Cluster-by-cluster phase unwrapping: Based on the corrected cluster distribution, calculate the blur vector of each cluster using the closed-form formula or the sparse-TSPA method, and calculate the absolute phase of each pixel within the corresponding cluster based on the blur vector of each cluster.

[0021] According to the first aspect, in some possible implementations, when there are M interferograms corresponding to different vertical baselines, one of them is selected as the master image, and compared with other interferograms... Figure 1 A combination of calculations yields M-1 sets of intercept information, which, together with the row and column of each pixel, constitutes an M+1-dimensional clustering feature; 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, thus dividing the multidimensional feature space into several grids; then the number of pixel targets contained in each grid is counted.

[0025] Wavelet transform is applied to the quantized multidimensional feature space to obtain a new transform space, in which each grid has a corresponding wavelet coefficient;

[0026] Significant grids with wavelet coefficients greater than or equal to the threshold τ in the new transform space are selected, and then clustering is completed by finding k-connected significant grids;

[0027] The cluster labels corresponding to each grid in the new transformation space are mapped to the grids in the original multidimensional feature space. Then, the cluster label of each grid in the multidimensional feature space is assigned to all pixel targets in that grid.

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

[0029] Assume the grid size for the q-th intercept dimension is Δb q express;

[0030] but:

[0031]

[0032] In the formula, This indicates rounding down, and ε represents the neighborhood radius. This represents the minimum distance between adjacent cluster centerlines in the q-th intercept dimension; therefore, Δb q The maximum value should not exceed At the same time, when the entangled phase When the intercept b∈(-1,B) q / B1), B1, B q These are two vertical baselines; therefore, we have:

[0033]

[0034] In the formula, This represents the number of intervals into which the q-th intercept dimension is divided; Indicates rounding up;

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

[0036]

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

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

[0039]

[0040]

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

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

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

[0044] The proportion of insignificant grids in all non-zero wavelet coefficients is determined based on the estimated proportion of noise grids.

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

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

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

[0048] Noise clusters and pseudo-clusters with fewer than a threshold number of pixels within a cluster are collectively considered as noise clusters;

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

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

[0051] The clustering dataset construction module is used to obtain the intercept information of each pixel based on the interferogram 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 clustering target to obtain the clustering dataset.

[0052] The wavelet clustering processing module is used to perform clustering processing on the dataset to be clustered based on wavelet clustering.

[0053] The clustering result correction module is used to correct the cluster labels of pixels in each noisy 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 the closed-formula solution or the sparse-TSPA method, and to calculate the absolute phase of each pixel within the corresponding cluster based on the blur vector of each cluster.

[0055] Thirdly, an electronic device is provided, comprising:

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

[0057] A processor for loading and executing the computer program or instructions to implement the multi-baseline INSAR phase unwrapping method as described above.

[0058] Fourthly, a readable storage medium is provided on which a computer program or instructions are stored, which, when executed by a processor, implement the multi-baseline INSAR phase unwrapping method as described above.

[0059] This invention proposes a multi-baseline InSAR phase unwrapping method, system, device, and medium. It employs a wavelet clustering algorithm to cluster pixels with multi-dimensional clustering features. Through multi-scale analysis using wavelet transform, it avoids the erroneous clustering problem caused by dimensionality mismatch in existing clustering methods (such as DBSCAN), thereby improving clustering efficiency and accuracy. Based on the unique information of the InSAR dataset, appropriate initial grid and neighborhood parameters are selected to address the dimensionality mismatch between intercept and location features. By optimizing these parameters, the number of erroneous clusters caused by dimensionality mismatch is reduced. Clustering correction is introduced during post-processing, which effectively identifies and corrects noise clusters identified by wavelet clustering, further improving the accuracy of phase unwrapping. Therefore, when processing large-size interferograms, this invention has advantages over existing multi-baseline phase unwrapping methods in terms of efficiency, accuracy, and baseline ratio adaptability. Attached Figure Description

[0060] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

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

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

[0063] Figure 3This is a comparison diagram of the experimental results of the second experiment provided by the present invention; where (a) is the reference DEM, (b) is the short baseline simulated interferogram, (c) is the long baseline simulated interferogram, (d) is the intercept diagram obtained by equation (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 elevation inversion of (c) by the TSPA algorithm, (h) is the result of elevation inversion of (c) by the CANet-PU algorithm, (i) is the result of elevation inversion of (c) by the CANOPUS algorithm, (j) is the result of 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 This is a comparison chart 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 a short baseline interferogram of TanDEM-X, (d) is a 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 Implementation

[0065] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention will be described in detail below. Obviously, the described embodiments are merely some embodiments of this invention, and not all embodiments. Based on the embodiments of this invention, all other implementation methods obtained by those skilled in the art without creative effort are within the scope of protection of this invention.

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

[0067] S1: Construct the dataset to be clustered: Based on the interferograms corresponding to different vertical baselines, obtain the intercept information corresponding to each pixel and magnify 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 dataset to be clustered.

[0068] In some embodiments, when there are M interferograms corresponding to different vertical baselines, one of them is selected as the master image and compared with other interferograms. Figure 1 A combination of calculations yields M-1 sets of intercept information, which, together with the row and column of each pixel, constitutes an M+1-dimensional clustering feature; where M≥2.

[0069] To facilitate understanding, we will first use a dual-baseline InSAR system as an example. The viewing angle, wavelength, and vertical baseline of this system are determined by parameters θ, λ, and B, respectively. i (i = 1, 2) represents the terrain height h(s) of the s-th pixel. Therefore, the terrain height h(s) of the s-th pixel can be represented by this equation:

[0070]

[0071] Where r(s) represents the slope distance of the s-th pixel, ψ i (s) represents the absolute phase of the s-th pixel in the i-th interferogram. Since the InSAR system can only measure the wrapped phase, the wrapped phase corresponding to the s-th pixel in the i-th interferogram... Expressed as follows:

[0072]

[0073] Where, k i (s) represents the ambiguity number corresponding to the s-th pixel in the i-th interferogram. According to equations (1) and (2), the elevations corresponding to pixels at the same position in interferograms measured from different baselines are the same, therefore:

[0074]

[0075] Taking the two unknowns k1(s) and k2(s) in the above equation as the independent and dependent variables respectively, equation (3) can be expressed in the following form:

[0076]

[0077] According to equation (4), the fuzzy numbers k1(s) and k2(s) corresponding to the same pixel in different interferograms are linearly related, and the two together constitute the fuzzy vector. This is the intercept of the s-th pixel, and this intercept, together with the pixel's position information (including row and column), can constitute a three-dimensional clustering feature. Furthermore, pixels with the same intercept correspond to the same blur vector; otherwise, equation (4) would be contradictory. Therefore, pixels with the same intercept can be grouped into a cluster, and each pixel within 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, then equation (3) above can be rewritten as:

[0079]

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

[0081] S2: Wavelet clustering: Clustering the dataset to be clustered based on wavelet clustering.

[0082] First, let's introduce some basic concepts of wavelet clustering:

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

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

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

[0086] (4) ε-neighborhood: Suppose there are two salient grids (in the transform 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 ε-neighborhoods of each other, 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 transform space) or both are non-empty grids (in the quantization feature space), and c1 is one of the k pre-set ε-neighborhoods of c2, then grid c1 is called a k-ε-neighborhood of grid c2.

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

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

[0090] As can be seen from the above definitions, the k-connectivity between salient grids in the transformation space is defined by k-ε-neighborhoods. If there exist two salient grids c1 and c2 that are k-connected, then we can imagine that there is a connecting path between them, making them reachable from each other in the transformation space. These definitions can be used to describe the relationships between grids in the feature space, thereby helping to understand their interactions.

[0091] The wavelet clustering process is explained in detail below, and the process includes:

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

[0093] Based on the unique information in the InSAR dataset, select an appropriate initial grid size and neighborhood parameters. The initial grid size should consider the feature range of the intercept dimension, row dimension, and column dimension to address the dimensionality mismatch between intercept and location features. The neighborhood parameters (such as the neighborhood radius ε) should ensure reasonable grid connectivity in the quantization space, avoiding clustering errors caused by neighborhoods that are too large or too small.

[0094] In the grid partitioning scheme, finding the optimal grid size is crucial. An initial value can be set based on InSAR dataset information, and then the grid parameters can be adjusted according to the clustering results. Wavelet clustering does not require prior knowledge of the exact number of clusters, but estimating the desired number of clusters helps in selecting a suitable grid size, as grids that are too small or too large are detrimental to the clustering results. When the grid is too small, pixels that should belong to the same cluster may be assigned to different clusters, leading to overquantization. In this case, the number of selected intervals in each dimension can be gradually reduced. As this process continues, the grid size gradually increases. When the grid is too large, pixels from different clusters may be assigned to the same cluster, leading to underquantization. In this case, the number of selected intervals in each dimension can be gradually increased, and the grid size will decrease. Heuristic methods are used to visually observe the cluster distribution, checking whether the distribution of each cluster and its boundaries correspond to the intercept distribution, thereby determining whether the boundaries between each cluster are reasonable. This process iterates repeatedly until a suitable cluster distribution is obtained. Thanks to the multi-resolution characteristics of wavelet transform, clustering results at different scales can be observed with the same grid size, thus reducing the number of times grid parameters need to be adjusted. Based on the above analysis, the choice of initial grid size is crucial. Using proprietary information from the InSAR dataset, the grid size can be initialized along three dimensions: intercept, rows, and columns.

[0095] Wavelet clustering is a grid-based clustering algorithm. It requires dividing the feature space into grids, then using these grids as basic clustering units to determine grid connectivity in the quantized feature space or transform space. Taking a dual-baseline InSAR system as an example, assume the grid sizes in the intercept, row, and column dimensions are represented by Δb, Δr, and Δl, respectively. Grid connectivity is determined in the quantized feature space, and a sphere of radius ε is defined to delineate the neighborhood. It must be ensured that all pixels in all grids within this neighborhood belong to the same cluster. This requires that all pixels in any grid within the neighborhood 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 the minimum distance between the center lines of adjacent clusters in the intercept dimension, i.e., d. min / 2. Therefore, the following conditions must be met:

[0096]

[0097] In the formula, This indicates rounding down. Therefore, the division of the intercept dimension should not be too coarse during the quantization process, Δb q The maximum value should not exceed Meanwhile, it can be seen from equation (4) that when the winding phase When the intercept b∈(-1,B2 / B1), B1 and B2 are two vertical baselines; therefore, we have:

[0098]

[0099] In the formula, m1 represents the number of intervals into which the intercept dimension is divided; This indicates rounding up to the nearest integer.

[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. A 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 interferogram size is N2 rows and N3 columns, then:

[0101]

[0102] In the formula, n 2l and n 2s n represents the number of fringes in the row dimension of the long baseline and short baseline interferograms, respectively. 3l and n 3s These represent the number of fringes in the column dimension of the long baseline and short baseline interferograms, respectively.

[0103] The above describes the initialization of grid parameters 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 the grid size for the q-th intercept dimension is Δb q express;

[0105] but:

[0106]

[0107] In the formula, This indicates rounding down, and ε represents the neighborhood radius. This represents the minimum distance between the centerlines of adjacent clusters in the q-th intercept dimension; therefore, Δb q The maximum value should not exceed At the same time, when the entangled phase When the intercept b∈(-1,B) q / B1), B1, B q These are two vertical baselines, B q ∈B, where B is the set of all vertical baselines except for the reference vertical baseline B1; therefore, we have:

[0108]

[0109] In the formula, This indicates the number of intervals into which the q-th intercept dimension is divided; Indicates rounding up;

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

[0111]

[0112] In the formula, n2 max and n3 max These 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 the number of pixel targets contained in each grid is counted.

[0114] S2.3: Apply wavelet transform to the quantized multidimensional feature space to obtain a new transform space, in which each grid 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 is important to note 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 (Power Processing) field. Therefore, this embodiment uses SWT, sacrificing slightly in efficiency for a significant improvement in accuracy. Subsequently, each grid in the transform space has a corresponding value, called a wavelet coefficient, thus obtaining a grid with wavelet coefficient information.

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

[0117] The process involves selecting salient grids in the transformed space whose wavelet coefficients are greater than or equal to a threshold τ. The clustering process essentially involves finding a set of k-connected salient grids. Determining the salient grid threshold τ is crucial, as it determines the proportion of noisy grids among all non-empty grids. If the threshold is set too low, some noisy grids may be incorrectly merged into a cluster; conversely, if the threshold is set too high, some normal grids may be incorrectly identified as noisy grids. The initial threshold can be estimated by setting an initial proportion of insignificant grids for all grids with non-zero wavelet coefficients based on the distribution of wavelet coefficients and the estimated proportion of noisy grids. A heuristic approach is then used to find the optimal threshold; gradually adjusting the threshold and observing the cluster boundaries helps determine the optimal threshold. Once the threshold τ is determined, grids in the transformed space with wavelet coefficients greater than or equal to τ are labeled as salient grids, while other grids are labeled as noisy grids.

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

[0119] In the new transform space, based on the identified cluster distribution, a label corresponding to the corresponding cluster is assigned to each grid cell. Since the clusters are identified in the transform space, a lookup table should be created to map them to the grid cells in the original multidimensional feature space. Then, the cluster label of each grid cell in the multidimensional feature space is assigned to all pixel targets within that grid cell. Finally, we can obtain the 2D cluster distribution of multidimensional clustered feature pixel targets at multiple resolutions. If this 2D cluster distribution seems unsuitable (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: Correct the cluster label for each pixel in the noisy cluster.

[0121] The accuracy of phase unwrapping results is affected by noise clusters identified in the wavelet clustering algorithm, and these noise clusters need to be corrected. Furthermore, some small clusters (with fewer than a threshold number of target pixels within the cluster) are pseudo-clusters caused by phase noise, and therefore can also be considered noise clusters. The basic idea for correction is to select the cluster label with the highest frequency within a preset window centered on each pixel in the noise cluster 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 the closed 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] To ensure the applicability of the present invention to the processing of large-size interferograms, the clustering process should have high execution efficiency. Assume the interferogram has N (N is very large) pixels (targets), and the clustering features have d dimensions (d is very small). The wavelet clustering algorithm first identifies all data and then divides them into corresponding grids; therefore, the time complexity of this step is O(N). For convenience, assume each dimension is divided into m equal intervals, and the original feature space is divided into K = m... d The second step is to apply SWT (Standard Wavelet Transform) to these grids; its time complexity is O(ldK) = O(dK) = O(K), where l represents the filter length used in SWT, which is usually small. If we apply SWT to decompose the space in T layers, and each layer downsamples the space by half, the time required will be:

[0124]

[0125] The above equation shows that the cost of applying wavelet transform is at most O(4K / 3) = O(K). This demonstrates that wavelet clustering algorithms 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, thus taking O(N). In summary, if we disregard the time spent reading the dataset and assigning cluster labels to targets, the time complexity of wavelet clustering is O(K). Therefore, the intermediate processes depend only on the number of grid cells and are independent of the size of the interferogram. When the size of the interferogram is very large, K is typically much smaller than N, and the efficiency advantage of wavelet clustering becomes even more pronounced.

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

[0127] This invention also provides a multi-baseline INSAR phase unwrapping system, comprising:

[0128] The clustering dataset construction module is used to obtain the intercept information of each pixel based on the interferogram 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 clustering target to obtain the clustering dataset.

[0129] The wavelet clustering processing module is used to perform clustering processing on the dataset to be clustered based on wavelet clustering.

[0130] The clustering result correction module is used to correct the cluster labels of pixels in each noisy 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 the closed-formula solution or the sparse-TSPA method, and to calculate the absolute phase of each pixel within the corresponding cluster based on the blur vector of each cluster.

[0132] It should be understood that the functional unit modules in the 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 hardware or software.

[0133] This invention also provides an electronic device, comprising:

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

[0135] A processor for loading and executing the computer program or instructions to implement the multi-baseline INSAR phase unwrapping method as described above.

[0136] This invention also provides a readable storage medium storing a computer program or instructions thereon, which, when executed by a processor, implements the multi-baseline INSAR phase unwrapping method as described above.

[0137] It is understood that the same or similar parts in the above embodiments can be referred to each other, and the contents not described in detail in some embodiments can be referred 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 with reference to specific experiments.

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

[0140] In the first experiment, a simulated scenario was used for testing. This simulated dataset was chosen because the scenario features hilly terrain with significant undulations, making phase unwrapping challenging. Furthermore, the large size of the dataset (3000×3000 pixels) effectively validates the algorithm's efficiency. The scene's height range was set from 0 meters to 162 meters, with blur heights set to 60 meters and 36 meters, respectively. To test the noise robustness of WCRFPU, phase noise with coherence coefficients of 0.8 and 0.7 was added to the short-baseline and long-baseline simulated interferograms, respectively. (Refer to the DEM as follows...) Figure 2 As shown in (a), simulated interferograms corresponding to different baselines are displayed in [the image]. Figure 2 In (b) and (c), the intercept plots obtained according to equation (4) are as follows. Figure 2 As shown in (d).

[0141] For CANOPUS, with a radius of 4 and a minimum number of points in the neighborhood of 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 set the percentage of non-significant grids for all non-zero wavelet coefficients to a certain value. The final cluster distribution is as follows. Figure 2 As shown in (f). By comparing the cluster distributions obtained by the two algorithms, it can be observed that WCRFPU identifies significantly fewer clusters than CANOPUS, and the cluster distribution of WCRFPU is closer to the actual distribution.

[0142] To present the final results more fairly and objectively, all algorithms used closed-form solutions to obtain clustered blur vectors. The phase unwrapping results and corresponding error distributions of both algorithms were annotated using the same color scale range. To comprehensively evaluate the phase unwrapping results of different algorithms, the root mean square error (RMSE) and phase unwrapping success rate (PUSR) were calculated. PUSR is defined as the proportion of pixels that successfully recovered blur numbers to the total number of pixels in the long baseline interferogram. The elevation inversion results obtained using CANOPUS and WCRFPU are shown respectively. Figure 2 In (g) and (h), the error distributions of CANOPUS and WCRFPU are shown respectively. Figure 2 In (i) and (j), 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 clusters identified is also much lower than CANOPUS, WCRFPU has a significant advantage in efficiency.

[0146] The second experiment used a DEM (458×157 pixels) from a mountainous area to simulate an MB InSAR interferogram. The reference DEM had an altitude range of 0–136.7 meters, as shown below. Figure 3 As shown in (a), the ambiguity heights corresponding to the two different baselines are 73 meters and 43.8 meters, respectively, with a baseline ratio of B2 / B1 = 5 / 3. To test the noise robustness of WCRFPU, phase noise with coherence coefficients of 0.85 and 0.8 was added to the short-baseline and long-baseline interferograms, respectively. The simulated short-baseline and long-baseline interferograms are shown in Figures (a). Figure 3 As shown in (b) and (c). The intercept for each pixel can be calculated using equation (4), as follows. Figure 3 As shown in (d), the intercept value is much smaller than the number of rows and columns. Next, we will compare TSPA, CANet-PU (a clustering analysis phase unwrapping method based on deep learning), and CANOPUS to verify the performance of WCRFPU. For TSPA, we use a residual-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 as follows... Figure 3 As shown in (e). For WCRFPU, the three-dimensional feature space is divided into a 40×40×40 grid. For all grids with non-zero wavelet coefficients, the percentage of insignificant grids is set to 60%, and the final cluster distribution is as follows. Figure 3As shown in (f). The elevation inversion results using TSPA, CANet-PU, CANOPUS, and WCRFPU are respectively shown in (f). Figure 3 As shown in (g)–(j). The elevation inversion errors of TSPA, CANet-PU, CANOPUS, and WCRFPU are respectively as follows: 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, causing errors generated during phase unwrapping to gradually accumulate in space. CANOPUS and WCRFPU can overcome the baseline ratio limitation and thus handle interferograms with small baseline ratios. However, CANOPUS identifies too many clusters and is prone to calculating incorrect blur vectors at cluster boundaries. Due to dimensional mismatch, some spurious clusters are also generated, resulting in numerous errors in the phase unwrapping results. The WCRFPU algorithm significantly reduces errors by selecting appropriate grid and neighborhood parameters, matching the intercept dimension to the position dimension, and performing additional processing on noisy clusters. In terms of efficiency, all algorithms can quickly obtain the corresponding phase unwrapping results due to the small interferogram size. However, the WCRFPU algorithm still has a slight edge.

[0148] Table 2 Performance Comparison of Four Algorithms

[0149]

[0150] The third experiment used a real dual-baseline InSAR dataset acquired from the Weinan region of Shaanxi Province, collected by the TanDEM-X satellite. This dataset was chosen because the terrain is steep and does not satisfy the phase continuity assumption. Therefore, the SBPU algorithm is unsuitable for handling such terrain, posing a significant challenge to multi-baseline phase unwrapping techniques. Figure 4 (a) shows the unwrapping results of the long baseline reference phase obtained by SRTM DEM, from which a large number of phase jumps can be observed. Figure 4 Figures (b) and (c) show the interferograms for the short baseline and the long baseline, respectively. Figure 4 Tables (e)-(h) show the long baseline phase unwrapping results obtained using 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 ambiguity number gradient. For WCRFPU, the 3D feature space is divided into a 70×70×70 grid. The percentage of insignificant grids in all grids with non-zero wavelet coefficients is set to 40%. Figure 4Figures (i)-(l) show the error distributions for 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 baseline ratio of the TanDEM-X system is approximately 2.8, it meets the baseline ratio requirements for TSPA and CANet. Therefore, TSPA and CANet did not exhibit 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-size 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 intrinsic structure of the dataset, resulting in poor clustering performance. Furthermore, CANet requires longer training time and more computational resources. WCRFPU was the fastest in this experiment. However, the high number of clusters (457 clusters) not only affected the efficiency of the cluster-by-cluster PU but also inevitably led to large errors at cluster boundaries. Conversely, WCRFPU reduced the adverse effects of dimensionality mismatch between clustering features on the clustering results; the 107 clusters identified by WaveCluster were closer to the true spatial distribution of the data. Therefore, WCRFPU demonstrates superior efficiency and accuracy when processing large-size interferograms.

[0152] Table 3. Main Interference Parameters of the Real Dataset

[0153]

[0154] Table 4 Performance Comparison of Four Algorithms

[0155]

[0156] Those skilled in the art will understand that embodiments of this application can be provided as methods, systems, or computer program products. Therefore, this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, this application can take the form of a computer program product embodied on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0157] This application is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of this application. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart... Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0158] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0159] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0160] Although embodiments of the present invention have been shown and described above, it is understood that the above embodiments are exemplary and should not be construed as limiting the present invention. Those skilled in the art can make changes, modifications, substitutions and variations to the above embodiments within the scope of the present invention.

Claims

1. A multi-baseline INSAR phase unwrapping method, characterized in that, Includes the following steps: S1: Construct the dataset to be clustered: Based on the interferograms corresponding to different vertical baselines, obtain the intercept information corresponding to each pixel, 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. S2: Wavelet clustering processing: Clustering the dataset to be clustered based on wavelet clustering; During wavelet clustering, the following conditions must be met when initializing the grid size and neighborhood radius: Assume the grid size of the q-th intercept dimension is... express; but: ； In the formula, Indicates rounding down. Represents the neighborhood radius. This represents the minimum distance between adjacent cluster centerlines in the q-th intercept dimension; therefore, The maximum value should not exceed / 2 / ( +1); at the same time, when the winding phase φ i When ∈[0, 2π), the intercept b∈(-1, B q / B1), B1, B q These are two vertical baselines; therefore, we have: ； In the formula, This represents the number of intervals into which the q-th intercept dimension is divided; Indicates rounding up; Assuming the grid sizes in the row and column dimensions are represented by ∆r and ∆l respectively, for a dual-baseline INSAR system, the number of intervals m2 and m3 in the row and column dimensions satisfy the following conditions: ； ； In the formula, n 2l and n 2s n represents the number of fringes in the row dimension of the long baseline and short baseline interferograms, respectively. 3l and n 3s N1 and N2 are the number of fringes in the column dimension of the long baseline and short baseline interferograms, respectively; N3 and N2 are the number of rows and columns of the interferograms, respectively. When using a multi-baseline INSAR system, the number of intervals m2 and m3 in the row and column dimensions must satisfy the following conditions: ； ； In the formula, and These are the maximum number of fringes in the row dimension and the maximum number of fringes in the column dimension, respectively, in all baseline interferograms; S3: Clustering result correction: Correct the cluster label for each pixel in the noisy cluster; S4: Cluster-by-cluster phase unwrapping: Based on the corrected cluster distribution, calculate the blur vector of each cluster using the closed-form formula or the sparse-TSPA method, and calculate the absolute phase of each pixel within the corresponding cluster based on the blur vector of each cluster.

2. The multi-baseline INSAR phase unwrapping method according to claim 1, characterized in that, In step S1, when there are M interferograms corresponding to different vertical baselines, one of them is selected as the main image, and M-1 sets of intercept information are calculated by combining it with other interferograms one by one. This information is then 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, thus dividing the multidimensional feature space into several grids; then the number of pixel targets contained in each grid is counted. Wavelet transform is applied to the quantized multidimensional feature space to obtain a new transform space, in which each grid has a corresponding wavelet coefficient; Significant grids with wavelet coefficients greater than or equal to the threshold τ in the new transform space are selected, and then clustering is completed by finding k-connected significant grids; The cluster labels corresponding to each grid in the new transformation space are mapped to the grids in the original multidimensional feature space. Then, the cluster label of each grid in the multidimensional feature space is assigned to all pixel targets in that grid.

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

5. 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 used.

6. The multi-baseline INSAR phase unwrapping method according to claim 1, characterized in that, Step S3 specifically includes: Noise clusters and pseudo-clusters with fewer than a threshold number of pixels within a cluster are collectively considered as noise clusters; For each pixel in a noise cluster, the cluster label with the highest frequency within a preset window centered on it is selected as the corrected cluster label for that pixel.

7. A multi-baseline INSAR phase unwrapping system, characterized in that, The system for implementing the multi-baseline INSAR phase unwrapping method as described in any one of claims 1 to 6 includes: The clustering dataset construction module is used to obtain the intercept information of each pixel based on the interferogram 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 clustering target to obtain the clustering dataset. The wavelet clustering processing module is used to perform clustering processing on the dataset to be clustered based on wavelet clustering. The clustering result correction module is used to correct the cluster labels of pixels in each noisy 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 the closed-form formula or the sparse-TSPA method, and to calculate the absolute phase of each pixel within the corresponding cluster based on the blur vector of each cluster.

8. An electronic device, characterized in that, include: A memory on which computer programs or instructions are stored; A processor for loading and executing the computer program or instructions to implement the multi-baseline INSAR phase unwrapping method as described in any one of claims 1 to 6.

9. A readable storage medium having a computer program or instructions stored thereon, characterized in that, When the computer program or instructions are executed by the processor, they implement the multi-baseline INSAR phase unwrapping method as described in any one of claims 1 to 6.

Images