CN101753513A - Doppler frequency and phase estimation method based on polynomial forecasting model - Google Patents

Abstract

The invention belongs to the technical field of wireless communication and satellite navigation receivers, in particular provides a Doppler frequency and phase estimation method based on a polynomial forecasting model. The invention firstly provides a novel dynamic model, the polynomial forecasting model, for describing Doppler frequency and phase, combines with an odorless Kalman filter on the basis of the model, and finally provides a novel self-validating filtering algorithm for estimating the Doppler frequency and phase. No matter how the relative motion among the receivers changes, as long as the expression of the relative motion meets a sectioned polynomial form, the method provided in the invention can effectively estimate the Doppler frequency and phase. Numerical value simulation result and true GPS receiver experimental data all show that the method provided in the invention is superior to most of the Doppler frequency estimation methods reported in other literatures under the signal-to-noise ratio condition that the GPS receivers can work normally.

Description

A kind of Doppler frequency and phase estimation method based on polynomial forecast model

Technical field

The invention belongs to radio communication and satellite navigation receiver technical field, be specifically related to a kind of Doppler frequency and phase estimation method based on polynomial forecast model.

Background technology

Because the existence of Doppler effect makes the transmitting-receiving carrier wave of communication system may have bigger frequency deviation, causes receiver can't demodulate correct emission data ^[1]Be in the spread spectrum satellite communication of representative with global positioning system (GPS) particularly, carrier tracking loop and code phase track loop have constituted most important two parts of receiver together ^[2,3]

In the past few decades, estimation and tracking at Doppler frequency have proposed a large amount of methods in the document.Phase-lock-loop algorithm ^[2-5]Simple in structure owing to having, need the few characteristics of priori to be widely adopted, be a kind of maximal possibility estimation to doppler phase.Another kind of based on model, Doppler frequency algorithm for estimating optimum under the least mean-square error meaning comprises expanded Kalman filtration algorithm ^[6], the odorlessness Kalman filtering algorithm ^[7-9]And particle filter algorithm ^[10]Or the like.Under the known prerequisite of Doppler frequency prior information, will be based on the performance of accurate model algorithm significantly better than those non-algorithms based on model.But when the model that adopts is inaccurate, the performance of algorithm will worsen, even disperse ^[11,12]Therefore the selection of dynamic model is very important based on the algorithm of model to this class.

The present invention at first proposes a kind of brand-new description Doppler frequency and the dynamic model of phase place---polynomial forecast model, on this model based,, obtained a kind of state-space method---polynomial prediction-odorlessness kalman filter method of new estimating Doppler frequency in conjunction with the odorlessness Kalman filtering algorithm.Under the framework of state filtering, this algorithm also can utilize the detection of its new breath average and judge model is carried out from confirming.No matter how the relative motion between the transceiver changes, as long as the expression formula of relative motion satisfies the piecewise polynomial form, the algorithm that the present invention proposes can effectively be estimated Doppler frequency and phase place.

Summary of the invention

The objective of the invention is to propose a kind of Doppler frequency and phase estimation algorithm of the GPS of being applied to receiver.

1. the polynomial forecast model of Doppler frequency and phase place

The Doppler frequency deviation of receiver is because the relative motion of receiving-transmitting sides causes ^[1], its expression formula is:

f＝v/λ (1)

Wherein v sends out relative velocity between the machine for the receipts machine, and λ is a carrier wavelength.Usually, the carrier wavelength of system is definite known, and like this according to formula (1), the relative motion between the Receiver And Transmitter has just determined Doppler frequency fully.According to the newtonian motion formula, for the constant motion of acceleration, the expression formula of its speed is

v _t＝v ₀+at (2)

V wherein _tFor receiving-transmitting sides at time t relative velocity constantly, v ₀Be the relative velocity of initial time, a is the acceleration of relative motion.

Convolution (1) and formula (2), can obtain acceleration when constant the expression formula of Doppler frequency and phase place be:

f＝v ₀/λ+at/λ

θ _t＝θ ₀+2πtv ₀/λ+πat ²/λ (3)

F wherein _tAnd θ _tBe respectively t Doppler frequency and phase place constantly, θ ₀It is initial doppler phase.As can be seen from the above equation, under the constant situation of acceleration, Doppler frequency and phase place satisfy polynomial form.After formula (3) is by discretization, can use the polynomial prediction filter ^[13,14]Formula (3) is set up state equation.

The basic conception of polynomial prediction filter is: for the L rank multinomial signal of a discretization

$x\left(n\right)=\underset{l=0}{\overset{L}{\mathrm{\Σ}}}p\left(l\right)n^l---\left(4\right)$

(wherein p (l) (l=0,1 ..., L) be polynomial coefficient, n is the moment of dispersing), can use the x (n) and the value in front (K-1) the individual moment thereof [x (n-K+1) ..., x (n-1)] linear combination come the future value x (n+N) of prediction signal, promptly

