CN104316945A - Satellite interference source three-satellite positioning method based on high-order cumulants and unscented Kalman filtering - Google Patents

Abstract

The invention relates to a satellite interference source three-satellite positioning method based on high-order cumulants and unscented Kalman filtering. The effects of accurate estimation of time delay difference in a time-varying and non-Gaussian-noise environment and the stability and the real-time performance of a time difference positioning algorithm are mainly achieved. The positioning method includes the steps that (1), preprocessing such as zero mean and normalization is carried out on received satellite interference signals, and the high-order cumulants of the received satellite interference signals are calculated; (2), an improved time delay estimation method based on four-order cumulants minimum mean square errors is used for estimating the time delay difference generated by retransmitting the interference signals through three satellites; (3), according to estimation of the time delay difference, initial estimation of the position of the interference signals is obtained by using an analytical method, and then the positioning result of the previous moment is modified and updated by using an unscented Kalman filtering method. Compared with the prior art and method, a satellite interference source can be positioned in the time-varying and non-Gaussian-noise environment, and the positioning method is high in convergence rate, high in stability, easy to implement and small in calculation quantity.

Description

A kind of satellite interference source three star problem method based on Higher Order Cumulants and Unscented kalman filtering

One technical field

The present invention relates to a kind of Satellite tool kit method, particularly a kind of satellite interference source three star problem method based on Higher Order Cumulants and Unscented kalman filtering.

Two background technologies

Satellite communication has the features such as broad covered area, communication distance is long, communication line is stable, communication band is wide, capacity is large, has played indispensable effect in national defence and commercial communication field.But because satellite transit is on disclosed space orbit, even destroy so satellite system is easily interfered.Wherein, that have a mind to from ground or Radio frequency interference (RFI) is day by day serious to the threat of telstar unintentionally.No matter the elimination of which kind of interference, first all accurately must locate interference source, utilize interference source localization method, when not changing telstar duty, test, estimate the approximate location of ground interference source, coordinate other means more on this basis, effectively can determine the accurate location of ground interference source.So Satellite tool kit method is to ensureing the effective of satellite communication system and reliability service, play a very important role.

Double-Star Positioning System technology is the main method of the not clear interference source in location at present, and it is by realizing interference source location to the estimation of time of arrival poor (TDOA) and arrival rate difference (FDOA).But there is the Accurate Measurement difficulty of frequency shift parameters in Double-Star Positioning System technology, and affects seriously by factors such as attitude of satellite changes, generally needs to set up multiple reference station to carry out error correction, the positioning precision that just can be comparatively satisfied with.

Samsung time difference positioning method utilizes the satellite group of space three satellite compositions in a distance, the signal of measurement ground interference source propagates into the mistiming between two between satellite, calculate two hyperboloids formed thus, and utilize the intersection point of itself and interference source place earth sphere, draw interference source position.Different with Double-Star Positioning System algorithm, this is that one utilizes two with reference to satellites, the method for being located by twice calculation delay, only make use of delay parameter, relative to frequency shift parameters, the mensuration of delay parameter is comparatively easily with accurate, finally higher to the positioning precision of interference source.

The performance of Samsung time difference positioning method depends on the estimated accuracy to difference time of arrival.Mainly contain about delay time estimation method: generalized correlation method, maximum likelihood estimate, least-mean-square error algorithm and high-order statistic method etc.Wherein, generalized correlation method is early stage delay time estimation method, and its principle is simple, calculated amount is also very little, but the precision estimated is not high, so more in the applications that background condition meets and accuracy requirement is not too high.Maximum likelihood estimate generally needs the probability density understanding signal in advance, and in practice, and the probability distribution of signal is difficult to obtain often, and maximum likelihood estimate has very large limitation when applying.Least mean-square error Time Delay Estimation Method is widely applicable, and calculated amount is little, and its real-time is better, practical, but they can not suppress the impact of correlation noise at all.The Higher Order Cumulants perseverance of Gaussian process is zero, therefore high-order statistic method is very effective when estimating the signal time delay difference with Gaussian Correlative Noise, but in practical application, the noise being mingled in signal may be non-gaussian Complex Noise, and signal to noise ratio (S/N ratio) also may change in time, during high-order statistic method pair, when change, non-Gaussian noise, satisfactory performance can not be obtained.Therefore, time delay estimation problem under change when existing delay parameter method of estimation needs to solve, non-Gaussian noise situation.

