METHOD AND SYSTEM FOR MONITORING LAND DEFORMATION FIELD OF THE INVENTION

The invention relates to a method and system for monitoring land deformation, particularly for use in monitoring landslide movement patterns.

BACKGROUND 0 F T H E I N V E NTI 0 N

The monitoring of landslide movement patterns remains a challenge despite recent availability of techniques such as laser-based geodetic techniques, GPS-based systems and ground/satellite-based radar.

Most, if not all, existing techniques in use rely on the availability of accurate location information of sparsely deployed monitoring points in the landslide region together with data from other sensing devices. The various types of sensors in use, including extensometers, in-place inclinometers, tiltmeters, pressure transducers and rain gauges, all tend to be expensive. This high cost tends to limit their deployment to only very high- risk sites.

A further disadvantage of prior systems is that they require a sustained power source. In practice a sustained power source is not always available meaning that GPS data are not continuously available. This in turn has the potential to affect adversely the accuracy of a typical GPS-based positioning system.

It is an object of at least preferred embodiments of the present invention to address some of the aforementioned disadvantages. An additional or alternative object is to at least provide the public with a useful choice.

SUMMARY OF THE INVENTION

In accordance with an aspect of the invention, a method of determining a relative position associated to a first receiver and a second receiver comprises receiving first location data associated to a first receiver; receiving second location data associated to second receiver; determining a first estimate relative position associated to the first receiver and the second receiver; determining a second estimate relative position associated to the first receiver and the second receiver at least partly from the first estimate relative position; determining a third estimate relative position associated to the first receiver and the second receiver at least partly from the second estimate relative position; and determining a fourth estimate relative position associated to the first receiver and the second receiver at least partly from the third estimate relative position. The term 'comprising' as used in this specification means 'consisting at least in part of. When interpreting each statement in this specification that includes the term

'comprising', features other than that or those prefaced by the term may also be present. Related terms such as 'comprise' and 'comprises' are to be interpreted in the same manner. In an embodiment determining the first relative position comprises determining a difference between a first autonomous position associated to the first receiver and a second autonomous position associated to the second receiver.

In an embodiment the first location data and/or the second location data comprise one or more of a common second based epoch, a satellite identifier, a receiver identifier, a time, an averaged autonomous position.

In an embodiment determining the second estimate relative position comprises receiving respective double difference estimates associated to at least two pairs of satellites;

forming an observation matrix associated to the double difference estimates; and defining a solution estimate at least partly from the observation matrix. In an embodiment determining the third estimate relative position comprises receiving at least one ambiguity estimate; fixing the at least one ambiguity estimate to respective integer ambiguity estimate(s); and defining a solution at least partly from the integer ambiguity estimates(s).

In accordance with a further aspect of the invention, a method of determining a relative position associated to a first receiver and a second receiver comprises receiving respective double difference estimates associated to at least two pairs of satellites;

forming an observation matrix associated to the double difference estimates; and defining a solution estimate at least partly from the observation matrix.

In an embodiment the method further comprises defining at least one ambiguity estimate.

In an embodiment the method further comprises defining a covariance matrix. In accordance with a further aspect of the invention, a method of determining a relative position associated to a first receiver and a second receiver comprises receiving at least one ambiguity estimate; fixing the at least one ambiguity estimate to respective integer ambiguity estimate(s); and defining a solution at least partly from the integer ambiguity estimates(s).

In an embodiment defining the solution comprises obtaining a first solution estimate; rounding the first solution estimate; and obtaining a second solution estimate at least partly from the rounded first solution estimate.

In accordance with a further aspect of the invention, a system for determining a relative position associated to a first receiver and a second receiver comprises a code module configured to receive location data associated to the first receiver, receive location data associated to the second receiver, and determine a first estimate relative position associated to the first receiver and the second receiver; a float module configured to determine a second estimate relative position associated to the first receiver and the second receiver at least partly from the first estimate relative position; a fix module configured to determine a third estimate relative position associated to the first receiver and the second receiver at least partly from the second estimate relative position; and a sidereal module configured to determine a fourth estimate relative position associated to the first receiver and the second receiver at least partly from the third estimate relative position.

In an embodiment the code module is configured to determine the first relative position by determining a difference between a first autonomous position associated to the first receiver and a second autonomous position associated to the second receiver.

In an embodiment the first location data and/or the second location data comprise one or more of a common second based epoch, a satellite identifier, a receiver identifier, a time, an averaged autonomous position.

In an embodiment the float module is configured to receive respective double difference estimates associated to at least two pairs of satellites; form an observation matrix associated to the double difference estimates; and define a solution estimate at least partly from the observation matrix.

In an embodiment the fix module is configured to receive at least one ambiguity estimate; fix the at least one ambiguity estimate to respective integer ambiguity estimate(s); and define a solution at least partly from the integer ambiguity

estimates(s). In accordance with a further aspect of the invention, a system for determining a relative position associated to a first receiver and a second receiver comprises a float module configured to receive respective double difference estimates associated to at least two pairs of satellites; form an observation matrix associated to the double difference estimates; and define a solution estimate at least partly from the observation matrix.

In an embodiment the float module is further configured to define at least one ambiguity estimate.

In an embodiment the float module is further configured to define a covariance matrix.

In accordance with a further aspect of the invention, a system for determining a relative position associated to a first receiver and a second receiver comprises a fix module configured to receive at least one ambiguity estimate; fix the at least one ambiguity estimate to respective integer ambiguity estimate(s); and define a solution at least partly from the integer ambiguity estimates(s).

In an embodiment the fix module is further configured to define the solution by obtaining a first solution estimate; round the first solution estimate; and obtain a second solution estimate at least partly from the rounded first solution estimate.

In accordance with an aspect of the invention, a computer readable medium has stored thereon computer-executable instructions that, when executed on a computing device, cause the computing device to perform a method of determining a relative position associated to a first receiver and a second receiver, the method comprising receiving first location data associated to a first receiver; receiving second location data associated to a second receiver; determining a first estimate relative position associated to the first receiver and the second receiver; determining a second estimate relative position associated to the first receiver and the second receiver at least partly from the first estimate relative position; determining a third estimate relative position associated to the first receiver and the second receiver at least partly from the second estimate relative position; and determining a fourth estimate relative position associated to the first receiver and the second receiver at least partly from the third estimate relative position. In an embodiment determining the first relative position comprises determining a difference between a first autonomous position associated to the first receiver and a second autonomous position associated to the second receiver. In an embodiment the first location data and/or the second location data comprise one or more of a common second based epoch, a satellite identifier, a receiver identifier, a time, an averaged autonomous position.

In an embodiment determining the second estimate relative position comprises receiving respective double difference estimates associated to at least two pairs of satellites;

forming an observation matrix associated to the double difference estimates; and defining a solution estimate at least partly from the observation matrix.

In an embodiment determining the third estimate relative position comprises receiving at least one ambiguity estimate; fixing the at least one ambiguity estimate to respective integer ambiguity estimate(s); and defining a solution at least partly from the integer ambiguity estimates(s).

In accordance with a further aspect of the invention, a computer readable medium has stored thereon computer-executable instructions that, when executed on a computing device, cause the computing device to perform a method of determining a relative position associated to a first receiver and a second receiver, the method comprising receiving respective double difference estimates associated to at least two pairs of satellites; forming an observation matrix associated to the double difference estimates; and defining a solution estimate at least partly from the observation matrix.

In an embodiment the method further comprises defining at least one ambiguity estimate.

In an embodiment the method further comprises defining a covariance matrix.

In accordance with a further aspect of the invention, a computer readable medium has stored thereon computer-executable instructions that, when executed on a computing device, cause the computing device to perform a method of determining a relative position associated to a first receiver and a second receiver, the method comprising receiving at least one ambiguity estimate; fixing the at least one ambiguity estimate to respective integer ambiguity estimate(s); and defining a solution at least partly from the integer ambiguity estimates(s).

In an embodiment defining the solution comprises obtaining a first solution estimate; rounding the first solution estimate; and obtaining a second solution estimate at least partly from the rounded first solution estimate.

The invention in one aspect comprises several steps. The relation of one or more of such steps with respect to each of the others, the apparatus embodying features of construction, and combinations of elements and arrangement of parts that are adapted to affect such steps, are all exemplified in the following detailed disclosure.

To those skilled in the art to which the invention relates, many changes in construction and widely differing embodiments and applications of the invention will suggest themselves without departing from the scope of the invention as defined in the appended claims. The disclosures and the descriptions herein are purely illustrative and are not intended to be in any sense limiting. Where specific integers are mentioned herein which have known equivalents in the art to which this invention relates, such known equivalents are deemed to be incorporated herein as if individually set forth. In addition, where features or aspects of the invention are described in terms of Markush groups, those persons skilled in the art will appreciate that the invention is also thereby described in terms of any individual member or subgroup of members of the Markush group.

As used herein, \s)' following a noun means the plural and/or singular forms of the noun.

As used herein, the term 'and/or' means 'and' or 'or' or both.

It is intended that reference to a range of numbers disclosed herein (for example, 1 to 10) also incorporates reference to all rational numbers within that range (for example, 1, 1.1, 2, 3, 3.9, 4, 5, 6, 6.5, 7, 8, 9, and 10) and also any range of rational numbers within that range (for example, 2 to 8, 1.5 to 5.5, and 3.1 to 4.7) and, therefore, all sub-ranges of all ranges expressly disclosed herein are hereby expressly disclosed.

These are only examples of what is specifically intended and all possible combinations of numerical values between the lowest value and the highest value enumerated are to be considered to be expressly stated in this application in a similar manner. In this specification where reference has been made to patent specifications, other external documents, or other sources of information, this is generally for the purpose of providing a context for discussing the features of the invention. Unless specifically stated otherwise, reference to such external documents or such sources of information is not to be construed as an admission that such documents or such sources of information, in any jurisdiction, are prior art or form part of the common general knowledge in the art.

Although the present invention is broadly as defined above, those persons skilled in the art will appreciate that the invention is not limited thereto and that the invention also includes embodiments of which the following description gives examples. The term 'computer-readable medium' should be taken to include a single medium or multiple media. Examples of multiple media include a centralised or distributed database and/or associated caches. These multiple media store the one or more sets of computer executable instructions. The term 'computer readable medium' should also be taken to include any medium that is capable of storing, encoding or carrying a set of instructions for execution by a processor and that cause the processor to perform any one or more of the methods described above. The computer-readable medium is also capable of storing, encoding or carrying data structures used by or associated with these sets of

instructions. The term 'computer-readable medium' includes solid-state memories, optical media and magnetic media .

The terms 'component', 'module', 'system', 'interface', and/or the like as used in this specification in relation to a processor are generally intended to refer to a computer- related entity, either hardware, a combination of hardware and software, software, or software in execution. For example, a component may be, but is not limited to being, a process running on a processor, a processor, an object, an executable, a thread of execution, a program, and/or a computer. By way of illustration, both an application running on a controller and the controller can be a component. One or more components may reside within a process and/or thread of execution and a component may be localized on one computer and/or distributed between two or more computers. The term 'connected to' as used in this specification in relation to data or signal transfer includes all direct or indirect types of communication, including wired and wireless, via a cellular network, via a data bus, or any other computer structure. It is envisaged that they may be intervening elements between the connected integers. Variants such as 'in communication with', 'joined to', and 'attached to' are to be interpreted in a similar manner. Related terms such as 'connecting' and 'in connection with' are to be

interpreted in the same manner.

BRIEF DESCRIPTION OF THE DRAWINGS

Preferred forms of the method and system for monitoring land deformation will now be described by way of example only with reference to the accompanying figures in which : Figure 1 shows a hardware embodiment adapted to monitor land deformation.

Figure 2 shows an example of a node from figure 1.

Figure 3 shows an example of a power management block from figure 2. Figure 4 shows a functional diagram of an example of a system for land deformation monitoring.

Figure 5 shows an example of distribution of the method between multiple devices.

Figure 6 and Figure 7 show a detailed view of an embodiment of a method for land deformation monitoring .

Figure 8 shows an example system that is configured to minimise data flow between the components of the system.

Figure 9 shows an example of an algorithm implemented by an embodiment of a node of figure 4. Figure 10 shows an example of an algorithm implemented by the creation and storage module of figure 6.

Figure 11 shows an example of an algorithm implemented by the float module of figures 4 and/or 6.

Figure 12 shows an example of an algorithm implemented by the fix module of figures 4 and/or 7.

Figure 13 shows an example of a method implemented by the sidereal module of figures 4 and/or 7.

DETAILED DESCRIPTION Figure 1 shows a hardware embodiment adapted to monitor land deformation. A base station 100 and a main computer 102. In an embodiment the base station 100 includes a wireless modem 104, for example a 3G modem, and the main computer 102 includes or is connected to an internet modem 106. The wireless modem 104 and the internet modem 106 exchange data over a network or collection of networks for example the internet 108.

In an embodiment the base station 100 is configured to gather broadcasted 'node algorithm data' from a wireless sensor network (WSN) 110 in addition to creating its own 'node algorithm data' from a physically attached GPS receiver 112.

This data is then processed to remove any time ambiguities and passed on to a server 114 so as to allow the main computer 102 access to this data. The main computer 102 is configured to take this data along with downloading of broadcast navigation data and execute one or more sets of computer-executable instructions to calculate position solutions.

The wireless sensor network 110 comprises a plurality of nodes for example node 116A, node 116B, and node 116C. In an embodiment the nodes 116 are configured to harvest solar energy, control their energy usage so as to allow as many satellite observations as possible, implement a node algorithm described below, and broadcast location data and/or 'node algorithm data'.

In an embodiment the wireless sensor network 110 includes at least one repeater for example repeater 118. The repeater 118 is configured to rebroadcast broadcast data so as to allow larger wireless sensor networks and/or wireless sensor networks with poor line of sight characteristics to be formed.

In an embodiment the base station 100 further comprises one or more of a WSN gateway, a base station computer, a mains power source. Figure 2 shows an example of a node 116. In an embodiment the node 116 includes a Power Management Block (PMB) 200. The Power Management Block 200 configured to store solar energy for example energy generated by a solar panel 202. The Power Management Block 200 and solar panel 202 collectively form a power supply.

In an embodiment the node 116 includes a microprocessor unit (MPU) 204, an 802.15.4 radio 206, and a GPS receiver 208. The MPU 204, the radio 206 and the GPS receiver 208 collectively form a host.

In an embodiment the PMB 200 stores energy in a super capacitor (not shown) and delivers it when needed to the host. In an embodiment the PMB 200 is configured to ensure that the super capacitor does not overcharge, manage current needs by information it is supplied by the host, ensure there is enough energy before supplying any power to the host and inform the host how much energy there is at any one time.

In an embodiment the host section controls operation of the node 116 and enables communication to the outside world.

The MPU 204 is partly responsible for energy management along with the power management block 200. In an embodiment the MPU 204 is configured to turn off and on the radio 206 and/or the GPS receiver 208 to conserve energy. In an embodiment the GPS receiver 208 passes raw and processed data from satellite observations to the microprocessor 204 so as to be able to implement a node algorithm on the node itself. Embodiments of the node algorithm are further described below.

In an embodiment the node includes an antenna comprising one or more of a passive helical antenna, an active patch antenna.

In an embodiment the node 116 has one of more of the characteristics low cost, does not include a battery, is solar powered, is compact enough to be held in one hand, able to operate even on dark cloudy days.

In an embodiment the radio 206 comprises an 802.15.4 radio being relatively long range while being cheaper, lower power and simpler than other standards such as Wifi.

In an embodiment the MPU 204 comprises a nanopower item from microchip having low power consumption and low cost.

In an embodiment the GPS receiver 208 comprises a singleband Ublox LEA6T module. The receiver 208 tends to be a relatively expensive component when compared with other components of the node 116. The receiver 208 also uses the most power when compared with other components in the node 116 (in the order of 100 mW). This kind of energy consumption is typical of small devices containing GPS receivers.

In an embodiment the GPS receiver 208 is configured to be put in some kind of off state or extremely low power state while still retaining ephemerides data or at least be able to be injected with ephemerides data from the MPU 204 when switched back on again.

In an embodiment the node 116 is enclosed within a housing (not shown). The solar panel 202 comprises one or more polycrystalline solar panels covering substantially all the surface area of the top of the housing as possible. Polycrystalline and/or

monocrystalline solar panels have been found to be suitable due to cost and efficiency with usual daylight conditions.

In an embodiment the PMB 200 is configured to link the solar panels 202 and the host section together.

Figure 3 shows a block diagram of an example of a power management block 200. In an embodiment the PMB 200 is configured to minimize quiescent power for times of low power demands so as to be able to store charge on the super capacitor even at times of low light levels. In an embodiment the PMB 200 has no Maximum Power Point Tracking (MPPT). The energy from the solar panel 202 is simply fed into the super capacitor via clamping circuitry to avoid the voltage over the capacitor going too high.

Upon the initial application of power to the super capacitor the PMB 200 is in an "off" state. The PMB 200 includes a boost converter 300 initially set to off.

A voltage monitor 302 is connected to the solar power. In an embodiment voltage monitor 302 is set to 2.3 V to avoid overvoltage on the super capacitor in addition to enabling the boost converter 300.

A voltage monitor 304 is set in an embodiment to 1.7 V which when triggered enables the boost converter 300. Voltage monitor 304 once triggered enables the boost converter 300. Filling of the local reservoir is then started. This filling of the local reservoir causes voltage monitor 304 to be disconnected. The capacitor voltage line will be connected to the super capacitor so as to allow the host to monitor the super capacitor voltage and vdd will be 3.3 V. When the host initially boots up it will be due to voltage monitor 304 which shortly gets disconnected due to vdd going high thus turning off the boost converter 300. This means that the host has to take immediate action as it is not using the super capacitor reservoir and instead is just using the local reservoir which is tiny in comparison and will not be entirely full. Either, the host must enter an extremely low power state waking up occasionally to enable the boost converter to replenish the local reservoir, else it must enable the boost converter to keep the local reservoir full. Either way it is configured to take the responsibility of keeping its vdd powerline from collapsing.

During normal times of operation when the host is sleeping and the boost converter 300 is disabled quiescent power is extremely low. This allows the solar panel 202 to deposit charge on the super capacitor. When the host wishes to turn on either the radio 206 or GPS receiver 208 it must keep the boost enable line active.

In an embodiment the host is configured to maintain a supper capacitor voltage of between approximately 1.4 and 2.2V. In an embodiment the power management block 200 further comprises a charging limiter. The quiescent power of the node 116 comprises the average power used by the node 116 at the nominal voltage of 2 V to only keep the power lines from failing while the node tries to use as little energy as possible. This means the host periodically refills the local reservoir but nothing else. The host needs to remain in a sleeping state occasionally waking itself up to refill the local reservoir while all the peripherals that it can control remain switched off.

Making the assumption that capacitor self-discharge is proportional to capacitor size, the quiescent power with respect to super capacitor in an embodiment is calculated as :

where Cs is super capacitor size in farads and PQ is quiescent power in microwatts.

In an embodiment approximately half the quiescent power goes to the super capacitor self discha rge when using the 33 F capacitor while for the 10 F capacitor approximately 20% goes to the self discharge of the super capacitor. In an embodiment a slight effect due to the different capacitors of 90μ\Λ/ or 130 μ\Λ/ are acceptably small that either ca pacitor could be used .

In an embodiment the peripherals associated to the node 116 comprise the radio 206 and the GPS receiver 208. As the boost regulator must be on when using the

peripherals, the voltage supplied to the peripherals is 3.3 V regulated .

As energy consumption is power times time, the radio 206 has the potential to use an insignificant amount of energy compa red to that of the GPS receiver 208. Generally the nodes 116 do not go into the RX mode and are entirely disconnected altogether using near zero power. The period of time that a node 116 turns on its radio 206 to transmit an epoch packet is in the 2 or 3 ms. A large part of this time is not in the TX mode but rather in some other mode while the oscillators of the radio stabilize. For a worse case scenario of 3 ms for packet using 495 mW sent once a second equates to 1.5 mW on average. At the same time the GPS receiver 208 needs to be on continuously to supply the radio epoch packet a nd uses around 123 mW of power; two orders of magnitude greater than the power req uired by the receiver.

Therefore, to a first approximation only the GPS receiver 208 uses energy when the receiver is active. The energy consumption of the radio 206 can therefore be considered negligible. When the GPS receiver 208 is inactive the radio 206 is disconnected and the GPS receiver 208 is in a backup power mode using 70 μ\Λ

The node 116 is configured to obtain observations and send them back to the base station. Active node power use is defined as the a pproximate average power at the nominal voltage of 2 V used by the node 116 during times of taking observations and sending them every second . This is energy used by all devices of the node, including the boost converter in the PMB 200, capacitor leakage, M PU 204, radio 206, GPS receivers 208, a GPS Low Noise Amplifier (LNA), and other devices.

In an embodiment the node 116 is configured to maximise the time that the GPS receiver 208 is in an operative 'on' mode taking observations and returning them to the base station 100.

In an embodiment the solar panel 202 is reasonably effective at powering the node 116 at low light levels. A solar panel having an MPP voltage close to 2V when light levels are low is suitable. Sola r panel typical voltages are typically given for direct sunlight conditions and an irradiance of lOOOWrrr ² [9] . The MPP voltage of the solar panel decreases when irradiance decreases, therefore, a solar panel with a rated typical voltage of more than 2 V is suita ble. While lOOOWrrr ² may be obtained in summer, even on clear days in winter depending on the location irradiance can be more tha n half of this. Even at times of year where the theoretica l maximum light levels can be large, from day to day light levels can alter dramatically.

A solar panel can be modeled as a constant current source over a d iode. The following equation shows the ubiquitous way of modeling the current obtained from a solar panel using just one d iode and consists of the superposition of the photocurrent and the dark current.

To a first approximation assuming the series resistance of the solar panel is small the short circuit current is approximately equal to the photo current. The photo current in turn is proportional to the irradiance. Therefore the short circuit current Isc is

approximately proportiona l to the irradiance . This permits a way of measuring irradiance once the proportionality consta nt is determined . The solar panel modeling equation as given is an implicit function that can be made explicit using the Lambert W function but generally it is just left in its implicit form.

Dark current iDark is the current that flows through the solar panel 202 due to an applied external voltage when no light is applied to the solar panel.

In an embodiment a blocking diode such as the ubiquitous 1N4148 is added to the solar panel 202 to reduce this current to near zero. Such a blocking diode reduces the MPP voltage by approximately 0.6 V which brings the nominal voltage of 2 V closer to that of the MPP voltage. However, as 2 V is considerably lower than the MPP voltage even at very low light levels, in an embodiment there is no blocking diode at all.

For applications that require the nodes 116 to function continuously at night and hosts that use very little energy it is desirable to use a blocking diode as this would

considerably reduce the quiescent power used by the node 116 and the solar panel 202 at night.

For applications that do not require the nodes to function at night or ones where the host uses considerable energy it is desirable not to use a blocking diode as this would increase the power delivered to the node 116 from the solar panel 202 while the 200 μ\Λ/ caused by the dark current of the solar panel 202 would be insignificant compared to that of the energy used by the node 116.

In an embodiment the GPS receiver 208 uses relatively considerable energy therefore it is desirable not to use a blocking diode.

Super capacitor sizing depends on the application. If the node 116 is to function at night than it is desirable for the super capacitors to be as large as possible to facilitate as many observations at night as possible while there is no sun.

However, if the node 116 is to function only during daylight hours and the role of the super capacitor is solely to allow duty cycling during the daylight hours then the size is less critical as a smaller one is possible.

Super capacitors up to 5000 F are cheaply available and sufficiently large. A GPS receiver using 240 mW has the potential to drop the capacitor voltage from 2.2 V to 1.4 V within about eight hours or say approximately the length of the night and is therefore approximately the minimum capacitor size needed for continuous night operation.

If continuous night and continuous day operation is required, approximately 500 mW on average would be required every day during daylight hours. As even a 5000 F super capacitor is unable to store enough charge for more than one night of continuous operation, even one day where the average daily power during daylight hours is not greater than 500 mW would require duty cycling.

If this 500 mW was to be obtained even on days of light levels as low as 20Wm ^~2 this would require a monocrystalline/polycrystalline solar panel roughly half a meter by half a meter. This is unacceptably large. Therefore, duty cycling algorithms are desirable.

Such duty cycling algorithms are configured to carefully decide who gets the limited energy from the solar panel 202 and how should the energy that has been stored on the super capacitor be used such that the voltage on it does not drop below 1.4 V, or such that the observations are as evenly distributed as possible throughout the 24-hour day period.

In an embodiment the super capacitor provides a method of allowing duty cycling during daylight hours rather than obtaining observations at night thus making the super capacitor size less critical.

GPS receivers tend to be subject to initialization time or Time To First Fix (TTFF) after they are turned on before useful observations are obtained. The GPS receiver 208 is configured to acquire and track the satellites and then resolve the code observable ambiguity. This acquisition time is between 1 and 26 seconds depending on many factors such as whether or not ephemerides have been retained, what timing information the GPS receiver has, how long has it been off for, and so forth. It is assumed that this value of between 1 and 26 seconds excludes the time required before the code observable ambiguity has been resolved as this can be as much as 6 seconds as shown in the code observable subsection A.4.1.

This initialization energy is energy spent not usefully obtaining observations. Therefore duty cycles that have long periods can help mitigate the total energy going towards continuously initializing the GPS receiver every duty cycle period. Duty cycle periods can be increased using larger super capacitors. The approximate wasted percentage of on time can be calculated as follows: where tm is TTFF, PGPS power of the GPS receiver, D Duty cycle, C capacitance, and v, and Vf the initial a nd final voltages over the super capacitor respectively.

Using la rger capacitors has the potential to red uce wasted time a nd the waste becomes greater when using low duty cycles. Using larger capacitors causes the observations to become more clumpy as well as a general trend in greater self-discharge of the super ca pacitors neither of which we want. However, compared to GPS power consumption the quiescent power tends to be insignificant for both 10 F capacitor and 33 F capacitor measurements. Both quiescent powers tend to be reasonably similar. TTFF for an autonomous GPS receiver can be somewhat unpredictable. There is no gua rantee that the GPS receiver 208 will ever acquire the satellite. In a n embodiment, keeping the backup mode of the Ublox LEA6T active continuously generally causes TTFF to become more consistent and shorter while using minimal power. Values of 10 seconds are typical . As the wasted percentage of 'on' time is a lso proportional to TTFF it is desira ble to implement a very short TTFF that is consistent. A 10 second TTFF tends to produce a lesser wasted percentage of 'on' time thereby reducing the benefits of using a larger capacitor.

