Helicopter rotor feature extraction method based on scintillation correlation of external emitter radar

Abstract

The present invention discloses a method for extracting features of helicopter rotors from external radiation source radar based on flicker correlation. It includes the following steps: Step 1: Obtain the slow time dimension echo of the distance unit where the target is located; Step 2: Time-frequency analysis, pixel value accumulation and cross-correlation processing to determine the parity of the number of rotor blades; Step 3: Determine the reference frequency of the echo flicker; Step 4: Perform phase compensation on the echo; Step 5: Time-frequency analysis of the echo after phase compensation; Step 6: Save the center frequency of each flicker at each search node, calculate the Euclidean distance, find the minimum position, estimate the number of rotor blades and the initial phase of a single blade; Step 7: Verify the estimated value of the number of blades; Step 8: Estimate the blade length and the initial phase of all blades. The present invention overcomes the defects of the existing phase compensation method, which has a large amount of computation in the iterative traversal compensation process and is easy to fall into the local optimal solution; and has the advantage of being able to complete the feature extraction of helicopter targets.

Description

Technical Field

[0001] The invention relates to the field of external radiation source radar signal processing, and more specifically to a method for extracting helicopter rotor features of external radiation source radar based on scintillation correlation. Background Art

[0002] Exogenous radiation source radar is a new type of radar that uses electromagnetic signals emitted by a third-party non-cooperative radiation source to illuminate the target, and only passively receives the target's scattered signals to implement detection. Compared with traditional active radars, exogenous radiation source radars usually have natural excellent low-altitude coverage, strong slow and small target detection capabilities, and low network detection costs. Therefore, exogenous radiation source radars have broad application prospects in the military field.

[0003] In recent years, major countries in the world have successively formulated digital broadcasting and television standards with independent intellectual property rights, such as the Advanced Television System Committee (ATSC) of the United States, Digital Video Broadcasting-Terrestrial (DVB-T) of Europe, Integrated Services digital Broadcasting (ISDB) of Japan, and China Mobile Multimedia Broadcasting (CMMB) of China. Due to the differences in base station layout, transmission power and signal format of different broadcasting and television standards, the research on external radiation source radar must be combined with reality and adapted to local conditions. As a widely covered broadcasting and television signal in China, the basic structure of CMMB signal is signal frame, 1 second is 1 frame, divided into 40 time slots, each time slot contains 1 beacon and 53 orthogonal frequency division multiplexing (OFDM) data symbols. Specifically, the beacon contains 1 transmission identification signal and 2 synchronization signals, and each OFDM symbol has 4096 subcarriers. For the same transmitting station, the fixed components in its CMMB signal will not change over a long period of time, so the signal can be used as a stable external radiation source. The Doppler modulation phenomenon generated when the target undergoes non-rigid motion such as rotation and vibration is called the micro-Doppler effect, which is often used for the classification and identification of aerial targets. In the study of the micro-Doppler effect, external radiation source radar based on terrestrial wireless broadcasting and television has shown unique advantages: (1) The separation of transmission and reception can achieve spatial diversity and effectively avoid detection blind spots; (2) The third-party radiation source is mostly continuous wave, there is no close-range blind spot, and long-term coherent accumulation can observe multiple continuous echo flashes, which is conducive to improving the detection and identification capabilities of targets with low radar scattering cross-section; (3) The extraction of micro-Doppler features does not require high distance resolution, and parameter estimation is not limited by the bandwidth of the third-party radiation source.

[0004] Most of the existing methods for extracting features from external radiation source radar helicopters are completed using traditional Hough transform, orthogonal matching pursuit and phase compensation methods. The team led by Wan Xianrong of Wuhan University successfully extracted the parameters of the rotor from the echo of a dual-base radar helicopter using the orthogonal matching pursuit algorithm. However, the amount of computation for this method will explode as the estimated dimension increases, and there is a potential defect of grid mismatch. The team led by Zhang Qun of the Air Force Engineering University used time-frequency analysis and Hough transform to extract features from external radiation source radar targets based on long-term evolution signals, but this method fails when the sinusoidal modulation curve is not obvious in the time-frequency spectrum, and its scope of application is small. The team led by Yang Jun of the Air Force Early Warning Academy uses the phase compensation method to focus on the echo phase, iteratively searches for the number of target blades, and realizes the feature extraction of the rotor by calculating the offset between the center frequency of the target flicker and the reference frequency. However, this method is prone to fall into a local optimal solution, and the search process is slow.

[0005] Therefore, it is necessary to develop a method for extracting helicopter rotor features from external radiation source radar, which can overcome the defects of existing phase compensation methods, such as large amount of computation in the iterative ergodic compensation process and easy to fall into local optimal solutions. Summary of the invention

[0006] The purpose of the present invention is to provide a method for extracting features of helicopter rotors from an external radiation source radar based on flicker correlation. The method is a method for estimating and optimizing helicopter rotor parameters based on flicker correlation after phase compensation. The method can extract features of helicopter targets with small computational complexity, high estimation accuracy and high stability. The method overcomes the defects of existing phase compensation methods in that the iterative traversal compensation process has a large computational complexity and is prone to falling into local optimal solutions.

[0007] The method of the present invention first performs time-frequency analysis on the echo of the distance unit where the target is located, and accumulates pixel values ​​for the positive and negative frequency parts respectively, and determines the parity of the number of target blades through cross-correlation processing. Secondly, the reference frequency of the flicker in the target echo is estimated in the time-frequency diagram, and the search grid is divided according to the parity of the number of blades, the initial phase search range is set in segments, and the echo is phase compensated using the parameters of each search node. Then, by calculating the Euclidean distance of the compensated flicker center frequency, it is determined whether the flicker belongs to the same blade, so as to estimate the specific number of target blades. Finally, other parameters are estimated by combining the search node parameters and the reference frequency of the flicker in the original echo. The simulation experimental results show that this method can effectively complete the rotor feature extraction under flickering conditions.

