Time difference of arrival positioning method based on optimized station distribution and gradient search solution

Abstract

The invention relates to a time difference of arrival positioning method, in particular to a time difference of arrival positioning method based on optimized station distribution and gradient search solution. The objective of the invention is to solve the problem that the existing method only pays attention to the positioning of a certain target at a specific position and cannot give consideration to the space area of the whole unmanned aerial vehicle. The method comprises the following steps: step 1, setting a station distribution optimization model of an observation station; 2, solving the station distribution optimization model set in the step 1 to obtain an optimal coordinate set of the observation station; 3, setting a target positioning model based on the optimal coordinate set of the observation station obtained in the step 2; and 4, solving the target positioning model set in the step 3 to obtain an optimal target positioning result. The method is applied to the field of time difference of arrival positioning.

Description

Technical Field

[0001] The invention relates to a time difference of arrival positioning method. Background Art

[0002] With the development of the low-altitude economy, the number of drones has increased significantly. However, due to the lack of relevant regulations, drone interference incidents frequently occur in key areas such as airports. Therefore, it is necessary to crack down on illegal drones in key areas.

[0003] For the optimization of station layout in key areas, the geometric dilution of precision (GDOP) and Cramer-Rao lower bound (CRLB) of the target movement area are often used as fitness functions. However, since CRLB is the minimum variance that can be achieved by an unbiased estimator, there will be deviations in the actual environment. This problem is highly nonlinear and non-convex. Generally, it is iteratively updated based on a metaheuristic algorithm (MA) to achieve the optimal value of the fitness function. A metaheuristic algorithm is a type of optimization algorithm that uses relevant physical theories, swarm intelligence algorithms, or biological theories. Some classic metaheuristic algorithms are already well known, such as the genetic algorithm (GA) derived from Darwin's theory of evolution, and the particle swarm optimization (PSO) algorithm derived from the social behavior of bird flocks and fish schools. In recent years, some scholars have proposed some effective algorithms, such as Backtracking Search Optimization Algorithm (BSA), Black-winged Kite Algorithm (BKA), Secretary Bird Optimization Algorithm (SBOA), etc.

[0004] Positioning can be divided into active positioning and passive positioning. The observation station of passive positioning does not transmit signals itself, but passively receives electromagnetic waves emitted by the target, and uses the observation quantity between stations to achieve positioning. Based on the Time Difference of Arrival (TDOA), the positioning equipment is simple and has high positioning accuracy, which has been widely studied. Its solution method can be divided into analytical method, iterative method and search method. The analytical method is represented by the two-step weighted least squares method, which has small computational complexity but poor noise robustness; the iterative method is represented by the gradient descent method and Taylor expansion method, which starts from solving the nonlinearity of the positioning equation, but there may be non-convergence problems; the search method is represented by the swarm intelligence algorithm, and some studies have introduced the gradient method into the search method. For example, Salajegheh combined the quasi-Newton method with the particle swarm algorithm to obtain better performance.

[0005] Most existing studies only focus on the positioning of a target at a specific location and fail to take into account the entire key area. Summary of the invention

[0006] The purpose of the present invention is to solve the problem that the existing methods only focus on the positioning of a certain target at a specific location and cannot take into account the entire UAV space area, and propose an arrival time difference positioning method based on optimized station layout and gradient search solution.

[0007] The specific process of the arrival time difference positioning method based on optimized station layout and gradient search solution is as follows:

[0008] Step 1: Set up the optimization model of observation station layout;

[0009] Step 2: Solve the station optimization model set in step 1 to obtain the set of coordinates of the optimal observation stations;

[0010] Step 3: Based on the set of coordinates of the optimal observation station obtained in step 2, a target positioning model is set;

[0011] Step 4: Solve the target positioning model set in step 3 to obtain the optimal target positioning result.

[0012] The beneficial effects of the present invention are:

[0013] The present invention studies how to scientifically deploy stations in the UAV space area to improve the overall accuracy in the area. On this basis, a high-precision positioning solution method is studied to form a complete station deployment and positioning system.

[0014] Since CRLB is the minimum variance that can be achieved by an unbiased estimator, there may be deviations in the actual environment. Therefore, the present invention uses the average value of the root mean square error of the sampling point positioning as the fitness function.

[0015] The present invention studies the TDOA positioning system. In terms of station optimization, the average RMSE minimization of the sampling points in the UAV space area is used as the fitness function, and the SBOA algorithm is introduced for particle iteration; in terms of solving the TDOA positioning model, considering that the gradient method has a fast optimization speed but may have non-convergence problems, and the search method takes into account the global but slow convergence problem, the two are combined to form a gradient-based search algorithm, which has higher positioning accuracy than the classic analytical method and iterative method, and higher stability than the search algorithm.

[0016] In view of the possible blind spots in positioning due to regular station layout, the present invention introduces the snake heron optimization algorithm to optimize station layout in key areas, and the mean value of the root mean square error of positioning is one order of magnitude smaller than that of regular station layout. On the basis of reasonable station layout, a high-precision solution method based on the time difference of arrival positioning problem is studied. The present invention combines the advantages of fast optimization speed of gradient descent and the overall consideration of the search method, and proposes to use the gradient-based search method (Gradient-based Search, GBS) to optimize the target position. The positioning accuracy of this method is improved compared with the use of gradient and search methods alone, and the performance is more stable than the search method; this method can locate targets at all positions and can take into account the problems of the entire drone space area. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] Figure 1 It is a flow chart of the present invention;

[0018] Figure 2 It is a graph showing the parameter value changing with the number of iterations;

[0019] Figure 3 It is the average positioning RMSE iteration curve graph;

[0020] Figure 4 The average RMSE curve of different configurations as the TDOA error increases;

[0021] Figure 5 This is a comparison chart of the positioning accuracy of the five methods;

[0022] Figure 6 This is a comparison chart of the positioning accuracy of the five algorithms when the target is flying;

[0023] Figure 7 This is the correspondence diagram between the hunting behavior of snake herons and the SBOA steps. DETAILED DESCRIPTION

[0024] Specific implementation method 1: The specific process of the arrival time difference positioning method based on optimized station layout and gradient search solution in this implementation method is as follows:

[0025] Step 1: Set up the optimization model of observation station layout;

[0026] Step 2: Solve the station optimization model set in step 1 to obtain the set of coordinates of the optimal observation stations;

[0027] Step 3: Based on the set of coordinates of the optimal observation station obtained in step 2, a target positioning model is set;

[0028] Step 4: Solve the target positioning model set in step 3 to obtain the optimal target positioning result.

[0029] The targets of positioning are drones in key areas.

[0030] Specific implementation method 2: This implementation method is different from the specific implementation method 1 in that: in step 1, an optimization model for the layout of observation stations is set;

[0031] The specific process is:

[0032] Step 1. Set decision variables;

[0033] The specific process is:

[0034] The purpose of solving the optimal multi-station location is to determine the coordinates of the stations to obtain the highest average accuracy for the target area, so the coordinates of the stations are the decision variables of the mathematical model;

[0035] Assume that the number of observation stations required for the passive time difference (the time difference between the arrival of the same target signal at two observation stations) positioning system is M;

[0036] The coordinates of the i-th observation station in the three-dimensional space in the three-dimensional rectangular coordinate system are (x i ,y i ,z i ), then the decision variables of the station layout optimization model are θ=(x1,y1,z1,…,x i ,y i ,z i ,…,x M ,y M ,z M ), i=1,2,…,M, 4≤M≤8;

[0037] (x1, y1, z1) are the coordinates of the first observation station in the three-dimensional space in the three-dimensional rectangular coordinate system;

[0038] (x M ,y M ,z M ) is the coordinate of the Mth observation station in the three-dimensional space under the three-dimensional rectangular coordinate system;

[0039] Step 1 and 2: Set constraints; the specific process is:

[0040] When constructing the optimization model of passive TDOA station layout for radiation sources, it is also necessary to consider some physical factors and constraints on the optimal deployment of observation station nodes according to the actual situation of the optimization problem.

[0041] First, in order to ensure the communication distance requirements and the signal-to-noise ratio requirements of the TDOA estimation algorithm, the distance between observation stations will be subject to certain constraints; secondly, the positioning target is a drone, and the safety and performance limits of the drone during flight need to be considered, that is, the flight altitude needs to be constrained; finally, the key detection range needs to be delineated, which can be determined according to actual needs;

[0042] The constraints are shown in formula (1);

[0043]

[0044] In the formula,

[0045] R s The spatial location of the observation station determines the station layout area;

[0046] R u Space area for drones;

[0047] s i is the position of the i-th observation station;

[0048] u j o is the true position of the jth sampling point;