$x(n+N)=\underset{k=0}{\overset{K-1}{\mathrm{\Σ}}}h\left(k\right)x(n-k)---\left(5\right)$

Formula (5) be one with h (k) (k=0,1 ..., be the FIR filter of coefficient K-1), be called the polynomial prediction filter.Consider that the signal in the real system all is mixed with noise, if the noise gain of constraint filter

Minimum ^[13]Utilize Lagrangian number multiplication, convolution (4) and formula (5) make the optimal solution of the multinomial filter coefficient of noise gain minimum be ^[13]:

Work as N=1, during L=1,

$h\left(k\right)=\frac{4K-6k-4}{K(K-1)}---\left(6\right)$

Work as N=1, during L=2,

Pairing optimal solution can be referring to document [13] when N and L are worth for other.Be worth noting that to be that the exponent number of polynomial prediction filter need satisfy K 〉=L+1.From formula (6) and formula (7) as can be seen, the coefficient h (k) of polynomial prediction filter formula is only relevant with N, K and L, and with multinomial coefficient p (the l) (l=0 of concrete description signal, 1, ..., L) irrelevant, this means when the multinomial signal being predicted with the polynomial prediction filter, do not need to know prior informations such as the coefficient of multinomial signal own.System state equation based on formula (5) can be expressed as

x _t＝A _PPMx _t-1 (8)

A wherein _PPMFor

H (k) (k=0,1 ..., K-1) be the coefficient of polynomial prediction filter.

In the present invention, Doppler frequency and phase place are estimated simultaneously.The state equation of Doppler frequency can be with reference to the form of formula (8) and is directly obtained; For doppler phase, based on the hypothesis of constant acceleration, doppler phase can be written as

${\mathrm{\&theta;}}_t={\mathrm{\&theta;}}_{t-1}+w_{t-1}{T}_s+0.5w^{\&prime;}{\mathrm{<figure-callout\; id="2"\; label="T\; s"\; filenames="H2010100230630E00012.png,H2010100230630E00041.png"\; state="\{\{state\}\}">T</figure-callout>}}_s^{}$

$={\mathrm{\&theta;}}_{t-1}+w_{t-1}{T}_s+0.5T_s[w_{t-1}-w_{t-2}]---\left(9\right)$

$={\mathrm{\&theta;}}_{t-1}+1.5w_{t-1}{T}_s-0.5w_{t-2}{T}_s$

θ wherein _tBe t doppler phase constantly, w _tBe t doppler angle frequency constantly, T _sBe the sampling period, w ' is the angular frequency rate of change.

To sum up, based on the Doppler frequency of polynomial forecast model and the state equation of phase place be

x _t＝A _PPMx _t-1， (10)

X wherein _t=[θ _tw _tw _T-1] ^T.A _PPMExpression formula be

H (l) (l=0,1) is to work as N=1, L=1, the coefficient of polynomial prediction filter during K=2..

It is irrelevant to describe multi-form concrete coefficient in parameter h (l) in the formula (11) and the formula (3).When this shows with the discrete form of formula (5) described polynomial prediction filter formula or the next equivalently represented formula (3) of the described state-space model formula of formula (8), parameter informations such as the acceleration of the relative motion that need not to know that receiving-transmitting sides is definite, initial velocity and initial doppler phase, and do not have any equivalent error.That is: under the definite known situation of multinomial exponent number, polynomial forecast model formula (8) does not need to introduce the process noise of descriptive model inaccuracy, also do not need to know the multinomial coefficient of definite description multinomial signal, this is different with the middle multinomial model of describing of document [11].

2 odorlessness Kalman filtering algorithms based on polynomial forecast model

In the channel of additive Gaussian noise, suppose that the signal (observational equation) that receives can be expressed as:

Wherein, A _tFor receiving the amplitude of signal, θ _tBe t doppler phase constantly, n _tBe white Gaussian noise, its covariance matrix is R.

Convolution (10) (11) and formula (12) can obtain a kind of new state-space model and come Doppler frequency and phase place are described.Because observational equation is non-linear, any nonlinear filtering algorithm based on state-space model can be used for iterative, from taking all factors into consideration of confirming, the present invention selects for use the odorlessness Kalman filtering algorithm to carry out iterative based on performance, operand and model.

Owing to do not comprise process noise in the polynomial forecast model, therefore when using the odorlessness Kalman filtering algorithm to estimate, the covariance matrix of assignment procedure noise is Q _PPM=0.

In the model of this paper, because the supposition observation noise is an additive white Gaussian noise,, reduce amount of calculation in order to reduce the quantity that Sigma is ordered, the present invention has adopted the odorlessness Kalman filtering algorithm of non-extend type ^[10,15,16], the concrete steps of polynomial prediction-odorlessness Kalman filtering algorithm are as follows:

1. initialization:

${\stackrel{\&OverBar;}{x}}_0=E\left[x_0\right]$ ${\mathrm{<figure-callout\; id="0"\; label="P"\; filenames="H2010100230630E00011.png,H2010100230630E00012.png"\; state="\{\{state\}\}">P</figure-callout>}}_0=E\left[(x_0-{\stackrel{\&OverBar;}{x}}_0){(x_0-{\stackrel{\&OverBar;}{x}}_0)}^T\right]---\left(13\right)$

X wherein ₀Be the initial value of state vector,

Be the expectation of initial state vector, P ₀It is the covariance matrix of initial state vector.

To t ∈ 1,2 .... ∞,

(a) produce the Sigma point: $X_{t-1}=\left[\begin{array}{cc}{\stackrel{\&OverBar;}{X}}_{t-1}& {\stackrel{\&OverBar;}{X}}_{t-1}\&PlusMinus;\sqrt{(\mathrm{nx}+\mathrm{\&lambda;})P_{t-1}}\end{array}\right]$

X _T-1Be the state vector when the t, Be state vector X _T-1Expectation.

Below add horizontal line above similar mark all be expressed as expectation.

(b). the time upgrades:

x _t|t-1＝A _PPMx _t-1 (14)

${\stackrel{\&OverBar;}{X}}_{t|t-1}=\underset{i=0}{\overset{2\mathrm{nx}}{\mathrm{\&Sigma;}}}{W}_i^{\left(m\right)}{X}_{i,t|t-1}---\left(15\right)$

(c). secondary produces the Sigma point,

$y_{t|t-1}=h(x_{t|t-1}^{\&prime;},t)---\left(18\right)$

${\stackrel{\&OverBar;}{y}}_{t|t-1}=\underset{i=0}{\overset{2\mathrm{nx}}{\mathrm{\&Sigma;}}}{w}_i^{\left(m\right)}{y}_{i,t|t-1}---\left(19\right)$

(d) observation is upgraded:

$K_t=P_{x_t{y}_t}{P}_{v_t{v}_t}^{-1}---\left(22\right)$

${\stackrel{\&OverBar;}{x}}_t={\stackrel{\&OverBar;}{x}}_{t|t-1}+K_t(z_t-{\stackrel{\&OverBar;}{y}}_{t|t-1})---\left(23\right)$

$P_t=P_{t|t-1}-K_t{P}_{v_t{v}_t}{K}_t^T---\left(24\right)$

In the superincumbent step, function h ( ^*) be the nonlinear function in the observational equation, z _tBe observation vector, λ is compound scale factor, and nx is the dimension of state vector, w _i ^(m)The weight that corresponding Sigma is ordered during for computation of mean values, w _i ^(c)The weight that corresponding Sigma is ordered during for the calculating covariance.

Certainly the affirmation of 3 polynomial forecast models

In the 1st joint, suppose that the constant acceleration form is satisfied in the relative motion of receiving-transmitting sides.But in practice, become when the acceleration of receiving-transmitting sides relative motion may be or the time constant, this means and cannot directly be applied to polynomial forecast model in the actual signal.According to Weierstrass ^[17]Approximation theorem, any continuous function in the closed interval can approach with arbitrary accuracy with a multinomial.Therefore represent that with multinomial relative motion continuous between the transceiver is rational, particularly when choosing a suitable observation time window, the relative motion of describing between the transceiver with the multinomial or the piecewise polynomial of a low order can obtain satisfactory accuracy, that is to say in actual applications, the time fluctuating acceleration relative motion also can be similar to a multinomial or piecewise polynomial.But, for piecewise polynomial, there is discrete point between each section multinomial, at these discontinuity poinies, acceleration is worth to another from a value mutation, the rate of change that is Doppler frequency is worth to another from a value mutation, and signal enters into another kind of polynomial form from a kind of polynomial form.If use polynomial forecast model at these discontinuity poinies, can not predict accurately signal, promptly predict that with the sampled point before the discontinuity point signal behind the discontinuity point is inaccurate.To become Doppler frequency or phase place in order following the tracks of, must to try every possible means to detect these discontinuity poinies and revised.

