A star-machine bistatic SAR imaging method in a bidirectional sliding spotlight mode

Abstract

The application discloses a kind of star-machine bistatic SAR frequency domain imaging methods under two-way sliding bunching mode, for the low PRF star-machine bistatic SAR system, first, the echo signal is preprocessed, combined with the idea of spectrum analysis technique, using equivalent convolution transform and compensation to remove the azimuth spectrum aliasing of signal;Subsequently, the method of Legendre polynomial expansion is used to approximate the slant range, and a high-precision slant range model is established;Then, the idea of series inversion is used to obtain the accurate two-dimensional spectrum phase analytical expression of point target;Finally, the algorithm is designed in frequency domain to obtain high-resolution image.The implementation process is: (1) establish the slant range expression of star-machine bistatic SAR;(2) obtain point target echo data;(3) carry out dealiasing preprocessing to point target echo;(4) use Legendre polynomial to approximate slant range;(5) use Legendre polynomial expansion to obtain spectrum phase analytical expression;(6) design frequency domain imaging algorithm to obtain imaging result.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of radar signal processing methods, in particular to a star-machine bistatic SAR imaging method under a bidirectional sliding beamforming mode, which has an improving effect on the imaging quality of scene edge points. BACKGROUND

[0002] Due to the separation of the transmitting and receiving platforms, the bistatic SAR system can be applied to various configurations, and thus has higher flexibility and practicality than the monostatic SAR system, and can obtain more abundant information. The star-machine bistatic SAR system has the characteristics of long range and wide coverage, and the advantages of strong anti-interference ability and high security. However, due to the separation of the transmitting and receiving platforms, the bistatic SAR system also has corresponding technical problems.

[0003] Firstly, the problem of slant range approximation. Unlike the traditional monostatic SAR, the total slant range of the bistatic SAR is the sum of the corresponding slant ranges of the transmitting and receiving platforms to the target point, which presents the form of the sum of double square roots, and thus has higher complexity in system and echo characteristics. In view of the approximation of the two-way slant range, scholars at home and abroad have proposed various methods to solve this problem. In document [1], Xiong Tao, Li Yachao, Li Qi, etc. proposed to approximate the two-way slant range by combining the monostatic equivalent component and the bistatic compensation component, wherein the monostatic equivalent component includes three equivalent parameters of equivalent slant range, velocity and squint angle. In document [2], Li Dong, Liao Guisheng, etc. approximated the slant range by using the Taylor series expansion method. In document [3], Hua Zhong, Guangyong Zheng, Ronghua Zhao, etc. proposed a derivation scheme based on the distance equivalent ellipse model of the high squint BiSAR focusing. In document [4], Li Menghui, Tan Gewei, etc. proposed to approximate the slant range by using the fourth-order Chebyshev polynomial, but only simulated and verified the airborne bistatic SAR, and did not consider the possibility and applicability under more configurations (such as the star-machine bistatic SAR, etc.).

[0004] Secondly, the design of imaging algorithm is a problem. The signal processing of bistatic SAR is more complex, and the imaging method of monostatic SAR is no longer directly applicable. The FBP algorithm mentioned in the literature [5] is a typical representative of time-domain imaging algorithm, which is considered as the preferred algorithm for bistatic SAR imaging. This algorithm can be used for any bistatic SAR, and there is no restriction on the configuration and flight trajectory in theory. However, the disadvantages are large amount of calculation, low computational efficiency, and great influence of platform motion error and inertial navigation measurement error. In the frequency domain imaging, the first task is to derive the two-dimensional spectrum of the target. The current mainstream spectrum solving method is the method of series reversion (MSR) method, which uses Taylor series high-order approximation method to solve the high-precision two-dimensional spectrum of any configuration bistatic SAR point target, so it is the most widely used. On the basis of obtaining the two-dimensional spectrum, the selection of RD, NCS, OMEGA-K and other classic algorithms can complete the imaging task of the target. In the star-machine bistatic SAR system, due to the huge difference in the speed of the transmitting and receiving platforms, there may be spectrum blurring problem at the edge of the imaging scene, which affects the imaging quality. Therefore, the step of eliminating spectrum blurring must be considered in the design of frequency domain imaging algorithm. The current solution mainly includes sub-aperture segmentation in the literature [6] and full-aperture processing based on upsampling in the literature [7], [8].

[0005] The following gives the relevant literature for retrieval in the background art:

[0006] [1] Xiong Tao, Li Yachao, Li Qi, et al. Using an Equivalence-Based Approach to Derive 2-D Spectrum of BiSAR Data and Implementation Into an RDA Processor [J]. IEEE Transactions on Geoscience and Remote Sensing, 2021, 59(6): 4765-4774.

[0007] [2] Li Dong, Liao Guisheng, Wang Wei, et al. Extended azimuth nonlinear chirp scaling algorithm for bistatic SAR processing in high-resolution highly squinted mode[J]. IEEE Geoscience and Remote Sensing Letters, 2014, 11(6): 1134-1138.

[0008] [3] Hua Zhong, Guangyong Zheng, Ronghua Zhao, et al. Focus Improvement for Highly Squinted One-Stationary BISAR Imaging Based On A Range Equivalent Model[C]. IGARSS 2019.

[0009] [4] Li Menghui, Tan Ewei, Yang Jingjing, et al. Bistatic SAR imaging algorithm based on motion compensation and orthogonal decoupling[J]. Signal Processing, 2021, 37(1): 75-85.

[0010] [5] Shi Tianyue, Mao Xinhua, Andreas Jakobsson, Liu Yanqi. Parametric Model-Based 2-D Autofocus Approach for General BiSAR Filtered Backprojection Imagery[J]. IEEE Transactions on Geoscience and Remote Sensing, 2022, 60: 1-14.