[0008] In order to achieve the above object, the technical solution of the present invention is: a method for extracting helicopter rotor features using external radiation source radar based on scintillation correlation, characterized in that: by using an optimized phase compensation algorithm, the characteristic parameters of the rotor target can be quickly extracted;

[0009] The specific extraction method comprises the following steps:

[0010] Step 1: Obtain the slow-time echo of the distance unit where the target is located;

[0011] Step 2: Time-frequency analysis, pixel value accumulation and cross-correlation processing to determine the parity of the number of rotor blades;

[0012] Step 3: Determine the reference frequency of the echo flicker;

[0013] Step 4: Perform phase compensation on the echo;

[0014] Step 5: echo time-frequency analysis after phase compensation;

[0015] Step 6: Save the center frequency of each flash at each search node, calculate the Euclidean distance, find the minimum position, estimate the number of rotor blades and the initial phase of a single blade;

[0016] Step 7: Verify the estimated number of leaves;

[0017] Step 8: Estimate the leaf length and initial phase of all leaves.

[0018] In the above technical solution, in step one, pulse compression is performed on the original echo of the target to obtain the slow time dimension echo of the distance unit where the target is located, and the time when the flashing occurs is recorded.

[0019] Since the helicopter rotor is composed of multiple blades, and each blade contains multiple scattering points, referring to the rotor modeling method in the single-base radar system, the echo of a single blade can be obtained by integrating the scattering points, and the complete echo of the target rotor can be obtained by summing the echoes of multiple blades. Assuming that the target has I blades, the complete echo of the target rotor can be written as

[0020]

[0021] In formula (1), j is an imaginary unit, ρ r is the amplitude envelope of the received signal, dimensionless; λ n is the wavelength of the nth subcarrier frequency, in meters; For fast time, t m is the slow time, unit: second; θ i is the initial phase of the i-th blade, unit: radian; is the bistatic factor (the bistatic angle is β, and the angle between the bistatic angle bisector and the horizontal plane is the azimuth ), under the condition that the observation position and the target position are determined, this item is a constant and dimensionless; L is the equivalent length of the blade in the reference coordinate system, unit: meter; w is the equivalent rotation speed of the rotor in the reference coordinate system, unit: revolutions per second.

[0022] From formula (1), we can know that the micro-Doppler frequency (unit: Hertz) caused by the rotation of the i-th blade is:

[0023]

[0024] In formula (2), w is the equivalent rotation speed of the rotor in the reference coordinate system, unit: revolutions per second; L is the equivalent length of the blade in the reference coordinate system, unit: meter; λ n is the wavelength of the nth subcarrier frequency, in meters; t m is the slow time, unit: second; θ i is the initial phase of the i-th blade, unit: radian; is the bistatic factor (the bistatic angle is β, and the angle between the bistatic angle bisector and the horizontal plane is the azimuth ), under the condition that the observation position and target position are determined, this item is a constant and dimensionless;

[0025] In the above technical solution, in step 2, the Gabor transform is used to perform time-frequency analysis on the original echo of the target, and the pixel values ​​of the positive and negative frequency parts of the time-frequency diagram are accumulated, and the parity of the number of rotor blades is judged by the time delay value after cross-correlation processing.

[0026] The pixel values ​​of the positive and negative frequency parts of the time-frequency graph are accumulated along the frequency axis, and the positive frequency accumulation result is recorded as TFD. P (represents the amplitude matrix of the positive frequency accumulation result), and the negative frequency accumulation result is recorded as TFD N (representing the amplitude matrix of the negative frequency accumulation result), the cross-correlation processing of the accumulation result is:

[0027]

[0028] In formula (3): is the convolution operator symbol.

[0029] The parity of the target blade number can be determined by the delay value processed by the cross-correlation. If the delay value is 0, the target has an even number of blades; if the delay value is not 0, the target has an odd number of blades.

[0030] In the above technical solution, in step three, the reference frequency of the target original echo flicker is determined according to the parity of the number of rotor blades.

[0031] When the target rotor contains an odd number of blades, due to the asymmetry of the blades, positive and negative frequency flickers appear alternately in the time-frequency domain. Define κ as the amplitude factor of the echo phase, and select a flicker center in the original echo as the reference frequency (unit: Hz).

[0032]

[0033] In formula (4), w is the equivalent rotation speed of the rotor in the reference coordinate system, unit: revolutions per second; L is the equivalent length of the blade in the reference coordinate system, unit: meter; λ n is the wavelength of the nth subcarrier frequency, in meters; is the bistatic factor (the bistatic angle is β, and the angle between the bistatic angle bisector and the horizontal plane is the azimuth ), under the condition that the observation position and target position are determined, this item is a constant and dimensionless; f b is the reference frequency of the scintillation band, unit: Hertz.

[0034] When the target rotor contains an even number of blades, due to the symmetrical distribution of the blades, the positive and negative frequency flickers in the time-frequency domain appear simultaneously, forming a flicker pair. The flicker center of a single blade is selected as the reference frequency, which is one-fourth of the bandwidth of the entire flicker pair.

[0035]

[0036] In formula (5), w is the equivalent rotation speed of the rotor in the reference coordinate system, unit: revolutions per second; L is the equivalent length of the blade in the reference coordinate system, unit: meter; λ n is the wavelength of the nth subcarrier frequency, in meters; is the bistatic factor (the bistatic angle is β, and the angle between the bistatic angle bisector and the horizontal plane is the azimuth ), under the condition that the observation position and target position are determined, this item is a constant and dimensionless; f b is the reference frequency of the scintillation band, unit: Hertz.

