Transcribed TextTranscribed Text

Reading Objectives In this project you’ll perform a navigation solution based on pseudorange data that was downloaded from a CORS receiver. Your solution will implement the least-squares techniques we discussed in lecture. You will account for ionospheric, tropospheric, earth rotation, and relativistic effects. Modeled Pseudorange Calculation Recall that in lecture we derived the following model for the pseudorange measurement: ρ1/2 = R + c(δtu − δts) + T + Iρ1/2 + Mρ1/2 + ερ1/2 (1) where • ρ1/2 pseudorange measured on L1/2 frequency based on code at the user’s signal receive time tR • R is the geometrical range from satellite s to user u based upon where the user is at tR and the satellite at (tR−δtu(tR)−tOWLT) with tOWLT being the one way light time defined as R = p (xs − xu) 2 + (ys − yu) 2 + (zs − zu) 2 • δtu is the user/receiver clock error (bias) at tR • δts is the satellite clock error at (tR − δtu(tR) − tOWLT) • T is the tropospheric delay • Iρ1/2 is the ionospheric delay in code measurement on L1/2 • Mρ1/2 is the multipath delay in code measurement on L1/2 • ερ1/2 other delay/errors in code measurement on L1/2 at tR 1 For signal transmission through a vacuum, the one way light time tOWLT is related to R by c · tOWLT = R + T + Iρ1/2 + Mρ1/2 = ||ru (tR) − rs (tR − δtR(tR) − tOWLT)|| + T + Iρ1/2 + Mρ1/2 (2) We call this an implicit definition for tOWLT because tOWLT shows up on both sides of the equality. At any particular tR, a GNSS receiver will produce multiple pseudorange measurements—one for each unique signal tracked. Each measurement can be modeled by (1). In preparation for least-squares estimation, we can cast (1) into following standard-form model for the ith pseudorange measurement: z 0 i = h 0 i (x) + w 0 i (3) where z 0 i = ρ, h 0 i (x) = R + c(δtu − δts) + T + Iρ + Mρ, w 0 i = ερ, and x = r T u , c · δtu T is the 4-by-1 unknown parameter vector that we wish to estimate. The primes in (3) indicate that we have not yet performed normalization to make the measurement noise have unity variance—we’ll do this later. Suppose we have nz different pseudorange measurements modeled as in (3), all taken at user’s receiver time tR but corresponding to different satellites. We can stack these together to obtain a vector measurement equation of the form z 0 = h 0 (x) + w 0 (4) Recall from lecture that we can approximate h 0 (x) by a first-order Taylor series approximation h 0 (x) ≈ h 0 (¯x) + H 0 (x − x¯) (5) where H 0 is the nz-by-4 Jacobian matrix for nz independent pseudorange measurements and ¯x = ¯r T R, c ¯δtu T is an approximate value (a guess) for x. Equation (5) is a better approximation the closer ¯x is to x. Recall from lecture that the ith row in H 0 is the 1-by-4 Jacobian ∂h0 i ∂x x=¯x =  ∂h0 i ∂ru , ∂h0 i ∂cδtR  x=¯x (6) with ∂h0 i ∂ru x=¯x ≈ r¯ T u (tR) − r T s (tR − ¯δtu − t¯OWLT) ||r¯u (tR) − rS (tR − ¯δtu − t¯OWLT)|| (7) which is just the transpose of the 3-by-1 unit vector pointing from the satellite (at the transmission event) to the receiver (at the receipt event), and ∂h0 i ∂cδtu x=x¯ ≈ 1 (8) In (7), t¯OWLT is also an approximation because it depends on ¯ru and ¯δtu through (2). Your first analysis task in the Analysis Tasks section below will be to derive (7) and (8) by taking partial derivatives (this was actually done in lecture). Our first step toward a full navigation solution will be to write a Matlab function named satpr which, given r¯u and c · ¯δtu, generates modeled pseudorange ¯ρi = h 0 i (¯x) and its corresponding 1-by-4 Jacobian H0 i for the ith SV. You can generate ¯ρi by application of (1) and (2) above. Your satpr function will need to account for each of the effects discussed previously. You will not be required to develop functions for modeling the ionospheric and neutral atmospheric delay (troposphere). Rather, you will invoke the functions getIonoDelay and getTropoDelay, which are found on the course website, within your satpr implementation. The getTropoDelay function will need a set of model coefficients stored in NMFcoeffs.mat, also found on the course website. 2 The getIonoDelay and getTropoDelay functions will require a conversion from ECEF coordinates to latitude, longitude, and altitude for which you can use the ecef2lla function you wrote in homework 1. An essential step in your satpr function will be to calculate t¯OWLT from ¯ru and ¯δtu. If (2) were a simple explicit expression with t¯OWLT only on one side, this would be easy. But t¯OWLT depends on rS (tR − ¯δtu(tR) − t¯OWLT), which in turn depends on t¯OWLT, so solving for t¯OWLT is a bit more tricky. We’ll take an iterative approach, as follows: Given tR, ¯ru (tR), and ¯δtu(tR), find t¯OWLT: 1. Start out with a rough guess of t¯OWLT, say t¯OWLT = 0.075 seconds. That’s about how long a signal takes to reach the earth from a satellite in GPS orbit. 2. Calculate a guess of the true time of the transmission event: t¯= tR − ¯δtu(tR) − t¯OWLT. 3. Plug t¯ into your satloc function to get the satellite position rS (t¯). 4. Obtain a better guess for t¯OWLT by applying (2): t¯OWLT = 1 c ||r¯u (tR) − rS (t¯)|| 5. Repeat steps 2-4 until rS (t¯) changes by less than a centimeter from one iteration to the next. No more than 10 iterations should be required. On the last iteration of this procedure, t¯ will be a very good guess of t, the true time of the transmission event. To complete your calculation of ¯ρ, you’ll need to evaluate the satellite clock error δS at t¯. Recall that the GPS control segment broadcasts a model for calculating δS(t). Relativistic Effects and Group Delay Differential The satellite clock delay δtS(t) can be modeled in simplified form as δtS(t) = af0 + af1 (t − tc) + af2 (t − tc) 2 (9) where t represents true time (equivalent to GPS system time for our purposes), the afi , i = 0, 1, 2, are polynomial clock model coefficients, and tc is the time of clock epoch. In this homework we’ll augment the model for δtS to include the group delay differntial tGD and a residual relativistic correction δtrel: δtS(t) = af0 + af1 (t − tc) + af2 (t − tc) 2 + δtrel − tGD (10) You can get tGD directly from satdata.TGD. Calculation of δtrel was discussed in lecture. Here, tc is the time of clock epoch, which is stored in the wc and tc fields of satdata. Having calculated δS at the transmission event, you now have all the ingredients necessary to generate a modeled pseudorange: ρ¯ = h 0 (¯x) = ||¯ru (tR) − rS (t¯)|| + c ¯δtu(tR) − δtS(t¯) (11) where t¯= tR − ¯δtu(tR) − t¯OWLT comes from the last iteration of the procedure for finding t¯OWLT. With this preparation, proceed to writing a satpr function that conforms to the following interface: function [rhoBar, Hrho] = ... 3 satpr(gpsWeekRx, gpsSecRx, cdtRx, rRxEcef, sd, ionodata, tiFlags) % satpr : Calculate the modeled pseudorange between the satellite specified in % sd and the receiver located at rRxEcef at the given time of signal % reception. % % % INPUTS % % gpsWeekRx -- GPS week of signal reception event in receiver time. % % gpsSecRx --- GPS seconds of week of signal reception event in receiver % time. % % cdtRx ------ Receiver clock error scaled by the speed of light: cdtRx = % c*dtRx. True time t is related to receiver time tRx by t = tRx % - dtRx. % % rRxEcef ---- 3-by-1 ECEF coordinates of the receiver antenna’s phase center, % in meters. % % sd --------- Ephemeris structure array for a single SV. Let ii be the % numerical identifier (PRN identifier) for the SV whose location % is sought. Then sd = satdata(ii) has the following fields: % % SVID - satellite number % health - satellite health flag (0 = healthy; otherwise unhealthy) % we - week of ephemeris epoch (GPS week, unambiguous) % te - time of ephemeris epoch (GPS seconds of week) % wc - week of clock epoch (GPS week) % tc - time of clock epoch (GPS seconds of week) % e - eccentricity (unitless) % sqrta - sqrt of orbit semi-major axis (m^1/2) % omega0 - argument of perigee (rad.) % M0 - mean anomaly at epoch (rad.) % L0 - longitude of ascending node at beginning of week (rad.) % i0 - inclination angle at epoch (rad.) % dOdt - longitude rate (rad / sec.) % dn - mean motion difference (rad / sec.) % didt - inclination rate (rad / sec.) % Cuc - cosine correction to argument of perigee (rad.) % Cus - sine correction to argument of perigee (rad.) % Crc - cosine correction to orbital radius (m) % Crs - sine correction to orbital radius (m) % Cic - cosine correction to inclination (rad.) % Cis - sine correction to inclination (rad.) % af0 - 0th order satellite clock correction (s) % af1 - 1st order satellite clock correction (s / s) % af2 - 2nd order satellite clock correction (s / s^2) % TGD - group delay time for the satellite (s) % % ionodata -- Ionospheric data structure array with the following fields: % % alpha0, alpha1, alpha2, alpha3 - power series expansion coefficients % for amplitude of ionospheric TEC % % beta0, beta1, beta2, beta3 - power series expansion coefficients for 4 % period of ionospheric plasma density cycle % % tiFlags --- 2-by-1 vector of flags that indicate whether to enable % tropospheric delay correction (tiFlag(1)) or ionospheric % delay correction (tiFlags(2)). % % % OUTPUTS % % rhoBar ---- Modeled pseudorange between the receiver at time of the signal % reception event (the time specified by gpsWeekRx and gpsSecRx) % and the satellite at time of the signal transmission event, in % meters. % % Hrho ------ 1 x 4 Jacobian matrix containing partials of pseudorange with % respect to rRxEcef and c*dtRx % %+------------------------------------------------------------------------------+ % References: % % % Author: %+==============================================================================+ Earth Rotation During Signal Propagation Interval In this homework we will calculate tOWLT taking into account the rotation of the earth during the interval of time between the signal transmission and reception events. You’ll do this by transforming both events into the same reference frame, i.e. the ECEF frame at the receipt event, as discussed in lecture. To test satpr and for further analysis, modify the viewsats script that you wrote in the previous lab as follows. After retrieving navigation data and plotting the SVs above the elevation mask angle at a specified time and location, call satpr for each SV above the elevation mask angle. Form a complete H 0 matrix by stacking the individual 4-by-1 Jacobian row vectors that satpr calculates for each visible SV. Verify your satpr function by running your modified viewsats script with the following inputs, assuming that gpsWeek and gpsSec refer to receiver time. You’ll find the getAntLoc.m function on Canvas. Pseudorange Model Verification Verify the accuracy of your satpr function by running your viewsats.m script with the following inputs, assuming that gpsWeek and gpsSec refer to receiver time. gpsWeek = 1712; gpsSec = 143800; elMask = 15; rRx = getAntLoc(’WRW0’); cdtRx = 0; You should find that the visible SVs are [2, 4, 8, 9, 17, 20, 24, 28]. The corresponding modeled pseudorange vector 5 ρ¯ = h 0 (x¯) should be within a centimeter or so of rhoBarVec = 23563949.9389694 20489373.1818216 23420865.0801816 22397464.3850917 21107992.5368601 23961491.4049471 22821307.9488751 21242700.8789888 and the H0 matrix should be nearly equivalent to Hprime = 0.57487 0.65411 0.49159 1 0.22441 0.97208 -0.068582 1 -0.11912 0.77048 0.62624 1 0.68363 0.043497 -0.72853 1 -0.040289 0.47493 -0.8791 1 -0.88517 0.19064 -0.42441 1 0.13052 0.87905 0.45852 1 -0.49349 0.76644 -0.41114 1 To help you diagnose problems you might have in your pseudorange modeling code, compare your intermediate values for the calculation of the modeled pseudorange corresponding to SV ID (PRN) 2 with the following values. Your values should not differ from the ones below by more than a centimeter or so. 1. If you’ve retrieved the right ephemeris data, your satdata structure for SV ID (PRN) 2 should be equal to satdata(2) = SVID: 2 health: 0 we: 1712 te: 144000 wc: 1712 tc: 143999.999996647 e: 0.011546980706 sqrta: 5153.7800827 omega0: -2.7277551532 M0: 2.33821189593 L0: 2.27738442636 i0: 0.938581804728 dOdt: -8.17462622028e-009 dn: 5.09949812885e-009 didt: -6.03953728526e-010 Cuc: -2.14762985706e-006 Cus: 8.45082104206e-006 6 Crc: 207.46875 Crs: -40.28125 Cic: -2.77534127235e-007 Cis: -3.91155481338e-008 af0: 0.000407380517572 af1: 1.13686837722e-012 af2: 0 TGD: -1.76951289177e-008 2. The receiver position ru expressed in latitude (rad), longitude (rad), and altitude (meters) should be lat = 0.528616927215839 lon = -1.70581017303855 alt = 164.932940838858 3. The ionospheric delay from the broadcast ionospheric model should be Iρ = 10.6903434747677 meters. 4. The tropospheric delay from getTropoDelay should be T = 6.44739205828992 meters. 5. The total time of flight multiplied by the speed of light should be c · tOWLT = 23686078.9649843 meters. 6. The satellite position at the time of transmission expressed in the ECEF frame of the receipt epoch should be rS = [−14358424.8331188, −20955469.9934282, −8445935.18315543]T meters. 7. The initial guess of receiver clock error should be c · δtu = 0 meters. 8. The relativistic SV clock correction should be c · δtrel = −5.8173636919967 meters. 9. The full SV clock correction (including the relativistic correction) should be c · δts = 122129.026014883 meters. Navigation Solution Now that we have a function for estimating ¯ρi = h 0 i (¯x) and generating the 1-by-4 Jacobian H0 i for the ith SV, we can proceed to the full least-squares solution of (4). Implement the steps outlined in the lectures on least squares estimation in a function called performNavigationSolution having the following interface: function [solutionMat,tRVecSolution,tTrueVecSolution,... residualVec,badSvIdVec,thetaNominalMat] = ... performNavigationSolution(tRVec,svIdVec,obsValidMat,prMat,fDMat,thetaMat,... iiStart,iiStop,rRxAppx, ... satdata,ionodata,svIdAllow,lambda,tiFlags) % performNavigationSolution : Perform the navigation solution using nonlinear % least squares techniques. % % % INPUTS % % tRVec --------------- Structure of two Nt-by-1 vectors that jointly define % the unique receiver time instants corresponding to all % received observables; tRVec.w contains the GPS week % and tRVec.s contains the GPS seconds of week. % % svIdVec ------------- Nsv-by-1 vector of unique SVIDs corresponding to 7 % SVs that were tracked at some point during the % measurement interval spanned by tRVec. % % obsValidMat --------- Nt-by-Nsv matrix whose (ii,jj)th element indicates % whether (1) or not (0) observables were valid for SVID % svIdVec(jj) at time tRVec.s(ii). % % prMat --------------- Nt-by-Nsv matrix whose (ii,jj)th element is the % measured pseudorange value for SVID svIdVec(jj) at % time tRVec.s(ii), in meters. % % fDMat --------------- Nt-by-Nsv matrix whose (ii,jj)th element is the % measured Doppler value for SVID svIdVec(jj) at time % tRVec.s(ii), in Hz. % % thetaMat ------------ Nt-by-Nsv matrix whose (ii,jj)th element is the % measured beat carrier phase value for SVID svIdVec(jj) % at time tRVec.s(ii), in cycles. % % iiStart,iiStop ------ Start and stop indices into tRVec between which a % navigation solution will be computed. These indices % are used to isolate a specific section of data. The % number of solutions Ns will nominally be Ns = iiStop - % iiStart + 1 but will be less if for some epochs a % solution could not be computed (e.g., too few SVs). % % rRxAppx ------------- 3-by-1 approximate location of receiver antenna, in % ECEF meters. % % satdata, ionodata --- SV and ionospheric navigation data parameters. % See getephem.m for details. % % svIdAllow ----------- Vector of SVIDs for SVs that are allowed to % participate in the navigation solution. % % lambda -------------- Nominal wavelength of GNSS signal, in meters. % % tiFlags ------------- 2-by-1 vector of flags that indicate whether to enable % tropospheric delay correction (tiFlag(1)) or % ionospheric delay correction (tiFlags(2)). % % % OUTPUTS % % solutionMat --------- Ns-by-4 matrix of navigation solutions whose rows % are formatted as [xEcef,yEcef,zEcef,c*deltR]. % % tRVecSolution ------- Structure of two Ns-by-1 vectors that jointly define % the unique receiver time instants corresponding to % navigation solutions in solutionMat; tRVecSolution.w % contains the GPS week and tRVecSolution.s contains the % GPS seconds of week. % % tTrueVecSolution ---- Structure of two Ns-by-1 vectors that jointly define % the unique receiver time instants corresponding to % navigation solutions in solutionMat. This vector is 8 % an estimate of ’true’ GPS time for each % solution. tTrueVecSolution.w contains the GPS week and % tTrueVecSolution.s contains the GPS seconds of week. % % residualVec --------- Ns-by-1 vector of maximum least-squares residuals % magnitudes from each solution epoch, in meters. % % badSvIdVec ---------- Ns-by-1 vector of SVIDs corresponding to the SV with % the worst residual at each solution epoch. % % thetaNominalMat ----- Ns-by-Nsv vector of nominal beat carrier phase values, % in the same arrangement as those in thetaMat, based % solely on the SV-to-rRxAppx range evaluated at % receiver time, in cycles. % %+------------------------------------------------------------------------------+ % References: % % % Author: %+==============================================================================+ At each measurement epoch (each unique receiver time tR at which pseudorange measurements are available), you’ll take the following steps within this function: 1. Start with a guess for x = r T u , c · δtu T . Call this ¯x. Usually it’s sufficient to let your guess for ru be anywhere on the correct continent and your guess for c · δtu be zero. 2. Stack all pseudorange measurements that are valid at epoch tR into the vector z 0 . You’ll draw these measurements from the appropriate row of prMat. Keep track of which SV corresponds to each element of z 0 . Be sure to exclude pseudoranges corresponding to entries marked 0 in obsValidMat. 3. For the ith element (ith pseudorange measurement) in z 0 , i = 1, . . . , nz, find the corresponding predicted pseudorange ¯ρi and the 1-by-4 Jacobian H0 i by making a call to satpr. Stack the ¯ρi , i = 1, 2, . . . , nz into the nz-by-1 vector ¯z 0 and the Jacobians H0 i , i = 1, 2, . . . , nz into the nz-by-nx matrix H0 , where nx = 4 is the dimension of the parameter vector x. Your pseudorange measurements can now be approximately modeled as z 0 = ¯z 0 + H 0 (x − x¯) + w 0 which can be rearranged as zˇ 0 = H 0 x + w 0 where ˇz 0 , z 0 − z¯ 0 + H 0 x¯. Note that the rearranged equation is in the standard form that is the starting point for linear least squares. 4. Assume that all the pseudorange measurements are independent and identically distributed (i.i.d.), with wρ ∼ N (0, σ2 ρ ), so that w 0 can be modeled as w 0 ∼ N (0, P w0 ), with P w0 = σ 2 ρ I. (Note that this makes for an especially simple Cholesky-factorized weighting matrix R −T a , where R T a Ra = P w0 .) For now, set σρ = 2 meters. 5. Normalize the measurement equation to obtain z = Hx + w where z = R −T a zˇ 0 , H = R −T a H 0 , and w = R −T a w 0 . 9 6. Proceed as discussed in lecture to solve for the x that would minimize the cost function J(x) = kHx−zk 2 . This involves QR factorization of H. Call the solution you obtain ¯xnew. If kx¯new − x¯k is very small (say, less than 10−4 meters), then stop, declaring ¯xnew to be your solution at time tR; if not, set ¯x = ¯xnew and go to step 3. Note that, within this function, you will call satpr for each allowed SV for which you have data at each epoch between tRxVec(iiStart) and tRxVec(iiStop), so your new satpr function will see plenty of use. For this lab, you’ll only perform a pseudorange-based navigation solution so you won’t need to use the inputs fDMat, thetaMat, ionodata, and lambda, but keep them in the interface as placeholders. More sophisticated solutions that you may implement may make use of these quantities. Also, just return an empty matrix [] for thetaNominalMat. All of the other return arguments should be properly populated. The jth element of the vector residualVec equals the element with the maximum magnitude (absolute value) in the vector ∆z 0 , z 0 − zˆ 0 on the last iteration of the jth solution epoch. The SV corresponding to this value is stored in badSvIdVec(j). Top-level Script It will be convenient to have a top-level script called topsol that performs the following steps. You can either write your own function or use the example topsol.m on Canvas. Even if you write your own topsol function, make sure your performNavigationSolution works correctly within the example topsol function, since this is how your code will be evaluated by the TA. Note that example topsol.m is slightly more sophisticated than the simple procedure below; it divides the input data into small processing chunks and refreshes the ephemeris and ionospheric data, and performs elevation masking, before processing each chunk. 1. Set an initial guess for the receiver position. 2. Draw in a channel.mat file and format the data as described in the interface to performNavigationSolution. Extract only the data corresponding to GPS L1 C/A signals and having valid pseudorange measurements. This data formatting task is best implemented as a separate function. You can either write this yourself or use the function prepareTimeHistory.m posted on the course website. Note that the vector tRVec in performNavigationSolution should be made up of offset receiver time (ORT) measurement time stamps. 3. Pull the first ORT from tRVec. This first ORT is very close to the true GPS time of the first valid pseudorange measurement. 4. Plot the satellites at the first ORT by a call to satmap assuming some elevation mask angle elMaskDeg. 5. Set up other input arguments necessary for a call to performNavigationSolution such as iiStart, iiStop, and lambda. If your initial guess for the receiver position is good to, say, 100 km or so, then set svIdAllow equal to the list of SVs returned by your call to satmap. Otherwise, if your initial guess for the receiver position is way off, then satmap can’t be trusted to provide you a reliable list of SVs above your mask angle. In this case, just set svIdAllow = [1:32];. 6. Exclude from svIdAllow any unhealthy SVs, as indicated by the satdata(i).health flag (0 = healthy; otherwise unhealthy). 7. Call performNavigationSolution. 8. Solve for the average receiver antenna position in ECEF coordinates assuming a stationary antenna. 9. Plot the horizontal offset of each individual solution from the average solution as a single dot with East along the x axis and North along the y axis, in meters. The processing for this type of horizontal displacement plot was described in Homework 1. 10 10. Plot as a function of time (or solution index) the vertical displacement from the mean, in meters. 11. Plot as a function of time the worst-case residual at each solution epoch as contained in the residualVec output of performNavigationSolution. 12. Plot as a function of time the SV identification number for the SV with the largest residual at each solution epoch as contained in the badSvIdVec output of performNavigationSolution. Data Collection You’ll be provided with the data already in a channel.mat format. My guess is that it will take a very long time to proces the entire data set. Limit the data to the first 4 hours of data for the anlysis tasks below. Analysis Tasks 1. By taking partial derivatives, show that (7) and (8) are true. What is neglected in the approximation in (8)? 2. Run your topsol script on the data set you collected, assuming an elevation mask angle of 0 degrees. Set the initial guess of cδtu to 0 and the initial guess of ru to rRxAppx = getAntLoc(’Fiction_Island’); which is the approximate ECEF location of somewhere far away. Include all plots generated by topsol in your report. Pro tip: Work at first with a short set of data—perhaps just a single epoch of data—until your solver begins functioning correctly. You can also set epochStride to N > 1 during debugging to process only every Nth epoch, which will speed up your processing. 3. The nonlinear least squares solution is surprisingly forgiving as regards your initial choice of rRxAppx. Investigate what happens to successive position estimates as you iterate within the nonlinear least squares solver (i.e., as you successively set ˆxnew = x ∗ after each iteration of your solver). Where does Fiction Island lie on a map? Print out the position component of first 6 iterations of your navigation solution at some solution epoch and include these in your report. Then plot the position components of the first 6 iterations in a 3-d plot in Matlab, showing how your solution moves with each iteration from Fiction Island to the correct location. Experiment with other values of rRxAppx. See if you can find a value of rRxAppx that causes the least squares solver to fail. Draw some conclusions about the “attractor region” of the nonlinear solver. Note that because you don’t have even an approximate guess for the true receiver location for this mystery data set, you’ll need to disable elevation masking on your first run by setting elMaskDeg = -90. Having obtained an approximate position guess from your first run, you can perform a second run with elMaskDeg ≥ 0 if you’d like to eliminate some low-elevation signals. 4. Calculate HDOP and VDOP for a single epoch of data within the 10-minute interval. These values should be representative of the HDOP and VDOP values over the whole interval, since the SV geometry doesn’t change much over 10 minutes. Calculate the standard deviation in the offset-from-mean of your solutions in the East, North, and Vertical directions. Compare these standard deviations to HDOP·σρ and VDOP·σρ, with σρ = 2 meters. What value of σρ best fits the data? 11 5. Run your topsol script again on the data set you collected, but this time impose an elevation mask angle high enough that you only see 4 SVs. Again calculate HDOP and VDOP for a single epoch of data within the 10-minute interval. Compare the DOP quantities from these runs against the corresponding quantities from the previous run with elevation mask angle of 0 degrees. 6. Isolate the pseudorange measurements for a particular measurement epoch. Using these, solve for x ∗ . Now add 1000 meters to every one of the pseudorange values and solve again for x ∗ . What is the difference in the two values of x ∗ ? Explain. 7. Perform a full navigation solution on the first 100 seconds of data from your data set. Examine the average size, as an absolute value, expressed in meters, of (1) the modeled ionospheric delay, (2) the modeled neutral atmospheric delay, (3) the relativistic correction δtrel, and (4) the change in modeled pseudorange with and without accounting for earth rotation during signal transit. Average over all pseudoranges at each solution epoch and over all solution epochs in the first 100 seconds of data. Make a bar plot showing the average size of each effect. Rank these four effects in terms of size, and therefore, importance. 8. Choose one of the L2C-capable SV that was visible throughout your data set. Extract the beat carrier phase and the pseudorange data corresponding to this SV for both the GPS L1 C/A signal (GenericType = 0) and whichever GPS L2C variant is being tracked, whether CM (GenericType = 1), CL (GenericType = 2), or CLM (GenericType = 3). From these data, form the pseudorange-based and carrier-phase-based delay estimates ˆIρ and ˆIφ at the L1 frequency as we discussed in lecture. Negate ˆIφ so that it trends in the same direction as ˆIρ. Then take the mean of the difference between ˆIρ and the negated ˆIφ to estimate the constant bias in the negated ˆIφ. Add this bias to the negated ˆIφ and plot the rough ˆIρ and the smooth (negated and unbiased) ˆIφ together on the same graph. The smooth line should trace through the mean of the rough line. Approximate the RMS noise on the smooth and the rough plots. In terms of RMS noise (units of meters), how much more precise is the carrier-phase-based delay estimate than the pseudorange-based delay estimate? Now repeat these steps for another L2C-capable SV that was visible throughout your data set. 9. Your ionospheric delay plots from Task 8 may be significantly biased from the true delay. You may even end up with a negative delay estimate, which as you know is impossible for a group delay (it would violate special relativity). Whatever bias is present in your delay plots is due to (1) a slight misalignment of the L1 C/A code and the L2C code at the transmitter (transmitter differential code bias, δρS) and (2) an excess delay in the receiver’s L1 signal path vs. its L2 signal path, or a sample train misalignment between the L1 and L2 portionas of the receiver’s RF front end (collectively called receiver inter-frequency code bias, δρR). Both biases are expressed in meters. If you knew the exact values of these biases, you could eliminate them by calculating ˆIρ as ˆIρ = γ1 (ρL2 − ρL1 − δρS − δρR) One way to estimate these biases is to compare the raw value of ˆIρ (without attempting to eliminate the biases) with a model-based estimate of Iρ. Pretend that the broadcast model embodied in getIonoDelay is perfectly accurate over the 4-hour interval of your data set. Plot the ionospheric delay predicted by this model for your chosen SVs on a copy of your plots from Task 8. Do the trends in this estimate match those of your dual-frequency-based estimates? Assuming the getIonoDelay model’s predictions are perfectly accurate, estimate a value for the sum δρS + δρR, one for each of the two L2C-capable satellites. Hint: You can get the ionospheric delay predicted by getIonoDelay for a chosen SV over the interval by selectively saving the getIonoDelay outputs during topsol processing. But processing CORS data for the whole data set interval takes a long time when the epoch spacing is small. You can shorten your wait by skipping every Nstride epochs in the data output by prepareTimeHistory.m. In the example topsol found on the course website, this is done by setting epochStride to Nstride. The modeled ionosphere will still look smooth even if you only process epochs spaced by, say, 30 seconds. 12 Wrap Up Prepare a formal laboratory report for this lab. As an appendix to this report, attach a hard copy of your Matlab satpr and performNavigationSolution functions and your topsol scr