[0049] and are the coordinate components of the ith observation station in the three-dimensional space on the x-axis, y-axis and z-axis in the three-dimensional rectangular coordinate system; the superscript T indicates the transposition;

[0050] e1 is the x-axis unit vector of the ith observation station in the three-dimensional space under the three-dimensional rectangular coordinate system, e2 is the y-axis unit vector of the ith observation station in the three-dimensional space under the three-dimensional rectangular coordinate system, and e3 is the z-axis unit vector of the ith observation station in the three-dimensional space under the three-dimensional rectangular coordinate system;

[0051] x min is the minimum value in the x-axis direction within the station layout range, x max It is the maximum value in the x-axis direction within the station layout range;

[0052] y min is the minimum value in the y-axis direction within the station layout range, y max It is the maximum value of the y-axis direction within the range of the station layout;

[0053] z min is the minimum value in the z-axis direction within the station layout range, z max It is the maximum value in the z-axis direction within the station layout range;

[0054] Step 13: Based on decision variables and constraints, build an optimization model for the layout of observation stations. The specific process is as follows:

[0055] The optimization problem of TDOA station layout for radiation sources is actually to take the positioning accuracy obtained by the time difference passive model as the optimization objective function, and use the intelligent algorithm to infer the optimal deployment method that can improve the accuracy. Therefore, the selection of the optimization objective function directly affects the final optimization effect.

[0056] Based on UAV space area R u The root mean square error of the position calculation of all sampling points in no-fly zones (such as airports, core cities, etc.);

[0057] The average value of the root mean square error is calculated based on all the root mean square errors, and the average value of the root mean square error is used as the fitness function; then the station layout optimization model of the observation station is:

[0058]

[0059] in,

[0060] g(θ) is the fitness function of the overall positioning accuracy;

[0061] N is the UAV space area R u The total number of sampling points in , j is the UAV space area R u The jth sampling point in ;

[0062] RMSE j is the root mean square error of the positioning result of the jth sampling point, as shown in formula (3);

[0063]

[0064] in,

[0065] K is the number of repeated experiments at each sampling point;

[0066] u jk It is the position estimate obtained by using a set of TDOA observations through weighted least squares method under a specific station layout;

[0067] TDOA observation (TDOA observation is the distance difference vector d = [r 21 ,r 31 ,…,r M1 ] T ) is through the real location and the location of the cloth station;

[0068] u jk is the drone space area R uThe positioning result of the kth experiment at the jth sampling point in ;

[0069] is the true position of the jth sampling point.

[0070] K Monte Carlo repeated experiments are performed at each sampling point.

[0071] The ultimate goal of the station layout is to optimize the station layout. However, each iteration of the particles represents a set of station layout positions, which means that a set of u jk , through each iteration selection, the particle that minimizes formula (2) is the optimal station location finally determined.

[0072] The other steps and parameters are the same as those in the first embodiment.

[0073] Specific implementation method three: This implementation method is different from specific implementation methods one or two in that: in step two, the station layout optimization model set in step one is solved to obtain the set of coordinates of the optimal observation station; the specific process is:

[0074] In order to solve the high-dimensional, non-convex problem of formula (2), meta-heuristic algorithms are often used to solve the optimization problem through specific global search rules.

[0075] Step 21: Population initialization:

[0076] X i,j =lb j +ran×(ub j -lb j ),i=1,2,...,Num,j=1,2,...,D (4)

[0077] in,

[0078] X i,j represents the initialization value of the j-dimensional variable population of the i-th particle;

[0079] ub j and lb j are the upper and lower bounds of the j-th dimension variable respectively;

[0080] ran is a random number uniformly distributed between (0,1) (excluding 0 and 1);

[0081] Num is the number of populations, and D is the dimension of decision variables;

[0082] X i represents the i-th particle, each X i All have D dimensions;

[0083] All X i The Num×D-dimensional population matrix is ​​composed of, denoted as population X;

[0084] Initialize the maximum number of iterations to Time;

[0085] Initialize the global optimal value g best , the global optimal particle position p best ;

[0086]

[0087] Among them, g(X) represents the Num fitness function values ​​corresponding to the initialized Num particles, and Num represents the number of particles in each iteration;

[0088] best means the global optimum;

[0089] g is the fitness function, as shown in formula (2);

[0090] Decision variables θ=(x1,y1,z1,…,x i ,y i ,z i ,…,x M ,y M ,z M ), X includes many X i , X i =(x i ,y i ,z i ), g(X) represents the i Find the value of fitness function;

[0091] Initialize the historical global optimal value g gbest , historical global optimal particle position p gbest :

[0092]

[0093] For the particle update method, the present invention introduces a snake egret optimization algorithm to perform particle update. Each iteration of the algorithm is divided into two processes: exploration and development. Each individual of the snake egret population is updated repeatedly through the iterations of these two stages, and the optimal solution to the problem is finally found. The state of each snake egret represents the position of a set of observation stations, and the behavior of the snake egret corresponds to the update method of the particle. Let the number of particle iterations be Time.

[0094] Exploration is that the algorithm conducts a comprehensive search of the entire solution space to locate the possible optimal region, and exploitation is to guide the algorithm to conduct a fine search locally. The two need to be balanced. SBOA simulates the hunting process of secretarybirds. In the exploration stage, referring to the hunting strategy of secretarybirds, it first uses the differential evolution algorithm for global search to find potential optimal solutions, and then comprehensively utilizes the randomness of Brownian motion and the position of the current optimal solution to balance local and global. Finally, it uses the Levy flight strategy to accelerate the algorithm convergence; in the exploitation stage, referring to the escape strategy of secretarybirds, when the C1 condition is met, it conducts a fine search near the optimal solution until convergence, otherwise it randomly updates the particles by difference, taking into account exploration. The algorithm process is as Figure 7 shown.

[0095] The exploration process is the first process, simulating the three periods when secretarybirds hunt, and updating each particle X in the first process i : searching for prey, consuming prey, and attacking prey. These three stages divide the exploration process into three equal time periods; while the exploitation process is carried out in each iteration of exploration.

[0096] Step 22: When the iteration number t < Time / 3, for the i-th particle X i Adopt differential evolution global search to search for prey (differential evolution global search) and update the population position;

[0097] Step 23: When the iteration number Time / 3 ≤ t < 2Time / 3, for the updated i-th particle X obtained in Step 22 i Evolve global exploration, consume prey and update the population position;

[0098] Step 24: When the iteration number t ≥ 2Time / 3, for the updated i-th particle X obtained in Step 23 i Evolve global exploration, attack prey and update the population position;

[0099] The historical global optimal particle position p at t = Time gbest is the decision variable θ to be solved in the station layout optimization model, θ = (x1, y1, z1, …, x i , y i , z i , …, x M , y M , z M ).

[0100] Other steps and parameters are the same as those in the first or second specific implementation manner.

[0101] Specific implementation manner four: The difference between this implementation manner and one of the first to third specific implementation manners is that: in Step 22 when the iteration number t < Time / 3, for the i-th particle X iUse differential evolution global search to find prey (differential evolution global search) to update the population position; the specific process is:

[0102] A differential evolution strategy is adopted to generate new populations using the differences between individuals to enhance the diversity and global search capability of the population, and differential mutation operations are introduced to avoid falling into local optimality.

[0103] Step 221: Set the number of iterations t = 1;

[0104] Step 222: For each particle X in the population matrix X i Explore and obtain the particle X after each position update in the population matrix X i ; The specific process is:

[0105] The location update strategy is shown in formula (7):

[0106]

[0107] in,

[0108] t is the current iteration number, Time is the maximum iteration number;

[0109] X i is the value of the i-th particle at the current iteration;

[0110] represents the newly generated i-th particle;

[0111] and are two candidate solutions randomly selected from population X (representing two particles X randomly selected from population X i );

[0112] R1 is a 1×D dimensional vector, and the elements in R1 satisfy U(0,1), where U(0,1) means uniform distribution between 0 and 1 (excluding 0 and 1);

[0113] The i-th particle newly generated in stage 1 The corresponding fitness function value is the fitness function value shown in formula (2);

[0114] g(X i ) represents the i-th particle X i The corresponding fitness function value is the fitness function value shown in formula (2);

[0115] g is the fitness function, as shown in formula (2);

[0116] Development process:

[0117] Simulate the behavior of snake herons escaping from natural enemies after hunting. According to the environmental threat, there are two escape strategies: environmental camouflage and escape, corresponding to strategies C1 and C2. A dynamic disturbance factor (1-t / Time) is introduced in the camouflage behavior. 2 , so that the updated particles are at the optimal solution p best A fine search is performed nearby, and as the number of iterations increases, the search range becomes smaller until convergence; the escape strategy is similar to the idea of ​​differential evolution, which makes the algorithm also have a certain exploration ability.