For an optimum Kalman filter, enter stable state after, its innovation sequence is the white Gaussian noise sequence of a zero-mean ^[18], but when filter be not when being operated in stable state, this statistical property of innovation sequence will not exist, and that is to say that innovation sequence will no longer be a white Gaussian noise sequence, its average is not 0 yet.According to analysis, state equation (8) can only be described a multinomial signal with certain exponent number, and when describing a piecewise polynomial signal with formula (8), these discontinuity poinies will cause algorithm to transfer to transition state from stable state, destroy the statistical property of innovation sequence.Document has provided the method that different being used to detects new breath statistical property in [18,19].Adopt χ among the present invention ²Whether the average that the detection rule that distributes detects innovation sequence is 0, thereby detects these discontinuity poinies.The new breath of polynomial prediction-odorlessness Kalman filtering can be expressed as $S_t=(z_t-{\stackrel{\&OverBar;}{y}}_{t|t-1}),$ The covariance matrix of new breath is P _Vv. therefore, S _t ^TP _Vv ^-1S _t ^TBe a χ that m the degree of freedom arranged ²Variable, wherein m is the dimension of new breath.Under given level of confidence α, S _t ^TP _Vv ^-1S _t ^TWith χ _α ²The relative size average that embodied innovation sequence whether be zero.When ${S_t}^T{P_{\mathrm{vv}}}^{-1}{S_t}^T>{{\mathrm{\&chi;}}_{\mathrm{\&alpha;}}}^2$ The time, the average of innovation sequence is non-vanishing, means that filtering algorithm does not enter stable state or stable state is broken; When ${S_t}^T{P_{\mathrm{vv}}}^{-1}{S_t}^T\&le;{{\mathrm{\&chi;}}_{\mathrm{\&alpha;}}}^2,$ The time, prove that filtering algorithm has converged to stable state.

Document [20] points out, when all about the hypothesis of model when all correct, the estimation error covariance matrix can converge to a constant matrix gradually, and this moment, filtering algorithm entered stable state.In case stable state is destroyed, the status predication that the user mode equation obtains is no longer accurate, and the information of a large amount of " accurately " is included in the observation sequence, promptly is included in the new breath.Be unlikely to disperse for algorithm can be restrained fast once more, newly ceasing gain matrix need be increased.Under the situation that state equation and observation noise variance are determined, new breath gain increases [20] along with the increase of process noise variance.On the other hand, the variance of process noise has been described the inaccuracy of model, and when model was inaccurate, the variance of process noise also should be increased.Therefore at these discontinuity poinies, when stable state is broken, promptly ${S_t}^T{P_{\mathrm{vv}}}^{-1}{S_t}^T>{{\mathrm{\&chi;}}_{\mathrm{\&alpha;}}}^2,$ State equation is no longer accurate, this seasonal Q _PPM=β I σ _n ², (σ n wherein ²Be the variance of observation noise, β be one greater than 1 positive number), can quicken the convergence of filtering algorithm, make it enter next stable state fast.And after filtering algorithm enters stable state once more, promptly ${S_t}^T{P_{\mathrm{vv}}}^{-1}{S_t}^T\&le;{{\mathrm{\&chi;}}_{\mathrm{\&alpha;}}}^2,$ Again make Q _PPM=0, carry out normal iteration filtering.Therefore complete polynomial prediction-odorlessness Kalman filtering algorithm should comprise from the step of confirming and handles discontinuity point in the permanent accelerated motion of segmentation.After observation is upgraded, calculate $S_t=(z_t-{\stackrel{\&OverBar;}{y}}_{t|t-1})$ And under certain level of confidence α S relatively _t ^TP _Vv ^-1S _t ^TWith χ ²Relation.When ${S_t}^T{P_{\mathrm{vv}}}^{-1}{S_t}^T>{{\mathrm{\&chi;}}_{\mathrm{\&alpha;}}}^2$ The time, make Q _PPM=β I σ _n ²When ${S_t}^T{P_{\mathrm{vv}}}^{-1}{S_t}^T\&le;{{\mathrm{\&chi;}}_{\mathrm{\&alpha;}}}^2$ The time, make QPPM=0.

It should be noted that time interval between adjacent two discontinuity poinies less than convergence of algorithm during the time, the method that the present invention proposes will be in transition state always and can't enter stable state.This moment algorithm to deteriorate to state model inaccurate and comprise the odorlessness Kalman filtering algorithm of big process noise.State estimation result mainly depends on the renewal of new breath to the status predication result, and algorithm can not dispersed, but the optimum also can't reach stable state the time.

4.GPS the processing of modulation signal

In the gps signal of reality, navigation data is modulated on the carrier wave by the BPSK mode, like this generation ± 180 ° the phase overturn of meeting not timing on carrier wave.Therefore need carry out certain processing to data, overcome the influence of phase overturn Frequency Estimation.Ignoring under the prerequisite of noise, the signal of establishing the I that receives, Q two-way satisfies

Wherein, A is the amplitude of carrier wave, supposes that it is a constant, and D (t) is a modulation intelligence, and value is ± 1, and w is a Doppler frequency deviation.The algorithm that proposes for this method in order to overcome ± influence of 180 ° of phase overturns, can be converted to observation

D ²(t)A ²cos[2θ(t)]＝D ²(t)A ²[cos ²(θ(t))-sin ²(θ(t))]＝I ²-Q ² (27)

D ²(t)A ²sin[2θ(t)]＝2D ²(t)A ²sin(θ(t))·cos(θ(t))＝2IQ (28)

D in the following formula ²(t)=1, so just avoided the influence of modulating data to carrier phase.Observational equation can be expressed as:

It should be noted that the result of direct estimation is 2 times of actual frequency deviation, need carry out reprocessing when employing formula (29) during as observational equation.For fear of the influence of signal amplitude, need carry out normalized to I, Q two paths of signals, i.e. I=I/ (I simultaneously to estimated result ²+ Q ²), Q=Q/ (I ²+ Q ²).

Technique effect

The present invention for a kind of based on polynomial forecast model, be applied to the estimating Doppler frequency in the GPS receiver and the algorithm method of phase place.No matter how the relative motion between the transceiver changes, as long as the expression formula of relative motion satisfies the piecewise polynomial form, the algorithm that the present invention proposes can effectively be estimated Doppler frequency and phase place.

Description of drawings

The time dependent track of Fig. 1 emulation Doppler frequency signal.

Fig. 2 algorithm of the present invention and the performance of all contrast algorithms under the simulate signal situation compare, and wherein Fig. 2 a is when 17dB, the estimation for mean frequency error of 100 Monte Carlo simulations; Fig. 2 b is under the different signal to noise ratios, the comparison of various algorithm square root mean square errors.

Fig. 3 is the track of high Dynamic Doppler Frequency signal.Wherein, 3a is a Doppler frequency signal, Fig. 3 b is frequency one Jie's rate of change, and Fig. 3 c is frequency two Jie's rates of change.

Fig. 4 algorithm of the present invention compares with the performance of contrast algorithm under high Dynamic Doppler Frequency RST, when wherein Fig. 4 a is 17dB, and the estimation for mean frequency error of 100 Monte Carlo simulations; Fig. 4 b is under the different signal to noise ratios, the comparison of various algorithm square root mean square errors.

Fig. 5 the inventive method and the Frequency Estimation experimental result of contrast algorithm to actual gps signal, wherein Fig. 5 a is the frequency-tracking result of whole experiment; Fig. 5 b is in experiment starting stage frequency-tracking result.

Embodiment

1. set up the polynomial forecast model of Doppler frequency and phase place according to formula (10) (11);

2. if do not have modulation intelligence on the carrier wave, then adopt formula (12) as observational equation; If have modulation intelligence on the carrier wave, then utilize formula (27), (28) to carry out preliminary treatment to received signal, employing formula (29) is as observational equation.

3. make Q _PPM=0, set certain level of confidence α and constant β;

4. the odorlessness Kalman filtering algorithm shown in the use formula (14-24) carries out iteration filtering, and last in each iteration calculates the variance S of new breath _t ^TP _Vv ^-1S _t ^TIf, greater than threshold value χ _α ², then make Q _PPM=β I σ _n ²If, then make Q less than threshold value _PPM=0.

5. as if no modulation intelligence, phase place that then estimates and frequency are exactly actual doppler phase and frequency, if modulation intelligence is arranged, Doppler frequency that then estimates and phase place are actual phase and frequency 2 times, need carry out reprocessing.

Simulation result:

1. simulated conditions:

In order to assess the performance of method proposed by the invention, method that the present invention is proposed and phase-lock-loop algorithm, expanded Kalman filtration algorithm, odorlessness Kalman filtering algorithm and particle filter algorithm based on traditional white-noise excitation model compare.In the emulation below, in order to represent that conveniently, the amplitude of received signal is a constant in the assumption (12), is without loss of generality, and makes A in the observational equation _t=1.

The algorithm that proposes for the present invention, based on a large amount of numerical simulation results, when doing new breath average and detect, setting level of confidence α is 0.25, constant β=10.

For phase-lock-loop algorithm, in emulation and experiment, all adopted second-order PLL structure commonly used ^[21], the parameter setting of phase-locked loop is middle consistent with document [5], noise bandwidth B _n=50Hz, damping ratio is 0.707.

In the maneuvering target motion model, white-noise excitation model commonly used comprises white noise acceleration (WNA) ^[11]Model and white noise acceleration (WNJ) ^[11]Model.According to formula (3) as can be known, the speed of maneuvering target and Doppler frequency exist the Linear Mapping relation, and therefore the dynamic model with WNA and the pairing Doppler frequency of WNJ is called white noise Doppler frequency rate of change model and white noise Doppler frequency second order rate of change model.

The state equation of white noise Doppler frequency rate of change model can be expressed as:

x(t)＝Ax(t-1)+q(t-1) (30)

State vector x (t)=[θ (t) w (t)] wherein ^T, the covariance matrix of q (t-1) is

T _sBe the sampling time of system, σ _v ²Variance for white noise Doppler frequency rate of change.

The state equation of white noise Doppler frequency second order rate of change can be expressed as

x(t)＝Ax(t-1)+q(t-1) (27)