In an embodiment, this method of initialization results in a typical time of 5 seconds before useful observations are obtained . A 5 second TTFF produces a lesser wasted percentage of 'on' thereby reducing the benefits of using a larger capacitor. However, ephemeris a nd time injections add to the complexity of a system.

In an embodiment a 33 F capacitor is used . This allows a reasonable amount of flexibility when testing types systems. It allows for a reasonably small wasted percentage of 'on' time due to initialization even without the use of ephemeris and time injections by permanently supplying GPS backup power instead ; generally less than 5% compa red to less than 20% for a 10 F capacitor.

Figure 4 shows a functional diagram of an exa mple of a system 400 for land

deformation monitoring . The system 400 is particularly suited to processing location data obta ined by nodes 116.

In an embodiment the system 400 includes a code module 402, a float module 404, a fix module 406, and a sidereal module 408. Each of the modules is designed to produce a position solution more accurate than the last. In an embodiment the float module 404, the fix module 406, and the sidereal module 408 require the previous module position solution as their input.

In an embodiment, each of the four modules produce solutions. For example, the code module 402 outputs a code solution, the float module 404 outputs a float solution, the fix module 406 outputs a fix solution, and the sidereal module 408 outputs a sidereal solution.

In an embodiment the float module 404 takes as input the location data from the nodes 116 or from another module rather than the code solution from the code module 402. In an embodiment the fix solution from the fix module 406 is output as the final solution from the system 400, thereby bypassing the sidereal module 408.

In an embodiment the code module 402 is configured to subtract course autonomous position solutions to produce an initial estimate of the final solution.

In an embodiment involving the float module 404, a real number N is given to each unique pair of satellites. Least squares (LS) solves for all Ns, the position solution, and the covariance information of both. For this the double differences are unwrapped using an approximate position estimate. If the code solution estimate is reasonably accurate (less than say 20 m) the wrapped double differences can be unwrapped. The wrapped double differences when corrected using the approximate position estimate don't change too rapidly even when the GPS receivers are separated by less than a few kilometers.

The float module 404 in an embodiment is configured to operate iteratively.

In an embodiment the fix module 406 is configured to take the real valued Ns and their covariance and perform Integer Least Squares (ILS) using the LAMBDA algorithm to fix the Ns to integers. Using LS, the fixed Ns are used to produce a fixed solution. The double differences are then rounded to the nearest solution plane and another solution obtained. Solution planes are the position solutions given one double difference for a particular integer ambiguity N.

In an embodiment, the sidereal module 408 is configured to perform sidereal filtering to improve daily solution precision. Absolute wrapped residual stacking with respect to a biased location is used in addition to fix solution stacking. Arithmetic averaging of both when possible results in daily solution precision improvements due to at least in part the random nature with which observations are obtained. Due to the random nature of obtaining observations, sidereal filtering has the potential to improve the precision of daily solutions. In an embodiment the improvement is at least partly due to the random and intermittent nature of the observations as well as an inability to obtain high rate observations continuously. The system 400 does not need code pseudoranges and the integer component of the phase observable to be sent over the radio link. This has the potential to provide a considerable reduction of bandwidth by removing these items before sending these things over the radio link. This results in smaller packets which in turn results in lower epoch losses thus mitigating the effect of a lossy channel. Navigation data is still required and need not be of high accuracy. Broadcast navigation is suitable and can be obtained from the satellites directly or downloaded from the Internet.

In an embodiment, system 400 is configured to perform a distributed algorithm referred to as the Code Float Fix Sidereal (CFFS) algorithm. The CFFS algorithm includes a node section and a main section. The main section comprises five main modules, namely a creation and storage module (not shown), in addition to a code module 402, a float module 404, a fix module 406, and a sidereal module 408.

The required data for the CFFS algorithm location data. In an embodiment the location data includes one or more of phase and code observations, time of reception as given by the receiver, broadcast data, navigation data.

The CFFS algorithm is configured to reduce or minimise use of bandwidth. This means some processing of this data has to be done before being sent over the network. In an embodiment the algorithm is split between multiple hardware devices reducing the data transmitted from one device to another. Figure 5 shows an example of distribution of the algorithm between multiple devices and the data input to the algorithm. The algorithm is configured to find a relative positioning solution for receiver 500 and receiver 502.

The receivers 500 and 502 can be regarded as data sources supplying the algorithm with its data. This data contains for example code, phase, time, and navigation information. The nodes 116A and 116B take the raw data, process it, and then send this processed data to the main section 504. In an embodiment the main section 504 comprises one or more of the modules from figure 4, for example the code module 402, the float module 404, the fix module 406 and/or the sidereal module 408. In an embodiment, each node 116 contains different phase, code and time data while navigation data is common to all.

Figure 6 and Figure 7 show a more detailed view of the CFFS algorithm. For example, the functions of the nodes 116, the creation and storage module 600, the float module 404, the fix module 406 and the sidereal module 408 are set out in more particular detail.

Figure 8 shows an example that is configured to minimise data flow between the node sections and the main section of the algorithm. Navigation data and code observations are not passed between the nodes and the main section. In an embodiment, what is passed comprise wrapped phase observations extrapolated to a common second based epoch 0|(t), a satellite identifier for example a satellite number S, a receiver identifier for example a receiver number B and the second that these observation estimates are for t.

In addition, occasionally averaged autonomous code based solutions x _B are passed to the main section along with the receiving number B. In an embodiment, code solutions generated by the code module 402 are required to be passed only once and hence contribute an insignificant fraction of data that flows between the node and the main sections.

Generally GPS receivers try to take their observations on the second and it has become somewhat of a de facto standard for measurements. For example the Ublox LEA-6T takes its observations on the second ±0.5 ms. Receiver Independent Exchange Format (RINEX) observation files are typically aligned to the second as well. In an embodiment the CFFS algorithm in an embodiment is configured to work on observations that are approximately aligned to the second.

Figure 9 shows an example of an algorithm implemented by an embodiment of the node section. The method 900 includes receiving 902 data from the receiver 208.

In an embodiment, the location data received from the receiver 208 comprise wrapped phase observations extrapolated to a common second based epoch 0|(t), a satellite identifier for example a satellite number S, a receiver identifier for example a receiver number B and the second that these observation estimates are for t. In addition, the location data occasionally comprises averaged autonomous code based solutions x _B , along with the receiving number s. The method includes determining 904, at least partly from the received data, an autonomous position solution XB. In an embodiment the method further includes determing a clock bias Af _B for a given reception time TB as given by the receiver.

GPS signals and carrier frequencies a re derived from a single time source and

synchronised with one a nother thus causing the signals to be bit phased with one another in addition to the signals being synchronised with one another at a higher level .

In an embodiment the instantaneous LI wavefront a satellite sends Wrx(t) at a time t, ca n be written as follows assuming the satellite clock keeps perfect freq uency:

W _TX (t) = A _TX 3l {s(t +Δ r _ra (t))e ² -W™ ⁺⁰ o ^X )}

In an embodiment the a bove equation is used to calculate XB and/or A† _B .

In an embodiment the method further comprises determining 906 a running average of the autonomous position solution XB. This running average is referred to as X _B .

In an embodiment the method further comprises determining 908 a value for t representing the second that the determined observation estimates are for. An estimate for reception time te is calculated as t _B = T _B — A† _B .

An iteration is used to determine satellite range at transmission time to reception time. In an embodiment the satellite ra nge is calculated for example by p _B (t _B ) = \\ s (t _B - The satellite tra nsmission time is calculated for example by tf = t _B - p _B ^s (t _B )/ c. The satellite's radial velocity is calculated for example by:

In an embodiment the value of t is initially set to te. The value of t is checked 910 for acceptability. If the absolute difference between t and te \t - t _B \ is above a threshold then the estimate determined for this second is discarded . The method returns to step 902 waiting for data from the receiver.

In an embodiment the threshold comprises approximately 200ms greater than which the estimate determined for this second is discarded . In an embodiment the threshold comprises approximately 500ms greater than which the estimate determined for this second is discarded . In an embodiment the signal strength of the satellite is also checked for acceptability. For example if it is too weak, unhealthy, and/or phase is not being tracked, then the estimate determined for this second is discarded. The method returns to step 902 waiting for data from the receiver. In an embodiment the method includes determining 912 an extrapolated wrapped phase ple equation for calculating this observable

In an embodiment the method includes transmitting 914 one or more of the wrapped phase observations extrapolated to a common second based epoch 0|(t), a satellite identifier for example a satellite number S, a receiver identifier for example a receiver number B and the second that these observation estimates are for t. In an embodiment these values are transmitted regularly.

In an embodiment the method further includes occasionally transmitting the averaged autonomous code based solutions x _B . In an embodiment the x _B values are not transmitted as frequently as the values transmitted in step 914.

In an embodiment the frequencies of the values transmitted regularly and occasionally are determined at least partly by the available amount of energy able to be harvested. In an embodiment a priority is to send 0f (t), S, B, and/or t values regularly as soon as they become available and if there is sufficient energy to send them. In an embodiment values for X _B are only transmitted after the 0f (t), S, B, and/or t values if there is sufficient energy remaining after transmitting the 0f (t), S, B, and/or t values. This means that the X _B values are transmitted occasionally.

In an embodiment the method pauses or stops when a node runs out of energy.

In an embodiment the main section of the algorithm takes the data that is given by the node sections and calculates a relative position solution.

Figure 10 shows an example of an algorithm implemented by the creation and storage module 600. In an embodiment the creation and storage module is configured to calculate a relative position solution at least partly from the data received from the node sections. In an embodiment the received data includes one or more of the wrapped phase observations extrapolated to a common second based epoch 0|(t), a satellite identifier for example a satellite number S, a receiver identifier for example a receiver number B and the second that these observation estimates are for t.

In an embodiment the received data includes the averaged autonomous code based solutions X _B . The method 1000 shown in figure 10 determines a wrapped double difference from wrapped phase observations. The creation and storage module receives 1002 node data, for example from the node sections.

In an embodiment the method includes extrapolating 1004 to a common epoch one or more of the wrapped phase observations ¾ (t), ¾ ² (t), 0B ¹ (t), ^(t). In an embodiment the method includes determining 1006 at least one possible wrapped double difference 0¾ ⁵¹ (t) from the wrapped phase observations. In an embodiment the method includes determining all such possible wrapped double differences.

In an embodiment the method includes determining 1008 a delta norm A _2S1 from navigation data. In an embodiment the creation and storage module stores in a memory one or more of the wrapped double difference ¾ ^si (t) the delta norm An _S2S1 , the common epoch t.

The code module 402 in an embodiment calculates the first relative position solution and is simply the difference of the two autonomous position solutions calculated in the previous node section. One example equation comprises letting P _1BA = X _B -X _A .

An embodiment of the float module 404 is configured to obtain a rough solution with an accuracy of a few centimeters along with the ambiguity N estimates and their covariance information for the fix stage.

Figure 11 shows an example of an algorithm implemented by the float module 404. An example model of a double difference is below:

¾¾T Si Λ";ϊί.· _; :ίί * ^:i * ^* * ¾ ^* ::i; ^: r * -iff ·> _: :? ^: ·<: ^* .

A wrapped double difference is given by the wrapped double difference creation algorithm, which can be unwrapped to re-create double difference estimates by the following equation: In an embodiment this unwrapping produces an integer N that is fixed for a ll observations to these two satellites for a ll times. It is assumed that iV j ^sl is fixed for all observations to these two satellites for all epochs tk. In an embodiment the unwrapped double differences given a position estimate P _BA where X = λ/Cf can be determined as follows:

In matrix form with many observations to many satellites the following is obtained where rBA is the residual vector:

The above set of equations has just one unknown integer N for each pair of satellites along with three unknown position variables. The number of equations can be much greater than the number of unknowns hence this form is suitable for being solved using LS which minimises the mea n squa re error of the residual . Being a linear method the integer N solutions end up being "floating" numbers hence the algorithm name. The unwrapping algorithm requires all possible double differences while the set of double differences required for the LS is a linearly independent set.

The independent set for the LS is created using a reference satellite, where one satellite for each epoch is chosen as the reference satellite and all double differences within that epoch must contain this reference satellite. The reference satellite is selected as the satellite with the highest elevation on an epoch by epoch basis. The satellite with the highest elevation is assumed to have the least noise and hence brings to the LS solution the least amount of error. This removes linear dependence a mongst double differences but however does not remove stochastic dependence amongst the double differences in each epoch . We ignore the stochastic dependence.

The above matrix equation has the form : A LS solution would result in the following estimates:

Solving for the solution using LS depends on the initial position estimate of P _BA for upwrapping. It is assumed that this procedure of first unwrapping the data given an approximate solution then solving for the solution produces a more accurate solution. The data can be unwrapped again given a better approximation of the solution.

This process can be repeated until the solution does not change much. Therefore, an iterative method is used feeding the output solution back into the initial position estimate to produce a new solution. This iteration is continued until convergence or divergence. The solution is assumed to be divergent when the iteration count becomes too large or clearly distinct solutions are repeating. Convergence is assumed when the input solution estimate and the output solution are very similar to one another. To help avoid convergence to incorrect solutions, the float algorithm is run twice. The first time it is iterated until convergence or divergence with the code solution as the initial solution, the second time it is iterated only once and yesterday's best solution is used as the initial solution for the float algorithm. The solution with the least mean square error of the wrapped residuals is chosen as the correct solution.

Referring to figure 11, in an embodiment the method 1100 first checks 1102 that there are sufficient wrapped double differences and that an initial solution P _IBA exists. If either of these conditions fails then the method fails.

In an embodiment the method includes obtaining 1104 double difference estimates. In an embodiment the unwrapping algorithm described above is used to unwrap at least one double difference given solution estimate P _BA to obtain double difference estimates <§T ^u (t _k ,P _BA ) for at least two pairs of satellites. In an embodiment all double differences are unwrapped. In an embodiment double difference estimates are obtained for all pairs of satellites.

In an embodiment the method includes selecting 1106 a reference satellite. The independent set of double differences is formed such that for each epoch the satellite with the highest elevation is selected as the reference satellite. In an embodiment the method includes forming 1108 an observation matrix using only the independent set of double differences formed above. In an embodiment the observation matrix is as follows:

In an embodiment the method includes defining 1110 an LS solution estimate and a covariance matrix. LS is used to find an estimate for the unknowns along with the covariance estimate.

An example of the resulting LS solution estimate is BA An example of the covariance

a

matrix is

The solution is checked 1112 for convergence or divergence. If there is not convergence or divergence the method returns to step 1104. Otherwise the solution is output 1114.

In an embodiment yesterday's best solution estimate is used as the initial solution for the float algorithm. The float algorithm is performed again with just one iteration. If the mean square error of the wrapped residuals of this solution is less than that of the iterative float solution then this solution is output.

In an embodiment the fix module 406 is configured to improve the accuracy of the float solution using the position estimate, the ambiguity estimate, and the covariance of the ambiguity estimate from the float section.

The fix solution fixes the integer N _s , also referred to as fixing the integer ambiguities, or just fixing the ambiguities. The float solution method is a linear one and hence the N _s when solved by the LS method are not integers as they should be.

If the Ns can be fixed to the correct integers then the number of unknowns is reduced from being the position vector plus all the N _s to just the position vector. This has the potential effect of producing more accurate solutions if the N _s is set to the correct values. For fixing the ambiguities, the covariance as well is the ambiguity estimates themselves obtained from the float solution a and Q _& are fixed to integers via ILS using a LAMBDA algorithm. The fixed ambiguities a are then used to correct the unwrapped double differences and hence remove the integer ambiguities from the above equation to produce the following :

LS can then be used on the above equation to obtain a fix solution estimate P _BA . The solution P _BA is then used to obtain another fix solution by rounding to the nearest solution plane a nd obtaining a nother LS estimate. In this step an integer ambiguity is given to each double difference and rounded to the nearest solution plane given the solution P _BA as follows where f () is the sawtooth difference function .

From the above equation, LS can be used to obtain a solution P _BA using the following equation :

This solution is also a fixed solution but in doing this rounding, any cycle slips introduced whilst unwrapping a re removed as this rounding is invariant to cycle slips.

It is expected that the mean square error of the residuals given the actual true solution would result in the lowest va lue for all given solutions. Because of this it is desirable to include one final step in the fix algorithm whereby the mean square error of the residuals for today's best estimate is compared with yesterday's best estimate.

If the mean squa re error of the residua ls using yesterday's best estimate results in a lower value than using today's, this means it is likely that for whatever reason that yesterday's best estimate is a better estimate for today than the current estimate for today. In this case it is desirable to round the double differences to the nea rest solution planes using yesterday's best estimate and then use least squares to find a new estimate for today. Figure 12 shows an example of an algorithm implemented by the fix module 406. The method 1200 checks 1202 for the existence of P _BA , Q _a , and . If one of these variables does not exist then the method 1200 fails.

In an embodiment the method uses ILS to find 1204 integer estimates a of a given ambiguity covariance Q _a -

In an embodiment the method corrects 1206 double differences based on integer ambiguity estimates a. An example of an equation to perform this correction is as follows:

In an embodiment the method uses LS to obtain 1208 a first solution estimate P _BA to the following equation :

In an embodiment the method rounds 1210 double differences to nearest solution planes given the solution P _BA . In an embodiment the equation is used : λΦ _Β ³ _Α ²³¹ {^, Ρ _ΒΑ ) = An _S2S1 (t _k ). P _BA + rf(C _f $ _B ^s _A ^2sl (t _k ) - ^An _S2S1 (t _k ). P _BA )

In an embodiment the method uses LS to obtain 1212 a second solution estimate the following equation :

In an embodiment the method compares the mean square error of the wrapped residuals using yesterday's best solution estimate Pyesterday with the mean square error of the wrapped residuals using P _BA . If the mean square error is less using yesterday's solution then P _BA is assigned the value of Pyesterday- The double differences are rounded to the nearest solution plane as set out in step 1210 to obtain a new estimate for P _BA as set out in step 1212. In an embodiment a daily sidereal filtered equation is set out below, where the subscript n represents day n, the bar notation is average over days of matching double

differences, P _{M n} is the average multipath bias error, P _n is the average of a priori known solutions and r _n is the double difference residual. All notation for intraday time (time within the same day), satellites and receiver nomenclature have been omitted for clarity. λ >4 444 : ^' :·.« «' :1Vf> - » 4 :>.« i's; .* 4

:4 ¾·;··■¾: -Si ¾·:···:.¾ :.:;.··.?;

The above equation is derived and explored more fully below.

In an embodiment the sidereal module 408 is configured to perform sidereal filtering where the satellite constellation approximately repeats every 24h -246 seconds also erroneously referred to as the GPS sidereal day. It is erroneous because the truth sidereal day length is approximately 24h -236 seconds on Earth.

Double differences that belong to a pair of satellites and a pair of receivers are considered to match if they are separated in time by an integer number of nominal GPS satellite orbital repeat periods and are on different sidereal days.

As multipath is a function of the GPS receiver's surroundings and satellite geometry, as long as the GPS receiver's surroundings do not change with respect to the GPS receiver then the multipath error should be near identical from one sidereal day to the next.

The bar notation is the arithmetic mean calculated when possible including the current day n. For rk at day n this is _n = me n(n u {r _n }) where Ω contains a set of matching residuals (the residuals of matching double differences); in this particular instance {r _n } represents the set containing only r _n and ensures _n is sensibly defined. This means there may be m days to average over but Ω may contain far fewer entries.

For calculating the sidereal filtered solution there is a need to calculate P _n , f _n , r _n , and (r _n - f _n ) . Incorrectly estimating the double difference integer ambiguity has the potential to introduce an integer ambiguity offset of the double difference residual.

Therefore, instead of calculating (r _n - f _n ) directly it is desirable to use an integer ambiguity invariant method using absolute wrapped residual stacking with respect to a biased location in addition to expected observation stacking using fix solutions from m consecutive days of the past including day n but excluding days k where ||P _fe - P _n || is excessive so as to prevent outliers from affecting the sidereal solution. P _n values are calculated from solutions obtained from the fix module. f _n are calculated considering only the absolute value of wrapped residuals, while for r _n only the wrapped residual is considered . To calculate (r _n - _n ) the sign of _n is chosen such to minimise the absolute magnitude of (r _n - f _n ).

In this way calculation of (r _n - _n ) is invariant with respect to the estimated double difference integer ambiguities. The above equation represents one double difference, many equations are needed to use LS therefore it is desirable to use matrix notation and describe the implementation of the algorithm as follows where a _P is the approximate satellite constellation repeat period (approximately 24h -246s as already stated) .

Figure 13 shows an exa mple of a method 1300 implemented by the sidereal module 408. The value of PB is first initialised by selecting it as the first ever fix solution obtained from the fix module.

In an embodiment the method includes shifting 1302 wrapped residuals rk temporarily by an integer multiple of a _P to a lign them with today, day n.

In an embodiment the method 1300 determines 1304 biased wrapped residuals. An example of an equation to perform this function is:

In an embodiment the method excludes 1306 outliers. One technique for excluding outliers is to exclude all days k from the m days where ||P _fe - P _n || is above a threshold . This has the potential to prevent outliers affecting the sidereal solution.

In an embodiment the threshold is approximately 1cm. In an embodiment the threshold is approximately 2cm.

In an embodiment the method determines or obtains 1308 averages. One example of an average is the average of a ll \r _k \ over all matching double differences from the m days to produce A further exa mple of an average is the average of all N. P _k over all matching double differences from the m days to produce N. P _n .

In an embodiment the method determines 1310 an estimate. For example _n =

magmin(r _n - |rjj, r _n + jr^j). LS is then used to solve (xf _n + N. P^j = N. (p _n + P _{M n} ) + e _n and hence obtain an estimate of Pn. This estimate is referred to as P _n . The function magmin (x, y) selects the value to return of the argument with the least absolute magnitude. It is acknowledged that each satellite repeat period differs from one satellite to another. The value of a _P is chosen to be a fixed value for simplicity, for example 24h -246s. For both the float and the fix stages, sets of linear equations are solved using LS. As each double difference can be created as the addition or subtraction of two other double differences Φ¾ ⁵¹ = Β ^{3 2} = Φ¾ ⁵¹ if all double differences were used in the equations, biased solutions could result as each single difference contributes a different amount of noise. In an embodiment, an independent set is used to avoid this. To decide which double differences stay and which ones a re removed, a reference satellite is used . A reference satellite is the satellite in each epoch that must appear in each double difference. This produces a star pattern of double differences where each double difference has been ta ken with respect to this reference satellite. Which reference satellite is used may make a large difference to the position solution; this is because the reference satellite brings its noise into each equation in the set of equations.

Therefore selecting the reference satellite to be the one with the lowest noise has the potential to result in the most accurate solution . A typical method of selecting the reference satellite is to choose the one with the highest elevation due to the assumption that satellites with high elevation have low noise and low multipath.

The selection of the reference satellite can alter the accuracy of the final solution and is critical for precise difference based solutions. Both the fix and float stages of the algorithm use the idea of a reference satellite and use the same set of independent double differences. Choosing a satellite with the highest elevation as the reference satellite has the potential to result in a smaller variance in all directions than when randomly choosing a reference satellite.

Appendix A titled 'GPS basics' provides backg round reading on the basics of GPS signals.

Appendix B titled 'Differencing the phase observable' expla ins more specific aspects of solution processing .

Appendix C titled 'Node functionality at night' provides details of experiments used to obtain observations at night. Appendix D titled 'Network' describes the structure and the algorithms used to create an embodiment of a communications network.

Appendix E titled 'Chapter 3 - Multipath and sidereal filtering' describes the form that the solutions take in an embodiment both before and after sidereal filtering has been implemented.

Appendix F titled 'Chapter 4 - Algorithm comparison' compares the solutions obtained using an embodiment disclosed above with solutions obtained from other solution processing software.

Appendix G titled 'Chapter 6 - Testbed results and analysis' describes results and analysis of two experiments.

The foregoing description of the invention includes preferred forms thereof. Modifications may be made thereto without departing from the scope of the invention, as defined by the accompanying claims

In this appendix a general introduction is given to GPS operation and in particular, aspects that relate to acquiring observations for solving position solutions using singleband receivers.

As of writing GPS is currently undergoing a modernization to improve both civilian and military use. Between 1990 and 2004 legacy satellites were launched while from 2005 modernized satellites have been launched. According to the National Coordination Office for Space-Based Positioning, Navigation, and Timing [16] this is in an effort to upgrade the features and performance of GPS. Currently GPS transmits on three different RF links from the satellites to end-users. The RF links are called LI, L2 and L5 and are named after the bands that they transmit in. Code Division Multiple Access (CDMA) is used as the channel access method so all satellites used the same carrier frequencies. LI has a nominal frequency of 1575.42Mhz as seen from Earth, L2 1227.6QMhz and 1.5 1176.45Mhz. These nominal frequencies are modulated with various signals to aid navigation . The current GPS modernization consists of generally improving the hardware as well as adding more signals that are sent over the RF links. The GPS modernization currently underway will take many years and satellites producing signals such as LlC on LI axe not expected to be launched until 2016 with 24 satellites expected by around 2026 [31] . As we are APPENDIX A. GPS BASICS interested in only singleband receivers because they are cheaper than multi- band receivers and that no new modernized signals are currently in operation on LI, we restrict our investi ation here to that of legacy LI signals.

Figure A.l : Legacy Li signal generation block diagram

GPS legacy signals on LI consist of a coarse acquisition Pseudo Random Number (PRINT) code called C/A sent a,t 1.023Mb/s, a precise PRN code P sent at 10.23Mb/s which is called Y if encrypted, and navigation data called NAV sent at 50b/s. The two PRN codes are unique to each satellite and each spread the navigation signal. Bits of the PRN codes are also called chips. LI consists of in-phase and quadrature components. Each component is separately modulated using Binary Phase Shift Keying (BPSK) as their modulation technique. One modulator is supplied with a bit train from modulo-2 addition of P(Y) and NAV while the other is supplied with a bit train from modulo-2 addition of C/A and NAV. The signals and the carrier frequencies are derived from a single time source and synchronized with one another thus causing the signals to be bit phased with one another in addition to the signals being synchronized with one another at a higher level. Figure A. l depicts a block diagram of the generation of such signals. The instantaneous LI wavefront a satellite sends W ^' X if) at a time £, can be written as follows assuming the satellite clock keeps perfect frequency. APPENDIX A, GPS BASICS

Where S (t) = NAV it) (P (Y) (t) e ⁱ⁰ - C (t) e ^i7r/'2 ) is the composite signal, C (t) C/A code, P (Y) (t) P(Y)-code, NAV (t) navigation data, φξ ^χ satellite oscillator phase at time zero, Α χ transmission amplitude, φ ^χ satellite oscillator phase at time zero, frx nominal satellite oscillator frequency as seen from Earth so as to account for general reiativistic effects, JQ 1575.42Mhz as seen from Earth, and TJ _} ^X is the satellite composite signal offset at time zero. C, P (Y) ; NAV £ {— 1, 1} and are functions so as to produce the correct composite signal where a mapping of 0 --- I and 1 --> — 1 has been applied. At the receiver the instantaneous LI wavefront a receiver receives WRX (t) at a time t, can then be written as follows assuming no hindrance by the atmosphere.

WRX (t) A _RX {S [t + AT _J X (t) - At it)) ₆ « ^~ «Μτ _Χ+ Φ ) ] (A.2)

