CN103529461A - Receiver quick positioning method based on strong tracking filtering and Hermite interpolation method - Google Patents

Abstract

The invention discloses a receiver quick positioning method based on strong tracking filtering and a Hermite interpolation method. According to the receiver quick positioning method, based on a known satellite running state, the position of a satellite can be predicted in an interpolation manner, so that the computation complexity is greatly reduced while the positioning precision is ensured; by utilizing a strong tracking filtering device, the quick positioning of a receiver is realized. Specifically the receiver quick positioning method comprises the following steps of firstly extracting ephemeris parameters in a navigation message, and computing the position and the speed of the satellite according to an air interface document; secondly selecting a time interval of a sliding window, and estimating the position of the satellite orbit according to the Hermite interpolation method; finally adjusting a system state vector and a noise covariance matrix. According to the receiver quick positioning method based on the strong tracking filtering and the Hermite interpolation method, which is disclosed by the invention, the positioning resolution of the receiver is carried out by utilizing the strong tracking filtering; meanwhile, the time interval of the sliding window is adjusted, and the prediction and the filtering on the orbit are carried out repeatedly, so that the quick positioning of the receiver is realized; on the premise of ensuring the positioning precision, the computation complexity is greatly reduced, and the requirements of a user on the positioning precision, the real time performance and the robustness of the receiver can be met.

Description

A kind of receiver method for rapidly positioning based on strong tracking filter and Hermite interpolation method

Technical field

The invention belongs to GPS (Global Position System) receiver positioning calculation field, specifically a kind of receiver method for rapidly positioning based on nonlinear filtering and method of interpolation.

Background technology

GLONASS (Global Navigation Satellite System) (GNSS) mainly comprises the Galileo in the GPS of the U.S., Muscovite GLONASS, Europe and the large system of BDS tetra-of China at present.GNSS is extremely important status of performer in military and civilian field, and particularly, in military use, it can improve the targeting rate of guided missile; At civil area, the object that GNSS can be sea, land and sky every aspect positions, navigates, and comprises boats and ships oceangoing voyage and the diversion of approaching, Automobile automatic navigation and Waypoint guiding and the landing etc. of marching into the arena.In addition the most crucial technology that, various surface cars are followed the tracks of and municipal intelligent traffic is managed is exactly GNSS location; It can also be electric power, post and telecommunications and network system time service and the calibrating frequency such as communicate by letter, and comprises and produces lock in time, precision time service and calibrating frequency etc.; The various measurement tasks such as geodetic surveying, earth movement monitoring, engineering survey, engineering project deformation monitoring and resource exploration are inseparable with GNSS.

The ultimate principle that GNSS carries out receiver location is to utilize satellite ephemeris parametric solution to calculate the position of satellite in space coordinates, then according to pseudo range measurement value information, solves the first Nonlinear System of Equations of following n:

$\left\{\begin{array}{c}\sqrt{{(x_S^{\left(2\right)}-x)}^2+{(y_S^{\left(2\right)}-y)}^2+{(z_S^{\left(2\right)}-z)}^2}+{\mathrm{\δt}}_u={\mathrm{\ρ}}_c^{\left(1\right)}\\ \sqrt{{(x_S^{\left(2\right)}-x)}^2+{(y_S^{\left(2\right)}-y)}^2+{(z_S^{\left(2\right)}-z)}^2}+{\mathrm{\δt}}_u={\mathrm{\ρ}}_c^{\left(2\right)}\\ \sqrt{{(x_S^{\left(3\right)}-x)}^2+{(y_S^{\left(3\right)}-y)}^2+{(z_S^{\left(3\right)}-z)}^2}+{\mathrm{\δt}}_u={\mathrm{\ρ}}_c^{\left(3\right)}\\ ...\\ \sqrt{{(x_S^{\left(n\right)}-x)}^2+{(y_S^{\left(n\right)}-y)}^2+{(z_S^{\left(n\right)}-z)}^2}+{\mathrm{\δt}}_u={\mathrm{\ρ}}_c^{\left(n\right)}\end{array}\right.---\left(1\right)$

Wherein, X _recvr=[x, y, z] ^trepresent the position of receiver in WGS-84 coordinate system,

represent the position of n satellite in WGS-84 coordinate system, δ t _urepresent receiver clock correction, for the pseudo-range measurements after proofreading and correct, updating formula is:

ρ is the pseudorange information before proofreading and correct, δ t ^(s)for satellite clock clock correction correction term, computing formula is δ t ^(s)=Δ t ^(s)+ Δ t _r-T _gD, and Δ t ^(s)=a _f0+ a _f1(t-t _oc)+a _f2(t-t _oc) ²,

at ^(s)represent satellite clock correction, Δ t _rrepresent relativistic effect correcting value, T _gDrepresent the time delay of group's ripple.

Satellite transit track is approximately from earth surface 20200Km, accurately calculate the correct location of satellite spatial position and receiver closely related, it is the ephemeris parameter that navigation message is broadcast via satellite that classic method is obtained satellite spatial position coordinates, wherein comprise 6 Keplerian orbit parameters, 9 disturbed modify parameters, through series of complex computing, comprise and solve overdetermined equation, iteration convergence, inverse trigonometric function etc., calculate accurate satellite spatial position coordinates, but need larger computation complexity simultaneously.The term of validity of satellite ephemeris parameter is conventionally at ephemeris reference time t _oewithin former and later two hours, consider that the orbit of satellite is in WGS-84 coordinate system, under non-motor-driven and non-failure conditions, we can find out that satellite spatial orbit has good slickness, therefore the thought of available method of interpolation is carried out satellite transit orbital prediction, than tradition, according to ephemeris parameter, calculate satellite transit orbital position algorithm, method of interpolation can guarantee under the prerequisite of track valuation accuracy, reduce greatly computational complexity, for receiver lays the foundation location fast.Conventional interpolation method has Lagrangian linear interpolation method, Newton divided difference type quadratic interpolattion, Hermite interpolation method, segmentation Hermite interpolation method and splines method of interpolation etc., consider computational complexity and estimation accuracy, the present invention adopts Hermite interpolation method three times, requires interpolating function to equate with original function on node, and the first order derivative on node also equates.

In receiver positioning calculation, according to resolved data origin classification, conventionally there are snapshot algorithm and filtering algorithm.Snapshot algorithm comprises least square method, weighted least-squares method etc.; Filtering algorithm comprises Kalman filtering, strong tracking filter etc.Snapshot algorithm has only been considered each location time information epoch, does not merge the receiver information in the moment in the past, and generally receiver location is gradual change, so often adopt filtering with level and smooth positioning result.Kalman filtering is carried out optimal estimation with mathematical recursion formula to system state, makes the estimated value of system state have least mean-square error (MMSE), and its governing equation and observation equation generally represent in order to lower linear difference equation:

$\left\{\begin{array}{c}x\left(t_k\right)=A(t_k,t_{k-1})\·x\left(t_{k-1}\right)+B\left(t_{k-1}\right)\·u\left(t_{k-1}\right)+w\left(t_{k-1}\right)\\ y_k=C\·x_k+v_k\end{array}\right.---\left(2\right)$

Wherein, t _krepresent that k is measured constantly or k corresponding time of measurement epoch, the input vector of u representative system, the state vector of x representative system, w represents process noise vector, A (t _k, t _k-1) represent from t _k-1to t _kstate-transition matrix constantly, B (t _k-1) represent at t _k-1time relational matrix between etching system input vector and system state vector, the input quantity u of system is an option, GNSS system is generally considered as without any input, C represents the relational matrix between observed quantity and system state, v _knoise vector is measured in representative.Kalman filtering can be divided into prediction and trimming process, by adjusting kalman gain with level and smooth output vector.In GNSS, consider that positioning calculation equation (1) is non-linear overdetermined equations, conventionally utilize Jacobi's transformation, adopt EKF to position and solve.

For Kalman filtering, the reasons such as parameter change due to model simplification, noise statistics out of true, and real system inaccurate to the statistical property modeling of real system original state, make model have a large amount of uncertainties, the nonlinear model set up can not mate completely with the nonlinear system of its description, EKF is poor to the robustness of model uncertainty, therefore introduced strong tracking filter, it is compared with common wave filter, has following good characteristic:

(1) the stronger robustness about model uncertainty;

(2) the extremely strong tracking power about mutation status;

(3) moderate computation complexity.

Strong tracking filter is described nonlinear dynamic system and is conventionally had following structure:

$\left\{\begin{array}{c}{X}_k=f(X_{k-1},u_{k-1})+W_{k-1}\\ Z_k=h\left(X_k\right)+V_k\end{array}\right.---\left(3\right)$

Wherein, X _k∈ R ⁿrepresent system state vector, u _k-1∈ R ^lfor controlled quentity controlled variable, Z _k∈ R ^mrepresent observation vector, W _k-1and V _kfor mutual incoherent zero-mean white Gaussian noise, variance is respectively Q _kand R _k.Strong tracking filter result generally has following structure:

X _k＝X _k，k-1+K _k·γ _k (4)

In formula, X _{k, k-1}=f (X _k-1, u _k-1), γ _k=Z _k-h (X _{k, k-1}), X _{k, k-1}represent priori estimates, K _krepresent filter gain, γ _krepresent predicted value and observed reading residual error.The fading factor that strong tracking filter becomes while passing through adjust in real time the covariance matrix of status predication error and gain, and adopt accordingly fades to the data in past, weakens the impact of stale data on current filter value.

Summary of the invention

The object of the invention is to overcome the weak point in above-mentioned background, a kind of receiver method for rapidly positioning is provided.

The method via satellite ephemeris parameter is calculated satellite spatial position coordinates boundary value, in sliding window time interval, by Hermite interpolation method, carry out satellite orbit estimation, reduce the computational complexity that solves satellite spatial position, and adopt strong tracking filter to carry out receiver to locate fast, avoided the divergence problem causing due to model evaluated error in traditional EKF method, dynamically adjusted afterwards sliding window time interval, with the validity that guarantees that satellite orbit is estimated.

To achieve these goals, the invention provides a kind of receiver method for rapidly positioning based on strong tracking filter and Hermite interpolation method, the method comprises the following steps:

(1) extract satellite navigation message, and demodulate 16 ephemeris parameters according to ICD file (air interface document), XYZ coordinate shaft position and the speed component of initial calculation satellite in WGS-84 coordinate system, calculates satellite position and speed formula as follows:

${\mathrm{<figure-callout\; id="0"\; label="n"\; filenames="HSA0000096181630000021.png,HSA0000096181630000022.png"\; state="\{\{state\}\}">n</figure-callout>}}_0=\sqrt{\frac{\mathrm{GM}}{{\mathrm{<figure-callout\; id="3"\; label="a\; s"\; filenames="HSA0000096181630000021.png,HSA0000096181630000022.png"\; state="\{\{state\}\}">a</figure-callout>}}_s^{}}}$

t _k＝t-t _oe

n＝n ₀+Δn

M _k＝M ₀+nt _k ${\stackrel{\&CenterDot;}{M}}_k=n$

M _k＝E _k-esinE _k ${\stackrel{\&CenterDot;}{E}}_k=\frac{{\stackrel{\&CenterDot;}{M}}_k}{1-e_s\cos E_k}$

$v_k=a\tan \left(\frac{\sqrt{1-{\mathrm{<figure-callout\; id="2"\; label="e"\; filenames="HSA0000096181630000021.png,HSA0000096181630000022.png"\; state="\{\{state\}\}">e</figure-callout>}}^2}\sin E_k}{\cos E_k-e}\right),{\stackrel{\&CenterDot;}{v}}_k=\frac{\sqrt{1-e_s^2}{E}_k}{1-e_s\cos E_k}$

Ф _k＝v _k+w ${\stackrel{\&CenterDot;}{\mathrm{\&Phi;}}}_k={\stackrel{\&CenterDot;}{v}}_k$

δu _k＝Cussin(2Ф _k)+Cuccos(2Ф _k) $\mathrm{\&delta;}{\stackrel{\&CenterDot;}{u}}_k=2{\stackrel{\&CenterDot;}{\mathrm{\&Phi;}}}_k(\mathrm{Cus}\cos \left(2{\mathrm{\&Phi;}}_k\right)-\mathrm{Cuc}\sin \left(2{\mathrm{\&Phi;}}_k\right))$

δr _k＝Crssin(2Ф _k)+Crccos(2Ф _k) $\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_k=2{\stackrel{\&CenterDot;}{\mathrm{\&Phi;}}}_k(\mathrm{Crs}\cos \left(2{\mathrm{\&Phi;}}_k\right)-\mathrm{Crc}\sin \left(2{\mathrm{\&Phi;}}_k\right))$

δi _k＝Cissin(2Ф _k)+Ciccos(2Ф _k) $\mathrm{\&delta;}{\stackrel{\&CenterDot;}{i}}_k=2{\stackrel{\&CenterDot;}{\mathrm{\&Phi;}}}_k(\mathrm{Cis}\cos \left(2{\mathrm{\&Phi;}}_k\right)-\mathrm{Cic}\sin \left(2{\mathrm{\&Phi;}}_k\right))$

u _k＝Φ _k+δu _k ${\stackrel{\&CenterDot;}{u}}_k={\stackrel{\&CenterDot;}{\mathrm{\&Phi;}}}_k+\mathrm{\&delta;}{\stackrel{\&CenterDot;}{u}}_k$

r _k＝a _s(1-ecosE _k)+δr _k ${\stackrel{\&CenterDot;}{r}}_k=a_s{e}_s{\stackrel{\&CenterDot;}{E}}_k\sin E_k+\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_k$

i _k＝i ₀+δi _k+(IDOT)t _k ${\stackrel{\&CenterDot;}{i}}_k=\stackrel{\&CenterDot;}{i}+\mathrm{\&delta;}{\stackrel{\&CenterDot;}{i}}_k$

Ω _k＝Ω ₀+(Ω-Ω _e)t _k-Ω _et _oe ${\stackrel{\&CenterDot;}{\mathrm{\&Omega;}}}_k=\stackrel{\&CenterDot;}{\mathrm{\&Omega;}}-{\stackrel{\&CenterDot;}{\mathrm{\&Omega;}}}_e$

x′ _k＝r _kcosu _k ${\stackrel{\&CenterDot;}{x}}_k^{\&prime;}={\stackrel{\&CenterDot;}{r}}_k\cos u_k-r_k{\stackrel{\&CenterDot;}{u}}_k\sin u_k$

y′ _k＝r _ksinu _k ${\stackrel{\&CenterDot;}{y}}_k^{\&prime;}={\stackrel{\&CenterDot;}{r}}_k\sin u_k-r_k{\stackrel{\&CenterDot;}{u}}_k\cos u_k$

x _k＝x′ _kcosΩ _k-y′ _kcosi _ksinΩ _k ${\stackrel{\&CenterDot;}{x}}_k=-y_k{\stackrel{\&CenterDot;}{\mathrm{\&Omega;}}}_k-({\stackrel{\&CenterDot;}{y}}_k^{\&prime;}\cos i_k-z_k{\stackrel{\&CenterDot;}{i}}_k)\sin {\mathrm{\&Omega;}}_k+{\stackrel{\&CenterDot;}{x}}_k^{\&prime;}\cos {\mathrm{\&Omega;}}_k$

y _k＝x′ _ksinΩ _k+y′ _kcosi _kcosΩ _k ${\stackrel{\&CenterDot;}{y}}_k=x_k{\stackrel{\&CenterDot;}{\mathrm{\&Omega;}}}_k+({\stackrel{\&CenterDot;}{y}}_k^{\&prime;}\cos i_k-z_k{\stackrel{\&CenterDot;}{i}}_k)\cos {\mathrm{\&Omega;}}_k+{\stackrel{\&CenterDot;}{x}}_k^{\&prime;}\sin {\mathrm{\&Omega;}}_k$

z _k＝y′ _ksini _k ${\stackrel{\&CenterDot;}{z}}_k={\stackrel{\&CenterDot;}{y}}_k^{\&prime;}\sin i_k+y_k^{\&prime;}{\stackrel{\&CenterDot;}{i}}_k\cos i_k$

Wherein, t represents GNSS system satellite launch time, t _oerepresent the ephemeris parameter reference time, e _srepresent excentricity, i ₀represent the inclination angle of reference time, Ω ₀represent right ascension of ascending node, w represents the argument of perigee, M ₀represent the average abnormal of reference time, Δ n represents mean motion correction, and IDOT represents change of pitch angle rate, and Ω represents right ascension rate of change, and Cus, Cuc represent latitude correction, and Cis, Cic represent inclination correction, and Crs, Crc represent radius correction.

(2) satellite is carried out respectively under WGS-84 coordinate system to orbit interpolation calculating on XYZ component, take directions X as example, choose t ₁constantly as first interpolation knot, t ₂constantly as second interpolation knot, and there is t ₂-t ₁=Δ t, Δ t value is 20 to 100 seconds, according to ephemeris parameter, calculates this locus f (t of satellite constantly ₁)=x ₁, f (t ₂)=x ₂and satellite velocity f ' (t ₁)=v ₁, f ' (t ₂)=v ₂.

(3) according to three Hermite interpolation methods, to t ₁, t ₂between the X-axis coordinate of satellite carry out Interpolate estimation, interpolation formula is as follows:

H(t)＝h ₁(t)x ₁+h ₂(t)x ₂+H ₁(t)v ₁+H ₂(t)v ₂ (5)

Wherein, t ∈ [t ₁, t ₂], $h_1\left(t\right)=(1+2\frac{t-t_1}{t_2-t_1}){\left(\frac{t-t_2}{t_1-t_2}\right)}^2,h_2\left(t\right)=(1+2\frac{t-t_2}{t_1-t_2}){\left(\frac{t-t_1}{t_2-t_1}\right)}^2,$ $H_1\left(t\right)=(t-t_1){\left(\frac{t-t_2}{t_1-t_2}\right)}^2,H_2\left(t\right)=(t-t_2){\left(\frac{t-t_1}{t_2-t_1}\right)}^2$

(4) by strong tracking filter, carry out receiver positioning calculation.The state vector that makes system is Xstate=[X Vx Y Vy Z Vz b d] ^t, in GNSS system, control vector is 0, state transition function is f=[X+TsVx Vx Y+TsVy Vy Z+TsVz Vz b+Tsd d] ^t, process noise vector W _kcovariance matrix be made as $Q=\left[\begin{array}{cccc}\mathrm{<figure-callout\; id="0"\; label="Qxyz"\; filenames="HSA0000096181630000021.png,HSA0000096181630000022.png"\; state="\{\{state\}\}">Qxyz</figure-callout>}& & 0& \\ & \mathrm{Qxyz}& & \\ & & \mathrm{Qxyz}& \\ 0& & & Q_{\mathrm{bd}}\end{array}\right],$ Wherein

$\mathrm{Qxyz}={\mathrm{Sv}}_x\left[\begin{array}{cc}{T}_s^3/3& T_s^2/2\\ {\mathrm{<figure-callout\; id="2"\; label="T\; s"\; filenames="HSA0000096181630000021.png,HSA0000096181630000022.png"\; state="\{\{state\}\}">T</figure-callout>}}_s^{}/2& \mathrm{Ts}\end{array}\right],Q_{\mathrm{bd}}=\left[\begin{array}{cc}\mathrm{St}\&CenterDot;\mathrm{Ts}+\mathrm{Sf}\&CenterDot;T_s^3/3& \mathrm{Sf}\&CenterDot;T_s^2/2\\ \mathrm{Sf}\&CenterDot;T_s^2/2& \mathrm{Sf}\&CenterDot;\mathrm{Ts}\end{array}\right]$ And St is receiver clock correction noise figure, Sf is that receiver floats noise figure, Sv frequently _xfor the speed noise power spectral density value on X coordinate components, Ts is resolving time interval, observation noise V _kcovariance matrix be made as $R={\mathrm{\&sigma;}}_{\mathrm{<figure-callout\; id="2"\; label="p"\; filenames="HSA0000096181630000021.png,HSA0000096181630000022.png"\; state="\{\{state\}\}">p</figure-callout>}}^2\&CenterDot;\mathrm{eye}\left(4\right),$ σ _pfor receiver is measured noise, ${\mathrm{\&sigma;}}_{\mathrm{<figure-callout\; id="2"\; label="p"\; filenames="HSA0000096181630000021.png,HSA0000096181630000022.png"\; state="\{\{state\}\}">p</figure-callout>}}^2={\mathrm{\&sigma;}}_{\mathrm{<figure-callout\; id="2"\; label="CS"\; filenames="HSA0000096181630000021.png,HSA0000096181630000022.png"\; state="\{\{state\}\}">CS</figure-callout>}}^2+{\mathrm{\&sigma;}}_{\mathrm{<figure-callout\; id="2"\; label="d"\; filenames="HSA0000096181630000021.png,HSA0000096181630000022.png"\; state="\{\{state\}\}">d</figure-callout>}}^2+{\mathrm{\&sigma;}}_{\mathrm{<figure-callout\; id="2"\; label="RNM"\; filenames="HSA0000096181630000021.png,HSA0000096181630000022.png"\; state="\{\{state\}\}">RNM</figure-callout>}}^2+\mathrm{\&sigma;},$ σ _cSthe error to standard deviation of the satellite clock bias model that representative is partly produced by GNSS ground monitoring, σ _drepresentative is atmosphere time delay correction error standard deviation on route of transmission, σ _rNMthe measuring error standard deviation that representative is relevant with receiver and multipath, σ represents the error that other directly do not enumerate out, such as antenna phase noise, the time delay of group's ripple, receiver thermonoise etc.Prior estimate covariance matrix is made as P=eye (8), and observation equation is:

$h=\sqrt{{(\mathrm{Xs}-X)}^2+{(\mathrm{Ys}-Y)}^2+{(\mathrm{Zs}-Z)}^2}+b,$ Receiver original state is made as receiver approximate location and travelling speed.Utilize strong tracking filter to carry out state estimation algorithm as follows:

$X_k=X_{k,k-1}+{\tilde{y}}_k$

(6)

P _k＝[I-K _kH(X _k，k-1)]P _k，k-1

Wherein

X _k，k-1＝f(X _k-1，u _k-1)

P _k，k-1＝λ _k-1f(X _k-1，u _k-1)P _k-1f ^T(X _k-1，u _k-1)+Q _k-1

K _k＝P _k，k-1H ^T(X _k-1)[H(X _k-1)P _k，k-1H ^T(X _k-1)+R _k] ^-1

${\tilde{y}}_k=y_k-h\left(X_{k,k-1}\right)$

$P_{{\tilde{y}}_k}=E\left[{\tilde{y}}_k{\tilde{y}}_k^T\right]\&ap;H\left(X_{k,k-1}\right)P_{k,k-1}{H}^T\left(X_{k,k-1}\right)+R_k$

In formula, λ _kfor self-adaptation fading factor, can be determined by method below:

${\mathrm{\&lambda;}}_k=\left\{\begin{array}{cc}{\mathrm{\&lambda;}}_{0,k},& {\mathrm{\&lambda;}}_{0,k}\&GreaterEqual;1\\ 1,& {\mathrm{\&lambda;}}_{0,k}<1\end{array}\right.$

${\mathrm{\&lambda;}}_{0,k}=\frac{\mathrm{Tr}\left(N_k\right)}{\mathrm{Tr}\left(M_k\right)},\begin{array}{c}{N}_k={\mathrm{<figure-callout\; id="0"\; label="V"\; filenames="HSA0000096181630000021.png,HSA0000096181630000022.png"\; state="\{\{state\}\}">V</figure-callout>}}_{0,k}-H\left(X_{k,k-1}\right)Q_{k-1}{H}^T\left(X_{k,k-1}\right)-l_k{R}_k\\ M_k=H\left(X_{k,k-1}\right)F(X_{k,k-1},u_k)P_{k-1}{F}^T(X_{k,k-1},u_k)H^T\left(X_{k,k-1}\right)\end{array}$

$H\left(X_{k,k-1}\right)=\frac{\&PartialD;h\left(X_k\right)}{\&PartialD;X}{|}_{X_k=X_{k,k-1}},F(X_{k,k-1},u_k)=\frac{\&PartialD;f(X_k,u_k)}{\&PartialD;X}{|}_{X_k=X_{k,k-1}}$

$F(X_{k-1},u_{k-1})=\frac{\&PartialD;f(X_{k-1},u_{k-1})}{\&PartialD;X}{|}_{X_{k-1}=X_{k-1}},{\mathrm{<figure-callout\; id="0"\; label="V"\; filenames="HSA0000096181630000021.png,HSA0000096181630000022.png"\; state="\{\{state\}\}">V</figure-callout>}}_{0,k}=\left\{\begin{array}{cc}{\tilde{y}}_1{{\tilde{y}}_1}^T,& k=0\\ \frac{\mathrm{\&rho;}{\mathrm{<figure-callout\; id="0"\; label="V"\; filenames="HSA0000096181630000021.png,HSA0000096181630000022.png"\; state="\{\{state\}\}">V</figure-callout>}}_{0,k-1}+{\tilde{y}}_k{{\tilde{y}}_{\mathrm{<figure-callout\; id="1"\; label="k\; T"\; filenames="HSA0000096181630000021.png,HSA0000096181630000022.png"\; state="\{\{state\}\}">k</figure-callout>}}}^T}{1+\mathrm{\&rho;}},& k\&GreaterEqual;1\end{array}\right.$

And there is 0.95≤ρ≤0.995 for forgetting factor; l _k>=1 for weakening the factor, can get following equation:

l _k＝1-d _k， $d_k=\frac{1-\mathrm{\&rho;}}{1-{\mathrm{\&rho;}}^{k+1}}$

Parameter in above formula is to utilize orthogonality principle, by certain approximate and abbreviation, obtains.

(5) complete after above-mentioned steps, interval to [t with " sliding window " form traveling time ₂, t ₃], and have t ₃-t ₂=Δ t, returns to step (2), completes at [t ₂, t ₃] satellite orbital position Interpolate estimation and receiver positioning calculation in the period.

The receiver that the present invention has realized based on strong tracking filter and Hermite interpolation method is located fast.This invention has the characteristic of flatness according to Navsat orbit, use Hermite interpolation method to carry out satellite orbit Interpolate estimation, and channel estimation error is 10 in-orbit ^-6under meter Qian Ti, use strong tracking filter to carry out receiver positioning calculation, receiver location is resolved almost without any impact, and the computings such as iteration, inverse trigonometric function and overdetermined equation while as far as possible having avoided calculating satellite orbital position in classic method, greatly reduce computational complexity, meet the requirement of user to receiver location accuracy, real-time and robustness.

Accompanying drawing explanation

Fig. 1 is receiver method for rapidly positioning process flow diagram of the present invention;

Fig. 2 is position and the speed component figure of GNSS satellite transit track under WGS-84 coordinate system, and Data Source is the WAMC research station almanac data that 7:02 issued on March 14th, 2013 in IGS (international GNSS service);

Fig. 3 is satellite transit position and the physical location comparison diagram that Hermite interpolation method is estimated, wherein XYZ site error horizontal ordinate is an interpolation time, ordinate is the absolute value of the difference of Interpolate estimation value and satellite actual value, average error is 4.320618e-06 rice, and variance is 1.527342e-10 rice.

Fig. 4 adopts receiver that strong tracking filter the calculates location components variation diagram under WGS-84 coordinate system, at transverse axis 70 to 80 time point place receiver locations, changing is more greatly because now ionosphere time delay changes, and several meters of deviations appear in corresponding positioning result.

Embodiment

Below in conjunction with accompanying drawing, the present invention is further elaborated.

As shown in Figure 1, a kind of receiver method for rapidly positioning process flow diagram based on strong tracking filter and Hermite interpolation method of realizing to step 5 according to summary of the invention step 1.Receiver signal processing section is synchronized navigation text first, realize bit synchronization and frame synchronization, and extract ephemeris parameter information, signal resolves part and judges ephemeris validity, and initial calculation satellite orbital position and speed component, choose sliding window time interval frontier point, in interval, utilize Hermite interpolation method to carry out satellite orbital position interpolative prediction, co-ordinates of satellite carries out earth rotation rotation correction simultaneously, system state vector and process noise are set, observation noise covariance matrix, utilize the level and smooth positioning calculation result of strong tracking filter, final updating sliding window time interval, re-starting receiver locates fast.

As shown in Figure 2, emulated data source is the almanac data of the WAMC research station in IGS (intemational GNSS service) in 7:02 issue on March 14th, 2013, choose number one satellite as sample, by signal processing, extract navigation message, and change into standard RINEX form, net result is as follows:

2.10 NAVIGATION DATA RINEX VERSION/TYPE

SPIDER V4，3，0，4633 2013 03 14 07：02PGM/RUN BY/DATE

2.1420D-08 7.4506D-09-1.1921D-07 0.0000D+00 IONALPHA

1.2288D+05 0.0000D+00-2.6214D+05 1.9661D+05 ION BETA

6.519258022308D-09 1.598721155460D-14 589824 1731 DELTA-UTC：A0，A1，T，W

16 LEAP SECONDS

END OF HEADER1 13 03 14 05 59 44.0 7.309485226870D-06 3.069544618484D-12 0.000000000000D+006.000000000000D+00-5.234375000000D+01 4.516973906021D-09 1.060559374361D+00-2.568587660789D-06 1.647748751566D-03 1.029483973980D-05 5.153707595825D+033.671840000000D+05-1.862645149231D-09-8.735865195850D-02 1.303851604462D-089.606694685394D-01 1.785937500000D+02 3.614369992561D-01-8.019619635036D-09-5.268076530562D-10 1.000000000000D+00 1.731000000000D+03 0.000000000000D+002.000000000000D+00 0.000000000000D+00 8.381903171539D-09 6.000000000000D+003.599100000000D+05 0.000000000000D+00

The term of validity of ephemeris parameter is conventionally within former and later two hours of ephemeris reference time, when the present invention passes through to adjust GNSS system, carrying out satellite orbital position and speed calculates, result is presented under WGS-84 coordinate system the most at last, as can be seen from the figure the orbit of satellite is smoother, therefore can carry out the prediction of satellite transit orbital position by method of interpolation.

As shown in Figure 3, sliding window time interval is elected 20 seconds as, and the resolving time is spaced apart 1 second, and first interpolation point of take * axle of the present invention is example, and satellite interpolation reference position and speed under WGS-84 coordinate system are respectively:

[-1.046425875853372e+07，1.197094385364485e+07，2.105884511265871e+07]；

[1.669041327874326e+03 ,-2.159592586864055e+03,4.467912702351642e+02]; Interpolation final position and speed are respectively:

[-1.049766826649375e+07，1.192776156128369e+07，2.016768955414204e+07]；

[-1.671907245399891e+03，-2.158632823180821e+03，4.376523731653159e+02]；

Hermite interpolation polynomial is H (t)=h ₁(t) x ₁+ h ₂(t) x ₂+ H ₁(t) v ₁+ H ₂(t) v ₂, h wherein ₁=9.9275e-01, h ₂=7.2500e-03, H ₁=9.0250e-01, H ₂=-4.7500e-02, the satellite orbital position that final Interpolate estimation goes out is-1.046592787167068e+07.In an interpolation section, the difference of the satellite orbital position estimating and real satellite position is 1e-05 rice magnitude, causes micron order positioning result error, therefore in receiver positioning calculation process, the impact of satellite orbital position evaluated error can be ignored.

As shown in Figure 4, receiver initial state vector is [2336245.4900,0 ,-3846430.3478,0,4505277.2221,0 ,-3.150713910084174,0.05], receiver clock correction is that the result calculating by least square method is carried out initial assignment, St=36, Sf=0.01

σ _p=6.0, forgetting factor is set to ρ=0.995, system state transfer function, observation equation, process noise vector W _kwith observation noise vector V _kcovariance matrix be made as respectively f=[X+TsVx Vx Y+TsVy Vy Z+TsVz Vz b+Tsd d] ^t, $h=\sqrt{{(\mathrm{Xs}-X)}^2+{(\mathrm{Ys}-Y)}^2+{(\mathrm{Zs}-Z)}^2}+b,Q=\left[\begin{array}{cccc}\mathrm{<figure-callout\; id="0"\; label="Qxyz"\; filenames="HSA0000096181630000021.png,HSA0000096181630000022.png"\; state="\{\{state\}\}">Qxyz</figure-callout>}& & 0& \\ & \mathrm{Qxyz}& & \\ & & \mathrm{Qxyz}& \\ 0& & & Q_{\mathrm{bd}}\end{array}\right],R={\mathrm{\&sigma;}}_{\mathrm{<figure-callout\; id="2"\; label="p"\; filenames="HSA0000096181630000021.png,HSA0000096181630000022.png"\; state="\{\{state\}\}">p</figure-callout>}}^2\&CenterDot;\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right],$

Utilize strong tracking filter to carry out system state estimation, equation is as follows:

$X_k=X_{k,k-1}+{\tilde{y}}_k$

P _k＝[I-K _kH(X _k，k-1)]P _k，k-1

Wherein,

X _k，k-1＝f(X _k-1，u _k-1)

P _k，k-1＝λ _k-1f(X _k-1，u _k-1)P _k-1f ^T(X _k-1，u _k-1)+Q _k-1

K _k＝P _k，k-1H ^T(X _k-1)[H(X _k-1)P _k，k-1H ^T(X _k-1)+R _k] ^-1

${\tilde{y}}_k=y_k-h\left(X_{k,k-1}\right)$

$P_{{\tilde{y}}_k}=E\left[{\tilde{y}}_k{\tilde{y}}_k^T\right]\&ap;H\left(X_{k,k-1}\right)P_{k,k-1}{H}^T\left(X_{k,k-1}\right)+R_k$

In formula, λ _kfor self-adaptation fading factor, can be determined by method below:

${\mathrm{\&lambda;}}_k=\left\{\begin{array}{cc}{\mathrm{\&lambda;}}_{0,k},& {\mathrm{\&lambda;}}_{0,k}\&GreaterEqual;1\\ 1,& {\mathrm{\&lambda;}}_{0,k}<1\end{array}\right.$

${\mathrm{\&lambda;}}_{0,k}=\frac{\mathrm{Tr}\left(N_k\right)}{\mathrm{Tr}\left(M_k\right)},\begin{array}{c}{N}_k={\mathrm{<figure-callout\; id="0"\; label="V"\; filenames="HSA0000096181630000021.png,HSA0000096181630000022.png"\; state="\{\{state\}\}">V</figure-callout>}}_{0,k}-H\left(X_{k,k-1}\right)Q_{k-1}{H}^T\left(X_{k,k-1}\right)-l_k{R}_k\\ M_k=H\left(X_{k,k-1}\right)F(X_{k,k-1},u_k)P_{k-1}{F}^T(X_{k,k-1},u_k)H^T\left(X_{k,k-1}\right)\end{array}$

$H\left(X_{k,k-1}\right)=\frac{\&PartialD;h\left(X_k\right)}{\&PartialD;X}{|}_{X_k=X_{k,k-1}},F(X_{k,k-1},u_k)=\frac{\&PartialD;f(X_k,u_k)}{\&PartialD;X}{|}_{X_k=X_{k,k-1}}$

$F(X_{k-1},u_{k-1})=\frac{\&PartialD;f(X_{k-1},u_{k-1})}{\&PartialD;X}{|}_{X_{k-1}=X_{k-1}},{\mathrm{<figure-callout\; id="0"\; label="V"\; filenames="HSA0000096181630000021.png,HSA0000096181630000022.png"\; state="\{\{state\}\}">V</figure-callout>}}_{0,k}=\left\{\begin{array}{cc}{\tilde{y}}_1{{\tilde{y}}_1}^T,& k=0\\ \frac{\mathrm{\&rho;}{\mathrm{<figure-callout\; id="0"\; label="V"\; filenames="HSA0000096181630000021.png,HSA0000096181630000022.png"\; state="\{\{state\}\}">V</figure-callout>}}_{0,k-1}+{\tilde{y}}_k{{\tilde{y}}_{\mathrm{<figure-callout\; id="1"\; label="k\; T"\; filenames="HSA0000096181630000021.png,HSA0000096181630000022.png"\; state="\{\{state\}\}">k</figure-callout>}}}^T}{1+\mathrm{\&rho;}},& k\&GreaterEqual;1\end{array}\right.$

And there is 0.95≤ρ≤0.995 for forgetting factor; l _k>=1 for weakening the factor, gets following equation:

l _k＝1-d _k， $d_k=\frac{1-\mathrm{\&rho;}}{1-{\mathrm{\&rho;}}^{k+1}}$

Considering that under ionospheric corrections prerequisite, the final positioning calculation result of receiver as shown in Figure 4.

Claims (6)

1. the receiver method for rapidly positioning based on strong tracking filter and Hermite interpolation method, its feature comprises following five steps:

Steps A: according to ephemeris parameter initial calculation satellite orbital position and speed component;

Step B: extract sliding window time interval boundary satellite transit orbital position and speed component;

Step C: utilize Hermite interpolation method to carry out satellite transit orbital position in sliding window time interval and estimate;

Step D: adopt strong tracking filter to carry out receiver positioning calculation;

Step e: adjust sliding window time interval, repeating step B.

2. receiver method for rapidly positioning according to claim 1, is characterized in that: steps A is extracted 16 ephemeris parameters by navigation message, XYZ coordinate position and the speed component of initial calculation satellite under WGS-84 coordinate system.

3. receiver method for rapidly positioning according to claim 1, it is characterized in that: step B determines starting point and the terminal of interpolation knot by sliding window time interval, and extract corresponding satellite orbital position and speed component, the interpolation section duration is 20～100 seconds.

4. receiver method for rapidly positioning according to claim 1, is characterized in that: it is consistent with protosatellite orbit curve at Nodes that step C meets interpolation curve, and first order derivative also equates.

5. receiver method for rapidly positioning according to claim 1, it is characterized in that: in step D, system state vector comprises position, speed, receiver clock correction and frequently floats, process noise and observation noise are additive white Gaussian noise, and forgetting factor value is 0.995.

6. receiver method for rapidly positioning according to claim 1, is characterized in that: in step e, adjust sliding window time initial point position, the sliding window duration remains unchanged.

Images