An SAR interferogram generation method based on L 1 / 2 -norm regularization

Abstract

The present invention discloses a method for generating an SAR interferogram based on L<subgt;1 / 2< / subgt>-norm regularization, specifically as follows: Step 1: Use SAR to obtain the main and secondary SAR images of the same target scene, and obtain the phase of the secondary SAR image by taking the phase value of the complex matrix of the secondary image; Step 2: Based on the interferometric phase of the main and secondary SAR images, and establish a generation model of the SAR interferogram by solving the L<subgt;1 / 2< / subgt>-norm regularization problem; Step 3: Use the threshold iteration algorithm to solve the L<subgt;1 / 2< / subgt>-norm regularization problem to obtain a wavelet coefficient matrix related to the restored interferogram; Step 4: Perform an inverse discrete wavelet transform on the wavelet coefficient matrix obtained in Step 3, and take the transformed angle value as the restored interferogram. The interferogram reconstructed by the present invention has better performance and smaller phase noise; the quality of the SAR interferogram obtained by the present invention is higher.

Description

Technical Field

[0001] The present invention belongs to the technical fields of interferometric synthetic aperture radar imaging and sparse signal processing. Background Art

[0002] SAR is a high-resolution imaging radar. Compared with optical imaging methods, it can obtain images of large scenes and wide swaths all day and all weather, and is an important research direction of current microwave imaging technology. Traditional SAR imaging technology is limited to two-dimensional (2D) imaging, while interferometric SAR (InSAR) imaging technology breaks through this limitation. It uses the phase difference between two SAR complex images to restore the topographic map of the target scene, thereby obtaining the height information of the observed target and the digital elevation model (DEM). Therefore, it is widely used in fields such as topographic surveying and terrain exploration. High-quality DEM has high requirements for the processing accuracy of interferograms. Currently, typical phase noise filtering techniques, such as the Lee filter, Goldstein filter, etc., will all lead to a reduction in the spatial resolution of the image. Summary of the Invention

[0003] Object of the Invention: To solve the problems existing in the above-mentioned prior art, the present invention provides a method for generating an SAR interferogram based on L 1 / 2 -norm regularization.

[0004] Technical Solution: The present invention provides a method for generating an SAR interferogram based on L 1 / 2 -norm regularization, specifically including the following steps:

[0005] Step 1: Use SAR to obtain the SAR main image and the SAR secondary image of the same target scene, and obtain the phase of the SAR secondary image by taking the phase value of the secondary image complex matrix;

[0006] Step 2: Based on the interferometric phase of the SAR main and secondary images, and by solving the L 1 / 2 -norm regularization problem, establish a generation model for the SAR interferogram;

[0007] Step 3: Use the threshold iteration algorithm to solve the L 1 / 2 -norm regularization problem to obtain the wavelet coefficient matrix related to the restored SAR interferogram

[0008] Step 4: Perform an inverse discrete wavelet transform on the result obtained in Step 3 , and take the transformed angle value as the restored interferogram.

[0009] Further, the phase y of the SAR sub-image in step 1 is:

[0010]

[0011] θ m = exp{j(φ m - φ flat )}

[0012] u = exp(-jφ topo )

[0013] n = θ m u{exp(-jφ noise ) - 1}

[0014] where φ m is the phase of the SAR main image y m and φ flat is the flat-earth phase, and φ topo is the terrain phase; φ noise = φ s2 - φ s1 and φ s1 and φ s2 represent the scattering phases of y m and y s respectively, y s represents the SAR sub-image, exp{.} is the exponential power of e, and j represents the imaginary part.

[0015] Further, the generation model of the SAR interferogram in step 2 is:

[0016]

[0017] where represents rewriting the vector form of y into the matrix form Y, and Y is a two-dimensional SAR sub-image matrix with an amplitude of 1; represents rewriting the vector form of θ m into the matrix form Θ m , and Θ m is a matrix containing the phase of the SAR main image and the flat-earth phase; R(.) is the inverse discrete wavelet transform; represents the Hadamard product; X = W(U), represents rewriting the vector form of u into the matrix form U, and U represents a matrix containing the interferometric phase of the reconstructed SAR main and sub-images; represents rewriting the vector form of n into the matrix form N, and N represents the noise matrix, β represents the regularization parameter, represents the Fibonacci norm, argmin{.} represents taking the minimum value, and W(.) is the discrete wavelet transform.

[0018] Furthermore, step 3 is specifically as follows:

[0019] Step 3.1: Set the initial value X of the wavelet coefficient (0) = 0; the wavelet coefficient is the wavelet coefficient after the DWT transform of the interference phase of the SAR main and secondary images;