[0011] [6] Zhang Jingdong, Chen Jiarui, Zhu Daoyin, et al. A sub-aperture processing imaging method for squint sliding spotlight SAR[J]. Data Acquisition and Processing, 2017, 32(04): 776-784.

[0012] [7] Guan Yifu, Yang Yuan, Jia Wentong, et al. Research on airborne TOPS mode two-step imaging algorithm[J]. Journal of Air Force Early Warning Academy, 2020, 34(03): 167-172.

[0013] [8] LUO X L, DENG Y K, WANG R, et al. Image formation procfessing for sliding spotlight SAR with stepped frequency chirps[J]. IEEE Geoscience and Remote Sensing Letters, 2014, 11(10): 1692-1696. SUMMARY

[0014] In view of the problems in the prior art, i.e., the slant range approximation accuracy and the design of the space-airborne bistatic SAR imaging algorithm, the application provides a space-airborne bistatic SAR imaging method in a bidirectional sliding spotlight mode, so as to realize unfuzzy imaging of scene edge point targets under a low pulse repetition frequency (PRF) condition, and the effect is better than that of a traditional frequency domain imaging algorithm based on Taylor series expansion.

[0015] To achieve the above object, the technical scheme adopted by the application is as follows:

[0016] A space-airborne bistatic SAR imaging method in a bidirectional sliding spotlight mode, comprising the following steps:

[0017] Step 1: establishing a slant range expression of the bidirectional sliding spotlight mode according to a geometric model of the space-airborne bistatic SAR system;

[0018] Step 2: obtaining a ground point target echo signal generated when a low-orbit satellite platform transmits a linear frequency modulation signal;

[0019] Step 3: performing dealiasing preprocessing on the point target echo to obtain an unfuzzy spectrum;

[0020] Step 4: approximating the slant range expression established in step 1 by using a Legendre polynomial to obtain a Legendre expansion of the two-way slant range;

[0021] Step 5: combining a series inversion method, using the Legendre expansion of the two-way slant range in step 4 to obtain a phase ψ(f r ,f a ) of a point target original two-dimensional spectrum, expanding the phase ψ(f r ,f a ) into a power series form in a range frequency domain f r , and discarding high-order terms in the range frequency domain f r and an azimuth frequency domain f a to obtain a two-dimensional spectrum phase based on the Legendre polynomial expansion;

[0022] Step 6, combining each phase term of the two-dimensional frequency spectrum based on the Legendre polynomial expansion obtained in step 5, performing distance compression and secondary compression, distance migration correction, constant term and high-order term compensation and azimuth compression on the unambiguous spectrum obtained in step 3, and obtaining an imaging result.

[0023] Further, in step 3, the process of de-aliasing preprocessing of the point target echo is as follows:

[0024] 3a), constructing a reference function:

[0025] First, the reference function of the azimuth time domain τ is constructed as shown in formula (1):

[0026] H ref (τ)=exp(jπητ 2 ) (1),

[0027] Wherein, η is a frequency modulation coefficient, which is related to the actual configuration parameters of the star-machine bistatic SAR system, and satisfies the relationship λ is the wavelength of the carrier, v T and v R are the moving speeds of the transmitting platform and the receiving platform respectively, R tef and R ref are the distances from the transmitting platform and the receiving platform to the virtual rotation center point in the two-way sliding spotlight mode.

[0028] Step 2, performing distance FFT on the time domain echo signal S(t,τ) obtained in step 2 to convert it into the distance frequency domain to obtain S(f r ,τ), and multiplying S(f r ,τ) with the reference function H ref (τ);

[0029] 3b), equivalent convolution transform:

[0030] Performing azimuth FFT on the result of step 3a), and multiplying it with the reference function H ref (τ') under the new time domain coordinate τ' to obtain the unambiguous signal S(f r ,τ') in the azimuth frequency domain under the new time domain coordinate τ'. This step can also be regarded as equivalent convolution transform of the time domain echo S(f r ,τ) and the reference function H ref (τ) to obtain S(f r ,τ'), wherein:

[0031] The FFT transform is equivalent to the conversion between the time-frequency axes. Before and after the FFT transform, the original time domain coordinate τ of the echo signal, the new time domain coordinate τ', the original frequency coordinate f a of the echo signal, and the new frequency coordinate f′ a satisfy the relationship

[0032] 3c) Conjugate compensation:

[0033] For S(f) obtained in step 3b), r Perform an azimuth FFT again to obtain a new two-dimensional frequency domain signal S(f). r ,f′ a ), for S(f r ,f′ a Multiply by the frequency domain conjugate compensation function H ref * (f′ a The non-aliased spectrum is obtained as shown in formula (2):

[0034] S'(f r ,f′ a )=S(f r ,f′ a )·H ref * (f′ a (2),

[0035] In further step 4, the process of obtaining the two-dimensional spectrum expression based on Legendre polynomial expansion is as follows:

[0036] 4a) Normalize the slow time τ in the direction of orientation, that is, we have T a The time for synthesizing the aperture is [time]. Then, the slope distance R is [time]. bi (τ) is decomposed using Legendre orthogonal decomposition and expanded into a fourth-order power series of x as shown below:

[0037] R bi (x)=α0+α1x+α2x 2 +α3x 3 +α4x 4 (3),

[0038] Where the corresponding term coefficient α i satisfy:

[0039]

[0040] Where L n (x) is a Legendre polynomial that satisfies the formula

[0041] In equation (4), the operation within the curly braces {} is a definite integral operation. The result can be calculated by substituting the corresponding order n. The specific operation content satisfies the formula:

[0042] 4b) will In the case of return (3), the double-path slant range R is obtained by reorganizing bi (τ) is the fourth power series expansion of the slow time τ, i.e. the Legendre polynomial approximation of the slant range is:

[0043] R bi (τ) = k0+ k1τ+ k2τ 2 + k3τ 3 + k4τ 4 (5),

[0044] where, are the Legendre decomposition coefficients of the slant range.

[0045] Further, in step 5, the process of obtaining the two-dimensional spectrum expression based on the Legendre polynomial expansion is as follows:

[0046] 5a), first, the azimuth time-frequency transform of the de-linear phase is completed by using the stationary phase principle, and the result is shown in equation (6):

[0047]

[0048] Then, the corresponding relationship between the azimuth slow time and the frequency is obtained by series inversion method:

[0049]

[0050] In equations (6) and (7), k n (n = 2, 3, 4) are the Legendre decomposition coefficients of the slant range, f r is the range frequency, f a is the azimuth frequency, f c is the carrier frequency, and c is the speed of light.

[0051] 5b), the linear phase removed before is recovered, and the original two-dimensional spectrum can be obtained by using the frequency shift property of Fourier transform, as shown in equation (8):

[0052] S(f r ,f a ) = W r (f r ) W a (f a ) exp(jψ(f r ,f a )) (8),

[0053] where, W r represents the range envelope, W a represents the azimuth envelope, and ψ(f r ,f a ) represents the original two-dimensional spectrum phase, and the specific expression is shown in equation (9):

[0054]

[0055] where k n (n = 0, 1, 2, 3, 4) are the corresponding Legendre decomposition coefficients of slant range, f r is the range frequency, f a is the azimuth frequency, f c is the carrier frequency, and c is the speed of light.

[0056] The expression of (9) is expanded by using the Legendre polynomial. First, the normalization processing is performed, i.e., let B is the signal bandwidth. The three fractional expressions are respectively arranged into the third power series of y:

[0057]

[0058] where the coefficients in front of y of different power terms satisfy the relationship as shown in equation (11):

[0059]

[0060] where the operation process in the braces {} is the definite integral operation, and the corresponding order n is brought into the calculation result during the operation. The specific operation content satisfies the formula: Then, the expression of (10) is brought into the power series of y, and is arranged into the third power series of f r , to obtain equation (12):

[0061]

[0062] where is the corresponding Legendre decomposition coefficient of f .

[0063] 5c), the result of equation (12) is brought into equation (9) to arrange and combine, and the high-order terms of f r and f a are discarded, so that the expression of the Legendre spectrum phase is shown in equation (13):

[0064]

[0065] where the different power terms of the range frequency f r satisfy the following relationship:

[0066]

[0067]

[0068] where k​​n (n = 0, 1, 2, 3, 4) are the slant range corresponding Legendre decomposition coefficients, f r is the range frequency, f a is the azimuth frequency, f c is the carrier frequency, c is the speed of light.

[0069] Observing each phase term in equation (13), we find that, is independent of f r and f a , is a constant phase term; is the linear term coefficient of f r , corresponding to the range migration; is the quadratic term coefficient of f r , corresponding to the range compression; is the cubic term coefficient; is only related to f a , is the azimuth modulation term, corresponding to the azimuth compression.

[0070] Further, in step 6, the alias-free two-dimensional spectrum S'(f r ,f′ a ) is subjected to range compression and quadratic compression, migration correction, constant term and high-order term compensation, and azimuth compression, in sequence. The specific process is as follows:

[0071] 6a), first, construct the range compression and quadratic compression function as shown in equation (16), where γ is the frequency modulation:

[0072]

[0073] Secondly, perform the range migration correction, and construct the compensation function as shown in equation (17):

[0074]

[0075] Then, perform the constant term and high-order term compensation, and the compensation function is as shown in equation (18):

[0076]

[0077] 6b), after the above-mentioned range compression, migration correction, and constant term and high-order term compensation, the compensated two-dimensional frequency domain signal S'(f r ,f′ a ) is subjected to range IFFT, and the signal is converted back to the time domain to obtain S'(t,f′ a ). Then, combining the idea of Dechirp processing, construct the azimuth compression function to convert the phase of the IFFT signal, as shown in equation (19)

[0078]

[0079] For S'(t,f′) a ) and H az The result of multiplication is used to perform azimuth IFFT to obtain an unambiguous azimuth time domain signal, which is then multiplied by the Decirp function to obtain formula (20):

[0080]

[0081] In equations (19) and (20), ξ satisfies the relation η satisfies the relation λ is the wavelength of the carrier wave, v T and v R These represent the movement speed of the transceiver platform, R tef and R ref R represents the distance from the transceiver platform to the virtual rotation center point in bidirectional sliding beamforming mode. Tcen and R Rcen These are the minimum slant distances from the transceiver platform to the ground target.

[0082] 6c) Finally, the azimuth time-domain unambiguous signal after the above series of operations is transformed back to the frequency domain by azimuth FFT to obtain an unambiguous image.

[0083] Compared with the prior art, the advantages of the present invention are:

[0084] (1) The traditional method of approximating slant range using Taylor series approximates the slant range near the center of the synthetic aperture by using the value of the slant range model at that moment. When the azimuth time is far from the center time, the slant range approximation error of the target point will increase. In contrast, this invention uses Legendre polynomial expansion to perform a fourth-order approximation of the slant range expression of the satellite-aircraft bistatic SAR. This method uses the instantaneous slant range at the corresponding azimuth time for weighted calculation, which has higher accuracy and smaller error compared to the traditional Taylor series expansion approximation method.

