CN104833949A - Multiple-unmanned aerial vehicle cooperative passive location method based on improved distance parameterization - Google Patents

Abstract

The invention provides a multiple-unmanned aerial vehicle cooperative passive location method based on improved distance parameterization. First, an initial distance interval is divided, and the weight of the distance interval is initialized; then, squared root cubature information filtering is performed separately on the sub intervals, and the weight of the distance interval is updated; and weighted fusion is performed on filtering results of all the sub intervals to obtain a final location result. The problem that the covariance matrix is not positive definite due to round-off error of a computer is eliminated, and the robustness of the filtering method is improved. Compared with a traditional distance parameterization method, the amount of calculation of the adopted improved distance parameterization method is greatly reduced. The influence of the initial filtering value to the location performance is avoided effectively, and the real-time performance of the filtering method is improved. Information filtering under a fully distributed fusion estimation structure makes fusion estimation calculation easier and improves the location accuracy.

Description

A kind of multiple no-manned plane based on improving distance parameter works in coordination with passive location method

Technical field

The invention belongs to passive location field, relate to a kind of multiple no-manned plane and work in coordination with passive location method.

Background technology

Passive location refers to the radiation detection target utilizing target self, and determines a kind of technology of target state.Compared to active locations such as radars, passive location, due to the non-radiating signal of observer own, has the advantages such as disguise is strong, detection range is far away, is widely used in the fields such as electronic reconnaissance, electronic countermeasure and precise guidance.

Observed quantity due to passive sensor only has the angle information of target, cannot obtain the range information of target, when taking single unmanned plane passive location, in some cases, even if carrier aircraft self is motor-driven, still can not ensure the complete observability to target.Multiple no-manned plane is worked in coordination with passive location and is realized information sharing by communication between machine, adopts information fusion technology, can improve the observability to target, become the study hotspot of current passive location.

In passive location, the observed quantity due to sensor is all the nonlinear function of quantity of state, so need the matter of utmost importance solved to be that a kind of estimated accuracy of searching is high, the nonlinear filtering algorithm of better numerical value stability.Nonlinear function is carried out first order Taylor expansion by EKF (EKF), ignores higher order term, by nonlinear problem linearization.But for nonlinearity problem, use EKF can produce larger error, even cause filtering divergence, and EKF needs when linearization to calculate Jacobian matrix, and be difficult to obtain in practical problems.For the deficiency of EKF, Julier proposes Unscented kalman filtering (UKF), the weight sampling point determined by one group approaches the distribution function of stochastic variable, by the nonlinear transformation of sampled point, catches the statistical property of stochastic variable after nonlinear transformation.But UKF is carrying out needing in iterative filtering process to carry out matrix decomposition and inversion operation, and state estimation covariance battle array is difficult to keep positive definite.The volume Kalman filtering (CKF) proposed in recent years, uses the numerical integration method based on volume principle to calculate the stochastic variable average after nonlinear transformation and covariance.But CKF needs to decompose error co-variance matrix and invert in iterative filtering process, and the non-positive definite of error covariance matrix can cause the instability of filtering numerical value.

In addition, in passive location, filtering initial value is larger on locating effect impact.Because passive sensor cannot obtain the range information of target, make filtering initial distance and true initial distance there is comparatively big error, and then have a strong impact on positioning precision.Peach proposes a kind of method of distance parameter, namely at the initial time of filtering, by the supposition of target range information in some sub-ranges, and distribute initial weight to each sub-range, the independent parallel filtering of each sub-range also upgrades weights, finally using the filter result Weighted Fusion in each sub-range as final state estimation.But the method exists a large amount of computing when solving sub-range likelihood function, system real time is poor and easily cause wave filter to lose efficacy because of index exploding.

Summary of the invention

In order to overcome the deficiencies in the prior art, the invention provides a kind of multiple no-manned plane based on improving distance parameter and working in coordination with passive location method.

The technical solution adopted for the present invention to solve the technical problems is: first improve traditional distance parameter method, choose numerical stability good square root volume Kalman filtering (SRCKF) algorithm, the subduplicate mark of utilization state covariance matrix carries out adaptive updates to the weights in each sub-range, obtains distance parameter square root volume Kalman filtering (RP-SRCKF) algorithm.On this basis, combining information filtering, adopts complete distributed fusion to estimate structure, improves RP-SRCKF algorithm, obtains distance parameter square root volume information filtering (RP-SRCIF) algorithm under complete distributed fusion estimation structure.Specifically comprise the following steps:

Initial distance interval division is N by step 1, employing geometric ratio rule _findividual sub-range, a jth sub-range is [R _minρ ^j-1, R _minρ ^j], 1≤j≤N _f, R _minfor the lower limit in initial distance interval, R _maxfor the higher limit in initial distance interval, for scale factor, then the distance range in a jth sub-range is R _j=R _min(ρ ^j-ρ ^j-1);

Step 2, initialization distance interval right weight, distribute initial weight to a jth sub-range

Step 3, the independent operating square root volume information filtering of each sub-range, comprise the following steps:

The original state in 3.1 given each sub-ranges and error co-variance matrix obtain the square root of the initial information matrix in each sub-range and information state subscript j represents a jth sub-range, and the volume basic point ξ of calculating sampling _i, i=1 ..., 2n, n are the dimension of system state vector;

3.2 computational prediction states with the square root of predicting covariance matrix and computational prediction information state thus with the square root of predicted state information matrix subscript j represents a jth sub-range, and subscript k represents the kth moment;

3.3 computing information matrix contribution square roots contribute with information state and the square root of lastest imformation matrix with information state vector finally obtain system state vector subscript j represents a jth sub-range, and subscript d represents d sensor, and subscript k represents the kth moment;

Step 4, according to the subduplicate mark of each sub-range state covariance matrix, real-time adaptive ground regulates the weights in each sub-range, the weights in a kth moment jth sub-range for a jth subduplicate mark of sub-range error co-variance matrix;

