A method and system for extracting radar target motion characteristics

Abstract

The present invention relates to a method and system for extracting radar target motion characteristics. The method includes preprocessing radar one-dimensional range profile data to obtain the phase signal of the target; using empirical mode decomposition to separate the signal of the phase signal to obtain a plurality of intrinsic mode functions; and performing motion parameter estimation on the intrinsic mode functions respectively based on the spectrum analysis method of Fourier analysis to obtain the radar target motion characteristics. The present invention can adapt to the actual situations of more diverse, high-speed motion, and motion coupling of radar targets, and the extracted target motion characteristics have higher accuracy.

Description

Technical Field

[0001] The present invention relates to the technical field of signal processing, and particularly to a method and system for extracting radar target motion characteristics. Background Art

[0002] A radar is an important system in the field of information acquisition and perception, having unique advantages such as all-day, all-weather, and long operating range. By using a radar system to obtain the echo signals of targets in environments such as aerospace, ground, and ocean, and then through signal processing tools such as noise suppression, imaging, and feature analysis, key information such as the geometric shape, motion state, target type, and working state of the targets can be obtained. Radar target feature extraction can provide conditions for target identification and interpretation.

[0003] The motion characteristics of radar targets are key elements for identifying target attributes, including specific forms such as translational motion, spin motion, and precession motion, which reflect the dynamic characteristics of the targets. Research shows that the accuracy of radar target motion feature extraction not only depends on the radar system itself, but is also closely related to radar target characteristics and fine signal processing. Existing target motion feature methods are mainly developed and designed for low-complexity situations such as single target motion, low-speed motion, and no coupled motion. The new characteristics of radar target motion such as diversification, high speed, and coupling bring many problems and challenges to radar target motion feature extraction.

[0004] First, the motion of radar targets is more diverse, making it difficult for existing motion feature extraction methods to take into account multiple composite motions. For example, a satellite that has lost its attitude control ability and fails, due to reasons such as space debris impact or insufficient power supply, the satellite exhibits rotational and tumbling motion forms, and the rotation speed will gradually decrease, making the motion non-uniform. For some space targets, due to the lateral force generated during impact or separation of the target, precession may also occur. Generally speaking, in addition to orbital motion, space targets will also generate various motion forms such as spin, tumbling, precession, and non-uniform rotation. Second, the high speed of radar target motion leads to a decline in the quality of radar echo signals. For example, the high-speed motion of supersonic moving radar targets such as missiles leads to an increase in radar echo errors and a decrease in signal-to-noise ratio at the same time, directly affecting subsequent motion feature extraction. Finally, the motion coupling caused by multiple scattering points of radar targets leads to easy misjudgment in motion feature extraction. For example, when multiple scattering points of a target perform the same spin motion, it is easy to double the rotation period of motion feature extraction, thereby causing the failure of motion feature extraction.

[0005] Facing the above bottleneck problems, there is an urgent need to develop new methods for extracting radar target motion characteristics. Summary of the Invention

[0006] In order to overcome the deficiencies of the prior art, the purpose of the present invention is to provide a method and system for extracting radar target motion characteristics.

[0007] To achieve the above object, the present invention provides the following solutions:

[0008] A method for extracting radar target motion characteristics, comprising:

[0009] Preprocessing the radar one-dimensional range profile data to obtain the phase signal of the target;

[0010] Using empirical mode decomposition to separate the signal of the phase signal to obtain a plurality of intrinsic mode functions;

[0011] Based on the spectrum analysis method of Fourier analysis, performing motion parameter estimation on the intrinsic mode functions respectively to obtain radar target motion characteristics.

[0012] Preferably, the preprocessing of the radar one-dimensional range profile data to obtain the phase signal of the target includes:

[0013] For the radar one-dimensional range profile X ∈ C M×N Calculating the absolute value, and then summing the absolute value column by column to obtain a summation parameter where X is the radar one-dimensional range profile, Y(n) is the summation parameter, C is the complex domain, M is the length of the one-dimensional range profile time series, N is the number of range cells of the one-dimensional range profile, m is the index of the one-dimensional range profile time series, n is the index of the one-dimensional range profile range cell, n = 1, 2,..., N, and Y ∈ R N ;

[0014] Setting a threshold value T 1 ∈ R, for Y ∈ R N , statistically analyzing the index set of the summation parameter and extracting the time series signal of the strong scattering unit according to the index set; where Γ is the index set, and the elements in Γ are the columns where the strong scattering units are located in the radar one-dimensional range profile X ∈ C M×N and k is the index of the one-dimensional range profile range cell;

[0015] Let y = arg(x) represent the phase data of the time series signal; where y = arg(x) is an operator for extracting complex phase, and y ∈ R M ;

[0016] Calculating the phase gradient D(m) = y(m) - y(m + 1) according to the phase data of the time series signal; where m = 1, 2,..., M - 1;

[0017] Calculating the wrapped phase difference according to the phase data of the time series signal: E(m) = arctan(sin(y(m)) / cos(y(m))), where m = 1, 2,..., M - 1;

[0018] Let \(z(1) = y(1)\); \(z(1)\) is the initialization data;

[0019] Sum the wrapped phases according to the formula \(z(m)=y(m)+E(m - 1)\) to obtain the unwrapped phase signal; where \(z(m)\) is the \(m\)-th element of the unwrapped phase signal.

[0020] Preferably, using empirical mode decomposition to separate the phase signal to obtain a plurality of intrinsic mode functions, including: obtaining the first intrinsic mode function, obtaining other intrinsic mode functions, and setting a stop condition;

[0021] The specific steps for obtaining the first intrinsic mode function are as follows:

[0022] Calculate the set of coordinates \(I\) of local maximum points and the set of coordinates \(J\) of local minimum points of the phase signal \(z\in R\) M ;

