Three-dimensional millimeter wave imaging method based on laplacian prior sparse bayesian learning

Abstract

The application belongs to the technical field of millimeter wave radar, and discloses a three-dimensional millimeter wave imaging method based on Laplace prior sparse Bayesian learning, comprising the following steps: obtaining echo data s after distance focusing, constructing imaging operators and conjugate operators; when the termination condition is not met, iteratively performing the steps until the termination condition is met; and outputting a three-dimensional imaging result. The application avoids input of an additional regularization parameter, avoids construction, storage and large-scale matrix-vector multiplication in an optimization process of solving of an observation matrix, and uses an imaging operator constructed based on an MF imaging model to replace vectorization operation in a CS model, so that the calculation efficiency is improved.

Description

Technical Field

[0001] This invention belongs to the field of millimeter-wave radar technology, and particularly relates to a three-dimensional millimeter-wave imaging method based on Laplace prior sparse Bayesian learning. Background Technology

[0002] Millimeter-wave (MMW) imaging technology plays a crucial role in numerous fields, including security detection, driver assistance systems, medical imaging, security monitoring, and concealed target detection, gradually becoming one of the key technologies for intelligent and efficient perception and monitoring [Millimeter-Wave SAR Imaging With Radar Networks Based on Radar Self-Localization][Millimeter-Wave SAR Sparse Imaging With 2-D Spatially Pseudorandom Spiral-Sampling Pattern]. Compared to traditional radar imaging technology, three-dimensional millimeter-wave imaging technology achieves three-dimensional resolution of the imaging space by expanding the aperture dimension and combining it with pulse compression technology. However, with the improvement of imaging performance, imaging methods based on matched filtering theory (MF), such as back projection algorithm (BPA) and range migration algorithm (RMA), are often limited by the Nyquist sampling criterion, leading to an increased burden on the acquisition, storage, and transmission of large amounts of echo data. In cases of insufficient or undersampled data, these methods may produce problems such as high sidelobes, aliasing, and artifacts, affecting tasks such as image interpretation and applications [A 3D Sparse SAR Imaging Method Based on Plug-and-Play][RMIST-Net: Joint Range Migration and Sparse Reconstruction Network for 3D mmWImaging].

[0003] Sparse imaging algorithms based on compressed sensing (CS) have attracted widespread attention for reducing sampling data while maintaining image quality [SparseArray Optimization Using Simulated Annealing and Compressed Sensing for Near-Field Millimeter Wave Imaging]. These methods typically transform the radar imaging problem into an underdetermined linear equation system, combining prior information such as sparsity to construct a joint constraint optimization model, and then optimizing the imaging results through iterative updates. For example, Zeng et al., in their research [Sparse Regularization: Convergence of Iterative Jumping Thresholding Algorithm], explored the convergence and rate of a non-convex iterative jump thresholding algorithm in solving sparse regularization optimization problems. This research proposed an iterative jump thresholding algorithm (IJT), which demonstrated good imaging performance in sparse imaging applications of ship targets. Yang et al. proposed a near-end separation algorithm based on a nonconvex Cauchy penalty function, demonstrating good convergence characteristics and SAR image reconstruction performance [On Solving SAR Imaging Inverse Problems Using Nonconvex Regularization With a Cauchy-Based Penalty][Convergence Guarantees for Non-Convex Optimisation With Cauchy-Based Penalties]. Yang et al. studied the CS method's effect on sparse measurement image reconstruction in broadband microwave SAR systems, verifying its effectiveness at low sampling rates. Wang et al. used the latest denoising techniques to replace the thresholding step in the Alternating Direction Multiplier Method (ADMM), improving the quality of 3D MMW images [Plug-and-play priors for model-based reconstruction][A 3D Sparse SAR Imaging Method Based on Plug-and-Play]. Chen et al., based on the ADMM framework, optimized the ADMM iteration steps according to the updated augmented Lagrange function, avoiding large-scale matrix inversion and effectively reducing computational complexity [Efficient MMW Image Reconstruction Algorithm Based on ADMM Framework for Near-Field MIMO-SAR]. Ge et al. modeled the single-bit SAR imaging problem as a sparse logistic regression optimization problem, and conducted research on SAR sparse imaging using single-bit measurements [SparseLogisticRegressionBasedOne-BitSARImaging]. Furthermore, Wang et al., based on the theoretical framework of sparse learning imaging, proposed several network models focusing on 3D imaging model representation, network construction, and error compensation [Efficient ADMM Framework Based on Functional Measurement Model for mmW 3D SAR Imaging][3D SAR Autofocusing With Learned Sparsity].

[0004] The aforementioned sparse reconstruction-based imaging methods typically require regularization parameters as algorithm inputs, and the selection of these parameters significantly impacts imaging performance. Furthermore, sparse learning-based imaging methods often rely on dataset construction, but obtaining a large number of real-world echo data for network training is a challenge. Sparse Bayesian Learning (SBL), a classic method in Computer Science, can model signals based on probabilistic models, more accurately capturing sparse structures within signals. It also possesses self-regularization properties, automatically adjusting model parameters to prevent overfitting and improve the model's generalization ability. With the deepening of SBL theoretical analysis, its advantages in radar imaging are gradually becoming apparent. For example, Zhao et al. studied the phase error correction problem in the context of sparse signal recovery [An Autofocus Technique for High-Resolution Inverse Synthetic Aperture Radar Imagery], using a multi-task Bayesian model for hierarchical modeling to obtain focused high-resolution radar images. Zhou et al. proposed a SBL method based on conjugate hierarchical priors for sparse subband imaging [High-Resolution Sparse Subband Imaging Based on Bayesian Learning With Hierarchical Priors]. Liu et al. proposed a fully automated ISAR imaging algorithm based on SBL. Experimental results based on simulation and measurement data show that this algorithm maintains a better balance between computational load and reconstructed signal than existing algorithms [Superresolution ISAR Imaging Based on Sparse Bayesian Learning]. Wang et al. proposed a fast SBL algorithm suitable for high-resolution imaging with sparse apertures and achieved better noise robustness by using the minimum Tsallis entropy algorithm to achieve autofocusing [Sparse Aperture Autofocusing and Imaging Based on FastSparse Bayesian Learning From Gapped Data]. Xu proposed a locally structured SBL algorithm by utilizing the joint sparse modes of adjacent scatterers. Structured priors are utilized by modeling adjacent correlations or dependencies to encode joint sparse modes. Meanwhile, a parameter dictionary with unknown rotation parameters was constructed to represent the target's maneuverability, further improving imaging performance [Enhanced ISAR Imaging and Motion Estimation With Parametric and Dynamic Sparse Bayesian Learning].However, current sparse imaging methods based on SBL are mainly concentrated in the field of ISAR imaging. In contrast, in 3D MMW imaging, there is usually a large amount of echo data, and due to the large number of pixels in the imaging scene, the problem of constructing a large-scale observation matrix when using the SBL method is more prominent. Therefore, research on 3D MMW sparse imaging methods based on SBL is relatively limited. Summary of the Invention