Step 5, fusion is weighted to the filter result in each sub-range, obtains final positioning result ${\hat{x}}_{k|k}=\underset{j=1}{\overset{N_F}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^j{\hat{x}}_{k|k}^j;$

Step 6, compare sub-range weights with the metric-threshold ω of setting _thsize, if then stop wave filter corresponding to a jth sub-range, otherwise directly enter step 7;

Step 7, judge whether filtering algorithm reaches the execution number of times of setting, if so, terminate algorithm; Otherwise, return step 3.2.

The invention has the beneficial effects as follows: eliminate by introducing square root filtering technology the covariance matrix not positive definite problem brought by computing machine round-off error, improve the robustness of filtering method, the relatively traditional distance parameter method of improvement distance parameter method adopted significantly decreases operand, avoiding filtering initial value to the impact of positioning performance improves the real-time of filtering method, complete distributed fusion estimates that the information filter under structure makes fusion estimation calculating become simple, and improves positioning precision.

Accompanying drawing explanation

Fig. 1 is the process flow diagram of implementation method of the present invention;

Fig. 2 is the movement locus figure of target and unmanned plane;

Fig. 3 is the accuracy comparison figure of the present invention and other localization methods, and wherein, (a) is the location estimation relative error comparison diagram of the improvement RP-SRCKF algorithm that proposes of the present invention and traditional RP-SRCKF algorithm.B () is the location estimation relative error comparison diagram of the multi-machine collaborative passive location method that proposes of the present invention and unit passive location method.

Embodiment

Below in conjunction with drawings and Examples, the present invention is further described, the present invention includes but be not limited only to following embodiment.

Setting non-linear system status spatial model:

$\left\{\begin{array}{c}{X}_k=f\left(X_{k-1}\right)+{\mathrm{\ω}}_{k-1}\\ Z_k=h\left(X_k\right)+v_k\end{array}\right.$

In formula, X _k∈ R ⁿit is the system state in k moment; Z _k∈ R ^mit is the measurement in k moment; ω _k-1~ N (0, Q _k-1) be the process noise in k-1 moment; v _k~ N (0, R _k) be the measurement noise in k-1 moment; F and h is respectively state model function and measurement model function.

Based on above-mentioned nonlinear system model, under complete distributed fusion estimation structure, the detailed step of RP-SRCIF algorithm is as follows:

Step 1 divides initial distance interval

Geometric ratio rule is adopted to be N by initial distance interval division _findividual sub-range, a jth sub-range is [R _minρ ^j-1, R _minρ ^j], 1≤j≤N _f, R _minfor the lower limit in initial distance interval, R _maxfor the higher limit in initial distance interval, for scale factor, then the distance range in a jth sub-range is R _j=R _min(ρ ^j-ρ ^j-1);

Step 2 initialization distance interval right weight

Distribute initial weight to each sub-range by zone distance coverage:

${\mathrm{\ω}}_0^j=\frac{R_j}{R_{\max }-R_{\min }}$

The independent operating square root volume information filtering of step 3 each sub-range

3.1 sub-range filtering initialization

3.1.1 the original state in given sub-range and error co-variance matrix obtain initial information matrix and the information state in sub-range:

$\left\{\begin{array}{c}{Y}_{0|0}^j={\left(P_{0|0}^j\right)}^{-1}\\ {\hat{y}}_{0|0}^j={\left(P_{0|0}^j\right)}^{-1}{\hat{X}}_{0|0}^j\end{array}\right.$

In formula, subscript j represents a jth sub-range.

3.1.2 right carry out matrix decomposition, obtain the square root of information matrix:

$Y_{0|0}^j=S_{y,0|0}^j{\left(S_{y,0|0}^j\right)}^T$

3.1.3 the volume basic point ξ of calculating sampling _i:

${\mathrm{\ξ}}_i=\sqrt{n}{\left[1\right]}_i$

In formula, i=1 ..., 2n, n are the dimension of system state, [1] _ibe defined as follows:

${\left[1\right]}_i=\left\{\begin{array}{cccccc}\left(\begin{array}{c}1\\ 0\\ .\\ .\\ .\\ 0\end{array}\right),& \left(\begin{array}{c}0\\ 1\\ .\\ .\\ .\\ 0\end{array}\right),...,& \left(\begin{array}{c}0\\ 0\\ .\\ .\\ .\\ .\end{array}\right),& \left(\begin{array}{c}-1\\ 0\\ .\\ .\\ .\\ 0\end{array}\right),& \left(\begin{array}{c}0\\ -1\\ .\\ .\\ .\\ 0\end{array}\right),...,& \left(\begin{array}{c}0\\ 0\\ .\\ .\\ .\\ -1\end{array}\right)\end{array}\right\}$

3.2 times upgraded

3.2.1 computing system state covariance matrix square root:

$S_{k-1|k-1}^j={\left(S_{y,k-1|k-1}^j\right)}^{-T}$

3.2.2 computing mode volume point:

${\mathrm{\χ}}_{i,k-1|k-1}^j=S_{k-1|k-1}^j\·{\mathrm{\ξ}}_i+{\hat{X}}_{k-1|k-1}^j$

3.2.3 spread state volume point:

${\mathrm{\χ}}_{i,k|k-1}^j=f\left({\mathrm{\χ}}_{i,k-1|k-1}^j\right)$

3.2.4 computational prediction state:

${\widehat{\mathrm{\χ}}}_{k|k-1}^j=\frac{1}{2n}\underset{i=1}{\overset{2n}{\mathrm{\Σ}}}{\mathrm{\χ}}_{i,k|k-1}^j$

3.2.5 the square root of computational prediction error co-variance matrix:

$S_{k|k-1}^j=\mathrm{Tria}([\begin{array}{cc}{\mathrm{\χ}}_{k|k-1}^j& S_{Q,k-1}\end{array}])$

In formula, and S _{q, k-1}be defined as follows:

${\mathrm{\χ}}_{k|k-1}^j=\frac{1}{\sqrt{2n}}[{\mathrm{\χ}}_{1,k|k-1}^j-{\hat{X}}_{k|k-1}^j,{\mathrm{\χ}}_{2,k|k-1}^j-{\hat{X}}_{k|k-1}^j,...,{\mathrm{\χ}}_{2n,k|k-1}^j-{\hat{X}}_{k|k-1}^j]$

$Q_{k-1}=S_{Q,k-1}{S}_{Q,k-1}^T$

3.2.6 the square root of computational prediction state information matrix:

$S_{y,k|k-1}^j={\left(S_{k|k-1}^j\right)}^{-T}$

3.2.7 computational prediction information state:

${\hat{y}}_{k|k-1}^j=S_{y,k|k-1}^j{\left(S_{y,k|k-1}^j\right)}^T{\hat{X}}_{k|k-1}^j$

3.3 measure renewal

3.3.1 computational prediction state volume point, and upgrade

${\mathrm{\χ}}_{i,k|k-1}^j=S_{k|k-1}^j{\mathrm{\ξ}}_i+{\hat{X}}_{k|k-1}^j$

${\mathrm{\χ}}_{k|k-1}^j=\frac{1}{\sqrt{2n}}[{\mathrm{\χ}}_{1,k|k-1}^j-{\hat{X}}_{k|k-1}^j,{\mathrm{\χ}}_{2,k|k-1}^j-{\hat{X}}_{k|k-1}^j,...,{\mathrm{\χ}}_{2n,k|k-1}^j-{\hat{X}}_{k|k-1}^j]$

3.3.2 propagation forecast state volume point, the prediction obtaining each sensor measures volume point:

$Z_{d,i,k|k-1}^j=h_d\left({\mathrm{\χ}}_{i,k|k-1}^j\right)$

H in formula _dbe the measurement model of d sensor, subscript d represents d sensor.D=1 ..., D, D are the number of sensor.

3.3.3 computational prediction measures:

$Z_{d,k|k-1}^j=\frac{1}{2n}\underset{i=1}{\overset{2n}{\mathrm{\Σ}}}{Z}_{d,i,k|k-1}^j$

3.3.4 computational prediction state and prediction measure Cross-covariance:

$P_{d,\mathrm{xz},k|k-1}^j={\mathrm{\χ}}_{k|k-1}^j{\left[Z_{d,k|k-1}^{j*}\right]}^T$

In formula, be defined as follows:

$Z_{d,k|k-1}^{j*}=\frac{1}{\sqrt{2n}}[Z_{d,1,k|k-1}^j-Z_{d,k|k-1}^j,Z_{d,2,k|k-1}^j-Z_{d,k|k-1}^j,...,Z_{d,2n,k|k-1}^j-Z_{d,k|k-1}^j]$

3.3.5 computing information matrix contribution square root:

$S_{d,i,k}^j=S_{y,k|k-1}^j{\left(S_{y,k|k-1}^j\right)}^T{P}_{d,\mathrm{xz},k|k-1}^j{\stackrel{\‾}{S}}_{d,R,k}$

In formula, for

${\stackrel{\‾}{S}}_{d,R,k}=S_{d,R,k}^{-T}$

Wherein S _{d, R, k}it is the measurement noise covariance matrix square root of d sensor.

3.3.6 computing information state contribution:

$i_{d,k}^j=S_{d,i,k}^j{\stackrel{\‾}{S}}_{d,R,k}^T(Z_{d,k}-Z_{d,k|k-1}^j)+S_{d,i,k}^j{\left(S_{d,i,k}^j\right)}^T{\hat{X}}_{k|k-1}^j$

3.3.7 lastest imformation state vector:

${\hat{y}}_{k|k}^j={\hat{y}}_{k|k-1}^j+\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{i}_{d,k}^j$

3.3.8 the square root of lastest imformation matrix:

$S_{y,k|k}^j=\mathrm{Tria}([\begin{array}{ccccccc}{S}_{y,k|k-1}^j& S_{1,i,k}^j& S_{2,i,k}^j& .& .& .& S_{D,i,k}^j\end{array}])$

3.3.9 system state vector is upgraded:

${\hat{X}}_{k|k}^j={\left(S_{y,k|k}^j\right)}^{-T}{\left(S_{y,k|k}^j\right)}^{-1}{\hat{y}}_{k|k}^j$

Right value update between step 4 distance regions

According to the subduplicate mark of each sub-range state covariance matrix, real-time adaptive ground regulates the weights in each sub-range, then the weights in a kth moment jth sub-range are:

${\mathrm{\ω}}_k^j=\frac{\frac{1}{{\mathrm{trS}}_{k|k}^j}{\mathrm{\ω}}_{k-1}^j}{\underset{i=1}{\overset{N_F}{\mathrm{\Σ}}}\frac{1}{{\mathrm{trS}}_{k|k}^j}{\mathrm{\ω}}_{k-1}^j}$

In formula, for a jth subduplicate mark of sub-range error co-variance matrix;

The filter result of step 5 to each sub-range is weighted fusion, obtains final positioning result:

${\hat{X}}_{k|k}=\underset{j=1}{\overset{N_F}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^j{\hat{X}}_{k|k}^j$

Step 6 compares sub-range weights with the metric-threshold ω of setting _thsize, if then stop wave filter corresponding to a jth sub-range;

Step 7 judges whether filtering algorithm reaches the execution number of times of setting, if so, terminates algorithm; Otherwise, return step 3.2.

Effect of the present invention can be further illustrated by following emulation experiment:

Example for the same terrain object of two frame unmanned plane colocated carries out simulation analysis, suppose that two frame unmanned planes spiral at level altitude, UAV1 initial position is [0m, 4000m, 0m], initial velocity is 30m/s, and angle, initial heading is 0.25 π, continues to do at the uniform velocity turning motion at level altitude with angular velocity-5mrad/s; UAV2 initial position is [0m ,-4000m, 0m], and initial velocity is 30m/s, and angle, initial heading is-0.25 π, continues to do at the uniform velocity turning motion at level altitude with angular velocity 5mrad/s.Target initial position is [0m, 0m ,-1000m], and speed is 30m/s, makees S type motor-driven on ground.Process noise is zero mean Gaussian white noise, and each change in coordinate axis direction location criteria difference is 10m, and velocity standard difference is 0.1m/s.Target and unmanned plane movement locus are as shown in Figure 2.

Suppose that two frame unmanned planes all carry infrared sensor, its detection range interval is [1000m, 32000m]; get scale factor ρ=2, then can will be divided into 5 independently sub-ranges between distance regions, i.e. [1000m; 2000m], and [2000m, 4000m]; [4000m, 8000m], [8000m; 16000m]; [16000m, 32000m], setting metric-threshold ω _th=0.001.The unmanned plane tracking target time is 1000s, and unmanned plane observation cycle is T=1s, and observation noise is zero mean Gaussian white noise, and position angle and angle of pitch standard deviation are 1mrad.

In order to verify the performance of improvement distance parameter method proposed by the invention, for above-mentioned simulating scenes, choose filtering original state for [10000 ,-10000 ,-5000,0,20,0] ^t, through 100 Monte Carlo simulations, obtain location estimation relative error curve as shown in Fig. 3 (a).As can be seen from Fig. 3 (a), the filtering accuracy of the improvement RP-SRCKF algorithm that the present invention proposes and traditional RP-SRCKF algorithm quite but be all significantly higher than SRCKF algorithm.

Contrast the working time of these three kinds of filtering algorithms again, as shown in table 1.

Table 1 filtering algorithm working time

As can be seen from Table 1, adopt distance parameter to add working time, and the improvement RP-SRCKF algorithm that the present invention proposes is owing to greatly reducing calculated amount, has shorter working time, improve the real-time of system compared to traditional RP-SRCKF algorithm.

In order to verify the performance of Fusion Estimation Algorithm proposed by the invention, for above-mentioned simulating scenes, respectively simulation comparison being carried out to three kinds of algorithms, being respectively: the SRCKF algorithm of the RP-SRCIF algorithm of multi-machine collaborative, the SRCIF algorithm of multi-machine collaborative and single unmanned plane.Through 100 Monte Carlo simulations, obtain location estimation relative error curve as shown in Fig. 3 (b).As can be seen from Fig. 3 (b), multi-machine collaborative location has higher filtering accuracy compared to unit location, and when filtering initial distance differs larger with actual distance, because the RP-SRCIF algorithm of the present invention's proposition is by the impact of filtering initial distance, its filtering accuracy is significantly higher than all the other two kinds of algorithms.

It should be noted that in 400s ~ 600s time period, position estimation error is larger in Fig. 3 (a) and Fig. 3 (b).Analyzing its reason is because this time period internal object motion model changes, thus causes filtering accuracy poor relative to other times section.

Claims (1)

1. work in coordination with a passive location method based on the multiple no-manned plane improving distance parameter, it is characterized in that comprising the steps:

Initial distance interval division is N by step 1, employing geometric ratio rule _findividual sub-range, a jth sub-range is [R _minρ ^j-1, R _minρ ^j], 1≤j≤N _f, R _minfor the lower limit in initial distance interval, R _maxfor the higher limit in initial distance interval, for scale factor, then the distance range in a jth sub-range is R _j=R _min(ρ ^j-ρ ^j-1);

Step 2, initialization distance interval right weight, distribute initial weight to a jth sub-range

Step 3, the independent operating square root volume information filtering of each sub-range, comprise the following steps:

The original state in 3.1 given each sub-ranges and error co-variance matrix obtain the square root of the initial information matrix in each sub-range and information state subscript j represents a jth sub-range, and the volume basic point ξ of calculating sampling _i, i=1 ..., 2n, n are the dimension of system state vector;

3.2 computational prediction states with the square root of predicting covariance matrix and computational prediction information state thus with the square root of predicted state information matrix subscript j represents a jth sub-range, and subscript k represents the kth moment;

3.3 computing information matrix contribution square roots contribute with information state and the square root of lastest imformation matrix with information state vector finally obtain system state vector subscript j represents a jth sub-range, and subscript d represents d sensor, and subscript k represents the kth moment;

Step 4, according to the subduplicate mark of each sub-range state covariance matrix, real-time adaptive ground regulates the weights in each sub-range, the weights in a kth moment jth sub-range for a jth subduplicate mark of sub-range error co-variance matrix;

Step 5, fusion is weighted to the filter result in each sub-range, obtains final positioning result ${\hat{X}}_{k|k}=\underset{j=1}{\overset{N_F}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^j{\hat{X}}_{k|k}^j;$

Step 6, compare sub-range weights with the metric-threshold ω of setting _thsize, if then stop wave filter corresponding to a jth sub-range, otherwise directly enter step 7;

Step 7, judge whether filtering algorithm reaches the execution number of times of setting, if so, terminate algorithm; Otherwise, return step 3.2.