On the other hand, Localization Estimate Algorithm of TDOA can be summed up as the Solve problems of Nonlinear System of Equations.It solves positioning equation according to one group of time difference of measuring, and solves position of interference source exactly from geometrically saying by two crossing time difference hyperboloids and ground spheres intersect, its essence is to combine three ternary quadratic equations to solve.At present, that commonly uses in Localization Estimate Algorithm of TDOA has analytical method, Newton iteration method, least-squares iteration method, particle filter etc.Analytical method a kind ofly unknown number is directly calculated to the method solved, and generally has two solutions, can produce location ambiguity.Newton iterative relates to the calculating of trigonometric function, and calculated amount is larger.Least-squares iteration legal position stability is poor, and not easily restrains.Particle filter needs to adopt a large amount of particle samplers to go the error of approaching to reality to distribute, and calculated amount is large, and therefore real-time is poor.Therefore, need to improve existing Localization Estimate Algorithm of TDOA to solve the problem such as stability, real-time.

Three summary of the invention

The object of the invention is the problem such as stability, real-time in order to solve time-vary delay system tracking and Localization Estimate Algorithm of TDOA in Satellite tool kit delay parameter method of estimation, estimate to combine with Unscented kalman filtering by based on the fourth order cumulant least mean-square error time delay improved, propose a kind of satellite interference source three star problem method based on Higher Order Cumulants and Unscented kalman filtering.

Realizing step of the present invention is:

The first, the pre-service such as zero-mean, normalization is carried out to the satellite interference signal received in ground receiving station, and calculates its Higher Order Cumulants;

The second, utilize improve based on fourth order cumulant least mean-square error delay time estimation method, estimate undesired signal and forward the delay inequality produced through three satellites;

3rd, according to the estimation of delay inequality, utilize analytical method to obtain the initial estimation of target location, then utilized the positioning result of Unscented kalman filtering method to a upper moment to carry out correction and upgrade.

Beneficial effect of the present invention is: 1, have employed improvement based on fourth order cumulant least mean-square error delay time estimation method, the method effectively can suppress Gaussian noise and can accurately estimate time-vary delay system parameter; 2, utilize Unscented kalman filtering method constantly to revise renewal to positioning result, the method fast convergence rate, stability are high, and realize simple, and calculated amount is little.

Four accompanying drawing explanations

Fig. 1 is satellite interference source Samsung positioning using TDOA schematic diagram;

Fig. 2 be improve based on fourth order cumulant least mean-square error delay time estimation method theory diagram;

Fig. 3 is the positioning result correction renewal process flow diagram based on Unscented kalman filtering.

Five embodiments

Below in conjunction with accompanying drawing, the invention will be further described.

1, Fig. 1 gives satellite interference source Samsung positioning using TDOA system schematic.The undesired signal that ground interference source is launched, through the forwarding of different satellite, receives by processing enter is unified.First processing enter carries out the pre-service such as zero-mean, normalization to signal, then the mistiming parameter that three tunnel undesired signals arrive is estimated, recycling analytical method obtains the initial estimation of target location, finally utilizes the positioning result of Unscented kalman filtering method to a upper moment to carry out correction and upgrades.

2, Fig. 2 give improvement based on fourth order cumulant least mean-square error delay time estimation method theory diagram.

The model that undesired signal time delay is estimated can be expressed as follows:

$\left.\begin{array}{c}x\left(k\right)=s\left(k\right)+w_1\left(k\right)\\ y\left(k\right)=\mathrm{As}(k-D)+w_2\left(k\right)\end{array}\right\}---\left(1\right)$

Wherein k is discrete-time variable, and D is time delay, and A is decay factor, and x (k), y (k) are respectively to be had different delay and has superposed the undesired signal of noise, and s (k) is undesired signal, w ₁(k) and w ₂k () is respectively ground unrest, and suppose that it is uncorrelated, zero-mean independent with undesired signal mutually.

According to convolution theory, inhibit signal s (k-D) can be represented as

$s(k-D)=\underset{n=-\∞}{\overset{\∞}{\mathrm{\Σ}}}\sin c(n-D)s(k-n)---\left(2\right)$

Wherein sinc function is the even function being axis of symmetry with it to the maximum, and meets

