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

Abstract

The present invention discloses a phase unwrapping method, apparatus, terminal device, and storage medium based on Poisson correction, relating to the field of interferometric measurement. The method comprises: obtaining an interferometric phase map to be processed at a target position, inputting the map into a winding number prediction model to generate a predicted winding number; calculating a winding number gradient based on the predicted winding number, and performing a lower-edge filtering correction to obtain a corrected winding number gradient; performing a discrete Fourier transform on the horizontal component of the corrected winding number gradient to obtain a first frequency domain transform result; performing a Laplace operator calculation on the corrected winding number gradient and then performing a discrete Fourier transform to obtain a second frequency domain transform result; performing an inverse discrete Fourier transform on the sum of the two frequency domain transform results divided by a preset frequency domain adjustment factor to obtain an error winding number; determining the true phase based on the predicted winding number and the error winding number, and then calculating the actual elevation of the target position. By implementing the present invention, the accuracy of the unwrapped phase can be improved.

Description

Technical Field

[0001] The present invention relates to the field of interferometric measurement, and in particular to a phase unwrapping method, apparatus, terminal equipment 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). It aims to restore the actual phase value from the wrapped phase and phase difference through specific mathematical methods. In the InSAR process, the accuracy of the phase unwrapping algorithm directly determines the accuracy of the generated Digital Elevation Model (DEM).

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

[0004] However, in the traditional least squares unwrapping algorithm, since the wrapped phase is affected by factors such as noise, 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, causing the unwrapping result to deviate from the true value. Summary of the Invention

[0005] The embodiments of the present invention provide a phase unwrapping method, apparatus, terminal device and storage medium based on Poisson correction, which can improve the accuracy of the unwrapped phase and solve the problem that the error of the traditional least squares unwrapping algorithm is evenly distributed to the global unwrapping result, causing the unwrapping result to deviate from the true value.

[0006] An embodiment of the present invention provides a phase unwrapping method based on Poisson correction, comprising:

[0007] Obtaining an interferometric phase image to be processed at a target position;

[0008] Inputting the interference phase image 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 image to be processed;

[0009] According to the predicted winding number, the winding number gradient is calculated, and the winding number gradient is corrected by lower edge filtering to obtain the corrected winding number gradient;

[0010] Performing discrete Fourier transform on the horizontal component of the corrected winding number gradient to obtain a first frequency domain transform result;

[0011] Performing Laplace operator calculation on the corrected winding number gradient and then performing discrete Fourier transform to obtain a second frequency domain transform result;

[0012] Performing an inverse discrete Fourier transform on the sum of the first frequency domain transform result and the second frequency domain transform result, divided by the preset frequency domain adjustment factor, to obtain an error winding number;

[0013] The true phase is determined according to the predicted winding number and the error winding number, and the actual elevation of the target position is calculated according to the true phase.

[0014] Furthermore, a predicted winding number is generated according to the interference phase image to be processed, including:

[0015] Perform image division on the interference phase image to be processed to obtain a number of block images;

[0016] For each block image, multi-scale pyramid pooling is performed on the block image to obtain a multi-scale feature map;

[0017] After performing frequency domain transformation on the multi-scale feature map, spectrum feature fusion is performed to obtain the frequency domain feature map;

[0018] Decode the frequency domain feature map to obtain the predicted winding number of the corresponding block image;

[0019] According to the predicted winding numbers of all the block images, the predicted winding number corresponding to the interference phase image to be processed is determined.

[0020] Furthermore, based on the predicted winding number, the winding number gradient is calculated, including:

[0021] According to the predicted winding number, the winding number gradient is calculated by the following formula:

[0022] ;

[0023] in, represents the winding number gradient, represents the gradient operator, , represents the predicted winding number, Indicates rounding.

[0024] Furthermore, the preset frequency domain adjustment factors include:

[0025] ;

[0026] in, Indicates the preset frequency domain adjustment factor, represents the length of the interference phase pattern, represents the width of the interference phase pattern, and Indicates the horizontal and vertical coordinate values ​​on a two-dimensional plane of the interference phase pattern.

[0027] Furthermore, the true phase is determined based on the predicted winding number and the error winding number, including:

[0028] determining the interference phase according to the interference phase image to be processed;

[0029] According to the interference phase, the predicted winding number and the error winding number, the true phase is calculated by the following formula:

[0030] ;

[0031] in, represents the true phase, represents the interference phase, represents the predicted winding number, Indicates the error winding number.

