CN106384152A - PF space non-cooperative target track prediction method based on firefly group optimization - Google Patents

Abstract

The invention discloses an PF space non-cooperative target track prediction method, and the method comprises the steps: an initial particle swarm{x0<1>, x0<2>,...x0<N>} is constructed according to an initial state value of a system, a weight of an initial particle x0<i> is set as 1/N, and i=1, 2...N; a particle swarm {xk<1>, xk<2>,..., xk<N>} in k time can be calculated according to a system state equation, and a normalized weight wk<i> of a particle xk<i> can be calculated; the particle swarm in the k time can be optimized according to a firefly algorithm, the process comprises the steps: in each update process, attraction degree [beta]i, j from each particle xk<i> to other partcles can be calculated, the maximum attraction degree [beta]i, max of the partical xk<i> can be selected from the [beta]i, j; when the [beta]i, max is greater than a first threshold [beta]th, position update and weight calculate can be performed on a particle which is corresponding to the [beta]i, max and has small weight; when preset time m is completed, an optimized particle swarm can be obtained; and a state mean value xk^ of system state parameters in k time can be calculated according to the optimized particle swarm. According to the invention, the firefly algorithm can be adopted to optimize a resampling process of particle filtering, the particles cannot be got rid of, system information contained by a particle with low weight can be reserved, and a particle impoverishment problem can be solved.

Description

The PF space non-cooperative target orbital prediction method being optimized based on firefly group

Technical field

The present invention relates to orbital prediction technical field, more particularly, to a kind of PF space non-cooperative being optimized based on firefly group Target track Forecasting Methodology.

Background technology

Noncooperative target refers to a class can not provide the space object of effective cooperation information, including fault satellites, lost efficacy and defend Star, space junk, enemy's spacecraft etc..Orbital prediction for noncooperative target refers to known noncooperative target at a time State on the premise of, according to dynamics of orbits model prediction its afterwards a period of time in track condition.

In general, the orbital prediction technology for noncooperative target is the relative motion shape based on C-W equation analytic solutions State equation, in conjunction with estimation theories such as EKF, Unscented kalman filtering, particle filters, to Space borne detection platform Metrical information is tracked filtering process, to obtain the orbit information of noncooperative target.Specifically, non-based on particle filter Cooperative target orbit prediction algorithm mainly includes the following steps that：Build primary group and assignment is carried out to the weights of particle；Profit Estimate the population of subsequent time with state equation and calculate the weights of corresponding particle；Resampling, carries out evaluating it to particle Afterwards, particle big for weights is replicated, to substitute some sample points little to prediction contribution.Although this resampling process is permissible Alleviate sample degenerate problem, but after multiple circulation, may lead to only remain one or several big weights particles in population.So One, although the method ensure that number of particles, lose particle diversity, and then have impact on orbital prediction precision.

The particle depletion issues existing for the existing noncooperative target orbit prediction algorithm based on particle filter, need badly A kind of not only can guarantee that number of particles but also can guarantee that the multifarious new noncooperative target orbit prediction algorithm of particle, to improve rail Road precision of prediction.

Content of the invention

It is an object of the invention to proposing a kind of new noncooperative target orbit prediction algorithm based on particle filter, with Do not give up particle, retain solution particle depletion issues while comprised system information in low weights particle.

The invention provides a kind of PF space non-cooperative target orbital prediction method being optimized based on firefly group, including：

S1, the initial state value construction primary group according to systemAnd by each primaryWeights be set toI=1,2 ... N；

S2, the state equation according to system's relative motionCalculate the population in k momentAnd each particle is calculated according to the measured value in this momentNormalization weightsI=1,2 ... N；

S3, according to glowworm swarm algorithm, the population in k moment is optimized, including：In each renewal process, calculate every Individual particleAttraction Degree β to other particles_i,j, j=1,2 ... N and j ≠ i, and from β_i,jMiddle selection particleMaximum Attraction Degree β_i,max；When maximum Attraction Degree β_i,maxMore than default first threshold β_thWhen, to β_i,maxCorresponding little weights particle carries out position Update and calculate its weights；After completing default update times m, obtain the population after optimizing；

S4, according to optimize after population estimating system state parameter the k moment state average

Preferably, in step s3, β is calculated according to formula 1_i,j；

In formula, β_i,jFor particleTo particleAttraction Degree,For particle corresponding state variable norm.

Preferably, in step s3, location updating is carried out to particle according to formula 2；

In formula,It is to particleProduce maximum Attraction Degree β_i,maxParticle, ε_iIt is the random number being produced by Gaussian Profile Vector, K is that often step updates step-length, and K=β_i,max, α is the arbitrary width factor.

Preferably, in step s 2, each particle is calculated according to formula 3,4,5Normalization weights

