Multi-residence observation 2D super-resolution ISAR imaging method based on fast IWF

Abstract

The invention relates to a multi-resident observation 2D super-resolution ISAR imaging method based on a fast IWF. The method comprises the following steps: S1, modeling 2D super-resolution ISAR imaging of a multi-resident observation signal; s2, initializing the precision of a signal in the rapid IWF method, and setting the maximum number of iterations of the rapid IWF method; s3, calculating an LC decomposition factor of an inverse matrix G-1 of the auxiliary observation covariance matrix; s4, estimating a target rotating speed for compensating a distance space-variant phase error, and obtaining an observation signal after error compensation; s5, based on the LC decomposition factor of G-1 and the observation signal after distance space-variant phase error compensation, calculating a signal estimation value by using 2D-FFT, and updating the precision of the signal; and S6, repeating the steps S3-S5 to carry out loop iteration, stopping after convergence is achieved, obtaining a super-resolution ISAR image based on the signal estimation value, and realizing transverse calibration based on the target rotation speed estimation value. The method is low in calculation complexity, strong in robustness, high in reconstruction precision and fast in convergence, and can efficiently obtain a super-resolution ISAR image with a good focusing effect and an accurate target size estimation value.

Description

Technical Field

[0001] The present invention belongs to the technical field of computer simulation and method optimization, and particularly relates to a multi-dwell observation 2D (two-dimensional) super-resolution ISAR (inverse synthetic aperture radar) imaging method based on fast IWF (i.e., iterative Wiener filter, also known as iterative Wiener filtering). Fast IWF, that is, Fast IWF, can also be simply referred to as FIWF. This method is applicable to two-dimensional joint super-resolution ISAR imaging of non-cooperative targets under multi-dwell observation. Background Art

[0002] In the prior art, super-resolution images of non-cooperative targets such as aircraft, ships, missiles, and satellites can be obtained through inverse synthetic aperture radar (ISAR) technology, providing strong support for subsequent accurate target recognition and three-dimensional attitude estimation. With the progress and development of the times, ISAR observation targets are also continuously developing towards miniaturization, densification, high speed, high mobility, etc., which puts forward higher requirements for super-resolution ISAR imaging. According to the radar imaging principle, the range resolution and azimuth resolution of the target are respectively proportional to the radar transmit signal bandwidth and the angle of rotation of the target relative to the radar line of sight during the pulse accumulation time. However, in practical applications, a large bandwidth makes the radar sampling rate too high and the data volume too large, increasing the difficulty of radar hardware implementation. For multifunctional phased array radars, ISAR imaging usually jointly processes multi-dwell observation signals coherently to obtain sufficient azimuth resolution. However, the traditional Doppler processing method based on fast Fourier transform (FFT) will generate high grating lobes and side lobes in the azimuth dimension. Therefore, it is of great practical significance to study a robust and efficient two-dimensional (2D) joint super-resolution ISAR imaging method for multi-dwell observation signals.

[0003] At present, considering the sparse electromagnetic scattering characteristics of ISAR observation targets, sparse signal reconstruction methods have become the most effective methods for super-resolution ISAR imaging of sparse aperture (SA) signals. Sparse reconstruction methods can be divided into two categories. The first category is non-compressive sensing methods, and representative methods include Sequential-CLEAN, Relaxation (RELAX), and Gapped-data Amplitude and Phase Estimation (GAPES) methods, etc. These methods reconstruct missing data by obtaining model parameters. The second category is compressive sensing methods. Compressive sensing technology breaks through the limitation of the Nyquist sampling theorem and can achieve super-resolution ISAR imaging when echo data is limited and data is missing. Compressive sensing-based sparse signal reconstruction methods are divided into three categories: greedy methods, p l-norm regularization methods, and Bayesian methods.

[0004] For existing SA-ISAR imaging methods, although non-compressive sensing methods can reconstruct the target ISAR image, they are extremely vulnerable to model errors and noise, and also have a high computational complexity. Although compressive sensing methods can obtain super-resolution ISAR images, they all have certain limitations. Greedy pursuit methods have a low computational complexity but poor reconstruction accuracy and are difficult to meet the requirements of high-precision imaging; regularization methods require manual setting of regularization coefficients and have limited adaptability; Bayesian methods have higher reconstruction accuracy and are more robust than the first two types of methods, but have a high computational complexity and are difficult to perform real-time imaging. Some existing fast sparse Bayesian learning (SBL) methods use approximate methods to replace the time-consuming matrix inversion operation, which improves the computational efficiency, but the reconstruction accuracy will be lost to a certain extent, so the imaging effect is poor when the effective echo data is small. Summary of the Invention

[0005] To solve the above technical problems existing in the prior art, the present invention proposes a multi-dwell observation 2D super-resolution ISAR imaging method based on fast IWF, specifically including the following steps:

[0006] S1. Model the 2D super-resolution ISAR imaging of multi-dwell observation signals;

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

[0008] S3. Calculate the elements within the auxiliary observation covariance matrix using 2D-FFT based on the accuracy of the reconstructed signal and the noise variance estimated using the auxiliary data, and then calculate the inverse matrix G of the auxiliary observation covariance matrix using the 2D Levinson-Durbin method -1 for the LC decomposition factors of;

[0009] S4. Estimate the target rotation speed using the minimum entropy criterion and the trust region method in the image reconstruction iteration to compensate for the range-variant phase error, and obtain the observation signal after error compensation

[0010] S5. Based on the inverse matrix G of the observation covariance matrix -1 for the LC decomposition factors and the observation signal after range-variant phase error compensation, calculate the signal estimate value using 2D-FFT and update the signal accuracy value

[0011] S6. Repeat steps S3 to S5 for cyclic iteration, stop after convergence, obtain the super-resolution ISAR image based on the signal estimate value, and achieve transverse calibration based on the target rotation speed estimate value

[0012] Furthermore, S1 specifically includes:

[0013] First, assume that the radar emits a chirp signal. After de-chirping processing and removing the residual phase term, the target echo signal received by the radar is expressed in the range frequency domain and the azimuth time domain as:

[0014]

[0015] where f ∈ [-B / 2, B / 2], representing the range frequency, B represents the bandwidth, t m represents the slow time, c is the speed of light, f c is the center carrier frequency of the radar. Assume there are I scattering points on the target, δ i is the backscattering coefficient of the i-th scattering point, R i (t m ) is the instantaneous distance between the i-th scattering point and the radar, R i (t m ) is expressed 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 motion of the target relative to the radar is decomposed into a translational component and a rotational component, r(t m) represents the instantaneous slant range change caused by the translational component, x i sinθ(t m ) + y i cosθ(t m ) represents the instantaneous slant range change caused by the rotational component, (x i , y i ) represents the coordinates of the i-th scatterer, θ(t m ) represents the angle of the target relative to the radar, θ(t m ) = ωt m , ω is the rotational speed of the target. Perform a second-order Taylor series expansion on the trigonometric functions in formula (2), and R i (t m ) is re-expressed as:

[0018]

[0019] Assume that the translational error has been compensated and the range linear walk caused by rotation has been corrected by the Keystone transform. By approximation and ignoring the envelope curvature caused by rotation, the echo signal with range-variant phase error is expressed in the range frequency domain and azimuth time domain as:

[0020]

[0021] Assume that the number of effective pulses within the coherent integration 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-variant phase error is obtained as:

[0022]

[0023] Among them, and represent the echo signal, the super-resolution ISAR image, and the complex Gaussian white noise matrix respectively, represents the complete Fourier dictionary matrix, and represent the range-dimensional dictionary matrix and the azimuth-dimensional dictionary matrix respectively, both of which are over-complete Fourier dictionaries. K1 and K2 represent the number of range cells and Doppler cells after super-resolution respectively. If the full-aperture signal contains M pulses, then the super-resolution multiples in the range dimension and azimuth dimension are K1 / N and K2 / M respectively. ⊙ represents the Hadamard product, represents the range-variant phase error caused by target rotation, and the elements in E are Among them, is the coordinate of the n-th range cell, (·) H represents the conjugate transpose operation of the matrix, (·) T represents the transpose operation of the matrix,

[0024] Vectorize formula (6) to obtain the 2D sparse signal reconstruction model under the SA signal:

[0025]

[0026] Among them, represents the 2D Fourier dictionary matrix, represents the Kronecker product, represents the vectorized form of the SA echo signal matrix, represents the vectorized form of the noise matrix. Vec(·) represents vectorizing the matrix column by column, represents the vectorized form of the super-resolution ISAR image, represents a diagonal matrix composed of the elements in E. Diag(·) represents constructing a diagonal matrix with the vector as the diagonal elements, is a matrix composed of F R as diagonal elements, and its form is:

[0027]

[0028] Furthermore, S2 specifically includes:

[0029] Let each element value in the vector constituting the signal precision be 1, and set the maximum number of iterations of the fast IWF method to J max , where the superscript (0) represents the initial value, and K1 and K2 are the dimensions of the super-resolution ISAR image, represents the set of complex numbers.

[0030] Furthermore, S3 specifically includes:

[0031] For j = 0, 1, 2,..., J max , based on the signal precision γ (j) in the j-th iteration and the noise variance estimate value δ obtained from the auxiliary data 2 , according to Q = (HΛH H +δ 2 I) = S T GS, use 2D-FFT to quickly calculate the elements in the auxiliary observation covariance matrix G, and then use the 2D levinson-Durbin method to calculate the LC decomposition factor of G -1 , where Q is the observation covariance matrix, H is the 2D dictionary matrix, Λ is the signal autocorrelation matrix, S is a permutation matrix composed of 0 and 1 elements, I is the identity matrix, the superscript (j) represents the j-th iteration, (·) H represents the conjugate transpose operation of the matrix, (·) T represents the transpose operation of the matrix, (·) -1Represents matrix inversion operation.

[0032] Further, calculate the auxiliary observation covariance matrices G and G in S3 -1 and their LC decomposition factors, which specifically include the following steps:

[0033] The observation model of Wiener filtering is:

[0034] y = z + n (9)

[0035] where y = [y(0) y(1)... y(M - 1)] T is the observed data, y(m) represents the element in y, m is the index of the element in y, z = [z(0) z(1)... z(M - 1)] T is the desired signal, z(m) represents the element in z, m is the index of the element in z, n is the observation noise, and the optimal weight of Wiener filtering and the estimation of the desired signal are respectively:

[0036]

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

[0038] Substitute the 2D sparse signal reconstruction model represented by formula (7) into the observation model of the Wiener filtering. Assuming that the range - variant phase error has been compensated, let the echo after range - variant phase error compensation be where (·) * represents the conjugate operation of the matrix. Then the forms of the echo model, the autocorrelation matrix of the observed signal, the cross - correlation matrix, and the desired signal are respectively:

[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 represents the m-th row element of H, obtaining R nn =δ 2 I, δ 2 is the noise variance. The noise variance of the echo signal in radar imaging is estimated with the aid of auxiliary data. Λ = E[xx H is the autocorrelation matrix of the reconstructed signal. The autocorrelation matrix of the signal is constructed using the signal power value, that is, Λ is a diagonal matrix with diagonal elements γ k =|x k | 2 where x k represents the k-th element in x, and γ k represents the precision of the signal x k ;

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

[0045]

[0046] Among them, the observation covariance matrix Q is:

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

[0048] Among them, I represents an identity matrix;

[0049] Using the permutation matrix S composed of 1 and 0 elements, Q is written as:

[0050] Q = S T GS(26)

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

[0052]

[0053] Among them, the inner sub-matrix of G has the form:

[0054]

[0055] Among them, the superscript <·> represents the element in Q. Since G has a fourth-order Toeplitz tensor matrix structure, G is written in two different forms as follows:

[0056]

[0057] Among them:

[0058]

[0059] I Nq is a reverse matrix, the elements on its secondary diagonal are 1, and other elements are 0, that is

[0060]

[0061] Using the block matrix inversion formula for formulas (29) and (30) respectively, the inverse matrix of G is:

[0062]

[0063] where 0 represents a zero matrix of appropriate dimension, and

[0064]

[0065] Define a lower triangular Toeplitz matrix of dimension M gp ×M gp and a circulant matrix and a circulant matrix

[0066]

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

[0068]

[0069] where and respectively represent a lower triangular Toeplitz-block permutation matrix composed of 1 and 0 elements and a circulant-block permutation matrix composed of 1 and 0 elements, and their forms are:

[0070]

[0071] Let:

[0072]

[0073] where t l (l = 0, 1,..., Nq - 1) is the l-th column vector in T, p l is the l-th column vector in P, is the l-th column vector in;

[0074] Based on G -1 is written as:

[0075]

[0076] Among them, T, P, and are the LC decomposition factors of G, and T, P, and -1 are calculated by the 2D Levinson - Durbin method. and are 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 represents the echo signal after phase error compensation, is a matrix composed of the Fourier dictionary matrix, represents the range - variant phase error, s g represents the echo signal, and (·) * represents the conjugate operation of the matrix;

[0080] Use the Tsallis entropy of the image obtained by the fast IWF method as the cost function to estimate the target rotational speed as:

[0081]

[0082] Among them, is the estimated value of the rotational speed ω, T(X) is the Tsallis entropy of the image, and the form of T(X) is:

[0083]

[0084] Among them, the total energy of the image is p is the exponent of the Tsallis entropy, represents the elements within the 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 rotational speed is:

[0086]

[0087] Among them, the first - order partial derivative of the range - variant phase error with respect to the rotational speed is:

[0088]

[0089] Re{·} represents the operation of taking the real part, and solve the nonlinear equation to obtain the estimated rotational speed of the target.

[0090] Further, S5 specifically includes:

[0091] Based on the LC decomposition factors of G -1 and the echo after phase error compensation, according to Use 2D-FFT to quickly calculate the estimated value of the sparse signal and update the precision vector γ of the sparse signal according to γ k = |x k | 2 where k = 0, 1, …, K-1 is the index value of the signal precision vector; (j+1)

[0092] As known from formula (16), the calculation of is divided into three steps: φ = Q -1 y g = S T G -1 Sy g , and Substitute formula (47) into φ. The right side of the equation is the product of some TB matrices and vectors and the product of CB matrices and vectors, which are respectively transformed into linear convolution and circular convolution, and calculated quickly using 2D-FFT to transform the φ matrix into Φ = [φ0 φ1 … φ q-1 , and the dimension of φ i is N×M gp , construct a matrix using Φ where then get Finally, calculate using dot product

[0093] Further, S6 specifically includes:

[0094] Based on and calculate the iterative error and judge whether the iteration converges based on the iterative convergence threshold η. If Er > η and j ≤ J max , then perform the next loop iteration; if Er < η, or j > J max , then stop the iteration, and the final reconstruction result is And perform azimuth calibration based on the finally estimated target rotation speed.

[0095] The present invention can be used for the detection, imaging and recognition of mobile targets in complex electromagnetic environments, and has the advantages of low computational complexity, strong robustness, high reconstruction accuracy and fast convergence, and can efficiently obtain a super-resolution ISAR image with good focusing effect and a relatively accurate target size estimated value. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0097] Figure 2 is the diagram of the 2D multi-dwell echo signal model;

[0098] Figure 3 are the 2D sparse signal reconstruction results of Model A and multi-dwell observation data. Among them, (a) is the scattering point model of letter A, (b) is the 2D multi-dwell 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 are the performance curve diagrams of the method under different q, MR, and SNR. Among them, (a) is the calculation time of each method under different q values, (b) is the reconstruction error of each method under different q values, (c) is the calculation time of each method under different MR values, (d) is the reconstruction error of each method under different MR values, (e) is the calculation time of each method under different SNR values, and (f) is the calculation time and reconstruction error of each method under different SNR values;

[0100] Figure 5 are the scattering point aircraft model, its 2D echo signal, and the RD imaging results. Among them, (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-dwell echo signal, and (e) is the RD imaging result of the multi-dwell echo signal;

[0101] Figure 6 are different types of multi-dwell echo signals and the ISAR imaging results of OMP and FIWF. Among them, (a)-(c) are the echo signal when the number of dwells is 4, the OMP imaging result of the ideal echo, and the imaging result of FIWF-ME-TR, and (d)-(f) are the echo signal when the number of dwells is 8, the OMP imaging result of the ideal echo, and the imaging result of FIWF-ME-TR;

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

[0103] Figure 8Echo signals and ISAR images obtained by OMP and FIWF at different SNRs. Among them, (a)-(c) are the echo signal, the OMP imaging result of the ideal echo, and the imaging result of FIWF-ME-TR when the SNR is 15 dB, respectively; (d)-(f) are the echo signal, the OMP imaging result of the ideal echo, and the imaging result of FIWF-ME-TR when the SNR is 5 dB, respectively;

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

[0105] Figure 10 Optical image, HRRPs and full-aperture RD imaging results of Boeing 787. Among them, (a) is the optical image of 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 Multi-dwell HRRPs and ISAR images of Boeing-787 obtained by GP-FIWF, SAMP and FIWF-ME-TR. Among them, (a)-(d) are the HRRPs with MR of 70% and q of 2, and the corresponding GP-FIWF imaging result, SAMP imaging result and IWF-ME-TR imaging result, respectively; (e)-(h) are the HRRPs with 75% and q of 4, and the corresponding GP-FIWF imaging result, SAMP imaging result and FIWF-ME-TR imaging result, respectively. Detailed implementation manners

[0107] To make the objectives, contents and advantages of the present invention clearer, the following further describes in detail the specific implementation manners of the present invention in conjunction with the accompanying drawings and embodiments.

[0108] In view of the fact that the traditional imaging method has failed for the multi-dwell echo signals used for ISAR imaging received by a multi-functional radar, and in order to relieve the pressure of radar hardware design, the present invention proposes a multi-dwell observation 2D super-resolution ISAR imaging method based on the iterative Wiener filter (IWF).

[0109] Based on the fact that the noise variance based on radar observations can be estimated by auxiliary data and the amplitudes of target scatterers can be regarded as unknown deterministic variables, the present invention proposes another implementation method of the minimum mean square error (MMSE) criterion estimator, namely the IWF method. Compared with SBL, this method has the advantages of small computational complexity, fast convergence speed, and high reconstruction accuracy.

[0110] Moreover, the present invention utilizes the fact that in the IWF reconstruction of 2D signals for multi-dwell observation signals, the covariance matrix in the iteration has an expansion structure of a fourth-order Toeplitz tensor, and proposes a fast IWF method, which can efficiently and accurately realize the 2D sparse signal reconstruction of multi-dwell observation data.

[0111] In addition, in order to make the fast IWF method compatible with the range-variant phase error compensation method, the present invention designs a new range-variant autofocus method. In the fast IWF-based image reconstruction, the target rotation speed is estimated by using image minimum entropy and trust region methods, and 2D super-resolution ISAR imaging and azimuth calibration can be jointly realized. The present invention has the advantages of low computational complexity, high reconstruction accuracy, small occupied computing memory, and good noise suppression ability.

[0112] The present invention utilizes the fact that the noise variance of radar observations can be estimated by auxiliary data and the amplitudes of target scatterers can be regarded as unknown deterministic variables, and proposes the IWF method. And for multi-dwell observation signals, a 2D super-resolution ISAR imaging method based on the fast IWF method is proposed. In the proposed fast IWF method (abbreviated as FIWF), by using the fact that the covariance matrix in the IWF reconstruction of 2D signals has a fourth-order Toeplitz tensor expansion structure, the present invention designs a lower-triangular-Toeplitz-cyclic (LC) decomposition method to avoid time-consuming matrix inversion operations. Based on the LC decomposition factors, the signal can be quickly reconstructed using FFT in each iteration, greatly reducing the computational complexity of the method.

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

[0114] S1. Model the 2D super-resolution ISAR imaging of multi-dwell observation signals;

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

[0116] Let each element value in the vector constituting the signal accuracy be 1, and set the maximum number of iterations of the FIWF method to Jmax Among them, 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 the auxiliary data, calculate the elements within the auxiliary observation covariance matrix using 2D-FFT, and then calculate the LC decomposition factors of the inverse matrix of the auxiliary observation covariance matrix 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) and the noise variance estimate δ obtained from the auxiliary data 2 , according to Q=(HΛH H +δ 2 I)=S T GS, quickly calculate the elements within the auxiliary covariance matrix G using 2D-FFT, and then calculate the LC decomposition factors of G -1 . Among them, Q is called the observation covariance matrix, H is the 2D dictionary matrix, Λ is the signal autocorrelation matrix, S is the permutation matrix composed of 0 and 1 elements, I represents the identity matrix, the superscript (j) represents the j-th iteration, (·) H represents the conjugate transpose operation of the matrix, (·) T represents the transpose operation of the matrix, (·) -1 represents the matrix inverse operation;

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

[0120] Obtain the echo signal after phase error compensation Among them, y g represents the echo signal after phase error compensation, is the matrix composed of the Fourier dictionary matrix, represents the range-variant phase error, s g represents the echo signal, (·) * represents the conjugate operation of the matrix;

[0121] S5. Based on the LC decomposition factors of the inverse matrix G -1 of the observation covariance matrix and the observation signal after range-variant phase error compensation, calculate the signal estimate using 2D-FFT and update the accuracy value of the signal.

[0122] Based on G -1The echo after LC decomposition factor and phase error compensation, according to Use 2D-FFT to quickly calculate the estimated value of the sparse signal And according to γ k = |x k | 2 Update the precision vector γ of the signal (j+1) . Where k = 0, 1, …, K-1 is the index value of the signal precision vector;

[0123] S6. Repeat the above steps S3-S5 for cyclic iteration and stop after convergence. Based on the signal estimated value, obtain the super-resolution ISAR image and achieve transverse calibration based on the target rotation speed estimated value.

[0124] Based on and Calculate the iteration error And judge whether the iteration converges based on the iteration convergence threshold η. If Er > η and j ≤ J max , then perform the next cyclic iteration; if Er < η, or j > J max , then stop the iteration, and the final reconstruction result is And based on the finally estimated target rotation speed, perform azimuth calibration.

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

[0126] First, assume that the radar emits a chirp signal. After de-chirping processing and removing the residual phase term, the target echo signal received by the radar can be expressed in the range frequency domain and azimuth time domain as

[0127]

[0128] where f ∈ [-B / 2, B / 2] represents the range frequency, B represents the bandwidth, t m represents the slow time, c is the speed of light, f c is the center carrier frequency of the radar. Assume that there are I scattering points on the target, δ i is the backscattering coefficient of the i-th scattering point, R i (t m ) is 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] Among them, R0 represents the initial distance. The motion of the target relative to the radar can be decomposed into a translational component and a rotational component. r(t m ) represents the instantaneous slant-range change caused by the translational component. The latter two terms represent the instantaneous slant-range change caused by the rotational component. (x i , y i ) represents the coordinates of the i-th scatterer, and θ(t m ) represents the angle of rotation of the target relative to the radar. Since the coherent integration angle of ISAR imaging is relatively small (about 5°), during the imaging coherent integration time, the motion of the maneuvering target can be approximately regarded as uniform rotation, that is, θ(t m ) = ωt m , where ω is the rotational speed of the target. To represent the instantaneous distance of the target more precisely, the trigonometric function is expanded by the second-order Taylor series, and R i (t m ) can be re-expressed as

[0131]

[0132] Substitute formula (6) into formula (4), and assume that the translational error has been compensated and the range linear migration has been corrected by the Keystone transform. After approximation, the echo signal is expressed as After approximation, the echo signal is expressed as

[0133]

[0134] It can be seen from formula (7) that the quadratic term will cause range curvature and range-variant phase error in the azimuth dimension. Since the rotational speed of the maneuvering target is small, the range curvature can be ignored. However, the phase error compensation accuracy is relative to the wavelength level, so the range-variant phase error compensation must be carried out before imaging. After eliminating the constant term, the echo signal with range-variant phase error can be expressed in the range frequency domain and azimuth time domain as

[0135]

[0136] Assume that the number of effective pulses within the coherent integration time of the sparse aperture signal is L, and the number of range frequency points is N. The discretized model of 2D-SA (two-dimensional sparse aperture) super-resolution ISAR imaging with range-variant phase error can be obtained as follows:

[0137]

[0138] Among them, and represent the echo signal, the super-resolution ISAR image, and the complex Gaussian white noise matrix respectively, represents the complete Fourier dictionary matrix, and They represent the range - dimension dictionary matrix and the azimuth - dimension dictionary matrix respectively, both of which are over - complete Fourier dictionaries. K1 and K2 represent the number of range cells and Doppler cells after super - resolution respectively. If the full - aperture signal contains M pulses, then the super - resolution multiples in the range dimension and azimuth dimension are K1 / N and K2 / M respectively. ⊙ represents the Hadamard product. represents the range - variant phase error caused by target rotation. The elements in E are Among them, is the coordinate of the nth range cell and is measurable. Therefore, only by solving the target rotation speed ω can E be calculated.

[0139] Vectorize formula (9), and the 2D sparse signal reconstruction model under the SA (sparse aperture) signal can be obtained:

[0140]

[0141] Among them, represents the 2D Fourier dictionary matrix, represents the Kronecker product, represents the vectorized form of the SA echo signal matrix, represents the vectorized form of the noise matrix. Vec(·) represents vectorizing the matrix column - by - column. represents the vectorized form of the super - resolution ISAR image, and its specific form can be expressed as Among them, the internal elements of x are k2 = 0, 1, …, K2 - 1. represents a diagonal matrix composed of the elements in E at the SA time period. Diag(·) represents constructing a diagonal matrix with the vector as the diagonal elements. is a matrix composed of F R as the diagonal elements, and its form is

[0142]

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

[0144] The observation model of Wiener filtering is:

[0145] y = z + n (12)

[0146] Among them, y = [y(0) y(1) … y(M - 1)] TIt is the observed data. y(m) represents the element in y, and m is the index of the element in y. z = [z(0) z(1) … z(M-1)] T is the desired signal. z(m) represents the element in z, and m is the index of the element in z. n is the observed noise. The essence of Wiener filtering is to use the observed data for the linear estimation of the desired signal.

[0147] The optimal weight of Wiener filtering and the estimation of the desired signal are respectively

[0148]

[0149] where, R yy 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 are respectively the autocorrelation matrix of the desired signal and the autocorrelation matrix of the noise. r yz is the cross-correlation matrix between the input signal and the desired signal. E[·] represents the expectation operation.

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

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

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

[0153]

[0154] where, H m represents the m-th row element of H. Since the noise is a complex Gaussian white noise, we can obtain R nn = δ 2 I, δ 2 is the noise variance. In radar imaging, the noise variance of the echo signal can be estimated with the help of auxiliary data. Λ = E[xx His the autocorrelation matrix of the reconstructed signal. It can be seen from Equation (18) that Wiener filtering requires the known 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 an initial value and iteratively solving according to known parameters. This is the key to realizing super-resolution imaging using IWF. Since the amplitude of the target scatterer can be regarded as an unknown deterministic variable, the autocorrelation matrix of the signal can be constructed using the signal power value. Therefore, Λ is a diagonal matrix with diagonal elements γ k = |x k | 2 , where x k represents the k-th element in x, and γ k represents the precision of x k . The Iterative Adaptive Approach (IAA) method is also the same.

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

[0156]

[0157] where the observation covariance matrix Q is

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

[0159] where I represents an identity matrix.

[0160] Observing Equation (19), it can be found that the signal estimation formula in the IWF method is the same as that in the SBL method. Different from SBL, the calculation methods of the noise variance and the 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 there are, the closer the estimated value of the noise variance 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 is that, aiming at the time-consuming matrix inversion operation in the iteration of the IWF method, in order to reduce the computational complexity of the method, taking advantage of the fact that Q has a fourth-order Toeplitz tensor matrix structure under multi-dwelling echo signals, a fast IWF method is designed to obtain the elements in the auxiliary observation covariance matrix G calculated quickly using FFT in S3 and the LC decomposition factors of G calculated using the 2D Levinson-Durbin (LD) method. -1

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

[0163] First, assume that the geometric structure of the 2D multi-dwell echo signal is as Figure 2 shown. Assume that the number of dwells is q, one dwell contains M gp pulses, and there are M mp pulses between two dwells. Then, we can get M sp = M gp + M mp , and L = qM gp . The data missing rate (MR) is defined as MR = M mp / M sp .

[0164] For multi-dwell echo signals, when using IWF to reconstruct the 2D signal, the azimuthal dictionary matrix H A and the range dictionary matrix H R are a partial Fourier dictionary matrix and a complete Fourier dictionary matrix respectively. The form of H A is:

[0165]

[0166] where H A(i) represents the dictionary matrix corresponding to the i-th segment of dwell data. The forms of the column elements in H R and H A(i) are respectively

[0167]

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

[0169] Let the observation covariance matrix Q = R + δ 2 I, where R is the signal covariance matrix. Substituting the dictionary H, we can get that R is a matrix with a fourth-order Toeplitz tensor expansion form:

[0170]

[0171] where, the sub-matrix in R is in the form of a third-order Toeplitz tensor expansion, that is, a Toeplitz-block-Toeplitz matrix, and its form is:

[0172]

[0173] where, Inner sub - matrix is a Toeplitz matrix, and its form is:

[0174]

[0175]

[0176] Here represents the row vector of the (j + 1)-th row of H A(i+1) , and the elements in can be expressed as

[0177]

[0178] Obviously , the elements in

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

[0180] A further technical solution is that with the help of a permutation matrix, an auxiliary observation covariance matrix G can be constructed. The specific steps include:

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

[0182] Q = S T GS (29)

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

[0184]

[0185] where the inner sub - matrix of G has the form of:

[0186]

[0187] Here the superscript <·> represents the elements in Q. For example represents the element in Obviously, the elements in G can also be quickly calculated by 2D - FFT.

[0188] A further technical solution is that using the fourth - order Toeplitz tensor matrix structure of G, the LC decomposition can be used to calculate G -1 , and obtain the LC decomposition factors of G -1 . Calculate G -1 and G -1The steps of 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] Where:

[0192]

[0193] is a reverse matrix, with elements on the sub-diagonal being 1 and other elements being 0, that is

[0194]

[0195] Using the block matrix inversion formula for equations (32) and (33) respectively, the inverse matrix of G can be obtained as

[0196]

[0197] Where:

[0198]

[0199]

[0200] Define a lower triangular Toeplitz matrix of dimension M gp ×M gp and a circulant matrix and a circulant matrix

[0201]

[0202] Using equations (38) and (39), the displacement matrix of G -1 is:

[0203]

[0204] Where, and represent a lower triangular Toeplitz-block (TB) permutation matrix composed of 1 and 0 elements and a cyclic-block (CB) permutation matrix composed of 1 and 0 elements respectively, and their forms are:

[0205]

[0206] Let

[0207]

[0208] where \(t\) l (\(l = 0, 1, \ldots, N_q - 1\)) is the \(l\)-th column vector in \(T\) and can be written as \(t\) i represents the \((i + 1)\)-th sub-vector in \(t\) with dimension \(N_q\), which can be written as l \(t\) \(t\) i,j represents the \((j + 1)\)-th sub-vector in \(t\) with dimension \(N\), and its form is \(t\) i \( = [t_0\ t_1\ \ldots\ t\) i,j N-1 T . \(p\) l (\(l = 0, 1, \ldots, N_q - 1\)) is the \(l\)-th column vector in \(P\), is the \(l\)-th column vector in, and its form is the same as \(t\) l .

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

[0210]

[0211] where \(T\), \(P\) and are called the LC decomposition factors of \(G\) -1 and can be calculated by the 2D levinson - Durbin (LD) method. and are respectively defined as a TB matrix and a CB matrix:

[0212]

[0213]

[0214] From formula (50), it can be seen that \(G\) -1 can be obtained by LC decomposition, and the right side of \(G\) -1 is the product operation of the TB matrix and the CB matrix. Therefore, formula (50) is called the LC decomposition formula of \(G\) -1 , and the process of calculating \(G\) -1 based on the \(G\) structure is called LC decomposition. The displacement rank of \(G\) -1 is \(2N_q\). Using \(Q\) -1 = \(S\) T G -1 S, it can be known that the displacement rank of \(Q\) -1 is also \(N_q\).

[0215] ​​A further technical solution is to estimate the target rotation speed in the image reconstruction iteration based on the rough estimate of the target rotation speed, and then compensate the range-variant phase error to obtain the echo signal after phase error compensation in S4. The specific steps include:

[0216] As can be seen from the super-resolution imaging model shown in formula (10), the range-variant 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 cell is measurable, the target rotation speed can be obtained, and then the accurate compensation of the range-variant phase error can be realized. In imaging, image entropy is the most commonly used index to measure the imaging focusing effect. Based on the minimum entropy criterion, this paper transforms the target rotation speed estimation into a non-linear optimization problem as shown in formula (50).

[0217] Using the Tsallis entropy of the image obtained by fast IWF as the cost function to estimate the target rotation speed as

[0218]

[0219] Where is the estimated value of the rotation speed ω, T(X) is the Tsallis entropy of the image, and the form of T(X) is:

[0220]

[0221] Where the total image energy is p is the exponent of the Tsallis entropy, usually taking the value of 1.6. represents the elements in the 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 rotation speed is:

[0223]

[0224] Where the first-order partial derivative of the range-variant phase error with respect to the rotation speed is:

[0225]

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

[0227] Using the trust region algorithm to solve the non-linear equation Then the solution of this non-linear equation is the estimated rotation speed of the target. The parameter estimation method based on the minimum entropy criterion and the trust region method is simply referred to as ME-TR.

[0228] A further technical solution is to use the LC decomposition factor of G -1 and the echo signal y after range-variant phase error compensation g, obtain the estimated value of the sparse signal calculated quickly using FFT in S5 Quick calculation The steps of

[0229] As can be seen from formula (19), The calculation of -1 y g = S T G -1 Sy g , and Substitute formula (50) into φ. The right side of the equation is the product of some TB matrices and vectors and the product of CB matrices and vectors, which can be respectively transformed into linear convolution and circular convolution. Therefore, 2D-FFT can be used for quick calculation. Transform the φ matrix into Φ = [φ0 φ1 … φ q-1 , the dimension of φ i is N×M gp . Use Φ to construct the matrix where Then we can get Finally, use dot product to calculate

[0230] The FIWF method uses a fourth-order Toeplitz tensor expansion structure, uses LC decomposition to find the inverse matrix, and almost all operations can be quickly calculated using FFT except for the solution of LC decomposition factors. It can be seen that FIWF can efficiently and accurately achieve 2D super-resolution ISAR imaging. The displacement rank of the inverse matrix of the auxiliary observation covariance matrix in the FIWF method is 2Nq. The larger the displacement rank, the higher 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 azimuth calibration based on the finally estimated target rotation speed, which provides strong support for subsequent target feature extraction and recognition.

[0232] To verify the effectiveness of the multi-dwell observation super-resolution ISAR imaging method based on the fast IWF method proposed by the present invention, the following simulation and actual measurement comparison experiments are set up. Since there are few 2D super-resolution imaging methods proposed at present, only OMP and SBL are used as comparison methods.

[0233] Example 1:

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

[0235]

[0236] where and x represent the reconstructed signal and the true signal respectively.

[0237] A multi-dwell observation data is constructed by missing 50% of the complete observation pulse period stage by stage, as shown in Figure 3 (b). Figure 3 (c)-(f) are the images reconstructed by 2D-FFT, OMP, SBL, and IWF / FIWF respectively. The calculation time and reconstruction error of the methods are also marked on the images. It can be seen from Figure 3 that under multi-dwell observation data, the grating lobes of the 2D-FFT reconstruction result are obvious and ghosting occurs. The reconstruction accuracy of OMP is low, and the reconstruction effect is poor in the places where the target scattering points are dense. There are some false points and missing points. However, the reconstruction effects of SBL and IWF are particularly good, and clear target images can be obtained. Compared with SBL, the reconstruction accuracy of IWF is higher, and the calculation efficiency of FIWF is two orders of magnitude higher than that of IWF.

[0238] Example 2:

[0239] In this example, in order to quantitatively analyze the performance of OMP, SBL, and the proposed IWF / FIWF method in the present invention, 100 Monte Carlo experiments are carried out using the 2D multi-dwell observation data of the A letter scattering point model, and the performance curves of the methods under different q, MR, and SNR are obtained, as shown in Figure 4 . Figure 4 (a) and Figure 4 (b) respectively show the reconstruction time and reconstruction error of the methods under different types of multi-dwell data. In this experiment, the number of range frequency points and pulses of the 2D observation are 32 and 64 respectively, 50% of the azimuth pulse period is missing, and the signal-to-noise ratio is 10 dB. From Figure 4 (a) and Figure 4(b) It can be seen that the OMP has high computational efficiency but poor reconstruction accuracy. The SBL has high reconstruction accuracy but long computational time. The IWF has higher computational efficiency than the SBL and higher reconstruction accuracy than the SBL. The computational time of the FIWF is much less than that of the IWF, but the reconstruction accuracy is the same. The fewer the number of dwells of the observed data, the higher the computational efficiency of the 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 the FIWF. Figure 4 (c) and Figure 4 (d) show the computational time and reconstruction error of the method with different azimuth apertures MR. The number of range frequency points and pulses of the 2D observation are 32 and 64 respectively, and the signal-to-noise ratio is 10 dB. From Figure 4 (c) and Figure 4 (d) it can be seen that the more pulses are missing, the shorter the computational time and the larger the reconstruction error. When less data is missing, the reconstruction error of the SBL is very small. When more data is missing, the reconstruction error is also acceptable. Since the noise variance in the IWF is closer to the true value, the reconstruction error of the IWF is smaller than that of the SBL. Due to the relatively dense model scatterers, the reconstruction errors of the OMP are all large. The computational time and reconstruction error of the method under different SNRs are shown in Figure 4 (e) and Figure 4 (f). The MR is 50%, and the number of range sampling points and pulses of the 2D observation 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 the IWF is less than that of the SBL. However, when the signal-to-noise ratio is 0 dB, the reconstruction errors of the SBL and IWF / FIWF are relatively large. This is because the 2D observation data in the experiment contains noise. The more sampling points, the more energy accumulation of the reconstructed signal. However, the number of 2D sampling points in the experiment is small, so the reconstruction error is relatively large when the signal-to-noise ratio is low.

[0240] Example 3:

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

[0242]

[0243] Among them, and ω respectively represent the estimated value and the true value of the target rotational speed.

[0244] In this example, Figure 5 (b) and Figure 5 (c) give the complete echo signal with range-variant phase error and the RD imaging result at a signal-to-noise ratio of 10 dB, which are used for comparison with the SA signal imaging result. Figure 5 (d) and Figure 5 (e) give the multi-dwell echo signal with 75% pulse missing and the RD imaging result. It can be seen from the figure that due to the uncompensated range-variant phase error, the RD image is significantly defocused, 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-dwell echo data, different degrees of azimuth aperture missing, and different signal-to-noise ratios.

[0245] Table 1 Radar and target motion parameters

[0246]

[0247]

[0248] Figure 6 Shows different types of multi-dwell echo signals and super-resolution ISAR images. SNR and MR are set to 10 dB and 75% respectively. Figure 6 (a)-(c) are respectively the echo signal with 4 dwells, the OMP imaging result of the ideal echo signal (the ideal echo signal refers to the echo signal without range-variant error), and the FIWF-ME-TR imaging result (the parameter estimation method based on the minimum entropy and trust region method is called ME-TR, and the method using the FIWF method and the parameter estimation method based on the minimum entropy and trust region method is abbreviated as FIWF-ME-TR). Figure 6 (d)-(f) are respectively the echo signal with 8 dwells, the OMP imaging result of the ideal echo signal, and the FIWF-ME-TR imaging result. From Figure 6 it can be seen that for the incomplete ideal echo signal, OMP has missing points and false points. However, the fast IWF method can obtain a well-focused high-quality ISAR image, and the target rotational speed estimated using ME-TR is very close to the true value. Since q has no influence on the image entropy and the target rotational speed estimation, Figure 9(a) shows only the computation times of OMP and FIWF-ME-TR for different q. The analysis of the computation times indicates that OMP is computationally efficient, and FIWF has higher computational efficiency in processing echo signals with fewer dwell times, which is consistent with the theoretical analysis of the method in Section III.

[0249] The imaging results under different MR are as Figure 7 shown. Combining 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 contour of the target cannot be seen in the imaging result of OMP, while only a few scattering points are reconstructed incorrectly in the IWF imaging result. Figure 9 (b) illustrates the influence of MR on the method. The more data is missing, the shorter the computation time of the method, the larger the image entropy, and the greater the deviation between the target rotation speed estimated by ME-TR and the true value. However, when more data is missing, the estimation error of ME-TR is acceptable. The image entropy reconstructed by IWF is smaller than that of OMP, indicating that the target image obtained by IWF has a better focusing effect. The target rotation speed estimated in the image reconstruction based on IWF is also closer to the true value.

[0250] Figure 8 (a)-(c) are the echo signal, the OMP imaging result, and the FIWF-ME-TR imaging result at SNR = 15 dB, respectively. Figure 8 (d)-(f) show the echo and the ISAR image at SNR = 5 dB, respectively. From Figure 6 (a)-(c) and Figure 8 it can be seen that under the same conditions, the higher the SNR, the better the imaging effect of the method. However, even at high SNR, there will still be missing points and false points in the target image reconstructed by OMP due to more missing data. While IWF can obtain a clear target image, and only a very small number of scattering points cannot be accurately reconstructed at low SNR. Figure 9 (c) details the imaging time, the image entropy, and the REE of the target rotation speed estimated by ME-TR for various methods at different SNRs. From Figure 9 (c), it can be seen that the larger the SNR, the smaller 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 iterations of FIWF to 40 times, the computation times are basically the same at different SNRs.

[0251] Example 4:

[0252] In this example, the processing results of the measured data of the Boeing-787 aircraft will be given to verify the effectiveness of the proposed method. The main parameters of the radar system are as follows: the carrier frequency is 10 GHz, the bandwidth is 400 MHz, and the PRF is 100 Hz. The full-aperture data includes 1024 range sampling points and 256 pulses. The Boeing 787 optical image, high-resolution range profile sequence (HRRPs), and RD imaging results at full aperture are as Figure 10 shown. The one-dimensional super-resolution ISAR imaging method based on fast IWF, abbreviated as GP-FIWF, and SAMP are used as comparison methods. To avoid the influence of range-variant phase error phase, both SAMP and GP-FIWF process the echoes after phase error correction according to the target rotation speed estimated by the method proposed in the present invention. Based on the full-aperture data, two different types of multi-dwell observation data are simulated. Figure 11 (a)-(d) show the HRRPs with 2 dwells and 70% pulse missing, as well as the imaging results of GP-FIWF, SAMP, and FIWF-ME-TR. Figure 11 (e)-(h) show the HRRPs with 4 dwells and 75% pulse missing, and the imaging results of GP-FIWF, SAMP, and FIWF-ME-TR. It can be seen from Figure 11 the figure that although the general shape of the target can be seen in the image obtained by SAMP, 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 after 2D super-resolution is clearer, especially at the wing of the 787 aircraft marked by the circle.

[0253] The calculation time of the above measurement data experiment is shown in Table 2. Combining Table 2 and the above ISAR images of Boeing-787, it can be concluded that the fast IWF method proposed in the present 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] In view of the multi-dwell observation signals received by a multi-functional radar for imaging and to reduce the pressure on radar hardware design, the present invention proposes a fast IWF method for efficiently and highly accurately realizing 2D super-resolution ISAR imaging of multi-dwell observation data. In the proposed fast IWF method, taking advantage of the fact that the observed covariance matrix has a fourth-order Toeplitz tensor expansion form in the IWF iteration, a LC decomposition method is designed to find the inverse matrix of the auxiliary covariance matrix, avoiding time-consuming matrix inversion, and based on the LC decomposition factors, almost all the remaining operations in the IWF method can be quickly calculated using FFT. It can be seen that the proposed fast IWF method has the advantages of low computational complexity, high reconstruction accuracy, strong robustness and fast convergence speed. In addition, in image reconstruction, the target rotation speed is obtained by minimizing the image entropy, and then the range-variant phase error is compensated and the azimuth calibration of the target is achieved.

[0257] The present invention proposes an IWF method and a fast IWF method for multi-dwell 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, the present invention designs a LC decomposition to find the inverse matrix of a matrix with a fourth-order Toeplitz tensor expansion form, avoiding time-consuming matrix inversion operations.

[0258] In order to make the proposed fast IWF method compatible with the range-variant autofocus method, the present invention designs a new autofocus method. In the image reconstruction iteration, the target rotation speed is obtained by minimizing the image entropy, and then the range-variant phase error is compensated and the azimuth calibration is achieved.

[0259] The above are only the preferred embodiments of the present invention. It should be noted that for those of ordinary skill in the art, without departing from the technical principle of the present invention, several improvements and modifications can be made, and these improvements and modifications should also be regarded as the protection scope of the present invention.

Claims

1. A multi-resident observation 2D super-resolution ISAR imaging method based on fast IWF, characterized in that It includes the following steps: S1. Model the 2D super-resolution ISAR imaging of multi-dwelling observation signals; S2. Initialize the precision of the signals in the fast IWF method and set the maximum number of iterations of the fast IWF method; S3. Based on the accuracy of the reconstructed signal and the noise variance estimated using the auxiliary data, calculate the elements within the auxiliary observation covariance matrix using 2D-FFT, and then calculate the inverse matrix G of the auxiliary observation covariance matrix using the 2D Levinson-Durbin method -1 of the LC decomposition factor; S4. In the image reconstruction iteration, use the minimum entropy criterion and the trust region method to estimate the target rotation speed, compensate for the range-variant phase error, and obtain the observation signal after error compensation; S5. Inverse matrix G of the observation covariance matrix -1 Using the LC decomposition factor and the observation signal after distance-variant phase error compensation, calculate the signal estimate value by 2D-FFT, and update the precision value of the signal; S6. Repeat steps S3 to S5 for cyclic iteration, stop after convergence, obtain the super-resolution ISAR image based on the signal estimation value, and achieve transverse calibration based on the target rotation speed estimation value.

2. A multi-dwelling observation 2D super-resolution ISAR imaging method based on fast IWF according to claim 1, characterized in that, S1 specifically includes: First, assume that the radar emits a chirp signal. After de-chirping processing and removing the residual phase term, the target echo signal received by the radar is expressed in the range frequency domain and azimuth time domain as: Among them, f ∈ [-B / 2, B / 2], representing the range frequency, B represents the bandwidth, t m represents the slow time, c is the speed of light, f c is the center carrier frequency of the radar. Assuming there are I scattering points on the target, δ i is the backscattering coefficient of the i-th scattering point, R i (t m ) is the instantaneous distance between the i-th scattering point and the radar, R i (t m ) is expressed as: R i (t m ) = R0 + r(t m ) + x i sinθ(t m ) + y i cosθ(t m ) (2) Among them, R0 represents the initial distance. The movement of the target relative to the radar is decomposed into a translational component and a rotational component. r(t m ) represents the instantaneous slant range change caused by the translational component. x i sinθ(t m ) + y i cosθ(t m ) represents the instantaneous slant range change caused by the rotational component. (x i , y i ) represents the coordinates of the i-th scattering point. θ(t m ) represents the angle of rotation of the target relative to the radar. θ(t m ) = ωt m , where ω is the rotational speed of the target. Perform a second-order Taylor series expansion on the trigonometric function in formula (2), and R i (t m ) is re-expressed as: Assume that the translational error has been compensated and the range linear walk caused by rotation has been corrected by the Keystone transform. By Approximating and neglecting the envelope curvature caused by rotation, the echo signal with range-variant phase error is expressed in the range-frequency domain and azimuth-time domain as follows: Assume that the number of effective pulses within the coherent integration 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-variant phase error is obtained as: Among them, and represent the echo signal, the super-resolution ISAR image, and the complex Gaussian white noise matrix respectively, represents the complete Fourier dictionary matrix, and represent the range-dimensional dictionary matrix and the azimuth-dimensional dictionary matrix respectively, both of which are over-complete Fourier dictionaries. K1 and K2 represent the number of range cells and Doppler cells after super-resolution respectively. If the full-aperture signal contains M pulses, then the super-resolution multiples in the range dimension and the azimuth dimension are K1 / N and K2 / M respectively. ⊙ represents the Hadamard product, represents the range-variant phase error caused by target rotation. The elements in E are Among them, is the coordinate of the nth range cell, (·) H represents the conjugate transpose operation of the matrix, (·) T represents the transpose operation of the matrix, Vectorize formula (6) to obtain the 2D sparse signal reconstruction model under SA signals: Among them, represents the 2D Fourier dictionary matrix, represents the Kronecker product, represents the vectorized form of the SA echo signal matrix, represents the vectorized form of the noise matrix, and Vec(·) represents vectorizing the matrix by columns, represents the vectorized form of the super-resolution ISAR image, represents a diagonal matrix composed of the elements in E. Diag represents constructing a diagonal matrix with the vector as the diagonal elements, is a matrix composed of F R as the diagonal elements, and its form is:

3. A method for multi-resident observation 2D super-resolution ISAR imaging based on fast IWF according to claim 2, characterized in that, S2 specifically includes: Let each element value in the vector composed of signal precision be 1, and set the maximum number of iterations of the fast IWF method to J max , where the superscript (0) represents the initial value, K1 and K2 are the dimensions of the super-resolution ISAR image, represents the set of complex numbers.

4. A method for multi-dwelling observation 2D super-resolution ISAR imaging based on fast IWF according to claim 3, characterized in that, S3 specifically includes: For j = 0, 1, 2, …, J max , based on the precision γ of the signal in the j-th iteration (j) and the estimated noise variance δ obtained from the auxiliary data 2 , according to Q = (HΛH H +δ 2 I)=S T GS, use 2D-FFT to quickly calculate the elements in the auxiliary observation covariance matrix G, and then use the 2D levinson-Durbin method to calculate the LC decomposition factors of G -1 , where Q is the observation covariance matrix, H is the 2D dictionary matrix, Λ is the signal autocorrelation matrix, S is the permutation matrix composed of 0 and 1 elements, I is the identity matrix, the superscript (j) represents the j-th iteration, (·) H represents the conjugate transpose operation of the matrix, (·) T represents the transpose operation of the matrix, (·) -1 represents the matrix inversion operation.

5. A multi-resident observation 2D super-resolution ISAR imaging method based on fast IWF according to claim 4, characterized in that, Calculate the auxiliary observation covariance matrices G and G in S3 -1 and the LC decomposition factors thereof, specifically including the following steps: The observation model of Wiener filtering is: y = z + n (9) where \(y = [y(0) y(1) \ldots y(M - 1)]\) T is the observed data, \(y(m)\) represents an element in \(y\), \(m\) is the index of the element in \(y\), \(z = [z(0) z(1) \ldots z(M - 1)]\) T is the desired signal, \(z(m)\) represents an element in \(z\), \(m\) is the index of the element in \(z\), \(n\) is the observation noise, and the optimal weights of Wiener filtering and the estimate of the desired signal are respectively: where R yy is the autocorrelation matrix of the filter input signal, and R yy = E[yy H = E[(z + n)(z + n) H = R zz + R nn , where R zz and R nn are the autocorrelation matrix of the desired signal and the autocorrelation matrix of the noise respectively, r yz is the cross-correlation matrix between the input signal and the desired signal, and E[·] represents the expectation operation; Substitute the 2D sparse signal reconstruction model represented by formula (7) into the observation model of the Wiener filter. Assuming that the range-variant phase error has been compensated, let the echo after the range-variant phase error compensation be where (·) * represents the conjugate operation of the matrix. Then the forms of the echo model, the autocorrelation matrix of the observation signal, the cross-correlation matrix, and the desired signal are respectively: y g = Hx + n (12) R yy = R zz + R nn = HΛH H + δ 2 I (13) Among them, H m represents the m-th row element of H, obtaining R nn = δ 2 I, where δ 2 is the noise variance. The noise variance of the echo signal in radar imaging is estimated with the aid of auxiliary data. Λ = E[xx H is the autocorrelation matrix of the reconstructed signal. The autocorrelation matrix of the signal is constructed using the signal power value. That is, Λ is a diagonal matrix with diagonal elements γ k = |x k | 2 , where x k represents the k-th element in x, and γ k represents the precision of the signal x k ; Because z(m) = H m x, the signal estimation formula is as follows: where the observation covariance matrix Q is: Q = (H Λ H H + δ 2 I) (17) where I represents an identity matrix; Using the permutation matrix S composed of 1 and 0 elements, Q is written as: Q = S T GS (26) where G is the auxiliary observation covariance matrix, which is a matrix with a fourth-order Toeplitz tensor expansion form: Among them, the sub-matrix within G is in the form of: where the superscript <·> represents the elements in Q. Since G has a fourth-order Toeplitz tensor matrix structure, G is written in two different forms as follows: where: is an inverse matrix, where the elements on the secondary diagonal are 1 and the other elements are 0, that is Use the block matrix inversion formula for formula (29) and formula (30) respectively to obtain the inverse matrix of G as: where 0 represents a zero matrix of appropriate dimension, and Define a lower triangular Toeplitz matrix of dimension M gp ×M gp and a circulant matrix and a circulant matrix Using formulas (35) and (36), the displacement matrix of G -1 is as follows: Among them, and respectively represent a lower triangular Toeplitz-block permutation matrix composed of 1 and 0 elements and a cyclic-block permutation matrix composed of 1 and 0 elements, and their forms are as follows: Let: where t l (l = 0, 1, …, Nq - 1) is the l-th column vector in T, p l is the l-th column vector in P, is the l-th column vector in; Based on ▽G -1 , G -1 is written as: where T, P, and are the LC decomposition factors of G -1 and are calculated by the 2D Levinson-Durbin method. T, P, and are respectively defined as a TB matrix and a CB matrix:

6. A method for multi-dwelling observation 2D super-resolution ISAR imaging based on fast IWF according to claim 5, characterized in that, S4 specifically includes: Obtain the echo signal after phase error compensation where y g represents the echo signal after phase error compensation, is a matrix composed of the Fourier dictionary matrix, represents the range-variant phase error, s g represents the echo signal, (·) * represents the conjugate operation of the matrix; Use the Tsallis entropy of the image obtained by the fast IWF method as the cost function to estimate the target rotation speed as: wherein, is the estimated value of the rotational speed ω, T(X) is the Tsallis entropy of the image, and the form of T(X) is: where the total image energy is p is the exponent of the Tsallis entropy, representing the elements within the 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: where the first-order partial derivative of the range-variant phase error with respect to the rotation speed is: Re{·} represents the operation of taking the real part and solving the non - linear equation Obtain the estimated rotational speed of the target.

7. A multi-dwelling observation 2D super-resolution ISAR imaging method based on fast IWF according to claim 6, characterized in that, S5 specifically includes: Based on G -1 of the LC decomposition factor and the echo after phase error compensation, according to Use 2D-FFT to quickly calculate the estimated value of the sparse signal and according to γ k = |x k | 2 Update the precision vector γ of the sparse signal (j+1) , where k = 0, 1, …, K - 1 is the index value of the signal precision vector; As can be seen from formula (16), The calculation of is divided into three steps: φ = Q -1 y g = S T G -1 Sy g , and Substitute formula (47) into φ. The right side of the equation is the product of some TB matrices and vectors and the product of CB matrices and vectors, which are respectively transformed into linear convolution and circular convolution, and calculated quickly using 2D-FFT. The φ matrix is transformed into Φ = [φ0 φ1 … φ q-1 , and the dimension of φ i is N×M gp , and use Φ to construct a matrix where Then obtain Finally, calculate using dot product 8. A multi-resident 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 iterative error and determine whether the iteration converges based on the iterative convergence threshold η. If Er > η and j ≤ J max , then perform the next loop iteration; if Er < η, or j > J max , then stop the iteration, and the final reconstruction result is And perform azimuth calibration based on the finally estimated target rotation speed.

Images