[0118] Step 223: Update the particle X at each position in the population matrix X obtained in step 222 i Develop and obtain the particle X after each position update in the population matrix X i ; The specific process is:

[0119] The position update formula of the development process is shown in formula (8):

[0120]

[0121] in,

[0122] C1 represents the environmental camouflage strategy, and C2 represents the escape strategy;

[0123] represents the i-th particle newly generated in the development process;

[0124] rand is a uniformly distributed random number in the range (0,1);

[0125] R2 is a 1×D dimensional vector, and the elements in R2 satisfy U(0,1);

[0126] K represents a random selection of an integer 1 or 2;

[0127] X random represents a randomly selected particle;

[0128] Represents the i-th particle newly generated according to the development process The corresponding fitness function value is the fitness function value shown in formula (2);

[0129] The exploration process to the development process is called an iteration. For each particle X in the population matrix X i After completing the above exploration and development processes, the population X update is completed;

[0130] Step 224: Particle X after each position is updated in the population matrix X obtained in step 223 i , update the global optimal fitness function value g of the tth iteration according to formula (9) bestand the global optimal particle position p of the tth iteration best ;

[0131]

[0132] Among them, g best represents the global optimal fitness function value of the t-th iteration, g(X) represents the Num fitness function values ​​corresponding to the Num particles of the t-th iteration, Num represents the number of particles of the t-th iteration, and best represents the global optimal of the t-th iteration;

[0133] Formula 8 gives the particle X after all positions in the population matrix X are updated. i , for each particle X after the position is updated i Calculate g(X i ), all g(X i ) to form g(X), in which g(X) is selected g(X i ) minimum value as g best , g(X i ) The minimum value of X i As p best ;

[0134] Step 225

[0135] The global optimal fitness function value g of the tth iteration obtained in step 224 is best and the historical global optimal fitness function value g gbest Make comparisons;

[0136] The global optimal particle position p of the tth iteration obtained in step 224 is best and the historical global optimal particle position p gbest Make comparisons;

[0137] Update the historical global optimal fitness function value g gbest and the historical global optimal particle position p searched by the entire population so far gbest ;

[0138] It is expressed as:

[0139]

[0140] Step 226: Let the number of iterations t = t + 1, and repeat steps 222 to 226 until t = (Time / 3) - 1, and obtain the historical global optimal fitness function value g at t = (Time / 3) - 1. gbest and the historical global optimal particle position p searched by the entire population so far gbest .

[0141] The other steps and parameters are the same as those in Specific Embodiments 1 to 3.

[0142] Specific implementation mode 5: This implementation mode is different from the specific implementation modes 1 to 4 in that: in step 23, when the number of iterations Time / 3≤t<2Time / 3, the updated i-th particle X obtained in step 22 is i Evolutionary global exploration, consuming prey to update the population position; the specific process is:

[0143] Step 231: Let the number of iterations t = Time / 3;

[0144] Step 232: For each particle X in the population matrix X obtained in step 22 i Explore and obtain the particle X after each position update in the population matrix X i ; The specific process is:

[0145] Exploration process stage 2 consumes prey (Brownian motion balances exploration and exploitation):

[0146] After finding the prey, the snake heron will not hunt directly, but will hover around the prey to consume the opponent's endurance. The Brownian Motion (BM) shown in Equation (11) is introduced to simulate the random movement of the snake heron around the prey, that is, the D-dimensional Gaussian vector with a mean of 0 and a variance of 1 is best Local random search better explores the solution space near the optimal solution. At the same time, particles search based on global information and their own historical best positions, thereby increasing the chance of finding the global optimum and avoiding premature convergence to the local optimum. The position update method of stage 2 of the exploration process is shown in formula (12), adding the disturbance factor exp((t / Time) 4 ) is to flexibly control the impact of the difference between the optimal solution and the current position. As the iteration proceeds, the impact will gradually increase, and a wider range of exploration can be used to avoid falling into the local optimum too early.

[0147] The location update method is shown in formula (12):

[0148] RB=randn(1,D) (11)

[0149]

[0150] Where RB represents Brownian motion; D represents the dimension of each particle; randn(1,D) represents a row vector composed of D-dimensional Gaussian distributed random numbers;

[0151] exp((t / Time) 4 ) represents the disturbance factor; p bestrepresents the global optimal particle position of the tth iteration;

[0152] Step 233: Update the particle X at each position in the population matrix X obtained in step 232 i Develop and obtain the particle X after each position update in the population matrix X i ; The specific process is:

[0153] The development process is shown in formula (13):

[0154]

[0155] Step 234: Particle X after each position is updated in the population matrix X obtained in step 233 i , update the global optimal fitness function value g of the tth iteration according to formula (14) best and the global optimal particle position p of the tth iteration best ;

[0156]

[0157] Among them, g best represents the global optimal fitness function value of the t-th iteration, g(X) represents the Num fitness function values ​​corresponding to the Num particles of the t-th iteration, Num represents the number of particles in each iteration, and best represents the global optimal value of this iteration;

[0158] Step 235: The global optimal fitness function value g of the tth iteration obtained in step 234 is best and the historical global optimal fitness function value g gbest Make comparisons;

[0159] The global optimal particle position p of the tth iteration obtained in steps 2, 3, and 4 is best and the historical global optimal particle position p gbest Make comparisons;

[0160] Update the historical global optimal fitness function value g gbest and the historical global optimal particle position p searched by the entire population so far gbest ;

[0161] It is expressed as:

[0162]

[0163] Step 236: Let the number of iterations t = t + 1, repeat steps 232 to 236 until t = (2Time / 3) - 1, and obtain the historical global optimal fitness function value g at t = (2Time / 3) - 1 gbestand the historical global optimal particle position p searched by the entire population so far gbest .

[0164] The other steps and parameters are the same as those in Specific Implementation 1 to 4-1.

[0165] Specific implementation method 6: This implementation method is different from the specific implementation methods 1 to 5 in that: in step 24, when the number of iterations t≥2Time / 3, the updated i-th particle X obtained in step 23 i Evolutionary global exploration, attacking prey to update population position;

[0166] The historical global optimal particle position p at t = Time gbest That is, the decision variable θ to be determined in the optimization station layout modeling, θ=(x1,y1,z1,…,x i ,y i ,z i ,…,x M ,y M ,z M );

[0167] The specific process is:

[0168] Exploration process stage 3 attacking prey (Levy flight accelerated convergence):

[0169] After exhausting the prey's physical strength, the snake heron quickly attacks and kills the prey by kicking continuously. These behaviors correspond to the Levy flight strategy introduced at this stage. It is a random movement pattern characterized by short and continuous movements and occasional long-distance jumps in a short period of time near the optimal solution. Long-distance jumps enhance the ability to explore the search space, while small step sizes help improve optimization accuracy. In addition, in order to make the algorithm more dynamic and flexible and better balance the exploration and development capabilities, a nonlinear perturbation factor (1-t / Time) is introduced. (2t / Time) , its role is to reduce the influence of the current individual position and random motion when the number of iterations t increases, so the global optimal position p best The influence of will increase accordingly, which is conducive to accelerating the convergence speed at the end of the iteration.

[0170] Step 241: Let the number of iterations t = 2Time / 3;

[0171] Step 242: For each particle X in the population matrix X obtained in step 23 i Explore and obtain each updated particle X in the population matrix X i ; The specific process is:

[0172] The position update formula is shown in formula (16):

[0173]

[0174] Among them, RL is the weighted Levy flight formula, expressed as:

[0175]

[0176] in,

[0177] Levy (D) represents the D-dimensional row vector corresponding to the Levy flight, where D represents the dimension of each particle;

[0178] s and η are constant factors;

[0179] u and v are D-dimensional standard normal distribution random vectors;

[0180] σ represents the scaling factor of the step size, and the formula for σ is:

[0181]

[0182] Where Γ is the gamma function;

[0183] Step 243: Update the particle X at each position in the population matrix X obtained in step 242 i Develop and obtain the particle X after each position update in the population matrix X i ; The specific process is:

[0184] The position update formula is shown in formula (19):

[0185]

[0186] Step 244: Based on the updated particle X at each position in the population matrix X obtained in step 243 i , update the global optimal fitness function value g of the tth iteration according to formula (20) best and the global optimal particle position p of the tth iteration best ;

[0187]

[0188] in,

[0189] g(X) represents the Num fitness function values ​​corresponding to the Num particles in the tth iteration;

