CN101839987B - 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 tracing channel as much as possible.Realize the GPS software receiver of low signal-to-noise ratio, the adjustment of high dynamic self-adapting.

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 the length of correlation integral to some extent, 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 after the Doppler frequency skew (below be called carrier frequency) 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 through integration.If the energy of signal obviously surpasses noise energy, explain 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 need search for altogether 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,, C/A sign indicating number start address adopt FFT to be optimized when searching for, 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 the 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 DFTs 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

is the complex conjugate of H (k).Calculation procedure when catching is following:

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 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 through 1 millisecond of correlation time.

Traditional algorithm flow process to a certain satellite caught is 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 be subject to noise effect.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 up of ring retard (DLL) and phaselocked loop (PLL), and peeling off at phaselocked loop of carrier wave accomplished.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 _iBy the carrier frequency of tracking satellite, θ _fBe local carrier frequency, NCO is a digital controlled oscillator.In tracing process, θ _iCan not directly obtain, in fact, θ _iAnd θ _fDifference estimate through 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 following,

$\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="HSA00000134415200024.png,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="HSA00000134415200024.png,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="HSA00000134415200024.png,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; Get

and make system satisfy the critical damping condition, loop noise bandwidth is by computes

$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 be as shown in Figure 4, and 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 like Fig. 5 and shown in 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, like 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 and the step response of ring retard 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 following,

