A multi-dwelling observation 2D super-resolution ISAR imaging method based on fast IWF

Abstract

The application relates to a kind of multi-resident observation 2D super-resolution ISAR imaging methods based on fast IWF, comprising S1 modeling 2D super-resolution ISAR imaging of multi-resident observation signal;S2 initialize the accuracy of signal in fast IWF method, and set the maximum iteration number of fast IWF method;S3 calculate the LC decomposition factor of the inverse matrix G ‑1 of auxiliary observation covariance matrix;S4 estimate target rotation speed, for compensating distance space variable phase error, and obtain error compensation after observation signal;S5 based on the LC decomposition factor of G ‑1 And distance space variable phase error compensation after observation signal, signal estimation value is calculated using 2D-FFT, and the accuracy of signal is updated;S6 repeat steps S3-S5 to carry out cyclic iteration, stop after reaching convergence, obtain super-resolution ISAR image based on signal estimation value and realize transverse calibration based on target rotation speed estimation value.The application has low computational complexity, strong robustness, high reconstruction accuracy and fast convergence, and can efficiently obtain super-resolution ISAR image with good focusing effect and accurate target size estimation value.

Description

Technical Field

[0001] This invention belongs to the field of computer simulation and method optimization technology, specifically relating to a multi-stationary observation 2D (two-dimensional) super-resolution ISAR (inverse synthetic aperture radar) imaging method based on fast IWF (iterative Wiener filter). Fast IWF, also known as FIWF, is applicable to two-dimensional joint super-resolution ISAR imaging of non-cooperative targets under multi-stationary observation. Background Technology

[0002] In existing technologies, inverse synthetic aperture radar (ISAR) can obtain super-resolution images of non-cooperative targets such as aircraft, ships, missiles, and satellites, providing strong support for subsequent accurate target identification and three-dimensional attitude estimation. With advancements in technology, ISAR observation targets are increasingly becoming smaller, denser, faster, and more maneuverable, placing higher demands on super-resolution ISAR imaging. According to radar imaging principles, the range resolution and azimuth resolution of a target are proportional to the bandwidth of the radar's transmitted signal and the angle of rotation of the target relative to the radar's line of sight during the pulse accumulation time, respectively. However, in practical applications, large bandwidths result in excessively high radar sampling rates and large data volumes, increasing the difficulty of radar hardware implementation. For multi-functional phased array radars, ISAR imaging typically involves coherent processing of multiple stationary observation signals to obtain sufficient azimuth resolution. However, traditional Doppler processing methods based on fast Fourier transform (FFT) produce high grating lobes and sidelobes in the azimuth dimension. Therefore, researching a robust and efficient two-dimensional (2D) joint super-resolution ISAR imaging method for multi-resident observation signals is of great practical significance.

[0003] Currently, given the sparse electromagnetic scattering characteristics of ISAR observation targets, sparse signal reconstruction methods have become the most effective methods for super-resolution ISAR imaging using sparse aperture (SA) signals. Sparse reconstruction methods can be divided into two categories: the first is uncompressed sensing methods, represented by methods such as sequence-CLEAN, relaxation (RELAX), and gapped-data amplitude and phase estimation (GAPES). These methods reconstruct missing data by obtaining model parameters. The second category is compressed sensing methods. Compressed sensing technology breaks the limitations of the Nyquist sampling theorem, enabling super-resolution ISAR imaging with limited echo data and data gaps. Compressed sensing sparse signal reconstruction methods are further divided into three categories: greedy methods, l p Norm regularization methods and Bayesian methods.

[0004] For existing SA-ISAR imaging methods, while uncompressed sensing methods can reconstruct target ISAR images, they are highly susceptible to model errors and noise, and also have high computational complexity. Compressed sensing methods, although capable of producing super-resolution ISAR images, all have certain limitations. Greedy tracking methods have low computational complexity but poor reconstruction accuracy, making them unsuitable for high-precision imaging; regularization methods require manual setting of regularization coefficients, limiting their adaptability; Bayesian methods offer higher reconstruction accuracy and are more robust than the previous two types, but their high computational complexity makes real-time imaging difficult. Some proposed fast sparse Bayesian learning (SBL) methods use approximation methods instead of time-consuming matrix inversion operations, improving computational efficiency, but at the cost of some loss in reconstruction accuracy, resulting in poor imaging performance when effective echo data is limited. Summary of the Invention

[0005] To address the aforementioned technical problems in existing technologies, this invention proposes a multi-stationary observation 2D super-resolution ISAR imaging method based on fast IWF, which specifically includes the following steps:

[0006] S1. Modeling 2D super-resolution ISAR imaging of multi-stationary observation signals;

[0007] S2. Initialize the precision of the signal in the fast IWF method and set the maximum number of iterations for the fast IWF method;

[0008] S3. Based on the accuracy of the reconstructed signal and the noise variance estimated using auxiliary data, the elements within the auxiliary observation covariance matrix are calculated using 2D-FFT. Then, the inverse matrix G of the auxiliary observation covariance matrix is ​​calculated using the 2D Levinson-Durbin method. -1 LC decomposition factor;

[0009] S4. In the image reconstruction iteration, the minimum entropy criterion and the trust region method are used to estimate the target rotation speed to compensate for the range-space-varying phase error and obtain the error-compensated observation signal.

[0010] S5. Based on the inverse matrix G of the observation covariance matrix -1 The LC decomposition factor and the observed signal after range-space-varying phase error compensation are used to calculate the signal estimate using 2D-FFT and update the signal accuracy value.

[0011] S6. Repeat steps S3 to S5 for iterative processing until convergence is achieved. Obtain super-resolution ISAR images based on signal estimation values ​​and achieve lateral calibration based on target rotation speed estimation values.

[0012] Furthermore, S1 specifically includes:

[0013] First, assume the radar transmits a linear frequency modulated (LFM) signal. After LFM demodulation and removal of residual phase terms, the target echo signal received by the radar can be represented in the range frequency domain and azimuth time domain as follows:

[0014]

