A method for multi-UAV collaborative dynamic direction of arrival estimation

Abstract

The present invention discloses a method for collaborative dynamic direction of arrival estimation of multiple unmanned aerial vehicles (UAVs). The method uses multiple UAVs, each equipped with a receiving antenna, to form a signal receiving array with a variable structure. According to the number of signal sources to be measured, a suitable array structure is selected, and a fractional low-order covariance matrix or a virtual fractional low-order covariance array of the received data under the structure is derived. Then, a maximum likelihood equation is derived, and the equation is solved through the cultural migratory bird mechanism to obtain a snapshot direction of arrival estimation result. The weighted sample fractional low-order covariance array is updated through the update equation, and the above solution process is repeated to obtain the final dynamic direction of arrival estimation result. This method solves the problem that existing direction of arrival estimation methods cannot take into account UAV maneuverability, dynamic direction of arrival estimation, and direction of arrival estimation in an impact noise environment, while breaking through the limitations of the fixed array structure.

Description

Technical Field

[0001] The present invention belongs to the field of array signal processing and relates to a method for estimating the dynamic direction of arrival of multiple unmanned aerial vehicles in a coordinated manner, in particular to a method for estimating the dynamic direction of arrival of multiple unmanned aerial vehicles based on a cultural migratory bird mechanism in an impact noise environment. Background Art

[0002] Signal direction of arrival (DOA) estimation, also known as direction finding, is an important research area in array signal processing. However, in practical applications, a relatively small number of antennas are required to measure a large number of signal sources. Therefore, special non-uniform linear arrays based on nested arrays have attracted considerable attention. Furthermore, complex noise environments exist in reality, and the targets being measured are often in motion. Therefore, it is crucial to design a dynamic DOA estimation method based on special non-uniform arrays and to generalize it to impulsive noise environments.

[0003] Coprime arrays are non-uniform arrays with closed-form solutions for their array structure and degrees of freedom. Furthermore, the large spacing between array elements reduces the impact of inter-element coupling on DOA estimation accuracy. However, traditional DOA estimation methods based on coprime arrays are mostly applicable only in Gaussian noise environments and cannot dynamically track signal sources. Furthermore, the fixed array structure limits DOA estimation performance. Therefore, it is necessary to design methods with flexible array structures that can perform dynamic DOA estimation in impulsive noise environments.

[0004] Using the maximum likelihood principle for direction-of-arrival (DOA) estimation can theoretically yield superior results and directly resolve coherent sources. However, this process requires a global minimum search for a multidimensional nonlinear optimization problem. Rapid and accurate search results are a key limitation to the development of maximum likelihood DOA estimation methods. Currently, intelligent algorithms are often used to solve these problems. However, existing intelligent algorithms suffer from various drawbacks, such as slow convergence and a tendency to fall into local optima. Therefore, a solution that overcomes these drawbacks is needed.

[0005] Through searching the existing technical literature, it is found that Zhang Yankui et al. published a paper titled "High-dimensional Direction of Arrival Estimation Algorithm Based on Coprime Array Reconstruction" in the Journal of Electronics (2018, 46(12): 2923-2929). By performing column vectorization on the covariance matrix of the received signal, a virtual array model was established, and then the virtual array flow pattern was reconstructed, thus breaking through the limitation of the existing algorithm in estimating the small number of signal sources. However, this method cannot perform accurate estimation in an impulsive noise environment. Yao Linhong et al. published a paper titled "Research on Improved TLS-ESPRIT Algorithm under Impulsive Noise" in the Journal of Changchun University (2011, 21(02): 33-36). They applied spatial smoothing theory to the TLS-ESPRIT algorithm and designed a Direction of Arrival Estimation algorithm that does not require spatial spectrum search by combining the knowledge of fractional low-order matrices. The algorithm works well in an impulsive noise environment and can estimate coherent signal sources. However, the array aperture is lost, and the number of estimable signal sources is more severely limited by the number of physical array elements.

[0006] The search results of existing literature show that the existing direction of arrival estimation methods have a narrow scope of application, high computational complexity, and fixed array structure. There is a lack of a direction of arrival estimation method that can use multiple UAVs to achieve variable array structure and dynamically track multiple signal sources using a multi-UAV array in an impact noise environment. Summary of the Invention

[0007] In response to the above-mentioned existing technologies, the technical problem to be solved by the present invention is to provide a multi-UAV dynamic direction of arrival estimation method based on the cultural migratory bird mechanism in an impact noise environment. By updating the weighted fractional low-order covariance matrix and constructing the weighted fractional low-order covariance matrix of the virtual uniform linear array, a maximum likelihood equation based on the weighted fractional low-order covariance matrix is ​​designed, and a cultural migratory bird mechanism is designed to solve it. The dynamic direction of arrival estimation results can be output quickly and accurately, and different array structures can be constructed by changing the position of the UAV to meet the direction of arrival estimation requirements under different conditions, breaking through the application limitations and performance limitations of existing direction of arrival estimation methods.

[0008] To solve the above technical problems, the present invention provides a multi-UAV collaborative dynamic direction of arrival estimation method, comprising:

