CN104007440B - One accelerated decomposition rear orientation projection spot beam SAR formation method - Google Patents

Abstract

The invention belongs to synthetic aperture radar image-forming technical field, particularly a kind of accelerated decomposition rear orientation projection spot beam SAR formation method. This accelerated decomposition rear orientation projection spot beam SAR formation method comprises the following steps: utilize platform motion to form synthetic aperture, set up polar coordinate system (r, θ), and full aperture is divided into N taking aperture center as initial point₀Individual isometric sub-aperture; Make α=sin θ, by u (r corresponding to sub-aperture_p, θ) and the impulse response function of locating pixel is expressed as I^(u)(α)＝I^(u)(r_p, α); Define two-dimentional wave number variable; Draw u two-dimentional wave-number spectrum corresponding to sub-aperture; For each two-dimentional wave-number spectrum corresponding apart from wave number, to each two-dimentional wave-number spectrum corresponding to sub-aperture spliced, draw each one dimension wave-number spectrum corresponding apart from wave number along orientation; Then by each apart from one dimension wave-number spectrum corresponding to wave number along distance to splicing, draw full aperture two dimension wave-number spectrum; Described full aperture two dimension wave-number spectrum is carried out to two-dimensional Fourier transform, draw total space resolution image under polar coordinate system.

Description

Accelerated decomposition backward projection bunching synthetic aperture radar imaging method

Technical Field

The invention belongs to the technical field of Synthetic Aperture Radar (SAR) imaging, and particularly relates to an accelerated decomposition backward projection bunching synthetic aperture radar imaging method which can be used for an airborne or spaceborne SAR imaging processing platform and is suitable for complex imaging occasions such as nonlinear flight paths, large squints and the like.

Background

The Back-Projection (BP) method realizes the continuous accumulation of energy of each pixel point on an imaging grid through integration along the slope course (called BP integration), and is an accurate bunching SAR imaging algorithm. The time domain back projection method ideally solves the coupling problem of distance and direction, is suitable for image reconstruction under complex imaging geometries such as large coherent accumulation angle or nonlinear flight path, and the focused image has no geometric distortion. However, the computational burden of the time domain backprojection method limits its application to large data scale imaging scenarios.

In order to improve the efficiency of the time domain back projection method, some fast back projection methods are successively proposed, among which fast decomposition back projection (FFBP) method is most representative. The fast decomposition back projection method divides the BP integral on the full aperture into step and segment integrals on the sub-apertures, namely aperture decomposition and recursive fusion. The aperture decomposition establishes a butterfly algorithm structure, which is a precondition for the rapidity of a quick decomposition back projection method; the recursive fusion realizes the formation of the image resolution from low to high, and is an important way for reconstructing a two-dimensional full-resolution image. Based on the thought, the fast decomposition back projection method greatly reduces the operation burden brought by BP integral, and obtains the operation efficiency close to the frequency domain algorithm. In the initial stage and the processing stage, the fast decomposition back projection method establishes a local polar coordinate system by taking the center of each sub-aperture as an origin. For a certain point target in an imaging scene, the point corresponds to different distance and angular domain coordinates under different coordinate systems. Therefore, when the sub-images corresponding to different apertures are fused, a mapping relation from the corresponding coordinate system to the new coordinate system needs to be established, and the process is essentially a two-dimensional processing process depending on distance and angular domain interpolation. However, two-dimensional (point-by-point) interpolation to achieve image fusion causes two main problems: 1) the interpolation operation inevitably introduces errors, and the errors are continuously accumulated and amplified along with the progress of recursive fusion, and finally the loss of image quality is caused; 2) two-dimensional interpolation is a very time-consuming part of the fast decomposition backprojection method. When the number of the apertures divided in the initial stage is less, the interpolation error accumulated by image recursive fusion is smaller, the image quality is improved therewith, but the corresponding operation efficiency is reduced; when the number of the apertures divided in the initial stage is large, the number of times of image recursive fusion is increased, the accumulated interpolation error is large, and the image quality is reduced accordingly.

Disclosure of Invention

The invention aims to provide an accelerated decomposition backward projection bunching synthetic aperture radar imaging method, which is an accelerated fast backward projection method (AFBP), and is abbreviated as an AFBP method. The AFBP method avoids the accumulation and amplification of interpolation errors and accurately keeps the original form of a wave number spectrum because two-dimensional interpolation and recursive fusion are not needed. Therefore, the AFBP method has higher operation efficiency and more ideal focusing performance than the fast decomposition back projection method. The invention is realized by the following steps:

analyzing the frame configuration and processing links of the existing fast decomposition backward projection method, and finding out the part with high operation complexity in the algorithm; the AFBP method optimizes and improves the part with high operation complexity in the existing fast decomposition back projection method, and the algorithm efficiency is accelerated. The method comprises the following specific steps:

in order to achieve the technical purpose, the invention is realized by adopting the following technical scheme.

An accelerated decomposition backward projection bunching synthetic aperture radar imaging method comprises the following steps:

s1: receiving an echo signal by using a synthetic aperture radar on a motion platform, and performing matched filtering on the echo signal to obtain a range pulse compression signal, wherein the motion platform is an airplane or a satellite;

s2: establishing a polar coordinate system with the synthetic aperture center O as an origin, wherein the coordinate is (r)_pTheta) is represented by (r)_pTheta) pixel point, dividing the full aperture into N₀Sub-apertures of equal length, N₀Setting α to sin theta, and making the u-th sub-aperture correspond to (r)_pTheta) the impulse response function of the pixel point is represented as I^(u)(r_pα), u is 1 to N₀(ii) a Let I^(u)(α)＝I^(u)(r_p,α)；

S3: defining a distance wavenumber variable K_rAnd the angular wavenumber variable K_α(ii) a By the pair I^(u)(α) performing inverse Fourier transform to obtain a two-dimensional wavenumber spectrum H corresponding to the u-th sub-aperture^(u)(K_r,K_α) (ii) a Splicing the two-dimensional wave number spectrums corresponding to the sub-apertures in the two-dimensional wave number spectrums corresponding to the sub-apertures along the azimuth direction aiming at the two-dimensional wave number spectrums corresponding to the sub-apertures to obtain one-dimensional wave number spectrums corresponding to the sub-apertures; then, the one-dimensional wave number spectrum corresponding to each distance wave number is spliced along the distance direction to obtain a full-aperture two-dimensional wave number spectrum; and performing two-dimensional Fourier transform on the full-aperture two-dimensional wave number spectrum to obtain a full-space resolution image under a polar coordinate system.

The invention is characterized by further improvement:

in step S1, when the moving platform flies at a constant speed, a synthetic aperture with a length L is formed, and the center of the synthetic aperture is represented as O; in the flying process of the moving platform, the synthetic aperture radar transmits signals outwards and receives echo signals to obtain demodulated echo signals;

in step S1, the process of obtaining the demodulated echo signal is: establishing a plane rectangular coordinate system by taking the center O of the synthetic aperture as an origin, and forming a rectangular coordinate on the planeIn the system, the X-axis direction is the motion direction of the motion platform, and the absolute value of the X-axis coordinate value at the antenna phase center is X; then, a polar coordinate system is established by taking the center O of the synthetic aperture as the origin, and the phase center of the antenna reaches a point target P (r)_p,θ_p) Is determined by the instantaneous slope distance R (X; r is_p,θ_p) Comprises the following steps:

$R(X;r_p,{\mathrm{\&theta;}}_p)=\sqrt{r_p^2+X^2-2r_{p}X\sin {\mathrm{\&theta;}}_p},-\frac{L}{2}\&le;X<\frac{L}{2}$

wherein r is_pAnd theta_pTwo coordinates of the point target P under a polar coordinate system respectively, and then an angular sine coordinate α is defined_p，α_p＝sinθ_pThen, the above formula is rewritten as:

$R(X;r_p,{\mathrm{\&alpha;}}_p)=\sqrt{r_p^2+X^2-2r_{p}X{\mathrm{\&alpha;}}_p},-\frac{L}{2}\&le;X<\frac{L}{2}$

the demodulated echo signal s (τ, X) is then:

$s(\mathrm{\&tau;},X)=\mathrm{rect}\left[\frac{X}{L}\right]\&CenterDot;s_T(\mathrm{\&tau;}-{\mathrm{\&Delta;t}}_p)$

wherein s is_T(. to) is the emission signal waveform, τ is the fast time, Δ t_p＝2R(X；r_p,α_p) C, c is the speed of light, rect [. C]Representing a rectangular window function;

in step S1, after the demodulated echo signal is obtained, the demodulated echo signal is subjected to matched filtering to obtain a range pulse compressed signal S_M(τ,X)：

$s_M(\mathrm{\&tau;},X)=\sin c[B\&CenterDot;(\mathrm{\&tau;}-\mathrm{\&Delta;}{t}_p)]\&CenterDot;\mathrm{rect}\left[\frac{X}{L}\right]\&CenterDot;\exp [-j{K}_{\mathrm{rc}}R(X;r_p,{\mathrm{\&alpha;}}_p)]$

Wherein B is the transmission signal bandwidth of the synthetic aperture radar, K_rc＝4π/λ_mid，λ_midThe wavelength corresponding to the carrier frequency when the synthetic aperture radar transmits signals.

In step S1, when the moving platform flies at a constant speed, a synthetic aperture with a length L is formed, and the center of the synthetic aperture is represented as O; in the flying process of the moving platform, the synthetic aperture radar transmits signals outwards and receives echo signals to obtain demodulated echo signals;