State vector x (t)=[θ (t) w (t) w ' (t)] wherein ^T, the covariance matrix of q (t-1) is

σ _j ²It is the variance of white noise Doppler frequency second order rate of change.Above-mentioned two kinds of models have been used respectively in below first group and second group of emulation.

At normal temperatures, the power spectral density of thermal noise is-204dB W/Hz that for GPS L1 signal, the bandwidth of signal is 2M Hz ^[2], thermal noise power is about-141dB W like this.For the ordinary GPS receiver, the sensitivity minimization of requirement is-160dB W ^[4], under sensitivity minimization, corresponding signal to noise ratio is approximately-19dB, and carrier-to-noise ratio is 44dB-Hz.Because gps signal adopts band spectrum modulation, the spreading rate of spreading code is 1.023M Hz, and the correlation intergal time of adding up is 1ms, so just exists the storage gain of about 33dB.Therefore in every simulation result below, all provided the Frequency Estimation result of signal to noise ratio from 14dB to 20dB, with the requirement of realistic system, and sampling time Ts all is set at 1ms.

1.1 Doppler frequency simulate signal

In order to verify the performance of algorithm of the present invention when signal satisfies piecewise polynomial form and continuous variation of the satisfied sine of Doppler frequency, provide following one section Doppler frequency simulate signal:

Comprised the Doppler frequency rate of change in this segment signal and be constant, stochastic variable, changed and the situation when undergoing mutation continuously.Formula (30) medium frequency rate of change f '=25Hz/s, v is a stochastic variable that satisfies normal distribution, its variances sigma ²=(25Hz/s) ²The track of corresponding Doppler frequency as shown in Figure 1.

In the emulation of this group, because Doppler frequency rate of change time to time change, therefore the contrast algorithm based on the white-noise excitation model all adopts white noise Doppler frequency rate of change model, the variances sigma of Doppler frequency rate of change _v ²=σ ²=(25Hz/s) ²

1.2 high Dynamic Doppler Frequency signal

In order to verify the performance of algorithm of the present invention under high dynamic model, the high dynamic model that provides in the document [5] has been used in emulation, and initial velocity is-40m/s, in tracing process, one constant acceleration-25g (g is an acceleration of gravity) is arranged, middle two sections each lasting 0.5s in addition, size is 100g/s ²Acceleration.According to formula (3), under the high dynamic motion model of GPS, the corresponding high dynamic model of Doppler frequency as shown in Figure 3.

Because in the emulation of this group, the Doppler frequency rate of change all is one and is not equal to zero constant in the overwhelming majority times, therefore in this emulation, adopt white noise Doppler frequency second order rate of change model, the standard deviation of Doppler frequency second order rate of change is taken as 1/3rd of actual maximum doppler frequency second order rate of change, i.e. 3 σ _j=4900Hz/s ², σ _j≈ 1633Hz/s ²

1.3 actual gps signal

In order to verify the validity of algorithm of the present invention in real system, adopt actual reception to gps data come algorithm is verified.Used actual gps signal is that document [4] is with the data in the subsidiary CD of book in this group emulation.These data receive sampling by the AFE (analog front end) R30 of NordNav company to gps signal, and the sampling position is positioned at gondola Turin.The sample frequency of R30 AFE (analog front end) is 16.3676M Hz, and IF-FRE is 4.1304MHz, and quantified precision is 4bit.The first processing of catching of used data through document [4] software receiver that provides, and through thick frequency offset estimating and compensation, make the absolute value of residual doppler frequency deviation less than 250Hz, and data are carried out preliminary treatment, eliminate the influence of modulating data according to formula (27) (28).To through pretreated data, adopt method among the present invention and other four kinds of contrast algorithms Doppler frequency is estimated and to be followed the tracks of respectively.In experiment.We have intercepted the actual gps data of 2.8s, and coherent accumulation time of integration is 1ms.Algorithm based on the white-noise excitation model all adopts white noise Doppler frequency rate of change model, the variances sigma of Doppler frequency rate of change _v ²=(25Hz/s) ²

2. experimental result:

In first group of emulation, in each stage of Doppler emulating true signal, the algorithm estimated performance that the present invention proposes all is better than other contrast algorithm, and simulation result as shown in Figure 2.When t=600, the Doppler frequency rate of change is undergone mutation, and for the model that the present invention proposes, this moment, discontinuity point appearred in piecewise polynomial.According to the analysis of a last joint, this moment, method of the present invention can detect discontinuity point, and increased the process noise variance in a period of time, restrained again up to algorithm.Therefore when t=600, moving party's algorithm of the present invention has the little peak value of an evaluated error, descends rapidly thereafter, and algorithm is restrained once more, thereby has illustrated in this aspect method from the validity of confirming part.When Doppler frequency satisfied sinusoidal variations, owing to be not that polynomial form is satisfied in strictness, the performance of algorithm of the present invention descended to some extent, but significantly better than other contrast algorithm.

Simulation result under the high Dynamic Doppler Frequency RST as shown in Figure 4.As can be seen, the algorithm that the present invention proposes is 0 o'clock at Doppler frequency second order rate of change from the way, significantly better than other contrast algorithm.And when the absolute value of Doppler frequency second order rate of change is a very big value, because the Doppler frequency rate of change constantly changes, the quick situation of change of Doppler frequency rate of change will can accurately not described with the model that the present invention proposes based on segmentation constant acceleration hypothesis.Therefore at interval t ∈ [2000,2500] and t ∈ [4500,5000], the method that the present invention proposes is in transition state always, the value of process noise is bigger, and algorithm deteriorates to model odorlessness Kalman filtering algorithm when inaccurate, and its estimated performance is poorer than the method that adopts the white-noise excitation model, but when stable state, the algorithm performance that the present invention proposes will contrast algorithm significantly better than other.

The frequency-tracking result of actual signal as shown in Figure 5.As can be seen from the figure, various algorithms can comparatively fast converge to a more stable value for the Frequency Estimation of actual carrier signal, and can not disperse, and the variance that the present invention proposes the estimated result of algorithm obviously is less than other contrast algorithm, thereby verified the validity of this paper algorithm, can in the GPS of reality receiver, be applied.

List of references

[1]Proakis?J?G，Digital?Communications[M]，3rd?ed.New?York：McGrawHill，1995.

[2]E?D?Kaplan?and?C.J.Hegatry.Understanding?GPS?Principles?and?Applications，2 ^nd?ed.Norwood，MA：Artech?House，2006.

[3]James?Bao-Yen?and?TSUI.Fundamentals?of?Global?Positioning?System?Receivers-Asoftware?Approach，2 ^nd?ed.John?Wiley?&?Sons，Inc.Hoboken，New?Jersey，2005.

[4]Kai?Borre，Dennis?M.Akos，Nicolaj?Bertelsen?et?al.A?Software-defined?GPS?andGalileo?Receiver.Boston：Birkhauser，2006.

[5]V?A?Vilnrotter，S?Hinedi，R?Kumar.Frequency?Estimation?Techniques?for?High?DynamicTrajectories.IEEE?Trans.on?Aerospace?and?Electronic?Systems.25(4)：1989.

[6]Steven?M.Kay.Fundamentals?of?Statistical?Signal?Processing?Volume?I：EstimationTheory.New?Jersey：Prentice?Hall，1993.

[7]Simon?J.Julier，Jeffrey?K.Uhlmann.A?New?Extension?of?the?Kalman?Filterto?NonlinearSystems，Proc.Of?AeroSense：The?11 ^th?International?Symposium?on?Aerospace/Defence?Sensing，Simulation?and?Controls，Orlando，Florida，Vol.Multi?Sensor?Fusion，Tracking?and?ResourceManagement?II，1997.

[8]Rudolph?van?der?Merwe，Arnaud?doucet，nando?de?Freitas?and?Eric?Wan.The?unscentedparticle?filter，Aug.2000，Oregon?Graduate?Institute

[9]Julier，Simon?J，Jeffrey?K.Uhlman.Unscented?filtering?and?nonlinear?estimation.Proceedings?of?the?IEEE，vol.92(3)：401-422，2004.

[10]M.Sanjeev?Arulampalam，Simon?Maskell，Neil?Gordon?and?Tim?Clap.A?tutorial?onparticle?filters?for?online?nonlinear?non-Gaussian?Bayesian?tracking.IEEE?Trans.On?SignalProcessing，vol.50(2)：174-188，2002.

[11]X.Rong?Li?and?Vesselin?P.Jilkov，Survey?of?Maneuvering?Target?Tracking.Part?I：Dynamic?Models.IEEE?Trans.on?Aerospace?and?Electronic?Systems，2003，39(4)：1333-1364.

[12]Bar-Shalom，Y.，Li，X.R.，and?Kirubarajan，T.(2001)Estimation?with?Appliations?toTracking?and?Navigation：Theory，Algorithm，and?Software.New?York：Wiley，2001.

[13]P.Heinonen?and?Y.Neuvo，FIR-Median?Hybrid?Filters?with?Predictive?FIRSubstructures[J]，IEEE?Trans.on?Acousitics，Speech?and?Signal?Processing，1998，36(6)：892-899.

[14]Sami?Valiviita，Seppo?J.Ovaska?and?Olli?Vainio.Polynomial?Predictive?Filtering?inControl?Instrumentation：A?Review.IEEE.Trans.On?Industrial?Electronics，vol.46(5)，1999.

[15]Jouni?HartiKainen，Simo?Sarkka.Optimal?filtering?with?Kalman?filters?and?smoothers-a?Manual?for?Matlab?toolbox?EKU/UKF.

[16]Wu，Y.，Hu，D，Wu?M?and?Hu?X.Unscented?Kalmen?Filtering?for?additive?noise?case：augmented?versus?nonaugmented.IEEE?Signal?Processing?Letters?12(5)：357-360，2005.

[17]Eric?W?Weisstein.″Weierstrass?Approximation?Theorem.″From?MathWorld--AWolfram Web Resource.http://mathworld.Wolfram.com/WeierstrassApproximationTheorem.html

[18]R.K.Melrra，On?the?identification?of?variance?and?adaptive?Kalman?filtering.IEEETrans.OnAutomatic?Control，Vol?15(2)：175-184，1970.

[19]E.Parzen.An?approach?to?time?series?analysis.Ann.Math，Statist，Vol.32(5)：951-989，1961.

[20]Fitzgerald?R.J.，Divergence?of?the?Kalman?Filter，IEEE.Trans.on?Automation?Control，1971，16(6)：736-747.

[21]Roland?E.Best.Phase-Locked?Loop?Design，Simulation?and?Application，4 ^th?edition.Columbus，OH：McGraw-Hill，1999.

Claims (4)

1. Doppler frequency and phase estimation method based on a polynomial forecast model, it is characterized in that adopting dynamic model---the polynomial forecast model of describing Doppler frequency, under the assumed condition of constant Doppler frequency rate of change, be expressed as based on the Doppler frequency of polynomial forecast model and the state equation of phase place:

x _t＝A _PPMx _t-1，

X wherein _t=[θ _tw _tw _T-1] ^T, θ _tBe t doppler phase constantly, w _tBe t doppler angle frequency constantly, A _PPMExpression formula be

H (l), l=0,1, be to work as N=1, L=1, K=2, the time polynomial prediction filter coefficient, T _sBe the sampling period, L is the exponent number of the multinomial signal of polynomial prediction filter.

2. Doppler frequency according to claim 1 and phase estimation method is characterized in that in the channel of additive Gaussian noise, suppose that observational equation is expressed as:

$z_t=\left[\begin{array}{c}{A}_t\cos \left({\mathrm{\&theta;}}_t\right)\\ A_t\sin \left({\mathrm{\&theta;}}_t\right)\end{array}\right]+n_t---\left(12\right)$

Wherein, A _tFor receiving the amplitude of signal, θ _tBe t doppler phase constantly, n _tBe white Gaussian noise, its covariance matrix is R; Polynomial forecast model combined with the odorlessness Kalman filtering algorithm carry out iterative, concrete steps are as follows:

1) initialization:

${\stackrel{\&OverBar;}{x}}_0=E\left[x_0\right]$ $P_0=E\left[(x_0-{\stackrel{\&OverBar;}{x}}_0){(x_0-{\stackrel{\&OverBar;}{x}}_0)}^T\right]---\left(13\right)$

X wherein ₀Be the initial value of state vector,

Be the expectation of initial state vector, P ₀It is the covariance matrix of initial state vector;

2) to t ∈ 1,2 .... ∞,

(a) produce the Sigma point:

x _T-1Be the state vector when the t,

Be state vector x _T-1Expectation;

Below add horizontal line above similar mark all be expressed as expectation;

(b). the time upgrades:

x _t|t-1＝A _PPMx _t-1 (14)

(c) secondary produces the Sigma point,

y _t|t-1＝h(x′ _t|t-1，t) (18)

(d) observation is upgraded:

$K_t=P_{x_t{y}_t}{P}_{v_t{v}_t}^{-1}---\left(22\right)$

${\stackrel{\&OverBar;}{x}}_t={\stackrel{\&OverBar;}{x}}_{t|t-1}+K_t(z_t-{\stackrel{\&OverBar;}{y}}_{t|t-1})---\left(23\right)$

$P_t=P_{t|t-1}-K_t{P}_{v_t{v}_t}{K}_t^T---\left(24\right)$