[0009] Step 1: Control the position of the drone group according to the number of signal sources to be measured, place one array element on each drone, and construct the required array structure. If the number of signal sources is less than the number of array elements, proceed to step 2; otherwise, proceed to step 3.

[0010] Step 2: Establish a uniform array direction finding model and derive the fractional low-order covariance matrix and its corresponding steering matrix under the model;

[0011] Step 3: Establish a co-prime array direction-finding model, and deduce the fractional lower-order covariance matrix and its corresponding steering matrix under this model;

[0012] Step 4: Initialize the upper bound and lower bound of the direction-of-arrival search, and initialize the maximum number of iterations, search factor, and genetic factor;

[0013] Step 5: Update the fractional lower-order covariance matrix of the array γ is the genetic factor, k is the snapshot number label. When k = 1,

[0014] Step 6: Initialize the space of the migratory bird flock, calculate the fitness values of all migratory bird individuals, and initialize the belief space;

[0015] Step 7: Search the neighborhood of the bird flock and complete the update of the bird flock;

[0016] Step 8: Determine whether to replace the leading bird;

[0017] Step 9: Update the belief space, and set the belief space to be updated at a ratio of the selection probability β; <00001​​​​​​​​​​​​​​​​​​​​​​​​​​​​​ n=1,2,...,N, d is the minimum distance between two array elements, θ=[θ1,θ2,…,θ N ] is the arrival angle vector, λ is the wavelength of the signal source, X(k)=[x1(k),x2(k),…,x M (k)] T is an M×1-dimensional array snapshot vector, where k represents the snapshot number, s(k)=[s1(k),s2(k),…,s N (k)] T is the N×1 dimensional signal vector of the k-th snapshot sampling, n(k) is the M×1 dimensional independent and identically distributed additive SaS noise vector, j is the complex unit, and T is the matrix transpose; then the weighted fractional low-order covariance matrix of the k-th snapshot data is in It means multiplying the corresponding dimension elements of two vectors, abs() means taking the absolute value of each dimension of the vector in the brackets, p1 is the fractional low-order covariance parameter, and 0<p1<1.

[0022] Furthermore, in step 3, a coprime array direction finding model is established, and the fractional low-order covariance matrix and its corresponding steering matrix under the model are derived, including:

[0023] N far-field narrowband signals are incident on a coprime array, which consists of two subarrays, B and C, each of which is a uniform linear array. Assume that subarray B has 2×b1 elements and subarray C has c1 elements. The element spacing of subarray B is c1d and the element spacing of subarray C is b1d. b1 and c1 are coprime numbers. Assume λ is the wavelength of the signal source and the minimum element spacing is d. In the coprime array, the first element positions of subarrays B and C coincide, so the actual number of physical elements is The position of the array element is in is an integer, Get integer combination Taking the continuous part, we can get a set of natural numbers. Let the number of continuous natural numbers in D be M-1, then the coprime array virtualizes a uniform linear array of M elements. Let N narrowband signals be incident on the array in the far field, then the k-th snapshot data received by the coprime array is In the formula for dimensional steering vector, where the nth steering vector is n=1,2,…,N, θ=[θ1,θ2,…,θ N ] is the arrival angle vector, for dimensional array snapshot vector, where k is the number of snapshots, s(k) = [s1(k), s2(k), ..., s N (k)] Tis an N×1 dimensional signal vector, n(k) is dimensional independent and identically distributed additive SaS noise vector, then the weighted fractional low-order covariance matrix of the k-th snapshot data is Where "⊙" means multiplying the corresponding dimension elements of the two vectors, p1 is the fractional low-order covariance parameter, and 0<p1<1; The weighted fractional low-order covariance matrix C expanded to a virtual uniform array k ; in make Where E(·) represents the mathematical expectation, 1≤ρ,τ≤M, then the fractional low-order covariance matrix C of the virtual uniform linear array is k =[c1(k),c2(k),…,c M (k)], where c m (k)=[c 1m (k),c 2m (k),…,c Mm (k)] T , m=1,2,…,M, then the steering matrix corresponding to the virtual uniform linear array is A(θ)=[a(θ1),a(θ2),…,a(θ N )], the nth extended steering vector is 1≤n≤N.

[0024] Furthermore, in step 4, the initialization of the upper bound and lower bound of the direction of arrival search, the initialization of the maximum number of iterations, the search factor, and the genetic factor includes:

[0025] The upper bound of the direction of arrival search for the kth snapshot is Search lower bound for and Represent the search upper and lower bounds of the n-th dimension variable, n = 1, 2, ..., N, to determine the initial search upper bound Define the upper bound of the domain for each search variable and determine the initial search lower bound Define the lower bound of the domain for each search variable; the maximum number of iterations corresponding to the k-th snapshot is And when hour, Where χ is the search factor, [·] is the rounding function, and max(·) is the function for finding the maximum value. Initialize the genetic factor γ.

[0026] Furthermore, in step 6, the migratory bird flock space is initialized, and the fitness values ​​of all migratory birds are calculated. Initializing the belief space includes:

[0027] Initialize each dimension variable of each position of the migratory bird flock to take a random value within the search upper bound and search lower bound, the migratory bird position set Q is the total number of migratory birds in the flock, where the position of the qth migratory bird is t is the number of generations; calculate the fitness values ​​of all migratory birds, and according to the maximum likelihood principle, the fitness function of the qth migratory bird position is obtained as follows: Where trace(·) is the matrix trace operator, and the superscript H is the conjugate transpose of the matrix; t The situational knowledge represented by the position of the migratory bird with the largest fitness value in the t-th generation of bird flock is equal to the initial leading bird of this generation. The canonical knowledge of all dimensional variables of the t-generation bird flock is expressed as The n-th dimension canonical knowledge is expressed as in is the upper bound of the n-dimensional canonical knowledge of the t-th generation, is the lower bound of the n-dimensional canonical knowledge of the t-th generation, is the fitness value corresponding to the upper bound of the n-dimensional canonical knowledge of the t-th generation, is the fitness value corresponding to the lower bound of the n-th dimension normative knowledge of the t-th generation; initialize the belief space, o 1 is the initial situation knowledge represented by the position of migratory birds with the best fitness value of the first generation of bird flocks; the initial normative knowledge is in is the fitness value corresponding to the upper limit of the initial n-th dimension variable, is the fitness value corresponding to the lower limit of the initial n-th dimension variable.

[0028] Furthermore, the search for the flock neighborhood in step 7 to complete the flock update includes:

[0029] The leading bird of the tth iteration The neighborhood location set of To lead the neighborhood of birds positions, The qth one is only in the position of the flying bird The neighborhood location is The qth is only connected to the first neighbor of the flying bird Locations is the number of neighbors explored by the leader bird, o lf is the number of neighbors passed to the next bird, let the exploration step of the n-th dimension variable be n=1,2,…,N, the upper bound of the n-th dimension variable in the neighborhood Lower bound of the n-th dimension variable in the neighborhood Among them no b is the number of birds in the flock, b p is the number of neighborhoods explored by the current bird; the number of neighborhoods explored by the qth migratory bird is The evolution formula of the nth dimension of a position is:

[0030]

[0031] and are all random numbers generated in the range [0,1]. η is the weight of normative knowledge, and its value is a positive number less than 1. is a random number from the standard normal distribution, ζ* is a constant between [0,1]; then a boundary check is performed to adjust the neighboring positions that are not in the search area to the feasible region:

[0032]

[0033] To determine whether the leader bird can evolve, calculate the fitness value of the leader bird's neighborhood, select the neighborhood position with the best fitness value, compare its fitness value with that of the leader bird, and if it is better than the leader bird, then the neighborhood position is used as the new leader bird position, and the best 2×o in the remaining neighborhood of the leader bird is selected. lf The neighbors are evenly distributed to the two following birds; otherwise, the position of the leading bird is not updated, and the best 2×o in the neighborhood of the leading bird is lf The neighborhood is evenly distributed to the two following birds; judge whether the following bird can evolve, calculate the number of birds that the following bird searches for (no n -o lf ) neighbors and the last bird assigned o lf The fitness value of the neighborhood is selected, and the neighborhood position with the best fitness value is compared with the fitness value of the following bird. If it is better than the following bird, the neighborhood position is used as the new following bird position, and the best o in the remaining neighborhood of the following bird is selected. lf The neighborhood is assigned to the following bird; otherwise, the position of the following bird is not updated, and the optimal o in the neighborhood of the following bird is lf The neighborhoods are assigned to the following birds; all the migratory birds in the flock are judged in turn whether they can evolve, and one iteration of the flock is completed. The positions of the migratory birds retained after the t+1th iteration are recorded in order according to the size of their fitness.

[0034] Furthermore, the step 8 of determining whether to replace the leading bird includes:

[0035] Let n of The number of iterations of the flock after each replacement of the leading bird, that is, the number of tours, is used to determine the tour parameters. Is it greater than the number of rounds, where ξ is the maximum number of initialization cycles, is a positive integer greater than 1; if it is greater than, the leading bird is replaced with the optimal migratory bird position other than the leading bird, and the patrol parameters are reset Otherwise, the leader bird is not replaced and the tour parameter is increased by 1.

[0036] Furthermore, the update rule for updating the belief space in step nine is:

[0037] 1) If the optimal position of the t+1th iteration is The situation knowledge update rules are as follows:

[0038]

[0039] 2) Standard knowledge update: before selection A migratory bird location to update the normative knowledge, is a number between [0,1], if The migratory birds affect the lower bound of the n-th dimension normative knowledge and its fitness value, The migratory birds affect the upper bound of the n-th dimension normative knowledge and its fitness value, then:

[0040]

[0041]

[0042]

[0043]

[0044] The beneficial effects of the present invention are as follows: the present invention proposes a multi-UAV dynamic direction of arrival estimation method suitable for an impact noise environment, which specifically uses multiple UAVs, each UAV is equipped with a receiving antenna to form a signal receiving array with a variable structure, and selects a suitable array structure according to the number of signal sources to be measured, and derives the fractional low-order covariance matrix or virtual fractional low-order covariance array of the received data under the structure, and then derives the maximum likelihood equation, and solves the equation through the cultural migratory bird mechanism to obtain a snapshot direction of arrival estimation result, and updates the weighted sample fractional low-order covariance array through the update equation. Repeating the above solution process can obtain the final dynamic direction of arrival estimation result, which solves the problem that the existing direction of arrival estimation method cannot take into account the maneuverability of UAVs, dynamic direction of arrival estimation, and direction of arrival estimation in an impact noise environment, and at the same time breaks through the limitations of the fixed array structure.

