A Method for InSAR Sparse Imaging with Embedded Backprojection Algorithm

Abstract

The present invention discloses an InSAR sparse imaging method by embedding a back-projection algorithm. First, a conjugate multiplication inverse matrix, an azimuth phase compensation and coherent accumulation inverse matrix, a range interpolation inverse matrix, and a range matching filtering inverse matrix are respectively constructed according to the processing flow of the back-projection imaging algorithm. Then, based on the above inverse matrices, an observation matrix embedded with the back-projection algorithm is constructed. Secondly, according to the constructed observation matrix, a frequency-domain sparse regularization imaging equation is established. Finally, the standard least absolute shrinkage and selection operator is used to solve the equation to obtain the final InSAR interference phase diagram. This method utilizes the frequency-domain sparsity of the interference phase of the imaging scene and the advantages of the back-projection imaging algorithm without registration and post-processing steps for removing flat-earth phase, realizing the direct acquisition of a high-quality interference phase diagram from the echo and simplifying the processing flow.

Description

Technical Field

[0001] The present invention belongs to the technical field of synthetic aperture radar, and particularly relates to the field of interferometric synthetic aperture radar (InSAR) imaging. Background Art

[0002] Synthetic Aperture Radar (SAR) is a radar system with all-weather and all-day imaging capabilities, capable of imaging at any time, whether it is day or night, sunny or rainy or snowy. Conventional SAR obtains high-resolution capabilities in the azimuth direction by the azimuth movement of the platform, which equivalently synthesizes a large aperture antenna in the azimuth direction, and further combines with a large bandwidth transmitted signal to obtain high-resolution capabilities in the range direction. However, conventional SAR can only obtain a two-dimensional image of the imaging scene and it is difficult to obtain the height information of the scene.

[0003] InSAR is a SAR system that can obtain a three-dimensional image of the imaging scene. Compared with the conventional SAR system, it additionally adds a secondary antenna along the height direction on the basis of the main antenna to simultaneously receive the scene echo signal. Therefore, when the platform moves in the azimuth direction, the secondary antenna can equivalently synthesize an additional large aperture antenna. The echoes received by the main and secondary antennas are imaged respectively to obtain two two-dimensional imaging results from different perspectives of the scene, and the interference phase map is obtained by interfering the main and secondary antenna images. Since the interference phase map linearly corresponds to the scene height information, the acquisition of the scene height information is realized, and a three-dimensional image of the imaging scene is obtained. InSAR has wide application values in many fields such as environmental terrain mapping, geological disaster emergency response, building settlement monitoring, etc.

[0004] Currently, the InSAR imaging method includes multiple processing steps, such as separately imaging the echoes of the main and secondary antennas, registering the main and secondary antenna images, conjugate multiplication, removing the flat earth phase, interfering phase filtering, etc., which are relatively cumbersome and complicated, and the quality of the finally obtained interference phase map is affected by multiple steps; in addition, the echo imaging step usually adopts an imaging algorithm based on matched filtering, such as the range-Doppler method, and this kind of method belongs to a non-parametric method and does not consider the characteristics of the imaging scene itself, which results in interference from adverse factors such as sidelobes, clutter, and noise in the imaging result, reducing the quality of the subsequent obtained interference phase map; moreover, since the slant range-azimuth plane is used as the imaging plane, there is a flat earth phase in the interference phase map, resulting in a linear deviation in the range direction of the obtained scene height, and it must be compensated by the step of removing the flat earth phase.

[0005] Therefore, in order to simplify the processing steps of the imaging method and at the same time obtain a high-quality InSAR interferometric phase map, taking advantage of the sparsity of the interferometric phase in the frequency domain brought about by the gentle terrain variation of the imaging scene and the advantages of the back-projection imaging algorithm, such as no need for registration and post-processing steps of removing flat-earth phase, the present invention proposes an InSAR sparse imaging method incorporating the back-projection algorithm. Summary of the Invention

[0006] The present invention belongs to the field of Interferometric Synthetic Aperture Radar (InSAR) imaging, and discloses an InSAR sparse imaging method incorporating the back-projection algorithm, which is used to simplify the steps of the existing InSAR imaging method and at the same time obtain a high-quality InSAR interferometric phase map. This method first constructs a conjugate multiplication inverse matrix, an azimuth phase compensation and coherent accumulation inverse matrix, a range interpolation inverse matrix, and a range matching filtering inverse matrix respectively according to the processing flow of the back-projection imaging algorithm; then, based on the above inverse matrices, constructs an observation matrix embedded with the back-projection algorithm; secondly, according to the constructed observation matrix, establishes a frequency-domain sparse regularization imaging equation; finally, uses the standard least absolute shrinkage and selection operator to solve this equation to obtain the final InSAR interferometric phase map. This method takes advantage of the sparsity of the interferometric phase in the frequency domain of the imaging scene and the advantages of the back-projection imaging algorithm, such as no need for registration and post-processing steps of removing flat-earth phase, realizes directly obtaining a high-quality interferometric phase map from the echo, and greatly simplifies the processing flow.

