Phase unwrapping method and device based on Poisson correction, terminal equipment and storage medium

Abstract

The invention discloses a phase unwrapping method and device based on Poisson correction, terminal equipment and a storage medium, and relates to the field of interferometry, and the method comprises the steps: obtaining a to-be-processed interferometric phase diagram of a target position, and inputting the to-be-processed interferometric phase diagram into a winding number prediction model to generate a predicted winding number; according to the predicted winding number, calculating a winding number gradient, and performing lower edge filtering correction to obtain a corrected winding number gradient; performing discrete Fourier transform on the corrected horizontal component of the winding number gradient to obtain a first frequency domain transform result; carrying out Laplace operator calculation on the corrected winding number gradient and then carrying out discrete Fourier transform to obtain a second frequency domain transform result; performing inverse discrete Fourier transform on a result obtained by dividing the sum of the two frequency domain transform results by a preset frequency domain adjustment factor to obtain an error winding number; and determining a real phase according to the predicted winding number and the error winding number, and further calculating the actual elevation of the target position. According to the invention, the accuracy of phase unwrapping can be improved.

Description

Technical Field

[0001] The present invention relates to the field of interferometry, and in particular, to a phase unwrapping method, device, terminal device, and storage medium based on Poisson correction. Background Art

[0002] Phase Unwrapping is a key technology in the field of Interferometric Synthetic Aperture Radar (InSAR), aiming to restore the wrapped phase and phase difference to the actual phase value through specific mathematical methods. In the InSAR process, the accuracy of the phase unwrapping algorithm directly determines the accuracy of generating a Digital Elevation Model (DEM).

[0003] Ideally, the true phase value can be completed simply through integration. However, during the measurement process, the obtained phase information is usually affected by interference such as phase noise, undersampling, and fringe aliasing. Therefore, the phase unwrapping process needs to be cleverly designed to obtain more accurate calculation results.

[0004] However, in the traditional least squares unwrapping algorithm, due to the influence of factors such as noise on the wrapped phase, there is a large deviation between the wrapped phase gradient and the absolute phase gradient. These errors will be evenly distributed to the global unwrapping result, resulting in the situation where the unwrapping result deviates from the true value. Summary of the Invention

[0005] Embodiments of the present invention provide a phase unwrapping method, device, terminal device, and storage medium based on Poisson correction, which can improve the accuracy of the unwrapped phase and solve the problem that the errors of the traditional least squares unwrapping algorithm are evenly distributed to the global unwrapping result, resulting in the deviation of the unwrapping result from the true value.

[0006] An embodiment of the present invention provides a phase unwrapping method based on Poisson correction, including: Obtain an interferometric phase map to be processed at a target position; Input the interferometric phase map to be processed into a trained wrapped number prediction model, so that the wrapped number prediction model generates a predicted wrapped number according to the interferometric phase map to be processed; Calculate the wrapped number gradient according to the predicted wrapped number, and perform down-edge filtering correction on the wrapped number gradient to obtain a corrected wrapped number gradient; Perform a discrete Fourier transform on the horizontal component of the corrected wrapped number gradient to obtain a first frequency-domain transform result; Perform a Laplacian operator calculation on the corrected wrapped number gradient and then perform a discrete Fourier transform to obtain a second frequency-domain transform result; Perform an inverse discrete Fourier transform on the result of dividing the sum of the first frequency-domain transformation result and the second frequency-domain transformation result by a preset frequency-domain adjustment factor to obtain an error winding number; Determine the true phase based on the predicted winding number and the error winding number, and calculate the actual elevation of the target position based on the true phase.

[0007] Further, generate a predicted winding number based on the interferometric phase map to be processed, including: Perform image partitioning on the interferometric phase map to be processed to obtain a number of sub-block images; For each sub-block image, perform multi-scale pyramid pooling on the sub-block image to obtain a multi-scale feature map; After performing a frequency-domain transformation on the multi-scale feature map, perform spectral feature fusion to obtain a frequency-domain feature map; Decode the frequency-domain feature map to obtain the predicted winding number corresponding to the sub-block image; Determine the predicted winding number corresponding to the interferometric phase map to be processed based on the predicted winding numbers of all sub-block images.

[0008] Further, calculate the winding number gradient based on the predicted winding number, including: Calculate the winding number gradient based on the predicted winding number through the following formula: ; Wherein, represents the winding number gradient, represents the gradient operator, , represents the predicted winding number, represents the rounding process.

[0009] Further, preset the frequency-domain adjustment factor, including: ; Wherein, represents the preset frequency-domain adjustment factor, represents the length of the interferometric phase map, represents the width of the interferometric phase map, and represent the coordinate values in the horizontal direction and the vertical direction on the two-dimensional plane of the interferometric phase map.

[0010] Further, determine the true phase based on the predicted winding number and the error winding number, including: Determine the interferometric phase based on the interferometric phase map to be processed; Calculate the true phase based on the interferometric phase, the predicted winding number, and the error winding number through the following formula: ; Among them, represents the true phase, represents the interference phase, represents the predicted winding number, represents the error winding number.

