CN101839987A - Implementation method for self-adaptive GPS software receiver - Google Patents

Abstract

The invention designs a GPS software receiver, which is suitable for various working environments with low signal to noise ratio, high dynamic, assistance of an inertial navigation system (INS) or independent working and the like, automatically regulates performance as required to adapt to different working environments and achieves the optimal information utilization to the fullest extent as possible. An implementation method for the GPS software receiver comprises the following steps of: independently regulating the tracking bandwidth of each tracking channel according to signal strength; combining the advantages of long-term noncorrelation integral and short-term correlation integral to design a tracking loop with high-sensitivity real-time output; realizing real-time estimation of a receiver crystal frequency error model by using a method for estimating fixed intervals; and finally, based on the deep combination, improving capture velocity and tracking dynamic performance of a weak signal. The implementation method for the GPS software receiver can fully use all detectable satellite signals, can improve the tracking dynamic performance of the weak signal by utilizing auxiliary information of INS information, and also accelerates new star search when INS information assistance exists.

Description

A kind of implementation method of self-adaptive GPS software receiver

Technical field

The invention provides catching and tracking of a kind of self-adaptive GPS software receiver, this method is utilized supplementary as far as possible and is carried out the real-time estimation of crystal oscillator error, realize catching fast and tenacious tracking, and under the condition that signal intensity allows, improve the dynamic property of following the tracks of passage as much as possible.Realize the GPS software receiver that low signal-to-noise ratio, high dynamic self-adapting are adjusted.

Background technology

GPS receiver high sensitivity and high dynamic performance are mutual contradiction.General highly sensitive GPS receiver need carry out long irrelevant or correlation product assigns to catch and tracking weak signal, requires carrier motor-driven unsuitable excessive.High dynamic GPS receiver then limits to some extent to the length of correlation integral, and this just makes weak signal to be hunted down and to follow the tracks of, and under the situation that signal is disturbed easy losing lock.

High dynamic GPS software receiver requires the correlation time of signal trace enough short, with 1 millisecond catch with tracking scheme be example, maximum exportable navigation information frequency is 1kHz.

Signal capture

If the maximum speed of related movement of carrier and satellite is 2500m/s, Doppler frequency offset variation scope is-14K～14KHz, the input intermediate-freuqncy signal is divided into per 1 millisecond of signal frame that sampling number is a unit, each signal frame comprises 1023 pseudorandom code elements (C/A code element), just in time is a C/A sign indicating number cycle (ignoring the cycle shift that the Doppler frequency skew causes).But the C/A sign indicating number start offset of these signal frames is unknown, and the IF-FRE (hereinafter referred to as carrier frequency) after the Doppler frequency skew also is unknown.The C/A sign indicating number start offset and the carrier frequency variation of adjacent two signal frames are very little, the drift speed of C/A sign indicating number start address is 0～7.5Hz, carrier frequency variation is no more than 14Hz between adjacent two milliseconds, therefore two adjacent signal frame C/A code element drifts are no more than 0.0075, and carrier frequency variation is less than 0.088rad/ms.The purpose of signal capture is to find out the rough C/A sign indicating number start address and the carrier frequency of a certain signal frame for the signal trace link.

The intermediate-freuqncy signal of receiving at receiving end can be expressed as:

S _IF＝∑S _i+N (1)

S wherein _iBe the intermediate-freuqncy signal of each visible star, N is a noise.

S _i＝A _iD _iC _i?cos(ω _it) (2)

A _iBe signal amplitude, D _iBeing the navigation data sign indicating number, is 1 or-1, C _iIt is C/A sign indicating number sequence.In signal capture, utilize local C/A sign indicating number, digital carrier and the S that generates earlier _IFMultiply each other, obtain the energy of navigation data by integration.If the energy of signal obviously surpasses noise energy, illustrate that the C/A sign indicating number start address and the carrier frequency of local generation C/A sign indicating number start address and carrier frequency and signal is close, then this acquiring satellite success.If signal energy surpasses peak value, then local C/A sign indicating number start address of possibility or carrier frequency are incorrect, perhaps do not have this satellite-signal.Because the Doppler frequency offset variation scope of signal is-14K～14KHz that step-size in search can be made as 1KHz when catching, and needs the search of 29 secondary frequencies altogether.Suppose that the signal sampling frequency is 5MHz, 1 millisecond signal frame has 5000 sampled points, in order to aim at C/A sign indicating number start address, need carry out 5000 times C/A sign indicating number address searching, the number of times that needs altogether to search for is 29 * 5000 times, each search needs 5000 multiplication of execution and addition, and calculated amount is 29 * 5000 * 5000.In order to accelerate search speed, when searching for, C/A sign indicating number start address adopt fast fourier transform to be optimized, and can make search complexity drop to the computation complexity of 29 * 0 (5000log5000).Here introduce the ABC of cyclic convolution earlier.If the impulse response of continuous system is h (t), input signal is x (t), and output signal is y (t), then

$y\left(t\right)=x\left(t\right)*h\left(t\right)={\∫}_{-\∞}^{\∞}x\left(\mathrm{\τ}\right)h(t-\mathrm{\τ})\mathrm{d\τ}---\left(3\right)$

Behind Fourier transform, Y (f)=X (f) H (f), wherein Y (f), X (f), H (f) they are respectively the Fourier transforms of y (t), x (t), h (t).Similar theorem is also arranged in discrete system

$y\left(n\right)=x\left(n\right)*h\left(n\right)=\underset{m=0}{\overset{N-1}{\mathrm{\Σ}}}x\left(m\right)h(n-m)=L^{-1}\left[X\left(k\right)H\left(k\right)\right]---\left(4\right)$

Wherein X (k) and H (k) are respectively the discrete Fourier transform (DFT) of x (n) and h (n), L ^-1It is the inverse discrete Fourier transformer inverse-discrete operator.Character and formula (4) by Fourier transform are not difficult to derive

$\underset{m=0}{\overset{N-1}{\mathrm{\Σ}}}x\left(m\right)h(n+m)=L^{-1}\left[X\left(k\right)\stackrel{\‾}{H}\left(k\right)\right]---\left(5\right)$

Wherein