[0005] Three-dimensional millimeter-wave (MMW) sparse imaging technology based on compressed sensing (CS) can achieve high-quality scene reconstruction using undersampled echo data. Sparse Bayesian learning (SBL), a classic signal recovery method in CS, effectively promotes the sparsity of solutions by modeling sparse signals based on probabilistic models. Considering that the Laplace distribution is a peaked distribution, it exhibits excellent performance in promoting sparse solutions, reducing overfitting, and demonstrating strong robustness in handling outliers.

[0006] Compared to traditional CS imaging methods, SBL methods provide posterior probability estimates for solutions and utilize higher-order statistical knowledge to improve the reconstruction error propagation problem during model iteration, demonstrating significant application value in the field of sparse imaging [The Variational Approximation for Bayesian Inference]. Therefore, for the 3D MMW sparse imaging problem, this application constructs a hierarchical Bayesian framework with an adaptive Laplacian prior and proposes a novel SBL method (LSBL). Theoretical analysis shows that the proposed model is equivalent to the reweighted l1 norm problem, but utilizes the self-regularization property of SBL to avoid additional regularization parameters. Then, a variational inference strategy is used to iteratively optimize the proposed method. To avoid the construction, storage, and large-scale matrix-vector multiplication in the optimization process of the observation matrix, this application also explores imaging operators based on MF imaging model construction to replace the vectorization operations in the CS model and introduces them into the proposed SBL method framework. This improvement effectively integrates the advantages of SBL and MF methods in imaging, significantly reducing the computational and storage costs of sparse imaging methods that directly apply SBL. Finally, through three-dimensional MMW sparse imaging experiments using simulated and measured data, the ability of the proposed method to achieve efficient image reconstruction using a small number of echoes was verified.

[0007] To achieve the above objectives, this application discloses a three-dimensional millimeter-wave imaging method based on Laplace prior sparse Bayesian learning, comprising the following steps:

[0008] Acquire the echo data s after range focusing, and construct the imaging operator. and conjugate operator

[0009] in,

[0010]

[0011] Where z represents the distance from the array plane to the target along the range direction, and k z For the corresponding wavenumber; IFFT 2D Represents the two-dimensional inverse fast Fourier transform, FFT 2D Θ represents the two-dimensional Fast Fourier Transform; e is the natural constant, j is the imaginary unit, ⊙ represents the Hadamard product; Θ represents the three-dimensional imaging result;

[0012] If the termination condition is not met, perform the following steps:

[0013] Calculate ∑ Θ =(LγI+diag(α) -1 )) -1 α is the variance (α is the variance when Gaussian modeling Θ, and can also be regarded as a non-negative hyperparameter controlling the sparsity of Θ in the 3D imaging result), L is the Lipschitz constant, γ is the noise variance, and ∑ Θ It is the covariance matrix of the posterior distribution;

[0014] Calculate the estimated value of the three-dimensional imaging result Θ δ is an intermediate variable in the iteration, and its estimated value is μ. Θ Initialize to And it changes with iteration;

[0015] The estimated value μ of the three-dimensional imaging result Θ Θ Assign the value to Θ;

[0016] calculate α follows a gamma distribution, where β is the scale parameter;

[0017] calculate β follows a gamma prior distribution, where a and b are the shape and scale parameters of the β gamma prior distribution;

[0018] calculate M represents the amount of echo data, c and d are hyperparameters of γ, and g(Θ,δ) is a relaxation function with respect to Θ and δ.

[0019] The estimated value μ of the three-dimensional imaging result Θ Θ Assign the value to δ;

[0020] Until the termination condition is met;

[0021] Output three-dimensional imaging result Θ=μ Θ .

[0022] Preferably, within the Bayesian modeling framework, all unknown variables are treated as random variables and assigned a prior distribution. A three-level hierarchical model is adopted, wherein the first level models the elements in the imaging result Θ according to the Gaussian distribution with a mean of 0 and a variance of α; the second level models the elements of α using the gamma distribution with a scale parameter of β; and the third level assigns a gamma prior to β with shape and scale parameters a and b, respectively, and models the imaging result Θ to obtain the Laplace prior form.

[0023] Specifically, this includes: imaging models considering distance focusing:

[0024] s=ΦΘ+w

[0025] Where Φ is the observation matrix constructed based on the imaging geometry, and w is noise;

[0026] First order: The elements of Θ follow an independent Gaussian distribution:

[0027]

[0028] Where α = [α1, α2, ..., α N ] T ;α n It is the nth element of α, N is the number of elements in α, and p() represents the conditional probability density function;

[0029] Second order: For the nth element α n Described using a gamma distribution:

[0030]

[0031] Where β = [β1, β2, ..., β 2N ] T ,β n It is the nth element of β, and:

[0032]

[0033] Third order: β is endowed with a gamma prior, that is

[0034]

[0035] The posterior Laplace distribution is:

[0036]

[0037] Where, Θ n Let Θ be the nth element of vector Θ;

[0038] The synthesized Laplace prior form is as follows:

[0039]

[0040] The noise w follows a complex mean of zero and a variance of γ. -1 The Gaussian distribution of is given by:

[0041]

[0042] In practical radar systems, the noise variance γ -1 It is unknown, and the probability density function of γ is described using a gamma distribution:

[0043] p(γ)=Γ(γ∣c,d)=Γ(c) -1 d c γ c-1 e -dγ

[0044] Based on the imaging model and the probability density function of the noise, the likelihood is described as follows:

[0045]