[0020] Step 3.2: Calculate the residual estimate value during the t-th iterative calculation:

[0021]

[0022] X (t) represents the wavelet coefficient during the t-th iterative calculation;

[0023] Step 3.3: Update the wavelet coefficient during the t-th iterative calculation:

[0024] S (t) = X (t) + ΔX (t)

[0025] S (t) represents the updated wavelet coefficient during the t-th iterative calculation;

[0026] Step 3.4: Update the regularization parameter β during the t-th iterative calculation (t) :

[0027] β (t) = |S (t) | K+1 / μ

[0028] where μ is a preset parameter, |S (t) | represents the amplitude component of S (t) | (t) represents the (K + 1)-th component after arranging the amplitude components of S K+1 in descending order; (t)

[0029] Step 3.5: Update the wavelet coefficient X for the next iterative calculation according to the value of S (t) : (t+1) :

[0030] F(S (t) , μβ) = f(g, μβ)

[0031]

[0032] g = S (t)

[0033] where μ is an iterative parameter, and the value of f(g, μβ) is used as X ( t+1) value;

[0034] Step 3.6: Calculate the error of the (t + 1)-th iteration calculation:

[0035] Resi = ||X (t+1) - X (t) || F

[0036] Step 3.7: If Resi is greater than the preset threshold ε and the number of iterations t is less than or equal to the preset maximum number of iterations, then increment the number of iterations by t + 1 and go to Step 3.2; otherwise, stop the calculation.

[0037] Beneficial effects: The present invention can perform interferogram reconstruction based on MF images. The interferogram reconstructed using the present invention has better performance and smaller phase noise; the SAR interferogram obtained using the present invention has higher quality. Description of the Drawings

[0038] Figure 1 is the implementation flowchart of the SAR interferogram generation method based on L 1 / 2 norm regularization of the present invention;

[0039] Figure 2 is the geometric position relationship diagram between the target scene and the antenna;

[0040] Figure 3 is the iterative implementation flowchart of the SAR interferogram generation method based on L 1 / 2 norm regularization of the present invention;

[0041] Figure 4 is the simulation result diagram based on MF images, based on L1 norm regularization, and based on L 1 / 2 norm regularization; where (a) is the simulation result diagram of SAR interferogram generation based on MF images when the noise standard deviation is 1.0; (b) is the simulation result diagram of SAR interferogram generation based on L1 norm regularization; (c) is the simulation result diagram of SAR interferogram generation based on L 1 / 2 norm regularization. Detailed Embodiments

[0042] The drawings constituting a part of the present invention are used to provide a further understanding of the present invention. The schematic embodiments of the present invention and their descriptions are used to explain the present invention and do not constitute an improper limitation of the present invention.

[0043] As Figure 1 shown, this embodiment proposes a SAR interferogram generation method based on L 1 / 2 norm regularization, and the specific implementation steps are as follows:

[0044] Step S1: SAR Interferometric Phase Analysis

[0045] Let the SAR platform obtain the primary and secondary SAR images of the target scene as y m and y s respectively, and their phases be φ m and φ s respectively. According to the InSAR imaging geometric relationship, the phases of the primary and secondary images can be expressed as

[0046]

[0047] where R1 and R2 respectively represent the slant ranges of the two antennas from the center of the target scene, and φ s1 and φ s2 represent the scattering phases of y m and y s respectively, and λ c represents the wavelength of the received signal. The geometric position relationship between the target scene and the antennas is as shown in Figure 2 , where B is the distance between the two antennas and H is the height.

[0048] According to the InSAR theory, the phase difference between the primary and secondary images can be expressed as

[0049]

[0050] where ΔR = R2 - R1, and the noise phase φ noise is the difference between the scattering phases φ s1 and φ s2 , and can be expressed as

[0051] φ noise = φ s2 - φ s1 (4)

[0052] Since the expression on the right side of equation (3) can be composed of the flat-earth phase φ flat and the terrain phase φ topo , then the phase difference between the primary and secondary images can be expressed as the sum of three types of phases

[0053] φ m - φ s = φ flat + φ topo + φ noise (5)

[0054] It can be seen from equation (5) that the terrain phase φ topo can be obtained by removing the flat-earth phase φ flat and the noise phase φ noise . Therefore, to obtain a high-quality interferometric phase map, it is necessary to suppress the flat-earth phase and the noise phase as much as possible.

[0055] Step S2: Based on the L 1 / 2 Construction of the SAR interferogram generation model based on the norm regularization