[0085] (2) This invention designs an echo dealiasing preprocessing scheme under low PRF conditions. Due to the huge difference in flight speed between the low-orbit satellite launch platform and the UAV receiving platform, there is a significant difference in the actual PRF required by the two. For the satellite-UAV dual-base SAR system under low PRF conditions, the Doppler bandwidth of the target point located at the edge of the imaging scene is significantly greater than the PRF. If conventional frequency domain imaging methods are used directly for processing, spectral aliasing will inevitably occur. However, this invention first preprocesses the echo, and through equivalent convolution transformation and compensation, makes the preprocessed PRF greater than the Doppler bandwidth. This scheme not only eliminates the azimuth spectral ambiguity, but also avoids the excessive echo data volume and range ambiguity caused by directly expanding the PRF, thereby enabling efficient frequency domain ambiguity-free imaging in subsequent algorithm design.

[0086] (3) The application expands the frequency domain imaging application range of using Legendre polynomial to expand the target spectrum. The application background of the research of the Legendre polynomial expansion proposed in the prior art is all airborne bistatic SAR systems. The application expands the bidirectional sliding spotlight working mode under the more complex space-airborne bistatic SAR system, and overcomes the spectrum ambiguity problem of the space-airborne bistatic SAR system. The working mode can realize the coordination of the mapping width and the resolution by setting a reasonable sliding coefficient.

[0087] (4) The application designs a set of processing flow for the space-airborne bistatic SAR system imaging of low PRF. Firstly, the echo signal is preprocessed, the equivalent convolution transform and compensation are used in combination with the idea of spectrum analysis technology to eliminate the azimuth aliasing of the signal; then the two-way slant range is approximated by using the Legendre polynomial expansion method to establish a high-precision slant range model; subsequently, the accurate two-dimensional spectrum phase expression of the point target is obtained in combination with the idea of series inversion; finally, the imaging result is obtained by designing the algorithm in the frequency domain. BRIEF DESCRIPTION OF DRAWINGS

[0088] Figure 1 The flow chart of the frequency domain imaging algorithm proposed in the embodiment.

[0089] Figure 2 The schematic diagram of the bistatic SAR geometric model.

[0090] Figure 3 The schematic diagram of the bidirectional sliding spotlight working mode of the space-airborne bistatic SAR, wherein (a) is the satellite forward sliding spotlight, and (b) is the unmanned aerial vehicle reverse sliding spotlight.

[0091] Figure 4 The simulation schematic diagram of the ground 5*5 uniformly distributed point target array.

[0092] Figure 5 The slant range error comparison diagram under different slant range approximation methods.

[0093] Figure 6 The two-dimensional spectrum comparison diagram before and after echo de-aliasing preprocessing, wherein (a) is the two-dimensional spectrum before de-aliasing, and (b) is the two-dimensional spectrum after de-aliasing.

[0094] Figure 7 The scene different position imaging point height diagram comparison, wherein (a) is the imaging height diagram of the center point P1 under the algorithm of the application, (b) is the imaging height diagram of the edge point P3 under the algorithm of the application, (c) is the imaging height diagram of the center point P3 under the Taylor series expansion method, and (d) is the imaging height diagram of the edge point P3 under the Taylor series expansion method.

[0095] Figure 8Distance and azimuth profile of scene edge imaging point P3 under different imaging algorithms, wherein: (a) is the distance profile of scene edge point P3 under the algorithm of the application; (b) is the azimuth profile of scene edge point P3 under the algorithm of the application; (c) is the distance profile of edge point P3 under the Taylor series expansion frequency domain algorithm; (d) is the azimuth profile of edge point P3 under the Taylor series expansion frequency domain algorithm. DETAILED DESCRIPTION

[0096] The application will be further described below in conjunction with the accompanying drawings and examples.

[0097] Example 1

[0098] This embodiment discloses a star-machine dual-basis SAR two-dimensional frequency domain imaging method combining echo preprocessing and Legendre polynomials, as shown in Figure 1 , comprising the following steps:

[0099] Step 1: Establishing the slant range expression of the bidirectional sliding spotlight mode according to the geometric model of the star-machine dual-basis SAR system.

[0100] The schematic diagram of the bidirectional sliding spotlight mode is shown in Figure 3 , wherein W and U are the virtual rotation centers of the transmitter and receiver respectively, R tef and R ref are the distances from the transmitting and receiving platforms to the rotation center, R tcen and R rcen are the shortest action distances from the transmitter and receiver to the ground.

[0101] The geometric model of the dual-basis SAR system is shown in Figure 2 , wherein point P represents an arbitrary scattering point in the scene, with coordinates (x0, y0, 0), and the three-dimensional real-time positions of the transmitter and receiver radar platforms are (x t , y t , z t ) and (x r , y r , z r ) respectively. The slant range expression R bi (τ) under the bidirectional sliding spotlight working mode is established according to the geometric model of the star-machine dual-basis SAR system as follows:

[0102]

[0103] Step 2: The ground point target echo signal generated when the low-orbit satellite platform transmits the linear frequency modulation signal is shown as follows:

[0104]

[0105] where σ is the backscattering coefficient of the point target, t and τ are the range fast time and azimuth slow time, respectively, W r and W a are the range window function and azimuth window function, respectively, c is the light speed, λ is the carrier wavelength, and γ is the transmitted signal modulation frequency.

[0106] Step 3, de-aliasing preprocessing is performed on the point target echo to obtain a non-ambiguous spectrum. The specific process is as follows:

[0107] 3a), a reference function is constructed:

[0108] First, the reference function of the azimuth time domain τ is constructed as follows:

[0109] H ref (τ)=exp(jπητ 2 ) (1),

[0110] where η is the modulation frequency coefficient, which is related to the actual configuration parameters of the space-aircraft bistatic SAR system and satisfies the relationship λ is the wavelength of the carrier, v T and v R are the moving speeds of the transmitting and receiving platforms, R tef and R ref are the distances from the transmitting and receiving platforms to the virtual rotation center point in the two-way sliding spotlight mode.

[0111] The time domain echo signal S(t,τ) obtained in step 2 is subjected to a range FFT to convert it into the range frequency domain to obtain S(f r ,τ), and S(f r ,τ) is multiplied by the reference function H ref (τ);

[0112] 3b), equivalent convolution transform:

[0113] The result of step 3a) is subjected to an azimuth FFT, and multiplied by the reference function H ref (τ') in the new time domain coordinate τ', to obtain the non-ambiguous signal S(f r ,τ') in the azimuth frequency domain under the new time domain coordinate τ'. This step can also be regarded as an equivalent convolution transform of the time domain echo S(f r ,τ) and the reference function H ref (τ) to obtain S(f r ,τ'), where:

[0114] The FFT transform is equivalent to the conversion between the time-frequency axes. Before and after the FFT transform, the original time domain coordinate τ of the echo signal, the new time domain coordinate τ', the original frequency coordinate f a of the echo signal, and the new frequency coordinate f′ a satisfy the relationship

[0115] 3c) Conjugate compensation:

[0116] For S(f) obtained in step 3b), r Perform an azimuth FFT again to obtain a new two-dimensional frequency domain signal S(f). r ,f′ a ), for S(f r ,f′ a Multiply by the frequency domain conjugate compensation function H ref * (f′ a The non-aliased spectrum is obtained as shown in formula (2):

[0117] S'(f r ,f′ a )=S(f r ,f′ a )·H ref * (f′ a (2),

[0118] Step 4: Approximate the slope distance expression established in Step 1 using Legendre polynomials to obtain the Legendre expansion of the two-way slope distance. The specific process is as follows:

[0119] In this embodiment, the slope distance R mentioned in step 1 bi (τ) When using Legendre polynomials for fourth-order approximation, arrange the power series of the azimuth slow time.

[0120] 4a) Normalize the slow time τ in the direction of orientation, that is, we have Ta is the time for synthesizing the aperture. Then, the slope distance R is... bi (τ) is decomposed into Legendre orthogonal form and expanded into a fourth power series of x as shown in formula (3):

[0121] R bi (x)=α0+α1x+α2x 2 +α3x 3 +α4x 4 (3),

[0122] Where the corresponding term coefficient α i satisfy:

[0123]

[0124] Where L n (x) is a Legendre polynomial that satisfies the formula

[0125] And in formula (4), the operation process in the braces {} is definite integral operation, and the calculation result can be obtained by bringing in the corresponding order n during operation. The specific operation content meets the formula:

[0126] 4b), the substitute back into formula (5), and rearrange to get the two-way slant range R bi (τ) is the fourth-order power series expansion of slow time τ, that is, the Legendre polynomial approximation result of slant range, as shown in formula (5):

[0127] R bi (τ) = k0+k1τ+k2τ 2 +k3τ 3 +k4τ 4 (5),

[0128] Among them, is the Legendre decomposition coefficient corresponding to the slant range.

[0129] Step 5, combined with the series inversion method, the Legendre expansion of the two-way slant range in step 4 is used to obtain the phase ψ(f r ,f a ) of the original two-dimensional spectrum of the point target. The phase ψ(f r ,f a ) is expanded into a power series form in the range frequency f r , and the high-order terms in the range frequency f r and the azimuth frequency f a are discarded, to obtain the two-dimensional spectrum phase based on the Legendre polynomial expansion. The specific process is as follows:

[0130] 5a), first, the azimuth time-frequency transform of the de-linear phase is completed by using the stationary phase principle, and the result is shown in formula (6):

[0131]

[0132] Then, the corresponding relationship between the azimuth slow time and the frequency is obtained by the series inversion method:

[0133]

[0134] In formulas (6) and (7), k n (n=2, 3, 4) is the Legendre decomposition coefficient corresponding to the slant range, f r is the range frequency, f a is the azimuth frequency, f c is the carrier frequency, and c is the speed of light.

[0135] 5b), the linear phase removed before is supplemented back, and the frequency shift property of the Fourier transform is used to obtain the original two-dimensional spectrum as shown in formula (8):

[0136] S(f r ,f a )=W r (f r )W a (f a )exp(jψ(f r ,f a )) (8),

[0137] where W r represents the range envelope, W a represents the azimuth envelope, and ψ(f r ,f a ) represents the original two-dimensional spectrum phase, and the specific expression is shown in equation (9):

[0138]

[0139] where k n (n = 0, 1, 2, 3, 4) is the corresponding Legendre decomposition coefficient of slant range, f r is the range frequency, f a is the azimuth frequency, f c is the carrier frequency, and c is the speed of light.

[0140] The Legendre polynomial is used to expand in equation (11). First, normalization processing is performed, i.e., let B is the signal bandwidth. The three fractional expressions are respectively arranged into the third power series of y as shown in equation (10):

[0141]

[0142] where the coefficients in front of y of different power terms satisfy the relationship shown in equation (11):

[0143]

[0144] The operation process in the curly braces {} is the definite integral operation, and the calculation result can be obtained by bringing in the corresponding order n during operation, and the specific operation content satisfies the formula: Then, the is brought into the power series of y, and is arranged into the third power series of f r , to obtain equation (12):

[0145]