[0046] Maximum a posteriori estimation corresponds to the adaptive reweighted l1 norm problem:

[0047]

[0048] Solve the following optimization problem to obtain the maximum a posteriori estimate:

[0049]

[0050] in, Assign different weights related to the data to the N distinct elements in Θ.

[0051] Preferably, the regularization parameter λ is estimated. n It is estimated through a variational inference iterative process.

[0052] Preferably, the variational inference iterative process includes:

[0053] use Representing variables in a hierarchical model, the goal of Bayesian inference is to find hidden variables in echo data s. posterior distribution Approximately:

[0054]

[0055] Where q Θ (Θ), q α (α), q β (β) and q γ(γ) denote the posterior distribution functions in factorized form with respect to variables Θ, α, β, and γ, respectively;

[0056] Maximizing a lower bound L(q) of the log-likelihood lnp(s) by using Each variable in the process is implemented alternately, which leads to:

[0057]

[0058] Where const represents a constant. express A subset, express China regarding The complement;

[0059] For f(Θ) = ||s - ΦΘ|| 2 And for any Θ, δ, we have:

[0060] f(Θ)≤g(Θ,δ)=||s-Φδ|| 2 +2(Θ-δ) T Φ T (Φδ-s)+L||Θ-δ|| 2

[0061]

[0062] and

[0063] Where g(Θ,δ) is a relaxation function constructed with respect to Θ and δ given s, Φ and L;

[0064] Let g(Θ,δ) represent the likelihood function constructed based on g(Θ,δ), and φ represent the observation matrix constructed based on the imaging geometry.

[0065] p(s|Θ,γ) represents the likelihood description of s if and only if Θ=δ;

[0066] Each hidden variable is calculated using an alternating optimization approach. The posterior distribution q Θ (Θ)q α (α)q β (β)q γ The approximate value of (γ).

[0067] Preferably, each hidden variable is calculated using an alternating optimization method. The posterior distribution q Θ (Θ)q α (α)q β (β)q γ Approximate values ​​for (γ) include:

[0068] q Θ Update of (Θ): Ignore terms irrelevant to Θ, approximate posterior distribution q Θ (Θ) is calculated in the following way:

[0069]

[0070] Where E(.) is the expectation operator, q Θ (Θ) follows a Gaussian distribution, and its mean and covariance matrices are as follows:

[0071] μ Θ =γ∑ Θ (Lδ+Φ T s-Φ T Φδ)

[0072] ∑ Θ =(LγI+diag(α) -1 )) -1

[0073] The update formula for Θ is:

[0074] Θ=μ Θ

[0075] q α (α) update: approximate posterior q α (α) yielded the following results:

[0076]

[0077] in Represents ∑ Θ The nth diagonal element is calculated using the block coordinate descent method:

[0078]

[0079] Equaling this expression to zero results in:

[0080]

[0081] Update q β (β): The approximate posterior of β is:

[0082]

[0083] This indicates that β follows the parameters a+1 and Gamma distribution:

[0084]

[0085] Therefore, β nThe update formula is:

[0086]

[0087] q γ Update (γ): for q γ (γ) is subjected to variational optimization, resulting in:

[0088]

[0089] Therefore, γ follows the parameter and Gamma distribution:

[0090]

[0091] Update of δ: Find an estimate of δ through the following optimization:

[0092] Taking the derivative of g(Θ,δ) with respect to δ and setting it to 0, we get:

[0093] E((LI-Φ T Φ(δ-Θ))=0

[0094] And: δ=Θ

[0095] The imaging result Θ of the range plane z = z0 is obtained through the above iterative optimization process.

[0096] Preferably, based on the imaging result Θ of the distance to the z = z0 plane, all reconstructed N... z Stacking 2D slices along the range direction yields a 3D imaging result; then, on N... z When performing 2D imaging of a range slice, this is achieved by replacing the exact observation function in the CS-based sparse imaging framework with approximate observations obtained from the MF-based imaging process:

[0097] After distance focusing, s represents the corresponding echo data, and the reconstructed 3D image is denoted as Θ. Zero-padding is applied to both s and Θ, and then a 2D Fourier transform of the appropriate size is performed to obtain:

[0098]

[0099] Where, k z It corresponds to the wavenumber component in the distance direction, i.e., the z-axis direction, k z The values ​​of z and z vary with different imaging slices in the distance direction. The imaging operator M(·) and its inverse operator corresponding to the above equation are used. It approximates the large-scale matrix-vector multiplication operation, thus directly performing three-dimensional imaging on the range-focused echo data s.

[0100] Preferably, the termination condition is that the maximum number of iterations is reached or that the condition is met during the t-th iteration. ε = 10 -1 This represents a minimum value, with a maximum number of iterations of 20. Attached Figure Description

[0101] Figure 1 The geometric configuration of the imaging system is shown, where black dots represent active acquisition elements and white dots represent missing acquisition elements.

[0102] Figure 2 This is a graphical model of the hierarchical Bayesian framework of this application;

[0103] Figure 3 The left side shows the simulated complex target shape, and the right side shows the imaging scene of the simulated experiment;

[0104] Figure 4 The images show the imaging results of a simulated scene with SR=40%, from left to right: the results of RMA, IJT, Cauchy, G-SBL, and LSBL algorithms, respectively.

[0105] Figure 5 The images show the results of imaging a wrench target with an SR of 40%, from left to right: the results of the RMA, IJT, Cauchy, GSBL, and LSBL algorithms, respectively. Detailed Implementation

[0106] The present invention will be further described below with reference to the accompanying drawings, but this is not intended to limit the present invention in any way. Any modifications or substitutions made based on the teachings of the present invention shall fall within the protection scope of the present invention.

[0107] Before introducing the embodiments of this application, the symbols involved in this application will be explained.

[0108] Table 1 List of Important Symbols

[0109]

[0110]

[0111] Consider as Figure 1 The millimeter-wave imaging system shown is based on a single-state array. The transmitting and receiving antennas are located at the same position and scanned using electrical or mechanical control techniques to form a two-dimensional planar aperture. Here, X, Y, and Z represent the azimuth, elevation, and range directions, respectively. The number of grid points in the horizontal and vertical directions of this planar aperture are Nv and Nh, respectively, and the number of sampling frequency points at each location is N. k Let P a = (x', y', 0) is the position coordinate of the antenna, and the scatterer is at the target position P. tThe complex reflection coefficient at (x,y,z) is denoted as Θ(x,y,z).

[0112] For 3D imaging tasks, the number of pixels in the reconstructed image along the x-, y-, and z- directions is set to N, respectively. x N y N z By stacking all reconstructed N values ​​along the distance axis z 2-DN x ×N y Slicing allows for the acquisition of three-dimensional imaging results. Assuming the antenna transmits a frequency-modulated continuous wave signal, and considering round-trip amplitude attenuation, after range focusing, the echo data received by the planar target at a distance z = z0 can be represented as:

[0113]

[0114] Where k = (2πf) / c is the wave number, f is the frequency, c represents the speed of light, and R represents the distance from the antenna to the point scatterer, defined as:

[0115]

[0116] The amplitude coefficient R in equation (1) -2 Approximately R -1 The resulting approximation error is negligible in near-field applications. 2D image reconstruction will depend on the phase term e. -j2kR and amplitude coefficient R -1 The combination of .

[0117]

[0118] According to the fixed-phase method (MSP), a spherical wave can be represented as a superposition of plane waves. Neglecting constant terms, we can obtain:

[0119]

[0120] Where k corresponds to the wavenumber components x, y, and z, respectively. x ,k y and k z satisfy:

[0121]

[0122] Substituting (4) into (3) and rearranging the order of the integrals, we can write:

[0123]

[0124] Among them, FT 2D (·) and IFT 2D(·) represent the two-dimensional Fourier transform and the inverse Fourier transform, respectively. Considering the extremely short pulse duration, k can be regarded as a constant term of the pulse duration. Therefore, for simplicity, the integrals related to k are omitted, and the imaging formula for the distance plane z0 is obtained as follows:

[0125]

[0126] Sparse imaging techniques can reduce system complexity and decrease the amount of echo data by efficiently reconstructing scenes under limited measurement conditions. Assuming sparse sampling, κ (0 < κ < 1) represents the sparse sampling rate. Figure 1 This sparse sampling pattern is illustrated, where black dots represent antennas transmitting and receiving signals at corresponding locations, while white dots indicate that the antenna does not transmit or receive signals at those locations. However, using a sparse sampling pattern may not satisfy the Nyquist sampling criterion, significantly reducing the imaging quality of imaging methods based on MF theory. Representing the 2-D image reconstructed from the range plane z = z0 as Θ, after vectorizing Θ and s, the corresponding linear model can be constructed as follows:

[0127] s=ΦΘ+w (8)

[0128] Where s, Θ, and Φ are observation matrices determined by the imaging geometry. According to (1), the elements of Φ depend on the wavenumber k and the distance R from the sampling point to the point scatterer, i.e., they are composed of exp(-j2kR). w represents the noise vector. For equation (8), a l0 regularization term can be introduced to promote the sparse solution of the model:

[0129]

[0130] Here, λ>0 is an adjustable regularization parameter used to control the trade-off between sparsity and reconstruction fidelity. Directly solving for the l0 norm is NP-hard, and it is often transformed into an l1 norm optimization problem.

[0131]

[0132] Furthermore, by using a series of data-related weights instead of a single ordinary weight λ, we obtain the adaptive reweighted L1 problem, as shown below:

[0133]

[0134] Where [λ1,λ2,...,λ] N ] T Assigning different weights to the N distinct elements in Θ can achieve more accurate recovery performance. However, these models typically require regularization parameters as input to the algorithm. The choice of parameters has a significant impact on the algorithm's performance, and obtaining accurate parameter inputs is often challenging in practical applications.

[0135] Since its inception, the SBL method has been widely applied to sparse signal recovery and other fields. By introducing the sparsity of the prior distribution, the model can automatically select important features, preventing overfitting while maintaining interpretability of the observed data. In particular, compared with the common Gaussian distribution, the Laplace distribution, as a peaked distribution, shows great potential in promoting the generation of sparse solutions, reducing overfitting, and exhibiting strong robustness in handling outliers. Based on these advantages, this application proposes a hierarchical SBL model (LSBL) based on adaptive Laplace prior for 3D MMW imaging.

[0136] Within the Bayesian modeling framework, all unknown variables are treated as random variables and assigned appropriate prior distributions. This application employs a three-layer hierarchical model to model the range-focused echo data. For example... Figure 2 The graphical model shown depicts the model structure in detail.

[0137] Imaging Result Θ: To construct a Bayesian model with a Laplace prior, a three-layer framework was established for the unknown signal Θ. The first layer assumes that the elements of Θ follow independent Gaussian distributions:

[0138]

[0139] Where α = [α1, α2, ..., α N ] T Θ represents variance, which is a non-negative hyperparameter that controls the sparsity of Θ.

[0140] Second order, for each element α of variance α n This can be described using a gamma distribution:

[0141]

[0142] Where β = [β1, β2, ..., β 2N ] T ,and

[0143]

[0144] For the third level of the hierarchy, β is assigned a gamma prior: that is

[0145]

[0146] Where a and b are the shape and scale parameters of the β-gamma distribution, in this embodiment, a = b = 10. -6 .

[0147] Combining equations (13)-(15), we obtain the prior Laplace distribution:

[0148]

[0149] Noise (w): Assume w follows a complex mean of zero and a variance of γ. -1 The Gaussian distribution of has a probability density function (PDF) as follows:

[0150]

[0151] Noise variance (γ) -1 In practical radar systems, the noise variance γ -1 It is unknown; this application describes the probability density function (PDF) of γ using a gamma distribution:

[0152] p(γ)=Γ(γ∣c,d)=Γ(c) -1 d c γ c-1 e -dγ (18)

[0153] c and d are hyperparameters of γ. In this embodiment, we set c = 10. -6 and d=10 -6 .

[0154] Based on the observation model and the noise model (17), the likelihood can be described as:

[0155]

[0156] Where M represents the amount of echo data, by assigning the prior (8) to (19), it can be observed that the maximum a posteriori (MAP) estimation corresponds to the adaptive reweighted l1 norm problem.

[0157]

[0158] Solve the following optimization problem to obtain the maximum a posteriori estimate:

[0159]