$\mathrm{\&Delta;c}=c_i-c_f=\frac{1}{2}\frac{\sqrt{{\mathrm{<figure-callout\; id="2"\; label="I\; E"\; filenames="HSA00000134415200024.png,HSA00000134415200042.png"\; state="\{\{state\}\}">I</figure-callout>}}_E^{}+Q_E^2}-\sqrt{I_{\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HSA00000134415200024.png,HSA00000134415200042.png"\; state="\{\{state\}\}">L</figure-callout>}}^2+Q_L^2}}{\sqrt{I_{\mathrm{<figure-callout\; id="2"\; label="E"\; filenames="HSA00000134415200024.png,HSA00000134415200042.png"\; state="\{\{state\}\}">E</figure-callout>}}^2+{\mathrm{<figure-callout\; id="2"\; label="Q\; E"\; filenames="HSA00000134415200024.png,HSA00000134415200042.png"\; state="\{\{state\}\}">Q</figure-callout>}}_E^{}}+\sqrt{I_{\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HSA00000134415200024.png,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 generates two paths of signals after the quadrature reproduction carrier wave that digital controlled oscillator produces reappears carrier multiplication with intersection, the sequence dot product through the local C/A that generates and time-code sequence, the sequence that shifts to an earlier date half code element and delay half code element generates I again _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, possibly influence the stability of tracking.And, sign indicating number (D sign indicating number) border of in tracing process, need constantly judging and Update Information.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 with 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, automatic track loop bandwidth according to each passage of signal intensity adjustment; 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 tracing channel can be according to signal intensity independent regulation tracking bandwidth; The 2nd, the advantage of irrelevant integration and short-term correlation integral designs the track loop that a kind of high-sensitivity real-time is exported during in conjunction with length; 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 through 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 squints 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.By finding out among the figure that 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.

Description of drawings

Fig. 1 is the calculation flow chart of traditional acquisition algorithm, and wherein: T1 representes 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 representes the number of search rate point, initial value m=29; J representes 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 representes the y that each correlation integral obtains _i(n) maximal value; The Max that Sum representes to obtain through irrelevant integration with; 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 representes iterations, initial value k=0;

Figure 15 is the flow process that algorithm is confirmed on the numeric data code border, and wherein: j representes 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 representes that receiver is started working, along with the 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.

is 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 representes 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 possibly change state after finishing, and the transfer of state is 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, ready state and tracking mode correspondence tracing process.

1, catches

Based on irrelevant integration acquisition mode, because under this mode, 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 joint 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; Adopted two optimisation strategy, 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 to catch accelerating algorithm 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 squint and variance thereby calculate Doppler frequency in the hope of satellite and user's relative velocity by above-mentioned three conditions.Need to confirm 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 need not search for all 29 components.The main flow process of tradition acquisition algorithm is as shown in Figure 1, when supplementary is arranged to catch flow process shown in figure 13.Can know with Figure 13 through comparison diagram 1,, only need calculate the frequency component that is close with quenching frequency having under the auxiliary condition of information.

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); When the carrier dynamic range is big this hunting zone increase to (14KHz ,+14KHz), even 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,, possibly 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.It is maximum once to obtain in the aforementioned calculation integration absolute value, and 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.The algorithm of refining of capture frequency is 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, and need be on the formal border of specified data sign indicating number before of following the tracks of.Secondly, in order to satisfy high dynamic requirement, result's (I, Q value) of each millisecond correlation intergal is carried out filtering, filtering bandwidth is decided based on signal strength signal intensity.Therefore to solve two problems, the one, specified data sign indicating number border, the 2nd, confirm 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 confirming 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 to begin 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 need not 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, like 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.0～19 punishment section that in these 1000 continuous datas, squints successively asks absolute value to sue for peace again, chooses the corresponding skew of maximum one-component, can specified data sign indicating number border.Through 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 is shown in figure 15.

2.2 confirm signal intensity and noise bandwidth

When specified data sign indicating number border, the intensity of statistical signal in passing is to confirm the bandwidth of signal.About the calculating of signal intensity, many algorithms have been proposed, do not enumerate here.Confirm 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 through to the crystal oscillator error modeling, estimates, thereby reaches the purpose of round-off error.But, research shows that the factor that influences clock correction is a lot, comprises carrier vibration and temperature etc., causes being difficult to set up practical accurate model, so uses 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 correction 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$

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

Receiver crystal oscillation frequency error according to obtaining in a period of time is set up the Real-time Error estimation model; 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 amount of calculation 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 based on 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="HSA00000134415200024.png,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)

To 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 representes that receiver is started working, along with the 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.

is 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 representes 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 the fixed interval L unit interval section that moves forward 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 so crystal oscillation frequency error constantly is taken as up-to-date observation frequency error constantly; 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 consider the situation interval apart from the current time short time, if some observation datas are arranged, because the correlativity of these data and next moment error is more intense in short interval near current time; 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 estimates that the crystal oscillation frequency error algorithm flow chart is 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 converts the observational variance battle array R in the observation equation into.INS also estimates state variance in combined filter, through relatively generating the weight of junction filter feedback subcarrier frequency after conversion and the coordinate axis mapping and by the variance of phase detector statistics.

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

Hi-pass filter to phase detector designs like 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="HSA00000134415200024.png,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 confirm 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="HSA00000134415200024.png,HSA00000134415200042.png"\; state="\{\{state\}\}">w</figure-callout>}}_2}+{\mathrm{<figure-callout\; id="2"\; label="C"\; filenames="HSA00000134415200024.png,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, shown in figure 17.w ₂Be that the junction filter velocity error is mapped to the component on satellite-user's line and converts frequency error into.Get C ₁=C ₂=0.01.

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 maybe be because of not obtaining auxiliary fully losing lock.When combined filter auxiliary lock phase round trip degree 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.So should draw

and

two control auxiliary or corrected weights.

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 characteristic 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 and progressively dwindle the interval quick acquisition algorithm of irrelevant integration of search rate; If the maximum stack of this integral result peak value is then successfully caught this satellite greater than threshold value, get into 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 satellite signal acquisition, 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 gets into 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 through to search rate reducing one by one at interval; If frequency accuracy reaches the requirement of ready link, the link of refining is accomplished, and gets into step (C);

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

(D) follow the tracks of link: after signal gets into 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 accomplished to a frame signal, get into 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 skew 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, according to the stack peak value size ordering of 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 are that position and speed via satellite, position and speed, receiver clock correction and the corresponding variance of 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 search rate is spaced apart 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, it is maximum once to obtain in the aforementioned calculation integration absolute value, 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, get into ready link.

7. 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 (C); The method on specified data sign indicating number border is: be all numeric data codes of T to a period of time length; The numeric data code of individual segmentation summation by per 20 milliseconds of summations obtain approximately

; Ask absolute value to sue for peace again to this

individual value again, obtain an output; 0～19 punishment section that in the continuous data that this time length is T, squints successively asks absolute value to sue for peace again, chooses the corresponding skew of maximum one-component, gets final product specified data sign indicating number border; Through 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 adjustment channel noise bandwidth when carrying out the essence tracking phase

B _N＝P-Cons

B wherein _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, for the intermediate frequency data of 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 confirming to be revised by GPS INS or integrated navigation velocity information according to stabilized property; Promptly utilize Hi-pass filter to confirm the phaselocked loop tracking error, and to INS subsystem or junction filter evaluated error is provided as the Doppler frequency offset error; The estimation variance of final Doppler frequency skew does

$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-user's line and converts frequency error into, 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}$

A 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 is characterized in that, in the step (F), under the condition that has the INS signal, utilizes 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_I\&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, and

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 then next constantly crystal oscillation frequency error 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