A Three-Dimensional Sparse Imaging Method for Millimeter Waves Based on Generalized MCP Constraint

Abstract

The present invention discloses a method for millimeter-wave three-dimensional sparse imaging based on generalized MCP constraint. Firstly, range pulse compression processing of the echo signal is completed by using pulse compression and frequency upsampling techniques; then, the standard range migration correction algorithm is used to process the signal after range pulse compression to complete range migration correction; then, for each range cell, an imaging equation approximating the RMA imaging operator based on generalized MCP non-convex constraint is constructed respectively, and optimized to solve to obtain the corresponding azimuth-height imaging result; finally, the overall three-dimensional sparse imaging result of the target scene is obtained by stacking the imaging results of each range cell along the range direction. Compared with the existing millimeter-wave three-dimensional sparse imaging method based on L1 norm constraint, the method of the present invention is more suitable for the situation of unbalanced energy distribution of scene targets, has small estimation deviation, more complete retention of geometric contours, and high image quality.

Description

Technical Field

[0001] The present invention belongs to the technical field of radar imaging, and particularly relates to the technical field of millimeter-wave imaging. Background Art

[0002] In recent years, three-dimensional millimeter-wave imaging technology has received great attention. Compared with microwaves, millimeter waves have a certain penetration ability through clothing; compared with X-ray rays, millimeter waves do not cause ionizing radiation damage and have high biological safety. Therefore, millimeter-wave imaging technology has a wide range of applications in many fields such as concealed weapon inspection, biomedical diagnosis, and non-contact security inspection. Millimeter-wave imaging has a similar imaging principle to traditional microwave imaging. However, compared with microwave imaging, the emitted signal has a higher carrier frequency and bandwidth, so higher-resolution imaging results can be obtained. In addition, the system can form a two-dimensional sampling aperture on the two-dimensional azimuth-height plane through a multi-channel two-dimensional planar array or a single-channel scanning synthesis planar array, so that a three-dimensional imaging result of the scene can be obtained, and the resolution on the azimuth-height plane can reach the millimeter level, which is very beneficial for the detection and recognition of targets in the subsequent imaging results.

[0003] However, limited by the Nyquist sampling theorem, the adjacent spatial sampling intervals on the two-dimensional aperture of the azimuth-height plane need to meet the requirement of not exceeding half a wavelength. This will result in a large number of sampling channels required by the system, significantly increasing the cost of the system. In recent years, with the development of sparse imaging technology, when the imaging scene has sparsity, the system sampling can break through the limitation of the Nyquist sampling theorem, and the required number of two-dimensional aperture samplings can be greatly reduced. The sparse imaging method forms an imaging model with multiple constraints by adding an additional sparse constraint term about the observed scene to the traditional imaging model. During imaging, through the method of optimal iterative solution, an estimation of the observed scene is obtained, which is the imaging result of the observed scene. In addition to reducing the requirement for the number of two-dimensional aperture samplings, compared with the imaging results of traditional three-dimensional imaging methods, the sparse imaging method can also obtain higher-quality imaging results, manifested in higher target resolution, lower sidelobes and clutter, etc.

[0004] However, the current methods mainly use constraint terms based on the L1 norm, which are applicable to the case where the energy distributions of individual point scatterers are relatively balanced. When it is unbalanced, such as for targets composed of different materials, the weights of point scatterers with stronger energy are too large, resulting in unbalanced constraints, and thus estimation deviation. In the imaging result, it is manifested as energy loss of some scatterers, usually corresponding to the geometric contour of the target, which has an adverse impact on the subsequent detection and recognition of the target. Therefore, in order to solve the problem of estimation deviation and thus retain more geometric contours of the target, the present invention proposes a millimeter-wave three-dimensional sparse imaging method based on the generalized MCP constraint. Summary of the Invention

[0005] The present invention proposes a millimeter-wave three-dimensional sparse imaging method based on a generalized MCP constraint. First, it uses pulse compression and frequency upsampling techniques to complete the range compression processing of the echo signal; then it uses a standard range migration correction algorithm to process the signal after range compression to complete range migration correction; then, for each range cell, an imaging equation approximating the RMA imaging operator based on the generalized MCP non-convex constraint is constructed respectively, and optimized and solved to obtain the corresponding azimuth-height imaging result; finally, the overall three-dimensional sparse imaging result of the target scene is obtained by stacking the imaging results of each range cell along the range direction. Compared with the millimeter-wave three-dimensional sparse imaging method based on the L1-norm constraint, the method of the present invention is more applicable to the situation where the energy distribution of the scene target is uneven, the estimation deviation is smaller, the geometric contour is more completely retained, and the image quality is higher.

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

[0007] Definition 1. Azimuth, height, and range

[0008] The direction of the lateral movement of the radar platform is called the azimuth, the direction of the longitudinal movement is called the height, and the direction perpendicular to the former two is called the range.

[0009] Definition 2. L1-norm constraint term