[0007] To facilitate the description of the content of the present invention, the following term definitions are made first:

[0008] Definition 1. Synthetic Aperture Radar

[0009] Synthetic Aperture Radar (SAR) is a high-resolution microwave imaging radar, which has the advantages of all-weather and all-day operation and has been widely used in various fields, such as topographic mapping, guidance, environmental remote sensing, and resource exploration. An important prerequisite for SAR applications and the main goal of signal processing is to obtain high-resolution and high-precision microwave images through imaging algorithms. For details, see "Pi Yiming, Yang Jianyu, Fu Yusheng, Yang Xiaobo. Principles of Synthetic Aperture Radar Imaging [M]. University of Electronic Science and Technology Press. 2007".

[0010] Definition 2. Traditional Back-projection Algorithm

[0011] The traditional back-projection algorithm uses the position information of the radar platform to calculate the range history between the platform and the scene pixels. Then, by traversing the range history, the corresponding echoes in the data after pulse compression interpolation of the echoes are found, and their phase compensation for the corresponding range is performed. Then, they are coherently accumulated, and the accumulated result is projected into the image space to complete the imaging process. The main steps of the back-projection algorithm include: range matching filtering, range zero-padding interpolation, range echo indexing, azimuth phase compensation, and azimuth coherent accumulation. For details of the traditional back-projection algorithm, see "Shi Jun. Research on the Principles and Imaging Technologies of Bistatic SAR and Linear Array SAR [D]. Doctoral Thesis of University of Electronic Science and Technology of China. 2009".

[0012] Definition 3. Range-matching filtering reference signal

[0013] Range-matching filtering is the range imaging step in the SAR imaging algorithm. By multiplying the spectrum of the reference signal with the spectrum of the echo signal, the phase compensation of the quadratic term in the echo signal is achieved, so that the echo signal can perform pulse compression and obtain the range imaging result. The reference signal is usually the replicated signal of the transmitted signal after time reversal. For details, see "Lan G. Cumming Frank H. Wong. Synthetic Aperture Radar Imaging: Algorithms and Implementation [M]. Publishing House of Electronics Industry, 2012".

[0014] Definition 4. Traditional vector-matrix diagonal operator method

[0015] The vector-matrix diagonal operator diag(a) method refers to the operation method of generating matrix A from the input vector a, and the diagonal elements of A are composed of the vector a. Specifically, the vector-matrix diagonal operator forms the diagonal elements of matrix A from top to bottom in the order of the elements of a column vector a of dimension n×1 from top to bottom, and the other elements of the matrix are 0. diag(a) is expressed as:

[0016]

[0017] where a i , i = 1, 2, 3,... N represents the i-th element of a.

[0018] Definition 5 Traditional matrix-vectorization operator method

[0019] The traditional matrix-vectorization operator vec(A) method refers to the operation method of arranging the input matrix A in columns to form a column vector. Specifically, the matrix-vectorization operator arranges each column of a matrix A of dimension m×n from left to right end to end to form a column vector vec(A):