[0015] Where f∈[-B / 2,B / 2], representing the distance frequency, B represents the bandwidth, and t m Represents slow time, c is the speed of light, f c Let δ be the center carrier frequency of the radar, and assume there are I scattering points on the target. i R is the backscattering coefficient at the i-th scattering point. i (t m R is the instantaneous distance between the i-th scattering point and the radar. i (t m ) is represented as:

[0016] R i (t m )=R0+r(t m )+x i sinθ(t m )+y i cosθ(t m (2)

[0017] Where R0 represents the initial distance, the target's motion relative to the radar is decomposed into translational and rotational components, r(t) m() represents the instantaneous change in slope distance caused by the translational component, x i sinθ(t m )+y i cosθ(t m (x) represents the instantaneous change in slant distance caused by the rotational component. i ,y i ) represents the coordinates of the i-th scattering point, θ(t) m θ(t) represents the target's rotation angle relative to the radar. m )=ωt m ω is the target rotational speed. A second-order Taylor series expansion is performed on the trigonometric functions in formula (2), R... i (t m Represented as:

[0018]

[0019] Assuming translational errors have been compensated and linear distance travel due to rotation has been corrected via Keystone transformation, then... Approximately ignoring the envelope curvature caused by rotation, the echo signal with range-space-varying phase error is expressed in the range frequency domain and azimuth time domain as follows:

[0020]

[0021] Assuming the number of effective pulses within the coherent accumulation time of the sparse aperture signal is L, and the number of range frequency points is N, the discretized model of 2D-SA super-resolution ISAR imaging with range-varying phase errors is obtained as follows:

[0022]

[0023] in, and These represent the echo signal, the super-resolution ISAR image, and the complex Gaussian white noise matrix, respectively. Represents a complete Fourier dictionary matrix. and K1 and K2 represent the range and azimuth dictionary matrices, respectively, both being overcomplete Fourier dictionaries. K1 and K2 represent the number of range and Doppler elements after super-resolution, respectively. If the full-aperture signal contains M pulses, then the super-resolution factors for the range and azimuth dimensions are K1 / N and K2 / M, respectively. ⊙ represents the Hadamard product. The element in E represents the range-varying phase error caused by the target's rotation. in, It is the coordinate of the nth distance unit, (·) H The conjugate transpose operation of a matrix is ​​represented by (·). T This represents the transpose operation of a matrix.

[0024] Vectorizing equation (6) yields the 2D sparse signal reconstruction model for SA signals:

[0025]

[0026] in, Represents a 2D Fourier dictionary matrix. Representing Kronecker, This represents the vectorized form of the SA echo signal matrix. Vec(·) represents the vectorized form of the noise matrix, indicating that the matrix is ​​vectorized by columns. This represents the vectorized form of a super-resolution ISAR image. Diag(·) represents a diagonal matrix composed of elements from E. Diag(·) indicates constructing a diagonal matrix using vectors as diagonal elements. It is by F R The matrix formed by the diagonal elements has the following form:

[0027]

[0028] Furthermore, S2 specifically includes:

[0029] Let the vector of signal precision Each element in the array has a value of 1, and the maximum number of iterations for the fast IWF method is set to J. max Where the superscript (0) represents the initial value, and K1 and K2 are the dimensions of the super-resolution ISAR image. It represents the set of complex numbers.

[0030] Furthermore, S3 specifically includes:

[0031] For j = 0, 1, 2, ..., J max Based on the accuracy γ of the signal in the j-th iteration (j) The noise variance estimate δ obtained from the auxiliary data 2 According to Q=(HΛH) H +δ 2 I)=S T GS uses 2D-FFT to quickly compute the elements in the auxiliary observation covariance matrix G, and then uses the 2D Levinson-Durbin method to compute G. -1 The LC decomposition factor is given by: where Q is the observation covariance matrix, H is the 2D dictionary matrix, Λ is the signal autocorrelation matrix, S is the permutation matrix consisting of 0 and 1 elements, I is the identity matrix, and the superscript (j) represents the j-th iteration. H The conjugate transpose operation of a matrix is ​​represented by (·). T The transpose operation represents a matrix, (·). -1This represents the matrix inversion operation.

[0032] Furthermore, calculate the auxiliary observation covariance matrices G and G in S3. -1 The LC decomposition factor is determined by the following steps:

[0033] The observation model for Wiener filtering is:

[0034] y = z + n (9)

[0035] Among them, y=[y(0) y(1) … y(M-1)] T These are observation data, where y(m) represents an element in y, m is the index of an element in y, and z = [z(0) z(1) … z(M-1)]. T Let z(m) be the desired signal, m be the index of an element in z, and n be the observation noise. The optimal weights and the estimate of the desired signal for Wiener filtering are as follows:

[0036]

[0037] Among them, R yy R is the autocorrelation matrix of the filter input signal. yy =E[yy H ] = E[(z+n)(z+n) H ] = R zz +R nn , where R zz and R nn These are the autocorrelation matrix of the desired signal and the autocorrelation matrix of the noise, respectively. yz It is the cross-correlation matrix between the input signal and the desired signal, and E[·] represents the expectation operation;

[0038] Substituting the 2D sparse signal reconstruction model represented by formula (7) into the Wiener filter observation model, assuming that the range-space-varying phase error has been compensated, let the echo after range-space-varying phase error compensation be... in(·) * Representing the conjugate operation of matrices, the echo model, the autocorrelation matrix of the observed signal, the cross-correlation matrix, and the desired signal are respectively in the following forms:

[0039] y g =Hx+n (12)

[0040] R yy =R zz +R nn =HΛH H +δ 2 I (13)

[0041]

[0042]

[0043] Among them, H m Representing the m-th row element of H, we obtain R. nn =δ 2 I, δ 2 This is the noise variance, which is estimated using auxiliary data in radar imaging. Λ=E[xx H ] is the autocorrelation matrix of the reconstructed signal. It is constructed using the signal power values; that is, Λ is a diagonal matrix with elements γ. k =|x k | 2 A diagonal matrix, where x k γ represents the k-th element in x. k Representing signal x k The accuracy;

[0044] Because z(m) = H m x, so the signal estimation formula is:

[0045]

[0046] The observation covariance matrix Q is:

[0047] Q=(HΛH H +δ 2 I) (17)

[0048] Where I represents an identity matrix;

[0049] Using the permutation matrix S consisting of elements of 1 and 0, Q can be written as:

[0050] Q = S T GS (26)

[0051] Where G is the auxiliary observation covariance matrix, which is a matrix with a fourth-order Toeplitz tensor expansion:

[0052]

[0053] Wherein, G inner submatrix The form is:

[0054]

[0055] In this context, the superscript <·> represents an element in Q. Since G has a fourth-order Toeplitz tensor matrix structure, G can be written in two different forms as follows:

[0056]

[0057] in:

[0058]

[0059] I Nq It is a reversed matrix, where the elements on the secondary diagonal are 1s and the other elements are 0s.

[0060]

[0061] Applying the block matrix inverse formula to equations (29) and (30) respectively, the inverse matrix of G is obtained as follows:

[0062]

[0063] Where 0 represents a zero matrix of appropriate dimensions, and

[0064]

[0065] Define a dimension M gp ×M gp Lower triangular Toeplitz matrix and a circular matrix

[0066]

[0067] Using formulas (35) and (36), G -1 The displacement matrix is:

[0068]

[0069] in, and Let represent a lower triangular Toeplitz-block permutation matrix consisting of 1s and 0s, and a cyclic-block permutation matrix consisting of 1s and 0s, respectively. Their forms are:

[0070]

[0071] make:

[0072]

[0073] Among them, t l (l=0,1,…,Nq-1) is the l-th column vector in T, p l It is the l-th column vector in P. yes The l-th column vector;

[0074] based on G -1 Written as:

[0075]

[0076] Among them, T, P and For G -1 The LC decomposition factor was calculated using the 2D Levinson-Durbin method, yielding T, P, and [other parameters]. and Defined as a TB matrix and a CB matrix respectively:

[0077]

[0078] Furthermore, S4 specifically includes:

[0079] Obtain the echo signal after phase error compensation Among them, y g This represents the echo signal after phase error compensation. It is a matrix composed of Fourier dictionary matrices. Represents the distance-space phase error, s g Represents the echo signal, (·) * Represents the conjugate operation of a matrix;

[0080] The target rotation speed is estimated using the Tsallis entropy of the image obtained by the fast IWF method as the cost function:

[0081]

[0082] in, This is an estimate of the rotational speed ω, and T(X) is the Tsallis entropy of the image. T(X) has the following form:

[0083]

[0084] The total energy of the image is p is the exponent of Tsallis's entropy. Represents the elements within signal x.

[0085] The first-order partial derivative of the Tsallis entropy of the image reconstructed by the fast IWF method with respect to the target rotation speed is:

[0086]

[0087] The first-order partial derivative of the distance-space phase error with respect to rotational speed is:

[0088]

[0089] Re{·} represents the operation of taking the real part and solving nonlinear equations. Obtain the estimated rotational speed of the target.

[0090] Furthermore, S5 specifically includes:

[0091] Based on G -1 The LC decomposition factor and the echo after phase error compensation, according to Using 2D-FFT to quickly compute estimates of sparse signals And according to γ k =|x k | 2 Update the precision vector γ of the sparse signal (j+1) Where k = 0, 1, ..., K-1 are the index values ​​of the signal precision vector;

[0092] From formula (16), we know that The calculation consists of three steps: φ = Q -1 y g =S T G -1 Sy g , and Substituting formula (47) into φ, the right side of the equation contains products of TB matrices and vectors and products of CB matrices and vectors. These are converted into linear convolution and circular convolution, respectively. Using 2D-FFT for fast calculation, the φ matrix is ​​transformed into Φ = [φ0φ1…φ]. q-1 ], φ i The dimension is N×M gp Construct a matrix using Φ in Then get Finally, the dot product calculation is used.

[0093] Furthermore, S6 specifically includes:

[0094] based on and Calculate iterative error The iteration convergence is determined based on the convergence threshold η. If Er > η and j ≤ J, the convergence is determined. max If Er < η, or j > J, then proceed to the next iteration. max If the iteration stops, the final reconstruction result is... Furthermore, based on the final estimated target rotational speed, azimuth calibration is performed.