[0032] Furthermore, the actual elevation of the target position is calculated based on the true phase, including:

[0033] Obtain the antenna position height of the measurement 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;

[0034] The actual elevation of the target position is calculated based on 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.

[0035] Furthermore, the entanglement number prediction model is determined by:

[0036] Acquire a plurality of training samples; each training sample includes: an interference phase training graph and a corresponding actual winding number;

[0037] Inputting a number of training samples into a winding number prediction model to be trained, so that the winding number prediction model takes the interference phase training image as input and the predicted winding number corresponding to the interference phase training image as output for training, and during the training process, calculating a loss function according to the predicted winding number corresponding to the interference phase training image and the corresponding actual winding number, and adjusting the parameters of the winding number prediction model according to the loss function until the loss function converges, thereby obtaining a trained winding number prediction model;

[0038] Each training sample is determined by:

[0039] Get Gaussian random matrix;

[0040] Expand the Gaussian random matrix to obtain an expanded matrix;

[0041] Adjust the expanded matrix to the preset size and introduce Gaussian noise to generate a true phase training map;

[0042] Performing twisting processing on the real phase training image to obtain the interference phase training image;

[0043] According to the real phase training graph and the interference phase training graph, the actual winding number corresponding to the interference phase training graph is determined.

[0044] Based on the above method embodiment, the present invention provides a corresponding device embodiment, 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;

[0045] An interference phase image acquisition module, used to acquire an interference phase image to be processed at a target position;

[0046] A winding number prediction module is used to input the interferometric phase image to be processed into a trained winding number prediction model, so that the winding number prediction model generates a predicted winding number based on the interferometric phase image to be processed;

[0047] A winding number gradient correction module is used to 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 a corrected winding number gradient;

[0048] A first frequency domain transformation module is used to perform discrete Fourier transform on the horizontal component of the corrected winding number gradient to obtain a first frequency domain transformation result;

[0049] A 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 a second frequency domain transformation result;

[0050] an error winding number determining module, configured to perform an inverse discrete Fourier transform on the sum of the first frequency domain transform result and the second frequency domain transform result divided by a preset frequency domain adjustment factor to obtain an error winding number;

[0051] 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.

[0052] Based on the above-mentioned method embodiment, the present invention provides a corresponding terminal device 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.

[0053] Based on the above-mentioned method embodiment, the present invention provides a corresponding computer-readable storage medium embodiment, including: a stored computer program, which 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 when the computer program is running.

[0054] Compared with the prior art, the beneficial effects of the embodiment of this solution are:

[0055] The present invention obtains a pending interferometric phase image at a target location and inputs it into a trained winding number prediction model. The model, trained with extensive data, learns the complex relationship between phase and winding number and generates a predicted winding number. Based on the predicted winding number, a winding number gradient is calculated, reflecting the change in winding number between adjacent pixels. This winding number gradient is then corrected by applying a lower-edge filtering correction. This filtering suppresses minor gradient variations caused by noise while preserving gradients due to prediction errors, resulting in a corrected winding number gradient. Next, a discrete Fourier transform is performed on the horizontal component of the corrected winding number gradient to obtain the first frequency domain transform result. Low-frequency components correspond to slowly changing gradients, while high-frequency components correspond to locally abnormal gradients. The corrected winding number gradient is then subjected to a Laplace operator calculation and then a discrete Fourier transform to obtain the second frequency domain transform result. The Laplace operator calculation can highlight areas where the gradient changes are particularly fast or slow. A further discrete Fourier transform is performed to separate these characteristic information by frequency. The sum of the first and second frequency domain transform results divided by the preset frequency domain adjustment factor is then inversely discretely transformed to obtain the error winding number. These two sets of information at different frequencies are combined to form a multi-dimensional description of the error. Finally, the true phase is determined based on the predicted winding number and the error winding number. Based on the true phase, the actual elevation of the target location is calculated.

[0056] In summary, the present invention is based on the entanglement number predicted by the model, and decomposes the error entanglement number by performing a two-dimensional frequency domain analysis of the horizontal component discrete Fourier transform and the Laplace operator combined with the discrete Fourier transform on the entanglement number gradient. This effectively solves the problem that the error of the traditional least squares unwrapping algorithm is evenly distributed to the global unwrapping result, causing the unwrapping result to deviate from the true value, thereby improving the accuracy of the unwrapping phase. BRIEF DESCRIPTION OF THE DRAWINGS

[0057] Figure 11 is a flow chart of a phase unwrapping method based on Poisson correction provided by one embodiment of the present invention;