[0146] where is the corresponding Legendre decomposition coefficient of .

[0147] 5c) Substitute the result of equation (12) into equation (9) and rearrange and merge, discarding f. r and f a After the higher-order terms, the Legendre spectrum phase expression can be obtained as shown in formula (13):

[0148]

[0149] Where the distance frequency f r The coefficients of different power terms satisfy the following relationship:

[0150]

[0151]

[0152] Where, k n (n = 0, 1, 2, 3, 4) are the Legendre decomposition coefficients corresponding to the slant distance, f r f is the distance frequency. a f is the azimuth frequency. c denoted as carrier frequency, and c as the speed of light.

[0153] Observing the phase terms in equation (13), we find that with f r and f a They are all irrelevant and are constant phase terms; It is f r The coefficients of the linear term correspond to distance migration; It is f r The coefficient of the squared term corresponds to distance compression; It is the coefficient of the cubic term; Then only with f a It is related to the azimuth modulation term, which corresponds to azimuth compression.

[0154] Step 6: Combining the phase terms of the two-dimensional spectrum based on Legendre polynomial expansion obtained in Step 5, perform range compression and secondary compression, range migration correction, constant and higher-order term compensation, and azimuth compression on the unambiguous spectrum obtained in Step 3 to obtain the imaging result. The specific process is as follows:

[0155] 6a) First, construct the distance compression and quadratic compression functions as shown in formula (16), where γ is the modulation frequency:

[0156]

[0157] Next, distance migration correction is performed, and the compensation function is constructed as shown in formula (17):

[0158]

[0159] Constant term and high order term compensation is then performed, and the compensation function is shown in equation (18):

[0160]

[0161] 6b) After the above distance compression, migration correction, and constant term and high order term compensation, the compensated two-dimensional frequency domain signal S'(f r ,f′ a ) is subjected to distance IFFT to convert the signal back to the time domain to obtain S'(t,f′ a ), and then the idea of Dechirp processing is combined to construct an azimuth compression function to convert the phase of the IFFT signal, as shown in equation (19)

[0162]

[0163] The result of multiplying S'(t,f′ a ) and H az is subjected to azimuth IFFT to obtain an azimuth time domain unambiguous signal, and then multiplied by a Dechirp function to obtain equation (20):

[0164]

[0165] In equations (19) and (20), ξ satisfies the relationship η satisfies the relationship λ is the wavelength of the carrier, v T and v R are the moving speeds of the transmitting and receiving platforms, R tef and R ref are the distances from the transmitting and receiving platforms to the virtual rotation center point in the two-way sliding spotlight mode, R Tcen and R Rcen are the minimum values of the slant ranges from the transmitting and receiving platforms to the ground point targets.

[0166] 6c) Finally, the azimuth time domain unambiguous signal after the above series of operations is subjected to azimuth FFT to convert it back to the frequency domain, and the unambiguous image is obtained.

[0167] Embodiment Two

[0168] This embodiment provides a simulation experiment of the star-machine bistatic SAR imaging method in the two-way sliding spotlight mode described in embodiment one, as follows:

[0169] (I) Simulation conditions:

[0170] The simulation parameters in Table 1 are used for ground point target array simulation experiments to verify the focusing effect of the star-machine bistatic SAR frequency domain imaging method based on echo preprocessing combined with Legendre polynomials proposed in embodiment one.

[0171] Table 1 Star-Plane Bistatic SAR simulation parameters

[0172]

[0173] Using the parameters in Table 1, the size of the star-plane bistatic SAR ground imaging scene is set to 3km x 3km, and 5 x 5 uniformly distributed point targets are arranged in the imaging scene, as shown in FIG. 1, wherein P1 is a center point of the scene, and P3 is an edge point of the scene. Figure 4

[0174] (II) Simulation content:

[0175] (1) Comparison of slant range approximation errors of two expansion methods

[0176] In this embodiment, the slant range approximation errors of the slant range based on the Legendre polynomial expansion and the slant range based on the traditional Taylor series expansion are compared, and the results are shown in FIG. 2. It can be seen that when the azimuth time is small, the slant range model errors of the two expansion methods are similar; and when the azimuth time increases, the approximation accuracy of the fourth-order Legendre expansion of the slant range model is higher than that of the fourth-order Taylor series, and the error is smaller. Figure 5

[0177] This is because the Taylor series expansion method is to approximate the slant range near the center time of the slant range model at the center time of the synthetic aperture. When the azimuth time is far away from the center time, the slant range error of the target point will increase. If this problem is to be solved, the slant range needs to be expanded to a higher order expression, which further increases the complexity of formula derivation and calculation amount. The Legendre polynomial expansion method is to use the instantaneous slant range at the corresponding azimuth time for weighted operation, which reduces the approximation slant range error of the target point far away from the center time and ensures the stable approximation of the slant range model.

[0178] (2) Effect display of echo dealiasing preprocessing

[0179] According to the echo preprocessing scheme proposed in embodiment one, the two-dimensional frequency domain results of the echo before and after dealiasing can be obtained, as shown in FIG. 3. Wherein Figure 6 (a) is the two-dimensional frequency spectrum before dealiasing, and (b) is the two-dimensional frequency spectrum after dealiasing. It can be seen that after the echo signal is preprocessed by the dealiasing process, the two-dimensional frequency spectrum width is contracted, and aliasing no longer occurs. Figure 6

[0180] (3) Comparison of point target imaging contour maps

[0181] In order to evaluate the focusing effect of the imaging algorithm proposed in embodiment one on the point targets at different positions in the scene, the contour maps of the center point P1 and the edge point P3 in the scene are selected for analysis, and the contour results are shown in FIG. 4. Figure 7 ​​​The energy concentration degrees of the center point P1 and the edge point P3 are higher, thereby verifying that the frequency domain imaging method proposed in the embodiment one has better focusing performance.