[0037] It can be seen that by selecting the flicker center of a single blade as the reference frequency, the amplitude factor in the phase compensation operator under different blade parity conditions can be uniformly described, and the parameter κ (amplitude factor of the echo phase) is associated with w (equivalent rotation speed of the rotor in the reference coordinate system), so the phase compensation operator Υ can be written as

[0038]

[0039] In formula (6), j is an imaginary unit; f b is the reference frequency of the scintillation band, unit: Hertz; t m is the slow time, unit: seconds; μ is the angular velocity factor, ξ is the phase factor. At this time, the rotor feature extraction problem is transformed into a parameter estimation problem, and the parameter search dimension is reduced from three dimensions (κ, μ, ξ) to two dimensions (μ, ξ).

[0040] In the above technical solution, in step 4, the number of target search blades I is set in reverse order, the search grid (μ, ξ) is segmented to construct a two-dimensional search space, and the parameters on each search node are used to construct a phase compensation operator Υ to perform phase compensation on the echo.

[0041] For the parameter search process, if the entire two-dimensional search space is directly searched in an ergodic manner, the amount of calculation is large. However, since the length and rotation speed of all blades of the same helicopter rotor are the same, the only difference is the initial phase of each blade. If the parameters of one of the blades can be obtained, the relevant parameters of the entire rotor can be obtained through the relationship with the number of blades. Therefore, the algorithm described in the present invention considers using the phase compensation of a single blade to achieve the feature extraction of the entire rotor target. In the two-dimensional search space, the grid division of the angular velocity factor μ can only be ergodic, but the grid division of the phase factor ξ can be divided into sections by analyzing the relationship between the number of blades and the initial phase. When the number of blades on the rotor is more, the initial phase range of a single blade is smaller, and the search nodes are fewer. Therefore, the present invention searches the number of blades of the rotor in reverse order. This setting method can ensure that the number of ξ nodes in the grid is minimized during each search process, thereby reducing the amount of calculation.

[0042] The search range of rotation speed w is set to 3-7 rpm, with a step of 0.1 rpm. When the number of blades is an even number, the search range of blade number N is set to [8, 6, 4, 2]; when the number of blades is an odd number, the search range of blade number N is set to [7, 5, 3]. At the same time, the initial phase range is set to [0, 2π / N) with a step of π / 180. The phase compensation of the slow time dimension echo is performed by searching the grid nodes. The compensated echo is expressed as

[0043]

[0044] If and only if μ=w,ξ=θ i When the single blade echo phase is fully compensated,

[0045]

[0046] In formula (7) and formula (8), j is an imaginary unit; f b is the reference frequency of the scintillation band, unit: Hertz; μ is the angular velocity factor, ξ is the phase factor; t m is the slow time, unit: second; L is the equivalent length of the blade in the reference coordinate system, unit: meter; λ n is the wavelength of the nth subcarrier frequency, in meters; is the bistatic factor (the bistatic angle is β, and the angle between the bistatic angle bisector and the horizontal plane is the azimuth ), under the condition that the observation position and the target position are determined, this item is a constant and dimensionless; w is the equivalent rotation speed of the rotor in the reference coordinate system, unit: revolutions per second; θ i is the initial phase of the i-th blade, unit: radian;

[0047] Given that sinc(t m ) and rect(f) are a Fourier transform pair. The Fourier transform of the echo envelope after compensation is a rectangular strip occupying a certain bandwidth, that is, the time-frequency distribution of the single-blade echo after phase compensation is a series of zero-frequency flickers with a center frequency of zero.

[0048] In the above technical solution, in step five, a time-frequency analysis is performed on the compensated echo to obtain a time-frequency spectrum, and the center frequency of each flicker in the time-frequency spectrum is recorded.

[0049] In the above technical solution, in step 6, the Euclidean distance of the center frequency between each flash under different search nodes is obtained to determine whether the flash belongs to the same blade. At the same time, after summing in the ξ direction, the position with the smallest Euclidean distance of the three groups is found, and the number of rotor blades is estimated by finding the greatest common divisor.

[0050] When performing phase compensation on the rotor echo with an odd number of blades, it is necessary to determine whether a certain flash is focused to the zero frequency position at the search node. When performing phase compensation on the rotor echo with an even number of blades, it is necessary to observe whether the center of the flash is shifted by the reference frequency f. b To determine whether a certain flash is focused to the zero frequency position at the search node, and to prevent the algorithm from falling into the local optimal solution due to the unstable focusing effect of a single flash, the central frequency of all flashes at different nodes within the echo duration is estimated, which is recorded as C f If the flickers belong to the same leaf, then when one flicker is focused to zero frequency after phase compensation, the other flickers of the same leaf should also be focused to zero frequency at the same time. The correlation between the flickers can be determined by calculating the Euclidean distance of the center frequencies of different flickers in the two-dimensional search space.

[0051]

[0052] In formula (9): D(k) is the Euclidean distance; C f is a matrix used to record the center frequency of each flicker at different nodes (where C f {k,n,m} is a matrix used to record the center frequency of the kth flicker at different nodes; C f {k+1,n,m} is a matrix used to record the center frequency of the k+1th flicker at different nodes);

[0053] If the flicker belongs to the same blade, then the Euclidean distance D(k) takes the minimum value. The same blade may flicker multiple times during the observation time. Therefore, it is only necessary to find several flicker positions where the Euclidean distances simultaneously take the minimum value. By finding the greatest common divisor of the flicker positions, it is possible to verify whether this group of flickers occurs periodically, thereby estimating the number of rotor blades. The corresponding grid nodes are the estimated values ​​of the rotational speed and initial phase.