in step S1, the process of obtaining the demodulated echo signal is: establishing a plane rectangular coordinate system by taking a synthetic aperture center O as an origin, wherein in the plane rectangular coordinate system, the X-axis direction is the motion direction of the motion platform, and the absolute value of the X-axis coordinate value at the antenna phase center is X; then, a polar coordinate system is established by taking the center O of the synthetic aperture as the origin, and the phase center of the antenna reaches a point target P (r)_p,θ_p) Is determined by the instantaneous slope distance R (X; r is_p,θ_p) Comprises the following steps:

$R(X;r_p,{\mathrm{\&theta;}}_p)=\sqrt{r_p^2+X^2-2r_{p}X\sin {\mathrm{\&theta;}}_p},-\frac{L}{2}\&le;X<\frac{L}{2}$

wherein r is_pAnd theta_pTwo coordinates of the point target P under a polar coordinate system respectively, and then an angular sine coordinate α is defined_p，α_p＝sinθ_pThen, the above formula is rewritten as:

$R(X;r_p,{\mathrm{\&alpha;}}_p)=\sqrt{r_p^2+X^2-2r_{p}X{\mathrm{\&alpha;}}_p},-\frac{L}{2}\&le;X<\frac{L}{2}$

the demodulated echo signal s (τ, X) is then:

$s(\mathrm{\&tau;},X)=\mathrm{rect}\left[\frac{X}{L}\right]\&CenterDot;s_T(\mathrm{\&tau;}-{\mathrm{\&Delta;t}}_p)$

wherein s is_T(. to) is the emission signal waveform, τ is the fast time, Δ t_p＝2R(X；r_p,α_p) C, c is the speed of light, rect [. C]Representing a rectangular window function;

in step S1, after obtaining the demodulated echo signal, calculating a pitch error caused by the motion error according to platform parameters, where the platform parameters include a position, an attitude, a speed, and an acceleration of the motion platform; according to the slant range error caused by the motion error, performing motion compensation on the demodulated echo signal to obtain a motion-compensated echo signal; then, the echo signals after the motion compensation are subjected to matched filtering to obtain a range pulse compression signal s_M(τ,X)：

$s_M(\mathrm{\&tau;},X)=\sin c[B\&CenterDot;(\mathrm{\&tau;}-\mathrm{\&Delta;}{t}_p)]\&CenterDot;\mathrm{rect}\left[\frac{X}{L}\right]\&CenterDot;\exp [-j{K}_{\mathrm{rc}}R(X;r_p,{\mathrm{\&alpha;}}_p)]$

Wherein B is the transmission signal bandwidth of the synthetic aperture radar, K_rc＝4π/λ_mid，λ_midThe wavelength corresponding to the carrier frequency when the synthetic aperture radar transmits signals.

In step S2, the full aperture is first divided into N₀Sub-apertures of equal length, N₀Is a natural number greater than 1, and the length of each sub-aperture is l;

in the polar coordinate system, the coordinate is (r)_pTheta) is represented by (r)_pTheta), let α be sin theta, then in the polar coordinate system, the u-th sub-aperture corresponds to (r_pTheta) impulse response function I of the pixel^(u)(r_pα) is:

$I^{\left(u\right)}(r_p,\mathrm{\&alpha;})={\&Integral;}_{-l/2}^{l/2}{s}_M[\frac{2R(X+x_u;r_p,\mathrm{\&alpha;})}{c},X+x_u]\&CenterDot;\exp [j{K}_{\mathrm{rc}}\&CenterDot;R(X+x_u;r_p,\mathrm{\&alpha;})]\mathrm{dX}$

wherein x is_u＝(u-N₀L, u is 1 to N₀And c is the speed of light; using X to X + X in the above formula_uIf the variables are replaced, the above formula is rewritten as:

$\begin{array}{c}{I}^{\left(u\right)}(r_p,\mathrm{\&alpha;})={\&Integral;}_{x_u-l/2}^{x_u+l/2}{s}_M[\frac{2R(X;r_p,\mathrm{\&alpha;})}{c},X]\&CenterDot;\exp \left[j{K}_{\mathrm{rc}}R(X;r_p,\mathrm{\&alpha;})\right]\mathrm{dX}\\ ={\&Integral;}_{x_u-l/2}^{x_u+l/2}\exp [{\mathrm{jK}}_{\mathrm{rc}}\&CenterDot;\mathrm{\&Delta;R}(X;r_p,\mathrm{\&alpha;})]\mathrm{dX}\end{array}$

wherein, Delta R (X; R)_p,α)＝R(X；r_p,α)-R(X；r_p,α_p) (ii) a Then for Δ R (X; R)_pα) performing a second order taylor series expansion, resulting in:

$\begin{array}{c}\mathrm{\&Delta;R}(X;r_p,\mathrm{\&alpha;})=R(X;r_p,\mathrm{\&alpha;})-R(X;r_p,{\mathrm{\&alpha;}}_p)\\ =-(\mathrm{\&alpha;}-{\mathrm{\&alpha;}}_p)X-\frac{{\mathrm{\&alpha;}}^2-{\mathrm{\&alpha;}}_p^2}{2r_p}{X}^2\end{array}$

ignoring X in the above equation²Item, then I^(u)(r_pα) to:

$I^{\left(u\right)}(r_p,\mathrm{\&alpha;})={\&Integral;}_{x_u-l/2}^{x_u+l/2}\exp [-{\mathrm{jK}}_{\mathrm{rc}}\&CenterDot;(\mathrm{\&alpha;}-{\mathrm{\&alpha;}}_p)X]\mathrm{dX}$

then, the angular wave number domain variables are defined as follows:

$\left\{\begin{array}{c}{K}_{\mathrm{\&alpha;}}=K_{\mathrm{rc}}\&CenterDot;X\\ \mathrm{\&Delta;}{K}_{\mathrm{\&alpha;}}=K_{\mathrm{rc}}\&CenterDot;l\\ K_{\mathrm{\&alpha;c}}^{\left(u\right)}=K_{\mathrm{rc}}\&CenterDot;x_u\end{array}\right.$

let I^(u)(α)＝I^(u)(r_pα), mixing K_αSubstitution into I^(u)(α), then I^(u)(α) is:

$I^{\left(u\right)}\left(\mathrm{\&alpha;}\right)=I^{\left(u\right)}(r_p,\mathrm{\&alpha;})={\&Integral;}_{K_{\mathrm{\&alpha;c}}^{\left(u\right)}-\mathrm{\&Delta;}{K}_{\mathrm{\&alpha;}}/2}^{K_{\mathrm{\&alpha;c}}^{\left(u\right)}+\mathrm{\&Delta;}{K}_{\mathrm{\&alpha;}}/2}\exp [-j{K}_{\mathrm{\&alpha;}}\&CenterDot;(\mathrm{\&alpha;}-{\mathrm{\&alpha;}}_p)]{\mathrm{dK}}_{\mathrm{\&alpha;}}$

wherein, Δ α is the angular imaging range for each sub-aperture for the sine of the view angle from the center of the u-th sub-aperture to the center of the scene.

In step S3, the two-dimensional wavenumber variables are first defined as follows:

$\left\{\begin{array}{cc}{K}_{\mathrm{\&alpha;}}=K_r\&CenterDot;X& \\ \mathrm{\&Delta;}{K}_{\mathrm{\&alpha;}}=K_r\&CenterDot;l& K_r\&Element;[\frac{4\mathrm{\&pi;}}{{\mathrm{\&lambda;}}_{\max }},\frac{4\mathrm{\&pi;}}{{\mathrm{\&lambda;}}_{\min }}]\\ K_{\mathrm{\&alpha;c}}^{\left(u\right)}=K_r\&CenterDot;x_u& \end{array}\right.$

wherein λ is_maxFor the wavelength, λ, corresponding to the minimum frequency at which the synthetic aperture radar transmits signals_minThe wavelength corresponding to the maximum frequency when the synthetic aperture radar transmits signals; k_r4 pi/lambda, wherein lambda is the wavelength of the synthetic aperture radar when transmitting signals;

by the pair I^(u)(α) performing inverse Fourier transform to obtain a two-dimensional wave number spectrum H corresponding to the u-th sub-aperture^(u)(K_r,K_α)，H^(u)(K_r,K_α) Comprises the following steps:

$\begin{array}{c}{H}^{\left(u\right)}\left(K_{\mathrm{\&alpha;}}\right)={\&Integral;}_{{\mathrm{\&alpha;}}_c^{\left(u\right)}-\mathrm{\&Delta;\&alpha;}/2}^{{\mathrm{\&alpha;}}_c^{\left(u\right)}+\mathrm{\&Delta;\&alpha;}/2}{I}^{\left(u\right)}\left(\mathrm{\&alpha;}\right)\&CenterDot;\exp \left(j{K}_{\mathrm{\&alpha;}}\mathrm{\&alpha;}\right)\mathrm{d\&alpha;}\\ =\mathrm{rect}\left[\frac{K_{\mathrm{\&alpha;}}-K_{\mathrm{\&alpha;c}}^{\left(u\right)}}{\mathrm{\&Delta;}{K}_{\mathrm{\&alpha;}}}\right]\&CenterDot;\exp (-j{K}_{\mathrm{\&alpha;}}{\mathrm{\&alpha;}}_p)\end{array}$

$\left\{\begin{array}{c}{K}_{\mathrm{\&alpha;}}\&Element;[-\frac{\mathrm{\&Delta;}{K}_{\mathrm{\&alpha;}}}{2}+K_{\mathrm{\&alpha;c}}^{\left(u\right)},\frac{\mathrm{\&Delta;}{K}_{\mathrm{\&alpha;}}}{2}+K_{\mathrm{\&alpha;c}}^{\left(u\right)}]\\ K_r\&Element;[-\frac{\mathrm{\&Delta;}{K}_r}{2}+K_{\mathrm{rc}},\frac{\mathrm{\&Delta;}{K}_r}{2}+K_{\mathrm{rc}}]\end{array}\right.$

wherein, Δ K_rAnd c is the speed of light 4 pi B/c.

After step S3, converting the full-space-resolved image in the polar coordinate system into the planar rectangular coordinate system by two-dimensional interpolation, so as to obtain the full-space-resolved image in the planar rectangular coordinate system.

The invention has the beneficial effects that:

(1) according to the method, the imaging plane is selected in a r-sin theta coordinate system, the position of the point to be projected on a sin theta axis is directly calculated according to the sine value of the projection visual angle, and the projection visual angle theta is not required to be calculated. The method avoids the pixel-by-pixel operation of the arcsine of the projection visual angle, and effectively reduces the calculation burden of the algorithm.

(2) Using the global r-sin θ coordinates, the invention establishes the aperture position X and the angular wavenumber K_αFourier transform pair relationship between them, i.e. K_r＝K_αX, wherein K_rIs the distance wavenumber. Due to K_αThe above formula can provide a one-to-one correspondence relationship between the phase history and the azimuth time corresponding to the distance compression phase history domain, which is also a premise that the invention can realize aperture fusion through the wave number spectrum.

(3) Based on Fourier transform and cyclic shift, the invention can realize image fusion through one-dimensional processing of the wave number spectrum without two-dimensional interpolation and recursive fusion, thereby avoiding the damage of interpolation errors to the wave number spectrum form. Therefore, the method has the operation efficiency close to a frequency domain algorithm and the focusing performance comparable to a time domain back projection method.

Drawings

FIG. 1 is a flow chart of an accelerated decomposition backward projection beaming synthetic aperture radar imaging method of the present invention;

FIG. 2 is a geometrical diagram of sub-aperture imaging of the synthetic aperture radar of the present invention in a polar coordinate system;

FIG. 3 is a schematic diagram of image fusion achieved by one-dimensional angular wavenumber spectrum stitching according to the present invention;

FIG. 4 is a schematic diagram of image fusion achieved by two-dimensional spectral stitching according to the present invention;

FIG. 5 is a geometric diagram of the distribution of imaging regions and point targets in a simulation experiment of the present invention;

FIG. 6a is a schematic diagram of an azimuth impulse response function of three points P1, P2 and P3 reconstructed by three methods when the number of aperture positions of each sub-aperture is 32 in a simulation experiment; FIG. 6b is a schematic diagram of an azimuth impulse response function of three points P1, P2 and P3 reconstructed by three methods when the number of aperture positions of each sub-aperture is 16 in a simulation experiment; fig. 6c is a schematic diagram of an azimuth impulse response function of three points P1, P2 and P3 reconstructed by using three methods when the number of aperture positions of each sub-aperture is 8 in a simulation experiment;

FIG. 7a is a full-space resolution image in a polar coordinate system reconstructed by a fast decomposition back-projection method in an actually measured data test; FIG. 7b is a full-space resolution image in a polar coordinate system reconstructed by the present invention in a measured data test;

FIG. 8a is a full-space resolution image under a polar coordinate system of a sub-scene A obtained by a fast decomposition back projection method in an actually measured data test; FIG. 8b is a full-space resolution image of a polar coordinate system of a sub-scene A obtained by using the present invention in an actual measurement data test; FIG. 8c is a full-space resolution image of a polar coordinate system of a sub-scene B obtained by a fast decomposition back-projection method in an actual measurement data test; FIG. 8d is a full-space resolution image of a polar coordinate system of a sub-scene B obtained by using the present invention in an actual measurement data test; FIG. 8e is a full-space resolution image under a polar coordinate system of a sub-scene C obtained by a fast decomposition back-projection method in an actually measured data test; FIG. 8f is a full-space resolution image under a polar coordinate system of a sub-scene C obtained by the present invention in an actual measurement data test;

FIG. 9a is a schematic diagram of an imaging result after two-dimensional interpolation, which is obtained by two point targets in an actual measurement data experiment by using a fast decomposition back projection method; FIG. 9b is a schematic diagram of an imaging result of two-dimensional interpolation of two point targets obtained by the present invention in an actual measurement data experiment; FIG. 9c is a schematic diagram of an azimuth impulse response function obtained by two methods for two point targets in an actual measurement data experiment;

FIG. 10a is a schematic diagram of an imaging result of a point target 3 after two-dimensional interpolation, which is obtained by a fast decomposition back projection method in an actual measurement data experiment; FIG. 10b is a schematic diagram of an imaging result of the point target 3 after two-dimensional interpolation in the measured data experiment; fig. 10c is a schematic diagram of the azimuth impulse response function of the point target 3 obtained by two methods in the measured data experiment.

Detailed Description

The invention will be further described with reference to the accompanying drawings in which:

referring to fig. 1, a flow chart of an accelerated decomposition backward projection beaming synthetic aperture radar imaging method of the present invention is shown. The accelerated decomposition backward projection bunching synthetic aperture radar imaging method comprises the following steps:

s1: and receiving the echo signal by using a synthetic aperture radar on a motion platform, and performing matched filtering on the echo signal to obtain a range pulse compression signal, wherein the motion platform is an airplane or a satellite. The concrete description is as follows:

and a synthetic aperture radar is arranged on the moving platform. When the moving platform flies at a constant speed v (the flying direction of the moving platform is represented as the x direction), a synthetic aperture with the length of L is formed, the center of the synthetic aperture is represented as O, in the flying process of the moving platform, the synthetic aperture radar transmits signals outwards, and the beam center of an antenna (namely a receiving antenna and a transmitting antenna) of the synthetic aperture radar always points to the center C of a scene. Referring to fig. 2, it is a geometrical diagram of sub-aperture imaging of the synthetic aperture radar of the present invention in a polar coordinate system. Establishing a polar coordinate system with the synthetic aperture center O as the origin, and then the antenna phase center to the point target P (r)_p,θ_p) Has an instantaneous slope of

$R(X;r_p,{\mathrm{\&theta;}}_p)=\sqrt{r_p^2+X^2-2r_{p}X\sin {\mathrm{\&theta;}}_p},-\frac{L}{2}\&le;X<\frac{L}{2}$

Wherein r is_pAnd theta_pRespectively two coordinates of the point target P under a polar coordinate system; the meaning of X is explained below: establishing a plane rectangular coordinate system by taking the center O of the synthetic aperture as an origin, wherein in the plane rectangular coordinate system, the direction of an X axis is the movement direction of the movement platform, and the X axis is the skyThe absolute value of the x-axis coordinate value of the line phase center.

Defining angular sine coordinates α_p，α_p＝sinθ_pThen the above formula can be rewritten as

$R(X;r_p,{\mathrm{\&alpha;}}_p)=\sqrt{r_p^2+X^2-2r_{p}X{\mathrm{\&alpha;}}_p},-\frac{L}{2}\&le;X<\frac{L}{2}$

The echo signal after demodulation is

$s(\mathrm{\&tau;},X)=\mathrm{rect}\left[\frac{X}{L}\right]\&CenterDot;s_T(\mathrm{\&tau;}-{\mathrm{\&Delta;t}}_p)$

Wherein s is_TIs a transmission signal waveForm, τ is fast time, Δ t_pRepresenting two-way delay, Δ t_p＝2R(X；r_p,α_p) C, c is the speed of light, rect [. C]A rectangular window function is represented.

And then, calculating the slope distance error caused by the motion error according to platform parameters recorded by a GPS (satellite positioning system) or an INS (inertial navigation system), wherein the platform parameters comprise the position, the attitude, the speed and the acceleration of the motion platform. And performing motion compensation on the echo signal according to the slant range error caused by the motion error to obtain the echo signal after the motion compensation.

Performing matched filtering on the echo signal after motion compensation to obtain a range pulse compression signal s_M(τ,X)：

$s_M(\mathrm{\&tau;},X)=\sin c[B\&CenterDot;(\mathrm{\&tau;}-\mathrm{\&Delta;}{t}_p)]\&CenterDot;\mathrm{rect}\left[\frac{X}{L}\right]\&CenterDot;\exp [-j{K}_{\mathrm{rc}}R(X;r_p,{\mathrm{\&alpha;}}_p)]$

Wherein B is the transmission signal bandwidth of the synthetic aperture radar, K_rcIs a distance from the center of the wave number, K_rc＝4π/λ_mid，λ_midWhich is the wavelength corresponding to the center frequency (i.e., carrier frequency) at which the synthetic aperture radar transmits signals.

S2: establishing a polar coordinate system with the synthetic aperture center O as an origin, wherein the coordinate is (r)_pTheta) is represented by (r)_pθ) pixel points; dividing full aperture into N₀Sub-apertures of equal length, N₀Setting α to sin theta, and making the u-th sub-aperture correspond to (r)_pTheta) the impulse response function of the pixel point is represented as I^(u)(r_pα), let I^(u)(α)＝I^(u)(r_pα) as follows:

in step S2, the full aperture is first divided into N₀Sub-apertures of equal length, N₀Is a natural number greater than 1, and the length of each sub-aperture is L, L ═ L/N₀。

The invention projects the pulse compression signals corresponding to each sub-aperture backwards to a polar coordinate system with the center of the synthetic aperture as the origin. In the polar coordinate system, the coordinate is (r)_pTheta) is represented by (r)_pα is made to be sin theta, so as to establish a r-sin theta coordinate system, and then (r) is converted into a new vector_pTheta) pixel point is converted into r-sin theta coordinate system_pTheta) the coordinate of the pixel point in the r-sin theta coordinate system is (r)_pα), will have a coordinate of (r) in the r-sin θ coordinate system_pα) is represented as (r)_pα).

Then (r) corresponds to the u-th sub-aperture in the polar coordinate system_pTheta) impulse response function I of the pixel^(u)(r_pα) is:

$I^{\left(u\right)}(r_p,\mathrm{\&alpha;})={\&Integral;}_{-l/2}^{l/2}{s}_M[\frac{2R(X+x_u;r_p,\mathrm{\&alpha;})}{c},X+x_u]\&CenterDot;\exp [j{K}_{\mathrm{rc}}\&CenterDot;R(X+x_u;r_p,\mathrm{\&alpha;})]\mathrm{dX}$

wherein x is_u＝(u-N₀L/2-1/2). From I^(u)(r_pα) is given by the formula I^(u)(r_pα) is obtained by following the course of the slope distance R (X + X)_u；r_pα), corresponding integral variable X ∈ [ -l/2, l/2)]. By X to I^(u)(r_pα) of the calculation formula_uBy carrying out a variable substitution, then I^(u)(r_pα) is further rewritten as:

$\begin{array}{c}{I}^{\left(u\right)}(r_p,\mathrm{\&alpha;})={\&Integral;}_{x_u-l/2}^{x_u+l/2}{s}_M[\frac{2R(X;r_p,\mathrm{\&alpha;})}{c},X]\&CenterDot;\exp \left[j{K}_{\mathrm{rc}}R(X;r_p,\mathrm{\&alpha;})\right]\mathrm{dX}\\ ={\&Integral;}_{x_u-l/2}^{x_u+l/2}\exp [{\mathrm{jK}}_{\mathrm{rc}}\&CenterDot;\mathrm{\&Delta;R}(X;r_p,\mathrm{\&alpha;})]\mathrm{dX}\end{array}$

wherein, Delta R (X; R)_pα) is the projection slope distance, Δ R (X; R)_p,α)＝R(X；r_p,α)-R(X；r_p,α_p). For Δ R (X; R)_pα) performing a second order taylor series expansion, resulting in:

$\begin{array}{c}\mathrm{\&Delta;R}(X;r_p,\mathrm{\&alpha;})=R(X;r_p,\mathrm{\&alpha;})-R(X;r_p,{\mathrm{\&alpha;}}_p)\\ \&ap;-(\mathrm{\&alpha;}-{\mathrm{\&alpha;}}_p)X-\frac{{\mathrm{\&alpha;}}^2-{\mathrm{\&alpha;}}_p^2}{2r_p}{X}^2\end{array}$

if the above formula can ignore X²Item, then I^(u)(r_pα) can be simplified to:

$I^{\left(u\right)}(r_p,\mathrm{\&alpha;})\&ap;{\&Integral;}_{x_u-l/2}^{x_u+l/2}\exp [-{\mathrm{jK}}_{\mathrm{rc}}\&CenterDot;(\mathrm{\&alpha;}-{\mathrm{\&alpha;}}_p)X]\mathrm{dX}$

then, the angular wave number domain variables are defined as follows

$\left\{\begin{array}{c}{K}_{\mathrm{\&alpha;}}=K_{\mathrm{rc}}\&CenterDot;X\\ \mathrm{\&Delta;}{K}_{\mathrm{\&alpha;}}=K_{\mathrm{rc}}\&CenterDot;l\\ K_{\mathrm{\&alpha;c}}^{\left(u\right)}=K_{\mathrm{rc}}\&CenterDot;x_u\end{array}\right.$

Wherein K_αIs the angular wave number, Δ K_αIn order to be a bandwidth of the angular wave number,is the wave number center corresponding to the u-th sub-aperture. Let I^(u)(α)＝I^(u)(r_pα), mixing K_αSubstitution into I^(u)(α), then I^(u)(α) can be further rewritten as to the angular wavenumber K_αFunction of (c):

$I^{\left(u\right)}\left(\mathrm{\&alpha;}\right)=I^{\left(u\right)}(r_p,\mathrm{\&alpha;})\&ap;{\&Integral;}_{K_{\mathrm{\&alpha;c}}^{\left(u\right)}-\mathrm{\&Delta;}{K}_{\mathrm{\&alpha;}}/2}^{K_{\mathrm{\&alpha;c}}^{\left(u\right)}+\mathrm{\&Delta;}{K}_{\mathrm{\&alpha;}}/2}\exp [-j{K}_{\mathrm{\&alpha;}}\&CenterDot;(\mathrm{\&alpha;}-{\mathrm{\&alpha;}}_p)]{\mathrm{dK}}_{\mathrm{\&alpha;}}$

wherein,which is the sine of the view angle from the uth sub-aperture center to the scene center, which can be directly calculated from the geometric relationship, and Δ α is the angular imaging range of each sub-aperture (the angular imaging range of each sub-aperture is equal).By introducing the concept of angular wavenumber, the above equation establishes K_αAnd α, while producing a sinc-like azimuth impulse response function I^(u)(α) it is worth noting that the invention utilizes the angular sine coordinates to establish the Fourier transform relationship between the backward projection image and the angular wave number domain, breaks through the barrier of the wave number spectrum operation under the time domain algorithm, and injects some frequency domain algorithm characteristics into the invention.

Based on Δ R (X; R)_pα) can ignore X²Assumption of terms, we focus below on deriving neglect X²For convenience of analysis, let α be α_p+ Δ ρ, then Δ R (X; R)_pα) of X²The items can be simplified into

$-\frac{{\mathrm{\&alpha;}}^2-{\mathrm{\&alpha;}}_p^2}{2r_p}{X}^2\&ap;-\frac{{\mathrm{\&alpha;}}_p\&CenterDot;\mathrm{\&Delta;\&rho;}}{r_p}{X}^2$

The quadratic phase error QPE corresponding to the above equation is

$\mathrm{QPE}=-\frac{4\mathrm{\&pi;}{\mathrm{\&alpha;}}_p\&CenterDot;\mathrm{\&Delta;\&rho;}}{\mathrm{\&lambda;}{r}_p}{X}^2$

It is generally believed that when QPE ≦ π/8, the energy of the impulse response function is primarily concentrated at α_pCentered, width angular resolution ρ_αThe size of the adjacent area. Angular resolution ρ according to the beam forming principle_αComprises the following steps:

${\mathrm{\&rho;}}_{\mathrm{\&alpha;}}=\frac{{\mathrm{\&lambda;}}_{\min }}{2L}$

wherein λ is_minThe wavelength corresponding to the maximum frequency at which the synthetic aperture radar transmits signals (i.e. the minimum value of the wavelength of the synthetic aperture radar transmitted signals).

Because QPE takes the maximum value when the antenna is positioned at two ends of the synthetic aperture, replacing Delta rho with rho_αThen QPE satisfies:

$\mathrm{QPE}\&le;\frac{\mathrm{\&pi;L}}{2r_p}\&CenterDot;\left|{\mathrm{\&alpha;}}_p\right|\&le;\frac{\mathrm{\&pi;L}}{2r_p}\&le;\frac{\mathrm{\&pi;}}{8}$

further simplifying the above formula, then:

$L\&le;\frac{r_p}{4}$

i.e. the synthetic aperture length is not more than r_pIs one fourth of (a), Δ R (X; R)_pα) about X is sufficient to guarantee that the approximation accuracy of the first order term：I^(u)(r_pα) can be further rewritten as to the angular wavenumber K_αAs a function of (c). And this constraint is satisfied in most imaging scenarios.

S3: defining a distance wavenumber variable K_rAnd the angular wavenumber variable K_α(ii) a By the pair I^(u)(α) performing inverse Fourier transform to obtain a two-dimensional wavenumber spectrum H corresponding to the u-th sub-aperture^(u)(K_r,K_α) (ii) a Splicing the two-dimensional wave number spectrums corresponding to the sub-apertures in the two-dimensional wave number spectrums corresponding to the sub-apertures along the azimuth direction aiming at the two-dimensional wave number spectrums corresponding to the sub-apertures to obtain one-dimensional wave number spectrums corresponding to the sub-apertures; then, the one-dimensional wave number spectrum corresponding to each distance wave number is spliced along the distance direction to obtain a full-aperture two-dimensional wave number spectrum; and performing two-dimensional Fourier transform on the full-aperture two-dimensional wave number spectrum to obtain a full-space resolution image under a polar coordinate system. The concrete description is as follows:

first to I^(u)(α) performing inverse Fourier transform to obtain corresponding angular wavenumber spectrum H^(u)(K_α)，H^(u)(K_α) Comprises the following steps:

$\begin{array}{c}{H}^{\left(u\right)}\left(K_{\mathrm{\&alpha;}}\right)={\&Integral;}_{{\mathrm{\&alpha;}}_c^{\left(u\right)}-\mathrm{\&Delta;\&alpha;}/2}^{{\mathrm{\&alpha;}}_c^{\left(u\right)}+\mathrm{\&Delta;\&alpha;}/2}{I}^{\left(u\right)}\left(\mathrm{\&alpha;}\right)\&CenterDot;\exp \left(j{K}_{\mathrm{\&alpha;}}\mathrm{\&alpha;}\right)\mathrm{d\&alpha;}\\ =\mathrm{rect}\left[\frac{K_{\mathrm{\&alpha;}}-K_{\mathrm{\&alpha;c}}^{\left(u\right)}}{\mathrm{\&Delta;}{K}_{\mathrm{\&alpha;}}}\right]\&CenterDot;\exp (-j{K}_{\mathrm{\&alpha;}}{\mathrm{\&alpha;}}_p)\end{array}$

wherein, the rectangular window function represents the width and position of the angular wave number spectrum under the u-th sub-aperture, and the exponential term corresponds to the angular sine position α of the P point under the polar coordinate system_p。

Referring to fig. 3, a schematic diagram of image fusion realized by one-dimensional angle wave number spectrum splicing of the present invention is shown. Order toThe angular resolution unit size corresponding to the u-th sub-apertureOrWhen the spectrum is folded, the angular wavenumber spectrum is folded. However, this phenomenon does not present insurmountable difficulties for the inventive corner-wavenumber spectrum stitching as long as the corner-domain sampling satisfies the Nyquist sampling theorem. According to the Nyquist sampling theorem,is selected to meet