Solution PreviewSolution Preview

These solutions may offer step-by-step problem-solving explanations or good writing examples that include modern styles of formatting and construction of bibliographies out of text citations and references. Students may use these solutions for personal skill-building and practice. Unethical use is strictly forbidden.

antId = upper(antId);

switch antId
case 'TXAU'
    % NGS CORS antenna. Type: TRM57971.00. Source of coordinates:
    rEcef = [-743774.390; -5460644.512; 3200347.710];
    vEcef = [-0.0120; -0.0003; -0.0031];
    epochYear = 2005.0;
case 'SAM2'
    % NGS CORS antenna. Type: TRM57971.00. Source of coordinates:
    rEcef = [-751970.424; -5463638.752; 3193399.042];
    vEcef = [-0.0123; -0.0003; -0.0033];
    epochYear = 2005.0;
case 'TXTA'
    % NGS CORS antenna. Type: TRM57971.00. Source of coordinates:
    rEcef = [-712254.169; -5450502.920; 3224462.604];
    vEcef = [-0.0114; -0.0024; -0.0027];
    epochYear = 2005.0;
case 'WRW0'
    % WRW primary antenna. Type: TRM57971.00. Source of coordinates: L1-only CDGPS
    % solution on the WRW0-TXAU baseline using the above TXAU coordinates.
    % Velocity copied from TXAU velocity.
    rEcef = [-741990.378; -5462227.718; 3198019.532];
    vEcef = [-0.0120; -0.0003; -0.0031];
    epochYear = 2005.0;
case 'WRW1'
    error('Not yet defined');
case 'AUT0'
    % CSR primary antenna. Type: TRM57971.00. Source of coordinates: Integrated
    % MGEX network solution.
    rEcef = [-741133.4180; -5456322.9775; 3208395.3863];
    vEcef = [-0.0120; -0.0003; -0.0031];
    epochYear = 2005.0...

By purchasing this solution you'll be able to access the following files:

for this solution

PayPal, G Pay, ApplePay, Amazon Pay, and all major credit cards accepted.

Find A Tutor

View available Aerospace Engineering Tutors

Get College Homework Help.

Are you sure you don't want to upload any files?

Fast tutor response requires as much info as possible.

Upload a file
Continue without uploading

We couldn't find that subject.
Please select the best match from the list below.

We'll send you an email right away. If it's not in your inbox, check your spam folder.

  • 1
  • 2
  • 3
Live Chats