[0011] Furthermore, according to the true phase, the actual elevation of the target position is calculated, including: Obtain the antenna position height at the measured target position, the local incident angle of the antenna position, the interference baseline, the wavelength of the radar wave, the baseline tilt angle, and the distance from the antenna position to the target position; According to the antenna position height, the local incident angle of the antenna position, the interference baseline, the wavelength of the radar wave, the baseline tilt angle, the distance from the antenna position to the target position, and the true phase, calculate the actual elevation of the target position.

[0012] Furthermore, the winding number prediction model is determined in the following manner: Obtain a number of training samples; each training sample includes: an interference phase training map and the corresponding actual winding number; Input a number of training samples into the winding number prediction model to be trained, so that the winding number prediction model is trained with the interference phase training map as the input and the predicted winding number corresponding to the interference phase training map as the output, and during the training process, according to the predicted winding number corresponding to the interference phase training map and its corresponding actual winding number, calculate the loss function, and adjust the parameters of the winding number prediction model according to the loss function until the loss function converges to obtain the trained winding number prediction model; Among them, each training sample is determined in the following manner: Obtain a Gaussian random matrix; Expand the Gaussian random matrix to obtain an expanded matrix; Adjust the expanded matrix to a preset size and introduce Gaussian noise to generate a true phase training map; Perform winding processing on the true phase training map to obtain an interference phase training map; According to the true phase training map and the interference phase training map, determine the actual winding number corresponding to the interference phase training map.

[0013] Based on the above method item embodiments, the present invention correspondingly provides device item embodiments, including: an interference phase map acquisition module, a winding number prediction module, a winding number gradient correction module, a first frequency domain transformation module, a second frequency domain transformation module, an error winding number determination module, and a phase unwrapping application module; The interference phase map acquisition module is used to acquire the interference phase map to be processed at the target position; The winding number prediction module is used to input the interference phase map to be processed into the trained winding number prediction model, so that the winding number prediction model generates a predicted winding number according to the interference phase map to be processed; The winding number gradient correction module is used to calculate the winding number gradient according to the predicted winding number, and perform a down-edge filtering correction on the winding number gradient to obtain the corrected winding number gradient; The first frequency domain transformation module is used to perform a discrete Fourier transform on the horizontal component of the corrected winding number gradient to obtain a first frequency domain transformation result; The second frequency domain transformation module is used to perform a Laplacian operator calculation on the corrected winding number gradient and then perform a discrete Fourier transform to obtain a second frequency domain transformation result; The error winding number determination module is used to perform an inverse discrete Fourier transform on the result of dividing the sum of the first frequency domain transformation result and the second frequency domain transformation result by a preset frequency domain adjustment factor to obtain an error winding number; The phase unwrapping application module is used to determine the true phase according to the predicted winding number and the error winding number, and calculate the actual elevation of the target position according to the true phase.

[0014] Based on the above method item embodiments, the present invention correspondingly provides a terminal device item embodiment, including: a processor, a memory, and a computer program stored in the memory and configured to be executed by the processor. When the processor executes the computer program, the steps of the phase unwrapping method based on Poisson correction as described in the present invention are implemented.

[0015] Based on the above method item embodiments, the present invention correspondingly provides a computer-readable storage medium item embodiment, including: a stored computer program. When the computer program runs, it controls the device where the computer-readable storage medium is located to execute the steps of the phase unwrapping method based on Poisson correction as described in the present invention.

[0016] Compared with the prior art, the beneficial effects of the embodiments of the present solution are as follows: The present invention obtains an interferometric phase diagram to be processed at a target location, inputs the interferometric phase diagram to be processed into a trained winding number prediction model. The model is trained with a large amount of data and can learn the complex relationship between the phase and the winding number to generate a predicted winding number. According to the predicted winding number, a winding number gradient is calculated, which reflects the change amount of the winding number between adjacent pixels. The lower edge filtering correction is performed on the winding number gradient. Through the filtering correction, the small gradient changes caused by noise will be suppressed, and the gradients caused by prediction errors will be retained, obtaining a corrected winding number gradient. Then, the horizontal component of the corrected winding number gradient is subjected to a discrete Fourier transform to obtain a first frequency-domain transform result. The low-frequency components correspond to the slowly changing gradients, and the high-frequency components correspond to the locally abnormal gradients. After calculating the Laplacian operator on the corrected winding number gradient and then performing a discrete Fourier transform, a second frequency-domain transform result is obtained. The Laplacian operator calculation can highlight the regions where the gradient changes particularly fast or particularly slow, and then through the discrete Fourier transform, these characteristic information is separated according to the frequency. The result of dividing the sum of the first frequency-domain transform result and the second frequency-domain transform result by a preset frequency-domain adjustment factor is subjected to an inverse discrete Fourier transform to obtain an error winding number. These two sets of information with different frequencies are combined to form a multi-dimensional description of the error. Finally, according to the predicted winding number and the error winding number, the true phase is determined, and according to the true phase, the actual elevation of the target location is calculated.