[0023] Use nearest neighbor interpolation to interpolate the local maximum points \(\{z(i)|i\in I\}\) to obtain the upper envelope signal \(a\in R\) M and use nearest neighbor interpolation to interpolate the local minimum points \(\{z(j)|j\in J\}\) to obtain the lower envelope signal \(b\in R\) M ;

[0024] Calculate the mean envelope signal \(m\) of the upper envelope signal and the lower envelope signal 1 \(=(a + b) / 2\);

[0025] Let \(h\) 1 \(=z - m\) 1 , and determine whether \(h\) 1 is an intrinsic mode function. If not, let \(h\) 1 replace \(z\), and jump to the step "Calculate the set of coordinates \(I\) of local maximum points and the set of coordinates \(J\) of local minimum points of the phase signal \(z\in R\)" M until the first intrinsic mode function is obtained;

[0026] The specific steps for obtaining other intrinsic mode functions are as follows:

[0027] Let \(r\) k \(=r\) k-1 \(-h\) k , and let \(r\) k replace \(z\), repeat the step "Obtain the first intrinsic mode function" until the \((k + 1)\)-th intrinsic mode function \(h\) k+1 is obtained, where \(k = 2,3,\cdots\);

[0028] The specific steps for setting the stop condition are as follows:

[0029] Set a stop threshold \(T\) 2 \(>0\), calculate the discriminant

[0030] If SD ≤ T 2 , stop the empirical mode decomposition, and record K = k - 1 at this time.

[0031] Preferably, for the spectral analysis method based on Fourier analysis, the motion parameters of the intrinsic mode functions are estimated respectively to obtain the radar target motion characteristics, including:

[0032] Perform Fourier transform and absolute value processing on the intrinsic mode function in sequence to obtain spectral data f k = |Fh k |; where F is the discrete Fourier matrix, and the elements of the discrete Fourier matrix are F ij = exp(-2πj 0 (i - 1)(j - 1) / M), where i, j ∈ {1, 2,..., M} and j 0 is the imaginary unit, f k is the spectral data; h k is the intrinsic mode function;

[0033] Take the frequency value corresponding to the maximum value of the spectral data as the estimation of the motion parameter to obtain the radar target motion characteristic w k = argmaxf k , k = 1, 2,..., K; where w k is the radar target motion characteristic.

[0034] A radar target motion characteristic extraction system includes:

[0035] A preprocessing unit for preprocessing the radar one-dimensional range profile data to obtain the phase signal of the target;

[0036] A signal separation unit for separating the phase signal by using empirical mode decomposition to obtain a plurality of intrinsic mode functions;

[0037] A parameter estimation unit for estimating the motion parameters of the intrinsic mode functions respectively based on the spectral analysis method of Fourier analysis to obtain the radar target motion characteristics.

[0038] Preferably, the preprocessing unit includes:

[0039] A summation module for calculating the absolute value of the radar one-dimensional range profile X ∈ C M×N and then summing the absolute values column by column to obtain the summation parameter Wherein, X is the one-dimensional range profile of the radar, Y(n) is the summation parameter, C is the complex domain, M is the length of the one-dimensional range profile time series, N is the number of range cells of the one-dimensional range profile, m is the index of the one-dimensional range profile time series, n is the index of the one-dimensional range profile range cell, n = 1, 2, …, N, and Y ∈ R N ;

[0040] An extraction module, configured to set a threshold value T 1 ∈R, for Y ∈ R N , and count the index set of the summation parameter and extract the time series signal of the strong scattering unit according to the index set; wherein Γ is the index set, and the elements in Γ are the columns where the strong scattering units in the radar one-dimensional range profile X ∈ C M×N are located, and k is the index of the one-dimensional range profile range cell;

[0041] A phase representation module, configured to make y = arg(x) represent the phase data of the time series signal; wherein y = arg(x) is an operator for extracting the complex phase, and y ∈ R M ;

[0042] A first calculation module, configured to calculate the phase gradient D(m) = y(m) - y(m + 1) according to the phase data of the time series signal; wherein m = 1, 2, …, M - 1;

[0043] A second calculation module, configured to calculate the wrapped phase difference according to the phase data of the time series signal: E(m) = arctan(sin(y(m)) / cos(y(m))), wherein m = 1, 2, …, M - 1;

[0044] An initialization module, configured to make z(1) = y(1); z(1) is the initialization data;

[0045] A third calculation module, configured to sum the wrapped phases according to the formula z(m) = y(m) + E(m - 1) to obtain the unwrapped phase signal; wherein z(m) is the m-th element of the unwrapped phase signal.

[0046] Preferably, the signal separation unit includes:

[0047] A first acquisition module, configured to acquire the first intrinsic mode function;

[0048] A second acquisition module, configured to acquire other intrinsic mode functions;

[0049] A stop module, configured to set a stop condition;

[0050] The first acquisition module includes:

[0051] An extreme value calculation sub-module, configured to calculate the phase signal z ∈ RM The local maximum point coordinate set I and the local minimum point coordinate set J;

[0052] An interpolation sub-module, configured to interpolate the local maximum points {z(i)|i∈I} using nearest neighbor interpolation to obtain an upper envelope signal a∈R M and interpolate the local minimum points {z(j)|j∈J} using nearest neighbor interpolation to obtain a lower envelope signal b∈R M ;

[0053] A mean calculation sub-module, configured to calculate the mean envelope signal m of the upper envelope signal and the lower envelope signal 1 =(a + b) / 2;

[0054] A first acquisition sub-module, configured to set h 1 =z - m 1 , and determine whether h 1 is an intrinsic mode function. If not, let h 1 replace z, and jump to the step "calculate the local maximum point coordinate set I and the local minimum point coordinate set J of the phase signal z∈R M " until the first intrinsic mode function is obtained;