[0045] Compared with the existing technology, the present invention addresses the problems of narrow applicability, high computational complexity, and fixed array structure of existing direction of arrival estimation methods. A dynamic direction of arrival estimation method for a multi-UAV cooperative coprime array based on the cultural migratory bird mechanism is designed. By utilizing the advantages of the maximum likelihood method, coherent signal sources can be directly distinguished. At the same time, according to different direction of arrival estimation requirements, the array structure can be changed by remotely controlling the UAV to quickly and accurately obtain the required direction of arrival estimation results.

[0046] The cultural bird migration algorithm designed in this paper has strong global search capabilities, good convergence and robustness, is less likely to fall into local optimal solutions when solving the maximum likelihood function, and has a short convergence time. Simulation experiments demonstrate the effectiveness of the dynamic direction-of-arrival estimation method for a multi-UAV cooperative coprime array based on the cultural bird migration mechanism. This method overcomes the limitation of traditional direction-of-arrival estimation methods, which limit the number of measured signal sources to the number of physical array elements. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] Figure 1 This is a schematic diagram of a dynamic direction of arrival estimation method for a multi-UAV cooperative mutual prime array based on the cultural migratory bird mechanism;

[0048] Figure 2 is a schematic diagram of the coprime array structure;

[0049] Figure 3 This is the flow chart of the Chinese migratory bird algorithm of the present invention;

[0050] Figure 4 This is the simulation diagram of estimating three dynamic signal sources using a 6-element coprime array in Example 1;

[0051] Figure 5 This is a simulation diagram of the 6-element uniform linear array in Example 1 using the PSO-FLOM-ML algorithm to estimate three dynamic signal sources;

[0052] Figure 6 This is the simulation diagram of estimating five dynamic signal sources using a 6-element coprime array in Example 2;

[0053] Figure 7 This is a simulation diagram of the 6-element uniform linear array in Example 2 using the PSO-FLOM-ML algorithm to estimate 5 dynamic signal sources;

[0054] Figure 8 Simulation diagram of estimating 7 dynamic signal sources using a 6-element coprime array in Example 3. DETAILED DESCRIPTION

[0055] The present invention will be further described below with reference to the accompanying drawings and embodiments.

[0056] Combine Figure 1 、 Figure 2 and Figure 3 , the present invention comprises the following steps:

[0057] Step 1: Control the position of the drone group according to the number of signal sources to be measured, place an array element on each drone, and construct the required array structure. If the number of signal sources is less than the number of array elements, execute step 2; otherwise, execute step 3.

[0058] Step 2: Establish a uniform array direction finding model and derive the fractional low-order covariance matrix and its corresponding steering matrix under the model.

[0059] Assuming that N far-field narrowband signals are incident on a uniform linear array of M elements, the received snapshot data for the kth time is X(k) = A(θ)s(k) + n(k), where A(θ) = [a(θ1), a(θ2), …, a(θ N )] is an M×N dimensional steering vector, where the nth steering vector is n=1,2,…,N, d is the minimum distance between two array elements, θ=[θ1,θ2,…,θ N ] is the arrival angle vector, λ is the wavelength of the signal source, X(k)=[x1(k),x2(k),…,x M (k)] T is an M×1-dimensional array snapshot vector, where k represents the snapshot number, s(k)=[s1(k),s2(k),…,s N (k)] T is the N×1 dimensional signal vector of the k-th snapshot sampling, n(k) is the M×1 dimensional independent and identically distributed additive SaS noise vector, j is the complex unit, and T is the matrix transpose; then the weighted fractional low-order covariance matrix of the k-th snapshot data is in It means multiplying the corresponding dimension elements of two vectors, abs() means taking the absolute value of each dimension of the vector in the brackets, p1 is the fractional low-order covariance parameter, and 0<p1<1.

[0060] Step three, establish a coprime array direction finding model and derive the fractional low-order covariance matrix and its corresponding steering matrix under the model.