[0054] In the above technical solution, in step 7, verify the estimated value of the number of leaves Is it equal to the number of leaves I set in the search?

[0055] If they are equal, the estimate of the number of leaves is correct;

[0056] If they are not equal, return to step 4 and update the target search blade number and search grid.

[0057] In the above technical solution, in step eight, the rotation speed of a single blade on the target is estimated based on the search grid node by combining the prior information such as the target posture and the target position. and the initial phase of the first leaf Then calculate the blade length At the same time, based on the relationship with the number of leaves, the initial phase of all leaves is estimated.

[0058]

[0059] In formula (10): is the estimate of the reference frequency of the scintillation band, in Hertz; n is the wavelength of the nth subcarrier frequency, in meters; is the bistatic factor (the bistatic angle is β, and the angle between the bistatic angle bisector and the horizontal plane is the azimuth ), under the condition that the observation position and target position are determined, this item is a constant and dimensionless.

[0060] The present invention has the following advantages:

[0061] (1) The advantage of the method of the present invention is that the phase compensation process is optimized. Compared with the iterative traversal search process, the amount of computation is greatly reduced. The amount of computation in the entire feature extraction process mainly comes from the need to perform a time-frequency analysis for each search. Therefore, the amount of computation can be compared by comparing the number of time-frequency analyses used before and after the optimization of the search method (specifically, the number of times the traversal search process needs to be performed). time-frequency analysis, k is the number of iterative searches, N μ is the number of grid nodes for the angular velocity factor, N ξis the number of grid nodes of the phase factor. The actual amount of computation of the optimized algorithm is closely related to the target attribute. The more the number of target leaves, the smaller the amount of computation. When the number of target leaves I>2, the search process only needs It can be seen that the search amount after optimization of the present invention is only 2 / kI of the iterative traversal search;

[0062] (2) The present invention determines whether the flashes belong to the same blade by calculating the Euclidean distance between the flashes under the search node. At least three flashes belonging to the same blade need to be found. The number of rotor blades can be estimated by the greatest common divisor of the flash positions. The estimated result is tested in step seven, which can effectively prevent the situation of local optimal solution and avoid wrong estimation. Therefore, the present invention has high stability.

[0063] (3) Under the condition that the signal-to-noise ratio after pulse compression is greater than 0 dB, the blade number estimation of the present invention is completely accurate, the average error of the rotation speed estimation is less than 2%, and the average error of the blade length estimation is less than 3%. BRIEF DESCRIPTION OF THE DRAWINGS

[0064] Figure 1 This is the processing procedure of the algorithm in the present invention.

[0065] Figure 2 It is a flow chart of the present invention.

[0066] Figure 3 This is a diagram of the AH-64 echo phase compensation pre-processing result in simulation 1 of the present invention.

[0067] Figure 4 This is the time-frequency spectrum of the AH-64 echo after phase compensation in the simulation one of the present invention.

[0068] Figure 5 This is a diagram of the AH-64 parameter estimation results in simulation 1 of the present invention.

[0069] Figure 6 This is the Mi-28N echo phase compensation pre-processing result diagram in the second simulation of the present invention.

[0070] Figure 7 This is the time-frequency spectrum of the AH-64 echo after phase compensation in the second simulation of the present invention.

[0071] Figure 8 This is a diagram of Mi-28N parameter estimation results in simulation 2 of the present invention.

[0072] exist Figure 1 In the equation, ω represents the rotation speed, in units of r / s; represents the initial phase, unit is rad; k represents the kth flash, dimensionless; C fA matrix used to record the center frequency of each flicker at different nodes; D represents the Euclidean distance between the center frequency matrices of the first flicker and the kth flicker, dimensionless; loc represents the position, dimensionless; δ represents the greatest common divisor, dimensionless.

[0073] exist Figure 3 In the figure, (a) represents slow time dimension echo; (b) represents time-frequency spectrum; (c) represents pixel value accumulation; (d) represents cross-correlation processing.

[0074] exist Figure 5 In the figure, (a) shows the estimation results of the number of blades and the rotation speed; (b) shows the initial phase estimation results.

[0075] exist Figure 6 In the figure, (a) represents slow time dimension echo; (b) represents time-frequency spectrum; (c) represents pixel value accumulation; (d) represents cross-correlation processing.

[0076] Figure 8 In the figure, (a) shows the estimation results of the number of blades and the rotation speed; (b) shows the initial phase estimation results. DETAILED DESCRIPTION

[0077] The following is a detailed description of the implementation of the present invention in conjunction with the accompanying drawings, but they do not constitute a limitation of the present invention and are only given as examples. At the same time, the advantages of the present invention are made clearer and easier to understand through the description.

[0078] The present invention establishes a micro-motion model of the helicopter rotor of an external radiation source radar for CMMB signals, and proposes a helicopter rotor parameter estimation optimization method based on the flicker correlation after phase compensation in combination with the rotor structure characteristics and the flicker phenomenon in the time-frequency domain. First, the Gabor transform is used to perform time-frequency analysis on the target slow time dimension echo to obtain the time-frequency spectrum and flicker reference frequency of the echo, and then the pixel values ​​of the time-frequency graph are accumulated and the accumulated curve is cross-correlated to determine the parity of the target. Finally, the echo signal is phase compensated by dividing the two-dimensional search grid and constructing a phase compensation operator. The number of blades, rotation speed, initial phase and length of the target are estimated according to the phase compensation result. During the algorithm execution process, in order to obtain a smaller amount of calculation, the search is performed by setting the phase search range in segments. At the same time, in order to obtain higher algorithm stability, the correlation is determined by calculating the Euclidean distance between each flicker, so that the algorithm is more adaptable to targets with different rotor structures. Finally, through steps six and eight, simulation experiments are used to show that the algorithm in the present invention has a higher parameter estimation accuracy, which provides a solution for subsequent research on helicopter target recognition.