Where ARX is reception amplitude and At it) is the transmission flight time from the satellite at transmission time to the receiver at reception time t. It's important to realize tha,t Δί (t) is not the flight time between the satellite and receiver at time t, but rather the flight time based as the receiver sees it. It is similar to w ^T hen a airplane passes by and one hears the sound of the plane lagging where the plane actually is. The flight time of the sound from the plane as determined by the listener is different from the flight time one would get from calculating where the plane actually is to the user with respect to the same reception time.

APPENDIX A. GPS BASICS

Figure A.2: Simplified version of a possible GPS satellite acquisition and tracking scheme

Figure A.2 shows a simplified version of what a GPS satellite receiver can do to acquire the satellite and track the signals that it produces. Upon reception the receiver tries to separate the signals from LI and track the C /A code and the carrier phase. This can effectively be performed by mixing WRX (i) with a, complex local oscillator LO (t)— where f _RX is the frequency of the oscillator and φ ^χ the phase of the oscillator at time zero, filtering using a low pass filter LPF, then phase and/or frequency tracking to stop rotation along with correct C/A timeing to match the C/A code as sent by the satellite so as to acquire access to the data being sent by the satellite.

A, 2, 0.1 Obtaining the baseband signal R(t )

After mixing with the local oscillator and filtering using the low pass filter, the receiver obtains the following baseband, where, Δ/ = fax — J _T X and

R (t LPF {W _TX (t) LO ) }

A _RX

S it + AT, (t) ~~ At it)) _E 2* ^i(t f+W ⁾ frx- APPENDIX A, GPS BASICS

We define tAf -I- Δί ( /;) frx + Φο as the received beat carrier phase (carrier phase).

A.2.0.2 Stopping rotation

We see that this is a constellation of four points that rotates due to the frequency difference between the receiver's local oscillator and the satellite's oscillator, and also rotates due to the radial motion of the satellite itself with respect to the receiver.

If we let Φ (t) — tAf + At (t) frx + Αφο then we can correct for rotation by multiplying i? (t) as follows.

S (t + AT _TX (ί) - Δί (ί)) β 2 ^: -τ¾(Φ(ί) ) -27τ?:(Φ(ί))

R it) e

S (t + AT _TX (έ) - Δί (ί)) (A.4)

This stops the constellation from rotating and removes any constant constellation rotation offset. This then resolves the composite signal. Hence we define Φ (t) as the estimated received beat carrier phase by the receiver (estimated carrier phase). The constellation's phase (or equivalently the constellation's rotation offset) is defined as the difference between the carrier phase and the estimated carrier phase tAf + At (t) frx + Αφ ₀ — Φ (t). More generally as long as the estimated carrier phase is in phase with the carrier phase the constellation stops rotating and there is no constant constellation rotation offset.

A.2.0.3 C/A Code alignment

After the constellation rotation has been stopped by letting the estimated carrier phase be in phase with the carrier phase, a local replica of the C/A code has to be mixed with the composite signal and phase shifted in time until the local replica of the C/A code is in phase with the one that is in the received composite signal. We then define r (t) as the C/A code alignment offset and is how much the incoming C/A is misaligned with the local replica. We define the local replica of the C/A code as LCA (t)— C (t + AT _RX ({ )) APPENDIX A. GPS BASICS where Δ¾ (/. j is the receiver's clock offset. First we notice that if we offset the local C/A replica by r (t) and let r (t) = Δί (t) + AT _RX (t) - ΔΤ _Τ χ (t) and multiply this offseted C/A with R (t) we expect the following for a random time t.

- E [LCA [t - T (i) ) S (t + ΑΤτχ (i) - Δί (i))]

= E [C (t + AT I (ί) - Δί (ί)) 5 (ί + ΔΓ _ΤΧ (t) - At (t))] fa %E [NAV (t + A rx (t) - At (t))j (A.5)

This is due to the fact that P(Y) and C are not well cross-correlated, while of course C is perfectly correlated with itself.

LCA (t— T (t)) R (t) β ^"·2 ™^* ⁾⁾ is the input to the acquisition and tracking block in figure A.2. Visually LCA (t - r (t)) R (t) e ^{' 2} ^^ is a constellation without any rotation of four points. The two points lying on the imaginary axis move slowly at no more than 50 times a second and contain the navigation data, while, the two points that lie on the real axis move very rapidly at up to 10.23 million times a second in a seemingly random way with a mean value of zero. More generally when this happens the incoming C/A code is aligned with the local replica and the carrier phase is in phase with the estimated carrier phase.

A.3.0.1 Extraction of NAV data, using filtering

When the incoming C/A code is aligned with the local replica and the carrier phase is in phase with the estimated carrier phase, navigation data can easily be obtained simply by using a low pass filter APPENDIX A, GPS BASICS

NAV (t + AT _TX (t) - At (t)) = -iLPF \ LCA (t - τ (i)) R (i) e ^" - ^w( * ^(f))

(A.6)

This allows the GPS receiver access to the satellite's navigation data which includes a wealth of data including the satellite's estimate of ΑΤ χ (t) and time of transmission of specific navigation data transitions. Along with data to calculate the satellites oositions.

A.3.0.2 A metric for C/A code alignment and rotation

We restrict t to be within a small period of time 2δ which is less than NAV's period while still ensuring a long enough period such that P(Y) and C are still not well correlated over that period. To find such a time period is possible as NAVs period is 204600 times longer than P(F )'s period. We assume the C/A codes are perfectly alinged and the the carrier phase is in phase with the estimated carrier phase. We then expect the following.

E \ LCA (t - r (<)) R (t) e ^~2wi W ^l »] « iNAV (t + AT _TX (t) - At (i)) (A. *7\

k=i+6

LCA (k (k) ) R (k) e ^ ⁿ ^>dk (A

2δ k=t-S

As the navigation data has a constant magnitude of 1, any imperfections in the correlation between LCA (t - r (t)) and C (t + AT _TX (t)— At (t)) due to incorrectly estimating r ( ) will affect the magnitude of the acquired navigation data. In addition, incorrectly estimating Φ (k) causing a rotating constellation motion, will, after integration also negatively affect the magnitude of the acquired navigation data. Therefore, we can say the following where is defined as the correlation coefficient between a local time shifted replica of the C/A code and the one being received that may be rotating. APPENDIX A. GPS BASICS

This coefficient more generally can be calculated more accurately with averaging; therefore a GPS receiver could calculate it as follows.

Squaring eliminates the sign of the BPSK. The magnitude of 7 is related to how well the C/A codes are aligned and how well the rotating motion of the constellation has been stopped, while the angle of 7 is related to the constellation's constant rotation offset with a half cycle ambiguity. A constant constellation rotation offset does not affect the magnitude of 7 as a constant rotation offset is just a constant that can be taken out of the integral. Therefore, gamma will be maximized when the C /A codes are perfectly aligned and the constellation is not rotating while it is invariant for constant rotation offset of the constellation.

• /jis maximized when C/A codes are aligned and the constellation is not rotating.

Treating Φ and τ as variables a GPS receiver can vary Φ and τ to maximize 7. If 7 is above a certain threshold the GPS receiver can assume that the satellite is acquired and to commence tracking τ, Φ, and decoding NAV data.

We are interested in maximizing 7 because when it is a maximum, with the addition of some ambiguity both Φ and r are good estimates for tl j + At (t) fox + φο and At (t) + ATRX (t) — ΑΤτχ (t) respectively which turns out to be useful in finding positions solutions. In addition when 7 is maximized we are able to obtain navigation data which is also useful for finding position solutions. APPENDIX A. GPS BASICS

A .3.0.3 First-order linear approximations of unknown fnnctions r

and

We wish to maximize 7 , the correlation coefficient. Due to the surface shape of the magnitude of 7 with respect to and τ initial estimates for both Φ and r are required. Without good initial estimates. 7 is dominated by noise making standard tracking schemes such as Phase Locked Loop (PLL) , Frequency Locked Loop FLL) , Delay Locked Loop (DLL) and early/late time useless. It's like trying to track an ant crawling in long grass; you have to find it first before you can track it as the grass makes it difficult to see the ant from afar.

First we create a first order linear approximation model of how Φ and r change. We have already seen that when Φ (t) = tAf + At it) f _TX + Δ<¾ and τ (t) = At (t) + ATux (t) — ΑΤτχ (i) we are able to stop rotation and align our local C /A code replica with the incoming one. Therefore, these are the Φ and r that we are looking for. Linear approximations of these two equations are written below where F (t _m ) = (Δ/ + Af _TX (t _m )) , Af _TX (t _m ) is the change of frequency of f _TX due to Doppler at time t _m , where a, positive value is for the satellite moving away from the receiver and Θ (t _m ) and Ξ (t _m ) are some constants. Proofs can be found in A.5.1 and A.5.2.

; ) « tF (t _m ) + Θ (A. l l

(t) ?a it - t _r F <7 (A.12)

These approximations are only valid if F (t _m ) does not change to rapidly around time t _m . The maximum rate at which velocity will change is about 0.1178 nis ^'"2 and is when the satellite is directly overhead (A.5.4) . On the LI band this implies that the constellation rotation speed will change by less than about 0.9 Hzs ^'1 if F U _m ) is left unchanged in equation A.13 ( see A.5.3 and A.5.4) . Now compare this to the range of F (t _m ). A fox {t _m ) can be as large as about d::5 kHz (A.5.4) , and depending on the receiver clock accuracyA could be out by another 5 kHz if we assume a receiver clock accuracy of 3.5 ppm. This means the range of F (t _m ) is in the order of 20 kHz. If we restrict our time of interest to 1 ms, then, F (t _m ) will change by less than 0.0009 Hz which APPENDIX A. GPS BASICS is far less than the range of the 20 kHz of F _m ) . Therefore, for short periods of time this is a valid approximation.

Equations A.11 and A.12 form an approximate model of how Φ and r will change in the short term.

A.3.0.4 Further simplifications to the first-order linear approximations of unknown functions τ and Φ in regard to acquisition

In figure A.2 the input of the acquisition and tracking block as we have already seen is written as in equation A.13.

LCA it '2πί(Φ(ί)) (A.13)

To acquire Φ and r initially we can make further simplifications to equations A.11 and A.12. Θ (t _m ) in equation A.11 when placed in equation A.13 has no effect on changing constellation rotation with respect to time. A.13 will still be a non-rotating four-point constellation but just with a constant rotation offset of Θ (t _m ) cycles. As we have already mentioned in equation A.10, any constant rotation offset has no effect on the magnitude of the correlation coefficient. Therefore we can ignore Θ (t _m ) when initially acquiring Φ.

The carrier wave frequency is 1540 times greater than the bit rate of the C/A code. As frequency times time is phase, if F = 10, 000 Hz then it takes 0.025 vcis for Φ to change by a quarter of the cycle, while it takes 38.5 ms for the C/A code to change by a quarter of a chip. If we then assume a digitalization of R at the rate of 4.092 Mb/s and sampling 1 ms worth of R, more often than not we couldn't even detect the difference between F ----- 0 and F ----- 10, 000 in the C/A code directly while it would be easy to detect in the carrier wave.

Due to these two points we make the following two approximations when considering initial acquisition of Φ and r. Here we acknowledge that the constellation will have an arbitrary constant rotation offset, are only valid for periods of time of a few milliseconds, and baseband sampling rate is no more than a few times per chip.

Φ it) « tF (t _m ) (A.14) APPENDIX A, GPS BASICS

A.3.0.5 Acquisition

With these two approximations A, 14 and A.15 the receiver can do a two- dimensional search, F with a frequency dimension and the other Ξ with a, time dimension to find the point that maximizes the correlation coefficient in equation A.10. Assuming a receiver clock accuracy of 3.5 ρρπι, the receiver would have to search from—10 kHz to 10 kHz. As the C/A. code is a periodic function with a period of 1 ms, the receiver would have to search from 0 ms to 1 ms.

Estimating F and Ξ by trying to maximize the correlation between the local C/A code replica with the incoming one turns out to be computationally demanding using more energy than tracking, and is a major concern for GPS receivers. Because of this much research has been directed towards this problem to reduce the computational effort to estimate these two parameters [63]. A parallelized 2D search by using Fast Fourier Transform (FFT) is a conventional method currently used in software defined receivers [63] .

As an example we cross-correlated a local C /'A code replica with an incoming one on LI . We used a sample rate of 4.092 Mb/s being four times the nominal frequency of the chip rate and searched by varying the frequency term by ±10 kHz in steps of 125 Hz and then using cross correlation varying the time term by no more than ±511.5 chips (±511.5 chips covers the entire 1 ms). Navigation data was simulated by using random data, while no P(Y) data was added. An offset of 5000 Hz was applied to the 1.57542 GHz carrier frequency and a phase offset of 250.15 chips was applied to the C/A code of the incoming signal. Using 1 ms worth of sequential incoming data, figure A.3 was obtained. APPENDIX A. GPS BASICS

Figure A.3: Correlation coefficient versus phase and frequency offsets. 1 ms of sequential data of C/A PRN 16 cross-correlation. 4.092Mb/s sampling rate. 16 kb of data, 658,651 evaluation points.

As can be seen there is a clear peak in the graph representing the estimates of F and Ξ . The estimated phase offset using interpolation around the highest peak was 250.18 chips while that of the frequency was 5000 H 2; this matches well with the exact values. From this figure it is clear to see why tracking will not work without good initial estimations of τ and . It's interesting to note in passing that the 0.03 chips that the interpolation value was out equates to the time light takes to travel 9m. This is roughly the correct order of accuracy cheap consumer grade GPS receivers have.

The number of points needed to be evaluated in figure A.3 was 658, 651. While less points could have been used and still be able to obtain reasonable estimates for F and Ξ , by the 2-D search method there are inherently always going to be a large number of points needed to be evaluated. While other methods such as the parallelized 2D search by using FFT exist and are less computationally demanding and hence less energy demanding, acquisition to the best of our knowledge still uses more energy than tracking. As an example the Ublox NEO-7N which is a modern consumer grade GPS receiver uses 23% more energy during its acquisition state than it's tracking state [17]. APPENDIX A. GPS BASICS

A .3.0.6 Tracking

Once the acquisition has been performed and the point of the maximum correlation coefficient has been found, the receiver can then track r and using standard techniques such as PLL, FLL, DLL and earl /late tracking methods. The linear approximations do not stop the maximum correlation coefficient point from moving but it slows it down sufficiently that one can treat it as a stationary point until it has been found then one can simply track it as it moves.

From acquisition, estimates for F and Ξ are obtained, which via equation A.12 gives an estimate for r (' ) . To track r (£) then early/late tracking can be used. Such a method usually consists of three locally produced C/A replicas, one slightly ahead of what is expected from the satellite, one as expected from the satellite, and one slightly behind what is expected from the satellite.

Earlv LCA (t ~ r it) + ξ)

Prompt LCA (t— τ (t))

Late LCA (t - r (t) - ξ)

These three C/A replicas are then each correlated with R (t) e ^~2W H*i* ⁾ ) to produce three correlation coefficients ( E ,Ί ^' Ρ, 7Z. ) and using interpolation a new estimate for r (/;) can be obtained to keep the code aligned. Keeping the code aligned is the one of the two requirements for maximizing 7.

LCA (i - r (t)) R (t) has the effect of removing the C/A code from R. This removal is called wiping the code. Once it is removed a carrier tracking scheme such as a PLL can be used on LCA \ i— r (t) ) R {t) because in one direction it appears as a standard BPSK signal. A costas PLL could be performed on LCA (t— T (t)) R (i) to estimate Φ (t) as it is invariant to the navigation transitions. Assume that we designed the costas loop to align on the imaginary axis. Then, the costas loop will align NAV's BPSK signal along the imaginary axis with an ambiguity as to which way around it is aligned. The costas APPENDIX A. GPS BASICS loop will also stop the constellation from rotating which is one of the two requirements for maximizing 7.

So, using code tracking and carrier tracking simultaneously the maximum point of correlation can be continuously tracked. It is not sufficient for a receiver solely to track only one of r (t) or Φ (i) ; both need to be tracked simultaneously.

We have already seen that in A.3.0.2 r (t) = At (t) + AT _RX (t) - AT _TX (t) implies maximum correlation. Due to the 1 ms C/A ambiguity the converse is not true. Therefore we can say that maximum correlation implies r it) —

At (t) + AT _RX (t )— ΑΤτ' χ (t) + / 10(3(3 where M is some fixed integer.

Likewise Φ (/. j = t/Af + At (t) frx + φο implies maximum correlation but due to carrier phase cycle ambiguity and that gamma is maximized for any constant rotation offset the converse is not true. The costas loop removes the constant rotation offset with an ambiguity of half a cycle, and if the PLL accumulates its phase offset rather than resetting it as it passes through an angle of zero, then maximum correlation implies Φ (t) — tAf + At (t) f _TX + φο + N/2 for some fixed integer N. r and Φ when these ambiguities are considered become the two observables used by almost all low end consumer grade GPS receivers for position solution calculations.

A, 4 Observables

Observables are measurements taken by the GPS receiver of quantities that the GPS receiver can directly measure. Observables do not directly tell you where the GPS receiver is situated but with using various techniques will allow you to calculate position solutions that do tell you where the GPS is situated. The two observables we consider are the code observable and the phase observable.

We have seen by tracking the maximum point of 7 using early late timing and a costas PLL that we have found r (t) = At (t) + AT _RX (t) - AT _TX (t) + /1000 for some fixed integer M and Φ it ) = t.Af + At (t) f _TX + Αφ ₀ + N/2 for some fixed integer N. These are the code observable and the phase observable respectively so far. However, there are some added complications and the form they take can differ. APPENDIX A, GPS BASICS

A.4.1 The code observable

Once is tracked the receiver has access to the navigation data. The navigation data is sent as 30 bits per word. There are 10 words in a subframe taking 6 seconds to transmit. Each subframe contains a, Hand Over Word (HOW) word that indicates the exact time when the leading-edge of the first bit of the navigation data was transmitted from the satellite. In the satellite the first bit of every NAV data transition is aligned to first chip of the the C/A code. This is possible as they are derived from the same oscillator (see A. l) . Figure A.4 shows the C/A NAV timing relationship.

Figure A.4: C/A NAV timing relationship

Because of this unique time stamp every 6 seconds and that the receiver is continuously tracking the C/A code of the satellite, each chip of a C/A code can be uniquely identified with an exact time of transmission. Therefore, the receiver can resolve the ambiguity M in the code observable and hence can estimate τ (t) such that r (t) = At (t) + ATRX (t) - AT _TX (t).

Usually the code observable is in units of meters and is called the pseudorange. Converting r (t) into meters by multiplying by the speed of light c results in the following pseudorange equation where the variable time has been removed for brevity and p is the range from the transmitter at transmission time to receive at reception time. APPENDIX A. GPS BASICS

p = + c (A¾ - ΑΤτχ) (A.16 )

A.4.1.1 Calculating code based solutions

The satellite's current clock hiasATrx is transmitted in the navigation data and therefore is a known value. The unknown values are therefore the receiver's position and clock bias: together these are P = [x, y, z, ΑΤρχ j . This means the pseudorange is a function of these unknown variables. Obtaining one such pseudoranges for a satellite results in the following nonlinear equation.

Because of the ease of solving linear equations, linearization of this equation using a first order Taylor expansion is sensible for deriving a generalized method of solving sets of pseudoranges of arbitrary sizes. We let an estimated solution be x Cj, z, AT' _f tx l for time of reception. Then, a first order Taylor expansion for p _n around P is as follows.

a a

Pn (P) « Pn (?) + {X - X) P; y) TPn

dx

_^ O \ d

[ Z - Z)—p _n (A.18) dAT R _f1 X _>

Given a satellite's position S _n — [x _n , y _n , z _n \ at time of transmission for a pseudorange p _n , then, the partial derivatives can be calculated given p _n {x, y, z) = \! (x— x _n + (y— y _n ) ² + {z— z _n y. Therefore, upon evaluation, equation A.18 can be written as follows.

{y - y + (AT _RX - T _RX ) c (A.19) APPENDIX A, GPS BASICS

Given a good estimate P, this equation has four independent unknowns, hence at least four equation just like it are needed to solve the unknowns. These four equations require the receiver's clock bias to be the same for all equations and the receiver's position to be the same for all equations. Therefore, the receiver has to obtain four pseudoranges simultaneously. A set of m such pseudoranges for different satellites obtained simultaneously are given below.

Pi (P) = pi + c (ATRX— Δ ^ p _m (P) = p _m + c {AT _R x - T _j x' ) (A.20)

Because equation A.19 is linear this set of pseudo ranges can be written in matrix form as follows.

(A.21)

Using bold tvpe notation vectors or matrices, this has the form of A ( P— P ) ss b where A, b. and P are known. Using LS and rearranging for the unknown yields P f» (A ¹ A] A ¹ b + P. The right-hand side of this equation when

from an old solution estimate. This process can be reiterated as the following algorithm describ es . APPENDIX A. GPS BASICS

Algorithm A. I Iterative LS solution using code observable

L Po = [0, 0, 0, 0] , n - 0

2. increment n

3. estimate reception time t _rx — TRX - TRX

'-ί, calculate S _m at time t _rx

-5. calculate S _m at time t _rx — \\ S _m — [x _n , ΰη, z _n \ jc and reitterate a few times

6. if [ ¹ A) stop wih error

P., - f Λ ' Λ ) ^' A ' b · P.,

8. if ( < less than desired error) stop

9. if ( n > to big ) stop with error

10. goto 2

The correction made to the time of reception as believed by the receiver Τ#χ- to produce an estimate of the reception time in step 3, for most receivers would be small in the order of less than a millisecond and could be conceivably ignored. The correction needed for calculating an estimate of the time of transmission as performed in step 5 is generally comparatively large, and is in the order of around 60 ms being the approximate flight time from the satellite to the receiver; this step can't be ignored.

Figure A.5 shows an example of algorithm A. l converging for a set of six satellites and code observations taken of them . As can be seen using the center of the Earth as the initial estimate within six iterations the algorithm has converged. APPENDIX A, GPS BASICS

2 4 6 8 10 12 14 16 18

Iteration

Figure A, 5: E Example of convergence of algorithm A. l

The output of this algorithm for the example in figure A, 5 resulted in a receiver clock bias estimate of—116.644 μ-s. This compares to an estimate as computed by the receiver itself of—116.738 μα with a spatial solution discrepancy of 32 ra between the two.

A.4.1.2 Final code observable model

One reason why the discrepancy of 32 ra between our solution and the solution as calculated by the GPS receiver itself is we have neglected some things in our modeling of equation A.16. The model in equation A.16 can be extended by adding tropospheric delays T, ionospheric delays 1 , multipath M _r and miscellaneous errors e _r . Therefore, a more exact model of the phase observable can be written as in the following equation.

p = p + c (AT _RX - AT _TX ) + T + I + M _r + e _r { A..∑ 2 j

Of the sources of error ionospheric are typically the greatest. The ionosphere stretches from 50 km to 1000 km above the Earth consisting mainly of charged particles, charged atoms and charged molecules. A large part of the ionization is caused by UltraViolet (UV) rays from the sun and hence there is a large APPENDIX A. GPS BASICS diurnal change in the Total Electron Count (TEC) which in turn effects the ionospheric correction term /. The ionosphere can produce a satellite range error as little as 1 rn to as much as 100 m | 13] , The ionospheric correction term is frequency dependent and with dual band receivers ionospheric free combinations of observables are possible. For singleband receivers no such combination is possible, instead Klobuchar ionospheric model is used for singleband receivers GPS. Klobuchar coefficients are transmitted in the navigation message so that the receiver can then estimate ionospheric correction terms. The Klobuchar algorithm corrects about 50% of the ionospheric errors [40 j.

The important thing for us in this section is that we are able to obtain an approximate solution, the spatial component with some sort of accuracy less than 100 m, and a time accuracy of some sort less than a 1 .,s. This is all we are concerned about regarding the code observable.

A.4.2 The phase observable

The phase observable is a measured quantity taken by the GPS receiver for a particular satellite for a particular time. Phase observables allows higher accuracy GPS measurements to be made than compared to that of ones solely using code observables. This is due to the much shorter wavelength of GPS carrier than compared to the chip length of the code observable. The wavelength of Li is approximately 20 cm compared to approximately 300 m length for the code chip of the C/A signal and can result in a correspondingly large increase in accuracy. The measurement comes from monitoring the phase difference between the received satellite carrier and a reference oscillator on the GPS receiver. The receiver accumulates this instantaneous phase difference by tracking and outputs this to the user as the phase observable. Figure A.6 shows a block diagram of what the GPS receiver is doing when observing a satellite for the phase observable neglecting all signals sent on the carrier such as C/A, navigation and P(Y) code. APPENDIX A, GPS BASICS satellite channel receiver

P ase observable

Figure A.6: Simplified block diagram of phase measurement

For GPS. phase is customary in units of cycles rather than radians or degrees for GPS work. Phase is the argument inside a trigonometric function that accepts units of cycles. Φ ^ώ is the phase of the carrier of the satellite while Φ ^υ is the receiver's reference phase. Both Φ ⁸ and Φ ^υ can become arbitrarily large.

Inherent in accumulation of phase is an ambiguity N that depends on when you started accumulating phase. In addition to the ambiguity there is the possibility of missing some rotations. Counting the number of times a car tire rotates depends on when you started counting its rotations and also depends on whether you missed any rotations. While GPS receivers try to continuously track the phase, this is not always possible. Due to noise, loss of signal or turning the GPS receiver off and on again, the tracking of the phase can be lost resulting in an integer change in the value of N. This produces what is called a "cycle slip". So, ideally this phase ambiguity should be fixed while the satellite is being tracked but due to cycle slips occasionally it will change.

As we have seen satellites don't send out continuous waves, they are modulated with two BPSK signals, one in the quadrature phase and the other in the in- phase. Code tracking has the effect of wiping the C/A code from one of the BPSK signals but still leaves the navigation adding a level of complexity when trying to track it. When the BPSK data is not used to regenerate the original carrier wave a half cycle ambiguity in the carrier phase is introduced into the APPENDIX A. GPS BASICS phase observable and it is said that the phase observable has a code factor of two Cf — 2. When the BPSK data is used to regenerate the original carrier wave, the original carrier wave can be fully regenerated with an ambiguity of one cycle and it is said that the phase observable has a code factor of two Cf — 1. Thus the ambiguity of the phase observable can be reduced when a code factor of one is used.

A.4.2.1 Final phase observable model

As we have seen tracking the carrier phase maximizes which in turn means the accumulated phase while tracking is Φ + At fox + Δ ο + N/2 for some fixed intesrer N when using a costas PLL. When considering ^' the code factor this can be written as follows.

Φ - tAf + At f _TX + Δ ₀ - N/Cf (A.23)

Assuming the receiver's clock is based around its local oscillator and its frequency keeps perfect time, true GPS time can be converted into the time as determined by the receiver as ¾— {tjRx + Φ§ ^Χ ) / fo- By definition true time plus clock bias is also the time as determined by the receiver TRX t + ATRX . Equating the two and rearranging yields the clock bias in term of true GPS time. Likewise this can be done for the satellite's clock.

Δ¾ - (t f _RX + φ ) /f _Q - t (A.24) TTX - (ifox + 4 ^Χ ) /fo - t (A.25) Subtracting the two and multiplying by ₀ results in the following.

fo (Δ¾ ^~ AT _TX ) - <Δ/ + Αφ ₀ ( A.26)

Equating this with equation A.23 we see that we can write equation A.23 as follows APPENDIX A, GPS BASICS

Φ = f ₀ { Tnx - AT _TX ) + -r-pfrx/fo + N/C, (A.27) frx is the oscillator of the GPS transmitter and is an atomic clock being extremely close to / ₀ . Therefore frx/fo = 1 for our purposes and we can rewrite equation A.27 as follows. = /ο (ATRX - ATTX) + ρ/λο + N/Cf (A.28)

This model in equation A.28 can be extended by adding tropospheric delays T, ionospheric delays J, multipath Μφ and miscellaneous errors βψ as was done with the code observable. However, the ionospheric correction for the phase while being of the same magnitude of that of the code observable is of the opposite sign. Therefore, a more exact model of the phase observable can be written as in the following equation.

Φ = fo (ΑΤηχ - ATTX) + (p + T - I + Μ _φ + β _φ ) /λ ₀ + N/C _f (A.29)

The range term p in equation A.29 is for the receiver at reception time ί _¾ and the satellite at transmission time of ίτχ- So the distance p is a measure of where you are to where the satellite was a short period of time ago because trx is an earlier time than the current time of ίχχ; the difference between these two values is typically around the 60 ras mark and a satellite can move a few? hundred meters in this time.

The phase observable was measured at GPS time tux, this variable itself has to be solved for, as you are not going to know exactly w ^T hat time the measurement was performed; you know you performed a measurement now, but you don't know when now is. This can be obtained using the code observable as previously shown.

The receiver's clock bias is at t x while the satellite's clock bias is at txx, but these aren't so critical as these don't change rapidly over 60 ms and can safely be assumed to be constant over the short time periods. APPENDIX A. GPS BASICS liquation A.29 is our final model for the phase observable. The left-hand side is what the receiver gives us, while the right-hand side is what we interpret it as. Multiplying it by the satellite's nominal wavelength is still classified as the single difference but rather than units of cycles for the phase observable, the phase observables units become meters.

APPENDIX A. GPS BASICS

A.5 Selected proofs

A.5.1 Received phase using flight time approximation.

Ta lor expantion of flight time

At \i) J r^ {t - i _m )

At (i) - At (t _m ) + {t- t _m ) + tmf Ι^ + '··

First order linear approximation

At (t) » At (t _m ) + (t- t _m ) ^ ¹

Define received phase

Φ (t) - tAf - At (t) frx + Αφο

→ Φ i t) « ίΔ/ (Δέ (t _m ) -I- (i - f, _H )∑¾si) _TX - Δ ο

→ Φ ~ t (Δ/ - + i _w ^/r - ^■ A* (i _m ) frx + ₀

→ Φ (t) « t (Δ/ - Af _TX ) + Θ (t _m ) □

A.5.2 T First-order linear approximation

At (*) « ^ ^{h itm)} it - tm) + constantx

AT _TX (t) « /rx(tm)- o ^ _ _im ) ₊ constmi

ATRX (?) « " ^x(f r ^{)" /o} ίί -- t _m ) + constants