${\mathrm{\&rho;}}_{\mathrm{\&alpha;}}^{\left(u\right)}\&le;\frac{{\mathrm{\&lambda;}}_{\min }}{2l}$

According to the foregoing description, the synthetic aperture is first divided into N₀A sub-aperture. Each sub-aperture corresponds to a sub-image and a separate angular wavenumber spectrum, wherein the wavenumber center and the wavenumber width of the angular wavenumber spectrum correspond to the sub-aperture center and the sub-aperture length, respectively. Under a polar coordinate system, ordered combination of the angular wave number spectrum, namely image fusion, can be realized through a series of simple wave number spectrum shifts. When the angular wave number spectrum is spliced, the center of the angular wave number spectrum of the u-th sub-aperture needs to be shifted from zero frequency to zero frequency through cyclic shift. The image domain is multiplied by the linear phase to realize the shift of the angular wave number spectrum by considering the Fourier transform characteristic of the signal. The shifting of the angular wavenumber spectrum can eliminate the folding phenomenon of the spectrum. And splicing the aperture wave number spectrums into one-dimensional vectors along the azimuth direction, thereby completing the ordered recombination of the wave number spectrums.

One-dimensional wave number spectrum H^(u)(K_α) Generalizing to the two-dimensional wavenumber spectrum, the two-dimensional wavenumber variables are first defined as follows:

$\left\{\begin{array}{cc}{K}_{\mathrm{\&alpha;}}=K_r\&CenterDot;X& \\ \mathrm{\&Delta;}{K}_{\mathrm{\&alpha;}}=K_r\&CenterDot;l& K_r\&Element;[\frac{4\mathrm{\&pi;}}{{\mathrm{\&lambda;}}_{\max }},\frac{4\mathrm{\&pi;}}{{\mathrm{\&lambda;}}_{\min }}]\\ K_{\mathrm{\&alpha;c}}^{\left(u\right)}=K_r\&CenterDot;x_u& \end{array}\right.$

wherein λ is_maxFor the wavelength, K, corresponding to the minimum frequency at which the synthetic aperture radar transmits signals_rIs the distance wave number, K_rAnd 4 pi/lambda, and lambda is the wavelength (lambda is a variable) when the synthetic aperture radar transmits signals. The two-dimensional wavenumber spectrum H corresponding to the u-th sub-aperture^(u)(K_r,K_α) Comprises the following steps:

$H^{\left(u\right)}(K_r,K_{\mathrm{\&alpha;}})=\mathrm{rect}\left[\frac{K_r-K_{\mathrm{rc}}}{\mathrm{\&Delta;}{K}_r}\right]\&CenterDot;\exp (-j\&CenterDot;K_r{r}_p)\&CenterDot;\mathrm{rect}\left[\frac{K_{\mathrm{\&alpha;}}-K_{\mathrm{\&alpha;c}}^{\left(u\right)}}{\mathrm{\&Delta;}{K}_{\mathrm{\&alpha;}}}\right]\&CenterDot;\exp (-j\&CenterDot;K_{\mathrm{\&alpha;}}{\mathrm{\&alpha;}}_p)$

$\left\{\begin{array}{c}{K}_{\mathrm{\&alpha;}}\&Element;[-\frac{\mathrm{\&Delta;}{K}_{\mathrm{\&alpha;}}}{2}+K_{\mathrm{\&alpha;c}}^{\left(u\right)},\frac{\mathrm{\&Delta;}{K}_{\mathrm{\&alpha;}}}{2}+K_{\mathrm{\&alpha;c}}^{\left(u\right)}]\\ K_r\&Element;[-\frac{\mathrm{\&Delta;}{K}_r}{2}+K_{\mathrm{rc}},\frac{\mathrm{\&Delta;}{K}_r}{2}+K_{\mathrm{rc}}]\end{array}\right.$

wherein Δ K_r4 π B/c, c for the speed of light, Δ K_rRepresenting the distance wavenumber width.

Under the backward projection framework, step S2 is performed by ignoring Δ R (X; R)_pα) X²Constraint of itemsUnder the condition, the Fourier transform relation between the image domain and the angular wave number domain is established, so that the analysis of the two-dimensional wave number spectrum expresses H^(u)(K_r,K_α) Is an approximate, resolved rectangular wavenumber spectrum. It should be noted that neither fourier transform nor cyclic shift in the actual processing procedure introduces an approximation, and thus the two-dimensional wavenumber spectrum after stitching is an exact spectrum with a trapezoidal shape, as shown in fig. 4. Referring to FIG. 4, the two-dimensional spectrum splicing of the present invention is a schematic diagram of image fusion, the angular wavenumber centerIs the distance wave number K_rIs a linear function of (a). By referring to the angular domain wave number spectrum splicing method shown in fig. 3, the two-dimensional wave number spectrum splicing method is shown in fig. 4. Specifically, a two-dimensional wave number spectrum H corresponding to the u-th sub-aperture is obtained^(u)(K_r,K_α) And then, splicing the two-dimensional wave number spectrums corresponding to the sub-apertures in the two-dimensional wave number spectrums corresponding to all the sub-apertures along the azimuth direction to obtain one-dimensional wave number spectrums corresponding to all the sub-apertures (expressed in the form of one-dimensional vectors). Then, the one-dimensional wave number spectrum corresponding to each distance wave number is spliced along the distance direction to obtain a full-aperture two-dimensional wave number spectrum; and performing two-dimensional Fourier transform on the full-aperture two-dimensional wave number spectrum to obtain a full-space resolution image (angular domain full resolution image) under a polar coordinate system.

The operation complexity of step S2 and step S3 is analyzed as follows:

assuming that the number of aperture positions (the number of azimuth pulses and the number of azimuth sampling points) in the full aperture is M, and the number of pixel points of the imaging grid is M × M, the number of aperture positions per sub-aperture is M, where M is M/N₀In the two-dimensional wavenumber spectrum stitching stage of step S3, the cyclic shift operation only occupies a very small portion of the overall computation, so that the portion is ignored.

a) The amount of computation required for step S2 is (M/N)₀×M/N₀×M)×N₀；

b) In step S3, the calculation amount required for obtaining the two-dimensional wave number spectrum corresponding to each sub-aperture is: $[M\&times;M/N_0\&times;({\log }_2^{M/N_0}+{\log }_2^M)]\&times;N_0;$

c) in step S3, the operation amount required for obtaining the full aperture two-dimensional wave number spectrum is M × N₀The method comprises the steps of firstly creating a full-aperture wave number spectrum matrix (memory space) when obtaining a full-aperture two-dimensional wave number spectrum, calculating an angular wave number center of the wave number spectrum, realizing ordered splicing of the sub-aperture wave number spectrum along the azimuth direction through cyclic shift, repeating the process along the distance direction, finally realizing splicing of the two-dimensional full-aperture wave number spectrum, obtaining the full-aperture two-dimensional wave number spectrum, and if the operation amount of cyclic shift operation is ignored, calculating the operation amount required by the angular wave number center to be M × N₀。

d) In step S3, the amount of computation required for performing the two-dimensional fourier transform on the full-aperture two-dimensional wave number spectrum is

In summary, the total operation amount OP of step S2 and step S3_AFBPComprises the following steps:

${\mathrm{OP}}_{\mathrm{AFBP}}=M^3/N_0+{3M}^2{\log }_2^M+M^2{\log }_2^{M/N_0}+{\mathrm{MN}}_0$

observe the above formula, OP_AFBPIs about N₀Is a monotonically decreasing function of (a). Therefore, when N is₀The larger the calculation amount of step S2 and step S3, the smaller the total calculation amount, the higher the calculation efficiency. When the length of the sub-aperture is very small, such as m 4 or m 8, the angular wavenumber spectrum bandwidth is very narrow, blurring of the full-aperture wavenumber spectrum is easily caused at the time of stitching of the wavenumber spectrum, and degradation of image quality is caused. When in actual data processing, the wave number spectrum can be ensured to have good precision by taking m as 16 or m as 32, and the calculation efficiency of the algorithm can be considered.

After step S3, converting the full-space-resolved image in the polar coordinate system into the planar rectangular coordinate system by two-dimensional interpolation, so as to obtain the full-space-resolved image in the planar rectangular coordinate system. And then displaying the full-space resolution image under the plane rectangular coordinate system on a display.

The effects of the invention can be further illustrated by the following simulation experiments and actual measurement data experiments:

simulation conditions are as follows:

the simulation experiment uses an X-band radar, and the parameters of the simulation experiment are shown in the table I.

table-X band radar parameters

Center frequency	5GHz	Center slope distance	1000m
				Two dimensional resolution	0.175m	Synthetic pore size length	192m

9 point targets are equally spaced within the 100m × 100m imaging region, as shown in fig. 5. Fig. 5 is a geometric schematic diagram of the distribution of the imaging region and the point target in the simulation experiment of the present invention. One synthetic aperture contains 2048 columns of pulse echoes and the pixel distribution of the imaging grid is 2048 x 2048.

The main parameters of the computer simulation platform are as follows: a 2.53GHz processor, 4GB memory (RAM), a 64-bit operating system, and 64-bit MATLAB.7.10.0 software.

Simulation content and simulation result analysis:

in the simulation experiment, the present invention, a Back-Projection (BP) method, and a fast decomposition Back-Projection (FFBP) method were used to perform imaging simulation, respectively. When the imaging simulation is carried out by adopting a fast decomposition back projection method, a local polar coordinate system subimage is reconstructed through the sub-aperture BP integral, but in the invention, the subimage is reconstructed in a global polar coordinate system. The invention and the fast decomposition back projection method only need to carry out distance interpolation at the corresponding stage, and do not need to carry out angular domain interpolation. In order to ensure the fairness and the comparability of a simulation experiment, the method and the fast decomposition back projection method are realized by up-sampling when distance interpolation is carried out, and the method comprises the following basic steps: zero filling at two ends of a distance frequency domain, inverse Fourier transform and low-order interpolation, wherein the up-sampling multiple is 8 to ensure the interpolation precision.

When obtaining a full-space resolution image under a polar coordinate system, the fast decomposition back projection method realizes the recursive fusion of the image by means of two-dimensional distance interpolation, wherein the distance interpolation is still realized by 8 times of upsampling, and the angular domain interpolation is realized by truncated sinc interpolation. And obtaining a full-space resolution image under a polar coordinate system after the recursive fusion is finished. When a full-space resolution image under a polar coordinate system is obtained, the full-aperture wave number spectrum is recombined through Fourier transform and cyclic shift, and interpolation processing is not needed in the process. And performing two-dimensional inverse Fourier transform on the full-aperture wave number spectrum to obtain a full-space resolution image under a polar coordinate system.

For the method and the fast decomposition back projection method, different sub-aperture lengths are set (at the moment, the number m of the aperture positions of each sub-aperture is respectively set to be 32, 16 and 8), three groups of different imaging results are obtained through three different experiments, and the reconstructed images are all full-space resolution images under a polar coordinate system. In order to illustrate the image focusing quality of the present invention and the fast decomposition back projection method, three point targets of P1, P2 and P3 in fig. 5 were selected for image quality analysis. Limited to space limitations and emphasis, only the azimuthal impulse response function of the three point targets at different sub-aperture lengths l is given here, as shown in fig. 6. Referring to fig. 6a, an azimuth impulse response function diagram of three points P1, P2 and P3 is reconstructed by using three methods when the number of aperture positions of each sub-aperture in the simulation experiment is 32. Referring to fig. 6b, an azimuth impulse response function diagram of three points P1, P2 and P3 is reconstructed by using three methods when the number of aperture positions of each sub-aperture in the simulation experiment is 16. Referring to fig. 6c, an azimuth impulse response function diagram of three points P1, P2 and P3 is reconstructed by using three methods when the number of aperture positions of each sub-aperture in the simulation experiment is 8. In fig. 6a, 6b and 6c, the horizontal axis represents angular domain interpolation position and the vertical axis represents normalized amplitude in dB.

It can be seen from fig. 6 that even though the values of the sub-aperture length l are different, the present invention and the fast decomposition back projection method can both achieve good focusing performance. The quality of the reconstructed image of the present invention and the fast decomposition backprojection method improves as the sub-aperture length l increases. When the length L of the sub-aperture is equal to the length L of the synthetic aperture, namely when the sub-aperture is not divided, the quality of the reconstructed image of the method disclosed by the invention and the method for fast decomposing the back projection is equivalent to a time domain back projection method. And as the subaperture length is decreased, the notch depth of the main lobe of the IRF (impulse response function) in the fast decomposition back projection method is continuously changed to be shallow, the side lobe of the IRF is continuously increased, and the response form of the side lobe of the IRF in the fourth and the later stages is not ideal. When the aperture length is equal to the initial stage aperture length of the fast decomposition back projection method, the impulse response function curve (the azimuth profile of the three points P1, P2 and P3) of the invention is closer to the time domain back projection method. At this time, the zero depth of the main lobe and other positions of the invention is lower than that of the fast decomposition back projection method. As shown in fig. 6a, when the number of initial sub-aperture positions is 32, the impulse response function curves of the three points P1, P2 and P3 obtained by the present invention are substantially consistent with the time domain back projection method, and there is no undesirable situation of high order side lobe (IRF side lobe) when the far-field point P3 is focused by the fast decomposition back projection method. From the above analysis of fig. 6, the present invention intuitively has a more desirable focusing performance than the fast decomposition backprojection method. In order to provide quantitative evaluation of image quality and computational efficiency, parameters such as corresponding Peak Side Lobe Ratio (PSLR), Integrated Side Lobe Ratio (ISLR), processing time, and the like are obtained as shown in table two. As seen from the table II, the peak-to-side lobe ratio performance of the method is slightly better than that of the fast decomposition back projection method, but the method has obvious improvement on the integral-to-side lobe ratio. In the case of close peak sidelobe performance, the lower the integral sidelobe ratio, the higher the contrast of the image, and this phenomenon will be better exemplified in the processing of the measured data. From the aspect of processing time, under the condition that the sub-aperture lengths are equal, the method has higher operation efficiency than a fast decomposition back projection method. The aperture decomposition brings good operation efficiency for the fast decomposition back projection method, but the two-dimensional interpolation realizes recursive fusion and brings extra operation burden to the algorithm. For the fast decomposition backprojection method, the shorter the sub-aperture length, the larger the amount of computation required for recursive fusion. Corresponding to the above three different sub-aperture lengths, the time consumed by the present invention in performing step S2 is 73.8S, 36.6S, and 18.7S, respectively. As can be seen from the total processing time in the second comparison table, the calculation amount of the present invention is mostly concentrated on the step S2, and the time consumption of the spectrum stitching stage is very small. In summary, the invention and the fast decomposition back projection method both greatly improve the operation efficiency of the BP integration, but the invention is superior to the fast decomposition back projection method in terms of imaging accuracy and operation efficiency.

TABLE II image quality evaluation and processing time comparison

Actually measured data experiment:

the performance of the invention will be verified by X-band high resolution beamforming SAR imaging processing. In the actual measurement data experiment, 16384 pulse echoes are contained in the corresponding echo signal, and the synthetic aperture time is 7.8 seconds. The transmission signal bandwidth is 1.16GHz, and the center slant distance is 10.5 Km. The number of distance sampling points used in data processing is 16384, and the imaging range is 1.6 × 1.2 km. Considering the influence of two-dimensional windowing on the broadening of the main lobe of the IRF, the theoretical values of the distance and azimuth resolutions are about 0.15 m.

In order to explain the difference of image fusion realized by two-dimensional interpolation processing and two-dimensional wave number spectrum splicing and avoid interpolation errors introduced by subsequent coordinate system transformation, the given processing results are all located in a polar coordinate system, and the reconstructed image is a full-space resolution image in the polar coordinate system. In the actual measurement data experiment, the image reconstruction is respectively carried out by adopting the method and the fast decomposition back projection method, and the obtained images are respectively shown in fig. 7a and fig. 7 b. Referring to fig. 7a, a full-space resolution image in a polar coordinate system is reconstructed by using a fast decomposition back projection method in a measured data test. Referring to fig. 7b, the full-space resolution image in the polar coordinate system reconstructed by the present invention in the measured data test is shown.

It can be seen from fig. 7a and 7b that the accuracy of the fast decomposition backprojection method is inherited, and both the fast decomposition backprojection method and the present invention achieve good focusing of the image. During data processing, 1024 sub-apertures are divided, each sub-aperture containing 16 aperture locations. The fast decomposition back projection method needs 10 times of recursion processing to realize image fusion, and the operation time is 5751.4 s; the invention only needs one-time spectrum splicing to realize image fusion, the operation time is 1723.4s, and the operation time for realizing image fusion is shorter than that of the fast decomposition back projection method.

As can be seen from fig. 7a and 7b, both the fast decomposition backprojection method and the present invention achieve good focusing of the image. However, from the imaging of field texture and dense building clusters in the lower part of fig. 7, the present invention contributes slightly better to image contrast than the fast decomposition back projection method. In order to evaluate the image quality of the present invention and the fast decomposition backprojection method in only one step, three different sub-scenes (indicated as A, B and C in fig. 7a and 7b, respectively) extracted from fig. 7a and 7b are subjected to an enlargement process to obtain three sub-scene images, as shown in fig. 8. Referring to fig. 8a, the full-space resolution image under the polar coordinate system of the sub-scene a obtained by adopting the fast decomposition back projection method in the actual measurement data test is shown, and referring to fig. 8b, the full-space resolution image under the polar coordinate system of the sub-scene a obtained by adopting the method in the actual measurement data test is shown; referring to fig. 8c, the full-space resolution image under the polar coordinate system of the sub-scene B obtained by adopting the fast decomposition back projection method in the actual measurement data test is referred to fig. 8d, the full-space resolution image under the polar coordinate system of the sub-scene B obtained by adopting the method in the actual measurement data test is referred to; referring to fig. 8e, the full-space resolution image of the sub-scene C in the polar coordinate system obtained by the fast decomposition back projection method in the measured data test is shown, and referring to fig. 8f, the full-space resolution image of the sub-scene C in the polar coordinate system obtained by the fast decomposition back projection method in the measured data test is shown. From fig. 8a to fig. 8f, it is seen that the point target is well focused by the present invention and the fast decomposition back projection method, but the focused point target has more concentrated energy, lower IRF side lobe level and better image contrast, and this phenomenon is attributed to the present invention that reasonably avoids the influence of interpolation error and its accumulation effect on the wave number spectrum precision. In sub-scene a, there are two azimuthally adjacent 0.15m apart corner reflectors (labeled point target 1 and point target 2 in fig. 8 b) and an independently placed corner reflector (point target 3), which are identified by a circular ring in fig. 8 a. To illustrate the resolution performance of the present invention, point targets 1 and 2 were imaged and analyzed for the azimuthal impulse response function, as shown in FIG. 9. Referring to fig. 9a, a schematic diagram of an imaging result after two-dimensional interpolation obtained by using a fast decomposition back projection method for two point targets in an actually measured data experiment is shown. Referring to fig. 9b, a schematic diagram of an imaging result after two-dimensional interpolation is obtained for two point targets in an actually measured data experiment by using the method of the invention. Referring to fig. 9c, schematic diagrams of azimuth impulse response functions obtained by two methods for two point targets in the measured data experiment are shown. In fig. 9a and 9b, the horizontal axis represents the angular interpolation position, and the vertical axis represents the distance interpolation position. In fig. 9c, the horizontal axis represents angular domain interpolation position and the vertical axis represents normalized amplitude. AFBP refers to the present invention, and FFBP refers to the fast decomposition backprojection method. As can be seen in fig. 9c, the IRF main lobes for the two targets have been separated, indicating that adjacently placed corner reflectors can be well resolved.