[0079] Referring to the attached drawings, it can be seen that the method for extracting helicopter rotor features based on scintillation correlation using external radiation source radar is characterized in that: by using an optimized phase compensation algorithm, the characteristic parameters of the rotor target can be quickly extracted;

[0080] The specific extraction method comprises the following steps:

[0081] Step 1: Obtain the slow-time echo of the distance unit where the target is located;

[0082] Step 2: Time-frequency analysis, pixel value accumulation and cross-correlation processing to determine the parity of the number of rotor blades;

[0083] Step 3: Determine the reference frequency of the echo flicker;

[0084] Step 4: Perform phase compensation on the echo;

[0085] Step 5: echo time-frequency analysis after phase compensation;

[0086] Step 6: Save the center frequency of each flash at each search node, calculate the Euclidean distance, find the minimum position, estimate the number of rotor blades and the initial phase of a single blade;

[0087] Step 7: Verify the estimated number of leaves;

[0088] Step 8: Estimate the leaf length and initial phase of all leaves.

[0089] Furthermore, in step one, pulse compression is performed on the original echo of the target to obtain the slow time dimension echo of the distance unit where the target is located, and the moment when the flashing occurs is recorded to provide an object for algorithm processing.

[0090] Furthermore, in step 2, the Gabor transform is used to perform time-frequency analysis on the echo, and the pixel values ​​of the positive and negative frequency parts of the time-frequency graph are accumulated, and the parity of the number of rotor blades is determined by the time delay value after cross-correlation processing; the above three methods are used simultaneously to determine the parity of the number of target blades; the specific method is:

[0091] The cross-correlation processing of the pixel value accumulation results is:

[0092]

[0093] In formula (3): is the convolution operator symbol; TFD P The amplitude matrix representing the accumulation result of positive frequencies; TFD N The magnitude matrix representing the accumulation result of negative frequencies;

[0094] The parity of the target blade number is determined by the delay value processed by the cross-correlation. If the delay value is 0, the target has an even number of blades; if the delay value is not 0, the target has an odd number of blades.

[0095] Furthermore, in step 3, the reference frequency of the echo flicker is determined according to the parity of the number of rotor blades, providing a reference standard for the phase compensation algorithm; when the target rotor contains an odd number of blades, due to the asymmetry of the blades, positive and negative frequency flickers appear alternately in the time-frequency domain; κ is defined as the amplitude factor of the echo phase, and a flicker center in the original echo is selected as the reference frequency.

[0096]

[0097] In formula (4), w is the equivalent rotation speed of the rotor in the reference coordinate system, unit: revolutions per second; L is the equivalent length of the blade in the reference coordinate system, unit: meter; λ n is the wavelength of the nth subcarrier frequency, in meters; is the dual-base factor, which is a constant when the observation position and target position are determined; f b is the reference frequency of the scintillation band, in Hertz; f (i) max is the maximum micro-Doppler frequency of the i-th glint or glint pair.

[0098] When the target rotor contains an even number of blades, the flicker center of a single blade is selected as the reference frequency, which is one-fourth of the entire flicker pair bandwidth.

[0099]

[0100] In formula (5), w is the equivalent rotation speed of the rotor in the reference coordinate system, unit: revolutions per second; L is the equivalent length of the blade in the reference coordinate system, unit: meter; λ n is the wavelength of the nth subcarrier frequency, in meters; is the dual-base factor, which is a constant when the observation position and target position are determined; f b is the reference frequency of the scintillation band, in Hertz; f (i) max is the maximum micro-Doppler frequency of the ith glint or glint pair; f (i) min is the minimum micro-Doppler frequency of the i-th scintillation or scintillation pair.

[0101] By selecting the flicker center of a single leaf as the reference frequency, the amplitude factor in the phase compensation operator under different leaf parity conditions is uniformly described, and the parameter κ is associated with w, so the phase compensation operator Υ is written as

[0102]

[0103] In formula (6), j is an imaginary unit; f b is the reference frequency of the scintillation band, unit: Hertz; tm is the slow time, unit: second; μ is the angular velocity factor, ξ is the phase factor.

[0104] Furthermore, in step 4, the number of target search blades I is set in reverse order, the search grid (μ, ξ) is segmented to construct a two-dimensional search space, and the phase compensation operator Υ is constructed using the parameters on each search node to perform phase compensation on the echo; the potential problem of explosive growth of the amount of calculation and grid mismatch in the prior art with the increase of the estimation dimension is solved;

[0105] The echo after compensation is expressed as

[0106]

[0107] If and only if μ=w,ξ=θ i When the single blade echo phase is fully compensated,

[0108]

[0109] In formula (7) and formula (8), j is an imaginary unit; f b is the reference frequency of the scintillation band, unit: Hertz; μ is the angular velocity factor, ξ is the phase factor; t m is the slow time, unit: second; L is the equivalent length of the blade in the reference coordinate system, unit: meter; λ n is the wavelength of the nth subcarrier frequency, in meters; is the dual-base factor. Under the condition that the observation position and the target position are determined, this item is a constant and dimensionless; w is the equivalent rotation speed of the rotor in the reference coordinate system, unit: revolutions per second; θ i is the initial phase of the ith blade, in radians.

[0110] Furthermore, in step five, a time-frequency analysis is performed on the compensated echo to obtain a time-frequency spectrum, and the center frequency of each flicker in the time-frequency spectrum is recorded; the time-frequency spectrum after phase compensation is obtained to facilitate the determination of the center frequency of each flicker; the problem that this method is ineffective and has a small scope of application when the sinusoidal modulation curve is not obvious in the time-frequency spectrum is solved.