-- r w (i - t _m ) (Δ/ - Δ _, χχ (t _m )) -jJ- + E (i _m ) as / _rx ¾ / ₀ APPENDIX A. GPS BASICS

A.5.3 Radial velocity with constant radial velocity offset e:

What happens with constant velocity offset. If radial velocity is out by e at all times then lim Φ it) = tAf - (At (t _m ) + (t - t _m ) JTX + + δ + Αφ ₀ lim Φ (t)

Where θ' (t _m ) and δ are some constants

A.5.4 Maximum radial velocity and acceleration of the satellite with respect to the receiver

Figure A.7 is a simplified model of satellite orbiting the ^' Earth while transmitting to a receiver. No relativistic effects are considered and it is assumes that the satellite's orbit is perfectly circular with constant tangential velocity and when the satellite is closest to the receiver the satellite is directly overhead.

APPENDIX A, GPS BASICS

Figure A.7: Simplified model of a satellite transmitting to a receiver while orbiting

From figure A.7 we obtain the following equations to describe the system.

2πί

x {t) r cos

p [t) = ix* it) + y ² [t]

The following approximate generally recognized values for GPS satellite orbital period T , radius of the earth and radius of GPS satellite orbits were used. De ^¬ differentiating the range equation A.30 figure A.8 was obtained. APPENDIX A. GPS BASICS

r _e = 6371000 m

r _s = 26600000 m (A.31) T - 43080 s

Modeled radial velocity and acceleration with respect receiver versus time after satellite rise

Figure A.8: Modeled raxiial velocity and acceleration

The maximum radial velocity derived from the model was 929 ms ^"1 at satellite rise time and set time. Maximum acceleration was 0.1178 ms ^{~ 2} and was when the satellite was directly overhead. [57] (pg 91) states that GPS satellites can have a radial velocity of up to 800 ms ^"1 with respect to a, stationary receiver on earth. This is consistent with our simplified model. On the LI band an acceleration of 0.1178 ms ^{" 2} is approximately a 0.9 Hzs ^{'" 1} doppler shift rate while 929 ms ^"1 on the LI band is approximately a 5 kHz doppler shift, therefore /o ± 5 kHz must be searched for the carrier frequency.

Using multiple phase observations and taking the difference of them produces various differences that can be used finding position solutions. In this appendix we look at these differences and what form the position solutions take.

A single difference is the difference between two phase observables at the same epoch for two different GPS receivers (two different users) . Performing a single difference eliminates systematic errors when observing one satellite. Such errors include the satellite clock bias and to a large extent atmospheric errors.

We now use two GPS receivers, A and B to observe satellite S. The phase observables are measured at approximately the same time by both receivers. User A reports a time of as reported by the system clock on the receiver itself, while user B reports that the clock on it's GPS receiver when taking it's measurement was 7¾. The relationship between actual GPS time and times according to the clocks on the receivers themselves are as follows.

(B. l APPENDIX B. DIFFERENCING THE PHASE OBSERVABLE

LB — TB - AT _B (B.2)

Using the pseudo-range observables it is possible to obtain the receivers clock biases AT A and Δ¾ . A good estimate of the clock bias of less than 1 μ,3 is assumed as given from now on.

Two observable models can be written down as follows where it is understood that the multipath errors and miscellaneous errors refer to the phase observable and not code.

λ | (t _B ) = i (t _B ) + c (AT _B - ATs) + T - I + f+4 + XNi/Cf ΑΦ¾ U A ) = p ^s _A (t _A ) + c (AT _A - ATs) + T - / - M + + XN¾/C _f Two receivers one satellite near simultaneous observations (B.3)

Where ATs is the satellite's clock bias, T the tropospheric correction, 1 the ionospheric correction, (7/ the code factor, iV and N integer ambiguities, Μ _β and multipath correction, ef and e _A miscellaneous errors, c speed of light, pg is) and 0% {IA ) receiver to satellite ranges, Φ (is ) and Φ¾ (ί _/Λ ) receiver phase measurements, and λ is the satellite's nominally wavelength. The subscript on Φ indicates which receiver has performed the measurement and the superscript indicates to which satellite, f (i) is the range from the satellite S at transmission time to the receiver A at reception time of i . N' indicates that the unknown integer N is for receiver A monitoring satellite S. M f is the phase multipath affecting the signal from satellite S to receiver A and e ^s are the miscellaneous errors affecting the phase observable that the receiver measures for satellite S that are not modeled. These two phase observations ( λ ' and Α ¾ ) are given in meters rather than cycles by multiplying by the satellite's nominal wavelength.

While the tropospheric and ionosphere corrections are not exactly the same for both receivers we make an assumption that they are. Also as measurements were not entirely done simultaneously, the satellite clock bias will be different when when the two receivers take measurements, but here we assume the satellite clock bias ATs is the same for both receivers. APPENDIX B. DIFFERENCING THE PHASE OBSERVABLE

As the observations were not simultaneous the two users see the satellite in a different place; the following figure describes the situation.

As seen from Earth

Figure B.l: One satellite as seen from two different receivers when the observations are not simultaneous

We therefore extrapolate to a common epoch on the second t. We now have to adjust the phase observables as if we had actually taken the measurement at this common epoch. To do this we can use a first order approximation using a Taylor expansion on the two range terms and rearrange which results in the following.

ρ' _β (tjg) is the range of the satellite calculated from where the satellite was at transmission time to where the receiver is at reception time given a reception time of tfj . This can be written as ¾ = jj S iiif) — Β (ίβ) || where S is the satellite's position, B is the receiver's position, t _B is reception time and ' _/ f ^::::: - Δ/; is the transmission time.

Substituting this back into the observable model in equation B.3 we get the following. APPENDIX B. DIFFERENCING THE PHASE OBSERVABLE

λΦ| (ts)+ - (t - is) « PB (t)+c ( T _B - &Ts)+T-I+M§+e ^s _B +XN§/Cf

B.5)

The range rate ^ - can be calculated either by using an approximate receiver's position and using navigation data, or, by measureing the Doppler. To find the range rate using the receiver's approximate position, differentiating ρ|. (i) using the chain rule produces the following (see B.4.1).

at

Range rate calculation

If we assume the receiver is more or less stationary and the range rate is much less than the speed of light then we obtain the following.

Approximate range rate calculation

The maximum range rate as we have already seen is less than 800 ms ^{" 1} which is much less than the speed of light so it's a valid assumption. Generally t— ts will be much less than 1 ms which would produce a correction error of less than 10 μπι in equation B.5 which is an acceptable error.

The terms on the left-hand side of equation B.5 are all obtainable. Φ is the measured phase observable at the common epoch t we have chosen, the actual time of the measurement ίβ is known and the range rate we also can estimate. The left-hand side is a corrected version of the phase observable to account for the fact that the phase observables were not taken at the same time. We therefore rewrite the given phase observables as follows. APPENDIX B. DIFFERENCING THE PHASE OBSERVABLE

λΦ| (t) « λΦ| it _B ) it

The hat on top of the symbol is to show that we have not measured the phase at this time, but instead it is an estimate of it at that time. The time difference between say t and t _A is likely to be less than a millisecond which is a relatively short period of time for what we're talking about here and hence is an acceptable approximation. Upon substituting equations B.8 into equations B.3 we get the following phase observables extrapolated to the common epoch t.

ΛΦ| (*) = pf (t) + c (AT _B - AT _S ) + Γ - / + | + e| + AJVf/C) λ¾ if) - pi (t ) + c (AT _A - AT _S ) + T - I + M _A ^S + e ^s _A + AN£/C> Two receivers one satellite extrapolated to a common epoch (B.9)

Taking the difference between the two observables in equation B.9 we get the following where N§ _A = Nf - NI, Φ _ΒΑ (ί) - H i t) - (ί), Δ | _Λ (t) - p| (t) - P¾ it) , Mp _)A = - f , m.cl = β| -

ΑΦΐ ₄ (*) = PBA (*) + c (Δ7¾ - A7k) + M _B ^S _A + 4A + N _B ^S _A /C _F Single difference model extrapolated to a common epoch (B.10)

This single difference has the effect of removing errors that are common to both receivers. These errors include atmospheric ones and ones that relate to the satellite itself. It assumes that the two receivers are close enough such that atmospheric conditions are identical for both users and the time between these two observations to be as small as possible. However it also doesn't account for multi-path and miscellaneous errors such as noise. APPENDIX B. DIFFERENCING THE PHASE OBSERVABLE

B.l.l Single difference interpretation

We now make the approximation that the satellite is much further away than the two receivers are to one another. This means the satellite appears almost at infinity to the receivers meaning that angle SAB is equal to the angle ABS. as is depicted in the following figure.

As seen from Earth

Figure B.2: Satellite almost appears at infinity for two receivers close to one another

GPS satellites are at least around 20,000 km away so for any baseline of a few tens of kilometers this is a valid assumption.

This means that there is a common vector pointing towards the satellite for both receivers. This we call the unit normal vector ή§. Below is a figure that depicts the current situation. We can see that there is a vector between the two receivers which we now call the position vector P = B— A and we also are able to visualize the range difference term in the right-hand side of equation APPENDIX B. DIFFERENCING THE PHASE OBSERVABLE

Satellite at infinity

Figure B.3: Unit norm in direction of satellite at infinity

Taking the dot product of the position vector (displacement vector) with the unit normal vector we can rewrite the range difference term as follows.

(*) = *s (t) · P (B.i l)

This allows us to rewrite equation B.10 as follows.

λΦΐ. ₄ (t) = ri _s (t) ^■ P + c (AT _B - AT _A ) + | _A + e% _A + XN§ _A /C _f Model of a single difference extrapolated to a common epoch (B.12)

The position vector is what we want to solve for. The position vector tells how user B's position is related to user A.

Using the code observable it is possible to obtain autonomous code based solutions which also consist of the receiver clock bias. However, these are not accurate enough to allow an improvement in position accuracy using Equation B.12 compared the autonomous code based solutions themselves. The Ubiox APPENDIX B. DIFFERENCING THE PHASE OBSERVABLE

LEA-6T which is a singleband receiver designed for precision timing applications claims "an accuracy of up to 15 ns" |15| which equates to theintroducing an error in equation B.12 of around 4.5 m. This gives us motivation for performing a double difference which eliminates systematic errors of the receivers themselves performing simultaneous observation

B.2 Double difference

Double differences are the differences of single differences. A double difference eliminates systematic errors of the receivers themselves when performing simultaneous observations. Such errors include receiver clock biases.

We take two single differences to two satellites S2 and SI still using receiver A and B. we also require that these single differences are taken simultaneously. This simultaneity is possible as GPS uses GDMA which means the GPS receiver sees every satellites in the sky simultaneously. Equation B.13 shows the pair of single differences.

λ ¾ it) - n _S2 (t) ^■ P + c (AT _B - AT _A ) + l + e& + XN§ _A /C _f

(t) = Asi (*) · P + c (AT _B - AT _A ) + M§A + e ^S BA + ANJi/C/ (B.13)

Taking the difference of the two we get the following where xis ₂ Si — ³ 2 ^~ ¾'i ;

][†S2S1 .... A, S2 Λ f Sl „S251 _¾ S2 ~S1 I TS2S1 ..... T TS2 ATSI „ J &S2S1 .... ^M BA ^"" BA ^{"" 1 W} BA ; ^e BA ' BA ^{"" e} BA ; ¹ BA. ~~ ^V S . ^~ ^ BA ^dna ' BA —

λ ¾ f ¹ (*) - Δη _525χ it) · P - Mi ¹ + ef + XN _A ^sl /C _f Model of double difference extrapolated to a common epoch (B.14)

Phase multipath can be as high as a quarter of a cycle for an observation and is not well correlated from one site to another |28] therefore the multipath component in the double difference can be as much as one cycle. The miscellaneous APPENDIX B. DIFFERENCING THE PHASE OBSERVABLE error component depends on many things but if we assume it is solely due to the ability to accurately measure the phase and we assume we can measure the phase to 1% of a cycle we would expect the miscellaneous error component to be 4% of the cycle. Indeed because of this multipath can be the dominant source of positioning error when using double differences.

B.2,1 Double difference interpretation

Ans ₂ siwe call the delta norm, while any vector that is dotted with a position vector is called a direction vector. The delta norm is theefore also a direction vector. We notice that the direction vector for the double difference is no longer a unit vector. The interpretation of this can be visualized with the following diagram.

Figure B.4: Visual interpretation of the Delta norm

For brevity let m = <Sg ' (t) C _f - C _f ( ¹ + ef ^] ) /Λ - JV ¹ and call it a corrected measurement, the double difference can be modeled as follows using our model in eauation B.14.

Simplified model of the double difference (B.15) APPENDIX B. DIFFERENCING THE PHASE OBSERVABLE

P is the unknown we wish to solve for so we treat it as a variable. As a variable many values of P exits that satisfy equation B.15, not just the one we're looking for of P— B— A as hence we can now no longer say P = B— A. We now look for what form the solutions to equation B.15 take thus finding the solutions to our double difference model. Equating equation B.15 when P = B— A— AB which we know satisfies equation B.15 and is the actual answer, with the general case when P satisfies equation B.21 but need not equal AB we obtain the following

Δ ^η ,5· ₂ ,5Ί · P = Δ ^η ₅ · ₂ 5ΐ · AB (B.16)

Upon rearranging we see

If we treat P and AB as point vectors then this has the form of a plane. The end of the AB vector is a point on this plane and the plane is normal to the delta norm. This can be visualized in the following figure.

Figure B.5: Double difference solution plane

This means the solution to equation B.14 forms a plane sitting on receiver B. We say that this plane is the solution plane of the double difference. APPENDIX B. DIFFERENCING THE PHASE OBSERVABLE

B.2,2 Double difference cycle slips

Ideally the integer JV¾ in the double difference in equation B.14 would be a constant, but as we have already mentioned in subsection A.4.2 on page 188 the GPS receiver can lose lock causing what is called a cycle slip. Figure B.6 shows an example of how these cycle slips effect double differences.

Double difference cycle slips

for SV14 and SV22

Time (s)

Figure B.6: Cycle slips in double differences

As can be seen both half cycle and full cycle ambiguities have been detected therefore the receiver (a Ublox LEA-6T) is not using the navigation data to reconstruct the carrier wave when this plot was produced. In this case from the plot is easy to see where the cycle slips have occurred and what is needed to correct for them but this is not always so especially when the receivers are widely separated. Various ways have been devised to detect and/or correct double difference cycle slips. Uncorrected and/or undetected cycle slips severely affect position solution accuracy [46] .

B.2,3 Double difference residual

A double difference residual (or simply a residual) is the difference between the double difference measured and the one expected given a position P; in APPENDIX B. DIFFERENCING THE PHASE OBSERVABLE cycles it is written as follows

Residuals are useful as a way of obtaining an idea of how correct a solution is or ho much error is in a, double difference.

We define the following notation for the positive fractional component of a double difference with respect to satellite SI and S2 and receivers A and B and call it the the wrapped double difference where {} is the sawtooth function.

Wrapped double difference (B.19)

Equating this with our double difference model we see the following.

? _/ ¾Sf ¹ (*) - { An _S2S1 (t) · P + Cj ( j ¹ + eg ) ,/λ}

Model of wrapped double difference (B.20)

This is our model for the wrapped double difference. The wrapped double difference removes the unknown integer N removing cycle slips.

B.3.1 Wrapped double difference interpretation

The wrapped double difference removes the unknown integer N by wrapping the double difference. When the wrapped double differences from one pair of satellites is plotted against time we get something that looks like the following figure. The data appears wrapped, hence its name. APPENDIX B. DIFFERENCING THE PHASE OBSERVABLE

Figure Β.7: Wrapped double difference data for one pair of satellites using two GPS receivers spaced a few centimeters apart.

Figure B.7 is for two receivers placed only a few centimeters apart. For receivers that are at distances greater than this the plot one gets when plotting wrapped double differences becomes visually meaningless and looks like nothing but a homogeneous mass of points.

For brevity let m _w - {C/<¾f (t) - C _j ( ¹ + egf ¹ ) /λ} and call it a corrected wrapped measurement then our model for the wrapped double difference can be written as follows.

Simplified model of the wrapped double difference

The form of the the solutions to our model of the wrapped double difference can be found in a similar way as was done for the double difference model in subsection B.2.1. We look for what form the solutions to equation B.21 take. Equating equation B.21 when P = B— A = AB which we know satisfies equation B.21 and as the actual answer, with the general case when P satisfies equation B.21 but need, not equal AB we obtain the following. APPENDIX B. DIFFERENCING THE PHASE OBSERVABLE

Taking the difference of the two.

251 i t )

. Λ λ

This is similar to the form the solutions take for the double differences except this time due to the sawtooth function there are an infinite number of planes. To remove the sawtooth function we note that there exists a k £ Z such that

Which means

This can be rewritten as follows where us2Si ^~ A/ (Cf il) ·

Wrapped double difference solution form (B.26 )

Thus the solutions form an infinite set of planes. One plane is the same as the solution plane for the double difference and the others are parallel to this APPENDIX B. DIFFERENCING THE PHASE OBSERVABLE plane where each plane has a fixed spacing of fis2Si meters between each other. Therefore we call ^.5251 the plane spaceing. The figure below visualizes three of the set of infinite solution planes.

Figure B.8: Three of the infinite set of wrapped double difference solution planes

B.3,2 Wrapped phase observable

We defined the wrapped phase observable as follows.

C _f ≠ _B (t _B ) = {θ _} Φ _Β (t _B ) \

Wrapped phase observable (B.27

This can be extrapolated to common epoch t as follows using equation B.8.

APPENDIX B. DIFFERENCING THE PHASE OBSERVABLE

^Γ ·' ^φ "' ''· ^· '/ ^· T ^' 7/Γ (t - t _B )

|6' ¾ (ij I ^* by equation B.8

Cf(j) _B (t) by definition

(B.28)

A single difference can then be formed given two receivers A and B that we call a wrapped single difference φ _ΒΑ (/;).

(B.29)

And finally given two wrapped single differences to two different satellites Si and S2 a wrapped double difference can be formed.

(i)}

C }fiΨjΒΑ \ t)

(B.30)

Thus it is possible to form wrapped double differences extrapolated to a com- APPENDIX B. DIFFERENCING THE PHASE OBSERVABLE mon epoch given wrapped phase observables. The advantages of doing it this way is that fewer bits are needed to represent a wrapped phase observable compared to a normal phase observable for a given accuracy; this saves network bandwidth.

Due to the properties of the saw tooth function, when creating a wrapped double difference from four wrapped observations, only one sawtooth operation is needed: this can be written as follows.

Cfi i ¹ ) - { {Cfi it) - c _f i (t)) - (CfH ¹ (t) - c _f (*)) } (B.3i)

Unwrapping is a method of re-creating the double difference from the wrapped double difference.

First off we define the sawtooth difference function as follows.

/ (x) = {x + 0.5} - 0.5 (B.32)

Here we define an algorithm that does unwrapping and is used in the CFF algorithm, G |0. f | and determines how sensitive the algorithm is to sudden changes in input.

APPENDIX B. DIFFERENCING THE PHASE OBSERVABLE

Algorithm B.l Unwrapping algorithm

1. Let /· U

2. ! : ./ ^' a¾)

3. ar _fc = (h) ^~ ¾ ^£ Δ ^η ₅₂ 5 _ΐ ( ) · P

4. Increment ;

5. yk— f {xk - <¾-i ) + jfc_i

6. Let i¾ = «¾¾ ·-}-- ( 1 --- a) _¾ „i

7. If not end of wrapped double differences goto 4

8- C/df ² / ¹ (t _k , P ) = y _fe -t- Ali52Sl · P

This unwrapping algorithm needs as an input an approximate position estimate P. This is because the phase of the double difference for a pair of satellites can change very rapidly with respect to time, using a position estimate in effect slows the wrapping effect so as to increase the likelihood that the unwrapping algorithm will be successful and not introduce cycle slips. The figure below shows an example of this algorithm on one pair of satellites with a position estimate half a meter different from reality for the residuals before and after the unwrapping has been performed.

Epoch

Figure B.9: Unwrapping algorithm given a position estimate half a meter out from realitv APPENDIX B. DIFFERENCING THE PHASE OBSERVABLE

The unwrapping algorithm can make mistakes by introducing cycle slips, particularly when the data is intermittent as is often the case for our scenario. However, the unwrapping algorithm, as we shall see, performs surprisingly well for static scenarios such as ours.

B.3.4 Wrapped double difference residual

A wrapped double difference residual (or simply a wrapped residual) is the sawtooth difference between the wrapped or unwrapped double difference measured and the one expected given a position P; in cycles it is written as follows.

Wrapped residuals in contrast to double difference residuals are cycle slip invariant. Wrapped residuals like double difference residuals express the error between what is measured in what is expected.

APPENDIX B. DIFFERENCING THE PHASE OBSERVABLE

Given p| (i _B ) S (if) - B (t _B )\\ and i - t _B - p (t _B ) /c then derive for 2D.

V (¾ - ¾ (ίβ)) ² + (Sy (if) -- ¾ (t _B ))

d

X ({S _x i ) ~ ¾ (¾)) ² + (S _y if) - ¾ ( _S )) ² di

1 d

(,¾ (if) - B _X (i _s )) ² + (S _y (if) - J3 _y (is))'

, A ,

