Bistatic inverse synthetic aperture radar imaging and calibration method based on motion parameter estimation

Abstract

The invention discloses a bistatic inverse synthetic aperture radar imaging and calibration method based on motion parameter estimation, and belongs to the technical field of radar signal processing and radar imaging. The method solves the problems that an existing method is poor in anti-noise performance, related scattering points are difficult to find, and the azimuth defocusing effect is not eliminated. According to the method, an equation set is constructed by using a bistatic ISAR echo phase coefficient, the ratio of a first-order coefficient to a constant coefficient of a bistatic coefficient and an estimated value of an effective rotation angular velocity are obtained by solving the equation set, and then the constant coefficient of the bistatic coefficient is estimated by using a receiver beam center direction and target distance information; and constructing a compensation phase based on the estimated value to effectively compensate image distortion and range defocus, thereby realizing calibration of the inverse synthetic aperture radar image. The method can be applied to bistatic inverse synthetic aperture radar imaging and calibration.

Description

Technical Field

[0001] The invention belongs to the technical field of radar signal processing and radar imaging, and in particular relates to a bistatic inverse synthetic aperture radar (ISAR) imaging and calibration method based on motion parameter estimation. Background Art

[0002] Due to its characteristics of all-day, all-weather, long-distance monitoring and two-dimensional high-resolution imaging, ISAR has been widely used in target recognition and classification in the past few decades. Among them, most research work focuses on monostatic ISAR. However, when the target moves along the radar line of sight (LOS) direction, the monostatic ISAR cannot obtain a two-dimensional image of the target, that is, the blind angle phenomenon of the monostatic ISAR. Bistatic ISAR effectively overcomes the blind angle phenomenon of the monostatic ISAR by deploying the transmitter and the receiver separately. It also has the advantages of concealment, anti-interference, long detection distance, and rich target scattering characteristics. Therefore, in recent years, bistatic ISAR has received more and more attention.

[0003] The main task of bistatic inverse synthetic aperture radar is to obtain high-quality and high-precision calibration images. The main process includes two stages: imaging and calibration. During the imaging process, image distortion and defocusing effects may occur due to the influence of the time-varying bistatic coefficient. Once the bistatic factor is estimated, these negative effects can be eliminated. The purpose of the calibration process is to extract the size information of the target and improve the accuracy of target recognition and classification. In bistatic inverse synthetic aperture radar, the range and lateral scaling factors depend not only on the bandwidth of the transmitted signal and the effective rotation speed of the target, but also on the bistatic coefficient. Therefore, the core of bistatic inverse synthetic aperture radar imaging and calibration lies in estimating the bistatic coefficient and the effective rotation speed.

[0004] Kang MS et al. proposed a bistatic inverse synthetic aperture radar imaging and calibration algorithm based on tracking information in the document "Bistatic-ISAR Cross-Range Scaling. (IEEE Transactions on Aerospace and Electronic Systems, vol. 53, no. 4, pp. 1962-1973, Aug. 2017)". This algorithm estimates the bistatic coefficients by using the moving target trajectory obtained by the tracking device, and estimates the effective rotational angular velocity by calculating the phase coefficient of the selected prominent point. However, this method cannot be implemented without a tracking device.

[0005] Later, Kang MS et al. proposed a bistatic inverse synthetic aperture radar distortion correction and calibration algorithm without tracking equipment in the document "Bistatic-ISAR Distortion Correction and Rangeand Cross-Range Scaling. (IEEE Sensors Journal, vol. 17, no. 16, pp. 5068-5078, 15Aug. 15, 2017)". This algorithm predicts the phase coefficient of the extracted scatterer by implementing phase unpacking and least squares fitting, and then applies the particle swarm optimization algorithm to maximize the established cost function to estimate the effective rotation angular velocity and bistatic coefficient. However, the anti-noise performance of this algorithm is weak.

[0006] Sun Sibo et al. proposed an algorithm for eliminating image distortion based on scattering point association in the document "Bistatic ISAR distortion and defocusing analysis. (IEEE Transactions on Aerospace and Electronic Systems, vol. 52, no. 3, pp. 1168-1182, June 2016)". This method requires finding associated scattering points in bistatic and monostatic ISAR images. However, since the image projection planes of bistatic ISAR and monostatic ISAR are usually different, it makes it difficult to find associated scattering points in the two images. In addition, all the above methods do not consider eliminating the azimuth defocus effect.

[0007] In summary, the existing methods for estimating the dual-base coefficients and effective rotational angular velocity without the need for tracking equipment have the problems of poor anti-noise performance, difficulty in finding associated scattering points, and failure to eliminate the azimuth defocus effect. Summary of the invention

[0008] The purpose of the present invention is to solve the problems of poor anti-noise performance, difficulty in finding associated scattering points and failure to eliminate azimuth defocus effect in existing methods, and to propose a bistatic inverse synthetic aperture radar imaging and calibration method based on motion parameter estimation.

[0009] The technical solution adopted by the present invention to solve the above technical problems is: a bistatic inverse synthetic aperture radar imaging and calibration method based on motion parameter estimation, the method specifically comprising the following steps:

[0010] Step 1: The bistatic inverse synthetic aperture radar obtains echo data reflected by scattering points on the uniformly rotating target, and then performs distance compression and motion compensation processing on the echo data to obtain echo data reflected by each scattering point on the uniformly rotating target after distance compression and motion compensation processing;

[0011] Step 2: Based on the echo data after range compression and motion compensation processing of each scattering point on the uniformly rotating target, obtain the estimated value of the ratio of the first-order coefficient to the constant coefficient in the dual-base coefficient and an estimate of the square of the effective angular velocity

[0012] Step 3: Calculate the bistatic angle at the initial moment, and then estimate the constant coefficient K0 in the bistatic coefficient based on the bistatic angle at the initial moment;

[0013] Step 4: Based on the estimated value and Construct the phase compensation term s of the i-th distance unit com,i , and then multiplying the phase compensation term of the i-th distance unit with the echo data after range compression and motion compensation corresponding to the i-th distance unit to obtain the echo data of the i-th distance unit after phase compensation;

[0014] Step 5, compressing the echo data of all range units after phase compensation in step 4 in azimuth direction to obtain echo data after azimuth compression;

[0015] Step 6: Obtain a scattering point echo map based on the echo data compressed in azimuth;

[0016] Step 7: Obtain the target azimuth dimension y according to the scattering point echo map obtained in step 6 max , and then according to the constant coefficient K0 and the estimated value Calculate K1 and determine the target orientation size y max Is it satisfied? Among them, κ w is the widening factor of the window function;

[0017] If it is not satisfied, the scattering point echo map is post-compensated and the post-compensated image is used as the inverse synthetic aperture radar calibration image;

[0018] If it is satisfied, the scattering point echo image obtained in step 6 is the inverse synthetic aperture radar calibration image.

[0019] Furthermore, the echo data after distance compression and motion compensation processing reflected by each scattering point on the uniformly rotating target is specifically:

[0020]

[0021] Where σ is the bistatic scattering coefficient of the scattering point, B is the signal bandwidth, and T obs is the observation duration, rect(·) is the rectangular window function, sinc(·) is the sinc function, j is the imaginary unit, λ is the carrier wavelength, t is the slow time, τ is the fast time delay variable, c is the speed of light, and f c is the center frequency, ξ, α, β and γ are the phase coefficients;

[0022]

[0023] Among them, K0 and K1 are the constant coefficient and the first-order coefficient in the dual base coefficient, respectively, (x, y) is the coordinate of the scattering point relative to the rotation center, Ω eff is the effective angular velocity.

[0024] Furthermore, the specific process of step 2 is as follows:

[0025] Step 21: selecting the distance unit based on the minimum normalized amplitude variance criterion and the processed echo data reflected from each scattering point;

[0026] Step 2: Determine the scattering point with the largest echo energy in each selected distance unit, and then use the polynomial phase estimation method to estimate the phase coefficient corresponding to each determined scattering point, and record the phase coefficient corresponding to the scattering point with the largest echo energy in the kth distance unit as ξ k ,α k ,β k ,γ k , k=1,2,…,K;

[0027] Step 23: Use the phase coefficients estimated in step 22 to construct the equation system:

[0028]

[0029] Among them, x k '=K0x k ,y k '=K0y k ,(x k ,y k ) is the coordinate value of the scattering point with the largest echo energy in the kth range unit;

[0030] Step 24: Convert the equation system of formula (3) into the following form:

[0031]

[0032] make but:

[0033]

[0034] Step 25: Substitute any set of phase coefficients obtained in step 22 into equation (5), and then solve equation (5) to obtain the three roots corresponding to the set of phase coefficients; similarly, solve the three roots corresponding to each set of phase coefficients obtained in step 22 respectively;

[0035] Then, according to all the positive roots obtained, we can get Estimated value of

[0036] Step 26: Estimated value Substitute into the equation group of formula (3), calculate a ratio according to each set of phase coefficients, and then calculate the mean of each ratio, and use the calculated mean as the estimation result

[0037] Furthermore, the specific process of step 21 is as follows:

[0038] Normalized amplitude variance NAV of the i-th range unit i for:

[0039]

[0040] Among them, s i (n) is the processed echo data corresponding to the nth scattering point in the ith range unit, E n (·) is the mean operation, |s i (n)| means s i (n) the magnitude;

[0041] Similarly, after calculating the normalized amplitude variance of each distance unit, the normalized amplitude variances of each distance unit are sorted in order from small to large, and the distance units corresponding to the normalized amplitude variances in the first K positions are selected.

[0042] Further, in the step 25, according to all the positive roots obtained, Estimated value of Specifically:

[0043] Cluster all positive roots, select the cluster with the largest density from the clustering results, calculate the mean of all positive roots in the selected cluster, and use the calculated mean as the estimated value

[0044] Furthermore, the specific process of step three is:

[0045] Let T, R and O represent the positions of radar transmitter Tx, receiver Rx and target rotation center respectively, and let the line connecting T and R be represents the baseline of the bistatic inverse synthetic aperture radar, and the straight-line distance between T and R is L, and the straight-line distance between T and O is R T , let the line connecting R and O be The straight-line distance between R and O is R R ;

[0046] Will and The angle between Indicates the angle between the LOS direction of the receiving radar and the baseline direction, extended To point A, The length is equal to Then connect point A to T and make a line from R to The perpendicular line of , and the foot of the perpendicular is C;