[0182] By using the same idea, the two-dimensional frequency domain imaging of the points P1 and P3 is performed by using the traditional frequency domain imaging algorithm based on Taylor series expansion, and the imaging results are as shown in Figure 7 The imaging results are as shown in (c) and (d). Figure 7 By comparing (a) and (c), for the center point P1, the imaging effects of the two imaging algorithms are similar. Figure 7 By comparing (b) and (d), it can be found that, for the edge point, the traditional frequency domain algorithm based on Taylor series expansion has obvious sidelobe lengthening in the azimuth direction of the imaging result, and the focusing effect is poor. Figure 7

[0183] (4) Point target imaging contour effect comparison

[0184] The performance analysis of the imaging point results of the method in the embodiment one and the traditional frequency domain algorithm based on Taylor expansion is further performed, the distance and azimuth profile graphs of the contour of the edge point P3 are obtained, and the results are as shown in (a) to (d). Figure 8 The peak sidelobe ratio (PSLR) and integral sidelobe ratio (ISLR) of the distance and azimuth profile graphs corresponding to the center point P1 and the edge point P3 are calculated and recorded in Table 2.

[0185] In (a) and (b), (c) and (d), respectively, the distance and azimuth profile graphs of the edge point P3 are shown. Figure 8 By comparing (a) and (c), it can be found that, for the edge point P3, the performances of the two imaging algorithms in the distance direction are similar. By comparing (b) and (d), it can be obviously found that, at the edge point P3, the Taylor expansion imaging algorithm has the phenomena of main lobe widening and sidelobe rising in the azimuth profile graph, and the imaging quality is reduced. It can be seen that the profile result of the method proposed in the embodiment one at the edge is closer to the result of the center point.

[0186] Table 2: Imaging quality evaluation of the two methods

[0187]

[0188] ​It can be seen from the data in Table 2 that, in terms of the performance parameters of the range direction, the effect of the method proposed in Embodiment 1 is slightly improved relative to the MSR algorithm of the traditional Taylor expansion, all below-13.2 dB; and in terms of the performance parameters of the azimuth direction, at the edge point P3, the PSLR of the method proposed in Embodiment 1 is improved from-12.773 dB to-13.311 dB, and the ISLR is improved from-8.557 dB to-10.728 dB, compared with the algorithm of the Taylor expansion, which has a greater improvement. In summary, the two-dimensional frequency domain imaging method of the space-airborne bistatic SAR proposed in Embodiment 1 makes the point target generally obtain better performance parameters, and can effectively improve the imaging quality of the scene edge point.

[0189] The preferred embodiments of the present application are described in detail above with reference to the accompanying drawings, and the embodiments described in the present application are merely descriptions of the preferred embodiments of the present application, and do not limit the concept and scope of the present application. In the above specific embodiments, various specific technical features described can be combined in any appropriate manner without contradiction, and such combination should also be considered as disclosed by the present disclosure, as long as it does not deviate from the concept of the present application. In order to avoid unnecessary repetition, the present application will not further describe various possible combinations.

[0190] The present application is not limited to the specific details described in the above embodiments, and various modifications and improvements of the technical solutions of the present application made by those skilled in the art within the scope of the technical concept of the present application and without departing from the design concept of the present application should fall within the protection scope of the present application. The technical content claimed by the present application has been fully recorded in the claims.

Claims

1. A star-aircraft bistatic SAR imaging method in bidirectional sliding spotting mode, characterized in that, Includes the following steps: Step 1: Establish the slant range expression for the bidirectional sliding spotting mode based on the geometric model of the star-aircraft dual-base SAR system; Step 2: Obtain the ground point target echo signal generated when the low-orbit satellite platform transmits a linear frequency modulated signal; Step 3: Perform dealiasing preprocessing on the point target echo to obtain an unambiguous spectrum; Step 4: Approximate the slope distance expression established in Step 1 using Legendre polynomials to obtain the Legendre expansion of the two-way slope distance. Step 5: Using the series inversion method, obtain the phase of the original two-dimensional spectrum of the point target using the Legendre expansion of the two-way slant range described in Step 4. Phase Expanded into the distance frequency domain The power series form, and discarding the distance-frequency domain. and azimuth frequency domain The higher-order terms yield a two-dimensional spectral phase based on Legendre polynomial expansion; Step 6: Combining the phase terms of the two-dimensional spectrum based on Legendre polynomial expansion obtained in Step 5, perform range compression and secondary compression, range migration correction, constant term and higher-order term compensation and azimuth compression on the unambiguous spectrum obtained in Step 3 to obtain the imaging result. In step 5, the two-dimensional spectral phase process based on Legendre polynomial expansion is obtained as follows: 5a) First, the azimuth-direction time-frequency transformation of the delinear phase is completed using the stationary phase principle, and the result is shown in formula (6): (6)， Then, the correspondence between azimuth slow time and frequency is obtained by series inversion: (7)， In formulas (6) and (7), These are the Legendre decomposition coefficients corresponding to the slant distance. For distance frequency, For azimuth frequency, Where c is the carrier frequency and c is the speed of light; 5b) By replenishing the previously removed linear phase and utilizing the frequency shift property of the Fourier transform, the original two-dimensional spectrum can be obtained as shown in formula (8): (8) in, Represents distance to the envelope, Represents the direction of the envelope. The phase of the original two-dimensional spectrum is represented by the formula (9): (9) in, These are the Legendre decomposition coefficients corresponding to the slant distance. For distance frequency, For azimuth frequency, Where c is the carrier frequency and c is the speed of light; Using Legendre polynomials to modify equation (9) To expand, first perform normalization, that is, let B is the signal bandwidth; the three fractional expressions are rearranged into the third power series form of y as shown in formula (10): (10)， The coefficients of y for different power terms satisfy the relationship shown in formula (11): (11)， The operations within the curly braces {} are definite integral operations; the corresponding order should be substituted during the calculation. The calculation result is sufficient; the specific operations must satisfy the formula: Then Substituting into the power series of y, we can simplify to: The third power series form is given by formula (12): (12)， in, yes The corresponding Legendre decomposition coefficients; 5c) Substitute the result of equation (12) into equation (9), rearrange and merge, and discard. and After the higher-order terms, the Legendre spectrum phase expression can be obtained as shown in formula (13): (13)， Among them, distance frequency The coefficients of different power terms satisfy the following relationship: (14)， (15)， in, These are the Legendre decomposition coefficients corresponding to the slant distance. For distance frequency, For azimuth frequency, Where c is the carrier frequency and c is the speed of light; Observing the phase terms in equation (13), we find that and and They are all irrelevant and are constant phase terms; yes The coefficients of the linear term correspond to distance migration; yes The coefficient of the squared term corresponds to distance compression; It is the coefficient of the cubic term; Then only with It is related to the azimuth modulation term, which corresponds to azimuth compression.