[0061] N far-field narrowband signals are incident on a coprime array, which consists of two subarrays, B and C, each of which is a uniform linear array. Assume that subarray B has 2×b1 elements and subarray C has c1 elements. The element spacing of subarray B is c1d and the element spacing of subarray C is b1d. b1 and c1 are coprime numbers. Assume λ is the wavelength of the signal source and the minimum element spacing is d. In the coprime array, the first element positions of subarrays B and C coincide, so the actual number of physical elements is The position of the array element is in is an integer, Get integer combination Taking the continuous part, we can get a set of natural numbers. Let the number of continuous natural numbers in D be M-1, then the coprime array virtualizes a uniform linear array of M elements. Let N narrowband signals be incident on the array in the far field, then the k-th snapshot data received by the coprime array is In the formula for dimensional steering vector, where the nth steering vector is n=1,2,…,N, θ=[θ1,θ2,…,θ N ] is the arrival angle vector, for dimensional array snapshot vector, where k is the number of snapshots, s(k) = [s1(k), s2(k), ..., s N (k)] T is an N×1 dimensional signal vector, n(k) is dimensional independent and identically distributed additive SaS noise vector, then the weighted fractional low-order covariance matrix of the k-th snapshot data is in Indicates multiplying the corresponding dimension elements of the two vectors, p1 is the fractional low-order covariance parameter, and 0<p1<1; Weighted fractional low-order covariance matrix expanded to a virtual uniform array in make Where E(·) represents the mathematical expectation, 1≤ρ,τ≤M, then the fractional low-order covariance matrix of the virtual uniform linear array is where c m (k)=[c 1m (k),c 2m (k),…,c Mm (k)] T , m=1,2,…,M, then the steering matrix corresponding to the virtual uniform linear array is A(θ)=[a(θ1),a(θ2),…,a(θ N )], the nth extended steering vector is 1≤n≤N.

[0062] Step 4: Initialize the upper bound and lower bound of the direction of arrival search, and initialize the maximum number of iterations, search factor, and genetic factor.

[0063] The upper bound of the direction of arrival search for the kth snapshot is Search lower bound for and Represent the search upper and lower bounds of the n-th dimension variable, n = 1, 2, ..., N, to determine the initial search upper bound Define the upper bound of the domain for each search variable and determine the initial search lower bound Define the lower bound of the domain for each search variable. The maximum number of iterations corresponding to the k-th snapshot is And when hour, Where χ is the search factor, [·] is the rounding function, and max(·) is the function for finding the maximum value. Initialize the genetic factor γ.

[0064] Step 5: Update the fractional low-order covariance matrix of the array γ is the genetic factor, when k=1,

[0065] Step 6: Initialize the migratory bird flock space, calculate the fitness values ​​of all migratory bird individuals, and initialize the belief space.

[0066] Initialize each dimension variable of each position of the migratory bird flock to take a random value within the search upper bound and search lower bound, the migratory bird position set Q is the total number of migratory birds in the flock, where the position of the qth migratory bird is t is the number of generations; calculate the fitness values ​​of all migratory birds, and according to the maximum likelihood principle, the fitness function of the qth migratory bird position is obtained as follows: Where trace(·) is the matrix trace operator, and the superscript H is the conjugate transpose of the matrix. t The situational knowledge represented by the position of the migratory bird with the largest fitness value in the t-th generation of bird flock is equal to the initial leading bird of this generation. The canonical knowledge of all dimensional variables of the t-generation bird flock is expressed as n=1,2,…,N, the n-th dimension canonical knowledge is expressed as in is the upper bound of the n-dimensional canonical knowledge of the t-th generation, is the lower bound of the n-dimensional canonical knowledge of the t-th generation, is the fitness value corresponding to the upper bound of the n-dimensional canonical knowledge of the t-th generation, is the fitness value corresponding to the lower bound of the n-th dimension normative knowledge of the t-th generation. Initialize the belief space, o 1 is the initial situation knowledge represented by the position of migratory birds with the best fitness value of the first generation of bird flocks; the initial normative knowledge is in is the fitness value corresponding to the upper limit of the initial n-th dimension variable, is the fitness value corresponding to the lower limit of the initial n-th dimension variable.

[0067] Step 7: Search the neighborhood of the bird flock and complete the bird flock update.

[0068] The leading bird of the tth iteration The neighborhood location set of To lead the neighborhood of birds positions, The qth one is only in the position of the flying bird The neighborhood location is The qth one is only connected to the first one in the neighborhood of the flying bird. Locations is the number of neighbors explored by the leader bird, o lfis the number of neighbors passed to the next bird, and the exploration step of the n-th dimension variable is n=1,2,…,N, the upper bound of the n-th dimension variable in the neighborhood Lower bound of the n-th dimension variable in the neighborhood where n ob is the number of birds in the flock, b p is the number of neighborhoods explored by the current bird. The evolution formula of the nth dimension of a position is

[0069] and are all random numbers generated in the range [0,1]. η is the weight of normative knowledge, and its value is a positive number less than 1. is a random number from standard normal distribution, ζ * is a constant between [0,1]. Then a boundary check is performed to adjust the neighboring positions that are not in the search area to the feasible region:

[0070]

[0071] To determine whether the leader bird can evolve, calculate the fitness value of the leader bird's neighborhood, select the neighborhood position with the best fitness value, compare its fitness value with that of the leader bird, and if it is better than the leader bird, then the neighborhood position is used as the new leader bird position, and the best 2×o in the remaining neighborhood of the leader bird is selected. lf The neighbors are evenly distributed to the two following birds; otherwise, the position of the leading bird is not updated, and the best 2×o in the neighborhood of the leading bird is lf The neighborhoods are evenly distributed to the two following birds. To determine whether the following bird can evolve, calculate the number of times the following bird searches (n on -o lf ) neighbors and the last bird assigned o lf The fitness value of the neighborhood is selected, and the neighborhood position with the best fitness value is compared with the fitness value of the following bird. If it is better than the following bird, the neighborhood position is used as the new following bird position, and the best o in the remaining neighborhood of the following bird is selected. lf The neighborhood is assigned to the following bird; otherwise, the position of the following bird is not updated, and the optimal o in the neighborhood of the following bird is lf The neighborhood is assigned to the following birds. All migratory birds in the flock are judged in turn whether they can evolve, and one iteration of the flock is completed. The positions of the migratory birds retained after the t+1th iteration are recorded in order according to the fitness size.