[0047] Calculate θ0 using the law of sines for a triangle:

[0048]

[0049] Among them, θ0 is the bistatic angle at the initial moment, R RT =R R +R T , that is, R RT Represents R T and R R The distance and;

[0050] Then the constant coefficient K0 of the double base coefficient is:

[0051] K0=cos(θ0 / 2) (10).

[0052] Furthermore, the specific process of step 4 is as follows:

[0053]

[0054] Among them, ξ i =τ i c / 2, i represents the distance unit index value, τ i is the fast time delay corresponding to the i-th range unit.

[0055] Furthermore, the specific process of step five is as follows:

[0056]

[0057] Among them, s rc (τ,f d ) is the echo data after azimuth compression, f d is the Doppler variable.

[0058] Furthermore, the specific process of step six is ​​as follows:

[0059]

[0060] Among them, ΔX and ΔY are the range resolution and azimuth resolution of the bistatic inverse synthetic aperture radar image, respectively, and X and Y are the coordinates of the range dimension and azimuth dimension, respectively;

[0061] Substituting equation (13) into equation (12) and taking the absolute value of the right side of equation (12) we can obtain the scattering point echo map of the bistatic inverse synthetic aperture radar:

[0062]

[0063] Among them, I(X,Y) is the scattering point echo pattern of the bistatic inverse synthetic aperture radar.

[0064] Furthermore, the specific process of the post-compensation processing is as follows:

[0065] Step 7: Segment the scattering point echo map obtained in step 6, that is, define the segmentation step length. Then, the scattering point echo map is divided into various regions along the azimuth according to the segmentation step size, and the azimuth intervals corresponding to each region are recorded as [-Ry p , -(R-1)y p ],(-(R-1)y p , -(R-2)y p ], ..., (-y p ,0],(0,y p ],…,((R-2)y p ,(R-1)y p ],((R-1)y p ,Ry p ], 2R is the total number of azimuth intervals;

[0066] The azimuth interval is (ry p ,(r+1)y p ] is taken as the rth positive defocus area, and the azimuth interval is (-(r-1)y p , -(r-2)y p ] is taken as the rth negative defocus area, where r is a positive integer; the azimuth interval is (-y p ,0] and (0,y p ] is the non-defocused area;

[0067] Step 72: extract the image of each defocused area, perform inverse Fourier transform on the image of each defocused area in the azimuth dimension, and convert the image from the range-azimuth domain to the range-slow time domain;

[0068] The rth positive defocused area image in the range-slow time domain is combined with the post-compensation phase signal Multiply the rth negative defocus area image in the distance-slow time domain with the post-compensation phase signal Multiply, perform Fourier transform on each multiplied image along the slow time dimension, transform each defocused area image from the range-slow time domain to the range-azimuth domain, synthesize the transformed image of each defocused area with the image of the non-defocused area, and obtain the focused and calibrated inverse synthetic aperture radar image.

[0069] The beneficial effects of the present invention are:

[0070] The present invention does not need to rely on external tracking and angle measurement equipment. It can accurately estimate the time-varying bistatic coefficient and the effective rotation angular velocity only through the phase coefficient of the echo signal, the center direction of the receiver beam and the measured distance information, and eliminate the image distortion and defocusing effect based on the estimation result, so as to realize the calibration of the inverse synthetic aperture radar image without the operation of finding the associated scattering points. Moreover, the method of the present invention can also perform well under low signal-to-noise ratio conditions and has high robustness. BRIEF DESCRIPTION OF THE DRAWINGS

[0071] Figure 1 It is a flow chart of a bistatic inverse synthetic aperture radar imaging and calibration method based on motion parameter estimation of the present invention;

[0072] Figure 2 This is a geometric model diagram of bistatic inverse synthetic aperture radar and space targets;

[0073] Figure 3 Schematic diagram of the target view estimation process;

[0074] Figure 4 is the scattering point model of the simulation target;

[0075] Figure 5(a) is the imaging result of the range Doppler algorithm at a -5dB signal-to-noise ratio;

[0076] FIG5( b ) is an imaging result diagram of the range Doppler algorithm under a 0 dB signal-to-noise ratio;

[0077] Figure 5(c) is the imaging result of the range Doppler algorithm at a 5dB signal-to-noise ratio;

[0078] Figure 5(d) is the imaging result of the range Doppler algorithm at a signal-to-noise ratio of 10 dB;

[0079] FIG6( a ) is a diagram showing the imaging result of the method of the present invention at a signal-to-noise ratio of -5 dB;

[0080] FIG6( b ) is a diagram showing the imaging result of the method of the present invention at a signal-to-noise ratio of 0 dB;

[0081] FIG6( c ) is a diagram showing the imaging result of the method of the present invention at a signal-to-noise ratio of 5 dB;

[0082] FIG6( d ) is a diagram showing the imaging result of the method of the present invention at a signal-to-noise ratio of 10 dB;

[0083] Figure 7(a) shows the scattering point model of Boeing 737;

[0084] FIG7( b ) is an image of the Boeing 737 model without post-compensation processing;

[0085] FIG7( c ) is an image showing the result of post-compensation processing of the Boeing 737 model. DETAILED DESCRIPTION