[0160] in, Different weights related to the data are assigned to each element of Θ. In the CS model of equations (10)-(12), it is necessary to estimate the regularization parameter λ, which is difficult in practical applications. In the SBL method, these variables can be estimated through an iterative process. Next, this application calculates all variables in the model based on variational inference.

[0161] let The variables in the hierarchical model described above are used to find hidden variables in the observed data s. posterior distribution However This is typically not possible to calculate analytically. Considering that the marginal probability of s can be decomposed into:

[0162] lnp(s)=F(q)+KL(q||p) (22)

[0163] in

[0164]

[0165] and

[0166]

[0167] in It is an arbitrary probability density function. KL(q||p) is and The Kullback-Leibler divergence between them, since KL(q||p)≥0, therefore we have:

[0168] lnp(s)≥F(q)

[0169] That is, F(q) is a lower bound of the log-likelihood lnp(s), and the equality holds if and only if KL(q||p)=0, in which case:

[0170]

[0171] Therefore, maximizing F(q) is equivalent to minimizing KL(q||p), and the posterior distribution... F(q) can be approximated by the variational distribution q(θ). Its factorized form is:

[0172]

[0173] Maximizing F(q) can be done on Each variable in the process is implemented alternately, which leads to:

[0174]

[0175] Where const represents a constant. express A subset, express China regarding The supplement to .

[0176] Using the smoothing function f(Θ)=||s-ΦΘ|| 2 The continuous differentiability of F(q) can be used to relax the lower bound F(q) according to the following lemma, which can improve computational efficiency.

[0177] Lemma 1 [A fast iterative shrinkage-thresholding algorithm for linear inverse problems]: For f(Θ) to be a continuously differentiable function with a Lipschitz continuous gradient and a Lipschitz constant of L, then for any u, υ:

[0178]

[0179] Where L=λ max (φ T φ)+τ, τ>0 is a constant. Using Lemma 1, for any... have:

[0180] f(Θ)≤g(Θ,δ)=||s-Φδ|| 2 +2(Θ-δ) T Φ T (Φδ-s)+L||Θ-δ|| 2 (29)

[0181] Combined with equation (19), let

[0182]

[0183] have

[0184]

[0185] The equality holds if and only if Θ = δ;

[0186] g(Θ,δ) is a relaxation function constructed with respect to Θ and δ, given s, Φ, and L.

[0187] Let g(Θ,δ) represent the likelihood function constructed based on g(Θ,δ), and Φ represent the observation matrix constructed based on the imaging geometry.

[0188] make

[0189]

[0190] The lower bound F(q) can be further relaxed as follows:

[0191]

[0192] in, It is a normalization term, making It becomes a strict distribution. Therefore, the posterior distribution By maximizing It is approximated by the variational distribution q(θ).

[0193] Next, each hidden variable is computed using an alternating optimization approach. Approximate posterior distribution:

[0194]

[0195] 1)q Θ Update of (Θ): Ignore terms irrelevant to Θ, approximate posterior distribution q Θ (Θ) can be calculated in the following way:

[0196]

[0197] That is, q Θ (Θ) follows a Gaussian distribution, and its mean and covariance matrices are as follows:

[0198] μ Θ =γ∑ Θ (Lδ+Φ T s-Φ T Φδ) (34)

[0199] ∑ Θ =(LγI+diag(α) -1 )) -1 =1 / (LγI+diag(1 / α)) (35)

[0200] Therefore, ∑ Θ Since it is a diagonal matrix, ∑ is calculated during the iterative update process. Θ It has low complexity, and the update formula for Θ is:

[0201] Θ=μ Θ (36)

[0202] 2)q α (α) update: approximate posterior q α (α) yields the following result:

[0203]

[0204] in Represents ∑ Θ The nth diagonal element. Since the specific distribution is not satisfied in equation (37), this application uses the block coordinate descent (BCD) method for calculation:

[0205]

[0206] Equaling this expression to zero will result in:

[0207]

[0208] 3) Update q β (β): The approximate posterior of β is:

[0209]

[0210] This indicates that β follows the parameters a+1 and Gamma distribution:

[0211]

[0212] Therefore, β n The update formula is:

[0213]

[0214] 4)q γ Update (γ): Similarly, for q γ (γ) By performing variational optimization, we can obtain:

[0215]

[0216] Therefore, γ follows the parameter and Gamma distribution:

[0217]

[0218] The updated formula is:

[0219]

[0220] in

[0221]

[0222] 5) δ update: Find an estimate of δ through the following optimization:

[0223]

[0224] Taking the derivative of g(Θ,δ) with respect to δ and setting it to 0, we get:

[0225] E((LI-Φ T Φ)(δ-Θ))=0

[0226] as well as:

[0227] δ=Θ

[0228] The above iterative optimization process can obtain the imaging result Θ of the range direction z = z0 plane. Then, all the reconstructed 2-D slices in the range direction are stacked to obtain the three-dimensional imaging result.

[0229] However, 2D imaging requires constructing an observation matrix. This operation is not only time-consuming but also consumes a significant amount of memory. In contrast, traditional MF-based imaging processes avoid the construction and storage of Φ. Therefore, effectively integrating the advantages of CS and MF methods in imaging can significantly reduce the computational and storage costs of directly applying CS-based sparse imaging methods. This idea is achieved by replacing the exact observation function in the CS-based sparse imaging framework with approximate observations obtained from traditional MF-based imaging processes.

[0230] Taking the RMA algorithm as an example, assuming that after distance focusing, s is the echo data, and the reconstructed 3D image is represented as Θ, s and Θ are zero-paddingd and then used to perform a 2D Fourier transform of the corresponding size. Therefore, (7) and (8) can be expressed as:

[0231]

[0232] and

[0233]

[0234] Where k z And z varies with distance to different slices. This application adopts the imaging operators corresponding to equations (38) and (39). and inverse operators Approximate replacement of large-scale matrix-vector multiplication operations Φ H The operations (·) and Φ(·), where Replace Φ H (·) operation, conjugate operator Instead of the Φ(·) operation, this substitution improves computational efficiency and reduces memory usage during the imaging process, enabling direct 3D imaging of the range-focused 3D complex-valued echo data s without the need to separately construct the observation matrix and process N. z A 2-D slice along the range direction. Given that using the imaging operator instead of the observation matrix is ​​an approximation strategy, this application omits Φ in equation (1) during the iteration process. T The (Φ(·)) operation is performed, and L is initialized to 1, and Φ T The calculation of Φ(δ) is simplified to δ. This application summarizes the proposed LSBL method in Algorithm 1. The iteration termination condition is either reaching the maximum number of iterations or satisfying the condition at the t-th iteration. ε = 10 -1 This represents a minimum value.