[0190] Num represents the number of particles in each iteration, and best represents the global optimality of this iteration;

[0191] Step 245: The global optimal fitness function value g of the tth iteration obtained in step 244 is best and the historical global optimal fitness function value g gbestMake comparisons;

[0192] The global optimal particle position p of the tth iteration obtained in step 244 is best and the historical global optimal particle position p gbest Make comparisons;

[0193] Update the historical global optimal fitness function value g gbest and the historical global optimal particle position p searched by the entire population so far gbest ;

[0194] It is expressed as:

[0195]

[0196] Step 246: Let the number of iterations t = t + 1, repeat steps 242 to 246 until t = Time, and obtain the historical global optimal fitness function value g at t = Time. gbest and the historical global optimal particle position p searched by the entire population so far gbest , the historical global optimal particle position p at t = Time gbest That is, the decision variable θ to be determined in the optimization station layout modeling, θ=(x1,y1,z1,…,x i ,y i ,z i ,…,x M ,y M ,z M ).

[0197] The other steps and parameters are the same as those in Specific Implementation Methods 1 to 5-1.

[0198] Specific implementation method 7: This implementation method is different from any one of the specific implementation methods 1 to 6 in that: in step 3, based on the set of coordinates of the optimal observation station obtained in step 2, a target positioning model (Formula (26)) is set; the specific process is:

[0199] TDOA positioning uses the distance difference between each observation station and the target to determine the target position through the intersection of spatial hyperboloids. Since there are certain errors in the observed quantity and the known observation station positions, multiple hyperboloids generally intersect in a certain area.

[0200] Step 31: Assume that the real position of the target in the three-dimensional rectangular coordinate system is u o =(x u ,y u ,z u ) T ;

[0201] There are M observation stations, located at s1=(x1,y1,z1)T , s2=(x2,y2,z2) T ,…,s i =(x i ,y i ,z i ) T ,…,s M =(x M ,y M ,z M ) T ;

[0202] Then the actual distance between the target and the i-th observation station is r i o =|u o -s i |, i = 1, ..., M, there are M distances in total, which are integrated into a matrix to obtain the true distance vector between the M observation stations and the target.

[0203] Step 3.2: Take the first observation station as the main station, then r1 o With r i o The actual distance difference between for:

[0204]

[0205] in,

[0206] c is the propagation speed of electromagnetic waves;

[0207] r1o is the actual distance between the target and the first observation station;

[0208] is the actual arrival time difference between the i-th observation station and the master station to the target;

[0209] For r1 o With r i o The actual distance difference between

[0210] Step 33: Because the TDOA measurement parameters have errors, based on the r1 obtained in step 32 o With r i o The actual distance difference between Set the observation vector with error; expressed as:

[0211] d=d o +n (23)

[0212] in,

[0213] d=[r 21 ,r 31 ,…,r M1 ] T is the distance difference vector with error, r 21 is the observed distance difference with error between the second observation station and the main station to the target, r 31 is the observed distance difference with error between the third observation station and the main station to the target, r M1 is the observed distance difference with error between the Mth observation station and the main station to the target; the superscript T indicates the transposition;

[0214] is the true distance difference vector, for With r1 o The actual distance difference between for With r1 o The actual distance difference between for With r1 o The actual distance difference between is the true distance between the target and the second observation station; is the true distance between the target and the third observation station; is the actual distance between the target and the Mth observation station;

[0215] n=[n 21 ,n 31 ,…,n M1 ] T is the observation error vector, which satisfies the mean of 0 and the covariance matrix is ​​E[nn T ] = Q; n 21 is the distance difference measurement error between the second observation station and the main station to the target, n 31 is the distance difference measurement error between the third observation station and the main station to the target, n M1 is the distance difference measurement error between the Mth observation station and the master station to the target, Q is the covariance of the distance difference measurement error, and E[] is the expected sign;

[0216] Step 34, formula (24) is the relationship between the observed value and the actual value of the distance difference. The observation vector formula (23) with error in step 33 is rewritten as the distance difference vector with error d = [r 21 ,r 31 ,…,r M1 ] T With r o (r oEach element in is related to the target real position u o There is an equal relationship)

[0217] d=Gr o +n (24)

[0218] in,

[0219] G is a (M-1)×M dimensional vector, G=[-1 M-1 I M-1 ], 1 M-1 represents the (M-1)-dimensional unit column vector, I M-1 represents the (M-1)-dimensional unit matrix; r o Represents the true distance vector between M observation stations and the target;

[0220] Each surface that determines the position u of the radiation source in space is determined by a set of distance difference measurements. When there is no measurement error, multiple surfaces intersect at a point, which means that the position of the target can be determined;

[0221] On the contrary, the surfaces determined by the measurement equations intersect to form a certain area, and the maximum likelihood problem needs to be derived based on the probability density function of the error;

[0222] Step 35: Based on the distance difference vector with error in step 34 and r o The relationship (24) is used to obtain the joint probability density function of the target position; it is expressed as:

[0223]

[0224] in,

[0225] p(d;u) represents the joint probability density function of the target position;

[0226] r represents the distance vector between M observation stations and the target. Each element in r has an equal relationship with the target position u to be determined. o is not available, so r is used here instead of r o ;

[0227] r is unknown. r is just a symbol containing the target position u. The objective function (26) contains r. By minimizing the objective function and solving r, the target position u can be obtained indirectly.

[0228] r=[r1,r2,…,r i ,…,r M ] T ;

[0229]

[0230] u is the target position to be determined, u=uo The real position of the target is obtained at that time;

[0231] Step 36: Based on the joint probability density function of the target location (25), the positioning problem is expressed as the following maximum likelihood problem:

[0232] f(u)=(d-Gr) T Q -1 (d-Gr) (26)

[0233] Where f(u) represents the objective function of the maximum likelihood positioning problem; the target position u = u o The real position of the target is obtained at that time.

[0234] The other steps and parameters are the same as those in Specific Embodiments 1 to 6.

[0235] Specific implementation eight: This implementation differs from any one of specific implementations one to seven in that: in step four, the target positioning model set in step three is solved to obtain the optimal target positioning result;

[0236] The specific process is:

[0237] The TDOA positioning problem is a maximum likelihood problem. The gradient descent method has the advantages of simple principle and small amount of calculation. The search algorithm has the characteristics of taking the whole situation into consideration, but generally has the disadvantage of slow convergence speed. The TDOA optimization objective function is a non-convex function. It is easy to fall into the local minimum when using only the gradient method. At this time, it is necessary to introduce a search algorithm to enhance the ability to explore the global space, which is called the gradient-based search method.

[0238] Step 41: Initialize the gradient vector; the specific process is:

[0239] Assume that there are N gradient vectors in the three-dimensional space under the three-dimensional rectangular coordinate system, and each gradient vector is recorded as

[0240] Y n =[Y n,1 ,Y n,2 ,Y n,3 ] T ,n=1,2,…,N (27)

[0241] Among them, Y n,1 Represents the value of the x direction of the nth gradient vector, Y n,2 Represents the value of the nth gradient vector in the y direction, Y n,3 represents the value of the nth gradient vector in the z direction, and the superscript T indicates the transpose; Y n represents the nth gradient vector;

[0242] Initialize the gradient matrix Y = [Y1, Y2, ..., YN ] T ;

[0243] Among them, Y1 represents the first gradient vector; Y2 represents the second gradient vector; Y N represents the Nth gradient vector;

[0244] Assume that the upper and lower bounds of each gradient vector are Y max and Y min ;

[0245] The total number of initialization iterations is Total;

[0246] Step 42: Set the number of iterations t′=1;

[0247] Initialize the global optimal fitness function value f best 、Global optimal solution Y best , the global worst fitness function value f worst , the global worst solution Y worst :

[0248]

[0249] Where f(Y) is the N fitness function values ​​corresponding to the N gradient vectors of the t′th iteration, and N represents the total number of gradient vectors;

[0250] Best means the global best at the t′th iteration, and Worst means the global worst at the t′th iteration;

[0251] f is the objective function of the maximum likelihood positioning problem, as shown in formula (26);

[0252] The nth gradient vector Y n represents the coordinate value of u at the t′th iteration;

[0253] Gradient matrix Y = [Y1, Y2, ..., Y N ] T , each Y n =u, calculate each Y according to formula (26) n The corresponding f(Y n ), all f(Y n ) to form f(Y), select f(Y) from f(Y) n ) minimum value as f best , f(Y n )The minimum value of Y n As Y best ;

[0254] Gradient matrix Y = [Y1, Y2, ..., Y N ] T , each Y n=u, calculate each Y according to formula (26) n The corresponding f(Y n ), all f(Y n ) to form f(Y), select f(Y) from f(Y) n ) is the maximum value of f worst , f(Y n )The maximum value of Y n As Y worst ;