It is the complex conjugate of H (k).Calculation procedure when catching is as follows:

The local complex carrier signal sequence that generates 29 frequency components

If centre carrier frequency f _C=1.3MHz, then these frequency components are represented 1300KHz, 1300 ± 1KHz ... 1300 ± 14KHz;

(1) the local complex carrier signal sequence that generates 29 frequency components If centre carrier frequency f _C=1.3MHz, then these frequency components are represented 1300KHz, 1300 ± 1KHz ... 1300 ± 14KHz;

(2) with 29 L _Fi(n) multiply each other with intermediate-freuqncy signal x (n) respectively, obtain 29 burst w _i(n), to w _i(n) carry out discrete Fourier transform (DFT) and obtain W _i(k);

(3) with local C/A sign indicating number Fourier transform and get complex conjugate and obtain sequence

(4) 29 W _i(k) respectively with

Multiply each other and obtain Y _i(k);

(5) to each Y _i(k) carry out inverse fourier transform and obtain sequence of complex numbers y _i(n), find out | y _i(n) | in maximum norm, if maximum norm surpasses threshold value, acquisition success then, the pairing subscript n of maximum norm is the start address of C/A sign indicating number, the pairing frequency f of i _iBe initial tracking frequencies;

(6) if (5) step was not found the mould value above threshold value, a little less than then explanation does not exist this satellite-signal or this satellite-signal excessively, can't catch by 1 millisecond of correlation time.

The traditional algorithm flow process that a certain satellite caught as shown in Figure 1.

Signal trace

Since adopt 1 millisecond of time span tracking scheme, the signal more weak to intensity, the track loop potentially unstable, and easily affected by noise.At first introduce the transport function of track loop.Similar with signal capture, signal trace equally need be to the intermediate-freuqncy signal despreading, promptly carries out carrier wave and the C/A sign indicating number is peeled off.Track loop is made of ring retard (DLL) and phaselocked loop (PLL), and peeling off at phaselocked loop of carrier wave finished.In the signal trace starting stage, owing to the C/A sign indicating number start address and carrier frequency and the out of true that obtain from the signal capture device, need to use the carrier tracking loop of big bandwidth, be called FLL (FLL), the design frame chart of FLL such as Fig. 2.

θ among Fig. 2 _iBe the carrier frequency of tracked satellite, θ _fBe local carrier frequency, NCO is a digital controlled oscillator.In tracing process, θ _iCan not directly obtain, in fact, θ _iAnd θ _fDifference estimate by phase detector.Phase detector utilizes the orthogonal code integrated value I of current time and the integrated value I of cross code integrated value Q and previous moment _-1And Q _-1The phase differential of estimation carrier wave.The phase detector of FLL is an approximately linear, and computing formula is as follows,

$\mathrm{\&Delta;\&theta;}={\mathrm{\&theta;}}_i-{\mathrm{\&theta;}}_f=\frac{I_{-1}Q-{\mathrm{IQ}}_{-1}}{N(I_{-1}^2+Q_{-1}^2+{\mathrm{<figure-callout\; id="2"\; label="I"\; filenames="HSA00000134415200042.png"\; state="\{\{state\}\}">I</figure-callout>}}^2+Q^2)}---\left(6\right)$

Wherein N is the sampling number of every frame.F (s) is the kernel function of FLL, adopts proportional-plus-integral controller

$F\left(s\right)=a{\mathrm{\&omega;}}_0+\frac{{\mathrm{\&omega;}}_0^2}{s}---\left(7\right)$

Closed loop transfer function, then

$H\left(s\right)=\frac{k_0{k}_{f}F\left(s\right)}{s+k_0{k}_{f}F\left(s\right)}=\frac{k_0{k}_{f}a{\mathrm{\&omega;}}_{0}s+k_0{k}_f{\mathrm{\&omega;}}_0^2}{{\mathrm{<figure-callout\; id="2"\; label="s"\; filenames="HSA00000134415200042.png"\; state="\{\{state\}\}">s</figure-callout>}}^2+k_0{k}_{f}a{\mathrm{\&omega;}}_{0}s+k_0{k}_f{\mathrm{\&omega;}}_0^2}---\left(8\right)$

Get k ₀=k _f=1, then

$H\left(s\right)=\frac{a{\mathrm{\&omega;}}_{0}s+{\mathrm{\&omega;}}_0^2}{{\mathrm{<figure-callout\; id="2"\; label="s"\; filenames="HSA00000134415200042.png"\; state="\{\{state\}\}">s</figure-callout>}}^2+a{\mathrm{\&omega;}}_{0}s+{\mathrm{\&omega;}}_0^2}---\left(9\right)$

This is a second-order system, gets

Make system satisfy the critical damping condition, loop noise bandwidth is calculated by following formula

$B_n=\frac{{\mathrm{\&omega;}}_0(1+a^2)}{4a}=\frac{{3\mathrm{\&omega;}}_0}{4\sqrt{2}}---\left(10\right)$

In the FLL tracing process, local carrier frequency is accurate inadequately, needs bigger noise bandwidth in order to avoid losing lock makes B _n=300, then Its frequency characteristic and unit-step response such as Fig. 3.

When FLL is followed the tracks of, because the local carrier phase fluctuation is bigger, be not suitable for exporting the navigation data sign indicating number, when local carrier was followed the tracks of accurately, Phase Tracking switched to phaselocked loop by FLL.

The design frame chart of phaselocked loop and FLL are similarly, and as shown in Figure 4, different is that phaselocked loop is except following the tracks of carrier frequency Δ ω _f, also followed the tracks of carrier phase Δ θ simultaneously.Kernel function still adopts proportional-plus-integral controller, and the peak acceleration of supposing carrier is 400m/s ², along tangential movement, the average height angle of satellite is 60 degree, then average Doppler frequency offset variation rate is 1kHz/s=2 π k rad/s ²≈ 6.28rad/s/ms, acquiescence kernel function noise bandwidth is got B _n=30, equivalent phase noise bandwidth is 26, and equivalent frequency noise bandwidth is 14.Can guarantee enough bandwidth nargin.

The frequency of phaselocked loop and the frequency characteristic of Phase Tracking and unit-step response are respectively as shown in Figure 5 and Figure 6.The phaselocked loop phase detector as shown in the formula