[0086] Specific implementation method 1: Combination Figure 1 This embodiment describes a bistatic inverse synthetic aperture radar imaging and calibration method based on motion parameter estimation, wherein the motion parameter estimation includes estimation of bistatic coefficients and effective rotation angular velocity, and the method specifically includes the following steps:

[0087] Step 1: The bistatic inverse synthetic aperture radar obtains echo data reflected by scattering points on the uniformly rotating target, and then performs distance compression and motion compensation processing on the echo data to obtain echo data reflected by each scattering point on the uniformly rotating target after distance compression and motion compensation processing;

[0088] Step 2: Based on the echo data after range compression and motion compensation processing of each scattering point on the uniformly rotating target, obtain the estimated value of the ratio of the first-order coefficient to the constant coefficient in the dual-base coefficient and an estimate of the square of the effective angular velocity

[0089] Step 3: Calculate the bistatic angle at the initial moment, and then estimate the constant coefficient K0 in the bistatic coefficient based on the bistatic angle at the initial moment;

[0090] Step 4: Based on the estimated value and Construct the phase compensation term s of the i-th distance unit com,i , and then multiplying the phase compensation term of the i-th distance unit with the echo data after range compression and motion compensation corresponding to the i-th distance unit to obtain the echo data of the i-th distance unit after phase compensation;

[0091] Step 5, compressing the echo data of all range units after phase compensation in step 4 in azimuth direction to obtain echo data after azimuth compression;

[0092] Step 6: Obtain a scattering point echo map based on the echo data compressed in azimuth;

[0093] Step 7: Obtain the target azimuth dimension y according to the scattering point echo map obtained in step 6 max , and then according to the constant coefficient K0 and the estimated value Calculate K1 and determine the target orientation size y max Is it satisfied? Among them, κ w is the widening factor of the window function;

[0094] If it is not satisfied, the scattering point echo map is post-compensated and the post-compensated image is used as the inverse synthetic aperture radar calibration image;

[0095] If it is satisfied, the scattering point echo image obtained in step 6 is the inverse synthetic aperture radar calibration image.

[0096] The present invention uses the bistatic ISAR echo phase coefficient to construct an equation group, obtains the ratio of the first-order coefficient to the constant coefficient of the bistatic coefficient and the estimated value of the effective rotation angular velocity by solving the equation group, and then estimates the constant coefficient of the bistatic coefficient using the receiver beam center direction and target distance information. On this basis, a compensation phase is constructed to effectively compensate for image distortion and distance defocus, and calibration processing is performed on the image. For targets of larger size, azimuth defocus effect may occur. Therefore, after image calibration, it is necessary to determine whether the target needs post-compensation processing. If necessary, the image is post-compensated to obtain a focused calibrated bistatic ISAR image.

[0097] Specific implementation method 2: This implementation method is different from the specific implementation method 1 in that the echo data after distance compression and motion compensation processing of each scattering point reflected on the uniformly rotating target is specifically:

[0098]

[0099] Where σ is the bistatic scattering coefficient of the scattering point, B is the signal bandwidth, and T obs is the observation duration, rect(·) is the rectangular window function, sinc(·) is the sinc function, j is the imaginary unit, λ is the carrier wavelength, t is the slow time, τ is the fast time delay variable, c is the speed of light, and f c is the center frequency, ξ, α, β and γ are the phase coefficients;

[0100]

[0101] Among them, K0 and K1 are the constant coefficient and the first-order coefficient in the dual base coefficient, respectively, (x, y) is the coordinate of the scattering point relative to the rotation center, Ω eff is the effective angular velocity.

[0102] The other steps and parameters are the same as those in the first embodiment.

[0103] Specific implementation method three: This implementation method is different from specific implementation methods one or two in that the specific process of step two is as follows:

[0104] Step 21: selecting the distance unit based on the minimum normalized amplitude variance criterion and the processed echo data reflected from each scattering point;

[0105] Step 2: Determine the scattering point with the largest echo energy in each selected distance unit, and then use the polynomial phase estimation method to estimate the phase coefficient corresponding to each determined scattering point, and record the phase coefficient corresponding to the scattering point with the largest echo energy in the kth distance unit as ξ k ,α k ,β k ,γ k , k=1,2,…,K;

[0106] Step 23: Use the phase coefficients estimated in step 22 to construct the equation system:

[0107]

[0108] Among them, x k '=K0x k ,y k '=K0y k ,(x k ,y k ) is the coordinate value of the scattering point with the largest echo energy in the kth range unit;

[0109] Step 24: Convert the equation group of formula (3) into the following form:

[0110]

[0111] make but:

[0112]

[0113] Step 25: Substitute any set of phase coefficients obtained in step 22 into equation (5), and then solve equation (5) to obtain the three roots corresponding to the set of phase coefficients; similarly, solve the three roots corresponding to each set of phase coefficients obtained in step 22 respectively;

[0114] Then, according to all the positive roots obtained, we can get Estimated value of

[0115] Step 26: Estimated value Substitute into the equation group of formula (3), calculate a ratio according to each set of phase coefficients, and then calculate the mean of each ratio, and use the calculated mean as the estimation result

[0116] The other steps and parameters are the same as those in the first or second embodiment.