[0255] Step 43: Based on the global optimal fitness function value f best 、Global optimal solution Y best , the global worst fitness function value f worst , the global worst solution Y worst , use the Gradient Search Rule (GSR) to control the movement of the gradient vector and search for possible solutions in the feasible domain;

[0256] Step 44: Set the number of iterations t′=t′+1, and repeat step 43 until t′=Total. The last iteration updates Y best This is the final target position u.

[0257] The other steps and parameters are the same as those in Specific Implementation 1 to 7-1.

[0258] Specific implementation method 9: This implementation method is different from any one of the specific implementation methods 1 to 8 in that: in step 43, based on the global optimal fitness function value f best 、Global optimal solution Y best , the global worst fitness function value f worst , the global worst solution Y worst , use the Gradient Search Rule (GSR) to control the movement of the gradient vector and search for possible solutions in the feasible domain; the specific process is:

[0259] GSR uses a gradient-based approach to enhance the search tendency and speed up the convergence, thereby obtaining a better position in the search space;

[0260] Step 431: In the three-dimensional space under the three-dimensional rectangular coordinate system, based on the initialization of the global optimal solution Y best , the global worst solution Y worst , calculate the nth gradient vector Y n The GSR operator of is:

[0261]

[0262] Among them, randn represents a standard normal distribution random number;

[0263] ΔY n is the jump step size of the gradient vector,

[0264] Among them, rand(1:N) represents N-dimensional random numbers. Represents a random vector in the current gradient matrix;

[0265] Y worst represents the global worst solution of the t′th iteration, Y best represents the global optimal solution of the t′th iteration;

[0266] ε represents a small positive number to avoid the denominator being zero;

[0267] ρ1 is an adaptive parameter introduced to improve the search capability and randomness of the GBS algorithm, expressed as

[0268]

[0269] in,

[0270] Total represents the total number of iterations; rand is a uniformly distributed random number in the range of (0,1);

[0271] β represents the parameter that changes with the number of iterations t′;

[0272] α represents a parameter;

[0273] α changes with the number of iterations. When the total number of iterations is 500, the α value changes with the number of iterations as follows: Figure 2 As shown in the figure, in the initial iteration, a larger α value is conducive to global exploration and optimization. As the number of iterations increases, α gradually decreases to promote convergence. Between 250 and 325 times, in order to get rid of the local optimum, the parameter value increases briefly, and finally the α value continues to decrease to search around the optimal solution.