[0017] In summary, based on the winding number predicted by the model, the present invention decomposes the error winding number through a two-dimensional frequency-domain analysis combining the discrete Fourier transform of the horizontal component of the winding number gradient and the discrete Fourier transform combined with the Laplacian operator, effectively solving the problem that the errors of the traditional least squares unwrapping algorithm are evenly distributed to the global unwrapping result, resulting in the unwrapping result deviating from the true value, and improving the accuracy of the unwrapped phase. BRIEF DESCRIPTION OF THE DRAWINGS

[0018] Figure 1 is a schematic flowchart of a phase unwrapping method based on Poisson correction provided by an embodiment of the present invention; Figure 2 is a schematic geometric structure diagram of spaceborne SAR interferometric altimetry provided by an embodiment of the present invention; Figure 3 is a schematic network structure diagram of a winding number prediction model provided by an embodiment of the present invention; Figure 4 is a schematic flowchart of the training process of a winding number prediction model provided by an embodiment of the present invention; Figure 5 is a schematic structure diagram of a phase unwrapping device based on Poisson correction provided by an embodiment of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0019] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without making creative efforts belong to the scope of protection of the present invention.

[0020] In the description of the present invention, it should be understood that the terms "first" and "second" are only used for descriptive purposes and cannot be construed as indicating or implying relative importance or implicitly specifying the quantity of the indicated technical features.

[0021] As Figure 1 shown, in order to solve the problem that the errors of the traditional least squares phase unwrapping algorithm are evenly distributed to the global unwrapping result, resulting in the unwrapping result deviating from the true value, an embodiment of the present invention provides a phase unwrapping method based on Poisson correction, which at least includes the following steps: Step S1: Obtain the interferometric phase map to be processed at the target position; For step S1, as Figure 2 shown is a schematic structural diagram of the spaceborne SAR interferometric altimetry geometric model. The satellite emits and receives radar waves at different positions S1 and S2. The radar wave emitted by the satellite from position S1 is reflected by the target position P and then received. The radar wave emitted from position S2 is also reflected by P and then received. Due to the spatial position difference of the target position P relative to S1 and S2, and the different propagation path lengths of the radar waves, a phase difference is generated when these two radar waves are received. These phase difference information constitutes the interferometric phase map of the target position.

[0022] It should be noted that during the process of extracting phase information from the interferometric phase map, the absolute phase value will be wrapped within the range of [-π, π) or [0, 2π), which is called the wrapped phase value. And the task of phase unwrapping is to determine the wrapping number from the wrapped phase, so as to restore the true phase value. The wrapped phase value and the true phase value can be expressed by the following formula: ; where represents the true phase value, represents the wrapped phase value, represents the wrapping number, as an integer function, representing the integer multiple of 2π corresponding to the true phase value.

[0023] Step S2: Input the interferometric phase map to be processed into the trained wrapping number prediction model, so that the wrapping number prediction model generates a predicted wrapping number according to the interferometric phase map to be processed; In a preferred embodiment, generating a predicted winding number according to the interference phase map to be processed includes: Performing image partitioning on the interference phase map to be processed to obtain a plurality of sub-block images; For each sub-block image, performing multi-scale pyramid pooling on the sub-block image to obtain a multi-scale feature map; After performing frequency domain transformation on the multi-scale feature map, performing spectral feature fusion to obtain a frequency domain feature map; Decoding the frequency domain feature map to obtain the predicted winding number corresponding to the sub-block image; Determining the predicted winding number corresponding to the interference phase map to be processed according to the predicted winding numbers of all sub-block images.

[0024] For step S2, inputting the interference phase map to be processed obtained in step S1 into the trained winding number prediction model. The winding number prediction model is trained to reflect the complex relationship between the phase and the winding number, and generates a predicted winding number according to the interference phase map to be processed.

[0025] Specifically, as Figure 3 shown, the winding number prediction model adopts a PUformer network structure. First, before encoding, the interference phase map to be processed is divided into images, and it is divided into a plurality of non-overlapping sub-block images (patches) of 4×4 pixels. Each sub-block image contains the phase values of 16 pixel points. If the size of the original interference phase map is H×W, then a total of (H / 4)×(W / 4) sub-block images are obtained after division. Each sub-block image is used as the smallest feature unit for the subsequent encoding process of the winding number prediction model, effectively retaining local structural information while reducing the amount of calculation.

[0026] The phase unwrapping based on PUFormer is divided into two stages. In the first stage, for each sub-block image, feature extraction and encoding are performed through the MiT (Mix Vision Transformer) encoder set inside the winding number prediction model. Specifically, the MiT encoder contains four stages (Stage 1 to Stage 4). Each stage is a transformer module (Transformerblock) composed of multiple Transformer blocks, and spatial compression and channel expansion are achieved through downsampling operations between stages. The feature maps output by each stage respectively correspond to 1 / 4, 1 / 8, 1 / 16, and 1 / 32 resolutions of the interference phase map, forming a multi-scale feature map to capture multi-level phase information from details to the global.

