CN109358328A - Polar Coordinate Format Imaging Method for Bistatic Forward Looking SAR on Mobile Platforms - Google Patents

Abstract

The invention discloses a polar coordinate imaging method for a bistatic forward-looking beamforming synthetic aperture radar on a mobile platform. The method of affine transformation of the wavenumber spectrum support domain is adopted to transform the parallelogram wavenumber spectrum support domain of MP-BFSAR into a horizontal quasi-rectangle. , and then realize the resampling of the wavenumber spectrum from polar coordinates to Cartesian coordinates through distance chirp‑z transformation, and then perform distance compression and scene center motion compensation, and then use the method of upsampling mapping accumulation in the azimuth direction to realize the wavenumber spectrum azimuth polar Resampling from coordinates to Cartesian coordinates; then perform two-dimensional inverse Fourier transform to obtain rough imaging results; finally obtain final imaging results through defocus correction; the utilization rate of the collected wavenumber spectrum in the present invention is compared with the existing method is significantly improved.

Description

Polar Coordinate Format Imaging Method for Bistatic Forward Looking SAR on Mobile Platforms

technical field

The invention belongs to the field of bistatic forward-looking synthetic aperture radar (Bistatic Forward-looking Synthetic Aperture Radar, BFSAR) imaging, and particularly relates to a polar coordinate format algorithm of a mobile platform bistatic forward-looking beamforming synthetic aperture radar.

Background technique

Due to the spatial separation of the transmitting station and the receiving station, bistatic SAR breaks through the forward-looking constraint of SAR, enabling the receiving platform to image the forward-looking terrain in the flight direction. This unique feature makes BFSAR a promising sensing technology. BFSAR can provide all-day, all-weather forward-looking high-resolution images.

Due to its bistatic forward-looking configuration, the echoes of MP-BFSAR exhibit severe two-dimensional spatial errors, and its wavenumber spectrum support region is usually a parallelogram rather than a rectangle, making it difficult to efficiently use the collected wavenumber spectrum; due to MP -BFSAR beams are often operated in mode to consistently obtain high-resolution images of the region of interest in many applications, which leads to complex time-frequency properties of the echoes such as Doppler aliasing.

The PFA algorithm (Parabola fitting algorithm: parabola fitting algorithm) deals with echoes in k-space (ie, wavenumber domain). PFA is suitable for processing SAR data due to its unique properties such as modeless Doppler aliasing, low platform trajectory constraints, and low computational complexity. At present, PFA has been widely used in monostatic SAR and bistatic SAR.

Since the PFA adopts the plane wave assumption, but the actual wavefront is curved, when the imaging scene is large, the wavefront curvature correction needs to be performed, which may be invalid for the large scene.

Documents "Polar format algorithm for bistatic SAR, IEEETrans.Aerosp.Electron.Syst., vol.40, no.4, pp.1147–1159, Oct.2004" and "Space-variant filtering for wavefront curvature correction in polar formatted bistatic SAR image, IEEE Trans.Aerosp.Electron.Syst., vol.48, no.2, pp.940–950, Apr.2012" proposed a k-space coordinate axis rotation method to support the wavenumber spectrum of the inclined parallelogram The domain is transformed into a horizontal parallelogram, however, the k-space axis rotation is only valid if the support region is rectangular or quasi-rectangular. Therefore, it is not suitable for MP-BFSAR where the support region is a parallelogram; the paper "Polar format algorithm wavefrontcurvature compensation under arbitrary radar flight path, IEEE Geosci.RemoteSens.Lett., vol.9, no.3, pp.526–530, May 2012", in order to compensate for the phase error caused by wavefront bending, a space-varying wavefront bending compensation filter is proposed, which exploits the internal homogeneity of the polar-to-rectangular resampling operation, however, the space-varying wavefront The forward curvature compensation filter is implicit and its calculation involves two-dimensional interpolation, which means it is inconvenient and inefficient in practice; the paper "Space-variant filtering for wavefront curvaturecorrection in polar formatted bistatic SAR image, IEEETrans.Aerosp. Electron.Syst.,vol.48,no.2,pp.940–950,Apr.2012", an analytical space-varying wavefront bending compensation filter is proposed, however, it is made by assuming that the platform moves in a straight line with uniform velocity Therefore, it is not suitable for MP-BFSAR because the motion trajectory of MP-BFSAR is highly nonlinear.