[0111] Furthermore, in step 6, the Euclidean distance of the center frequency between each flash under different search nodes is calculated to determine whether the flash belongs to the same blade. At the same time, after summing in the ξ direction, the positions with the smallest Euclidean distances of the three groups are found, and the number of rotor blades is estimated by finding the greatest common divisor. The present invention determines the correlation of flickers by minimizing the Euclidean distance, that is, finding out which flickers belong to the same leaf, which is convenient for estimating the number of leaves and the initial phase of a single leaf. Compared with the existing methods, the present invention has the following characteristics: 1) the search method is improved (two-dimensional reverse search is adopted, and the phase search range varies with the number of leaves), which reduces the amount of calculation; 2) the correlation between flickers is used (by judging the Euclidean distance of the flicker center frequency of each node) to estimate the target leaf;

[0112] The calculation method of Euclidean distance is as follows:

[0113]

[0114] In formula (9): D(k) is the Euclidean distance; C f {k,n,m} is a matrix used to record the center frequency of the kth flicker at different nodes; C f {k+1,n,m} is a matrix used to record the center frequency of the k+1th flicker at different nodes.

[0115] Further, in step 7, the estimated number of leaves is verified Is it equal to the number of leaves I set in the search?

[0116] If they are equal, the estimate of the number of leaves is correct;

[0117] If they are not equal, return to the fourth step and update the target search blade quantity and search grid; this step prevents inaccurate estimation and achieves high stability of the present invention; and solves the problem that the prior art is prone to fall into local optimal solutions and the search process is slow.

[0118] Furthermore, in step eight, the target attitude and target position and other prior information are combined to estimate the rotation speed of a single blade on the target based on the search grid node. And the initial phase Then calculate the blade length At the same time, based on the relationship with the number of leaves, the initial phase of all leaves is estimated. Complete the estimation of all characteristic parameters of helicopter rotors using external radiation source radar;

[0119] The blade length The calculation method is as follows:

[0120]

[0121] In formula (10): is the estimate of the reference frequency of the scintillation band, in Hertz; n is the wavelength of the nth subcarrier frequency, in meters; is the dual-base factor. Under the condition that the observation position and target position are determined, this item is a constant and dimensionless.

[0122] Simulation Analysis

[0123] In order to verify the effectiveness of the feature extraction method proposed in this invention, the simulation parameters are set as follows: a broadcast signal transmitting station is located at the origin of the radar reference system, the receiving station is located at (8000, 0, 0), the coordinates of the target are (10000, 10000, 2000), and the initial Euler angle between the target coordinate system and the reference coordinate system is assumed to be (0, 0, 0), that is, the target attitude is horizontal. The target parameters select the real parameters of the AH-64 helicopter and the Mi-28N helicopter. At the same time, in order to simulate the weak characteristics of the external radiation source radar echo, the signal-to-noise ratio after Dechirp processing is set to 0dB.

[0124] Table 1 Simulation parameter settings

[0125] Radar parameters value Target parameters value Carrier frequency f(MHz) 658 Number of blades 4、5 Number of subcarrier frequencies 4096 Blade length L(m) 7.3、8.6 Subcarrier frequency spacing Δf (kHz) 1.95 <![CDATA[Rotational angular velocity w (r·s -1 )]]> 4.8、4 Subcarrier frequency coefficient 1 <![CDATA[Translational velocity v (m·s -1 )]]> 80 Repetition period PRF(Hz) 2160 Scattering point coefficient 1

[0126] Simulation 1:

[0127] The simulation target of this part is the AH-64 helicopter. The slow time dimension echo processing results before phase compensation are as follows: Figure 3 As shown, at this time, the initial phases of the four leaves are set to [0°, 90°, 180°, 270°]. Figure 3 (a) is the slow time dimension echo of the target range unit. It can be observed that the target has 9 complete flashes. Due to the initial phase setting of the rotor blade, the first flash was not completely observed. Figure 3 (b) is the time-frequency spectrum of the slow-time echo after Gabor transformation. Due to the low signal-to-noise ratio of the echo, the sinusoidal modulation curve can no longer be observed in the time-frequency spectrum. It can be seen that it is not feasible to use the sinusoidal modulation characteristics in the echo to extract the features of the helicopter rotor. However, the occurrence of time-frequency domain flicker can still be observed in the time-frequency spectrum. The number and position are consistent with the time-domain flicker. The reference frequency for extracting flicker is 365.3Hz. The pixel values ​​of the positive and negative frequency parts of the time-frequency spectrum are accumulated along the frequency axis. The accumulation results are shown in Figure 2. Figure 3 (c) is shown. Figure 3 (c) Whether the positive frequency pixel value accumulation curve and the negative frequency pixel value accumulation curve appear at the same time, cross-correlation processing is performed on them. The processing results are as follows: Figure 3 (d) Figure 3 (d) shows that the time delay of the two pixel value accumulation curves is 0, that is, the symmetry of the positive frequency part and the negative frequency part in the time-frequency diagram is correctly judged, thereby determining that the number of target rotor blades is an even number of blades.

