Micro-motion feature extraction method for multiple UAVs based on synchronous compression transformation and ridge reorganization

Abstract

The present invention discloses a method for extracting micro-Doppler of a UAV based on synchronous compression transformation and ridge path recombination, and belongs to the field of signal detection technology. The present invention improves the resolution of the separation process by integrating synchronous compression transformation into the signal separation strategy of ridge path recombination, establishes a multi-UAV rotor echo model, and uses the algorithm in the signal separation strategy of the present invention to extract the micro-Doppler features of the rotors of multiple UAVs in the constructed echo model, thereby realizing the identification of multiple UAVs. The present invention integrates synchronous compression transformation into the "ridge path recombination" method for separating overlapping non-stationary signals, thereby improving the resolution of the ridge path recombination method. This improved signal separation strategy is used to extract the micro-motion features of the rotors of multiple UAVs, extracts parameter information such as the rotor speed and blade length of each UAV, and realizes the identification of multiple UAVs.

Description

Technical Field

[0001] The present invention belongs to the technical field of signal detection, and in particular relates to a method for extracting micro-Doppler of an unmanned aerial vehicle based on synchronous compression transformation and ridge path reconstruction. Background Art

[0002] With the rapid growth in the number of drones, their widespread use has brought convenience to people's lives, but has also posed a great threat to public safety and citizens' privacy. Therefore, research on drone detection, target recognition, etc. is currently a hot topic.

[0003] When detecting drones, in addition to the movement of the drone body, the rotation of the drone rotor blades will also cause sinusoidal modulation of the radar echo phase, which VCChen calls the micro-Doppler effect, also known as micro-motion features. Since drones have this characteristic, the micro-motion features of drone rotor blades are an effective method for detecting drones. And through the extracted micro-motion features, the parameters such as the rotor speed and blade length of the drone can also be calculated, which can be used to identify different drones. The commonly used method for extracting micro-motion features is time-frequency analysis, such as wavelet transform and short-time Fourier transform. The time-frequency representation often contains multiple components. In order to extract the information of each component, it is necessary to separate the components in the time-frequency domain. Commonly used signal separation methods include empirical mode decomposition and intrinsic chirp component signal separation. For overlapping signals in the time-frequency domain, these commonly used separation methods are not effective. Shiqian Chen et al. proposed a ridge path recombination method to separate overlapping signals in the time-frequency domain. However, this method has not been used in the specific scenario of extracting micro-motion features of multiple drone rotors. When extracting the micro-motion features of drone rotors, it often contains information about more than one rotor blade. For micro-motion features in multiple drone scenes, it contains multiple components. In order to analyze each component, it is necessary to separate each component. Moreover, the resolution of this method is not high enough, and the error needs to be improved when extracting the rotor speed and blade length of drones. Summary of the invention

[0004] The present invention provides a multi-UAV micro-Doppler feature extraction method based on synchronous compression transformation and ridge path recombination, so as to extract the micro-Doppler features of each UAV in a multi-UAV scenario, and then identify the rotor speed and blade length of each UAV.

[0005] The technical solution adopted by the present invention is as follows:

[0006] A method for extracting micro-motion features of multiple UAVs based on synchronous compression transformation and ridge reorganization.

[0007] Step 1, modeling the multi-UAV rotor echo;

[0008] The one-dimensional time domain rotor echo signal s0(t) is expressed as:

[0009]

[0010] Wherein, t represents the sampling point, i represents the UAV ID, and k represents the UAV component ID. The UAV components include the UAV body and the UAV rotor blades. represents the radar cross section (RCS) of the components of the drone, represents the distance between the UAV component and the radar, λ represents the wavelength of the radar emission signal, and e represents the natural base;

[0011] Furthermore, when k represents the drone body (for example, when k=0, it represents the drone body), then represents the RCS of the main part of UAV i, which is a preset value, and the preferred value can be set to 3 times the RCS of any UAV rotor blade of UAV i;

