A satellite-aircraft dual-base synthetic aperture radar combined mode imaging method

Abstract

The present invention discloses a satellite-aircraft bistatic synthetic aperture radar combined mode imaging method, belonging to the field of signal processing technology. The phase error introduced by azimuth spectrum aliasing in the range processing process is avoided by improving the range processing flow; the azimuth processing flow is improved, while avoiding aliasing in the time-frequency domain of azimuth processing, more degrees of freedom are introduced, and the high-order space variation in azimuth is balanced; an accurate calculation method for Doppler space variation parameters of combined mode satellite-aircraft bistatic SAR is proposed, which can obtain high-precision de-skew factors and NLCS equalization factors, and realize efficient and high-precision processing of combined mode satellite-aircraft bistatic SAR. The method of the present invention solves the problem of imaging processing difficulties caused by the complex range migration characteristics and Doppler characteristics of bistatic SAR echo signals, and solves the problem of decreased keystone interpolation accuracy caused by azimuth spectrum aliasing. Simulation results show that the method can effectively realize accurate imaging processing in the satellite-aircraft bistatic SAR combined imaging mode.

Description

Technical Field

[0001] The invention belongs to the technical field of signal processing, and in particular relates to a satellite-machine dual-base synthetic aperture radar combined mode imaging method. Background Art

[0002] Spaceborne synthetic aperture radar (SAR) has all-day and all-weather imaging capabilities and can be applied to terrain mapping, resource exploration, military reconnaissance and other fields. Compared with spaceborne monostatic SAR, spaceborne bistatic SAR has a more flexible configuration, can obtain the non-backscattering coefficient of ground objects, has stronger concealment, better anti-interference and survivability, and has broad application prospects. With the rapid development of SAR technology, the demand for SAR image size and resolution has also continued to deepen. Compressed sensing technology, multi-channel technology and multiple working modes are used to improve the size and resolution of SAR images. Spaceborne SAR, such as TerraSAR-X and GF-3, have multiple imaging modes to obtain ground observation images of different resolutions and sizes, such as strip, scanning, TOPS, sliding beam, etc. For spaceborne bistatic SAR, there is a large speed difference between platforms, and it is difficult to coordinate the transceiver beam. The transceiver platform adopts different imaging modes to achieve a trade-off between spatial resolution and imaging width, which is conducive to the coordination of transceiver beams, ensure the quality of echoes, and achieve the improvement of resolution and width. In addition, several research institutions have conducted imaging tests on the combined working mode of satellite-aircraft bistatic SAR. The German Space Agency (DLR) conducted a satellite-aircraft bistatic SAR imaging test, with the satellite working in sliding beamforming mode and the aircraft working in reverse sliding beamforming mode or strip mode, which compensated for the speed difference between beam footprints, improved the beam coordination time, and increased the imaging width. The University of Electronic Science and Technology of China used a similar imaging mode in the oblique satellite-aircraft bistatic SAR imaging test, which increased the data acquisition time within a shorter startup time window, and improved the imaging resolution and imaging width.

[0003] Due to the large speed difference and spatial scale difference between satellite and airborne platform, the satellite-aircraft bistatic SAR presents heterogeneous and shift-variant characteristics in geometric structure, that is, the geometric parameters of the satellite-aircraft bistatic SAR, such as squint angle, pitch angle, and bistatic angle, will change with the azimuth slow time. In SAR echo data, this shift-variant characteristic is manifested as different range cell migration (RCM) and Doppler parameters for different target points in azimuth, and will lead to two problems, namely how to efficiently solve the coupling of echo signals and how to achieve consistent processing of azimuth signals. In response to these two problems, some scholars proposed to achieve range migration balance through CSA and balance the azimuth third-order phase space variability through azimuth nonlinear scaling (NLCS, Nonlinear Chirp Scaling). However, CSA cannot solve more serious range migration space variability, such as squint. Some scholars proposed to solve the azimuth space variability through the subaperture method, but in the case of more serious space variability, the efficiency will be sharply reduced. Some scholars have combined linear range migration correction with the NLCS method to achieve azimuth space-variation correction for bistatic SAR. Some scholars have solved the two-dimensional coupling problem by using the keystone transform and high-order RCM compensation method, and achieved azimuth Doppler phase balance by introducing multiple NLCS factors. However, the keystone transform needs to be implemented by interpolating the azimuth signal. When the PRF is smaller than the total width of the azimuth spectrum, it will cause the azimuth spectrum to fold, resulting in ambiguity. Then for the combined mode satellite-machine bistatic SAR, the above method is no longer applicable.