[0027] In the second stage, for each scale of feature map, through the pointwise convolutional decoder in the winding number prediction model, a frequency domain transformation is performed on it to obtain the representation of each scale of feature map in the frequency domain. Then, the spectrum is divided into low-frequency components and high-frequency components. Low-frequency features retain the large-scale structural information of the image, such as background regions, contour lines, etc., and have good global semantic consistency; high-frequency features correspond to the mutations and detailed changes at the boundaries of the winding number, manifested as local features such as phase jumps, and are suitable for boundary perception and sharpening processing. The light blue arrows in the figure represent the step-by-step upsampling operation, which upsamples the low-resolution features (1 / 32, 1 / 16, 1 / 8) layer by layer (Upsample) and finally unifies them to the 1 / 4 resolution. After the feature maps of each scale are consistent in the spatial dimension, they are fused through a frequency fusion block in a concatenated (Concat) manner, and the fused frequency domain feature map is input into the internal convolutional fusion layer (Conv fusion Layer) for feature integration, normalization processing, and introduction of a non-linear transformation. The output processed and transformed feature map (Feature map) is input into the pixelwise convolutional layer (Pixelwise Conv Layer) for refined feature extraction, and finally the predicted winding number is output. 。

[0028] Finally, the predicted winding number of each sub-block image reflects the phase winding situation of the local area. By combining the predicted winding numbers of all sub-block images, the predicted winding number corresponding to the interference phase map to be processed is obtained.

[0029] The present invention introduces a frequency-aware pointwise convolutional decoder. This module can sense the gradient jumps existing in the interference phase and obtain relevant high-frequency and low-frequency information, improving the accuracy of the semantic segmentation model for predicting gradient jump regions. In addition, the convolutional fusion layer of the present invention is a lightweight decoder composed of a convolution (Conv), a batch normalization (BN), and an activation function (ReLU). It replaces the widely used decoder with a large number of parameters with a lightweight pointwise convolution. While reducing the model parameters, it can also maintain the prediction accuracy.

[0030] Next, the training of the winding number prediction model will be described in detail: As Figure 4 shown, the training of the above-mentioned winding number prediction model includes the following steps: Step S21: Obtain a number of training samples; each training sample includes: an interference phase training map and the corresponding actual winding number; Among them, each training sample is determined by the following method: Obtain a Gaussian random matrix; Expand the Gaussian random matrix to obtain the expanded matrix; Adjust the expanded matrix to a preset size and introduce Gaussian noise to generate a true phase training map; Perform wrapping processing on the true phase training map to obtain an interferometric phase training map; Determine the actual wrapping number corresponding to the interferometric phase training map according to the true phase training map and the interferometric phase training map; For step S21, obtain training samples for training the wrapping number prediction model. Each training sample includes an interferometric phase training map and the corresponding actual wrapping number. When training a deep learning algorithm, a large amount of data is required to determine the mapping relationship between the input wrapped phase and the output wrapping number. However, it is difficult to obtain a large-scale real wrapped phase map and its corresponding accurate wrapping number. Therefore, training samples are generated in the following way.

[0031] First, obtain a Gaussian random matrix with a matrix size of 2*2 to 10*10 and a value range of 10 to 150. Then, use the mathematical fractal algorithm to expand the Gaussian random matrix. In this embodiment, the diamond-square algorithm (The diamondsquare method) is used to process and augment the random matrix 2 to 5 times to generate a more complex true phase matrix. Adjust the expanded matrix to a preset size. In this embodiment, it is adjusted to 256*256. Interpolation is performed on a matrix with a smaller size to obtain a smoother true phase, and downsampling is performed on a matrix with a larger size to obtain a rougher true phase, so as to obtain phase data with both smooth and rough characteristics. Then, introduce Gaussian noise with a standard deviation between 0 and 0.3 to simulate the interference to the phase in the real world and generate a true phase training map.

[0032] Finally, perform wrapping processing on the generated true phase training map through the following formula to obtain an interferometric phase training map: ; Then, according to the true phase training map and the interferometric phase training map, determine the actual wrapping number corresponding to the interferometric phase training map through the following formula: ; where, represents the interferometric phase training map, represents the true phase training map, represents the modulo operation, represents the actual wrapping number corresponding to the interferometric phase training map.

[0033] Step S22: Input a number of training samples into the winding number prediction model to be trained, so that the winding number prediction model takes the interference phase training diagram as the input and the predicted winding number corresponding to the interference phase training diagram as the output for training. During the training process, calculate the loss function based on the predicted winding number corresponding to the interference phase training diagram and its corresponding actual winding number, and adjust the parameters of the winding number prediction model according to the loss function until the loss function converges, obtaining the trained winding number prediction model; For step S22, input the number of training samples obtained in step S21 into the winding number prediction model to be trained. When the winding number prediction model is being trained, with the interference phase training diagram as the input, a series of feature extraction, encoding, and decoding operations will be performed on the input interference phase training diagram inside the model, and the predicted winding number for this interference phase training diagram will be output. During the training process, compare the predicted winding number output by the model with the actual winding number corresponding to it in the training samples, calculate the average value (MSE) of the square of the difference between the predicted winding number and the actual winding number, and use it as the loss function value. The loss function is an index used to evaluate the difference between the model's prediction result and the true result. Then, according to the calculated loss function value, adjust the parameters of the winding number prediction model through gradient descent. Continuously repeat the above training process. As the training progresses, the loss function value will gradually decrease. When the loss function value no longer decreases significantly and reaches a relatively stable state, that is, the so-called convergence, it means that the model has learned enough knowledge. At this time, the obtained is the trained winding number prediction model.