[0072] Step 8: Determine whether to replace the leading bird.

[0073] Let n ofLet \(n\) be the number of iterations (rounds) of the bird flock after each replacement of the leading bird, and judge the round parameter whether it is greater than the number of rounds, where \(\xi\) is the initialized maximum number of rounds, and \(m\) is a positive integer greater than 1; if it is greater, then replace the leading bird with the optimal migratory bird position except the leading bird, and reset the round parameter Otherwise, do not replace the leading bird and increment the round parameter by 1.

[0074] Step Nine, update the belief space. Let the belief space be updated at a ratio of the selection probability \(\beta\), and the update rule is as follows:

[0075] 1) If the optimal position in the \((t + 1)\)-th iteration is then the situation knowledge update rule is as follows:

[0076]

[0077] 2) Update the normative knowledge. Select the first migratory bird positions to update the normative knowledge, where \(\alpha\) is a number between \([0, 1]\). If the \(i\)-th migratory bird affects the lower bound of the \(n\)-th dimensional normative knowledge and its fitness value, and the \(j\)-th migratory bird affects the upper bound of the \(n\)-th dimensional normative knowledge and its fitness value, then:

[0078] [[ID=,36]]

[0079]

[0080] Step Ten, judge whether the maximum number of iterations \(T\) k is reached. If so, output the position of the leading bird as the estimated direction of arrival of the \(k\)-th snapshot Otherwise, set \(t=t + 1\) and execute Step Seven.

[0081] Step Eleven, judge whether the maximum number of snapshots is reached. If so, output the estimated results of all snapshot dynamic directions of arrival; otherwise, update the upper and lower bounds of the search for the \(n\)-th dimensional variable, \(n = 1, 2, \cdots, N\), and the rule is: where when \(k = 1\), \(r\) is the search radius, \(z\) is the convergence coefficient, and \(0 < z < 1\); update the maximum number of iterations If<00004"21>then receive the data of the \((k + 1)\)-th snapshot, set \(k = k + 1\), \(t = 1\), and execute Step Five.

[0082] At Figure 4 、 Figure 5 、 Figure 6 and Figure 7 In the figure, the tracking performance of the algorithm designed by the present invention and the PSO-FLOM-ML algorithm when tracking 3 dynamic information sources and 5 dynamic information sources respectively are compared.

[0083] In the simulation experiment, a 6-element coprime array is set, and the minimum element spacing is The array element position vector is Where λ is the wavelength of the signal source. The number of snapshots is 1000, the generalized signal-to-noise ratio in the impact noise environment is 20dB, the fractional low-order covariance parameter p1 = 0.5, and the search radius η=0.06,convergence coefficient z=0.99,genetic factor γ=0.95, β=0.2,population size 31,search factor χ=1.5,impact noise characteristic index α=1.85, The number of neighbors n for each bird oi =3, the number of neighbors passed to the next bird o lf =1.

[0084] Example 1: The initial angle of the source is The k-th snapshot dynamic signal source arrival direction is n=1,2,3, each information source is independent of each other.

[0085] Example 2: The initial angle of the source is The k-th snapshot dynamic signal source arrival direction is n=1,2,3,4,5, each information source is independent of each other.

[0086] Example 3: The initial angle of the source is The k-th snapshot dynamic signal source arrival direction is n=1,2,3,4,5,6,7, each information source is independent of each other.

[0087] Depend on Figure 4 and Figure 6 It can be seen that under the impact noise environment, the method designed by the present invention can effectively track multiple dynamic targets; Figure 5 and Figure 7 It can be seen that compared with the PSO-FLOM-ML algorithm, the proposed method has better direction of arrival estimation effect. Figure 8 It can also be seen that the present invention can estimate the situation where the number of signal sources is greater than the number of array elements, while the PSO-FLOM-ML algorithm has failed in this case, which shows that the present invention can effectively expand the array aperture.

Claims