[0004] For the combined mode satellite-aircraft bistatic SAR, due to the beam rotation of the transceiver antenna, the observation angle of its integrated beam relative to the transceiver platform changes continuously with the azimuth slow time. Therefore, compared with the strip mode, the combined imaging mode of the satellite-aircraft bistatic SAR has a more complex azimuth space-variant characteristic. On the one hand, the azimuth space-variant characteristic will further broaden the azimuth spectrum, resulting in azimuth spectrum aliasing; on the other hand, the high-order azimuth space-variant characteristic will cause the equalization effect of the NLCS factor to deteriorate, that is, more NLCS factors need to be introduced to solve the high-order space-variant. Some scholars have achieved an equivalent improvement in the azimuth spectrum sampling rate through multi-channel technology, solved the problem of azimuth spectrum folding, and achieved azimuth space-variant equalization through two NLCS. However, multi-channel technology puts forward higher hardware requirements for the receiving system, and the phase error between channels will also affect the imaging quality. In addition, the two NLCS in this method do not make full use of all the scaling factors, but compensate for the azimuth phase after the first NLCS, so it cannot effectively solve the high-order azimuth high-order phase, and cannot be applied to combined mode bistatic SAR imaging. Some scholars have used azimuth preprocessing methods to equivalently improve the azimuth sampling rate, and used two-dimensional Stolt interpolation to remove RCM spatial variation and azimuth phase spatial variation. However, two-dimensional Stolt interpolation can still only solve low-order spatial variation, thus limiting the imaging depth of field. Some scholars have used azimuth de-skewing methods to solve the Doppler centroid spatial variation problem in the beamforming mode, but the accuracy of the distance model affects the imaging accuracy, and the de-skewing factor cannot be accurately calculated, which will not effectively solve the azimuth spectrum aliasing problem under certain circumstances. Summary of the invention

[0005] In order to solve the above technical problems, the technical solution adopted by the present invention is: a satellite-aircraft bistatic SAR combined mode imaging method, the specific steps are as follows:

[0006] Step S1: input echo signals, geometric parameters, and beam steering parameters to construct a satellite-machine bistatic synthetic aperture radar SAR combination mode geometric model and signal model;

[0007] Step S2: solving the Doppler space-variable coefficient;

[0008] Step S3: Perform a range Fourier transform to obtain the signal spectrum of the original echo, and then compensate the spectrum, specifically including pulse compression, spectrum tilt correction, and spectrum regularization;

[0009] Step S4: Performing Keystone transformation on the spectrum corrected in step S3 to obtain a transformed signal;

[0010] Step S5: After performing high-order range motion correction, secondary range compression, and additional phase compensation on the transformed signal, perform inverse Fourier transform in the range direction to obtain the signal ;

[0011] Step S6: Azimuth Fourier transformation is performed to obtain the azimuth phase; then the non-space-variant term of the azimuth space-variant coefficient is compensated; and then the azimuth degree of freedom is introduced by multiplying the variable scaling function;

[0012] Step S7: Perform an inverse Fourier transform in azimuth to obtain the azimuth time domain echo phase;

[0013] Step S8: Compensate the azimuth time domain echo phase to remove the non-space variable term in the time domain; then multiply it by the NLCS equalization function;

[0014] Step S9: perform Fourier transform in azimuth to obtain frequency domain phase in azimuth;

[0015] Step S10: Obtain and multiply the azimuth pulse pressure function according to the azimuth frequency domain phase ;

[0016] Step S11: Perform inverse Fourier transform on the azimuth signal to obtain an imaging result.

[0017] The beneficial effects of the present invention are as follows:

[0018] The present invention analyzes in detail the spatial variation law of RCM and Doppler parameters in the combined imaging mode, proposes a calculation method for each order of spatial variation rate, and proposes an improved NLCS to solve the spatial variation problems in the range and azimuth directions. In view of the problem that the keystone method is not applicable due to the spatial variation of the Doppler centroid, the present invention improves the range processing flow, and implements the keystone transformation after equalizing the azimuth centroid. On the one hand, it avoids the influence of azimuth spectrum folding on the keystone, and on the other hand, it equivalently improves the azimuth time domain sampling accuracy, thereby avoiding the phase error introduced by the keystone interpolation processing. In addition, based on the spatial variation law of the combined imaging mode, a calculation method for the combined imaging spatial variation factor is proposed, and the NLCS equalization factor is improved, while avoiding the aliasing problem in the azimuth time-frequency domain conversion process and realizing the spatial variation equalization of the azimuth Doppler parameters, thereby realizing the azimuth focusing of the combined mode echo. BRIEF DESCRIPTION OF THE DRAWINGS

[0019] Figure 1 The present invention is a flow chart of the satellite-aircraft bistatic SAR combined mode imaging method.

[0020] Figure 2 It is a point target arrangement diagram in the point target simulation of the present invention.

[0021] Figure 3 It is a geometric configuration diagram of the satellite-aircraft bistatic SAR combination mode in the point target simulation of the present invention.

[0022] Figure 4It is the imaging resolution of the satellite-aircraft bistatic SAR combination mode in the point target simulation of the present invention.

[0023] Figure 5 This is the imaging result corresponding to the de-skewing process in the point target simulation of the present invention.

[0024] Figure 6 It is the distance processing result in the point target simulation of the present invention.

[0025] Figure 7 It is the point target imaging result in the point target simulation of the present invention.