[0012] When k represents a rotor blade of the UAV, then represents the RCS of the kth blade of the rotor of UAV i, and its expression is:

[0013]

[0014]

[0015]

[0016] Among them, parameter a i =L i / 2, parameter b i =W i / 2,L i represents the rotor blade length of UAV i, W i is the blade width of UAV i, θ represents the pitch angle of the kth blade of the rotor of UAV i relative to the radar, represents the azimuth angle of the kth blade of the rotor of UAV i relative to the radar.

[0017] Step 2, perform synchronous compression transformation on the one-dimensional time domain rotor echo signal s0(t) to obtain its time-frequency representation, denoted as TF(t,f);

[0018] Step 3, extract all instantaneous frequency ridges in the time-frequency representation TF(t,f):

[0019] Find the time t corresponding to the maximum value in the time-frequency representation TF(t,f) max , frequency f max ;

[0020] Based on time t max , traverse the time axis from the left and right directions respectively, for each moment t k , find the time t k The time-frequency representation of TF(t,f) corresponds to the frequency value f at the maximum point k :

[0021] Based on the corresponding f at each moment k Get the frequency sequence [f1...,f k ,...,f l ] and as a frequency ridge IF i =[f1...,f k ,...,f l ] Where l represents the number of sampling points;

[0022] Frequency Ridge IF i The instantaneous frequency ridge is extracted by Fourier series smoothing fitting. And get the corresponding time domain signal s i (t), where represents the frequency value after fitting;

[0023] Subtract the time domain signal s from the one-dimensional time domain rotor echo signal s0(t) i (t), and obtain the processed one-dimensional time domain rotor echo signal s′0(t);

[0024] If ||s′0(t)||≥ε1, return to step 2; otherwise, end the extraction and obtain the instantaneous frequency ridge set based on all instantaneous frequency ridges Wherein, ε1 represents a preset threshold value. Preferably, ε1 can be set to 0.2 of the norm of the one-dimensional time domain rotor echo signal s0(t);

[0025] Step 4, find the intersection intervals of each frequency ridge and merge the intersection intervals with similar distances;

[0026] Based on the preset threshold ε2 (preferably the sampling frequency f s 30%), for the instantaneous frequency ridge set Any two instantaneous frequency ridges in and

[0027] Traverse each time t k , if the current time t k The two frequencies corresponding to and If satisfied Then store the time number k in the vector findindex;

[0028] For all the time numbers in the vector findindex, if the absolute difference between two adjacent numbers is less than the threshold ε3 and greater than the threshold ε4, then the time corresponding to these two numbers That is, the two instantaneous frequency ridges and The edge of a crossover interval of and The intersection interval of

[0029] Traversing the instantaneous frequency ridge set All any two instantaneous frequency ridge combinations in , obtain their corresponding crossover intervals, and obtain a crossover interval set;

[0030] For any two intersection intervals in the intersection interval set, if the Euclidean distance between their points is less than the threshold ε5, the current two intersection intervals are merged to obtain the intersection interval set after the intersection intervals are merged;

[0031] Step 5, reorganize the ridges of the crossover interval to obtain a reorganized instantaneous frequency ridge set;

[0032] For any two instantaneous frequency ridges in each crossover interval in the crossover interval set after the crossover interval is merged, and Compare the slopes at the edge of the crossover interval. If the absolute value of the slope difference is less than the threshold ε6, the two current instantaneous frequency ridges are and Reconnect the combination to obtain the reorganized instantaneous frequency ridge;

[0033] Step 6, the reorganized instantaneous frequency ridges are fitted with Fourier series to obtain their time domain signals, and the extracted component signals are obtained, thereby obtaining the UAV rotor micro-Doppler characteristics.

[0034] When the present invention is used for the identification of drones, the component signals extracted in step 6 can be subjected to synchronous compression transformation to obtain the time-frequency representation of the extracted component signals. Based on the time-frequency representation of each component signal, the rotation speed and blade length of the drone rotor represented by each component signal can be calculated, thereby realizing the identification of multiple drones.