In formula,For important density function,Particle for the k momentState value take observation y_kProbability,Weights for the k moment after updating.

Preferably, in step s 4, calculated according to formula 6

In formula,For system status parameters the k moment state average.

Preferably, step S4 also includes：According to optimize after population estimating system state parameter the k moment covariance P_k；

In formula,ForTransposition.

Preferably, first threshold β_thMeet：β_th=(rand (1)+1)/N；In formula, rand (1) is the random number within 1.

Preferably, N meets:100≤N≤500.

Preferably, default update times m meet：50≤m≤100.

As can be seen from the above technical solutions, the PF space non-cooperative target orbital prediction method being optimized based on firefly group Mainly include the following steps that：Initial state value according to system constructs primary group, and the weights of each primary are set For 1/N；Calculate the population in k moment according to system state equation, and calculate the normalization weights of each k moment particleRoot According to glowworm swarm algorithm, the population in k moment is optimized, including：In each renewal process, calculate each particle? Big Attraction Degree β_i,max；Work as β_i,maxMore than first threshold β_thWhen, to β_i,maxCorresponding little weights particle carry out location updating and Calculate its weights；When completing default update times m, obtain the population after optimizing；Estimated according to the population after optimizing System status parameters are in the state average in k momentThe present invention passes through the resampling to particle filter using glowworm swarm algorithm Journey is optimized, and solves particle depletion issues while not giving up particle, retain low weights particle comprised system information.

Brief description

By the specific embodiment part providing referring to the drawings, the features and advantages of the present invention will become more Easy to understand, in the accompanying drawings：

Fig. 1 is the PF space non-cooperative target orbital prediction method being optimized based on firefly group in the embodiment of the present invention Schematic flow sheet；

Fig. 2 is the site error schematic diagram according to existing Forecasting Methodology；

Fig. 3 is the site error schematic diagram according to Forecasting Methodology of the present invention；

Fig. 4 is the velocity error schematic diagram according to existing Forecasting Methodology；

Fig. 5 is the velocity error schematic diagram according to Forecasting Methodology of the present invention.

Specific embodiment

With reference to the accompanying drawings the illustrative embodiments of the present invention are described in detail.Illustrative embodiments are retouched State merely for the sake of demonstration purpose, and be definitely not the restriction to the present invention and its application or usage.

Due to there are particle depletion issues, significantly shadow in the existing noncooperative target orbit prediction algorithm based on particle filter Ring precision of prediction.In consideration of it, the present inventor proposes a kind of PF space non-cooperative mesh optimizing based on firefly group Mark orbital prediction method.The main thought of the present invention is：The particle analogy of selection is become to meet the firefly of the worm cluster characteristics of motion Fireworm is individual, using individual optimizing mode in population, so that the relatively low particle of weights is moved to higher particle.Such one Come, do not give up particle, retain low weights particle comprise system information while solve particle depletion issues.

With specific embodiment, technical scheme is described in detail below in conjunction with the accompanying drawings.Fig. 1 shows this The schematic flow sheet of the PF space non-cooperative target orbital prediction method being optimized based on firefly group in bright embodiment.From Fig. 1 It can be seen that, the method mainly includes step S1～S4.

It is necessary first to determine emulation system used before carrying out orbital prediction.In particular it is necessary to determine system phase State equation to motionAnd the observational equation of system.In this embodiment it is assumed that two star distances Much smaller than star ground distance and reference orbit is circular orbit, ignore the impact perturbed to satellite execution, and carry out first-order linear, The C-W equation of classics then can be obtained.If ignoring controling power further, C-W non trivial solution analysis solution can be released, refer to formula 9.

Wherein, x, y, z is relative position,For relative velocity, n is to follow the trail of star mean orbit angular speed.If will C-W non trivial solution analyses the form that solution is write as state matrix, then the state equation of system's relative motion can be expressed as：X (k+1)=Φ (T)X(k).Wherein, state-transition matrix Φ (T) is：

$\mathrm{\&Phi;}\left(T\right)=\left[\begin{array}{cccccc}4-3\cos nt& 0& 0& \frac{\sin nT}{n}& \frac{2(1-\cos nT)}{n}& 0\\ 6(\sin nT-nT)& 1& 0& \frac{2(\cos nT-1)}{n}& \frac{4\sin nT}{n}-3T& 0\\ 0& 0& \cos nT& 0& 0& \frac{\sin nT}{nT}\\ 3n\sin nT& 0& 0& \cos nT& 2\sin nT& 0\\ 6n(\cos nT-1)& 0& 0& -2\sin nT& 4\cos nT-3& 0\\ 0& 0& -n\sin nT& 0& 0& \cos nT\end{array}\right]$