[0117] Specific implementation method 4: This implementation method is different from specific implementation methods 1 to 3 in that the specific process of step 21 is as follows:

[0118] Normalized amplitude variance NAV of the i-th range unit i for:

[0119]

[0120] Among them, s i (n) is the processed echo data corresponding to the nth scattering point in the ith range unit, E n (·) is the mean operation, |s i (n)| means s i (n) magnitude;

[0121] It should be noted that distance unit i corresponds to the fast time dimension, scattering point n corresponds to the slow time dimension, and one distance unit corresponds to multiple scattering points;

[0122] Similarly, after calculating the normalized amplitude variance of each distance unit, the normalized amplitude variances of each distance unit are sorted in order from small to large, and the distance units corresponding to the normalized amplitude variances in the first K positions are selected.

[0123] The other steps and parameters are the same as those in Specific Embodiments 1 to 3.

[0124] In this implementation manner, the value of K is 5-8.

[0125] Specific implementation method 5: This implementation method is different from the specific implementation methods 1 to 4 in that in step 25, all positive roots are obtained. Estimated value of Specifically:

[0126] Cluster all positive roots, select the cluster with the largest density from the clustering results, calculate the mean of all positive roots in the selected cluster, and use the calculated mean as the estimated value

[0127] The other steps and parameters are the same as those in Specific Embodiments 1 to 4.

[0128] The cubic equation of formula (5) has three roots. Obviously, for this cubic equation, there exists an equal The positive root of . Since the cubic coefficient of equation (5) and the constant term 4γ2 are all positive, and according to the graphical characteristics of a cubic function, we know that the equation of formula (5) must have a negative root (invalid solution). As for the third root, if it is negative, then the only positive root among the three roots is the true solution ( ), but in fact, the true solutions obtained according to each set of phase coefficients are not necessarily exactly the same. If the third root is also a positive root, it is difficult to determine which positive root is the true solution. In this case, the three roots of the equation (5) are defined as Λ1, Λ2, and Λ3, where Λ2 represents a negative root, and Λ3 represents another positive root. The relationship between the roots and coefficients of a cubic equation is:

[0129]

[0130] Will Substituting the equation of formula (3) into the equation of formula (7), the expression of the third root Λ3 is transformed into:

[0131]

[0132] It is difficult to distinguish the magnitudes of Λ1 and Λ3 from equation (8). However, the value of Λ3 is spatially variable because it depends on the distance of the scatterer and the crossover distance position. In contrast, Λ1 is a constant (i.e., spatially invariant) because it is equal to Based on this, the present invention can use each set of phase coefficients to solve the equation to obtain the three roots of the equation. Then use the clustering algorithm to find the positive roots that are clustered together from these roots, and use the mean of the clustered positive roots as Ω eff 2 The estimated results.

[0133] It should be noted that due to Ω eff The sign of is unknown, and the present invention can only estimate Ω eff The absolute value of the target is not used to estimate its true value. This may only lead to the mirror symmetry of the final calibrated ISAR image relative to the cross-range axis, without affecting the estimation results of the target size and shape. Therefore, the present invention does not strictly distinguish between absolute values ​​and actual values.

[0134] Specific implementation method six: Combination Figure 2 This embodiment is different from the first to fifth embodiments in that the specific process of step three is as follows:

[0135] Let T, R and O represent the positions of radar transmitter Tx, receiver Rx and target rotation center respectively, and let the line connecting T and R be represents the baseline of the bistatic inverse synthetic aperture radar, and the straight-line distance between T and R is L, and the straight-line distance between T and O is R T , let the line connecting R and O be The straight-line distance between R and O is R R ;

[0136] Will and The angle between Indicates the angle between the LOS direction of the receiving radar and the baseline direction (i.e., the target viewing angle), extended To point A, The length is equal to Then connect point A to T and make a line from R to The perpendicular line of , and the foot of the perpendicular is C;

[0137] Calculate θ0 using the law of sines for a triangle:

[0138]

[0139] Among them, θ0 is the bistatic angle at the initial moment, R RT =R R +R T , that is, R RT Represents R T and R R The distance and;

[0140] Then the constant coefficient K0 of the double base coefficient is:

[0141] K0=cos(θ0 / 2) (10)

[0142] The other steps and parameters are the same as those in Specific Implementation Methods 1 to 5.

[0143] From the one-dimensional distance image at the initial moment, we can get R RT -L is an estimate with an estimated error less than half of the radar range resolution. The positions of the receiver Rx and the transmitter Tx can be determined by the Global Positioning System or the BeiDou Navigation Satellite System. Then, the distance between the positions of the receiver Rx and the transmitter Tx can be calculated to estimate the length of the baseline L. Since the precise point positioning (PPP) accuracy of GPS or BDS can reach the centimeter level, the estimated error of L can be ignored. Then, by replacing R RT -Calculate R by adding the estimated values ​​of L and L RT The estimated value of R RT -Estimation error of L and L, R RT The estimated error will be in the meter level.

[0144] When the signal received by the receiver Rx exceeds the detection threshold, it indicates that the target is currently in the main beam of the receiver antenna, such as Figure 3According to the directional information of the receiver's beam center (provided by the antenna control system), the beam center angle, that is, the angle between the receiver's beam center and the baseline, can be determined and used as an estimate of the target viewing angle.