[0035] Furthermore, based on the time-frequency representation of each component signal, the rotation speed and blade length of the drone rotor represented by each component signal can be calculated as follows: The maximum frequency deviation Δf obtained from the time-frequency representation of each component signal k , period T k , we can calculate the rotation speed of the UAV rotor as The rotor tip speed is v k =Δf k ×λ, and then the length of the UAV rotor blade is obtained as

[0036] The technical solution provided by the present invention brings at least the following beneficial effects:

[0037] The present invention establishes a multi-UAV rotor echo model, applies the signal separation method of ridge path recombination to the extraction of UAV rotor micro-Doppler features, and adds synchronous compression transformation to the processing of ridge path recombination to improve the resolution of signal separation. The blade length and rotor speed are calculated from the extracted multi-UAV rotor micro-Doppler features to identify the information of each UAV. BRIEF DESCRIPTION OF THE DRAWINGS

[0038] In order to more clearly illustrate the technical solutions in the embodiments of the present invention, the drawings required for use in the description of the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work.

[0039] Figure 1 It is a processing flow chart of a method for extracting micro-Doppler of a drone based on synchronous compression transformation and ridge path reconstruction provided by an embodiment of the present invention;

[0040] Figure 2 is a schematic diagram of the geometric relationship between one rotor of two UAVs and a radar in an embodiment of the present invention;

[0041] Figure 3 is a schematic diagram of the intersection interval of two instantaneous frequency ridges in an embodiment of the present invention;

[0042] Figure 4 In the embodiment of the present invention, Figure 2 The schematic diagram of the cross-interval merging is shown;

[0043] Figure 5 It is a schematic diagram of the reorganization and combination connection of the same component in an embodiment of the present invention. DETAILED DESCRIPTION

[0044] In order to make the objectives, technical solutions and advantages of the present invention more clear, the embodiments of the present invention will be further described in detail below with reference to the accompanying drawings.

[0045] The UAV micro-Doppler extraction method based on synchronous compression transformation and ridge path recombination provided by the embodiment of the present invention first establishes a rotor echo model in a multi-UAV scenario, integrates synchronous compression transformation into the signal separation method of ridge path recombination to improve its resolution, and uses this method of separating signals based on ridge path recombination of synchronous compression transformation to extract micro-Doppler of multi-UAV rotors, thereby identifying parameters such as rotor speed and rotor blade length of multiple UAVs.

[0046] like Figure 1 As shown, the specific implementation steps of the UAV micro-Doppler extraction method based on synchronous compression transformation and ridge path reconstruction provided by the embodiment of the present invention include:

[0047] Step 1: Model the multi-UAV rotor echoes.

[0048] In this embodiment, a rotor echo model is established in a two-UAV scenario. Consider two UAVs A and B. Each UAV has two rotor blades. The geometric relationship between one rotor of the two UAVs and the radar is as follows: Figure 2 In the figure, (X1, Y1, Z1) are the coordinates of the radar. is the coordinate of the center of one rotor of UAV A, is the coordinate of the center of one rotor of UAV B, θ A is the pitch angle of the rotor blade of UAV A relative to the radar, is the azimuth angle of the rotor blade of UAV A relative to the radar, θ B is the pitch angle of the rotor blades of UAV B relative to the radar, The azimuth of the rotor blades of UAV B relative to the radar, Ω A ,Ω B are the rotor speeds of drones A and B, L A , L B are the rotor blade lengths of UAV A and B, respectively, W A , W B are the blade widths of UAV A and UAV B respectively. The one-dimensional time domain rotor echo signal s0(t) is obtained:

[0049]