2. The method for star-aircraft bi-static SAR imaging in bi-directional sliding spotting mode according to claim 1, characterized in that, In step 3, the specific process of dealiasing preprocessing for the point target echo is as follows: 3a) Constructing a reference function: First, construct the orientation time domain. The reference function is shown in formula (1): (1)， in, The frequency tuning coefficient is related to the actual configuration parameters of the satellite-aircraft bistatic SAR system and satisfies the following relationship. , The wavelength of the carrier wave. and These represent the movement speed of the sending and receiving platforms, and These represent the distances from the transceiver platform to the virtual rotation center point in the bidirectional sliding beam-gathering mode. The time-domain echo signal obtained in step 2 Perform a distance FFT to transform it into the distance frequency domain. and will With reference function Multiply; 3b) Equivalent convolution transformation: Perform an azimuth FFT on the result of step 3a) and compare it with the new time-domain coordinates. Reference functions below Multiply to obtain new time-domain coordinates Unambiguous signals in the lower bit frequency domain This step can also be viewed as a process for analyzing time-domain echoes. With reference function Perform an equivalent convolution transformation to obtain ,in: The FFT transform is equivalent to a transformation between the time and frequency axes. Before and after the FFT transform, the original time-domain coordinates of the echo signal change. New time-domain coordinates The original frequency coordinates of the echo signal New frequency coordinates The relationship between them is satisfied ; 3c) Conjugate compensation: For the result obtained in step 3b) Perform azimuth FFT again to obtain a new two-dimensional frequency domain signal. ,right Multiply by the frequency domain conjugate compensation function The resulting non-aliasing spectrum is shown in equation (2): (2)。 3. The method for star-aircraft bi-static SAR imaging in bidirectional sliding spotting mode according to claim 2, characterized in that, In step 4, the process of approximating the slope distance using Legendre polynomials is as follows: 4a) Opposite direction slow time After normalization, we have , For the synthesis aperture time; then for the slant distance Perform Legendre orthogonal decomposition and expand it into The fourth-order power series form is shown below: (3)， Where the corresponding term coefficient satisfy: (4)， in For Legendre polynomials, satisfying the formula ; In equation (4), the operation within the curly braces {} is a definite integral operation, and the corresponding order is substituted during the operation. The calculation result is sufficient; the specific operations must satisfy the formula: ; 4b) will Substituting back into equation (3), we can rearrange the equation to obtain the double-run slant distance. Slow time The fourth-order power series expansion, that is, the Legendre polynomial approximation of the slant distance, is: (5)， in, is the Legendre decomposition coefficient corresponding to the slant distance.

4. The method for star-aircraft bi-static SAR imaging in bi-directional sliding spotting mode according to claim 1, characterized in that, In step 6, the non-aliased two-dimensional spectrum is processed. The following operations are performed sequentially: range compression and secondary compression, migration correction, constant and higher-order term compensation, and azimuth compression. The specific process is as follows: 6a) First, construct the distance compression and quadratic compression functions as shown in formula (16), where To adjust the frequency: (16)， Next, distance migration correction is performed, and the compensation function is constructed as shown in formula (17): (17)， Then, constant terms and higher-order terms are compensated, and the compensation function is shown in formula (18): (18)； 6b) After the above-mentioned distance compression, migration correction, and compensation for constant and higher-order terms, the compensated two-dimensional frequency domain signal... Perform a distance IFFT to convert the signal back to the time domain. Furthermore, combining the ideas of Decirp processing, an azimuth compensation function is constructed to transform the phase of the signal after IFFT, as shown in formula (19). （19）， right and The result of multiplication is used to perform azimuth IFFT to obtain an unambiguous azimuth time domain signal, which is then multiplied by the Decirp function to obtain formula (20): (20)， In equations (19) and (20), Satisfying the relation , Satisfying the relation , The wavelength of the carrier wave. and These represent the movement speed of the sending and receiving platforms, and These represent the distances from the transceiver platform to the virtual rotation center point in bidirectional sliding beamforming mode. and These are the minimum slant distances from the transceiver platform to the ground target point; 6c) Finally, the azimuth time-domain unambiguous signal after the above series of operations is transformed back to the frequency domain by azimuth FFT to obtain an unambiguous image.

Images