[0058] Figure 2 Schematic diagram of the geometric structure of spaceborne SAR interferometric altimetry provided by one embodiment of the present invention;

[0059] Figure 3 Schematic diagram of the network structure of the entanglement number prediction model provided by one embodiment of the present invention;

[0060] Figure 4 1 is a flow chart of a training process of an entanglement number prediction model provided by one embodiment of the present invention;

[0061] Figure 5 3 is a schematic structural diagram of a phase unwrapping device based on Poisson correction provided in one embodiment of the present invention. DETAILED DESCRIPTION

[0062] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. All other embodiments obtained by ordinary technicians in this field based on the embodiments of the present invention without making any creative efforts shall fall within the scope of protection of the present invention.

[0063] 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 understood as indicating or implying relative importance or implicitly indicating the number of the indicated technical features.

[0064] like Figure 1 As shown, in order to solve the problem that the error of the traditional least squares unwrapping algorithm is evenly distributed to the global unwrapping result, causing the unwrapping result to deviate from the true value, an embodiment of the present invention provides a phase unwrapping method based on Poisson correction, which includes at least the following steps:

[0065] Step S1: obtaining an interferometric phase image to be processed at a target position;

[0066] For step S1, Figure 2 Figure 2 shows the schematic diagram of the geometric model of spaceborne SAR interferometric altimetry. The satellite transmits and receives radar waves at different locations S1 and S2. The radar wave transmitted from S1 is reflected by the target location P and then received. The radar wave transmitted from S2 is also reflected by P and then received. Due to the spatial position difference of the target location P relative to S1 and S2, as well as the different propagation path lengths of the radar waves, the two radar waves produce a phase difference when they are received. This phase difference information constitutes the interferometric phase pattern of the target location.

[0067] It should be noted that in the process of extracting phase information from the interferometric phase pattern, the absolute phase value will be wrapped in the range of [-π,π) or [0, 2π), which is called the wrapped phase value. The task of phase unwrapping is to determine the wrapping number from the wrapped phase and thus recover the true phase value. Wrapping phase value The true phase value The relationship can be expressed by the following formula:

[0068] ;

[0069] in, represents the true phase value, represents the winding phase value, represents the winding number, which is an integer function representing the integer multiple of 2π corresponding to the true phase value.

[0070] Step S2: inputting the interference phase image 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 image to be processed;

[0071] In a preferred embodiment, generating a predicted winding number according to the interference phase image to be processed includes:

[0072] Perform image division on the interference phase image to be processed to obtain a number of block images;

[0073] For each block image, multi-scale pyramid pooling is performed on the block image to obtain a multi-scale feature map;

[0074] After performing frequency domain transformation on the multi-scale feature map, spectrum feature fusion is performed to obtain the frequency domain feature map;

[0075] Decode the frequency domain feature map to obtain the predicted winding number of the corresponding block image;

[0076] According to the predicted winding numbers of all the block images, the predicted winding number corresponding to the interference phase image to be processed is determined.

[0077] For step S2, the interference phase image to be processed obtained in step S1 is input into the trained winding number prediction model. The winding number prediction model is trained to reflect the complex relationship between phase and winding number, and generates a predicted winding number based on the interference phase image to be processed.

[0078] Specifically, if Figure 3As shown, the entanglement number prediction model uses a PUformer network structure. First, before encoding, the interferometric phase image is segmented into several non-overlapping 4×4 pixel patches, each containing 16 pixel phase values. If the original interferometric phase image is of size H×W, the resulting image patches total (H / 4)×(W / 4). Each patch serves as the minimum feature unit for subsequent encoding in the entanglement number prediction model, effectively preserving local structural information while reducing computational effort.

[0079] PUFormer-based phase unwrapping is divided into two stages. In the first stage, for each block of the image, the MiT (Mix Vision Transformer) encoder, built into the entanglement number prediction model, extracts and encodes features. Specifically, the MiT encoder consists of four stages (Stage 1 to Stage 4). Each stage is a transformer block composed of multiple Transformer blocks. Downsampling operations are performed between stages to achieve spatial compression and channel expansion. The feature maps output by each stage correspond to 1 / 4, 1 / 8, 1 / 16, and 1 / 32 resolutions of the interference phase image, respectively, forming multi-scale feature maps that capture multi-level phase information from detailed to global.

