CN104391295A - Compressive sensing SAR sparse self-focusing imaging method with optimum image entropy - Google Patents

Abstract

The invention discloses a compressive sensing SAR sparse self-focusing imaging method with optimum image entropy. The method is targeted for the influence of azimuth phase errors in a SAR echo signal measurement model on compressive sensing SAR imaging, and the problem of estimation and compensation of unknown phase errors in a compressive sensing SAR imaging model; azimuth phase error characteristics and sparse target characteristics in the imaging model are utilized, and compressive sensing SAR image entropy is adopted as an evaluation criteria; relationship of SAR azimuth echoes and observed objects are utilized to estimate the azimuth phase errors in each iteration processing process; then, phase error compensation is carried out on the compressive sensing imaging model; and next, compressive sensing SAR imaging is carried out, and a successive iteration method is utilized to enable the image entropy of the compressive sensing SAR imaging to be optimum, thereby improving compressive sensing SAR imaging quality.

Description

A kind of compressed sensing SAR sparse self-focusing formation method of image entropy optimum

Technical field:

This technological invention belongs to Radar Technology field, and it is in particular to synthetic-aperture radar (SAR) technical field of imaging.

Background technology:

Owing to having the advantages such as round-the-clock, round-the-clock and large scene observation, synthetic-aperture radar (SAR) has become an important remote sensing technology of current large-scale terrain mapping, plays a greater and greater role in fields such as topographic mapping, Natural calamity monitoring and survey of natural resources.The sparse reconstruct of compressed sensing is as a kind of new signal processing theory proposed in recent years, breach the constraint of traditional Nyquist sampling thheorem, can utilize far below Nyquist sampling rate Accurate Reconstruction original sparse signal (referring to list of references " D.L.Donoho.Compressed sensing.IEEE Transactions on Information Theory; 2006; 52 (4): 1289-1306 "), in reduction SAR system sampling rate and raising image quality etc., have huge application potential.Therefore, compressed sensing SAR imaging has become the emerging hot subject in SAR field.Compressed sensing reconstructing method is very high to the precise requirements of signal measurement model, if signal measurement model exists error or out of true, the reconstruct degree of accuracy of compressed sensing reconstructing method will serious degradation, sometimes even there will be the result of mistake.For the imaging of compressed sensing SAR real data, the key factor affecting its imaging effect is motion error extraction and compensation (referring to list of references " Jihua Tian; Sun Jinping; Han Xiao; and Zhang Bingchen.Motion compensation for Compressive Sensing SAR imaging with autofocus.IEEE Conference on Industrial Electronics and Applications (ICIEA); 1564-1567,2011 ").Need when compressed sensing SAR imaging to utilize each orientation to set up SAR signal observation model to position of platform.But because the factors such as Platform movement exist measuring error, make each orientation be mixed into phase error in echoed signal, causing compressed sensing SAR image quality to decline even cannot imaging.Therefore in compressed sensing SAR actual imaging process, must first to each orientation of SAR data to unknown phase error carry out effectively estimating and compensating, just can carry out the imaging of high precision compressed sensing to SAR raw radar data.

Entropy, as the probabilistic yardstick of descriptor, has a wide range of applications in signal transacting field.For SAR imaging, out-of-focus image contrast is lower, shows that image each pixel probability is close, and uncertainty is comparatively large, and focusedimage contrast is higher, illustrates that the determinacy of each pixel value increases.Therefore, image entropy can be used as the key criteria weighing SAR image quality good or not.To based in the SAR data process of conventional imaging method, image entropy optimization has been used to the important means (referring to list of references " X.Li; G.S.Liu; and J.L.Ni.Autofocusing of ISAR images based on entropy minimization.IEEE Transactions on Aerospace and Electronic Systems; 1999,35 (4): 1240-1252 ") realizing self-focusing imaging.But, on signal transacting, this qualitative difference is had based on the SAR formation method of compressed sensing and traditional SAR formation method, based on compressed sensing SAR formation method can from lack sampling echo data Accurate Reconstruction observed object, thus realize high precision imaging, but traditional SAR formation method is difficult to realize the imaging of lack sampling echo data high precision.At present, the traditional auto-focus method based on image entropy optimum processes for the SAR echo data under full sampling condition, and in lack sampling situation, this class methods self-focusing performance can even lose efficacy by degradation.Therefore, in order to utilize image entropy to carry out self-focusing imaging in compressed sensing SAR imaging, directly can not utilize traditional images entropy self-focusing formation method, must process in conjunction with compressed sensing SAR signal reconstruction model and reconstruction result.

Summary of the invention:

In order to solve unknown orientation in compressed sensing SAR sparse imaging process to the estimation of phase error and compensation problem, the present invention according to orientation in compressed sensing SAR imaging model to phase error properties and sparse target signature, in conjunction with the sparse reconstructing method of compressed sensing and SAR image entropy, propose a kind of compressed sensing SAR sparse self-focusing formation method of image entropy optimum, the present invention to utilize in imaging model orientation to phase error properties and sparse target signature, utilize compressed sensing SAR image entropy as interpretational criteria, propose a kind of compressed sensing SAR sparse self-focusing formation method of image entropy optimum, the method utilizes SAR orientation to estimate that orientation is to phase error to the relation of echo and observed object in iterative process each time, then phase error compensation is carried out to compressed sensing imaging model, then compressed sensing SAR imaging is carried out again, and adopt iterative technique to make the image entropy of compressed sensing SAR imaging optimum, thus improve compressed sensing SAR image quality.

Content of the present invention for convenience of description, first make following term definition:

Definition 1, sparse signal

If the number of nonzero value is much smaller than the length of signal itself in a discrete signal, then this signal can be thought sparse.If X=is [x ₁, x ₂..., x _n] ^tfor the column vector that N number of discrete signal forms, wherein x ₁represent the 1st element in vectorial X, x ₂represent the 2nd element in vectorial X, x _nrepresent the N number of element in vectorial X, upper right corner roman symbol T is transpose operation symbol.If only have K in vectorial X ₀individual element non-zero or much larger than zero, then vectorial X is defined as K ₀sparse vector.Refer to document " S.Mallat.A Wavelet Tour of Signal Processing:The Sparse Way.Access Online via Elsevier, 2008 ".

Definition 2, norm

If X is number field linear Space, represent complex field, if it meets following character: || X||>=0, and || X||=0 only has X=0, || aX||=|a|||X||, a are arbitrary constant, || X ₁+ X ₂||≤|| X ₁||+|| X ₂||, then claiming || X|| is X norm spatially, || || represent norm sign, wherein X ₁and X ₂for X any two values spatially.Discrete signal vector X=[x is tieed up for N × 1 in definition 1 ₁, x ₂..., x _n] ^t, the LP norm expression formula of vectorial X is wherein x _ifor i-th element of vectorial X, || represent absolute value sign, Σ || represent absolute value summation symbol, the L1 norm expression formula of vectorial X is the L2 norm expression formula of vector X is the L0 norm expression formula of vector X is and x _i≠ 0.Refer to document " matrix theory ", Huang Tingzhu etc. write, and Higher Education Publishing House publishes.

Definition 3, linearly measurement model

For a digital Signal Measurement System, suppose that discrete signal vector X=[x is tieed up in N × 1 in definition 1 ₁, x ₂..., x _n] ^tfor this digital signal measuring system needs the signal of measurement, vectorial Y=[y ₁, y ₂..., y _m] ^tdiscrete signal vector is tieed up, wherein y in M × 1 exported for this digital signal measuring system ₁represent the 1st element in vectorial Y, y ₂represent the 2nd element in vectorial Y, y _mrepresent M element in vectorial Y, upper right corner T is transpose operation symbol.The linearly measurement model of this measuring system refers to that the relation of measuring-signal Y and measured signal X can be expressed as Y=AX, and wherein A is M × N matrix, and matrix A is called the calculation matrix of signal X in measuring system.

Definition 4, synthetic-aperture radar slow moment and fast moment

SAR motion platform flies over an orientation and is called the slow time to the time required for length of synthetic aperture, radar system is with the repetition period launch and accept pulse of certain hour length, therefore the slow moment can be expressed as the time discretization variable that take pulse repetition time as step-length, and wherein each pulse repetition time discrete-time variable value is a slow moment.The synthetic-aperture radar fast moment refers to that, within a pulse repetition time, distance is to the time interval variable of sampled echo signals.Refer to document " synthetic aperture radar image-forming principle ", Pi Yiming etc. write, and publishing house of University of Electronic Science and Technology publishes.

Definition 5, the compression of synthetic-aperture radar gauged distance

The compression of synthetic-aperture radar gauged distance refers to and utilizes synthetic-aperture radar transmission signal parameters, adopts the distance of matched filtering technique Technologies Against Synthetic Aperture Radar to carry out the process of filtering to signal.Refer to document " radar imaging technology ", protect polished etc. writing, Electronic Industry Press publishes.