1. A method for multi-UAV collaborative dynamic direction of arrival estimation, characterized in that: include: Step 1: Control the position of the drone group according to the number of signal sources to be measured, place one array element on each drone, and construct the required array structure. If the number of signal sources is less than the number of array elements, proceed to step 2. Otherwise, go to step 3; Step 2: Establish a uniform array direction finding model and derive the fractional low-order covariance matrix and its corresponding steering matrix under the model; Step 3: Establish a coprime array direction finding model and derive the fractional low-order covariance matrix and its corresponding steering matrix under the model; Step 4: Initialize the upper bound and lower bound of the direction of arrival search, the maximum number of iterations, the search factor, and the genetic factor; Step 5: Update the fractional low-order covariance matrix of the array γ is the genetic factor, k is the snapshot number, when k=1, Step 6: Initialize the migratory bird flock space, calculate the fitness values ​​of all migratory birds, and initialize the belief space; Step 7: Search the neighborhood of the bird flock and complete the bird flock update; Step 8: Determine whether to replace the leading bird; Step 9: Update the belief space. Assume that the belief space is updated in proportion to the selection probability β. Step 10: Determine whether the maximum number of iterations T has been reached k If yes, the output is the direction of arrival of the k-th snapshot estimate. N is the number of signal sources; otherwise, let the flock algebra t = t + 1 and execute step 7; Step 11: Determine whether the maximum number of snapshots has been reached. If so, output the dynamic direction of arrival estimation results of all snapshots. Otherwise, update the n-th dimension variable to search for the upper and lower bounds, where n = 1, 2, ..., N. The rules are: in When k=1, is the search radius, z is the convergence coefficient, and 0 <z<1; Update the maximum number of iterations like but Receive the k+1th snapshot data, set k=k+1, t=1, and execute step five.

2. The method for multi-UAV collaborative dynamic direction of arrival estimation according to claim 1, characterized in that: In step 2, a uniform array direction finding model is established, and the fractional low-order covariance matrix and its corresponding steering matrix under the model are derived as follows: Assuming that N far-field narrowband signals are incident on a uniform linear array of M elements, the received snapshot data for the kth time is X(k) = A(θ)s(k) + n(k), where A(θ) = [a(θ1), a(θ2), …, a(θ N )] is an M×N dimensional steering vector, where the nth steering vector is d is the minimum distance between two array elements, θ=[θ1,θ2,…,θ N ] is the arrival angle vector, λ is the wavelength of the signal source, X(k)=[x1(k),x2(k),…,x M (k)] T is an M×1-dimensional array snapshot vector, where k represents the snapshot number, s(k)=[s1(k),s2(k),…,s N (k)] T is the N×1 dimensional signal vector of the k-th snapshot sampling, n(k) is the M×1 dimensional independent and identically distributed additive SaS noise vector, j is the complex unit, and T is the matrix transpose; then the weighted fractional low-order covariance matrix of the k-th snapshot data is Where "⊙" means multiplying the corresponding dimensional elements of the two vectors, abs() means taking the absolute value of each dimension of the vector in the brackets, p1 is the fractional low-order covariance parameter, and 0<p1<1.

3. The method for multi-UAV collaborative dynamic direction of arrival estimation according to claim 1, characterized in that: Step 3 establishes a coprime array direction finding model, and derives the fractional low-order covariance matrix and its corresponding steering matrix under the model, including: N far-field narrowband signals are incident on a coprime array, which consists of two subarrays, B and C, each of which is a uniform linear array. Assume that subarray B has 2×b1 elements and subarray C has c1 elements. The element spacing of subarray B is c1d and the element spacing of subarray C is b1d. b1 and c1 are coprime numbers. Assume λ is the wavelength of the signal source and the minimum element spacing is d. In the coprime array, the first element positions of subarrays B and C coincide, so the actual number of physical elements is The position of the array element is in is an integer, Get integer combination Taking the continuous part, we can get a set of natural numbers. Let the number of continuous natural numbers in D be M-1, then the coprime array virtualizes a uniform linear array of M elements. Let N narrowband signals be incident on the array in the far field, then the k-th snapshot data received by the coprime array is In the formula for dimensional steering vector, where the nth steering vector is θ=[θ1,θ2,…,θ N ] is the arrival angle vector, for dimensional array snapshot vector, where k is the number of snapshots, s(k) = [s1(k), s2(k), ..., s N (k)] T is an N×1 dimensional signal vector, n(k) is dimensional independent and identically distributed additive SaS noise vector, then the weighted fractional low-order covariance matrix of the k-th snapshot data is Where "⊙" means multiplying the corresponding dimension elements of the two vectors, p1 is the fractional low-order covariance parameter, and 0<p1<1; The weighted fractional low-order covariance matrix C expanded to a virtual uniform array k ; in make Where E(·) represents the mathematical expectation, 1≤ρ,τ≤M, then the fractional low-order covariance matrix C of the virtual uniform linear array is k =[c1(k),c2(k),…,c M (k)], where c m (k)=[c 1m (k),c 2m (k),…,c Mm (k)] T , m=1,2,…,M, then the steering matrix corresponding to the virtual uniform linear array is A(θ)=[a(θ1),a(θ2),…,a(θ N )], the nth extended steering vector is 1≤n≤N.

4. The method for multi-UAV collaborative dynamic direction of arrival estimation according to claim 1, characterized in that: In step 4, the initialization of the upper bound and lower bound of the direction of arrival search, the initialization of the maximum number of iterations, the search factor, and the genetic factor includes: The upper bound of the direction of arrival search for the kth snapshot is Search lower bound for and Represent the search upper and lower bounds of the n-th dimension variable, n = 1, 2, ..., N, to determine the initial search upper bound Define the upper bound of the domain for each search variable and determine the initial search lower bound Define the lower bound of the domain for each search variable; the maximum number of iterations corresponding to the k-th snapshot is And when hour, Where χ is the search factor, [·] is the rounding function, and max(·) is the function for finding the maximum value. Initialize the genetic factor γ.