[0080] In the second stage, the point-by-point convolutional decoder built into the entanglement prediction model performs a frequency domain transform on the feature maps at each scale, obtaining a frequency domain representation of the feature maps at each scale. The spectrum is then divided into low-frequency and high-frequency components. Low-frequency features preserve large-scale structural information in the image, such as background areas and contour lines, and have good global semantic consistency. High-frequency features correspond to sudden changes and detailed changes at the boundaries of the entanglement number, manifesting as local features such as phase jumps, making them suitable for edge perception and sharpening. The light blue arrows in the figure represent the gradual upsampling operation, which upsamples the low-resolution features (1 / 32, 1 / 16, and 1 / 8) layer by layer until they are finally unified to 1 / 4 resolution. After the spatial dimensions of the feature maps of each scale are consistent, they are fused through the frequency fusion block in a concatenated manner, and the fused frequency domain feature maps are input into the built-in convolution fusion layer for feature integration, normalization, and introduction of nonlinear transformation. The output processed and transformed feature maps are input into the pixelwise convolution layer for refined feature extraction, and finally the predicted entanglement number is output. .

[0081] Finally, the predicted winding number of each block image reflects the phase winding situation of the local area. The predicted winding numbers of all block images are combined to obtain the predicted winding number corresponding to the interference phase image to be processed.

[0082] This paper introduces a frequency-aware point-by-point convolutional decoder. This module can perceive gradient jumps in the interference phase and obtain relevant high- and low-frequency information, improving the semantic segmentation model's accuracy in predicting gradient jump regions. Furthermore, the convolutional fusion layer of this paper is a lightweight decoder composed of a convolution (Conv), a batch normalization (BN), and an activation function (ReLU). This lightweight point-by-point convolution replaces the widely used decoder with larger parameters, reducing model parameters while maintaining prediction accuracy.

[0083] Next, the training of the entanglement number prediction model is described in detail:

[0084] like Figure 4 As shown, the training of the above-mentioned entanglement number prediction model includes the following steps:

[0085] Step S21: obtaining a plurality of training samples; each training sample includes: an interference phase training graph and a corresponding actual winding number;

[0086] Each training sample is determined by:

[0087] Get Gaussian random matrix;

[0088] Expand the Gaussian random matrix to obtain an expanded matrix;

[0089] Adjust the expanded matrix to the preset size and introduce Gaussian noise to generate a true phase training map;

[0090] Performing twisting processing on the real phase training image to obtain the interference phase training image;

[0091] Determining the actual winding number corresponding to the interference phase training map according to the real phase training map and the interference phase training map;

[0092] In step S21, training samples are obtained for training the winding number prediction model. Each training sample includes an interferometric phase training image and a corresponding actual winding number. Deep learning algorithm training requires a large amount of data to determine the mapping between the input winding phase and the output winding number. However, obtaining large-scale, realistic winding phase images and their corresponding precise winding numbers is difficult. Therefore, training samples are generated using the following method.

[0093] First, a Gaussian random matrix is ​​obtained, with a matrix size of 2*2~10*10 and a value range of 10~150. Next, the Gaussian random matrix is ​​expanded using a mathematical fractal algorithm. In this embodiment, the diamond-square method is used to process and augment the random matrix 2~5 times to generate a more complex true phase matrix. The expanded matrix is ​​adjusted to a preset size, in this embodiment, it is adjusted to 256*256. The smaller matrix is ​​interpolated to obtain a smoother true phase, and the larger matrix is ​​downsampled to obtain a rougher true phase, thereby obtaining phase data with both smooth and rough characteristics. Then, Gaussian noise with a standard deviation between 0 and 0.3 is introduced to simulate the interference of the phase in the real world and generate a true phase training map.

[0094] Finally, the generated real phase training map is twisted using the following formula to obtain the interference phase training map:

[0095] ;

[0096] Then, based on the real phase training graph and the interference phase training graph, the actual winding number corresponding to the interference phase training graph is determined by the following formula:

[0097] ;

[0098] in, represents the interferometric phase training graph, represents the true phase training graph, Remaining operation, Indicates the actual winding number corresponding to the interferometric phase training diagram.

[0099] Step S22: inputting a plurality of training samples into a winding number prediction model to be trained, so that the winding number prediction model takes the interference phase training map as input and the predicted winding number corresponding to the interference phase training map as output for training, and during the training process, calculating a loss function based on the predicted winding number corresponding to the interference phase training map and the corresponding actual winding number, and adjusting the parameters of the winding number prediction model based on the loss function until the loss function converges, thereby obtaining a trained winding number prediction model;