SUMMARY OF THE INVENTION

In order to solve the above technical problems, the present invention proposes a polar coordinate imaging method for bistatic forward-looking SAR on a mobile platform. The KAM method is used to transform the parallelogram wavenumber spectrum support domain into a horizontal quasi-rectangle, so that the collected wavenumber spectrum has The utilization rate has increased significantly.

The technical scheme adopted by the present invention is as follows: a polar coordinate format imaging method of a bistatic forward-looking SAR on a mobile platform, and a method of affine transformation (k-set affine mapping, KAM) in the support domain of the wave number spectrum is adopted to convert the parallelogram wave number of the MP-BFSAR The spectral support domain is transformed into a horizontal quasi-rectangle, and then the wavenumber spectrum is resampled from polar coordinates to rectangular coordinates by distance chirp-z transform (Chirp-Z Transform, CZT), and then distance compression and scene center motion compensation are performed. The Upsampling mapping accumulation (UMA) method is used in the azimuth direction to realize the resampling of the wavenumber spectrum from polar coordinates to rectangular coordinates; The imaging result; finally, the final imaging result is obtained through defocus correction.

Wavenumber Spectral Affine Transform (KAM), the specific process is as follows:

A1. Calculate echo delay gradient;

A2. Calculate echo Doppler gradient;

A3. Calculate the affine matrix and its inverse matrix according to the echo delay gradient and the echo Doppler gradient;

A4. Calculate the wavenumber in the x-coordinate and y-coordinate directions;

A5. Calculate the wavenumber distribution range after KAM according to the wavenumbers in the x-coordinate and y-coordinate directions and the affine matrix;

A6. Calculate the coordinate range of the image according to the wavenumber distribution range and the inverse matrix of the affine matrix.

The distance CZT realizes the resampling of the wavenumber spectrum from polar coordinates to rectangular coordinates; specifically:

B1, echo signal modeling;

B2. Perform distance CZT on the echo signal modeled in step B1;

Azimuth upsampling map accumulation, including the following processes:

C1, upsampling the azimuth direction of the signal obtained in step B2 by _Nu times;

C2. Perform azimuth signal accumulation on the signal after step C1.

Azimuth signal accumulation

C2. Calculate the mapping relationship from the n-axis to the k and _q -axis;

C3. Accumulate the signals obtained in step C1 according to the mapping relationship in step C2.

The final imaging result is obtained through defocus correction, which includes the following processes:

D1. Calculate the phase error caused by wavefront bending;

D2, according to the phase error of step D1, calculate and calculate the secondary phase error caused by wavefront bending;

D3. Construct a defocus correction filter according to the secondary phase error;

D4. Divide the coarse imaging result into several sub-images;

D5. Multiply the sub-images by the defocus correction filter constructed in step D3 to perform defocus correction;

D6. Combine all the sub-images after defocus correction obtained in step D5 to obtain a final imaging result.

Beneficial effects of the present invention: the method of the present invention adopts the KAM method to transform the support domain of the wavenumber spectrum of the parallelogram into a horizontal quasi-rectangle, so that the utilization rate of the collected wavenumber spectrum is increased by about 60%; It can be efficiently realized through Fourier transform (Fast fourier transform, FFT), IFFT and phase multiplication, so the PFA proposed by the present invention is more efficient; the analytical defocus correction filter proposed by the present invention is compared with the existing numerical filter. It is more convenient and efficient, and the defocus correction of the present invention enables the PFA algorithm to be suitable for large imaging scenes.

Description of drawings

Fig. 1 is the imaging processing flow chart of the method of the present invention;

Fig. 2 is the geometric model of the bistatic SAR echo recording of the mobile platform according to the present invention;

Figure 3 shows the wavenumber spectrum support domain before and after KAM;