Δθ＝arc?tan(Q/I) (11)

The design loop and the FLL of ring retard are similar, as Fig. 7.C wherein _iAnd c _fThe C/A sign indicating number start address of representing this locality of input signal respectively.

The drift velocity of C/A code element is much smaller than the drift velocity of frequency, at carrier acceleration 400m/s ²Situation under, the drift speed of C/A code element is 1.3Hz/s to the maximum.Under the phaselocked loop tracking mode, get k ₀=32, k _f=1, B _n=1, equivalent noise bandwidth is 21.6667.Under the FLL state, get B _n=20, equivalent noise bandwidth is 433.3333.The frequency domain characteristic of ring retard and step response are as shown in Figure 8.

Similar with phaselocked loop, the C/A sign indicating number start address c of input signal _iCan not directly obtain, need utilize the estimation of code element phase detector.The formula of phase detector calculating symbol error is as follows,

$\mathrm{\&Delta;c}=c_i-c_f=\frac{1}{2}\frac{\sqrt{{\mathrm{<figure-callout\; id="2"\; label="I\; E"\; filenames="HSA00000134415200042.png"\; state="\{\{state\}\}">I</figure-callout>}}_E^{}+Q_E^2}-\sqrt{I_{\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HSA00000134415200042.png"\; state="\{\{state\}\}">L</figure-callout>}}^2+Q_L^2}}{\sqrt{I_{\mathrm{<figure-callout\; id="2"\; label="E"\; filenames="HSA00000134415200042.png"\; state="\{\{state\}\}">E</figure-callout>}}^2+{\mathrm{<figure-callout\; id="2"\; label="Q\; E"\; filenames="HSA00000134415200042.png"\; state="\{\{state\}\}">Q</figure-callout>}}_E^{}}+\sqrt{I_{\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HSA00000134415200042.png"\; state="\{\{state\}\}">L</figure-callout>}}^2+Q_L^2}}---\left(12\right)$

I _E, Q _EAnd I _L, Q _LShift to an earlier date and postpone half C/A code element correlation integral by local code respectively and obtain.Fig. 9 has represented the calculation process of these correlation integrals.Intermediate-freuqncy signal is through quadrature reproduction carrier wave that digital controlled oscillator produces and intersect and generate two paths of signals after the reproduction carrier multiplication, again C/A that generates through this locality and time-code sequence, the sequence of half code element and the sequence dot product that postpones half code element generate I in advance _E, Q _E, I, Q and I _L, Q _LSix integrated values, sign indicating number ring phase detector utilizes these six integrated values to estimate the poor of actual C/A sign indicating number and local C/A sign indicating number start address.The phase place of NCO and frequency are followed the tracks of by phaselocked loop among Fig. 9, that is to say, ring retard and phaselocked loop are not independent.I and Q signal are the input signals of phaselocked loop or FLL phase detector.

High dynamic GPS software receiver requirement signal C/N generally more than 36dB-Hz, if the signal frame of current tracking is in 1 millisecond boundary, may influence the stability of tracking.And, in tracing process, need constantly to judge and upgrade numeric data code (D sign indicating number) border.High sensitive receiver then requires long enough integral time, has therefore limited the maneuvering range of carrier.The principle of muting sensitivity GPS software receiver is similar to high dynamic GPS software receiver, does not carefully state here.No matter dynamically high or high sensitive receiver, its design is underaction all, can not satisfy high dynamic, high sensitivity applications requirement simultaneously.

Summary of the invention

The technical matters that the present invention solves mainly contains three, the one, the automatic track loop bandwidth of adjusting each passage according to signal intensity, in allowed band, improve the dynamic property of each passage as far as possible, make receiver under the strong and weak signals and the situation of depositing, can make full use of all detectable satellite-signals.The 2nd, when existing INS information to assist, can utilize its supplementary to improve the tracking dynamic property of weak signal, quicken new star search simultaneously.The 3rd, carried out the estimation of receiver crystal oscillator error real-time model with the fixed interval estimation technique.

The self-adaptive GPS receiver of the present invention's design has carried out 4 innovations: the one, and each follows the tracks of passage can be according to signal intensity independent regulation tracking bandwidth; The 2nd, the advantage of irrelevant integration and short-term correlation integral during in conjunction with length designs a kind of high sensitivity track loop of output in real time; The 3rd, set up the real-time model of receiver crystal oscillator error with the fixed interval estimation technique, be used for the subcarrier tracking loop, improve the tracking stability of receiver; The 4th, on dark combination foundation, improve acquisition speed and weak signal tracking dynamic performance.

In order to meet design requirement, the problem that scheme need be solved be summarized as following some: in loop design, adopt the method for acquisition and tracking stage by stage, reduce the bit error rate and also improve loop stability.In tracing process, promptly by C/A code phase regulate the frame starting point of input intermediate frequency data stream to C/A sign indicating number boundary alignment integral time, solved the problem of numeric data code phase hit.Method with regretional analysis is estimated the crystal oscillation frequency error model, further obtains next Doppler frequency constantly and is offset the subcarrier tracking loop.Receive the high dynamic requirement of functional adaptation in order to make, capture time can not be long, otherwise catch the information that obtains with out-of-date.

Figure 10 is the overall flow figure of the simple form of adaptable software receiver.The present invention is mainly by catching link, the link of refining, ready link, following the tracks of link, combination link (being junction filter), auxiliary link (aided capture, the auxiliary tracking) and form as can be seen from Figure.

Description of drawings

Fig. 1 is the calculation flow chart of traditional acquisition algorithm, and wherein: T1 represents acquisition threshold;

Fig. 2 is a FLL;

Fig. 3 a is that FLL frequency characteristic, Fig. 3 b are unit-step responses;

Fig. 4 is the phase-locked loop systems block diagram;

Fig. 5 a is that phase-locked loop frequency tracking frequencies characteristic, Fig. 5 b are step response;

Fig. 6 a is that phaselocked loop Phase Tracking frequency characteristic, Fig. 6 b are step response;

Fig. 7 is the ring retard system chart;