[0026] Figure 8 It is a point target contour map in the point target simulation of the present invention.

[0027] Fig. 9 It is an azimuth profile diagram of a point target in the point target simulation of the present invention. DETAILED DESCRIPTION

[0028] The present invention is mainly verified by simulation experiment, and all steps and conclusions are verified to be correct on Matlab 2012. The method of the present invention is further described below in conjunction with the accompanying drawings and embodiments.

[0029] In order to facilitate the description of the content of the method of the present invention, the following terms are first explained:

[0030] Term 1: Bistatic SAR;

[0031] Bistatic SAR refers to a SAR system in which the transmitting station and the receiving station are placed on different platforms, at least one of which is a moving platform. Conceptually, it is a bistatic radar.

[0032] like Figure 1 As shown, a flow chart of a satellite-aircraft bistatic SAR combined mode imaging method provided by the present invention, the specific steps are as follows:

[0033] Step S1: input echo signals, geometric parameters, and beam steering parameters to construct a satellite-aircraft bistatic SAR combination mode geometric model and signal model;

[0034] Input echo signal, geometric parameters, and beam control parameters for parameter calculation in subsequent steps. This simulation experiment adopts point target simulation verification. In the simulation, the illumination source is a low-orbit satellite, operating in the C band, with an altitude of 700km and a speed of about 7km / s. It flies along the azimuth direction and illuminates the imaging target area in a positive side view. The illumination source works in sliding beam mode. The airborne receiving platform flies at an altitude of 10km, with a slant angle of about 4°, and works in progressive scanning imaging mode (TOPS, Terrain Observation with Progressive Scans). In addition, the simulation data acquisition time is 8.25s. Since the azimuth beam of the airborne platform is narrow, the synthetic aperture time is determined by its azimuth beam width, so the synthetic aperture time is 0.6s. It should be noted that the data window of the satellite-based bistatic SAR is short, so this simulation only images the data within the illumination time of a single sub-band in the TOPS mode. In addition, the point target scenes used for simulation experiments are as follows Figure 2 The point target scene is 1.5km×6km in size, with 5×5 points evenly distributed, and the lower left corner point is taken as P1, the center point is O, and the upper right corner point is P2.

[0035] The simulation parameters are shown in Table 1.

[0036] Table 1 System parameters of satellite-aircraft bistatic SAR combination mode

[0037] parameter Numeric parameter Numeric Carrier frequency 5.4 GHz Signal bandwidth 80 MHz Pulse Width 1μs Sampling bandwidth 120 MHz Pulse Repetition Frequency (PRF) 4000Hz Data collection time 8.25 s Synthetic Aperture Time 0.6 s Satellite altitude 700 km Satellite inclination 98° Satellite view 6° Satellite working mode Sliding bunching Beam Steering Factor 0.12 Aircraft operating mode TOPS Beam Steering Factor 3.3 Aircraft location [-15, -1, 10]km Aircraft speed [0, 240, 0]m / s

[0038] Geometric model construction: The bistatic SAR transceiver platform works in different working modes to meet different imaging mission requirements. It can also integrate the advantages of the transceiver platform combination to build a satellite-aircraft bistatic SAR earth observation system and achieve performance complementarity of the transceiver platform. Figure 3 As shown. In the Cartesian coordinate system, the origin represents the center of the imaging scene, and the y-axis is along the direction of the receiving platform velocity. The satellite platform operates in the sliding beamforming mode, and the airborne receiving platform operates in the TOPS mode. In the TOPS mode, as the airborne receiving platform moves, three sub-strips are scanned in sequence, first scanning sub-strip 1, and then scanning the next sub-strip (sub-strip 2) in the same way, and finally scanning the third sub-strip (sub-strip 3). The virtual rotation center of the transmitting station represents the virtual rotation center of the satellite platform beam, and the virtual rotation center of the receiving station represents the virtual rotation center of the airborne receiving platform beam.

[0039] Signal model construction: The beam movement speed of the transceiver platform in the satellite-aircraft bistatic SAR combined imaging mode is different. For the echo model, it is mainly reflected in the azimuth envelope and phase. For the satellite-aircraft bistatic SAR, the point target echo model of each imaging mode can be expressed as:

[0040] s r ( τ , t ; x , y ) = σ ( x , y ) w a [ t − t c ( x , y ) T a ] exp { j π B r T r [ τ − R b ( t ; x , y ) / c ] 2 } exp { − j 2 π R b ( t ; x , y ) λ }

[0041] in, represents the scattering coefficient of the target point, w a Indicates the azimuth signal envelope; represents bandwidth, Indicates the pulse width, Indicates distance, represents the wavelength, represents the speed of light, represents the double base distance history, j represents the sign of the imaginary number; Indicates the direction and time. represents the synthetic aperture time, represents the central irradiation time. This time variable is related to the bistatic SAR integrated beam, that is, the combined beam of the transceiver platform. Without loss of generality, it is assumed here that the azimuth width of the combined beam of the satellite bistatic SAR is ,therefore and It can be expressed as follows:

[0042]

[0043]

[0044] in, represents the ground projection of the receiving station's oblique angle, represents the ground moving speed of the beam, which is determined by the beam control factor Determine. And can be expressed as

[0045]

[0046] in, Indicates the speed of the transmitting station or receiving station. The specific formula is as follows:

[0047]

[0048] in, , , and , , Represent the location coordinates of the transmitting station and the receiving station respectively. , , is the beam center position. Indicates the direction and time. is the spatial variation coefficient of each order of Doppler, i=1,2…. And it can be expressed as:

[0049]

[0050] in, , and represents the position of the satellite at the reference zero time, , and represents the satellite speed, , and Represents the satellite acceleration. The receiving station coordinates can be expressed in the same way.

[0051] Step S2: solving the Doppler space-variable coefficient;

[0052] Doppler spatial variation coefficient Performing Taylor expansion, we can get:

[0053]

[0054] in, , , and Respectively represent the first, second, third, and fourth order Doppler parameter space-varying coefficients, i=1,2…, is the azimuth time, for hour At this time,

[0055]

[0056] in, is the symbol of partial derivative, represents the central irradiation time, Indicates the target point Axis position, Indicates the target point Axis position, . and satisfy:

[0057]

[0058] in, represents the central irradiation time, Indicates the time when different target points are illuminated at their azimuth center The corresponding double base distances and can be expressed as:

[0059]

[0060] also, The spatial variation coefficients of each order can be expressed as:

[0061]

[0062] Similarly, the other order space-varying coefficients , and It can be found in a similar way.

[0063] Step S3: Range Fourier transform (FFT) to obtain the spectrum of the original echo , then pulse compress the spectrum and compensate the phase;

[0064] The spectrum of the original echo is expressed as:

[0065] ss ( f τ , t ) = σ ( x , y ) W r [ f τ B r ] w a [ t − t c ( x , y ) T a ] exp { − j 2 π f c + f τ c R b ( t ; x , y ) } exp ( j π f τ 2 k r )

[0066] in, W r [ · ] represents the range signal envelope, represents the distance frequency, is the carrier center frequency, Indicates the range frequency modulation slope;

[0067] Then the spectrum is multiplied by the compensation phase function to achieve pulse compression and compensate the phase. The compensation phase function is expressed as follows:

[0068]

[0069] Step S4: construct a linear range movement (LRCM) correction function, multiply it with the spectrum, and correct the spectrum tilt;

[0070] For satellite-based bistatic SAR, the curved trajectory of the satellite platform, the speed difference between the transmitting and receiving platforms, and the change in the azimuth beam angle will lead to severe coupling characteristics between range and azimuth. Therefore, a linear range movement (LRCM) correction function can be constructed. , to compensate for the linear component of the range walk (RCM) of the reference point in the center of the imaging scene and to correct for the spectral tilt:

[0071]

[0072] Step S5: construct a spectrum regularization function and multiply it with the spectrum to achieve spectrum regularization;

[0073] Due to the spatial variation of the Doppler centroid, the azimuth spectrum is aliased. If the signal is converted to the azimuth frequency domain, on the one hand, it will affect the accuracy of the subsequent Keystone transformation, and on the other hand, it will cause errors in the subsequent NLCS method. At this time, the centroid needs to be equalized to achieve spectrum regularization. The spectrum regularization function is as follows:

[0074] H de ( t ) = exp [ − j π K Rot t 2 ]

[0075] in, represents the deskewing coefficient, which can be expressed as

[0076]

[0077] Finally, the corrected signal is obtained, which is expressed as follows:

[0078] S d ( f τ , t ) = σ ( x , y ) W r [ f τ B r ] w a [ t − t c ( x , y ) T a ] exp { − j 2 π f c + f τ c [ R b ( t ; x , y ) − k 1 ref t ] } exp { j π a 1 λ t 2 }

[0079] Step S6: Keystone transformation to obtain the transformed signal;

[0080] Perform keystone transformation, the specific interpolation method is as follows:

[0081]

[0082] It is important to note that after the keystone transformation, the sampling rate of the azimuth direction changes, and the new sampling rate is recorded as , and by The transformed signal can be expressed as

[0083] S d ( f τ , t ) = W r [ f τ B r ] w a [ f c t m f c + f τ − f c t c , m ( x , y ) f c + f τ T a , m ] exp { − j 2 π f c + f τ c k 0 − j 2 π f c c ( k 1 − k 1 ref ) t m } exp { − j 2 π f c + f τ c [ k 2 ( f c t m f c + f τ ) 2 + k 3 ( f c t m f c + f τ ) 3 + k 4 ( f c t m f c + f τ ) 4 ] } exp { j π a 1 λ ( f c t m f c + f τ ) 2 }