Among them, Figure 3(a) is the wavenumber spectrum support domain before KAM; Figure 3(b) is the wavenumber spectrum support domain after KAM;

Fig. 4 is the imaging result of point target;

Fig. 5 is the enlarged contour map of the imaging result of the point target in Fig. 4;

Figure 6 shows the imaging results of a surface target.

Detailed ways

In order to facilitate those skilled in the art to understand the technical content of the present invention, the content of the present invention will be further explained below with reference to the accompanying drawings.

As shown in FIG. 1 is the scheme flow chart of the present invention, including: adopting the method of wavenumber spectrum support domain affine transformation (k-setaffine mapping, KAM), the parallelogram wavenumber spectrum support domain of MP-BFSAR is transformed into a horizontal quasi-rectangle , and then realize the resampling of the wavenumber spectrum from polar coordinates to rectangular coordinates by distance chirp-z transform (Chirp-Z Transform, CZT), and then perform distance compression and scene center motion compensation, and then use upsampling map accumulation in the azimuth direction. Upsampling mapping accumulation (UMA) method realizes resampling of wavenumber spectrum from azimuth polar coordinates to rectangular coordinates; then performs two-dimensional inverse fast fourier transform (IFFT) to obtain coarse imaging results; finally defocusing Correction to get the final imaging result.

Specifically, it includes the following 5 steps:

1. Wavenumber Spectral Affine Transform (KAM):

A1. The echo delay gradient is

where c is the speed of light, Respectively represent the unit vector of the line-of-sight direction from the transmitting station and the receiving station to the center of the scene at azimuth zero time, r _T0 = (x _T , y _T , z _T ), r _R0 = (x _R , y _R , z _R ) represent the initial moment of transmission, respectively The position vectors of the station and the receiving station, R _T0 =|r _T0 |, R _R0 =|r _R0 | respectively represent the distance from the zero time azimuth of the transmitting station and the receiving station to the center of the scene.

The projection of the delay gradient on the ground is

G _gτ =PG _τ

in, Represents the ground projection matrix.

A2. Calculate the Doppler gradient

The echo Doppler gradient is

Among them, λ represents the carrier wavelength, v _T = (v _Tx , v _Ty , v _Tz ), v _R = (v _Rx , v _Ry , v _Rz ) represent the speed of the transmitting station and the receiving station at zero azimuth time, respectively, and I represent identity matrix.

The projection of the Doppler gradient on the ground is

G _gd =PG _d

A3. Calculate affine matrix and its inverse

in, represents the angle between the vector direction and the positive direction of the horizontal axis, and φ _c = |θ _gd - θ _gτ | represents the coupling angle of the wavenumber spectrum.

Invert the affine matrix to get

A4. Calculate the wavenumber in the x-coordinate and y-coordinate directions

ΔR(η; x, y)=R _bi (η; x, y)-R _bi (0; x, y) represents the differential distance, and a first-order Taylor approximation is made for the differential distance:

ΔR(n; x, y) = R _SAG (n; x, y) - R _SAG (0; x, y) ≈ C ₁₀ x + C ₀₁ y,

in

The projections of the wavenumber k on the x-axis and the y-axis are obtained as k _x =kC ₁₀ , k _y =kC ₀₁ respectively, where

k=2π(f _r +f _c )/c

That is, the wavenumbers in the x-coordinate and y-coordinate directions are respectively

k _x =kC ₁₀

k _y =kC ₀₁

R _bi (n; x, y) = R _T (n; x, y) + R _R (n; x, y) is the bistatic distance sum, n is the slow time,

R _T (n; x, y), R _R (n; x, y) represent the distance history of the transmitting station and the receiving station, respectively, a _T = (a _Tx , a _Ty , a _Tz ), a _R = (a _Rx , a _Ry , a _Rz ) represent the acceleration of the transmitting station and the receiving station respectively, and r _P =(x, y, 0) represents the position coordinates of any point target in the scene.

A5. The wavenumber distribution range after calculating KAM

Among them, k _xmin and k _xmax represent the minimum and maximum values of k _x , respectively, and k _ymin and k _ymax represent the minimum and maximum values of _ky , respectively;