[0010] The L1-norm regularization term is a kind of target prior information constraint term, which is composed of the L1-norm of the target scene image X and can be expressed as ‖X‖1 = ∑|x|. It is a convex function and is specifically calculated by accumulating the absolute values of each image pixel x. This regularization term is a common constraint term for existing three-dimensional SAR sparse imaging methods and is used to enhance the sparsity of the obtained target scene image. For details, see "Wang Y, Zhang X, Zhan X, et al. An RCS Measurement Method Using Sparse ImagingBased 3D SAR Complex Image[J]. IEEE Antennas and Wireless Propagation Letters, 2021".

[0011] Definition 3. Generalized MCP non-convex constraint term

[0012] The generalized minimax-concave penalty (GMC penalty) is a non-convex constraint term of target prior information and can be expressed as where int represents the infimum and ||v||1 represents the 1-norm of the vector v. Denotes the square of the 2-norm of the vector B(x - v). Matrix B is a hyperparameter matrix. For details, see Selesnick I. Sparse regularization via convex analysis[J]. IEEE Transactions on Signal Processing, 2017, 65(17): 4481 - 4494.

[0013] Definition 4. Standard pulse compression method

[0014] Pulse compression is a modern radar signal processing method. Briefly speaking, the radar emits wide pulses and then "compresses" them into narrow pulses at the receiving end, thereby improving two radar performances: operating range and range resolution.

[0015] For details of the standard pulse compression method, see "Pi Yiming, Yang Jianyu, Fu Yusheng, Yang Xiaobo. Principles of Synthetic Aperture Radar Imaging [M]. University of Electronic Science and Technology Press. 2007".

[0016] Definition 5. Traditional upsampling method

[0017] The traditional upsampling method is a method to increase the signal sampling rate in the discrete signal domain. There are two implementation methods: time-domain upsampling and frequency-domain upsampling. For example, frequency-domain upsampling is completed by padding zeros in the frequency domain. For details of the traditional upsampling method, see "Shi Jun. Research on the Principles and Imaging Technologies of Bistatic SAR and Linear Array SAR [D]. PhD Thesis of University of Electronic Science and Technology. 2009".

[0018] Definition 6. Traditional fast Fourier transform pair (FFT / IFFT) method

[0019] The traditional fast Fourier transform pair method is a fast algorithm for calculating the discrete Fourier transform pair, which can be divided into fast Fourier transform (FFT) and inverse fast Fourier transform (IFFT). Using this algorithm can greatly reduce the number of multiplications required by the computer to calculate the discrete Fourier transform. Especially when the number of sampling points to be transformed is larger, the saving of the FFT / IFFT algorithm calculation amount is more significant. For details of the traditional fast Fourier transform pair method, see "Cheng Qiansheng. Digital Signal Processing [M]. Peking University Press. 2003".

[0020] Definition 7. Traditional range migration correction method

[0021] The traditional range migration correction method is used to correct the range migration in the range direction imaging result caused by the different distances between the radar and the target at different sampling points, so that the imaging results of the target at different sampling points are located in the same range cell after range migration correction. The range migration correction can be completed by time-domain interpolation or frequency linear phase compensation. For details of the traditional range migration correction method, see "Lan G.Cumming, Frank H.Wong. Synthetic Aperture Radar Imaging: Algorithms and Implementation [M]. Publishing House of Electronics Industry, 2012".

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

[0023] The discrete Fourier transform matrix is an expression that represents the discrete Fourier transform in 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, that is, 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 of the following:

[0024]

[0025] where i is the imaginary unit. For details of the traditional 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 9. Traditional RMA method

[0027] The traditional Range Migration Algorithm (RMA) method is a three-dimensional SAR imaging algorithm in the wavenumber domain. The algorithm performs imaging processing on the echo after range direction pulse compression, upsampling, and range migration correction for each range cell. The main processing steps include azimuth-height frequency domain transformation, frequency phase compensation, and inverse frequency domain transformation. For details of the traditional RMA method, see "Wang M, Wei S, Liang J, et al. RMIST-net: Joint range migration and sparse reconstruction network for 3-D mmW imaging [J]. IEEE Transactions on Geoscience and Remote Sensing, 2021, 60: 1-17".

[0028] Definition 10. Traditional vector matrix diagonal operator method

[0029] The traditional vector matrix diagonal operator diag(a) refers to an operation method that generates a 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. dian(a) is expressed as:

[0030]

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

[0032] Definition 11. Traditional matrix vectorization operator method

[0033] The traditional matrix vectorization operator vec(A) method refers to an operation method that arranges the input matrix A by 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):