[0034] Step S3: Calculate the winding number gradient based on the predicted winding number, and perform a lower edge filtering correction on the winding number gradient to obtain the corrected winding number gradient; In a preferred embodiment, calculating the winding number gradient based on the predicted winding number includes: Calculate the winding number gradient according to the predicted winding number through the following formula: ; where represents the winding number gradient, represents the gradient operator, , represents the predicted winding number, represents the rounding process.

[0035] For step S3, obtain the predicted winding number through step S2 After that, in order to detect the rate of change of the phase, calculate the winding number gradient through the formula , which reflects the change of the winding number in space. The rounding operation means converting the continuous phase change into a discrete integer jump, that is, determining the increment of the winding number at each position.

[0036] When the phase gradient changes drastically, small gradient changes may cause errors in the estimation of the gradient field of the number of wrong windings. To meet the needs of correction, a phase filter is used to perform a lower-edge filtering correction on the gradient of the winding number to obtain the corrected gradient of the winding number : ; Through the phase filter, small gradient changes caused by noise will be suppressed, and the gradients caused by prediction errors will be retained.

[0037] Step S4: Perform a discrete Fourier transform on the horizontal component of the corrected gradient of the winding number to obtain a first frequency-domain transformation result; For step S4, the corrected gradient of the winding number is divided into the horizontal component of the corrected gradient of the winding number and the vertical component of the corrected gradient of the winding number , and the horizontal component reflects the change of the gradient of the winding number in the horizontal direction of the two-dimensional interference phase diagram, and the vertical component reflects the change of the gradient of the winding number in the vertical direction of the two-dimensional interference phase diagram.

[0038] Then, perform a discrete Fourier transform (DCT) on the horizontal component of the corrected gradient of the winding number to obtain a first frequency-domain transformation result: ; where, represents the first frequency-domain transformation result, represents the horizontal component of the corrected gradient of the winding number, represents the discrete Fourier transform process.

[0039] By performing a discrete Fourier transform, the interference phase diagram is transformed from the spatial domain to the frequency domain for analysis to understand the information of different frequency components in the horizontal component .

[0040] Step S5: Perform a Laplacian operator calculation on the corrected gradient of the winding number and then perform a discrete Fourier transform to obtain a second frequency-domain transformation result; For step S5, similarly, perform a Laplacian operator calculation on the corrected gradient of the winding number and then perform a discrete Fourier transform to obtain a second frequency-domain transformation result: ; where, represents the second frequency-domain transformation result, represents the Laplacian operator, Represents the corrected winding number gradient.

[0041] It should be noted that in the present invention, partial frequency domain information is obtained by performing DCT transformation on the horizontal component in step S4. After the winding number gradient processed by the Laplace operator in step S5 is subjected to DCT transformation, frequency domain information is obtained from another perspective. The two different frequency domain transformation results complement each other and can more comprehensively describe the frequency characteristics of the winding number gradient.

[0042] Step S6: Perform inverse discrete Fourier transform on the result of dividing the sum of the first frequency domain transformation result and the second frequency domain transformation result by a preset frequency domain adjustment factor to obtain an error winding number; In a preferred embodiment, the preset frequency domain adjustment factor includes: ; Wherein, Represents the preset frequency domain adjustment factor, Represents the length of the interference phase diagram, Represents the width of the interference phase diagram, And Represents the coordinate values in the horizontal direction and the vertical direction on the two-dimensional plane of the interference phase diagram.

[0043] For step S6, perform inverse discrete Fourier transform on the result of dividing the sum of the first frequency domain transformation result and the second frequency domain transformation result by a preset frequency domain adjustment factor through the following formula to obtain an error winding number: ; ; ; Wherein, Represents the result of dividing the sum of the first frequency domain transformation result and the second frequency domain transformation result by a preset frequency domain adjustment factor, Represents the preset frequency domain adjustment factor, Represents the length of the interference phase diagram, Represents the width of the interference phase diagram, And Represents the coordinate values in the horizontal direction and the vertical direction on the two-dimensional plane of the interference phase diagram, Represents the error winding number, Represents the inverse discrete Fourier transform process.

[0044] It should be noted that in the present invention, first, the first frequency domain transformation result obtained in step S4 Add them up, fuse the frequency-domain information obtained from different perspectives, and comprehensively consider the gradient characteristics in the horizontal direction and the frequency-domain characteristics after enhancing the features with the Laplacian operator. Moreover, through the Fourier transform, the Laplacian operator in the spatial domain can be converted into an algebraic operation in the frequency domain, thus simplifying the solution of the equation. After adding the two, the numerical ranges and magnitudes of different frequency components may vary greatly. Dividing the sum of the two by the preset frequency-domain adjustment factor param can normalize the fused frequency-domain information, making different frequency components in a relatively unified and appropriate numerical range. Finally, perform the inverse discrete Fourier transform (IDCT) to convert the frequency-domain information back to the spatial domain to obtain the error winding number.