In addition, in this embodiment, the observational equation of system can be set to Y_k=h_k(x_k,t)+v_k, refer to formula 10, system Noise may be set to：w_k～(0, Q_k), v_k～(0, R_k).

After determining emulation system used, below for this system, step S1～S4 is illustrated.

Step S1, the initial state value construction primary group according to systemAnd each is initial ParticleWeights be set to 1/N, i=1,2 ... N.

Specifically, in step sl, the initial state value of system refer to through pretreatment input data, such as position, Speed data etc..In actual track forecasting research, the measuring cell that space-based platform observation is taken is angle measurement camera and laser Rangefinder.The angle information for observation platform for the extraterrestrial target can be obtained by angle measurement camera, permissible by laser range finder Obtain the range information that extraterrestrial target is with respect to observation platform.Characteristic according to this two measuring cells sets up measurement mould respectively Type, can obtain input data.After obtaining input data, need to take some algorithms that input data is pre-processed, To improve efficiency and the precision of tracking filter algorithm, such as improved Laplace method.When constructing primary group, particle Number N can be chosen in COMPREHENSIVE CALCULATING amount and required precision.Such as, N can take more than or equal to 100 and whole less than or equal to 500 Number.

Step S2, the state equation according to system's relative motionCalculate the population in k momentAnd each particle is calculated according to the measured value in this momentNormalization weightsI=1,2 ... N.

Specifically, in step s 2, the particle in k-1 moment is substituted into system's relative motion state equation, you can when drawing k The particle carved.In addition, the measured value in this moment can be drawn according to observational equation, and then each grain can be calculated according to formula 3,4,5 SonNormalization weights

In formula,For important density function,Particle for the k momentState value take observation y_kProbability,Weights for the k moment after updating.

Can be seen that the priori value that particle substitution state equation can be obtained particle from step S2, through observational equation To this particle evaluation can be obtained again.So, just information of forecasting is merged in the distribution of particle, and observation is believed Breath has been integrated in the weight of each particle.

Step S3, according to glowworm swarm algorithm, the population in k moment is optimized, including：In each renewal process, meter Calculate each particleAttraction Degree β to other particles_i,j, j=1,2 ... N and j ≠ i, and from β_i,jMiddle selection particleMaximum suction Degree of drawing β_i,max；When maximum Attraction Degree β_i,maxMore than default first threshold β_thWhen, to β_i,maxCorresponding little weights particle is carried out Location updating simultaneously calculates its weights；After completing default update times m, obtain the population after optimizing.Wherein, renewal time Number m can be determined according to the actual requirements, and such as, m can take more than or equal to the 50, integer less than or equal to 100.

Before step S3 is discussed in detail, the glowworm swarm algorithm (Firefly Algorithm, i.e. FA) first to standard is entered Row is introduced.In glowworm swarm algorithm, the Attraction Degree that individuality is subject to determines moving direction and the distance of single individuality.By to this The renewal of two key elements, thus realizing seeking of optimal value is taken, its major parameter is defined as follows：Firefly is subject to Attraction Degree fixed Justice isThe location updating formula that each iteration is adopted is defined as：s_i=s_i×(1-β)+β×s_j+αε_i. Wherein, γ approximately represents the fixing absorption coefficient of light, its typical value scope (0.1,10), r_ijRepresent firefly away from another each and every one The distance of body, andBehalf firefly position coordinates.

During research glowworm swarm algorithm and particle filter, inventors have seen that：Can will be selected The firefly that particle analogy becomes to meet the worm cluster characteristics of motion is individual, using individual optimizing mode in population, makes weights relatively Low particle can move to higher particle, to improve the weight of itself, it is to avoid degradation phenomena.In consideration of it, the present invention Inventor with reference to the feature of particle filter algorithm, glowworm swarm algorithm is carried out with some corrections, and by revised glowworm swarm algorithm During particle filtering resampling.

Specifically, in this embodiment, the Attraction Degree formula redefining is as follows：

In formula, β_i,jFor particleTo particleAttraction Degree,For particle corresponding state variable norm.β_i,j Bigger, illustrate that the weights difference of two particles is bigger.Work as β_i,jDuring less than 0, particle in observation scope is describedWeights be less than ParticleNeed not be to particleIt is updated.

In addition, in step s3, the mobile formula of particle redefining is as follows:

In formula,It is to particleProduce maximum Attraction Degree β_{I, max}Particle, K is that often step updates step-length, and K=β_i,max. For avoiding glowworm swarm algorithm to be absorbed in local extremum phenomenon as other intelligent optimization algorithms, introduce α ε_i, wherein, α is random Step factor, general span is [0,1]；ε_iBe usually by Gaussian Profile, be uniformly distributed or other distribution produce random Number vector.

