Position estimation using a network of a global-positioning receivers

[0001] This application claims priority under 35 U.S.C. 119(e) to Provisional U.S. Patent Application Serial No. 60/413,205, filed Sep. 23, 2002, the entire contents of which are incorporated herein by reference.

[0002] The present invention relates to estimating precise position of a stationary or moving object using multiple satellite signals and a network of multiple receivers. The present invention is particularly suited to position estimation in real-time kinetic environments where it is desirable to take into account the spatial distribution of the ionosphere delay.

[0003] Satellite navigation systems, such as GPS (USA) and GLONASS (Russia), are intended for accuracy self-positioning of different users possessing special navigation receivers. A navigation receiver receives and processes radio signals broadcast by satellites located within line-of-sight distance, and from this, computes the position of the receiver within a pre-defined coordinate system. However, for military reasons, the most accurate parts of these satellite signals are encrypted with codes only known to military users. Civilian users cannot access the most accurate parts of the satellite signals, which makes it difficult for civilian users to achieve accurate results. In addition, there are sources of noise and error that degrade the accuracy of the satellite signals, and consequently reduce the accuracy of computed values of position. Such sources include carrier ambiguities, receiver time offsets, and atmospheric effects on the satellite signals.

[0004] The present invention is directed to increasing the accuracy of estimating the position of a rover station in view of carrier ambiguities, receiver time offsets, and atmospheric effects.

[0005] The present invention relates to a method and apparatus for determining the position of a receiver (e.g., rover) with respect to the positions of at least two other receivers (e.g., base receivers) which are located at known positions. The knowledge of the precise locations of the at least two other receivers (located at known positions within the coordinate system) makes it possible to better account for one or all of carrier ambiguities, receiver time offsets, and atmospheric effects encountered by the rover receiver, and to thereby increase the accuracy of the estimated receiver position of the rover (e.g., rover position).

[0006] In a first aspect of the present invention, an exemplary apparatus/method comprises receiving the known locations of a first base station and a second base station, obtaining a time offset representative of the time difference between the clocks of the first and second base stations, and obtaining measured satellite data as received by the rover, the first base station, and the second base station. The measured satellite data comprises pseudo-range information. The exemplary apparatus/method generates a first set of residuals of differential navigation equations associated with a set of measured pseudo-ranges related to a first baseline (R-B1) between the rover and the first base station. The residuals are related to the measured satellite data received by the rover station and the first base station, the locations of the satellites, and the locations of the rover station and the first base station. The exemplary apparatus/method also generates a second set of residuals of differential navigation equations associated with a set of measure pseudo-ranges related to a second baseline (R-B2) between the rover and the second base station. These residuals are related to the measured satellite data received by the rover station and the second base station, the locations of the satellites, and the locations of the rover station and the second base station. The exemplary apparatus/method estimates the rover's location from the first set of residuals, the second set of residuals, the time offset between the clocks of the first and second base stations.

[0007] In a second aspect of the present invention, the position information and measured satellite data from a third base station are received. Also, a time offset representative of the time difference between the clocks of the first and third base stations is obtained. Thereafter, the exemplary apparatus/method also generates a third set of residuals of differential navigation equations associated with a set of measured pseudo-ranges related to a third baseline (R-B3) between the rover and the third base station, and estimates the rover's location further with the third set of residuals and the time offset between the clocks of the first and third base stations.

[0008] In a third aspect of the present invention, which may be applied with any of the aspects of the present invention described herein, the exemplary apparatus/method generates the time offset(s) for the base stations from the measured satellite data of the base stations provided to it.

[0009] In a fourth aspect of the present invention, which may be applied with any of the aspects of the present invention described herein, the exemplary apparatus/method obtains measured satellite carrier phase data as received by the rover and the first base station (the first base line). The exemplary apparatus/method generates a fourth set of residuals of differential navigation equations associated with a set of measured carrier phase data related to the first base line, and resolves the cycle ambiguities in the carrier phase data from the fourth residual and one or more of the first, second, and third residuals. The resolved cycle ambiguities may take the form of floating ambiguities, fixed-integer ambiguities, and/or integer ambiguities. The exemplary apparatus/method estimates the rover's location further with the fourth set of residuals and the resolved cycle ambiguities associated with the first base line.

[0010] In a fifth aspect of the present invention, which may be applied with any of the aspects of the present invention described herein, the exemplary apparatus/method obtains measured satellite carrier phase data as received by the second base station, and further obtains a set of cycle ambiguities related to a set of satellite phase measurements associated with the baseline between the first and second base stations (B1-B2). The exemplary apparatus/method generates a fifth set of residuals of differential navigation equations associated with a set of measured carrier phase data related to the second base line between the rover and the second base station (R-B2), and resolves the cycle ambiguities in the carrier phase data from the fifth residual, the set of cycle ambiguities associated with the baseline between the first and second base stations, and one or more of the first, second, third, and fourth residuals. The resolved cycle ambiguities may take the form of floating ambiguities, fixed-integer ambiguities, and/or integer ambiguities. The exemplary apparatus/method estimates the rover's location further with the fifth set of residuals and the resolved cycle ambiguities associated with the second base line (R-B2).

[0011] In a sixth aspect of the present invention, which may be applied with any of the aspects of the present invention described herein, the exemplary apparatus/method receives measured satellite carrier phase data as received by the third base station, and further obtains a set of cycle ambiguities related to a set of satellite phase measurements associated with the baseline between the first and third base stations (B1-B3). The exemplary apparatus/method generates a sixth set of residuals of differential navigation equations associated with a set of measured carrier phase data related to the third base line between the rover and the third base station (R-B3), and resolves the cycle ambiguities in the carrier phase data from the sixth residual, the set of cycle ambiguities associated with the baseline between the first and third base stations, and one or more of the first, second, third, fourth, and fifth residuals. The resolved cycle ambiguities may take the form of floating ambiguities, fixed-integer ambiguities, and/or integer ambiguities. The exemplary apparatus/method estimates the rover's location further with the sixth set of residuals and the resolved cycle ambiguities associated with the third base line (R-B3).

[0012] In a seventh aspect of the present invention, which may be applied with any of the aspects of the present invention described herein, the exemplary apparatus/method obtains a first set of first ionosphere delay differentials associated with the satellite signals received along the base line formed by the first and second base stations, and generates ionosphere delay corrections to one or more of the above-described first through sixth residuals from the first set of first ionosphere delay differentials, the locations of the first and second base stations, and an estimated location of the rover station. As a further option, the exemplary apparatus/method obtains set of second ionosphere delay differentials associated with the satellite signals received along the base line formed by the first and third base stations (or the base line formed by the second and third base stations), and generates the ionosphere delay corrections to the one or more of the above-described first through sixth residuals further from the second set of ionosphere delay differentials and the location of the third base station.

[0013] In an eight aspect of the present invention, which may be applied to the seventh and further aspects of the present invention described herein, the exemplary apparatus/method forms one or more of the first through sixth residuals to account for second order effects in the ionosphere delay corrections applied to the baselines associated with the rover, and further generates an estimate of the second order effects.

[0014] In a ninth aspect of the present invention, which may be applied with any of the aspects of the present invention described herein, the exemplary apparatus/method generates the time offset representative of the time difference between the clocks of the first and second base stations, and optionally the time offset representative of the time difference between the clocks of the first and third base stations. As a further option, the apparatus/method generates the time offset representative of the time difference between the clocks of the second and third base stations, and performs a consistency check of the time offsets.

[0015] In an tenth aspect of the present invention, which may be applied with several of the aspects of the present invention described herein, the exemplary apparatus/method generates the resolved cycle ambiguities associated with the baseline between the first and second base stations (B1-B2), and optionally generates the resolved cycle ambiguities associated with the baseline between the first and third base stations (B1-B3). As a further option, the apparatus/method generates the resolved cycle ambiguities associated with the baseline between the second and third base stations (B2-B3), and performs a consistency check of the cycle ambiguities associated with the baselines between the base stations.

[0016] In an eleventh aspect of the present invention, which may be applied with any aspects of the present invention described herein which account for ionosphere delays, the exemplary apparatus/method generates the first set of first ionosphere delay differentials, and optionally the second set of first ionosphere delay differentials. As a further option, the apparatus/method generates a third set of ionosphere delay differentials associated with the satellite signals received along the base line formed by the second and third base stations, and generates therefrom three sets of ionosphere delay differentials which are self consistent.

[0017] Accordingly, it is an object of the present invention to increase the accuracy of estimating the rover's position using information from one or more baselines associated with two or more base stations.

[0018] It is a further object of the present invention to enable the spacing between the rover and the base stations to increase while maintaining or improving the accuracy of the estimation of the rover's position.

[0019] These and other objects of the present invention will become apparent to those skilled in the art from the following detailed description of the invention, the accompanying drawings, and the appended claims.

[0020] FIG. 1 is a perspective view of a rover station (R) and three base stations (B1, B2, B3) in an exemplary network according to exemplary embodiments of the present invention.

[0021] FIG. 2 is a top-plan schematic drawing the rover station (R) and three base stations (B1, B2, B3) in the exemplary network shown in FIG. 1 according to the present invention.

[0022] FIG. 3 is a perspective view of the ionosphere delay differentials to selected exemplary methods according to the present invention.

[0023] FIG. 4 is a top-plan schematic view of a road application where a reduced level of interpolation of ionosphere delays may be used according to the present invention.

[0024] FIG. 5 is a schematic diagram of an exemplary rover station according to the present invention.

[0025] FIG. 6 is a general flow diagram of embodiments of the present invention.

[0026] FIG. 7 is a schematic diagram of an exemplary computer program product according to the present invention.