[0055] The second acquisition module includes:

[0056] A second acquisition sub-module, configured to set r k =r k-1 -h k , and let r k replace z, and repeat the step "obtain the first intrinsic mode function" until the (k + 1)-th intrinsic mode function h k+1 is obtained, where k = 2, 3,...;

[0057] The stop module includes:

[0058] A discriminant calculation sub-module, configured to set a stop threshold T 2 >0, and calculate the discriminant

[0059] A judgment sub-module, configured to stop continuing empirical mode decomposition when SD≤T 2 , and record K = k - 1 at this time.

[0060] Preferably, the parameter estimation unit includes:

[0061] A Fourier processing module, configured to perform Fourier transform and absolute value processing on the intrinsic mode function in sequence to obtain spectrum data f k =|Fh k |; where F is a discrete Fourier matrix, and the elements of the discrete Fourier matrix are Fij = exp(-2πj 0 (i - 1)(j - 1) / M), where i, j ∈ {1, 2, …, M} and j 0 is the imaginary unit, f k is the said spectral data; h k is the said intrinsic mode function;

[0062] An estimation module, configured to take the frequency value corresponding to the maximum value of the said spectral data as the estimation of the motion parameter, and obtain the radar target motion feature w k = argmax f k , k = 1, 2, …, K; where w k is the said radar target motion feature.

[0063] According to the specific embodiments provided by the present invention, the following technical effects are disclosed by the present invention:

[0064] The present invention provides a method and system for extracting radar target motion features. The method includes preprocessing radar one-dimensional range profile data to obtain the phase signal of the target; using empirical mode decomposition to separate the signal of the phase signal to obtain a plurality of intrinsic mode functions; based on the spectral analysis method of Fourier analysis, respectively performing motion parameter estimation on the said intrinsic mode functions to obtain radar target motion features. The present invention can adapt to the actual situations of more diverse, high-speed moving, and motion-coupled radar targets, and the extracted target motion features have higher accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0065] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required to be used in the embodiments. Obviously, the drawings in the following description are only some embodiments of the present invention. For those of ordinary skill in the art, other drawings can be obtained based on these drawings without creative efforts.

[0066] Figure 1 is the flowchart of the method in the embodiment provided by the present invention;

[0067] Figure 2 is the flowchart of radar target motion feature extraction in the embodiment provided by the present invention;

[0068] Figure 3 is the flowchart of obtaining the intrinsic mode function in the embodiment provided by the present invention;

[0069] Figure 4 is the radar target structure diagram in the embodiment provided by the present invention;

[0070] Figure 5The first phase signal diagram before and after phase unwrapping in the embodiment provided by the present invention;

[0071] Figure 6 The second phase signal diagram before and after phase unwrapping in the embodiment provided by the present invention;

[0072] Figure 7 The empirical mode decomposition result in the embodiment provided by the present invention. Detailed implementation manners

[0073] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0074] Referring to "embodiments" herein means that the specific features, structures, or characteristics described in conjunction with the embodiments may be included in at least one embodiment of the present application. The phrase appears in various places in the specification does not necessarily refer to the same embodiment, nor is it an independent or alternative embodiment mutually exclusive with other embodiments. Those skilled in the art explicitly and implicitly understand that the embodiments described herein may be combined with other embodiments.

[0075] The terms "first", "second", "third", and "fourth", etc. in the specification and claims of the present application and the accompanying drawings are used to distinguish different objects, rather than to describe a specific order. In addition, the terms "comprising" and "having" and any variations thereof are intended to cover non-exclusive inclusion. For example, a series of steps, processes, methods, etc. included do not limit to the listed steps, but optionally further include steps not listed, or optionally further include other step elements inherent to these processes, methods, products, or devices.

[0076] The purpose of the present invention is to provide a method and system for extracting radar target motion features, which can adapt to the actual situations of more diverse, high-speed, and motion-coupled radar target motions, and the extracted target motion features have higher accuracy.

[0077] To make the above objects, features, and advantages of the present invention more obvious and understandable, the present invention will be further described in detail below in conjunction with the accompanying drawings and specific implementation manners.

[0078] Figure 1 and Figure 2 respectively are the method flow chart and the radar target motion feature extraction flow chart in the embodiment provided by the present invention, as Figure 1 and Figure 2As shown, the present invention provides a method for extracting radar target motion characteristics, including:

[0079] Step 100: Preprocess the radar one-dimensional range profile data to obtain the phase signal of the target;

[0080] Step 200: Use empirical mode decomposition to separate the signal of the phase signal to obtain multiple intrinsic mode functions;

[0081] Step 300: Based on the spectral analysis method of Fourier analysis, perform motion parameter estimation on the intrinsic mode functions respectively to obtain the radar target motion characteristics.

[0082] Preferably, the step 100 includes:

[0083] For the radar one-dimensional range profile X ∈ C M×N Calculate the absolute value, and then sum the absolute values column by column to obtain the summation parameter where X is the radar one-dimensional range profile, Y(n) is the summation parameter, C is the complex domain, M is the length of the one-dimensional range profile time series, N is the number of range cells of the one-dimensional range profile, m is the index of the one-dimensional range profile time series, n is the index of the one-dimensional range profile range cell, n = 1, 2,..., N, and Y ∈ R N ;

[0084] Set the threshold Set the threshold T 1 ∈ R, for Y ∈ R N , statistically analyze the index set of the summation parameter and extract the time series signal of the strong scattering unit according to the index set; where Γ is the index set, and the elements in Γ are the columns where the strong scattering units are located in the radar one-dimensional range profile X ∈ C M×N , and k is the index of the one-dimensional range profile range cell;