[0034] vec(A) = [a 1,1 , …, a m,1 , a 1,2 , …, a m,2 , … a 1,n , …, a m,n T

[0035] 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 details of the traditional matrix vectorization operator method, see "Zhang Xianda. Matrix Analysis and Applications [M]. Tsinghua University Press Co., Ltd., 2004".

[0036] Definition 12. Standard distance-wise stacking operation method

[0037] ​The standard range-stack operation method refers to a common operation method in the SAR frequency-domain 3D imaging algorithm. During the imaging process, 2D imaging is performed on the azimuth-height direction corresponding to each range cell to obtain 2D imaging results. Stacking arranges the 2D imaging results of each range cell in the order of range cells to form the final 3D imaging result. For details of the standard range-stack operation method, see "Wang M, Wei S, Liang J, et al. RMIST-net: Joint range migration and sparse reconstruction network for 3-D mmW imaging[J]. IEEE Transactions on Geoscience and Remote Sensing, 2021, 60: 1-17".

[0038] A millimeter-wave three-dimensional sparse imaging method based on the generalized MCP constraint provided by the present invention is characterized in that it includes the following steps:

[0039] Step 1. Initialize relevant parameters

[0040] The speed of light in air is denoted as c; the natural exponential function is denoted as exp(·); the imaginary unit is denoted as j; the pi is denoted as π; the wavelength of the transmitted signal is denoted as λ; the azimuth sampling sequence number is denoted as l = 1, 2, …, L, where L represents the total number of azimuth samplings; the azimuth sampling interval is denoted as d l ; the height sampling sequence number is denoted as m = 1, 2, …, M, where M represents the total number of height samplings; the height sampling interval is denoted as d m ; the total number of range samplings is denoted as N; the upsampling factor is denoted as K; the reference range is denoted as R0; the range sampling interval is denoted as d r ; the original target echo is denoted as S L×M×N ; the weight coefficient of the generalized non-convex MCP constraint term is denoted as β;

[0041] Step 2. Perform range pulse compression processing on the original target echo to obtain the range pulse compression result

[0042] Using the original target echo S in Step 1 L×M×N as the input, and adopting the pulse compression method defined in Standard 4 to compress the range signals corresponding to each pair of azimuth-height sampling points in S L×M×N to obtain the result P after range pulse compression L×M×N .

[0043] Step 3. Perform upsampling on the result after pulse compression

[0044] Using the result P after range pulse compression obtained in Step 2L×M×N Taking the upsampling multiple K initialized in Step 1 as the input, perform K-fold upsampling processing. After range compression, the result is P L×M×N For each pair of azimuth-height sampling points in the result P lm , the result after range compression is denoted as p lm , where lm represents the l-th azimuth - m-th height sampling point pair, l = 1, 2, …, L, m = 1, 2, …, M, and L and M respectively represent the total number of azimuth samplings and the total number of height samplings. For each result p L×M×(KN) Perform the following operations to obtain the range upsampling pulse compression result P′

[0045] Step 3.1.

[0046] Process the vector p using the traditional fast Fourier transform (FFT) described in Definition 6 lm to obtain the vector f lm .

[0047] Step 3.2.

[0048] Starting from the lm position of the vector f , insert (K - 1)N zero elements to obtain where represents the first lm elements in f , represents the lm -th to the last element in f , 0 (K-1)·N represents the inserted (K - 1)·N zero elements, represents rounding down, f lm is the vector obtained after fast Fourier transform processing in Step 3.1, N is the total number of range samplings initialized in Step 1, and K is the upsampling multiple.

[0049] Step 3.3.

[0050] Process the vector f′ using the traditional inverse fast Fourier transform (IFFT) described in Definition 6 lm to obtain the vector p′ lm .

[0051] Step 4. Range migration correction

[0052] Perform range migration correction on the range upsampling pulse compression result P′ obtained in Step 3 using the traditional range migration correction method described in Definition 7 L×M×(KN) to obtain the range upsampling pulse compression result Pc after range migration correction L×M×(KN), K, N, L, and M are the upsampling factor, the total number of range samples, the total number of azimuth samples, and the total number of height samples initialized in step 1, respectively.

[0053] Step 5. Construct an imaging solution equation for approximating the RMA imaging operator based on the generalized MCP non-convex constraint

[0054] Using the upsampled pulse compression result Pc after range migration correction obtained in step 4 L×M×(KN) as the input, construct an imaging solution equation for approximating the RMA imaging operator based on the generalized MCP non-convex constraint. The imaging solution equation corresponding to each range cell is constructed in the following step-by-step manner.

[0055] Step 5.1. Calculate the inverse matrix Φ of the RMA operator phase compensation corresponding to the current range cell using the following formula i :

[0056]

[0057] where exp, j, π, λ, R0, K, N, L, M, d r , d l and d m are the natural exponential, imaginary unit, pi, transmit signal wavelength, reference distance, upsampling factor, total number of range samples, total number of azimuth samples, total number of height samples, range sampling interval, azimuth sampling interval, and height sampling interval initialized in step 1, respectively.

[0058] Step 5.2. Construct the imaging solution equation corresponding to the current range cell as follows,

[0059]

[0060] where represents the X that minimizes , Pc(i) is the upsampled pulse compression result after range migration correction corresponding to the i-th range cell in Pc i , i = 1, 2,..., KN, and Pc L×M×(KN) is the upsampled pulse compression result after range migration correction obtained in step 4; F L×M×(KN) and F L and F M are the L-point and M-point discrete Fourier transform matrices described in Definition 8, respectively, represents the transpose of F M , represents the conjugate transpose of F L , represents F MThe conjugate of; ⊙ represents the matrix Hadamard product; represents the square of the matrix Fibonacci norm; β is the weight coefficient of the generalized MCP non-convex constraint term initialized in Step 1; is the weight of the generalized MCP non-convex constraint term described in Definition 3, x i = vec(X i ), Φ i is the inverse matrix of the phase compensation of the RMA operator corresponding to the current distance cell calculated in Step 5.1; represents the matrix Kronecker product, diag(·) is the vector matrix diagonal operator described in Definition 10, and vec(·) is the matrix vectorization operator described in Definition 11. X i is the imaging solution equation corresponding to the i-th distance cell constructed. K, N, L, and M are the upsampling multiple, the total number of range samples, the total number of azimuth samples, and the total number of height samples initialized in Step 1, respectively.

[0061] Step 6. Solve the imaging solution equation for the approximation of the RMA imaging operator based on the generalized MCP non-convex constraint

[0062] Solve the imaging solution equation corresponding to each distance cell obtained in Step 5 by an iterative method to obtain the azimuth-height imaging result corresponding to each distance cell.

[0063] Step 6.1. Initialize the azimuth-height imaging result X i (0) = 0 L×M , 0 L×M represents a zero matrix of dimension L×M, the iterative relative error ε, i = 1, 2,..., KN, and K, N, L, and M are the upsampling multiple, the total number of range samples, the total number of azimuth samples, and the total number of height samples initialized in Step 1, respectively.

[0064] Step 6.2. Perform the following steps in each iteration:

[0065] (1) Calculate Xf obtained by forward propagation using the following formula i (k),

[0066]

[0067] where Pc(i) is the upsampled pulse compression result after range migration correction corresponding to the i-th distance cell in Pc L×M×(KN) , i = 1, 2,..., KN, and Pc L×M×(KN) is the upsampled pulse compression result after range migration correction obtained in Step 4; F L and F MThey are the discrete Fourier transform matrices of points L and M described in Definition 8, denote F M transpose, denote F L conjugate transpose, denote F M conjugate; ⊙ denotes the matrix Hadamard product; X i (k - 1) is the result obtained from the (k - 1)-th iteration; Φ i is the inverse matrix of the phase compensation of the RMA operator corresponding to the current range cell calculated in Step 5.1; K, N, L, and M are the upsampling factor, the total number of range samples, the total number of azimuth samples, and the total number of height samples initialized in Step 1, respectively.

[0068] (2) Use the following formula to calculate the backpropagation Xb i (k)

[0069]

[0070] where denotes the conjugate of Φ i Φ i is the inverse matrix of the phase compensation of the RMA operator corresponding to the current range cell calculated in Step 5.1; F L and F M are the discrete Fourier transform matrices of points L and M described in Definition 8, denote F M transpose, denote F L conjugate transpose, denote F M conjugate; ⊙ denotes the matrix Hadamard product; Xf i (k) is obtained from the forward propagation calculation of the k-th iteration. i = 1, 2, …, KN, and K, N, L, and M are the upsampling factor, the total number of range samples, the total number of azimuth samples, and the total number of height samples initialized in Step 1, respectively.

[0071] (3) Use the following formula to calculate the gradient descent Xg i (k),

[0072] Xg i (k) = X i (k - 1) + Xb i (k)

[0073] where X i (k - 1) is the result obtained from the (k - 1)-th iteration, and Xb i (k) is obtained from the backpropagation calculation of the k-th iteration. i = 1, 2, …, KN, and K and N are the upsampling factor and the total number of range samples initialized in Step 1, respectively.

[0074] (4) Calculate the threshold shrinkage using the following formula

[0075] X i (k) = [firm(xg i (k; l, m); β, 2β)]

[0076]

[0077] where xg i (k; l, m) represents the (l, m)-th element of Xg i (k), Xg i (k) is obtained by calculating the gradient descent at the k-th iteration, β is the weight coefficient of the generalized MCP non-convex constraint term, l = 1, 2, …, L, m = 1, 2, …, M, |xg i (k; l, m)| represents the absolute value of xg i (k; l, m), and L and M represent the total number of azimuth samples and the total number of height samples, respectively. X i (k) is the result obtained at the k-th iteration. i = 1, 2, …, KN, where K and N are the upsampling factor and the total number of range samples initialized in step 1, respectively.

[0078] (5) Calculate the iteration error dX using the following formula i ,

[0079]

[0080] where X i (k - 1) is the result obtained at the (k - 1)-th iteration, and X i (k) is the result obtained at the k-th iteration, represents the square of the matrix Fibonacci norm. i = 1, 2, …, KN, where K and N are the upsampling factor and the total number of range samples initialized in step 1, respectively.

[0081] (6) Judge the iteration termination condition. If dX i ≤ ε, the iteration stops. dX i is the iteration error calculated in the previous step, and the azimuth-height imaging result X i corresponding to the i-th range cell is the result X i (k) obtained at the k-th iteration. Otherwise, continue the iteration. ε is the iteration relative error, i = 1, 2, …, KN, where K and N are the upsampling factor and the total number of range samples initialized in step 1, respectively.

[0082] Step 7. Stack along the range direction to obtain the overall 3D sparse imaging result of the target scene

[0083] For each azimuth-height imaging result X corresponding to the distance cell obtained in step 6 i , where i = 1, 2, …, KN, perform stacking operation along the range direction using the stacking operation method along the range direction defined in Definition 12 to obtain the final overall three-dimensional sparse imaging result X of the target; X is the final imaging result of the proposed millimeter-wave three-dimensional sparse imaging method based on the generalized MCP non-convex constraint. K and N are the upsampling factor and the total number of range samples obtained by initializing in step 1, respectively.

[0084] The innovation and advantages of the present invention are as follows: Different from the existing millimeter-wave three-dimensional sparse imaging method based on L1-norm constraint, a generalized MCP non-convex constraint is used to construct a sparse imaging equation, effectively avoiding the estimation deviation of strong scattering sources and the energy loss of weak scattering sources caused by the L1-norm convex constraint, and being able to obtain a more abundant imaging result of the target contour, which is beneficial to subsequent target detection and recognition in image aggregation; an RMA imaging operator is used to construct a sparse imaging sensing model, and the efficiency of the imaging operator is utilized to reduce the computational complexity of equation solving. Description of the Drawings

[0085] Figure 1 is the flowchart of the present invention.

[0086] Figure 2 is the verification result of the simulation experiment of the present invention. Detailed Embodiment

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

[0088] Step 1. Initialize relevant parameters

[0089] The speed of light c in air = 3×10 8 m / s; take the natural exponential function exp(·); the imaginary unit The pi π = 3.14; the wavelength of the transmitted signal λ = 0.0039 m; the azimuth sampling sequence number, denoted as l = 1, 2, …, 512; the azimuth sampling interval d l = 0.01 m; the elevation sampling sequence number, denoted as m = 1, 2, …, 512; the elevation sampling interval d m = 0.01 m; the total number of range samples N = 128; the upsampling factor K = 8; the reference distance R0 = 0.4 m; the range sampling interval d r = 0.03 m; the original target echo, denoted as S 512×512×128 ; the weight coefficient β of the generalized non-convex MCP constraint term = 0.5;

[0090] Step 2. Perform range-direction pulse compression processing on the target original echo to obtain the range-direction pulse compression result

[0091] Using the target original echo S in Step 1 512×512×128 as the input, adopt the pulse compression method defined by Standard 4 to compress the range-direction signals corresponding to each pair of azimuth-height-direction sampling points in S 512×512×128 to obtain the post-range-direction pulse compression result P 512×512×128 .

[0092] Step 3. Perform upsampling on the result after pulse compression

[0093] Using the post-range-direction pulse compression result P obtained in Step 2 512×512×128 and the upsampling multiple K initialized in Step 1 as the input, perform 8-fold upsampling processing. The post-range-direction pulse compression result corresponding to each pair of azimuth-height-direction sampling points in the post-pulse compression result P 512×512×128 is denoted as P lm , where lm represents the l-th azimuth direction - the m-th height direction sampling point pair, l = 1, 2,..., 512, m = 1, 2,..., 512. For each result p lm , perform the following operations to obtain the range-direction upsampled pulse compression result P' 512×512×1024 .

[0094] Step 3.1.

[0095] Adopt the fast Fourier transform (FFT) processing vector p described in Definition 6 lm to obtain the vector f lm .

[0096] Step 3.2.

[0097] Insert 896 zero elements starting from the 257th position of the vector f lm to obtain f' lm = [f lm (1, 2,... 256) 0 896 f lm (257,..., 512)], where f lm (1, 2,... 256) represents the first 256 elements in f lm , f lm (257,..., 512) represents the 257th element to the last element in f lm , 0 896 represents the inserted 896 zero elements, and f lm is the vector obtained by the fast Fourier transform processing in Step 3.1

[0098] Step 3.3.

[0099] Process the vector f′ using the inverse fast Fourier transform (IFFT) described in Definition 6 lm , to obtain the vector p′ lm .

[0100] Step 4. Range migration correction

[0101] Perform range migration correction on the range up-sampled pulse compression result P′ obtained in Step 3 using the range migration correction described in Definition 7 512×512×1024 to obtain the range up-sampled pulse compression result Pc after range migration correction 512×512×1024 .

[0102] Step 5. Construct the imaging solution equation for approximating the RMA imaging operator based on the generalized MCP non-convex constraint

[0103] Using the range up-sampled pulse compression result Pc after range migration correction obtained in Step 4 512×512×1024 as the input, construct the imaging solution equation for approximating the RMA imaging operator based on the generalized MCP non-convex constraint. The imaging solution equation corresponding to each range cell is constructed in the following step manner.

[0104] Step 5.1. Calculate the inverse matrix Φ of the RMA operator phase compensation corresponding to the current range cell i as follows

[0105]

[0106] where exp is the natural exponential initialized in Step 1.

[0107] Step 5.2. Construct the imaging solution equation corresponding to the current range cell as follows

[0108]

[0109] where denotes the X that minimizes , Pc(i) is the range migration correction up-sampled pulse compression result corresponding to the i-th range cell in Pc i , i = 1, 2,..., 1024, Pc 512×512×1024 is the range up-sampled pulse compression result after range migration correction obtained in Step 4; F 512×512×1024 and F 512 and F 512 are the 512-point and 512-point discrete Fourier transform matrices described in Definition 8 respectively, denotes the transpose of F 512 , denotes the conjugate transpose of F 512 . denotes F 512 the conjugate of; ⊙ denotes the matrix Hadamard product; denotes the square of the matrix Fibonacci norm; is the weight of the generalized MCP non-convex constraint term described in Definition 3, x i = vec(X i ), Φ i is the inverse matrix of the RMA operator phase compensation corresponding to the current range cell calculated in Step 5.1; denotes the matrix Kronecker product, diag(·) is the vector matrix diagonal operator described in Definition 10, and vec(·) is the matrix vectorization operator described in Definition 11. X i is the imaging solution equation corresponding to the i-th range cell constructed.

[0110] Step 6. Solve the imaging solution equation for the approximation of the RMA imaging operator based on the generalized MCP non-convex constraint

[0111] Taking the imaging solution equation corresponding to each range cell obtained in Step 5 as the input, solve the equation through an iterative method to obtain the azimuth-height imaging result corresponding to each range cell. Initialize the azimuth-height imaging result X i (0) = 0 512×512 , 0 512×512 denotes a all-zero matrix with dimensions of 512×512, the iterative relative error ε = 0.001, i = 1, 2,..., 1024. In each iteration, perform the following steps:

[0112] (1) Calculate the forward propagation to obtain Xf i (k), as follows,

[0113]

[0114] where Pc(i) is the range migration corrected and upsampled pulse compression result corresponding to the i-th range cell in Pc 512×512×1024 , i = 1, 2,..., 1024, Pc 512×512×1024 is the range migration corrected and upsampled pulse compression result obtained in Step 4; F 512 and F 512 are the 512-point and 512-point discrete Fourier transform matrices described in Definition 8, denotes the transpose of F 512 , denotes the conjugate transpose of F 512 , denotes the conjugate of F 512 ; ⊙ denotes the matrix Hadamard product; X i(k - 1) is the result obtained from the (k - 1)-th iteration; Φ i is the RMA operator phase compensation inverse matrix corresponding to the current distance cell calculated in step 5.1, where i = 1, 2, …, 1024.

[0115] (2) Calculate the backward propagation Xb i (k), as follows,

[0116]

[0117] where denotes the conjugate of Φ i ; Φ i is the RMA operator phase compensation inverse matrix corresponding to the current distance cell calculated in step 5.1; F 512 and F 512 are respectively the 512-point and 512-point discrete Fourier transform matrices described in Definition 8, denotes the transpose of F 512 ; denotes the conjugate transpose of F 512 ; denotes the conjugate of F 512 ; ⊙ represents the matrix Hadamard product; Xf i (k) is the result obtained from the forward propagation calculation in the k-th iteration, where i = 1, 2, …, 1024.

[0118] (3) Calculate the gradient descent Xg i (k), as follows,

[0119] Xg i (k) = X i (k - 1) + Xb i (k)

[0120] where X i (k - 1) is the result obtained from the k-th iteration, and Xb i (k) is the result obtained from the backward propagation calculation in the k-th iteration, where i = 1, 2, …, 1024.

[0121] (4) Calculate the threshold shrinkage as follows,

[0122] X i (k) = [firm(xg i (k; l, m); 0.5, 1)]

[0123]

[0124] where xg i (k; l, m) represents the (l, m)-th element of Xg i (k), and Xg i(k) is obtained from the k-th calculation of gradient descent, β is the weight coefficient of the generalized MCP non-convex constraint term, l = 1, 2, …, 512, m = 1, 2, …, 512, |xg i (k; l, m)| represents the absolute value of xg i (k; l, m), and i = 1, 2, …, 1024.

[0125] (5) Calculate the iteration error dX i , as follows,

[0126]

[0127] where X i (k - 1) is the result obtained from the (k - 1)-th iteration, and X i (k) is the result obtained from the k-th iteration. represents the square of the matrix Fibonacci norm, and i = 1, 2, …, 1024.

[0128] (6) Determine the iteration termination condition. If dX i ≤0.001, the iteration stops. dX i is the iteration error obtained from the previous step. The azimuth-height imaging result X i corresponding to the i-th range cell is the result X i (k) obtained from the k-th iteration. Otherwise, continue the iteration, and i = 1, 2, …, 1024.

[0129] Step 7. Stack along the range direction to obtain the overall 3D sparse imaging result of the target scene

[0130] Using the azimuth-height imaging results X i corresponding to each range cell obtained in Step 6, i = 1, 2, …, 1024 as the input, perform the stacking operation along the range direction according to the standard shown in Definition 12 to form the final overall 3D sparse imaging result X of the target. X is the final imaging result of the proposed millimeter-wave 3D sparse imaging method based on the generalized MCP non-convex constraint.

[0131] The results of computer simulation show that compared with the millimeter-wave 3D sparse imaging method based on the L1-norm constraint, the proposed method has a smaller estimation deviation, a more complete geometric contour retention, and a higher image quality.

Claims

1. A millimeter-wave three-dimensional sparse imaging method based on generalized MCP constraints, characterized in that it It includes the following steps: Step 1. Initialize relevant parameters The propagation speed of light in air is denoted as c; the natural exponential function is denoted as exp(·); the imaginary unit is denoted as j; the pi is denoted as π; the wavelength of the transmitted signal is denoted as λ; the azimuth sampling sequence number is denoted as l = 1, 2, …, L, where L represents the total number of azimuth samplings; the azimuth sampling interval is denoted as d l ; the elevation sampling sequence number is denoted as m = 1, 2, …, M, where M represents the total number of elevation samplings; the elevation sampling interval is denoted as d m ; the total number of range samplings is denoted as N; the upsampling multiple is denoted as K; the reference range is denoted as R0; the range sampling interval is denoted as d r ; the original target echo is denoted as S L×M×N ; the weight coefficient of the generalized non-convex MCP constraint term is denoted as β; Step 2. Perform range-direction pulse compression processing on the target original echo to obtain the range-direction pulse compression result Using the target original echo S in step 1 L×M×N as the input, the standard pulse compression method is used to compress the range signals corresponding to each pair of azimuth-height sampling points in S L×M×N to obtain the result P after range pulse compression L×M×N ; Step 3. Upsample the result after pulse compression Using the distance obtained in step 2 as the input together with the upsampling multiple K initialized in step 1, perform K-fold upsampling processing on the result P after range compression; for each pair of azimuth-height sampling points corresponding to the result P after range compression L×M×N in the azimuth direction, denote it as P L×M×N , where lm represents the l-th azimuth - m-th height sampling point pair, l = 1, 2, …, L, m = 1, 2, …, M, and L and M respectively represent the total number of azimuth samplings and the total number of height samplings; for each result p lm , perform the following operations to obtain the range upsampled pulse compression result P′ lm ; L×M×(Kn) ; Step 3.

1. Process vector f using the traditional Fast Fourier Transform (FFT) lm to obtain vector f lm ; Step 3.

2. Starting from the position of vector f lm Insert (K - 1)·N zero elements to obtain where denotes the first elements in f lm , denotes the th element to the last element in f lm , 0 denotes the (K - 1)·N zero elements inserted (K-1)·N , denotes rounding down, f lm is the vector obtained after the fast Fourier transform processing in step 3.1, N is the total number of range samples initialized in step 1, and K is the upsampling factor;​ Step 3.

3. Process the vector f′ using the traditional inverse fast Fourier transform (IFFT) lm to obtain the vector p′ lm ; Step 4. Range migration correction Perform range migration correction on the range upsampling pulse compression result P′ obtained in step 3 by using the traditional range migration correction method L×M×(KN) to obtain the range upsampling pulse compression result Pc after range migration correction L×M×(KN) , where K, N, L, and M are respectively the upsampling factor, the total number of range samples, the total number of azimuth samples, and the total number of height samples initialized in step 1; Step 5. Construct an imaging solution equation for approximating the RMA imaging operator based on the generalized MCP non-convex constraint Using the distance obtained in step 4 as the input to the upsampled pulse compression result Pc after migration correction, an imaging solution equation for approximating the RMA imaging operator based on the generalized MCP non-convex constraint is constructed; the imaging solution equation corresponding to each range cell is constructed in the following step-by-step manner; L×M×(KN) ​ Step 5.

1. Calculate the RMA operator phase compensation inverse matrix Φ corresponding to the current distance cell using the following formula i :[[]]END]] where i = 0, 1, 2, … (KN - 1), f l = 0, 1, 2, …, (L - 1), f m = 0, 1, 2, …, (M - 1), l = 1, 2, …, L, m = 1, 2, …, M; exp, j, π, λ, R0, K, N, L, M, d r 、d l and d m are respectively the natural exponent, imaginary unit, pi, wavelength of the transmitted signal, reference distance, upsampling factor, total number of range samples, total number of azimuth samples, total number of height samples, range sampling interval, azimuth sampling interval, and height sampling interval obtained by initialization in step 1; Step 5.