[0100] In step S22, the training samples obtained in step S21 are input into the winding number prediction model to be trained. During training, the winding number prediction model uses the interferometer phase training image as input. The model performs a series of feature extraction, encoding, and decoding on the input interferometer phase training image, and outputs a predicted winding number for the interferometer phase training image. During training, the predicted winding number output by the model is compared with the actual winding number corresponding to the training sample. The mean squared difference (MSE) of the squares of the differences between the predicted and actual winding numbers is calculated and used as the loss function value. The loss function is used to evaluate the difference between the model's predicted results and the actual results. Based on the calculated loss function value, the parameters of the winding number prediction model are adjusted using gradient descent. This training process is repeated continuously, and the loss function value gradually decreases as training progresses. When the loss function value no longer decreases significantly and reaches a relatively stable state, which is called convergence, it indicates that the model has learned enough knowledge and the trained winding number prediction model is obtained.

[0101] Step S3: 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 a corrected winding number gradient;

[0102] In a preferred embodiment, the winding number gradient is calculated based on the predicted winding number, including:

[0103] According to the predicted winding number, the winding number gradient is calculated by the following formula:

[0104] ;

[0105] in, represents the winding number gradient, represents the gradient operator, , represents the predicted winding number, Indicates rounding.

[0106] For step S3, the predicted winding number is obtained through step S2. Then, in order to detect the rate of change of the phase, the winding number gradient is calculated by the formula , which reflects the change in the winding number in space. The rounding operation means converting the continuous phase change into discrete integer jumps, that is, determining the winding number increment at each position.

[0107] When the phase gradient changes drastically, a small gradient change may cause an error in the estimation of the winding number gradient field. In order to meet the correction needs, a phase filter is used. Perform down-edge filtering correction on the winding number gradient to obtain the corrected winding number gradient :

[0108] ;

[0109] By using the phase filter, small gradient changes due to noise will be suppressed, and the gradient caused by the prediction error will be preserved.

[0110] Step S4: performing discrete Fourier transform on the horizontal component of the corrected winding number gradient to obtain a first frequency domain transform result;

[0111] For step S4, the corrected winding number gradient Partition into the horizontal component of the modified winding number gradient and the vertical component of the modified winding number gradient , horizontal component It reflects the change of winding number gradient in the horizontal direction of the two-dimensional interference phase pattern. It reflects the change of the winding number gradient in the vertical direction of the two-dimensional interference phase pattern.

[0112] Next, the horizontal component of the modified winding number gradient is Perform discrete Fourier transform (DCT) to obtain the first frequency domain transform result:

[0113] ;

[0114] in, represents the first frequency domain transformation result, represents the horizontal component of the modified winding number gradient, Denotes discrete Fourier transform processing.

[0115] The interference phase pattern is converted from the spatial domain to the frequency domain through discrete Fourier transform for analysis to understand the horizontal component The information of different frequency components in the .

[0116] Step S5: performing Laplace operator calculation on the corrected winding number gradient and then performing discrete Fourier transform to obtain a second frequency domain transform result;

[0117] In step S5, similarly, the modified winding number gradient is subjected to Laplace operator calculation and then discrete Fourier transform to obtain the second frequency domain transform result:

[0118] ;

[0119] in, represents the second frequency domain transform result, represents the Laplace operator, represents the corrected winding number gradient.

[0120] It should be noted that the present invention obtains part of the frequency domain information by performing DCT transform on the horizontal component in step S4. The winding number gradient after Laplace operator processing in step S5 is then subjected to DCT transform to obtain frequency domain information from another perspective. The two different frequency domain transform results complement each other and can more comprehensively describe the frequency characteristics of the winding number gradient.

[0121] Step S6: performing an inverse discrete Fourier transform on the sum of the first frequency domain transform result and the second frequency domain transform result divided by the preset frequency domain adjustment factor to obtain an error winding number;

[0122] In a preferred embodiment, the preset frequency domain adjustment factor includes:

[0123] ;

[0124] in, Indicates the preset frequency domain adjustment factor, represents the length of the interference phase pattern, represents the width of the interference phase pattern, and Indicates the horizontal and vertical coordinate values ​​on a two-dimensional plane of the interference phase pattern.

[0125] In step S6, the sum of the first frequency domain transform result and the second frequency domain transform result divided by the preset frequency domain adjustment factor is subjected to inverse discrete Fourier transform to obtain the error winding number using the following formula:

[0126] ;

[0127] ;

[0128] ;