[0095] This invention can be used for the detection, imaging and identification of maneuvering targets in complex electromagnetic environments. It has the advantages of low computational complexity, strong robustness, high reconstruction accuracy and fast convergence. It can efficiently acquire super-resolution ISAR images with good focusing effect and relatively accurate target size estimates. Attached Figure Description

[0096] Figure 1 This is a flowchart of the present invention;

[0097] Figure 2 A diagram of a 2D multi-stationary echo signal model;

[0098] Figure 3 The 2D sparse signal reconstruction results are shown for the A model and the multi-station observation data. Among them, (a) is the A letter scattering point model, (b) is the 2D multi-station observation data, (c) is the reconstruction result using 2D-FFT, (d) is the reconstruction result using OMP, (e) is the reconstruction result using SBL, and (f) is the reconstruction result using IWF / FIWF.

[0099] Figure 4 The graphs show the performance of the methods under different q, MR, and SNR values. (a) represents the computation time of each method under different q values, (b) represents the reconstruction error of each method under different q values, (c) represents the computation time of each method under different MR values, (d) represents the reconstruction error of each method under different MR values, (e) represents the computation time of each method under different SNR values, and (f) represents the computation time and reconstruction error of each method under different SNR values.

[0100] Figure 5 The images show the scattering point aircraft model and its 2D echo signal and RD imaging results, where (a) is the scattering point aircraft model, (b) is the 2D complete echo signal, (c) is the RD imaging result of the complete data, (d) is the multi-stationary echo signal, and (e) is the RD imaging result of the multi-stationary echo signal.

[0101] Figure 6 The images show ISAR imaging results for different types of multi-stationary echo signals, as well as OMP and FIWF. (a)-(c) show the echo signal with a station number of 4, the OMP imaging results of the ideal echo, and the FIWF-ME-TR imaging results. (d)-(f) show the echo signal with a station number of 8, the OMP imaging results of the ideal echo, and the FIWF-ME-TR imaging results.

[0102] Figure 7 The images show the echo signals under different MR values ​​and the ISAR images obtained by OMP and FIWF. Among them, (a)-(c) are the echo signals, the OMP imaging results of ideal echo, and the FIWF-ME-TR imaging results when the azimuth aperture MR is 50%, respectively; and (d)-(f) are the echo signals, the OMP imaging results of ideal echo, and the FIWF-ME-TR imaging results when the azimuth aperture MR is 81%, respectively.

[0103] Figure 8The images show the echo signals and ISAR images obtained by OMP and FIWF under different SNR conditions. (a)-(c) are the echo signals, the OMP imaging results of the ideal echo, and the FIWF-ME-TR imaging results, respectively, when the SNR is 15dB. (d)-(f) are the echo signals, the OMP imaging results of the ideal echo, and the FIWF-ME-TR imaging results, respectively, when the SNR is 5dB.

[0104] Figure 9 The imaging time, image entropy, and REE of the target rotation estimated by ME-TR for OMP and FIWF-ME-TR under different imaging conditions are given. Among them, (a) is the imaging time of OMP and FIWF-ME-TR under different q values; (b) and (c) are the imaging time, image entropy, and REE of the target rotation estimated by ME-TR for OMP and FIWF-ME-TR under different MR and SNR conditions, respectively.

[0105] Figure 10 The images show the optical images, HRRPs, and full-aperture RD imaging results of the Boeing 787. (a) is the optical image of the Boeing 787, (b) is the HRRPs after envelope alignment and initial phase error compensation, and (c) is the full-aperture RD imaging result.

[0106] Figure 11 The images show ISAR images of the Boeing 787 obtained by multi-stationary HRRPs, GP-FIWF, SAMP, and FIWF-ME-TR. (a)-(d) show HRRPs with MR of 70% and q of 2, and the corresponding GP-FIWF, SAMP, and IWF-ME-TR imaging results, respectively. (e)-(h) show HRRPs with MR of 75% and q of 4, and the corresponding GP-FIWF, SAMP, and FIWF-ME-TR imaging results, respectively. Detailed Implementation

[0107] To make the objectives, contents, and advantages of the present invention clearer, the specific embodiments of the present invention will be described in further detail below with reference to the accompanying drawings and examples.

[0108] This invention addresses the ineffectiveness of traditional imaging methods for multi-stationary echo signals used in ISAR imaging by multi-functional radar, and proposes a multi-stationary observation 2D super-resolution ISAR imaging method based on iterative Wiener filter (IWF) to reduce the pressure on radar hardware design.

[0109] Based on the fact that the noise variance of radar observations can be estimated using auxiliary data and that the amplitude of the target scatterer can be regarded as an unknown deterministic variable, this invention proposes another implementation method for the minimum mean square error (MMSE) criterion estimator, namely the IWF method. Compared with SBL, this method has the advantages of low computational cost, fast convergence speed and high reconstruction accuracy.

[0110] Furthermore, this invention utilizes the fact that the covariance matrix in the 2D signal reconstruction iteration of IWF for multi-resident observation signals has an expansion structure of a fourth-order Toeplitz tensor, and proposes a fast IWF method that can achieve high-efficiency and high-precision reconstruction of 2D sparse signals from multi-resident observation data.

[0111] Furthermore, to ensure compatibility between the fast IWF method and the range-space-variable phase error compensation method, this invention designs a novel range-space-variable autofocusing method. In fast IWF-based image reconstruction, it utilizes image minimum entropy and trust region methods to estimate the target rotational speed, enabling joint 2D super-resolution ISAR imaging and azimuth calibration. This invention offers advantages such as low computational complexity, high reconstruction accuracy, small memory footprint, and good noise suppression capabilities.

[0112] This invention utilizes the fact that the noise variance of radar observations can be estimated through auxiliary data, and that the amplitude of the target scatterer can be considered as an unknown deterministic variable, to propose the IWF method. Furthermore, for multi-stationary observation signals, a 2D super-resolution ISAR imaging method based on the fast IWF method is proposed. In the proposed fast IWF method (FIWF for short), the covariance matrix in the 2D signal reconstruction iteration has a fourth-order Toeplitz tensor expansion structure. This invention designs a lower-triangular-Toeplitz-cyclic (LC) decomposition method to avoid time-consuming matrix inversion operations. Based on the LC decomposition factor, the signal can be quickly reconstructed using FFT in each iteration, greatly reducing the computational complexity of the method.

[0113] like Figure 1 As shown, the present invention includes the following steps:

[0114] S1. Modeling 2D super-resolution ISAR imaging of multi-stationary observation signals;

[0115] S2. Initialize the precision of the signal in the fast IWF method and set the maximum number of iterations for the fast IWF method.

[0116] Let the vector of signal precision Each element in the array has a value of 1, and the maximum number of iterations for the FIWF method is set to J.max Wherein, the superscript (0) represents the initial value, and the dimension of the super-resolution ISAR image is K1×K2. Represents the set of complex numbers;

[0117] S3. Based on the accuracy of the reconstructed signal and the noise variance estimated using auxiliary data, the elements within the auxiliary observation covariance matrix are calculated using 2D-FFT, and then the LC decomposition factor of the inverse matrix of the auxiliary observation covariance matrix is ​​calculated using the 2D Levinson-Durbin method.

[0118] For j = 0, 1, 2, ..., J max Based on the accuracy γ of the signal in the j-th iteration (j) The noise variance estimate δ obtained from the auxiliary data 2 According to Q=(HΛH) H +δ 2 I)=S T GS uses 2D-FFT to quickly compute the elements within the auxiliary covariance matrix G, and then uses the 2D Levinson-Durbin method to compute G. -1 The LC decomposition factor. Where Q is called the observation covariance matrix, H is the 2D dictionary matrix, Λ is the signal autocorrelation matrix, S is the permutation matrix consisting of 0 and 1 elements, I represents the identity matrix, and the superscript (j) represents the j-th iteration, (·) H The conjugate transpose operation of a matrix is ​​represented by (·). T The transpose operation represents a matrix, (·). -1 This represents the matrix inversion operation;