2. Construct the imaging solution equation corresponding to the current range cell as follows Among them denotes the X that makes the smallest i , Pc(i) is the result of upsampled pulse compression after range migration correction corresponding to the i-th range cell in Pc L×M×(KN) , i = 1, 2, …, KN, Pc L×M×(KN) is the result of upsampled pulse compression after range migration correction obtained in step 4; F L and F M are the discrete Fourier transform matrices of L points and M points respectively, denotes the transpose of F M , denotes the conjugate transpose of F L , denotes the conjugate of F M ; ⊙ represents the matrix Hadamard product; denotes the square of the matrix Fibonacci norm; β is the weight coefficient of the generalized MCP non-convex constraint term initialized in step 1; is the weight of the generalized MCP non-convex constraint term, x i = vec(X i ), Φ i is the inverse matrix of the RMA operator phase compensation corresponding to the current range cell calculated in step 5.1; denotes the matrix Kronecker product, diag(·) is the vector matrix diagonalization operator, and vec(·) is the matrix vectorization operator; X i is the imaging solution equation corresponding to the i-th range cell constructed; K, N, L, and M are the upsampling multiple, total number of range samples, total number of azimuth samples, and total number of height samples initialized in step 1 respectively; Step 6. Solve the imaging solution equation for approximating the RMA imaging operator based on the generalized MCP non-convex constraint Solve the imaging solution equation corresponding to each range cell obtained in Step 5 in an iterative manner to obtain the azimuth-height direction imaging result corresponding to each range cell; Step 6.