The distribution range of the range-directed wavenumber variable k _p after KAM is [k _pmin ,k _pmax ], where

k _pmin =min{k _p1 ,k _p2 ,k _p3 ,k _p4 }

k _pmax =max{k _p1 ,k _p2 ,k _p3 ,k _p4 }

The distribution range of the azimuthal wavenumber variable k _q behind the KAM is [k _qmin ,k _qmax ], where

k _qmin =min{k _q1 ,k _q2 ,k _q3 ,k _q4 }

k _qmax =max{k _q1 ,k _q2 ,k _q3 ,k _q4 }

A6. Calculate the coordinate range of the output image, as follows:

Wherein, x _min and x _max represent the minimum and maximum values of the abscissa of the original scene, respectively, and y _min and y _max respectively represent the minimum and maximum values of the ordinate of the original scene;

The range of distance coordinates of the output image is [p _min ,p _max ], where

p _min =min{p ₁ ,p ₂ ,p ₃ ,p ₄ }

p _max =max{p ₁ ,p ₂ ,p ₃ ,p ₄ }

The azimuth coordinate range of the output image is [q _min , q _max ], where

q _min =min{q ₁ ,q ₂ ,q ₃ ,q ₄ }

q _max =max{q ₁ ,q ₂ ,q ₃ ,q ₄ }

2. Distance to CZT:

B1, echo signal modeling

The range-frequency-domain azimuth baseband echo can be modeled as follows

Among them, K _r represents the frequency of the emission pulse modulation, τ represents the fast time, c represents the speed of light, and λ represents the carrier wavelength.

B2, distance to CZT

For the echo signal, find the spectrum of the echo at the following position on the complex plane unit circle from the CZT through the distance

z _i =exp{j(φ ₀ +iΔφ)},i=0,1,2,...,N _p

in,

where F _s is the sampling rate of fast time t, Represents the size of the central wave number, k _pmin and k _pmax represent the maximum and minimum values of k _p , respectively, and N _p is the number of sampling points in the distance direction of the output image, which must satisfy the following conditions:

The signal after distance to CZT can be expressed as

Among them, k _p represents the distance wavenumber variable after KAM.

3. Match filter (MF, Match filter) and scene center motion compensation (MoCo to SC), that is, multiply the result obtained by the distance CZT by the following phase function:

The signal after MF and MoCo to SC can be approximated as

Among them, p and q represent the range and azimuth coordinate variables of the output image, respectively.

4. Azimuth upsampling mapping accumulation:

C1, upsampling _Nu times the azimuth direction of the signal obtained in step 3, and record the obtained signal as S ₃ (k _p ,n(m); p,q), where m represents the sampling sequence number of n after the upsampling

Wherein, Na represents the number of sampling points of _n before _upsampling , and Nu generally takes the value of 4 or 8;

C2. Calculate the mapping relationship from the η axis to the k _q axis:

Among them, round{·} represents the rounding operation, Represents the sampling interval of the k and _q axes;

C3, azimuth signal accumulation

Accumulate the signal obtained by C1 in the azimuth direction according to the mapping relationship obtained by C2:

Among them, n and N _q respectively represent the sampling point number and sampling point number of the signal on the k _q axis.

So far, the signal in the (k _p ,k _q ) domain is obtained, which can be approximately expressed as

S ₄ (k _p ,k _q ; p,q)≈exp{-j(k _p p+k _q q)}

C4. Perform a two-dimensional inverse Fourier transform (IFFT) on the signal obtained by C3 to obtain a rough imaging result.

5. For large scenes, the defocus effect of edge targets caused by wavefront bending cannot be ignored, and defocus correction needs to be performed on the coarse imaging results:

D1. Calculate the phase error caused by wavefront bending, as follows:

The phase error due to wavefront bending is

Φ _WC (k, η; x, y)=-kΔR(η; x, y)+k _x x+k _y y

D2. Calculate the quadratic phase error caused by wavefront bending