[0145] Finally, using L, R RT The estimated value of θ0 can be solved by solving the equation.

[0146] Specific implementation method 7: This implementation method is different from the specific implementation methods 1 to 6 in that the specific process of step 4 is as follows:

[0147]

[0148] Among them, ξ i =τ i c / 2, i represents the distance unit index value, τ i is the fast time delay corresponding to the i-th range unit.

[0149] The other steps and parameters are the same as those in Specific Embodiments 1 to 6.

[0150] Specific implementation eight: This implementation differs from specific implementations one to seven in that the specific process of step five is as follows:

[0151]

[0152] Among them, s rc (τ,f d ) is the echo data after azimuth compression, f d is the Doppler variable.

[0153] The other steps and parameters are the same as those in Specific Embodiments 1 to 7.

[0154] Specific implementation method 9: This implementation method is different from specific implementation methods 1 to 8 in that the specific process of step 6 is as follows:

[0155]

[0156] Among them, ΔX and ΔY are the range resolution and azimuth resolution of the bistatic inverse synthetic aperture radar image, respectively, and X and Y are the coordinates of the range dimension and azimuth dimension, respectively;

[0157] Substituting equation (13) into equation (12) and taking the absolute value of the right side of equation (12) we can obtain the scattering point echo map of the bistatic inverse synthetic aperture radar:

[0158]

[0159] Among them, I(X,Y) is the scattering point echo pattern of the bistatic inverse synthetic aperture radar.

[0160] The other steps and parameters are the same as those in Specific Embodiments 1 to 8.

[0161] Specific implementation method 10: This implementation method is different from the specific implementation methods 1 to 9 in that the specific process of the post-compensation processing is as follows:

[0162] Step 7: Segment the scattering point echo map obtained in step 6, that is, define the segmentation step length. Then, the scattering point echo map is divided into various regions along the azimuth according to the segmentation step size, and the azimuth intervals corresponding to each region are recorded as [-Ry p , -(R-1)y p ],(-(R-1)y p , -(R-2)y p ], ..., (-y p ,0],(0,y p ],…,((R-2)y p ,(R-1)y p ],((R-1)y p ,Ry p ], 2R is the total number of azimuth intervals;

[0163] The azimuth interval is (ry p ,(r+1)y p ] is taken as the rth positive defocus area, and the azimuth interval is (-(r-1)y p , -(r-2)y p ] is taken as the rth negative defocus area, where r is a positive integer; the azimuth interval is (-y p ,0] and (0,y p ] is the non-defocused area;

[0164] Step 72: extract the image of each defocused area, perform inverse Fourier transform on the image of each defocused area in the azimuth dimension, and convert the image from the range-azimuth domain to the range-slow time domain;

[0165] The rth positive defocused area image in the range-slow time domain is combined with the post-compensation phase signal Multiply the rth negative defocus area image in the distance-slow time domain with the post-compensation phase signal Multiply, perform Fourier transform on each multiplied image along the slow time dimension, transform each defocused area image from the range-slow time domain to the range-azimuth domain, synthesize the transformed image of each defocused area with the image of the non-defocused area, and obtain the focused and calibrated inverse synthetic aperture radar image.

[0166] The other steps and parameters are the same as those in Specific Embodiments 1 to 9.

[0167] The simulation parameters are as follows:

[0168] The target scattering point model used in the simulation is as follows: Figure 4 As shown in the figure, in the simulation scenario, the transmitter is located at (5,0,0)Km, the receiver is located at (-5,0,0)Km, the target is located at (0,8.6,0)Km, and the bistatic angle is 60°. The speed of the target is (0,-300,347)m / s. The bistatic angular velocity of the target is 0.03rad / s, and the effective rotation angular velocity is 0.04rad / s. The transmission parameters of the BISAR system are shown in Table 1. The beam width of the receiver is 10°, and the maximum angular error between the beam center and the target is half of the beam width. The echo signal-to-noise ratio (after pulse compression) is set to -5dB, 0dB, 5dB, and 10dB respectively.

[0169] Table 1 System simulation parameters

[0170]

[0171] Figure 5(a)-Figure 5(d) The results of imaging the motion compensated data directly using the range Doppler algorithm under different signal-to-noise ratios show that the image has large distortion and the part close to the wing is defocused. Table 2 is the comparison result of the estimated value and the true value of the dual-base coefficient and the effective rotation angular velocity of the target by the present invention. It can be seen that the estimated value of the method of the present invention is very close to the true value.

[0172] Table 2 Estimation results of bistatic coefficients and effective rotation angular velocity

[0173]

[0174] Figure 6(a)-Figure 6(d) In order to obtain imaging results at different signal-to-noise ratios using the method of the present invention, it can be seen that the method of the present invention effectively compensates for the distortion and defocusing effect of the image at different signal-to-noise ratios and realizes the calibration of the image.

[0175] Figure 4 Six scattering points (ABCDEF) were selected, and the positions of these scattering points under different signal-to-noise ratios and the actual positions of these scattering points are listed in Table 3. By comparing the results, it can be found that the positions of the scattering points on the map are very close to the actual positions, which shows that the method of the present invention has high calibration accuracy under different signal-to-noise ratios.

