\f4.5 Nonlinear Least Squares | 239 At a given instant, the receiver collects the synchronized signal from the ith satellite and determines its transmission time ti , the difference between the times the signal was transmitted and received. The nominal speed of the signal is the speed of light, c 299792.458 km/sec. Multiplying transmission time by c gives the distance of the satellite from the receiver, putting the receiver on the surface of a sphere centered at the satellite position and with radius cti . If three satellites are available, then three spheres are known, whose intersection consists of two points, as shown in Figure 4.16. One intersection point is the location of the receiver. The other is normally far from the earth's surface and can be safely disregarded. In theory, the problem is reduced to computing this intersection, the common solution of three sphere equations. Figure 4.16 Three Intersecting Spheres. Generically, only two points lie on all three spheres. However, there is a major problem with this analysis. First, although the transmissions from the satellites are timed nearly to the nanosecond by onboard atomic clocks, the clock in the typical low-cost receiver on earth has relatively poor accuracy. If we solve the three equations with slightly inaccurate timing, the calculated position could be wrong by several kilometers. Fortunately, there is a way to x this problem. The price to pay is one extra satellite. Dene d to be the difference between the synchronized time on the (now four) satellite clocks and the earth-bound receiver clock. Denote the location of satellite i by (Ai , Bi , Ci ). Then the true intersection point (x, y, z) satises r1 (x, y, z, d) = (x A1 )2 + (y B1 )2 + (z C1 )2 c(t1 d) = 0 r2 (x, y, z, d) = (x A2 )2 + (y B2 )2 + (z C2 )2 c(t2 d) = 0 r3 (x, y, z, d) = (x A3 )2 + (y B3 )2 + (z C3 )2 c(t3 d) = 0 r4 (x, y, z, d) = (x A4 )2 + (y B4 )2 + (z C4 )2 c(t4 d) = 0 (4.37) to be solved for the unknowns x, y, z, d. Solving the system reveals not only the receiver location, but also the correct time from the satellite clocks, due to knowing d. Therefore, the inaccuracy in the GPS receiver clock can be xed by using one extra satellite. Geometrically speaking, four spheres may not have a common intersection point, but they will if the radii are expanded or contracted by the right common amount. The 240 | CHAPTER 4 Least Squares system (4.37) representing the intersection of four spheres is the three-dimensional analogue of (4.35), representing the intersection point of three circles in the plane. The system (4.37) can be seen to have two solutions (x, y, z, d). The equations can be equivalently written (x A1 )2 + (y B1 )2 + (z C1 )2 = [c(t1 d)]2 (x A2 )2 + (y B2 )2 + (z C2 )2 = [c(t2 d)]2 (x A3 )2 + (y B3 )2 + (z C3 )2 = [c(t3 d)]2 (x A4 )2 + (y B4 )2 + (z C4 )2 = [c(t4 d)]2 . (4.38) Note that by subtracting the last three equations from the rst, three linear equations are obtained. Each linear equation can be used to eliminate a variable x, y, z, and by substituting into any of the original equations, a quadratic equation in the single variable d results. Therefore, system (4.37) has at most two real solutions, and they can be found by the quadratic formula. Two further problems emerge when GPS is deployed. First is the conditioning of the system of equations (4.37). We will nd that solving for (x, y, z, d) is ill-conditioned when the satellites are bunched closely in the sky. The second difculty is that the transmission speed of the signals is not precisely c. The signals pass through 100 km of ionosphere and 10 km of troposphere, whose electromagnetic properties may affect the transmission speed. Furthermore, the signals may encounter obstacles on earth before reaching the receiver, an effect called multipath interference. To the extent that these obstacles have an equal impact on each satellite path, introducing the time correction d on the right side of (4.37) helps. In general, however, this assumption is not viable and will lead us to add information from more satellites and consider applying Gauss-Newton to solve a least squares problem. Consider a three-dimensional coordinate system whose origin is the center of the earth (radius 6370 km). GPS receivers convert these coordinates into latitude, longitude, and elevation data for readout and more sophisticated mapping applications using global information systems (GIS), a process we will not consider here. Suggested activities: 1. Solve the system (4.37) by using Multivariate Newtons Method. Find the receiver position (x, y, z) near earth and time correction d for known, simultaneous satellite positions (15600, 7540, 20140), (18760, 2750, 18610), (17610, 14630, 13480), (19170, 610, 18390) in km, and measured time intervals 0.07074, 0.07220, 0.07690, 0.07242 in seconds, respectively. Set the initial vector to be (x0 , y0 , z0 , d0 ) = (0, 0, 6370, 0). As a check, the answers are approximately (x, y, z) = (41.77271, 16.78919, 6370.0596), and d = 3.201566 103 seconds. 2. Write a Matlab program to carry out the solution via the quadratic formula. Hint: Subtracting the last three equations of (4.37) from the rst yields three linear equations in the four unknowns x ux + y uy + zuz + d ud + w = 0, expressed in vector form. A formula for x in terms of d can be obtained from 0 = det[uy |uz |x ux + y uy + zuz + d ud + w], noting that the determinant is linear in its columns and that a matrix with a repeated column has determinant zero. Similarly, we can arrive at formulas for y and z, respectively, in terms of d, that can be substituted in the rst quadratic equation of (4.37), to make it an equation in one variable. 4.5 Nonlinear Least Squares | 241 3. If the Matlab Symbolic Toolbox is available (or a symbolic package such as Maple or Mathematica), an alternative to Step 2 is possible. Dene symbolic variables by using the syms command and solve the simultaneous equations with the Symbolic Toolbox command solve. Use subs to evaluate the symbolic result as a oating point number. 4. Now set up a test of the conditioning of the GPS problem. Dene satellite positions (Ai , Bi , Ci ) from spherical coordinates (, i , i ) as Ai = cos i cos i Bi = cos i sin i Ci = sin i , where = 26570 km is xed, while 0 i /2 and 0 i 2 for i = 1, . . . , 4 are chosen arbitrarily. The coordinate is restricted so that the four satellites are in the upper hemisphere. Set x = 0, y = 0, z = 6370, d = 0.0001, and calculate the corresponding satellite ranges Ri = A2 + Bi2 + (Ci 6370)2 and travel times ti = d + Ri /c. i We will dene an error magnication factor specially tailored to the situation. The atomic clocks aboard the satellites are correct up to about 10 nanoseconds, or 108 second. Therefore, it is important to study the effect of changes in the transmission time of this magnitude. Let the backward, or input error be the input change in meters. At the speed of light, ti = 108 second corresponds to 108 c 3 meters. Let the forward, or output error be the change in position ||( x, y, z)|| , caused by such a change in ti , also in meters. Then we can dene the dimensionless error magnication factor = ||( x, y, z)|| , c||( t1 , . . . , tm )|| and the condition number of the problem to be the maximum error magnication factor for all small ti (say, 108 or less). Change each ti dened in the foregoing by ti = +108 or 108 , not all the same. Denote the new solution of the equations (4.37) by (x, y, z, d), and compute the difference in position ||( x, y, z)|| and the error magnication factor. Try different variations of the ti 's. What is the maximum position error found, in meters? Estimate the condition number of the problem, on the basis of the error magnication factors you have computed. 5. Now repeat Step 4 with a more tightly grouped set of satellites. Choose all i within 5 percent of one another and all i within 5 percent of one another. Solve with and without the same input error as in Step 4. Find the maximum position error and error magnication factor. Compare the conditioning of the GPS problem when the satellites are tightly or loosely bunched. 6. Decide whether the GPS error and condition number can be reduced by adding satellites. Return to the unbunched satellite conguration of Step 4, and add four more. (At all times and at every position on earth, 5 to 12 GPS satellites are visible.) Design a Gauss-Newton iteration to solve the least squares system of eight equations in four variables (x, y, z, d). What is a good initial vector? Find the maximum GPS position error, and estimate the condition number. Summarize your results from four unbunched, four bunched, and eight unbunched satellites. What conguration is best, and what is the maximum GPS error, in meters, that you should expect solely on the basis of satellite signals