$\underset{n\→\∞}{\lim }\sin c\left(n\right)=0---\left(3\right)$

Namely, when time variable n is tending towards infinite, sinc function decays to 0.In actual applications, usually n is limited in a suitable scope, namely can is inhibit signal approximate representation

$s(k-D)=\underset{n=-p}{\overset{p}{\mathrm{\Σ}}}\sin c(n-D)s(k-n)---\left(4\right)$

Replace infinite with an abundant large integer p, and make p > | D|.

So inhibit signal model can be changed into

$\begin{array}{c}y\left(k\right)=\mathrm{As}(k-D)+w_2\left(k\right)\\ =\underset{i=-p}{\overset{p}{\mathrm{\Σ}}}A\sin c(i-D)s(k-i)+w_2\left(k\right)\\ =\underset{i=-p}{\overset{p}{\mathrm{\Σ}}}{b}_i[x(k-i)-w_1(k-i)]+w_2\left(k\right)\end{array}---\left(5\right)$

Wherein,

b _i＝Asinc(i-D) (6)

For zero-mean stationary stochastic process, the as follows from fourth order cumulant of stochastic variable x (k) can be defined:

C _xxxx(τ，0，0)＝cum[x(k)，x(k+τ)，x(k)，x(k)] (7)

Definition stochastic variable x (k), the Cross fourth order cumulant of y (k) is as follows:

C _xyxx(τ，0，0)＝cum[x(k)，y(k+τ)，x(k)，x(k)] (8)

Formula (5) is substituted into (8), can obtain:

$\begin{array}{c}{C}_{\mathrm{xyxx}}(\mathrm{\τ},0,0)=\mathrm{cum}[x\left(k\right),\underset{i=-p}{\overset{p}{\mathrm{\Σ}}}{b}_i[x(k+\mathrm{\τ}-i)-w_1(k+\mathrm{\τ}-i)]+w_2(k+\mathrm{\τ}),x\left(k\right),x\left(k\right)]\\ =\mathrm{cum}[x\left(k\right),\underset{i=-p}{\overset{p}{\mathrm{\Σ}}}{b}_i[x(k+\mathrm{\τ}-i)],x\left(k\right),x\left(k\right)]\\ +\mathrm{cum}[x\left(k\right),\underset{i=-p}{\overset{p}{\mathrm{\Σ}}}{b}_i[-w_1(k+\mathrm{\τ}-i)],x\left(k\right),x\left(k\right)]\\ +\mathrm{cum}[x\left(k\right),w_2(k+\mathrm{\τ}),x\left(k\right),x\left(k\right)]\end{array}---\left(9\right)$

And the fourth order cumulant of Gaussian noise is zero, then

$\begin{array}{c}{C}_{\mathrm{xyxx}}(\mathrm{\τ},\mathrm{0,0})=\mathrm{cum}[x\left(k\right),\underset{i=-p}{\overset{p}{\mathrm{\Σ}}}{b}_i[x(k+\mathrm{\τ}-i)],x\left(k\right),x\left(k\right)]\\ =\underset{i=-p}{\overset{p}{\mathrm{\Σ}}}{b}_i\mathrm{cum}[x\left(k\right),x(k+\mathrm{\τ}-i),x\left(k\right),x\left(k\right)]\\ =\underset{i=-p}{\overset{p}{\mathrm{\Σ}}}{b}_i{C}_{\mathrm{xxxx}}(\mathrm{\τ}-i,\mathrm{0,0})\end{array}---\left(10\right)$

Weight coefficient vector is made to be:

W＝[b _-p，b _-p+1，…，b _p] ^T (11)

Cross fourth order cumulant vector is:

C _xy＝[C _xyxx(-p，0，0)，C _xyxx(-p+1，0，0)，…，C _xyxx(p，0，0)] ^T (12)

From fourth order cumulant vector be:

Then formula (10) can be expressed as with vector form:

C _xy＝C _4x·W (14)

Wherein C _4xfor (2p+1) × (2p+1) dimension, W is that (2p+1) × 1 is tieed up, C _xyfor (2p+1) × 1 is tieed up.

Input signal x (k) is through calculating C from fourth order cumulant _xxxx(τ, 0,0), then obtain outputing signal C by auto-adaptive fir filter _xzxx(τ, 0,0), and the C calculated with the Cross fourth order cumulant of x (k) and y (k) _xyxx(τ, 0,0) subtracts each other, and obtains error e (k).