[0129] in, represents 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, Indicates the preset frequency domain adjustment factor, represents the length of the interference phase pattern, represents the width of the interference phase pattern, and Represents the horizontal and vertical coordinate values ​​on the two-dimensional plane of the interference phase pattern, represents the error winding number, Indicates the inverse discrete Fourier transform process.

[0130] It should be noted that the present invention first converts the first frequency domain transformation result obtained in step S4 into and the second frequency domain transformation result obtained in step S5 The frequency domain information acquired from different angles is added together to fuse the information, taking into account both the horizontal gradient characteristics and the frequency domain characteristics after Laplace operator enhancement. Furthermore, the Fourier transform converts the spatial Laplace operator into an algebraic operation in the frequency domain, simplifying the solution of the equation. After adding the two, the numerical ranges and magnitudes of the different frequency components may vary significantly. Dividing the sum by the preset frequency domain adjustment factor param normalizes the fused frequency domain information, bringing the different frequency components into a relatively uniform and appropriate numerical range. Finally, an inverse discrete Fourier transform (IDCT) is performed to convert the frequency domain information back to the spatial domain, yielding the error winding number.

[0131] This process is similar to solving the Poisson equation to recover phase information. In the phase unwrapping problem, the Poisson equation describes the relationship between the phase gradient and the true phase. By performing an inverse transform on the frequency-domain processed results, 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.

[0132] 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.

[0133] In a preferred embodiment, determining the true phase according to the predicted winding number and the error winding number includes:

[0134] determining the interference phase according to the interference phase image to be processed;

[0135] According to the interference phase, the predicted winding number and the error winding number, the true phase is calculated by the following formula:

[0136] ;

[0137] in, represents the true phase, represents the interference phase, represents the predicted winding number, Indicates the error winding number.

[0138] In a preferred embodiment, calculating the actual elevation of the target position according to the true phase includes:

[0139] Obtain the antenna position height of the measurement 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;

[0140] The actual elevation of the target position is calculated based on 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.

[0141] In step S7, since the interference phase image is obtained using interferometry technology, based on the principle of wave interference, when two coherent waves meet, interference occurs due to factors such as path difference. The interference phase can be obtained by directly reading the interference fringe distribution in the interference phase image to be processed.

[0142] The present invention has obtained the predicted winding number through the winding number prediction model of step S2 and determine the error winding number through step S6 Then, by calculating the formula Compensate for missing or incorrect phase values ​​caused by phase wrapping, thereby restoring the true phase.

[0143] The unwrapped true phase is applied to spaceborne SAR interferometric altimetry, as Figure 2 As shown in Figure 1, the antenna positions include the main image antenna position and the secondary image antenna position. Both transmit radar waves to the same target position P. The two radar waves propagate to the target from different paths and are reflected back to the satellite. Due to the existence of the propagation path, interference occurs when the two coherent radar waves meet, and the interferometric phase diagram to be processed is obtained. 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 inclination 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, and P0 is a point on the reference plane. P0 has the same main image slant distance as P, and the distance from the main image antenna position S1 to P is equal to that to P0. Its elevation is 0, and θ0 is the local incident angle of P0.

[0144] Obtain the main image antenna position height H, the local incident angle θ of the main image antenna position, the interference 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 real phase , the actual elevation h of the target position can be calculated using the following formula:

[0145] ;

[0146] in, Indicates the actual elevation of the target location, represents the wavelength of the radar wave, Indicates the distance from the main image antenna position to the target position, represents the local angle of incidence at the main image antenna position, represents the interference baseline, represents the baseline tilt angle, Indicates the true phase.

[0147] The present invention uses a deep learning model to preliminarily predict the winding number. Based on winding number gradient analysis, the error winding number is determined by using a relationship between the phase gradient and the true phase similar to that in the Poisson equation, thereby correcting the winding number and restoring the true phase. Compared with traditional methods, the untangling success rate is higher when faced with extreme conditions such as large terrain undulations and dense stripes. At the same time, when encountering data with a distribution that differs from the training set, general deep learning methods are susceptible to the limitations of the training data, resulting in inaccurate predictions of the winding number. Based on this, the present invention improves the untangling effect by correcting the winding number, making the untangling effect better than traditional methods and deep learning methods.

[0148] like Figure 5 As shown, based on the above method embodiment, a corresponding device embodiment is provided;

[0149] An embodiment of the present invention provides a phase unwrapping device based on Poisson correction, comprising: 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;