[0176] Table 3 Comparison results of the coordinates of the six scattering points on the bistatic ISAR image and the true coordinates at different signal-to-noise ratios

[0177]

[0178] The post-compensation processing effect of the present invention is further described below with reference to FIG. 7 .

[0179] The target scattering point model is shown in Figure 7(a). In the simulation scenario, the transmitter is located at (5,0,0)Km, the receiver is located at (-5,0,0)Km, the target is located at (0,8.6,0)Km, and the bistatic angle is 60°. The speed of the target is (0,-300,347)m / s. The bistatic angular velocity of the target is 0.03rad / s, and the effective rotation angular velocity is 0.04rad / s. The transmission parameters of the BISAR system are shown in Table 1. The beam width of the receiver is 10°, the maximum angular error between the beam center and the target is half of the beam width, and the echo signal-to-noise ratio (after pulse compression) is 5dB.

[0180] FIG7( b ) shows the imaging result obtained without post-compensation processing. From the local magnified image of the nose, an obvious defocusing effect can be seen.

[0181] Figure 7(c) shows the imaging result after post-compensation processing. From the local enlarged view of the nose, it can be seen that the defocus effect has been eliminated and the scattering points have achieved a better focusing effect. The comparison of the results shows that post-compensation processing can effectively deal with the azimuth defocus effect of large-sized targets.

[0182] The above calculation examples of the present invention are only used to explain the calculation model and calculation process of the present invention in detail, and are not intended to limit the implementation methods of the present invention. For ordinary technicians in the relevant field, other different forms of changes or modifications can be made based on the above description. It is impossible to list all the implementation methods here. All obvious changes or modifications derived from the technical solution of the present invention are still within the scope of protection of the present invention.

Claims

1. A bistatic inverse synthetic aperture radar imaging and calibration method based on motion parameter estimation, characterized in that: The method specifically comprises the following steps: Step 1: The bistatic inverse synthetic aperture radar obtains echo data reflected by scattering points on the uniformly rotating target, and then performs distance compression and motion compensation processing on the echo data to obtain echo data reflected by each scattering point on the uniformly rotating target after distance compression and motion compensation processing; Step 2: Based on the echo data after range compression and motion compensation processing of each scattering point on the uniformly rotating target, obtain the estimated value of the ratio of the first-order coefficient to the constant coefficient in the dual-base coefficient and an estimate of the square of the effective angular velocity Step 3: Calculate the bistatic angle at the initial moment, and then estimate the constant coefficient K0 in the bistatic coefficient based on the bistatic angle at the initial moment; Step 4: Based on the estimated value and Construct the phase compensation term s of the i-th distance unit com,i , and then multiplying the phase compensation term of the i-th distance unit with the echo data after range compression and motion compensation corresponding to the i-th distance unit to obtain the echo data of the i-th distance unit after phase compensation; Step 5, compressing the echo data of all range units after phase compensation in step 4 in azimuth direction to obtain echo data compressed in azimuth direction; Step 6: Obtain a scattering point echo map based on the echo data compressed in azimuth; Step 7: Obtain the target azimuth dimension y according to the scattering point echo map obtained in step 6 max , and then according to the constant coefficient K0 and the estimated value Calculate K1 and determine the target orientation size y max Is it satisfied? Among them, κ w is the widening factor of the window function; If it is not satisfied, the scattering point echo map is post-compensated and the post-compensated image is used as the inverse synthetic aperture radar calibration image; If it is satisfied, the scattering point echo image obtained in step 6 is the inverse synthetic aperture radar calibration image.

2. The method for bistatic inverse synthetic aperture radar imaging and calibration based on motion parameter estimation according to claim 1, characterized in that: The echo data after distance compression and motion compensation processing reflected by each scattering point on the uniformly rotating target is specifically: Where σ is the bistatic scattering coefficient of the scattering point, B is the signal bandwidth, and T obs is the observation duration, rect(·) is the rectangular window function, sinc(·) is the sinc function, j is the imaginary unit, λ is the carrier wavelength, t is the slow time, τ is the fast time delay variable, c is the speed of light, and f c is the center frequency, ξ, α, β and γ are the phase coefficients; Among them, K0 and K1 are the constant coefficient and the first-order coefficient in the dual base coefficient, respectively, (x, y) is the coordinate of the scattering point relative to the rotation center, Ω eff is the effective angular velocity.