[0020] vec(A) = [a 1，1 , …, a m，1 , a 1，2 , …, am，2 , …a 1，n , …, a m，n T

[0021] where a i，j represents A(i, j), the element in the i-th row and j-th column of matrix A; the superscript T represents the matrix transpose operation. For the traditional matrix vectorization operator method, see "Zhang Xianda. Matrix Analysis and Applications [M]. Tsinghua University Press Co., Ltd., 2004".

[0022] Definition 6. N-point discrete Fourier transform matrix

[0023] The discrete Fourier transform matrix is an expression that represents the discrete Fourier transform in terms of matrix multiplication. The N-point discrete Fourier transform matrix can implement the N-point discrete Fourier transform. Specifically, the N-point discrete Fourier transform can be represented by an N×N matrix multiplication, i.e., x F = Wx, where x is the original input signal, and x F is the output signal obtained after the discrete Fourier transform. The N-point discrete Fourier transform matrix W is composed as follows:

[0024]

[0025] where i is the imaginary unit. For the N-point discrete Fourier transform matrix, see "He Zishu, Xia Wei. Modern Digital Signal Processing and Its Applications [M]. Beijing: Tsinghua University Press, 2009".

[0026] Definition 7. Standard least absolute value convergence and selection operator method

[0027] The standard least absolute value convergence and selection operator l(y; A, γ) is used to solve the inverse problem based on the vector L1-norm regularization, which can be expressed as the equation where γ is the regularization term weight parameter, is the square of the vector 2-norm. The operator is the solution of this equation, i.e., l(y; A, γ) = x. For the standard least absolute value convergence and selection operator method, see "Zhang Xianda. Matrix Analysis and Applications [M]. Tsinghua University Press Co., Ltd., 2004".

[0028] Definition 8. Traditional vector matrixization operator method

[0029] ​The traditional vector matrix operator unfold(a, m, n) method refers to rearranging the elements of the input vector a to form a matrix with dimensions m×n. When rearranging, it is arranged column by column. When each column reaches the required number of rows, the current column arrangement stops and the next column arrangement begins. For details of the traditional vector matrix operator method, see "Zhang Xianda. Matrix Analysis and Applications [M]. Tsinghua University Press Co., Ltd., 2004".

[0030] An InSAR sparse imaging method using an embedded backprojection algorithm provided by the present invention is characterized by including the following steps:

[0031] Step 1. Initialize relevant parameters

[0032] Initialize the following parameters:

[0033] The preliminary imaging result of the secondary antenna obtained by processing with the backprojection algorithm as defined in Definition 2 is denoted as X s ; The range sampling frequency is denoted as F s ; The speed of light in air is denoted as c; The number of range pixels in the image is denoted as I r ; The number of azimuth pixels in the image is denoted as I a ; The number of azimuth points of the echo is denoted as N a ; The set of scene space coordinates corresponding to the image pixels is denoted as is the scene space coordinate corresponding to the image pixel (i r , i a ); The set of main antenna space coordinates is denoted as is the space coordinate of the main antenna at time n a ; The range reference time delay is denoted as T ref ; The number of range points of the echo is denoted as N r ; The range interpolation multiple is denoted as N i ; The wavelength of the system center frequency is denoted as λ; The number of range matched filtering reference signal points is denoted as L r ; The range matched filtering reference signal as defined in Definition 3 is denoted as h; The echo received by the main antenna is denoted as y m ; The regularization weight coefficient is denoted as ω

[0034] Step 2. Construct the observation matrix embedded with the backprojection algorithm

[0035] According to the preliminary imaging result X of the secondary antenna initialized in Step 1 s , the range interpolation multiple N i , the range sampling frequency F s , the speed of light c in air, the number of range pixels I in the image r , the number of azimuth pixels I in the imagea 、The number of points N in the echo azimuth direction a 、The set P of the corresponding scene space coordinates of the image pixel points I 、The set P of the main antenna space coordinates S 、The range-direction reference time delay T ref 、The number of points N in the echo range direction r 、The system center frequency wavelength λ, the number of points L of the range-direction matched filtering reference signal r and the range-direction matched filtering reference signal h to construct the observation matrix embedded in the back-projection algorithm.

[0036] Step 2.1 According to the preliminary imaging result X of the secondary antenna obtained by initialization in Step 1 s , construct the conjugate multiplication inverse matrix A1 as follows:

[0037] A1 = diag(vec(X s ))

[0038] where diag(·) is the vector matrix diagonal operator defined in Definition 4, vec(·) is the matrix vectorization operator defined in Definition 5, and X s is the preliminary imaging result of the secondary antenna obtained by initialization in Step 1.

[0039] Step 2.2 Use the following formula to calculate the echo range-direction index index(i r , i a , n a ) of each pixel point in the image at different azimuth times:

[0040]

[0041] where ||·||2 represents the vector 2-norm, represents the floor operation, N i 、N r 、F s 、c、I r 、I a 、N a 、 and T ref are the range-direction interpolation multiple, the number of echo range-direction points, the range-direction sampling frequency, the speed of light in air, the number of image range-direction pixels, the number of image azimuth-direction pixels, the number of echo azimuth-direction points, the corresponding scene space coordinates of the image pixel point (i r , i a ), the space coordinates of the main antenna at time n a and the range-direction reference time delay obtained by initialization in Step 1. i r = 1, 2,... I r , i a = 1, 2,... Ia , n a = 1, 2, … N a .

[0042] Step 2.3 uses the following formula to calculate the azimuth phase compensation vector h for each pixel point in the image a (i r , i a ):

[0043]

[0044] where h as (i r , i a , n a ) is calculated as follows:

[0045]

[0046] where exp(·) represents the exponential operation, i represents the imaginary unit, ||·||2 represents the vector 2-norm, the symbol T represents the vector transpose operation, I r , I a , N a , N r , N i , λ,[[]] and are respectively the number of range pixels in the image obtained by initialization in Step 1, the number of azimuth pixels in the image, the number of azimuth points of the echo, the number of range points of the echo, the range interpolation factor, the wavelength of the system center frequency, the scene space coordinates corresponding to the image pixel point (i r , i a ), and the space coordinates of the main antenna at time n a . index(i r , i A , n a ) is the range echo index of each pixel point in the image calculated in Step 2.2 at different azimuth times. i r = 1, 2, … I r , i a = 1, 2, … I a , n a = 1, 2, … N a .

[0047] Step 2.4 constructs the azimuth phase compensation and coherent accumulation inverse matrix A2 of the image as follows:

[0048]

[0049] where hs a (i a ) is calculated as follows:

[0050]

[0051] where \(i\) a = 1, 2, … \(I\) a , the matrix denotes the conjugate transpose of the matrix , \(I\) r , \(I\) a and \(N\) a are respectively the number of pixels in the range direction of the image, the number of pixels in the azimuth direction of the image, and the number of points in the azimuth direction of the echo obtained by initialization in step 1, \(h\) a (\(i\) r , \(i\) a ) is the azimuth phase compensation vector of each pixel point in the image calculated in step 2.3.

[0052] Step 2.5 constructs the range interpolation matrix \(P\) d as follows:

[0053]

[0054] where denotes the first r × \(N\) r rows of the identity matrix of dimension \(N\) ; denotes a matrix with all elements equal to 0 of dimension \([(N\) i - 1)\(N\) r × \(N\) r ; denotes the last r × \(N\) r rows of the identity matrix of dimension \(N\) ; denotes the floor operation, \(N\) r and \(N\) i are the number of range points of the echo and the range interpolation multiple obtained by initialization in step 1.

[0055] Step 2.6 constructs the inverse range interpolation matrix \(A3\) as follows:

[0056]

[0057] where the matrix is the \(N\) i × \(N\) r point discrete Fourier transform matrix defined by Definition 6; the matrix denotes the conjugate transpose of the matrix \(P\) d ; denotes the identity matrix of dimension \(N\) a × \(N\) a , denotes the Kronecker product of matrices. \(N\)a , N r and N i are the number of echo azimuth points, the number of echo range points, and the range interpolation multiple obtained by initializing in Step 1, respectively. P d is the interpolation matrix constructed in Step 2.5.

[0058] Step 2.7 constructs the range matching filter matrix H r as follows:

[0059] H r = repmat(hf pd , N a )

[0060] where hf pd is calculated as follows:

[0061]

[0062] where represents a vector of all 0 elements with dimension , represents a vector of all 0 elements with dimension ; is the N r -point discrete Fourier transform matrix defined by Definition 6; repmat(hf pd , N a ) means replicating the column vector hf pd N a times and arranging them to form a matrix with N a columns; represents the floor operation; N a , N r , L r and h are the number of echo azimuth points, the number of echo range points, the number of range matching filter reference signal points, and the range matching filter reference signal obtained by initializing in Step 1, respectively.

[0063] Step 2.8 constructs the range matching filter inverse matrix A4 as follows:

[0064]

[0065] where the matrix represents the conjugate of the matrix H r ; the matrix represents the conjugate transpose of the matrix , is the N r -point discrete Fourier transform matrix defined by Definition 6; diag(·) is the vector matrix diagonal operator defined by Definition 4; represents a vector with dimension Na ×N a identity matrix of; denotes the matrix Kronecker product; N r 、N a and H r are respectively the number of range bins of the echo obtained by initialization in step 1, the number of azimuth bins of the echo, and the range-matching filtering matrix constructed in step 2.7.

[0066] In step 2.9, based on the conjugate multiplication inverse matrix A1, azimuth-phase compensation and coherent accumulation inverse matrix A2, range interpolation inverse matrix A3, and range-matching filtering inverse matrix A4 constructed respectively in steps 2.1, 2.4, 2.6, and 2.8, the observation matrix A is constructed as follows:

[0067] A = A4A3A2A1

[0068] Step 3. Construct the frequency-domain sparse regularization imaging equation

[0069] The frequency-domain sparse regularization imaging equation is constructed as follows:

[0070]

[0071] where denotes the θ that makes take the minimum value f ; denotes the matrix Kronecker product; the matrix denotes the matrix conjugate transpose of, the matrix denotes the matrix conjugate transpose of, is the I defined in Definition 6 a point discrete Fourier transform matrix, is the I defined in Definition 6 r point discrete Fourier transform matrix; denotes the square of the vector 2-norm, |·| denotes the norm of the vector; I r 、I a 、y m and ω are respectively the number of range pixels of the image, the number of azimuth pixels of the image, the echo received by the main antenna, and the regularization weight coefficient obtained by initialization in step 1, A is the observation matrix embedded in the back-projection algorithm constructed in step 2; θ f is the frequency-domain interferometric phase map to be solved;

[0072] Step 4. Solve the frequency-domain sparse regularization imaging equation

[0073] For the frequency-domain sparse regularization imaging equation constructed in step 3, the standard least absolute value convergence and selection operator method defined in Definition 7 is used to solve it, and the frequency-domain interferometric phase diagram θ to be solved is obtained. f , as follows:

[0074] θ f = l(y m ; A′, ω)

[0075] where A is the observation matrix constructed in step 2, denotes the matrix Kronecker product, the matrix denotes the conjugate transpose of the matrix , the matrix denotes the conjugate transpose of the matrix . is the I a point discrete Fourier transform matrix defined in Definition 6, is the I r point discrete Fourier transform matrix defined in Definition 6, I r , I a , y m and ω are the number of pixels in the range direction of the image, the number of pixels in the azimuth direction of the image, the echo received by the main antenna, and the regularization weight coefficient initialized in step 1 respectively. l(y m ; A′, ω) is the standard least absolute value convergence and selection operator defined in Definition 7.

[0076] Step 5. Conversion of the frequency-domain interferometric phase diagram

[0077] According to the number of pixels I r in the range direction of the image and the number of pixels I a in the azimuth direction of the image initialized in step 1, and the frequency-domain interferometric phase diagram θ f solved in step 4, the time-domain interferometric phase diagram Θ is calculated using the formula . Among them, the matrix denotes the conjugate transpose of the matrix , the matrix denotes the conjugate of the matrix , is the I a point discrete Fourier transform matrix defined in Definition 6, is the I r point discrete Fourier transform matrix defined in Definition 6; unfold(·) is the vector matrix operator defined in Definition 8, and Θ is the final interferometric phase diagram to be obtained.

[0078] The innovation and advantages of the present invention are as follows: An InSAR sparse imaging method incorporating the back-projection algorithm is proposed. The advantage lies in that the terrain slow-variation characteristics commonly possessed by the imaging scene are considered in the method design, and regularization imaging is carried out through the sparsity of the interferometric phase in the frequency domain, enabling the acquisition of a high-quality interferometric phase map. The innovation is that the observation matrix required for regularization imaging is constructed based on the processing flow of the back-projection algorithm, which can avoid numerous post-processing steps required by existing imaging methods, realizing the direct obtaining of a high-quality interferometric phase map from the echo and greatly simplifying the processing flow. Description of the Drawings

[0079] Figure 1 It is the flowchart of the present invention.

[0080] Figure 2 It is the verification result of the simulation experiment of the present invention. Detailed Implementation Manner

[0081] The present invention is mainly verified by the method of simulation experiment, and all steps and conclusions are verified correctly on the mathematical calculation software Matlab2019b. The specific implementation steps are as follows:

[0082] Step 1. Initialize relevant parameters

[0083] Initialize the following parameters:

[0084] Define the preliminary imaging result X of the secondary antenna obtained by processing with the back-projection algorithm as defined in 2 s ; Range sampling frequency F s = 6×10 8 Hz; The speed of light c in air = 3×10 8 m / s; The number of range pixels of the image I r = 512; The number of azimuth pixels of the image I a = 512; The number of azimuth points of the echo N a = 400; The set of scene space coordinates corresponding to the image pixels, denoted as The set of main antenna space coordinates Range reference time delay T ref = 0s; The number of range points of the echo N r = 1024; Range interpolation multiple N i = 8; The wavelength λ of the system center frequency = 0.0125m; The number of range matched filtering reference signal points L r = 300; The range matched filtering reference signal h defined in 3; The echo y received by the main antenna m ; Regularization weight coefficient ω = 0.5

[0085] Step 2. Construct the observation matrix embedded with the back projection algorithm

[0086] According to the preliminary imaging result X of the secondary antenna obtained by initialization in Step 1 s , the range interpolation multiple N i = 8, the range sampling frequency F s = 6×10 8 Hz, the speed of light c in air = 3×10 8 m / s, the number of range pixels I of the image r = 512, the number of azimuth pixels I of the image a = 512, the number of azimuth points N of the echo a = 400, the set P of the spatial coordinates of the scene corresponding to the image pixels I , the set P of the spatial coordinates of the main antenna S , the reference time delay T in the range direction ref = 0s, the number of range points N of the echo r = 1024, the wavelength λ of the system center frequency = 0.0125m, the number of reference signals L for range matched filtering r = 300 and the range matched filtering reference signal h, construct the observation matrix embedded with the back projection algorithm.

[0087] Step 2.1 According to the preliminary imaging result X of the secondary antenna obtained by initialization in Step 1 s , construct the conjugate multiplication inverse matrix A1 as follows:

[0088] A1 = diag(vec(X s ))

[0089] where diag(·) is the vector matrix diagonal operator defined in Definition 4, and vec(·) is the matrix vectorization operator defined in Definition 5, and X s is the preliminary imaging result of the secondary antenna obtained by initialization in Step 1.

[0090] Step 2.2 Calculate the range index index(i r , i A , n A ) of the echo at each pixel point in the image at different azimuth times as follows:

[0091]

[0092] where ||·||2 represents the vector 2-norm, represents the floor operation, and are the spatial coordinates of the scene corresponding to the image pixel (i r , i a ) obtained by initialization in Step 1, and the main antenna at na Spatial coordinates of the moment, i r = 1, 2, … 512, i a = 1, 2, … 512, n a = 1, 2, … 400.

[0093] Step 2.3 Calculate the azimuth phase compensation vector h for each pixel point in the image a (i r , i a ) is as follows:

[0094]

[0095] where h as (i r , i A , n a ) is calculated as follows:

[0096]

[0097] where exp(·) represents the exponential operation, i represents the imaginary unit, ||·||2 represents the vector 2-norm, the symbol T represents the vector transpose operation, and are respectively the scene spatial coordinates corresponding to the image pixel point (i r , i a ) initialized in Step 1 and the spatial coordinates of the main antenna at the moment n a , i r = 1, 2, … 512, i a = 1, 2, … 512, n a = 1, 2, … 400.

[0098] Step 2.4, construct the azimuth phase compensation and coherent accumulation inverse matrix A2 of the image as follows:

[0099]

[0100] where hs a (i a ) is calculated as follows:

[0101]

[0102] where the matrix represents the conjugate transpose of the matrix , h a (i r , i a ) is the azimuth phase compensation vector for each pixel point in the image calculated in Step 2.3, i a = 1, 2, … 512, i r= 1, 2, … 512.

[0103] Step 2.5 Construct the range interpolation matrix P d as follows:

[0104]

[0105] where Ib 512×1024 represents the first 512×1024 rows of the 1024×1024 identity matrix; 0 [(8-1)×1024]×1024 represents a matrix of all 0 elements with dimensions [(8 - 1)×1024]×1024; Ie (1024-512)×1024 represents the last (1024 - 512)×1024 rows of the 1024×1024 identity matrix.

[0106] Step 2.6 Construct the inverse range interpolation matrix A3 as follows:

[0107]

[0108] where the matrix F 8×1024 is the 8×1024-point discrete Fourier transform matrix defined by Definition 6; the matrix represents the conjugate transpose of the matrix P d ; I 400×400 represents the 400×400 identity matrix, represents the Kronecker product of matrices.

[0109] Step 2.7 Construct the range matched filtering matrix H r as follows:

[0110]

[0111] hf pd = F 1024 h pd

[0112] H r = repmat(hf pd , 400)

[0113] where represents a vector of all 0 elements with dimensions ; represents a vector of all 0 elements with dimensions ; F 1024 is the 1024-point discrete Fourier transform matrix defined by Definition 6; repmat(hf pd , 400) means replicating the column vector hf pd 400 times and arranging them to form a matrix with 400 columns; Denotes the floor operation. h is the range-direction matched filtering reference signal obtained by initialization in Step 1.

[0114] Step 2.8 constructs the range-direction matched filtering inverse matrix A4 as follows:

[0115]

[0116] where the matrix denotes the conjugate of matrix H r ; the matrix denotes the conjugate transpose of matrix F 1024 , F 1024 is the 1024-point discrete Fourier transform matrix defined by Definition 6; diag(·) is the vector matrix diagonal operator defined by Definition 4; I 400×400 denotes the identity matrix of dimension 400×400; denotes the Kronecker product of matrices. H r is the range-direction matched filtering matrix obtained by initialization in Step 1.

[0117] Step 2.9 constructs the observation matrix A as follows according to the conjugate multiplication inverse matrix A1, azimuth-phase compensation and coherent accumulation inverse matrix A2, range-direction interpolation inverse matrix A3, and range-direction matched filtering inverse matrix A4 constructed in Steps 2.1, 2.4, 2.6, and 2.8 respectively:

[0118] A = A4A3A2A1

[0119] Step 3. Construct the frequency-domain sparse regularization imaging equation

[0120] Construct the frequency-domain sparse regularization imaging equation as follows:

[0121]

[0122] where denotes the θ that makes take the minimum value f ; denotes the Kronecker product of matrices; the matrix denotes the conjugate transpose F of matrix F 512 is the 512-point discrete Fourier transform matrix defined by Definition 6, F 512 is the 512-point discrete Fourier transform matrix defined by Definition 6; 512 is the 512-point discrete Fourier transform matrix defined by Definition 6; denotes the square of the vector 2-norm, |·| denotes the modulus of the vector; θ f is the frequency-domain interferometric phase map to be solved.

[0123] Step 4. Solve the frequency-domain sparse regularization imaging equation

[0124] According to the frequency-domain sparse regularization imaging equation constructed in step 3, solve the equation using the standard least absolute shrinkage and selection operator defined in Definition 7 to obtain the frequency-domain interference phase diagram θ to be solved. f , as follows:

[0125] θ f = l(y m ; A′, 0.5)

[0126] where A is the observation matrix constructed in step 2, represents the matrix Kronecker product, and the matrix represents the conjugate transpose of matrix F 512 of, and the matrix represents the conjugate transpose of matrix F 512 of, F 512 is the 512-point discrete Fourier transform matrix defined in Definition 6, and F 512 is the 512-point discrete Fourier transform matrix defined in Definition 6, and l(y m ; A′, 0.5) is the standard least absolute shrinkage and selection operator defined in Definition 7.

[0127] Step 5. Conversion of the frequency-domain interference phase diagram

[0128] According to the number of pixels I r = 512 in the range direction of the image and the number of pixels I a = 512 in the azimuth direction of the image obtained by initialization in step 1, the frequency-domain interference phase diagram θ f solved in step 4 is converted into the time-domain interference phase diagram Θ through the formula , where the matrix represents the conjugate transpose of matrix F 512 of, and the matrix represents the conjugate of matrix F 512 of, F 512 is the 512-point discrete Fourier transform matrix defined in Definition 6, and F 512 is the 512-point discrete Fourier transform matrix defined in Definition 6; unfold(·) is the vector matrix operator defined in Definition 8, and Θ is the final interference phase diagram to be obtained. The computer simulation results are as Figure 2 shown. The computer simulation results show that the present invention realizes directly obtaining a high-quality interference phase diagram from the echo, and the processing flow is greatly simplified.

Claims

1. An InSAR sparse imaging method using an embedded back-projection algorithm, characterized in that Including the following steps: Step 1. Initialize relevant parameters Initialize the following parameters: The preliminary imaging result of the secondary antenna obtained by the back-projection algorithm is denoted as X s ; The range sampling frequency is denoted as F s ; The speed of light in air is denoted as c; The number of range pixels in the image is denoted as I r ; The number of azimuth pixels in the image is denoted as I a ; The number of azimuth points of the echo is denoted as N a ; The set of scene space coordinates corresponding to the image pixel points is denoted as For the scene space coordinates corresponding to the image pixel point (i r , i a ); The set of main antenna space coordinates is denoted as For the space coordinates of the main antenna at time n a ; The range - direction reference time delay is denoted as T ref ; The number of range - direction echo points is denoted as N r ; The range - direction interpolation multiple is denoted as N i ; The wavelength of the system center frequency is denoted as λ; The number of range - direction matched - filtering reference signal points is denoted as L r ; The range - direction matched - filtering reference signal is denoted as h; The echo received by the main antenna is denoted as y m ; The regularization weight coefficient is denoted as ω Step 2. Construct the observation matrix embedded with the back-projection algorithm The preliminary imaging result X of the secondary antenna obtained by initializing according to Step 1 s , the interpolation multiple N in the range direction i , the sampling frequency F in the range direction s , the speed c of light propagating in air, the number of pixels I in the range direction of the image r , the number of pixels J in the azimuth direction of the image a , the number of points N in the azimuth direction of the echo a , the set P of scene space coordinates corresponding to the image pixels I , the set P of main antenna space coordinates S , the reference time delay T in the range direction ref , the number of points N in the range direction of the echo r , the wavelength λ of the system center frequency, the number of points L of the reference signal for range direction matched filtering r and the reference signal h for range direction matched filtering, and construct the observation matrix embedded in the back projection algorithm; Step 2.1 Based on the preliminary imaging result X of the secondary antenna obtained by initialization in Step 1 s , construct the conjugate multiplication inverse matrix A1 as follows: A1 = diag(vec(X s )) where diag(·) is the diagonal operator for vectorizing a matrix, vec(·) is the operator for vectorizing a matrix, and X s is the preliminary imaging result of the secondary antenna obtained by the initialization in step 1; Step 2.2 uses the following formula to calculate the echo range index index(i r , i a , n a ) for each pixel point in the image at different azimuth times: where $\|\cdot\|_2$ represents the vector 2-norm, represents the floor operation, $N$ i , $N$ r , $F$ s , $c$, $I$ r , $I$ a , $N$ a , and $T$ ref are respectively the range interpolation multiple obtained by initialization in step 1, the number of range bins of the echo, the range sampling frequency, the speed of light in air, the number of range pixels of the image, the number of azimuth pixels of the image, the number of azimuth bins of the echo, the scene space coordinates corresponding to the image pixel $(i$ r , $i$ a ), the space coordinates of the main antenna at time $n$ a and the range reference time delay; $i$ r = 1, 2, … $I$ r , $i$ a = 1, 2, … $I$ a , $n$ a = 1, 2, … $N$ a ; Step 2.3 uses the following formula to calculate the azimuth phase compensation vector h for each pixel in the image a (i r ,i a ): where h as (i r , i a , n a ) is calculated as follows: where exp(·) represents the exponential operation, i represents the imaginary unit, ‖·‖2 represents the vector 2-norm, the symbol T represents the vector transpose operation, I r 、I a 、N a 、N r 、N i 、λ、 and are respectively the number of pixels in the range direction of the image, the number of pixels in the azimuth direction of the image, the number of points in the azimuth direction of the echo, the number of points in the range direction of the echo, the interpolation multiple in the range direction, the wavelength of the system center frequency, the scene space coordinates corresponding to the image pixel (i r ,i a ), and the space coordinates of the main antenna at time n a . index(i r ,i a ,n a ) is the range direction index of the echo of each pixel in the image calculated in step 2.2 at different azimuth times; i r = 1, 2, … I r ,i a = 1, 2, … I a ,n a = 1, 2, … N a ; Step 2.4 Construct the azimuth phase compensation and coherent accumulation inverse matrix A2 of the image as follows: Among which hs a (i a ) is calculated as follows: where i a = 1, 2, … I a , the matrix represents the conjugate transpose of the matrix , I r , I a and N a are respectively the number of range pixels, the number of azimuth pixels, and the number of echo azimuth points obtained by initialization in step 1, h a (i r , i a ) is the azimuth phase compensation vector of each pixel point in the image calculated in step 2.3; Step 2.5 Construct the range interpolation matrix P d As follows: Among them represents the first r ×N r rows of the identity matrix of dimension N×N; rows; represents a matrix of all 0 elements of dimension [(N i -1)N r ×N; r represents the last r ×N r rows of the identity matrix of dimension N×N; rows; represents the floor operation, where N r and N i are the number of points in the echo range direction and the range direction interpolation multiple obtained by initialization in step 1;​ Step 2.6 Construct the range interpolation inverse matrix A3 as follows: where the matrix is the N i × r point discrete Fourier transform matrix; the matrix represents the conjugate transpose of matrix P d ; represents the identity matrix of dimension N a × a ×N, represents the matrix Kronecker product; N a , N r and N i are the number of points in the azimuth direction of the echo, the number of points in the range direction of the echo, and the range interpolation factor obtained by the initialization in step 1 respectively, and P d is the interpolation matrix constructed in step 2.5; Step 2.7 Construct the range matching filter matrix H r As follows: H r = repmat(hf pd , N a ) where hf pd is calculated as follows: Among them represents a vector composed of all 0 elements with a dimension of ; represents a vector composed of all 0 elements with a dimension of ; is the N r point discrete Fourier transform matrix; repmat(hf pd , N a ) means replicating the column vector hf pd N a times and arranging them to form a matrix with the number of columns being N a ; represents the floor operation; N a , N r , L r and h are the number of echo azimuth points, the number of echo range points, the number of range - direction matched - filtering reference signals, and the range - direction matched - filtering reference signal obtained by initializing in step 1 respectively; Step 2.8 Construct the range matched filtering inverse matrix A4 as follows: where the matrix represents the conjugate of matrix H r ; the matrix represents the conjugate transpose of matrix ; is the N r -point discrete Fourier transform matrix; diag(·) is the vector-to-matrix diagonal operator; represents the identity matrix of dimension N a ×N a ; represents the Kronecker product of matrices; N r , N a and H r are the number of range bins of the echo obtained by initialization in step 1, the number of azimuth bins of the echo, and the range-direction matched filtering matrix constructed in step 2.7, respectively; Step 2.9 According to the conjugate multiplication inverse matrix A1, azimuth phase compensation and coherent accumulation inverse matrix A2, range interpolation inverse matrix A3, and range matched filtering inverse matrix A4 constructed in Steps 2.1, 2.4, 2.6, and 2.8 respectively, construct the observation matrix A as follows: A = A4A3A2A1 Step 3. Construct the frequency-domain sparse regularization imaging equation Construct the frequency-domain sparse regularization imaging equation as follows: wherein denotes the θ that takes the minimum value f ; denotes the matrix Kronecker product; the matrix denotes the matrix 's conjugate transpose, and the matrix denotes the matrix 's conjugate transpose, is I a point discrete Fourier transform matrix, is I r point discrete Fourier transform matrix; denotes the square of the vector 2-norm, |·| denotes the vector norm; I r 、I a 、y m and ω are respectively the number of pixels in the range direction of the image, the number of pixels in the azimuth direction of the image, the echo received by the main antenna, and the regularization weight coefficient obtained by initialization in step 1, A is the observation matrix embedded in the back-projection algorithm constructed in step 2; θ f is the interferometric phase diagram in the frequency domain to be solved; Step 4. Solve the frequency-domain sparse regularization imaging equation For the frequency-domain sparse regularization imaging equation constructed in step 3, the standard least absolute shrinkage and selection operator method is used to solve it, and the frequency-domain interference phase map θ to be solved is obtained. f , as follows: θ f = l(y m ; A′, ω) Among them A is the observation matrix constructed in step 2, represents the Kronecker product of matrices, and the matrix represents the matrix is the conjugate transpose of the matrix, and the matrix represents the matrix is the conjugate transpose of the matrix, is the I a point discrete Fourier transform matrix, is the I r point discrete Fourier transform matrix, I r 、I a 、y m and ω are respectively the number of pixels in the range direction of the image, the number of pixels in the azimuth direction of the image, the echo received by the main antenna, and the regularization weight coefficient obtained by initialization in step 1. l(y m ; A′, ω) is the standard least absolute value convergence and selection operator; Step 5. Frequency-domain interferometric phase diagram conversion According to the number of pixels I in the image distance obtained by initialization in step 1 r and the number of pixels I in the image azimuth direction a , the frequency-domain interference phase diagram θ obtained by solving in step 4 f , using the formula to calculate the time-domain interference phase diagram Θ, where the matrix represents the conjugate transpose of the matrix , and the matrix represents the conjugate of the matrix ; is the I a point discrete Fourier transform matrix, is the I r point discrete Fourier transform matrix; unfold(·) is the vector matrix operator, and Θ is the finally required interference phase diagram.

Images