[0235] Algorithm 1:

[0236] Input: Range-focused echo data s, construct imaging operator and conjugate operator Initialization: a, b, c, d, τ = 10 -6 L = 1, hyperparameters α = 1, β = 1, γ = 1, and If the termination condition is not met, execute:

[0237] 1)∑ Θ =(LγI+diag(α) -1 )) -1 ;

[0238] 2)

[0239] 3)

[0240] 4)

[0241] 4)

[0242]

[0243] 6) δ=μ Θ

[0244] 5) Until the termination condition is met.

[0245] Output: 3D imaging results

[0246] To verify the effectiveness of the proposed algorithm in 3D-MMW imaging, this application provides imaging experimental results based on simulation and measured data. This application compares the proposed method with novel Iterative Jump Thresholding (IJT), Non-convex Cauchy Regularization (Cauchy), and efficient SBL method with Gamma super-prior (GSBL). Furthermore, the imaging performance of the RMA method based on MF theory is demonstrated. It is worth noting that both the IJT and GSBL methods require the calculation of the Lipschitz constant L, and the accurate estimation of L has a significant impact on the performance and convergence of these methods. Considering the long computation time required for eigenvalue decomposition of large imaging scenes, this application employs the power method for fast computation [Numerical methods for scientists and engineers]. For the RMA and Fast LSBL methods, sparse sampling is achieved by setting the echo data at randomly selected missing elements to zero. To quantitatively evaluate and compare the imaging performance of the proposed method and contrast algorithms, image entropy (IE) and contrast (IC) are introduced for evaluation. The definition of IE is as follows:

[0247]

[0248] A smaller IE value reflects better image focusing quality. IC is calculated as follows [SequentialFrequency-Domain Imaging Algorithm for Near-Field MIMO-SAR With Arbitrary Scanning Paths]:

[0249]

[0250] Here, E(·) represents the averaging operator. A large value of IC indicates that the dominant coefficients are clustered in the target region, while artifacts and clutter are effectively suppressed. Furthermore, since CS-based imaging methods require only a small number of measurements to achieve high-precision imaging of sparse scenes, SR (sampling ratio) is used to measure the proportion of missing elements, defined as:

[0251]

[0252] This application uses the Terahertz Imaging Toolbox [Terahertz Imaging Toolbox with Interactive User Interface] to generate simulated data and conduct imaging experiments to verify the proposed method and comparative methods. The center frequency and bandwidth of the antenna are set to 79 GHz and 4 GHz, respectively. For the array element configuration, this application sets the number of horizontal and vertical sampling positions to 192, and sets 50 sampling points for each sampling position, thus shaping the echo data into a cube with a size of 192 × 192 × 50. Figure 3 As shown, to simulate the imaging performance of the proposed method for complex targets, cuts of different shapes were added inside the rectangular target, which consists of 11,748 scattering points and is 0.5m away from the array plane. For the IJT, Cauchy, and GSBL methods, the imaging scene is divided into a 200*200 imaging grid to construct the observation matrix.

[0253] Table 2 shows the quantitative evaluation and time consumption of all methods for simulation experiment imaging results.

[0254]

[0255] This application considers simulation experiments for cases with SR of 40% and 60% to evaluate the imaging performance of the sparse imaging method under undersampling conditions. Figure 4The imaging results of various methods under different SR conditions are compared. Observations show that in undersampling conditions, the RMA method can delineate the main outline of the target, but has a large number of side lobes. When SR = 40%, the Cauchy method still has some artifacts; however, when SR increases to 60%, the imaging quality of the Cauchy method improves, and the suppression of artifacts is significantly enhanced. The IJT method performs well in focusing the target, but some artifacts still exist at different shaped cuts within the rectangular target, and this phenomenon does not improve significantly with increasing SR. In contrast, the SBL-based GSBL method and the proposed method exhibit a clearer background, effectively delineate the different shaped cuts within the rectangular target, and have a more significant effect on artifact suppression.

[0256] Table 2 presents the quantitative evaluation of imaging results and time consumption for all methods. Due to effective suppression of artifacts and sidelobes, the quantitative evaluation results of GSBL and the proposed method are superior to RMA, IJT, and Cauchy methods, especially showing a significant advantage in eigenvalue decomposition. Furthermore, IJT and GSBL require the use of power methods to calculate the Lipschitz constant L, thus incurring high time consumption. It is worth noting that the computational efficiency of the power method introduced in this application is significantly lower than that of direct eigenvalue decomposition. In contrast, the proposed LSBL method effectively combines the advantages of CS and MF methods, achieving superior imaging results while maintaining low time consumption and memory requirements.

[0257] This application primarily uses the publicly available dataset [3D RIED: A High-Resolution 3D Millimeter-Wave Radar Dataset Dedicated to Imaging and Evaluation] to verify the performance of the proposed method and the comparative CS method in imaging. This dataset uses FMCW signals transmitted in the W band, with a frequency range of 77-81 GHz. The sensor's scanning trajectory forms an array plane with a horizontal sampling interval of 1 mm and a vertical sampling interval of 2 mm, and N... v =407, N h =200 and N k =256. The synthetic aperture size is 0.4m × 0.4m, resulting in 81,400 equivalent array elements on the array plane. The original echo data size is 407 × 200 × 256. Based on real measurement data, this application considers scenarios with a single wrench target. Furthermore, the imaging grid for 2-D imaging slice-by-slice using the IJT, Cauchy, and GSBL methods is divided into 150 × 150 grids.