[0119] S4. In the image reconstruction iteration, the minimum entropy criterion and trust region method are used to estimate the target rotation speed to compensate for the range-space phase error and obtain the error-compensated observation signal.

[0120] Obtain the echo signal after phase error compensation Among them, y g This represents the echo signal after phase error compensation. It is a matrix composed of Fourier dictionary matrices. Represents the distance-space phase error, s g Represents the echo signal, (·) * Represents the conjugate operation of a matrix;

[0121] S5. Based on the inverse matrix G of the observation covariance matrix -1 The LC decomposition factor and the observed signal after distance-space-varying phase error compensation are used to calculate the signal estimate using 2D-FFT and update the signal accuracy value.

[0122] Based on G -1The LC decomposition factor and the echo after phase error compensation, according to Using 2D-FFT to quickly compute estimates of sparse signals And according to γ k =|x k | 2 Update the precision vector γ of the signal (j+1) Where k = 0, 1, ..., K-1 are the index values ​​of the signal precision vector;

[0123] S6. Repeat steps S3-S5 for iterative processing until convergence is achieved. A super-resolution ISAR image is obtained based on the signal estimation value, and lateral calibration is performed based on the target rotation speed estimation value.

[0124] based on and Calculate iterative error The iteration convergence is determined based on the convergence threshold η. If Er > η and j ≤ J, the convergence is determined. max If Er < η, or j > J, then proceed to the next iteration. max If the iteration stops, the final reconstruction result is... Furthermore, based on the final estimated target rotational speed, azimuth calibration is performed.

[0125] A further technical solution is to model the 2D super-resolution ISAR imaging before initializing the S2 parameters. The specific steps are as follows:

[0126] First, assuming the radar transmits a linear frequency modulated (LFM) signal, after LFM demodulation and removal of residual phase terms, the target echo signal received by the radar can be expressed in the range frequency domain and azimuth time domain as follows:

[0127]