After Attraction Degree, particle movement are redefined, based on step S3, the population in k moment can be optimized. In the specific implementation, first threshold β_thCan be determined according to actual conditions.Such as, when initial population is 100, and finally effective When population is 10～20, can be by β_thIt is defined as：β_th=(rand (1)+1)/N.In formula, rand (1) is the random number within 1. Can be seen that the actual numerical value magnitude of Attraction Degree by observing Attraction Degree computing formula 10^-2, even more little.By above formula definition β_th, ensure that most particles move to more excellent position.In addition, the weights of the particle after to location updating are counted During calculation, the probability density function that this particle can be directly substituted into existing white noise is calculated.

Step S4, according to optimize after population estimating system state parameter the k moment state averageAnd, root According to optimize after population estimating system state parameter the k moment covariance P_k；

Specifically, in step s 4, calculated according to formula 6P is calculated according to formula 7_k；

In formula,For system status parameters the k moment state average.

In formula,ForTransposition.

In the specific implementation, step S2～S4 can repeatedly be circulated according to demand.After primary particle filtering terminates, under The population that the primary particle filtering last time was processed is directly inputted in state equation, enters and circulates again.In the present invention In embodiment, by being optimized to the resampling process of particle filter using glowworm swarm algorithm, not only ensure that number of particles, And ensure that the diversity of particle.Further, improve the precision of orbital position, prediction of speed.

With reference to Fig. 2 to Fig. 5 and a concrete simulation process, the technique effect of the embodiment of the present invention is carried out specifically Bright.In this simulation process, emulation data is as follows：In platform track parameter, semi-major axis of orbit is 17178.1km, centrifugation Rate is 0.001, and orbit inclination angle is 5 °, and the argument of perigee is 0 °, and the red footpath of ascending node is 0 °, and true anomaly is 100 °；In target track In road parameter, semi-major axis of orbit is 17128.1km, and eccentricity is 0.012, and orbit inclination angle is 10 °, and the argument of perigee is 30 °, rises The red footpath of intersection point is 0 °, and true anomaly is 70 °.Simulation time chooses 200S, and initial measured error variance matrix isSimulation result is referring to Fig. 2 to Fig. 5.From Fig. 2 to Fig. 5, and existing Particle filter orbital prediction method compare, the present invention improves the precision of position, prediction of speed to a great extent.

Although with reference to illustrative embodiments, invention has been described but it is to be understood that the present invention does not limit to The specific embodiment that Yu Wenzhong describes in detail and illustrates, in the case of without departing from claims limited range, this Skilled person can make various changes to described illustrative embodiments.

Claims (9)

1. a kind of PF space non-cooperative target orbital prediction method being optimized based on firefly group is it is characterised in that methods described Including：

S1, the initial state value construction primary group according to systemAnd by each primary's Weights are set to

S2, the state equation according to system's relative motionCalculate the population in k momentAnd each particle is calculated according to the measured value in this momentNormalization weights

S3, according to glowworm swarm algorithm, the population in k moment is optimized, including：In each renewal process, calculate each grain SonAttraction Degree β to other particles_i,j, j=1,2 ... N and j ≠ i, and from β_i,jMiddle selection particleMaximum Attraction Degree β_i,max；When maximum Attraction Degree β_i,maxMore than default first threshold β_thWhen, this particle is carried out with location updating and calculates its power Value；After completing default update times m, obtain the population after optimizing；

S4, according to optimize after population estimating system state parameter the k moment state average

2. the method for claim 1, wherein in step s3, β is calculated according to formula 1_i,j；

In formula, β_i,jFor particleBy particleAttraction Degree,For particle corresponding state variable norm.

3. the method for claim 1, wherein in step s3, location updating is carried out to particle according to formula 2；

In formula,It is to particleProduce maximum Attraction Degree β_i,maxParticle, ε_iIt is the random number vector being produced by Gaussian Profile, K is that often step updates step-length, and K=β_i,max, α is the arbitrary width factor.

4. the method for claim 1, wherein in step s 2, each particle is calculated according to formula 3,4,5Normalizing Change weights

In formula,For important density function,Particle for the k momentState value take observation y_k's Probability,Weights for the k moment after updating.

5. the method for claim 1, wherein in step s 4, calculated according to formula 6

In formula,For system status parameters the k moment state average.

6. the method for claim 1, wherein step S4 also includes：According to the population estimating system state after optimizing Parameter is in covariance P in k moment_k；

In formula,ForTransposition.

7. described method as arbitrary in claim 1-6, wherein, first threshold β_thMeet：

β_th=(rand (1)+1)/N；

In formula, rand (1) is the random number within 1.

8. described method as arbitrary in claim 1-6, wherein, N meets:100≤N≤500.

9. described method as arbitrary in claim 1-6, wherein, default update times m meet：50≤m≤100.