1. Initialize the azimuth-height imaging result X corresponding to the distance cell to be solved i (0) = 0 L×M , 0 L×M denotes a zero matrix of dimension L×M, the iterative relative error ε, i = 1, 2, …, KN, where K, N, L, and M are the upsampling factor, the total number of range samples, the total number of azimuth samples, and the total number of height samples obtained by initializing in Step 1, respectively; Step 6.

2. Perform the following steps in each iteration: (1) Use the following formula to calculate the forward propagation to obtain Xf i (k), where Pc(i) is the upsampled pulse compression result after range migration correction corresponding to the i-th range cell in Pc L×M×(KN) and i = 1, 2, …, KN. Pc L×M×(KN) is the upsampled pulse compression result after range migration correction obtained in step 4; F L and F M are the discrete Fourier transform matrices of L points and M points respectively, denotes the transpose of F M , denotes the conjugate transpose of F L , denotes the conjugate of F M ; ⊙ represents the matrix Hadamard product; X i (k - 1) is the result obtained in the (k - 1)-th iteration; Φ i is the inverse matrix of the RMA operator phase compensation corresponding to the current range cell calculated in step 5.1; K, N, L, and M are the upsampling factor, the total number of range samples, the total number of azimuth samples, and the total number of height samples initialized in step 1 respectively; (2) Calculate the backward propagation Xb using the following formula i (k) Among them represents the conjugate of Φ i , and Φ i is the inverse matrix of the RMA operator phase compensation corresponding to the currently calculated distance cell in step 5.1; F L and F M are the discrete Fourier transform matrices of L points and M points respectively, represents the transpose of F M , represents the conjugate transpose of F L , represents the conjugate of F M ; ⊙ represents the matrix Hadamard product; Xf i (k) is obtained from the forward propagation of the k-th iteration calculation; i = 1, 2, …, KN, and K, N, L, and M are the upsampling multiples, the total number of range samples, the total number of azimuth samples, and the total number of height samples initialized in step 1 respectively; (3) Calculate the gradient descent Xg i (k) using the following formula Xg i X(k) = i X(k - 1)+Xb i (k) Among them, X i (k - 1) is the result obtained from the first iteration, and Xb i (k) is obtained by backpropagation after the k-th iteration calculation; i = 1, 2,..., KN, where K and N are the upsampling factor and the total number of range samples initialized in step 1, respectively; (4) Calculate the threshold shrinkage using the following formula X i (k) = [firm(xg i (k; l, m); β, 2β)] where \(x_g\) i \((k; l, m)\) represents the \((l, m)\)-th element of \(x_g\) i \((k)\), \(x_g\) i \((k)\) is obtained by calculating the gradient descent at the \(k\)-th iteration, \(\beta\) is the weight coefficient of the generalized MCP non-convex constraint term, \(l = 1, 2, \ldots, L\), \(m = 1, 2, \ldots, M\), \(|x_g\) i \((k; l, m)|\) represents the absolute value of \(x_g\) i \((k; l, m)\), \(L\) and \(M\) respectively represent the total number of azimuthal samplings and the total number of elevation samplings; \(X\) i \((k)\) is the result obtained at the \(k\)-th iteration; \(i = 1, 2, \ldots, KN\), \(K\) and \(N\) are respectively the upsampling factor and the total number of range samplings obtained by initializing in step 1; (5) Calculate the iterative error dX using the following formula i , where X i (k - 1) is the result obtained from the (k - 1)-th iteration, X i (k) is the result obtained from the k-th iteration, represents the square of the matrix Fibonacci norm; i = 1, 2, …, KN, where K and N are the upsampling factor and the total number of range samples obtained by initializing in step 1, respectively; (6) Judgment of iteration termination condition. If dX i ≤ ε, the iteration stops. dX i is the iteration error obtained from the previous calculation. The azimuth-height imaging result X i corresponding to the i-th range cell is the result X i (k) obtained from the k-th iteration. Otherwise, continue the iteration; ε is the relative iteration error, i = 1, 2, …, KN, where K and N are the upsampling factor and the total number of range samples obtained from the initialization in step 1 respectively; Step 7. Stack along the range direction to obtain the overall three-dimensional sparse imaging result of the target scene For the azimuth-height imaging result X corresponding to each range cell obtained in step 6 i , where i = 1, 2, …, KN, perform a range stacking operation using the quasi-range stacking operation method to obtain the final overall three-dimensional sparse imaging result X of the target; X is the final imaging result of the proposed millimeter-wave three-dimensional sparse imaging method based on the generalized MCP non-convex constraint; K and N are the upsampling factor and the total number of range samples obtained by initializing in step 1, respectively.

Images