[0258] In the single-target imaging experiment, an imaging scenario of a wrench was performed, where the vertical distance between the target and the array plane was approximately 0.57 m. During 3D imaging, the number of slices selected in the distance-up direction was set to 6. Figure 5 The imaging results of all methods are shown in Table 3, with the imaging results of the third slice also presented as 2D results. It was observed that the RMA method's imaging results all exhibited a considerable degree of sidelobes and artifacts. At an SR of 40%, the IJT method's imaging results still showed some sidelobes and artifacts; however, its imaging quality improved accordingly with increasing SR. Although the Cauchy method based on non-convex regularization suppressed artifacts to some extent using the Cauchy penalty term, it also suppressed the target itself, leading to some distortion in the imaging results. In contrast, the proposed method based on SBL and the GSBL method showed good suppression of artifacts and clutter, presenting a cleaner background and better reconstructing the 3D image of the target. Table 3 lists the quantitative numerical evaluation and time consumption of the imaging results of all methods. Due to the presence of artifacts and the suppression of the target, the numerical evaluation of the IJT and Cauchy methods was worse than that of the GSBL and the proposed method. The proposed method, equivalent to an adaptive reweighted l1 optimization problem, played a better role in estimating the imaging target, thus achieving higher IC and lower IE. Comparing time consumption reveals that CS-based imaging methods require constructing the observation matrix slice by slice during 3D imaging, resulting in high time consumption and memory requirements. In contrast, the proposed LSBL method exhibits a significant advantage in computational efficiency.

[0259] Table 3 shows the quantitative evaluation and time consumption of imaging results for the wrench target using all methods.

[0260]

[0261] The beneficial effects of this invention are as follows:

[0262] This invention proposes a novel sparse imaging method, LSBL, for 3D MMW imaging. This method constructs a hierarchical Bayesian framework with an adaptive Laplacian prior, avoiding additional regularization parameter inputs, and employs a variational inference strategy for optimization. To further improve computational efficiency and reduce memory requirements, an imaging method based on MF theory is introduced, constructing an imaging operator to replace the vectorization operations in the CS model. This strategy effectively combines the advantages of CS and MF methods in imaging, resulting in a significant improvement in computational efficiency. Finally, extensive simulation and experimental data verification demonstrate the superior performance of this invention in achieving efficient 3D MMW image reconstruction from undersampled echoes.

[0263] As used herein, the term "preferred" is meant as an example, illustration, or illustration. Any aspect or design described herein as "preferred" need not be construed as being more advantageous than other aspects or designs. Rather, the use of the term "preferred" is intended to present the concept in a specific manner. As used in this application, the term "or" is intended to mean an inclusive "or" rather than an exclusionary "or." That is, unless otherwise specified or clear from the context, "X uses A or B" naturally includes either of the permutations. That is, if X uses A; X uses B; or X uses both A and B, then "X uses A or B" is satisfied in any of the foregoing examples.

[0264] Furthermore, although this disclosure has been shown and described with respect to one or more implementations, equivalent variations and modifications will occur to those skilled in the art based on a reading and understanding of this specification and the accompanying drawings. This disclosure includes all such modifications and variations and is limited only by the scope of the appended claims. In particular, with respect to the various functions performed by the aforementioned components (e.g., elements, etc.), the terminology used to describe such components is intended to correspond to any component (unless otherwise indicated) that performs the specified function of said component (e.g., is functionally equivalent to it), even if structurally not equivalent to the disclosed structure performing the functions in the exemplary implementations of this disclosure shown herein. Moreover, although specific features of this disclosure have been disclosed with respect to only one of several implementations, such features may be combined with one or more features of other implementations that may be desirable and advantageous for a given or particular application. Furthermore, with regard to the use of the terms “comprising,” “having,” “containing,” or variations thereof in the Detailed Description or claims, such terms are intended to be included in a manner similar to the term “including.”

[0265] The functional units in this invention embodiment can be integrated into a processing module, or each unit can exist physically separately, or multiple units can be integrated into a module. The integrated module can be implemented in hardware or as a software functional module. If the integrated module is implemented as a software functional module and sold or used as an independent product, it can also be stored in a computer-readable storage medium. The storage medium mentioned above can be a read-only memory, a disk, or an optical disk, etc. The aforementioned devices or systems can execute the storage methods in the corresponding method embodiments.

[0266] In summary, the above embodiments are one implementation of the present invention, but the implementation of the present invention is not limited to the embodiments described above. Any changes, modifications, substitutions, combinations, or simplifications made that deviate from the spirit and principle of the present invention should be considered equivalent substitutions and are included within the protection scope of the present invention.

Claims

1. A three-dimensional millimeter-wave imaging method based on Laplace prior sparse Bayesian learning, characterized in that, Includes the following steps: Acquire the echo data s after range focusing, and construct the imaging operator. and conjugate operator in, Where z represents the distance from the array plane to the target along the range direction, and k z For the corresponding wavenumber; IFFT 2D Represents the two-dimensional inverse fast Fourier transform, FFT 2D Θ represents the two-dimensional Fast Fourier Transform; e is the natural constant, j is the imaginary unit, ⊙ represents the Hadamard product; Θ represents the three-dimensional imaging result; Within the Bayesian modeling framework, all unknown variables are treated as random variables and assigned a prior distribution. A three-level hierarchical model is adopted, where the first level models the elements in the imaging result Θ according to a Gaussian distribution with a mean of 0 and a variance of α; the second level models the elements of α using a gamma distribution with a scale parameter of β; and the third level assigns a gamma prior to β with shape and scale parameters a and b, respectively, to model the imaging result Θ, resulting in a Laplace prior form. If the termination condition is not met, perform the following steps: Calculate ∑ Θ =(LγI+diag(α) -1 )) -1 α is the variance, L is the Lipschitz constant, γ is the noise variance, I is the identity matrix, and ∑ Θ It is the covariance matrix of the posterior distribution; Calculate the estimated value of the three-dimensional imaging result Θ δ is an intermediate variable in the iteration, and its estimated value is μ. Θ Initialize to And it changes with iteration; The estimated value μ of the three-dimensional imaging result Θ Θ Assign the value to Θ; calculate α follows a gamma distribution, where β is the scale parameter; calculate β follows a gamma prior distribution, where a and b are the shape and scale parameters of the β gamma prior distribution; calculate M represents the amount of echo data, c and d are hyperparameters of γ, and g(Θ,δ) is a relaxation function with respect to Θ and δ. The estimated value μ of the three-dimensional imaging result Θ Θ Assign the value to δ; Until the termination condition is met; Output three-dimensional imaging result Θ=μ Θ .