To illustrate the focusing performance of the present invention, the point target 3 is imaged and analyzed for the azimuth impulse response function, as shown in fig. 10. Referring to fig. 10a, a schematic diagram of a two-dimensional interpolated imaging result obtained by a fast decomposition back projection method for a point target 3 in an actual measurement data experiment is shown. Referring to fig. 10b, a schematic diagram of an imaging result of the point target 3 after two-dimensional interpolation obtained by the present invention in an experiment of measured data is shown. Referring to fig. 10c, the point target 3 is shown in the experimental data by two methods to obtain the azimuth impulse response function. In fig. 10a and 10b, the horizontal axis represents the angular interpolation position, and the vertical axis represents the distance interpolation position. In fig. 10c, the horizontal axis represents angular domain interpolation position and the vertical axis represents normalized amplitude. AFBP refers to the present invention, and FFBP refers to the fast decomposition backprojection method. As can be seen from FIG. 10c, the IRF first minor lobe of the present invention is close to-25 dB, whereas the first minor lobe on the left side of the IRF is about-15 dB in the fast decomposition backprojection method, but is already absorbed by the IRF main lobe. It can be found from fig. 9c and fig. 10c that the IRF side lobe performance of the present invention is still better than that of the fast decomposition back projection method when processing the measured data, which is consistent with the conclusion obtained by the simulation experiment. Of course, the fast decomposition back projection method can improve the side lobe performance of the azimuth impulse response function by properly increasing the aperture length of the initial stage, and the method necessarily increases the processing time of the algorithm. Compared with a fast decomposition back projection method, the method can obtain better image quality and higher operation efficiency under the condition of relatively shorter sub-aperture.

It will be apparent to those skilled in the art that various changes and modifications may be made in the present invention without departing from the spirit and scope of the invention. Thus, if such modifications and variations of the present invention fall within the scope of the claims of the present invention and their equivalents, the present invention is also intended to include such modifications and variations.

Claims (6)

1. An accelerated decomposition backward projection bunching synthetic aperture radar imaging method is characterized by comprising the following steps:

s1: receiving an echo signal by using a synthetic aperture radar on a motion platform, and performing matched filtering on the echo signal to obtain a range pulse compression signal, wherein the motion platform is an airplane or a satellite;

s2: establishing a polar coordinate system with the synthetic aperture center O as an origin, wherein the coordinate is (r)_pTheta) is represented by (r)_pTheta) pixel point, dividing the full aperture into N₀Sub-apertures of equal length, N₀Setting α to sin theta, and making the u-th sub-aperture correspond to (r)_pTheta) the impulse response function of the pixel point is represented as I^(u)(r_pα), u is 1 to N₀(ii) a Let I^(u)(α)＝I^(u)(r_p,α)；

S3: defining a distance wavenumber variable K_rAnd the angular wavenumber variable K_αThe following were used:

$\left\{\begin{array}{c}{K}_{\mathrm{\&alpha;}}=K_r\&CenterDot;X\\ {\mathrm{\&Delta;K}}_{\mathrm{\&alpha;}}=K_r\&CenterDot;l\\ K_{\mathrm{\&alpha;}c}^{\left(u\right)}=K_r\&CenterDot;x_u\end{array}\right.,K_r\&Element;\&lsqb;\frac{4\mathrm{\&pi;}}{{\mathrm{\&lambda;}}_{max}},\frac{4\mathrm{\&pi;}}{{\mathrm{\&lambda;}}_{min}}\&rsqb;$

wherein λ is_maxFor the wavelength, λ, corresponding to the minimum frequency at which the synthetic aperture radar transmits signals_minThe wavelength corresponding to the maximum frequency when the synthetic aperture radar transmits signals; k_r4 pi/lambda, wherein lambda is the wavelength of the synthetic aperture radar when transmitting signals; full aperture division into N₀Sub-apertures of equal length, N₀Is a natural number greater than 1, and the length of each sub-aperture is l, x_u＝(u-N₀L, u is 1 to N₀(ii) a The X-axis direction is the motion direction of the motion platform, and the absolute value of the X-axis coordinate value at the antenna phase center is X;

by the pair I^(u)(α) performing inverse Fourier transform to obtain a two-dimensional wavenumber spectrum H corresponding to the u-th sub-aperture^(u)(K_r,K_α)；

Splicing the two-dimensional wave number spectrums corresponding to the sub-apertures in the two-dimensional wave number spectrums corresponding to the sub-apertures along the azimuth direction aiming at the two-dimensional wave number spectrums corresponding to the sub-apertures to obtain one-dimensional wave number spectrums corresponding to the sub-apertures; then, the one-dimensional wave number spectrum corresponding to each distance wave number is spliced along the distance direction to obtain a full-aperture two-dimensional wave number spectrum; and performing two-dimensional Fourier transform on the full-aperture two-dimensional wave number spectrum to obtain a full-space resolution image under a polar coordinate system.

2. The method as claimed in claim 1, wherein in step S1, the moving platform forms a synthetic aperture with a length L when flying at a constant speed, and the center of the synthetic aperture is denoted as O; in the flying process of the moving platform, the synthetic aperture radar transmits signals outwards and receives echo signals to obtain demodulated echo signals;

in step S1, the process of obtaining the demodulated echo signal is: establishing a plane rectangular coordinate system by taking the synthetic aperture center O as an origin, wherein in the plane rectangular coordinate system, the X-axis direction is the motion direction of the motion platform, and the absolute value of the X-axis coordinate value at the antenna phase center is X; then, a polar coordinate system is established by taking the center O of the synthetic aperture as the origin, and the phase center of the antenna reaches a point target P (r)_p,θ_p) Is determined by the instantaneous slope distance R (X; r is_p,θ_p) Comprises the following steps:

$R(X;r_p,{\mathrm{\&theta;}}_p)=\sqrt{r_p^2+X^2-2r_{p}X{\mathrm{sin\&theta;}}_p},-\frac{L}{2}\&le;X<\frac{L}{2}$

wherein r is_pAnd theta_pTwo coordinates of the point target P under a polar coordinate system respectively, and then an angular sine coordinate α is defined_p，α_p＝sinθ_pThen, the above formula is rewritten as:

$R(X;r_p,{\mathrm{\&alpha;}}_p)=\sqrt{r_p^2+X^2-2r_p{\mathrm{X\&alpha;}}_p},-\frac{L}{2}\&le;X<\frac{L}{2}$

the demodulated echo signal s (τ, X) is then:

$s(\mathrm{\&tau;},X)=rect\&lsqb;\frac{X}{L}\&rsqb;\&CenterDot;s_T(\mathrm{\&tau;}-{\mathrm{\&Delta;t}}_p)$

wherein s is_T(. to) is the emission signal waveform, τ is the fast time, △ t_p＝2R(X；r_p,α_p) C, c is the speed of light, rect [. C]Representing a rectangular window function;

in step S1, after the demodulated echo signal is obtained, the demodulated echo signal is subjected to matched filtering to obtain a range pulse compressed signal S_M(τ,X)：

$s_M(\mathrm{\&tau;},X)=\sin c\&lsqb;B\&CenterDot;(\mathrm{\&tau;}-{\mathrm{\&Delta;t}}_p)\&rsqb;\&CenterDot;rect\&lsqb;\frac{X}{L}\&rsqb;\&CenterDot;\exp \&lsqb;-{\mathrm{jK}}_{rc}R(X;r_p,{\mathrm{\&alpha;}}_p)\&rsqb;$

Wherein B is the transmission signal bandwidth of the synthetic aperture radar, K_rc＝4π/λ_mid，λ_midThe wavelength corresponding to the carrier frequency when the synthetic aperture radar transmits signals.

3. The method as claimed in claim 1, wherein in step S1, the moving platform forms a synthetic aperture with a length L when flying at a constant speed, and the center of the synthetic aperture is denoted as O; in the flying process of the moving platform, the synthetic aperture radar transmits signals outwards and receives echo signals to obtain demodulated echo signals;

in step S1, the process of obtaining the demodulated echo signal is: establishing a plane rectangular coordinate system by taking the synthetic aperture center O as an origin, wherein in the plane rectangular coordinate system, the X-axis direction is the motion direction of the motion platform, and the absolute value of the X-axis coordinate value at the antenna phase center is X; then, a polar coordinate system is established by taking the center O of the synthetic aperture as the origin, and the phase center of the antenna reaches a point target P (r)_p,θ_p) Is determined by the instantaneous slope distance R (X; r is_p,θ_p) Comprises the following steps:

$R(X;r_p,{\mathrm{\&theta;}}_p)=\sqrt{r_p^2+X^2-2r_{p}X{\mathrm{sin\&theta;}}_p},-\frac{L}{2}\&le;X<\frac{L}{2}$

wherein r is_pAnd theta_pTwo coordinates of the point target P under a polar coordinate system respectively, and then an angular sine coordinate α is defined_p，α_p＝sinθ_pThen, the above formula is rewritten as:

$R(X;r_p,{\mathrm{\&alpha;}}_p)=\sqrt{r_p^2+X^2-2r_p{\mathrm{X\&alpha;}}_p},-\frac{L}{2}\&le;X<\frac{L}{2}$