Fig. 8 a is that ring retard frequency characteristic, Fig. 8 b are step response;

Fig. 9 is the ring retard schematic diagram;

Figure 10 is software receiver overall flow figure;

Figure 11 is the channel status transition diagram, and the channel status transition diagram has shown the condition of the state and the state exchange of receiver channel;

Figure 12 be no supplementary catch the accelerating algorithm process flow diagram,

Wherein: m represents the number of search rate point, initial value m=29; J represents to carry out the number of times of irrelevant integration, initial value j=1; T ₁Represent irrelevant integrated value threshold value; T ₂Represent irrelevant integral number of times threshold value; Max represents the y that each correlation integral obtains _i(n) maximal value; Sum represent to obtain through irrelevant integration Max's and; Int (*) expression rounding operation symbol:

Figure 13 is the flow process that the supplementary acquisition algorithm is arranged;

Figure 14 is the flow process of algorithm of refining, and wherein: k represents iterations, initial value k=0;

Figure 15 is the flow process that algorithm is determined on the numeric data code border, and wherein: j represents iterations, initial value j=0;

Figure 16 is that the fixed interval method is estimated the crystal oscillation frequency error algorithm flow chart,

Wherein: L and L* represent the length of fixed interval time, and L*＜L.The physical length of current interval when length represents that receiver is started working, along with entering of the new intermediate frequency data of every frame, the value of length is increased to till the L gradually.N _SAnd N _S* the number of representing current slot and legacy data time period visible satellite among the fixed interval L respectively.E _iWith

Represent that respectively the may observe number of satellite is greater than four or less than four unit interval section.

Be illustrated in the crystal oscillation frequency error observed reading that number of satellite obtains greater than 4 unit interval section.S _EWith

Be illustrated respectively in may observe number of satellite in the interval of consideration greater than four with less than the number of four unit interval section.N ₁And N ₂S when carrying out regretional analysis when representing burst length respectively for L _EThe upper and lower bound of value.The quantitative value of the frequency error observation data of actual usefulness when n represents to carry out regretional analysis.When the burst length of considering is L, if S _E＞=N ₁, then make n=N ₁If N ₁＞S _E＞=N ₂, then make n=S _EWhen the burst length of considering is L*, if S _E＞=2, then make n=S _E

Figure 17 Hi-pass filter schematic diagram,

Wherein: w ₁Be phase detector difference value evaluated error (absolute value); T ₁, T ₂Be time constant, be taken as 0.99; Abs is an operator, and expression takes absolute value.

Embodiment

Each passage of GPS software receiver has four kinds of possible states, be respectively out-of-lock condition, the state of refining, ready state, tracking mode, passage will be according to the state of passage to the processing of input intermediate-freuqncy signal, and processing may change state after finishing, and the transfer of state as shown in figure 11.Respectively the processing mode under every kind of state is introduced below.Wherein out-of-lock condition and the corresponding acquisition procedure of the state of refining, the corresponding tracing process of ready state and tracking mode.

1, catches

Based on irrelevant integration acquisition mode, because in this manner, integration lengths can be adjusted as required to the acquisition mode of intermediate-freuqncy signal.As everyone knows, the acquisition procedure calculated amount is excessive, and this section utilizes supplementary and some strategies to reduce its calculated amount.

1.1 there is not the accelerating algorithm of supplementary

For a single passage, traditional correlation capturing process prescription front was told about.Wherein the first two step is can be precalculated for all passages, and the 3rd step only need calculate once for all passages, and two steps of back can not be unifiedly calculated because the C/A sign indicating number of each passage has nothing in common with each other.If the if sampling rate is 5MHz, then one millisecond of needed calculated amount of correlation capturing of all passages is:

29×fs+29×fs×log(fs)+cn×29×fs×2+cn×29×fs×log(fs)

Wherein fs is a frame length 5000, and cn is the largest passages number.

In order further to catch weak signal, the result of the correlation integral that need repeatedly superpose, calculated amount is bigger.In order to reduce calculated amount, two optimisation strategy have been adopted, first is that threshold value is set, after each irrelevant integration, judge whether (5) the maximum norm mean value that obtains of step surpasses threshold value, if surpass, think acquisition success and end to surpass threshold value or acquisition success up to integral number of times otherwise continue irrelevant integration; Second is to remove half search rate behind each irrelevant integration, for example searches for for the first time 29 frequency components, searches for 15 for the second time, and the component that removes is the maximum norm of trying to achieve in (5) step half frequency component less than intermediate value.Search for 8 for the third time, the rest may be inferred, and up to remaining two components, this mode is captured as power with reduction and exchanges capture time for, in fact, is captured as power and reduces also not obviously, lower because the frequency component that eliminates comprises the likelihood ratio of actual frequency.Use after the accelerating algorithm, first millisecond to catch calculated amount the same, from second millisecond, for catching unsuccessful passage, after once carry out irrelevant calculated amount of catching and once reduce half than preceding.No supplementary catch accelerating algorithm as shown in figure 12.

1.2 the accelerating algorithm of supplementary is arranged

Under the condition that obtains the relative velocity supplementary, can remove the process of search rate from, reduce calculated amount in a large number.Supplementary must satisfy following three conditions simultaneously:

(1) satellite of certain passage has obtained complete ephemeris

(2) carrier positions, speed and variance are known

(3) receiver crystal oscillator error and variance thereof are successfully estimated

Can be offset and variance thereby calculate Doppler frequency in the hope of satellite and user's relative velocity by above-mentioned three conditions.Need to determine the frequency range of consideration by the variance of Doppler frequency.With immediate Frequency point in 29 Frequency points that generated be the center, about respectively select n Frequency point to generate the frequency of carrier wave as this locality, and do not need all 29 components are searched for.The main flow process of tradition acquisition algorithm is caught flow process as shown in figure 13 when supplementary is arranged as shown in Figure 1.By comparison diagram 1 and Figure 13 as can be known, having under the auxiliary condition of information, only need calculate the frequency component that is close with quenching frequency.