[0050] Among them, the sampling point t = t1, ..., t l , l represents the number of sampling points, i represents the UAV ID, k represents the UAV component ID (including the UAV body and the UAV rotor blades), represents the distance between the UAV component and the radar. In this embodiment, when k=0, that is, It represents the distance between the main part of UAV i (i is A or B) and the radar, which can be calculated based on the movement of the UAV relative to the radar. When k = 1 or 2, is the distance between the kth (k is 1 or 2)th blade of the rotor blade of UAV i (i is A or B) and the radar, which can be calculated from the motion of the UAV relative to the radar and the parameters of the rotor blades. λ represents the wavelength of the radar transmission signal, j represents the imaginary unit, The RCS (radar cross section) of the kth (k is 1 or 2) blade of the rotor of UAV i (i is A or B) is calculated as follows:

[0051]

[0052]

[0053]

[0054] Among them, a i =L i / 2, b i =W i / 2,L i is the blade length of UAV i, W i is the blade width of UAV i, θ is the pitch angle of the kth (k is 1 or 2)th blade of the rotor of UAV i (i is A or B) relative to the radar, is the azimuth angle of the kth (k is 1 or 2)th blade of the rotor of UAV i (i is A or B) relative to the radar. is the RCS of the main part of drone i (i is A or B), which can be set to It can also be set to other values ​​according to specific application scenarios.

[0055] Step 2: Use synchronous compression transformation on the time domain callback signal to obtain its time-frequency representation TFR.

[0056] That is, the one-dimensional time-domain rotor echo signal s0(t) is subjected to synchronous compression transformation to obtain its time-frequency representation, which is recorded as TF(t, f).

[0057] Step 3: Extract all instantaneous frequency ridges in the time-frequency representation.

[0058] (3-1) Find the time t corresponding to the maximum value in the time-frequency representation max , frequency f max ,Right now

[0059]

[0060] (3-2) From this moment, traverse the time axis from the left and right directions to find the moment t kThe time-frequency representation of TF(t, f) corresponds to the frequency value f at the maximum point k ,Right now

[0061]

[0062] Traverse all moments and find the f corresponding to each moment k , the obtained frequency sequence [f1...,f k , ..., f l ], which is a frequency ridge IF represented by the time-frequency i =[f1...,f k , ..., f l ].

[0063] (3-3) Use Fourier series to smooth the extracted instantaneous frequency ridge And get the corresponding time domain signal s i (t). Among them, Represents the frequency value after fitting.

[0064] Then, the extracted signal s is subtracted from the original echo time domain signal s0(t). i (t), and the processed one-dimensional time domain rotor echo signal s′0(t)=s0(t)-s i (t), if ||s′0(t)||≥ε1, where ε1 represents a preset threshold, preferably, ε1 can be set to 0.2 of the norm of the original time domain signal, then return to step 2; otherwise, end the extraction process and obtain all N instantaneous frequency ridges in the original time-frequency representation.

[0065] Step 4: Find the intersection intervals of each frequency ridge and merge the intersection intervals with similar distances.

[0066] Given a threshold ε2, this threshold is usually set to the sampling frequency f s 30%, find the intersection interval of any two instantaneous frequency ridges in the time-frequency diagram, as shown in the schematic diagram Figure 3 As shown, for any two instantaneous frequency ridges If for time t k The two frequencies corresponding to have Then store the number k in the findindex vector and find all the matching When this condition is met, the findindex vector is findindex=[k1,...,k f ]. For all labels in findindex, if two adjacent numbers k p , k qThe difference (absolute difference) is less than the preset threshold ε3 (ε3 can be 1 / 200) and greater than the threshold ε4 (ε4 can be 1 / 500), that is, the difference between two adjacent numbers ∈(ε4, ε3), then the time corresponding to these two numbers is The edge of a crossing interval of In this way, the set of all intersecting intervals is found.