— !— i (. _r ί ^' /y;) -B _t (t _R )) dt S _x [t%) -B _x {t _B )

+ (Sy (i ) By (t _B )) (Sy (i ) - By (t _B )) ) APPENDIX B. DIFFERENCING THE PHASE OBSERVABLE t(S, (i5) - ¾.(ie)) , (¾ (t'i) - B _y (t _B ))

PB (i ^j )

S(i¾ ^' ) Ui _: ) d(S(t%) -- B(i _B ))

PB tB) di

S(i¾ ^' ) -B(t _s ) (dSit ) dB(t _B )

PB M di dt j

B (i _B ) _ ( J+tX

B ! 2' i xlx

5 \ ^L B , i )

PB ^B)

_Ϊ 1 | (is)

^L B /

PB 'BJ

APPENDIX B. DIFFERENCING THE PHASE OBSERVABLE

B (t[

(s' .$'(>,.· ⁾ )

P ^S B (t B!

(t ⁱ B ^x J (t B\ 1

/ ¹ + s ^' ! ^/ y;

PB

(S(i¾f) rfS d

dt

ri?

dS

P' (t _B ) + ^■ (S (if) D the extension is obvious and results in same formula.

K ^L B) (t _B )) f dB

\ at dt

dt

c dt B ) UB),

As a proof of concept we designed and implemented an energy distribution algorithm based on the eircadian rhythm in addition to using an Assisted Global Positioning System (AGPS) system that injected timing, location and ephemerides into the GPS receivers. A 33 F capacitor was used for this test.

To spread the observations obtained reasonably uniformly at night using such an extremely small energy storage devices as a 33 F capacitor, an extremely low duty cycle with an extremely low on time was needed (Chapter 5 of [44] ); this motivated us to design a AGPS system. The AGPS system supplied the ephemerides, the almanac, and an approximate time and position to the GPS module via an alternative link rather than the normal SObaud communication channel from the satellites. This allowed the GPS receivers to save time not having to listen to the satellites just to obtain the ephemerides and almanac data.

Our implementation of assisted GPS comprised of obtaining the ephemerides, the almanac, and approximate time from the Assistnow server provided by Ublox. This data was saved on a local PC and synchronized periodically with the Assistnow server. This PC continuously sent out this data at a rate of approximately 2-50kb/s on a dedicated channel with a distinct Personal Area Network (PAN) IDentification (ID) (the AGPS channel) using an 802.15.4 APPENDIX C. NODE FUNCTIONALITY AT NIGHT radio. Another channel with a different PAN ID (the main channel) was used for retrieving the data from the GPS node and sending it back to the Personal Computer (PC). Figure C.l shows this communication implementation.

Figure C. l: AGPS network structure

Acknowledgments were required for the main channel due to cross talk from the AGPS channel to avoid packet losses on the main channel.

The amount of assisted GPS data at any particular time is not constant but is around the 3 kB amount. The GPS receiver could be initialized with GPS data within around 300 ms this compares to around 15 minutes to download this data directly from satellites.

We performed 10 trials with a clear view of the sky with and without our implementation of assisted GPS data without backup power being supplied to the GPS receiver. The TTFF was recorded for each trial and the results plotted in the figure below?. APPENDIX C. NODE FUNCTIONALITY AT NIGHT

Time to first fix (s)

Without AGPS without GPS backu power.

55

50 ■3- -□- -a-

45 □

ω 40

LL LL 35

h-

1—

30

□

25

20

15

6 10 12 Trial

Figure C.2: Time to first fix without GPS backup power

As can be seen the TTFF when using the AGPS system is considerably lower than when not using it. The averasre time when not usinsr AGPS was 42 seconds compared to 18 seconds when using the AGPS system. This TTFF test was then performed again supplying the GPS receiver with backup power using approximately 30 according measurements; a period of 60 seconds was used between each trial. The results are pictured in the following figure.

APPENDIX C. NODE FUNCTIONALITY AT NIGHT

Time to first fix with GPS backup power.

o With AG PS

□ With out AG PS

10

Figure C.3: Time to first fix with GPS backup power

The reason the first GPS trial TTFF is much greater than the remaining trials is due to the fact that the GPS backup power had not been given power until that trial. Excluding the first trial our AGPS implementation produced a TTFF of 4 seconds compared to 10 seconds when not using it.

Both with and without GPS backup power TTFF is approximately half in both cases.

The period of time between trials is important when using GPS backup power as whatever information the GPS receiver holds when the backup power is applied will age; given a, long enough period of time whatever information it holds will become useless. For example ephemerides are generally not used beyond four hours of age. Therefore a longer interval of 30 minutes was performed when using the AGPS and supplying GPS backup power. Seven sequential trials of this type was taken and the results are plotted in the figure below. APPENDIX C. NODE FUNCTIONALITY AT NIGHT

Time to first fix for AGPS

Trial

Figure C.4: Time to first fix for AGPS with GPS backup power and 30 minutes between trials

This test, excluding the first trial, the average TTFF was 5 seconds. This shows that periods of node inactivity of up to 30 minutes affects the TTFF when using AGPS and GPS backup power minimally.

The duty cycle at night can be calculated as follows

As can be seen the duty cycle is a function of night length. Therefore assuming a night length of 12 hours we see the duty cycle approximately 0.5% which needless to say is an extremely low duty cycle.

The wasted percentage of on time is given as follows where t _on is the on time for each active part of the duty cycle.

Hnit

(C.2) τ on

While the number of active periods throughout the night is given as follows. APPENDIX C. NODE FUNCTIONALITY AT NIGHT

Therefore we take a wasted percentage of on time of 50% which implies, given a TTFF of 5 seconds using our AGPS, an on time period of 10 seconds. This in turn implies approximately 20 active periods throughout the night. Each active period lasts 10 seconds and half of this energy used for each of these active periods is wasted due to initialization. While the duty cycle is 0.5%, as half of this percentage is wasted due to initialization we expect only 0.25% of the night to obtain observations; this we call effective duty cycle.

To perform measurements at night energy stored on a capacitor has to be distributed throughout the night to perform this task. In addition to this there must be enough energy reserves throughout the night to run other such things like the microprocessor. To allow this to happen we designed and implemented an energy distribution algorithm based on the idea of the circadian rhythm. The following describes the algorithm.

The first assumption is that it is always possible to fully charge the capacitor during one day of daylight hours. If not, then we say the capacitor is too big for the solar panel or the solar panel was not big enough. This assumption has to be checked for every solar panel and capacitor combination.

The second assumption is that all periods of time used to obtain observations take the same period of time and energy. Each one of these periods of times is called an event.

There are two types of events triggered and spaced.

Triggered events can happen at any time. They happen when the voltage of the capacitor exceeds the fixed voltage of 2.2 V.

The other type of event is the spaced event. These events happen evenly spaced between the last triggered event and the first trigger event as pictured in the following figure. APPENDIX C. NODE FUNCTIONALITY AT NIGHT

First triggered event

T T T

3

3 Not to scale

Figure C.5: Circadian distribution of stored energy

The sum of the energy of spaced events between the last and the first triggered events during the spaced event phase is roughly equal to the usable energy of the capacitor.

The first spaced event happens a/ (n + 1) after the last triggered event and happens approximately every a/ (n + 1) until the first trigger event, n is how many events a fully charged capacitor can do without going flat and without getting any new energy input, while a is the period of time between the last triggered event and the first triggered event during the spaced event phase. APPENDIX C. NODE FUNCTIONALITY AT NIGHT

When the node starts doing spaced events it will not do triggered events until at least a third of the spaced events have occurred. Under voltage is monitored during the spaced events phase and if encountered will stop doing spaced events. An event will be extended if the voltage is still above a fixed voltage say 2.2 V. The spacing of e vents is updated every time a spaced event happens. The number of spaced events is calculated by taking into consideration the current capacitor voltage after an event, the amount of energy each spaced event uses and how long away the next triggered event could possibly be. The node keeps a Real Time Calendar Clock (RTCC) active at all times. The RTCC is not set any particular time in the world it only cycles around every 24 hours. It is not necessary for it to know what the time is but it is necessary for it to have an idea on how long a rotation of the earth takes around its own axis. It performs the same function as our circadian rhythm.

The circadian energy distribution algorithm along with the AGPS system as described was then implemented on hardware and was tested. The following describes the results obtained from such a system.

At 5:00 PM a node was placed outside for approximately 18 hours and it's super capacitor voltage was recorded; figure C.6 shows the results of these measurements; as can be seen this graph looks similar to the one in figure C.5.

APPENDIX C. NODE FUNCTIONALITY AT NIGHT

Capacitor voltage night test.

4 8 12 16

Time Cap measured (hours)

¾ure C.6: Voltage over super capacitor at night

The short sharp drops throughout the night that can be seen in figure C.6 are caused by the GPS going active and taking measurements. The reason that there are only approximately 14 measurements throughout the night rather than closer to the 20 measurements as expected throughout the entire period between the last triggered event and the first triggered event is because the nodes circadian rhythm is not synchronized until a full 24 hours has past.

Performing a test lasting approximately 75 hours the autonomous code based position solutions of GPS module were recorded, the mean of this set of data was calculated and then for each position solution the offset from this mean was calculated and plotted in figure C.7.

APPENDIX C. NODE FUNCTIONALITY AT NIGHT

Position Offset from mean Versus time

10 20 30 40 50 60 70 80

Time (h) onomous code based solution offset from the mean code based ion with respect to time

The lower density places in figure C.7 are the position solutions obtained at night, the heavy density places are the daytime solutions. For this test one standard deviation equated 2.9 m. Approximately 7.7 position solutions were calculated for each spaced event. Nighttime approximately lasted 13 hours at the particular time of year at that location meaning a solution density at night of 0.003 position solutions a second, or an effective duty cycle of 0.3%; this compares very well with the approximate expected effective duty cycle of

During the day, to ail intents and purposes we can say the duty cycle was 100% by comparison, therefore 99.7% of the solutions were obtained during the day. However, the nighttime data was relatively evenly spaced throughout the night which was the goal of the energy distribution. Therefore our energy distribution algorithm is effective at operating the GPS receiver at night in such a way as not to cluster the observations it gets to greatly. A possible advantage of obtaining observations at night might be as a way of mitigating solution biases. In addition the GPS receiver is able to operate at times at night and at least obtain autonomous code based solution for those times of night increasing temporal information. APPENDIX C. NODE FUNCTIONALITY AT NIGHT

Due to time constraint, added complexity and that we where primarily interested in static solutions where many observations usually from a 24 hour period form one solution, distributing observations throughout the night was not taken further. However we have shown that it is indeed possible to obtain solutions at night using extremely low duty cycles and or very large capacitors. Clearly many other new algorithms could be created with greater complexity to allow such techniques as sidereal filtering and obtaining epoch based solutions.

APPENDIX C. NODE FUNCTIONALITY AT NIGHT

*

A i x

Data has to transfer between the nodes and the main computer that performs final processing. We divided this into two sections; the WSN section and the Internet section. The Internet section consists of the data flow from when it arrives at base station after any processing has been done such as time ambiguity resolution, duplicate packet removal, non-sequential packet removal, etc. The WSN section consists of the flow of data, and processing between the processed data of the nodes untill it is handed over to the Internet section.

D.l Internet section

The Internet structure of the network is a castor, server, client structure. This is the same idea that Networked Transport of RTCM via Internet Protocol (NTRIP) uses for accessing GPS receivers over the Internet. Out own implementation of a custom caster, server, client structure was implemented that we called JAC packet. The following figure shows the structure for two testbeds. APPENDIX D. NETWORK

Tesibed 1 base station

Figure D. l i Caster server client structure over the Internet

The device server acts as a driver for both the WSN gateway devices and physically connected GPS receivers as well as a supplier of data to the JAC packet caster. As such some of its tasks axe for the internet section while others are for the WSN section. In the previous figure the dotted line on the device server represents the junction between its two tasks. The device server acts as a server that the JAC packet caster can connect to to obtain data wanting to be logged at the JAC packet server. The JAC packet caster than attaches to the JAC packet server via a TCP connection that is kept open and the JAC packet server logs the incoming data from each test bed. JAC packet clients can then connect to the JAC packet server and either obtain the logged data from the test a bed or a live stream from one of the test beds. For a JAC packet caster to connect to a JAC packet server authentication is via a username and password; no authentication requirement was implemented for JAC packet clients. No logging of data was performed on the base station APPENDIX D. NETWORK meaning if the Internet connectio went down data would be lost.

During the design phase initially a single hop network was envisaged as the process of obtaining solutions was more important than creating difficult mul- tihop networks. However, in a real life testbed in LuShan Taiwan the node transmission range were found not to be sufficient and a way to extend the range was needed. Due to both limited time and limited resources on the nodes, a dedicated repeater node was designed and implemented to extend the range.

The following picture shows the nodes, repeaters and the part of the base station that is dedicated to creating the data that is handed over to the Internet section.

WSN

Testbed base station

Figure D.2: WSN and interface

D.2.1 Repeater

To use the repeaters two extra bytes are added to each data packet. One byte is a sequence number while the other byte is 8 flags where each flag corresponds to the modulo 8 of the repeater node number. The data packet sequence number is not the same as the sequence number in the header of each 802.15.4 packet. APPENDIX D. NETWORK

Every time a packet is transmitted by the node the data packet sequence number will increase. At the same time the flag byte transmitted will be set to zero. Therefore anyone listening to dat packets will be able to tell the order of such packets, have any of them been missed, in addition to knowing that it comes from a node and not a repeater as the flag byte has been set to zero.

A repeater contains a small list of the packets it has recently sent by storing the 802.15.4 header packet sequence number and source node number. When a repeater hears a packet if the repeater finds that it has already sent the packet it will not send it again: this helps prevent infinite transmission loops from occurring and repeating data the repeater has already sent. In addition when a repeater hears a packet, if the repeater finds that its bit of the repeater flag byte is already set it will not send this packet again; while maybe not strictly necessary this totally prevents infinite transmission loops from occurring.

When a repeater sends a packet it spoofs the 802.15.4 header so as to match the originating node. This spoofed information means both the source MAC address and the sequence number in the 802.15.4 header match that of the originating node. The repeater sets the bit of the repeater flag byte that corresponds to the repeater's node number. In addition it updates the list of transmitted packets using node number and the 802.15.4 header packet sequence number. The reason for using the 802,15.4 header packet sequence number rather than the dat packet sequence number is to allow the nodes the opportunity to send duplicate packets to reduce packet loss.

We are aware that the possibility still remains that two repeaters transmit simultaneously due to hearing the same packet from a node. However, if the repeaters can hear one another then the csma/ca algorithm should at least try to mitigate this effect so that the two resulting packets from the repeaters are sent sequentially,

D.2.2 Device server

The other half of the device server's role is to resolve time ambiguity from the "node algorithm data" that it hears, remove duplicate packets, and remove non sequential packets. APPENDIX D. NETWORK

First the the device server resolves the time ambiguity. To remove duplicate packets and non-sequential packets the device server does the following. It keeps track of the last 10 unique data packets received from each node in a cycle buffer. If the packet is found in the cycle buffer it is taken as a duplicate and dropped. If it's not a duplicate, the packet is put into the cycle buffer. If the GPS time of the packet is older than any of the other 10 packets in the cycle buffer the packet is dropped. Else, it is passed to the Internet section. This prevents both duplicate packets and non-sequential order of the packets being passed to the Internet section.

The repeaters simply repeats packets they hear; they have no knowledge as to when the packet will arrive. Because of this for the repeaters to function as repeaters, their radios have to be continuously on, which uses approximately 70% of the power that having the GPS receiver on uses. This is considerable power meaning duty cycling is required.

While there may not be enough irradiance to keep any one node active continuously, it is possible for there to always be one node active during daylight hours if the nodes all take a turn being active thus forcing the repeater to be on continuously. Of course if there is not enough irradiance, one repeater cannot be in an active state continuously. If the nodes and the repeaters were all somehow synchronized to become active at the same time and both nodes and repeaters were obtaining approximately the same amount of irradiance, the nodes would then all become inactive before the repeater had to, as the nodes use somewhat more power than the repeater; thus the repeater would be fully serving the nodes needs.

Another problem arises when a repeater is active and all nodes are inactive. In such a case when a node finally does become active the repeater may have to become inactive to save energy unable to function as a repeater for the active node. Therefore it makes sense for a repeater to become inactive when node activity ceases.

Prom these two problems the idea of phasing all nodes and repeaters in the WSN network was formed. The idea was to make the repeaters and nodes APPENDIX D. NETWORK become active at the same time and become inactive at the same time; such as to be in phase with one another.

As double differences can only be formed from observations from two GPS receivers that are on at the same time, phasing the nodes to become active at the same time means more double differences can be formed from nodes themselves than if the nodes were not in phase with one another. While an interesting observation, due to time constraints this effect was not investigated.

D.2.3.1 Phasing algorithm

We designed and implemented a phasing algorithm on our nodes and repeaters. The main assumption of this algorithm is that all nodes and repeaters obtain the same irradiance.

With the assumption that the irradiance is the same for all nodes and repeaters, all nodes can have the same duty cycle and on time, at the same time as the repeaters use slightly less power when active than the nodes, the repeaters can at least match the duty cycle and on time of the nodes.

Whether or not a node is obtaining useful observations it will send data packets to inform both nodes and repeaters that the node has its GPS receiver on; these will be empty when no useful observations are obtainable else they will contain the usual "node algorithm data ^" . If they are empty they are dropped at the device server.

When a node's super capacitor reaches 2.2 V it is considered to be full and the node will start a phasing period that times out after 10 minutes if no GPS packet are heard. The phasing algorithm as run on the node is pictured in the following figure. APPENDIX D. NETWORK

Figure D.3: Phasing algorithm on node

In our implementation during the node phasing a comparatively small power is required as the node only turns on its radio for 2 seconds and then off again for 18 seconds using approximately 10 mW or 4% of what the GPS receiver would use if switched on. In our implementation this phasing period may last up to 10 minutes. During this phasing period if any data packets are heard that signify other nodes have their GPS receivers on, the node will exit the phasing period and commence usual gathering of observational GPS data. If APPENDIX D. NETWORK on the other hand, after 10 minutes still no data packets are heard that signify other nodes have their GPS receivers on, the node will exit the phasing period and commence usual gathering of observational GPS data.

The phasing algorithm as run on the repeaters is depicted in the following figure.

APPENDIX D. NETWORK

Figure D.4: Phasing algorithm on repeater

During times of no node activity the repeater runs a phasing period, the repeater will turn its radio on for 2 seconds to listen for data packets and then turn off again. If data packets are heard that signify other nodes have their GPS receivers on, the repeater then becomes active and a timeout commences APPENDIX D. NETWORK that will get triggered if no new data packets are heard that signify other nodes have their GPS receivers on. When the timeout gets triggered the repeater goes back into the phasing period once again. Wherever the repeater is in its algorithm, if voltage level of the capacitor becomes too low, the radio will get switched off and the repeater will wait for its super capacitor voltage to rise once again to 2.2 V.

An epoch based solution is one where a solution is obtained using observations from one epoch at the time. A daily based solution is one where a solution is obtained using all available observations from one day at a time. We primarily focus on daily static based solutions in this thesis.

The sidereal part of the algorithm implements sidereal filtering using a stacking technique of both fix solutions and absolute wrapped residuals stacking with respect to a biased location; this is designed to improve solution accuracy and precision for our particular scenario where observations are intermittent. In this chapter we start by investigating multipath which turns out to be the largest source of error in our double differences. From this, assuming ambiguity resolution is a solved problem, we rev/rite our double difference model to include only multipath and additive white Gaussian noise and determine the distribution of our solutions we obtain when intermittently taking observations. We then derive our sidereal filtering formula and show that both in simulations and in real life it does indeed improve our solutions. We also show that due to the random intermittent nature that we obtain observations, we obtain a greater improvement in solution precision than if we were using continuously powered GPS receivers output ting constant high rate observations.

The multipath error of the phase observable is due to the superposition of the carrier wave reaching the receiver via different paths. This has the effect of appearing to delay the carrier wave reaching the receiver. For a single reflected signal in addition to the direct line of sight signal the multipath error in the CHAPTER, 3. MULTIPATH AND SIDEREAL FILTERING phase observable can be modeled as follows where a is a damping factor, λ is the nominal wavelength, and φ is the phase difference in cycles between the two signals [60] .

a sin (2πό)

^J 6B arctan

This produces an absolute maximum error of e _B — λ/4 when the damping factor is equal to one. Therefore the theoretical maximum double difference error can be one cycle or approximately 19 cm as four such phase observables are used to produce a double difference. Double differencing removes most systematic errors but leaves errors that are related to the geometry of the observed satellites such as multipath and PCV. Figure 3.1 shows double difference residuals for all pairs of satellites obtained using our hardware at the LuShan testbed site in Taiwan for one day.

All double difference residuals

Time (Hours)

Figure 3.1: All double difference residuals for one day of observations at the LuShan site for one pair of receivers

PCV appears very similar to multipath except for a dependence on antenna rotation. Because of this we lump PCV together with multipath errors. As can be seen from Figure 3.1 there are deviations of up to about one cycle which is the maximum value we would expect from the multipath error. However, the unwrapping algorithm can fail and introduce cycle slips which may also be a cause of some of the extremely high readings obtained in the previous figure, CHAPTER 3. MUUTIPATH AND SIDEREAL FILTERING

Figure 3.1. In processing our data we used a code factor of two so cycle slips in Figure 3.1 would results in an offset of 0.5 of a cycle. We shall see that multipath turns out to be the most major component of the variance in double difference residuals. This also supported by references such as |55| claiming that "multipath is the dominant error source in differential carrier phase and pseudorange measurements" and [45] saying "multipath is often considered the most limiting factor in precise GPS positioning". The double differences from Figure 3.1 are not the ones used to obtain the solution. Instead, the independent set of double differences using the satellite with the highest elevation as the reference satellite is used; the residuals of these double differences can be seen in Figure 3.2. independent set of double difference residuals used for solution using maximum elevation satellite

Time (Hours)

Figure 3.2: Independent set of double difference residuals for one day of observations at the LuShan site for one pair of receivers using the highest elevation satellite as the reference satellite

From Figure 3.2 it is clear that the variance is considerably less but as we shall see is still a great deal larger than the average measurement noise.

From the residuals in Figure 3.1 there appears to a large component that varies slowly over a relatively long time of a few hours. In addition to a slow component, a fast component that appears periodic and has a period of typically around 5 to 10 minutes can also be seen and has somewhat of triangular nature. The fast periodic component can be seen more clearly if we CHAPTER, 3. MULTIPATH AND SIDEREAL FILTERING magnify just a small amount of one of the traces from the figure; this can be seen in Figure 3.3.

694 696 698 700 702 704 706 708 710 712 714

Time (Minutes)

Figure 3.3: Fast periodic component of residuals

This highly time correlated multipath is called specular multipath and results from reflections of smooth surfaces (presumably the ground) and leads to positioning errors both when solutions are calculated every epoch and also as in our case when one solution is calculated for an entire day [41, 20]. In high multipath environments solution errors can be substantial and be in the order of a few centimeters [20] .

Due to the GPS satellite orbital period, the satellites repeat their locations approximately every 24 hours -246 seconds with respect to a ground based user, also erroneously referred to as the GPS sidereal period (erroneously, as a sidereal day is not 24 hours -246 seconds but instead 24 hours -236 seconds) . While the satellites have an orbit period of slightly less than 12 hours, due to the rotation of the earth the satellites appear approximately in the same place in the sky with respect to an earth based observer approximately every 24 hours rather than 12 hours. The GPS satellites perform manoeuvres occasionally to ensure they approximately maintain fixed ground tracks with respect to the earth. According to literature the earth oblateness causes a westward drift in the ascending node which is compensated for by decreasing the GPS orbital period by approximately four seconds [19]. There are two orbital periods in CHAPTER 3. MUUTIPATH AND SIDEREAL FILTERING a day which causes approximately eight seconds difference from the sidereal day. This value is only approximate and the satellite repeat times have been seen to vary by at least 9 seconds [24] . As the satellite repeat period is only approximate we simply use one fixed repeat time of 24 hours -246 seconds which is 1 second difference from the mean repeat time found by for a set of observations by [24] and matches the mean repeat period given by [34] in a set of data. In our discussion we refer to the GPS satellite repeat time of 24 hours -246 seconds erroneously as the sidereal period and take it as given that it refers to our given nominal satellite repeat period and not the actual sidereal day length, we do this as this filtering process is commonly called sidereal filtering not GPS orbital repeat period filtering.

Multipath is a function of the GPS receiver's surroundings and satellite geometry. Therefore, as long as the GPS receiver's surroundings do not change with respect to the GPS receiver then the multipath error should be similar from one sidereal day to the next. This allows a way of determining whether or not the large errors we see in Figure 3.1 are due to multipath or some other form of error. In the two Figures 3.4 and 3.5 that follow we plot the double difference residuals for one pair of satellites for two distinct sidereal days for fixed GPS receivers obtained from our hardware; one for the whole time the satellite pair were visible in a day and the other one for a small close up portion of them.

Time (Hours)

Figure 3.4: Double difference residuals for one SV pair over two different days: full view CHAPTER, 3. MULTIPATH AND SIDEREAL FILTERING days close up

9.55 9.6 9.65 9.7 9.75 9.8 9.85 9.9 9.95

Time (Hours)

Figure 3.5: Double difference residuals for one SV pair over two different days; close up view

As can be seen there is an extraordinary close match between days even for the small fast periodic triangular components. Therefore we conclude that most of the errors visible in Figure 3.1 are due to multipath and PCV.

Assuming integer ambiguity resolution has been successful, we can rewrite our double difference model as follows more simply where all subscripts, superscripts and indexing have been removed.

ΛΨ An · P + M + e (3.1)

We have divided the remaining errors into M consisting of errors that are due to geometry (Satellite location, receiver surroundings, and PCV) which we treat as multipath, and miscellaneous errors e. The miscellaneous errors will consist of measurement noise and other errors. As we only use the independent set of double differences where the satellite with the maximum elevation is used as the reference satellite for a solution, we are only concerned with the remaining errors in only these independent double differences not the others. The double CHAPTER 3. MUUTIPATH AND SIDEREAL FILTERING differences we're interested other ones as seen in Figure 3.2. We only use these double differences when estimating the remaining errors.

3.2.1 Measurement noise

To obtain an approximate value of the measurement noise, we took the one second time differenced double difference residuals (r (t)—r (t— 1)) for 60 pairs of satellites given an accurate position estimate, then measured the average variance. As the satellites appear to move relatively slowly from and earth- bound observer and assuming other errors change relatively slowly, the one second difference will consist mainly of measurement noise. Figure 3.6 below shows the one second time differenced double difference residuals for one of the pairs of satellites used to obtain the estimate. ware

Time (Hours)

^" Figure 3.6: One second time differenced double difference residuals for one pair of satellites using our receivers

As this is a difference of additive white gaussian noise, for a, double difference this needs to be divided by the square root of two. From this we obtain the average standard deviation of the double difference due to measurement noise of 0.013 cycles or a variance of var (noise)— 0.000169 cycles ² . As measurement error depends on signal strength, an error as low as 0.009 cycles was obtained for one pair of satellites while a maximum of 0.019 cycles was obtained for CHAPTER, 3. MULTIPATH AND SIDEREAL FILTERING another pair of satellites. As a double difference consists of two doublings of error, to obtain an estimate of the measurement error of one of our GPS receivers these values have to be divided by the square root of four. Therefore, on average our GPS receivers should have a 95% confidence that it will correctly measure the phase observable to better than 0.013 cycles of the true value or about 1% a cycle.

3,2.2 Multipath error

The multipath error is highly dependent on location hence no exact value can be obtained for the variance of it. An upper bound as we have already seen is one cycle for a double difference and quarter of a cycle for one receiver itself. In addition, measuring the multipath from actual data can also be tricky, as a sidereal differencing of double difference residuals does not entirely cancel multipath from one day to the next. This is because double differences needs two satellites and as the repeat times of these two satellites will not be identical, strictly speaking there is no exact time when the satellites will be in the exact same locations with respect to one another. Others have realised the same problem and techniques to transform double differences into single differences and then align each satellite with a much more exact repeat time then transform these back into a double difference have been designed |62j. However, we do not use this method and is beyond the scope of this thesis; therefore we make the approximation that the alignment of the satellites are exactly the same from one nominal satellite repeat period to the next. With this approximation using sidereal differencing we can obtain an approximate idea of the magnitude of multipath error and other errors that we cannot correct with our particular method of sidereal filtering.

Figure 3.7 shows the sidereal difference of two consecutive days along with the two days residuals.

Ill CHAPTER 3. MUUTIPATH AND SIDEREAL FILTERING

Time (Hours)

Figure 3.7: Sidereal differenced double difference residuals for one SV pair over two distinct days

Performing this sidereal differencing calculation for all pairs of satellites used for the final daily solution for one day at the LuShan testbed site, an average value of the sidereal variance was determined to be 0.0014 while that of day two was 0.0096. The sidereal variance can be written as var (sidereal difference) = 2 x var (noise) + 2 x var (other) while the variance of day 2 can be written as var (day 2) — var (noise) + var (other) + var (multipath that we can correct) Which in turn implies

var (ski ereal difference) var (multipath that we can correct) var (day 2) ¹

y ' 2

and var (other) = var (day 2) — var (noise)— var (multipath that we can correct) Hence var (multipath that we can correct) = 0.0089 cycles ² and var (other)— 0.000531 cycle

Making the approximation that the alignment of the satellites are exactly the same from one nominal satellite repeat period to the next, we say that the CHAPTER, 3. MULTIPATH AND SIDEREAL FILTERING average multipath in the signal at the LuShan testbed for two receivers used had a standard deviation of 0.094 cycles, while that of other errors had a standard deviation 0.023 cycles.

So the standard deviation of 0.094 cycles for the multipath is a lower bound for this particular instance while 1 cycle is an upper bound. Therefore it appears that measurement noise is at least an order of magnitude less than multipath error hence fairly insignificant. The other errors undoubtedly have some component of multipath but are hard to determine, we combine the other errors with the measurement noise which we call the miscellaneous errors, while that of the multipath errors that we can correct, we call the multipath errors. Hence in this particular circumstance we obtain the following standard deviations for these two errors.

cr _e - 0.026 Cycles

σ _Μ - 0.094 Cycles

Qualitative standard deviation values

While of course this is dependent on the GPS receivers surroundings it gives an idea of the magnitude of both multipath errors and ones that are not multipath. The multipath error here is a lower bound so is conservative. Even so, the repeating nature of the errors we have examined here which we lump together into multipath are still many times greater than that of the rest of the transient errors combined. Therefore some form of sidereal filtering seems advantageous to use.

An interday sidereal filtered solution (an interday solution) is a daily solution using all available observations in other days in addition to the all observations available in the day in question itself; the other days may not necessarily be consecutive days.

The satellite repeating nature allows us to treat double differences obtained on different sidereal days as if they were almost the same. When the satellites CHAPTER 3. MULTIPATH AND SIDEREAL FILTERING location axe the same, the delta norm and the multipath M are also the same if both of the receiver locations and their surroundings are unchanged. In this case if the double differences used today to obtain a position solution are also the same as the double differences obtained tomorrow to obtain a solution, then the two solutions should be very similar.

We define subscripts k and n to denote the sidereal day. For simplicity we do not give a notation that represents time within a day. For two double difference to match one another we define that they must be on different sidereal days separated by no more than a second with respect to an integer multiple of £¾ _> , and the double differences belong to the same pair of satellites and receivers, where a _v is nominal GPS satellite orbital repeat period. When double differences match one another we make the approximation that the satellite geometry is exactly the same for both double differences. We define N as containing every possible delta norm for the day, while N _¾ contains only the delta norms used on day k, similar notation is used for multipath, M for every possible multipath during a day and for multipath used on day k.

We can factorize multipath both for all double differences in the day and for just the ones sampled on day k; namely M — NP^ + AMNL, and M¾ = NfcP _¾ f.fc + M,¾,vi,i. This allows us to rewrite equation 3.1 as follows which means our solution is biased due to multipath.

λ _¾ = An _k ^■ [P _k + P _M> *) + M _MNL ,k + ¾ (3.2)

Therefore we can say for day k the solutions take the distribution as follows assuming multipath is normally distributed with zero mean, where <¾ is the percentage of daily double differences used for day k.

Pfc ~ N (p _k + P _M>K ,— σ? ( JNFNV ¹ ) (3.3) Daily solution distribution for day k

If the same double differences are used every day then Ρ · becomes fixed and the variance of the solution is due to solely the number of double differences CHAPTER, 3. MULTIPATH AND SIDEREAL FILTERING used and the miscellaneous errors. Whilst the solution would be a biased solution, for land deformation monitoring where relative movements is the merit of interest, this fixed biased would be irrelevant hence the solution would be as good as one could hope.

If we assume we cannot obtain the same double differences every day and instead we randomly sample a percent of all the daily double differences with a uniform distribution, then we have only a probability that given a double difference obtained today there will exist a matching one tomorrow, then our solutions should be distributed as follows given enough days such that the average value of Ρ , _/C more closely resembles that of Ρ _Λί as we expect Ρ , _/C ~

P _FC ~ N ( P, - ! ^> ,,- . ((I - a) σ + <ή) - ( V ) ) (3.4) Daily solution distribution given uniformly random double differences

In this case over enough days we expect a fixed bias but with a variance larger than before. We are now faced with two problems the first is obviously the increase in variance, while the other is how the daily multipath bias - changes from one day to another. The GPS satellite repeat time is about four minutes lagged every day with respect to 24 hours, therefore if we could obtain observations throughout all daylight hours we would expect a slow change in the daily multipath bias. The solutions in this case could appear similar as if the actual location of the node was moving. As an example we performed a simulation where we generated normally distributed multipath, delta norms, and miscellaneous errors where the miscellaneous error was significantly smaller than the multipath error. We set a fixed node location of [1, 1, l j . calculated the observable for all epochs, used 50% of consecutive epochs everyday and shifted this by 0.3% every day, thus simulating obtaining observations only during daylight hours. Figure 3.8 below shows the plot for the x coordinate of the solutions using our simulation. "Phat" is for P _¾ , "Pm offset" for P + P _¾ and "Actual" for P _¾ . As can be seen there is a slow movement where the deviations from P _M can be great. The distribution of the solutions matched the distribution as given in equation 3.4. Visually it is hard to tell whether or not the node is stationary or moving. Also, as can be seen from the graph the CHAPTER 3. MUUTIPATH AND SIDEREAL FILTERING solutions have a repeat period of approximately one year which is due to the four minute shift every day.

Figure 3.8: GPS model simulation of solutions for 50% of day observations and 4 minute daily shift

In comparison the same simulation was run again with the same parameters except exactly the same double differences were used everyday instead. The x coordinate of the solutions are plotted in Figure 3.9. As can be seen the remaining variance is due to miscellaneous errors as given in equation 3.3 when P\f,k is fixed.

CHAPTER, 3. MULTIPATH AND SIDEREAL FILTERING

Mode! simulation 50% of observations per day, same used everyday

1 .035

Phat

1 .03 Actual

Pm offset

1 .025

1 .02

1 .015

1 .01

1 .005

1

0.995

0.99

100 200 300 400 500 600

Day

igure 3.9: GPS model simulation 50% of observations and the same matching ouble differences used evervdav

Therefore our scenario the day to day the distribution may follow equation 3.3 but in the longer term the distribution may follow equation 3.4.

Therefore to reduce solution variance over the short and long terms, observations should be taken once every second 24 hours a day and data processed using the same double differences everyday. With permanently powered GPS receivers it is very easy to obtain the same double differences each day and obtai one every second 24 hours a day as there is enough power to do this. However, in our scenario we cannot simply turn our GPS receiver on in the middle of the night just because there is a double difference there that we need to match the double difference we used yesterday to calculate a solution. Therefore we would expect for our scenario a rather small value causing an increase in the solution variance, while that of permanently powered high rate GPS receivers sampling at once a second would have a large value producing solutions with less variance. We have a biased estimator with a large amount of variance. We now take advantage of the fact that over the long term M* is a random variable.

Using our model from equation 3.1 we see that the residual in cycles for sidereal day k, evaluated at P _# with an actual location P _¾ , can be written as follows.. CHAPTER 3. MULTIPATH AND SIDEREAL FILTERING

r _k = (Δ ^η · (P _fc - P _B ) + M + e _k ) (3.5)

For sidereal differencing using double differences, today's residual minus yesterday's residual results in the displacement since yesterday's observation; this is written as follows.

(r _n - r _n→ ) = - ^» (An · ι Δ Ρ,.. ) + 2e) (3.6)

A

The additive noise is double, the multipath is removed, and using LS results in ΔΡ _η the estimate of the displacement between yesterday and today. It's a very precise solution but the solution is for displacement and not for position, hence very good for precise epoch based solutions but gives nothing for daily solution improvements. Therefore, sidereal differencing removes the position information. If one wished to obtain a position solution one would have to add yesterday's position to the displacement which would once again add the multipath error back into the equations.

A (r _n — r _n _i) + Δη - P., = Δη ^• [ \ , + ΔΡ _η j + 2e

= Δη ^• (p _n _ _! -{-- ΔΡ„) + e _n - e„_i

- An ^• (P _n ... _! + P _M ,n-i ) + Δ ^η · (ΔΡ _η ) + e _n

^■ P„ + Δ ^η■ (P ,n-i) + e„ (3.7)

Averaging over many such days we obtain the following.

ave (λ (r _n — r _k ) + Δη · P _¾ = ave (Δη · P _n + Δη · + e _n

= Δη ^■ P _n ·-[-- Δη ^■ ave (P^ _¾ ) + e _n (3.8)

Over enough days we expect ave ( MA) = Ρ _λ ί then

ave (A (r _n - r _k ) + An · P _k ) - Δη · P + Δη · P _M + _n (3.9) CHAPTER, 3. MULTIPATH AND SIDEREAL FILTERING

Therefore using LS to obtain an estimate P _n for P _n using this equation we would then expect following distribution.

Which is identical to the distribution when we had used the same double differences every day. This of course is an upper bound to accuracy and precision. If we consider arithmetic mean averaging over many such days things become more complicated as there is only a probability that we can form the residual difference r _n — r _k in the first place. Because of this, we define the bar notation as the arithmetic mean calculated when possible including the current day n. For r _¾ at day n this is f _n = mean (Ω U {r _rt }) where Ω contains a set of matching residuals; in this particular instance {r _n } represents the set containing only r _n and ensures f _n is sensibly defined. Whilst the bar notation is valid for nonconsecutive days we only consider rn consecutive days less than or equal today n here for simplicity; this we take as read from here on. Therefore, more generally we write the following.

X (r _r r _n ) -]-- An ^■ P _n ------- An · P -I- An ^■ Pj. _/ .„ - e (3.11)

Sidereal filtering equation

The distribution of ΡΜ, _Π will be the same as M but with a scaled variance due to the days included in the average. For small values of m, PM, _« can approximately be written as follows.

Therefore, the LS estimate P for P _n + P _M> using equation 3.11 has a distri- bution as follows. CHAPTER 3. MUUTIPATH AND SIDEREAL FILTERING

In our scenario we cannot randomly saxnple a observations throughout the day as the nodes are only able to operate during daylight hours. If we define the percentage of daylight hours with respect to an entire day as β, then we are only able to randomly saxnple a percent of the observations that axe within daylight hours (i.e we obtain β percent of all the possible observations each day and all of these are in daylight hours). Every day the observations that β refer to are shifted by four minutes, therefore the probability of being able to form a, residual falls the further away in time day k is away from day n. Hence the more days the average is run over the smaller the increase in the number of items in r _n until given enough days there is no overlap and there is zero probability that adding that day to the average will make any difference whatsoever. Also, we note that f _n is the average of fewer things than if there was no four minute shift every day no matter how many days we averaged over. In this case we have some function Q (a, β, πι, δ) where we can write the distribution of the solutions as follows.

P« ~ N fp„ + PM, (g (cr, β, m, 6) σ% -\- f) J- (N ^T N) ^{" 1} ) (3.14)

We do not obtain a explicit formula for Q (a, β, πι, δ) instead we simply note that is it is positively correlated with β and rn.

During daylight hours the nodes obtain observations but do not obtain them uniformly during daylight hours. Instead they obtain many more around the middle of the day and fewer in the morning and evenings. However, to keep things simple we do not consider this extra complication and instead continue to assume uniform acquisition of observations during daylight hours. In the Paekakariki testbed in New Zealand over a two week period we measured the average percentage of epochs obtained during daylight hours as 30%. Therefore, for a simulation we use a— 0,3, β— 0.5 and δ— 0,003. In addition we move the node by [0.050.05 0.05] halfway through the simulation. We then use 15 days to perform the sidereal averaging (m = 15). Using the same simulation settings as before we obtain the following simulation result in Figure 3.10. CHAPTER, 3. MULTIPATH AND SIDEREAL FILTERING

Day

Figure 3, 10: Si.mulati.on of sidereal filtering using 15 days with intermittent daylight observations (a— 0.3)

As can be seen on average there is a slight solution bias and the sidereal filtered solution perform significantly better. As a comparison we performed the same simulation again but this time using = 1 signifying that all double difference during daylight hours were used.

Day

gure 3.11: Simulation of sidereal filtering using 15 days with continuou aylight observations (a — \ . — 0.5)

From Figures 3.10 and. 3.11 we see a low frequency noise component and a CHAPTER 3. MUUTIPATH AND SIDEREAL FILTERING high frequency noise component. The low frequency component is clearly the same, only the high frequency component has been reduced when using interday sidereal filtering. The low frequency component appears as a slow moving solution bias that repeats every year and is due to the four minute daily shift. The lo frequency component is problematic because it looks like a GPS receiver node slowly moving in time. From these figures we also see that in the long term interday sidereal filtering is only effective when the high frequency noise component has a magnitude significantly greater than the magnitude of the low-frequency component. This is equivalent to saying interday sidereal filtering is only really effective at improving daily solution accuracy when the magnitude of the usual solution precision is significantly greater than the magnitude of the slow moving bias that repeats every year.

To examine interday sidereal filtering on the the high frequency noise we look at just the noise obtained by using one day differencing. To be explicit P _k — P _f c_i and P _k — P _k -i- The the difference between consecutive days produces approximately double the variance of the high frequency noise that exists in the solutions and eliminates the low frequency component. From this we can then see the reduction of high frequency noise and hence the improvement in the short term solution accuracy. This time we run the simulation without the addition of the movement halfway through the simulation to avoid an outlier. Figure 3.12 shows simulations of P _k — P _k -i and P _k — P _k -i for an alpha of 1 and 0.3 while β = 0.5 and δ = 0.003 and rn = 15 as before; these correspond- to the settings of Figures 3.11 and 3.10 respectively.

Day Day

Figure 3.12: Day differenced solutions for = 1.0 and a = 0.3 while β = 0.5 and δ = 0.003 CHAPTER, 3. MULTIPATH AND SIDEREAL FILTERING

Figure 3.12 shows only the high frequency noise component once the low frequency noise comrjonent has been removed. As can be seen even when all double differences during daylight hours are used (a— 1.0) there is a reduction in the high frequency noise component when using sidereal filtering. However, the improvement in the high frequency noise comrjonent is more significant when only using a fraction (30% in this case) of the double differences during daylight hours. In this case, using all double differences during daylight hours caused a reduction in the high frequency noise component of approximately 50%, while when using only 30% of double differences during daylight hours produced a reduction in the high frequency noise component of approximately 80%. Therefore, interday sidereal filtering is effective at reducing the high frequency noise component when using any fraction of observations during daylight hours. A reduction in the high frequency component is of course an improvement in solution precision.

We finally examine the effect interday sidereal filtering has on solutions where every double difference in a day is used to compute the solutions; in this case a ------ 1 β ----- 1 and δ is irrelevant. Running this through the simulation we obtain Figure 3.13.

Day

Figure 3.13: Simulation of sidereal filtering using all double differences in a day ( = 1, 3 = 1)

As expected the solutions before and after sidereal filtering are identical and there is no improvement in accuracy or precision. In addition there is no slow moving solution bias that repeats every year. CHAPTER 3. MUUTIPATH AND SIDEREAL FILTERING

It appears that interday sidereal filtering improves precision when obtaining observations solely through daylight hours. The improvement in precision is most marked when only a few observations are taken during daylight hours. Interday sidereal filtering may or may not significantly improve accuracy depending on the magnitude of the slow moving bias that repeats every year compared to the magnitude of the usual solution precision. The magnitude of the usual solution precision of course at least in part depends on the number of observations obtained during daylight hours.

In summary, we see added benefit from using sidereal filtering for our unique intermittent observation scenario that would otherwise not be available to us had it not be intermittent. In addition because our nodes are duty cycled even during daylight hours we obtain more benefit from interday sidereal filtering than we would if we were continuously obtaining observations throughout daylight hours. However, we also observe a slow moving solution bias that repeats once a year when using only daylight hour observations that cannot be reduced using our form of interday sidereal filtering. This slow moving bias affects accuracy while for short periods of time it appears as a fixed bias.

We now proceed to verify that interday sidereal filtering does indeed improve solutions when intermittently obtaining observations with the use of real life observations.

3.3.1 Simulation validation

We have seen for simulations based on our model of muitipath, that for our particular unique hardware where observations are only available during the day, and even then only intermittently, that we predict our solutions should take the form which is seen in Figure 3.10. Firstly we predict the sidereal filter should respond instantly to a change in node location and the solutions should match that of the physical displacement performed. Secondly we predict that we do obtain an improvement in daily solution precision using the interday sidereal filtering described in this chapter when intermittently obtaining observations and that no improvement can be got when using continuous high rate observations. Thirdly, that there should be slow change in solution bias that will repeat once every year that our interday sidereal filtering cannot reduce. The third prediction we could not verify as it would require more than CHAPTER, 3. MULTIPATH AND SIDEREAL FILTERING a year of observations of a static node. Fourthly that any improvement is more pronounced when using fewer daylight observations than many daylight observations. The fourth prediction we do not validate here.

The first two predictions of response and improvement we validate as follows.

.3.1.1 Sidereal response

An experiment was run with a node as described in the hardware chapter 5 but modified to be run off a continuous power supply. The antennas were placed on a roof and the receivers were mounted beneath the eaves aw ^T ay from direct sunlight and the elements, the baseline was approximately 10 m. With continuous power it is possible to obtain the highest possible accuracy and to determine whether or not the sidereal solution moves instantaneously when the underlying movement of the node changes, in addition, that the change of the sidereal solution is consistent with the underlying movement. The node's antenna was moved in a NNW direction and slightly down; firstly by 20 ± 2mm and then by 5 ± 1mm. Using every epoch of the day, the primary principal component of the solutions were plotted with respect to the day the solution was calculated for. The Figure 3.14 shows this plot where the sidereal filter had a moving window period of 15 days.

Day

Figure 3.14: Primary principal component of daily solutions given near continuous 24-hour observations CHAPTER 3. MULTIPATH AND SIDEREAL FILTERING

The node's antenna was moved approximately at midday on the 357th day; this can be seen in Figure 3.14 above, as a point that is neither at -0.01 nor 0.01 as some double differences were obtained while at the former location and others at the latter location causing the final solution to be neither one nor the other. The 5 mm movement was performed at 7 p.m. on the 368th day. As can be seen the sidereal filter responds instantaneously to movements as expected, but unlike simulations, for larger movements like the 2 cm movement, the sidereal solution underestimates the movement and slowly adjusts to the new location. This makes sense because multipath is a function of location, hence, the previous days used for the sidereal solution have slightly incorrect multi- path information. For simulation we assumed that multipath was independent of location, hence why we do not see the slow linear rise of the sidereal solution that we see in Figure 3.14 in Figure 3.10. For the 5 mm movement the sidereal filter did not have sufficient time to discard ail the observations before the 2 cm movement, 4 days we used before the 2 cm movement was performed when calculating the first sidereal solution after the 5 mm movement. From Figure 3.14 we estimate the magnitude of the movements of the sidereal solutions by examining the sidereal solution before the movements were performed and obtaining an estimates for the positions before the movements were performed; these values were—9.9 ± 0.2m?n,8.0 ± 0.4mm and 12.6 ± 0.2mm in chronological order with a 95% confidence. From this we estimate the first movement was 17.8 ± 0.6mm while the second movement was 4.6 ± 0.6mm both with 95% confidence; therefore, both estimates are consistent with the estimates obtained from using a ruler. This shows for relatively small movements over the short timeframes at this particular test site at least, that sidereal solutions as far as we can determine are consistent with the underlying movement of of the GPS antennas themselves.

3.3.1.2 Sidereal solution improvement

All our static real life observational tests were relatively short lasting no more than a couple of months. This means we were able to see the high frequencies of the solutions but not the low frequencies. This means we could only observe solution precision improvement not accuracy. This of course we predict is where interday sidereal filtering will show a clear improvement.

According to our multipath model the solutions obtained when using the same CHAPTER, 3. MULTIPATH AND SIDEREAL FILTERING epochs every day to obtain solutions such as when using a permanently powered GPS receiver taking high rate observations interday sidereal filtering should have no effect. However, from Figure 3.14 this appears to not be so. If we consider just days before the 2 cm movement and after the 5 mm movement we see an improvement in 2DRMS of approximately 30%. We theorise that this improvement that we did not predict in our simulations maybe due to the fact that the satellite repeat period is not exactly 24 hours -246 seconds but instead varies for each satellite causing the multipath of the double differences not to entirely repeated each day. Plotting the first two principal components of the solutions when using all observations from the permanently powered receiver Figure 3.15 was obtained where this improvement can be seen more clearly.

x 10 ^~3 Without sidereal

PC1 (m)

Figure 3.15: First two principal components of daily solutions given all obser vat ions from a permanently powered receiver

To investigate the improvement when observations are not continuous we used the observations obtained from the continuously powered node and selected only ones that had a time in common with a solar powered GPS node placed nearby. Therefore, the distribution of activity of the permanently powered observations would match that of the observations obtained from the solar pow ^T ered GPS node itself. In this way w ^T e eliminate the possible dependence of location on improvement. Once again, considering just days before the 2 cm movement and after the 5 mm movement we obtain an improvement CHAPTER 3. MUUTIPATH AND SIDEREAL FILTERING in 2D RMS of approximately 60% in this ease. This improvement is double that of the unexpected improvement obtained when obtaining observations continuously once a second. Plotting the first two principal components of the solutions when only using observations from the permanently powered receiver when the solar powered GPS receiver node obtains an observation we obtain the following plot where this improvement can be seen. 10 ^"3 With sidereal

PC1 (m)

Figure 3.16: First two principal components of daily solutions given only observation times that match a solar powered GPS node

As can be seen from both of these figures (3.1 -5 and 3.16) using sidereal filtering produces solutions reasonably similar to one another. The 5 mm step without using sidereal filtering in the second figure (3.16) cannot be easily determined while it is easy to determine when using sidereal filtering. As desired the improvement when using intermittent observations during daylight hours is considerably greater than that obtained when using a permanently powered GPS receiver.

^" Finally we compare sidereal solutions with non-sidereal solutions for an actual solar powered GPS receiver placed at the Paekakariki testbed and left for 22 days. This GPS receiver node is one as described in Chapter Chapter 5. A birds eye view of the solutions can be seen in Figure Figure 3.17 on page 68. CHAPTER, 3. MULTIPATH AND SIDEREAL FILTERING

East (m)

Figure 3.17: Effect of interday sidereal filtering from observations obtained from intermittent solar powered hardware on solutions (one square centimeter)

As can be seen sidereal filtering has a significant effect on solution precision. An improvement of approximately -50% in 2DRMS from 5.6 mm to 2.8 mm was obtained. This improvement percentage in 2DRMS is comparable to the improvement we saw when using the permanently powered GPS receiver using only observations taken when a solar powered GPS receiver node was active.

An intraday sidereal filtered solution (an intraday solution) is a epoch based solution using observations in other days in addition to the observations in the second in question itself; the other days may not necessarily be consecutive days.

Intraday sidereal filtering is not a primary concern of ours as our system is not competing with real time systems such as RTK. In addition, the decision was made to concentrate on long term land deformation measurements using daily solutions and doing this well. However, there is no reason why our system can't be extended to epoch based solutions but more work would be needed to investigate this type of solution. In the section we briefly examine intraday solutions using past information and the present epoch to calculate a solution. CHAPTER 3. MUUTIPATH AND SIDEREAL FILTERING

Calculating the muitipath estimate for the intraday algorithm is performed by using previous days where the location of the node is unchanged to which we have interday sidereal solutions. These interday solutions are very accurate with the residuals mainly consisting of muitipath and hence we can calculate the muitipath for those previous days. This muitipath information we have then gets combined using simple geometric averaging to create a muitipath map to which we can correct any double difference for the present day. Rather than correcting the double differences themselves, we turn the double differences into residuals using a biased location, correct them, then convert these back into double differences to which least squares can be used. This is done so as to avoid using unwrapped double differences as was done in the interday sidereal filtering algorithm. Assuming we have an estimate to the true muitipath M as M' this procedure is summarised as follows initially for unwraped double differences.

• r _n ----- (ΑΦ _η — An ^■ P R ) definition of of residual for location P _B for day n

• Xr _n — „ _j = λΦ _η — Δη · Ρβ— Μ _η ' ___-, use muitipath map to correct residuals

• Xr _n - ML_ _X = An ^■ P„ + M _n + e - Αη · Ρ _β - Μ _η '_ _χ substituting model for Φ

• Ar _n — ^_ _] +Δη·Ρβ — Δη·Ρ ' _n +e+ AM _n defining AM _n as muitipath map error between today's actual and yesterday's estimate of the muitipath

• Xr _n — M _n '„ _x + An ^■ Pp, f¾ An ^■ P„ approximation if good muitipath map and noise is small

To obtain the muitipath estimate M' geometric averaging was used due to the simple implementation and is summarised as follows.

• m„ = ί λΦ — Δη - P j definition of of residual for location P„. This residual is close to muitipath if the receiver did not move on day n.

• M _n ' ----- (1 — ) Μ _η ^! _ + \m _n geometric averaging for predicting the multi- path of tomorrow. CHAPTER, 3. MULTIPATH AND SIDEREAL FILTERING

Therefore, two equations are used for our intraday sidereal filtering; a corrected equation based on yesterday's multipath estimate and an update equation of the multipath estimate, these two equations are equations (3.15) and (3.16)is not based.

Xr _n - Μ _η ' _ + Δ ^η - P _B « Δ ^η · P

Intraday sidereal filtering equations

Rather than using unwrapped residuals we chose instead to use wrapped residuals so as to avoid cycle slips as we did in the interday sidereal filtering. Upon using wrapped residuals our algorithm can be described as listed in the following algorithm.

CHAPTER 3. MUUTIPATH AND SIDEREAL FILTERING

Algorithm 3.1 Intraday sideral algorithm

Correct

* Shift multipath map m temporarily by an integer multiple of a _p to align them with today, day n

* Calculate todays wrapped residuals for epoch z ,r , _i2 = / (^· (λ _η , _ζ — N _z ^■ P _n -i ) ) /Cf based on yesterday's interday sidereal solution P _i

* Let ϊ _η>ζ — magmin (r„ _S — ηι„,__. _{1 ζ} , r _z + π½_.ι _¾ ) be the corrected residuals for this epoch

Solve

* Use LS to solve Xv _n ,z + N _¾ · P _n _i ~ N _z · P _UiZ if sufficient independent equations and PDOP is low enough, hence obtain an estimate of P _n . _z

* Use Kalman smoothing on the LS solution where the model is stationary spherical and the measurement covariance is taken from the LS solution

* Call the solution after smoothing P _n , _z Update

* At the end of the day update the multipath map in as follows if no movement has happened throughout the day

» Calculate todays wrapped residuals r _n = / ( ··- ( ΛΦ,, · I *,, ) ) /( ', based on today's interday sidereal solution P _n

* Update map _n = (1 — a) m _n _ _! + \ r _n \

3.4.1 Intraday sidereal filtering improvement

We compare epoch based solutions, with and without intraday sidereal filtering. Figure 3.18 below shows epoch based solutions without intraday sidereal filtering on observations obtained from the LuShan testbed by our wireless solar powered GPS receiver nodes as described in the hardware chapter Hardware. Epoch's with insufficient equations and/or excessive Position Dilution Of Precision (PDOP) were not processed. In addition Kalman smoothing was CHAPTER, 3. MULTIPATH AND SIDEREAL FILTERING performed on the solutions where the model was stationary spherical and the measurement covariance taken from the LS solution as (N{N ₇ ) . As the hardware was intermittent only a fraction of solutions during the day could be calculated producing many seconds where no solution could be calculated. Because of this these the x axis only represents the epochs where solutions could be calculated, Figure show the results 3.18.

Figure 3.18: Epoch based North solutions without intraday sidereal filtering

Figure 3.19 below shows the same observations as was used for Figure 3.18 but with intraday sidereal filtering instead. The algorithm had 10 days beforehand to adjust its multipath map. CHAPTER 3. MULTIPATH AND SIDEREAL FILTERING

As can be seen the effect the multipath has been considerably reduced on a per second basis. For the solutions themselves an improvement of approximately 60% was obser ved in 2DRMS accuracy, while an improvement of approximately 80% was observed once the solutions themselves had been run through Kalman filter. For this test the 2DRMS precision measurements are tabulated in the following table, Table 3.1.

Table 3.1: Epoch based accuracies with and without intraday sidereal filtering CHAPTER, 3. MULTIPATH AND SIDEREAL FILTERING

3,4.2 Intraday sidereal filtering response

3.4.2.1 5 mm permanent power

The response of the intraday sidereal filtering with respect to movement was first examined using the 5 mm movement of the permanently powered GPS receiver as seen in section 3.3.1.1. Firstly the primary horizontal component of the movement without intraday sidereal filtering over four days was plotted in Figure 3.20 with a one hour moving average filter applied to the solutions.

5mm movement for permanently powered receiver. No intraday sidereal filteri ng.

Figure 3.20: 5 mm movement no intraday sidereal filtering with a permanently powered GPS receiver

The same data was processed once again identically but with the addition of intraday sidereal filtering with an of 0.25; the results can be seen in the Figure 3.21 with the same scale as the previous figure, Figure 3.20. CHAPTER 3. MUUTIPATH AND SIDEREAL FILTERING

5mm movement for permanently powered receiver, = 0.25

Day

^" Figure 3.21 : -5 mm movement intraday sidereal filtering with a permanently powered GPS receiver a— 0.25

As can be seen, without intraday sidereal filtering it is very difficult to determine that movement has taken place, while, with intraday sidereal filtering, the movement can be detected more easily. The time of the movement was determine by Figure 3.21 as being between 6 p.m. and 8 p.m. on the 3 ^rd of January while that of the actual movement was at 7:10 p.m. on the 3 ^rd of January.

Performing Kernel Density Estimation (KDE) on Figure 3.21 Figure 3.22 was obtained. CHAPTER, 3. MULTIPATH AND SIDEREAL FILTERING

KDE plot i n primary honzonta! componentof 5 mm movement second based processing 400

X: 0.001343

350 Y: 382.7

; X: -0.001855

i Y: 358.2

300 i \ 250 ί \ g 200

o

150 100 50 V 0

0

Separati on (m)

x 1 0

Figure 3.22: KDE plot of the 5 mm movement using intraday sidereal filtering with the permanently powered GPS receiver

Figure 3.22 shows that the movement was approximately 3.2 mm. The slope of the roof as measured by an inclinometer was 15° therefore the total movement should have been 3.3 mm which is 1.7 mm out from the expected value of 5 mm.

Finally the intraday sidereal filtering was performed again but this time using an a value of one. The results of this can be seen in Figure 3.23.

CHAPTER 3. MUUTIPATH AND SIDEREAL FILTERING

5mm movement for permanently powered receiver. Intraday sidereal filteri ng a =

Day

^" Figure 3.23: -5 mm movement intraday sidereal filtering with a permanently powered GPS receiver a = 1

The difference between Figures 3.21 and 3.23 is the impulse response. It can be seen in ^" Figure 3.23 that there is a sharp dip exactly one day after the the movement that was not noticeable in Figure 3.21. This is expected as when = 1 only the previous day is used for obtaining a multipath estimate and the requirement when obtaining the multipath estimate is that the node stationery for that day which is clearly not the case when the node moves. When a— 0.25 less of the multipath estimate is based on the previous day and hence why we do not notice any dip one day after the actual movement in Figure 3.21.

3.4.2.2 2 cm solar power

A step response test using a solar powered node over four days was performed at the Paekakariki testbed in New Zealand. The first day was to record a multipath map, the second day to make sure the multipath mitigation was successful and to obtain a control, a third day where the node was moved, and a fourth day to monitor the response of the intraday sidereal filter. An a of CHAPTER, 3. MULTIPATH AND SIDEREAL FILTERING

1 was used for the intraday sidereal filtering algorithm which meant a finite impulse response and only one day after the step would be needed to see any type of response.

The movement of the node was in approximately a north direction and the movement as measured by a ruler was approximately 2 cm. The movement was done at 11:06 AM on the 5th of December 2014. Using the intraday sidereal filtering algorithm, a linear movement was seen starting at approximately 11 AM and finishing at 12 PM but it was difficult to determine as the solutions were only intermittent. As the motion was instantaneous, we could say that the first sign of the linear movement should be at or after the actual instantaneous motion as caused by us was performed. Therefore we determined from the intraday sidereal filtering algorithm that the instantaneous motion happened at approximately 11 AM which is consistent with the actual movement as performed by us.

The Figure 3.24 shows the density of the solutions over the second to the fourth days inclusive. As can be seen there is a clear movement in almost a due north direction; the displacement between the two locations as measured from Figure 3.24 below is approximately 1.8 cm. This displacement is slightly less than the 2 cm as measured using the ruler.

intraday sidereal filtering

Figure 3.25 below shows the primary principal component of the horizontal CHAPTER 3. MUUTIPATH AND SIDEREAL FILTERING solutions around the time when the movement was performed.

Figure 3.25: Epoch based intraday sidereal filtering around time of the movement of a solar oowered GPS receiver node

As can be seen the intermittent nature of using solar powered GPS receivers means if a movement happens when the GPS receiver is not on or actively obtaining useful observations no movement will be measured until after the GPS receiver becomes active again and obtains a sufficient number of observations.

Figure 3.26 below shows the primary principal component of the horizontal solutions for all four days.

CHAPTER, 3. MULTIPATH AND SIDEREAL FILTERING

Figure 3.26: Epoch based intraday sidereal filtering over four days of a solar powered GPS receiver node

As a comparison, without using the intraday sidereal filtering, the density of the solutions of days 2 to 4 inclusive were plotted as before in Figure 3.24 but without intraday sidereal filtering; this can be seen in Figure 3.27 where the scale of the plot is the same as in Figure 3.24. As can be seen in this case it is not possible to determine the two locations. without intraday sidereal filtering

-1 .935 -1 .93 -1 .925 -1 .92 -1 .915

East (m)

Figure 3.27: Birds eye view of the density of the epoch based solutions without intraday sidereal filtering

GPS receivers for precision applications output the following: LI pseudorange observable, LI phase observable, LI Doppler observable, LI signal strength observable, loss of lock indicator, quality indicator, satellite number, and the observation time as given by the receiver's clock. Usually a great deal more is given out. While not all are needed for high accuracy solutions, the pseudo range, phase, satellite number and time would be regarded as the minimum needed for some sort of high accuracy solution.

A typical breakdown of data size would be, the pseudorange, phase and time are all eight byte values, Doppler is four bytes, signal strength one byte, loss of lock indicator one byte, quality indicator one byte, and satellite number one byte.

For an epoch with n many observations to different satellites usually the time would be only given out once. Table 4.1 shows a breakdown of the number of bytes sent in this situation. CHAPTER 4. ALGORITHM COMPARISON

Table 4.1 : Typical breakdown of the minimum number of bytes sent per epoch from the receiver to the user for standard algorithms.

While n could be anywhere from 0 to 32 for one epoch a typical value would be between 4 and 10. Assuming an upper bound of 10 this would produce up to 248 bytes a second per node.

Low power wireless network solutions have relatively slow data rates restricting channel capacity. This is partly due to the fact that slower data rates mean larger transmission ranges for the same amount power and the same method of modulation. In addition, less bandwidth requirement for nodes results in lower epoch data loss thus mitigating the effect of a lossy channel. Because of channel capacity, power usage, and epoch loss, minimizing the amount of data that the nodes actually sends is advantageous. Another advantage of less bandwidth requirement is more nodes can occupy a channel increasing the number of nodes in a network on a given channel.

To reduce the amount of data from the GPS receivers themselves, some of the CFFS algorithm is implemented on the nodes themselves thus changing the data sent over the RF link. Algorithm 2.1 is the algorithm run on the nodes.

The extrapolated wrapped phase observable, satellite number and the second aligned time are the main variables returned. In addition, occasionally an approximate code based solution is returned; this can be done as infrequently as once or twice a day. This is so infrequent that effectively it contributes nothing to the needed bandwidth and can be ignored.

The main part of the algorithm has an approximate knowledge of the time via the Internet or some other mechanism and the node section has already CHAPTER 4. ALGORITHM COMPARISON aligned the time to on the second. If we assume that the main part of the algorithm knows the time to within one minute and the latency between the node section and the main section of the algorithm is also less than a minute, then it is possible for the main part of the algorithm to calculate the time of an observation if the node section sends just a modulo 256 version of the time. As these timing requirements are very easy to achieve, only one byte per epoch is needed to be sent between the node and the main section of the algorithm.

The wrapped phase observable removes the integer part of the phase observable and leaves a number between zero and one. This transforms a number that was in the order of 20,000 needing an accuracy of a fraction of one, to a number in the order of one needing a fraction of one; this means we only need a couple of bytes to encode the remaining information just as accurately and for our purposes 2 bytes should be sufficient. 2 bytes means we have the phase resolution of about one in 65,000.

The following table shows a breakdown of the bare minimum number of bytes that a node could send for the CFFS algorithm.

Table 4.2: Breakdown of the minimum number of bytes a node could send per epoch for the CFFS algorithm

As before, assuming an upper bound of 10 satellites in one epoch this would produce up to 31 bytes a second per node. This equates to an 87.5% reduction in required bandwidth for the same given epoch rate.

4.1.1 Real life reduction of bandwidth comparison

Low cost GPS receivers nodes powered by small solar panels and batteryiess using low power radios were designed and implemented. The radio links had no packet loss correction, the packets were sent using a broadcast method but if they did not arrive they would not be sent again. The GPS receivers CHAPTER 4. ALGORITHM COMPARISON themselves were singleband Ublox LEA6T while the antennas used were 25 mm patch and 20 mm helical antennas. The receivers were placed on the ground outside with a high amount of surrounding trees and buildings in winter. The nodes received no direct sunlight due to the time of year and the tali trees. The Sun's maximum elevation was 27° at noon and the day length was 9.5 hours. During these 9.5 hours of the day the nodes on average were able to have their GPS receivers on approximately 20% of the time (about 2 hours a day). A receiver was placed approximately 10 m away from the nodes and data gathered. A permanently powered reference receiver was used to maximize overlapping epochs of the nodes; this was placed on a roof.

Two experiments were performed, the first one where all the typical data as in Table 4.1 was returned ("typical data") for approximately one week and another one for approximately two weeks where the output from the node algorithm as in Table 4.2 in addition to signal strength and some debugging information was returned ("node algorithm data"). The epoch loss as well as the bytes per epoch sent over the radio link where compared when the nodes returned the "node algorithm data" and when the nodes returned "typical data"; the results of which are summarized in the table below.

Table 4.3: Comparison of channel capacity requirement and data loss for typical algorithm compared to CFFS algorithm

The reason for a reduction of only 74% of data compared to the estimated 87.5%; is due to the added signal strength and debugging information added to the "node algorithm data".

There is approximately a 74% reduction in bandwidth and a 46% reduction in epoch loss when using the "node algorithm data" compared to what is normally outputted. As the "node algorithm data" is from observations that have been partly processed through the CFFS algorithm, this data wouldn't be able to be used by typical GNSS packages. The "node algorithm data" is intermediate data designed for the CFFS algorithm, it's an integral part of the algorithm itself. CHAPTER 4. ALGORITHM COMPARISON

Accuracy and precision can be confusing. In particular accuracy is difficult to measure and a subtle subject made worse by common vague usage of the word in everyday speech. The figure below shows the distinction between precision and accuracy for a set of solutions colored blue with respect to a true solution colored red. The black circle signifies the level of precision for the set of measurements while the distance between the center of the black circle in the center of the red dot represents the accuracy of the set of measurements.

Precise solutions but

inaccurate Imprecise solutions but

accurate

Precise solutions and

Imprecise solutions and accurate

inaccurate

Figure 4.1: Precision and accuracy (red dot is true location)

The precision can be easily estimated from a set of solutions, the accuracy is an entirely different matter. As pointed out in [48] , it is not uncommon for accuracy to be calculated internally where solutions are obtained over a, long period such that they take on a normally distributed form and it is assumed that the mean of this distribution matches the true solution. This assumption is not necessarily valid. As we have seen in the multipath chapter 3, multipath causes a fixed solution bias that is present even if observations are taken continuously throughout the day and such a bias does not change at least in simulations from one day to the next. Unavoidably systematic errors CHAPTER 4. ALGORITHM COMPARISON causing biases creep into solutions. For many applications this would not be a problem as a fixed bias that does not change no matter what just means you're consistently measuring to a different point of your apparatus. However, multipath causes a location dependent bias, if the GPS receiver gets moved this bias changes with respect to the mean location of the solutions, hence invalidating the assumption that the mean of normally distributed solutions taken over a large enough period matches the true solution. To measure the true location requires calculating accuracy externally which means using some other device with a better accuracy. This is beyond the scope of this thesis so we only consider approximate internal accuracy assessment. As said by [48] "Keep in mind that any given GPS measurement can be represented by the following equation: measurement— exact value + bias 4- random error. The random error component presents roughly the same problem for both internal and external assessments. The bias however, requires external truth for detection. There is no easy way to detect a constant shift from truth in a dataset by studying only the shifted dataset.".

While probably difficult to measure, the bias caused due to multipath using 24 hours of observations in a low multipath environment with a high quality GPS antenna most likely produces a small bias partly due to the environment and equipment, and also partly due to the 24 hours observations obtaining an average of the instantaneous multipath bias. In this case obtaining accuracy internally is probably reasonably close to the actual accuracy. However, for our application we have poor uncalibrated antennas with no choke rings, a high multipath environment and a predicted yearly cycle of a bias caused clue to multipath and only being able to obtain daylight our observations. This means reliably calculating accuracy internally is far more difficult and potentially more inaccurate. The yearly solution bias cycle as predicted in Chapter 3 particularly makes an accuracy estimate very difficult as at least a year worth of observations would be needed to take into account the yearly solution bias cycle.

Accuracy and precision are meaningless without some type of measure. Accuracy and precision measures are generally the probability that a solution will be within some shape of some given size. Generally the shapes used are circles or spheres and the probabilities go from 50% and up. 2D horizontal (northing easting) seem to be the most commonly given probably because GPS solutions tend to be more accurate in this plane compared to the up direction and are CHAPTER 4. ALGORITHM COMPARISON reasonably symmetrical in this plane. The following table defines three very common measures where CEP DRMS and 2DRMS are for the 2D horizontal plane and SEP is for all three dimensions. For the one dimensional up direction RMS is usually used.

Table 4.4: Measures of accuracy and precision

The following table shows the approximate calculations we used to calculate the these measures where σχ, OE and συ are standard deviations of the solutions in the North East and up directions respectively. These formulas were obtained from [26] ,

Table 4.5: Approximate internal accuracy and precision calculations from a set of solutions

We only consider internal accuracies in this thesis and when we talk about the accuracy of a solution, we mean that we have obtained an estimate for the standard deviations of the population from a sample of solutions in sufficiently large quantity that the mean of these solutions we hope is the same as the true solution and have calculated the accuracies as in the previous table. For solution precision as we typically deal with daily solutions and only able to calculate one solution a day, the precision of the solution means the precision CHAPTER 4. ALGORITHM COMPARISON as measured from the previous table using a set of solutions taken over a period of time no more than a few months or over longer periods of time when the low frequency components of the solution have been removed leaving just the high frequency components.

4.2,1 Comparison methodology

We designed two programs called JAC-MM and J AC. JAC-MM implemented the CFFS algorithm while JAC implemented the CFFS algorithm bar the sidereal filtering stage. If we assume a remote GNSS receiver setup isolated from the power grid and connected to the user via a wireless link, then, if we reduce the cost and power requirements of the hardware and deployment, the user will most likely obtain fewer types of observations, fewer and more intermittent observations, increased muitipath, increased PCV, and generally poorer quality data. The question then remains is it possible to still obtain accurate / precise static solutions and if so how accurate/precise will they likely be. One would not expect current GNSS software to do a particularly good job in this respect as the people designing the software maybe assuming continuous high quality data is available. Therefore, we ran different types of modified observations through JAC and JAC-MM and compared these staic relitive solutions with the same modified observations run through the open source GNSS packages RTKLib [11] and C PS ^' ! k [8]. RTKLib and GPSTk were used because these were the only two open source GNSS packages freely available that supported singleband phase based relative position solutions that we found. Static solutions were obtained on a day by day basis.

RTKLib performs static solutions in addition to a myriad of other solution types such as RTK, Precise Point Positioning PPP), Differential GPS DGPS), Single, Fixed, and so on. Both singleband and dual band receivers are supported. For static solutions on short base lines less than 10 km it uses some form of an extended kalman filter and double differences [52] (pg 161).

GPSTk has in it an application called DDBase that is capable of calculating static solutions using only singleband receivers and also uses double differences

[35] .

JAC implements the CFFS algorithm but without sidereal filtering. Therefore, all daily solutions that JAC calculates is solely based on observations obtained CHAPTER 4. ALGORITHM COMPARISON for that day just as RTKLib and GPSTk are. This allows for meaningful performance comparisons between the three programs.

Sidereal filtering uses information from previous days to calculate solutions and hence is not a fair comparison with JAG GPSTk and RTKLib. Sidereal filtering algorithms could be appended to anyone of these three programs to obtain solutions of higher accuracies / precision. However, we do append an implementation of sidereal filtering to the JAG program so as to see the effect that sidereal filtering has compared to the other three programs; this JAC program with the appended sidereal filtering is called JAC-MM and is a full implementation of our CFFS algorithm as described in Chapter 2.

4.2.2 Continuous high quality data

The purpose of this subsection is to show that our implementation of the CFFS algorithm produces solutions that are comparable to pre-existing solution processing programs under good quality conditions. For initial performance testing we used two months of high quality data obtained from GPS receivers, the first located at Te Papa museum and the second at Wellington Airport. This data was obtained from the RINEX data archive of Land and information New Zealand (LINZ) [5]. The data rate was one epoch every 30 seconds. The baseline between these two receivers was around 4 km. Both receivers were Trimble NetRS [12] while the antennas were different, one being an LEIAT504 and the other a TRM41249.00. The data was supplied to us by the GeoNet project. These types of receivers and antennas are high-end dualband. The static solution accuracy of such a receiver as given by the data sheet is 5 mm - 0.5 ppm DRM S for the horizontal, and the up direction is 5 mm + 1 ppm RMS, where ppm refers to the baseline length. This means the datasheet implies that the static solutions as calculated by the receivers themselves would be 14 mm 1DRM S for the horizontal. Power consumption is between 3 and 4 W.

For a fair comparison, as our algorithm is intended to only deal with singleband GPS receivers, observations were modified first by stripping of precise P codes as well as L2 band observations, so as to be left with just Cl, Dl, SI , and LI observation types and broadcast navigation data; these types of observations are typical for low cost low power singleband GPS receivers. CHAPTER 4. ALGORITHM COMPARISON

DDBase failed on six days due to not being able to read the navigation files correctly despite these files being obtained from the IGS products [27|. Both GPSTk and JAC obtained solutions for all days processed. Figure 4.2 shows a bird's eye view of the solutions while table 4.6 shows some statistics of the solutions where NEU is an acronym for North East Up. The JAC-MM solutions are also included.

WGTN-WGTT day 1 to da 59 of 2010

01.,Ι 81 _> 01

-2035.365 -2035.355 -2035.345 -2035.335

East (m)

Figure 4.2: Birds eye view using high quality data CHAPTER 4. ALGORITHM COMPARISON

Table 4.6: Statistics using high quality data

As can be seen all programs produced solutions with comparable precision. The standard deviations of all solutions are similar. In the horizontal plane DDBase and RTKLib agree with each other on average to within a millimeter, in the vertical RTKLib is 1.6 cm higher than DDBase. JAC agrees with DDBase in the vertical direction to within 8mm while in the horizontal plane JAC differs from both DDBase and RTKLib by approximately 1.8 cm. The JAC to DDBase difference is 20mm, JAC to RTKLib difference is 20mm and the DDBase to RTKLib difference is 16mm. Three measures of accuracy are summarized in the table below.

Table 4.7: Solution precision using high quality data

As already mentioned, JAC-MM computes solutions using data from previous days. Therefore the results JAC-MM's solutions are unfair when directly compared to to the other three programs solutions. However, as expected JAC-MM's solutions are more precise than the other three programs. It's also interesting to note that sidereal filtering improves precision despite the observations only having an epoch rate of once every 30 seconds. This is interesting because the sidereal shift required means there is no mathching observation of a pair of satellites today as there was yesterday, as 246/30 is not an integer. In fact it takes five days of separation for a matching observation of a pair of satellites to occur. Examining the numbers, JAC-MM is about 20% more precise than JAC for this data. CHAPTER 4. ALGORITHM COMPARISON

All programs had significant solution biases when compared with one another ranging from 16 to 20 mm. For our purposes we are not concerned with solution biases that do not change with location or over time. Whether or not these biases caused by using the three different programs change or not with location or extended time we do not know. The figure below shows the mean solutions with 95% confidence spheres over these two months worth of processing. As can be seen none of them are consistent with one another. JAC-MM's sphere was not plotted as this was within JAC's shphere. It's interesting to note that all programs used the same observations but all programs produced different results. All programs produced very precise results.

Mean with 95% confi dence

Figure 4.3: Mean solutions calculated over two months

4,2.3 Simulated intermittent high quality data

Simulated intermittent high quality data was created by using the same high quality Te Papa Wellington Airport data as before but with randomly removing 50% of satellite observations. As before the data was run through J AC, DDBase and RTKLib. DDBase and RTKLib failed to return any solutions. Figure 4.4 shows the bird's eye view of the solutions while Table 4.8 shows CHAPTER 4. ALGORITHM COMPARISON some statistics of the solutions and Table 4.9 shows measures of precision of the solutions.

WGTN-WGTT day 1 to day 59 of 2010

C1,L1,S1,D1, 50% of satellite observations removed

3666.695

3666.690

3666.685

Δ _^

Ά ..Δ

t 3666.680

o

z:

3666.675

3666.670

3666.665

-2035.365 -2035.355 -2035.345 -2035.335

East (m)

Figure 4.4: Birds eye view using high tv data with. 50% random removal of satellite observations

Mean NEU Standard deviation NEU

(m,m,m) (mm, mm, mm)

[ JAC (3666.685 ,-2035.340 ,15.530) (2.51,2.38,4.61)

JAC-MM (3666.685 ,-2035.340 , 15.530 ) (2.21 ,2.14,4.01 )

Table 4.8: Statistics of simulated intermittent high quality data CHAPTER 4. ALGORITHM COMPARISON

Table 4.9: Solution precision using simulated intermittent high quality data

J AG's solutions are almost identical to the situation before where no satellite observations were removed: the mean location is 4 mm difference between the two. As can be seen the effective randomly removing 50% of the satellite observations makes little difference to the precision of JAC's solutions.

4,2.4 Algorithm performance using real life intermittent poor quality data

The purpose of this subsection is to determine whether or not any of the three programs (JAC, DDBase and RTKLib) could obtain solutions of reasonable precision using real life intermittent poor quality data obtained from the GPS WSN nodes that we designed. The same two experiments used in subsection 4.1.1 used to obtain epoch loss were also used to obtain real life intermittent poor quality data.

The CFFS algorithm is designed to be a distributed algorithm where some of the algorithm is run on the nodes themselves. This has the advantage of a reduction in epoch losses as seen in subsection 4.1.1. It would be unfair to compare the CFFS algorithm with DDBase and RTKLib without considering this. If the nodes just return the "typical data" and the node section of the CFFS algorithm was to be run in a non-distributed manner on the base station, the CFFS algorithm could conceivably be incurring a 26% loss of epochs compared to a 14% loss of the epochs as seen in Table 4.3. Therefore, two node setups were configured. The first set up where the "typical data" was returned by the nodes, and another set up where the "node algorithm data" was returned by the nodes. The first set up was designed to supply data to DDBase and RTKLib, whilst the second set up was designed to supply data to the main part of the CFFS algorithm.

Approximately 1 week of "typical data" was obtained for DDBase and RTKLib, while approximately 2 weeks of "node algorithm data" was obtained for the CHAPTER 4. ALGORITHM COMPARISON main section of the CFFS algorithm to process. The observations were not taken simultaneously, and due to the need to reprogram the nodes to output "node algorithm data" rather than "typical data", the location of the nodes where unavoidably shifted slightly during the two data gathering periods. Such movements were minimal and do not change our results. The following figure shows a birds eye view of the solution results. JAC calculated solutions for all days. DDBase could calculate solutions for all days but they were inaccurate, RTKLib did not successfully calculate any solutions.

Node on ground in shade of trees

8/ a /2014 to 15/ ay/2014 (DDBase & RTKLib) and 16/May/2014 to 27/May/2014 (JAC & JAC-MM)

18.900 : : l : :

18.850

□ DDBase

18.800 -> RTKLib

Δ JAC-MM

18.750 o 18.700

18.650

18.600

18.550

-3.000 -2.800 -2.600 -2.400 -2.200 -2.000 -1 .800

East (m)

Figure 4.5: Bird's eye view using intermittent poor quality da

data" for DDBase and RTKLib while "node algorithm data" fo

DDBase solution not plotted due to large variation)

To compare the effect on the ability to obtain solutions due to the reduction in epoch loss when transmitting "node algorithm data" rather than "typical data" between the nodes in the base statio when using JAC, the node algorithm was implemented on the base station and the one week worth of data that DDBase and RTKLib used to calculate solutions above was run through JAC. JAC calculated solutions for all days; the birds eye view of the obtained solution are shown in the following Figure. CHAPTER 4. ALGORITHM COMPARISON

Node on ground in shade trees

-3.000 -2.800 -2.600 -2.400 -2.200 -2.000 -1.800

East (m)

Figure 4.6; Bird's eye view using intermittent poor quality data, "typical data" for DDBase and RTKLib and JAC. (One DDBase solution not plotted due to large variation)

As can be seen despite the increase in loss of epochs, JAC performs very similar to when transmitting "node algorithm data" between the nodes and the base station. The following two tables summarize some statistics of the solutions.

Fable 4.10: Statistics of intermittent poor quality CHAPTER 4. ALGORITHM COMPARISON j CEP (mm) 2D RMS (mm) SEP (mm) j j JAC (average) 5 11 8 j j JAC-MM (average) 3 1 8 j 6 j

DDBase 1800 5400 1700 j

Table 4.11: Solution precision using intermittent poor quality data

As can be seen only JAC was able to obtain solutions of any real value. JAC's measures of precision are approximately double of that of the high quality Te Papa Wellington Airport data. Due to the small sample sizes the sample standard deviations of JAC when transmitting "typical data" compared to "node algorithm data" are different than one would expect. The higher epoch loss incurred when "typical data" was transmitted between the node and the base station had no noticeable effect on being able to obtain valid solutions than when "node algorithm data" was transmitted.

4.2.5 Solution comparison summary

A qualitative summary is pictured in the following figure.

^" Figure 4.7: Qualitative GNSS software solution results CHAPTER 4. ALGORITHM COMPARISON

All three programs obtained comparable precisions whe high quality GPS data over lossless radio links that were rjermanently powered were used, while, at the same time produced significantly different biases with respect to each other ranging from 16 mm to 20 mm. Both DDBase and RTKLib were extremely sensitive to missing satellite observations while JAC was insensitive to such events. Due to JAC's insensitivity to missing satellite observations and small sample sizes, no noticeable improvement in precisio was detected when transmitting "node algorithm data ^" compared to transmitting "typical data". As expected JAC's solutions when using our low end hardware in a high multipath environment were worse than when using high end equipment in a low multipath environment. Only JAC and JAC-MM could produce solutions of any meaningful accuracy when using our GPS WSN nodes. Therefore, as JAC-MM is a implementation of our CFFS algorithm, the CFFS algorithm is suited for our application of land deformation monitoring using our hardware.

While accuracy may not be entirely known, it is possible to measure small sudden movements of just a few millimeters due to the high precision of our system.

The Paekakariki testbed shown in the figure below. Nodes were placed relatively close to one another as depicted by the mauve circle. In addition permanently powered GPS receiver nodes were placed on the roof of the building at the bottom right hand corner of the figure. The left of the figure is the north direction.

Figure 6.1: Paekakariki testbed site CHAPTER 6. TESTBED RESULTS AND ANALYSIS

As can be seen there is poor visibility of the sky the Paekakariki testbed for the nodes on the ground while the receivers on the roof have somewhat better visibility. Observations obtained from this testbed and results using them are used throughout this thesis.

First we look into the number of hours the GPS receivers tend to be active for in a day and how this varies. The figure below shows the number of hours node 272's GPS was active for per day and successfully returned epochs when placed at the Paekakariki test site for approximately two weeks in summer and then in winter.

Summer Winter

0

Day Day

Figure 6.2: Number of hours a day GPS a receiver is active for

In winter we observed a minimum of less than 0.5 hours, while in summer we observed a maximum of more than 8 hours. On average in winter there was 1.7 hours of activity while in summer there was 6.2 hours of activity: this means 3.6 times as many epoch are used for processing summer solutions than winter ones. The horizontal solutions using JAC-MM to implement the GFFS algorithm are pictured in the following figure.

CHAPTER 6. TESTBED RESULTS AND ANALYSIS

Summer winter accuracy comparison

-6 -4 -2 0 2 4 6

East (mm)

Figure 6.3: Qualitative summer winter solution precision comparison

As can be seen there appears to be better precision during summer than winter as one would expect as there are 3.6 times as many epochs used for the summer solutions than the winter ones.

j j 2DRMS (mm)

JAC-MM (Summer) 3

j JAC-MM (Winter) 8

Table 6.1: Solution precision for summer winter comparison test

The picture below shows a deployment at the Paekakariki testbed site run in summer for one month where north is in the up direction of the picture and the nodes are placed directly on the ground; these three nodes can just be seen in figure 6.1. As can be seen in this figure a shadow of a tree is passing the nodes. CHAPTER 6. TESTBED RESULTS AND ANALYSIS

Figure 6.4: Three nodes at the testbed site

The measurements to the approximate physical horizontal center of each of the antennas were measured using a tape measure with an accuracy of approximately 1 cm. The results of these measurements is pictured in the following figure.

Measured with tape measure

■Ac accuracy

Figure 6.5: Tape measure measurements of the nodes to the approximate physical center of the top of the antennas

After one month the daily solutions were processed, and the horizontal solutions all plotted in the figure below; the mean range between each node is also CHAPTER 6. TESTBED RESULTS AND ANALYSIS

Lured in the following figure with 95% uncertainties as derived from izontal components of the CFFS solutions.

1 month Paekakriki static test

3

East (m) figure 6.6: Birds eye view of CFFS solutions for one month Paekakariki static est

As the phase center of each of the three antennas were not known as they were uncalibrated antennas, the difference between the relative distance measurements using the tape measure between the approximate approximate physical centers of the antennas and the difference between the mean values of the horizontal CFFS solutions, one should expect these to only be approximate. However, despite this, the difference in these range values are still consistent with the ones measured using the tape measure.

The average 2DRMS precision calculated for the three nodes for this one month static test was 5 mm and appears typical for solution precision for high mul- tipath environments.

The LuShan site is a known high risk landslide site on a hill above the town of Lushan in Taiwan's mountainous Nantou region. It is a site studied using CHAPTER 6. TESTBED RESULTS AND ANALYSIS satellite imagery, aerial imagery, laser survey and GPS. The active movement zone approximately covers 1 km ² and resembles the shape of an upsidedown V starting at the top of the hill.

A testbed of 4 small batteryless GPS wireless sensor nodes were deployed at the site. The nodes were deployed alongside retro reflectors used for laser based land deformation monitoring using a total station by National Chi Nan University (NCNU).

Due to dense and tall vegetation, both network coverage and adequate visibility of GPS satellites were problematic. For adequate network coverage two batteryless solar powered radio repeaters were used. One placed on the hill of interest while the other on the opposite hill. For adequate visibility of GPS satellites, node sites were found by visual inspection and deciding whether or not there was enough of a clear view of the sky. Node 274 was placed at a knowingly stationary location at the bottom of the hill as a point of reference. The other three nodes or placed on active movement sites.

A permanently powered base station (also known as the router) that functioned as a sink for the data obtained by the nodes, a router to route the data via a 3G modem to us via the Internet for solution processing, and a GPS node itself to obtain satellite observations was placed on the hill called node 1. The base station was mounted in a large pre-existing weatherproof metal cabinet box with mains outlet sockets inside.

The figure below shows the locations of the nodes, repeaters and the base station; the size of the testbed is approximately 300m by 700m.

CHAPTER 6. TESTBED RESULTS AND ANALYSIS

Figure 6.7: LuShan testbed site (Courtesy of NCNIJ)

The direction of land slippage on the hill of interest is in a direction approximately southeast. It is uncertain whether repeater 104 was situated on stationary land or whether or not it too is sliding, as repeater 104 was not on the hill of interest. Node 274 w ^T as the only known stationary location. Node 276 and 277 were on the very west side of the sliding area while node 275 was on the very east of the sliding area. The base station which was permanently powered also functioned as a GPS receiver node called node 1. Node 1 was also on the westside of the sliding area.

LI GPS phase observabies were the principal observations used to calculate relative node positions. These phase observabies along with autonomous solutions calculated by the GPS receiver itself were combined to modify the phase CHAPTER 6. TESTBED RESULTS AND ANALYSIS observables on the nodes before sending them back to the sink using 2.4Ghz 802.15.4 radios either directly or via a series of repeaters. The sink then passed this data onto us via a 3G Internet connection for processing. Despite the lack of a battery the nodes were designed to operate in extremely low light conditions as expected on dark days during heavy rain, the nodes however would not work at night. Because of this the nodes only returned data during the day.

Temperature measurements using a mock node installation at NCNU near the test site in summer were performed and a maximum temperature of 58.1°C was measured within the node. The minimum rated value of any component inside the node was 70°C so the maximum recorded temperature was within component tolerances .

The first node was installed on 4/7/2014 while the last node installed on the

15/7/2014. The mounting method of the nodes can be seen in the following figure where a metal flange is bolted to the side of one of the walls of the road. The node or repeater is then glued to the top using RTV silicone.

Figure 6.8: Node and repeater mounting method at LuShan

Taiwan is prone to typhoons with a season lasting approximately from June to October. It was hoped that heavy rain due to tropical cyclones might cause land deformation that we could measure. Only two tropical cyclones hit Taiwan during the 2014 typhoon season; a category 2 typhoon called Mat mo approximately on 23/7/2014 (day 204) and a tropical storm called Fung-wong approximately on 21/9/2014 (day 264). CHAPTER 6. TESTBED RESULTS AND ANALYSIS

(a) Matmo 22/7/2014 [7] (b) Fimg-wong 20/9/2014 [6|

Figure 6.9: 2014 tropical cyclones for Taiwan

Node failure was a significant issue. Node 274 ultimately failed on the 11/8/2014, node 275 on 13/8/2014, and node 277 on 21 /12/2014. Node 274 and 275 performed normally until they ultimately failed. Node 277 failed during typhoon Matmo and returned no signal for two days after that. Node 277 became increasingly intermittent ultimately the last signal being heard on 20/12/2014 while the second to last signal was heard on the 16/11/2014. Node 276 functioned almost flawlessly and only on three days during tropical cyclone Fung- wong did it fail to transmit any signals.

Both the base station hardware and the repeaters performed flawlessly and not one day was observed where either of them malfunctioned.

The reason for the node failures was presumed to be caused by water entering electronic circuitry and causing corrosion and short circuits. However as we were based in New Zealand we could not easily go back and examine the reason for these node failures. The figure below shows node 275 some time after it failed. Originally the GPS antenna was a shiny silver colour while clearly in the figure below it is been tarnished. The antenna is under RTV silicone which was designed to protect it against the elements. A small amount of brown colour at the bottom right corner of the antenna is also presumably corrosion. Therefore it seems likely that this particular node failed due to water damage. CHAPTER 6. TESTBED RESULTS AND ANALYSIS

Figure 6.10: Suspected water related corrosion in node 275 (Circa 28/8/2014. Courtesy of NCNU)

Figure 6.11 : Vegetative overgrowth on node 277 (Circa 28/8/2014, Courtesy of NCNU)

Previous attempts by ourselves to create nodes impervious to water damage proved difficult. At the Paekakariki testbed after six months outside the following picture was obtained upon opening up the GPS receiver node. As can be seen considerable amount of water has entered the housing and has also managed to get in past the conformal coating that covers the PCB. In fact, because of the conformal coating on the PCB the water can't drain away making the situation worse. CHAPTER 6. TESTBED RESULTS AND ANALYSIS

Figure 6.12: Corrosion at the Paekakariki testbed

At the LuShan testbed at the beginning of 2014, we sent initial GPS receiver nodes to be tested on the hill. By the time we arrived in late June, upon opening one GPS receiver node up, the housing contained hot water that would not drain away as can be seen in the following figure, figure 6.13. The housing of the new nodes we put on the hill in July had bee redesigned and the PC boards inside the nodes themselves were entirely covered in epoxy resin except for the programming pins. The GPS antenna was covered in epoxy. This we hoped would prevent node failure due to water damage. However, we did not have enough time to rigorously do testing beforehand.

Figure 6.13: Water inside an initial testing node at the LuShan site

We were wanting to use node 274 as a fixed reference so as to be able to measure the relative locations to the other nodes. However, as node 274 failed approximately a month after installation this could not be done for most of the time testbed was running. This left for the majority of the time only node 1 and 276 usable. Node 1 was mounted on a metal I-beam that once upon CHAPTER 6. TESTBED RESULTS AND ANALYSIS a time presumably was vertical but due to land deformation was on a slight south tilted slope. Both these nodes were on the active sliding zone meaning over the long term all we could hope to measure as the long term relationship between a wall along a road and the top of an I-beam.

6.2,1 Movements relative to node 274

In the figures that follow we show the NEU time series of the wireless solar powered GPS receivers with respect to node 274. This was done by taking the difference of the final solutions between each node with respect to node 274. When taking the difference between solutions additive noise increases. However, it seems that there is a considerable stochastic dependence between solutions obtained from two distinct nodes with respect to a common reference station, presumably as the similar sunlight conditions causes observations to be taken reasonably similarly between the two nodes. Therefore the precision of the solutions obtain this way is not as bad as one might first anticipate. This method of calculating solutions has the advantage of being able to use a permanently powered reference station that is on an active sliding zone as a reference station such to maximise the probability that there exists an epoch for an observation on a given solar powered GPS receiver, hence more double differences can be used in the final solution. This contrasts to the method by using the two wireless GPS receiver nodes themselves to obtain double differences and hence solve for position; in this case the two GPS receiver nodes must be on at the same time to obtain a double difference which does not necessarily happen. In addition it is far easier subtracting final solutions than creating double difference is between the two GPS receivers themselves.

For the figures that follow, we plot the NEU solutions produced with JAC-MM which implements the CFFS algorithm as described in chapter 2. The time series of all the wireless GPS receiver nodes with respect to node 274 were plotted as we were told that no 274 was on a non- sliding zone. As node 274 failed shortly after installation only approximately one month of solutions are plotted on each of the figures. A black vertical line just before the 204th day as to signify the approximate time observations would be affected by typhoon Matmo. In the 277-274 NEU timeseries plot two outliers were removed by hand as they were significantly inconsistent with the rest of the solutions. CHAPTER 6. TESTBED RESULTS AND ANALYSIS

Ignoring the first 14 solutions to allow the sidereal filtering to stabilise the average 2DRMS precision was 3.3 mm using this indirect way of obtaining relative solutions relative to node 274.

275-274 NEU time series

Dav

Figure 6.14: 27-5-274 NEU timeseries

276-274 NEU time series

415.945

f! 415.940

415.935

415.930

E _. 95.935

in 95.930

210.22

210.20

195 200 205 210 215

Day

Figure 6.15: 276-274 NEU timeserie CHAPTER 6. TESTBED RESULTS AND ANALYSIS

277-274 NEU time series

Day

Figure 6.16: 276-274 NEU timeseries

As can be seen there is no noticeable movement caused due to typhoon Matmo. Any solution aberrations around the few days after the typhoon may be due to the lack of double differences because of cloudy days or swelling of the soil or some other temporary phenomena. All nodes had a southward trend between 0.12 mm/day to 0.22 mm/day. Node 277 had a slightly eastward component while 275 and 276 had very slightly westward components all nodes had an up component of between 0.3 mm per day to 0.7 mm per day. There is no compelling evidence that this very slight trend over such a short period of time is related to any land deformation and is entirely conceivable that it is instead the slow bias variation predicted in simulations in chapter 3 caused due to the four minute sidereal shift with respect to the solar day or the nodes settling over time.

6,2.2 277 and 276 relative to node 1

After day 222 node 274 was never heard of again. Node 277 was heard of less and less, while node 276 still continues today as of writing. Therefore at this point we compare solutions relative to node 1 and in particular concentrate on 276 relative to node 1 as this pair was most consistently active. CHAPTER 6. TESTBED RESULTS AND ANALYSIS

The figure below shows the NEU timeseries solutions for node 277 relative to node 1 (277-1) . The second black vertical line centered just before day 264 is to signify the time tropical storm Fung-wong might affect solutions.

0

Day

Figure 6.17: 277-1 NEU timeseries

As can be seen a significant anomaly occurs centered around date 230 and extends at least 10 days to either side of the date. There is a drop in the north direction of about 5 mm while there is a rise in the east direction of approximately the same magnitude at the same time; what this as we do not know. As can be seen there is no significant change around the day tropical storm Fung-wong hit. As with all the solutions shown so far for LuShan. there is a slow approximately linear increase in the up direction; in this case about 0.03 mm/day to 0.05 mm/day.

The following figure shows north east up time series solutions for 276-1 obtained using JAC-MM. CHAPTER 6. TEST BED RESULTS AND ANALYSIS

276- 1 NEU time series

Day

Figure 6.18: 276-1 NEU timeseries. 2D RMS precision 3.6mm. 2DRMS internal accuracy between 3.6 to 9.12mm

Between approximately days 340 to 360 there was an unusual increase in noise, this was during the month of December. We believe this noise to be as a result of roadworks, as during this time some of the asphalt of road was being replaced, the heavy machinery passing by the sensors presumably vibrating them. However, apart from the month with the roadworks, generally the solutions were very precise and slowly changed over the seven months they were obtained. As this was our test that lasted the greatest period of time, we considered seven months close enough to 12 months that we could at least obtain an idea of the expected accuracy a system such as ours would have in the typical setting. If node 276 did not move and/or its surroundings did not change then ignoring the first 14 days to allow the sidereal filtering to settle the internal solution accuracy was 9.12 mm 2DRMS; this of course is only valid if there was no movement or the surroundings did not change. This can also be considered approximately as a worst case scenario for internal accuracy. If on the other hand the slow change of the solutions where caused by physical movement of the node itself, then using 10 dgree polynomial fits to remove the CHAPTER 6. TESTBED RESULTS AND ANALYSIS slow change of the solutions, a 2DRMS internal accuracy of approximately 3.6 mm was obtained; this can be considered the best case for what the solution accuracy was and also considered as the solution precision.

First off we note that the up direction unlike node 277 in not linear, the east solutions undulate smoothly while the north solutions are fairly flat until day 310 where they rapidly drop while at the same time the up direction rapidly increases. The most significant thing we noticed was that the up solutions obtained just prior to tropical storm Fung-wong where significantly higher than the ones obtain just after the storm. There was also a reasonably significant quick change in the north solutions at the same time. In addition the rate of change of the solutions changed abruptly before and after the storm. These rapid changes can be seen better if we obtain the primary principal component direction vector between days 242 and 287 inclusive and plot these solutions as can be seen in the following figure.

Day

Figure 6.19: Primary principal component of 276-1 during tropical cyclone Fung- Wong

This figure shows approximately 20 days before the cyclone and 20 days after the cyclone. The green and the mauve lines are linear fits for the time before CHAPTER 6. TESTBED RESULTS AND ANALYSIS and after the cyclone respectively. The solutions are changing at a different rate before the cyclone and after the cyclone. In addition the two linear fits mismatch at the time of the cyclone by 7.45 mm. The solutions also appear to be more or less normally distributed around their linear fits. We then estimate the accuracy of this displacement by removing the linear trends from these solutions to obtain the following figure with standard deviation of dotted lines.

Day

Figure 6.20: Primary principal component of 276-1 during tropical cyclone Fung- Wong with linear trends removed

As can be seen it appears to be reasonably normally distributed with a standard deviation of 1.758 mm. Therefore we should obtain a 95% confidence that the movement measured was ±4 x 1.758/ y 20 mm— 1.6 mm around our estimate. In other words 7.5 ± 1.6 mm is our estimate for the change in the primary principal component before and after the cyclone. From chapter 3, the slow solution bias change due to the four minute daily shift, we did not see that in our simulations there should be a discontinuity of the solution derivative or a discontinuity in the solutions themselves. As we expect from simulations that the solution bias change should have a repeating period of one year we expect over a short period such as 20 days that the solution bias should not alter significantly over such a short period of time; this justifies what we have done by using linear fits on both sides of the cyclone. In addition it seems that a linear fit does indeed fit very well as can be seen in the previous figure. While it appears that there are other discontinuities of the solution derivative in figure 6.18 that may or may not be related to movement, that this disconti- CHAPTER 6. TESTBED RESULTS AND ANALYSIS nuity at day 264 appears exactly when cyclone Fung- Wong hit, it seems more than probable that this discontinuity was caused either due to one of the nodes moving or the surroundings of the nodes changing such as trees falling down changing the multipath environment. Therefore, ignoring the possibility that the nodes surroundings changed we say that it's quite possible that one of the nodes did indeed move by 7.5 ± 1.6 mm. Node 276 relative to node 1 moved in a direction of travel near vertical with an elevation of approximately -75° and a small component in the horizontal plane with an azimuth of approximately 280°; that is to say either node 276 moved down or node 1 moved up.

The discontinuities in the solution derivative are are not entirely understood by us as they did not appear in our simulations. One significant discontinuity in figure 6.18 is in the up direction around day 315 which is also accompanied by a slow drop on the north direction. This just happens to coincide with a small shallow earthquake of 3.6 magnitude and 4.9 km in depth very close by but could just be coincidence.

6.2.3 Comparison with retro reflector observations

The following figure shows the locations of the retro reflectors at the site. The green coloured retro reflectors are ones that are outside the active slipping zone and are used for the relative positioning of the reflectors in the slipping zone.

CHAPTER 6. TESTBED RESULTS AND ANALYSIS

Figure 6.21: Retro reflector locations (Courtesy of NCNU)

As can be seen there is no retro reflector fitted at the site of our base station's GPS receiver node 1. Retro reflector 8 was situated next to our node 276 while retro reflector 9 was next to our node 274. The total station used by NCNU was a Trimble S8; according to NCNU, the report they received from Trimble said that range accuracy was 2mm ± Ippm and an angle accuracy of 0.5", though we were not informed of what type of accuracies these were. From the Trimble S8 datasheet [56] there are many different accuracies given for different situations; where retroreflectors were used all accuracies were around the same order and are given as RMSE for range and standard deviation for angle. For performance under "DR HP" a range accuracy lmm + Vppm RMSE and an angle accuracy of 0.5" standard deviation. Therefore we assume a range accuracy 2mm + lppro RMSE and an angle accuracy of 0.5" standard deviation which approximately matches the "DR HP" performance as given in the datasheet and the numbers are very similar to the ones given to us by NCNU. The total station was placed approximately a kilometer due south of the site; the angle accuracy and the distance accuracy along with orientation CHAPTER 6. TESTBED RESULTS AND ANALYSIS can be seen in the following figure.

Fisnire 6.22: Trimble S8 orientation and accuracies

To obtain an idea of what 3D accuracy should be expected from a total station, a very approximate estimate was derived. For the east and up accuracies we used the small angle approximation assuming rotational symmetry around the north axis, while that of the north accuracy simply to be the distance accuracy of the total station. Doing this means the standard deviation in the east and up directions is about 1000 | ₍ ~ [ ^ ] ^¾ 2mm while that of the north direction is about 2 x 10 ^"" + J .LO,,O _r .O,LOO„O,, 1000 = 3mm. Therefore, ' a 95% confidence interval produces east or up ±4mm and north ±&mm and a 2DRMS accuracy of 7mm.

The following three figures show the total displacement since δ/Νον/2013 of the retro reflectors with respect to retro reflector 2 and 9 to the date as indicated by the colour. For example, between 2/7/2014 and 28/8/2014 retro reflector 6 moved 10 mm south in 57 days.

We are interested in retro reflector 8 as this corresponds to our node 276, and we are interested in retroreflectors 3 and 4 as these are close to our node 1 and hence we would expect that movements of the land are more likely to be similar than other retro reflectors. CHAPTER 6. TESTBED RESULTS AND ANALYSIS

CHAPTER 6. TESTBED RESULTS AND ANALYSIS

Up totai s¾sis¾is

Figure 6.25; Up total station measurements (Courtesy of NCNU)

As can be seen there are some considerable changes in the solutions obtained using total station; both third up solutions and the fourth east solutions are considerably out from what one would expect. There also appears to be considerable correlation of the solutions for each time the total station was used to obtain measurements as can be seen when one solution moves south all the other solutions tend to move at the same time. The sparseness of the measurements are problematic as this makes it hard to know how reliable the measurements are. Of all the solutions the north solutions seem to be the most consistent and shows a slow general trend south as expected. Both the east and up solutions do not seem to say a great deal even after removing the measurements performed in July and August. Clearly there are other errors introduced into these solutions that are considerably greater than the ones we estimated. Due to these other errors how the actual movement of the hill relates to total station measurements is unknown. While unsatisfactory we really can't compare the GPS solutions with the total station solutions.

The total station measurements are time sparse and variable making conclusions difficult. The GPS solutions for instantaneous movements (ie. ones that are not slow and continuous) are clearly very precise as we saw for the solution change due to tropical storm Fung- Wong. In contrast it is difficult to make conclusions about the long-term movement of the GPS receivers over CHAPTER 6. TESTBED RESULTS AND ANALYSIS the six months from the GPS solutions. This is because from Multipath and sidereal filtering we predict a solution component due to the multipath bias that is predicted to change over the year making it difficult to distinguish true movement and merely a slowly changing solution bias. In addition to GPS movements can only be compared for for the month of July before node 274 failed. However, in this month the GPS receivers all showed a, slight southward trend relative to node 274 which is given as being stationary, this matches the general trend of that of the north solutions obtained from the total station and what one would expect.