[0085] Let y = arg(x) represent the phase data of the time series signal; where y = arg(x) is the operator for extracting the complex phase, and y ∈ R M ;

[0086] Calculate the phase gradient D(m) = y(m) - y(m + 1) according to the phase data of the time series signal; where m = 1, 2,..., M - 1;

[0087] Calculate the wrapped phase difference according to the phase data of the time series signal: E(m) = arctan(sin(y(m)) / cos(y(m))), where m = 1, 2,..., M - 1;

[0088] Let z(1) = y(1); z(1) is the initialization data;

[0089] Sum the wrapped phase according to the formula z(m) = y(m) + E(m - 1) to obtain the unwrapped phase signal; where z(m) is the m-th element of the unwrapped phase signal.

[0090] The first step in this embodiment is the "data preprocessing" step. Let X ∈ C M×N be the one-dimensional range image of the radar after range alignment. Data preprocessing mainly realizes obtaining the phase signal of the target from the radar one-dimensional range image data for subsequent empirical mode decomposition and motion feature extraction.

[0091] Specifically, it includes:

[0092] S1.1: Extract strong scattering units

[0093] For the radar one-dimensional range image X ∈ C M×N take the absolute value and then sum column by column, that is:

[0094]

[0095] where n = 1, 2, …, N, Y ∈ R N . Set the threshold T 1 ∈ R. For Y ∈ R N , count its index set:

[0096]

[0097] Then the elements in Γ are the columns where the strong scattering units in the radar one-dimensional range image X ∈ C M×N are located. For the convenience of discussion, assume that there is only one element n 0 in Γ. At this time, let x(m) = X(m, n 0 ) represent the time series signal corresponding to the strong scattering unit, where m = 1, 2, …, M. If Γ is an empty set, adjust the threshold T; if the number of elements in Γ is greater than 1, then repeat the same processing steps when there is only 1 element for each element respectively.

[0098] S1.2: Phase unwrapping

[0099] For the time series signal x ∈ C M of the extracted strong scattering unit, let y = arg(x) represent the phase signal of this signal, where arg(·) is the operator for extracting the complex phase, and y ∈ R M . Since directly extracting the phase, its phase value will be in the range of (-π, π], it is necessary to use phase unwrapping to solve the complete phase signal. The main steps of phase unwrapping are as follows:

[0100] 1) Calculate the phase gradient: D(m) = y(m) - y(m + 1), where m = 1, 2, …, M - 1;

[0101] 2) Calculate the wrapped phase difference: E(m) = arctan(sin(y(m)) / cos(y(m))), where m = 1, 2, …, M - 1;

[0102] 3) Initialize: Let z(1) = y(1);

[0103] 4) Sum the wrapped phases to obtain the unwrapped phase signal: z(m) = y(m) + E(m - 1), where m = 1, 2, …, M - 1;

[0104] The phase signal z ∈ R after phase unwrapping M can be used for subsequent signal separation and motion feature extraction.

[0105] Preferably, the step 200 includes: obtaining the first intrinsic mode function, obtaining other intrinsic mode functions, and setting a stop condition;

[0106] The specific steps for obtaining the first intrinsic mode function are:

[0107] Calculate the local maximum point coordinate set I and the local minimum point coordinate set J of the phase signal z ∈ R M ;

[0108] Interpolate the local maximum points {z(i)|i ∈ I} using nearest-neighbor interpolation to obtain the upper envelope signal a ∈ R M , and interpolate the local minimum points {z(j)|j ∈ J} using nearest-neighbor interpolation to obtain the lower envelope signal b ∈ R M ;

[0109] Calculate the mean envelope signal m 1 =(a + b) / 2;

[0110] Let h 1 = z - m 1 , and determine whether h 1 is an intrinsic mode function. If not, let h 1 replace z, and jump to the step "Calculate the local maximum point coordinate set I and the local minimum point coordinate set J of the phase signal z ∈ R M " until the first intrinsic mode function is obtained;

[0111] The specific steps for obtaining other intrinsic mode functions are:

[0112] Let r k = r k-1 - h k , and let r k replace z, and repeat the step "Obtain the first intrinsic mode function" until the (k + 1)-th intrinsic mode function h is obtainedk+1 , where k = 2, 3, …;

[0113] The specific steps for setting the stop condition are as follows:

[0114] Set a stop threshold T 2 > 0, and calculate the discriminant

[0115] If SD ≤ T 2 , then stop the empirical mode decomposition, and record K = k - 1 at this time.

[0116] The second step in this embodiment is "empirical mode decomposition". For the phase signal z ∈ R M , first, use empirical mode decomposition to separate the signals corresponding to different motion forms. The specific steps are as follows:

[0117] S2.1: Obtain the first intrinsic mode function

[0118] The flowchart for obtaining the first intrinsic mode function is shown in Figure 3 , and the specific steps are as follows:

[0119] 1) Calculate the set of local maximum point coordinates I and the set of local minimum point coordinates J of the signal z ∈ R M ;

[0120] 2) Use nearest-neighbor interpolation to interpolate the local maximum points {z(i)|i ∈ I} to obtain the upper envelope signal a ∈ R M ; at the same time, use nearest-neighbor interpolation to interpolate the local minimum points {z(j)|j ∈ J} to obtain the lower envelope signal b ∈ R M ;

[0121] 3) Calculate the mean envelope signal m 1 = (a + b) / 2;

[0122] 4) Let h 1 = z - m 1 ; if h 1 is an intrinsic mode function, the algorithm stops; otherwise, let h 1 replace z, and repeat the above steps 1)-4) until the first intrinsic mode function is obtained.

