CN104251990B - Synthetic aperture radar self-focusing method - Google Patents

Abstract

The invention discloses a kind of synthetic aperture radar self-focusing method；Specifically include following steps: distance to pulse compression, range migration correction, orientation is gone tiltedly, coordinate declines multidimensional phase error estimation and phase error, estimate single localizer unit phase error, output multidimensional phase error estimation and phase error value and compensate phase error.The synthetic aperture radar self-focusing method of the present invention uses coordinate descent that each orientation moment phase error is carried out multivariate joint probability estimation, linear search is solved the problem of maximum-contrast and is converted into the problem solving quartic polynomial by the method using multidimensional projection, substantially reduces operand；The present invention has the feature limited without exponent number to orientation to the estimation of phase error simultaneously.

Description

Synthetic aperture radar self-focusing method

Technical Field

The invention belongs to the technical field of radars, and particularly relates to a synthetic aperture radar self-focusing method.

Background

Synthetic Aperture Radar (SAR) is a full-time, all-weather, high-resolution microwave remote sensing imaging Radar. In the fields of military reconnaissance, topographic mapping, marine and hydrological observation, environmental and disaster monitoring, resource exploration, crustal micro-variation detection and the like, the SAR plays an increasingly important role. In SAR imaging under ideal conditions, the platform is assumed to move linearly at a constant speed, however, in practice, the platform often deviates from the actual motion track due to the influence of wind and turbulence, so that phase errors are introduced, and the imaging quality is reduced. Thus, motion compensation is a key technique for SAR high resolution imaging. At present, SAR motion compensation is mainly divided into two types, one is motion compensation based on sensor information, and the other is motion compensation based on echo data, so-called auto-focusing. In recent years, with the requirement of imaging resolution being continuously improved, the requirement on motion compensation is higher and higher, the motion information provided by the sensor often cannot meet the requirement on compensation accuracy, and particularly under certain conditions, such as a small unmanned aerial vehicle, even no inertial navigation information is available, so that motion compensation based on echo data, namely self-focusing, becomes more important. Among the existing auto-focusing methods, documents "basic auto focus for Synthetic Aperture Radar, T.M. Calloway, and G.W.Dohoe, IEEE Transaction on Aero space and electronic Systems, Vol.30, No.2, pp.617-621, April 1994" propose a sub-view correlation (MD) method, which estimates Doppler modulation frequency by dividing aperture and inter-aperture correlation operation, and realizes auto-focusing. However, this method is based on a second order approximation model, and cannot compensate for the influence of a high order phase error. According to the concept of Phase gradient, a Phase gradient auto-focusing algorithm (PGA) built on an order-free limited Phase error model is provided, which can estimate each Phase error and make up for the defect that an MD algorithm can only estimate low-order Phase errors. However, the PGA algorithm has a high requirement for scene contrast in the field, and under the conditions of a uniform scene and a low signal-to-noise ratio, the focusing effect is often not good because the phase history of the salient point cannot be successfully extracted. In documents "Autofocusing of ISAR images based on elementary minimization, l.xi and j.ni, IEEE Transactions on atomic and electronic Systems, vol.35, No.4, pp.1240-1252,1999", a step-by-step approach autofocus algorithm is proposed with an image quality evaluation function as a metric, and phase errors in a high-dimensional space are searched by gradually reducing a step size, so that Autofocusing is realized, and the method is suitable for various types of scenes, but the amount of operation is large due to the involvement of high-dimensional search.

Disclosure of Invention

The purpose of the invention is: aiming at the defects in the background art, a self-focusing algorithm based on the maximum contrast of an image is researched and designed so as to solve the problem that the existing self-focusing algorithm only aims at low-order phase errors and does not have the order to limit the large operation amount of a high-dimensional space searching process in phase error estimation.

For the convenience of describing the contents of the present invention, the following terms are first explained:

the term 1: SAR self-focusing algorithm

The SAR self-focusing algorithm refers to an azimuth phase error estimation and compensation algorithm based on SAR echo data.

The term 2: coordinate descent method

The coordinate descent method is a multivariable non-gradient optimization algorithm. In each iteration of the algorithm, one-dimensional search is carried out at the current point along a coordinate direction to obtain a local minimum value of a function, and different coordinate directions are circularly used in the whole process.

The term 3: multidimensional projection method

The multidimensional projection method refers to a method of projecting a vector in a multidimensional space onto a specific two-dimensional plane.

The technical scheme of the invention is as follows: a synthetic aperture radar self-focusing method comprises the following steps:

s1, acquiring two-dimensional echo data, utilizing a matched filtering method to realize range direction pulse compression, and obtaining point target echo data s after pulse compression₀(τ, η) is expressed as:

$s_0(\mathrm{\τ},\mathrm{\η})=\sin c\[\mathrm{\τ}-\frac{2R\left(\mathrm{\η}\right)}{c}\]w_{az}\left(\mathrm{\η}\right)\exp \[-j\frac{4\mathrm{\π}}{\mathrm{\λ}}R\left(\mathrm{\η}\right)\]$

where c is the speed of light, η is the azimuth slow time, τ is the distance fast time, w_az(η) is azimuth time domain envelope, lambda is wavelength, and R (η) is scene center distance history;

s2, point target echo data S in step S1₀(tau, η) performing two-dimensional Fourier transform to obtain two-dimensional frequency domain echo signal S₀(f_τ,f_η) Then two-dimensional frequency domain echo signal S₀(f_τ,f_η) Multiply-migrate corrected phaseH_RCMC(f_τ,f_η) And performing Fourier inversion on the result obtained by the multiplication to a two-dimensional time domain to obtain a two-dimensional time domain signal s after range migration correction₁(τ, η), expressed as:

$s_1(\mathrm{\τ},\mathrm{\η})=\sin c\[\mathrm{\τ}-\frac{2R_0}{c}\]w_{az}\left(\mathrm{\η}\right)\exp \[-j\frac{4\mathrm{\π}}{\mathrm{\λ}}R\left(\mathrm{\η}\right)\]$

wherein the migration corrects the phasef_cIs the carrier frequency, V is the platform motion velocity, R₀Instantaneous slope distance, f, at the moment of Doppler center crossing_τIs a distance to a reference frequency, f_ηIs the azimuth reference frequency;

s3, constructing an azimuth declivity function H by using the ideal speed and action distance information of the platform_dechirp(η), expressed as:

$H_{dechirp}\left(\mathrm{\η}\right)=\exp \left(j\frac{2{\mathrm{\πV}}^2}{{\mathrm{\λR}}_0}{\mathrm{\η}}^2\right)$

declivity function H_dechirp(η) and the distance migration corrected two-dimensional time domain signal S obtained in step S2₁(τ, η) to obtain y (τ, η) and discretely represent y (τ, η) as y_m,nWherein M is 1,2, 3., M, N is 1,2, 3., N, M is the number of sampling points in the direction of distance, and N is the number of sampling points in the direction of direction; the phase error of the echo data in the azimuth direction is phi ═ phi₁,φ₂,φ₃,…,φ_N) Wherein phi_n(N ═ 1,2, 3.., N) denotes the phase error of the nth azimuth unit in the echo data, then y will be_m,nExpressed as:

$y_{m,n}={\tilde{y}}_{m,n}\·e^{{\mathrm{j\φ}}_n}$

wherein,the method is a result of discrete representation after multiplication of a distance migration corrected two-dimensional time domain signal and an azimuth deskew function under an ideal condition; will y_m,nPerforming orientation discrete Fourier transform to obtain a two-dimensional image z without self-focusing processing_m,nExpressed as:

$z_{m,n}=1/N\underset{k=1}{\overset{N}{\Σ}}\exp (j2\mathrm{\π}kn/N){\tilde{y}}_{m,k}{e}^{{\mathrm{j\φ}}_k}$

wherein k is 1,2.. N is the number of azimuth units; for two-dimensional image z without self-focusing processing_m,nCalculating the image contrast C thereof₀Expressed as:

$C_0=\underset{q=1}{\overset{N\×M}{\Σ}}{(z_q\×z_q^{*})}^2$

wherein z is_qAs a two-dimensional image z_m,nQ is 1,2, and N × M is the number of matrix elements, [ · q]^*Is a conjugate operation;

s4, setting a phase error estimation vectorWhen i is 1, is fixedEstimate phi₁When i is>At 1 time, fixEstimate phi_iWherein i is the number of azimuth units;

s5, when the number i of azimuth units of the phase error is equal to 1, utilizing the phase error phi to be estimated_iAnd the phase error estimates of the remaining azimuth elementsConstructing a phase error compensation vectorWhen i ≠ 1, the phase error phi to be estimated is utilized_iAnd the phase error estimates of the remaining azimuth elementsConstructing a phase error compensation vectorWill y in step S3_m,nMultiplying the phase error compensation vector, and then performing azimuth Fourier transform on the result obtained after multiplication to obtain

$z_{m,n}^{\′}=\frac{1}{N}\exp (j2\mathrm{\π}n/N)y_{m,i}{e}^{-{\mathrm{j\φ}}_i}+\underset{k\≠i}{\Σ}\frac{1}{N}\exp (j2\mathrm{\π}kn/N)y_{m,k}{e}^{-j{\widehat{\mathrm{\φ}}}_k}$

To z 'obtained'_m,nAccording to the image contrast C in step S3₀The expression (c) leads to an analytical expression of the contrast ratio, and is simplified to obtain:

C(φ_i)＝||F₀+A cosφ_i+B sinφ_i||₂

wherein | · | purple sweet₂Is a vector two norm, F₀Is a constant vector quantity, satisfies

F₀＝[(f₀)_1,1,(f₀)_1,2,...,(f₀)_1,N,(f₀)_2,1,(f₀)_2,2,...,(f₀)_2,N,...,(f₀)_M,1,(f₀)_M,2,...,(f₀)_M,N]

${\left(f_0\right)}_{m,n}={\left|\underset{k\≠i}{\Σ}\frac{1}{N}\exp (j2\mathrm{\π}kn/N)y_{m,k}{e}^{-j{\widehat{\mathrm{\φ}}}_k}\right|}^2+{\left|\frac{1}{N}\exp (j2\mathrm{\π}n/N)y_{m,i}\right|}^2$

A. B is an elliptical parameter vector satisfying

A＝[(a)_1,1,(a)_1,2,...,(a)_1,N,(a)_2,1,(a)_2,2,...,(a)_2,N,...,(a)_M,1,(a)_M,2,...,(a)_M,N]

B＝[(b)_1,1,(b)_1,2,...,(b)_1,N,(b)_2,1,(b)_2,2,...,(b)_2,N,...,(b)_M,1,(b)_M,2,...,(b)_M,N]

$a_{m,n}=2\mathrm{Re}\{\frac{1}{N}\exp (-j2\mathrm{\π}n/N)y_{m,i}^{*}\[\underset{k\≠i}{\Σ}\frac{1}{N}\exp (j2\mathrm{\π}kn/N)y_{m,k}{e}^{-j{\widehat{\mathrm{\φ}}}_k}\]\}$