Definition adaptive error function for

$\begin{array}{c}J\left(\hat{W}\right)=\underset{\mathrm{\τ}=-p}{\overset{p}{\mathrm{\Σ}}}{[C_{\mathrm{xzxx}}(\mathrm{\τ},\mathrm{0,0})-C_{\mathrm{xyxx}}(\mathrm{\τ},\mathrm{0,0})]}^2\\ =\underset{\mathrm{\τ}=-p}{\overset{p}{\mathrm{\Σ}}}{[\underset{i=-p}{\overset{p}{\mathrm{\Σ}}}{\hat{b}}_i{C}_{\mathrm{xxxx}}(\mathrm{\τ}-i,\mathrm{0,0})-C_{\mathrm{xyxx}}(\mathrm{\τ},\mathrm{0,0})]}^2\\ ={(C_{4x}\·\hat{W}-C_{\mathrm{xy}})}^T(C_{4x}\·\hat{W}-C_{\mathrm{xy}})\end{array}---\left(15\right)$

Wherein, weight coefficient vector is $\hat{W}={[{\hat{b}}_{-p},{\hat{b}}_{-p+1},...,{\hat{b}}_p]}^T.$

The gradient of error function is:

$\frac{\∂J\left(\hat{W}\right)}{\∂\hat{W}}=2C_{4x}^T(C_{4x}\·\hat{W}-C_{\mathrm{xy}})---\left(16\right)$

Adopt steepest descent method to realize the iterative process of optimized parameter vector in adaptive process.Make initial weight vector each numerical value changing current vector with the increment being proportional to negative gradient vector, namely

$\hat{W}(k+1)=\hat{W}\left(k\right)-\mathrm{\μ}{\left[\frac{\∂J\left(\hat{W}\right)}{\∂\hat{W}}\right]}_{\hat{W}=\hat{W}\left(k\right)}---\left(17\right)$

Wherein, μ is converging factor, is used for regulating adaptive stability and speed.

Formula (17) represents with scalar, then

${\hat{b}}_i(k+1)={\hat{b}}_i\left(k\right)-2\mathrm{\μ}\underset{i=-p}{\overset{p}{\mathrm{\Σ}}}{C}_{\mathrm{xxxx}}(\mathrm{\τ}-i,\mathrm{0,0})e\left(k\right)---\left(18\right)$

Wherein

$e\left(k\right)=\underset{i=-p}{\overset{p}{\mathrm{\Σ}}}{\hat{b}}_i\left(k\right)C_{\mathrm{xxxx}}(\mathrm{\τ}-i,\mathrm{0,0})-C_{\mathrm{xyxx}}(\mathrm{\τ}-i,\mathrm{0,0})---\left(19\right)$

In adaptive process, make FIR filter coefficient close to the time delays parameter of channel, namely make error function minimum, in last hunting time delay parameter estimated value maximum value, thus indirect estimation time delay D.

3, Fig. 3 gives and carries out correction renewal process flow diagram based on Unscented kalman filtering method to positioning result

Unscented kalman filtering is a kind of filtering based on deterministic sampling.It avoids by the method for deterministic sampling carries out linearization by nonlinear equation, and can obtain and be accurate to the average suitable with three rank taylor series expansions and variance, performance exceedes EKF, and calculated amount increases few.

Its step is as follows:

1. initialization

$\hat{X}\left(0\right)=E\left[X\left(0\right)\right]---\left(20\right)$

$P_0=E\left[(X\left(0\right)-\hat{X}\left(0\right))\right){(X\left(0\right)-\widehat{X\left(0\right)})}^T]---\left(21\right)$

In formula, be respectively interference position with X (0) and estimate initial value and actual position initial value, P ₀for error initial value.

2. Sigma value and weighting coefficient is calculated

$\left.\begin{array}{cc}{\mathrm{\χ}}_{0,k}=\stackrel{\‾}{X}\left(k\right)& \\ {\mathrm{\χ}}_{i,k}=\stackrel{\‾}{X}\left(k\right)+{\left(\sqrt{(L+\mathrm{\λ})P_k}\right)}_i& i=\mathrm{1,2},...,L\\ {\mathrm{\χ}}_{i,k}=\stackrel{\‾}{X}\left(k\right)-{\left(\sqrt{(L+\mathrm{\λ})P_k}\right)}_{i-L}& i=L+1,L+2,...,2L\end{array}\right\}---\left(22\right)$