3. The method for bistatic inverse synthetic aperture radar imaging and calibration based on motion parameter estimation according to claim 2, characterized in that: The specific process of step 2 is as follows: Step 21: selecting the distance unit based on the minimum normalized amplitude variance criterion and the processed echo data reflected from each scattering point; Step 2: Determine the scattering point with the largest echo energy in each selected distance unit, and then use the polynomial phase estimation method to estimate the phase coefficient corresponding to each determined scattering point, and record the phase coefficient corresponding to the scattering point with the largest echo energy in the kth distance unit as ξ k ,α k ,β k ,γ k , k=1,2,…,K; Step 23: Use the phase coefficients estimated in step 22 to construct the equation system: Among them, x k '=K0x k ,y k '=K0y k ,(x k ,y k ) is the coordinate value of the scattering point with the largest echo energy in the kth range unit; Step 24: Convert the equation system of formula (3) into the following form: make but: Step 25: Substitute any set of phase coefficients obtained in step 22 into equation (5), and then solve equation (5) to obtain the three roots corresponding to the set of phase coefficients; similarly, solve the three roots corresponding to each set of phase coefficients obtained in step 22 respectively; Then, according to all the positive roots obtained, we can get Estimated value of Step 26: Estimated value Substitute into the equation group of formula (3), calculate a ratio according to each set of phase coefficients, and then calculate the mean of each ratio, and use the calculated mean as the estimation result 4. The method for bistatic inverse synthetic aperture radar imaging and calibration based on motion parameter estimation according to claim 3, characterized in that: The specific process of step 21 is as follows: Normalized amplitude variance NAV of the i-th range unit i for: Among them, s i (n) is the processed echo data corresponding to the nth scattering point in the ith range unit, E n (·) is the mean operation, |s i (n)| means s i (n) magnitude; Similarly, after calculating the normalized amplitude variance of each distance unit, the normalized amplitude variances of each distance unit are sorted in order from small to large, and the distance units corresponding to the normalized amplitude variances in the first K positions are selected.

5. The method for bistatic inverse synthetic aperture radar imaging and calibration based on motion parameter estimation according to claim 4, characterized in that: In the step 25, according to all the positive roots obtained, Estimated value of Specifically: Cluster all positive roots, select the cluster with the largest density from the clustering results, calculate the mean of all positive roots in the selected cluster, and use the calculated mean as the estimated value 6. A bistatic inverse synthetic aperture radar imaging and calibration method based on motion parameter estimation according to claim 5, characterized in that: The specific process of step three is: Let T, R and O represent the positions of radar transmitter Tx, receiver Rx and target rotation center respectively, and let the line connecting T and R be represents the baseline of the bistatic inverse synthetic aperture radar, and the straight-line distance between T and R is L, and the straight-line distance between T and O is R T , let the line connecting R and O be The straight-line distance between R and O is R R ; Will and The angle between Indicates the angle between the LOS direction of the receiving radar and the baseline direction, extended To point A, The length is equal to Then connect point A to T and make a line from R to The perpendicular line of , and the foot of the perpendicular is C; Calculate θ0 using the law of sines for a triangle: Among them, θ0 is the bistatic angle at the initial moment, R RT =R R +R T , that is, R RT Represents R T and R R The distance and; Then the constant coefficient K0 of the double base coefficient is: K0=cos(θ0 / 2) (10).

7. A bistatic inverse synthetic aperture radar imaging and calibration method based on motion parameter estimation according to claim 6, characterized in that: The specific process of step 4 is as follows: Among them, ξ i =τ i c / 2, i represents the distance unit index value, τ i is the fast time delay corresponding to the i-th range unit.

8. The method for bistatic inverse synthetic aperture radar imaging and calibration based on motion parameter estimation according to claim 7, characterized in that: The specific process of step five is as follows: Among them, s rc (τ,f d ) is the echo data after azimuth compression, f d is the Doppler variable.

9. A bistatic inverse synthetic aperture radar imaging and calibration method based on motion parameter estimation according to claim 8, characterized in that: The specific process of step six is ​​as follows: Among them, ΔX and ΔY are the range resolution and azimuth resolution of the bistatic inverse synthetic aperture radar image, respectively, and X and Y are the coordinates of the range dimension and azimuth dimension, respectively; Substituting equation (13) into equation (12) and taking the absolute value of the right side of equation (12) we can obtain the scattering point echo map of the bistatic inverse synthetic aperture radar: Among them, I(X,Y) is the scattering point echo pattern of the bistatic inverse synthetic aperture radar.

10. The method for bistatic inverse synthetic aperture radar imaging and calibration based on motion parameter estimation according to claim 9, characterized in that: The specific process of the post-compensation processing is as follows: Step 7: Segment the scattering point echo map obtained in step 6, that is, define the segmentation step length. Then, the scattering point echo map is divided into various regions along the azimuth according to the segmentation step size, and the azimuth intervals corresponding to each region are recorded as [-Ry p , -(R-1)y p ],(-(R-1)y p , -(R-2)y p ], ..., (-y p ,0],(0,y p ],…,((R-2)y p ,(R-1)y p ],((R-1)y p ,Ry p ], 2R is the total number of azimuth intervals; The azimuth interval is (ry p ,(r+1)y p ] is taken as the rth positive defocus area, and the azimuth interval is (-(r-1)y p , -(r-2)y p ] is taken as the rth negative defocus area, where r is a positive integer; The azimuth interval is (-y p ,0] and (0,y p ] is the non-defocused area; Step 72: extract the image of each defocused area, perform inverse Fourier transform on the image of each defocused area in the azimuth dimension, and convert the image from the range-azimuth domain to the range-slow time domain; The rth positive defocused area image in the range-slow time domain is combined with the post-compensation phase signal Multiply the rth negative defocus area image in the distance-slow time domain with the post-compensation phase signal Multiply, perform Fourier transform on each multiplied image along the slow time dimension, transform each defocused area image from the range-slow time domain to the range-azimuth domain, synthesize the transformed image of each defocused area with the image of the non-defocused area, and obtain the focused and calibrated inverse synthetic aperture radar image.

Images