[0027] FIG. 1 is a perspective view of a rover station (R) and three base stations (B1, B2, B3) in an exemplary network according exemplary embodiments of the present invention, and FIG. 2 is a top-plan schematic drawing thereof. The present invention pertains to estimating the position of the rover station with information provided from two or three of the base stations and with satellite measurements made by the rover station. Referring to FIG. 1, each station has a receiver that receives the satellite positioning signals with a satellite antenna (shown as a substantially flat disk). A plurality of satellites S1-S4 are depicted in FIG. 1, with the range between each satellite and each antenna being depicted by a respective dashed line. In the example shown and FIG. 1, the Rover station is operated by a human user, and has a positioning pole for positioning the Rover's satellite antenna over a location whose coordinates are to be determined. The Rover's satellite antenna is coupled to the receiver's processor, which may be disposed in the backpack and carried by the user. The user interacts with the receiver's processor through a keypad/display. The receiver also has a radio modem (more formally a demodulator) that can receive data from the base stations through a conventional RF antenna and relay it to the processor. In other embodiments, the Rover's satellite antenna may be mounted to a vehicle, and the receiver may operate independently and automatically without the need of a human operator. The present invention is applicable to these and other physical embodiments.

[0028] The satellite signals comprise carrier signals which are modulated by pseudo-random binary codes, which are then used to measure the delay relative to local reference clock or oscillator. These measurements enable one to determine the so-called pseudo-ranges between the receiver and the satellites. The pseudo-ranges are different from true ranges (distances) between the receiver and the satellites due to variations in the time scales of the satellites and receiver and various noise sources. To produce these time scales, each satellite has its own on-board atomic clock, and the receiver has its own on-board clock, which usually comprises a quartz crystal. If the number of satellites is large enough (four or more), then the measured pseudo-ranges can be processed to determine the user location (e.g., X, Y, and Z coordinates) and to reconcile the variations in the time scales. Finding the user location by this process is often referred to as solving a navigational problem or task.

[0029] More specifically, the GPS system employs a constellation of satellites in orbit around the earth at an altitude of approximately 26,000 km. Each GPS satellite transmits microwave radio signals in two frequency bands located around 1575.42 MHz and 1227.6 MHz, referred to as L1 band and L2 band, respectively. The GPS L1-band signal is modulated by a coarse/acquisition code (C/A) and a precision ranging code (P-code). The L2-band signal is binary modulated by the P-code. The GPS C/A code is a pseudo-random (PR) Gold code that is specific to each satellite and is used to identify the source of a received signal. The P-code is a pseudo-random code signals and is also specific to each satellite, having a symbol rate which is ten time more than C/A, which reduces the granularity by a factor of ten. The GPS satellite transmission standards are set in detail by the ICD-GPS-200, Revision C, ARINC Research Corporation, 10 Oct., 1993.

[0030] The satellites of the GLONASS system transmit signals in two frequency bands around 1602 MHz, and around 1246 MHZ, referred to also as L1 band and L2 band, respectively. The GLONASS L1-band signal is modulated by a pseudo-random C/A code and a pseudo-random P-code. The L2-band signal is modulated by the P-code. Unlike GPS, in which all of the satellites transmit on the same nominal frequency, the GLONASS satellites each transmit at a unique frequency in order to differentiate between the satellites. The GLONASS L1-band carrier frequency is equal to 1602 MHzl*0.5625 MHz, where l is an integer index ranging between 0 and 24 that identifies the satellites. The GLONASS L2-band carrier frequency is equal to 1246 MHzl*0.4375 MHz. Details of the GLONASS signals maybe found in GLONASS ICD, Version 4.0, 1998.

[0031] The distance between a receiver and a satellite (called the receiver-to-satellite range) is determined by measuring the time that it takes for the signal to pass from the satellite to the receiver, provided that the position of the satellite is known. The satellites and receivers have internal clocks that are synchronized to a single GPS time. For each satellite signal being tracked, the receiver generates a local version of the satellite's expected PR-code, and then retards that version in time until the local version correlates (i.e., matches) with the received satellite signal. Thereafter, the satellite signal is tracked by advancing or retarding the local version of the PR-code. The carrier phase of the satellite can also be tracked, which is usually done by tracking the Doppler shift of the satellite signal. The positions of the satellites are, except for minor variations, highly predictable as a function of time, and the receiver generally carriers a model of the satellite's position as a function of GPS time. In theory, by determining the ranges to three different satellites, the receiver can perform a three-dimensional triangulation to find its position. But because limitations in the accuracy of the receiver's clock, the internally generated time is offset somewhat from true GPS time. Thus, the ranges to at least four different satellites are simultaneously measured in order to be able to solve for four unknowns, namely the three coordinates of the position of the receiver location (e.g., x, y, and z) and an offset of the receiver clock time from the GPS time. The location is usually performed with respect the defined Cartesian coordinates frame.

[0032] In theory, a GPS receiver can track both the C/A code and the P-code of a satellite. The C/A code generally provides accuracy to within 20-50 meters, and the P-code generally provides greater accuracy to within 10 meters because of its higher repetition rate (less granularity) relative to the C/A code. However, knowledge of the P-code is restricted to military users and not available to civilian users. Nonetheless, some approaches for tracking P-codes have been developed.

[0033] The desire to guarantee the solution of navigational tasks with accuracy better than 10 meters, and the desire to raise the stability and reliability of measurements, have led to the development of the mode of differential navigation ranging, also called differential navigation (DN). In the DN mode, the task of finding the user position is performed relative to a Base station (Base), the coordinates of which are known with the high accuracy and precision. The Base station has a navigation receiver that receives the signals of the satellites and processes them to generate measurements. The results of these measurements enable one to calculate corrections, which are then transmitted to a roving GPS receiver, which the user has set up. We call this GPS receiver the Rover station, or Rover receiver. By using these corrections, the roving GPS receiver gains the ability to compensate for the major part of the strongly correlated errors in the measured pseudo-ranges, and to substantially improve the accuracy of the estimate of its position.

[0034] Usually, the Base station is immobile during measurements. The rover station may be either immobile or mobile. Depending on the navigational tasks to be solved, different modes of operation may be used in the DN mode. They differ in the way in which the measurement results are transmitted from the Base to the Rover. In the post-processing (PP) mode, these results are transmitted as digital recordings and go to the user after all the measurements have been finished. In the PP mode, the user reconstructs his or her location for definite moments in the past.

[0035] Another mode is the Real-Time Processing (RTP) mode, and it provides for the positioning of the Rover receiver just during the measurements. The RTP mode uses a communication link (such as the radio communication links shown in FIG. 1), through which all the necessary information is transmitted from the Base to the Rover receiver in digital form.

[0036] Further improvement of accuracy of differential navigation may be reached by supplementing the measurements of the pseudoranges with the measurements of the phases of the satellite carrier signals. If one measures the carrier phase of the signal received from a satellite in the Base receiver and compares it with the carrier phase of the same satellite measured in the Rover receiver, one can obtain measurement accuracy to within several percent of the carrier's wavelength, i.e., to within several centimeters.

[0037] The practical implementation of those advantages, which might be guaranteed by the measurement of the carrier phases, runs into the problem of there being ambiguities in the phase measurements.

[0038] The ambiguities are caused by two factors. First, the difference of distances D from any satellite to the Base and Rover is usually much greater than the carrier's wavelength . Therefore, the difference in the phase delays of a carrier signal D/ received by the Base and Rover receivers may substantially exceed one cycle. Second, it is not possible to measure the integer number of cycles in from the incoming satellite signals; one can only measure the fractional part of . Therefore, it is necessary to determine the integer part of , which is called the ambiguity. More precisely, we need to determine the set of all such integer parts for all the satellites being tracked, one integer part for each satellite. One has to determine this set along with other unknown values, which include the Rover's coordinates and the variations in the time scales.

[0039] In a general way, the task of generating highly-accurate navigation measurements may be formulated by defining a set of unknowns and system of relationships between the unknowns and measured navigation parameters. The vector of unknowns, denoted herein as n, include three Rover coordinates (usually along Cartesian axes X, Y, Z) in a given coordinate system (sometimes time derivatives of coordinates are added too); the variations of the time scales which is caused by the phase drift of the local main reference oscillator; and n integer unknown values associated with the ambiguities of the phase measurements of the carrier frequencies. The value of n is determined by the number of different carrier signals being processed, and accordingly coincides with the number of satellite channels actively functioning in the receiver. At least one satellite channel is used for each satellite whose broadcast signals are being received and processed by the receiver. Some satellites broadcast more than one code-modulated carrier signal, such as a GPS satellite, which broadcasts a carrier in the L₁ frequency band and a carrier in the L₂ frequency band. If the receiver processes the carrier signals in both of the L₁ and L₂ bands, the number of satellite channels (n) increases correspondingly.

[0040] Two sets of navigation parameters are measured by the Base and Rover receivers, respectively, and are used to determine the set of unknowns, which is sometimes called the state vector. Each set of navigation parameters includes the pseudo-range of each satellite to the receiver, and the full (complete) phase of each satellite carrier signal. Each pseudo-range is obtained by measuring the time delay of a code modulation signal of the corresponding satellite (C/A code or P-code). The code modulation signal is tracked by a delay-lock loop (DLL) circuit in each satellite-tracking channel. The full phase of a satellite's carrier signal is tracked by a phase-lock-loop (PLL) in the corresponding satellite tracking channel. (The DLL and PLL are, for example, provided by the demodulator 120 of an exemplary rover shown in FIG. 5). An observation vector is generated as the collection of the measured navigation parameters for specific (definite) moments of time.

[0041] The relationship between the state vector and the observation vector is defined by a well-known system of navigation equations. Given an observation vector, the system of equations may be solved to find the state vector if the number of equations equals or exceeds the number of unknowns in the state vector. In the latter case, conventional statistical methods are used to solve the system: the least squares method, the method of dynamic Kalman filtering, the method of least squares and various modifications of these methods.