[0084] in, Indicates the azimuth center illumination time corresponding to different target points in the azimuth time domain after keystone transformation; It represents the synthetic aperture time in the azimuth time domain after the keystone transformation. The first phase indicates the reference center point compensation phase, the second phase includes the high-order RCM, the secondary range frequency modulation phase and the azimuth phase, and the third phase includes the additional RCM introduced by the azimuth de-skewing and the azimuth spectrum equalization phase.

[0085] Step S7: Perform high-order range migration correction on the transformed signal, multiplying it by ;

[0086] The second phase of the transformed signal contains high-order RCM. At this time, the RCM in the range direction still has high-order null-variant terms, so it is necessary to compensate for the high-order RCM and construct a high-order range migration correction function. as follows:

[0087] H HRCMC ( f τ ) = exp [ − 2 π c f τ ( k 2 ref t m 2 + 2 k 3 ref t m 3 + 3 k 4 ref t m 4 ) ]

[0088] Step S8: Perform secondary distance compression, multiply by ;

[0089] At this time, the keystone transformation introduces an additional range modulation rate, see the second phase. Therefore, it is necessary to compensate for the range modulation rate and construct a quadratic range compression function as follows:

[0090] H SRC ( f τ ) = exp [ j 2 π cf c f τ 2 ( k 2 ref t m 2 + 2 k 3 ref t m 3 + 3 k 4 ref t m 4 ) ]

[0091] Step S9: Additional compensation phase, multiplied by ;

[0092] In addition, since the azimuth centroid balance factor is introduced before the keystone, there will be additional azimuth offset terms and high-order azimuth phases after the keystone. Therefore, they need to be compensated. The additional compensation phase as follows:

[0093] H com ( f τ ) = exp [ − j ( 2 π K Rot f c t m 2 f τ − 3 π K Rot f c 2 t m 2 f τ 2 + 4 π K Rot f c 3 t m 2 f τ 3 − 20 π K Rot f c 4 t m 2 f τ 4 ) ]

[0094] Step S10: Distance-directed inverse Fourier transform (IFFT) to obtain the signal ;

[0095] The purpose of azimuth processing is to introduce enough equalization factors to ensure high degrees of freedom of each order of space-variant phase, so as to achieve the purpose of equalization of each order of space-variant coefficients in azimuth. However, it should be noted that in the process of NLCS equalization, azimuth phase cannot be aliased in both time domain and frequency domain, that is, non-aliasing must be ensured before and after each Fourier transform and inverse transform. After step S7, RCM has been corrected. At this time, when converted to the distance time domain, it can be approximately seen that RCM has been eliminated. Let , the signal can be obtained as

[0096] S d ( τ , t m ) = sin c { B r [ τ − R bi ( t m ; x , y ) c ] } w a ( t m ) exp { − j 2 π [ k 0 + ( k 1 − k 1 ref ) t m + k 2 t m 2 + k 3 t m 3 + k 4 t m 4 ] λ }

[0097] Step S11: Azimuth Fourier transform to obtain azimuth phase ;

[0098] Fourier transform to frequency domain can get azimuth phase for

[0099]

[0100] in, represents the azimuth frequency variable, is the azimuth time, , , and represents the non-space-varying azimuth frequency modulation coefficient and can be expressed as follows:

[0101]

[0102] in, Represents the spatial variation coefficients of various orders in azimuth.

[0103] Step S12: Obtain and multiply the compensation filter according to the azimuth phase ;

[0104] To compensate The non-empty variables in the compensation filter are obtained as follows

[0105] H 1 ( f a ) = exp [ − j ( G 03 f a 3 + G 04 f a 4 ) ]

[0106] In addition, more azimuth degrees of freedom need to be introduced to achieve the processing of azimuth space-variable phase. Therefore, for the signal processed by the compensation filter, multiply it by the following variable scaling function ,as follows

[0107] H 2 ( f a ) = exp [ j π ( p 2 f a 2 + p 3 f a 3 + p 4 f a 4 ) ]

[0108] in, , and Represents the coefficient of the azimuth equalization function and needs to be solved by the subsequent equation. Since the azimuth spatial variable coefficient needs to be equalized, it is necessary to multiply the variable scaling function in the time domain. However, in the process of converting to the time domain, aliasing cannot occur in the time domain, otherwise when the new variable scaling function is introduced, additional phase errors will be caused, so the following constraints need to be met:

[0109]

[0110] in, is a constant, which determines the envelope center of the azimuth echo in the time domain, so it is necessary to ensure .in Indicates the number of azimuth signal points before keystone transformation. Indicates the azimuth sampling rate after keystone transformation; It can be expressed as

[0111]

[0112] Solving the equation, we get

[0113]

[0114] Step S13: Perform an inverse Fourier transform in azimuth to obtain the azimuth time domain echo phase ;

[0115] Convert the echo to the time domain to obtain the azimuth time domain echo phase , which is expressed as follows

[0116]

[0117] in, , , , represents the time-domain non-space-variant azimuth frequency modulation coefficient and can be expressed as follows

[0118]

[0119] also, Represents the spatial variation coefficients of various orders in azimuth.