[0067] Merge the intersection intervals with small distance, that is, the Euclidean distance between the midpoints of any two intersection intervals is less than the threshold ε5 (ε5 can be taken as 4(ε2) 2 ), the two intersecting intervals are merged. Figure 4 For example, for any crossover interval The time-frequency coordinate interval of the crossover interval in the time-frequency diagram is Take the midpoint of this interval If there is another time-frequency interval corresponding to the crossover interval The midpoint of the interval is And the distance between the time-frequency midpoints corresponding to these two crossover intervals is less than ε5, that is, The two intersecting intervals are merged into All intersecting intervals that can be merged are merged in this way.

[0068] Step 5, breaking the ridges in the crossover interval, reorganizing the ridges in the crossover interval according to the slope of the ridges, and obtaining a reorganized instantaneous frequency ridge set.

[0069] For each frequency ridge in the crossover interval, reconnect them at the crossover interval position, and reconnect and combine them according to the slope of the frequency ridge in the crossover interval, that is, for any two ridge components in the crossover interval, compare their slopes at the edge of the crossover interval. If the absolute value of the slope difference is less than the threshold ε6 (ε6 can be taken as l / 70), they are considered to belong to the same component and recombined and connected. Specifically, the schematic diagram is as follows Figure 5 As shown, for example, the pair is located in the crossover interval The instantaneous frequency ridge and Calculate the slope of several points on the edge of the crossing interval, such as Figure 5 As shown, the slopes of the four points are recorded as For these four slopes, if the absolute value of the slope difference between any two of them is less than ε6, such as Figure 5 middle, Then the reorganized new instantaneous frequency ridge is recorded as The instantaneous frequency ridges of all crossover intervals are reorganized to obtain a new instantaneous frequency ridge:

[0070] Step 6, fitting each reorganized instantaneous frequency ridge with Fourier series to obtain its time domain signal, and obtain each extracted component signal.

[0071] Each instantaneous frequency ridge obtained by recombining is expressed as For each instantaneous frequency ridge Apply Fourier series to fit the corresponding time domain signal Then we can get the time domain signal of all components extracted.

[0072] Step 7: extract the information of each component and calculate the rotor speed and blade length of each UAV.

[0073] That is, by performing synchronous compression transformation on the extracted time domain signal to obtain the time-frequency representation of each extracted component, the rotation speed and blade length of the drone rotor represented by each component can be calculated, thereby realizing the identification of multiple drones.

[0074] In the present invention, by integrating synchronous compression transformation into the signal separation strategy of ridge path recombination, the resolution of the separation process is improved, a multi-UAV rotor echo model is established, and the algorithm in the signal separation strategy of the present invention is used to extract the micro-Doppler characteristics of the multi-UAV rotors in the constructed echo model, thereby realizing the identification of multiple UAVs. The present invention integrates synchronous compression transformation into the "ridge path recombination" method for separating overlapping non-stationary signals, and improves the resolution of the ridge path recombination method. This improved signal separation strategy is used to extract the micro-motion characteristics of multiple UAV rotors, and the parameter information such as the rotor speed and blade length of each UAV is extracted to realize the identification of multiple UAVs.

[0075] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some of the technical features therein. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present invention.

[0076] The above are only some embodiments of the present invention. For those skilled in the art, several modifications and improvements can be made without departing from the creative concept of the present invention, which all belong to the protection scope of the present invention.

Claims