[0045] This process is similar to solving the Poisson equation to recover the phase information. In the phase unwrapping problem, the Poisson equation is used to describe the relationship between the phase gradient and the true phase. By performing the inverse transform on the result after frequency-domain processing, it is equivalent to constructing the distribution of the error winding number in the spatial domain based on the information obtained from the frequency-domain analysis. This error winding number reflects the difference between the predicted winding number and the true winding number.

[0046] Step S7: Determine the true phase according to the predicted winding number and the error winding number, and calculate the actual elevation of the target position according to the true phase.

[0047] In a preferred embodiment, determining the true phase according to the predicted winding number and the error winding number includes: Determine the interference phase according to the interference phase diagram to be processed; Calculate the true phase through the following formula according to the interference phase, the predicted winding number, and the error winding number: ; where, represents the true phase, represents the interference phase, represents the predicted winding number, represents the error winding number.

[0048] In a preferred embodiment, calculating the actual elevation of the target position according to the true phase includes: Obtain the antenna position height at the measured target position, the local incident angle at the antenna position, the interference baseline, the wavelength of the radar wave, the baseline tilt angle, and the distance from the antenna position to the target position; Calculate the actual elevation of the target position according to the antenna position height, the local incident angle at the antenna position, the interference baseline, the wavelength of the radar wave, the baseline tilt angle, the distance from the antenna position to the target position, and the true phase.

[0049] For step S7, since the interferometric phase diagram is obtained using interferometry technology, based on the principle of wave interference. When two coherent waves meet, an interference phenomenon will occur due to factors such as the path difference. By directly reading the interference fringe distribution in the interferometric phase diagram to be processed, the interference phase can be obtained.

[0050] The present invention has obtained the predicted winding number through the winding number prediction model in step S2 and determined the error winding number through step S6 Then, through the calculation formula compensate for the missing or incorrect phase values caused by phase winding, so as to recover the true phase.

[0051] Apply the unwrapped true phase to spaceborne SAR interferometric altimetry. As Figure 2 shown, the antenna positions include the main image antenna position and the secondary image antenna position. Both emit radar waves towards the same target position P. The two radar waves propagate to the target along different paths and are reflected back to the satellite. Due to the existence of the propagation paths, when these two coherent radar waves meet, an interference phenomenon will occur, obtaining the interferometric phase diagram to be processed. P is the target position, S1 is the main image antenna position, S2 is the secondary image antenna position, H is the height of the main image antenna position, θ is the local incident angle of the main image antenna position, B is the interferometric baseline, that is, the distance between the two antenna positions, a represents the baseline tilt angle, that is, the angle between the interferometric baseline and the horizontal direction, R1 is the distance from the main image antenna position to the target position P, R2 is the distance from the secondary image antenna position to the target position P, h is the elevation of the target position P, P0 is a point on the reference plane, the slant range of P0 and P in the main image is the same, and the distances from the main image antenna position S1 to P and to P0 are equal, and its elevation is 0, and θ0 is the local incident angle of P0.

[0052] Obtain the height H of the main image antenna position, the local incident angle θ of the main image antenna position, the interferometric baseline B, the wavelength λ of the radar wave, the baseline tilt angle a, and the distance R1 from the main image antenna position to the target position, combined with the true phase Through the following formula, the actual elevation h of the target position can be calculated: ; where represents the actual elevation of the target position, represents the wavelength of the radar wave, represents the distance from the main image antenna position to the target position, represents the local incident angle of the main image antenna position, represents the interferometric baseline, represents the baseline tilt angle, represents the true phase.

[0053] The present invention initially predicts the winding number through a deep learning model. Based on the analysis of the winding number gradient, by leveraging the relationship between the phase gradient and the true phase in a Poisson-like equation, the error winding number is determined, thereby realizing the correction of the winding number and restoring the true phase. In the face of extreme situations such as large terrain undulations and dense fringes, the success rate of phase unwrapping is higher than that of traditional methods. At the same time, when encountering data with a different distribution from the training set, general deep learning methods are vulnerable to the limitations of training data, resulting in inaccurate prediction of the winding number. On this basis, the present invention improves the phase unwrapping effect through winding number correction, making the phase unwrapping effect better than that of traditional methods and deep learning methods.