[0128] After completing the above processing, the rotor speed search range is divided into μ∈[3,7]r·s-1 , step Δμ=0.1, and the preset blade number range is I∈[8,6,4,2]. When searching for the preset value of the target blade I=8, the initial phase search range can be set to ξ∈[0,π / 4], step Δξ=π / 180. At this time, the search grid has a total of 1845 nodes. Different phase compensation operators are constructed using the parameters of each node to perform phase compensation on the target slow-time echo. In this search, the estimated value of the target rotor blade number is 4, which is inconsistent with the search preset value. At this time, only the determination that the target rotor blade number is not 8 is made. Therefore, the preset value of the target number of blades is updated to I=6, and the initial phase search range is updated to ξ∈[π / 4,π / 3]. At this time, the search grid has a total of 615 nodes. After completing the phase compensation of all search nodes, the algorithm has obtained the phase compensation result of ξ∈[0,π / 3]. According to the relationship between the number of blades and the initial phase, if the target number of rotor blades is 6, there must be an initial phase of a blade within this phase range. However, the estimated value of the target number of rotor blades is 4, which is still inconsistent with the search preset value. Therefore, it is determined that the target number of rotor blades is not 6. Continue to update the preset value of the target number of blades I=4, and update the initial phase search range to ξ∈[π / 3,π / 2]. At this time, the search grid has a total of 1230 nodes. Combined with the previous two search results, a set of parameters that can focus multiple flashes to zero frequency at the same time is successfully searched within the range of ξ∈[0,π / 2]. The time-frequency spectrum is shown as follows. Figure 4 As shown, at this time, the center frequencies of the 1st, 3rd, 5th, 7th, and 9th flash pairs are shifted by the reference frequency f b , that is, the flashes with a single blade in these flash pairs are focused to the zero frequency position.

[0129] from Figure 5 (a) It can be seen that after calculating the greatest common divisor of the minimum Euclidean distance between different flashes under the search grid, it is found that the greatest common divisor is 2, that is, the estimated value of the number of blades is 4, which is consistent with the preset value, so the number of target rotor blades is finally determined to be 4. It is worth pointing out that in the above three search processes, the estimated value of the number of rotor blades is 4, but the first two times were judged to be invalid. This is because the initial phase of a blade on the target rotor is 0, which is within the first phase search range. Therefore, the phase compensation of the blade was completed in the first search process. The reason why it is not used is that the estimated result has not been tested and its accuracy cannot be guaranteed. At this time, the speed search nodes are the 18th and 19th. The two different speed estimates appear because when the search grid is divided finely, the parameters adjacent to the true value can also better focus the flash to the zero frequency position. When this happens, the average value method can be used to obtain the final estimate. At this time, the speed estimate is 4.75r·s -1, the error is 1.04%. At the same time, under the above speed estimation value, a single flicker that can be focused to zero frequency is selected, and the phase search node that can achieve the best phase focusing effect is found. The parameter value of this node is the initial phase estimation value of the blade. Figure 5 As shown in (b), the estimated value of the initial phase of the blade is 2°. Combined with the estimated value of the target blade number, it can be estimated that the initial phase of each blade of the entire rotor is (2°, 92°, 182°, 272°). Finally, the target blade length is calculated to be 7.34m with an error of 0.55%.

[0130] Simulation 2:

[0131] The simulation target of this part is the Mi-28N helicopter. The slow time dimension echo processing results before phase compensation are as follows: Figure 6 To further illustrate the segmented search process of the algorithm, the initial phase of the blade is not within the first phase search range [0, 2π / 7]. Specifically, the initial phases of the five blades are set to [65°, 137°, 209°, 281°, 353°]. Figure 6 (a) is the slow time dimension echo of the target distance unit, where 20 complete flashes of the target can be observed. Figure 6 (b) is the corresponding time-frequency spectrum, from which the reference frequency of flicker is extracted as 361.3 Hz. The pixel value accumulation result is as follows: Figure 6 (c) shows the cross-correlation processing result. Figure 6 As shown in (d), since the time delay between the two pixel value accumulation curves is not 0, the target number of rotor blades is an odd number of blades.

[0132] The rotor speed search grid is the same as that in simulation 1, and the range of the number of blades is I∈[7,5,3]. When searching for the target blade preset value I=7, the initial phase search range is set to ξ∈[0,2π / 7], and the step Δξ=π / 180. At this time, the search grid has a total of 2091 nodes. By performing phase compensation on the target slow-time echo, it is estimated that the number of target rotor blades is 1. This indicates that no correlation is found between the flashes within this search range, and the number of blades is not 7. Update the preset value of the target number of blades I=5, and update the initial phase search range to ξ∈[2π / 7,2π / 5]. At this time, the search grid has a total of 861 nodes. At this time, the estimated value of the target rotor blade number is 5, which is consistent with the search preset value. Therefore, the number of rotor blades is successfully estimated. The time-frequency spectrum after phase compensation is as follows: Figure 7 As shown, at this time, the 1st, 6th, 11th, and 16th flashes are successfully focused to the zero frequency position.

[0133] Figure 8(a) It can be seen that there are four speed search nodes that can correctly estimate the number of blades, namely the 10th, 11th, 12th and 13th nodes. After averaging, the speed estimation value is 4.05 r·s -1 , with an error of 1.24%. Figure 8 (b) It can be seen that the estimated initial phase value of a blade on the rotor is 65°. Combined with the target number of blades being 5, it can be estimated that the initial phase of the entire rotor blade is [65°, 137°, 209°, 281°, 353°]. Finally, the target blade length is calculated to be 8.4m with an error of 2.33%.

[0134] It can be seen from the parameter estimation results of the two groups of simulations that the algorithm in the present invention can accurately estimate the number of rotor blades, and the estimation accuracy of the rotor speed, blade length and phase is relatively high. Under the same search accuracy, the algorithm has a small amount of computation. In the prior art, the estimation accuracy of the OMP algorithm can be controlled within 1% when the grid is not mismatched, but once there is a mismatch, the target parameters cannot be successfully estimated; the Hough transform and IRadon transform methods cannot use flicker information, and the estimation error is greater than 5% or even fails under the simulation conditions of the present invention; the traditional iterative phase compensation method has similar accuracy to the present method, but it may fall into a local optimal solution, resulting in estimation errors and the computational complexity is dozens of times that of the present invention.