The innovative point of this method is that general acquisition algorithm is in order to take into full account the Doppler shift that relative motion produced between satellite and the carrier interior, the frequency range of search is (10KHz, + 10KHz), this hunting zone increases to (14KHz when the carrier dynamic range is big, + 14KHz), in addition bigger.When INS, can record the movement velocity of carrier, utilize satellite almanac data to obtain the position and the movement velocity of satellite again, and under the known situation of crystal oscillator error, can calculate the speed of related movement between carrier and the satellite and produce Doppler shift thus, though the Doppler shift that calculates has certain error, can reduce frequency search range greatly, reduce operand and improve real-time.

1.3 refining of capture frequency

For catching of no supplementary,, may be discontented with the requirement of heel track because the capture frequency precision is 500Hz.After acquisition success, increase by one and refine the stage, its thought is: for the frame that newly advances, establishing and catching the frequency that obtains is f, because phase place do not know yet, utilizes f, f-f simultaneously _s, f+f _sThree components, and to 0, π/4, pi/2,3 π/4 four component carry out carrier wave and remove and integration, amount to ten second carrier waves and remove, and after the trigonometric function precomputation, carry out the calculated amount of above ten secondary calculating and little.Obtain in the aforementioned calculation integration absolute value maximum once, its correspondent frequency component is as new frequency, f subsequently _sReduce by half the frame iterative refinement that the next one is newly advanced.f _sInitial value is made as 500Hz.Through three iteration, frequency accuracy can be accurate to 62.5Hz, generally satisfies tracer request.Capture frequency refine algorithm as shown in figure 14.

2. follow the tracks of

In order to improve the accuracy of track loop stability and numeric data code, this paper is the result of 20 frame correlation integrals summation to the output of each numeric data code (D sign indicating number), in order to improve the correctness of numeric data code, need be before formal tracking the border of specified data sign indicating number.Secondly, in order to satisfy high dynamic requirement, result's (I, Q value) of each millisecond correlation integral is carried out filtering, filtering bandwidth is decided according to signal intensity.Therefore to solve two problems, the one, specified data sign indicating number border, the 2nd, determine signal intensity.

2.1 specified data sign indicating number border

The front was said, removal carrier wave and C/A sign indicating number are the starting points of determining every frame intermediate frequency data with the phase place of C/A sign indicating number in the tracing process, promptly for the intermediate frequency data of per 20 frames, are not only the numeric data code boundary alignment, and be from C/A sign indicating number starting point, so just thoroughly solved boundary problem.Facts have proved that with the tracking mode of C/A sign indicating number starting point alignment, calculated amount is slightly less than the mode with receiver time boundary alignment, and does not need to consider boundary problem, precision is higher.The method on specified data sign indicating number border is, to all numeric data codes of a period of time length, as 1000 milliseconds continuous data sign indicating number, obtains the numeric data code of about 50 segmentations summation by per 20 milliseconds of summations, asks absolute value to sue for peace to these 50 values more again, obtains an output.In these 1000 continuous datas, be offset 0～19 punishment section successively and ask absolute value to sue for peace again, choose the skew of maximum one-component correspondence, can specified data sign indicating number border.By increasing the time of summation, can improve the accuracy on numeric data code border, be cost certainly with time.Specified data sign indicating number border algorithm as shown in figure 15.

2.2 determine signal intensity and noise bandwidth

When specified data sign indicating number border, the intensity of statistical signal in passing is to determine the bandwidth of signal.About the calculating of signal intensity, many algorithms have been proposed, do not enumerate here.Determine noise bandwidth by signal intensity, adopt a kind of linear mapping relation:

B _N＝P-Cons

B wherein _NBe noise bandwidth; P is the power of signal, and unit is dB-Hz; Cons is a constant, such as being taken as 20.

2.3 the output of numeric data code

The correlation integral length of tracing process is 1 millisecond, and this is in order to adapt to high dynamic requirement.But each used time span of computational data sign indicating number is 20 milliseconds (satellite times), has so just guaranteed the high sensitivity of receiver.And the adjusting key of dynamic property and sensitivity is the setting of noise bandwidth, because so with regard to receiver generally speaking track loop, can guarantee the highest dynamic property to the different noise bandwidth of signal employing of different capacity, can receive weak signal again.And following assisting of dark combined information, even weak signal also is not easy losing lock.

2.4 crystal oscillation frequency error is estimated

Have frequency error because receiver front end produces the crystal oscillator of the local carrier that is used for down coversion, this will cause centre frequency to be drifted about in time, and the error of crystal oscillator is equivalent to the error of centre frequency.Thereby the auxiliary Doppler shift error of receiver crystal oscillation frequency error influence.Under the situation of compensated receiver crystal oscillator frequency error not, contain error by the speed of GPS observation, the frequency stability of crystal oscillator can reach 10 at present ^-11～10 ^-10, that is to say that the velocity error that causes p.s. can increase the magnitude of 0.01～0.001 meter per second.But the same with INS, the high precision clock price is also very expensive.Thereby more feasible method is by to the crystal oscillator error modeling, estimates, thereby reaches the purpose of round-off error.But, studies show that the factor that influences clock correction is a lot, comprise carrier vibration and temperature etc., cause being difficult to set up practical accurate model, so use the Real-time modeling set method to estimate the receiver crystal oscillation frequency error in this project.

The calculating of may observe number of satellite crystal oscillation frequency error under 4 and above condition.

If air index is ignored to the influence of pseudo range observed quantity, the pseudorange observation equation can be expressed as:

${\mathrm{\&rho;}}_i^j=R_i^j+\mathrm{c\&delta;}{t}_i-c{\mathrm{\&delta;t}}^j$

In the formula:

With

Observation pseudorange and the geometric distance of representing user and receiver respectively, c δ t _iWith c δ t ^jRepresent receiver clock correction and satellite clock correction respectively.

Notice satellite clock correction c δ t ^jCan related parameter revised be arranged according to what provide in the navigation message, can obtain the pseudorange rates observation equation so following formula is carried out differential:

${\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_i^j={\stackrel{\&CenterDot;}{R}}_i^j+c{\stackrel{\&CenterDot;}{\mathrm{\&delta;}}t}_i$

Following formula is carried out linearization to be had:

${\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_i^j=\left(\begin{array}{ccc}{l}_i^j& m_i^j& n_i^j\end{array}\right)\left\{\left[\begin{array}{c}{\stackrel{\&CenterDot;}{x}}^j\\ {\stackrel{\&CenterDot;}{y}}^j\\ {\stackrel{\&CenterDot;}{z}}^j\end{array}\right]-\left[\begin{array}{c}{\stackrel{\&CenterDot;}{x}}_i\\ {\stackrel{\&CenterDot;}{y}}_i\\ {\stackrel{\&CenterDot;}{z}}_i\end{array}\right]\right\}+c{\stackrel{\&CenterDot;}{\mathrm{\&delta;}}t}_i$

In the formula

Be radial vector the direction cosine in the earth's core body-fixed coordinate system of receiver to observation satellite.Because pseudorange rates and satellite velocities information obtain, if observe 4 or above satellite simultaneously, just can obtain the error of receiver crystal oscillator frequency.

Set up the Real-time Error estimation model according to the receiver crystal oscillation frequency error that obtains in a period of time, utilize model further estimate next constantly crystal oscillation frequency error and obtain corresponding carrier doppler frequency shift (FS), be used for the subcarrier tracking loop, realize the stability of carrier tracking loop.

The process of setting up of following mask body statement error model.

1) the low order model of error estimate is a unique selection in estimating in real time

The model that this method is set up is the low order multinomial model.Though higher order polynomial can be similar to arbitrary function, because calculated amount is big, can not satisfy the needs of real-time navigation, so just got rid of high-order model, the exponent number of lower-order model can be chosen according to the performance of crystal oscillator.

The crystal oscillation frequency error model that this algorithm uses is:

${\stackrel{\&CenterDot;}{\mathrm{\&delta;}}t}_i=x_0+x_{1}t+x_2{\mathrm{<figure-callout\; id="2"\; label="t"\; filenames="HSA00000134415200042.png"\; state="\{\{state\}\}">t</figure-callout>}}^2+...+x_n{t}^n$

T express time wherein, x ₀, x ₁, x ₂X _nBe model parameter, n is a little positive integer (such as getting 1,2 or 3)

Two problems mainly considering when 2) using lower-order model

(1) length of time period can not be long, just will solve the problem of time in the past availability of data.

Could adopt low order multinomial model (such as inearized model) to the estimation of model short time period, therefore when setting up model, the low order multinomial model has strict demand to the length of time period, if the time period is long, the low order multinomial model is probably no longer suitable so.

(2) the data point number can not be very few.

When burst length is not very hour, observation data point number can not be very few, otherwise it is excessive to obtain result's error with regretional analysis; If burst length reduces, regression analysis is to corresponding the loosening of requirement meeting of observation data number.

3) receiver crystal oscillator frequency error algorithm for estimating (fixed interval algorithm for estimating)

At two problems using lower-order model mainly to consider, the fixed interval algorithm for estimating is proposed.

(1) algorithm symbol description: L and L* represent the length of fixed interval time, and L*＜L.The physical length of current interval when length represents that receiver is started working, along with entering of the new intermediate frequency data of every frame, the value of length is increased to till the L gradually.N _SAnd N _S* the number of representing current slot and legacy data time period visible satellite among the fixed interval L respectively.E _iWith Represent that respectively the may observe number of satellite is greater than four or less than four unit interval section.

Be illustrated in the crystal oscillation frequency error observed reading that number of satellite obtains greater than 4 unit interval section.S _EWith

Be illustrated respectively in may observe number of satellite in the interval of consideration greater than four with less than the number of four unit interval section.N ₁And N ₂S when carrying out regretional analysis when representing burst length respectively for L _EThe upper and lower bound of value.The quantitative value of the frequency error observation data of actual usefulness when n represents to carry out regretional analysis.When the burst length of considering is L, if S _E＞=N ₁, then make n=N ₁If N ₁＞S _E＞=N ₂, then make n=S _EWhen the burst length of considering is L*, if S _E＞=2, then make n=S _E

(2) at first to adopt a time span be the fixed interval of L to this algorithm, and record is the number situation of each unit interval section visible satellite during this period of time.When new intermediate frequency data frame having been carried out after acquisition and tracking handles, judge the number of current unit interval section visible satellite, if the visible satellite number more than or equal to 4, S then _EAdd 1, otherwise S _EConstant.If the oldest unit interval section visible satellite number is greater than 4 in the interval, then S _ESubtract 1, otherwise S _EConstant.Then the current time state is added interval head end, will carve state the old times and remove from interval end.In brief, make fixed interval L move forward a unit interval section exactly.When fixed interval length is L, no matter E wherein _iWith

How to distribute, as long as S _EMore than or equal to threshold value N ₂, then utilize a up-to-date n observed reading

Set up the crystal oscillation frequency error real-time model, and estimate next crystal oscillation frequency error constantly, be used for the subcarrier tracking loop.

If S _ELess than threshold value N ₂, then fixed interval length is reduced to L* (can get L*=10ms).If S _E* more than or equal to 2, then utilize nearest S _EIndividual

Carry out linear regression and set up error model, and estimate next crystal oscillation frequency error constantly; If S _EEqual 1, next that estimate crystal oscillation frequency error constantly is taken as up-to-date observation frequency error constantly so; If S _E* equal 0, error model is not temporarily estimated.The thought here is, under the situation of long-time interval observation data quantity not sufficient, then consideration is apart from the situation in current time short time interval, if some observation datas are arranged in short interval near current time, because the correlativity of these data and next moment error is more intense, utilize these to carry out regretional analysis than new data so, the model accuracy that obtains also is an acceptable.Estimate next crystal oscillation frequency error constantly with the model that obtains then, be used for the subcarrier tracking loop.

The fixed interval method is estimated the crystal oscillation frequency error algorithm flow chart as shown in figure 16.

2.5 subcarrier ring

Intermediate-freuqncy signal generates tracking frequencies, phase place and pseudo-random code skew through generating I and Q integrated value after the decorrelation through code phase discriminator and frequency phase detector.Simultaneously, the difference value process Hi-pass filter by phase detector output counts frequency variance and pseudorange variance, and is converted to the observational variance battle array R in the observation equation.INS also estimates state variance in combined filter, through relatively generate the weight of junction filter feedback subcarrier frequency after conversion and the coordinate axis mapping with the variance by the phase detector statistics.