the demodulated echo signal s (τ, X) is then:

$s(\mathrm{\&tau;},X)=rect\&lsqb;\frac{X}{L}\&rsqb;\&CenterDot;s_T(\mathrm{\&tau;}-{\mathrm{\&Delta;t}}_p)$

wherein s is_T(. to) is the emission signal waveform, τ is the fast time, △ t_p＝2R(X；r_p,α_p) C, c is the speed of light, rect [. C]Representing a rectangular window function;

in step S1, the method obtainsAfter the demodulated echo signals are obtained, calculating the slope distance error caused by the motion error according to platform parameters, wherein the platform parameters comprise the position, the posture, the speed and the acceleration of the motion platform; according to the slant range error caused by the motion error, performing motion compensation on the demodulated echo signal to obtain a motion-compensated echo signal; then, the echo signals after the motion compensation are subjected to matched filtering to obtain a range pulse compression signal s_M(τ,X)：

$s_M(\mathrm{\&tau;},X)=\sin c\&lsqb;B\&CenterDot;(\mathrm{\&tau;}-{\mathrm{\&Delta;t}}_p)\&rsqb;\&CenterDot;rect\&lsqb;\frac{X}{L}\&rsqb;\&CenterDot;\exp \&lsqb;-{\mathrm{jK}}_{rc}R(X;r_p,{\mathrm{\&alpha;}}_p)\&rsqb;$

Wherein B is the transmission signal bandwidth of the synthetic aperture radar, K_rc＝4π/λ_mid，λ_midThe wavelength corresponding to the carrier frequency when the synthetic aperture radar transmits signals.

4. An acceleration component as claimed in claim 3A method for imaging a synthetic aperture radar by backward projection solution and beamforming, wherein in step S2, a full aperture is first divided into N₀Sub-apertures of equal length, N₀Is a natural number greater than 1, and the length of each sub-aperture is l;

in the polar coordinate system, the coordinate is (r)_pTheta) is represented by (r)_pTheta), let α be sin theta, then in the polar coordinate system, the u-th sub-aperture corresponds to (r_pTheta) impulse response function I of the pixel^(u)(r_pα) is:

$I^{\left(u\right)}\left(r_p,\mathrm{\&alpha;}\right)={\&Integral;}_{-l/2}^{l/2}{s}_M\&lsqb;\frac{2R\left(X+x_u;r_p,\mathrm{\&alpha;}\right)}{c},X+x_u\&rsqb;\&CenterDot;\exp \&lsqb;{\mathrm{jK}}_{rc}\&CenterDot;R\left(X+x_u;r_p,\mathrm{\&alpha;}\right)\&rsqb;dX$

wherein x is_u＝(u-N₀L, u is 1 to N₀And c is the speed of light; using X to X + X in the above formula_uIf the variables are replaced, the above formula is rewritten as:

$\begin{array}{c}{I}^{\left(u\right)}\left(r_p,\mathrm{\&alpha;}\right)={\&Integral;}_{x_u-l/2}^{x_u+l/2}{s}_M\&lsqb;\frac{2R\left(X;r_p,\mathrm{\&alpha;}\right)}{c},X\&rsqb;\&CenterDot;\exp \&lsqb;{\mathrm{jK}}_{rc}R\left(X;r_p,\mathrm{\&alpha;}\right)\&rsqb;dX\\ ={\&Integral;}_{x_u-l/2}^{x_u+l/2}\exp \&lsqb;{\mathrm{jK}}_{rc}\&CenterDot;\mathrm{\&Delta;}R\left(X;r_p,\mathrm{\&alpha;}\right)\&rsqb;dX\end{array}$

wherein △ R (X; R)_p,α)＝R(X；r_p,α)-R(X；r_p,α_p) Then pair △ R (X; R)_pα) performing a second order taylor series expansion, resulting in:

$\begin{array}{c}\mathrm{\&Delta;}R(X;r_p,\mathrm{\&alpha;})=R(X;r_p,\mathrm{\&alpha;})-R(X;r_p,{\mathrm{\&alpha;}}_p)\\ =-\left(\mathrm{\&alpha;}-{\mathrm{\&alpha;}}_p\right)X-\frac{{\mathrm{\&alpha;}}^2-{\mathrm{\&alpha;}}_p^2}{2r_p}{X}^2\end{array}$

ignoring X in the above equation²Item, then I^(u)(r_pα) to:

$I^{\left(u\right)}\left(r_p,\mathrm{\&alpha;}\right)={\&Integral;}_{x_u-l/2}^{x_u+l/2}\exp \&lsqb;-{\mathrm{jK}}_{rc}\&CenterDot;(\mathrm{\&alpha;}-{\mathrm{\&alpha;}}_p)X\&rsqb;dX$

then, the angular wave number domain variables are defined as follows:

$\left\{\begin{array}{c}{K}_{\mathrm{\&alpha;}}=K_{rc}\&CenterDot;X\\ {\mathrm{\&Delta;K}}_{\mathrm{\&alpha;}}=K_{rc}\&CenterDot;l\\ K_{\mathrm{\&alpha;}c}^{\left(u\right)}=K_{rc}\&CenterDot;x_u\end{array}\right.$

let I^(u)(α)＝I^(u)(r_pα), mixing K_αSubstitution into I^(u)(α), then I^(u)(α) is:

$I^{\left(u\right)}\left(\mathrm{\&alpha;}\right)=I^{\left(u\right)}(r_p,\mathrm{\&alpha;})={\&Integral;}_{K_{\mathrm{\&alpha;}c}^{\left(u\right)}-{\mathrm{\&Delta;K}}_{\mathrm{\&alpha;}}/2}^{K_{\mathrm{\&alpha;}c}^{\left(u\right)}+{\mathrm{\&Delta;K}}_{\mathrm{\&alpha;}}/2}\exp \&lsqb;-{\mathrm{jK}}_{\mathrm{\&alpha;}}\&CenterDot;\left(\mathrm{\&alpha;}-{\mathrm{\&alpha;}}_p\right)\&rsqb;{\mathrm{dK}}_{\mathrm{\&alpha;}}$

wherein, △α is the angular domain imaging range for each sub-aperture, being the sine of the view angle from the center of the u-th sub-aperture to the center of the scene.

5. The method of accelerated decomposition backward projection beaming synthetic aperture radar imaging of claim 4, wherein in step S3, the two-dimensional wavenumber variables are first defined as follows:

$\left\{\begin{array}{c}{K}_{\mathrm{\&alpha;}}=K_r\&CenterDot;X\\ {\mathrm{\&Delta;K}}_{\mathrm{\&alpha;}}=K_r\&CenterDot;l\\ K_{\mathrm{\&alpha;}c}^{\left(u\right)}=K_r\&CenterDot;x_u\end{array}\right.,K_r\&Element;\&lsqb;\frac{4\mathrm{\&pi;}}{{\mathrm{\&lambda;}}_{max}},\frac{4\mathrm{\&pi;}}{{\mathrm{\&lambda;}}_{min}}\&rsqb;$

wherein λ is_maxFor synthetic aperture radar transmissionWavelength, λ, corresponding to the minimum frequency of the signal_minThe wavelength corresponding to the maximum frequency when the synthetic aperture radar transmits signals; k_r4 pi/lambda, wherein lambda is the wavelength of the synthetic aperture radar when transmitting signals;

by the pair I^(u)(α) performing inverse Fourier transform to obtain a two-dimensional wave number spectrum H corresponding to the u-th sub-aperture^(u)(K_r,K_α)，H^(u)(K_r,K_α) Comprises the following steps:

$\begin{array}{c}{H}^{\left(u\right)}\left(K_{\mathrm{\&alpha;}}\right)={\&Integral;}_{{\mathrm{\&alpha;}}_c^{\left(u\right)}-\mathrm{\&Delta;}\mathrm{\&alpha;}/2}^{{\mathrm{\&alpha;}}_c^{\left(u\right)}+\mathrm{\&Delta;}\mathrm{\&alpha;}/2}{I}^{\left(u\right)}\left(\mathrm{\&alpha;}\right)\&CenterDot;\exp \left({\mathrm{jK}}_{\mathrm{\&alpha;}}\mathrm{\&alpha;}\right)d\mathrm{\&alpha;}\\ =rect\&lsqb;\frac{K_{\mathrm{\&alpha;}}-K_{\mathrm{\&alpha;}c}^{\left(u\right)}}{{\mathrm{\&Delta;K}}_{\mathrm{\&alpha;}}}\&rsqb;\&CenterDot;\exp \left(-{\mathrm{jK}}_{\mathrm{\&alpha;}}{\mathrm{\&alpha;}}_p\right)\end{array}$

$\left\{\begin{array}{c}{K}_{\mathrm{\&alpha;}}\&Element;\&lsqb;-\frac{{\mathrm{\&Delta;K}}_{\mathrm{\&alpha;}}}{2}+K_{\mathrm{\&alpha;}c}^{\left(u\right)},\frac{{\mathrm{\&Delta;K}}_{\mathrm{\&alpha;}}}{2}+K_{\mathrm{\&alpha;}c}^{\left(u\right)}\&rsqb;\\ K_r\&Element;\&lsqb;-\frac{{\mathrm{\&Delta;K}}_r}{2}+K_{rc},\frac{{\mathrm{\&Delta;K}}_r}{2}+K_{rc}\&rsqb;\end{array}\right.$

wherein, △ K_rAnd c is the speed of light 4 pi B/c.

6. The method as claimed in claim 1, wherein after step S3, the full-space-resolved image in the polar coordinate system is transformed into a rectangular planar coordinate system by two-dimensional interpolation to obtain the full-space-resolved image in the rectangular planar coordinate system.