Unmanned aerial vehicle monitoring method based on extended Kalman filter and MUSIC algorithm

Abstract

The invention discloses an unmanned aerial vehicle monitoring method based on an extended Kalman filter and a MUSIC algorithm, and the method comprises the steps: a data center processes an echo signal through the MUSIC algorithm, obtains a spectral function, uniformly takes a plurality of angle values in a search window, and substitutes the angle values into a spectral function curve, and the angle value spacing is a step length; finding out a candidate peak in the spectral function curve and matching the candidate peak with the unmanned aerial vehicle, wherein a direction angle of the candidate peak is a measurement azimuth angle of the matched unmanned aerial vehicle; and then, based on the EKF in combination with the measured azimuth angle, obtaining the state vector of the unmanned aerial vehicle, and fitting the flight path of the unmanned aerial vehicle according to the state vectors at multiple moments. The search window and the step length at each moment can be adjusted for spectral function curve fitting at the next moment. According to the search window and the step length dynamically adjusted according to the target state, the accuracy of the result is improved on the premise of controlling the calculation amount; according to the method, a plurality of maximum values in the curve can be processed, and track aliasing is effectively avoided; a result is obtained through EKF, noise influence is reduced, and authenticity is guaranteed.

Description

Technical Field

[0001] This invention belongs to the field of signal processing and array processing technology, specifically relating to a UAV monitoring method based on an extended Kalman filter and the MUSIC algorithm. Background Technology

[0002] In current UAV (Unmanned Aerial Vehicle) monitoring, the Multiple Signal Classification (MUSIC) algorithm is the first commonly used algorithm for echo signal processing. This algorithm constructs a spatial spectrum function based on the echo signal and determines the azimuth angle of the UAV by the direction of arrival corresponding to the peak value in the spectrum function curve. However, with the rapid development of technology, the performance of UAVs has been greatly improved, and the existing MUSIC algorithm can no longer adequately meet the needs of real-time UAV trajectory monitoring.

[0003] First, current drones possess considerable maneuverability, with frequent abrupt changes in speed and heading during flight. For example, when flying low over cities, a drone may suddenly turn to avoid obstacles, and during reconnaissance missions, it may also engage in irregular maneuvers. The accuracy of positioning using spectral function curves largely depends on the set azimuth search range. If the search range is set too small, sudden changes in target speed and heading may cause the actual position to deviate from the prediction, resulting in missed targets. If the search range is set too large, the computational complexity increases significantly, reducing real-time performance. Second, multiple drones often enter the monitoring range simultaneously, with adjacent or overlapping flight paths. For instance, when drones employ "swarm" tactics for attack or defense, multiple targets appear in dense formations. In such cases, the spectral function curves are prone to blurred or overlapping peaks, affecting not only the differentiation of different targets but also the continuity of flight paths.

[0004] To improve target trajectory continuity and positioning accuracy, the results obtained through the MUSIC algorithm require further processing using filtering algorithms. The Extended Kalman Filter (EKF) is a widely used filtering algorithm for dynamic target state estimation. This algorithm combines measured and predicted values ​​to produce an estimate, which is then used as the output to eliminate noise interference and make the output clearer. However, in complex environments, the predicted values ​​given by the EKF are prone to significant deviations from the measured values. For example, in urban environments, UAVs may experience drastic changes in their trajectories due to strong winds or electromagnetic interference, changes that are difficult to reflect in the predicted values. This disconnect between measurement and prediction leads to insufficient robustness of the monitoring system in dynamic and complex environments. Furthermore, in multi-target monitoring, the trajectory of each target needs to be calculated independently, and the features contained in the echo signal need to be matched with each target using joint probability data association or multi-hypothesis tracking algorithms. This not only results in high computational complexity but also makes it easy for errors in algorithm matching to cause aliasing or intersection of trajectories of different targets, affecting the accuracy of the obtained trajectories.

[0005] In summary, existing monitoring signal processing methods struggle to produce accurate and highly precise flight path results with limited computational resources under complex environments and multi-target conditions. Therefore, it is necessary to propose a UAV flight path monitoring algorithm that can meet the above requirements. Summary of the Invention

[0006] The purpose of this invention is to provide a UAV monitoring method based on the Extended Kalman Filter (EKF) and the MUSIC algorithm. This method can dynamically guide the adaptive changes of the search range and step size through the state estimation function of the EKF, adjust the range of the generated spatial spectrum function curve, reduce the possibility of missed detection under conditions such as high-speed maneuvering of a single UAV, multi-UAV formation flight, and complex interference, and achieve real-time and accurate monitoring of the target under the above conditions.

[0007] The technical solution adopted in the method of this invention is as follows:

[0008] A drone monitoring method based on extended Kalman filter and MUSIC algorithm specifically includes the following steps:

[0009] S1. Construct a drone monitoring system, which includes a transmitting antenna, a receiving device, and a data center. The receiving device includes multiple receiving antennas arranged in a uniform linear array. The transmitting antenna and all receiving antennas are connected to the data center for signal transmission. All receiving antennas transmit the received echo signals to the data center in real time. The receiving antennas and any drone being monitored meet the far-field condition.

[0010] S2. The data center constructs a signal matrix based on all echo signals received at the same time, and then processes the signal matrix using the MUSIC algorithm to obtain a spectral function. Multiple angle values ​​are taken within a certain directional angle range, with equal intervals between adjacent angle values. Each angle value is substituted into the spectral function to obtain the corresponding spectral function value, thus obtaining data points. The spectral function curve is then fitted based on the data points. The aforementioned directional angle range is the search window, and the interval between adjacent angle values ​​is the step size.

[0011] S3. Determine all maxima in the spectral function curve that can indicate the location of the UAV, and call the maxima candidate peaks. Then match all candidate peaks with all monitored UAVs. The matched candidate peaks correspond one-to-one with UAVs. Then use the direction angle corresponding to the candidate peak as the measurement azimuth angle of the matched UAV.

[0012] S4. Construct a state vector to describe the motion state of the UAV, and obtain the latest state vector of the UAV based on the extended Kalman filter (EKF) combined with the measured azimuth angle; for UAVs that did not match the candidate peak in step S3, their latest state vector is obtained directly through EKF.

[0013] S5. Establish an adaptive search scale model and derive the adaptive search scale at the current time using the model. Then, adjust the two endpoints and step size of the search window according to the adaptive search scale. The adjusted search window and step size are used to fit the spectral function curve at the next time step.

[0014] S6. Set the default value for the step size; take multiple consecutive time points, with a fixed interval between adjacent time points. The echo signals received at each time point are processed according to the process described in steps S2-S5. When processing the echo signals at the initial time point, the search window is the complete monitoring range of the system, and the default value is used for the step size. Then, the UAV trajectory is fitted based on the state vectors obtained from all time points.

[0015] The MUSIC algorithm processes echo signals using a received data model of the following form.

[0016] X(k)=A(θ)·S(k)+N(k) (1) where X(k) is the received data matrix, A(θ) is the angle matrix, S(k) is the echo signal matrix, N(k) is the noise matrix, and N(k)=[N1(k), N2(k), ..., N M (k)] T N1(k), N2(k), ..., N M (k) represents the noise during the reception of echo signals by the 1st to Mth receiving antennas, where A(θ) = [α(θ1), α(θ2), ..., α(θ)]. N ], column vector Let \(M\) be the number of array elements, i.e., the number of receiving antennas in the system, \(N\) be the number of signal sources, the signals reflected by different signal sources are uncorrelated, and \(N < M\), \(j\) be the imaginary unit, \(d\) be the element spacing, \(\lambda\) be the wavelength of the echo signal, \(T\) denote the transpose, \(\theta_1,\theta_2,\cdots,\theta\) N are the azimuth angles of the 1st, 2nd, \(\cdots\), \(N\)th signal sources respectively, \(\theta\in\{\theta_1,\theta_2,\cdots,\theta\}\) N .

[0017] Then construct the covariance matrix \(R\) xx , \(R\) xx is obtained through the following formula

[0018] \(R\) XX =E[X(k)X H (k)] (2) where \(X\) H (k) is the conjugate transpose of the received data matrix, and \(E(\cdot)\) is the covariance operation.

[0019] Substitute formula (1) into formula (2) to get

[0020] \(R\) XX =E\{[A(\theta)\cdot S(k)+N(k)][A(\theta)\cdot S(k)+N(k)] H \} (3) After simplification, we get

[0021]

[0022] where \(R\) S is the eigenmatrix, its number of rows and columns are both \(M\), and \(R\) S is a full-rank matrix,[[]] is the noise power.

[0023] Perform eigenvalue decomposition on \(R\) xx to get \(R\) =λ XX v p =λ p v p , \(λ\) p is the eigenvalue, \(v\) p is the eigenvector, \(p\) is the serial number of the eigenvalue and eigenvector, \(1\leq p\leq M\). Substitute formula (4) to get In addition, since \(A(\theta)R\) S A H (\theta)\) is a matrix with the number of rows and columns both being \(N\). When \(p\in[N - 1,M]\) Therefore, when \(p\in[N + 1,M]\), \(A(\theta)R\) S A H (\theta)v p = 0. Matrices \(A\) H A and \(R\) S both have inverse matrices, and we can get \(R\) S-1 [A H (θ)A(θ)] -1 A(θ)R S A H (θ)v p =0, after simplification we get A H (θ)v p =0.