[0123] The necessary condition that the intrinsic mode function in step 4) should satisfy is: the number of its extreme points p and the number of zero-crossing points q satisfy q - 1 ≤ p ≤ q + 1.

[0124] S2.2 Obtain other intrinsic mode functions

[0125] Let the first intrinsic mode function of the signal z ∈ R M be h1 , then let r 1 = z - h 1 , and let r 1 replace z, and repeat the above step S2.1 until the second intrinsic mode function h 2 is obtained.

[0126] Let r 2 = r 1 - h 2 , and let r 2 replace z, and repeat the above step S2.1 until the third intrinsic mode function h 3 is obtained.

[0127] And so on. Without loss of generality, let r k = r k-1 - h k , and let r k replace z, and repeat the above step S2.1 until the (k + 1)-th intrinsic mode function h k+1 is obtained, where k = 2, 3,....

[0128] S2.3 Stop condition

[0129] Set a stop threshold T 2 > 0, and calculate the discriminant:

[0130]

[0131] If SD ≤ T 2 , then stop the empirical mode decomposition, and at this time, record K = k - 1, that is, the phase signal z ∈ R M is decomposed into K intrinsic mode functions in total.

[0132] Preferably, the step 300 includes:

[0133] Perform Fourier transform and absolute value processing on the intrinsic mode functions in sequence to obtain spectral data f k = |Fh k |; where F is the discrete Fourier matrix, and the elements of the discrete Fourier matrix are F ij = exp(-2πj 0 (i - 1)(j - 1) / M), where i, j ∈ {1, 2,..., M} and j 0 is the imaginary unit, f k is the spectral data; h k is the intrinsic mode function;

[0134] Take the frequency value corresponding to the maximum value of the spectral data as the estimation of the motion parameter to obtain the radar target motion feature w k = argmax fk , k = 1, 2, …, K; where w k is the motion feature of the radar target.

[0135] Specifically, the third step in this embodiment is motion feature extraction. For the intrinsic mode functions obtained after empirical mode decomposition, motion parameter estimation is performed respectively, and then motion feature extraction is realized. Here, a spectral analysis method based on Fourier analysis is mainly used to estimate the motion parameters. The specific steps are as follows:

[0136] S3.1: Spectral analysis

[0137] Considering the intrinsic mode function h k ∈R M , perform Fourier transform on it and take the absolute value, that is:

[0138] f k = |Fh k |.

[0139] Among them, k = 1, 2, …, K, F is the discrete Fourier matrix, and its element is F ij = exp(-2πj 0 (i - 1)(j - 1) / M), where i, j ∈ {1, 2, …, M} and j 0 is the imaginary unit.

[0140] S3.2: Motion parameter estimation

[0141] Take the frequency value corresponding to the maximum value of f k as the estimation of the motion parameter:

[0142] w k = argmax f k , k = 1, 2, …, K.

[0143] At this time, the obtained motion feature of the radar target is {w 1 , w 2 , …, w K}.

[0144] This embodiment also provides another specific implementation manner, in which the radar parameters are set as follows: carrier frequency 9.5 GHz, bandwidth 1 GHz, pulse repetition frequency 300 Hz. Among them, the target is four ideal scattering points for simulation, as Figure 4 shown. Under the radar parameters, target echo data is simulated and pulse compression and range alignment are performed to generate a one-dimensional range profile X ∈ C 5000×1024 as the input signal of this method. The specific steps are as follows:

[0145] S1: Data preprocessing

[0146] Taking the radar one-dimensional range image \(X\in\mathbb{C}\) after distance alignment as the input, the data preprocessing mainly realizes obtaining the phase signal of the target from the radar one-dimensional range image data for subsequent empirical mode decomposition and motion feature extraction. 5000×1024

[0147] S1.1: Extracting strong scattering cells

[0148] Taking the absolute value of the radar one-dimensional range image \(X\in\mathbb{C}\) 5000×1024 and then summing up column by column, that is

[0149]

[0150] where \(n = 1,2,\cdots,1024\), \(Y\in\mathbb{R}\) 1024 Set the threshold \(T\) 1 \(= 0.3\), for \(Y\in\mathbb{R}\) 1024 and count its index set:

[0151]

[0152] Then the columns where the strong scattering cells in the radar one-dimensional range image \(X\in\mathbb{C}\) 5000×1024 are located are the 503rd column and the 518th column. For the convenience of discussion, only the 503rd strong scattering cell is considered in this example, and a similar treatment is adopted for the 518th strong scattering cell. At this time, let \(x(m)=X(m,503)\) represent the time series signal corresponding to the strong scattering cell, where \(m = 1,2,\cdots,5000\).

[0153] S1.2: Phase unwrapping

[0154] For the time series signal \(x\in\mathbb{C}\) of the extracted strong scattering cell 5000 , let \(y = \arg(x)\) represent the phase signal of this signal, where \(\arg(\cdot)\) is the operator for extracting the complex phase, and \(y\in\mathbb{R}\) 5000 . Since the directly extracted phase value will be between \((-\pi,\pi]\), it is necessary to use phase unwrapping to solve the complete phase signal. The main steps of phase unwrapping are as follows:

[0155] 1) Calculate the phase gradient: \(D(m)=y(m)-y(m + 1)\), where \(m = 1,2,\cdots,4999\);

[0156] 2) Calculate the wrapped phase difference: \(E(m)=\arctan(\sin(y(m)) / \cos(y(m)))\), where \(m = 1,2,\cdots,4999\);

[0157] 3) Initialize: Let \(z(1)=y(1)\);

[0158] 4) Sum the wrapped phases to obtain the unwrapped phase signal: z(m) = y(m) + E(m - 1), where m = 1, 2, …, 4999;