[0054] As Figure 5 shown, based on the above method item embodiments, corresponding apparatus item embodiments are provided; An embodiment of the present invention provides a phase unwrapping apparatus based on Poisson correction, including: an interference phase diagram acquisition module, a winding number prediction module, a winding number gradient correction module, a first frequency domain transformation module, a second frequency domain transformation module, an error winding number determination module, and a phase unwrapping application module; The interference phase diagram acquisition module is used to acquire the interference phase diagram to be processed at the target position; The winding number prediction module is used to input the interference phase diagram to be processed into the trained winding number prediction model, so that the winding number prediction model generates a predicted winding number according to the interference phase diagram to be processed; The winding number gradient correction module is used to calculate the winding number gradient according to the predicted winding number, and perform down-edge filtering correction on the winding number gradient to obtain the corrected winding number gradient; The first frequency domain transformation module is used to perform discrete Fourier transform on the horizontal component of the corrected winding number gradient to obtain the first frequency domain transformation result; The second frequency domain transformation module is used to perform Laplace operator calculation on the corrected winding number gradient and then perform discrete Fourier transform to obtain the second frequency domain transformation result; The error winding number determination module is used to perform inverse discrete Fourier transform on the result of dividing the sum of the first frequency domain transformation result and the second frequency domain transformation result by the preset frequency domain adjustment factor to obtain the error winding number; The phase unwrapping application module is used to determine the true phase according to the predicted winding number and the error winding number, and calculate the actual elevation of the target position according to the true phase.

[0055] It can be understood that the above apparatus item embodiments correspond to the method item embodiments of the present invention, and can implement the phase unwrapping method based on Poisson correction provided by any one of the above method item embodiments of the present invention.

[0056] It should be noted that the device embodiments described above are merely illustrative, and some or all of the modules can be selected according to actual needs to achieve the purpose of the solution of this embodiment. In addition, in the attached drawings of the device embodiments provided by the present invention, the connection relationships between the modules indicate that there is a communication connection between them, which can be specifically implemented as one or more communication buses or signal lines. Those of ordinary skill in the art can understand and implement it without creative work.

[0057] Based on the above embodiments of the phase unwrapping method based on Poisson correction, another embodiment of the present invention provides a terminal device, which includes a processor, a memory, and a computer program stored in the memory and configured to be executed by the processor. When the processor executes the computer program, the phase unwrapping method based on Poisson correction according to any embodiment of the present invention is implemented.

[0058] Exemplarily, in this embodiment, the computer program can be divided into one or more modules, and the one or more modules are stored in the memory and executed by the processor to complete the present invention. The one or more module elements can be a series of computer program instruction segments capable of completing specific functions, and the instruction segments are used to describe the execution process of the computer program in the terminal device.

[0059] The terminal device can be a computing device such as a desktop computer, a notebook, a palm computer, and a cloud server. The terminal device may include, but is not limited to, a processor and a memory.

[0060] The so-called processor may be a central processing unit (CPU), or may also be other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), off-the-shelf programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. The general-purpose processor may be a microprocessor or the processor may also be any conventional processor, etc. The processor is the control center of the terminal device, and connects various parts of the entire terminal device through various interfaces and lines.

[0061] Based on the above method embodiments, another embodiment is provided: A computer-readable storage medium provided by another embodiment of the present invention includes a stored computer program, wherein when the computer program runs, it controls the device where the computer-readable storage medium is located to execute the phase unwrapping method based on Poisson correction described in any one of the above method embodiments of the present invention.

[0062] Among them, the module / unit integrated in the phase unwrapping device / terminal device based on Poisson correction, if implemented in the form of a software functional unit and sold or used as an independent product, can be stored in a computer-readable storage medium. Based on such an understanding, to implement all or part of the processes in the above method embodiments of the present invention, it can also be completed by instructing relevant hardware through a computer program. The computer program can be stored in a computer-readable storage medium. When the computer program is executed by a processor, the steps of the above various method embodiments can be implemented. Among them, the computer program includes computer program code, and the computer program code can be in the form of source code, object code, executable file or some intermediate form, etc. The computer-readable medium can include: any entity or device capable of carrying the computer program code, recording medium, USB flash drive, mobile hard disk, magnetic disk, optical disc, computer memory, read-only memory (ROM, Read-Only Memory), random access memory (RAM, Random Access Memory), electrical carrier signal, telecommunication signal, and software distribution medium, etc.

[0063] The above are the preferred embodiments of the present invention. It should be noted that for those of ordinary skill in the art, without departing from the principle of the present invention, several improvements and refinements can be made, and these improvements and refinements are also regarded as the protection scope of the present invention.

Claims

1. A phase unwrapping method based on Poisson correction, characterized in that, Including: Obtain the interferometric phase diagram to be processed at the target position; Input the interferometric phase diagram to be processed into the trained winding number prediction model, so that the winding number prediction model generates a predicted winding number according to the interferometric phase diagram to be processed; Calculate the winding number gradient according to the predicted winding number, and perform a lower edge filtering correction on the winding number gradient to obtain the corrected winding number gradient; Perform a discrete Fourier transform on the horizontal component of the corrected winding number gradient to obtain a first frequency domain transformation result; Perform a Laplacian operator calculation on the corrected winding number gradient and then perform a discrete Fourier transform to obtain a second frequency domain transformation result; Perform an inverse discrete Fourier transform on the result of dividing the sum of the first frequency domain transformation result and the second frequency domain transformation result by a preset frequency domain adjustment factor to obtain an error winding number; Determine the true phase according to the predicted winding number and the error winding number, and calculate the actual elevation of the target position according to the true phase.