[0056] Write the phase of the sub-image y s in the following form:

[0057]

[0058] exp{.} is the exponential power of e, and j represents the imaginary part.

[0059] The amplitude of = 1, |y s | represents the amplitude of y s .

[0060] Rewrite Equation 6 in the following form:

[0061] y = θ m u + n (7)

[0062]

[0063] θ m = exp{j(φ m - φ flat )} (9)

[0064] u = exp(-jφ topo ) (10)

[0065] n = θ m u{exp(-jφ noise ) - 1} (11)

[0066] Let denote the vector form of y rewritten as the matrix form Y, where Y is a two-dimensional SAR sub-image matrix with an amplitude of 1; denote the vector form of θ m rewritten as the matrix form Θ m , and Θ m is a matrix containing the phase of the SAR main image and the flat-earth phase; denote the vector form of u rewritten as the matrix form U, where U represents a matrix containing the reconstructed interferometric phase, denote the vector form of n rewritten as the matrix form N, where N represents the noise matrix;

[0067] Therefore, Equation 7 can be expressed as

[0068]

[0069] where denotes the Hadamard product operator, and R(.) represents the inverse discrete wavelet transform.

[0070] For Equation (12), first, the matrix U is sparsely represented using DWT. Let W(·) and R(·) denote DWT (after discrete wavelet transform) and IDWT (inverse discrete wavelet transform, abbreviated as IDWT), respectively. Then, X = W(U), and the matrix X represents the wavelet coefficients after the discrete wavelet transform of the matrix U. Thus, Equation (12) can be expressed as

[0071]

[0072] By solving the L 1 / 2 -norm regularization problem to reconstruct the interference image.

[0073]

[0074] is the wavelet coefficient matrix related to the restored interference image, argmin{.} represents taking the minimum value, β represents regularization, denotes the Fibonacci norm.

[0075] Step S3: Threshold iteration

[0076] For the L 1 / 2 -norm regularization problem in Formula (14), this embodiment uses a threshold iteration algorithm to implement the reconstruction of the interference pattern. The threshold iteration algorithm inputs the secondary image Y and the phase matrix Θ containing the phase of the primary image and the flat-earth phase m . Let the initial value of the wavelet coefficient be X (0) = 0; the wavelet coefficient is the wavelet coefficient after the DWT transform of the interference phase of the SAR primary and secondary images, the iteration parameter is μ, usually set as a constant, the error parameter is ε, and the maximum number of iteration steps is T max . When the conditions t ≤ T max and the iteration error Resi > ε are satisfied, as shown in Figure 3 , the following steps are executed.

[0077] Step S31: Residual estimation

[0078]

[0079] The symbol (·) H represents the conjugate transpose operator, and N (t) represents the wavelet coefficient at the t-th iteration calculation.

[0080] Step S32: Update the wavelet coefficient at the t-th iteration calculation:

[0081] S(t) = X (t) + ΔX (t) (16)

[0082] where S (t) represents the wavelet coefficient at the t-th iteration calculation after update.

[0083] Step S33: Update the regularization parameter β at the t-th iteration calculation (t) :

[0084] β (t) = |S (t) | K+1 / μ (17)

[0085] where |S (t) | represents the amplitude component of S (t) ; |S (t) | k+1 represents the (K + 1)-th component after arranging the amplitude components of S (t) in descending order;

[0086] Step S34: Half threshold shrinkage

[0087] The threshold operator F(S (t) , μβ) is

[0088] F(S (t) , μβ) = f(g, μβ) (18)

[0089] The threshold function f(g, μβ) is expressed as

[0090]

[0091] Take the value of f(g, μβ) as the updated value of X (t+1) and take the value of this X (t+1) as the wavelet coefficient for the next iteration calculation; that is, f(g, μβ) = X (t+1) .

[0092] Step S35: Calculate the iteration error at step t + 1

[0093] Resi = ||X (t+1) - X (t) || F (21)

[0094] If the conditions t ≤ T max and Resi > ε are satisfied, continue the iteration, that is, t = t + 1, and repeat the above steps, where ε is a preset threshold. If the conditions are not satisfied, the iteration process ends and X (t) ,

[0095] Finally, through the IDWT and taking the angle value after the inverse transformation as the restored SAR interferogram:

[0096]

[0097] As the restored SAR interferogram, angle(.) represents taking the angle value.