[0150] An interference phase image acquisition module, used to acquire an interference phase image to be processed at a target position;

[0151] A winding number prediction module is used to input the interferometric phase image to be processed into a trained winding number prediction model, so that the winding number prediction model generates a predicted winding number based on the interferometric phase image to be processed;

[0152] A winding number gradient correction module is used to 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 a corrected winding number gradient;

[0153] A first frequency domain transformation module is used to perform discrete Fourier transform on the horizontal component of the corrected winding number gradient to obtain a first frequency domain transformation result;

[0154] A 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 a second frequency domain transformation result;

[0155] an error winding number determining module, configured to perform an inverse discrete Fourier transform on the sum of the first frequency domain transform result and the second frequency domain transform result divided by a preset frequency domain adjustment factor to obtain an error winding number;

[0156] 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.

[0157] It can be understood that the above-mentioned device embodiment corresponds to the method embodiment of the present invention, which can implement the phase unwrapping method based on Poisson correction provided by any of the above-mentioned method embodiments of the present invention.

[0158] It should be noted that the device embodiments described above are merely illustrative, and some or all of the modules may be selected according to actual needs to achieve the purpose of the present embodiment. Furthermore, in the drawings of the device embodiments provided by the present invention, the connection relationship between modules indicates that they have a communication connection, which may be implemented as one or more communication buses or signal lines. Those skilled in the art can understand and implement the present invention without inventive effort.

[0159] Based on the above-mentioned embodiment 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 of any embodiment of the present invention is implemented.

[0160] For example, in this embodiment, the computer program may be divided into one or more modules, which are stored in the memory and executed by the processor to implement the present invention. The one or more module elements may be a series of computer program instruction segments capable of performing specific functions, and the instruction segments are used to describe the execution process of the computer program in the terminal device.

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

[0162] The processor may be a central processing unit (CPU), other general-purpose processors, digital signal processors (DSP), application-specific integrated circuits (ASIC), field-programmable gate arrays (FPGA) or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. The general-purpose processor may be a microprocessor or any conventional processor, etc. The processor is the control center of the terminal device, connecting various parts of the entire terminal device using various interfaces and lines.

[0163] Based on the above method embodiment, another embodiment is provided: another embodiment of the present invention provides a computer-readable storage medium, including a stored computer program, wherein when the computer program is running, the device where the computer-readable storage medium is located is controlled to execute the phase unwrapping method based on Poisson correction described in any one of the above method embodiments of the present invention.

[0164] If the module / unit integrated into the Poisson correction-based phase unwrapping device / terminal equipment is implemented as a software functional unit and sold or used as a standalone product, it can be stored in a computer-readable storage medium. Based on this understanding, the present invention can also implement all or part of the process steps in the above-mentioned method embodiments by instructing the relevant hardware through a computer program. The computer program can be stored in a computer-readable storage medium. When executed by a processor, the computer program can implement the steps of each of the above-mentioned method embodiments. The computer program includes computer program code, which can be in source code form, object code form, executable file, or some intermediate form. The computer-readable medium can include any entity or device capable of carrying the computer program code, a recording medium, a USB flash drive, a mobile hard drive, a magnetic disk, an optical disk, a computer memory, a read-only memory (ROM), a random access memory (RAM), an electrical carrier signal, a telecommunications signal, and a software distribution medium.

[0165] The above is a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principles of the present invention. These improvements and modifications are also considered to be within the scope of protection of the present invention.

Claims

1. A phase unwrapping method based on Poisson correction, characterized in that: include: Obtaining an interferometric phase image to be processed at a target position; Inputting the interferometric phase image to be processed into a trained winding number prediction model so that the winding number prediction model performs image segmentation on the interferometric phase image to obtain a plurality of block images; performing multi-scale pyramid pooling on each block image to obtain a multi-scale feature map; performing frequency domain transformation on the multi-scale feature map and then performing spectrum feature fusion to obtain a frequency domain feature map; decoding the frequency domain feature map to obtain a predicted winding number of the corresponding block image; and determining the predicted winding number corresponding to the interferometric phase image to be processed based on the predicted winding numbers of all block images; According to the predicted winding number, a winding number gradient is calculated, and the winding number gradient is corrected by lower edge filtering to obtain a corrected winding number gradient; Performing discrete Fourier transform on the horizontal component of the corrected winding number gradient to obtain a first frequency domain transform result; Performing Laplace operator calculation on the corrected winding number gradient and then performing discrete Fourier transform to obtain a second frequency domain transform result; Performing an inverse discrete Fourier transform on the sum of the first frequency domain transform result and the second frequency domain transform result, divided by the preset frequency domain adjustment factor, to obtain an error winding number; The preset frequency domain adjustment factor includes: ； in, Indicates the preset frequency domain adjustment factor, represents the length of the interference phase pattern, represents the width of the interference phase pattern, and Represents the horizontal coordinate value and the vertical coordinate value on the two-dimensional plane of the interference phase pattern; The true phase is determined according to the predicted winding number and the error winding number, and the actual elevation of the target position is calculated according to the true phase.