2. The phase unwrapping method based on Poisson correction according to claim 1, wherein Generating a predicted winding number according to the interferometric phase diagram to be processed includes: Perform image partitioning on the interferometric phase diagram to be processed to obtain a number of sub-block images; For each sub-block image, perform multi-scale pyramid pooling on the sub-block image to obtain a multi-scale feature map; After performing a frequency domain transformation on the multi-scale feature map, perform spectral feature fusion to obtain a frequency domain feature map; Decode the frequency domain feature map to obtain the predicted winding number corresponding to the sub-block image; Determine the predicted winding number corresponding to the interferometric phase diagram to be processed according to the predicted winding numbers of all sub-block images.

3. The phase unwrapping method based on Poisson correction according to claim 1, wherein Calculating the winding number gradient according to the predicted winding number includes: Calculate the winding number gradient according to the predicted winding number through the following formula: ； Among them, represents the winding number gradient, represents the gradient operator, , represents the predicted winding number, represents the rounding process.

4. The phase unwrapping method based on Poisson correction according to claim 1, characterized in that The preset frequency domain adjustment factor includes: ； Among them, represents a preset frequency domain adjustment factor, represents the length of the interference phase diagram, represents the width of the interference phase diagram, and represent the coordinate values in the horizontal direction and the vertical direction on the two-dimensional plane of the interference phase diagram.

5. The phase unwrapping method based on Poisson correction according to claim 1, characterized in that, Determining the true phase according to the predicted winding number and the error winding number includes: Determine the interferometric phase according to the interferometric phase diagram to be processed; Calculate the true phase through the following formula according to the interferometric phase, the predicted winding number and the error winding number: ； Among them, represents the true phase, represents the interference phase, represents the predicted winding number, represents the error winding number.

6. The phase unwrapping method based on Poisson correction according to claim 1, wherein Calculating the actual elevation of the target position according to the true phase includes: Obtain the antenna position height for measuring the target position, the local incident angle of the antenna position, the interferometric baseline, the wavelength of the radar wave, the baseline tilt angle, and the distance from the antenna position to the target position; Calculate the actual elevation of the target position according to the antenna position height, the local incident angle of the antenna position, the interferometric baseline, the wavelength of the radar wave, the baseline tilt angle, the distance from the antenna position to the target position, and the true phase.

7. The phase unwrapping method based on Poisson correction according to claim 1, characterized in that, The winding number prediction model is determined by the following method: Obtain a number of training samples; each training sample includes: an interferometric phase training diagram and a corresponding actual winding number; Input a number of training samples into the winding number prediction model to be trained, so that the winding number prediction model is trained with the interferometric phase training diagram as the input and the predicted winding number corresponding to the interferometric phase training diagram as the output, and during the training process, calculate the loss function according to the predicted winding number corresponding to the interferometric phase training diagram and its corresponding actual winding number, and adjust the parameters of the winding number prediction model according to the loss function until the loss function converges to obtain the trained winding number prediction model; Among them, each training sample is determined in the following manner: Obtain a Gaussian random matrix; Expand the Gaussian random matrix to obtain an expanded matrix; Adjust the expanded matrix to a preset size and introduce Gaussian noise to generate a true phase training map; Perform wrapping processing on the true phase training map to obtain an interferometric phase training map; Determine the actual wrapping number corresponding to the interferometric phase training map according to the true phase training map and the interferometric phase training map.

8. A phase unwrapping device based on Poisson correction, characterized in that, Including: An interferometric phase map acquisition module, a wrapping number prediction module, a wrapping number gradient correction module, a first frequency domain transformation module, a second frequency domain transformation module, an error wrapping number determination module, and a phase unwrapping application module; The interferometric phase map acquisition module is used to acquire an interferometric phase map to be processed at a target position; The wrapping number prediction module is used to input the interferometric phase map to be processed into a trained wrapping number prediction model, so that the wrapping number prediction model generates a predicted wrapping number according to the interferometric phase map to be processed; The wrapping number gradient correction module is used to calculate a wrapping number gradient according to the predicted wrapping number and perform downward edge filtering correction on the wrapping number gradient to obtain a corrected wrapping number gradient; The first frequency domain transformation module is used to perform a discrete Fourier transform on the horizontal component of the corrected wrapping number gradient to obtain a first frequency domain transformation result; The second frequency domain transformation module is used to perform a Laplace operator calculation on the corrected wrapping number gradient and then perform a discrete Fourier transform to obtain a second frequency domain transformation result; The error wrapping number determination module is used to perform an inverse discrete Fourier transform on the result of dividing the sum of the first frequency domain transformation result and the second frequency domain transformation result by a preset frequency domain adjustment factor to obtain an error wrapping number; The phase unwrapping application module is used to determine a true phase according to the predicted wrapping number and the error wrapping number, and calculate the actual elevation of the target position according to the true phase.

9. A terminal device, characterized in that, Including a processor, a memory, and a computer program stored in the memory and configured to be executed by the processor. When the processor executes the computer program, it implements the Poisson correction-based phase unwrapping method according to any one of claims 1-7.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium includes a stored computer program. When the computer program runs, it controls the device where the computer-readable storage medium is located to execute the Poisson correction-based phase unwrapping method according to any one of claims 1-7.

Images