[0159] The phase signal z ∈ R after phase unwrapping 5000 As Figure 5 and Figure 6 shown, it can be used for subsequent signal separation and motion feature extraction.

[0160] S2: Empirical Mode Decomposition

[0161] For the phase signal z ∈ R 5000 , first use empirical mode decomposition to separate the signals corresponding to different motion forms. The specific steps are as follows:

[0162] S2.1: Obtain the first intrinsic mode function

[0163] According to Figure 3 , obtain the first intrinsic mode function;

[0164] S2.2 Obtain other intrinsic mode functions

[0165] Let the first intrinsic mode function of the signal z ∈ R M be h 1 , then let r 1 = z - h 1 , and let r 1 replace z, and repeat the above step S2.1 until the second intrinsic mode function h 2 is obtained.

[0166] And so on. Let's assume r k = r k-1 - h k , and let r k replace z, and repeat the above step S2.1 until the (k + 1)-th intrinsic mode function h k+1 is obtained, where k = 2, 3, ….

[0167] S2.3 Stopping condition

[0168] Set the stopping threshold T 2 = 0.2, and calculate the discriminant:

[0169]

[0170] If SD ≤ T 2 , then stop further empirical mode decomposition. At this time, record K = k - 1 = 2, that is, the phase signal z ∈ R 5000 is decomposed into 2 intrinsic mode functions in total. The result is shown in Figure 7 .

[0171] S3: Motion Feature Extraction

[0172] For the intrinsic mode functions after empirical mode decomposition, motion parameter estimation is performed respectively, and then motion feature extraction is realized. Here, the spectral analysis method based on Fourier analysis is mainly used to realize the estimation of motion parameters. Considering the above two intrinsic mode functions, Fourier transform is performed on them, and the absolute value is taken, and then the frequency value corresponding to the maximum value is taken as the estimation of the motion parameter. The results are shown in the following table. The results show that this method can extract two motion features, namely the target spin frequency and the precession frequency, and the estimated values of the motion parameters are very close to the true values, that is, the accuracy is relatively high.

[0173] Table 1

[0174] Estimated value True value Spin frequency 9.75 Hz 10 Hz Precession frequency 0.49 Hz 0.5 Hz

[0175] This embodiment also provides a radar target motion feature extraction system, including:

[0176] A preprocessing unit for preprocessing the radar one-dimensional range profile data to obtain the phase signal of the target;

[0177] A signal separation unit for separating the phase signal by using empirical mode decomposition to obtain a plurality of intrinsic mode functions;

[0178] A parameter estimation unit for respectively performing motion parameter estimation on the intrinsic mode functions based on the spectral analysis method of Fourier analysis to obtain the radar target motion features.

[0179] Preferably, the preprocessing unit includes:

[0180] A summation module for calculating the absolute value of the radar one-dimensional range profile X ∈ C M×N and then summing the absolute values column by column to obtain a summation parameter where X is the radar one-dimensional range profile, Y(n) is the summation parameter, C is the complex domain, M is the length of the one-dimensional range profile time series, N is the number of range cells of the one-dimensional range profile, m is the index of the one-dimensional range profile time series, n is the index of the range cell of the one-dimensional range profile, n = 1, 2,..., N, and Y ∈ R N ;

[0181] An extraction module for setting a threshold T 1 ∈ R, and for Y ∈ R N , statistically analyzing the index set of the summation parameter and extracting the time series signal of the strong scattering unit according to the index set; where Γ is the index set, and the elements in Γ are the columns where the strong scattering units are located in the radar one-dimensional range profile X ∈ C M×N and k is the index of the range cell of the one-dimensional range profile;

[0182] A phase representation module, configured to make y = arg(x) represent the phase data of the timing signal; where y = arg(x) is an operator for extracting complex phases, and y ∈ R M ;

[0183] A first calculation module, configured to calculate a phase gradient D(m) = y(m) - y(m + 1) according to the phase data of the timing signal; where m = 1, 2, …, M - 1;

[0184] A second calculation module, configured to calculate a wrapped phase difference according to the phase data of the timing signal: E(m) = arctan(sin(y(m)) / cos(y(m))), where m = 1, 2, …, M - 1;

[0185] An initialization module, configured to make z(1) = y(1); z(1) is initialization data;

[0186] A third calculation module, configured to sum the wrapped phases according to the formula z(m) = y(m) + E(m - 1) to obtain the unwrapped phase signal; where z(m) is the m-th element of the unwrapped phase signal.

[0187] Preferably, the signal separation unit includes:

[0188] A first acquisition module, configured to acquire the first intrinsic mode function;

[0189] A second acquisition module, configured to acquire other intrinsic mode functions;

[0190] A stop module, configured to set a stop condition;

[0191] The first acquisition module includes:

[0192] An extreme value calculation sub-module, configured to calculate a set of coordinates I of local maximum points and a set of coordinates J of local minimum points of the phase signal z ∈ R M ;

[0193] An interpolation sub-module, configured to perform interpolation on the local maximum points {z(i)|i ∈ I} using nearest neighbor interpolation to obtain an upper envelope signal a ∈ R M and perform interpolation on the local minimum points {z(j)|j ∈ J} using nearest neighbor interpolation to obtain a lower envelope signal b ∈ R M ;

[0194] An average value calculation sub-module, configured to calculate an average envelope signal m 1 =(a + b) / 2;

[0195] A first acquisition sub-module, configured to make h 1 = z - m1 , determine h 1 whether it is an intrinsic mode function. If not, let h 1 replace z and jump to the step "calculate the local maximum point coordinate set I and local minimum point coordinate set J of the phase signal z ∈ R M ", until the first intrinsic mode function is obtained;

[0196] The second acquisition module includes:

[0197] The second acquisition sub-module is used to set r k = r k-1 -h k , and let r k replace z, and repeat the step "obtain the first intrinsic mode function" until the (k + 1)-th intrinsic mode function h k+1 is obtained, where k = 2, 3,...;

[0198] The stop module includes:

[0199] The discriminant calculation sub-module is used to set the stop threshold T 2 > 0, and calculate the discriminant

[0200] The judgment sub-module is used to stop continuing the empirical mode decomposition if SD ≤ T 2 At this time, record K = k - 1.

[0201] Preferably, the parameter estimation unit includes:

[0202] The Fourier processing module is used to perform Fourier transform and absolute value processing on the intrinsic mode function in sequence to obtain the spectrum data f k = |Fh k |; where F is the discrete Fourier matrix, and the elements of the discrete Fourier matrix are F ij = exp(-2πj 0 (i - 1)(j - 1) / M), where i, j ∈ {1, 2,..., M} and j 0 is the imaginary unit, f k is the spectrum data; h k is the intrinsic mode function;

[0203] The estimation module is used to take the frequency value corresponding to the maximum value of the spectrum data as the estimation of the motion parameter to obtain the radar target motion feature w k = argmaxf k , k = 1, 2,..., K; where w k is the radar target motion feature.

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

[0205] (1) The present invention can extract two motion characteristics, namely the target spin frequency and the precession frequency, and the estimated values of the motion parameters are very close to the true values, that is, the accuracy is relatively high.

[0206] (2) The present invention can adapt to more diverse, high-speed and motion-coupled actual situations of radar targets, and the extracted target motion characteristics have higher accuracy.

[0207] The various embodiments in this specification are described in a progressive manner. Each embodiment focuses on the differences from other embodiments. For the same and similar parts among the various embodiments, reference can be made to each other. For the system disclosed in the embodiments, since it corresponds to the method disclosed in the embodiments, the description is relatively simple, and reference can be made to the description in the method part for related parts.

[0208] Specific examples are used in this article to elaborate on the principle and implementation manner of the present invention. The description of the above embodiments is only used to help understand the method of the present invention and its core idea; at the same time, for those of ordinary skill in the art, according to the idea of the present invention, there will be changes in the specific implementation manner and application scope. In summary, the content of this specification should not be construed as a limitation on the present invention.

Claims

1. A method for extracting radar target motion characteristics, characterized in that, comprising: Preprocessing the radar one-dimensional range image data to obtain the phase signal of the target; Using empirical mode decomposition to separate the signal of the phase signal to obtain a plurality of intrinsic mode functions; Based on the spectrum analysis method of Fourier analysis, performing motion parameter estimation on the intrinsic mode functions respectively to obtain radar target motion characteristics; The using empirical mode decomposition to separate the signal of the phase signal to obtain a plurality of intrinsic mode functions includes: obtaining the first intrinsic mode function, obtaining other intrinsic mode functions, and setting a stop condition; The specific steps for obtaining the first intrinsic mode function are: Calculate the set of coordinates I of local maximum points and the set of coordinates J of local minimum points of the phase signal z ∈ ℝ M ; Interpolate the local maximum points \(\{z(i)|i\in I\}\) using nearest-neighbor interpolation to obtain the upper envelope signal \(a\in\mathbb{R}\) M and interpolate the local minimum points \(\{z(j)|j\in J\}\) using nearest-neighbor interpolation to obtain the lower envelope signal \(b\in\mathbb{R}\) M ; Calculate the mean envelope signal m of the upper envelope signal and the lower envelope signal 1 =(a + b) / 2; Let h 1 = z - m 1 , and determine whether h 1 is an intrinsic mode function. If not, let h 1 replace z, and jump to the step "calculate the local maximum point coordinate set I and local minimum point coordinate set J of the phase signal z ∈ R M " until the first intrinsic mode function is obtained; The specific steps for obtaining other intrinsic mode functions are: Let r k = r k-1 - h k , and let r k replace z, and repeat the step of "obtaining the first intrinsic mode function" until the (k + 1)-th intrinsic mode function h k+1 is obtained, where k = 2, 3,...; The specific steps for setting the stop condition are: Set the stop threshold T 2 > 0, calculate the discriminant If SD ≤ T 2 , stop continuing with empirical mode decomposition, and at this time record K = k - 1; The based on the spectrum analysis method of Fourier analysis, performing motion parameter estimation on the intrinsic mode functions respectively to obtain radar target motion characteristics includes: The Fourier transform and absolute value processing are successively performed on the intrinsic mode function to obtain spectral data f k = |Fh k |; where F is a discrete Fourier matrix, and the elements of the discrete Fourier matrix are F ij = exp(-2πj 0 (i - 1)(j - 1) / M), where i, j ∈ {1, 2, …, M}, and j 0 is the imaginary unit, f k is the spectral data; h k is the intrinsic mode function; Take the frequency value corresponding to the maximum value of the said spectrum data as the estimation of the motion parameter, and obtain the radar target motion feature w k = argmax f k , k = 1, 2, …, K; where w k is the said radar target motion feature.