2. The phase unwrapping method based on Poisson correction according to claim 1, characterized in that: According to the predicted winding number, a winding number gradient is calculated, including: According to the predicted winding number, the winding number gradient is calculated by the following formula: ； in, represents the winding number gradient, represents the gradient operator, , represents the interference phase, represents the predicted winding number, Indicates rounding.

3. The phase unwrapping method based on Poisson correction according to claim 1, characterized in that: Determine the true phase based on the predicted winding number and the error winding number, including: determining the interference phase according to the interference phase image to be processed; According to the interference phase, the predicted winding number and the error winding number, the true phase is calculated by the following formula: ； in, represents the true phase, represents the interference phase, represents the predicted winding number, Indicates the error winding number.

4. The phase unwrapping method based on Poisson correction according to claim 1, characterized in that: Calculating the actual elevation of the target position according to the true phase includes: Obtain the antenna position height of the measurement 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; The actual elevation of the target position is calculated based on 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.

5. 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: Acquire a plurality of training samples; each training sample includes: an interference phase training graph and a corresponding actual winding number; Inputting a number of training samples into a winding number prediction model to be trained, so that the winding number prediction model takes the interference phase training image as input and the predicted winding number corresponding to the interference phase training image as output for training, and during the training process, calculating a loss function according to the predicted winding number corresponding to the interference phase training image and the corresponding actual winding number, and adjusting the parameters of the winding number prediction model according to the loss function until the loss function converges, thereby obtaining a trained winding number prediction model; Each training sample is determined by: Get Gaussian random matrix; Expanding the Gaussian random matrix to obtain an expanded matrix; Adjust the expanded matrix to the preset size and introduce Gaussian noise to generate a true phase training map; Performing twisting processing on the real phase training image to obtain the interference phase training image; According to the real phase training graph and the interference phase training graph, the actual winding number corresponding to the interference phase training graph is determined.

6. A phase unwrapping device based on Poisson correction, characterized in that: include: Interference phase image acquisition module, winding number prediction module, winding number gradient correction module, first frequency domain transformation module, second frequency domain transformation module, error winding number determination module and phase unwrapping application module; The interference phase image acquisition module is used to acquire the interference phase image to be processed at the target position; The winding number prediction module is configured to input the interference phase image to be processed into a trained winding number prediction model, so that the winding number prediction model performs image segmentation on the interference phase image to be processed to obtain a plurality of block images; for each block image, perform multi-scale pyramid pooling on the block image to obtain a multi-scale feature map; perform frequency domain transformation on the multi-scale feature map, and then perform spectrum feature fusion to obtain a frequency domain feature map; decode the frequency domain feature map to obtain a predicted winding number of the corresponding block image; and determine the predicted winding number corresponding to the interference phase image to be processed based on the predicted winding numbers of all block images; The winding number gradient correction module is used to 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 a 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 a 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 a second frequency domain transformation result; The error winding number determining module is configured to perform an inverse discrete Fourier transform on the sum of the first frequency domain transform result and the second frequency domain transform result divided by a preset frequency domain adjustment factor to obtain an error winding number; The preset frequency domain adjustment factor includes: ； in, Indicates the preset frequency domain adjustment factor, represents the length of the interference phase pattern, represents the width of the interference phase pattern, and Represents the horizontal coordinate value and the vertical coordinate value on the two-dimensional plane of the interference phase pattern; 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.

7. A terminal device, characterized in that: The method comprises a processor, a memory, and a computer program stored in the memory and configured to be executed by the processor, wherein when the processor executes the computer program, the phase unwrapping method based on Poisson correction according to any one of claims 1 to 5 is implemented.

8. A computer-readable storage medium, characterized in that The computer-readable storage medium includes a stored computer program, wherein when the computer program is executed, the device where the computer-readable storage medium is located is controlled to execute the phase unwrapping method based on Poisson correction according to any one of claims 1 to 5.

Images