The quadratic phase error caused by wavefront bending is calculated by the implicit function theorem and the composite function derivation rule as

D3. Build a defocus correction filter

D4, sub-image division

Divide the coarse imaging result into N _sub sub-images along the q-axis, and then transform each sub-image into the (p,k _q ) domain through the Fourier transform (Fast fourier transform, FFT) about q, where N _sub is satisfying a positive integer of the formula

where q _m is solved by the following equation

D5, sub-image defocus correction

Multiply each sub-image by the following defocus correction filter

H _RF (k _q ; p,q _c )=exp{-jk ² D(p,q _c )}

Among them, q _c is the q-axis coordinate of the center of the sub-image, and then IFFT about k _q is performed to obtain the sub-image after defocus correction in the (p, q) domain.

D6, sub-image stitching

Combine all the defocus corrected sub-images from D5 to get the final imaging result.

The technical effects of the present invention will be described below in conjunction with specific data.

Table 1 shows the system parameters used in this embodiment, and Fig. 2 is a schematic diagram of the geometric configuration of the mobile bistatic SAR in this embodiment. In this embodiment, it is assumed that there are 49 uniformly distributed point targets on the ground, and the x-axis and y-axis are The distances between the above point targets are all 800 meters.

Table 1 System parameters

Various calculations are performed based on the system parameter values in Table 1, specifically:

The echo delay gradient calculation result is:

The projection of the delay gradient on the ground is: G _gτ =PG _τ =1.0×10 ^-8 [-0.2711,0.2634] ^T

The echo Doppler gradient calculation result is:

The projection of the Doppler gradient on the ground is: G _gd =PG _d =[-0.0001,0.0458] ^T

The calculation result of the coupling angle is:

φ _c =|θ _gτ -θ _gd |=45.7°

The affine matrix is:

Inverting the affine matrix, we get:

Perform the following affine transformation on the wavenumber spectrum

The size of the central wave number

The maximum and minimum calculation results of k _p are:

k _pmin =28.38rad/m, k _pmax =30.99rad/m,

The number of sampling points in the distance direction of the output image is: N _p =3328.

Nu = 8, the sampling sequence number of n after up _- sampling is:

The pulse repetition interval is: PRI=1/PRF=4ms,

The number of sampling points of the signal on the k _q axis is: N _q =2224,

In addition, from k _qmax =-0.9044rad/m, k _qmax =0.8413rad/m, the sampling interval of k _q axis is obtained as

The sampling number of the signal on the k and _q axis is:

Divide the coarse focus map into N _sub =8 sub-images along the q-axis, and transform each sub-image into the (p,k _q ) domain by FFT on q;

The defocus correction filters are:

Wherein, D(p,q _c )=7.31×10 ^-6 p ² +1.02×10 ^-5 pq _c +1.03×10 ^-5 q ² ,

p∈[-4.0020×10 ³ ,3.9996×10 ³ ]m; q _c is the q-axis coordinate corresponding to the center of each sub-image, and its value set is [-3500,-2500,-1500,-500,500,1500,2500 ,3500]m;

According to the above calculation results, the simulation images shown in Figure 3, Figure 4, Figure 5, and Figure 6 are obtained; the simulation results of the wavenumber spectrum support domain before and after KAM are shown in Figure 3. After KAM, the use of the wavenumber spectrum support domain of the point target The efficiency is increased from 33% to 98%. It can be seen from this that the method proposed in the present invention can effectively improve the utilization efficiency of the wavenumber spectrum. The imaging result of the point target of the present invention is shown in Figure 4; the corresponding enlarged contour map of the imaging result of the point target is shown in Figure 5. It can be seen that the contour map of the imaging result of the point target is evenly and regularly distributed, and its peak sidelobe ratio ( PSLR) and integral side lobe ratio (ISLR) are very close to the theoretical values of -13.26dB and -9.96dB, indicating that the algorithm proposed in the present invention has a good focusing imaging effect on the point target echo; the imaging result of the surface target is shown in Figure 6, It can be seen in FIG. 6 that the surface target is clearly distinguishable, which indicates that the algorithm proposed in the present invention has a good focusing imaging effect of the echo of the surface target. 4, 5 and 6, it can be seen that the method proposed in the present invention can effectively focus and image the echoes of the bistatic forward-looking SAR target on the mobile platform.