[0098] Next, an experiment is carried out to verify a method for generating SAR interferograms based on L 1 / 2 -norm regularization provided in the present invention. Figure 4 For the simulated interferogram, Gaussian white noise with a standard deviation of 1.0 is added to it to obtain the interferogram shown in (a) in Figure 4 , which represents the interferogram obtained based on MF image processing. Figure 4 The (b) in Figure 4 is the SAR interferogram based on L1-norm regularization; 1 / 2 The (c) in 1 / 2 is the SAR interferogram based on L 1 / 2 -norm regularization. By comparing the quality of the SAR interferograms based on L1 and L 1 / 2 -norm regularization, it is not difficult to see that the methods based on L1 and L

[0099] -norm regularization can both effectively suppress phase noise. The L 1 / 2 -norm regularization method proposed in the present invention can obtain an interferogram with better quality, having fewer residual points and a smaller mean square error value than the L1-norm regularization method. In addition, it should be noted that, in the above specific embodiments, the various specific technical features described can be combined in any suitable manner without contradiction. To avoid unnecessary repetition, the present invention does not further explain various possible combination methods.

Claims

1. A method for generating SAR interferograms based on L 1 / 2 -norm regularization, characterized in that: Specifically, it includes the following steps: Step 1: Use SAR to obtain the SAR main image and the SAR sub-image of the same target scene, and obtain the phase of the SAR sub-image by taking the phase value of the complex matrix of the sub-image; Step 2: Based on the interferometric phase of the SAR master and slave images, and by solving the 1 / 2 ℓ0 -norm regularization problem, establish a generation model for the SAR interferogram; Step 3: Solve the L 1 / 2 -norm regularization problem to obtain the wavelet coefficient matrix related to the restored SAR interferogram Step 4: Perform an inverse discrete wavelet transform on the obtained in Step 3, and take the angle value after the transform as the restored interference image; The phase y of the SAR sub-image in Step 1 is: θ m = exp{j(φ m - φ flat )} u = exp(-jφ topo ) n = θ m u{exp(-jφ noise ) - 1} Among them, φ m is the phase of the SAR main image y m , φ flat is the flat-earth phase, and φ topo is the terrain phase; φ noise = φ s2 - φ s1 , where φ s1 and φ s2 represent the scattering phases of y m and y s respectively, y s represents the SAR sub-image, exp{.} is the exponential power of e, and j represents the imaginary part; The generation model of the SAR interferogram in Step 2 is: Among them, represents the reconstruction of the vector - form y into the matrix - form Y, where Y is a two - dimensional SAR sub - image matrix with an amplitude of 1; represents the vector - form of θ m reconstructed into the matrix - form Θ m , Θ m is a matrix containing the phase of the SAR main image and the flat - ground phase; R(.) is the inverse discrete wavelet transform; represents the Hadamard product; X = W(U), represents the reconstruction of the vector - form u into the matrix - form U, where U represents a matrix containing the interferometric phase of the reconstructed SAR main and sub - images; represents the reconstruction of the vector - form n into the matrix - form N, where N represents the noise matrix, and β represents the regularization parameter, represents the Fibonacci norm, argmin{.} represents taking the minimum value, and W(.) is the discrete wavelet transform.

2. A method for generating an SAR interferogram based on L 1 / 2 norm regularization, characterized in that Step 3 is specifically: Step 3.1: Set the initial value X of the wavelet coefficients (0) = 0; the wavelet coefficients are the wavelet coefficients after the DWT transformation of the interferometric phase of the SAR main and secondary images; Step 3.2: Calculate the residual estimate value during the t-th iteration calculation: X (t) represents the wavelet coefficient at the t-th iterative calculation; Step 3.3: Update the wavelet coefficients during the t-th iteration calculation: S (t) = X (t) + ΔX (t) S (t) Denotes the wavelet coefficient at the t-th iteration calculation after update; Step 3.4: Update the regularization parameter β during the t-th iteration calculation (t) : β (t) = |S (t) | K+1 / μ where μ is a preset parameter, and |S (t) | represents the magnitude component of S (t) , and |S (t) | K+1 represents the (K + 1)-th component after arranging the magnitude components of S (t) in descending order; Step 3.5: Update the wavelet coefficient X for the next iterative calculation according to the value of S (t) (t+1) :​ F(S (t) ,μβ) = f(g,μβ) g = S (t) where μ is the iteration parameter, and the value of f(g,μβ) is taken as the value of X (t+1) ; Step 3.6: Calculate the error of the (t + 1)-th iteration calculation: Resi = ||X (t+1) -X (t) || F Step 3.7: If Resi is greater than the preset threshold ε and the number of iterations t is less than or equal to the preset maximum number of iterations, then the number of iterations is t + 1, and go to Step 3.2; otherwise, stop the calculation.

Images