Ring retard and the design of phaselocked loop phase detector Hi-pass filter

Hi-pass filter to phase detector designs as Figure 17.

Junction filter is auxiliary to phase-locked loop frequency

$f_C=\frac{f_G\&CenterDot;{\mathrm{<figure-callout\; id="2"\; label="q"\; filenames="HSA00000134415200042.png"\; state="\{\{state\}\}">q</figure-callout>}}_2+f_I\&CenterDot;q_1}{{\mathrm{<figure-callout\; id="1"\; label="q"\; filenames="HSA00000134415200021.png,HSA00000134415200022.png"\; state="\{\{state\}\}">q</figure-callout>}}_1+q_2}---\left(13\right)$

F wherein _CBe the frequency of combination calculation, f _GBe the output frequency of phase detector, f _IIt is the frequency of deriving out by junction filter speed.q ₁And q ₂Ratio reflected f _IAuxiliary f _GDegree.Q is set ₂≡ 1, q ₁Change according to Doppler frequency offset variation rate, promptly determine according to acceleration and satellite sight line,

${\mathrm{<figure-callout\; id="1"\; label="q"\; filenames="HSA00000134415200021.png,HSA00000134415200022.png"\; state="\{\{state\}\}">q</figure-callout>}}_1={\mathrm{<figure-callout\; id="1"\; label="C"\; filenames="HSA00000134415200021.png,HSA00000134415200022.png"\; state="\{\{state\}\}">C</figure-callout>}}_1\left|\stackrel{\&CenterDot;}{f}\right|---\left(14\right)$

The estimation variance of Doppler frequency skew

$q_F=\frac{w_1{w}_2}{{\mathrm{<figure-callout\; id="1"\; label="w"\; filenames="HSA00000134415200021.png,HSA00000134415200022.png"\; state="\{\{state\}\}">w</figure-callout>}}_1+{\mathrm{<figure-callout\; id="2"\; label="w"\; filenames="HSA00000134415200042.png"\; state="\{\{state\}\}">w</figure-callout>}}_2}+{\mathrm{<figure-callout\; id="2"\; label="C"\; filenames="HSA00000134415200042.png"\; state="\{\{state\}\}">C</figure-callout>}}_2\left|\stackrel{\&CenterDot;}{f}\right|---\left(15\right)$

W wherein ₁Be the variance of the Doppler frequency skew of phaselocked loop estimation, as shown in figure 17.w ₂Be that the junction filter velocity error is mapped to the component on satellite-user's line and is converted to frequency error.Get C ₁=C ₂=0.01.

It is the rate of change of Doppler frequency skew.

$\stackrel{\&CenterDot;}{f}\&ap;(a+{\mathrm{\&omega;}}^2{r}_s)\&CenterDot;(r_s-r_u)\&times;\frac{f_{L1}}{c}---\left(16\right)$

α is a specific force, and ω is the angular velocity that satellite changes around the earth, r _sAnd r _uIt is respectively satellite and user's position.f _L1Be the L1 carrier frequency, c is the light velocity.

If the auxiliary degree of junction filter is bigger, phaselocked loop is just more little to the correcting action of junction filter, and this has formed serious positive feedback effect.If the auxiliary degree of junction filter is too small, then phaselocked loop may be because of not obtaining auxiliary fully losing lock.When the auxiliary phase-locked round trip degree of combined filter was big, for fear of positive feedback effect, the degree of phaselocked loop correction combined filter should reduce; When the degree of the auxiliary phaselocked loop of combined filter hour, should increase the correcting action of phaselocked loop, otherwise cause information waste junction filter.Therefore should by

With Two weights that control is assisted or revised.

Claims (12)

1. the method for work of a self-adaptive GPS software receiver, its can be used for high dynamically, low signal-to-noise ratio and have that inertial navigation system (INS) is auxiliary, no INS assists the navigation under four kinds of conditions, its feature comprises:

(A) catch link: according to whether auxiliary signal is arranged, the acquisition loop selected parts are selected and are adopted one of following two kinds of acquisition algorithms: the acquisition algorithm during (1) no auxiliary signal: in the search procedure of carrying out certain satellite signal acquisition, under the condition of no auxiliary signal, adopt the quick acquisition algorithm of irrelevant integration that progressively dwindles the search rate interval; If this integral result maximal value greater than threshold value, is then successfully caught this satellite, enter step (B); If this integral result maximal value is less than threshold value, repeating step (A) then; If irrelevant integral number of times reaches threshold value, think that then this satellite-signal does not exist or signal is crossed weak can not catching; Acquisition algorithm when (2) auxiliary signal being arranged: in the search procedure of carrying out certain acquiring satellite, exist under the condition of auxiliary signal,, adopt the quick acquisition algorithm of frequency search among a small circle by the carrier doppler frequency and the variance thereof that obtain; If acquisition success enters step (C); If catch unsuccessfully, think then that this satellite-signal does not exist or signal is crossed weak can not catching;

(B) link of refining: the link of refining improves the precision of search rate by to search rate reducing one by one at interval; If frequency accuracy reaches the requirement of ready link, the link of refining is finished, and enters step (C);

(C) ready link: the numeric data code to certain-length is handled, and the reference position of specified data sign indicating number enters step (D);

(D) follow the tracks of link: after signal enters the tracking link, according to the intensity adaptive design tracking bandwidth of signal; When the tracking link was carried out decorrelation to each signal frame, making each signal frame to be demodulated was starting point with C/A sign indicating number zero phase all; After tracking is finished to a frame signal, enter step (E);

(E) combination link: step (D) is handled information and the INS metrical information that obtains, and is admitted to junction filter together, carries out information fusion; The combination link obtains the carrier navigation information and with navigation information input step (F);