Definition 6, synthetic-aperture radar original echo emulation mode

Standard synthetic aperture radar original echo emulation mode refers to given simulation parameter, there is under going out certain systematic parameter condition based on synthetic aperture radar image-forming principles simulation the method for the original signal of Synthetic Aperture Radar Echo characteristic, detailed content can list of references: " InSAR echoed signal and system emulation are studied ", Zhang Jianqi, Harbin Institute of Technology's Master's thesis.

Definition 7, phase error vector and phase error matrix

Phase error vector refers to measures by SAR the vector that in echo data, each orientation forms to phase error, and phase error matrix refers to the diagonal matrix be made up of to phase error each orientation in SAR measurement data.

Definition 8, SAR observe scene object space

SAR observes scene object space refer to the geographic position set of all scene objects scattering points to be observed in realistic space.Observation scene object space has different expressions under different spaces coordinate system, but once it is unique that coordinate system establishes its expression later.Conveniently imaging, in the present invention, SAR observes scene object space be taken as earth axes.

Definition 9, SAR imaging space

SAR imaging space refer to the scattering point in scene objects space projected to orientation to-distance to space coordinates, this space is determined by two mutually orthogonal coordinate bases in SAR imaging space.

Definition 10, the sparse reconstructing method of compressed sensing

Higher-dimension original signal is mainly carried out non-self-adapting linear projection to lower dimensional space with the structural information of holding signal by compressed sensing, reconstruct the theory of original signal again by solving linear optimal solution, this theory mainly comprises sparse signal representation, sparseness measuring and sparse reconstruct three aspects.The basic thought of the sparse reconstructing method of compressed sensing is solve optimum solution under particular constraints condition or near-optimal solution, and main method has greedy tracing algorithm and convex relaxed algorithm etc.Detailed content can list of references " Donoho D L.Compressed sensing.IEEE Transactions on Information Theory, 2006,52 (4): 1289-1306. ".

Definition 11, diagonal matrix and unit matrix

First diagonal matrix is the square formation that line number is equal with columns, and the main diagonal element of matrix is not zero entirely simultaneously, and matrix non-master diagonal element is zero entirely.It is 1 entirely that unit matrix refers to main diagonal element, and non-master diagonal element is the diagonal matrix of 0 entirely.

Definition 12, image entropy

Image entropy is expressed as the bit average of image gray levels set, and also illustrate the average information of video source, image blur is larger, and image entropy is larger.Detailed content can list of references " Brink A D.Using spatial information as an aid to maximum entropy image threshold selection.Pattern Recognition Letters; 1996,17 (1): 29-36. "

Definition 13, Schmidt process

The basic thought of Schmidt process utilizes projection theory to construct a new orthogonal basis on the basis of existing orthogonal basis, one group of linearly independent vector become the method for a unit Orthogonal Vectors.Detailed content can with reference to " linear algebra ", and occupy remaining horse etc. and write, publishing house of Tsing-Hua University publishes.

The SAR sparse self-focusing formation method of a kind of image entropy optimum provided by the invention, it comprises the following steps:

Step 1, initialization SAR system parameter:

Initialization SAR system parameter comprises: platform speed vector, is denoted as V; Antenna initial position vector, namely orientation is to 0 slow moment aerial position, is denoted as P (0); Radar operating center frequency, is denoted as f _c; Radar carrier frequency wavelength, is denoted as λ; The signal bandwidth of radar emission baseband signal, is denoted as B _r; Radar emission signal pulse width, is denoted as T _p; The chirp rate of radar emission signal, is denoted as f _dr; Radar pulse repetition frequency, is denoted as PRF; Orientation, to antenna eliminator length, is denoted as D _a; The sample frequency of Radar Receiver System, is denoted as f _s; The aerial velocity of propagation of light, is denoted as C; Distance, to fast moment sequence, is denoted as t, t=1,2 ..., N _r, N _rfor distance is to fast moment sum; Orientation, to slow moment sequence, is denoted as l, l=1,2 ..., N _a, N _afor orientation is to slow moment sum.Above-mentioned parameter is SAR system canonical parameter, determines in SAR system design and observation process.According to SAR imaging system scheme and observation program, the initializes system parameters that SAR formation method needs is known.

Step 2, initialization SAR imaging space parameter and acquisition original echoed signals:

The imaging space parameter of initialization SAR, comprising: the volume coordinate formed using radar beam exposure field region ground level is as the imaging space of SAR, and this imaging space is designated as Ω; Platform, to the reference oblique distance at SAR imaging space center, is denoted as R _ref; Imaging space Ω is evenly divided into equal-sized flat unit lattice, and also referred to as resolution element, unit grid is designated as d respectively in horizontal cross, the longitudinal length of side of level _xand d _y, unit grid size is chosen as 1/2nd of SAR system traditional theory imaging resolution or SAR system traditional theory imaging resolution; In observation scene object space Ω, the coordinate vector of m cell, is denoted as P _m, m represents m cell in SAR imaging space Ω, m=1,2 ..., M, M are the cell sum in imaging space Ω; In SAR imaging space Ω, the scattering coefficient opsition dependent order of all cells rearranges vector, is denoted as α, and vectorial α is made up of capable 1 row of M; M element in scattering coefficient vector α, is denoted as α _m, m=1,2 ..., M.

According to platform speed vector V initialized in step 1, the pulse repetition rate PRF of antenna initial position vector P (0) and radar system, adopts formula P (l)=P (0)+Vl/PRF, l=1,2 ..., N _a, calculate the position vector of antenna in l orientation to the slow moment, be designated as P (l), N _afor the orientation of step 1 is to slow moment sum.Adopt formula R (P (l), P _m)=|| P (l)-P _m|| ₂, l=1,2 ..., N _a, m=1,2 ..., M, to calculate in l orientation to m cell in slow moment SAR imaging space Ω to the distance of antenna, is designated as R (P (l), P _m), wherein || || ₂represent the vectorial L2 norm in definition 2, P _mfor initialization obtains the coordinate vector of m cell in imaging space Ω, M is cell sum in initialized imaging space Ω.Adopt formula τ _m(l)=2R (P (l), P _m)/C, l=1,2 ..., N _a, m=1,2 ..., M, to calculate in l orientation to m cell in slow moment SAR imaging space Ω to the time delays of antenna, is designated as τ _ml (), wherein C is the aerial velocity of propagation of light that in step 1, initialization obtains.

L orientation is designated as s (t, l) to the raw radar data of SAR antenna in the fast moment, t=1,2 to slow moment and t distance ..., N _r, l=1,2 ..., N _a, wherein N _rfor distance initialized in step 1 is to fast moment sum.In SAR actual imaging, echo data s (t, l), t=1,2 ..., N _r, l=1,2 ..., N _a, can be provided by radar data receiver in SAR system.

Step 3, set up the linear measurement model of SAR echo signal:

Rearranging vector in order by obtaining all SAR original echoed signals s (t, l) in step 2, being designated as echoed signal vector S, echoed signal vector S is made up of capable 1 row of D, wherein D=N _an _r, N _rfor distance initialized in step 1 is to fast moment sum, N _afor orientation initialized in step 1 is to slow moment sum.

Adopt formula φ _d(m)=exp (-j2 π f _cτ _m(l)) exp (j π f _dr[t-τ _m(l)] ²), t=1,2 ..., N _r, l=1,2 ..., N _a, m=1,2 ..., M, d=1,2 ..., D, to calculate in imaging space Ω m cell at time delay function corresponding to echoed signal vector S d elemental signals, is designated as φ _dm (), wherein exp () represents e index sign of operation, f _cfor the radar operating center frequency that step 1 initialization obtains, f _drfor the chirp rate that transmits that step 1 initialization obtains, τ _m(l) for step 2 obtain in l orientation in slow moment SAR imaging space Ω m cell to the time delays of antenna, t be distance to t fast moment, j is imaginary unit's (namely-1 open root), and π is circular constant.

Calculation matrix in SAR original echoed signals vector S and imaging space Ω between all cell scattering coefficient vector α, is designated as A.Calculation matrix A is made up of the time delay function that all cells in SAR imaging space Ω are corresponding, and A is the two-dimensional matrix of the capable M row of D, and expression is