Those of ordinary skill in the art will appreciate that the embodiments described herein are intended to assist readers in understanding the principles of the present invention, and it should be understood that the scope of protection of the present invention is not limited to such specific statements and embodiments. Various modifications and variations of the present invention are possible for those skilled in the art. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included within the scope of the claims of the present invention.

Claims (7)

1. The polar coordinate format imaging method of the bistatic forward-looking beamforming synthetic aperture radar on the mobile platform is characterized in that, Firstly, the parallelogram wavenumber spectrum support domain of MP-BFSAR is transformed into a horizontal quasi-rectangle through affine transformation of wavenumber spectrum; secondly, the resampling of wavenumber spectrum from polar coordinates to rectangular coordinates is realized by distance chirp-z transformation; , performing distance compression and scene center motion compensation; resampling from the secondary, wavenumber spectrum azimuth to polar coordinates to rectangular coordinates; then, performing two-dimensional inverse Fourier transform to obtain coarse imaging results; finally, defocus correction to obtain the final imaging results. 2. the polar coordinate format imaging method of mobile platform bistatic forward-looking SAR according to claim 1, is characterized in that, wavenumber spectrum affine transformation, concrete process is as follows: A1. Calculate the echo delay gradient; A2. Calculate echo Doppler gradient; A3. Calculate the affine matrix and its inverse matrix according to the echo delay gradient and the echo Doppler gradient; A4. Calculate the wave number in the x-coordinate and y-coordinate directions; A5. Calculate the wavenumber distribution range after KAM according to the wavenumbers in the x-coordinate and y-coordinate directions and the affine matrix; A6. Calculate the coordinate range of the image according to the wavenumber distribution range and the inverse matrix of the affine matrix. 3. the polar coordinate format imaging method of mobile platform bistatic forward looking SAR according to claim 2, is characterized in that, the distance to CZT realizes the resampling of wavenumber spectrum distance to polar coordinate to Cartesian coordinate; Be specially: B1, echo signal modeling; B2. Perform range CZT on the echo signal modeled in step B1. 4. the polar coordinate format imaging method of the bistatic forward-looking SAR of the motorized platform according to claim 3, is characterized in that, the resampling of wavenumber spectrum azimuth direction polar coordinate to rectangular coordinate, is specially: adopt upsampling mapping in azimuth direction Accumulation; includes the following steps: C1, upsampling the azimuth direction of the signal obtained in step B2 by _Nu times; C2. Perform azimuth signal accumulation on the signal after step C1. 5. the polar coordinate format imaging method of the bistatic forward-looking SAR of the motorized platform according to claim 4, is characterized in that, the signal after step C1 adopts is: S ₃ (k _p ,n (m); p, q) ; Among them, k _p represents the range wavenumber variable after the affine transformation of the wavenumber spectrum support domain, η represents the slow time, p represents the range coordinate of the image, q represents the azimuth coordinate of the image, and m represents the sampling sequence number of η after upsampling. 6. the polar coordinate format imaging method of mobile platform bistatic forward looking SAR according to claim 5, is characterized in that, azimuth direction signal accumulation is specially: First, calculate the mapping relationship from the n-axis to the k _q -axis; Then, the signals obtained in step C1 are accumulated according to the mapping relationship. 7. The polar coordinate format imaging method of the bistatic forward-looking SAR motorized platform according to claim 6, is characterized in that, obtains final imaging result by defocus correction, specifically comprises following process: D1. Calculate the phase error caused by wavefront bending; D2, according to the phase error of step D1, calculate the secondary phase error caused by wavefront bending; D3. Construct a defocus correction filter according to the secondary phase error; D4. Divide the coarse imaging result into several sub-images; D5. Multiply the sub-images by the defocus correction filter constructed in step D3 to perform defocus correction; D6. Combine all the sub-images after defocus correction obtained in step D5 to obtain a final imaging result.