$b_{m,n}=-2\mathrm{Im}\{\frac{1}{N}\exp (-j2\mathrm{\π}n/N)y_{m,i}^{*}\[\underset{k\≠i}{\Σ}\frac{1}{N}\exp (j2\mathrm{\π}kn/N)y_{m,k}{e}^{-j{\widehat{\mathrm{\φ}}}_k}\]\}$

Vector F is projected by multi-dimensional projection method₀Projecting the vector A, B to form a two-dimensional plane, and solving the ellipse on the two-dimensional plane to x_oIs the farthest point from, where x_oIs a projected point with the origin of coordinates in a two-dimensional plane spanned by A, B and having coordinates of x₀＝-[E₁E₂]^T·F₀，E₁、E₂Is the unit orthogonal basis of vector A, B [ ·]^TPerforming matrix transposition operation; using the farthest point to x_oIs parallel to the geometric relationship of the ellipse normal vector, solving a fourth order polynomial about the parameter to be solved α, expressed as:

$\underset{p=0}{\overset{4}{\Σ}}{c}_p{\mathrm{\α}}^p=0$

where p is 0,1,2,3,4, which is a polynomial order, and the parameter α is related to the phase error phi to be solved_iThe relationship is as follows:

$\left[\begin{array}{c}cos{\mathrm{\φ}}_i\\ {\mathrm{sin\φ}}_i\end{array}\right]={\[AB\]}^{-1}{({\mathrm{\α}}_{min}R+I)}^{-1}{x}_0$

$\left\{\begin{array}{c}{c}_0={\mathrm{\λ}}_1{\mathrm{\β}}_1^2+{\mathrm{\λ}}_2{\mathrm{\β}}_2^2-1\\ c_1=2{\mathrm{\λ}}_1({\mathrm{\λ}}_2{\mathrm{\β}}_2^2-1)+2{\mathrm{\λ}}_2({\mathrm{\λ}}_1{\mathrm{\β}}_1^2-1)\\ c_2=({\mathrm{\λ}}_1{\mathrm{\β}}_1^2-1){\mathrm{\λ}}_2^2+({\mathrm{\λ}}_2{\mathrm{\β}}_2^2-1){\mathrm{\λ}}_1^2-4{\mathrm{\λ}}_1{\mathrm{\λ}}_2\\ c_3=-2{\mathrm{\λ}}_1{\mathrm{\λ}}_2({\mathrm{\λ}}_1+{\mathrm{\λ}}_2)\\ c_4=-{\left({\mathrm{\λ}}_1{\mathrm{\λ}}_2\right)}^2\end{array}\right.$

λ₁、λ₂is the eigenvalue, v, of an elliptic parameter matrix R₁、v₂For its corresponding feature vector, [ β ]₁β₂]^T＝[v₁v₂]^Tx₀Solving a quartic polynomial about the parameter to be solved α to obtain a minimum real number solution α of α_minα will be_minCarry-in parameter α and phase error phi to be solved_iIn the relation, the estimated value of the phase error is obtained

S6, updating the phase error estimation vector in the step S4The phase error estimated in the judgment step S5If the number i of azimuth cells in (b) is equal to N, if not, let i be i +1, and proceed to step S5 to estimate phi_i(ii) a If true, then y in step S3_m,nMultiplying by a phase error compensation vectorAnd performing azimuth Fourier transform to obtain a two-dimensional focusing image with compensated phase errorExpressed as:

${\tilde{z}}_{m,n}=1/N\underset{k=1}{\overset{N}{\Σ}}\exp (j2\mathrm{\π}kn/N)y_{m,k}{e}^{-j{\widehat{\mathrm{\φ}}}_k}$

according to the image contrast C in step S3₀Is calculated from the expressionContrast ratio C of₁Judging whether the image contrast at the moment meets the conditionIf yes, directly outputting the phase error estimation vectorIf not, updating the initial image contrast C₀Contrast C of the image at this time₁Is given to C₀Returning to step S4, repeating steps S5 and S6, and estimating phi₁,φ₂,φ₃,…,φ_NUntil the condition is satisfied, outputting the phase error estimation vector at the moment

S7. the phase error estimation vector output in the step S6Constructing a phase error compensation vectorWill y in step S3_m,nMultiplying the phase error compensation vector and performing azimuth Fourier transform to obtain a two-dimensional focusing image after phase error compensationExpressed as:

${\tilde{z}}_{m,n}=1/N\underset{k=1}{\overset{N}{\Σ}}\exp (j2\mathrm{\π}kn/N)y_{m,k}{e}^{-j{\widehat{\mathrm{\φ}}}_k}.$

further, the ellipse parameter matrix in step S5 is specifically defined as:

$R=\left[\begin{array}{cc}{r}_1& r_3\\ r_3& r_2\end{array}\right]$