5. The method for multi-UAV collaborative dynamic direction of arrival estimation according to claim 2 or 3, characterized in that: In step 6, the migratory bird flock space is initialized, and the fitness values ​​of all migratory birds are calculated. The initialization of the belief space includes: Initialize each dimension variable of each position of the migratory bird flock to take a random value within the search upper bound and search lower bound, the migratory bird position set Q is the total number of migratory birds in the flock, where the position of the qth migratory bird is t is the number of generations; calculate the fitness values ​​of all migratory birds, and according to the maximum likelihood principle, the fitness function of the qth migratory bird position is obtained as follows: Where trace(·) is the matrix trace operator, and the superscript H is the conjugate transpose of the matrix; t The situational knowledge represented by the position of the migratory bird with the largest fitness value in the t-th generation of bird flock is equal to the initial leading bird of this generation. The canonical knowledge of all dimensional variables of the t-generation bird flock is expressed as The n-th dimension canonical knowledge is expressed as in is the upper bound of the n-dimensional canonical knowledge of the t-th generation, is the lower bound of the n-dimensional canonical knowledge of the t-th generation, is the fitness value corresponding to the upper bound of the n-th dimension canonical knowledge of the t-th generation, is the fitness value corresponding to the lower bound of the n-th dimension normative knowledge of the t-th generation; initialize the belief space, o 1 is the initial situation knowledge represented by the position of migratory birds with the best fitness value of the first generation of bird flocks; the initial normative knowledge is in is the fitness value corresponding to the upper limit of the initial n-th dimension variable, is the fitness value corresponding to the lower limit of the initial n-th dimension variable.

6. The method for multi-UAV collaborative dynamic direction of arrival estimation according to claim 1, characterized in that: Step 7, searching for the flock neighborhood and completing the flock update, includes: The leading bird of the tth iteration The neighborhood location set of To lead the neighborhood of birds positions, The qth one is only in the position of the flying bird The neighborhood location of The qth one is only connected to the first one in the neighborhood of the flying bird. Locations is the number of neighbors explored by the leader bird, o lf is the number of neighbors passed to the next bird, and the exploration step of the n-th dimension variable is Upper bound of the n-th dimension variable in the neighborhood Lower bound of the n-th dimension variable in the neighborhood where n ob is the number of birds in the flock, b p is the number of neighborhoods explored by the current bird; the number of neighborhoods explored by the qth migratory bird is The evolution formula of the nth dimension of a position is: and are all random numbers generated in the range [0,1]. η is the weight of normative knowledge, and its value is a positive number less than 1. is a random number from a standard normal distribution, is a constant between [0,1]; then a boundary check is performed to adjust the neighboring positions that are not in the search area to the feasible region: To determine whether the leader bird can evolve, calculate the fitness value of the leader bird's neighborhood, select the neighborhood position with the best fitness value, compare its fitness value with that of the leader bird, and if it is better than the leader bird, then the neighborhood position is used as the new leader bird position, and the best 2×o in the remaining neighborhood of the leader bird is selected. lf The neighbors are evenly distributed to the two following birds; otherwise, the position of the leading bird is not updated, and the best 2×o in the neighborhood of the leading bird is lf The neighborhoods are evenly distributed to the two following birds; judge whether the following bird can evolve, calculate the number of neighbors searched by the following bird itself (n on -o lf ) neighbors and the last bird assigned o lf The fitness value of the neighborhood is selected, and the neighborhood position with the best fitness value is compared with the fitness value of the following bird. If it is better than the following bird, the neighborhood position is used as the new following bird position, and the best o in the remaining neighborhood of the following bird is selected. lf The neighborhood is assigned to the following bird; otherwise, the position of the following bird is not updated, and the optimal o in the neighborhood of the following bird is lf The neighborhoods are assigned to the following birds; all the migratory birds in the flock are judged in turn whether they can evolve, and one iteration of the flock is completed. The positions of the migratory birds retained after the t+1th iteration are recorded in order according to the size of their fitness.

7. The method for multi-UAV collaborative dynamic direction of arrival estimation according to claim 1, characterized in that: The step 8 of determining whether to replace the leading bird includes: Let n of The number of iterations of the flock after each replacement of the leading bird, that is, the number of tours, is used to determine the tour parameters. Is it greater than the number of rounds, where ξ is the maximum number of initialization cycles, is a positive integer greater than 1; if it is greater than, the leading bird is replaced with the optimal migratory bird position other than the leading bird, and the patrol parameters are reset Otherwise, the leader bird is not replaced and the tour parameter is increased by 1.

8. The method for multi-UAV collaborative dynamic direction of arrival estimation according to claim 1, characterized in that: The update rules for updating the belief space in step 9 are: 1) If the optimal position of the t+1th iteration is The situation knowledge update rules are as follows: 2) Standard knowledge update: before selection A migratory bird location to update normative knowledge, is a number between [0,1], if The migratory birds affect the lower bound of the n-th dimension normative knowledge and its fitness value, The migratory birds affect the upper bound of the n-th dimension normative knowledge and its fitness value, then:

Images