[0120] Step S14: Obtain and multiply the compensation function according to the azimuth time domain echo phase ;

[0121] In order to remove the non-empty variables in the time domain, the compensation function is obtained as follows

[0122] H 3 ( t m ) = exp [ − j ( G 03 ′ t m 3 + G 04 ′ t m 4 + G 05 ′ t m 5 ) ]

[0123] Step S15: Multiply by NLCS equalization function ;

[0124] At this point, an additional azimuth equalization factor is further introduced and multiplied by the NLCS equalization function ,as follows

[0125] H 4 ( t m ) = exp { j π [ m 2 t m 2 + m 3 t m 3 + m 4 t m 4 ] }

[0126] in , , They represent the equalization factors of each order of NLCS respectively.

[0127] Step S16: Azimuth FFT to obtain the azimuth frequency domain phase ;

[0128] After introducing enough equalization factors, the signal is converted to the frequency domain, and the azimuth frequency domain phase is obtained:

[0129]

[0130] in, represents the frequency domain non-space-varying azimuth frequency modulation coefficient and can be expressed as follows

[0131]

[0132] in, express , express In addition, in order to balance the high-order space-variant phase, the following equation needs to be satisfied:

[0133]

[0134] in, , , , and Represent the spatial variation coefficients of each order respectively. Solving the above equation, we can get

[0135]

[0136]

[0137] in, express .for , and , whose expression is relatively complex and can be derived through numerical derivation software. is a constant that determines the envelope center of the azimuth echo in the frequency domain.

[0138] Step S17: Obtain and multiply the azimuth pulse pressure function according to the azimuth frequency domain phase ;

[0139] The azimuth pulse pressure is realized by azimuthally consistent compression. The azimuth pulse pressure function is as follows:

[0140]

[0141] Step S18: Output the imaging result.

[0142] Perform inverse Fourier transform on the azimuth signal to obtain the imaging result.

[0143] The corresponding two-dimensional resolution for this simulation case is shown in Figure 4 As shown, Figure 4 (a) is the azimuth resolution, Figure 4 (b) in the figure is the range resolution. The azimuth resolution and range resolution are about 3m, and the imaging width is 6km×1.5km.

[0144] The pulse repetition frequency of low-orbit satellites is usually 3000~4000Hz, and 4000Hz is used in this simulation. According to the Doppler analysis of the satellite-aircraft bistatic SAR combination mode, the total azimuth bandwidth corresponding to the point target scene in this case is 9035Hz, which is much larger than the pulse repetition frequency of the transmitted signal. Therefore, azimuth preprocessing and other methods are needed to achieve azimuth signal balance. The specific processing results are shown in Figure 5 shown. Figure 5 (a) shows the azimuth spectrum without de-skewing. The Doppler centroid of the target point at [1.5km, 1.5km] has exceeded the azimuth sampling rate, resulting in aliasing problems. Figure 5(b) shows the result of the preprocessing method proposed in this paper. At this time, the azimuth spectrum width has been compressed to a large extent, about 250Hz. Figure 5 (c) shows the spectrum result in azimuth when the de-slant factor for the satellite-aircraft bistatic SAR combination mode is not accurately calculated. At this time, the azimuth spectrum cannot be completely compressed, resulting in a wide azimuth spectrum. Although in the current simulation case, the pulse repetition frequency is 4000Hz, the bandwidth after processing is about 1250Hz, and there is no aliasing problem, in other cases, there is a risk of azimuth aliasing, which affects the robustness of the imaging algorithm. Therefore, accurate calculation of the de-slant factor is a key factor in frequency domain imaging of the satellite-aircraft bistatic SAR combination imaging mode.

[0145] The processing method proposed in this paper also optimizes the range processing flow. Since the interpolation operation of the keystone processing method cannot be applied to undersampled signals, this method combines the azimuth preprocessing and the keystone method and compensates for the phase error introduced by the process. The comparison results of this processing flow and the unoptimized range processing are given. Figure 6 (a) is the distance processing result of the signal after keystone transformation; Figure 6 (b) The range processing result of the signal after azimuth preprocessing + keystone + preprocessing compensation; Figure 6 (c) The distance processing result of the signal after the keystone transformation without preprocessing; Figure 6 (d) The distance processing result of the signal after keystone+ preprocessing; Figure 6 (e) The distance processing result of the signal after keystone+preprocessing+HRCMC+SRC processing; Figure 6 (a) and (b) in the figure show the range processing results of the azimuth preprocessing, keystone, HRCMC, SRC and preprocessing compensation mentioned in this paper. The azimuth preprocessing eliminates the azimuth spectrum aliasing problem, see Figure 5 Therefore, after the keystone transformation, accurate data mapping results can be obtained. Figure 6 (a) in the figure. Figure 6 (c), (d), and (e) in the figure show the unoptimized distance processing results. Figure 6 From (a) and (c), we can see that the azimuth spectrum aliasing problem has a significant impact on the keystone interpolation effect and will introduce additional blur energy, which is divided into two aspects: range and azimuth. The azimuth blur energy is reflected in the envelope distortion and phase error of the target point corresponding to the spectrum folding in the time domain. The range blur energy is reflected in the false target generated by the target point along the range direction at the spectrum folding. Figure 5 (b) and Figure 6 In (d), the blur energy introduced by the keystone transform will further interfere with the results of the orientation preprocessing. Figure 6 It can be seen from (b) and (e) that compared with the unoptimized range processing flow, the optimized range processing flow proposed in this paper can achieve accurate range migration correction and azimuth spectrum regularization.