1. A UAV micro-Doppler extraction method based on synchronous compression transformation and ridge path reconstruction, characterized in that: The following steps are involved: Step 1: Modeling the multi-UAV rotor echoes; The one-dimensional time domain rotor echo signal s0(t) is expressed as: Wherein, t represents the sampling point, i represents the UAV ID, and k represents the UAV component ID. The UAV components include the UAV body and the UAV rotor blades. The radar cross section (RCS) of the components of the drone. represents the distance between the UAV component and the radar, λ represents the wavelength of the radar emission signal, and e represents the natural base; Step 2, perform synchronous compression transformation on the one-dimensional time domain rotor echo signal s0(t) to obtain its time-frequency representation, denoted as TF(t,f); Step 3, extract all instantaneous frequency ridges in the time-frequency representation TF(t,f): Find the time t corresponding to the maximum value in the time-frequency representation TF(t,f) max , frequency f max ; Based on time t max , traverse the time axis from the left and right directions respectively, for each moment t k , find the time t k The time-frequency representation of TF(t,f) corresponds to the frequency value of the maximum value point Based on the corresponding f at each moment k Get the frequency sequence [f1...,f k ,...,f l ] and as a frequency ridge IF i =[f1...,f k ,...,f l ], where l represents the number of sampling points; Frequency Ridge IF i The instantaneous frequency ridge is extracted by Fourier series smoothing fitting. And get the corresponding time domain signal s i (t), where represents the frequency value after fitting; Subtract the time domain signal s from the one-dimensional time domain rotor echo signal s0(t) i (t), and obtain the processed one-dimensional time domain rotor echo signal s′0(t); If ||s′0(t)||≥ε1, return to step 2; otherwise, end the extraction and obtain the instantaneous frequency ridge set based on all instantaneous frequency ridges Among them, ε1 represents the preset threshold; Step 4, find the intersection intervals of each frequency ridge and merge the intersection intervals with similar distances; Based on the preset threshold ε2, for the instantaneous frequency ridge set Any two instantaneous frequency ridges in and Traverse each time t k , if the current time t k The two frequencies corresponding to and If satisfied Then store the time number k in the vector findindex; For all the time numbers in the vector findindex, if the absolute difference between two adjacent numbers is less than the threshold ε3 and greater than the threshold ε4, then the time corresponding to these two numbers That is, the two instantaneous frequency ridges and The edge of a crossover interval of and The intersection interval of Traversing the instantaneous frequency ridge set All any two instantaneous frequency ridge combinations in , obtain their corresponding crossover intervals, and obtain a crossover interval set; For any two intersection intervals in the intersection interval set, if the Euclidean distance between their points is less than the threshold ε5, the two current intersection intervals are merged to obtain the intersection interval set after the intersection intervals are merged; Step 5, reorganize the ridges of the crossover interval to obtain a reorganized instantaneous frequency ridge set; For any two instantaneous frequency ridges in each crossover interval in the crossover interval set after the crossover interval is merged, and Compare the slopes at the edge of the crossover interval. If the absolute value of the slope difference is less than the threshold ε6, the two current instantaneous frequency ridges are and Reconnect the combination to obtain the reorganized instantaneous frequency ridge; Step 6, fitting each reorganized instantaneous frequency ridge with Fourier series to obtain its time domain signal, and obtain each extracted component signal.

2. The method according to claim 1, characterized in that When k represents the drone body, then Indicates the RCS of the main part of drone i, which is the preset value; When k represents a rotor blade of the UAV, then represents the RCS of the kth blade of the rotor of UAV i, and its expression is: Among them, parameter a i =L i / 2, parameter b i =W i / 2,L i represents the rotor blade length of UAV i, W i is the blade width of UAV i, θ represents the pitch angle of the kth blade of the rotor of UAV i relative to the radar, represents the azimuth angle of the kth blade of the rotor of UAV i relative to the radar.

3. The method according to claim 2, characterized in that When k represents the drone body, then The value of is set to 3 times the RCS of any UAV rotor blade of UAV i.

4. The method according to any one of claims 1 to 3, characterized in that The threshold ε1 is set to 0.2 of the norm of the one-dimensional time domain rotor echo signal s0(t).

5. The method according to any one of claims 1 to 3, characterized in that: The value of the threshold ε2 is set to the sampling frequency f s 30% of.

6. The method according to any one of claims 1 to 3, characterized in that: The values ​​of thresholds ε3, ε4, ε5 and ε6 are set as follows: ε3 is 1 / 200, ε4 is 1 / 500, and ε5 is 4(ε2) 2 , ε6 is taken as l / 70.

Images