(F) auxiliary link: auxiliary link utilizes satellite position and velocity information and carrier positions and velocity information to estimate to obtain carrier doppler frequency shift (FS) and variance thereof, thereby can estimate the frequency range of carrier wave; Be at current satellite-signal under the situation of losing lock, this information is used for step (A), be used for dwindling catch link search frequency range with aided capture, thereby improved receiver to height dynamically, the acquisition performance of low signal-to-noise ratio signal; Be at current satellite-signal under the situation of tracking, carrier doppler frequency shift (FS) and variance information thereof that this information and step (D) obtain merge, further obtain accurate more carrier doppler frequency shift (FS) and variance information thereof and be used to follow the tracks of link, receiver is dynamic to height, the tracking performance of low signal-to-noise ratio signal thereby improved;

(G) crystal oscillator estimation of error link: the real-time model that carries out the receiver crystal oscillation frequency error is estimated, obtains corresponding Doppler frequency and information is sent into step (F).

2. the method for work of a kind of self-adaptive GPS software receiver as claimed in claim 1, it is characterized in that, in step (A) method (1), if the maximum stack of this irrelevant integral result peak value is less than threshold value, stack peak value size ordering according to frequency to be searched, eliminate half less search rate of peak value, dwindle frequency component set to be searched next time with maximum acquisition probability.

3. the method for work of a kind of self-adaptive GPS software receiver as claimed in claim 1 is characterized in that, in step (A) method (1), can regulate the irrelevant integral time of catching link automatically according to the intensity of signal.

4. the method for work of a kind of self-adaptive GPS software receiver as claimed in claim 1, it is characterized in that, in step (A) method (2), Doppler frequency and variance thereof be via satellite the position and position and speed, receiver clock correction and the corresponding variance of speed, carrier obtain.

5. the method for work of a kind of self-adaptive GPS software receiver as claimed in claim 1, it is characterized in that, in step (A) method (2), owing to obtained carrier doppler frequency and variance thereof, therefore obtained carrier frequency and variance, caught the link frequency search range thereby dwindled greatly.

6. the method for work of a kind of self-adaptive GPS software receiver as claimed in claim 1 is characterized in that, after acquisition success, increase by one and refine the stage in the step (B): establishing and catching the frequency that obtains is f, supposes that frequency interval is f _s,, utilize f, f-f simultaneously for the frame that newly advances _s, f+f _sThree components, and to initial phase be 0, π/4, pi/2,3 π/4 four component carry out carrier wave and remove and integration, amount to ten second carrier waves and remove, and obtain in the aforementioned calculation integration absolute value maximum once, its correspondent frequency component is as new frequency, f subsequently _sReduce by half the frame iterative refinement that the next one is newly advanced; After the link of refining, enter ready link.

7. the method for work of a kind of self-adaptive GPS software receiver as claimed in claim 1 is characterized in that, in the step (C), the method on specified data sign indicating number border is: to a period of time length is all numeric data codes of T, obtains approximately by per 20 milliseconds of summations

The numeric data code of individual segmentation summation is again to this

Individual value asks absolute value to sue for peace again, obtains an output.Successively at this moment between length be to be offset 0～19 punishment section in the continuous data of T to ask absolute value to sue for peace again, choose the skew of maximum one-component correspondence, get final product specified data sign indicating number border; By increasing the time of summation, can improve the accuracy on numeric data code border, be cost certainly with time.

8. the method for work of a kind of self-adaptive GPS software receiver as claimed in claim 1 is characterized in that, in the step (D), signal utilizes following formula to adjust the channel noise bandwidth when carrying out the essence tracking phase

B _N＝P-Cons

B wherein in the formula _NBe noise bandwidth, P is the power of signal, and Cons is a constant.

9. the method for work of a kind of self-adaptive GPS software receiver as claimed in claim 1, it is characterized in that, step (C) and (D) in, when solving the numeric data code boundary problem, intermediate frequency data for per 20 frames, be not only the numeric data code boundary alignment, and every frame signal is a starting point with C/A sign indicating number zero phase all.

10. the method for work of a kind of self-adaptive GPS software receiver as claimed in claim 1 is characterized in that, in the step (E): the degree of determining to be revised by GPS INS or integrated navigation velocity information according to stabilized; Promptly utilize Hi-pass filter to determine the phaselocked loop tracking error, and provide evaluated error to INS subsystem or junction filter as the Doppler frequency offset error; The estimation variance of final Doppler frequency skew is

$q_F=\frac{w_1{w}_2}{w_1+w_2}+C\left|\stackrel{\&CenterDot;}{f}\right|$

W wherein ₁Be the Doppler frequency offset error that phaselocked loop is estimated, w ₂Be that the junction filter velocity error is mapped to the component on satellite one user's line and is converted to frequency error, C is normal value,

It is the rate of change of Doppler frequency skew;

$\stackrel{\&CenterDot;}{f}\&ap;(a+{\mathrm{\&omega;}}^2{r}_s)\&CenterDot;(r_s-r_u)\&times;\frac{f_{L1}}{c}$

α is a specific force, and ω is the angular velocity that satellite changes around the earth, r _sAnd r _uBe respectively satellite and user's position, f _L1Be the L1 carrier frequency, c is the light velocity.

11. the method for work of a kind of self-adaptive GPS software receiver as claimed in claim 1, it is characterized in that, in the step (F), under the condition that has the INS signal, utilize the dynamic property of the auxiliary Doppler frequency skew raising track loop of junction filter prediction; Junction filter to the phase-locked loop frequency side formula is:

$f_C=\frac{f_G\&CenterDot;q_2+f_1\&CenterDot;q_1}{q_1+q_2}$

F wherein _CBe the frequency of combination calculation, f _GBe the output frequency of phase detector, f _IBe the frequency of deriving out by junction filter speed, q ₁And q ₂Ratio reflected f _IAuxiliary f _GDegree, q wherein ₂Be set to normal value, q ₁Change according to Doppler frequency offset variation rate; q ₁Computing formula is:

$q_1=C\left|\stackrel{\&CenterDot;}{f}\right|$

C is a constant in the formula,

It is Doppler frequency offset variation rate.

12. the method for work of a kind of self-adaptive GPS software receiver as claimed in claim 1, it is characterized in that, utilize the fixed interval method to estimate the crystal oscillation frequency error model in real time in the step (G), further estimate the crystal oscillation frequency error in next moment then and obtain corresponding Doppler frequency skew, in order to the subcarrier tracking loop, thereby guaranteed to follow the tracks of the stability of link.

Images