[0042] One method may comprise the following steps. The measured values of the pseudo-ranges and full phases at specific (definite) moments of time, along with an indication of the satellites to which these measurements belong to and the time moments of the measurements, are transmitted from the Base to the Rover (such as through the communication link or as recordings). Corresponding values are measured in the Rover receiver. The processing includes the determination of the single differences of the pseudo-ranges and full phases between the Base and Rover measurements for each satellite. The strongly correlated errors are largely compensated (i.e., substantially cancelled) in the single differences. Then, the residuals of the single differences are calculated by subtraction of calculated values from the measured results. The processing of residuals allows one to linearize the initial system of navigation equations (sometimes several subsequent iterations are necessary for that), which makes possible the use of the well-developed body of mathematics for solving systems of linear equations. The subsequent iterative solution of the linearized system of navigation equations is equivalent to the minimization of the sum of squared non-linear residuals by the Gauss-Newton minimization method. The components of the state vector, with the n ambiguities included, are found as a result of the solution. But the calculated values of the ambiguities are not necessarily integer numbers, and are often floating point numbers. Because of this, they are called float ambiguities, or floating ambiguities, at this stage of the solution. To find true values of the integer ambiguities one uses the procedure of rounding off the float ambiguity vector to the nearest set of integers. This process is called the ambiguity resolution. After the ambiguity resolution has been done, is it possible to generate more accurate values of residuals and then, by solving the system of equation again, to find the coordinate values for the baseline connecting the Base and Rover, and consequently to more accurately estimate the coordinates of Rover and the correction to its clock drift.

[0043] There are fractional parts associated with each ambiguity due to the time offset in the receiver's carrier generator, and as a practical matter, it is more convenient to include this fractional offset into the floating ambiguities and integer ambiguities. In this case, one additionally solves for the unknown fractional offset of the receiver's carrier generator, and then resolves the integer ambiguities in the form of integers plus a common fractional part that is related to the time offset of the receiver's carrier generator. We call these resolved ambiguities fixed-integer ambiguities because, although these ambiguities have fractional parts, the difference between any two fixed-integer ambiguities for two corresponding satellites measured by the same pair of receivers is an integer.

[0044] Even with all of the-above described processing to account for clock offsets and carrier ambiguities, there are additional factors which affect the accuracy of measurements made with GPS and/or GLONASS signals. As one factor, the trajectory of each satellite (or its initial data called the ephemeris), is elliptical and is affected by natural causes such as solar winds. The accuracy of any measurement is dependent upon knowledge of the position of the satellites at certain time. An estimate of the ephemeris is calculated on earth for each of the satellites and is periodically uploaded to the satellite. The position information of a satellite is encoded onto a low frequency (50 Hz) signal which is modulated on to one of the carrier signals, and transmitted to the GPS receiver on earth.

[0045] Two additional and important factors that affect the accuracy of measurements made with GPS and/or GLONASS signals are the effects of the troposphere and ionosphere on propagation of signals from the satellites to the receivers. The troposphere is the lower part of the atmosphere and variations in the temperature, pressure, and humidity lead to spatial variations in the signal propagation. The ionosphere is at the upper part of the atmosphere and has a slice of ionized gas at the altitude of around 300 km. The density of ionized particles is sufficiently high to affect the propagation of electromagnetic signals, and has a spatial variation and a time variation. The ionosphere effect becomes even more important during years of high solar activity. These variations in the troposphere and ionosphere introduce variable delays in the propagation the satellite signals to the receivers. If the base and rover are widely separated, these delays will be significantly different for the base and rover stations, and will introduce error into the estimation of the rover's position.

[0046] With that general overview, we now describe the invention in greater detail. For the sake of simplifying the presentation, and without loss of generality, we assume N satellites identified by the index s, s1, . . . , N, and we assume that each receiver can track the satellite's L1 band signal and the satellites L2 band signal. The receivers extract timing information from the satellite signals, and report this information a predefined increments k of time, which we call epochs k. The time between epochs can be selected by the user, and generally ranges between 0.1 seconds to 2 seconds, with 1 second being typical. The clocks of all of the receivers are typically accurate to within several milliseconds of the true GPS time, and for practical purposes the receivers can determined the number k of the current epoch from their clocks. The following timing information can be extracted from each satellite s at each epoch k by each receiver r:

[0047] 1. The low-frequency (50 Hz) information signal which provides information on the orbit (position) of the satellites which enables the receiver to determine the satellite's position.

[0048] 2. the pseudo-range _r,k^L1,s derived from the L1-band C/A code (or optionally the L1-band P-code for military users), and having the units of length (e.g., meters m);

[0049] 3. the pseudo-range _r,k^L2,s derived from the L2-band C/A code (or optionally the L2-band P-code for military users), and having the units of length (e.g., meters m);