[0146] Finally, the corresponding point target imaging result is as follows Figure 7 In addition, Figure 8 A contour plot of a dot matrix target is given, where Figure 8 (a), (b), and (c) correspond to the contour maps of point P1, point O, and point P2, respectively. Fig. 9 The azimuth profile of the dot matrix target is given, where Fig. 9 (a), (b), and (c) in the figure correspond to the azimuth profiles of point P1, point O, and point P2, respectively. It can be seen from the contour map and the azimuth profile that the method proposed in this paper has good imaging performance at the center and edge of the imaging result, and the imaging effect is good.

[0147] In summary, the present invention proposes a satellite-aircraft bistatic SAR combined mode imaging method, which avoids the phase error introduced by azimuth spectrum aliasing in the range processing process by improving the range processing flow; improves the azimuth processing flow, while avoiding aliasing in the time-frequency domain of azimuth processing, introduces more degrees of freedom, and balances the azimuth high-order space-variation; proposes an accurate combined mode satellite-aircraft bistatic SAR Doppler space-variation parameter calculation method, which can obtain high-precision de-skew factor and NLCS equalization factor, and realizes efficient and high-precision processing of the combined mode satellite-aircraft bistatic SAR.

[0148] Those skilled in the art will appreciate that the above embodiments are intended to help readers understand the principles of the present invention, and should be understood that the scope of protection of the present invention is not limited to such specific statements and embodiments. For those skilled in the art, the present invention may have various changes and variations. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of the claims of the present invention.

Claims

1. A satellite-aircraft dual-base synthetic aperture radar combined mode imaging method, characterized in that: The following steps are involved: Step S1: input echo signals, geometric parameters, and beam steering parameters to construct a satellite-machine bistatic synthetic aperture radar SAR combination mode geometric model and signal model; Step S2: Solve the Doppler space-variable coefficient. The specific steps are as follows: For satellite-to-aircraft bistatic SAR, the point target echo model in imaging mode is expressed as: Among them, σ(x,y) represents the scattering coefficient of the target point, w a Indicates the azimuth signal envelope; B r represents bandwidth, T r represents pulse width, τ represents distance, λ represents wavelength, c represents speed of light, R b (t; x, y) represents the distance history of the two bases, j represents the imaginary number sign; t represents the azimuth time, T a represents the synthetic aperture time, t c Indicates the central irradiation time; the dual-base distance history R b (t; x, y) is obtained by the following formula: R b (t;x,y)=k0+k1t+k2t 2 +k3t 4 ; Among them, k i is the Doppler spatial variation coefficient of each order, i = 1, 2…; the Doppler spatial variation coefficient k i Performing Taylor expansion, we get: Among them, a i , b i 、c i and d i Represent the first, second, third, and fourth order Doppler parameter space-varying coefficients, i = 1, 2…, t R is the azimuth time, t R =0 when k i The value of ; at this time, Among them, t c represents the central irradiation time, x represents the x-axis position of the target point, y represents the y-axis position of the target point, and satisfy: Among them, t c represents the central irradiation time, R bi,0 Indicates the illumination time t of different target points at their azimuth center c The corresponding double base distance and is expressed as: R bi,0 =R b (t c (x,y);x,y); In addition, the spatial variation coefficients of k1 are expressed as: Similarly, the other order space-varying coefficients b i 、c i and d i Using the above solution a i Calculate it in the following way; Step S3: Perform a distance Fourier transform to obtain the signal spectrum of the original echo, and then compensate the spectrum, specifically including pulse compression, spectrum tilt correction, and spectrum regularization; the details are as follows: The signal spectrum of the original echo is expressed as: Among them, W r [·] represents the range signal envelope, w a Indicates the azimuth signal envelope; B r represents bandwidth, f τ Represents the distance frequency, f c is the carrier center frequency, k r represents the frequency modulation slope in the range direction; t represents the time in the azimuth direction, T a represents the synthetic aperture time, t c Indicates the central irradiation moment; Then the spectrum is multiplied by the compensation phase function to achieve pulse compression and compensate the phase. The compensation phase function is expressed as follows: Construct a linear distance movement correction function and multiply it with the spectrum to correct the spectrum tilt; the linear distance movement correction function is specifically expressed as follows: Construct a spectrum regularization function and multiply it with the spectrum to achieve spectrum regularization; the spectrum regularization function is as follows: H de (t)=exp[-jπK rot t 2 ]; Among them, K rot is the deskewing coefficient, expressed as Finally, the corrected signal is obtained, which is expressed as follows: Step S4: Performing Keystone transformation on the spectrum corrected in step S3 to obtain a transformed signal; Step S5: After performing high-order range motion correction, secondary range compression, and additional phase compensation on the transformed signal, perform inverse Fourier transform in the range direction to obtain signal S d (τ,t m ); Step S6: Azimuth Fourier transformation is performed to obtain the azimuth phase; then the non-space-variant term of the azimuth space-variant coefficient is compensated; and then the azimuth degree of freedom is introduced by multiplying the variable scaling function; Step S7: Perform an inverse Fourier transform in azimuth to obtain the azimuth time domain echo phase; Step S8: Compensate the azimuth time domain echo phase to remove the non-space variable term in the time domain; then multiply it by the NLCS equalization function; Step S9: perform Fourier transform in azimuth to obtain frequency domain phase in azimuth; Step S10: Obtain and multiply the azimuth pulse pressure function H according to the azimuth frequency domain phase az (f a ); Step S11: Perform inverse Fourier transform on the azimuth signal to obtain an imaging result.