$\left\{\begin{array}{c}{r}_1=(A_2^2+B_2^2)/(A_2{B}_1-A_1{B}_2)\\ r_2=(A_1^2+B_1^2)/(A_2{B}_1-A_1{B}_2)\\ r_3=-(A_1{A}_2+B_1{B}_2)/(A_2{B}_1-A_1{B}_2)\end{array}\right.$

[A₁A₂]＝A^T[E₁E₂]^T

[B₁B₂]＝B^T[E₁E₂]^T。

the invention has the beneficial effects that: the synthetic aperture radar self-focusing method carries out multi-dimensional joint estimation on phase errors at each azimuth moment by using a coordinate descent method, converts the problem of solving the maximum contrast by one-dimensional search into the problem of solving a quartic polynomial by using a multi-dimensional projection method, and greatly reduces the operation amount; meanwhile, the method has the characteristic of no order limitation on the estimation of the azimuth phase error.

Drawings

Fig. 1 is a schematic flow chart of a synthetic aperture radar self-focusing method of the present invention.

Fig. 2 is a schematic diagram of the bistatic forward-looking SAR geometry of the present invention.

FIG. 3 is a target scenario layout diagram of the present invention.

Fig. 4 is a schematic representation of the results of the pre-autofocus imaging of the invention.

FIG. 5 is a schematic of the azimuthal phase error of the present invention.

FIG. 6 is a graphical representation of the imaging results after autofocusing according to the invention.

Detailed Description

In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention is described in further detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.

The invention mainly adopts a simulation experiment method for verification, and all the steps and conclusions are verified to be correct on Matlab 2012. Fig. 1 is a schematic flow chart of a synthetic aperture radar self-focusing method according to the present invention. A synthetic aperture radar self-focusing method comprises the following steps:

s1, acquiring two-dimensional echo data, utilizing a matched filtering method to realize range direction pulse compression, and obtaining point target echo data s after pulse compression₀(τ, η) is expressed as:

$s_0(\mathrm{\τ},\mathrm{\η})=\sin c\[\mathrm{\τ}-\frac{2R\left(\mathrm{\η}\right)}{c}\]w_{az}\left(\mathrm{\η}\right)\exp \[-j\frac{4\mathrm{\π}}{\mathrm{\λ}}R\left(\mathrm{\η}\right)\]$

where c is the speed of light, η is the azimuth slow time, τ is the distance fast time, w_az(η) is the azimuth time domain envelope, λ is the wavelength, and R (η) is the scene center distance history.

As shown in fig. 2, is a schematic diagram of the bistatic forward-looking SAR geometry of the present invention. The parameters required for the simulation of the present invention are shown in the following table.

And calculating radar distance history of each point target in the central area of the imaging area to generate an SAR point target simulation echo matrix. Fig. 3 is a diagram showing the arrangement of the target scene according to the present invention. The distance direction Fourier transformation is carried out on the echo, and then the distance direction is multiplied by a matched filter matching function H₁(f_τ) And realizing the distance direction pulse compression.

Wherein,distance direction frequency f_τThe variation range is [ -24,24 [)]MHz, range of variation of distance time tau [6.13 × 10%^-5,7.20×10^-5]And second.

S2, point target echo data S in step S1₀(tau, η) performing two-dimensional Fourier transform to obtain two-dimensional frequency domain echo signal S₀(f_τ,f_η) Then two-dimensional frequency domain echo signal S₀(f_τ,f_η) Riding deviceCorrecting phase H by migration_RCMC(f_τ,f_η) And performing Fourier inversion on the result obtained by the multiplication to a two-dimensional time domain to obtain a two-dimensional time domain signal s after range migration correction₁(τ, η), expressed as:

$s_1(\mathrm{\τ},\mathrm{\η})=\sin c\[\mathrm{\τ}-\frac{2R_0}{c}\]w_{az}\left(\mathrm{\η}\right)\exp \[-j\frac{4\mathrm{\π}}{\mathrm{\λ}}R\left(\mathrm{\η}\right)\]$

wherein the migration corrects the phasef_cIs the carrier frequency, V is the platform motion velocity, R₀Instantaneous slope distance, f, at the moment of Doppler center crossing_τIs a distance to a reference frequency, f_ηFor the azimuth reference frequency, azimuth frequency f_ηThe variation range is [ -300,300 [ -300]Hz, azimuth time η range of [ -0.425,0.425]Second, η is 0 seconds when the beam center irradiates the target.

S3, constructing an azimuth declivity function H by using the ideal speed and action distance information of the platform_dechirp(η), expressed as:

$H_{dechirp}\left(\mathrm{\η}\right)=\exp \left(j\frac{2{\mathrm{\πV}}^2}{{\mathrm{\λR}}_0}{\mathrm{\η}}^2\right)$

declivity function H_dechirp(η) and the distance migration corrected two-dimensional time domain signal S obtained in step S2₁(τ, η) to obtain y (τ, η) and discretely represent y (τ, η) as y_m,nWherein M is 1,2, 3., M, N is 1,2, 3., N, M is the number of sampling points in the direction of distance, and N is the number of sampling points in the direction of direction; the phase error of the echo data in the azimuth direction is phi ═ phi₁,φ₂,φ₃,…,φ_N) Wherein phi_n(N ═ 1,2, 3.., N) denotes the phase error of the nth azimuth unit in the echo data, then y will be_m,nExpressed as:

$y_{m,n}={\tilde{y}}_{m,n}\·e^{{\mathrm{j\φ}}_n}$

wherein,the method is a result of discrete representation after multiplication of a distance migration corrected two-dimensional time domain signal and an azimuth deskew function under an ideal condition; will y_m,nPerforming orientation discrete Fourier transform to obtain a two-dimensional image z without self-focusing processing_m,nExpressed as:

$z_{m,n}=1/N\underset{k=1}{\overset{N}{\Σ}}\exp (j2\mathrm{\π}kn/N){\tilde{y}}_{m,k}{e}^{{\mathrm{j\φ}}_k}$

wherein k is 1,2.. N is the number of azimuth units; for two-dimensional image z without self-focusing processing_m,nCalculating the image contrast C thereof₀Expressed as:

$C_0=\underset{q=1}{\overset{N\×M}{\Σ}}{(z_q\×z_q^{*})}^2$

wherein z is_qAs a two-dimensional image z_m,nQ is 1,2, and N × M is the number of matrix elements, [ · q]^*Is a conjugate operation.

Here, we set the azimuth phase error in the echo signal to Φ (η) ═ 6 η²+3η³-2η⁴+η⁵Discretizing phi (η) into an N-dimensional phase error vector phi, adding the phase error vector phi to the echo signalObtaining the azimuth deskew function and s in the actual situation₁(τ, η) multiplied result y_m,k. For signal y_m,kCarrying out azimuth Fourier transform to obtain a pre-autofocus imaging result z_m,nAnd calculating the image contrast of the SAR imaging result before the self-focusing processing. Fig. 4 is a diagram illustrating the result of the imaging before autofocusing according to the present invention.The contrast of the image before the autofocus processing is C₀＝1.19×10¹⁸。

S4, setting a phase error estimation vectorWhen i is 1, is fixedEstimate phi₁When i is>At 1 time, fixEstimate phi_iAnd i is the number of azimuth units.

Due to estimating a multi-dimensional phase error phi₁,φ₂,φ₃,…,φ_NThe method adopts a coordinate descent method, namely, the phase errors of other azimuth units are kept fixed in each iteration, and one-dimensional search is carried out along the coordinate direction of the phase error of a specific azimuth unit.

S5, when the number i of azimuth units of the phase error is equal to 1, utilizing the phase error phi to be estimated_iAnd the phase error estimates of the remaining azimuth elementsConstructing a phase error compensation vectorWhen i ≠ 1, the phase error phi to be estimated is utilized_iAnd the phase error estimates of the remaining azimuth elementsConstructing a phase error compensation vectorWill y in step S3_m,nMultiplying the phase error compensation vector, and then performing azimuth Fourier transform on the result obtained after multiplication to obtain

$z_{m,n}^{\′}=\frac{1}{N}\exp (j2\mathrm{\π}n/N)y_{m,i}{e}^{-{\mathrm{j\φ}}_i}+\underset{k\≠i}{\Σ}\frac{1}{N}\exp (j2\mathrm{\π}kn/N)y_{m,k}{e}^{-j{\widehat{\mathrm{\φ}}}_k}$

To z 'obtained'_m,nAccording to the image contrast C in step S3₀The expression (c) leads to an analytical expression of the contrast ratio, and is simplified to obtain:

C(φ_i)＝||F₀+A cosφ_i+B sinφ_i||₂

wherein, C (phi)_i) Is a vector F₀+Acosφ_i+Bsinφ_iIs the length from the origin of coordinates in the N × M dimensional space to the ellipse Acos phi_i+Bsinφ_iDistance of (1) | · | | non-calculation₂Is a vector two norm, F₀Is a constant vector quantity, satisfies