In the superincumbent step, function h (*) is the nonlinear function in the observational equation, z _tBe observation vector, λ is compound scale factor, and nx is the dimension of state vector, w _i ^(m)The weight that corresponding Sigma is ordered during for computation of mean values, w _i ^(c)The weight that corresponding Sigma is ordered during for the calculating covariance.

3. Doppler frequency according to claim 2 and phase estimation method is characterized in that under the framework of state filtering, and be last in each iteration, utilizes the affirmation certainly that model is carried out in the detection and the judgement of new breath average; The process noise covariance matrix of promptly establishing the odorlessness Kalman filter is Q _PPM, new breath is S _t, the covariance matrix of new breath is P _Vv, under given level of confidence α, work as S _t ^TP _Vv ^-1S _t ^T＞χ _α ²The time, Q _PPM=β I σ _n ², work as S _t ^TP _Vv ^-1S _t ^T≤ χ _α ²The time, Q _PPM=0.

4. Doppler frequency according to claim 2 and phase estimation method is characterized in that having cycle information on the carrier wave, then signal is carried out preliminary treatment after, use following observational equation:

$z_t=\left[\begin{array}{c}A\cos \left(2{\mathrm{\&theta;}}_t\right)\\ A\sin \left(2{\mathrm{\&theta;}}_t\right)\end{array}\right]+n_t.$

Images