2. The radar target motion characteristic extraction method according to claim 1, characterized in that, The preprocessing the radar one-dimensional range image data to obtain the phase signal of the target includes: For the radar one-dimensional range profile \(X\in\mathbb{C}\) M×N Calculate the absolute value, and then sum the absolute values column by column to obtain the summation parameter where \(X\) is the radar one-dimensional range profile, \(Y(n)\) is the summation parameter, \(\mathbb{C}\) is the complex domain, \(M\) is the length of the one-dimensional range profile time series, \(N\) is the number of range cells of the one-dimensional range profile, \(m\) is the index of the one-dimensional range profile time series, \(n\) is the index of the one-dimensional range profile range cell, \(n = 1,2,\cdots,N\), and \(Y\in\mathbb{R}\) N ; Set the threshold T 1 ∈R, for Y∈R N , count the index set of the summation parameter and extract the timing signal of the strong scattering unit according to the index set; where Γ is the index set, and the elements in Γ are the one-dimensional radar range images X∈C M×N the column where the strong scattering unit is located in, k is the one-dimensional range image range cell index; Let \(y = \arg(x)\) represent the phase data of the timing signal; where \(y=\arg(x)\) is an operator for extracting the complex phase, and \(y\in R\). M ; Calculating the phase gradient D(m) = y(m) - y(m + 1) according to the phase data of the timing signal; where m = 1, 2,..., M - 1; Calculating the wrapped phase difference according to the phase data of the timing signal: E(m) = arctan(sin(y(m)) / cos(y(m))), where m = 1, 2,..., M - 1; Let z(1) = y(1); z(1) is the initialization data; Summing the wrapped phases according to the formula z(m) = y(m) + E(m - 1) to obtain the unwrapped phase signal; where z(m) is the m-th element of the unwrapped phase signal.

3. A radar target motion characteristic extraction system, characterized in that, comprising: A preprocessing unit for preprocessing the radar one-dimensional range image data to obtain the phase signal of the target; A signal separation unit for using empirical mode decomposition to separate the signal of the phase signal to obtain a plurality of intrinsic mode functions; A parameter estimation unit for performing motion parameter estimation on the intrinsic mode functions respectively based on the spectrum analysis method of Fourier analysis to obtain radar target motion characteristics; The signal separation unit includes: A first acquisition module for acquiring the first intrinsic mode function; A second acquisition module for acquiring other intrinsic mode functions; A stop module for setting a stop condition; The first acquisition module includes: An extreme value calculation sub-module, which is used to calculate the local maximum point coordinate set I and the local minimum point coordinate set J of the phase signal z ∈ R M ; An interpolation sub-module, which is used to interpolate the local maximum points {z(i)|i∈I} by nearest neighbor interpolation to obtain an upper envelope signal a∈R M and interpolate the local minimum points {z(j)|j∈J} by nearest neighbor interpolation to obtain a lower envelope signal b∈R M ; The mean calculation sub-module is used to calculate the mean envelope signal m of the upper envelope signal and the lower envelope signal 1 =(a + b) / 2; The first acquisition sub-module is used to make h 1 = z - m 1 , and determine whether h 1 is an intrinsic mode function. If not, then make h 1 replace z, and jump to the step "calculate the local maximum point coordinate set I and the local minimum point coordinate set J of the phase signal z ∈ R M ", until the first intrinsic mode function is obtained; The second acquisition module includes: The second acquisition sub-module is used to set r k = r k-1 - h k and let r k replace z, and repeat the step of "acquiring the first intrinsic mode function" until the (k + 1)-th intrinsic mode function h k+1 is obtained, where k = 2, 3,...; The stop module includes: Discriminant calculation sub-module, used to set the stop threshold T 2 > 0, calculate the discriminant A judgment sub-module, which is used to stop continuing the empirical mode decomposition if SD ≤ T 2 At this time, record K = k - 1 The parameter estimation unit includes: A Fourier processing module, configured to perform Fourier transform and absolute value processing on the intrinsic mode functions in sequence to obtain spectral data f k = |Fh k |; where F is a discrete Fourier matrix, and the elements of the discrete Fourier matrix are F ij = exp(-2πj 0 (i - 1)(j - 1) / M), where i, j ∈ {1, 2, …, M}, and j 0 is the imaginary unit, f k is the spectral data; h k is the intrinsic mode function; An estimation module, configured to use the frequency value corresponding to the maximum value of the spectrum data as an estimation of the motion parameter, so as to obtain a radar target motion feature w k = argmaxf k , k = 1, 2, …, K; where w k is the radar target motion feature.

4. The radar target motion characteristic extraction system according to claim 3, characterized in that, The preprocessing unit includes: A summation module for summing the radar one-dimensional range profile \(X\in\mathbb{C}\) M×N Calculate the absolute value, and then sum the absolute values column by column to obtain a summation parameter where \(X\) is the radar one-dimensional range profile, \(Y(n)\) is the summation parameter, \(\mathbb{C}\) is the complex domain, \(M\) is the length of the one-dimensional range profile time series, \(N\) is the number of range cells of the one-dimensional range profile, \(m\) is the index of the one-dimensional range profile time series, \(n\) is the index of the one-dimensional range profile range cell, \(n = 1,2,\cdots,N\), and \(Y\in\mathbb{R}\) N ; Extraction module, used to set threshold T 1 ∈R, for Y∈R N , statistically calculate the index set of the summation parameter and extract the timing signal of the strong scattering unit according to the index set; where Γ is the index set, and the elements in Γ are the one-dimensional range profiles of the radar X∈C M×N in the column where the strong scattering unit is located, and k is the range cell index of the one-dimensional range profile; A phase representation module, configured to use y = arg(x) to represent the phase data of the timing signal; where y = arg(x) is an operator for extracting complex phases, and y ∈ R M ; A first calculation module for calculating the phase gradient D(m) = y(m) - y(m + 1) according to the phase data of the timing signal; where m = 1, 2,..., M - 1; A second calculation module, configured to calculate the wrapped phase difference according to the phase data of the timing signal: E(m) = arctan(sin(y(m)) / cos(y(m))), where m = 1, 2, …, M - 1; An initialization module, configured to set z(1) = y(1); z(1) is the initialization data; A third calculation module, configured to sum the wrapped phases according to the formula z(m) = y(m) + E(m - 1) to obtain the unwrapped phase signal; where z(m) is the m-th element of the unwrapped phase signal.

Images