Wherein, φ ₁(1) for the 1st cell in imaging space Ω is at time delay function corresponding to echoed signal vector S the 1st elemental signals, φ ₁(2) for the 2nd cell in imaging space Ω is at time delay function corresponding to echoed signal vector S the 1st elemental signals, φ ₁(M) for M cell in imaging space Ω is at time delay function corresponding to echoed signal vector S the 1st elemental signals, φ ₂(1) for the 1st cell in imaging space Ω is at time delay function corresponding to echoed signal vector S the 2nd elemental signals, φ ₂(2) for the 2nd cell in imaging space Ω is at time delay function corresponding to echoed signal vector S the 2nd elemental signals, φ ₂(M) for M cell in imaging space Ω is at time delay function corresponding to echoed signal vector S the 2nd elemental signals, φ _d(1) for the 1st cell in imaging space Ω is at time delay function corresponding to echoed signal vector S D elemental signals, φ _d(2) for the 2nd cell in imaging space Ω is at time delay function corresponding to echoed signal vector S D elemental signals, φ _d(M) for M cell in imaging space Ω is at time delay function corresponding to echoed signal vector S D elemental signals, φ ₁(1), φ ₁(2) ..., φ ₁(M) the 1st, 2 are respectively in imaging space Ω ..., M cell at time delay functional vector corresponding to echoed signal vector S the 1st elemental signals, φ ₂(1), φ ₂(2) ..., φ ₂(M) the 1st, 2 are respectively in imaging space Ω ..., M cell at time delay functional vector corresponding to echoed signal vector S the 2nd elemental signals, φ _d(1), φ _d(2) ..., φ _d(M) the 1st, 2 are respectively in imaging space Ω ..., M cell is at time delay functional vector corresponding to echoed signal vector S D elemental signals.

Step 4, the scattering coefficient vector of SAR imaging space to be carried out according to a preliminary estimate:

Adopt formula calculate the initial estimate of scattering coefficient vector in SAR imaging space Ω, be denoted as wherein A ^hfor the conjugate transpose of matrix A, A is the SAR calculation matrix obtained in step 3, and subscript H is conjugate transpose operation symbol, and S is the SAR original echoed signals vector obtained in step 3.

Step 5, parameter needed for the sparse self-focusing imaging algorithm of initialization:

Initialization sparse self-focusing imaging algorithm desired parameters, comprising: the maximum iteration time of imaging processing iteration, is denoted as Maxiter; Phase error vector in imaging processing i-th iterative process in SAR imaging model, is designated as phase error vector Φ ⁽ⁱ⁾for the vector of capable 1 row of D, i is natural number, is expressed as i-th iteration of imaging processing, i=1,2 ..., Maxite, r element for vectorial Φ ⁽ⁱ⁾in the 1st element, element for vectorial Φ ⁽ⁱ⁾in the 2nd element, element for vectorial Φ ⁽ⁱ⁾in D element, upper right corner roman symbol T is transpose operation symbol; In imaging processing i-th iterative process, l orientation is to the phase error in slow moment, is designated as l=1,2 ..., N _a; The threshold value of imaging processing iteration, is denoted as τ; Adopt formula W ⁽ⁱ⁾=diag (exp (j Φ ⁽ⁱ⁾)), i=1,2 ..., Maxiter, calculates the initial phase error matrix value of SAR measurement model in imaging processing i-th iterative process, is denoted as W ⁽ⁱ⁾, i=1,2 ..., Maxiter, wherein phase error matrix W ⁽ⁱ⁾for the diagonal matrix of the capable D row of D, diag () is using the sign of operation of vector element as diagonal element in diagonal matrix; In imaging processing i-th iterative process, the image entropy of SAR image, is designated as L ⁽ⁱ⁾, i=1,2 ..., Maxiter.Before imaging iterative processing, in imaging processing the 0th iteration, the initial value of SAR measurement model phase error matrix, is designated as W ⁽⁰⁾, in imaging processing the 0th iteration, the initial value of SAR image entropy, is designated as L ⁽⁰⁾.

Step 6, estimate SAR orientation to phase error vector:

Comprise the following steps:

The response vector of step 6.1, acquisition target echo:

Adopt synthetic-aperture radar gauged distance compression method to SAR raw radar data s (t, l), t=1,2 ..., N _r, l=1,2 ..., N _a, carry out Range compress process, obtain the SAR echo data after Range compress, be denoted as s _r(t, l), t=1,2 ..., N _r, l=1,2 ..., N _a, wherein l orientation obtaining for step 2 of s (t, l) is to slow moment and t distance to the raw radar data of SAR antenna in the fast moment.

Adopt formula $\mathrm{\β}(l,m)=s_r(\mathrm{ceil}\left(\frac{R(P\left(l\right),P_m)-R_{\mathrm{ref}}}{2\·C\·f_s}\right),l)\·\exp (j\·\frac{4\·\mathrm{\π}\·R(P\left(l\right),P_m)}{\mathrm{\λ}}),$ L=1,2 ..., N _a, m=1,2 ..., M, and formula s _l=[β (l, 1), β (l, 2) ..., β (l, M)], l=1,2 ..., N _a, calculate the response vector of l orientation to target echo, be designated as s _lwherein β (l, 1) be that l orientation is to the 1st β that cell is corresponding (l in slow moment SAR imaging space Ω, m), β (l, 2) is l orientation to β (l, m) corresponding to the 2nd cell in slow moment SAR imaging space Ω, β (l, M) is ^lindividual orientation to β (l, m) corresponding to M cell in slow moment SAR imaging space Ω, rounding operation symbol in ceil () expression, R _refthe platform obtained for step 2 initialization is to SAR imaging space center reference oblique distance, and C is the aerial velocity of propagation of light that step 1 initialization obtains, f _sfor the sample frequency of the Radar Receiver System that step 1 initialization obtains, R (P (l), P _m) for l orientation obtaining in step 2, in slow moment SAR imaging space Ω, m cell is to the distance of antenna, λ is the wavelength that in step 1, initialization obtains.

Make q be natural number, q span is q=1,2 ..., N _a.As q=1, adopt formula calculate the data vector of q orientation to the slow moment, be designated as X.Work as q=N _atime, adopt formula calculate the data vector of q orientation to the slow moment, be designated as X; As 1 < q < N _atime, adopt formula calculate the data vector of q orientation to the slow moment, be designated as X, m the element of vectorial X is designated as x _m, wherein for l orientation in imaging processing i-th iterative process is to the phase error in slow moment, l=1,2 ..., N _a, s _lbe the response vector of l orientation to target echo.

Adopt formula Y=s _qcalculate the data vector of q orientation to the slow moment, be designated as Y, m the element of vectorial Y is designated as y _m, wherein s _qfor l orientation when l equals q is to the response vector of target echo, i.e. s _q=s _l.

Adopt formula m=1,2 ..., M, calculates element f _m; Adopt formula F=[f again ₁, f ₂..., f _m] data vector that calculates, be designated as F, wherein || ²for squared absolute value sign of operation, f ₁for element f corresponding during m=1 _m, f ₂for element f corresponding during m=2 _m, f _mfor element f corresponding during m=M _m, x _mfor m the element of vectorial X, y _mfor m the element of vectorial Y.

Step 6.2, Schmidt process is utilized to obtain unit orthogonal vector:

Adopt formula with m=1,2 ..., M, calculates element a _mand b _m, wherein x _mfor m element in the vectorial X that step 6.1 obtains, y _mfor m element in the vectorial Y that step 6.1 obtains; Adopt a=[a again ₁, a ₂..., a _m], b=[b ₁, b ₂..., b _m], m=1,2 ..., M, calculates vectorial a and vectorial b, a ₁for element a corresponding during m=1 _m, a ₂for element a corresponding during m=2 _m, a _mfor element a corresponding during m=M _m, b ₁for element b corresponding during m=1 _m, b ₂for element b corresponding during m=2 _m, b _mfor element b corresponding during m=M _m.

Utilize Schmidt process to carry out orthonormalization to vectorial a and vectorial b, obtain the unit orthogonal vector that the plane of vectorial a and vectorial b formation is corresponding wherein for vector the 1st element, for vector the 2nd element, for vector the 1st element, for vector the 2nd element.

Step 6.3, calculating quadratic form matrix and parameter estimation:

Adopt formula $r_1=\frac{{{\tilde{a}}_2}^2+{{\tilde{b}}_2}^2}{{\tilde{a}}_2\&CenterDot;{\tilde{b}}_1-{\tilde{a}}_1\&CenterDot;{\tilde{b}}_2},r_2=\frac{{{\tilde{a}}_1}^2+{{\tilde{b}}_1}^2}{{\tilde{a}}_2\&CenterDot;{\tilde{b}}_1-{\tilde{a}}_1\&CenterDot;{\tilde{b}}_2},r_3=-\frac{{\tilde{a}}_1\&CenterDot;{\tilde{a}}_2+{\tilde{b}}_1\&CenterDot;{\tilde{b}}_2}{{\tilde{a}}_2\&CenterDot;{\tilde{b}}_1-{\tilde{a}}_1\&CenterDot;{\tilde{b}}_2}$ With $R=\left[\begin{array}{cc}{r}_1& r_3\\ r_3& r_2\end{array}\right],$ Calculate quadratic form matrix R, wherein for vector the 1st element, for vector the 2nd element, for vector the 1st element, for vector the 2nd element, with be respectively the unit orthogonal vector that step 6.2 obtains.