2. The satellite-aircraft dual-base synthetic aperture radar combined mode imaging method according to claim 1, characterized in that: The Keystone transformation is as follows: Perform keystone transformation, the specific interpolation method is as follows: t m represents the azimuth time variable after keystone transformation; the transformed signal is expressed as: Among them, t c,m represents the azimuth center illumination time corresponding to different target points in the azimuth time domain after the keystone transformation; T a,m Represents the synthetic aperture time in the azimuth time domain after the keystone transform.

3. The satellite-aircraft dual-base synthetic aperture radar combined mode imaging method according to claim 2, characterized in that: The step S5 is specifically as follows: The transformed signal is multiplied by a high-order range migration correction function, a quadratic range compression function, and an additional compensation phase function; specifically, the high-order range migration correction function H HRCMC (f τ ) is expressed as follows: Quadratic distance compression function H SRC (f τ ) is expressed as follows: Additional compensation phase H com (f τ ) is expressed as follows: Then perform inverse Fourier transform to obtain the signal S d (τ,t m ); let f τ =0, the signal is 4. The satellite-aircraft dual-base synthetic aperture radar combined mode imaging method according to claim 3 is characterized in that: The step S6 is specifically as follows: Fourier transform to frequency domain to get azimuth phase for: Among them, f a represents the azimuth frequency variable, t R is the azimuth time, G 01 , G 02 , G 03 and G 04 represents the non-space-variable azimuth frequency modulation coefficient, G ij (i,j∈Z,i≠0) represents the spatially variable coefficients of various orders in azimuth; Then the non-space-variant term of the azimuth space-variant coefficient is compensated by multiplying the following compensation filter H1(f a ): H1(f a )=exp[-j(G 03 f a 3 +G 04 f a 4 )]; Then multiply by the following variable scaling function H2(f a ),as follows: Among them, p2, p3 and p4 represent the coefficients of the azimuth equalization function.

5. The satellite-aircraft dual-base synthetic aperture radar combined mode imaging method according to claim 4, characterized in that: The azimuth time domain echo phase is expressed as follows: Among them, G 02 ′, G 03 ′, G 04 ′, G 05 ′ represents the time domain non-space-variant azimuth frequency modulation coefficient, which is expressed as follows: Among them, G ij ′(i,j∈Z,i≠0) represents the spatially varying coefficients of various orders in azimuth.

6. The satellite-aircraft dual-base synthetic aperture radar combined mode imaging method according to claim 5, characterized in that: The step S8 is specifically as follows: When removing the non-empty variables in the time domain, multiply by the following compensation function H3(t m ): H3(t m )=exp[-j(G 03 ′t m 3 +G 04 ′t m 4 +G 05 ′t m 5 )]; The nonlinear variable-scale NLCS equalization function is expressed as follows: Among them, m2, m3, and m4 represent the equalization factors of each order of nonlinear scaling respectively.

7. The satellite-aircraft dual-base synthetic aperture radar combined mode imaging method according to claim 6 is characterized in that: The azimuth frequency domain phase is expressed as follows: Among them, G 01 ″、G 02 ″、G 03 ″、G 04 ″ represents the frequency domain non-space-variable azimuth frequency modulation coefficient, G 02 ″、G 03 ″、G 04 ″ means as follows Where ξ1 represents πp2m2-G 02 m2-π, ξ2 represents πp2+G 02 m2; In addition, the frequency domain spatially varying coefficients of each order need to satisfy the following equations: Among them, G 11 ″、G 21 ″、G 12 ″、G 22 ″ and G 13 ″ respectively represent the spatial variable coefficients of each order; solving the above equation, we get: Among them, ψ1 represents p3, m4 and p5 are obtained by numerical derivation software; γ is a constant; Then, according to the azimuth frequency domain phase, we can obtain and multiply it by the azimuth pulse pressure function H az (f a ):

Images