[0128] Where f∈[-B / 2,B / 2] represents the distance frequency, B represents the bandwidth, and t m Represents slow time, c is the speed of light, f c Let δ be the center carrier frequency of the radar. Assume there are I scattering points on the target, and δ... i R is the backscattering coefficient at the i-th scattering point. i (t m Let be the instantaneous distance between the i-th scattering point and the radar, which can be expressed as:

[0129] R i (t m )=R0+r(t m )+x i sinθ(t m )+y i cosθ(t m (5)

[0130] Where R0 represents the initial distance, the target's motion relative to the radar can be decomposed into translational and rotational components, r(t m The first two terms represent the instantaneous slope change caused by the translational component, and the last two terms represent the instantaneous slope change caused by the rotational component. i ,y i ) represents the coordinates of the i-th scattering point, θ(t) m θ(t) represents the target's rotation angle relative to the radar. Since the coherence accumulation angle of ISAR imaging is relatively small (around 5°), the motion of the maneuvering target during the imaging coherence accumulation time can be approximated as uniform rotation, i.e., θ(t) m )=ωt m ω represents the target's rotational speed. To more precisely represent the target's instantaneous distance, a second-order Taylor series expansion of the trigonometric function is performed, R... i (t m ) can be reformulated as

[0131]

[0132] Substituting equation (6) into equation (4), and assuming that the translational error has been compensated and the distance linearity has been corrected through the Keystone transformation, by... After approximation, the echo signal is represented as

[0133]

[0134] From formula (7), we can see that the quadratic term This leads to range curvature and range-space-varying phase errors in the azimuth dimension. Since the rotational speed of the maneuvering target is relatively small, the range curvature is negligible; however, the phase error compensation accuracy is on the order of wavelength, so range-space-varying phase error compensation must be performed before imaging. After eliminating the constant term, the echo signal with range-space-varying phase error can be expressed in the range frequency domain and azimuth time domain as follows:

[0135]

[0136] Assuming the number of effective pulses within the coherent accumulation time of the sparse aperture signal is L, and the number of range frequency points is N, the discretization model for 2D-SA (two-dimensional sparse aperture) super-resolution ISAR imaging with range-varying phase errors can be obtained as follows:

[0137]

[0138] in, and These represent the echo signal, the super-resolution ISAR image, and the complex Gaussian white noise matrix, respectively. Represents a complete Fourier dictionary matrix. and K1 and K2 represent the range and azimuth dictionary matrices, respectively, both being overcomplete Fourier dictionaries. K1 and K2 represent the number of range and Doppler elements after super-resolution, respectively. If the full-aperture signal contains M pulses, then the super-resolution factors for the range and azimuth dimensions are K1 / N and K2 / M, respectively. ⊙ represents the Hadamard product. The element in E represents the range-varying phase error caused by the target's rotation. in, It is the coordinate of the nth distance unit, which is measurable, so E can be calculated by solving for the target rotational speed ω.

[0139] Vectorizing equation (9) yields the 2D sparse signal reconstruction model for SA (sparse aperture) signals:

[0140]

[0141] in, Represents a 2D Fourier dictionary matrix. Representing Kronecker, This represents the vectorized form of the SA echo signal matrix. Vec(·) represents the vectorized form of the noise matrix, indicating that the matrix is ​​vectorized by columns. The vectorized form representing a super-resolution ISAR image can be expressed as follows: Where the internal elements of x are k2 = 0, 1, ..., K2-1. Diag(·) represents a diagonal matrix composed of elements from SA and E. Diag(·) means constructing a diagonal matrix with vectors as diagonal elements. It is by F R The matrix formed by the diagonal elements has the following form:

[0142]

[0143] A further technical solution is to utilize the sparsity of signal x to achieve 2D-SA super-resolution ISAR imaging using an iterative Wiener filtering method, obtaining the observation covariance matrix Q in S3, the auxiliary observation covariance matrix G, and the estimated value of the sparse signal in S5. Right now The specific steps of the iterative Wiener filtering method include:

[0144] The observation model for Wiener filtering is:

[0145] y = z + n (12)

[0146] Among them, y=[y(0) y(1) … y(M-1)] THere, y(m) represents the observed data, and m is the index of an element in y. z = [z(0) z(1) … z(M-1)] T Let z(m) be the desired signal, z(m) represent an element in z, and m be the index of the element in z. n is the observation noise. The essence of Wiener filtering is to use the observation data to linearly estimate the desired signal.

[0147] The optimal weights and the estimation of the desired signal for Wiener filtering are respectively...

[0148]

[0149] Among them, R yy R is the autocorrelation matrix of the filter input signal. Since the noise and the desired signal are independent, we can obtain R. yy =E[yy H ] = E[(z+n)(z+n) H ] = R zz +R nn , where R zz and R nn These are the autocorrelation matrices of the desired signal and the noise, respectively. yz It is the cross-correlation matrix between the input signal and the desired signal. E[·] denotes the expectation operation.

[0150] Substituting the 2D imaging model shown in formula (10) into the Wiener filter above, assuming that the range-space-varying phase error has been compensated, let... The echo model, the autocorrelation matrix of the observed signal, the cross-correlation matrix, and the desired signal are respectively in the following forms:

[0151] y g =Hx+n (15)

[0152] R yy =R zz +R nn =HΛH H +δ 2 I (16)

[0153]

[0154] Among them, H m This represents the m-th row element of H. Since the noise is complex Gaussian white noise, we can obtain R. nn =δ 2 I, δ 2 This refers to the noise variance. The noise variance of the echo signal in radar imaging can be estimated using auxiliary data. Λ=E[xx H] is the autocorrelation matrix of the reconstructed signal. From formula (18), it can be seen that Wiener filtering requires knowledge of the autocorrelation matrix of the reconstructed signal. However, in imaging, the autocorrelation matrix of the signal is unknown. The best solution is iteration, that is, setting initial values ​​and iteratively solving based on known parameters. This is the key to achieving super-resolution imaging using IWF. Since the amplitude of the target scatterer can be considered as an unknown deterministic variable, the autocorrelation matrix of the signal can be constructed using the signal power value. Therefore, Λ is a diagonal element with γ k =|x k | 2 A diagonal matrix, where x k Represents the k-th element in x, γ k Represents x k The accuracy is also similar to that of the Iterative Adaptive Approach (IAA) method.

[0155] Because z(m) = H m x, so the signal estimation formula is:

[0156]

[0157] Where the observation covariance matrix Q is

[0158] Q=(HΛH H +δ 2 I) (20)

[0159] Here, I represents an identity matrix.

[0160] Observing formula (19), it can be found that the signal estimation formula in the IWF method is the same as that in the SBL method. The difference between IWF and SBL is that the calculation methods for noise variance and signal autocorrelation matrix are different. SBL obtains them through the maximum expectation (EM) method, while IWF directly uses the signal power value to calculate the signal autocorrelation matrix and estimates the noise variance through auxiliary data. The more auxiliary data, the closer the noise variance estimate is to the true value. Therefore, the imaging accuracy of IWF is better than that of SBL, and the convergence speed is also faster than that of SBL.

[0161] A further technical solution addresses the time-consuming matrix inversion operation in the IWF method iteration. To reduce the computational complexity, a fast IWF method is designed, leveraging the fact that Q has a fourth-order Toeplitz tensor matrix structure under multi-resident echo signals. This method yields the elements of the auxiliary observation covariance matrix G calculated quickly using FFT in S3 and the elements of G calculated using the 2D Levinson-Durbin (LD) method. -1 LC decomposition factor.

[0162] A further technical solution is that Q has a fourth-order Toeplitz tensor matrix structure under multi-residual echo signals, and the elements within Q can be quickly calculated using FFT. The specific steps include:

[0163] First, assume the geometry of the 2D multi-residual echo signal is as follows: Figure 2 As shown. Assume the number of stays is q, and one stay contains M... gp A pulse, with M pulses between two dwell times. mp If there are pulses, then we can obtain M. sp =M gp +M mp L = qM gp The missing rate (MR) is defined as MR = M mp / M sp .

[0164] For multi-stationary echo signals, when reconstructing a 2D signal using IWF, the azimuth dimension dictionary matrix H... A And distance dimension dictionary matrix H R These are a partial Fourier dictionary matrix and a complete Fourier dictionary matrix, H. A The form is:

[0165]

[0166] Where H A(i) H represents the dictionary matrix corresponding to the i-th segment of resident data. R and H A(i) The forms of the elements in the middle column are as follows:

[0167]

[0168] in, k1=0,1,…,K1-1, k2=0,1,…,K2-1.

[0169] Let the observation covariance matrix Q = R + δ 2 Let I and R be the signal covariance matrices. Substituting the dictionary H, we obtain R as a matrix with a fourth-order Toeplitz tensor expansion:

[0170]

[0171] Wherein, the submatrix of R It is a third-order Toeplitz tensor expansion, namely a Toeplitz-block-Toeplitz matrix, with the following form:

[0172]

[0173] in, inner submatrix It is a Toeplitz matrix, with the following form:

[0174]

[0175]

[0176] here Represents H A(i+1) The row vector of the (j+1)th row, The element in the middle can be represented as

[0177]

[0178] Obviously, The elements in the middle can be quickly calculated using 2D-FFT.

[0179] Based on R, Q also has a fourth-order Toeplitz tensor matrix structure, and the structure is the same as that of R. Therefore, the elements in Q can be quickly calculated using 2D-FFT.

[0180] A further technical solution involves constructing an auxiliary observation covariance matrix G using a permutation matrix. The specific steps include:

[0181] Using the permutation matrix S consisting of elements of 1 and 0, Q can be written as:

[0182] Q = S T GS (29)

[0183] Here, G is called the auxiliary observation covariance matrix, which is also a matrix with a fourth-order Toeplitz tensor expansion:

[0184]

[0185] Wherein, the inner submatrix of G The form is:

[0186]

[0187] Here, the superscript <·> represents an element in Q, for example represent Middle elements Clearly, the elements in G can also be calculated quickly using 2D-FFT.

[0188] A further technical solution is to utilize the fourth-order Toeplitz tensor matrix structure of G, and then use LC decomposition to compute G. -1 , to obtain G -1 The LC decomposition factor. Calculate G. -1 and G -1The steps for determining the LC decomposition factor include:

[0189] Since G has a fourth-order Toeplitz tensor matrix structure, G can be written in two different forms:

[0190]

[0191] in:

[0192]

[0193] It is a reversed matrix, where the elements on the secondary diagonal are 1s and the other elements are 0s.

[0194]

[0195] Applying the block matrix inverse formula to equations (32) and (33) respectively, the inverse matrix of G can be obtained as follows:

[0196]

[0197] in:

[0198]

[0199]

[0200] Define a dimension M gp ×M gp Lower triangular Toeplitz matrix and a circular matrix

[0201]

[0202] Using formulas (38) and (39), G -1 The displacement matrix is:

[0203]

[0204] in, and Let represent a lower triangular Toeplitz-block (TB) permutation matrix consisting of 1s and 0s, and a cyclic-block (CB) permutation matrix consisting of 1s and 0s, respectively. Their forms are:

[0205]

[0206] make

[0207]

[0208] Among them, t l (l=0,1,…,Nq-1) is the l-th column vector in T, which can be written as t i Represents t l The (i+1)th subvector in the array, with dimension Nq, can be written as t i,j Represents t i The (j+1)th subvector in the array has dimension N and is of the form t. i,j =[t0t1 … t N-1 ] T p l (l=0,1,…,Nq-1) is the l-th column vector in P. yes The l-th column vector in the vector has the form t l Same.

[0209] based on G -1 It can be written as:

[0210]

[0211] Among them, T, P and Called G -1 The LC decomposition factor can be calculated using the 2D Levinson-Durbin (LD) method. and Defined as a TB matrix and a CB matrix respectively:

[0212]

[0213]

[0214] From formula (50), we can see that G -1 It can be obtained through LC decomposition, and G -1 The right side of the equation is the product of the TB matrix and the CB matrix; therefore, equation (50) is called G. -1 The LC decomposition formula, and the calculation of G based on the G structure. -1 The process is called LC decomposition. G -1 The displacement rank is 2Nq. Using Q... -1 =S T G -1 From S, we know that Q -1 The displacement rank is also Nq.

[0215] A further technical solution involves estimating the target rotational speed during image reconstruction iterations based on a coarse estimate of the target rotational speed, thereby compensating for the spatially varying phase error and obtaining the phase error-compensated echo signal in S4. The specific steps include:

[0216] As shown in the super-resolution imaging model (10), the spatially varying phase error caused by the target rotation component will lead to image defocusing. Observing the form of this phase error, since the position of the range unit is measurable, the target rotation speed can be obtained, thus achieving accurate compensation for the spatially varying phase error. In imaging, image entropy is the most commonly used indicator to measure the imaging focusing effect. Based on the minimum entropy criterion, this paper transforms the target rotation speed estimation into a nonlinear optimization problem as shown in formula (50).

[0217] The target rotation speed is estimated using the Tsallis entropy of the image obtained by fast IWF as the cost function.

[0218]

[0219] in, This is an estimate of the rotational speed ω, and T(X) is the Tsallis entropy of the image. T(X) has the following form:

[0220]

[0221] The total energy of the image is p is the exponent of Tsallis entropy, which is usually 1.6. Represents the elements within signal x.

[0222] Using formula (19), the first-order partial derivative of the Tsallis entropy of the image reconstructed by IWF with respect to the target rotational speed is:

[0223]

[0224] The first-order partial derivative of the distance-space phase error with respect to rotational speed is:

[0225]

[0226] Re{·} represents the operation of extracting the real part.

[0227] Solving nonlinear equations using the trust region algorithm The solution to this nonlinear equation is the estimated rotational speed of the target. The parameter estimation method based on the minimum entropy criterion and the trust region method is abbreviated as ME-TR.

[0228] A further technical solution is to utilize G -1 The LC decomposition factor and the echo signal y after distance space-varying phase error compensation gThis yields the estimated value of the sparse signal obtained in S5 using FFT for fast computation. Fast calculation The steps include:

[0229] From formula (19), we can see that The calculation can be divided into three steps: φ = Q -1 y g =S T G -1 Sy g , and Substituting formula (50) into φ, the right side of the equation consists of products of TB matrices and vectors, and products of CB matrices and vectors. These can be converted into linear convolution and circular convolution, respectively, so 2D-FFT can be used for fast calculation. The φ matrix is ​​transformed into Φ = [φ0 φ1 … φ q-1 ], φ i The dimension is N×M gp Construct a matrix using Φ. in Then we can obtain Finally, the dot product calculation is used.

[0230] The FIWF method utilizes a fourth-order Toeplitz tensor expansion structure and employs LC decomposition to obtain the inverse matrix. Furthermore, aside from solving for the LC decomposition factors, almost all operations can be rapidly computed using FFT. Therefore, FIWF can achieve efficient and high-precision 2D super-resolution ISAR imaging. In the FIWF method, the displacement rank of the inverse matrix of the auxiliary observation covariance matrix is ​​2Nq. A larger displacement rank increases the computational complexity of the method; therefore, FIWF is more efficient when processing observation data with fewer dwell times.

[0231] A further technical solution is to perform orientation calibration based on the final estimated target rotation speed, which provides strong support for subsequent target feature extraction and identification.

[0232] To verify the effectiveness of the multi-stationary observation super-resolution ISAR imaging method based on the fast IWF method proposed in this invention, the following simulation and field measurement comparison experiments were set up. Since there are few 2D super-resolution imaging methods proposed so far, only OMP and SBL are used as comparison methods.

[0233] Example 1:

[0234] In this example, based on such Figure 3(a) shows the scattering point model of letter A, comparing the reconstruction performance of OMP, SBL, and the IWF / FIWF method proposed in this invention. The signal-to-noise ratio of the 2D observed signal is set to 10 dB, with a range frequency of 32 and a pulse number of 64. To compare the sparse signal reconstruction accuracy of the methods, the normalized root mean square error (nRMSE) is defined.

[0235]

[0236] in, x and y represent the reconstructed signal and the real signal, respectively.

[0237] By periodically missing 50% of the complete observation pulse cycle, a multi-resident observation dataset is constructed, such as... Figure 3 As shown in (b). Figure 3 Images (c)-(f) are the reconstructed images using 2D-FFT, OMP, SBL, and IWF / FIWF, respectively. The computation time and reconstruction error of each method are also marked on the images. Figure 3 As can be seen, under multi-resident observation data, the 2D-FFT reconstruction results show obvious grating lobes and ghosting. OMP has lower reconstruction accuracy, performs poorly in areas with dense target scattering points, and contains some false points and missing points. However, SBL and IWF have particularly good reconstruction results, producing clear target images. Compared with SBL, IWF has higher reconstruction accuracy, and FIWF has a computational efficiency two orders of magnitude higher than IWF.

[0238] Example 2:

[0239] In this example, to quantitatively analyze the performance of OMP, SBL, and the IWF / FIWF method proposed in this invention, 100 Monte Carlo experiments were conducted using 2D multi-dwelling observation data of the A-letter scattering point model. The performance curves of the methods under different q, MR, and SNR values ​​were obtained, as shown below. Figure 4 As shown. Figure 4 (a) and Figure 4 (b) shows the reconstruction time and reconstruction error of the method under different types of multi-resident data. In this experiment, the number of range frequencies and pulses in the 2D observations were 32 and 64, respectively, with 50% missing in the azimuth pulse period phase, and the signal-to-noise ratio was 10 dB. From Figure 4 (a) and Figure 4(b) It can be seen that OMP has high computational efficiency but poor reconstruction accuracy. SBL has high reconstruction accuracy but long computation time. IWF has higher computational efficiency than SBL and higher reconstruction accuracy than SBL. FIWF has a much shorter computation time than IWF, but the same reconstruction accuracy. The fewer the number of dwells in the observation data, the higher the computational efficiency of FIWF. This is because the displacement rank of the auxiliary observation covariance inverse matrix in the FIWF method is proportional to the number of dwells. The larger the displacement rank, the higher the computational complexity. Therefore, the fewer the number of dwells in the 2D echo signal, the higher the computational efficiency of FIWF. Figure 4 (c) and Figure 4 (d) shows the computation time and reconstruction error for azimuth aperture MR methods that are not simultaneous, with 32 range frequency points and 64 pulses in 2D observations, and a signal-to-noise ratio of 10 dB. Figure 4 (c) and Figure 4 (d) It can be seen that the more missing pulses, the shorter the computation time, and the larger the reconstruction error. When the data is relatively small, the SBL reconstruction error is very small, and when the data is relatively large, the reconstruction error is still acceptable. Since the noise variance in IWF is closer to the true value, the reconstruction error of IWF is smaller than that of SBL. Due to the denser scattering points in the model, the OMP reconstruction error is large. The computation time and reconstruction error of the methods under different SNRs are as follows: Figure 4 (e) and Figure 4 As shown in (f). MR is 50%, and the distance sampling points and pulse number for 2D observations are 32 and 64, respectively. From Figure 4 (e) and Figure 4 (f) It can be seen that the larger the SNR, the shorter the reconstruction time and the smaller the reconstruction error. The reconstruction error of IWF is smaller than that of SBL. However, when the signal-to-noise ratio is 0dB, the reconstruction errors of SBL and IWF / FIWF are relatively large. This is because the 2D observation data in the experiment contains noise. The more sampling points there are, the more energy accumulates in the reconstructed signal. However, the number of 2D sampling points in the experiment is relatively small, so the reconstruction error is larger when the signal-to-noise ratio is low.

[0240] Example 3:

[0241] In this example, for example Figure 5 (a) shows the scattering point aircraft model, and super-resolution imaging is performed to verify the effectiveness of the fast IWF method and ME-TR. Complex Gaussian white noise is added to the radar echo signal. The simulation experimental system parameters are shown in Table 1. To prevent the optimization method from getting trapped in local optima, it is generally necessary to first use linear search to obtain a coarse estimate of the target rotational speed as the initial value for the next step of accurate target rotational speed estimation. Here, we first assume that the initial value of the target rotational speed is 0.088 rad / s. To analyze the accuracy of ME-TR in estimating the target rotational speed, its relative estimation error (REE) is defined as follows:

[0242]

[0243] in, ω and ω represent the estimated and actual values ​​of the target rotational speed, respectively.

[0244] In this example, Figure 5 (b) and Figure 5 (c) The complete echo signal and RD imaging results with range spatially varying phase error at a signal-to-noise ratio of 10 dB are given for comparison with the SA signal imaging results. Figure 5 (d) and Figure 5 (e) presents the multi-stationary echo signal with 75% pulse loss and the RD imaging results. As can be seen from the figure, the RD image is significantly out of focus due to the lack of compensation for range-varying phase errors, indicating that RD has failed under SA. Next, this example simulates various ISAR imaging environments to verify the effectiveness of OMP and IWF / FIWF under different types of multi-stationary echo data, varying degrees of azimuth aperture loss, and different signal-to-noise ratios.

[0245] Table 1 Radar and Target Motion Parameters

[0246]

[0247]

[0248] Figure 6 Different types of multi-stationary echo signals and super-resolution ISAR images are shown. SNR and MR are set to 10 dB and 75%, respectively. Figure 6 (a)-(c) are the OMP imaging results of the echo signal, the ideal echo signal (the ideal echo signal refers to the echo signal without distance spatial variation error) and the FIWF-ME-TR imaging results when the number of residences is 4, respectively. (The parameter estimation method based on the minimum entropy and trust region method is called ME-TR, and the parameter estimation method using the FIWF method and the method based on the minimum entropy and trust region method is simply referred to as FIWF-ME-TR.) Figure 6 (d)-(f) show the OMP imaging results of the echo signal with a dwell number of 8, the ideal echo signal, and the FIWF-ME-TR imaging results, respectively. Figure 6 It can be seen that OMP suffers from missing and spurious points for incomplete ideal echo signals. However, the fast IWF method can obtain well-focused, high-quality ISAR images, and the target rotation speed estimated using ME-TR is very close to the true value. Since q has no effect on image entropy and target rotation speed estimation, Figure 9(a) Only the computation time of OMP and FIWF-ME-TR at different q values ​​is shown. The analysis of the computation time shows that OMP has high computational efficiency, while FIWF has higher computational efficiency for echo signals with fewer dwell times, which is consistent with the theoretical analysis of the method in Section 3.

[0249] Imaging results under different MR conditions, such as Figure 7 As shown. Combined with Figure 6 (a)-(d) and Figure 7 It can be seen that the more data is missing, the worse the reconstruction effect of the method. When 80% of the data is missing, the basic outline of the target can no longer be seen in the OMP imaging results, while only a few scattering points are incorrectly reconstructed in the IWF imaging results. Figure 9 (b) illustrates the impact of MR on the method. The more missing data, the shorter the computation time, the higher the image entropy, and the greater the deviation between the ME-TR estimated target rotation speed and the true value. However, when the amount of missing data is large, the estimation error of ME-TR is acceptable. The image entropy reconstructed by IWF is smaller than that of OMP, indicating that IWF achieves better target image focusing. The estimated target rotation speed in IWF-based image reconstruction is also closer to the true value.

[0250] Figure 8 (a)-(c) show the echo signal, OMP imaging result, and FIWF-ME-TR imaging result when the SNR is 15dB, respectively. Figure 8 Images (d)-(f) show the echo and ISAR image at a signal-to-noise ratio of 5 dB, respectively. Figure 6 (a)-(c) and Figure 8 It can be seen that, under the same conditions, the higher the signal-to-noise ratio (SNR), the better the imaging effect of the method. However, even at a high SNR, the target image reconstructed by OMP still suffers from missing points and false points due to a large amount of missing data. In contrast, IWF can obtain a clear target image, and only a very small number of scattering points cannot be accurately reconstructed when the SNR is low. Figure 9 (c) Detailed representations of imaging time, image entropy, and REE for ME-TR-estimated target rotation speed for various methods at different signal-to-noise ratios. From Figure 9 As shown in (c), the higher the signal-to-noise ratio (SNR), the lower the image entropy, and the closer the target rotation speed estimated by the MT-TR method is to the true value. Since we set the number of FIWF iterations to 40, the computation time is basically the same under different SNRs.

[0251] Example 4:

[0252] This example presents the processed results of measured data from a Boeing 787 aircraft to verify the effectiveness of the proposed method. The main parameters of the radar system are as follows: carrier frequency of 10 GHz, bandwidth of 400 MHz, and PRF of 100 Hz. The full aperture data includes 1024 range sampling points and 256 pulses. The Boeing 787 optical images, high-resolution range profile sequences (HRRPs), and full aperture RD imaging results are shown below. Figure 10 As shown, a one-dimensional super-resolution ISAR imaging method based on Fast IWF, abbreviated as GP-FIWF, is compared with SAMP. To avoid the phase influence of range-varying phase errors, both SAMP and GP-FIWF process the echo after phase error correction of the target rotation speed estimated by the method proposed in this invention. Based on full aperture data, two different types of multi-dwelling observation data were simulated. Figure 11 (a)-(d) show HRRPs with a dwell number of 2 and a pulse loss of 70%, as well as imaging results from GP-FIWF, SAMP, and FIWF-ME-TR. Figure 11 (e)-(h) show the imaging results of HRRPs with a dwell number of 4 and a pulse loss of 75%, obtained using GP-FIWF, SAMP, and FIWF-ME-TR. Figure 11 It can be seen that although the image obtained by SAMP shows the general shape of the target, there are also some false points, and the more data is missing, the more false points there are. By comparing the one-dimensional super-resolution ISAR image obtained by IWF and the 2D super-resolution ISAR image obtained by IWF, it can be found that the target is clearer after 2D super-resolution, especially at the wing of the 787 aircraft marked with a circle.

[0253] The computation time for the above measurement data experiment is shown in Table 2. Combining Table 2 and the above ISAR images of the Boeing-787, it can be concluded that the fast IWF method proposed in this invention can achieve 2D joint super-resolution ISAR imaging with high efficiency and high precision.

[0254] Table 2 Comparison of Image Focusing Performance

[0255]

[0256] To address the challenges of receiving multi-stationary observation signals for imaging in multi-functional radars and to reduce the burden on radar hardware design, this invention proposes a fast IWF method for efficient and high-precision 2D super-resolution ISAR imaging of multi-stationary observation data. The proposed fast IWF method utilizes the fourth-order Toeplitz tensor expansion of the observation covariance matrix in the IWF iteration and designs an LC decomposition method to find the inverse of the auxiliary covariance matrix, avoiding time-consuming matrix inversion. Furthermore, based on the LC decomposition factor, almost all remaining operations in the IWF method can be rapidly calculated using FFT. Therefore, the proposed fast IWF method exhibits advantages such as low computational complexity, high reconstruction accuracy, strong robustness, and fast convergence speed. In addition, in image reconstruction, minimum entropy of the image is used to calculate the target rotation speed, thereby compensating for range-varying phase errors and achieving target azimuth calibration.

[0257] This invention proposes an IWF method and a fast IWF method for multi-resident observation signals for 2D super-resolution ISAR imaging, which has the advantages of low computational complexity, high reconstruction accuracy, strong robustness, and fast convergence speed. In the proposed fast IWF method, this invention designs an LC decomposition to find an inverse matrix with a fourth-order Toeplitz tensor expansion, avoiding time-consuming matrix inversion operations.

[0258] In order to make the proposed fast IWF method compatible with the range-space-variable autofocus method, this invention designs a new autofocus method. In the image reconstruction iteration, the target rotation speed is obtained by using the minimum entropy of the image, thereby compensating for the range-space-variable phase error and realizing the azimuth calibration.

[0259] The above description is only a preferred embodiment of the present invention. It should be noted that for those skilled in the art, several improvements and modifications can be made without departing from the technical principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A multi-station observation 2D super-resolution ISAR imaging method based on fast IWF, characterized in that, Includes the following steps: S1. Modeling 2D super-resolution ISAR imaging of multi-stationary observation signals; S2. Initialize the precision of the signal in the fast IWF method and set the maximum number of iterations for the fast IWF method; S3. Based on the accuracy of the reconstructed signal and the noise variance estimated using auxiliary data, the elements within the auxiliary observation covariance matrix are calculated using 2D-FFT. Then, the inverse matrix of the auxiliary observation covariance matrix is ​​calculated using the 2D Levinson-Durbin method. LC decomposition factor; S4. In the image reconstruction iteration, the minimum entropy criterion and the trust region method are used to estimate the target rotation speed to compensate for the range-space-varying phase error and obtain the error-compensated observation signal. S5. The inverse matrix based on the observation covariance matrix The LC decomposition factor and the observed signal after range-space-varying phase error compensation are used to calculate the signal estimate using 2D-FFT and update the signal accuracy value. S6. Repeat steps S3 to S5 for iterative processing until convergence is achieved. Obtain super-resolution ISAR images based on signal estimation values ​​and achieve lateral calibration based on target rotation speed estimation values.

2. The multi-station observation 2D super-resolution ISAR imaging method based on fast IWF according to claim 1, characterized in that, S1 specifically includes: First, assume the radar transmits a linear frequency modulated (LFM) signal. After LFM demodulation and removal of residual phase terms, the target echo signal received by the radar can be represented in the range frequency domain and azimuth time domain as follows: in, , where represents the distance frequency, and B represents the bandwidth. Represents slow time, where c is the speed of light. Let I be the center carrier frequency of the radar, and assume there are I scattering points on the target. Let be the backscattering coefficient at the i-th scattering point. Let be the instantaneous distance between the i-th scattering point and the radar. Represented as: in, Representing the initial distance, the target's motion relative to the radar is decomposed into translational and rotational components. This represents the instantaneous change in slant distance caused by the translational component. This represents the instantaneous change in slant distance caused by the rotational component. Represents the coordinates of the i-th scattering point. This indicates the angle of rotation of the target relative to the radar. , To determine the target rotational speed, a second-order Taylor series expansion is performed on the trigonometric functions in formula (2). Represented as: Assuming translational errors have been compensated and linear distance travel due to rotation has been corrected via Keystone transformation, then... Approximately ignoring the envelope curvature caused by rotation, the echo signal with range-space-varying phase error is expressed in the range frequency domain and azimuth time domain as follows: (5) Assuming the number of effective pulses within the coherent accumulation time of the sparse aperture signal is L, and the number of range frequency points is N, the discretized model of 2D-SA super-resolution ISAR imaging with range-varying phase errors is obtained as follows: (6) in, , and These represent the echo signal, the super-resolution ISAR image, and the complex Gaussian white noise matrix, respectively. Represents a complete Fourier dictionary matrix. and These represent the distance-dimensional dictionary matrix and the orientation-dimensional dictionary matrix, respectively, both of which are overcomplete Fourier dictionaries. and Let represent the number of range elements and the number of Doppler elements after super-resolution, respectively. If the full aperture signal contains M pulses, then the super-resolution factors in the range and azimuth dimensions are respectively... and , Represents the Hadamard product. This represents the spatial phase error in range caused by the target's rotation. The middle element is ,in, It is the coordinate of the nth distance unit. This represents the conjugate transpose operation of a matrix. This represents the transpose operation of a matrix. Vectorizing equation (6) yields the 2D sparse signal reconstruction model for SA signals: (7) in, Represents a 2D Fourier dictionary matrix. Representing Kronecker, This represents the vectorized form of the SA echo signal matrix. The vectorized form of the noise matrix. This indicates that the matrix is ​​processed as a column vector. This represents the vectorized form of a super-resolution ISAR image. This represents a diagonal matrix composed of elements from E. This means constructing a diagonal matrix using vectors as diagonal elements. It is by The matrix formed by the diagonal elements has the following form: (8)。 3. The multi-station observation 2D super-resolution ISAR imaging method based on fast IWF according to claim 2, characterized in that, S2 specifically includes: Let the vector of signal precision Each element in the array has a value of 1, and the maximum number of iterations for the fast IWF method is set to 1. , where superscript Represents the initial value. and The dimension of the super-resolution ISAR image. It represents the set of complex numbers.

4. The multi-station observation 2D super-resolution ISAR imaging method based on fast IWF according to claim 3, characterized in that, S3 specifically includes: For j=0,1,2,…, Based on the accuracy of the signal in the j-th iteration Noise variance estimate obtained from auxiliary data ,according to The elements within the auxiliary observation covariance matrix G are quickly computed using 2D-FFT, and then the 2D Levinson-Durbin method is used to calculate... The LC decomposition factor is given by where Q is the observation covariance matrix and H is the 2D dictionary matrix. Let S be the autocorrelation matrix of the signal, S be the permutation matrix consisting of 0 and 1 elements, and I be the identity matrix. The superscript... Represents the j-th iteration. This represents the conjugate transpose operation of a matrix. This represents the transpose operation of a matrix. This represents the matrix inversion operation.

5. The multi-station observation 2D super-resolution ISAR imaging method based on fast IWF according to claim 4, characterized in that, Calculate the auxiliary observation covariance matrix G in S3 and The LC decomposition factor is determined by the following steps: The observation model for Wiener filtering is: (9) in, It is observation data. represents an element in y, and m is the index of that element in y. It is the expected signal. Let represent an element in z, and m be the index of that element in z. Given the observation noise, the optimal weights and the estimated desired signal for Wiener filtering are as follows: (10) (11) in, It is the autocorrelation matrix of the filter input signal, which is obtained. ,in and These are the autocorrelation matrix of the desired signal and the autocorrelation matrix of the noise, respectively. It is the cross-correlation matrix between the input signal and the desired signal. This represents the expectation operation; Substituting the 2D sparse signal reconstruction model represented by formula (7) into the Wiener filter observation model, assuming that the range-space-varying phase error has been compensated, let the echo after range-space-varying phase error compensation be... ,in Representing the conjugate operation of matrices, the echo model, the autocorrelation matrix of the observed signal, the cross-correlation matrix, and the desired signal are respectively in the following forms: (12) (13) (14) (15) in, represent The element in the m-th row is obtained , This refers to the noise variance, which is estimated using auxiliary data in radar imaging echo signals. It is the autocorrelation matrix of the reconstructed signal, constructed using the signal power values, that is, It is a diagonal element a diagonal matrix, where Represents the k-th element in x. Representative signal The accuracy; because Therefore, the signal estimation formula is: (16) Among them, the observation covariance matrix for: (17) Where I represents an identity matrix; Using the permutation matrix S consisting of elements of 1 and 0, Written as: (26) in, The observation covariance matrix, used to assist observation, is a matrix with a fourth-order Toeplitz tensor expansion: (27) Wherein, the inner submatrix of G The form is: (28) Among them, superscript represent Since G has a fourth-order Toeplitz tensor matrix structure, G can be written in two different forms as follows: (29) (30) in: (31) (32) (33) It is a reversed matrix, where the elements on the secondary diagonal are 1s and the other elements are 0s. (34) Applying the block matrix inverse formula to equations (29) and (30) respectively, the inverse matrix of G is obtained as follows: (35) (36) Where 0 represents a zero matrix of appropriate dimensions, and (37) (38) Define a dimension as Lower triangular Toeplitz matrix and a circular matrix : (39) (40) Using formulas (35) and (36), The displacement matrix is: (41) in, and Let represent a lower triangular Toeplitz-block permutation matrix consisting of 1s and 0s, and a cyclic-block permutation matrix consisting of 1s and 0s, respectively. Their forms are: (42) (43) make: (44) (45) (46) in, yes The l-th column vector in the middle, , yes The l-th column vector in the middle, yes The l-th column vector; based on , Written as: (47) Among them, T, P and for The LC decomposition factor was calculated using the 2D Levinson-Durbin method, yielding T, P, and [other parameters]. , and Defined as a TB matrix and a CB matrix respectively: (48) (49)。 6. The multi-station observation 2D super-resolution ISAR imaging method based on fast IWF according to claim 5, characterized in that, S4 specifically includes: Obtain the echo signal after phase error compensation ,in, This represents the echo signal after phase error compensation. It is a matrix composed of Fourier dictionary matrices. Represents the distance-space phase error. Represents the echo signal. Represents the conjugate operation of a matrix; The target rotation speed is estimated using the Tsallis entropy of the image obtained by the fast IWF method as the cost function: (50) in, Rotational speed The estimated value, It is the Tsallis entropy of the image. The form is: (51) The total energy of the image is p is the exponent of the Tsallis entropy. Represents the elements within signal x; The first-order partial derivative of the Tsallis entropy of the image reconstructed by the fast IWF method with respect to the target rotation speed is: (52) The first-order partial derivative of the distance-space phase error with respect to rotational speed is: (53) Represents the operation of taking the real part and solving nonlinear equations. This allows us to obtain the estimated rotational speed of the target.

7. The multi-stationary observation 2D super-resolution ISAR imaging method based on fast IWF according to claim 6, characterized in that, S5 specifically includes: based on The LC decomposition factor and the echo after phase error compensation, according to Using 2D-FFT to quickly compute estimates of sparse signals and according to Update the precision vector of the sparse signal ,in, It is the index value of the signal precision vector; From formula (16), we know that The calculation consists of three steps: , and Substitute formula (47) into The right side of the equation contains products of TB matrices and vectors, and CB matrices and vectors. These are transformed into linear convolutions and circular convolutions, respectively, and then quickly calculated using 2D-FFT. Matrix transformation , The dimension is ,use Constructing a matrix ,in , and then get Finally, use dot product to calculate. .

8. The multi-station observation 2D super-resolution ISAR imaging method based on fast IWF according to claim 7, characterized in that, S6 specifically includes: based on and Calculate the iteration error And based on the iterative convergence threshold To determine if the iteration has converged, if Er > and If Er < ,or If the iteration stops, the final reconstruction result is... Furthermore, based on the final estimated target rotational speed, azimuth calibration is performed.

Images