[0024] Eigenvalue v p Constructing a noise subspace basis E n E n =[v N+1 ,v N+2 ,…,v M From the above inference, it can be concluded that when θ is the azimuth angle of the signal source, E n It is orthogonal to α(θ), meaning their product is 0. However, in practice, due to system sampling errors, E n The product of the two functions is not perfectly orthogonal to α(θ), but it decreases sharply as the azimuth angle approaches the signal source. The spatial spectrum function is constructed based on this property, and its form is as follows:

[0025]

[0026] Where P(θ) is the spectral function, E n H As the conjugate transpose of the noise subspace basis, the θ corresponding to the maximum value on the P(θ) curve is the azimuth angle of the signal source. In practice, by substituting multiple angle values ​​into formula (5) to obtain data points, and by fitting the P(θ) curve with the data points to find the maximum value, the value of the azimuth angle can be inferred.

[0027] As can be inferred from the above, the accuracy of the obtained azimuth angle depends on the density of data points near the azimuth angle; the more data points there are, the closer the P(θ) curve is to the actual situation. However, considering the limited computing resources in practical applications, the computational workload of substituting a large number of angle values ​​into formula (5) across the entire detection range is enormous, and it takes a long time to calculate the P(θ) curve, affecting the real-time performance of the monitoring results. Considering that the P(θ) curve only needs to reflect the maximum value, the method of this invention determines the search window based on the azimuth angle of the target at the previous moment, i.e., the reflected signal source, i.e., the range of angle values ​​at the current moment. This avoids a large amount of computation and makes the shape of the curve at the maximum value sufficiently detailed, ensuring the accuracy of the obtained azimuth angle.

[0028] Assume that n angle values ​​are selected within the search window, and the first angle value is ε. k Then the n angle values ​​are ε in sequence. k ε k +Δθ k 、…、ε k+(n-1)Δθ k And satisfy ε k ε k +Δθ k 、…、ε k +(n-1)Δθ k Substituting these values ​​into formula (2) yields P(ε) k ), P(ε k +Δθ k ), ..., P[ε k +(n-1)Δθ k Then, using θ as the abscissa and P(θ) as the ordinate, data points are determined based on the aforementioned angle values ​​and corresponding spectral function values, and the spectral function curve is obtained by fitting the data points.

[0029] Further optimization involves determining the candidate peak in step S3, specifically including the following steps:

[0030] S3.1. Let the spectral function be denoted as P. k (θ), where θ is the azimuth angle used to measure the direction of arrival, and k is the time. Assume P... k The curve (θ) has L maxima. All maxima are numbered, and the curvature of each maxima is calculated using the following formula.

[0031]

[0032] Among them, κ k,j Let P be the curvature, j be the index of the local maximum, 1≤j≤L and j is an integer. k (θ k,j Let θ be the j-th maximum value in the spectral function curve. k,j Let θk,0 be the direction angle corresponding to the j-th maximum value, Δθk,0 be the sampling interval, and P be the direction angle corresponding to the j-th maximum value. k (θ k,j -Δθ k,0 ) and P k (θ k,j +Δθ k,0 ) are the spectral function values ​​of the two sampling points used to measure the curvature of the j-th maximum.

[0033] Then, a curvature threshold T0 is chosen, and all κ values ​​are... k,j Compare the size with T0 one by one, if κ is satisfied k,j If the condition is ≤T0, then the j-th maximum is an effective peak value j; otherwise, it is not an effective peak value.

[0034] S3.2. Treat each monitored UAV as a target, with each target having a unique number, and then calculate the Mahalanobis distance between the effective peak value and the target using the following formula.

[0035]

[0036] Where i is the target index, d k,i_j Let the Mahalanobis distance be the distance between target i and effective peak j. Let σ be the predicted azimuth angle of target i at time k-1. θ,k-1 σ is the adaptive search scale for the echo signal at time k-1. j,k Let γ be the uncertainty factor of the effective peak value j, γ be the proportionality constant, ε be the boundary constant, both γ and ε are greater than 0, and max(·) be a function.

[0037] S3.3. Determine the threshold τ based on the noise variance obtained during the MUSIC algorithm processing. g,k Then, select an effective peak value j, and successively determine all d values ​​related to the effective peak value j. k,i_j Does d satisfy? k,i_j 2 ≤τ g,k If there exists a d that satisfies the condition... k,i_j Then, the effective peak value j is denoted as candidate peak j, and in all d k,i_j After all the judgments are completed, the valid peak j is marked as judged. Then, another unmarked valid peak is taken and the above process is repeated until all valid peaks are marked as judged.

[0038] When k=1, the above process only executes step S3.1 to obtain the effective peak value in the spectral function curve.

[0039] Ideally, when the targets are relatively dispersed, a P(θ) curve has only one maximum. However, considering that when multiple targets are close together, signal processing may treat them as a single entity. Therefore, to accurately identify each target, the curve's shape at each maximum must be sufficiently detailed. This invention details the processing flow for multiple maximum values ​​on the curve. In real-world environments, especially urban environments, there is significant interference. To avoid affecting subsequent calculations, maximum values ​​caused by interference must first be excluded.

[0040] Further optimization involves the following steps in step S3: Matching candidate peaks with the target peak.

[0041] S3.4. Select a candidate peak j, and calculate the likelihood cost between each objective and candidate peak j. The calculation formula is as follows:

[0042]

[0043] Among them, C k,i_j λ is the likelihood cost, λ is the peak height factor (λ≥0), and δ is a constant (δ>0).

[0044] Determine the cost threshold τ based on the noise variance. c,k Each C obtained k,i_j respectively with τ c,k Compare, if C is satisfied k,i_j ≤τ c,k The condition then applies to C. k,i_j The corresponding target i and candidate peak j are marked as pairable; when a candidate peak j completes C with all targets k,i_j After calculation, mark the candidate peak j as calculated, and then repeat the above process for another candidate peak that has not yet been marked, until all candidate peaks are marked as calculated.

[0045] S3.5. Use the Hungarian algorithm to find the case with the most pairings among all pairing cases, and take the pairing between the target and the candidate peak in this case as the matching result; there is a pairing relationship between any target and candidate peak that achieves pairing.

[0046] When k=1, the above process is not performed. Each valid peak appearing in the spectral function curve is taken as a candidate peak, and each candidate peak directly matches a new target, and the new target is numbered.

[0047] Since it's impossible to immediately determine whether maximum values ​​caused by interference have been completely eliminated during actual monitoring, and considering the potential for abrupt changes in the UAV's trajectory, to avoid mismatches between candidate peaks and targets, it's necessary to compare the calculated azimuth angle of the candidate peak with the azimuth angle of the target at the previous moment to determine the similarity of the azimuths. Matching candidate peaks with high similarity to the target effectively avoids trajectory matching errors, eliminates the influence of interference factors, and uses the scenario with the highest number of matches as the matching result. This ensures that each target is matched with a candidate peak at the current moment, reducing the probability of trajectory matching errors.

[0048] Further optimization involves the following steps in step S4, where the latest state vector is derived via EKF:

[0049] S4.1. Establish a Cartesian coordinate system to measure the entire measured airspace, with the plane containing this coordinate system being horizontal; then, through the state vector... Let p represent the motion state of target i. i k-x p i k-y v represents the x-axis and y-axis coordinates of target i at time k. i k-x v i k-y Let x and y be the velocity components of target i at time k, respectively, in the x and y directions.

[0050] S4.2. Establish the correspondence between the azimuth of the target and the state vector, the correspondence being expressed in the following form.

[0051]

[0052] in, Let h(·) and atan2(·) be the azimuth angle of target i, where h(·) and atan2(·) are both functions, and v k To observe the noise, v k ~N(0,R k ), R k To observe the noise covariance.

[0053] S4.3. Update the predicted state vector and predicted covariance matrix using the following formula.

[0054]

[0055] in, Let be the predicted state vector of target i at time k. Let F be the estimated state vector of target i at time k-1, and let F be the update matrix. T To update the transpose of the matrix, Q is the white noise matrix, q is the white noise spectral density, and w k For process noise, w k ~N(0,Q), P k|k-1 Let P be the predicted covariance matrix at time k. k-1|k-1 Estimate the covariance matrix at time k-1, where Δt is the time interval between time k and time k-1.

[0056] In addition, an initial estimated state vector is set. and the initial covariance matrix P 0|0 If the current time is the initial time, then... As the estimated state vector from the previous moment, P 0|0 The estimated covariance matrix from the previous moment is substituted into formula (10).

[0057] S4.4. The estimated state vector and estimated covariance matrix are obtained using the following formulas.

[0058]

[0059] Among them, S k Let S be the noise matrix. k -1 H is the inverse of the noise matrix. k H is the Jacobian matrix of the observation function. k T K is the transpose of the Jacobian matrix of the observation function. k For Kalman gain, Let be the estimated state vector of target i at time k. Let i be the predicted orientation angle of target i at time k. Let I be the measured azimuth angle of target i at time k, and let P be the identity matrix. k|k Let be the estimated covariance matrix at time k.

[0060] For targets that do not match a candidate peak at time k, their predicted azimuth at time k is substituted into formula (11) as the measured azimuth; the state vector is estimated. That is, the latest state vector at time k.

[0061] EKF is a filtering algorithm that provides predicted values ​​based on data from the previous time step and combines the predicted and measured values ​​to arrive at the final result. It reduces the impact of system noise on the results while ensuring the accuracy of the obtained results. In the method of this invention, the intermediate parameters generated by EKF during the estimation process—namely, the observation function Jacobian matrix and the prediction covariance matrix—can also be used to address the uncertainty of the results. In the method of this invention, the adaptive search scale can be derived based on the observation function Jacobian matrix and the prediction covariance matrix.

[0062] Before starting the EKF calculation, an initial covariance matrix P needs to be set. 0|0 Each target whose state vector is derived through EKF needs to have a corresponding initial estimated state vector set. The x-axis and y-axis coordinates are given by the azimuth data obtained from radar or optical angle measurement systems at the start of monitoring. The initial velocity value can be set to zero when external information is lacking, or given based on the previous target azimuth data and the assumption of uniform velocity. 0|0 The initial values ​​are used to characterize the uncertainty of the above estimation, where the element describing the position variance is taken as the square of the angular accuracy of the radar or optical angle measurement system, and the element describing the velocity variance is set to a large constant to reflect the low confidence level of the velocity estimate. The values ​​set in this way... With P 0|0 Using these as initial conditions for EKF recursive calculations helps ensure the effective operation of subsequent iterations.

[0063] Further optimization, step S5 specifically includes the following steps:

[0064] S5.1. Establish an adaptive search scaling model, which is described by the following formula.

[0065] σ θ,k 2 =H k P k|k-1 H k T (13)

[0066] Where, σ θ,kFor adaptive search scale.

[0067] S5.2. Denote the search window as Each target detected by the system corresponds to an independent search window, and the step size is denoted as Δθ. k , As the first endpoint, The second endpoint, All are angle values ​​and and Δθ k All are about σ θ,k The function, With Δθ k Expressed by the following formula

[0068]

[0069] Where α and β are both proportionality coefficients, and α and β > 0, Δθ max With Δθ min These are the upper and lower limits of the step size, respectively, and clip(·) is a function.

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

[0071] 1. The method of the present invention designs a search window and a step size, and the search window and step size can be adjusted according to the target's motion state at the previous moment. Under the premise of controlling the amount of computation, it can significantly improve the shape of the maximum value in the spectral function curve to be closer to the reality, and increase the accuracy of the obtained target azimuth angle.

[0072] 2. The search window and step size, which can be adjusted according to the target's previous motion state, can ensure effective monitoring of the UAV even when its trajectory changes abruptly;

[0073] 3. The method of the present invention can process multiple maxima in the same spectral function curve, thereby identifying multiple targets that are close to each other, and can effectively avoid track overlap.

[0074] 4. The track output results are obtained through EKF, which reduces the impact of noise generated during monitoring and data processing on the results, while also ensuring the authenticity of the results. Attached Figure Description

[0075] Figure 1 A schematic diagram of the data center structure of the UAV monitoring system constructed in this invention;

[0076] Figure 2 The spatial spectral function curve corresponding to Case 1;

[0077] Figure 3 The spatial spectral function curve corresponding to scenario 2;

[0078] Figure 4 A comparison chart of the monitoring methods of this invention and traditional monitoring methods;

[0079] Figure 5 The spatial spectral function curve corresponding to case 3;

[0080] Figure 6 The spatial spectral function curve corresponding to scenario 4. Detailed Implementation

[0081] To make the objectives, technical solutions, and advantages of the present invention clearer, the technical solutions of the present invention will be clearly and completely described below through specific embodiments. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0082] Example 1:

[0083] A drone monitoring method based on extended Kalman filter and MUSIC algorithm specifically includes the following steps:

[0084] S1. Construct a drone monitoring system, which includes a transmitting antenna, a receiving device, and a data center. The receiving device includes multiple receiving antennas arranged in a uniform linear array. The transmitting antenna and all receiving antennas are connected to the data center for signal transmission. All receiving antennas transmit the received echo signals to the data center in real time. The receiving antennas and any drone being monitored meet the far-field condition.

[0085] The data center includes a MUSIC search module, a feedback update module, a prediction module, and a collaborative control module. The MUSIC search module is signal-connected to both the feedback update module and the collaborative control module. The prediction module is also signal-connected to both the feedback update module and the collaborative control module. The echo signal is transmitted to the MUSIC search module in real time.

[0086] The S2.MUSIC search module constructs a signal matrix based on all echo signals received at the same time. Then, it processes the signal matrix using the MUSIC algorithm to obtain a spectral function. It takes multiple angle values ​​within a certain directional angle range, with equal intervals between adjacent angle values. Each angle value is substituted into the spectral function to obtain the corresponding spectral function value, thus generating data points. The spectral function curve is then fitted based on these data points. The aforementioned directional angle range is the search window, and the interval between adjacent angle values ​​is the step size. Afterward, the MUSIC search module transmits the obtained spectral function curve to the feedback update module.

[0087] S3. The feedback update module filters out all maxima in the spectral function curve that can indicate the azimuth of the UAV, and calls the maxima candidate peaks. Then, it matches all candidate peaks with all monitored UAVs. The matched candidate peaks correspond one-to-one with UAVs. Then, the direction angle corresponding to the candidate peak is used as the measured azimuth angle of the matched UAV. Then, the feedback update module transmits the measured azimuth angle data of the UAV to the prediction module.

[0088] S4. The prediction module constructs a state vector to describe the motion state of the UAV, and obtains the latest state vector of the UAV based on the extended Kalman filter (EKF) and the measured azimuth angle. For UAVs that did not match the candidate peak in step S3, their latest state vector is obtained directly through EKF. Then the prediction module transmits the latest state vector and the intermediate parameters of the process of obtaining the latest state vector through EKF to the cooperative control module.

[0089] S5. The collaborative control module establishes an adaptive search scale model and uses this model in conjunction with the intermediate parameters to obtain the adaptive search scale at the current moment. Then, it adjusts the two endpoints and step size of the search window according to the adaptive search scale. After that, the collaborative control module transmits the adjusted search window and step size to the MUSIC search module. The adjusted search window and step size are used for fitting the spectral function curve at the next moment.

[0090] S6. Set the default value for the step size; take multiple consecutive time points, with a fixed interval between adjacent time points. The echo signal received at each time point is processed according to the process described in steps S2-S5. The topology of each module in the data center is as follows: Figure 1 As shown, when processing the echo signal at the initial moment, the search window is the entire monitoring range of the system, and the step size uses the default value. Then, the UAV trajectory is fitted based on the state vector obtained at all moments.

[0091] In this embodiment, the system monitors only one drone in the airspace, and this drone remains stationary in the air. This drone is designated as target 1, and the scenario is designated as situation 1. During monitoring, the system narrows the search window based on the drone's motion state, and two consecutive moments during the monitoring process are designated as moment 1 and moment 2, respectively. The search windows for moment 1 and moment 2 are shown below. Figure 2 As shown in the figure, the horizontal axis represents the azimuth angle, and the vertical axis represents the MUSIC gain value, i.e., the spectral function value. There is only one candidate peak in the figure. In this case, the azimuth angle of target 1 is 0°. The search window at time 1 is -20° to 20°, with a span of 40°. At time 2, the search window is updated to -5° to 5°, with a span of only 10°. The step size corresponding to time 2 is 0.5°, and the adaptive search scale is 1.0.

[0092] Example 2:

[0093] In this embodiment, the monitoring method and system are the same as those described in Embodiment 1. Only one UAV is monitored in the airspace, and this UAV maintains uniform linear motion. This UAV is designated as Target 2, and the scenario is designated as Situation 2. The spectral function curve obtained at a certain moment during the monitoring process is shown below. Figure 3 As shown in the figure, the horizontal axis represents the direction angle, and the vertical axis represents the MUSIC gain value. There is only one candidate peak in the figure. At this moment, the azimuth angle of target 2 is 0°, the search window is -25° to 25°, the span is 50°, the step size is 1.0°, and the adaptive search scale is 3.0.

[0094] Furthermore, in this embodiment, another monitoring system is used to monitor target 2 using a traditional monitoring method, namely, obtaining the spectral function curve solely through the MUSIC algorithm and then directly identifying the maximum value in the curve. Both the track obtained using the traditional monitoring method and the track obtained using the monitoring method of this invention are plotted on [the map / plot]. Figure 4 In the graph, the horizontal axis represents time, and the vertical axis represents the target location. Figure 4 It can be seen that the tracking error obtained by the monitoring method of the present invention is relatively smaller, and therefore more accurate.

[0095] Example 3:

[0096] In this embodiment, the monitoring method and system are the same as those described in Embodiment 1. Two drones are monitored in the airspace, both maintaining uniform linear motion, but at different speeds. These two drones are designated as Target 3 and Target 4, respectively, and the scenario is referred to as Situation 3. The spectral function curve obtained at a certain moment during the monitoring process is shown below. Figure 5 As shown in the figure, the horizontal axis represents the direction angle, and the vertical axis represents the MUSIC gain value. There are two candidate peaks in the figure. At this moment, the azimuth angles of target 3 and target 4 are -10° and 10°, respectively. The search window corresponding to target 3 is -15° to -5°, with a span of 10°. The search window corresponding to target 4 is 5° to 15°, with a span of 10° as well. The step size is 0.8°. The adaptive search scale for target 3 is 1.2, and the adaptive search scale for target 4 is 2.5.

[0097] Example 4:

[0098] In this embodiment, the monitoring method and system are the same as those described in Embodiment 1. Two drones are monitored in the airspace, and their motion states are the same as in Embodiment 3. However, strong winds blow in the monitored airspace, causing both drones to move significantly away from the receiving device in a short period. These two drones are designated as Target 5 and Target 6, respectively, and this scenario is referred to as Situation 4. The spectral function curve obtained at a certain moment during the monitoring process is shown below. Figure 6As shown in the figure, the horizontal axis represents the azimuth angle, and the vertical axis represents the MUSIC gain value. There are two candidate peaks in the figure. At this moment, the azimuth angles of targets 5 and 6 are -15° and 15° respectively. The search window corresponding to target 5 is -35° to -15°, with a span of 20°, while the search window corresponding to target 6 is 15° to 35°, also with a span of 20° and a step size of 1.5°. The adaptive search scale for target 5 is 3.0, and for target 6 it is 2.8. (Comparison) Figure 5 and Figure 6 It can be seen that the search window setting needs to ensure that the range is moderate. Although a range that is too large can reduce the probability of losing the target, it will greatly increase the amount of calculation and reduce the real-time performance of the monitoring results. A range that is too small will result in insufficient ability to deal with emergencies, and the target may be lost when it suddenly moves significantly.

Claims

1. A method for monitoring unmanned aerial vehicles (UAVs) based on an extended Kalman filter and a MUSIC algorithm, characterized in that, Specifically, it includes the following steps: S1. Construct a drone monitoring system, which includes a transmitting antenna, a receiving device, and a data center. The receiving device includes multiple receiving antennas, which are distributed in a uniform linear array. The transmitting antenna and all receiving antennas are connected to the data center for signal transmission. All receiving antennas transmit the received echo signals to the data center in real time. The receiving antennas and any drone being monitored meet the far-field condition. S2. The data center constructs a signal matrix based on all echo signals received at the same time. Then, it processes the signal matrix using the MUSIC algorithm to obtain a spectral function. Multiple angle values ​​are taken within a certain directional angle range, with equal intervals between adjacent angle values. Each angle value is substituted into the spectral function to obtain the corresponding spectral function value, thus obtaining data points. A spectral function curve is then fitted based on the data points. The aforementioned directional angle range is the search window, and the interval between adjacent angle values ​​is the step size. S3. Determine all maxima in the spectral function curve that can indicate the location of the UAV, and call the maxima candidate peaks. Then match all candidate peaks with all monitored UAVs. The matched candidate peaks correspond one-to-one with the UAVs. Then take the direction angle corresponding to the candidate peak as the measured azimuth angle of the matched UAV. S4. Construct a state vector to describe the motion state of the UAV, and obtain the latest state vector of the UAV based on the extended Kalman filter (EKF) and the measured azimuth angle; for UAVs that did not match the candidate peak in step S3, their latest state vector is obtained directly through EKF. S5. Establish an adaptive search scale model and obtain the adaptive search scale at the current time through the model. Then, adjust the two endpoints and step size of the search window according to the adaptive search scale. The adjusted search window and step size are used to fit the spectral function curve at the next time step. S6. Set the default value for the step size; take multiple consecutive time points, with a fixed interval between adjacent time points. The echo signals received at each time point are processed according to the process described in steps S2-S5. When processing the echo signals at the initial time point, the search window is the complete monitoring range of the system, and the default value is used for the step size. Then, the UAV trajectory is fitted based on the state vectors obtained from all time points.

2. The UAV monitoring method based on extended Kalman filter and MUSIC algorithm as described in claim 1, characterized in that, The process of determining candidate peaks in step S3 specifically includes the following steps: S3.

1. Let the spectral function be denoted as P. k (θ), where θ is the azimuth angle used to measure the direction of arrival, and k is the time. Assume P... k The curve (θ) has L maxima. All maxima are numbered, and the curvature of each maxima is calculated using the following formula. Among them, κ k,j Let P be the curvature, j be the index of the local maximum, 1≤j≤L and j is an integer. k (θ k,j Let θ be the j-th maximum value in the spectral function curve. k,j Let θk,0 be the direction angle corresponding to the j-th maximum value, Δθk,0 be the sampling interval, and P be the direction angle corresponding to the j-th maximum value. k (θ k,j -Δθ k,0 ) and P k (θ k,j +Δθ k,0 ) are the spectral function values ​​of the two sampling points used to measure the curvature of the j-th maximum; Then, a curvature threshold T0 is chosen, and all κ values ​​are... k,j Compare the size with T0 one by one, if κ is satisfied k,j If the condition is ≤T0, then the j-th maximum is an effective peak value j; otherwise, it is not an effective peak value. S3.

2. Treat each monitored UAV as a target, with each target having a unique number, and then calculate the Mahalanobis distance between the effective peak value and the target using the following formula. Where i is the target index, d k,i_j Let the Mahalanobis distance be the distance between target i and effective peak j. Let σ be the predicted azimuth angle of target i at time k-1. θ,k-1 σ is the adaptive search scale for the echo signal at time k-1. j,k Let γ be the uncertainty factor of the effective peak value j, γ be the proportionality constant, ε be the boundary constant, both γ and ε are greater than 0, and max(·) be a function. S3.

3. Determine the threshold τ based on the noise variance obtained during the MUSIC algorithm processing. g,k Then, select an effective peak value j, and successively determine all d values ​​related to the effective peak value j. k,i_j Does d satisfy? k,i_j 2 ≤τ g,k If there exists a d that satisfies the condition... k,i_j Then, the effective peak value j is denoted as candidate peak j, and in all d k,i_j After all the judgments are completed, the valid peak j is marked as judged. Then, another unmarked valid peak is taken and the above process is repeated until all valid peaks are marked as judged. When k=1, the above process only executes step S3.1 to obtain the effective peak value in the spectral function curve.

3. The UAV monitoring method based on extended Kalman filter and MUSIC algorithm as described in claim 2, characterized in that, The matching process between the candidate peak and the target in step S3 includes the following steps: S3.

4. Select a candidate peak j, and calculate the likelihood cost between each objective and candidate peak j. The calculation formula is as follows: Among them, C k,i_j The likelihood cost is λ, the peak height factor is λ≥0, and δ is a constant, δ>0. Determine the cost threshold τ based on the noise variance. c,k Each C obtained k,i_j respectively with τ c,k Compare, if C is satisfied k,i_j ≤τ c,k The condition then applies to C. k,i_j The corresponding target i and candidate peak j are marked as pairable; when a candidate peak j completes C with all targets k,i_j After calculation, mark the candidate peak j as calculated, and then take another candidate peak that has not yet been marked and repeat the above process until all candidate peaks are marked as calculated. S3.

5. Use the Hungarian algorithm to find the case with the most pairings among all pairing cases, and take the pairing between the target and the candidate peak in this case as the matching result; there is a pairable relationship between any target and candidate peak that achieves pairing; When k=1, the above process is not performed. Each valid peak appearing in the spectral function curve is taken as a candidate peak, and each candidate peak directly matches a new target, and the new target is numbered.

4. The UAV monitoring method based on extended Kalman filter and MUSIC algorithm as described in claim 1, characterized in that, The specific process of obtaining the latest state vector through EKF in step S4 includes the following steps: S4.

1. Establish a Cartesian coordinate system to measure the entire measured airspace, with the plane containing this coordinate system being horizontal; then, through the state vector... Let represent the motion state of target i. Let x and y be the x-axis and y-axis coordinates of target i at time k, respectively. Let x and y be the velocity components of target i at time k, respectively; S4.

2. Establish the correspondence between the azimuth of the target and the state vector, the correspondence being expressed in the following form. in, Let h(·) and atan2(·) be the azimuth angle of target i, where h(·) and atan2(·) are both functions, and v k To observe the noise, v k ~N(0,R k ), R k To observe the noise covariance; S4.

3. Update the predicted state vector and predicted covariance matrix using the following formula. in, Let be the predicted state vector of target i at time k. Let F be the estimated state vector of target i at time k-1, and let F be the update matrix. T To update the transpose of the matrix, Q is the white noise matrix, q is the white noise spectral density, and w k For process noise, w k ~N(0,Q), P k|k-1 Let P be the predicted covariance matrix at time k. k-1|k-1 Estimate the covariance matrix at time k-1, where Δt is the time interval between time k and time k-1; In addition, an initial estimated state vector is set. and the initial covariance matrix P 0|0 If the current time is the initial time, then... As the estimated state vector from the previous moment, P 0|0 Substitute the estimated covariance matrix from the previous moment into formula (5); S4.

4. The estimated state vector and estimated covariance matrix are obtained using the following formulas. Among them, S k Let S be the noise matrix. k -1 H is the inverse of the noise matrix. k H is the Jacobian matrix of the observation function. k T K is the transpose of the Jacobian matrix of the observation function. k For Kalman gain, Let be the estimated state vector of target i at time k. Let i be the predicted orientation angle of target i at time k. Let I be the measured azimuth angle of target i at time k, and let P be the identity matrix. k|k Let be the estimated covariance matrix at time k; For targets that do not match a candidate peak at time k, their predicted azimuth at time k is substituted into formula (6) as the measured azimuth; the state vector is estimated. That is, the latest state vector at time k.

5. The UAV monitoring method based on extended Kalman filter and MUSIC algorithm as described in claim 1, characterized in that, Step S5 specifically includes the following steps: S5.

1. Establish an adaptive search scaling model, which is described by the following formula. σ θ,k 2 =H k P k|k-1 H k T (8) Where, σ θ,k To adapt the search scale; S5.

2. Denote the search window as Each target detected by the system corresponds to an independent search window, and the step size is denoted as Δθ. k , As the first endpoint, The second endpoint, All are angle values ​​and and Δθ k All are about σ θ,k The function, With Δθ k Expressed by the following formula Where α and β are both proportionality coefficients, and α and β > 0, Δθ max With Δθ min These are the upper and lower limits of the step size, respectively, and clip(·) is a function.

Images