F₀＝[(f₀)_1,1,(f₀)_1,2,...,(f₀)_1,N,(f₀)_2,1,(f₀)_2,2,...,(f₀)_2,N,...,(f₀)_M,1,(f₀)_M,2,...,(f₀)_M,N]

${\left(f_0\right)}_{m,n}={\left|\underset{k\≠i}{\Σ}\frac{1}{N}\exp (j2\mathrm{\π}kn/N)y_{m,k}{e}^{-j{\widehat{\mathrm{\φ}}}_k}\right|}^2+{\left|\frac{1}{N}\exp (j2\mathrm{\π}n/N)y_{m,i}\right|}^2$

A. B is an elliptical parameter vector satisfying

A＝[(a)_1,1,(a)_1,2,...,(a)_1,N,(a)_2,1,(a)_2,2,...,(a)_2,N,...,(a)_M,1,(a)_M,2,...,(a)_M,N]

B＝[(b)_1,1,(b)_1,2,...,(b)_1,N,(b)_2,1,(b)_2,2,...,(b)_2,N,...,(b)_M,1,(b)_M,2,...,(b)_M,N]

$a_{m,n}=2\mathrm{Re}\{\frac{1}{N}\exp (-j2\mathrm{\π}n/N)y_{m,i}^{*}\[\underset{k\≠i}{\Σ}\frac{1}{N}\exp (j2\mathrm{\π}kn/N)y_{m,k}{e}^{-j{\widehat{\mathrm{\φ}}}_k}\]\}$

$b_{m,n}=-2\mathrm{Im}\{\frac{1}{N}\exp (-j2\mathrm{\π}n/N)y_{m,i}^{*}\[\underset{k\≠i}{\Σ}\frac{1}{N}\exp (j2\mathrm{\π}kn/N)y_{m,k}{e}^{-j{\widehat{\mathrm{\φ}}}_k}\]\}$

The ellipse parameter matrix is defined as:

$R=\left[\begin{array}{cc}{r}_1& r_3\\ r_3& r_2\end{array}\right]$

