Transcription
Extended Kalman Filter Methods for TrackingWeak GPS SignalsMark L. Psiaki and Hee Jung, Cornell University, Ithaca, N.Y.BIOGRAPHIES(PRN) spreading code of a received signal so that thereceiver can determine navigation observables anddecode the navigation data message. One necessarytracking algorithm is a delay-locked loop (DLL), whichmaintains phase alignment between the received PRNcode and a replica in the receiver. A receiver uses eithera frequency-locked loop (FLL) or a phase-locked loop toalign a replica of the carrier signal with the receivedcarrier. Various designs exist for these tracking loops,e.g., see Refs. 1 and 2.Mark L. Psiaki is an Associate Professor of Mechanicaland Aerospace Engineering at Cornell University. Hereceived a B.A. in Physics and M.A. and Ph.D. degrees inMechanical and Aerospace Engineering from PrincetonUniversity. His research interests are in the areas ofestimation and filtering, spacecraft attitude and orbitdetermination, and GPS technology and applications.Hee Jung is a Ph.D. candidate in Mechanical andAerospace Engineering at Cornell University. Shereceived her BS and MS in Astronomy from the SeoulNational University in Korea and another MS inAerospace Engineering from Texas A&M University.Her main research interests are orbit and attitudedetermination of satellites with GPS applications.The goal of the present work is to design a combinedDLL/PLL code and carrier tracking loop that is effectiveat tracking GPS L1 Coarse/Acquisition signals (C/A) forvery low carrier-to-noise ratios.This is to beaccomplished without prior knowledge of the navigationdata bits. These tracking loops will be developed usingoptimal estimation techniques. Standard DLLs and PLLs,while robust, are not optimal. PLLs lose lock at lowcarrier-to-noise ratios because of the nonlinearities intheir discriminators and because of dynamic variations ofthe signal's phase 2. Uncertainty about the navigation databits greatly exacerbates the problems of PLLs at lowcarrier-to-noise ratios because bit error rates becomelarge and lead to a destabilizing feedback effect. Thegoal of using optimal techniques is to deal with theseproblems in the best possible way and thereby lower thecarrier-to-noise threshold at which loss of lock occursfor a given level of signal phase dynamics.ABSTRACTA combined phase-locked loop/delay-locked loop hasbeen developed for tracking weak GPS C/A signals. Thiswork enables the use of the weak side-lobe signals thatare available at geosynchronous altitudes. The trackingalgorithm is an extended Kalman filter (EKF) thatestimates code phase, carrier phase, Doppler shift, rateof change of Doppler shift, carrier amplitude and data bitsign. It forms a likelihood function that depends on theerrors between accumulations and their predicted values.It recursively minimizes this likelihood function in orderto track the signal. It deals with data bit uncertainty usinga Bayesian analysis that determines a posterioriprobabilities for each bit sign. A second filter is used toinitialize the EKF. This batch filter starts with coarsecarrier frequency and code phase estimates and refinesthem using maximum likelihood techniques whileestimating the carrier phase and the first PRN codeperiod of a navigation data bit. The resulting system canacquire and maintain lock on a signal as weak as 15 dBHz if the receiver clock is an ovenized crystal oscillatorand if the line-of-sight acceleration variations are as mildas those seen by a geostationary user vehicle.The motivation for this work comes from a desire to useGPS for high-altitude spacecraft navigation, above theconstellation. Typical off-the-shelf receivers can tracksignals down to a carrier-to-noise density, C/N0, of about35 dB Hz, which is marginally sufficient for trackingmain-lobe signals at geosynchronous altitudes using apatch antenna. Simulation studies of high-altitudenavigation performance predict that significant gains canbe made if signals can be used effectively down to 28 dBHz 3,4. This present work aims to track signals in the 1229 dB Hz range. This will enable the use of side-lobesignals received at geosynchronous altitudes using apatch antenna 5. Additional motivation for this workcomes from a desire to track transient weak signals thatoccur terrestrially during ionospheric scintillations 6.Such an ability would enable physicists to extract aI. INTRODUCTIONTracking algorithms allow a GPS receiver to maintainlock on the Doppler shift and the pseudo-random numberION GPS 2002, 24-27 September 2002, Portland, OR2539
greater amount of information from GPS soundings ofthe disturbed ionosphere.II. MODELS OF THE CORRELATION MEASUREMENTS AND THE SIGNAL DYNAMICSThe present work makes three main contributions. First,it develops a new fine acquisition algorithm. Thisprocedure is needed in order to provide the trackingalgorithm with accurate initial estimates of the carrierfrequency, carrier Doppler shift, code phase, and carrieramplitude. It starts with coarse estimates of the carrierDoppler shift and the code phase and uses a batchnonlinear filter to refine these estimates while solvingfor initial carrier phase and carrier amplitude. Thepaper’s second contribution is the development of acombined carrier- and code-tracking nonlinear Kalmanfilter. A Kalman filter is an optimal algorithm that isefficient for real-time implementation because of itsiterative-in-time nature. The third contribution is aBayesian adaptation of nonlinear Kalman filteringtechniques that deals effectively with the uncertainnavigation data bits when the carrier-to-noise density islow. The new fine acquisition and tracking algorithmsare tested using a high-fidelity simulation of the weakGPS signal that exits a receiver’s RF front end.The batch fine-acquisition algorithm and the signaltracking Kalman filter operate using dynamic models ofcarrier phase, code phase, and carrier amplitude andmeasurement models that give the relationship betweenthese signal quantities and the correlations that thereceiver uses to sense the signal. These models assumethat the sampled signal coming out of the receiver’s RFfront end takes the formyj A(τj)d[τj–ts(τj)] C[τj–ts(τj)] cos[ωIFτj-φ(τj)] nj(1)where yj is the measured RF front-end output at sampletime τj, A(τ) is the carrier amplitude, d(τ) is the 50 Hznavigation data bit stream of 1 values, C(τ) is the 1.023MHz C/A PRN bit stream of 1 values, ts(τ) is the PRNcode phase expressed as a relative code start time, ωIF isthe RF front-end’s intermediate image of the nominalGPS L1 carrier frequency, φ(τ) is the carrier phase thatresults from accumulated delta range, and nj is anelement of a zero-mean discrete-time Gaussian whitenoise sequence with a variance of σn2. The carrier-tonoise density of this sampled signal is C/N0 A2/(4σn2δτ), where δτ τj 1 - τj.Some of these algorithms require a special receiverarchitecture. The coarse and fine acquisitions must takeplace in a software receiver environment because theavailable algorithms for coarse acquisition of weaksignals 7 and the new algorithm for fine acquisition arebatch algorithms. The new tracking algorithm operates inan iterative causal manner using standard accumulations.It can be implemented in a software receiver or in astandard real-time receiver that uses a hardwarecorrelator. This paper envisions a system that acquiresan initial batch of data for coarse and fine acquisition.While it is processing this data, it continues to record astream of intermediate frequency data from its RF frontend for later use in tracking. After it has finished itscoarse and fine acquisition calculations, it initiates thenew weak signal tracking algorithm. It starts out trackingstored data and operates faster than real-time in order to“catch up”. Afterwards it continues to track a givensignal in real time. Systems that lack sufficient power todo batch-mode acquisition, interim RF bit storage, and“catch-up” tracking still could use the new trackingalgorithm. Its weak signal capabilities would be usefulfor a signal that suffered a power fade subsequent tohaving been acquired by standard techniques.Accumulation Measurement Models. The receiveraccumulates correlations between the yj data stream andreplicas of the code and carrier signals that it produces.These accumulations take the usual formsjk N kI k ( ) j j k 1y j C NCO (τ j t NCOk ) cos[ω IF τ j φ NCO (τ j )](2a)Q k ( ) jk N k j jk 1y j C NCO (τ j t NCOk ) sin [ω IF τ j φ NCO (τ j )](2b)where Ik( ) and Qk( ) are, respectively, the in-phase andquadrature accumulations with an “early” offset of ,which will be a late offset if 0. The function CNCO (τ)is the receiver's reproduction of the tracked PRN code.The subscript NCO is used because this functionsimulates a code numerically controlled oscillator. Thetime tNCOk is both the NCO's prompt code start time andthe start time of the accumulation interval. The initialsample jk 1 is the minimum value that respects the boundtNCOk τjk 1. Each accumulation spans one PRN codeperiod; therefore, the number of samples Nk is themaximum value that respects the limit τjk Nk tNCOk 1. Thefunction φNCO (τ) is the receiver's reproduction of thesignal’s carrier phase. The negative sign in front of itpresumes high-side mixing during one of the RF frontend’s down conversion stages.The remainder of this paper includes 4 main sectionsfollowed by conclusions. Section II presents models ofthe correlation measurements and the signal. Section IIIdescribes the fine acquisition algorithm. Section IVexplains the Kalman filter-based tracking algorithm.Section V presents simulation results for thesealgorithms. Section VI gives the conclusions.Equations (2a) and (2b) give a recipe for how theION GPS 2002, 24-27 September 2002, Portland, OR2540
receiver calculates its accumulations, but the fineacquisition algorithm and the Kalman filter need a modelof how these accumulations are related to the actualsignal parameters in eq. (1). This model takes the form:N AdI k ( ) k k m cos( φk ) R( tk ) n Ik(3a)2N A dQ k ( ) k k m sin( φ k )R ( t k ) nQk(3b)2where Ak is the average carrier amplitude over theaccumulation interval, dm is the navigation data bit, φk isthe interval average of the carrier phase errorφ(t)–φNCO (t), tk ts(tmidk)-tmidk is the code phase error at theinterval mid-point tmidk (tNCOk tNCOk 1)/2, R(t) is theautocorrelation function of the PRN code, and nIk and nQkare samples of zero-mean, uncorrelated Gaussiandiscrete-time white noise sequences, both with varianceequal to Nkσn2/2. The model of R(t) that gets used has theslope discontinuities of its triangular peak rounded off bycubic splines. These cause d2R/dt2 to be continuous,which avoids problems in the gradient-based numericaloptimizations that are part of the fine acquisitionalgorithm and the Kalman filter. This modification ofR(t) is reasonable because the limited bandwidth of theRF front end rounds off the actual correlation’s sharpcorners.ym cos( φm )R ( tm ) sin( φm )R ( tm )1ηeml 1η eml cos( φm )Reml ( tm) n ym sin( φm) Reml ( tm) dm hm( φm, tm, Am ) nym(5)where Am is the average carrier amplitude over the bitinterval, φm is the average carrier phase error over thebit interval, tm ts(tmidm)-tmidm is the code phase error atthe mid-point of the bit interval tmidm (tNCOkm tNCOkm 20)/2,Reml( t) R( t eml/2)-R( t- eml/2) is the early-minuslate correlation function, and nym is a sample from a zeromean, uncorrelated Gaussian discrete-time white noisevector sequence with covariance equal to the 4 4 identitymatrix. Equations (4) has been specifically designed inorder to achieve this normalization of the nymmeasurement error covariance.When performing optimal estimation it is best to modelthe raw measurements and their errors statistically and touse those models directly in the design of the estimator.One should avoid ad-hoc processing of measured dataprior to statistical analysis, which is why discriminatorsare not used here. This paper’s approach, however, is inslight violation of this principle. The rawest form of themeasurements is given in eq. (1). The new algorithms donot work directly with eq. (1) because their optimizationbased techniques would require many re-calculations ofcorrelations, which would be inefficient. Also, the PRNcode function C(t) is not everywhere differentiable,which would cause problems for these gradient-basedoptimization procedures. This "impure" approach doesnot cause significant loss of performance because of thefollowing facts: The measurements in eqs. (2a), (2b),and (4) capture almost all of the important informationabout the signal, the measurement models in eqs. (3a),(3b), and (5) are faithful, and the noise in thesemeasurements remains Gaussian.The Kalman filter make use of prompt and early-minuslate accumulations that are summed over entirenavigation data bit periods. These form the 4 1measurement vectork m 19 I k (0) k k m k m 19 Qk (0) 12 k k m ym σ n N m 1 k m 19 [ I k ( eml / 2) I k ( eml / 2)] ηeml k km 1 k m 19 [Qk ( eml / 2) Qk ( eml / 2)] η eml k km (4)where m is the index of the navigation data bit interval, Nm (Nkm Nkm 1 Nkm 19) is the number of RF front-endsamples in the data bit interval, km is the index of the firstPRN code period of the data bit interval, eml is the timeoffset between the early and late versions of thereceiver’s PRN code reconstruction, and ηeml 2[1-R( eml)] is a constant. Equation (4) gives a recipe forhow to compute ym from receiver data, but the Kalmanfilter also requires a model of how ym is related to signalparameters. It takes the form:ION GPS 2002, 24-27 September 2002, Portland, ORAm d mσn Nm 2 Carrier Phase, Code Phase, and Carrier AmplitudeDynamic Models. The carrier phase dynamics modeltakes the form of a discrete-time triple integrator drivenby discrete-time white noise: δtk2 1δt xk δ tk xφ 2 φ x 0 ω 01δtxk ω NCOk ω xα k 1 0 01 xα k 0 1000 0 1 0 0 wφk(6) 0 0 1 0 2541
where xφk φ(tNCOk) – φNCO (tNCOk) is the difference betweenthe true carrier phase and the receiver NCO’s carrierphase at the start of an accumulation interval, xωk is thetrue carrier Doppler shift at the start of the interval, xαk isthe rate of change of the carrier Doppler shift at the startof the interval, ωNCOk is the receiver NCO’s reconstructedcarrier Doppler shift for the interval, δtk tNCOk 1 - tNCOk isthe length of the interval, and wφk is a member of a zeromean, discrete time Gaussian white-noise vectorsequence. The states of this model can be used to deducethe average carrier phase error over the accumulationinterval: x δtk δtk2 φ δt φk 1 xω - k ωNCOk [0 0 0 1] wφk26 2 xα kThe dynamic model of the PRN code phase propagatesthe code start time from the beginning of one codeperiod to the next:ω L1δ tnom [1 0 0 0 ]wφk(9)tsk 1 t sk wtskω L1 xωk 0.5δ tnom xα kwhere tsk and tsk 1 are the true start and stop times of thePRN code period in question. This model includescarrier aiding via the second term on the right-hand sideof eq. (9). The time increment δtnom is the nominal codeperiod, 0.001 sec. The scalar wtsk is a white-noisesequence that effectively models code-carrierdivergence as a random walk. Its variance is E{wtsk 2} δtkqts, where qts is the random walk intensity. The PRNcode start/stop times in eq. (9) can be used to calculatethe code phase error tk that is used in eq. (3a) and (3b): tk (tsk 1 tsk )/2 - tmidkω L1δtnom [1 0 0 0]wφk tsk 12ωL1 xωk 0.5δtnom xαk(7)This quantity is needed in eqs. (3a) and (3b). A similarcarrier phase model applies when accumulations getsummed over bit intervals, as in eq. (4). In this case, thek index of the PRN code period in eqs. (6) and (7)changes to the m index of the data bit period, and thenominal accumulation interval increases from δtk 0.001 sec to δtm 0.020 sec. δtk3 3 δtk2 2 2δ tk2 δt k 2 SgωL1 00 32 δtk 8 δtk 6 δt k 0 SfωL12 0 δt k 20000Ak and Ak 1 are the carrier amplitudes at the times tNCOk andtNCOk 1. The scalar wAk is the zero-mean, white-noiseGaussian sequence that drives the random walk. Itsvariance is E{wAk2} δtkqA, where qA is the random walkintensity. The amplitudes in eq. (11) can be used tocompute the average amplitude that is needed in eqs. (3a)and (3b):Ak ( Ak 1 Ak ) / 2 Ak 0.5 w Ak(12)δtk3 6Equations (11) and (12) can be modified to propagate thecarrier amplitude from one navigation data bit start timeto the next. The only necessary change is to switch fromthe PRN code period index k to the data bit index m.0 δtk3 8 0 δ tk2 6 00 0 δ tk3 20 0 δ tk 2 00 00 0 δt k 3 III. BATCH FINE ACQUISITION ALGORITHMThe Kalman filter tracking algorithm of this paper needsan accurate initialization in order to function properly ina weak signal environment. The initialization procedurestarts with the usual search in code-phase/carrierDoppler space. For purposes of the present paper, thissearch is called the coarse acquisition. A suitable batchcoarse acquisition algorithm that works well in a weaksignal environment is described in Ref. 7. This algorithmgives Doppler shift to an accuracy of about 6 Hz andcode phase to an accuracy of 0.1 chips.(8)The quantity qLOS is the intensity of the accelerationrandom walk. Sg and Sf are the receiver clock's frequencyand phase random walk intensities, respectively 9, and ωL1is the nominal L1 carrier frequency.ION GPS 2002, 24-27 September 2002, Portland, ORwtsk -tmidk (10)A random walk is used to model the dynamics of thecarrier amplitude:Ak 1 Ak w Ak(11)δtk5 72 δtk3 3 δ tk2 2 δtk4 30 δ tk2 2δtkδtk3 24 43δ tk 30 δtk 24 δtk5 252 δtk4 812A model similar to eqs. (9) and (10) applies forpropagation of the navigation data bit start/stop timesfrom one bit to the next. In this case, the time index kgets replaced by the time index m, and the interval δtnomincreases to 0.020 sec.The covariance of the wφk white-noise sequence is acombination of three terms. One models a random walkacceleration of the line-of-sight (LOS) vector to the GPSsatellite, as in Ref. 8, and the other two model the effectsof receiver clock phase random walk and frequencyrandom walk, as in Ref. 9. The wφk covariance takes theform:TPwφwφk E{ wφk wφk } δtk5 20 4δt 8qLOS k3 δt 6 5k δtk 72 2542
J (k 0 ) N1 m 0 m The Kalman filter tracker needs the initial Doppler shiftto an accuracy about 0.25/τPLL Hz, where τPLL is thecharacteristic settling time of the filter. This time can beon the order of 1 sec when tracking a weak signal withslow dynamics. In addition, the Kalman filter needs toknow the PRN code period that corresponds to the starttime of a navigation data bit, it needs an initial carrierphase estimate that is accurate to within 45o, and itneeds initial estimates of the carrier amplitude and therate of change of Doppler shift. The Kalman filter alsoneeds a covariance matrix for the initial estimationerrors in its 5 states.M 2 (13)where M is the maximum integer that respects the boundM (Nacq-38)/20. M 1 is the number of full data bitintervals. The integer Nm is the total number of RF frontend samples during a proposed navigation data bitinterval, as in eq. (4). It normalizes the differentsummations.The bit start-stop determination is a brute-forceoptimization of J(k0). All 20 possible values of k0 aretried, and the value that yields the largest J(k0) getschosen as the correct bit start/stop index. If M is largeenough, then this method is very likely to yield thecorrect data bit start/stop time. If M is not very large,then one can increase M after the acquisition by using theaccumulations from the subsequent tracking interval.This approach further increases the likelihood of gettingthe correct value of k0, which is needed to calculate codepseudorange.A sequence of fine acquisition calculations is used todetermine suitable values for the Kalman filterinitialization. The first calculation determines thenavigation data bit start time. The second calculationestimates initial values for the three carrier phase states,xφ, xω, and xα, and it uses these estimates to get a firstestimate of the carrier amplitude, A.The finalcalculation refines the estimates of these 4 quantitiesalong with the estimate of the initial code phase, ts.Open-Loop Accumulation Data Used in Fine AcquisitionCalculations. The first step in the fine acquisition is tocalculate an “open-loop” time history of 1000 Hzaccumulations. This computation uses the Doppler shiftand code phase estimates from the coarse acquisition,ωNCO and tNCO0, to specify an “open-loop” carrier phasetime history via the formula φNCO (t) ωNCO t and an “openloop” code phase time history via iteration of thedifference equation tNCOk 1 tNCOk 0.001ωL1/(ωL1 ωNCO ).These phase time histories are used to generate mulations using the formulas in eqs. (2a) and (2b).Let the prompt accumulations be Ipk Ik(0) andQpk Qk(0), and let the early-minus-late accumulationsbe Iemlk Ik( eml/2) - Ik(- eml/2) and Qemlk Qk( eml/2) Qk(- eml/2). These quantities constitute the data that getused throughout the remainder of the fine acquisitioncalculations. They are calculated for k 0, , Nacq. Nacqis normally chosen to equal the number of PRN codeperiods used during the coarse acquisition.Batch Fine Acquisition Cost Functions. The batchacquisition calculations make use of several differentcost functions that are related to the joint probabilitydensity function of the noise terms in eqs. (3a) and (3b)or (5). The ultimate cost function that the initialestimates xφ0, xω0, xα0, A0, and ts0 must minimize isJa(xφ0,xω0,xα0,A0,ts0) M N k2 A02 2 k 0 20m 19 1 221 IQRt( ) pkpkkσ n2 m 0 k k 20 m N k 4 0k 20 m 191 0 2σ n k k 20 m0M- ln(2 cosh{m 0[ A0 R ( tk ){I pk cos ( φ k ) Q pk sin ( φ k )}] M k0 20 m 191 σ n2η eml m 0 k k 0 20 m{ Determination of the Start Time of a Navigation Data Bit.The objective of the bit start time calculation is todetermine the value of k0, the PRN code period indexwhose start time is also the start time of the firstnavigation data bit. The true value k0 is an element of theset {0,1,2, ,19}. It can be determined by an integeroptimization 1.Mm 0})122I emlk Q emlkNk[N k2 A02 2Reml ( t k )4]}- ln(2 cosh{k0 20m 19 σ n2ηeml k k 20 m01[ A0 R eml ( tk ){I emlk cos ( φ k ) Qemlk sin ( φ k )} ]})(14)This cost function is a negative log likelihood function.In other words, Cexp[-Ja(xφ0,xω0, xα0,A0,ts0)] is the jointprobability density function for the four accumulationtime histories Ipk, Qpk, Iemlk , and Qemlk for k k0 to(k0 20M 19) conditioned on the 5 unknown parametersxφ0, xω0, xα0, A0, and ts0, where C is a normalizing constant.The ln[2cosh()] terms are the results of probabilitysummations over the uncertain data bits dm under theThe merit function that gets optimized is the signalpower after summing over a pre-detection interval of 1navigation data bit:ION GPS 2002, 24-27 September 2002, Portland, OR2k0 20m 19 I pk N1 Q pkm k k 20 mk k0 20m0 k0 20 m 192543
assumption that the 1 and –1 values are equally likely.The dependence of Ja on xφ0, xω0, and xα0 comes partlythrough the carrier phase model in eq. (6), but with wφkset to 0 because the batch filter does not considerprocess noise. This yields the model φk(xφ0,xω0,xα0) [xφ0 xω0(tmidk-tNCOk0) 0.5xα0(tmidk-tNCOk0)2 -ωNCO tmidk]. Thedependence of Ja on ts0 comes through the tk terms.They can be computed as functions of xω0, xα0, and ts0 byiterating eq. (9) with wφk and wtsk both set zero. Thisiteration generates the time series tsk (xω0,xα0,ts0) for k k0t o (k0 20M 19), and eq. (10) yields tk(xω0,xα0,ts0) [tsk 1(xω0,xα0,ts0) tsk (xω0,xα0,ts0)]/2 - tmidk.Jd(xφ0, xω0, xα0) M1 Nm 0 mM k0 20 m 191 σ n2 m 0 k k 20 m0k0 20 m 19 k k 0 20mM k0 20 m 19 m 0 k k 0 20mk0 20 m 19 k k 0 20m[I pk cos( φ k xφ 0 ) Q pk sin( φk xφ 0 )](18a)[I pk sin( φk xφ 0) Q pk cos( φ k xφ 0 )](18b)where the r subscript of the accumulations stands for"rotated". Neither of these accumulations depends on xφ0because φk-xφ0 is independent of xφ0, and the value of xφ0that minimizes Jd for given values of xω0 and xα0 isxφ0opt(xω0, xα0) M I Q M [I 2 Q 2 ] - 12 atan2 2 rm rm , rm rm (19) Nm m 0 N m m 0 [I pk cos( φ k ) Q pk sin( φk )][I pk cos( φk ) Q pk sin( φk )]]Qrm(xω0, xα0) A Sequence of Fine Acquisition Optimizations. The firststep in the fine acquisition algorithm is an approximateglobal minimization of Jc(xφ0,xω0,xα0). This proceduresearches for the minimum on a grid in (xω0,xα0) space andcomputes xφ0 at each grid point by using eq. (19). Figure1 explains why the phase xφ0 that minimizes Jd will be agood approximation to the phase that minimizes Jc. Thefigure plots the rotated accumulations Ipsin( φ) Qpcos( φ) on the vertical axis versus Ipcos( φ)Qpsin( φ) on the horizontal axis for a case in which xω0and xα0 are nearly correct. The minimization of Jc withrespect to xφ0 performs an additional rotation in order tomaximize the spread of these points along the horizontalaxis, which would tend to align the solid line with thehorizontal axis. The minimization Jd, rotates xφ0 in orderto minimize the spread of these points along the verticalaxis, which would tend to align the dash-dotted line withthe vertical axis. Thus, these two minimizations tend toproduce nearly the same optimal xφ0 rotation, which iswhy the eq. (19) value of xφ0 can be used to approximatelyoptimize Jc(xφ0, xω0, xα0). constant(15)The above cost function is quadratic in A0. It can beminimized by first minimizing the following costfunctionJc(xφ0, xω0, xα0) - [This cost function is a weighted sum of the norm squaredof the projection of the first two elements of themeasurement error vector nym of eq. (5) perpendicular tothe signal direction. The importance of this costfunction is that it can be minimized with respect to xφ0analytically, which reduces the dimensionality of thesearch grid over which one must seek the globalminimum. In order to minimize Jd with respect to xφ0,defineIrm(xω0, xα0) An approximation to the cost function in eq. (14) isuseful for estimating the carrier phase parameters xφ0, xω0,and xα0 and the carrier amplitude A0. It presumes that thecode phase is correct, which implies that tk 0, Reml( tk) 0, and R( tk) 1. This approximation makes the lasttwo summations of Ja constant, and therefore irrelevantto any optimization.An additional approximationsubstitutes the absolute value function for the ln[2cosh()]function. Use of the absolute value function amounts to“hard” bit estimation, in which the bit sign is set equal tothe sum over the data bit interval of Ipkcos( φk) –Qpksin( φk). The cost function approximation isJb(xφ0, xω0, xα0, A0) M N k2 A02 k 0 20m 19 1 221 IQ pkpkσ n2 m 0 k k 20 m N k 4 0 - A02 k 0 20m 19 I pk sin( φ k ) Q pk cos( φ k ) k k0 20 m (17)(16)with respect to xφ0, xω0, and xα0. One then computes A0 -2Jcmin/Ntot to optimize Jb in eq. (15), where Ntot Nk0 Nk0 1 Nk0 20M 19.Yet another cost function needs to be defined in order todeal with the fact that the cost function in eq. (16) hasmany local minima. A search on a grid has to beperformed in order to find the global minimum, and analternate cost function allows a reduction from a 3dimensional grid to a 2-dimensional grid. The alternatecost function isThe extent and spacing of the search grid in (xω0,xα0) spacemust be chosen carefully in order to get a good solution.ION GPS 2002, 24-27 September 2002, Portland, OR2544
The xω0 grid should be centered at the coarse acquisitionDoppler shift estimate, and it should extend in eitherdirection as far as the possible uncertainty in the coarseDoppler shift. A good rule of thumb is to extend for 25% of the pre-detection bandwidth of the coarseacquisition or 100% of the Doppler grid spacing of thecoarse acquisition, whichever is greater. For most of theexamples of this paper a range of 12.5 Hz has beenused. The xα0 grid should be centered at zero. Recall thatxα0 models acceleration. The extent of this grid should bechosen to reflect the range of possible LOSaccelerations. For a geosynchronous receiver, the LOSacceleration comes mainly from gravitation and rangesup to 0.081 g, which translates into xα0 26.3 rad/sec2(4.2 Hz/sec). The actual search should expand to includea factor of safety which ensures that the global minimumdoes not fall outside of the grid. A range of 33.6rad/sec2 ( 5.3 Hz/sec) has been used in the present study.does this using Newton's method. It starts with the 10lowest isolated minima of Jc[xφ0opt(xω0,xα0), xω0, xα0] at thegrid points of the previous calculation and performs aniterative Newton search from each point in order tooptimize Jc with respect to xφ0, xω0, and xα0. The guardedNewton search includes precautions that ensure globalconvergence to a local minimum, such as modification ofthe cost Hessian and adaptation of the search step size 10.The point that achieves the lowest minimum usingNewton's method yields fine estimates for the carrierphase, xφ0, the Doppler shift, xω0, and the Doppler shiftrate, xα0. Ten different initial guesses are used in order toincrease the likelihood that one of the first guesses willreach the true global minimum of Jc.Results from the first two steps of a typical fineacquisition are shown in Fig. 2. It depicts a contour ofJb(xφ0,xω0,xα0,A0) at grid points in (xω0,xα0) space with xφ0given by eq. (19) and A0 given by the optimal formula thatappears just after eq. (16). The cost Jb(xφ0,xω0,xα0,A0),when optimized with respect to A0, is monotonicallyrelated to Jc(xφ0,xω
carrier. Various designs exist for these tracking loops, e.g., see Refs. 1 and 2. The goal of the present work is to design a combined DLL/PLL code and carrier tracking loop that is effective at tracking GPS L1 Coarse/Acquisition signals (C/A) for very low carrier-to-noise ratios. This is to be accomplished without prior knowledge of the navigation