2. The three-dimensional millimeter-wave imaging method based on Laplacian prior sparse Bayesian learning according to claim 1, characterized in that, Within the Bayesian modeling framework, all unknown variables are treated as random variables and assigned a prior distribution. A three-level hierarchical model is adopted, where the first level models the elements in the imaging result Θ according to a Gaussian distribution with a mean of 0 and a variance of α; the second level models the elements of α using a gamma distribution with a scale parameter of β; and the third level assigns a gamma prior to β with shape and scale parameters a and b, respectively, to model the imaging result Θ, resulting in a Laplace prior form. Specifically, this includes: imaging models considering distance focusing: s=ΦΘ+w Where Φ is the observation matrix constructed based on the imaging geometry, and w is noise; First order: The elements of Θ follow an independent Gaussian distribution: Where α = [α1, α2, ..., α N ] T ;α n It is the nth element of α, N is the number of elements in α, and p() represents the conditional probability density function; Second order: For the nth element α n Described using a gamma distribution: Where β = [β1, β2, ..., β 2N ] T ,β n It is the nth element of β, and Third order: β is endowed with a gamma prior, that is The synthesized Laplace prior form is as follows: Where, Θ n Let Θ be the nth element of vector Θ; The noise w follows a complex mean of zero and a variance of γ. -1 The Gaussian distribution of is given by: In practical radar systems, the noise variance γ -1 It is unknown, and the probability density function of γ is described using a gamma distribution: p(γ)=Γ(γ∣c,d)=Γ(c) -1 d c c c-1 e -dγ Based on the imaging model and the probability density function of the noise, the likelihood is described as follows: M represents the amount of echo data; Maximum a posteriori estimation corresponds to the adaptive reweighted l1 norm problem: Solve the following optimization problem to obtain the maximum a posteriori estimate: in, Assign different weights related to the data to the N distinct elements in Θ.

3. The three-dimensional millimeter-wave imaging method based on Laplacian prior sparse Bayesian learning according to claim 2, characterized in that, Estimate the regularization parameter λ n It is estimated through a variational inference iterative process.

4. The three-dimensional millimeter-wave imaging method based on Laplacian prior sparse Bayesian learning according to claim 3, characterized in that, The variational inference iterative process includes: use In a hierarchical model, the goal of Bayesian inference is to find the posterior distribution p(θ|s) of the hidden variable θ in the echo data s. p(θ|s) is approximately: p(θ∣s)≈q Θ (Θ)q α (a)q β (b)q γ (c), Where q Θ (Θ), q α (α), q β (β) and q γ (γ) denote the posterior distribution functions in factorized form for variables Θ, α, β, and γ, respectively; Maximizing a lower bound F(q) of the log-likelihood ln p(s) is achieved by alternating the operations on each variable in θ, which leads to: Where const represents a constant, θ i Denotes a subset of θ, θ\θ i Represents θ with respect to θ i The complement; For f(Θ) = ||s - ΦΘ|| 2 And for any Θ, δ, we have: f(Θ)≤g(Θ,δ)=||s-Φδ|| 2 +2(I-d) T F T (Φδ-s)+L||Θ-δ|| 2 and Where g(Θ,δ) is a relaxation function constructed with respect to Θ and δ given s, Φ and L; Let g(Θ,δ) represent the likelihood function constructed based on g(Θ,δ), and Φ represent the observation matrix constructed based on the imaging geometry. p(s|Θ,γ) represents the likelihood description of s if and only if Θ=δ; Each hidden variable is calculated using an alternating optimization approach. The posterior distribution q Θ (Θ)q α (α)q β (β)q γ The approximate value of (γ).

5. The three-dimensional millimeter-wave imaging method based on Laplace prior sparse Bayesian learning according to claim 4, characterized in that, The method of calculating each hidden variable through alternating optimization is described. The posterior distribution q Θ (Θ)q α (α)q β (β)q γ Approximate values ​​for (γ) include: q Θ Update of (Θ): Ignore terms irrelevant to Θ, approximate posterior distribution q Θ (Θ) is calculated in the following way: Where E(.) is the expectation operator, q Θ (Θ) follows a Gaussian distribution, and its mean and covariance matrices are as follows: m Θ =γ∑ Θ (Lδ+Φ T s-F T (F) ∑ Θ =(LγI+diag(α -1 )) -1 The update formula for Θ is: Θ=μ Θ q α (α) update: approximate posterior q α (α) yielded the following results: in Represents ∑ Θ The nth diagonal element is calculated using the block coordinate descent method: Equaling this expression to zero results in: Update q β (β): The approximate posterior of β is: This indicates that β follows the parameters a+1 and Gamma distribution: Therefore, β n The update formula is: q γ Update (γ): for q γ (γ) is subjected to variational optimization, resulting in: Therefore, γ follows the parameter and Gamma distribution: Update of δ: Find an estimate of δ through the following optimization: Taking the derivative of g(Θ,δ) with respect to δ and setting it to 0, we get: E((LI-Φ T Φ)(δ-Θ))=0 And: δ=Θ The imaging result Θ of the range plane z = z0 is obtained through the above iterative optimization process.

6. The three-dimensional millimeter-wave imaging method based on Laplacian prior sparse Bayesian learning according to claim 5, characterized in that, Based on the imaging results Θ of the distance plane z = z0, then all reconstructed N z Stacking 2D slices along the range direction yields a 3D imaging result; then, on N... z When performing 2D imaging of a range slice, this is achieved by replacing the exact observation function in the CS-based sparse imaging framework with approximate observations obtained from the MF-based imaging process: After distance focusing, s represents the corresponding echo data, and the reconstructed 3D image is denoted as Θ. Zero-padding is applied to both s and Θ, and then a 2D Fourier transform of the appropriate size is performed to obtain: Where, k z It corresponds to the wavenumber component in the distance direction, i.e., the z-axis direction, k z The values ​​of z and z vary with different imaging slices in the distance direction. The imaging operator M(·) and its inverse operator corresponding to the above equation are used. It approximates the large-scale matrix-vector multiplication operation, thus directly performing three-dimensional imaging on the range-focused echo data s.

7. The three-dimensional millimeter-wave imaging method based on Laplace prior sparse Bayesian learning according to claim 6, characterized in that, The termination condition is that the maximum number of iterations is reached or the condition is met on the t-th iteration. ε = 10 -1 This represents a minimum value, with a maximum number of iterations of 20.

Images