[0050] 4. the L1-band carrier phase _r,k^L1,s derived from the L1-band carrier, and typically having the dimensionless units of cycles (but is sometimes expressed in terms of meters by multiplying the cycles by the L1-band wavelength ^L1,s; and

[0051] 5. the L2-band carrier phase _r,k^L2,s derived from the L2-band carrier, and typically having the dimensionless units of cycles (but is sometimes expressed in terms of meters by multiplying the cycles by the L2-band wavelength ^L2,s.

[0052] Each of _r,k^L1,s, _r,k^L2,s, _r,k^L1,s, and _r,k^L2,s are known as observed quantities, or observables, because the are measured by the receiver from the satellite signals. The observables are related to the speed of light c, the wavelengths ^L1,s and ^L2,s of the L1-band and L2-band carriers of satellite s, the frequencies f^L1,Sc/^L1,s and f^L2,sc/^L2,s of the L1-band and L2-band carriers of satellite s, various variables that are to be determined, and various noise sources, in the following manner:

_r,k^L1,sR_r,k^sc(_k^s_r,k)T_r,k^sI_r,k^sn_r,k^L1,s1A $\begin{array}{cc}{}_{r,k}^{\mathrm{L2},s}=R_{r,k}^s+c\left({}_k^s-{}_{r,k}\right)+T_{r,k}^s+{\left(\frac{{}^{\mathrm{L2},s}}{{}^{\mathrm{L1},s}}\right)}^2{I}_{r,k}^s+n_{r,k}^{\mathrm{L2},s}& \text{[1B]}\\ {}_{r,k}^{\mathrm{L1},s}=\frac{1}{{}^{\mathrm{L1},s}}{R}_{r,k}^s+f^{\mathrm{L1},s}\left({}_k^s-{}_{r,k}\right)+\frac{1}{{}^{\mathrm{L1},s}}{T}_{r,k}^s-\frac{1}{{}^{\mathrm{L1},s}}{I}_{r,k}^s+N_r^{\mathrm{L1},s}+{}_r^{\mathrm{L1}}+v_{r,k}^{\mathrm{L1},s}& \text{[1C]}\\ {}_{r,k}^{\mathrm{L2},s}=\frac{1}{{}^{\mathrm{L2},s}}{R}_{r,k}^s+f^{\mathrm{L2},s}\left({}_k^s-{}_{r,k}\right)+\frac{1}{{}^{\mathrm{L2},s}}{T}_{r,k}^s-\frac{{}^{\mathrm{L2},s}}{{\left({}^{\mathrm{L1},s}\right)}^2}{I}_{r,k}^s+N_r^{\mathrm{L2},s}+{}_r^{\mathrm{L2}}+v_{r,k}^{\mathrm{L2},s}& \text{[1D]}\end{array}$

[0053] where variables (which are to be determined) are:

[0054] R_r,k^s is the true range between the satellite s and receiver r at the epoch k, having dimension of distance (usually expressed as meters m);

[0055] _k^s is the time offset from true GPS time of the clock of satellite s at epoch k;

[0056] _r,k is the time offset from true GPS time of the clock of receiver r at epoch k;

[0057] T_r,k^s is the delay due to tropospheric effects between satellite s and receiver r at epoch k;

[0058] I_r,k^s is the delay due to ionospheric effects between satellite s and receiver r at epoch k (the full amount of I_r,k^s affects the pseudo-range in the L1-band);

[0059] N_r^L1,s is the integer phase ambiguity in the L1-band observable _r,k^L1,s;

[0060] N_r^L2,s is the integer phase ambiguity in the L2-band observable _r,k^L2,s;

[0061] _r^L1 is the initial phase offset of the L1-band phase tracker of receiver r; and

[0062] _r^L2 is the initial phase offset of the L2-band phase tracker of receiver r,

[0063] and where various noise sources are:

[0064] n_r,k^L1,s is the total effect of the noise sources on the pseudo-range measurement at epoch k in receiver r when tracking the L1-band C/A-code signal (or P-code signal) of satellite s;

[0065] n_r,k^L2,s is the total effect of the noise sources on the pseudo-range measurement at epoch k in receiver r when tracking the L2-band C/A-code signal (or P-code signal) of satellite s;

[0066] v_r,k^L1,s is the total effect of the noise sources on the carrier-phase measurement at epoch k in receiver r when tracking the L1-band carrier of satellite s; and

[0067] v_r,k^L1,s is the total effect of the noise sources on the carrier-phase measurement at epoch k in receiver r when tracking the L2-band carrier of satellite s.

[0068] Forms 1A-1D are known as the satellite navigation equations.

[0069] In embodiments of the present invention, we will be forming differences between the observables of pairs of stations, and differences between instances of the satellite navigation equations. The vector, or straight line, between a pair of two such stations is called a baseline. The baselines between all of the stations shown in FIG. 1 are best seen in the top plan view of FIG. 2. We generally note the difference of an observable measured by two stations with the prefix notation _q,r, where each of q and r represents one of the stations B1, B2, B3, and R. For example,

[0070] Forms 2A-2D are often referred to as the between-station (or between-receiver) single differences of the signals of satellite s between stations (receivers) r and q. The between station operator _q,r can be applied to the variables and noise sources as well. For example, _q,rR_k^sR_r,k^sR_q,k^s, and _q,rn_k^L1,sn_r,k^L1,sn_q,k^L1,s. With this background, we can relate the between-station single differences of the observables to the differences in the variables as follows:

_q,rk^L1,s_q,rR_k^sc_q,rkq,rT_k^s_q,rI_k^s_q,rn_k^L1,s3A $\begin{array}{cc}{}_{q,r}{}_k^{\mathrm{L2},s}={}_{q,r}{R}_k^s-c{}_{q,r}{}_k+{}_{q,r}{T}_k^s+{\left(\frac{{}^{\mathrm{L2},s}}{{}^{\mathrm{L1},s}}\right)}^2{}_{q,r}{I}_k^s+{}_{q,r}{n}_k^{\mathrm{L2},s}& \text{[3B]}\\ {}_{q,r}{}_k^{\mathrm{L1},s}=\frac{1}{{}^{\mathrm{L1},s}}{}_{q,r}{R}_k^s-f^{\mathrm{L1},s}{}_{q,r}{}_k+\frac{1}{{}^{\mathrm{L1},s}}{}_{q,r}{T}_k^s-\frac{1}{{}^{\mathrm{L1},s}}{}_{q,r}{I}_k^s+{}_{q,r}{N}^{\mathrm{L1},s}+{}_{q,r}{}^{\mathrm{L1}}+{}_{q,r}{v}_k^{\mathrm{L1},s}& \text{[3C]}\\ {}_{q,r}{}_k^{\mathrm{L2},s}=\frac{1}{{}^{\mathrm{L2},s}}{}_{q,r}{R}_k^s-f^{\mathrm{L2},s}{}_{q,r}{}_k+\frac{1}{{}^{\mathrm{L2},s}}{}_{q,r}{T}_k^s-\frac{{}^{\mathrm{L2},s}}{{\left({}^{\mathrm{L1},s}\right)}^2}{}_{q,r}{I}_k^s+{}_{q,r}{N}^{\mathrm{L2},s}+{}_{q,r}{}^{\mathrm{L2}}+{}_{q,r}{v}_k^{\mathrm{L2},s}& \text{[3D]}\end{array}$

[0071] The difference form of 3A is generated by forming two instance of form 1A for receivers r and q, and then subtracting the two instances (the instance for receiver q is subtract from the instance for receiver r). Difference forms 3B, 3C, and 3D, are formed in a similar manner from corresponding instances of forms 1B, 1C, and 1D, respectively. Forms 3A-3D are the single-differences of the navigation equations. The benefit of forming the between-station single differences is that the error term representing the time offset of the satellite clock, _k^s, is cancelled out in the differences. We emphasize that the forms 3A-3D can be applied to each satellite that can be observed by receivers q and r. Higher-order differences of the navigation equations, such as the double-differences of the navigation equations, can be formed and are known to the art. For example, a common double-difference equation is the difference between two singles-difference equations associated with a common baseline, but with each single-difference equation being based on a different satellite. The present invention may also use these higher-order differences of the navigation equations, although the single differences of the navigation equations are currently preferred. Each of these difference forms is generically referred to as a differential navigation equation.

[0072] Between Base Station Processing, Part I.

[0073] Because the locations of the base stations are known, and because the locations of the satellites as a function of epoch time (k) are also known, the between station differences of each true range _q,rR_k^sR_r,k^sR_q,k^s associated with each satellite s can be generated from existing information. If vectors X_r, X_q, and X^s represent the locations of receiver r, receiver q, and satellite s, respectively, in a coordinate system (e.g., Cartesian coordinate system), and the operator represent the distance operator in that coordinate system, then _q,rR_k^sX^sX_rX^sX_q. This reduces by one the number of unknowns in forms 3A-3D. In addition, the Goad-Goodman model may be used to model the difference in troposphere effects (with an error less than a few percent), and thus the difference _q,rT_k^s may be estimated based on the positions of the receivers and the satellite (more specifically, the angles between the receivers and satellites). The noise sources _q,rn^L1,s, _q,rn^L2,s, _q,rv^L1,s, and _q,rv^L2,s can never be known, but they are generally zero-mean and their effect can be reduced by averaging. Thus, the number of solvable unknowns in forms 3A-3D may be reduced to the following six: (1) _q,rk, (2) _q,rI_k^s, (3) _q,rN^L1,s, (4) _q,r^L1, (5) _q,rN^L2,s (6) _q,r^L2. The first unknown varies with time and is common to all of the satellites being tracked by receivers q and r. The second unknown varies with time and is specific to the satellite s being tracked by the pair of receivers q and r. The third and fifth unknowns are each specific to the satellite s being tracked by the pair of receivers q and r, and each does not normally vary with time unless a cycle slip occurs in the receivers phase-lock loop. The forth and sixth unknowns are specific to the pair of receivers q and r, and do not normally vary with time.

[0074] As a practical matter, we note that the third and fourth unknowns are related to the L1-carrier phase measurement of satellite s and are time-independent (unless a cycle slip occurs). In theory, these unknowns can be represented as a combined unknown _q,r{circumflex over (N)}^L1,s_q,rN^L1,s_q,r^L1, which we call the between-station fixed-integer ambiguity for the L1 carrier of satellite s, as measured at stations q and r. We call these resolved ambiguities fixed-integer ambiguities because, although these ambiguities have fractional parts, the difference between any two fixed-integer ambiguities for two corresponding satellites measured by the same pair of receivers is an integer. In a similar manner, the fifth and sixth unknowns are related to the L2-carrier phase measurement of satellite s and are time-independent (unless a cycle slip occurs); they can be represented as a combined unknown _q,r{circumflex over (N)}^L2,s_q,rN^L2,s_q,r^L2, which we call the between-station fixed-integer ambiguity for the L2 carrier of satellite s, as measured at stations q and r.

[0075] As a further practical matter, one first generates floating point versions of the fixed-integer ambiguities, and then applies a discretization process to the floating point version to generate the fixed-integer ambiguities. We call these versions the floating-point ambiguities, and denote them as _q,r{overscore (N)}^L1,s_q,rN^L1,s_q,r^L1 and _q,r{overscore (N)}^L2,s_q,rN^L2,s_q,r^L2.

[0076] Thus, the solvable unknowns for each instance of form 3A-3D can be reduced to four: (1) _q,rk, (2) _q,rI_k^s, (3) _q,r{circumflex over (N)}^L1,s (_q,r{overscore (N)}^L1,s), and (4) _q,r{circumflex over (N)}^L2,s (_q,r{overscore (N)}^L2,s). The first unknown is common to all the satellites being tracked by the pair of receivers q and r, whereas the last three are specific to the satellite s. The first two unknowns vary with time, whereas the last two are time invariant (unless a cycle slip occurs). For the base stations, the known values of forms 3A-3D can be represented as:

_q,rb_k^L2,s_q,rk^L2,s_q,rR_k^s_q,rT_k^s4B $\begin{array}{cc}{}_{q,r}{p}_k^{\mathrm{L1},s}={}_{q,r}{}_k^{\mathrm{L1},s}-\frac{1}{{}^{\mathrm{L1},s}}{}_{q,r}{R}_k^s-\frac{1}{{}^{\mathrm{L1},s}}{}_{q,r}{T}_k^s& \text{[4C]}\\ {}_{q,r}{p}_k^{\mathrm{L2},s}={}_{q,r}{}_k^{\mathrm{L2},s}-\frac{1}{{}^{\mathrm{L2},s}}{}_{q,r}{R}_k^s-\frac{1}{{}^{\mathrm{L2},s}}{}_{q,r}{T}_k^s& \text{[4D]}\end{array}$

[0077] Forms 4A-4D comprise the known terms of the single-differences of the navigation equations, and we call them the residuals of the single differences, or more generally the residuals of a set of differential navigation equations. The magnitude of each residual is generally less than the magnitudes of one or more of the known terms that forms the residual.

[0078] Then, the unknowns can be estimated from the residuals using the following forms by any of number of solution methods known to the art:

66 _q,rb_k^L1,sc_q,rkq,rI_k^s_q,rn_k^L1,s5A $\begin{array}{cc}{}_{q,r}{b}_k^{\mathrm{L2},s}=-c{}_{q,r}{}_k+{\left(\frac{{}^{\mathrm{L2},s}}{{}^{\mathrm{L1},s}}\right)}^2{}_{q,r}{I}_k^s+{}_{q,r}{n}_k^{\mathrm{L2},s}& \text{[5B]}\\ {}_{q,r}{p}_k^{\mathrm{L1},s}=-f^{\mathrm{L1},s}{}_{q,r}{}_k-\frac{1}{{}^{\mathrm{L1},s}}{}_{q,r}{I}_k^s+{}_{q,r}{\hat{N}}^{\mathrm{L1},s}+{}_{q,r}{v}_k^{\mathrm{L1},s}& \text{[5C]}\\ {}_{q,r}{p}_k^{\mathrm{L2},s}=-f^{\mathrm{L2},s}{}_{q,r}{}_k-\frac{{}^{\mathrm{L2},s}}{{\left({}^{\mathrm{L1},s}\right)}^2}{}_{q,r}{I}_k^s+{}_{q,r}{\hat{N}}^{\mathrm{L2},s}+{}_{q,r}{v}_k^{\mathrm{L2},s}& \text{[5D]}\end{array}$

[0079] All of the residuals, forms, known values, and unknown variables are real-valued. The estimation process usually includes instances of forms 5A-5D for several satellites, and for several epochs of time.

[0080] As one aspect of the present invention, one or more of the unknown variables for several satellites are estimated for at least two base stations, and preferably for at least three base stations by a process that is performed at the rover. Using the numbers 1, 2, and 3 to represent the base stations B1, B2, and B3, respectively, in the difference operator _q,r, the estimated values of the between-base station unknowns for each satellite s at each epoch k are:

[0081] Base Stations 1 and 2

[0082] Base Stations 1 and 3

[0083] Base Stations 2 and 3

[0084] An exemplary estimation process is outlined in a section below entitled Between Base Station Processing, Part II. As one option, the estimation process may break down the fixed integer ambiguities into their components as follows:

[0085] Base Stations 1 and 2

[0086] Base Stations 1 and 3

[0087] Base Stations 2 and 3

[0088] In some applications of the present invention where lower estimate accuracy is acceptable, floating-point ambiguities may be used in place of the fixed-integer ambiguities, and the following set of between-base station unknowns may be generated:

[0089] Base Stations 1 and 2

[0090] Base Stations 1 and 3

[0091] Base Stations 2 and 3

[0092] Between Rover and Base Station Processing.

[0093] The processing of differences between the rover station and any one of the base stations is more complex because we do not know the position of the rover beforehand. Thus, the term _q,rR_k^sR_r,k^sR_q,k^s is not exactly known. However, given an initial estimated value of the rover's location, which we will denote as {overscore (X)}₀ and which can be generated from the 50 Hz information signals by means known to the art, we can generate an initial range estimate between satellite s and the rover, which we denote as {tilde over (R)}_0,k^s, and we can generate a linearized approximation of the term as follows, where we use r0 to denote the rover and use q to generally indicate any one of the base stations: $\begin{array}{c}{}_{q,0}{R}_k^s=R_{0,k}^s-R_{q,k}^s\left({\tilde{R}}_{0,k}^s+a_k^{s}X\right)-R_{q,k}^s\\ \left({\tilde{R}}_{0,k}^s-R_{q,k}^s\right)+a_k^s{X}_k\\ {}_{q,0}{\tilde{R}}_{0,k}^s+a_k^s{X}_k\end{array}$

[0094] where X_k is the vector difference between the actual location X₀ of the Rover and the initial estimated location {overscore (X)}₀ (X₀{overscore (X)}₀X_k), where a_k^s is a row vector comprised of the partial derivatives of the distance form R_0,k^sX^sX₀ with respect to the rover location X₀, as evaluated at {overscore (X)}₀, and where _q,0{tilde over (R)}_0,k^s has been used as a shorthand notation for ({tilde over (R)}_0,k^sR_q,k^s). Formally, row vector a_k^s is expressed as: ${{a_k^s=\frac{}{X_0}\left\{R_{0,k}^s\right\}}_{{\stackrel{\_}{X}}_0}=\frac{}{X_0}\left\{X^s-X_0\right\}}_{{\stackrel{\_}{X}}_0}$

[0095] The values of a_k^s depend upon the satellite location, which changes with each time moment k, and upon the initial location {overscore (X)}₀ of the Rover, which may be periodically updated, particularly if the Rover is moving. The row vectors a_k^s for all of the satellites may be collected together as a single s-by-3 matrix, which we denote as A_k. Matrix A_k is commonly known as the Jacobian matrix, the geometry matrix, and the matrix of directional cosines. It is commonly computed in the art and described in many tutorial textbooks on global positioning. The reader unfamiliar with the GPS is directed to these tutorial texts for more detailed explanations of how matrix A_k (which is oftentimes referred to matrix H or matrix G by these texts) is generated.

[0096] With _q,0R_k^s_q,0{tilde over (R)}_0,k^sa_k^sX, we may write the following versions of forms 4 (residuals) for the baselines between the Rover (0) and any one of the base stations (q):

_q,0b_k^L2,s_q,0k^L2,s_q,0{tilde over (R)}_k^s_q,0T_k^s8B $\begin{array}{cc}{}_{q,0}{p}_k^{\mathrm{L1},s}={}_{q,0}{}_k^{\mathrm{L1},s}-\frac{1}{{}^{\mathrm{L1},s}}{}_{q,0}{\tilde{R}}_k^s-\frac{1}{{}^{\mathrm{L1},s}}{}_{q,0}{T}_k^s& \text{[8C]}\\ {}_{q,0}{p}_k^{\mathrm{L2},s}={}_{q,0}{}_k^{\mathrm{L2},s}-\frac{1}{{}^{\mathrm{L2},s}}{}_{q,0}{\tilde{R}}_k^s-\frac{1}{{}^{\mathrm{L2},s}}{}_{q,0}{T}_k^s& \text{[8D]}\end{array}$

[0097] and the following versions of forms 5

_q,0b_k^L1,sa_k^sX_kc_q,0kq,0I_k^s_q,0n_k^L1,s9A $\begin{array}{cc}{}_{q,0}{b}_k^{\mathrm{L2},s}=a_k^s{X}_k-c{}_{q,0}{}_k+{\left(\frac{{}^{\mathrm{L2},s}}{{}^{\mathrm{L1},s}}\right)}^2{}_{q,0}{I}_k^s+{}_{q,0}{n}_k^{\mathrm{L2},s}& \text{[9B]}\\ {}_{q,0}{p}_k^{\mathrm{L1},s}=\frac{1}{{}^{\mathrm{L1},s}}{a}_k^s{X}_k-f^{\mathrm{L1},s}{}_{q,0}{}_k-\frac{1}{{}^{\mathrm{L1},s}}{}_{q,0}{I}_k^s+{}_{q,0}{\hat{N}}^{\mathrm{L1},s}+{}_{q,0}{v}_k^{\mathrm{L1},s}& \text{[9C]}\\ {}_{q,0}{p}_k^{\mathrm{L2},s}=\frac{1}{{}^{\mathrm{L2},s}}{a}_k^s{X}_k-f^{\mathrm{L2},s}{}_{q,0}{}_k-\frac{{}^{\mathrm{L2},s}}{{\left({}^{\mathrm{L1},s}\right)}^2}{}_{q,0}{I}_k^s+{}_{q,0}{\hat{N}}^{\mathrm{L2},s}+{}_{q,0}{v}_k^{\mathrm{L2},s}& \text{[9D]}\end{array}$

[0098] The solvable unknowns for form 9 includes X_k, as well as _q,0k, _q,0I_k^s, _q,0{circumflex over (N)}^L1,s, and _q,0{circumflex over (N)}^L2,s. As described below in greater detail, estimation processes for these unknowns generally include instances of forms 9A-9D for several satellites, and for several epochs of time.

[0099] While the presentation of this section has been made with respect to a specific rover station (0), it may be appreciated that the discussion and above forms generally apply to any pair of receivers (q, r), where r substitutes for 0 in the above.

[0100] Development of the Present Invention

[0101] As part of making their invention, the inventors have recognized that the base station time-offset unknowns, when accurately estimated, should satisfy the following relationship:

[0102] An underlying meaning of this relationship can be seen by noting that _3,1k1,3k, and by rewriting the relationship in the equivalent form of _2,1k3,2k1,3k0. That is, the sum of the time offsets around a loop of base stations is equal to zero. The other unknowns, when accurately estimated, should satisfy the following similar relationships:

[0103] As a further part of making their invention, the inventors have recognized that the above relationships 10A-10J should be satisfied by any set of three receivers, including sets where one of the receivers is the receiver of the rover station. We illustrate this point by replacing station B3 (3) with the rover station, which we will identify by the subscript number 0:

[0104] Still further, the inventors have recognized that relationships 10 and 11 can be generalized to apply to any loop of four or more receivers. Nonetheless, we provide exemplary embodiments of the present invention using loops of three receivers in order to simplify the presentation of the invention.

[0105] According to general aspects of the present invention, the position of the Rover is estimated by forming a primary baseline between the Rover and one of the base stations, usually the closest base station, and then forming one or more secondary baselines from the Rover to other base stations. Then, one or more of the above relationships are applied to the configuration of stations to relate the measured data associated with the secondary baseline(s) with the measured data associated primary baseline using one or more of the unknowns associated with the baselines between base stations through forms 11. To simplify the presentation of the present invention, and without loss of generality, we will use base station B1 to form the primary baseline with the Rover, and we use base stations B2 and B3 to form the secondary baselines with the Rover. Varying amounts measured data from the secondary baselines may be related with the primary baseline. In general, the greater amount of information so related increases the accuracy of the estimated position of the rover, and/or enables a greater spacing between base stations.

[0106] First General Group of Embodiments

[0107] In a first exemplary embodiment, the measured pseudo-range data from the secondary baselines are related with the measured pseudo-range data of the primary baseline. At first, only the L1-band data is used (forms 4A, 5A, 8A, and 9A). We take the case where the distances between the stations are small enough that the difference terms associated with the ionosphere terms (_q,rI_k^s) are negligible. For the primary baseline, forms 8A and 9A reduce to:

[0108] For the secondary baseline to the second base station, we have:

[0109] For the secondary baseline to the third base station, we have:

[0110] The terms in forms 12A, 14A, and 16A are the residuals of differential navigation equations and are known. The unknowns are contained in right-hand sides of forms 13A, 15A, and 17A. There are a total of six solvable unknowns: _1,0k, _2,0k, _3,0k, and the three components of X_k. However, using form 11A, _2,0k may be related to _0,1k through the estimated base-station time offset _2,1k as follows: _{2,0k2,1k0,1k2,1k1,0k}. In a similar manner, _3,0k may be related to _0,1k through the estimated base-station time offset _3,1k as follows: _{3,0k3,1k0,1k3,1k1,0k}. These two forms reduce number of true solvable unknowns to four, and the following modified set of forms may be used to estimate the true solvable unknowns:

[0111] As an equivalent, one may view forms _2,0k2,1k1,0k and _3,0k3,1k1,0k as increasing the total number of forms by two, and may then estimate the six unknowns the expanded form set:

[0112] The rover's location (X₀{overscore (X)}₀X_k) may be estimated from each of the above form sets using any one of several know methods. Here, we demonstrate estimating the rover's location with a least squares fitting approach based on the first form set. After generating the residuals _1,0b_k^L1,s for N satellites (s1 to sN) for the baseline between the rover and the first base station (R-B1), we collect them together into a first vector:

[0113] We do the same for the other two baselines, but modify their residuals according to the between-base station time offsets c_2,1k and c_3,1k:

[0114] where the superscript ^* is used to indicate the modification. The noise terms can be similarly grouped as: _1,0n_k^L1, _2,0n_k^L1, and _3,0n_k^L1. The first form set of 13A, 15A* and 17A* may then be written in matrix form as: $\begin{array}{cc}\left[\begin{array}{c}{}_{1,0}{B}_k^{\mathrm{L1}}\\ {}_{2,0}{B}_k^{{}^{*}\mathrm{L1}}\\ {}_{3,0}{B}_k^{{}^{*}\mathrm{L1}}\end{array}\right]=\left[\begin{array}{cc}{A}_k& -\stackrel{}{1}\\ A_k& -\stackrel{}{1}\\ A_k& -\stackrel{}{1}\end{array}\right]\left[\begin{array}{c}{X}_k\\ c{}_{1,0}{}_k\end{array}\right]+\left[\begin{array}{c}{}_{1,0}{n}_k^{\mathrm{L1}}\\ {}_{2,0}{n}_k^{\mathrm{L1}}\\ {}_{3,0}{n}_k^{\mathrm{L1}}\end{array}\right]& \text{[18]}\end{array}$

[0115] The solvable unknowns $\left[\begin{array}{c}{X}_k\\ c{}_{1,0}{}_k\end{array}\right]$

[0116] may be estimated by a least squares process as follows: $\begin{array}{cc}\left[\begin{array}{c}{X}_k\\ c{}_{1,0}{}_k\end{array}\right]={\left[{\left(A_k^{*}\right)}^T{C}_{n,k}^{-1}{A}_k^{*}\right]}^{-1}{\left(A_k^{*}\right)}^T{C}_{n,k}^{-1}\left[\begin{array}{c}{}_{1,0}{B}_k^{\mathrm{L1}}\\ {}_{2,0}{B}_k^{{}^{*}\mathrm{L1}}\\ {}_{3,0}{B}_k^{{}^{*}\mathrm{L1}}\end{array}\right]& \text{[19]}\end{array}$

[0117] where: $A_k^{*}=\left[\begin{array}{cc}{A}_k& -\stackrel{}{1}\\ A_k& -\stackrel{}{1}\\ A_k& -\stackrel{}{1}\end{array}\right],$

[0118] and where C_n,k is a 3N-by-3N covariance matrix for the noise terms _1,0n_k^L1, _2,0n_k^L1, and _3,0n_k^L1. The generation of covariance matrix C_n,k is well known to the GPS art and described in the technical and patent literature (see for example U.S. Pat. No. 6,268,824, incorporated herein by reference, where the covariance matrices are denoted as matrices R), and a description thereof is not needed for one of ordinary skill in the GPS art to make and use the present invention. (In a subsequent section, we describe an estimation process that we prefer to use, and which may be applied to these embodiments.)

[0119] For the reader not familiar with the GPS art, we briefly note that matrix C_n,k generally comprises a diagonal matrix, with each diagonal element being related to the noise sources in the two receivers that define the baseline (in this case the rover and one of the base stations). A covariance factor is usually associated with each satellite signal received by each receiver, and this covariance factor is usually related to the signal-to-noise ratio of the signal (as received by the receiver) and the elevation angle of the satellite (the multipath error has a strong correlation with the elevation angle). Each diagonal entry of matrix C_n,k usually comprises an addition of the two covariance factors associated with the two receivers which contributed to the underlying noise quantity, e.g., base station B1 and the rover for noise quantity _1,0n_k^L1,s. For more details and information on the generation of the covariance matrices, the reader is direct to A. Leick, GPS Satellite Surveying, John Wiley & Sons, (1995). Matrix A_k^* is known as the observation matrix, and it relates the solvable unknowns to the residuals. The solvable unknowns may also be estimated by other processes, such as various Kalman filtering processes.

[0120] The above may be carried out by omitting the rows related to one of the secondary base stations B2 and B3, such as the last N rows associated with the third base station. The above may also be carried out by adding rows related to a fourth base station (and even more base stations).

[0121] The second form set comprising 13A, 15A, 17A, 11A*, and 11A** may be written as: $\begin{array}{cc}\left[\begin{array}{c}{}_{1,0}{B}_k^{\mathrm{L1}}\\ {}_{2,0}{B}_k^{\mathrm{L1}}\\ {}_{3,0}{B}_k^{\mathrm{L1}}\\ {}_{2,1}{}_k\\ {}_{3,1}{}_k\end{array}\right]=\left[\begin{array}{cccc}{A}_k& -\stackrel{}{1}& \stackrel{}{0}& \stackrel{}{0}\\ A_k& \stackrel{}{0}& -\stackrel{}{1}& \stackrel{}{0}\\ A_k& \stackrel{}{0}& \stackrel{}{0}& -\stackrel{}{1}\\ 0& -1& 1& 0\\ 0& -1& 0& 1\end{array}\right]\left[\begin{array}{c}{X}_k\\ c{}_{1,0}{}_k\\ c{}_{2,0}{}_k\\ c{}_{3,0}{}_k\end{array}\right]+\left[\begin{array}{c}{}_{1,0}{n}_k^{\mathrm{L1}}\\ {}_{2,0}{n}_k^{\mathrm{L1}}\\ {}_{3,0}{n}_k^{\mathrm{L1}}\\ 0\\ 0\end{array}\right]& \text{[20]}\end{array}$

[0122] and the solvable unknowns may similarly be estimated by various least squares processes, and Kalman filtering processes. The observation matrix in the above form is denoted as A_k^**.

[0123] Thus, the above exemplary embodiments provide methods of estimating the location of the rover station (R) with the use of a first base station (B1) a second base station (B2), and optionally a third base station (B3) or more base stations. In summary, the locations of the base stations were obtained, and one of the base stations (e.g., B1) was selected to form a primary baseline with the Rover station. We refer to this base station as the primary base station, and the other base stations as the secondary base stations. Additionally, for each of the secondary base stations, the time offset representative of the time difference between the clocks of the primary station and secondary base station is obtained. Also, measured satellite data as received by the rover, the primary base station, and the secondary base station(s) is obtained. From this, the set of residuals _1,0B_k^L1 associated with the primary baseline (R-B1) is generated, and is related to the measured satellite data received by the rover station and the first base station, the locations of the satellites, and the locations of the rover station and the first base station. Similarly, the set(s) of residuals _2,0B_k^L1, _3,0B_k^L1, etc. associated with the secondary baseline(s) are generated, each set of residuals related to the measured satellite data received by the rover station and the secondary base station, the locations of the satellites, and the locations of the rover station and the secondary base station. Thereafter, the rover's location is estimated from the above sets of residuals, the time offset between the clocks of the base stations, and typically an observation matrix.

[0124] The above embodiments can be equally implemented using the L2-band signals. The forms for the primary baseline (neglecting ionosphere terms) are:

[0125] For the secondary baseline to the second base station, we have:

[0126] For the secondary baseline to the third base station, we have:

[0127] Not only may these forms be used in place of forms 12A-17B, but they may be used in addition to forms 12A-17B. In the latter case, the first form set 18 expands to $\begin{array}{cc}\left[\begin{array}{c}{}_{1,0}{B}_k^{\mathrm{L1}}\\ {}_{2,0}{B}_k^{{}^{*}\mathrm{L1}}\\ {}_{3,0}{B}_k^{{}^{*}\mathrm{L1}}\\ {}_{1,0}{B}_k^{\mathrm{L2}}\\ {}_{2,0}{B}_k^{{}^{*}\mathrm{L2}}\\ {}_{3,0}{B}_k^{{}^{*}\mathrm{L2}}\end{array}\right]=\left[\begin{array}{cc}{A}_k& -\stackrel{}{1}\\ A_k& -\stackrel{}{1}\\ A_k& -\stackrel{}{1}\\ A_k& -\stackrel{}{1}\\ A_k& -\stackrel{}{1}\\ A_k& -\stackrel{}{1}\end{array}\right]\left[\begin{array}{c}{X}_k\\ c{}_{1,0}{}_k\end{array}\right]+\left[\begin{array}{c}{}_{1,0}{n}_k^{\mathrm{L1}}\\ {}_{2,0}{n}_k^{\mathrm{L1}}\\ {}_{3,0}{n}_k^{\mathrm{L1}}\\ {}_{1,0}{n}_k^{\mathrm{L2}}\\ {}_{2,0}{n}_k^{\mathrm{L2}}\\ {}_{3,0}{n}_k^{\mathrm{L2}}\end{array}\right]& \left[{18}^{+}\right]\end{array}$

[0128] where:

[0129] From form 18, the solvable unknowns may be estimated by a least squares process.

[0130] One might envision a different approach whereby the Rover's position is estimated by using each of the three baselines separately to generate three separate estimates of the Rover's location, and then averaging the three separate estimates to generate a final estimate. However, the present invention provides more accurate results than this possible approach because the present invention uses the additional information provided by forms 11 (specifically, 11A, 11A* and 11A**), which correlates the values of the unknowns associated with the baselines. The achievement of more accurate results by the present invention is also true for the embodiments described below.

[0131] Second General Group of Embodiments

[0132] The second group of embodiments builds upon the first group of embodiments by generating and utilizing the residuals derived from the phase measurements at the receivers. We take the case where the distances between the stations are small enough that the difference terms associated with the ionosphere terms (_q,rI_k^s) are negligible. In later embodiments, we will take these ionosphere terms into account. With this assumption, forms 8C, 9C, 8D, and 9D for the primary baseline become the following forms 12C, 13C, 12D, and 13D, respectively: $\begin{array}{cc}{}_{1,0}{p}_k^{\mathrm{L1},s}={}_{1,0}{}_k^{\mathrm{L1},s}-\frac{1}{{}^{\mathrm{L1},s}}{}_{1,0}{\tilde{R}}_k^s-\frac{1}{{}^{\mathrm{L1},s}}{}_{1,0}{T}_k^s& \text{[12C]}\\ {}_{1,0}{p}_k^{\mathrm{L1},s}=\frac{1}{{}^{\mathrm{L1},s}}{a}_k^s{X}_k-f^{\mathrm{L1},s}{}_{1,0}{}_k+{}_{1,0}{\hat{N}}^{\mathrm{L1},s}+{}_{1,0}{v}_k^{\mathrm{L1},s}& \text{[13C]}\\ {}_{1,0}{p}_k^{\mathrm{L2},s}={}_{1,0}{}_k^{\mathrm{L2},s}-\frac{1}{{}^{\mathrm{L2},s}}{}_{1,0}{\tilde{R}}_k^s-\frac{1}{{}^{\mathrm{L2},s}}{}_{1,0}{T}_k^s& \text{[12D]}\\ {}_{1,0}{p}_k^{\mathrm{L2},s}=\frac{1}{{}^{\mathrm{L2},s}}{a}_k^s{X}_k-f^{\mathrm{L2},s}{}_{1,0}{}_k+{}_{1,0}{\hat{N}}^{\mathrm{L2},s}+{}_{1,0}{v}_k^{\mathrm{L2},s}& \text{[13D]}\end{array}$

[0133] For the secondary baseline to the second base station, we have: $\begin{array}{cc}{}_{2,0}{p}_k^{\mathrm{L1},s}={}_{2,0}{}_k^{\mathrm{L1},s}-\frac{1}{{}^{\mathrm{L1},s}}{}_{2,0}{\tilde{R}}_k^s-\frac{1}{{}^{\mathrm{L1},s}}{}_{2,0}{T}_k^s& \text{[14C]}\\ {}_{2,0}{p}_k^{\mathrm{L1},s}=\frac{1}{{}^{\mathrm{L1},s}}{a}_k^s{X}_k-f^{\mathrm{L1},s}{}_{2,0}{}_k+{}_{2,0}{\hat{N}}^{\mathrm{L1},s}+{}_{2,0}{v}_k^{\mathrm{L1},s}& \text{[15C]}\\ {}_{2,0}{p}_k^{\mathrm{L2},s}={}_{2,0}{}_k^{\mathrm{L2},s}-\frac{1}{{}^{\mathrm{L2},s}}{}_{2,0}{\tilde{R}}_k^s-\frac{1}{{}^{\mathrm{L2},s}}{}_{2,0}{T}_k^s& \text{[14D]}\\ {}_{2,0}{p}_k^{\mathrm{L2},s}=\frac{1}{{}^{\mathrm{L2},s}}{a}_k^s{X}_k-f^{\mathrm{L2},s}{}_{2,0}{}_k+{}_{2,0}{\hat{N}}^{\mathrm{L2},s}+{}_{2,0}{v}_k^{\mathrm{L2},s}& \text{[15D]}\end{array}$

[0134] For the secondary baseline to the third base station, we have: $\begin{array}{cc}{}_{3,0}{p}_k^{\mathrm{L1},s}={}_{3,0}{}_k^{\mathrm{L1},s}-\frac{1}{{}^{\mathrm{L1},s}}{}_{3,0}{\tilde{R}}_k^s-\frac{1}{{}^{\mathrm{L1},s}}{}_{3,0}{T}_k^s& \text{[16C]}\\ {}_{3,0}{p}_k^{\mathrm{L1},s}=\frac{1}{{}^{\mathrm{L1},s}}{a}_k^s{X}_k-f^{\mathrm{L1},s}{}_{3,0}{}_k+{}_{3,0}{\hat{N}}^{\mathrm{L1},s}+{}_{3,0}{v}_k^{\mathrm{L1},s}& \text{[17C]}\\ {}_{3,0}{p}_k^{\mathrm{L2},s}={}_{3,0}{}_k^{\mathrm{L2},s}-\frac{1}{{}^{\mathrm{L2},s}}{}_{3,0}{\tilde{R}}_k^s-\frac{1}{{}^{\mathrm{L2},s}}{}_{3,0}{T}_k^s& \text{[16D]}\\ {}_{3,0}{p}_k^{\mathrm{L2},s}=\frac{1}{{}^{\mathrm{L2},s}}{a}_k^s{X}_k-f^{\mathrm{L2},s}{}_{3,0}{}_k+{}_{3,0}{\hat{N}}^{\mathrm{L2},s}+{}_{3,0}{v}_k^{\mathrm{L2},s}& \text{[17D]}\end{array}$

[0135] The terms in forms 12C, D, 14C, D, and 16C, D are the residuals of the differential navigation equations and are known. The unknowns are contained on the right-hand sides of forms 13C, D, 15C, D, and 17C, D. For N satellites, there are 6*N equations and a total of (66*N) solvable unknowns: _1,0k, _2,0k, _3,0k, the three components of X_k, and N instances of each of _1,0{circumflex over (N)}^L1,s, _2,0{circumflex over (N)}^L1,s, _3,0{circumflex over (N)}^L1,s, _1,0{circumflex over (N)}^L2,s, _2,0{circumflex over (N)}^L2,s, and _3,0{circumflex over (N)}^L2,s for s1 to sN satellites. However, using forms 11A, 11C, and 11D, the number of unknowns maybe reduced to 42*N. As we saw above, _2,0k may be related to _0,1k through the estimated base-station time offset _2,1k using form 11A as follows: _{2,0k2,1k0,1k2,1k1,0k}. In a similar manner, _3,0k was related to _0,1k through the estimated base-station time offset _3,1k as follows: _{3,0k3,1k0,1k3,1k1,0k}. Form 11C may be used to related each of the L1 band ambiguities of the secondary baselines to the primary baseline as follows:

[0136] Similarly, Form 11D may be used to relate each of the L2 band ambiguities of the secondary baselines to the primary baseline as follows:

[0137] Forms 15C and 15D for the second base station B2 may then be modified as follows: $\begin{array}{cc}\left({}_{2,0}{p}_k^{\mathrm{L1},s}+f^{\mathrm{L1},s}{}_{2,1}{}_k-{}_{2,1}{\hat{N}}^{\mathrm{L1},s}\right)=\frac{1}{{}^{\mathrm{L1},s}}{a}_k^{s}X-f^{\mathrm{L1},s}{}_{1,0}{}_k+{}_{1,0}{\hat{N}}^{\mathrm{L1},s}+{}_{2,0}{v}_k^{\mathrm{L1},s}& \text{[15C*]}\\ \left({}_{2,0}{p}_k^{\mathrm{L2},s}+f^{\mathrm{L2},s}{}_{2,1}{}_k-{}_{2,1}{\hat{N}}^{\mathrm{L2},s}\right)=\frac{1}{{}^{\mathrm{L2},s}}{a}_k^{s}X-f^{\mathrm{L2},s}{}_{1,0}{}_k+{}_{1,0}{\hat{N}}^{\mathrm{L2},s}+{}_{2,0}{v}_k^{\mathrm{L2},s}& \text{[15D*]}\end{array}$

[0138] Forms 17C and 17D for the third base station B3 may be similarly modified: $\begin{array}{cc}\left({}_{3,0}{p}_k^{\mathrm{L1},s}+f^{\mathrm{L1},s}{}_{3,1}{}_k-{}_{3,1}{\hat{N}}^{\mathrm{L1},s}\right)=\frac{1}{{}^{\mathrm{L1},s}}{a}_k^{s}X-f^{\mathrm{L1},s}{}_{1,0}{}_k+{}_{1,0}{\hat{N}}^{\mathrm{L1},s}+{}_{3,0}{v}_k^{\mathrm{L1},s}& \text{[17C*]}\\ \left({}_{3,0}{p}_k^{\mathrm{L2},s}+f^{\mathrm{L2},s}{}_{3,1}{}_k-{}_{3,1}{\hat{N}}^{\mathrm{L2},s}\right)=\frac{1}{{}^{\mathrm{L2},s}}{a}_k^{s}X-f^{\mathrm{L2},s}{}_{1,0}{}_k+{}_{1,0}{\hat{N}}^{\mathrm{L2},s}+{}_{3,0}{v}_k^{\mathrm{L2},s}.& \text{[17D*]}\end{array}$

[0139] Similar to form 18, forms 13C, D and modified forms 15C*, D* and 17C*, D* can be written in matrix form: $\begin{array}{cc}\left[\begin{array}{c}{}_{1,0}{p}_k^{\mathrm{L1}}\\ {}_{2,0}{p}_k^{{}^{*}\mathrm{L1}}\\ {}_{3,0}{p}_k^{{}^{*}\mathrm{L1}}\\ {}_{1,0}{p}_k^{\mathrm{L2}}\\ {}_{2,0}{p}_k^{{}^{*}\mathrm{L2}}\\ {}_{3,0}{p}_k^{{}^{*}\mathrm{L2}}\end{array}\right]=\left[\begin{array}{cccc}{}_{\mathrm{L1}}^{-1}{A}_k& -f^{\mathrm{L1}}& I_{NN}& 0_{NN}\\ {}_{\mathrm{L1}}^{-1}{A}_k& -f^{\mathrm{L1}}& I_{NN}& 0_{NN}\\ {}_{\mathrm{L1}}^{-1}{A}_k& -f^{\mathrm{L1}}& I_{NN}& 0_{NN}\\ {}_{\mathrm{L2}}^{-1}{A}_k& -f^{\mathrm{L2}}& 0_{NN}& I_{NN}\\ {}_{\mathrm{L2}}^{-1}{A}_k& -f^{\mathrm{L2}}& 0_{NN}& I_{NN}\\ {}_{\mathrm{L2}}^{-1}{A}_k& -f^{\mathrm{L2}}& 0_{NN}& I_{NN}\end{array}\right]\left[\begin{array}{c}{X}_k\\ {}_{1,0}{}_k\\ {}_{1,0}{\hat{N}}^{\mathrm{L1}}\\ {}_{1,0}{\hat{N}}^{\mathrm{L2}}\end{array}\right]+\left[\begin{array}{c}{}_{1,0}{v}_k^{\mathrm{L1}}\\ {}_{2,0}{v}_k^{\mathrm{L1}}\\ {}_{3,0}{v}_k^{\mathrm{L1}}\\ {}_{1,0}{v}_k^{\mathrm{L2}}\\ {}_{2,0}{v}_k^{\mathrm{L2}}\\ {}_{3,0}{v}_k^{\mathrm{L2}}\end{array}\right]& \text{[21]}\end{array}$

[0140] where:

[0141] _1,0p_k^L1 is an N1 column vector of the L1-band phase residuals _1,0p_k^L1,s associated with the primary baseline,

[0142] _2,0p_k^*L1 is an N1 column vector of the modified L1-band phase residuals (_2,0p_k^L1,sf^L1,s_2,1k2,1{circumflex over (N)}^L1,s) associated with the secondary baseline between the rover and base station B2,

[0143] _3,0p_k^*L1 is an N1 column vector of the modified L1-band phase residuals (_3,0p_k^L1,sf^L1,s_3,1k3,1{circumflex over (N)}^L1,s) associated with the secondary baseline between the rover and base station B3,

[0144] _1,0p_k^L2 is an N1 column vector of the L2-band phase residuals _1,0p_k^L2,s associated with the primary baseline,

[0145] _2,0p_k^*L2 is an N1 column vector of the modified L2-band phase residuals (_2,0p_k^L2,sf^L2,s_2,1k2,1{circumflex over (N)}^L2,s) associated with the secondary baseline between the rover and base station B2,

[0146] _3,0p_k^*L2 is an N1 column vector of the modified L2-band phase residuals (_3,0p_k^L2,sf^L2,s_3,1k3,1{circumflex over (N)}^L2,s) associated with the secondary baseline between the rover and base station B3,

[0147] A_k is the derivative matrix previously described,

[0148] _L1 is a diagonal matrix of the carrier wavelengths of the satellites in the L1-band, and _L1¹ is its inverse matrix,

[0149] _L2 is a diagonal matrix of the carrier wavelengths of the satellites in the L2-band, and _L2¹ is its inverse matrix,

[0150] f^L1 is a column vector of the carrier frequencies of the satellites in the L1-band, where _L1f^L1c{right arrow over (1)}, where c is the speed of light and {right arrow over (1)} is a column vector of 1's,

[0151] f^L2 is a column vector of the carrier frequencies of the satellites in the L2-band, where _L2f^L2c{right arrow over (1)},

[0152] I_NN is an NN identity matrix, and

[0153] 0_NN is an NN matrix of zeros.

[0154] The solvable unknowns X_k, _1,0k, _1,0{circumflex over (N)}^L1, _1,0{circumflex over (N)}^L2T may be estimated by applying a number of processes to the combination of forms 21 and 18. Further below, in a subsequent section, we described a preferred process that can be used on the combination of forms 21 and 18. Here, we describe how a least squares process may be applied to combined forms of 18 and 21.

[0155] To combine form 21 with 18, we rewrite form 22 to substitute c_1,0k for _1,0k as one of the unknowns, using the facts that _L1f^L1c{right arrow over (1)} and _L2f^L2c{right arrow over (1)}: $\begin{array}{cc}\left[\begin{array}{c}{}_{1,0}{p}_k^{\mathrm{L1}}\\ {}_{2,0}{p}_k^{{}^{*}\mathrm{L1}}\\ {}_{3,0}{p}_k^{{}^{*}\mathrm{L1}}\\ {}_{1,0}{p}_k^{\mathrm{L2}}\\ {}_{2,0}{p}_k^{{}^{*}\mathrm{L2}}\\ {}_{3,0}{p}_k^{{}^{*}\mathrm{L2}}\end{array}\right]=\left[\begin{array}{cccc}{}_{\mathrm{L1}}^{-1}{A}_k& -{}_{\mathrm{L1}}^{-1}\stackrel{}{1}& I_{NN}& 0_{NN}\\ {}_{\mathrm{L1}}^{-1}{A}_k& -{}_{\mathrm{L1}}^{-1}\stackrel{}{1}& I_{NN}& 0_{NN}\\ {}_{\mathrm{L1}}^{-1}{A}_k& -{}_{\mathrm{L1}}^{-1}\stackrel{}{1}& I_{NN}& 0_{NN}\\ {}_{\mathrm{L2}}^{-1}{A}_k& -{}_{\mathrm{L2}}^{-1}\stackrel{}{1}& 0_{NN}& I_{NN}\\ {}_{\mathrm{L2}}^{-1}{A}_k& -{}_{\mathrm{L2}}^{-1}\stackrel{}{1}& 0_{NN}& I_{NN}\\ {}_{\mathrm{L2}}^{-1}{A}_k& -{}_{\mathrm{L2}}^{-1}\stackrel{}{1}& 0_{NN}& I_{NN}\end{array}\right]\left[\begin{array}{c}{X}_k\\ c{}_{1,0}{}_k\\ {}_{1,0}{\hat{N}}^{\mathrm{L1}}\\ {}_{1,0}{\hat{N}}^{\mathrm{L2}}\end{array}\right]+\left[\begin{array}{c}{}_{1,0}{v}_k^{\mathrm{L1}}\\ {}_{2,0}{v}_k^{\mathrm{L1}}\\ {}_{3,0}{v}_k^{\mathrm{L1}}\\ {}_{1,0}{v}_k^{\mathrm{L2}}\\ {}_{2,0}{v}_k^{\mathrm{L2}}\\ {}_{3,0}{v}_k^{\mathrm{L2}}\end{array}\right]& \text{[22]}\end{array}$

[0156] This facilitates the combination of forms 18 and 22 as: $\begin{array}{cc}\left[\begin{array}{c}{}_{1,0}{B}_k^{\mathrm{L1}}\\ {}_{2,0}{B}_k^{{}^{*}\mathrm{L1}}\\ {}_{3,0}{B}_k^{{}^{*}\mathrm{L1}}\\ {}_{1,0}{B}_k^{\mathrm{L2}}\\ {}_{2,0}{B}_k^{{}^{*}\mathrm{L2}}\\ {}_{3,0}{B}_k^{{}^{*}\mathrm{L2}}\\ {}_{1,0}{p}_k^{\mathrm{L1}}\\ {}_{2,0}{p}_k^{{}^{*}\mathrm{L1}}\\ {}_{3,0}{p}_k^{{}^{*}\mathrm{L1}}\\ {}_{1,0}{p}_k^{\mathrm{L2}}\\ {}_{2,0}{p}_k^{{}^{*}\mathrm{L2}}\\ {}_{3,0}{p}_k^{{}^{*}\mathrm{L2}}\end{array}\right]=\left[\begin{array}{cccc}{A}_k& -\stackrel{}{1}& & \\ A_k& -\stackrel{}{1}& 0_{3NN}& 0_{3NN}\\ A_k& -\stackrel{}{1}& & \\ A_k& -\stackrel{}{1}& & \\ A_k& -\stackrel{}{1}& 0_{3NN}& 0_{3NN}\\ A_k& -\stackrel{}{1}& & \\ {}_{\mathrm{L1}}^{-1}{A}_k& -{}_{\mathrm{L1}}^{-1}\stackrel{}{1}& & \\ {}_{\mathrm{L1}}^{-1}{A}_k& -{}_{\mathrm{L1}}^{-1}\stackrel{}{1}& I_{3NN}& 0_{3NN}\\ {}_{\mathrm{L1}}^{-1}{A}_k& -{}_{\mathrm{L1}}^{-1}\stackrel{}{1}& & \\ {}_{\mathrm{L2}}^{-1}{A}_k& -{}_{\mathrm{L2}}^{-1}\stackrel{}{1}& & \\ {}_{\mathrm{L2}}^{-1}{A}_k& -{}_{\mathrm{L2}}^{-1}\stackrel{}{1}& 0_{3NN}& I_{3NN}\\ {}_{\mathrm{L2}}^{-1}{A}_k& -{}_{\mathrm{L2}}^{-1}\stackrel{}{1}& & \end{array}\right]\left[\begin{array}{c}{X}_k\\ c{}_{1,0}{}_k\\ {}_{1,0}{\hat{N}}^{\mathrm{L1}}\\ {}_{1,0}{\hat{N}}^{\mathrm{L2}}\end{array}\right]+\left[\begin{array}{c}{}_{1,0}{n}_k^{\mathrm{L1}}\\ {}_{2,0}{n}_k^{\mathrm{L1}}\\ {}_{3,0}{n}_k^{\mathrm{L1}}\\ {}_{1,0}{n}_k^{\mathrm{L2}}\\ {}_{2,0}{n}_k^{\mathrm{L2}}\\ {}_{3,0}{n}_k^{\mathrm{L2}}\\ {}_{1,0}{v}_k^{\mathrm{L1}}\\ {}_{2,0}{v}_k^{\mathrm{L1}}\\ {}_{3,0}{v}_k^{\mathrm{L1}}\\ {}_{1,0}{v}_k^{\mathrm{L2}}\\ {}_{2,0}{v}_k^{\mathrm{L2}}\\ {}_{3,0}{v}_k^{\mathrm{L2}}\end{array}\right]& \text{[23]}\end{array}$