$\left\{\begin{array}{c}{r}_1=(A_2^2+B_2^2)/(A_2{B}_1-A_1{B}_2)\\ r_2=(A_1^2+B_1^2)/(A_2{B}_1-A_1{B}_2)\\ r_3=-(A_1{A}_2+B_1{B}_2)/(A_2{B}_1-A_1{B}_2)\end{array}\right.$

[A₁A₂]＝A^T[E₁E₂]^T

[B₁B₂]＝B^T[E₁E₂]^T

vector F is projected by multi-dimensional projection method₀Projected onto a two-dimensional plane spanned by vector A, B, the phase error φ corresponding to the maximum image contrast will be resolved₁Conversion to solving an ellipse to x on a two-dimensional plane_oIs the farthest point from, where x_oIs a projected point with the origin of coordinates in a two-dimensional plane spanned by A, B and having coordinates of x₀＝-[E₁E₂]^T·F₀，E₁、E₂Is the unit orthogonal basis of vector A, B [ ·]^TPerforming matrix transposition operation; using the farthest point to x_oThe vector of (2) is parallel to the geometric relation of the normal vector of the ellipse, the maximum distance from a point outside the ellipse to be solved is converted into a fourth-order polynomial about the parameter to be solved α, and the polynomial is expressed as:

$\underset{p=0}{\overset{4}{\Σ}}{c}_p{\mathrm{\α}}^p=0$

where p is 0,1,2,3,4, which is a polynomial order, and the parameter α is related to the phase error phi to be solved_iThe relationship is as follows:

$\left[\begin{array}{c}cos{\mathrm{\φ}}_i\\ {\mathrm{sin\φ}}_i\end{array}\right]={\[AB\]}^{-1}{({\mathrm{\α}}_{min}R+I)}^{-1}{x}_0$

$\left\{\begin{array}{c}{c}_0={\mathrm{\λ}}_1{\mathrm{\β}}_1^2+{\mathrm{\λ}}_2{\mathrm{\β}}_2^2-1\\ c_1=2{\mathrm{\λ}}_1({\mathrm{\λ}}_2{\mathrm{\β}}_2^2-1)+2{\mathrm{\λ}}_2({\mathrm{\λ}}_1{\mathrm{\β}}_1^2-1)\\ c_2=({\mathrm{\λ}}_1{\mathrm{\β}}_1^2-1){\mathrm{\λ}}_2^2+({\mathrm{\λ}}_2{\mathrm{\β}}_2^2-1){\mathrm{\λ}}_1^2-4{\mathrm{\λ}}_1{\mathrm{\λ}}_2\\ c_3=-2{\mathrm{\λ}}_1{\mathrm{\λ}}_2({\mathrm{\λ}}_1+{\mathrm{\λ}}_2)\\ c_4=-{\left({\mathrm{\λ}}_1{\mathrm{\λ}}_2\right)}^2\end{array}\right.$

λ₁、λ₂is the eigenvalue, v, of an elliptic parameter matrix R₁、v₂For its corresponding feature vector, [ β ]₁β₂]^T＝[v₁v₂]^Tx₀Solving a quartic polynomial about the parameter to be solved α to obtain a minimum real number solution α of α_minα will be_minCarry-in parameter α and phase error phi to be solved_iIn the relation, the estimated value of the phase error is obtained

S6, updating the phase error estimation vector in the step S4The phase error estimated in the judgment step S5If the number i of azimuth cells in (b) is equal to N, if not, let i be i +1, and proceed to step S5 to estimate phi_i(ii) a If true, then y in step S3_m,nMultiplying by a phase error compensation vectorAnd performing azimuth Fourier transform to obtain a two-dimensional focusing image with compensated phase errorExpressed as:

${\tilde{z}}_{m,n}=1/N\underset{k=1}{\overset{N}{\Σ}}\exp (j2\mathrm{\π}kn/N)y_{m,k}{e}^{-j{\widehat{\mathrm{\φ}}}_k}$

according to the image contrast C in step S3₀Is calculated from the expressionContrast ratio C of₁Judging whether the image contrast at the moment meets the conditionIf yes, directly outputting the phase error estimation vectorIf not, updating the initial image contrast C₀Contrast C of the image at this time₁Is given to C₀Returning to step S4, repeating steps S5 and S6, and estimating phi₁,φ₂,φ₃,…,φ_NUntil the condition is satisfied, outputting the phase error estimation vector at the momentParameter D here₀＝10^-3After three cycles, the inequality in equation (25) holds, at which time the image contrast is C₁＝7.37×10¹⁸. Actual value phi and estimated value of error phaseAs shown in fig. 5.

S7, the method comprises the following stepsPhase error estimate vector output in S6Constructing a phase error compensation vectorWill y in step S3_m,nMultiplying the phase error compensation vector and performing azimuth Fourier transform to obtain a two-dimensional focusing image after phase error compensationExpressed as:

${\tilde{z}}_{m,n}=1/N\underset{k=1}{\overset{N}{\Σ}}\exp (j2\mathrm{\π}kn/N)y_{m,k}{e}^{-j{\widehat{\mathrm{\φ}}}_k}.$

fig. 6 is a schematic diagram showing the imaging result after self-focusing according to the present invention. The invention estimates the azimuth phase error without order limitation and has accurate estimation, and the estimated mean square error is only 0.061. Under the same conditions, the method requires 61.42 seconds for the estimation of the azimuth phase error, and 487.89 seconds for processing the specific embodiment of the invention by adopting a successive approximation self-focusing method; compared with a progressive approach self-focusing method, the method has the advantage that the efficiency is improved by 7.9 times.

It will be appreciated by those of ordinary skill in the art that the embodiments described herein are intended to assist the reader in understanding the principles of the invention and are to be construed as being without limitation to such specifically recited embodiments and examples. Those skilled in the art can make various other specific changes and combinations based on the teachings of the present invention without departing from the spirit of the invention, and these changes and combinations are within the scope of the invention.

Claims (2)

1. A synthetic aperture radar self-focusing method, comprising the steps of:

s1, acquiring two-dimensional echo data, utilizing a matched filtering method to realize range direction pulse compression, and obtaining point target echo data s after pulse compression₀(τ, η) is expressed as:

$s_0(\mathrm{\τ},\mathrm{\η})=\sin c\[\mathrm{\τ}-\frac{2R\left(\mathrm{\η}\right)}{c}\]w_{az}\left(\mathrm{\η}\right)\exp \[-j\frac{4\mathrm{\π}}{\mathrm{\λ}}R\left(\mathrm{\η}\right)\]$

where c is the speed of light, η is the azimuth slow time, τ is the distance fast time, w_az(η) is azimuth time domain envelope, lambda is wavelength, and R (η) is scene center distance history;

s2, point target echo data S in step S1₀(tau, η) performing two-dimensional Fourier transform to obtain two-dimensional frequency domain echo signal S₀(f_τ,f_η) Then two-dimensional frequency domain echo signal S₀(f_τ,f_η) Multiplied by a migration-corrected phase H_RCMC(f_τ,f_η) And performing Fourier inversion on the result obtained by the multiplication to a two-dimensional time domain to obtain a two-dimensional time domain signal s after range migration correction₁(τ, η), expressed as:

$s_1(\mathrm{\τ},\mathrm{\η})=\sin c\[\mathrm{\τ}-\frac{2R_0}{c}\]w_{az}\left(\mathrm{\η}\right)\exp \[-j\frac{4\mathrm{\π}}{\mathrm{\λ}}R\left(\mathrm{\η}\right)\]$

wherein the migration corrects the phasef_cIs the carrier frequency, V is the platform motion velocity, R₀Instantaneous slope distance, f, at the moment of Doppler center crossing_τIs a distance to a reference frequency, f_ηIs the azimuth reference frequency;

s3, constructing an azimuth declivity function H by using the ideal speed and action distance information of the platform_dechirp(η), expressed as:

$H_{dechirp}\left(\mathrm{\η}\right)=\exp \left(j\frac{2{\mathrm{\πV}}^2}{{\mathrm{\λR}}_0}{\mathrm{\η}}^2\right)$

declivity function H_dechirp(η) and the distance migration corrected two-dimensional time domain signal S obtained in step S2₁(τ, η) to obtain y (τ, η) and discretely represent y (τ, η) as y_m,nWherein M is 1,2, 3., M, N is 1,2, 3., N, M is the number of sampling points in the direction of distance, and N is the number of sampling points in the direction of direction; the phase error of the echo data in the azimuth direction is phi ═ phi₁,φ₂,φ₃,…,φ_N) Wherein phi_n(N ═ 1,2, 3.., N) denotes the phase error of the nth azimuth unit in the echo data, then y will be_m,nExpressed as:

$y_{m,n}={\tilde{y}}_{m,n}\·e^{{\mathrm{j\φ}}_n}$

wherein,the method is a result of discrete representation after multiplication of a distance migration corrected two-dimensional time domain signal and an azimuth deskew function under an ideal condition; will y_m,nPerforming orientation discrete Fourier transform to obtain a two-dimensional image z without self-focusing processing_m,nExpressed as:

$z_{m,n}=1/N\underset{k=1}{\overset{N}{\Σ}}\exp (j2\mathrm{\π}kn/N){\tilde{y}}_{m,k}{e}^{{\mathrm{j\φ}}_k}$

wherein k is 1,2.. N is the number of azimuth units; for two-dimensional image z without self-focusing processing_m,nCalculating the image contrast C thereof₀Expressed as:

$C_0=\underset{q=1}{\overset{N\×M}{\Σ}}{(z_q\×z_q^{*})}^2$

wherein z is_qAs a two-dimensional image z_m,nQ is 1,2, and N × M is the number of matrix elements, [ · q]^*Is a conjugate operation;

s4, setting a phase error estimation vectorWhen i is 1, is fixedEstimate phi₁When i is>At 1 time, fixEstimate phi_iWherein i is the number of azimuth units;

s5, when the number i of azimuth units of the phase error is equal to 1, utilizing the phase error phi to be estimated_iAnd the phase error estimates of the remaining azimuth elementsConstructing a phase error compensation vectorWhen i ≠ 1, the phase error phi to be estimated is utilized_iAnd the phase error estimates of the remaining azimuth elementsConstructing a phase error compensation vectorWill y in step S3_m,nMultiplying the phase error compensation vector, and then performing azimuth Fourier transform on the result obtained after multiplication to obtain

$z_{m,n}^{\′}=\frac{1}{N}\exp (j2\mathrm{\π}n/N)y_{m,i}{e}^{-{\mathrm{j\φ}}_i}+\underset{k\≠i}{\Σ}\frac{1}{N}\exp (j2\mathrm{\π}kn/N)y_{m,k}{e}^{-j{\widehat{\mathrm{\φ}}}_k}$

To z 'obtained'_m,nAccording to the image contrast C in step S3₀The expression (c) leads to an analytical expression of the contrast ratio, and is simplified to obtain:

C(φ_i)＝||F₀+Acosφ_i+Bsinφ_i||₂

wherein | · | purple sweet₂Is a vector two norm, F₀Is a constant vector quantity, satisfies

F₀＝[(f₀)_1,1,(f₀)_1,2,...,(f₀)_1,N,(f₀)_2,1,(f₀)_2,2,...,(f₀)_2,N,...,(f₀)_M,1,(f₀)_M,2,...,(f₀)_M,N]

${\left(f_0\right)}_{m,n}=|\underset{k\≠i}{\Σ}\frac{1}{N}\exp (j2\mathrm{\π}kn/N)y_{m,k}{e}^{-j{\widehat{\mathrm{\φ}}}_k}{|}^2+|\frac{1}{N}\exp (j2\mathrm{\π}n/N)y_{m,i}{|}^2$

A. B is an elliptical parameter vector satisfying

A＝[(a)_1,1,(a)_1,2,...,(a)_1,N,(a)_2,1,(a)_2,2,...,(a)_2,N,...,(a)_M,1,(a)_M,2,...,(a)_M,N]

B＝[(b)_1,1,(b)_1,2,...,(b)_1,N,(b)_2,1,(b)_2,2,...,(b)_2,N,...,(b)_M,1,(b)_M,2,...,(b)_M,N]