$\left.\begin{array}{cc}{W_0}^{\left(m\right)}=\frac{\mathrm{\λ}}{L+\mathrm{\λ}}& \\ {W_0}^{\left(c\right)}=\frac{\mathrm{\λ}}{L+\mathrm{\λ}}+1-{\mathrm{\γ}}^2+\mathrm{\β}& i=\mathrm{1,2},...,2L\\ {W_i}^{\left(m\right)}={W_i}^{\left(c\right)}=\frac{1}{2(L+\mathrm{\λ})}& \end{array}\right\}---\left(23\right)$

In formula, χ is Sigma value, W weighting coefficient, for the average of location variable X (k), L is the quantity of Sigma point; λ=γ ²(L+ κ)-L is a compositely proportional factor; Constant γ determines that Sigma point is at mean value spreading range, it often get one very little on the occasion of, as 0≤γ≤10 ^-4; Constant κ is also a scale factor, and normal value is 3-L; β with distribution relevant, when for Gaussian distribution, β=2; for the i-th row of On Square-Rooting Matrices.

3. the time upgrades

$\left.\begin{array}{cc}{\mathrm{\χ}}_{i,k|k-1}=\mathrm{\Φ}{\mathrm{\χ}}_{i,k-1}& i=\mathrm{0,1,2},...,2L\\ {\hat{x}}_k^{-}=\underset{i=0}{\overset{2L}{\mathrm{\Σ}}}{W_i}^{\left(m\right)}{\mathrm{\χ}}_{i,k|k-1}& \end{array}\right\}---\left(24\right)$

$P_k^{-}=\underset{i=0}{\overset{2L}{\mathrm{\Σ}}}{W_i}^{\left(c\right)}({\mathrm{\χ}}_{i,k|k-1}-{\hat{x}}_k^{-}){({\mathrm{\χ}}_{i,k|k-1}-{\hat{x}}_k^{-})}^T+Q(k-1)---\left(25\right)$

In formula, Φ is state-transition matrix, for error matrix, Q is noise covariance matrix.

4. measurement updaue

In formula, h (.) is for measuring transforming function transformation function.

In formula, R is state-noise covariance matrix.

${\mathrm{\κ}}_k=P_{x_k{y}_k}{P}_{{\tilde{y}}_k{\tilde{y}}_k}^{-1}---\left(29\right)$

${\hat{x}}_k={\hat{x}}_k^{-}+{\mathrm{\κ}}_k(y_k-{\hat{y}}_y^{-})---\left(30\right)$

$P_k=P_k^{-}-{\mathrm{\κ}}_k{P}_{{\tilde{y}}_k{\tilde{y}}_k}^{-1}{{\mathrm{\κ}}_k}^T---\left(31\right)$

Claims (3)

1., based on a satellite interference source three star problem method for Higher Order Cumulants and Unscented kalman filtering, the steps include: the first, the pre-service such as zero-mean, normalization is carried out to the satellite interference signal received, and calculate its Higher Order Cumulants; The second, utilize improve based on fourth order cumulant least mean-square error delay time estimation method, estimate undesired signal and forward the delay inequality produced through three satellites; 3rd, according to the estimation of delay inequality, utilize analytical method to obtain the initial estimation of undesired signal position, then utilized the positioning result of Unscented kalman filtering method to a upper moment to carry out correction and upgrade.

2. method according to claim 1, it is characterized in that: propose improvement based on fourth order cumulant least mean-square error delay time estimation method, utilize two-way undesired signal from fourth order cumulant and Cross fourth order cumulant to build error function, utilize steepest descent method to find the minimum value of error function, make filter coefficient close to the delay parameter of channel, finally find the maximum value in delay parameter estimated value.

3. method according to claim 1, it is characterized in that: propose the positioning result correction update method based on Unscented kalman filtering device, the undesired signal position utilizing analytical method to calculate is as the initial value of Unscented kalman filtering, the time adopting the method for deterministic sampling to complete filtering upgrades and measurement updaue, thus constantly carries out correction renewal to the positioning result in a upper moment.