[0274] Step 432: In addition to the objective function gradient, the vector (Y best -Y n ) also guides the update direction and accelerates convergence. According to the vector (Y best -Y n 0, propose a Direction of Movement (DM) formula,

[0275] Based on the global optimal solution Y best and the nth gradient vector Y n , calculate the moving direction DM from the current position to the optimal position best-cur ; The expression is:

[0276] DM best-cur =rand×ρ2×(Y best -Y n ) (31)

[0277] in,

[0278] DM best-cur Indicates the moving direction from the current position to the optimal position;

[0279] Y best represents the global optimal solution of the t′th iteration;

[0280] rand is a uniformly distributed random number in the range (0,1);

[0281] ρ2 represents the adaptive parameter, and the calculation method is the same as that of formula ρ1; it is expressed as:

[0282]

[0283] Step 4: Based on the GSR operator obtained by formula (29) and the moving direction DM of the current position pointing to the optimal position obtained by formula (31), best-cur , calculate the phased update result of the nth particle in the t′th iteration It is expressed as:

[0284]

[0285] in, represents the phased update result of the nth gradient vector at the t′th iteration, represents the position of the nth gradient vector at the t′th iteration;

[0286] Step 4: Substitute the position of the nth gradient vector at the t′th iteration in equation (33) Replaced with the global optimal solution Y of the t′th iteration history best , get the staged update result of the nth gradient vector of the t′th iteration It is expressed as follows

[0287]

[0288] in,

[0289] r1 and r2 are two different integers randomly selected in the range [1, N], where N refers to the number of gradient vectors;

[0290] Represents the phased update result 2 of the nth gradient vector of the t′th iteration;

[0291] DM r1-r2Denote a random vector To the random vector Moving direction;

[0292] Denote a random vector, Denote a random vector;

[0293] Step Four Three Five, based on And the current vector Obtain the phased update result of the nth gradient vector at the t'-th iteration It is expressed as follows

[0294]

[0295] Obviously, Is more conducive to global search, while Is more conducive to local search near the current optimal value.

[0296] Step Four Three Six, in order to balance global search and local search, based on And Obtain the position of the nth gradient vector at the (t'+1)-th iteration It is expressed as:

[0297]

[0298] Where, r a And r b Are both random numbers within (0, 1);

[0299] Step Four Three Seven, Local Escape Operator:

[0300] Set a probability parameter rate. If the random number rand1 < rate, then enter the Local Escaping Operator (LEO) stage. In the LEO stage, according to the size of rand2, determine Value;

[0301] Where rand1 and rand2 are both uniformly distributed random numbers within the range of (0, 1), Is the intermediate value obtained by using the local escape operator.

[0302] Using LEO enables the GBS algorithm to escape from the local optimum and achieve a balance between global exploration and local exploitation. Based on the historical global optimum solution Y at the t'-th iteration best 、Position update vector Two random vectors And the newly generated solution randomly shown in Equation (36) Generate a solution using the local escape operator with stronger exploration and exploitation capabilities As shown in formula (37);

[0303] Randomly generated new solutions The generation strategy is as follows

[0304]

[0305] Among them, μ is a random number in the range of (0,1);

[0306] Y rand represents a randomly generated solution, Y rand =Y min +rand(0,1)×(Y max -Y min );

[0307] represents the solution of the nth gradient vector at the t′th iteration;

[0308] Y min Represents the lower bound of each gradient vector, Y max Represents the upper bound of each gradient vector, rand(0,1) represents a random number between (0,1);

[0309]

[0310] in,

[0311] Represents an intermediate parameter, represents the position of the nth gradient vector at the t′+1th iteration;

[0312] represents the solution obtained using the local escape operator;

[0313] l1 represents a uniformly distributed random number, l1~U(-1,1), l1~U(-1,1) represents a random number that satisfies a uniform distribution between -1 and 1 (excluding 1 and -1),

[0314] l2 represents a Gaussian distributed random number, l2~N(0,1), l2~N(0,1) represents a Gaussian distributed random number with a mean of 0 and a variance of 1;

[0315] U(-1,1) indicates uniform distribution between -1 and 1 (excluding 1 and -1), and N(0,1) indicates Gaussian distribution with mean 0 and variance 1.

[0316] Both rand1 and rand2 are random numbers uniformly distributed in the range (0,1);

[0317] u1,u2,u3 are three random numbers;

[0318] rand is a uniformly distributed random number in the range (0,1);

[0319] Update the historical global optimal fitness function value f best 、The historical global optimal solution Y best , the historical global worst fitness function value f worst 、The worst solution Y in history worst :

[0320]

[0321] The other steps and parameters are the same as those in Specific Implementation 1 to 8-1.

[0322] Specific implementation method 10: This implementation method is different from any one of the specific implementation methods 1 to 9 in that: the three random numbers u1, u2, u3 satisfy

[0323]

[0324] Among them, L1 takes 0 or 1;

[0325] rand is a uniformly distributed random number in the range (0,1).

[0326] The other steps and parameters are the same as those in Specific Implementation Methods 1 to 9.

[0327] The following examples are used to verify the beneficial effects of the present invention:

[0328] Embodiment 1:

[0329] The following examples are used to verify the beneficial effects of the present invention:

[0330] Embodiment 1:

[0331] (1) Experiment on site optimization based on SBOA

[0332] The number of observation stations is set to 5, and the coordinate origin is set to coincide with the position of the main station, then the decision variable is 12-dimensional; the constraint condition is to limit the station layout area R s ∈{x∈[-5,5],y∈[-5,5],z∈[0,0.5]}km, focusing on the target area R u ∈{x∈[-8,8],y∈[-8,8],z∈[0,1]}km; it is assumed that the time difference measurement error conforms to a Gaussian distribution with a mean of 0 and a variance of 10ns, the number of iterations is set to 500, and the number of single particles is 50.

[0333] In order to compare the performance with other search methods, the particle swarm optimization algorithm and the backtracking search algorithm are introduced to solve the station arrays of the target area under the same conditions. The iteration curve is as follows Figure 3 As shown, the final error at convergence is shown in Table 1. The coordinates of the five optimized observation stations are (0, 0, 0) km, (-5, -5, 0.0018) km, (-2.8, 0.1, 0.5) km, (-5, 5, 0.0032) km, and (5, -0.15, 0.5) km.

[0334] Table 1 Average positioning error when the algorithm converges

[0335]

[0336] The error of the SBOA algorithm drops most dramatically at the beginning of the iteration. The error is always the lowest, and the error is the smallest when it finally converges. The PSO algorithm always pursues the optimal solution during the iteration process, so the convergence speed is relatively fast, but it may fall into the local optimum; the BSA algorithm has a unique backtracking memory function, so the global optimization ability is stronger, but it also brings the problem of slow convergence speed. Both algorithms depend on the choice of parameters. The SBOA algorithm combines the advantages of the two algorithms, and the performance is not affected by the initial parameter settings.

[0337] In order to verify the performance of the algorithm under different noise conditions, the time difference measurement error was set to 0-50ns, and simulation tests were conducted on three optimized station arrays and two regular station arrays, square + center point and regular pentagon, in key observation areas. And compare, the results are as follows Figure 4 Two typical regular station layout methods All of them are maintained at a high level, and the three optimized station layout methods all show better average positioning performance. Among them, the average RMSE of the station layout method optimized by the PSO algorithm is higher than that of the other two algorithms. When the TDOA measurement error is 50ns, the error is 270m higher than that of the SBOA optimized station layout method. Although the performance of the BSA and SBOA algorithms in this experiment is similar, the latter converges significantly faster than the former.

[0338] Based on the above experiments, when using the meta-heuristic algorithm to design the station configuration and performing passive positioning in key areas, the average positioning accuracy of the sampling points in the area is significantly better than that using the typical regular station formation, and the performance is more stable, which helps to achieve scientific station layout.

[0339] (2) TDOA positioning solution experiment based on GBS

[0340] The station is arranged according to the results of the SBOA algorithm, and the target radiation source position is set to (8,8,1) km. The vector scale of the GBS algorithm is 50, the number of iterations is 500, the probability parameter rate = 0.5, and the time difference measurement error range is set to 0ns~100ns. 1000 Monte Carlo experiments are performed under each error to verify the effect of GBS in solving the TDOA problem, and the performance is compared with the Chan algorithm, Taylor algorithm, gradient descent (GD), and particle swarm algorithm. The gradient descent method and Taylor algorithm use the results of the Chan algorithm as the initial value, recorded as Chan-GD and Chan-Taylor respectively. The results are shown in Figure 5 shown.

[0341] Figure 5 It is shown in the figure that the positioning accuracy of the GBS algorithm is higher than that of the iterative and analytical algorithms under any TDOA error, and the performance improvement becomes more significant as the TDOA error increases. This is because the Chan algorithm and Taylor algorithm cannot use the prior information of the target position during the solution process, so the accuracy is low. Although the gradient descent method has a known target range R u , the positioning accuracy has been improved, but due to the lack of the ability to jump out of the local optimal value, the accuracy improvement is limited; the PSO algorithm also has good positioning accuracy, but the positioning accuracy is unstable at low TDOA errors; and because the GBS algorithm combines the ideas of gradient descent and search method, the gradient-based idea enables the algorithm to skip infeasible points and move to feasible areas, while at the same time being able to utilize the global search capability based on the search method, so the GBS-based solution method has the highest accuracy and high stability. When the time difference measurement error is 99ns, the accuracy of the GBS algorithm is 255m higher than that of the gradient descent-based algorithm under the same conditions.

[0342] In order to further study the stability of each algorithm, under the condition of TDOA measurement error of 30ns, the UAV is set to fly towards the observation station at a speed of (-40, -70, 0) m / s. The experiment verifies the changes of the target positioning accuracy of the five solution methods within the observation time of 0 to 30s. Figure 6 As shown in Figure 2, Chan-GD1 and Chan-GD are the positioning results of the gradient descent method with and without the target region prior condition.

[0343] In general, as the drone flies towards the observation station, the RMSE of the positioning results of the five algorithms shows a downward trend. When there is no prior knowledge of the target position, the Chan algorithm and the Chan-GD algorithm have similar accuracy, and the Chan-Taylor algorithm has slightly higher accuracy than the above two methods. When there is prior knowledge of the target position, the positioning accuracy of the Chan-GD1 algorithm is improved by about 26m. The positioning performance of the PSO algorithm is highly unstable as the target position changes, while the GBS solution method has the highest positioning accuracy and the most stable performance among all algorithms, which shows the effectiveness of the algorithm in solving TDOA-based positioning problems.

[0344] The present invention may also have many other embodiments. Without departing from the spirit and essence of the present invention, those skilled in the art may make various corresponding changes and modifications based on the present invention, but these corresponding changes and modifications should all fall within the scope of protection of the claims attached to the present invention.

Claims

1. The arrival time difference positioning method based on optimized station layout and gradient search solution is characterized by: The specific process of the method is: Step 1: Set up the optimization model of observation station layout; Step 2: Solve the station optimization model set in step 1 to obtain the set of coordinates of the optimal observation stations; Step 3: Based on the set of coordinates of the optimal observation station obtained in step 2, a target positioning model is set; Step 4: Solve the target positioning model set in step 3 to obtain the optimal target positioning result.

2. The arrival time difference positioning method based on optimized station layout and gradient search solution according to claim 1 is characterized in that: The station layout optimization model of the observation station is set in step 1; the specific process is: Step 1: Set decision variables; the specific process is: Assume that the number of observation stations required for the passive time difference positioning system is M; The coordinates of the i-th observation station in the three-dimensional space in the three-dimensional rectangular coordinate system are (x i ,y i ,z i ), then the decision variables of the station layout optimization model are θ=(x1,y1,z1,…,x i ,y i ,z i ,…,x M ,y M ,z M ), i=1,2,…,M, 4≤M≤8; (x1, y1, z1) are the coordinates of the first observation station in the three-dimensional space in the three-dimensional rectangular coordinate system; (x M ,y M ,z M ) is the coordinate of the Mth observation station in the three-dimensional space under the three-dimensional rectangular coordinate system; Step 1 and 2: Set constraints; The specific process is: In the formula, R s is the spatial location of the observation station; R u Space area for drones; s i is the position of the i-th observation station; u j o is the actual position of the jth sampling point; and are the coordinate components of the ith observation station in the three-dimensional space on the x-axis, y-axis and z-axis in the three-dimensional rectangular coordinate system; the superscript T indicates the transposition; e1 is the x-axis unit vector of the ith observation station in the three-dimensional space under the three-dimensional rectangular coordinate system, e2 is the y-axis unit vector of the ith observation station in the three-dimensional space under the three-dimensional rectangular coordinate system, and e3 is the z-axis unit vector of the ith observation station in the three-dimensional space under the three-dimensional rectangular coordinate system; x min is the minimum value in the x-axis direction within the station layout range, x max It is the maximum value in the x-axis direction within the station layout range; y min is the minimum value in the y-axis direction within the station layout range, y max It is the maximum value of the y-axis direction within the range of the station layout; z min is the minimum value in the z-axis direction within the station layout range, z max It is the maximum value in the z-axis direction within the station layout range; Step 13: Based on decision variables and constraints, build an optimization model for the layout of observation stations. The specific process is as follows: Based on UAV space area R u Calculate the root mean square error of the positions of all sampling points within; The average value of the root mean square error is calculated based on all the root mean square errors, and the average value of the root mean square error is used as the fitness function; The optimization model of observation station layout is: in, g(θ) is the fitness function; N is the UAV space area R u The total number of sampling points in , j is the UAV space area R u The jth sampling point in ; RMSE j is the root mean square error of the positioning result of the jth sampling point, as shown in formula (3); in, K is the number of repeated experiments at each sampling point; u jk is the drone space area R u The positioning result of the kth experiment at the jth sampling point in ; is the true position of the jth sampling point.

3. The arrival time difference positioning method based on optimized station layout and gradient search solution according to claim 2 is characterized in that: In the step 2, the station layout optimization model set in the step 1 is solved to obtain the set of coordinates of the optimal observation stations; the specific process is: Step 21: Population initialization: X i,j =lb j +ran×(ub j -lb j ),i=1,2,...,Num,j=1,2,…,D (4) in, X i,j represents the initialization value of the j-dimensional variable population of the i-th particle; ub j and lb j are the upper and lower bounds of the j-th dimension variable respectively; ran is a random number uniformly distributed between (0,1); Num is the number of populations, and D is the dimension of decision variables; X i represents the i-th particle, each X i All have D dimensions; All X i The Num×D-dimensional population matrix is ​​composed of, denoted as population X; Initialize the maximum number of iterations to Time; Initialize the global optimal value g best , the global optimal particle position p best ; Among them, g(X) represents the Num fitness function values ​​corresponding to the initialized Num particles; Num represents the number of particles in each iteration; best means the global optimum; g is the fitness function, as shown in formula (2); Initialize the historical global optimal value g gbest , historical global optimal particle position p gbest : Step 22: When the number of iterations t<Time / 3, for the i-th particle X i Use differential evolution global search to find prey and update the population position; Step 23: When the number of iterations Time / 3≤t<2Time / 3, update the i-th particle X obtained in step 22 i Evolutionary global exploration, consuming prey to update population position; Step 24: When the number of iterations t≥2Time / 3, update the i-th particle X obtained in step 23 i Evolutionary global exploration, attacking prey to update population position; The historical global optimal particle position p at t = Time gbest That is, the decision variable θ to be determined in the station layout optimization model, θ = (x1, y1, z1, …, x i ,y i ,z i ,…,x M ,y M ,z M ).

4. The arrival time difference positioning method based on optimized station layout and gradient search solution according to claim 3 is characterized in that: In step 22, when the number of iterations t<Time / 3, for the i-th particle X i Use differential evolution global search to find prey and update the population position; the specific process is: Step 221: Set the number of iterations t=1; Step 222: For each particle X in the population matrix X i Explore and obtain the particle X after each position update in the population matrix X i ; The specific process is: The location update strategy is shown in formula (7): in, t is the current iteration number, Time is the maximum iteration number; X i is the value of the i-th particle at the current iteration; represents the newly generated i-th particle; and are two candidate solutions randomly selected from population X; R1 is a 1×D dimensional vector, and the elements in R1 satisfy U(0,1), where U(0,1) means uniform distribution between 0 and 1; The i-th particle newly generated in stage 1 The corresponding fitness function value; g(X i ) represents the i-th particle X i The corresponding fitness function value; Step 223: Update the particle X at each position in the population matrix X obtained in step 222 i Develop and obtain the particle X after each position update in the population matrix X i ; The specific process is: The position update formula of the development process is shown in formula (8): in, C1 represents the environmental camouflage strategy, and C2 represents the escape strategy; represents the i-th particle newly generated in the development process; rand is a uniformly distributed random number in the range (0,1); R2 is a 1×D dimensional vector, and the elements in R2 satisfy U(0,1); K represents a random selection of an integer 1 or 2; X random represents a randomly selected particle; Represents the i-th particle newly generated according to the development process The corresponding fitness function value; Step 224: Particle X after each position is updated in the population matrix X obtained in step 223 i , update the global optimal fitness function value g of the tth iteration according to formula (9) best and the global optimal particle position p of the tth iteration best ; Among them, g best Represents the global optimal fitness function value of the tth iteration; g(X) represents the Num fitness function values ​​corresponding to the Num particles in the tth iteration; Num represents the number of particles in the tth iteration, and best represents the global optimum of the tth iteration; Step 225: The global optimal fitness function value g of the tth iteration obtained in step 224 is best and the historical global optimal fitness function value g gbest Make comparisons; The global optimal particle position p of the tth iteration obtained in step 224 is best and the historical global optimal particle position p gbest Make comparisons; Update the historical global optimal fitness function value g gbest and the historical global optimal particle position p gbest ; It is expressed as: Step 226: Let the number of iterations t = t + 1, and repeat steps 222 to 226 until t = (Time / 3) - 1, and obtain the historical global optimal fitness function value g at t = (Time / 3) - 1. gbest and the historical global optimal particle position p gbest .

5. The arrival time difference positioning method based on optimized station layout and gradient search solution according to claim 4 is characterized in that: In step 23, when the number of iterations Time / 3≤t<2Time / 3, the updated i-th particle X obtained in step 22 i Evolutionary global exploration, consuming prey to update the population position; the specific process is: Step 231: Let the number of iterations t = Time / 3; Step 232: For each particle X in the population matrix X obtained in step 22 i Explore and obtain the particle X after each position update in the population matrix X i ; The specific process is: The location update method is shown in formula (12): RB=randn(1,D) (11) Where RB represents Brownian motion; D represents the dimension of each particle; randn(1,D) represents a row vector composed of D-dimensional Gaussian distributed random numbers; exp((t / Time) 4 ) represents the disturbance factor; p best represents the global optimal particle position of the tth iteration; Step 233: Update the particle X at each position in the population matrix X obtained in step 232 i Develop and obtain the particle X after each position update in the population matrix X i ; The specific process is: Step 234: Particle X after each position is updated in the population matrix X obtained in step 233 i , update the global optimal fitness function value g of the tth iteration according to formula (14) best and the global optimal particle position pb of the tth iteration est ; Among them, g(X) represents the Num fitness function values ​​corresponding to the Num particles in the tth iteration; Num represents the number of particles in each iteration, and best represents the global optimality of this iteration; Step 235: The global optimal fitness function value g of the tth iteration obtained in step 234 is best and the historical global optimal fitness function value g gbest Make comparisons; The global optimal particle position p of the tth iteration obtained in steps 2, 3, and 4 is best and the historical global optimal particle position p gbest Make comparisons; Update the historical global optimal fitness function value g gbest and the historical global optimal particle position p gbest ; It is expressed as: Step 236: Let the number of iterations t = t + 1, repeat steps 232 to 236 until t = (2Time / 3) - 1, and obtain the historical global optimal fitness function value g at t = (2Time / 3) - 1 gbest and the historical global optimal particle position p gbest .

6. The arrival time difference positioning method based on optimized station layout and gradient search solution according to claim 5 is characterized in that: In step 24, when the number of iterations t≥2Time / 3, the updated i-th particle X obtained in step 23 is i Evolutionary global exploration, attacking prey to update population position; The historical global optimal particle position p at t = Time gbest That is, the decision variable θ to be determined in the optimization station layout modeling, θ=(x1,y1,z1,…,x i ,y i ,z i ,…,x M ,y M ,z M ); The specific process is: Step 241: Let the number of iterations t = 2Time / 3; Step 242: For each particle X in the population matrix X obtained in step 23 i Explore and obtain each updated particle X in the population matrix X i ; The specific process is: The position update formula is shown in formula (16): Among them, RL is the weighted Levy flight formula, expressed as: in, Levy(D) represents the D-dimensional row vector corresponding to the Levy flight, where D represents the dimension of each particle; s and η are constant factors; u and v are D-dimensional standard normal distribution random vectors; σ represents the scaling factor of the step size, and the formula for σ is: Where Γ is the gamma function; Step 243: Update the particle X at each position in the population matrix X obtained in step 242 i Develop and obtain the particle X after each position update in the population matrix X i ; The specific process is: The position update formula is shown in formula (19): Step 244: Based on the updated particle X at each position in the population matrix X obtained in step 243 i , update the global optimal fitness function value g of the tth iteration according to formula (20) best and the global optimal particle position pb of the tth iteration est ; in, g(X) represents the Num fitness function values ​​corresponding to the Num particles in the tth iteration; Num represents the number of particles in each iteration, and best represents the global optimality of this iteration; Step 245: The global optimal fitness function value g of the tth iteration obtained in step 244 is best and the historical global optimal fitness function value g gbest Make comparisons; The global optimal particle position p of the tth iteration obtained in step 244 is best and the historical global optimal particle position p gbest Make comparisons; Update the historical global optimal fitness function value g gbest and the historical global optimal particle position p gbest ; It is expressed as: Step 246: Let the number of iterations t = t + 1, repeat steps 242 to 246 until t = Time, and obtain the historical global optimal fitness function value g at t = Time. gbest and the historical global optimal particle position p gbest , the historical global optimal particle position p at t = Time gbest That is, the decision variable θ to be determined in the optimization station layout modeling, θ=(x1,y1,z1,…,x i ,y i ,z i ,…,x M ,y M ,z M ).

7. The arrival time difference positioning method based on optimized station layout and gradient search solution according to claim 6 is characterized in that: In step 3, the target positioning model is set based on the set of coordinates of the optimal observation station obtained in step 2. The specific process is as follows: Step 31: Assume that the real position of the target in the three-dimensional rectangular coordinate system is u o =(x u ,y u ,z u ) T ; There are M observation stations, located at s1=(x1,y1,z1) T , s2=(x2,y2,z2) T ,…,s i =(x i ,y i ,z i ) T ,…,s M =(x M ,y M ,z M ) T ; Then the actual distance between the target and the i-th observation station is There are M distances in total, which are integrated into a matrix to obtain the true distance vectors between the M observation stations and the target. Step 3.2: Take the first observation station as the main station. and The actual distance difference between for: in, c is the propagation speed of electromagnetic waves; is the actual distance between the target and the first observation station; is the actual arrival time difference between the i-th observation station and the master station to the target; for and The actual distance difference between Step 33: Based on step 32 and The actual distance difference between Set the observation vector with error; expressed as: d=do+n(23) Where d = [r 21 ,r 31 ,…,r M1 ] T is the distance difference vector with error, r 21 is the observed distance difference with error between the second observation station and the main station to the target, r 31 is the observed distance difference with error between the third observation station and the main station to the target, r M1 is the observed distance difference with error between the Mth observation station and the main station to the target; the superscript T indicates the transposition; is the true distance difference vector, for and The actual distance difference between for and The actual distance difference between for and The actual distance difference between is the true distance between the target and the second observation station; is the true distance between the target and the third observation station; is the actual distance between the target and the Mth observation station; n=[n 21 ,n 31 ,…,n M1 ] T is the observation error vector, which satisfies the mean of 0 and the covariance matrix is ​​E[nn T ] = Q; n 21 is the distance difference measurement error between the second observation station and the main station to the target, n 31 is the distance difference measurement error between the third observation station and the main station to the target, n M1 is the distance difference measurement error between the Mth observation station and the master station to the target, Q is the covariance of the distance difference measurement error, and E[] is the expected sign; Step 34: Rewrite the observation vector equation (23) with error in step 33 into the distance difference vector with error d = [r 21 ,r 31 ,…,r M1 ] T With r o The relationship d=Gro+n(24) Where G is a (M-1)×M dimensional vector, G=[-1 M-1 I M-1 ], 1 M-1 represents the (M-1)-dimensional unit column vector, I M-1 represents the (M-1)-dimensional unit matrix; r o Represents the true distance vector between M observation stations and the target; Step 35: Based on the distance difference vector with error in step 34 and r o The relationship (24) is used to obtain the joint probability density function of the target position; it is expressed as: in, p(d;u) represents the joint probability density function of the target position; r represents the distance vector between M observation stations and the target; u is the target position to be sought; Step 36: Based on the joint probability density function of the target location (25), the positioning problem is expressed as the following maximum likelihood problem: f(u)=(d-Gr) T Q -1 (d-Gr) (26) Among them, f(u) represents the objective function of the maximum likelihood positioning problem.

8. The arrival time difference positioning method based on optimized station layout and gradient search solution according to claim 7 is characterized in that: In step 4, the target positioning model set in step 3 is solved to obtain the optimal target positioning result; the specific process is: Step 41: Initialize the gradient vector; the specific process is: Suppose there are N gradient vectors in the three-dimensional space under the rectangular coordinate system of the three-dimensional space, and each gradient vector is recorded as Y n =[Y n,1 ,Y n,2 ,Y n,3 ] T ,n=1,2,…,N (27) Among them, Yn,1 represents the value of the nth gradient vector in the x direction, Yn,2 represents the value of the nth gradient vector in the y direction, and Yn,3 represents the value of the nth gradient vector in the z direction; The superscript T means to find the transpose; Y n represents the nth gradient vector; Initialize the gradient matrix Y = [Y1, Y2, ..., Y N ] T ; Among them, Y1 represents the first gradient vector; Y2 represents the second gradient vector; Y N represents the Nth gradient vector; Assume that the upper and lower bounds of each gradient vector are Y max and Y min ; The total number of initialization iterations is Total; Step 42: Set the number of iterations t′=1; Initialize the global optimal fitness function value f best 、Global optimal solution Y best , the global worst fitness function value f worst , the global worst solution Y worst : in, f(Y) is the N fitness function values ​​corresponding to the N gradient vectors of the t′th iteration, where N represents the total number of gradient vectors; Best means the global best at the t′th iteration, and Worst means the global worst at the t′th iteration; f is the objective function of the maximum likelihood positioning problem, as shown in formula (26); The nth gradient vector Y n represents the coordinate value of u at the t′th iteration; Step 43: Based on the global optimal fitness function value f best 、Global optimal solution Y best , the global worst fitness function value f worst , the global worst solution Y worst , using the gradient search rule to control the movement of the gradient vector and search for possible solutions in the feasible domain; Step 44: Set the number of iterations t′=t′+1, and repeat step 43 until t′=Total. The last iteration updates Y best This is the final target position u.

9. The arrival time difference positioning method based on optimized station layout and gradient search solution according to claim 8 is characterized in that: In step 43, based on the global optimal fitness function value f best 、Global optimal solution Y best , the global worst fitness function value f worst , the global worst solution Y worst , using the gradient search rule to control the movement of the gradient vector and search for possible solutions in the feasible domain; The specific process is: Step 431: In the three-dimensional space under the three-dimensional rectangular coordinate system, based on the initialization of the global optimal solution Y best , the global worst solution Y worst , calculate the nth gradient vector Y n The GSR operator of is: Among them, randn represents a standard normal distribution random number; ΔY n is the jump step size of the gradient vector, Among them, rand(1:N) represents N-dimensional random numbers. Represents a random vector in the current gradient matrix; Y worst represents the global worst solution of the t′th iteration, Y best represents the global optimal solution of the t′th iteration; ε represents a small positive number to avoid the denominator being zero; ρ1 is an adaptive parameter introduced and is expressed as in, Total represents the total number of iterations; rand is a uniformly distributed random number in the range of (0,1); β represents the parameter that changes with the number of iterations t′; α represents a parameter; Step 432: Based on the global optimal solution Y best and the nth gradient vector Y n , calculate the moving direction DM from the current position to the optimal position best-cur ; The expression is: DM best-cur =rand×ρ2×(Y best -AND n ) (31) in, DM best-cur Indicates the moving direction from the current position to the optimal position; Y best represents the global optimal solution of the t′th iteration; rand is a uniformly distributed random number in the range (0,1); ρ2 represents the adaptive parameter; it is expressed as: Step 4: Based on the GSR operator obtained by formula (29) and the moving direction DM of the current position pointing to the optimal position obtained by formula (31), best-cur , calculate the phased update result of the nth particle in the t′th iteration It is expressed as: in, represents the phased update result of the nth gradient vector at the t′th iteration, represents the position of the nth gradient vector at the t′th iteration; Step 4: Substitute the position of the nth gradient vector at the t′th iteration in equation (33) Replaced with the global optimal solution Y of the t′th iteration history best , get the staged update result of the nth gradient vector of the t′th iteration It is expressed as follows in, r1 and r2 are two different integers randomly selected in the range [1, N], where N refers to the number of gradient vectors; Represents the phased update result 2 of the nth gradient vector of the t′th iteration; DM r1-r2 Represents a random vector To random vector The direction of movement; represents a random vector, represents a random vector; Step 4, 3, 5, based on and the current vector Get the staged update result of the nth gradient vector of the t′th iteration It is expressed as follows Step 436: Based on and Get the position of the nth gradient vector at the t′+1th iteration It is expressed as: Among them, r a and r b All are random numbers within (0,1); Step 437: Global optimal solution Y based on the t′th iteration history best , position update vector Two random vectors And the randomly generated new solution shown in equation (36) Generate a solution using the local escape operator As shown in formula (37); Randomly generated new solutions The generation strategy is as follows Among them, μ is a random number in the range of (0,1); Y rand represents a randomly generated solution, Y rand =Y min +rand(0,1)×(Y max -Y min ); represents the solution of the nth gradient vector at the t′th iteration; Y min Represents the lower bound of each gradient vector, Y max Represents the upper bound of each gradient vector, rand(0,1) represents a random number between (0,1); in, Represents an intermediate parameter, represents the position of the nth gradient vector at the t′+1th iteration; represents the solution obtained using the local escape operator; l1 represents a uniformly distributed random number, l1~U(-1,1), l1~U(-1,1) represents a random number that satisfies a uniform distribution between -1 and 1, l2 represents a Gaussian distributed random number, l2~N(0,1), l2~N(0,1) represents a Gaussian distributed random number with a mean of 0 and a variance of 1; U(-1,1) means uniform distribution between -1 and 1, N(0,1) means Gaussian distribution with mean 0 and variance 1; Both rand1 and rand2 are random numbers uniformly distributed in the range (0,1); u1,u2,u3 are three random numbers; Update the historical global optimal fitness function value f best 、The historical global optimal solution Y best , the historical global worst fitness function value f worst 、The worst solution Y in history worst :

10. The arrival time difference positioning method based on optimized station layout and gradient search solution according to claim 9 is characterized in that: The three random numbers u1, u2, u3 satisfy Among them, L1 takes 0 or 1; rand is a uniformly distributed random number in the range (0,1).

Images