$a_{m,n}=2\mathrm{Re}\{\frac{1}{N}\exp (-j2\mathrm{\π}n/N)y_{m,i}^{*}\[\underset{k\≠i}{\Σ}\frac{1}{N}\exp (j2\mathrm{\π}kn/N)y_{m,k}{e}^{-j{\widehat{\mathrm{\φ}}}_k}\]\}$

$b_{m,n}=-2Im\{\frac{1}{N}\exp (-j2\mathrm{\π}n/N)y_{m,i}^{*}\[\underset{k\≠i}{\Σ}\frac{1}{N}\exp (j2\mathrm{\π}kn/N)y_{m,k}{e}^{-j{\widehat{\mathrm{\φ}}}_k}\]\}$

Vector F is projected by multi-dimensional projection method₀Projecting the vector A, B to form a two-dimensional plane, and solving the ellipse on the two-dimensional plane to x_oIs the farthest point from, where x_oIs a projected point with the origin of coordinates in a two-dimensional plane spanned by A, B and having coordinates of x₀＝-[E₁E₂]^T·F₀，E₁、E₂Is the unit orthogonal basis of vector A, B [ ·]^TPerforming matrix transposition operation; using the farthest point to x_oIs parallel to the geometric relationship of the ellipse normal vector, solving a fourth order polynomial about the parameter to be solved α, expressed as:

$\underset{p=0}{\overset{4}{\Σ}}{c}_p{\mathrm{\α}}^p=0$

where p is 0,1,2,3,4, which is a polynomial order, and the parameter α is related to the phase error phi to be solved_iThe relationship is as follows:

$\left[\begin{array}{c}cos{\mathrm{\φ}}_i\\ {\mathrm{sin\φ}}_i\end{array}\right]={\[AB\]}^{-1}{({\mathrm{\α}}_{min}R+I)}^{-1}{x}_0$

$\left\{\begin{array}{c}{c}_0={\mathrm{\λ}}_1{\mathrm{\β}}_1^2+{\mathrm{\λ}}_2{\mathrm{\β}}_2^2-1\\ c_1=2{\mathrm{\λ}}_1({\mathrm{\λ}}_2{\mathrm{\β}}_2^2-1)+2{\mathrm{\λ}}_2({\mathrm{\λ}}_1{\mathrm{\β}}_1^2-1)\\ c_2=({\mathrm{\λ}}_1{\mathrm{\β}}_1^2-1){\mathrm{\λ}}_2^2+({\mathrm{\λ}}_2{\mathrm{\β}}_2^2-1){\mathrm{\λ}}_1^2-4{\mathrm{\λ}}_1{\mathrm{\λ}}_2\\ c_3=-2{\mathrm{\λ}}_1{\mathrm{\λ}}_2({\mathrm{\λ}}_1+{\mathrm{\λ}}_2)\\ c_4=-{\left({\mathrm{\λ}}_1{\mathrm{\λ}}_2\right)}^2\end{array}\right.$

λ₁、λ₂is the eigenvalue, v, of an elliptic parameter matrix R₁、v₂For its corresponding feature vector, [ β ]₁β₂]^T＝[v₁v₂]^Tx₀Solving a quartic polynomial about the parameter to be solved α to obtain a minimum real number solution α of α_minα will be_minCarry-in parameter α and phase error phi to be solved_iIn the relation, the estimated value of the phase error is obtained

S6, updating the phase error estimation vector in the step S4The phase error estimated in the judgment step S5If the number i of azimuth cells in (b) is equal to N, if not, let i be i +1, and proceed to step S5 to estimate phi_i(ii) a If it is trueThen y in step S3_m,nMultiplying by a phase error compensation vectorAnd performing azimuth Fourier transform to obtain a two-dimensional focusing image with compensated phase errorExpressed as:

${\tilde{z}}_{m,n}=1/N\underset{k=1}{\overset{N}{\Σ}}\exp (j2\mathrm{\π}kn/N)y_{m,k}{e}^{-j{\widehat{\mathrm{\φ}}}_k}$

according to the image contrast C in step S3₀Is calculated from the expressionContrast ratio C of₁Judging whether the image contrast at the moment meets the conditionIf yes, directly outputting the phase error estimation vectorIf not, updating the initial image contrast C₀Contrast C of the image at this time₁Is given to C₀Returning to step S4, repeating steps S5 and S6, and estimating phi₁,φ₂,φ₃,…,φ_NUntil the condition is satisfied, outputting the phase error estimation vector at the moment

S7. the phase error estimation vector output in the step S6Constructing a phase error compensation vectorWill y in step S3_m,nMultiplying the phase error compensation vector and performing azimuth Fourier transform to obtain a two-dimensional focusing image after phase error compensationExpressed as:

${\tilde{z}}_{m,n}=1/N\underset{k=1}{\overset{N}{\Σ}}\exp (j2\mathrm{\π}kn/N)y_{m,k}{e}^{-j{\widehat{\mathrm{\φ}}}_k}.$

2. the method of self-focusing of synthetic aperture radar as claimed in claim 1, wherein the elliptical parameter matrix definition in step S5 is specifically as follows:

$R=\left[\begin{array}{cc}{r}_1& r_3\\ r_3& r_2\end{array}\right]$

$\left\{\begin{array}{c}{r}_1=(A_2^2+B_2^2)/(A_2{B}_1-A_1{B}_2)\\ r_2=(A_1^2+B_1^2)/(A_2{B}_1-A_1{B}_2)\\ r_3=-(A_1{A}_2+B_1{B}_2)/(A_2{B}_1-A_1{B}_2)\end{array}\right.$

[A₁A₂]＝A^T[E₁E₂]^T

[B₁B₂]＝B^T[E₁E₂]^T。