Adopt formula Γ Λ Γ ^t=R carries out Eigenvalues Decomposition to quadratic form matrix R, and the eigenvectors matrix obtaining quadratic form matrix R is designated as Γ, Γ ^tfor the transposition of matrix Γ, the eigenvalue matrix obtaining quadratic form matrix R is designated as Λ.

Adopt formula η=(Γ x ₀-I) R ^-1calculate the estimated value of parameter, be designated as η, wherein x ₀for vectorial F vectorial a and vectorial b the intersection point point of two dimensional surface opened, I is unit matrix, R ^-1for the inverse matrix of quadratic form matrix R.

Step 6.4, estimation orientation are to phase error:

Adopt formula σ=[a, b] ^-1(η R+I) ^-1x ₀the vector calculated, is designated as σ, and wherein vectorial σ dimension is 2 × 1, and in vectorial σ, the 1st element is designated as σ ₁, in vectorial σ, the 2nd element is designated as σ ₂, a and b is the vector that step 6.2 obtains, and η is the estimated value vector that step 6.3 obtains, and R is the quadratic form matrix that step 6.3 obtains, x ₀for vector f vectorial a with b intersection point point in the two dimensional surface Ω that opens, I is unit matrix, subscript ^-1representing matrix inversion operation symbol.

Adopt formula calculate and to upgrade in imaging processing i-th iterative process the phase error vector of q orientation to the slow moment, wherein atan () is for solving tan inverse function sign of operation.By phase error vector Φ ⁽ⁱ⁾in (q-1) N _r+ 1 to qN _rindividual element the whole assignment of value be wherein N _rfor distance is to fast moment sum, for phase error vector Φ ⁽ⁱ⁾in (q-1) N _r+ 1 element, for phase error vector Φ ⁽ⁱ⁾in qN _rindividual element.

Step 6.5, orientation are estimated one by one to phase error:

For all orientation to sequence number q, q=1,2 ..., N _a, adopt step 6.1 to step 6.4 one by one orientation to until estimate that all orientation are to phase error, finally obtain orientation to phase error vector Φ ⁽ⁱ⁾.

The sparse imaging of step 7, compressed sensing:

Adopt formula W ⁽ⁱ⁾=diag (exp (j Φ ⁽ⁱ⁾)) calculate phase error matrix in imaging processing i-th iterative process, be denoted as W ⁽ⁱ⁾, wherein Φ ⁽ⁱ⁾for the phase error vector in imaging processing i-th iterative process that step 6.5 obtains, diag () is using the sign of operation of vector element as diagonal matrix diagonal element.

Adopt formula calculate the scattering coefficient vector in SAR target imaging space in imaging processing i-th iterative process with the sparse reconstructing method of standard compression sensing, be denoted as wherein represent the optimal value asking for corresponding independent variable α when meeting minimum value in bracket, vectorial S is the SAR data echoed signal vector that step 3 obtains, represent the square operation symbol of vectorial L2 norm, || || ₁represent vectorial L1 norm sign of operation.

Step 8, calculating SAR image entropy:

Adopt formula m=1,2 ..., M, and formula calculate SAR image entropy in imaging processing i-th iterative process, be denoted as L ⁽ⁱ⁾, i=1,2 ..., Maxiter, wherein for the scattering coefficient vector obtained in step 8 m element, i is i-th iteration of imaging processing, i=1,2 ..., Maxiter, M are the cell sum divided in the SAR imaging space Ω that obtains of step 2 initialization, for the square operation symbol of vectorial L2 norm, Σ () represents vector element summation operation symbol, log ₂() represents that the truth of a matter is the logarithm operation symbol of 2.

Step 9, imaging algorithm iterated conditional judge:

If satisfy condition || L ⁽ⁱ⁾-L ^(i-1)|| ₂>=τ and i≤Maxiter, then make the value of imaging processing iterations i add 1, and then repeated execution of steps 6 is to step 9, wherein L ⁽ⁱ⁾for SAR image entropy in imaging processing i-th iterative process, L ^(i-1)for SAR image entropy in imaging processing the i-th-1 time iterative process, τ is the algorithm iteration threshold value that in step 5, initialization obtains, and Maxiter is the imaging processing maximum iteration time that in step 5, initialization obtains, || || ₂for vectorial L2 norm sign of operation; If satisfy condition || L ⁽ⁱ⁾-L ^(i-1)|| ₂< τ or i > Maxiter, then imaging processing i-th iterative process obtain scattering coefficient vector with phase error vector Φ ⁽ⁱ⁾for final SAR imaging scattering coefficient vector sum phase error estimation and phase error result.

Innovative point of the present invention: in order to obtain good compressed sensing SAR image quality, contemplated by the invention orientation in SAR echo signal measurement model to the impact of phase error on compressed sensing SAR imaging.For the unknown phase estimation of error in compressed sensing SAR imaging model and compensation problem, the present invention to utilize in imaging model orientation to phase error properties and sparse target signature, utilize compressed sensing SAR image entropy as interpretational criteria, propose a kind of compressed sensing SAR sparse self-focusing formation method of image entropy optimum, the method utilizes SAR orientation to estimate that orientation is to phase error to the relation of echo and observed object in iterative process each time, then phase error compensation is carried out to compressed sensing imaging model, then compressed sensing SAR imaging is carried out again, and adopt iterative technique to make the image entropy of compressed sensing SAR imaging optimum, thus improve compressed sensing SAR image quality.

The invention has the advantages that and utilize image entropy optiaml ciriterion to achieve unknown orientation in compressed sensing SAR imaging model to the estimation of phase error and compensation problem, improve compressed sensing SAR image quality and robustness.

Accompanying drawing illustrates:

Fig. 1 is the compressed sensing SAR sparse self-focusing formation method schematic flow sheet of a kind of image entropy optimum provided by the present invention

Fig. 2 estimates the schematic flow sheet of SAR orientation to phase error vector method in the present invention

Fig. 3 is the SAR system simulation parameter table that the specific embodiment of the invention adopts

Embodiment

The present invention mainly adopts the method for emulation experiment to verify, institute all verifies correct with conclusion in steps on MATLABR2008b software.Concrete implementation step is as follows:

Step 1, initialization SAR system parameter:

Initialization SAR system parameter comprises: platform speed vector V=[0,150,0] m/s; Antenna initial position vector P (0)=[0,0,6000] m; Radar operating center frequency f _c=10 × 10 ⁹hz; Radar carrier frequency wavelength X=0.03m; The signal bandwidth B of radar emission baseband signal _r=1.5 × 10 ⁸hz; Radar emission signal pulse width T _p=1 × 10 ^-6s; The chirp rate f of radar emission signal _dr=1.5 × 10 ¹⁴hz/s; The sample frequency f of Radar Receiver System _s=3 × 10 ⁸hz; Radar pulse repetition frequency PRF=500Hz; Orientation is D to antenna eliminator length _a=2m; Aerial velocity of propagation C=3 × 10 of light ⁸m/s; Distance is to fast moment sum N _r=128, distance to fast moment sequence t=1,2 ..., N _r; Orientation is to slow moment sum N _a=128, orientation to slow moment sequence l=1,2 ..., N _a.Above-mentioned parameter is SAR system canonical parameter, determines in SAR system design and observation process.According to SAR imaging system scheme and observation program, the initializes system parameters that SAR formation method needs is known.

Step 2, initialization SAR imaging space parameter and acquisition original echoed signals:

The imaging space parameter of initialization SAR, comprising: the volume coordinate formed using radar beam exposure field region ground level is as the imaging space of SAR, and this imaging space is designated as Ω; The size of initialization SAR imaging space Ω is 128 × 128 pixels, and the centre coordinate position of SAR imaging space Ω is positioned at [8000,0,0] m, and each unit grid is d in horizontal cross and the longitudinal length of side of level _x=d _y=0.5m, calculate total cell number M=16384 of SAR imaging space, in SAR imaging space Ω, the position of each cell is P _m=[(x '-64) 0.5, (y '-64) 0.5,0] m, wherein x '=1 ..., 2, y '=1,2 ..., 128, m=(x '-1) 128+y '.Scene center is with reference to oblique distance R _ref=10000m.P _mfor the position vector of m cell in SAR imaging space Ω, m represents m cell in SAR imaging space Ω, m=1,2 ..., M, M=16384.In SAR imaging space Ω, add simulated point target scattering body, the several number of point target scatterer is 5, and their scattering coefficient values are 1, coordinate position is respectively [0,0 ,], [20,20,0], [20 ,-21,0], [-20,20,0], [-20 ,-21,0], unit is m; The scattering coefficient not comprising point target cell in SAR imaging space Ω is set to 0.The target scattering coefficient of all cells in SAR imaging space Ω is rearranged scattering vector α by cell location order.After determining all unit scattering coefficients of SAR imaging space Ω, scattering coefficient vector α just determines in SAR imaging observation simulation process.Scene objects scattering coefficient vector α is made up of capable 1 row of M, α _mfor the scattering coefficient value of m cell in SAR imaging space Ω corresponding in vectorial α.In the Ω of this simulation imaging space, 5 the cell scattering coefficient value α only comprising point scattering target are set to 1, and the scattering coefficient of remaining element lattice is all 0.Traditional synthetic-aperture radar original echo emulation mode is utilized to produce the raw radar data of SAR.

According to platform speed vector V=[0 that initialization in step 1 obtains, 150,0] m/s, initial position vector P (0)=[0,0 of antenna, 6000] m and pulse repetition rate PRF=500Hz, adopt formula P (l)=P (0)+Vl/PRF, l=1,2,, N _a, calculate the position vector of antenna in l orientation to the slow moment, be designated as P (l), N _afor the orientation of step 1 is to slow moment sum N _a=128.Adopt formula R (P (l), P _m)=|| P (l)-P _m|| ₂, l=1,2 ..., N _a, m=1,2 ..., M, M=16384, to calculate in l orientation to m cell in slow moment SAR imaging space Ω to the distance of antenna, be designated as R (P (l), P _m), wherein || || ₂represent the vectorial L2 norm in definition 2, P _mfor initialization obtains the coordinate vector of m cell in imaging space Ω, M is cell sum M=16384 in initialized imaging space Ω.Adopt formula τ _m(l)=2R (P (l), P _m)/C, l=1,2 ..., N _a, m=1,2 ..., M, N _a=128, M=16384, to calculate in l orientation to m cell in slow moment SAR imaging space Ω to the time delays of antenna, is designated as τ _m(l), wherein C=3 × 10 ⁸m/s.

In simulation process or actual imaging process, obtain SAR raw radar data, be designated as s (t, l) to the raw radar data of SAR antenna in the fast moment, t=1,2 in l orientation to slow moment and t distance ..., N _r, l=1,2 ..., N _a, wherein N _rfor distance initialized in step 1 is to fast moment sum N _r=128, N _afor orientation initialized in step 1 is to slow moment sum N _a=128.In SAR actual imaging, s (t, l) can be provided by radar data receiver in SAR system; And in simulation process, s (t, l) adopts conventional synthesis aperture radar original echo emulation mode in definition 6 to produce and provides.

Step 3, set up the linear measurement model of SAR echo signal:

SAR original echoed signals s (t, l) obtained in step 2 is rearranged in order echoed signal vector S, echoed signal vector S is made up of capable 1 row of D, wherein D=N _an _r=16384, N _rfor distance initialized in step 1 is to fast moment sum N _r=128, N _afor orientation initialized in step 1 is to slow moment sum N _a=128.

Adopt formula φ _d(m)=exp (-j2 π f _cτ _m(l)) exp (j π f _dr[t-τ _m(l)] 2), t=1,2 ..., N _r, l=1,2 ..., N _a, m=1,2 ..., M, d=1,2 ..., D, N _r=128, N _a=128, M=16384, D=16384, to calculate in imaging space Ω m cell at time delay function corresponding to echoed signal vector S d elemental signals, be designated as φ _dm (), wherein exp () represents e index sign of operation, f _cfor the radar operating center frequency that step 1 initialization obtains, f _drfor the chirp rate that transmits that step 1 initialization obtains, τ _m(l) for step 2 obtain in l orientation in slow moment SAR imaging space Ω m cell to the time delays of antenna, t be distance to t fast moment, j is imaginary unit's (namely-1 open root), and π is circular constant, is about π=3.1416.

Adopt expression matrix formula

Calculate the calculation matrix A of SAR original echoed signals and all cells of SAR imaging space Ω, calculation matrix A is the two-dimensional matrix of the capable M row of D.

Step 4, the scattering coefficient vector of SAR imaging space to be carried out according to a preliminary estimate:

Adopt formula calculate the initial estimate of the scattering coefficient vector of SAR imaging space Ω, be denoted as wherein A ^hfor the conjugate transpose of matrix A, A is the SAR calculation matrix obtained in step 3, and subscript H is conjugate transpose operation symbol, and S is the SAR original echoed signals vector obtained in step 3.

Step 5, parameter needed for the sparse self-focusing imaging algorithm of initialization:

Initialization sparse self-focusing imaging algorithm desired parameters, comprising: the maximum iteration time Maxiter=15 of imaging iterative processing; Phase error vector in imaging processing i-th iterative process in SAR imaging model initial value be null vector, phase error vector Φ ⁽ⁱ⁾for the vector of capable 1 row of D, i is natural number, is expressed as i-th iteration of imaging processing, i=1,2 ..., Maxiter, element for vectorial Φ ⁽ⁱ⁾in the 1st element, element for vectorial Φ ⁽ⁱ⁾in the 2nd element, element for vectorial Φ ⁽ⁱ⁾in D element, upper right corner roman symbol T is transpose operation symbol; In imaging processing i-th iterative process, l orientation is to the phase error in slow moment, is designated as l=1,2 ..., N _a, N _a=128; Threshold tau=0.01 of imaging processing iteration; Adopt formula W ⁽ⁱ⁾=diag (exp (j Φ ⁽ⁱ⁾)), i=1,2 ..., Maxiter, calculates the initial phase error matrix value of SAR measurement model in imaging processing i-th iterative process, is denoted as W ⁽ⁱ⁾, i=1,2 ..., Maxiter, Maxiter=15, wherein phase error matrix W ⁽ⁱ⁾for the diagonal matrix of the capable D row of D, diag () is using the sign of operation of vector element as diagonal element in diagonal matrix; In imaging processing i-th iterative process, the image entropy of SAR image, is designated as L ⁽ⁱ⁾, i=1,2 ..., Maxiter, Maxiter=15.Before imaging processing iteration, SAR measurement model phase error matrix W in imaging processing the 0th iteration ⁽⁰⁾initial value be null matrix, the initial value L of SAR image entropy in imaging processing the 0th iteration ⁽⁰⁾=0.

Step 6, estimate SAR orientation to phase error vector:

Comprise the following steps:

The response vector of step 6.1, acquisition target echo:

Adopt standard synthetic aperture radar gauged distance compression method to SAR raw radar data s (t, l), t=1,2 ..., N _r, l=1,2 ..., N _a, carry out Range compress process, obtain the SAR echo data after Range compress, be denoted as s _r(t, l), t=1,2 ..., N _r, l=1,2 ..., N _a, N _r=128, N _a=128, wherein l orientation obtaining for step 2 of s (t, l) is to slow moment and t distance to the raw radar data of SAR antenna in the fast moment.

Adopt formula $\mathrm{\&beta;}(l,m)=s_r(\mathrm{ceil}\left(\frac{R(P\left(l\right),P_m)-R_{\mathrm{ref}}}{2\&CenterDot;C\&CenterDot;f_s}\right),l)\&CenterDot;\exp (j\&CenterDot;\frac{4\&CenterDot;\mathrm{\&pi;}\&CenterDot;R(P\left(l\right),P_m)}{\mathrm{\&lambda;}}),$ L=1,2 ..., N _a, m=1,2 ..., M, and formula s _l=[β (l, 1), β (l, 2) ..., β (l, M)], l=1,2 ..., N _a, calculate the response vector of l orientation to target echo, be designated as s _lwherein β (l, 1) is l orientation to β (l, m) corresponding to the 1st cell in slow moment SAR imaging space Ω, β (l, 2) be l orientation to β (l, m) corresponding to the 2nd cell in slow moment SAR imaging space Ω, β (l, M) be l orientation to β (l corresponding to M cell in slow moment SAR imaging space Ω, m), rounding operation symbol in ceil () expression, R _refthe platform obtained for step 2 initialization is to SAR imaging space center reference oblique distance, and C is the aerial velocity of propagation of light that step 1 initialization obtains, f _sfor the sample frequency of the Radar Receiver System that step 1 initialization obtains, R (P (l), P _m) for l orientation obtaining in step 2, in slow moment SAR imaging space Ω, m cell is to the distance of antenna, λ is the wavelength that in step 1, initialization obtains.

Make q be natural number, q span is q=1,2 ..., N _a, N _a=128.As q=1, adopt formula calculate the data vector of q orientation to the slow moment, be designated as X.Work as q=N _atime, adopt formula calculate the data vector of q orientation to the slow moment, be designated as X; As 1 < q < N _atime, adopt formula calculate the data vector of q orientation to the slow moment, be designated as X, m the element of vectorial X is designated as x _m, wherein for l orientation in imaging processing i-th iterative process is to the phase error in slow moment, l=1,2 ..., N _a, s _lbe the response vector of l orientation to target echo.

Adopt formula Y=s _qcalculate the data vector of q orientation to the slow moment, be designated as Y, m the element of vectorial Y is designated as ym, wherein s _qfor l orientation when l equals q is to the response vector of target echo, i.e. s _q=s _l.

Adopt formula m=1,2 ..., M, calculates element f _m; Adopt formula F=[f again ₁, f ₂..., f _m] data vector that calculates, be designated as F, wherein || ²for squared absolute value sign of operation, f ₁for element f corresponding during m=1 _m, f ₂for element f corresponding during m=2 _m, f _mfor element f corresponding during m=M _m, x _mfor m the element of vectorial X, y _mfor m the element of vectorial Y.

Step 6.2, Schmidt process is utilized to obtain unit orthogonal vector:

Adopt formula with m=1,2 ..., M, M=16384, calculate element a _mand b _m, wherein x _mfor m element in the vectorial X that step 6.1 obtains, y _mfor m element in the vectorial Y that step 6.1 obtains; Adopt a=[a again ₁, a ₂..., a _m], b=[b ₁, b ₂..., b _m], m=1,2 ..., M, M=16384, calculate vectorial a and vectorial b, a ₁for element a corresponding during m=1 _m, a ₂for element a corresponding during m=2 _m, a _mfor element a corresponding during m=M _m, b ₁for element b corresponding during m=1 _m, b ₂for element b corresponding during m=2 _m, b _mfor element b corresponding during m=M _m.

Utilize Schmidt process to carry out orthonormalization to vectorial a and vectorial b, obtain the unit orthogonal vector that the plane of vectorial a and vectorial b formation is corresponding wherein for vector the 1st element, for vector the 2nd element, for vector the 1st element, for vector the 2nd element.

Step 6.3, calculating quadratic form matrix and parameter estimation:

Adopt formula $r_1=\frac{{{\tilde{a}}_2}^2+{{\tilde{b}}_2}^2}{{\tilde{a}}_2\&CenterDot;{\tilde{b}}_1-{\tilde{a}}_1\&CenterDot;{\tilde{b}}_2},r_2=\frac{{{\tilde{a}}_1}^2+{{\tilde{b}}_1}^2}{{\tilde{a}}_2\&CenterDot;{\tilde{b}}_1-{\tilde{a}}_1\&CenterDot;{\tilde{b}}_2},r_3=-\frac{{\tilde{a}}_1\&CenterDot;{\tilde{a}}_2+{\tilde{b}}_1\&CenterDot;{\tilde{b}}_2}{{\tilde{a}}_2\&CenterDot;{\tilde{b}}_1-{\tilde{a}}_1\&CenterDot;{\tilde{b}}_2}$ With $R=\left[\begin{array}{cc}{r}_1& r_3\\ r_3& r_2\end{array}\right],$ Calculate quadratic form matrix R, wherein for vector the 1st element, for vector the 2nd element, for vector the 1st element, for vector the 2nd element, with be respectively the unit orthogonal vector that step 6.2 obtains.

Adopt formula Γ Λ Γ ^t=R carries out Eigenvalues Decomposition to quadratic form matrix R, and the eigenvectors matrix obtaining quadratic form matrix R is designated as Γ, Γ ^tfor the transposition of matrix Γ, the eigenvalue matrix obtaining quadratic form matrix R is designated as Λ.

Adopt formula η=(Γ x ₀-I) R ^-1calculate the estimated value of parameter, be designated as η, wherein x ₀for vectorial F vectorial a and vectorial b the intersection point point of two dimensional surface opened, I is unit matrix, R ^-1for the inverse matrix of quadratic form matrix R.

Step 6.4, estimation orientation are to phase error:

Adopt formula σ=[a, b] ^-1(η R+I) ^-1x ₀the vector calculated, is designated as σ, and wherein vectorial σ dimension is 2 × 1, and in vectorial σ, the 1st element is designated as σ ₁, in vectorial σ, the 2nd element is designated as σ ₂, a and b is the vector that step 6.2 obtains, and η is the estimated value vector that step 6.3 obtains, and R is the quadratic form matrix that step 6.3 obtains, x ₀for vector f vectorial a with b intersection point point in the two dimensional surface Ω that opens, I is unit matrix, subscript ^-1representing matrix inversion operation symbol.

Adopt formula calculate and to upgrade in imaging processing i-th iterative process the phase error vector of q orientation to the slow moment, wherein atan () is for solving tan inverse function sign of operation.By phase error vector Φ ⁽ⁱ⁾in (q-1) N _r+ 1 to qN _rindividual element the whole assignment of value be wherein N _rfor distance is to fast moment sum, for phase error vector Φ ⁽ⁱ⁾in (q-1) N _r+ 1 element, for phase error vector Φ ⁽ⁱ⁾in qN _rindividual element.

Step 6.5, orientation are estimated one by one to phase error:

For all orientation to sequence number q, q=1,2 ..., N _a, N _a=128, adopt step 6.1 to step 6.4 one by one orientation to until estimate that all orientation are to phase error, finally obtain orientation to phase error vector Φ ⁽ⁱ⁾.

The sparse imaging of step 7, compressed sensing:

Adopt formula W ⁽ⁱ⁾=diag (exp (j Φ ⁽ⁱ⁾)) calculate phase error matrix in imaging processing i-th iterative process, be denoted as W ⁽ⁱ⁾, wherein Φ ⁽ⁱ⁾for the phase error vector in imaging processing i-th iterative process that step 6.5 obtains, diag () is using the sign of operation of vector element as diagonal matrix diagonal element.

Adopt formula calculate the scattering coefficient vector in SAR target imaging space in imaging processing i-th iterative process with the sparse reconstructing method of standard compression sensing, be denoted as wherein represent the optimal value asking for corresponding independent variable α when meeting minimum value in bracket, vectorial S is the SAR data echoed signal vector that step 3 obtains, represent the square operation symbol of vectorial L2 norm, || || ₁represent vectorial L1 norm sign of operation.

Step 8, calculating SAR image entropy:

Adopt formula m=1,2 ..., M, M=16384, and formula calculate SAR image entropy in imaging processing i-th iterative process, be denoted as L ⁽ⁱ⁾, i=1,2 ..., Maxiter, Maxiter=15, wherein for the scattering coefficient vector obtained in step 8 m element, i is i-th iteration of algorithm, i=1,2 ..., 15, M is cell sum M=16384 divided in the SAR imaging space Ω that obtains of step 2 initialization, for the square operation symbol of vectorial L2 norm, Σ () represents vector element summation operation symbol, log ₂() represents that the truth of a matter is the logarithm operation symbol of 2.

Step 9, imaging processing iterated conditional judge:

If satisfy condition || L ⁽ⁱ⁾-L ^(i-1)|| ₂>=τ and i≤Maxiter, then make the value of imaging processing iterations i add 1, and then repeated execution of steps 6 is to step 9, wherein L ⁽ⁱ⁾for SAR image entropy in imaging processing i-th iterative process, L ^(i-1)for SAR image entropy in imaging processing the i-th-1 time iterative process, τ is algorithm iteration threshold tau=0.01 that in step 5, initialization obtains, and Maxiter is the imaging processing maximum iteration time Maxiter=15 that in step 5, initialization obtains, || || ₂for vectorial L2 norm sign of operation; If satisfy condition || L ⁽ⁱ⁾-L ^(i-1)|| ₂< τ or i > Maxiter, then imaging processing i-th iterative process obtain scattering coefficient vector with phase error vector Φ ⁽ⁱ⁾for final SAR imaging scattering coefficient vector sum phase error estimation and phase error result.

Claims (1)

1. a SAR sparse self-focusing formation method for image entropy optimum, is characterized in that it comprises the steps:

Step 1, initialization SAR system parameter:

Initialization SAR system parameter comprises: platform speed vector, is denoted as V; Antenna initial position vector, namely orientation is to 0 slow moment aerial position, is denoted as P (0); Radar operating center frequency, is denoted as f _c; Radar carrier frequency wavelength, is denoted as λ; The signal bandwidth of radar emission baseband signal, is denoted as B _r; Radar emission signal pulse width, is denoted as T _p; The chirp rate of radar emission signal, is denoted as f _dr; Radar pulse repetition frequency, is denoted as PRF; Orientation, to antenna eliminator length, is denoted as D _a; The sample frequency of Radar Receiver System, is denoted as f _s; The aerial velocity of propagation of light, is denoted as C; Distance, to fast moment sequence, is denoted as t, t=1,2 ..., N _r, N _rfor distance is to fast moment sum; Orientation, to slow moment sequence, is denoted as l, l=1,2 ..., N _a, N _afor orientation is to slow moment sum; Above-mentioned parameter is SAR system canonical parameter, determines in SAR system design and observation process; According to SAR imaging system scheme and observation program, the initializes system parameters that SAR formation method needs is known;

Step 2, initialization SAR imaging space parameter and acquisition original echoed signals:

The imaging space parameter of initialization SAR, comprising: the volume coordinate formed using radar beam exposure field region ground level is as the imaging space of SAR, and this imaging space is designated as Ω; Platform, to the reference oblique distance at SAR imaging space center, is denoted as R _ref; Imaging space Ω is evenly divided into equal-sized flat unit lattice, and also referred to as resolution element, unit grid is designated as d respectively in horizontal cross, the longitudinal length of side of level _xand d _y, unit grid size is chosen as 1/2nd of SAR system traditional theory imaging resolution or SAR system traditional theory imaging resolution; In observation scene object space Ω, the coordinate vector of m cell, is denoted as P _m, m represents m cell in SAR imaging space Ω, m=1,2 ..., M, M are the cell sum in imaging space Ω; In SAR imaging space Ω, the scattering coefficient opsition dependent order of all cells rearranges vector, is denoted as α, and vectorial α is made up of capable 1 row of M; M element in scattering coefficient vector α, is denoted as α _m, m=1,2 ..., M;

According to platform speed vector V initialized in step 1, the pulse repetition rate PRF of antenna initial position vector P (0) and radar system, adopts formula P (l)=P (0)+Vl/PRF, l=1,2 ..., N _a, calculate the position vector of antenna in l orientation to the slow moment, be designated as P (l), N _afor the orientation of step 1 is to slow moment sum; Adopt formula R (P (l), P _m)=|| P (l)-P _m|| ₂, l=1,2 ..., N _a, m=1,2 ..., M, to calculate in l orientation to m cell in slow moment SAR imaging space Ω to the distance of antenna, is designated as R (P (l), P _m), wherein || || ₂represent the vectorial L2 norm in definition 2, P _mfor initialization obtains the coordinate vector of m cell in imaging space Ω, M is cell sum in initialized imaging space Ω; Adopt formula τ _m(l)=2R (P (l), P _m)/C, l=1,2 ..., N _a, m=1,2 ..., M, to calculate in l orientation to m cell in slow moment SAR imaging space Ω to the time delays of antenna, is designated as τ _m(l), wherein C is the aerial velocity of propagation of light that in step 1, initialization obtains;

L orientation is designated as s (t, l) to the raw radar data of SAR antenna in the fast moment, t=1,2 to slow moment and t distance ..., N _r, l=1,2 ..., N _a, wherein N _rfor distance initialized in step 1 is to fast moment sum; In SAR actual imaging, echo data s (t, l), t=1,2 ..., N _r, l=1,2 ..., N _a, can be provided by radar data receiver in SAR system;

Step 3, set up the linear measurement model of SAR echo signal:

Rearranging vector in order by obtaining all SAR original echoed signals s (t, l) in step 2, being designated as echoed signal vector S, echoed signal vector S is made up of capable 1 row of D, wherein D=N _an _r, N _rfor distance initialized in step 1 is to fast moment sum, N _afor orientation initialized in step 1 is to slow moment sum;

Adopt formula t=1,2 ..., N _r, l=1,2 ..., N _a, m=1,2 ..., M, d=1,2 ..., D, to calculate in imaging space Ω m cell at time delay function corresponding to echoed signal vector S d elemental signals, is designated as φ _dm (), wherein exp () represents e index sign of operation, f _cfor the radar operating center frequency that step 1 initialization obtains, f _drfor the chirp rate that transmits that step 1 initialization obtains, τ _m(l) for step 2 obtain in l orientation in slow moment SAR imaging space Ω m cell to the time delays of antenna, t be distance to t fast moment, j is imaginary unit's (namely-1 open root), and π is circular constant;

Calculation matrix in SAR original echoed signals vector S and imaging space Ω between all cell scattering coefficient vector α, is designated as A; Calculation matrix A is made up of the time delay function that all cells in SAR imaging space Ω are corresponding, and A is the two-dimensional matrix of the capable M row of D, and expression is

Wherein, φ ₁(1) for the 1st cell in imaging space Ω is at time delay function corresponding to echoed signal vector S the 1st elemental signals, φ ₁(2) for the 2nd cell in imaging space Ω is at time delay function corresponding to echoed signal vector S the 1st elemental signals, φ ₁(M) for M cell in imaging space Ω is at time delay function corresponding to echoed signal vector S the 1st elemental signals, φ ₂(1) for the 1st cell in imaging space Ω is at time delay function corresponding to echoed signal vector S the 2nd elemental signals, φ ₂(2) for the 2nd cell in imaging space Ω is at time delay function corresponding to echoed signal vector S the 2nd elemental signals, φ ₂(M) for M cell in imaging space Ω is at time delay function corresponding to echoed signal vector S the 2nd elemental signals, φ _d(1) for the 1st cell in imaging space Ω is at time delay function corresponding to echoed signal vector S D elemental signals, φ _d(2) for the 2nd cell in imaging space Ω is at time delay function corresponding to echoed signal vector S D elemental signals, φ _d(M) for M cell in imaging space Ω is at time delay function corresponding to echoed signal vector S D elemental signals, φ ₁(1), φ ₁(2) ..., φ ₁(M) the 1st, 2 are respectively in imaging space Ω ..., M cell at time delay functional vector corresponding to echoed signal vector S the 1st elemental signals, φ ₂(1), φ ₂(2) ..., φ ₂(M) the 1st, 2 are respectively in imaging space Ω ..., M cell at time delay functional vector corresponding to echoed signal vector S the 2nd elemental signals, φ _d(1), φ _d(2) ..., φ _d(M) the 1st, 2 are respectively in imaging space Ω ..., M cell is at time delay functional vector corresponding to echoed signal vector S D elemental signals;

Step 4, the scattering coefficient vector of SAR imaging space to be carried out according to a preliminary estimate:

Adopt formula calculate the initial estimate of scattering coefficient vector in SAR imaging space Ω, be denoted as wherein A ^hfor the conjugate transpose of matrix A, A is the SAR calculation matrix obtained in step 3, and subscript H is conjugate transpose operation symbol, and S is the SAR original echoed signals vector obtained in step 3;

Step 5, parameter needed for the sparse self-focusing imaging algorithm of initialization:

Initialization sparse self-focusing imaging algorithm desired parameters, comprising: the maximum iteration time of imaging processing iteration, is denoted as Maxiter; Phase error vector in imaging processing i-th iterative process in SAR imaging model, is designated as phase error vector Φ ⁽ⁱ⁾for the vector of capable 1 row of D, i is natural number, is expressed as i-th iteration of imaging processing, i=1,2 ..., Maxiter, element for vectorial Φ ⁽ⁱ⁾in the 1st element, element for vectorial Φ ⁽ⁱ⁾in the 2nd element, element for vectorial Φ ⁽ⁱ⁾in D element, upper right corner roman symbol T is transpose operation symbol; In imaging processing i-th iterative process, l orientation is to the phase error in slow moment, is designated as l=1,2 ..., N _a; The threshold value of imaging processing iteration, is denoted as τ; Adopt formula W ⁽ⁱ⁾=diag (exp (j Φ ⁽ⁱ⁾)), i=1,2 ..., Maxiter, calculates the initial phase error matrix value of SAR measurement model in imaging processing i-th iterative process, is denoted as W ⁽ⁱ⁾, i=1,2 ..., Maxiter, wherein phase error matrix W ⁽ⁱ⁾for the diagonal matrix of the capable D row of D, diag () is using the sign of operation of vector element as diagonal element in diagonal matrix; In imaging processing i-th iterative process, the image entropy of SAR image, is designated as L ⁽ⁱ⁾, i=1,2 ..., Maxiter; Before imaging iterative processing, in imaging processing the 0th iteration, the initial value of SAR measurement model phase error matrix, is designated as W ⁽⁰⁾, in imaging processing the 0th iteration, the initial value of SAR image entropy, is designated as L ⁽⁰⁾;

Step 6, estimate SAR orientation to phase error vector:

Comprise the following steps:

The response vector of step 6.1, acquisition target echo:

Adopt synthetic-aperture radar gauged distance compression method to SAR raw radar data s (t, l), t=1,2 ..., N _r, l=1,2 ..., N _a, carry out Range compress process, obtain the SAR echo data after Range compress, be denoted as s _r(t, l), t=1,2 ..., N _r, l=1,2 ..., N _a, wherein l orientation obtaining for step 2 of s (t, l) is to slow moment and t distance to the raw radar data of SAR antenna in the fast moment;

Adopt formula $\mathrm{\&beta;}(l,m)=s_r(\mathrm{ceil}\left(\frac{R(P\left(l\right),P_m)-R_{\mathrm{ref}}}{2\&CenterDot;C\&CenterDot;f_s}\right),l)\&CenterDot;\exp (j\&CenterDot;\frac{4\&CenterDot;\mathrm{\&pi;}\&CenterDot;R(P\left(l\right),P_m)}{\mathrm{\&lambda;}}),$ L=1,2 ..., N _a, m=1,2 ..., M, and formula s _l=[β (l, 1), β (l, 2) ..., β (l, M)], l=1,2 ..., N _a, calculate the response vector of l orientation to target echo, be designated as s _lwherein β (l, 1) is l orientation to β (l, m) corresponding to the 1st cell in slow moment SAR imaging space Ω, β (l, 2) be l orientation to β (l, m) corresponding to the 2nd cell in slow moment SAR imaging space Ω, β (l, M) be l orientation to β (l corresponding to M cell in slow moment SAR imaging space Ω, m), rounding operation symbol in ceil () expression, R _refthe platform obtained for step 2 initialization is to SAR imaging space center reference oblique distance, and C is the aerial velocity of propagation of light that step 1 initialization obtains, f _sfor the sample frequency of the Radar Receiver System that step 1 initialization obtains, R (P (l), P _m) for l orientation obtaining in step 2, in slow moment SAR imaging space Ω, m cell is to the distance of antenna, λ is the wavelength that in step 1, initialization obtains;

Make q be natural number, q span is q=1,2 ..., N _a; As q=1, adopt formula calculate the data vector of q orientation to the slow moment, be designated as X; Work as q=N _atime, adopt formula calculate the data vector of q orientation to the slow moment, be designated as X; As 1 < q < N _atime, adopt formula calculate the data vector of q orientation to the slow moment, be designated as X, m the element of vectorial X is designated as x _m, wherein for l orientation in imaging processing i-th iterative process is to the phase error in slow moment, l=1,2 ..., N _a, s _lbe the response vector of l orientation to target echo;

Adopt formula Y=s _qcalculate the data vector of q orientation to the slow moment, be designated as Y, m the element of vectorial Y is designated as y _m, wherein s _qfor l orientation when l equals q is to the response vector of target echo, i.e. s _q=s _l;

Adopt formula m=1,2 ..., M, calculates element f _m; Adopt formula F=[f again ₁, f ₂..., f _m] data vector that calculates, be designated as F, wherein || ²for squared absolute value sign of operation, f ₁for element f corresponding during m=1 _m, f ₂for element f corresponding during m=2 _m, f _mfor element f corresponding during m=M _m, x _mfor m the element of vectorial X, y _mfor m the element of vectorial Y;

Step 6.2, Schmidt process is utilized to obtain unit orthogonal vector:

Adopt formula with m=1,2 ..., M, calculates element a _mand b _m, wherein x _mfor m element in the vectorial X that step 6.1 obtains, y _mfor m element in the vectorial Y that step 6.1 obtains; Adopt a=[a again ₁, a ₂..., a _m], b=[b ₁, b ₂..., b _m], m=1,2 ..., M, calculates vectorial a and vectorial b, a ₁for element a corresponding during m=1 _m, a ₂for element a corresponding during m=2 _m, a _mfor element a corresponding during m=M _m, b ₁for element b corresponding during m=1 _m, b ₂for element b corresponding during m=2 _m, b _mfor element b corresponding during m=M _m;

Utilize Schmidt process to carry out orthonormalization to vectorial a and vectorial b, obtain the unit orthogonal vector that the plane of vectorial a and vectorial b formation is corresponding wherein for vector the 1st element, for vector the 2nd element, for vector the 1st element, for vector the 2nd element;

Step 6.3, calculating quadratic form matrix and parameter estimation:

Adopt formula $r_1=\frac{{{\tilde{a}}_2}^2+{{\tilde{b}}_2}^2}{{\tilde{a}}_2\&CenterDot;{\tilde{b}}_1-{\tilde{a}}_1\&CenterDot;{\tilde{b}}_2},r_2=\frac{{{\tilde{a}}_1}^2+{{\tilde{b}}_1}^2}{{\tilde{a}}_2\&CenterDot;{\tilde{b}}_1-{\tilde{a}}_1\&CenterDot;{\tilde{b}}_2},r_3=-\frac{{\tilde{a}}_1\&CenterDot;{\tilde{a}}_2+{\tilde{b}}_1\&CenterDot;{\tilde{b}}_2}{{\tilde{a}}_2\&CenterDot;{\tilde{b}}_1-{\tilde{a}}_1\&CenterDot;{\tilde{b}}_2}$ With $R=\left[\begin{array}{cc}{r}_1& r_3\\ r_3& r_2\end{array}\right],$ Calculate quadratic form matrix R, wherein for vector the 1st element, for vector the 2nd element, for vector the 1st element, for vector the 2nd element, with be respectively the unit orthogonal vector that step 6.2 obtains;

Adopt formula Γ Λ Γ ^t=R carries out Eigenvalues Decomposition to quadratic form matrix R, and the eigenvectors matrix obtaining quadratic form matrix R is designated as Γ, Γ ^tfor the transposition of matrix Γ, the eigenvalue matrix obtaining quadratic form matrix R is designated as Λ;

Adopt formula η=(Γ x ₀-I) R ^-1calculate the estimated value of parameter, be designated as η, wherein x ₀for vectorial F vectorial a and vectorial b the intersection point point of two dimensional surface opened, I is unit matrix, R ^-1for the inverse matrix of quadratic form matrix R;

Step 6.4, estimation orientation are to phase error:

Adopt formula σ=[a, b] ^-1(η R+I) ^-1x ₀the vector calculated, is designated as σ, and wherein vectorial σ dimension is 2 × 1, and in vectorial σ, the 1st element is designated as σ ₁, in vectorial σ, the 2nd element is designated as σ ₂, a and b is the vector that step 6.2 obtains, and η is the estimated value vector that step 6.3 obtains, and R is the quadratic form matrix that step 6.3 obtains, x ₀for vector f vectorial a with b intersection point point in the two dimensional surface Ω that opens, I is unit matrix, subscript ^-1representing matrix inversion operation symbol;

Adopt formula calculate and to upgrade in imaging processing i-th iterative process the phase error vector of q orientation to the slow moment, wherein atan () is for solving tan inverse function sign of operation; By phase error vector Φ ⁽ⁱ⁾in (q-1) N _r+ 1 to qN _rindividual element the whole assignment of value be wherein N _rfor distance is to fast moment sum, for phase error vector Φ ⁽ⁱ⁾in (q-1) N _r+ 1 element, for phase error vector Φ ⁽ⁱ⁾in qN _rindividual element;

Step 6.5, orientation are estimated one by one to phase error:

For all orientation to sequence number q, q=1,2 ..., N _a, adopt step 6.1 to step 6.4 one by one orientation to until estimate that all orientation are to phase error, finally obtain orientation to phase error vector Φ ⁽ⁱ⁾;

The sparse imaging of step 7, compressed sensing:

Adopt formula W ⁽ⁱ⁾=diag (exp (j Φ ⁽ⁱ⁾)) calculate phase error matrix in imaging processing i-th iterative process, be denoted as W ⁽ⁱ⁾, wherein Φ ⁽ⁱ⁾for the phase error vector in imaging processing i-th iterative process that step 6.5 obtains, diag () is using the sign of operation of vector element as diagonal matrix diagonal element;

Adopt formula calculate the scattering coefficient vector in SAR target imaging space in imaging processing i-th iterative process with the sparse reconstructing method of standard compression sensing, be denoted as wherein represent the optimal value asking for corresponding independent variable α when meeting minimum value in bracket, vectorial S is the SAR data echoed signal vector that step 3 obtains, represent the square operation symbol of vectorial L2 norm, || || ₁represent vectorial L1 norm sign of operation;

Step 8, calculating SAR image entropy:

Adopt formula m=1,2 ..., M, and formula calculate SAR image entropy in imaging processing i-th iterative process, be denoted as L ⁽ⁱ⁾, i=1,2 ..., Maxiter, wherein for the scattering coefficient vector obtained in step 8 m element, i is i-th iteration of imaging processing, i=1,2 ..., Maxiter, M are the cell sum divided in the SAR imaging space Ω that obtains of step 2 initialization, for the square operation symbol of vectorial L2 norm, Σ () represents vector element summation operation symbol, log ₂() represents that the truth of a matter is the logarithm operation symbol of 2;

Step 9, imaging algorithm iterated conditional judge:

If satisfy condition || L ⁽ⁱ⁾-L ^(i-1)|| ₂>=τ and i≤Maxiter, then make the value of imaging processing iterations i add 1, and then repeated execution of steps 6 is to step 9, wherein L ⁽ⁱ⁾for SAR image entropy in imaging processing i-th iterative process, L ^(i-1)for SAR image entropy in imaging processing the i-th-1 time iterative process, τ is the algorithm iteration threshold value that in step 5, initialization obtains, and Maxiter is the imaging processing maximum iteration time that in step 5, initialization obtains, || || ₂for vectorial L2 norm sign of operation;

If satisfy condition || L ⁽ⁱ⁾-L ^(i-1)|| ₂< τ or i > Maxiter, then imaging processing i-th iterative process obtain scattering coefficient vector with phase error vector Φ ⁽ⁱ⁾for final SAR imaging scattering coefficient vector sum phase error estimation and phase error result.