[0135] Other parts not described belong to the prior art.

Claims

1. A method for extracting helicopter rotor features from external radiation source radar based on scintillation correlation. Features: The optimized phase compensation algorithm is used to quickly extract the characteristic parameters of the rotor target; The specific extraction method comprises the following steps: Step 1: Obtain the slow-time echo of the distance unit where the target is located; Step 2: Time-frequency analysis, pixel value accumulation and cross-correlation processing to determine the parity of the number of rotor blades; Step 3: Determine the reference frequency of the echo flicker; The reference frequency of the echo flash is determined by the odd or even number of rotor blades: When the target rotor contains an odd number of blades, due to the asymmetry of the blades, positive and negative frequency flickers appear alternately in the time-frequency domain; κ is defined as the amplitude factor of the echo phase, and a flicker center in the original echo is selected as the reference frequency of the flicker band. In formula (4), w is the equivalent rotation speed of the rotor in the reference coordinate system, unit: revolutions per second; L is the equivalent length of the blade in the reference coordinate system, unit: meters; λ n is the wavelength of the nth subcarrier frequency, in meters; is the bistatic factor, the bistatic angle is β, and the angle between the bistatic angle bisector and the horizontal plane is the azimuth Under the condition that the observation position and target position are determined, this item is a constant; f b is the reference frequency of the scintillation band, unit: Hertz; f (i) max is the maximum micro-Doppler frequency of the i-th scintillation or scintillation pair; When the target rotor contains an even number of blades, the flicker center of a single blade is selected as the reference frequency, which is one-fourth of the entire flicker pair bandwidth. In formula (5): f (i) min is the minimum micro-Doppler frequency of the i-th scintillation or scintillation pair; By selecting the flicker center of a single leaf as the reference frequency, the amplitude factor in the phase compensation operator under different leaf parity conditions is uniformly described, and the parameter κ is associated with w, so the phase compensation operator Υ is written as In formula (6), j is an imaginary unit; t m is the slow time, unit: second; μ is the angular velocity factor, ξ is the phase factor; Step 4: Perform phase compensation on the echo; The number of target search blades I is set in reverse order, the search grid (μ, ξ) is segmented to construct a two-dimensional search space, and the phase compensation operator Υ is constructed using the parameters on each search node to perform phase compensation on the echo. The compensated echo is expressed as If and only if μ=w,ξ=θ i When the single blade echo phase is fully compensated, In formula (7) and formula (8): n is the wavelength of the nth subcarrier frequency, in meters; θ i is the initial phase of the i-th blade, unit: radian; Step 5: echo time-frequency analysis after phase compensation; Step 6: Save the center frequency of each flash at each search node, calculate the Euclidean distance, find the minimum position, estimate the number of rotor blades and the initial phase of a single blade; The Euclidean distance of the center frequency between each flash under different search nodes is calculated to determine whether the flash belongs to the same blade. At the same time, the sum is taken in the ξ direction to find the position with the smallest Euclidean distance among the three groups. The number of rotor blades is estimated by finding the greatest common divisor. The calculation method of Euclidean distance is as follows: In formula (9): D(k) is the Euclidean distance; C f {k,n,m} is a matrix used to record the center frequency of the kth flicker at different nodes; C f {k+1,n,m} is a matrix used to record the center frequency of the k+1th flicker at different nodes; Step 7: Verify the estimated number of leaves; Step 8: Estimate the leaf length and initial phase of all leaves.

2. The method for extracting helicopter rotor features from external radiation source radar based on scintillation correlation according to claim 1, Features: In step 1, pulse compression is performed on the original echo of the target to obtain the slow-time echo of the distance unit where the target is located, and the moment when the flash occurs is recorded.

3. The method for extracting helicopter rotor features from external radiation source radar based on scintillation correlation according to claim 1 or 2, Features: In step 2, the Gabor transform is used to perform time-frequency analysis on the echo, and the pixel values ​​of the positive and negative frequency parts of the time-frequency graph are accumulated. The parity of the number of rotor blades is determined by the time delay value after cross-correlation processing. The specific method is as follows: The cross-correlation processing of the pixel value accumulation results is: In formula (3): is the convolution operator symbol; TFD P The amplitude matrix representing the accumulation result of positive frequencies; TFD N The magnitude matrix representing the accumulation result of negative frequencies; The parity of the target blade number is determined by the delay value processed by the cross-correlation. If the delay value is 0, the target has an even number of blades; if the delay value is not 0, the target has an odd number of blades.

4. The method for extracting helicopter rotor features from external radiation source radar based on scintillation correlation according to claim 1, Features: In step five, a time-frequency analysis is performed on the compensated echo to obtain a time-frequency spectrum, and the center frequency of each flicker in the time-frequency spectrum is recorded.

5. The method for extracting helicopter rotor features from external radiation source radar based on scintillation correlation according to claim 1, Features: In step 7, verify the estimated number of leaves Is it equal to the number of leaves I set in the search? If they are equal, the estimate of the number of leaves is correct; If they are not equal, return to step 4 and update the target search blade number and search grid.

6. The method for extracting helicopter rotor features from external radiation source radar based on scintillation correlation according to claim 5, Features: In step eight, the target attitude and the prior information of the target position are combined to estimate the rotation speed of a single blade on the target according to the search grid node. And the initial phase Then calculate the blade length At the same time, based on the relationship with the number of leaves, the initial phase of all leaves is estimated. The blade length The calculation method is as follows: In formula (10): is the estimate of the reference frequency of the scintillation band, in Hertz; n is the wavelength of the nth subcarrier frequency, in meters; is the dual-base factor. Under the condition that the observation position and target position are determined, this item is a constant and dimensionless.

Images