"Applied Science, Faculty of"@en . "Electrical and Computer Engineering, Department of"@en . "DSpace"@en . "UBCV"@en . "Roee, Diamant"@en . "2013-08-27T17:16:21Z"@en . "2013"@en . "Doctor of Philosophy - PhD"@en . "University of British Columbia"@en . "Ocean exploration, through the development of ocean-observation systems, is a key step towards a fuller understanding of life on Earth. Underwater acoustic communication networks (UWANs) will help to fulfill the needs of these ocean-observation systems, whose applications include gathering of scientific data, early warning systems, ecosystem monitoring and military surveillance. The data derived from UWANs is typically interpreted with reference to the location of a data collecting node, e.g. when reporting an event occurrence, or the location of an object itself is of interest, e.g. when tracking a moving underwater vehicle or diver. In this dissertation, we develop methods for localization and efficient data exchange in UWANs.\nIn the first part of this work, we focus on underwater localization (UWL). Since global positioning system signals do not propagate through water, UWL is often based on fusing information from acceleration-based sensors and ranging information to anchor nodes with known locations. We consider practical challenges of UWL. The propagation speed varies with depth and location, anchor and unlocalized nodes are not time-synchronized, nodes are moving due to ocean currents, propagation delay measurements for ranging of non-line-of-sight communication links are mistakenly identified as line-of-sight, and unpredictable changes in the ocean current makes it hard to determine motion models for tracking. Taking these features of UWL into account, we propose localization and tracking schemes that exploit the spatially correlated ocean current, nodes' constant motion, and the periodicity of ocean waves.\nIn the second part of this thesis, we use location information to develop medium access control scheduling algorithms and channel coding schemes. We focus on adaptive scheduling in which each node transmits based on timely network information. Specifically, our scheduling algorithms utilize the long propagation delay in the channel and the sparsity of the network topology to improve throughput, reliability and robustness to topology changes.\nTo evaluate performance, we have developed a simulator combining existing numerical models of ocean current and of power attenuation in the ocean. We have also verified simulation results in four sea experiments of different channel bathymetry structures, using both industry and self-developed underwater acoustic modems."@en . "https://circle.library.ubc.ca/rest/handle/2429/44893?expand=metadata"@en . "Spatial Reuse Scheduling andLocalization for UnderwaterAcoustic CommunicationNetworksbyRoee DiamantB.Sc., The Technion, Israel Institute of Technology, 2002M.Sc., The Technion, Israel Institute of Technology, 2007A THESIS SUBMITTED IN PARTIAL FULFILMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDoctor of PhilosophyinThe Faculty of Graduate and Postdoctoral Studies(Electrical and Computer Engineering)The University Of British Columbia(Vancouver)August 2013c? Roee Diamant 2013AbstractOcean exploration, through the development of ocean-observation systems, is a keystep towards a fuller understanding of life on Earth. Underwater acoustic commu-nication networks (UWANs) will help to fulfill the needs of these ocean-observationsystems, whose applications include gathering of scientific data, early warning sys-tems, ecosystem monitoring and military surveillance. The data derived from UWANsis typically interpreted with reference to the location of a data collecting node, e.g.when reporting an event occurrence, or the location of an object itself is of interest,e.g. when tracking a moving underwater vehicle or diver. In this dissertation, wedevelop methods for localization and efficient data exchange in UWANs.In the first part of this work, we focus on underwater localization (UWL). Sinceglobal positioning system signals do not propagate through water, UWL is oftenbased on fusing information from acceleration-based sensors and ranging informationto anchor nodes with known locations. We consider practical challenges of UWL. Thepropagation speed varies with depth and location, anchor and unlocalized nodes arenot time-synchronized, nodes are moving due to ocean currents, propagation delaymeasurements for ranging of non-line-of-sight communication links are mistakenlyidentified as line-of-sight, and unpredictable changes in the ocean current makesit hard to determine motion models for tracking. Taking these features of UWLinto account, we propose localization and tracking schemes that exploit the spatiallycorrelated ocean current, nodes? constant motion, and the periodicity of ocean waves.In the second part of this thesis, we use location information to develop mediumaccess control scheduling algorithms and channel coding schemes. We focus on adap-tive scheduling in which each node transmits based on timely network information.Specifically, our scheduling algorithms utilize the long propagation delay in the chan-nel and the sparsity of the network topology to improve throughput, reliability androbustness to topology changes. To evaluate performance, we have developed a simu-lator combining existing numerical models of ocean current and of power attenuationin the ocean. We have also verified simulation results in four sea experiments ofdifferent channel bathymetry structures, using both industry and self-developed un-derwater acoustic modems.iiPrefaceHereby, I declare that I am the first author of this thesis. The following publicationshave resulted from the thesis research.Journal Papers1. Roee Diamant, Ghasem Naddafzadeh Shirazi, and Lutz Lampe, ?Robust Spa-tial Reuse Scheduling in Underwater Acoustic Communication Networks,? IEEEJournal of Oceanic Engineering, Dec. 2012 (Included in Part II, Chapter 7).2. Roee Diamant, Hwee-Pink Tan, and Lutz Lampe, ?LOS and NLOS Classifi-cation for Underwater Acoustic Localization,? IEEE Transactions of MobileComputing, Nov. 2012 (Included in Part I, Chapter 4).3. Roee Diamant, Wenbo Shi, Wee-Seng Soh, and Lutz Lampe, ?Time and Spa-tial Reuse Handshake Protocol for Underwater Acoustic Communication Net-works,? IEEE Journal of Oceanic Engineering, Oct. 2012 (Included in Part II,Chapter 6).4. Roee Diamant and Lutz Lampe, ?Underwater Localization with Time Synchro-nization and Propagation Speed Uncertainties,? IEEE Transactions on MobileComputing, Oct. 2012 (Included in Part I, Chapter 2).5. Hwee-Pink Tan, Roee Diamant, Marc Waldmeyer, and Winston K.G. Seah, ?ASurvey of Techniques and Challenges in Underwater Localization,? ACM Jour-nal of Ocean Engineering,, vol. 38, no. 14, pp. 1663-1676, Oct. 2011 (Includedin Chapter 1).6. Roee Diamant and Lutz Lampe, ?Spatial Reuse TDMA for Broadcast Ad-HocUnderwater Acoustic Communication Networks,? IEEE Journal of Oceanic En-gineering,, vol. 36, no. 2, pp. 172-185, Apr. 2011 (Included in Part II, Chap-ter 7).7. Roee Diamant and Yunye Jin, ?AMachine Learning Approach for Dead Reckon-ing Navigation at Sea Using a Single Accelerometer,? IEEE Journal of OceanicEngineering, Apr. 2013 (Included in Part I, Chapter 5).iiiPreface8. Roee Diamant and Lutz Lampe, ?Adaptive Error-Correction Coding Schemefor Underwater Acoustic Communication Networks,? submitted (Included inPart II, Chapter 8).9. Roee Diamant, Lars Michael Wolff, and Lutz Lampe, ?Tracking for UnderwaterNavigation,? submitted (Included in Part I, Chapter 3).Conference Papers1. Roee Diamant, Lars Michael Wolff, and Lutz Lampe, ?Utilizing Ocean CurrentSpatial Correlation for Velocity Estimation of Underwater Drifting Nodes,? inProc. of IEEE Workshop on Advances in Network Localization and Navigation(ANLN), in ICC, Budapest, Hungary, Jun. 2013 (Included in Part I, Chap-ter 3).2. Lars Michael Wolff, Roee Diamant, and Lutz Lampe, ?Spatial and TemporalDependencies of Velocities of Underwater Drifting Nodes,? (extended abstract)in Proc. of ACM Conference on Underwater Networks and System (WUWNet),Los Angeles, USA, Nov. 2012 (Included in Part I, Chapter 3).3. Roee Diamant, Wenbo Shi, Wee-Seng Soh, and Lutz Lampe, ?Joint Time andSpatial Reuse Handshake Protocol for Underwater Acoustic CommunicationNetworks,? in Proc. of IEEE OCEANS Conference, Kona, Hawaii, Sep. 2011(Included in Part II, Chapter 6).4. Roee Diamant, Ghasem Naddafzadeh Shirazi, and Lutz Lampe, ?Robust Spa-tial Reuse Scheduling in Underwater Acoustic Communication Networks,? inProc. of IEEE Vehicular Technology Conference (VTC), San Francisco, USA,Sep. 2011, (Included in Part II, Chapter 7).5. Roee Diamant and Lutz Lampe, ?Underwater Localization with Time Synchro-nization and Propagation Speed Uncertainties,? in Proc. of IEEE Workshopon Positioning, Navigation and Communication (WPNC), Dresden, Germany,Apr. 2011 (Included in Part I, Chapter 2).6. Roee Diamant, Hwee-Pink Tan, and Lutz Lampe, ?NLOS Identification Using aHybrid ToA-Signal Strength Algorithm for Underwater Acoustic Localization,?in Proc. of IEEE OCEANS Conference, Seattle, USA, Sep. 2010 (Included inPart I, Chapter 4).7. Roee Diamant and Lutz Lampe, ?A Hybrid Spatial Reuse MAC Protocol forAd-Hoc Underwater Acoustic Communication Networks,? invited paper in Proc.of IEEE Workshop on Acoustic Underwater Networks, in ICC, Cape Town,South Africa, May. 2010 (Included in Part II, Chapter 7).ivPrefaceUnless stated differently, for all publications, I conducted the survey on related topics,identified the challenges, formalized the suggested solution, performed the analysis,and carried out all of the simulations and sea experiments. I also wrote all paperdrafts. My supervisor, Prof. Lutz Lampe, guided my research, validated analysisand methodology, and edited the manuscripts for papers co-authored by him. Partsof the thesis are a result of research collaboration with additional contributors. Theco-authors? contribution is listed below.1. Journal paper 1 and Conference paper 4: Ghasem Naddafzadeh Shirazi wroteroughly 10% of the simulation code and helped with editing the manuscript.2. Journal paper 2 and Conference paper 7: Hwee-Pink Tan organized the Singa-pore sea trial, and helped with editing the manuscript.3. Journal paper 3 and Conference paper 3: Wenbo Shi wrote the simulation codefor the slotted-FAMA benchmark method (roughly 5% of the code); Wee-SengSoh helped with editing the manuscript.4. Journal paper 5: For this survey publication I was a co-author. Hwee-Pink Tanwrote the paper drafts, and conducted half of the survey. Marc Waldmeyer andWinston K.G. Seah helped with editing the manuscript. I conducted and wrotehalf of the survey and identified the research challenges. The material includedin the thesis comprises only my contribution from this work.5. Journal paper 7: Yunye Jin organized the sea trial, and helped with editing themanuscript.6. Journal paper 9 and Conference papers 1 and 2: Lars Michael Wolff was co-supervised by me during this work. He wrote roughly 40% of the simulationcode, and formalized roughly 20% of the suggested solution.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiList of Abbreviations and Symbols . . . . . . . . . . . . . . . . . . . . . xvAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xixDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xx1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Underwater Localization and Tracking . . . . . . . . . . . . . . . . . 31.1.1 Underwater Localization . . . . . . . . . . . . . . . . . . . . . 41.1.2 Tracking the Location of Nodes . . . . . . . . . . . . . . . . . 91.1.3 Considerations of Measurement Errors . . . . . . . . . . . . . 111.2 Spatial Reuse Scheduling in UWANs . . . . . . . . . . . . . . . . . . 131.2.1 Scheduling Algorithms for UWANs . . . . . . . . . . . . . . . 131.2.2 Adaptive Transmissions . . . . . . . . . . . . . . . . . . . . . 171.3 Open Problems Addressed in this Thesis . . . . . . . . . . . . . . . . 171.3.1 Challenges for Underwater Localization . . . . . . . . . . . . 171.3.2 Challenges for Scheduling of UWANs . . . . . . . . . . . . . . 191.4 Thesis Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19I Underwater Localization and Tracking . . . . . . . . . . . 212 UWL with Time-Synchronization and Propagation Speed Uncer-tainties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222.1 Intuition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 232.2 System Setup and Assumptions . . . . . . . . . . . . . . . . . . . . . 23viTable of Contents2.3 The STSL Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . 252.3.1 Step 1: Time-Synchronization . . . . . . . . . . . . . . . . . . 262.3.2 Step 2: Localization . . . . . . . . . . . . . . . . . . . . . . . 282.3.3 Extensions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 302.3.4 Scalability . . . . . . . . . . . . . . . . . . . . . . . . . . . . 312.3.5 Pseudo-Code for STSL . . . . . . . . . . . . . . . . . . . . . 322.4 Crame?r-Rao Lower Bound . . . . . . . . . . . . . . . . . . . . . . . . 322.4.1 General Crame?r-Rao Lower Bound . . . . . . . . . . . . . . . 322.4.2 Application to STSL . . . . . . . . . . . . . . . . . . . . . . . 342.5 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 342.6 Sea Trial Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 392.6.1 Channel and System Characteristics . . . . . . . . . . . . . . 402.6.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 412.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 433 Spatially Dependent Underwater Navigation . . . . . . . . . . . . . 453.1 System Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 463.2 The SSM and Measurement Model . . . . . . . . . . . . . . . . . . . 473.2.1 State Space Model (SSM) . . . . . . . . . . . . . . . . . . . . 483.2.2 Measurement Model . . . . . . . . . . . . . . . . . . . . . . . 483.2.3 Crame?r Rao Lower Bound (CRLB) . . . . . . . . . . . . . . . 513.3 The DD-UT Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . 523.3.1 Tracking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 523.3.2 Drift Velocity Estimation . . . . . . . . . . . . . . . . . . . . 533.3.3 Confidence Index (CI) . . . . . . . . . . . . . . . . . . . . . . 553.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 573.4.1 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 583.4.2 Sea Trials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 633.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 664 LOS and NLOS Classification for UWL . . . . . . . . . . . . . . . . 674.1 System Setup and Assumptions . . . . . . . . . . . . . . . . . . . . . 684.1.1 RSS-Based Range Measurements . . . . . . . . . . . . . . . . 694.1.2 PDF for PD Measurements . . . . . . . . . . . . . . . . . . . 704.1.3 Remark on Algorithm Structure . . . . . . . . . . . . . . . . 714.2 Step One: Identifying ONLOS Links . . . . . . . . . . . . . . . . . . 714.2.1 Classification of non-ONLOS Links . . . . . . . . . . . . . . . 724.2.2 Classification of ONLOS Links . . . . . . . . . . . . . . . . . 734.3 Step 2: Classifying LOS and SNLOS Links . . . . . . . . . . . . . . 734.3.1 Basic Idea . . . . . . . . . . . . . . . . . . . . . . . . . . . . 734.3.2 Equivalence Constraints . . . . . . . . . . . . . . . . . . . . . 744.3.3 Formalizing the Log-Likelihood Function . . . . . . . . . . . . 74viiTable of Contents4.3.4 Estimating the Distribution Parameters ?1 and ?2 . . . . . . 754.3.5 Forming Initial Estimation ?0 . . . . . . . . . . . . . . . . . . 764.3.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 764.3.7 Summarizing the Operation of the Classifier . . . . . . . . . . 764.3.8 Deriving the HCRLB . . . . . . . . . . . . . . . . . . . . . . 784.4 Performance Evaluation . . . . . . . . . . . . . . . . . . . . . . . . . 794.4.1 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 794.4.2 Sea Trials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 824.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 875 A Machine Learning Approach for Underwater DR . . . . . . . . 895.1 System Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 905.2 The DR-A Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 935.2.1 Forming A?h: Estimation of the Heading Angle . . . . . . . . 935.2.2 Forming A?h,p: Estimating the Time-varying Pitch Angle . . . 945.2.3 Distance Estimation . . . . . . . . . . . . . . . . . . . . . . . 1005.2.4 Summarizing the Operation of the DR-A Method . . . . . . . 1005.2.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1025.3 Performance Evaluation . . . . . . . . . . . . . . . . . . . . . . . . . 1035.3.1 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . 1035.3.2 Sea-Trial Results . . . . . . . . . . . . . . . . . . . . . . . . . 1075.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110II Spatial Reuse MAC techniques for UWANs . . . . . . 1126 Time and Spatial Reuse Handshake Protocol for UWANs . . . . 1136.1 System Model and Objectives . . . . . . . . . . . . . . . . . . . . . . 1136.1.1 MAC Throughput . . . . . . . . . . . . . . . . . . . . . . . . 1156.1.2 Fairness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1166.1.3 Scheduling Delay . . . . . . . . . . . . . . . . . . . . . . . . . 1166.2 Maximizing Channel Utilization in Handshake Protocols . . . . . . . 1166.2.1 Types of Primary Conflicts . . . . . . . . . . . . . . . . . . . 1186.2.2 Formalizing Constraints . . . . . . . . . . . . . . . . . . . . . 1186.2.3 Channel Utilization Maximization Problem . . . . . . . . . . 1226.3 The TSR Protocol - A Sub-Optimal Approach . . . . . . . . . . . . 1236.3.1 Priority of Control Packets . . . . . . . . . . . . . . . . . . . 1266.3.2 Scheduling Control Packets . . . . . . . . . . . . . . . . . . . 1276.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1316.4.1 Simulation Setting . . . . . . . . . . . . . . . . . . . . . . . . 1336.4.2 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . 1346.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138viiiTable of Contents7 Robust Spatial Reuse Scheduling in UWANs . . . . . . . . . . . . . 1397.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1407.1.1 System Model . . . . . . . . . . . . . . . . . . . . . . . . . . 1407.1.2 Objectives of Resource Allocation . . . . . . . . . . . . . . . 1417.2 Formalizing the BSP . . . . . . . . . . . . . . . . . . . . . . . . . . . 1437.2.1 Basic Approach . . . . . . . . . . . . . . . . . . . . . . . . . 1437.2.2 Formalizing the T-BSP . . . . . . . . . . . . . . . . . . . . . 1447.2.3 Formalizing the R-BSP . . . . . . . . . . . . . . . . . . . . . 1467.3 Obtaining the Topology Matrix . . . . . . . . . . . . . . . . . . . . . 1547.3.1 Obtaining the Topology Matrix . . . . . . . . . . . . . . . . . 1547.3.2 Constructing the Conflict Graph . . . . . . . . . . . . . . . . 1567.4 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1577.4.1 Model-Based Topology . . . . . . . . . . . . . . . . . . . . . 1587.4.2 Sea-Trial-Based Topology . . . . . . . . . . . . . . . . . . . . 1627.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1658 Adaptive Error-Correction Coding Scheme for UWANs . . . . . . 1678.1 System Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1688.2 The Adaptive Coding Scheme . . . . . . . . . . . . . . . . . . . . . . 1688.2.1 Gain of Adaptive Coding . . . . . . . . . . . . . . . . . . . . 1698.2.2 Optimization of Nactualt,i . . . . . . . . . . . . . . . . . . . . . 1728.2.3 Extension to IR-HARQ . . . . . . . . . . . . . . . . . . . . . 1738.3 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1748.4 Performance Results and Discussion . . . . . . . . . . . . . . . . . . 1778.4.1 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1778.4.2 Sea Trial Results . . . . . . . . . . . . . . . . . . . . . . . . . 1848.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 185III Summary of Thesis and Further Research . . . . . . . 187Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .A Alternating Optimization Approach for Solving (4.21) . . . . . . . 206B Expressions for the HCRB . . . . . . . . . . . . . . . . . . . . . . . . 208ixList of Tables3.1 Drift velocity estimation results from the Israel and Singapore sea trials. 644.1 Israel sea trial: distribution of estimated values of ?Plastm , m = 1, 2. . . 844.2 Singapore sea trial: cx??E(di)E(di) . . . . . . . . . . . . . . . . . . . . . . . . 86xList of Figures1.1 Communications architecture for UWSNs [1]. . . . . . . . . . . . . . . 21.2 Mapping between the challenges and desirable properties of underwa-ter localization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51.3 Illustration of AUV-aided localization. . . . . . . . . . . . . . . . . . 71.4 Illustration of various types of communication links: LOS, SNLOS andONLOS links. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121.5 Illustration of the vessel?s wave-induced motion. . . . . . . . . . . . . 131.6 Structure of thesis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202.1 ?err from (2.31) as a function of 1/?2. Sound speed is known andall nodes are time-synchronized. Vertical bars show 95% confidenceintervals of the simulation results for STSL. . . . . . . . . . . . . . . 352.2 ?err from (2.31) as a function of esync. Sound speed is known and 1?2 =46 dB. Vertical bars show 95% confidence intervals of the simulationresults for STSL. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 362.3 ?err from (2.31) as a function of |c? c?|. All nodes are time-synchronizedand 1?2 = 46 dB. Vertical bars show 95% confidence intervals of thesimulation results for STSL. . . . . . . . . . . . . . . . . . . . . . . . 372.4 ?err from (2.31) for the STSL algorithm as a function of number ofiteration steps. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 382.5 Empirical PDF of cest for c = 1500 m/sec. . . . . . . . . . . . . . . . 392.6 ?err from (2.31) for time-synchronization and sound speed uncertainties. 402.7 Time-varying location of nodes in the sea trial. . . . . . . . . . . . . . 412.8 T? pddiff for all communication links as a function of time slots. . . . . . . 422.9 ??err from (2.32) as a function of |c? c?| for W = 30 time slots. R-STSLmethod. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 432.10 Probability that ??err ? x for STSL and MR-STSL. . . . . . . . . . . . 443.1 Simulated scenario: (a) simulated velocity field at depth 20 m, (b)simulated speeds of the anchor nodes and TN. . . . . . . . . . . . . . 583.2 (a) Mean of ev(k) from (3.32) as a function of number of anchor nodes,(b) PD vs. PF for different ?2anc and The values. . . . . . . . . . . . . 593.3 Simulated scenario: (a) speed of the TN, (b) tracking error edk from(3.40) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61xiList of Figures3.4 Tracking performance (simulations): (a) empirical C-CDF for MSNR=50dB,(b) edk as a function of MSNR (thick line shows results for the CRLB). 623.5 Node locations in Cartesian coordinates: (a) Singapore sea trial, (b)Israel sea trial. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 633.6 edk from (3.40) vs. time: (a) Singapore sea trial, (b) Israel sea trial. . 654.1 Classification results: (a) Prd,non?ONLOS and Prd,ONLOS vs. ?UB, (b)Empirical detection probabilities of LOS and SNLOS. . . . . . . . . . 804.2 Empirical C-CDF of ?err from (4.27). . . . . . . . . . . . . . . . . . . 814.3 Estimation error of LOS and SNLOS distribution parameters vs. EMiteration number. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 824.4 Israel harbor experiment: (a) location (picture taken from Googlemaps on Sep. 09), (b) results for ONLOS link classification. . . . . . 834.5 Empirical C-CDF of ?err from (4.28), T?c = 10Ts. . . . . . . . . . . . . 854.6 Histogram of ?erri from (4.29). Bin width 0.3 m, E(di) = 324.1 m. . . 864.7 Localization error, ?e, as a function of the number of received packets,W . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 875.1 Flow chart for the DR-A algorithm. . . . . . . . . . . . . . . . . . . . 905.2 Graph representation for positive and negative constraints. . . . . . . 965.3 Example of a modeled wave in the x? y plane. t = 10 seconds. . . . 1045.4 C-CDF results of ?angle. . . . . . . . . . . . . . . . . . . . . . . . . . . 1055.5 Average results of ?err as a function of (a) M? (for T?c = Tc), (b) T?c (forM? = 10). Tc = 6 seconds. . . . . . . . . . . . . . . . . . . . . . . . . 1065.6 Results of ?err as a function of 1? . M? = 10, T?c = Tc. . . . . . . . . . . 1075.7 C-CDF results of ?err. 1? = 20 dB, M? = 10, T?c = Tc. . . . . . . . . . . 1075.8 Sea trial: nodes location in Cartesian coordinates. (a) x-axis. (b) y-axis.1085.9 Sea trial: measured and projected acceleration along the x axis relativeto g for M? = 5 and T?c = 40 seconds. Boat 1. . . . . . . . . . . . . . . 1095.10 Sea trial: (a) C-CDF results of ?angle for the two boats, (b) ?angle vs.the accumulated change of the boat?s heading angle within the timeslot, Tslot = 300 seconds. . . . . . . . . . . . . . . . . . . . . . . . . . 1095.11 Sea trial: C-CDF results of ?err for the two boats. Tslot = 200 sec . . . 1106.1 Illustration of channel reservation process (note that the first datapacket of node j is smaller than the rest). . . . . . . . . . . . . . . . 1176.2 Illustration of different types of primary conflicts. . . . . . . . . . . . 1196.3 Illustration of packet exchange mechanism. . . . . . . . . . . . . . . . 1256.4 Structure of RTS, CTS, NT, and data packets. . . . . . . . . . . . . . 1316.5 Scheduling of RTS and CTS and determining of CS parameters in theTSR protocol. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132xiiList of Figures6.6 Average of ?through from (6.2) and ?delayTmsgj,j?from (6.4) as a function ofRmsg =Tmsgj?,jTmsgj,j?for TSR with and without DSSS. . . . . . . . . . . . . . 1346.7 Empirical C-CDF of ?fair,through in (6.3) for the TSR, TDMA, Slotted-FAMA, BiC-MAC, and APCAP protocols. |N | = 6 . . . . . . . . . . 1356.8 Empirical CDF of ?delay in (6.4) for the TSR, Slotted-FAMA and BiC-MAC protocols. |N | = 15. . . . . . . . . . . . . . . . . . . . . . . . . 1366.9 Empirical C-CDF of ?through from (6.2) for the TSR, TDMA, Slotted-FAMA, APCAP, and BiC-MAC protocols. |N | = 6. (a) without DSSSsignaling, (b) with DSSS signaling. . . . . . . . . . . . . . . . . . . . 1376.10 ?change(|N | = 6, 10, 15) from (6.37). . . . . . . . . . . . . . . . . . . . 1387.1 Example: (a) sample topology for a UWAN. (b) constructing matrixIskel using Algorithms 4 and 5 for the sample network. (c) maskingmatrix used in (7.14). . . . . . . . . . . . . . . . . . . . . . . . . . . . 1497.2 Flow chart to obtain the R-BSP scheduling matrix M . . . . . . . . . 1527.3 Probability (7.23) for k = u?1 and empirical probability of ? ? u?1when generating a Bernoulli random graph with N nodes and edgeprobability p. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1547.4 Convergence parameters b1 and b2 vs. transmission range . . . . . . . 1567.5 C-CDF of ?through from (7.2) for matched and mismatched OTT andstatic model-based topology. . . . . . . . . . . . . . . . . . . . . . . . 1597.6 C-CDF of ?through from (7.2), for static model-based topologies anddifferent MPR scenarios. . . . . . . . . . . . . . . . . . . . . . . . . . 1607.7 C-CDF of ?through from (7.2) for model-based topologies. (a) static (nooutdated topology). (b) dynamic (slowly propagating topology) . . . 1617.8 CDF of ?delay from (7.3) for static model-based topologies. . . . . . . 1627.9 Structure of Sea Trial: (a) satellite picture of the sea trial location(picture taken from Google maps on September 29, 2009). (b) recordednetwork topologies. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1637.10 Results of ?through from (7.2) for topologies recorded in the sea trial andT = 10 time-slots. Each vertical line represents the time a topologychange starts affecting results. . . . . . . . . . . . . . . . . . . . . . . 1647.11 ?delay from (7.3), for static and dynamic sea-trial-based topologies. . . 1658.1 Gain gcoding from (8.5) as a function of ?t,i from (8.2). . . . . . . . . 1698.2 Gain genergy from (8.7) as a function of Th and ?t,i from (8.2). Non-erasure channel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1718.3 Illustration of the procedure of updating channel conditions for multi-ple packet transmission. . . . . . . . . . . . . . . . . . . . . . . . . . 1748.4 Illustration of adaptive coding implementation for single packet trans-mission. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175xiiiList of Figures8.5 Transmission loss vs. range. Output from the Bellhop simulator [2]. . 1788.6 PER as a function of channel erasure rate. . . . . . . . . . . . . . . . 1798.7 Channel utilization as a function of channel erasure probability (for?t,i ? 650 m). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1808.8 Channel utilization as a function of ?t,i (for erasure probability ofroughly 0.21). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1818.9 Error ?p(m) as a function of packet number m. (Note that values inthe y-axis are reversed). . . . . . . . . . . . . . . . . . . . . . . . . . 1828.10 Histogram of network latency in terms of number of packets neededtill decoding. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1838.11 CDF of the number of symbols transmitted till successful decoding,Nfinal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1838.12 CDF of network goodput. . . . . . . . . . . . . . . . . . . . . . . . . 1848.13 Satellite picture of the experiment location (picture taken from GoogleEarth on July 23, 2012.) The three locations of the nodes are marked. 1858.14 Packet success rate for the three links from the sea trial. . . . . . . . 186xivList of Abbreviations and SymbolsAcronymsARQ automatic repeat requestAUV autonomous underwater vehicleBSP broadcast scheduling problemCA collision avoidenceC-CDF complementary cumulative density functionCI confidence indexCL connectivity listCRLB Cramer-Rao lower boundCS communication sessionCSI channel state informationCTS clear-to-sendCUMP channel-utilization-maximization problemDR-A DR navigation using a single accelerometerDR dead-reckoningDSSS direct sequence spread spectrumDVL Doppler-velocity-loggerEKF extended Kalman filterEM expectation-maximizationHCRLB hybrid Cramer-Rao lower boundHST-TDMA hybrid spatial-reuse time-division multiple accessINR interference-to-noise ratioINS inertial systemIR-HARQ incremental redundancy hybrid automatic repeat requestLDPC low-density parity checkLOS line-of-sightMAC medium access controlMIS maximal independent setMPR multiple packet receptionNLOS non-line-of-sightNT noticationONLOS object-related NLOSOTT orthogonal topology-transparentxvList of Abbreviations and SymbolsPCA principal component analysisPD propagation delayPDF probability density functionPER packet error rateR-BSP robust BSPRS Reed-SolomonRTS request-to-sendSD-UT spatially dependent underwater trackingSINR signal-to-interference-plus-noise ratioSNLOS sea-related NLOSSNR signal-to-noise ratioSSP sound speed profileSSM space state modelSSV state space vectorSTSL sequential time-synchronization and localizationT-BSP topology-dependent BSPTDMA time-division-multiple-accessTN tracked nodeToF time-of-flightTSR joint time and spatial reuseUAC underwater acoustic channelUKF unscented Kalman filterUL unlocalizedUT underwater trackingUWAC underwater acoustic communicationUWAN underwater acoustic communication networksUWL underwater localizationNotations and OperatorsBold upper case and lower case letters denote matrices and vectors, respectively.Accents x? and x? represents estimation and approximation of x, respectively. Theremaining notation and operators used in this thesis are listed below.L Number of anchor nodes directly connected to the UL nodeji 2-D UTM coordinates of the UL nodepi 2-D UTM coordinates of the lth anchor nodec Sound speed in water [m/sec]W Duration of the localization windowN Number of packets transmitted during the localization windowSl Clock skew of the UL node relative to the lth anchor nodexviList of Abbreviations and SymbolsOl Clock offset of the UL node relative to the lth anchor node [sec]T pdi Propagation delay of the ith packet [sec]Ti Transmission local time of the ith packet [sec]Ri Reception local time of the ith packet [sec]d?i,i? Self estimation of distance between locations ji and ji???i,i? Self estimation of angle between locations ji and ji? [rad]? Threshold for location quantization [m]?2 variance of ToA measurement noise [sec2]?li Ratio between Tpdi and the actual propagation delayak state space vectoryk measurement vectorhk measurement model vectorrk 3-D coordinates of the TNranck 3-D coordinates of the anchor nodedk distance vector rk ? ranckvk velocity of the TNv?drift drift velocity estimate of the TNvanc,driftk estimated drift speed and heading direction of the anchor nodec?k average sound speed between the anchor and TN? INS time interval between consecutive INS measurements? range time interval between consecutive range measurementsmINSk INS measurement vectormToFk ranging measurementmDopplerk Doppler shift measurementxi PD measurementsxLOS, dLOS delay and distance in the LOS link, respectivelyX vector of PD measurements xi of the same communication linkdRSS,LB lower bound of RSS-based range measurement? propagation loss factor? absorption loss factorM assumed number of classes in PD modelkm weight of the mth distribution in the mixture distribution model?m vector of parameters of the mth distribution? vector of parameters of the distribution of XTLIR upper bound on the length of the channel impulse responsec propagation speed in the channeldONLOS distance of the ONLOS linkRL reflection loss in ONLOS link?l group of xi measurements with the same distribution ?m?l classifier for group ?l%i classifier for xixviiList of Abbreviations and Symbolsdi,j distance traveled by the vessel at time period [ti, tj]N number of acceleration measurementsTc coherence time of acceleration along the horizontal plane?n the vessel orientation with respect to the reference system?n the vessel heading angle with respect to ?n?n the vessel pitch angle in time tnan the vessel 3-D acceleration in the horizontal plane coordinate systema?n acceleration measurementA? vector of acceleration measurementsA?h Matrix A? projected to match the vessel headingA?h,p Matrix A? projected into the horizontal plane along the vessel headingM? pre-defined number of pitch-states?l group of same pitch-stateL number of groups ?l?l classifier for group ?l??p the distribution parameters of A? at the pth EM iteration?m,x mean of the mth class for the x axis?m,x standard deviation of the mth class for the x axiskm prior probability of the mth classxviiiAcknowledgmentsFirst and foremost, I would like to express my utmost gratitude to my supervisor, Pro-fessor Lutz Lampe, for his assistance, support, and guidance throughout the course ofthis work. Your fine sense of constructive criticism was invaluable to me throughoutthe work.I would also like to thank Dr. Hwee-Pink Tan and Dr. Mandar Chitre for theirvaluable comments and assistance.I would like to thank my parents for believing in me and encouraging me through-out my life. I am indebted to my father and my mother for their endless love, care,and support.Last but not least, I would like to thank my lovely wife, Hadas, for raising ourfamily, my greatest accomplishment, joy, and pride. Without your support none ofthis could have happened.xixTo Hadas, Ronya, Gili, and Itamar?speak to the earth, and it shall teach thee?xxChapter 1IntroductionThe oceans with their diverse biology systems, vast energy resources, and climate andhistory records of our planet, have always attracted researchers and industries. In thelast decade, ocean exploration has considerably increased through the use of ocean-observation systems, autonomous underwater vehicles (AUVs), and fixed or mobilesensor networks. These submerged devices need to report the collective data backto base stations or to share information in the setting of a wireless communicationnetwork. Wireless communication underwater is usually established using acoustictransducers since radio frequency communication is only possible for very short dis-tances underwater. Underwater acoustic communication (UWAC) can fulfil the needsof a multitude of underwater applications, including: oceanographic data collection,warning systems for natural disasters (e.g., seismic and tsunami monitoring), ecolog-ical applications (e.g., pollution, water quality and biological monitoring), militaryunderwater surveillance, assisted navigation, industrial applications (offshore explo-ration), to name just a few [3]. For example, in offshore engineering applications,underwater sensors can measure and report parameters such as foundation strengthand mooring tensions to monitor the structural health of deepwater mooring sys-tems. In addition, underwater acoustic communication networks (UWANs) comprisecommunication between surface stations and AUVs. Two common communicationsarchitectures for UWANs are shown in Figure 1.1.Most of the UWAC research in the past has concentrated on the development ofmodels for the underwater acoustic channel (UAC) and the design of secure pointto point links. Only recently, networking aspects of UWAC consisting of a fixed butad-hoc infrastructure (alike base stations in a cellular network) and mobile AUVs(alike mobile phones) have started to attract significant interest in the research com-munity to enable fast, reliable, and cost effective UWAN [4]. One of the majorchallenges of UWANs is resource assignment which has become the bottleneck forUWAN applications. Scheduling for UWANs is governed by several factors such aslow sound speed (approximately 1500 m/sec [5]), half-duplex communication (due todesign constraints of the acoustic transducers [6]), large power attenuation, long delayspread, time varying propagation channel, and a very limited signal bandwidth dueto absorption loss (which increases with frequency) and noise level (which decreaseswith frequency) [7], as well as transducer constraints. These factors lead to generallypoor availability and reliability performance of UWAN systems and pose engineeringchallenges that are very different from those experienced in radio frequency wirelessnetworks [8, 9]. In particular, even though UWAC systems have been implemented1Chapter 1. Introduction(a) (b)Figure 1.1: Communications architecture for UWSNs [1].for many years, many design and implementation problems remain to be solved to-wards developing advanced communication network systems enabling applicationswith high quality of service requirements.The long propagation delay in the channel, as well as the often sparse nature ofUWANs, make location-dependent medium access control (MAC) protocols highlyattractive to improve latency and throughput performance of UWANs via spatial-reuse techniques. In this context, underwater localization (UWL) of unlocalized (UL)UWAN nodes (infrastructure elements and mobile devices) is a key element towardsefficient communication networks enabling, e.g., an ocean-observation system. (Wenote that the network has to provide its own localization, since global positioningsystem (GPS) systems do not work underwater.) Moreover, sensed data can only beinterpreted meaningfully when referenced to the location of the sensor. The followingare desirable properties of UWL:? AccuracyThe location of the sensor for which sensed data is derived should be accu-rate and unambiguous for meaningful interpretation of data. Localization al-gorithms usually minimize the distance between the estimated and the truelocation.? SpeedSince nodes constantly move, the localization and its tracking procedure shouldbe fast so that it reports the actual location when data is sensed.? Wide Coverage2Chapter 1. IntroductionThe localization scheme should ensure that most of the nodes in the networkcan be localized.? Low Communication CostsAccurate localization requires ranging estimation to anchor nodes whose loca-tion is known, usually via UWAC. Since the nodes are battery-powered and maybe deployed for long durations, communication overhead should be minimized.In addition to the above quantifiable properties, practical considerations such as easeand cost of deploying reference nodes and other required infrastructure should betaken into account too.The remainder of this chapter is organized as follows. In Section 1.1, we presenta survey of UWL and tracking schemes, and identify important challenges that needto be addressed. In Section 1.2, we describe the state-of-the-art in UWAN MACprotocols with a focus on spatial reuse scheduling that utilize location information.In Section 1.3, we list the open problems in UWL and UWAN MAC design addressedin this thesis. Finally, in Section 1.4 we present the structure of this thesis.1.1 Underwater Localization and TrackingAlthough localization has been widely studied for terrestrial wireless sensor networks,existing techniques cannot be directly applied to UWANs because of the followingunique characteristics:Deployment of Anchor Nodes Assuming depth sensors are used, to localize un-derwater nodes deployed in the 3D sea environment, reference locations of at leastthree anchor nodes are required. However, since due to restrictions on energy con-sumption for long deployment period and sparse network topology, localization cov-erage is limited, and a node may not always be in the communication range of atleast three anchor nodes.Node Mobility Underwater nodes will inevitably drift due to the water current,winds, shipping activity etc. [10]. While anchor nodes attached to surface buoyscan be precisely located through GPS updates, it is difficult to maintain submergedunderwater nodes at precise locations. This may affect localization accuracy, as somedistance measurements may have become obsolete by the time the node position isestimated. Furthermore, due to the unpredictable nature of the ocean current it ishard to track the location of drifting nodes using a predefined motion model.Inter-Node Time Synchronization Since GPS signals are severely attenuatedunderwater, it cannot be used to time-synchronize nodes deployed underwater to3Chapter 1. Introductioncompensate for clock drifts due to both offset and skew. Consequently, the accuracyof time-of-flight (ToF)-based range measurement may be affected. Furthermore, thespeed of sound underwater is five orders of magnitude lower than RF propagationover the air. Hence, both clock skew and offset should be compensated for.Signal Reflection In near-shore or harbor environments, where obstacles may existbetween nodes, non-line-of-sight (NLOS) signals reflected from object (e.g., vessels,harbor wall), or multipath (indirect) signals from the sea surface or bottom can bemistaken for line-of-sight (LOS) signals, and may significantly impact the accuracyof range measurement.Sound speed variation Unlike the speed of light which is constant, the speedof sound underwater varies with water temperature, pressure and salinity, givingrise to refraction. Without measuring the sound speed, the accuracy of distancemeasurements based on time-of-arrival approaches may be degraded.Asymmetric Power Consumption Unlike RF modems, acoustic modems typi-cally consume much more power (order of tens of Watts) in transmit mode comparedto receive mode (order of milliWatts). This asymmetry in transmission mode makesit preferable for ordinary nodes to be localized through passive / silent listening.Low bit rate Compared to RF communications, the bit rates achievable withacoustic communications is significantly lower. As a result, localization packets hold-ing anchors? location information are long and have high impact on network through-put.Figure 1.2 maps the above challenges to each desirable localization performancemetric.In the following, we start with a survey of the state-of-the-art in UWL. Next,we present current works on tracking the time-varying location of underwater nodes,and discuss approaches to mitigate localization measurement errors.1.1.1 Underwater LocalizationWe review both range-free and range-based UWL techniques. In range-free schemes,UL nodes may infer their proximity to anchor nodes (e.g., in terms of number of hops)so as to achieve coarse localization, e.g., in an area instead of a specific location.Range-based approaches rely on time and/or bearing information to evaluate thedistance to anchor nodes, usually using acoustic communications. The UL node thenutilizes multilateration/angulation to estimate its own location.4Chapter 1. IntroductionFigure 1.2: Mapping between the challenges and desirable properties of underwaterlocalization.Range-free LocalizationRange-free UWL schemes are designed for cases where range measurements sufferfrom large errors due to node mobility or strong attenuation. In [11], a UWL schemebased on an assumed attenuation model is proposed. When an anchor node i trans-mits at power Pi, the UL node can receive the transmission as long as it falls withinthe anchor node?s transmission range, which depends on the transmission power.Hence, by deploying several reference nodes that transmit beacons at multiple powerlevels, the plane is divided into many small sub-regions defined by intersecting circles.By receiving reports from each UL node of the minimum transmit power at which itreceived the respective anchor node beacon, a central sink can estimate the locationof the UL nodes. However, it is a centralized scheme where mobility is not consid-ered. A range-free UWL scheme for moving anchors is presented in [12]. Here, anAUV traverses a preprogrammed route and performs directional (vertical) beaconingperiodically. Assuming that the AUV moves with constant and known speed andknows its position underwater accurately, relying on the AUV periodic broadcasts,the UL node can estimate its position within a circle formed by the intersection ofthe transmitted beams with the horizontal plane.A different variant of range-free UWL schemes based on finger-printing have beenproposed in [13]. Such schemes involve an offline (or training) stage prior to theonline (or prediction) stage. The setup comprises an acoustic signal source capable oftransmitting at M different frequencies. During the offline stage, a receiver is placedat a reference location (with known position), and collects N samples of receivedpower at each frequency to constitute anM?N acoustic-signal map. This is repeatedat each reference locations. In the online stage, the receiver measures received powersamples from M different frequencies and compare these to the ones samples in the5Chapter 1. Introductionoffline stage using a likelihood function and a ?probabilistic-weighted? summation ofdifferent reference locations.Range-based LocalizationRange-based localization typically comprises the following steps:1. Range measurement: Each UL node estimates its distance from each anchornode within its communication range using either received signal strength indi-cator (RSSI), time difference of arrival (TDoA) or time of arrival (ToA). Sincethe path loss in UAC is usually time varying and multipath effect can result insignificant energy fading, the RSSI method is not the primary choice for under-water localization. Hence, most proposed range-based localization schemes useeither TDoA or ToA for ranging. The TDoA method involves computing thetime difference of arrival between beacons from different anchor nodes trans-mitted using acoustic signalling, and the ToA method performs ranging basedon the relationship among transmission time, speed and distance.2. Location estimation: Each UL node then estimates its position, typically,according to the intersection of various circles centered at each reference nodewith radii corresponding to the range measurements. In general, to localize anode in d-dimensional space, the number of independent range measurementsrequired should be at least d + 1.3. Tracking: The location estimate is refined e.g., using measurements from on-board sensors, measurement error models, mobility models, etc.Range-based UWL method can be classified into 1) single-stage schemes whichrely solely on message exchanges with the anchor nodes, and 2) multi-stage schemesin which newly localized nodes can serve as anchor nodes. The key innovation of thefirst type of UWL schemes lie in how they address localization inaccuracy due to time-synchronization and measurement errors and availability of anchor nodes. While itis usually assumed that clock offset is the main cause of time-synchronization errors,also clock skew cannot be neglected for UWL due to the long propagation delay in theUAC [14]. Furthermore, since anchor nodes are usually submerged we cannot assumethese nodes to be time synchronized. Regarding this problem, [14] suggested toestimate both skew and offset based on packet exchange with an already synchronizednode. Alternatively, the problem of time-synchronization can be avoided by usingTDoA techniques. In [15] ?silent positioning? is provided, i.e., UL nodes do nottransmit any beacon signal and just listen to the broadcasts of reference nodes toself-position, reducing the communication costs. The scheme relies on TDoA overmultiple beacon intervals, and thus does not require time synchronization amongstnodes. However, this kind of ?reactive beaconing? makes it susceptible to failure dueto transmission losses that are prevalent in harsh UACs. An improvement has been6Chapter 1. Introductionsuggested in [16], where a dynamic mechanism for leader reference node identificationand a time-out mechanism to trigger beaconing in the event of transmission loss ispresented. However, the scheme does not manage motion of nodes.AUV (t0)AUV (t2)AUV (t1)2-way message exchangeAUV routeFigure 1.3: Illustration of AUV-aided localization.Allowing node mobility is important for UWL, where the deployment of fixed ref-erence nodes such as surface buoys is time consuming, limits the localization coverageand may be infeasible or undesirable. Several algorithms have been suggested to com-pensate for node movements, either by regarding these movements as ToA noise [17]or by applying mobility prediction [10]. A few UWL schemes rely on a moving anchornode, where, as illustrated in Figure 1.3, a single node can localize the network. Inthe AUV-aided scheme proposed in [18], the AUV obtains position updates by risingto the surface to use GPS, and then dives to a predefined depth and periodicallyperforms a two-way message exchange with UL nodes. A 2D localization is achievedonce successful two-way message exchanges take place in at least three non-collinearAUV locations. However, UL nodes should remain static. Instead of AUVs, the?Dive?N?Rise? localization scheme [19] uses a weight/bladder mechanism to controlthe diving/rising of each mobile beacon. These beacons update their positions atthe surface, and broadcast them when they dive to a certain depth. However, theseapproaches consider node movements as an undesired phenomenon and do not utilizethe possibility of additional ToA measurements when coupled with self-localization.7Chapter 1. IntroductionWhile, ideally, three range measurements should be enough to localize a nodein a 2D plane through triangularization, ranging errors require the use of multi-ple range measurements through multilateration. Instead of the commonly-adoptedcircle-based least-squares multilateration, [20] suggested a hyperbola-based approach.The premise is that when range measurement errors due to imperfect time synchro-nization, or varying speed in acoustic transmission exist, two hyperbolas always in-tersect with each other with one cross point, or partial solution, while two circleswill likely intersect with either two or zero cross point(s). Several works suggestedmethods to compensate for location ambiguities such as flips and rotations that arisedue to NLOS-related range estimation errors. In [21], a three-phase algorithm issuggested, where first an ambiguity-free sub-tree of nodes is determined. Then, mul-tilateration is performed where the node is first assumed to be located in the center ofa rectangular area. Finally, a refinement phase is performed using a Kalman filter tomitigate remaining noise. A robust algorithm for mitigating localization ambiguitiesis suggested in [22] by rejecting measurements leading to ambiguities, e.g., when thereare insufficient anchor nodes or when the location of anchor nodes is almost collinear.In [23], an NLOS factor (i.e., the difference between the arrival times of the NLOSand LOS-based signals) is estimated using a maximum likelihood estimator basedon an attenuation model, and NLOS-based measurements are incorporated after afactor correction instead of being rejected. The problem of localization when all mea-surements are obtained from NLOS links is considered in [24], where the relationshipbetween anchor node distances and the NLOS factor is used to improve localization.Another significant challenge of UWL is the variability of the sound speed in wateras it depends on water temperature, depth and salinity [5]. Considering this problem,[25] suggested to first estimate the sound speed using packet exchange between float-ing buoys on the seabed and the sea surface. Alternatively, several works suggestedsound speed estimation based on measuring the channel characteristics and a soundspeed model, e.g., [26], [27]. Differently, [28] suggested to jointly estimate the nodelocation and the propagation speed in the channel by considering the propagationspeed as an additional variable. However, they made the assumptions that all nodesin the network are static, at least four anchor nodes are available and that all nodesare time synchronized, which do not hold true in most UWL applications. Whilechanges in the sound speed profile (SSP) with depth have been considered for UWL(e.g., [29]), for shallow water and short to medium range applications, a good approx-imation for the depth-related sound speed in water is a single parameter, c?, whichchanges as a function of the mean depth of the transmitter and receiver and can beconsidered as an average sound speed. In [28], c? is treated as a system parameterand is separately estimated during the initial localization process. However, due tochanges in the depth of the tracked node (TN), the evaluated parameter c? may shiftand its value should be tracked over time.The key innovations of proposed schemes within the category of multi-stage UWLlie in the trade-off between minimizing error propagation and delay while maximizing8Chapter 1. Introductioncoverage and energy efficiency. In [30], a two-phase algorithm is proposed where arelative coordinate system is built using the first three discovered UL nodes, followedby a selection of new nodes to localize based on their proximity from the new lo-calized nodes. However, the first-stage discovery algorithm requires high volume ofmessage exchange, node mobility is not considered, and only the relative coordinatesfrom the primary seed node is acquired. To combat nodes? mobility, the authors in[31] proposed a joint localization and synchronization scheme. The 3D network ispartitioned into cells, and localization is performed at the cell level. The authorsdetermined the required sensor node density, as well as cell partitioning in order tolocalize all nodes.Since multi-stage localization might lead to error-propagation in the network, self-evaluation of the localization accuracy is required to determine if the localized nodeis elidible to serve as a new anchor node. Several works suggested methods for suchself-evaluation, relying either on a node?s expected location in the network [30] orby assigning each node with a confidence index. In [32], a hard-decision confidenceindex is used based on detecting localization based on outdated range information.Alternatively, in [33], soft decision confidence value is obtained by normalizing theexpected position error with the sum of the Euclidean distance to anchor nodes. How-ever, in doing so these methods reuse the information used already for localization,which might cause biased self-evaluation of the localization accuracy.1.1.2 Tracking the Location of NodesAfter obtaining initial location information via UWL, in case of a moving node, atracking procedure begins where the time-varying location of the node is recursivelyestimated. This process is referred to as underwater tracking (UT). Since GPS isnot available underwater, UT has similarities to indoor navigation. However, due tothe difficulty of modeling motion at sea, UT includes additional challenges. First,the unpredictable water current with fast changes in speed and direction, and theexistence of water turbulences, cause irregularities in the motion of the TN [34].Second, the speed of acoustic signals changes with depth. Last, compared to RADARapplications, where directivity is applied and the emitter is fixed or slowly movingand the propagation speed in known, uncertainties in the sound speed as well as thecontinuous motion of nodes also make it more challenging to incorporate Dopplershift measurements in the UT scheme.The scenario adopted by most UT algorithms consists of several anchor nodesand a moving TN, usually an AUV, with an initial location estimate [34]. A varietyof on-board sensors are used. The most common are inertial sensors, providing ac-celeration measurements. Since acceleration is affected by the vessel pitch and roll,an inertial system (INS) also includes a gyrocompass, vibrating angular rate sensors,or Pendulum tilt sensors [34, 35]. During the UT process, occasionally, external in-formation through packet exchange with anchor node is also available. This includes9Chapter 1. Introductionranging information as well as Doppler shift measurements. The latter are commonlyused in RADAR and Doppler velocity log (DVL) systems, where angle-to-target ismeasured and propagation speed is assumed known, and it provides information ofthe relative velocity of a tracked object. In [36], Doppler shift is also used for UWL,assuming a fixed scenario with known sound speed. However, since in underwateracoustic communication omnidirectional transducers are used, the TN and the an-chor nodes constantly move, and the sound speed is unknown and may change overtime, incorporating Doppler shift measurements into UT is challenging.External and internal measurements for UT are often incorporated using an statespace model (SSM) to update a state space vector (SSV), and a key point is howto reliably determine the SSM. In [37], this problem is handled through a bank ofKalman filters (KFs), each of which uses a different motion model, and the outputof each branch is combined according to the accumulated filter errors. To combatmodel inaccuracies, [38] considered both regular and random motion, and the twocomponents are combined as inner-states of two Markov-like random states. Whilesuch an approach offers flexibility in determining the SSM, it is difficult to fit atracking filter to such a model. Indeed, a common SSM is of fixed speed with a randomGaussian acceleration [34]. To set the model parameters under motion irregularities,the velocity of the TN can directly be measured using a DVL [34]. Unfortunately,DVLs consume non-negligible energy and work well only close to the ocean bottomor surface. Alternatively, the velocity can be estimated by calculating the motor-induced thrust force [39]. However, such velocity measurement is only relative to thewater medium, and thus may not be accurate in the presence of high ocean current.Updating the SSM in time requires the implementation of a tracking filter. In[39], an extended Kalman filter (EKF)-based method is presented to fuse measure-ments from on-board sensors with occasional range measurements to anchor nodes.An interesting combination of the KF and the EKF is presented in [40], where theformer accounts for motion in the surge direction and the latter for angular motion.The scheme includes outlier mitigation using a probabilistic data association filterat the input to the KF. In [41], fusion of sensor information and ranging is formu-lated as a maximum likelihood optimization problem and the solution is found by anon-linear weighted least-squares optimization. Alternatively, [42] suggested incor-porating maximum likelihood data association in an EKF scheme. By using rangingand bearing measurements to several anchor nodes, compass sensor data, and DVL-related speed measurements, good tracking performance is obtained. In [43], it wassuggested to include previous locations of the TN in the SSV to combat delays in therange information due to the slow sound speed in the channel, and instead of an EKF,an extended information filter is used to reduce complexity in case the informationmatrix is sparse. A similar approach is presented in [35], where a post-processingcentralized EKF is used to incorporate both anchor and TN sensor information. Theused SSV is highly detailed and includes the TN pose, its depth, and its linear andangular velocities.10Chapter 1. Introduction1.1.3 Considerations of Measurement ErrorsLocalization and tracking accuracy is highly affected by measurement errors, wherethe dominant factors are ranging and acceleration measurement errors [34]. Theformer is affected by the long delay spread in the channel and strong multipath, andthe latter is caused by ambiguities of the node orientation, mostly due to ocean waves.In this section, we review current approaches to mitigate such errors.Propagation Delay ErrorsRanging information, required for both UWL and tracking, is obtained from esti-mating the ToF of a received signal, using either ToA or TDoA techniques. ToFmeasurements for range estimation can be obtained (i) from multiple impulse-typesignals transmitted in a short period of time, or (ii) from the symbols of a receiveddata packet. The former is a standard in many ultra short baseline systems for rang-ing underwater (e.g., [44]) and involves inspecting the output of an energy detector[45]. The latter involves inspecting the estimated channel impulse response by per-forming a matched filter operation at the receiver [46], or by performing a phase-onlycorrelation and using the kurtosis metric to mitigate channel enhanced noise [47].The ToF is then estimated by setting a detection threshold to identify the arrival ofthe first path. In [48], a fixed threshold is set based on the channel noise level anda target false alarm probability. In [49], an adaptive threshold is used based on theenergy level of the strongest path. A good overview of practical ToF estimators isgiven in [45].When estimating the ToF, one has to lock on to a certain arrival path, believed tobe the LOS path between sender and receiver. Existing UWL schemes, e.g., [15, 33],implicitly assume that PD measurements correspond to the LOS link between thetransmitter and receiver. However, signals can arrive from NLOS communicationlinks in several ways, as illustrated in Figure 1.4. For the node pairs (u, a2) and(u, a3), sea surface and bottom reflections links (referred to as sea-related NLOS(SNLOS )) exist, respectively, in addition to an LOS link. For (u, a1), the signalarrives from the reflection off a rock (referred to as object-related NLOS (ONLOS )).Lastly, between nodes u and a2, there is also an ONLOS link due to a ship. Whileit is expected that power attenuation in the LOS link is smaller than in NLOS links,it is common that the LOS signal is not the strongest. This is because, as shown inmultipath models [50, 2] as well as measurements [46], the UAC consists of groups ofNLOS links with small path delay, but significant phase differences, often resultingin negative superposition with the LOS link (if delay differences between the LOSand NLOS links are smaller than the system resolution for path separation) as wellas positive superposition between NLOS links. If PD measurements of NLOS linksare mistakenly treated as corresponding to delay in the LOS link, e.g., in node pairs(u, a2) and (u, a3), ranging accuracy will significantly be degraded. Clearly, rangingaccuracy affects localization performance. For example, using basic trilateration, the11Chapter 1. Introductiona1a2Unlocalized node, uAnchor nodea3Sea surfaceSea bottomONLOSSNLOSLOSd21dLOSd22ShiprocksuFigure 1.4: Illustration of various types of communication links: LOS, SNLOS andONLOS links.localization error grows quadratically with ranging offset, and a zero-mean Gaussiandistributed offset with a standard deviation of only 2 msec (which is quite commonin UWL [50, 2, 46]) would cause an average error of 6 m error.Acceleration ErrorsTracking the location of a node usually involves the use of acceleration measurementsproduced by an INS. By integrating INS measurements the distance traveled bythe TN can be estimated. This process is referred to as dead-reckoning (DR). Themain challenge in DR navigation is the possible drift and measurement noise ofinertial sensors, which may lead to errors in the order of 10% of the traveled distance,depending on the technology employed [51, 34]. For pedestrian applications, usingthe fact that velocity and orientation can be set to zero when the foot is on ground,it is customary to mount inertial sensors on foot or hip and estimate distance andorientation separately for each step, e.g., [52]. However, while at sea we can identifytime instances where the vessel pitch angle is zero, i.e., at the top or bottom ofthe ocean wave (see Figure 1.5), velocity of the vessel cannot be assumed zero atthese points. Instead, reference measurements are used, e.g., DVL, and measurementdrifts are controlled through fusion of large number of inertial sensors [53]. Forships, DR navigation involves dynamic positioning, heading autopilots, and thruster-assisted position mooring [54]. However, for small AUV-type TNs with strict energy-constraints, these options are not available.Another challenge in DR navigation is to determine the orientation of the inertial12Chapter 1. Introductionxz?zx?y?y?nz? x?z?x? z?x?x?z??nFigure 1.5: Illustration of the vessel?s wave-induced motion.sensor with respect to a reference coordinate system, usually using a gyrocompass [55].Only then the distance traveled by the tracked object can be estimated. However,when the vessel is located close to or on the sea surface such that its motion isaffected by the ocean waves, the vessel pitch angle is fast time-varying and orientationmeasurement may be too noisy to use directly [54]. For this reason, traditionally, DRnavigation at sea involves integrating a large number of inertial sensors and applyingBayesian filtering methods, e.g., EKF or particle filters (cf. [51, 56]), to mitigateoscillatory wave-dependent components and reduce measurement drifts [53].1.2 Spatial Reuse Scheduling in UWANsApart from navigation purposes, location information can greatly improve through-put in the setting of a wireless sensor network like a UWAN. For example, location-aware MAC protocols, for which scheduling of nodes? transmissions are set accordingto their geographical or relative location in the network (e.g., [57]), can greatly im-prove throughput and/or latency by utilizing network resources more efficiently whilemaintaining scheduling limitations. Another example are adaptive coding techniques,where due to the usually low reliability of communication in UWANs, the transmis-sion scheme changes adaptively as a function of the transmitter-receiver distance.In this section, we review current scheduling algorithms and adaptive transmissionstechniques for UWANs.1.2.1 Scheduling Algorithms for UWANsScheduling transmissions in UWANs is required for applications with wide range ofrequirements, e.g., latency, size of information packets, traffic rates, and reliability.Since UWANs are relatively small (usually in the order of tens of nodes), centralizedscheduling approaches are considered and both contention-based and contention-freescheduling algorithms are in use [58]. The need to develop reliable communicationlinks and the high cost of retransmissions due to the long transmission delay and large13Chapter 1. Introductionenergy consumption for transmission [59] make the handshake-based multiple accesswith collision avoidance (MACA) protocol the method of choice for contention-basedscheduling of long unicast transmissions in UWANs [60, 8]. MACA was inspired bythe carrier sense multiple access/collision avoidance (CSMA/CA) technique, stan-dardized in IEEE 802.11a, where a communication session (CS) is established byexchanging request-to-send (RTS) and clear-to-send (CTS) packets. Differently, dueto the narrow bandwidth available for UWAC, contention-free scheduling of UWANsrelies on time-division multiple-access (TDMA) algorithms, where each node is as-signed a unique time slot and each time slot includes guard interval to compensatefor the (long) propagation delay and (short) time-synchronization offset [61]. Suchscheduling is used for broadcast communication or short unicast transmissions inheavy load networks. By scheduling transmissions to avoid packet collisions, bothhandshake- or TDMA-based scheduling benefit from location information. In fact,in [62] it was shown that by carefully scheduling transmissions, unlike for terrestrialradio-frequency networks where network throughput decreases with the number ofnodes, the optimal network throughput achieved for UWANs is 12 .Handshake-based SchedulingDue to the long propagation delay in the UAC the exposure time to packet collisionsis long [60], and a modification to the basic handshake protocol is required. Con-sidering this problem, [60] suggested a slotted handshake protocol in which globallyestablished time slots, the size of which is comparable to the propagation delay, areused, and transmissions are restricted to the beginning of these time slots. The au-thors of [63] suggested improving channel utilization by employing separate time slotsfor control and data packets, and in [64] MAC throughput was further improved byallowing the receiver to warn the transmitter of expected interferences.In handshake-based scheduling, channel resources are allocated to the transmitterwhose RTS packet was the first to arrive at the receiver, whereas other nodes need totry again to gain channel access after waiting for a certain backoff period. This mightlead to an increased delay in packet transmission as the probability of successfullyreserving the channel is inversely proportional to the transmission distance [65]. In[60], this problem was resolved by using a fixed backoff-window size instead of anadaptive one as standardized in IEEE 802.11a. However, determining the size of thefixed backoff-window is difficult as it has opposite effects on throughput and transmis-sion delay [66]. For this reason, [66] suggested an algorithm to bring randomness tochannel reservation by allocating transmitters with time-varying random ranks andgiving priority to the transmitter with the currently highest rank. Unfortunately,they assumed control packets arriving almost simultaneously to the receiver, whichdoes not hold true in UWANs.In addition to the problem of increased delay, traditional handshake-based MACprotocols require nodes in the proximity of a CS to remain silent, which limits channel14Chapter 1. Introductionutilization. This effect is even more noticeable in UWANs, where longer silenceperiods are imposed by the long propagation delay [60]. This also leads to the exposedterminal problem1, which further decreases channel utilization. In [67], the longpropagation delay in the channel is utilized to allow simultaneous transmissions inexposed terminal links. Upon detecting a packet from node p directed to node p?, anode j schedules its transmissions such that its packets arrive at p before the responseof node p?. One way to increase channel utilization is to use timing-advance techniques,often called time reuse [67], such that more nodes can transmit. In UWAC, time reuseis related to the utilization of the long propagation delay in the UAC. In terrestrialradio-frequency networks, performance is inversely proportional to the network size.However, by utilizing the long propagation delay, in UWAC performance is potentiallyfixed for different network sizes [62]. In [68], time reuse is applied by allowing nodesto initiate another CS while waiting for a CTS response. [69] suggested a distance-aware protocol where channel utilization is improved by allowing both nodes involvedin a handshake CS to transmit simultaneously. Alternatively, in [70], time reuse wasapplied by allowing nodes to opportunistically transmit data packets to a node j,such that the packets arrive at j upon completion of its own transmissions. However,both in [68] and [69], nodes located within the interference range of a transmitteror receiver should remain silent, and in [70] simultaneous transmission in different(connected) CSs is not allowed, and thus the above mentioned limitation of traditionalhandshake-based protocols remains.Since network nodes at different locations experience different channel-access lim-itations, applying spatial reuse on top of time reuse can further improve channelutilization. Spatial reuse refers to simultaneous CSs in different parts of the network,and it is specifically applicable to UWANs, since low-power half-duplex transceivers,range and frequency dependent absorption loss [5], and acoustic NLOS scenarios leadto sparse network topologies. The spatial-reuse handshake protocol suggested in [71]identifies exposed terminal links when a node overhears RTS packets which are notfollowed by CTS responses. This node can then transmit simultaneously. [72] sug-gested using control gaps in each exposed terminal link, allowing nodes to scheduletheir transmissions via RTS/CTS packet exchange during those gaps. Alternatively,in [73] exposed terminal links are identified by building a conflict map using a trialand error procedure without the need to exchange RTS/CTS control packets.TDMA-based SchedulingWhen broadcast communication is required, or when transmitted messages are shortcompared to the propagation delay in the channel, the overhead of RTS/CTS pack-ets becomes significant due to the collision probability being comparable to that for1The exposed terminal problem occurs when a node, upon detecting other transmissions, decidesnot to transmit, even though these transmissions may not interfere with its own transmission, andvice versa [60].15Chapter 1. Introductionpayload packets. Furthermore, since nodes detecting an RTS or CTS packet shouldbe kept silent, channel utilization decreases. As shown in [74], the number of silencednodes grows quadratically with the communication range and the CSMA/CA algo-rithm becomes more and more conservative. Therefore, CSMA/CA is not suitable tomeet high network traffic demands for broadcast communication or for short unicastinformation packets [74]. Moreover, channel reservation cost using MACA is high inUWANs, where the long propagation delay necessitates longer silence periods [60].In fact, when considering high traffic networks, the conventional TDMA algorithmseems to outperform most of the existing random-access algorithms [74].In TDMA, besides packet duration, time slots include the expected propagationdelay at the maximal transmission range as well as a guard interval to compensate forpossible clock drifts between periodic time synchronization. To quantify the latter,consider a clock skew S, and guard-interval of ? sec. Then, time-synchronizationis required every k?S sec, where the value of k depends on the time it takes to re-synchronize. Since in UWAC propagation delay is long and message rate is low(transmission rates on the order of a few kbit/s are common [8]), the latter is notexpected to affect performance much. However, due to the long time-slots, end-to-end transmission delay in TDMA scheduling might be too large in practice even forsmall number of network nodes [75]. A different approach is to apply spatial-reusetechniques where significant improvement in channel utilization is possible. Spatialreuse allows several nodes to share network resources such as frequency bands (inmulticarrier systems, e.g. [76]) or time slots (in TDMA systems, e.g. [77]), increasingchannel utilization [78]. Spatial-reuse TDMA is an appealing technique in UWANswhere low power transceivers, range and frequency dependent absorption loss [5],and acoustic NLOS scenarios lead to sparse network graphs. The concept of spatialreuse in UWANs was first introduced in [75], where a scheduling algorithm that im-proves channel utilization by clustering the network was suggested. Assuming shortintra-cluster distances, the communication within clusters is based on TDMA, whileeach cluster is assigned a unique pseudo-random spreading sequence used for sig-nal modulation to reduce interferences between adjacent clusters. In [79] the longpropagation delay in the channel was utilized to allow staggered packet transmis-sion. Unfortunately, the scheduling algorithm is based on the propagation delayof specific node-to-node links, which cannot be exploited in node-to-multiple-nodestransmission. Assignment of network resources (i.e., time slots) to maximize chan-nel utilization is known as the broadcast scheduling problem (BSP) [80], [81]. TheBSP can be formulated as a graph-coloring problem, and various heuristics to solve(variants of) the BSP have been proposed in e.g., [77, 82, 80]. However, these BSPformulations do not consider transmission of broadcast packets, which requires packetflow control via routing.16Chapter 1. Introduction1.2.2 Adaptive TransmissionsSeveral approaches have been suggested to utilize location or propagation delay infor-mation by opportunistically transmitting more data in the channel (e.g., [83, 68, 70]).Another method to utilize location information is through transmitter-side adjust-ment of the transmission scheme or by receiver-initiated request of additional trans-missions, i.e., automatic repeat request (ARQ), to ensure successful data deliv-ery. Considering the benefit of channel-dependent adjustments of the transmissionscheme, an adaptive modulation scheme was implemented in [84] to optimize trans-mission rate for time-varying channel conditions. With regards to reliability of trans-missions, incremental redundancy hybrid ARQ (IR-HARQ) is particularly efficientas it does not suffer from a coding loss due to repetition of the same parity symbols[85]. Despite the coding efficiency of IR-HARQ, it suffers from high latency dueto retransmission requests and retransmissions. Due to the long propagation delayof sound transmission and the low link reliability, this disadvantage is particularlypronounced in UWANs. In part addressing this problem, several adaptive transmis-sion applications tailored to UWANs have been suggested. In [86], an HARQ usingrateless codes has been suggested for transferring large files underwater. A ratelesscoding scheme is also used in [87] and [88] to optimize throughput of UWANs forbroadcast communications. Recently, [89] offered to optimize the code rate for thecurrent channel conditions by forming transmitter and receiver collaboration.1.3 Open Problems Addressed in this ThesisIn the previous section, we have reviewed current approaches for UWL and UWANscheduling. While for the former, challenges associated with deployment of anchornode, time-synchronization, and mobility have been addressed to some extent in thereviewed schemes, and for the latter collision avoidance scheduling algorithms thatcombat the long propagation delay in the channel have been suggested, in our view,the following challenges should be, but have not been, fully addressed.1.3.1 Challenges for Underwater LocalizationSound Speed Variation While most range-based localization techniques assume aknown speed of sound underwater, the dependency of the speed of sound with depth,temperature, and salinity, makes it challenging to pre-estimate. In this thesis we willshow that a mismatch of roughly 10 m/sec in the assumed sound speed considerablyaffects localization accuracy. While some works suggested measuring the parametersaffecting the sound speed, it is not an easy task for small and relatively simplevehicles, and a model-based approach may induce localization errors. We thereforebelieve that the sound speed should be estimated and tracked over time.17Chapter 1. IntroductionInter-node Time Synchronization Localization schemes that rely on silent posi-tioning to minimize communication overhead assume that nodes are time-synchronized.However, unlike surface nodes that can be time-synchronized via GPS updates, theclocks of submerged nodes are subject to skew as well as offset. Although timesynchronization algorithms have been proposed for UWANs, to reduce latency andcommunication overhead they should be incorporated into localization schemes.Node mobility model Node mobility due to ocean current, which presents one ofthe greatest challenges for UWL, has only been accounted for up to various degrees.Although some schemes assumes a simple mobility model, either the anchor nodesor the UL nodes are always assumed fixed during the localization process. Further-more, node mobility exhibits different characteristics and irregularities which makesit hard to determine the motion model. By using inertial systems, often availablefor navigation, and the (possible) spatial correlation of the ocean current, we believethat more information can be available to account for such motion irregularities.Impact of MAC Delays Another important challenge that has not been fullyaddressed for UWL is MAC to resolve contention. MAC schemes inevitably introducedelays in transmission, and affect the accuracy of localization schemes that rely onimmediate or scheduled responses (e.g., two-way messaging). Due to this delay andthe constant motion of nodes in the channel, it is not possible to assume fixed nodesfor UWL.Impact of Channel Structure Since range is measured based on the ToF of thedirect path or its received power, it is essential to lock on the location of the directpath of the received signal. Existing algorithms assume that the direct path is thestrongest path and thus it?s location is easy to estimate. However, multipath fadingcan lead to destructive interference and as a result the energy of the direct path isnot always the strongest. Moreover, the presence of structures and obstacles in theUAC may result in the loss of the direct signal. Hence, designated mechanisms toclassify ToF measurements into LOS and NLOS are needed.Effect of Ocean Waves Last, current UWL schemes assume capability to es-timate or track the orientation angle of the vehicle, and thus project accelerationmeasurements to the horizontal plain. However, near the ocean surface the vehiclemay experience time-varying pitch and roll angles, which may be too irregular totrack and too rapid to directly measure. Therefore, a tailored solution to this specificcase is required.18Chapter 1. Introduction1.3.2 Challenges for Scheduling of UWANsLocation-Dependent Handshake-based Scheduling The problem of low chan-nel utilization when small-to-moderate information unicast packets are sent in handshake-based scheduling can be overcome by utilizing the long propagation delay in thechannel through location information. While some of the reviewed methods exploitexposed terminal links, this is mostly done opportunistically and performance are farfrom the optimal throughput for UWANs.Time-Varying Topology Changes While spatial-reuse TDMA scheduling algo-rithms can improve network throughput, the problem of time-varying network topol-ogy together with the slow propagation of topology information in the UWAN, cangreatly decrease performance of such scheduling algorithms. On the other hand, dueto the long duration of the time-slots in TDMA UWAN scheduling, the low through-put of topology-transparent scheduling algorithms (like pure TDMA) may not besufficient for network requirements. Considering the benefit of high throughput oftopology-dependent scheduling and reliability of topology-transparent scheduling al-gorithms, a scheme combining both approaches is the natural next step.Utilizing Location Information for Adaptive Scheduling The long propaga-tion delay in the channel together with transmitter-receiver distance information offergreat opportunities for improving performance by means of adaptive transmissions.Intuitively, by setting the transmission parameters according to the expected delayin the channel, throughput can be optimized. While some works consider adaptivecoding for UWANs, a pure distance-dependent adaptive channel coding scheme hasnot been suggested.1.4 Thesis StructureIn this thesis, we present UWL and spatial-reuse scheduling algorithms for UWANsthat address the open problems identified in Section 1.3. As illustrated in Figure 1.6,the former serves as a building tool for the latter. We divide the thesis into two parts:I) UWL, and II) spatial-reuse scheduling for UWANs. In Chapter 2, we present ascheme for joint time-synchronization and localization, which considers sound speeduncertainties and utilizes motion in the channel to allow localization even when onlya single anchor node is available. Using this location information, in Chapter 3 wepropose a location tracking scheme that combats the effect of motion irregularities,tracks the sound speed, and incorporates Doppler shift measurements. Next, weconsider the problem of ranging and INS measurement errors, which affect localizationand tracking accuracy. The former problem is considered in Chapter 4, where wepresent a classification approach of ToF information into classes of NLOS and LOS.The latter challenge is the focus of Chapter 5, where, for the case of a TN whose19Chapter 1. IntroductionFigure 1.6: Structure of thesis.motion is affected by ocean waves, we suggest a machine-learning approach to projectacceleration measurements into the horizontal plane and perform DR without usingorientation measurements. The connections between the above chapters is illustratedin Figure 1.6.In Part II, we present spatial-reuse MAC techniques that rely on either locationor topology information available through the UWL capability developed in Part I.In Chapter 6, we show how location information can be used in a handshake-basedscheduling algorithm to utilize all available network resources even in fully connectednetworks. For broadcast UWAC, in Chapter 7 we present an optimal spatial-reuseTDMA scheduling scheme to trade off robustness to topology changes and networkthroughput. Based on such TDMA scheme, in Chapter 8 we propose an adaptivechannel coding technique that utilizes location information to greatly increase net-work throughput. Finally, in Part III we summarize our contributions and suggesttopics for further research.20Part IUnderwater Localization andTracking21Chapter 2UWL with Time-Synchronizationand Propagation SpeedUncertaintiesConsidering the challenges identified in Section 1.3.1, in this chapter we propose anew algorithm for UWL. In particular, our algorithm takes into account anchor andUL node mobility as well as propagation speed uncertainties, can function with onlyone anchor node, and includes time-synchronization of nodes. These abilities areimportant to enable localization under varying conditions, such as static or mobilenodes, in shallow or deep water, and when nodes are submerged for short or longperiods of time. Our setting includes several anchor nodes at known locations andone or more UL node, whose location is estimated. We assume that UL nodes areequipped with means to self-evaluate their speed and direction such as accelerometerand compass. Since such inertial systems are relatively light weight and inexpensive,there is a large variety of applications that satisfy this assumption for UL nodes.These include AUV, remotely operated underwater vehicles (ROV), manned vehi-cles, and divers [34]. Operating in tandem with self-localization systems, our algo-rithm makes use of the permanent movements of underwater nodes. It also performsa self-evaluation of localization accuracy, by estimating the propagation speed andchecking its validity, relying on known model boundaries for it. According to thestructure of the proposed algorithm we refer to it as sequential time-synchronizationand localization (STSL) algorithm. We demonstrate the advantages of the STSL al-gorithm by simulation comparisons with two benchmark localization methods, whichreveal significant localization errors for the latter if nodes are not time-synchronizedor the propagation speed is not accurately estimated. Furthermore, we formalize theCrame?r-Rao lower bound for UWL and show that it is well approached by the STSLalgorithm. Considering the problem of accurately modeling the UAC in a simula-tion environment, we also conducted a sea trial in August 2010 in Haifa, Israel, andpresent results that confirm the performance of the proposed algorithm under realconditions.The remainder of this chapter is organized as follows. In Section 2.1, we brieflysummarize the general structure of our algorithm and the intuition behind it. InSection 2.2, we introduce the system model, followed by a detailed description anddiscussion of our STSL algorithm in Section 2.3. Crame?r-Rao lower bounds perti-22Chapter 2. UWL with Time-Synchronization and Propagation Speed Uncertaintiesnent to our problem are derived in Section 2.4. Simulation and sea trial results arepresented and discussed in Sections 2.5 and 2.6, respectively. Finally, conclusions aredrawn in Section 2.7.2.1 IntuitionThe intuition behind our approach is the use of relative speed and direction infor-mation available at the mobile UL node to compensate for node mobility. In doingso, three or more range measurements obtained at different times and locations canbe combined for 2-D localization. This approach allows us to readily include thelocalization procedure as part of the operation of a communication network. Morespecifically, instead of using designated localization packet exchange (which is nec-essary if node mobility is not compensated), we rely on periodic packet exchangebetween the network nodes. This characteristic renders our approach more flexibleand easy to integrate into a UWAC system and also reduces communication overhead.The STSL algorithm uses a two-step approach, in which first nodes are time-synchronized and then location is estimated. In both steps, the measured time offlight of packets exchanged between anchor and UL nodes and self-localization dataobtained at UL nodes are linked to the unknown location, synchronization (clock skewand offset), and propagation speed parameters through linearized matrix equations.Our algorithm is modular in the sense that both time-synchronization and localizationsteps can be readily replaced with alternative solutions (as we do in this chapter tobenchmark the STSL performance).Before describing the STSL algorithm in detail, we next present the system modeland assumptions used in this work.2.2 System Setup and AssumptionsOur setting includes one or more UL nodes directly connected to L ? 1 anchor nodes,which have means to accurately measure their time-varying 2-D location and transmitit to the UL node. Both UL and anchor nodes operate in a time-slotted UWACnetwork, where nodes transmit at the beginning of globally established time slots asin, for example TDMA, slotted handshake [60] or slotted Aloha [90] transmission2.Since usually UL nodes perform localization independently of each other, we considerlocalization of one UL node in the following.We are interested in estimating the 2-D location of the UL node in terms of theuniversal transverse mercator (UTM) coordinates jN = [jxN , jyN ]T (the subscript Nbecomes clear below) after a pre-defined localization window of duration W time-slots, which without loss of generality starts at the UL node local time tUL = 0. We2We note that a relaxation of this assumption is possible by having anchor nodes time-stamptheir packets, thereby informing the UL node of the transmission time.23Chapter 2. UWL with Time-Synchronization and Propagation Speed Uncertaintiesassume that nodes are not time-synchronized. Suppose that the local time at the ULnode tUL corresponds to the time tl according to the local clock of anchor node l,and let Sl and Ol denote the clock skew and offset of the UL node relative to node l,l = 1, . . . , L, which are constant within the localization window. Then,tUL = tl ? Sl +Ol . (2.1)We assume that the time-synchronization error is small relative to the globally es-tablished time-slot duration, such that a node can match a received packet with thetime slot it was transmitted in. Thus, local transmission times of received packetsare known and time stamps of packets are not required. We note that long time slotsare common in UWAC [8] due to the low propagation speed underwater, which ismodeled between 1420 m/sec and 1560 m/sec [5].For localization we rely on a two-way packet exchange between the anchor nodesand the UL node. Let N be the total number of packets exchanged between the ULnode and the L anchor nodes during the localization window W . For convenience,we define the sets N a and N b for enumerating the packets to and from the UL node,respectively, such that N a ? N b = N = {1, . . . , N}. Denote li the index of theanchor node that transmits (i ? N a) or receives (i ? N b) the ith packet, and Ti andRi the transmission and reception local times of this packet. Also denote T pdi thepropagation delay for packet i according to local clock of anchor node li. (Note thatthe propagation delay according to the clock of the UL node is T pdi ? Sli). Considerpacket i ? N a transmitted at local time Ti and detected by the UL node at anchornode li local time Ti + T pdi + ?i, where ?i is a propagation-delay-measurement-noisesample. Also, consider packet i ? N b received by anchor node li at local time Ri + ?iand transmitted at anchor node li local time Ri + ?i ? T pdi . Then, following therelation in (2.1), the above time variables are related byRi = Sli(Ti + Tpdi + ?i) +Oli , i ? N a (2.2a)Ti = Sli(Ri + ?i ? Tpdi ) +Oli , i ? N b . (2.2b)Since for i ? N a the UL node measures Ri and is aware of Ti via packet association(recall we assume transmissions in globally established time slots), and since fori ? N b the UL node records Ti and is informed of Ri through communication withanchor node li, the UL node is able to construct equations (2.2a) and (2.2b). Formathematical tractability and formulating a practical algorithm, we assume that thenoise ?i is a zero-mean i.i.d. Gaussian random variable with variance ?2. Since morecomplicated models, such as mixture models with one component having non-zeromean, would likely be a more faithful noise representations, we study the effect ofmodel mismatch in Section 2.5.Considering a dynamic scenario in which all nodes permanently move either byown means or by ocean current, we assume that the UL node uses an inertial systemto self-estimate its speed and direction. During the localization window, the inertial24Chapter 2. UWL with Time-Synchronization and Propagation Speed Uncertaintiessystem provides N position estimations j?i = [j?xi , j?yi ]T for the true locations ji of theUL node at the time of transmission (i ? N b) or reception (i ? N a) of the ith packet.These locations are translated into a series of motion vectors, ?i,i? = [d?i,i? , ??i,i? ]T ,where d?i,i? and ??i,i? are the distance and angle between two self-estimated locationsj?i and j?i? , respectively. More specifically, assuming depth differences to be small(extensions are straightforward but not included here for brevity), the elements of asingle motion vector ?i,i? ared?i,i? = ?j?i ? j?i??2 , tan(??i,i?) =j?yi ? j?yi?j?xi ? j?xi?. (2.3)While we do not directly use the self-estimated locations j?i, whose errors accumulatewith time, we rely on the accuracy of the motion vectors for all packet pairs i, i?transmitted or received by the UL node during the localization window. That is, weassume that for i, i? ? N the estimated distance d?i,i? equals the true distance di,i? andthat ??i,i? equals the true angle ?i,i? . We note that this assumption sets limits on thevalue of W , which is determined by the specifications of the inertial system in use.We are now ready to present the STSL algorithm for UWL.2.3 The STSL AlgorithmWe are interested in accurately estimating the position jN of the UL node at theend of the localization window. In this section we first formalize the optimizationproblem for estimating jN , using ToA measurements obtained from received packetsand taking into account inertial system information. Then, we derive a sub-optimalsolution, namely the STSL algorithm, in which first nodes are time-synchronized andthen localization is performed.According to the system description in the previous section, the location pi ofanchor node li when transmitting (i ? N a) or receiving(i ? N b)the ith packet, thetransmission and reception local times Ti and Ri, respectively, and the motion vector?i,i? between locations ji and ji? are available at the UL node. Hence, denoting thepropagation speed c, the UWL problem can be formulated asj?N = arg minjN?i?N(?liTpdi ?1c?ji ? pi?2)2(2.4a)s.t. 1420m/sec ? c ? 1560m/sec (2.4b)Ri = Sli(Ti + Tpdi ) +Oli , i ? N a (2.4c)Ti = Sli(Ri ? Tpdi ) +Oli , i ? N b (2.4d)di,i? = ?ji ? ji??2 , i, i? ? N (2.4e)tan(?i,i?) =jyi ? jyi?jxi ? jxi?, i, i? ? N , (2.4f)25Chapter 2. UWL with Time-Synchronization and Propagation Speed Uncertaintieswhere the factor ?li in (2.4a) accounts for the difference between Tpdi (according toanchor node li local clock) and the actual propagation delay of packet i. Note thatproblem (2.4) is different from conventional localization problems due to the time-synchronization constraints (2.4c) and (2.4d), and due to the unknown sound speedc. Since (2.4) is a non-convex problem, we device our STSL algorithm as a pragmaticsolution to the localization problem at hand, and we will compare its performance tothe Crame?r-Rao lower bound (CRLB) associated with estimating the desired locationjN as well as the unknown parameters Sl and Ol, l = 1, . . . , L, and c.In the following we describe the details of our STSL algorithm starting from thetime-synchronization step and followed by the localization step.2.3.1 Step 1: Time-SynchronizationThe objective of the time-synchronization step is to provide estimates of the prop-agation delays T pdi , ?i ? N . This is accomplished by two-way packet exchange,obtaining equations of type (2.4c) and (2.4d). However, due to the permanent mo-tion of nodes in the channel, propagation delays T pdi , i ? N a, and Tpdj , j ? N b,might not be equal, and thus (2.4c) and (2.4d) cannot be readily compared. Com-mon time-synchronization approaches for UWAC deal with this problem by lettingthe receiving node respond immediately to limit any possible movements (e.g., [14]).However, such a requirement limits the scheduling protocol. We choose a differentapproach and apply a quantization mechanism to allow for differences in the propa-gation delay of separate packets and to enable time-synchronization per anchor nodemaking use of the ongoing network communications.Quantized Representation of Node MovementsIn the quantization step, the locations of the UL node and the anchor nodesare quantized so that multiple ToA measurements from two-way communication areassociated with the same pair of quantized locations. More specifically, consider thetwo packets n,m, n ? N a, m ? N b. If the two sets of UL node locations jn, jmand anchor node locations pn, pm with ln = lm = l are associated with the samequantized location k? and ul,? of the UL node and anchor node l, respectively, weassume that T pdn = T pdm and (2.2a) and (2.2b) can be combined as we show furtherbelow. The variables ? and ? are used to enumerate quantized locations.To quantize the locations of anchor nodes, we introduce subsets Ul,? ? N includingall packets associated with the same anchor node l such that for each pair of packetsn,m ? Ul,? , ?pn ? pm?2 < ?, where ? is a fixed threshold. Next, we associatelocation pi, i ? Ul,? , with the quantized location ul,? . Similarly, to quantize locationsof the UL node we form subsets of packets K? ? N such that for each pair of packets26Chapter 2. UWL with Time-Synchronization and Propagation Speed Uncertaintiesn,m ? K?, dn,m < ?, and associate location j?i, i ? K?, with the quantized locationk?. We note that a single packet can be associated with multiple subsets Ul,? and K?.There is a tradeoff for choosing ?. If ? is too large, the assumption of identicalpropagation delay is notably flawed and thus the accuracy of the time-synchronizationprocess is low. If ? is too small, there might not be enough two-way ToA measure-ments associated with each pair of quantized locations ul,? and k?, and again accuracyof the time-synchronization process is degraded, as we further discuss below.Estimating the Clock Skews and OffsetsWe now use the quantized locations to estimate clock skews Sl and offsets Ol,l = 1, . . . , L. Let us define subsets N al ? N a and N bl ? N b, with cardinality Nal andN bl , respectively, including all packets associated with anchor node l. Consider thepair of packets n,m, n ? N al , m ? N bl , for which locations pn and pm are mappedonto the same quantized location ul,? , and locations j?n and j?m are mapped onto thesame quantized location k?. We assume that for each anchor node l, this mappingresults into Ml pairs of equations (2.2a) and (2.2b). Clearly, Ml increases with thequantization threshold, ?. As stated above, we neglect the differences between thepropagation delays T pdn and T pdm in (2.2a) and (2.2b) and thus obtain Ml equationsof the formRn + TmSl? 2OlSl= Tn +Rm + ?n + ?m, n ? N al , m ? N bl . (2.5)Note that equations of type (2.5) are introduced separately for each anchor node l.This is because the estimated clock skew and offset are different for each l.Introducing the variable vector ?l = [?l(1),?l(2)]T = [ 1Sl ,OlSl]T , we express (2.5)as the linear matrix equationBl?l = bl + \u000Fl (2.6)for each anchor node l, where Bl is an [Ml ? 2] matrix with rows [Rn + Tm,?2],and bl and \u000Fl are column vectors of appropriate length with elements Tn + Rm and?n + ?m, respectively, with n ? N al , m ? N bl . Next, we apply the LS estimator??l =(BTl Bl)?1BTl bl (2.7)for each anchor node l. By (2.7), the covariance matrix of ??l is [91]Q? = ?2(BTl Bl)?1 , (2.8)whose main diagonal elements are proportional to 1Ml and1M2l, respectively. Hence,for largeMl the estimates ??l(1) and ??l(2) are expected to have much smaller variancethan ?2.27Chapter 2. UWL with Time-Synchronization and Propagation Speed UncertaintiesEstimating Propagation DelaysAfter estimating ?l(1) and ?l(2), the quantized locations are no longer in use andwe return to our initial objective, which is to estimate the propagation delay. Thus,localization accuracy of the STSL algorithm is not limited to ?. Considering (2.2),for packets n ? N al and m ? N bl , the UL node estimates the propagation delay asT? pdn = Rn??l(1)? ??l(2)? Tn, n ? N alT? pdm = ??l(2)? Tm??l(1) +Rm, m ? N bl . (2.9)We observe from (2.9) that the propagation delay estimation error is a functionof both ToA measurement error, ?i, and clock skew and offset estimation errors.However, since during the localization window Rn and Tn are bounded by W , thevariances of Rn??l(1)? ??l(2), n ? N al , and ??l(2)? Tm??l(1), m ? N bl , are expected tobe much smaller than ?2. Thus, we use the approximation T? pdi = Tpdi + ?i in thefollowing.2.3.2 Step 2: LocalizationWe now introduce the localization step of the STSL algorithm. This step is performedimmediately after time-synchronization, using propagation delay estimations (2.9).The objective of the localization step of the STSL algorithm is to estimate the ULnode UTM coordinates jxN and jyN at the end of the localization window W . For thispurpose, we adopt the common approach to linearize the estimation problem [92],and first estimate the transformed variable vector ?N =[(jxN)2 + (jyN)2 , jxN , jyN]T.Define?i,i? =d?i,i??1 + tan(??i,i?)2?i,i? = ?i,i? tan(??i,i?) , (2.10)and assume d?i,i? and ??i,i? in (2.3) to be equal to di,i? and ?i,i? from (2.4), respectively(recall that we rely on the accuracy of the motion vectors during the localizationwindow). Thus,jxi? = jxi ? ?i,i? , jyi? = jyi ? ?i,i? , i, i? ? N . (2.11)Furthermore, we have from (2.4) thatT? pdi =1?lic?ji ? pi?2 + ?i , i ? N . (2.12)Since c is unknown, and assuming small differences between ?l such that ?l?l??= 1, l, l? =1, . . . , L (a relaxation of this assumption is given further below), we reduce the set ofN equations (2.12) to N ? 1 equations, which together with (2.11) can be written as?N,i ? ?N = aN,i + \u000FN,i, i = . . . , N ? 1 (2.13)28Chapter 2. UWL with Time-Synchronization and Propagation Speed Uncertaintieswith vector ?N,i = [?N,i(1), ?N,i(2), ?N,i(3)], whereaN,i =(T? pdi)2 ((pxN)2 + (pyN)2)?(T? pdN)2((pxi + ?N,i)2 + (pyi + ?N,i)2),?N,i(1) =(T? pdN)2?(T? pdi)2,?N,i(2) = 2(T? pdi)2pxN ? 2(T? pdN)2(pxi + ?N,i) ,?N,i(3) = 2(T? pdi)2pyN ? 2(T? pdN)2(pyi + ?N,i) , (2.14)and \u000FN,i is the noise component originating from the noisy estimations (2.9). For thelocalization window W , we construct an [(N ? 1)? 3] matrix A with rows ?N,i andvectors a and \u000F with elements aN,i and \u000FN,i, respectively. Then, the (N?1) equations(2.13) are arranged inA?N = a+ \u000F . (2.15)The elements of the error vector \u000F depend on the elements of ?N . Thus, directestimation of ?N from (2.15) will result in low accuracy. Hence, we follow [92] andoffer a two-step heuristic approach in which first we get a coarse estimate of ?N , andthen we perform a refinement step. The coarse estimate is given by??LSN =(ATA)?1Aa . (2.16)We note that \u000FN,i from (2.13) can be formalized as ?ifN,i, where fN,i is a function ofthe elements of ?N , not given here for brevity. Thus, \u000FN,i are i.i.d random variablesand the covariance matrix ?2QN of \u000F is a diagonal matrix whose ith diagonal elementequals ?2f 2N,i. Using ??LSN from (2.16) to estimate the elements of fN,i, i = 1, . . . , N?1,we estimate QN as Q?N . The refined estimate of ?N follows as??WLSN =(AT Q??1N A)?1AQ??1N a , (2.17)with the error covariance matrix [91]Q??N =(AT Q??1N A)?1. (2.18)Finally, we use the inner connection of the elements of ?N to estimate the locationvector jN . Defining GN =????WLSN (2) ??WLSN (3)1 00 1??, where ??WLSN (i) is the ith element of??WLSN , we obtainGNjN = ??WLSN + \u000FN , (2.19)29Chapter 2. UWL with Time-Synchronization and Propagation Speed Uncertaintieswhere \u000FN is a [3? 1] estimation noise vector of ??WLSN . Using (2.19), Q??N from (2.18)and ??WLSN from (2.17), the WLS estimator of jN isj?N =(GTNQ???1N GN)?1GNQ???1N ??WLSN , (2.20)whose elements j?N(1) = j?xN and j?N(2) = j?yN are the desired location coordinates.We would like to mention that if assumption ?l?l??= 1 used to obtain (2.13) does nothold, the localization process can be performed on a per-anchor-node basis. To thisend, packet index i in (2.13) is limited to packets transmitted or received by the sameanchor node, and the number of equations (2.13) is reduced to N?1L (assuming equalnumber of transmissions pre anchor node in the network). Since L is expected to besmall (we use L = 2 in our simulations and sea trial described below), the accuracyof the localization process is not expected to deteriorate much. Then, the UL nodelocation at the end of the localization window can be estimated by combining per-anchor-node based estimations j?N from (2.20) using data fusion techniques, cf., [93].Since per-anchor-node based estimations j?N are independent of ?l, such combinationis not affected by mismatch of ?l between anchor nodes.2.3.3 ExtensionsIn this section we introduce two extensions for the above location estimation. Thefirst is a refinement step in which we iteratively improve the location estimation(2.20). The second is a self-evaluation process to test the accuracy of the localizationprocess.Iterative RefinementThe accuracy of estimation (2.20) depends on the quality of the coarse estimate??LS from (2.16), used to construct the error covariance matrix, Q?N . We now follow[94] and propose an iterative refinement procedure in which the accuracy of Q?N isimproved.In the kth step of our iteration, vector j?N,k is estimated using (2.20) from whichthe vector ??N,k is constructed. Next, in the (k + 1)th step ??N,k replaces ??LSN in theconstruction of Q?N . As a stopping criterion, we use the covariance matrix of the kthestimation (2.20),Q???N,k =(GTNQ???1N GN)?1. (2.21)Since the determinant, |Q???N,k|, is directly proportional to the estimation accuracy [91],the iteration stops when the absolute value of |Q???N,k| ? |Q???N,k?1| is below some em-pirically chosen threshold, ?iter, or if the number of iterations exceeds its maximum,Niter. While we could not prove the convergence of this process, we demonstrate itby means of numerical simulations in Section 2.5.30Chapter 2. UWL with Time-Synchronization and Propagation Speed UncertaintiesSelf-Evaluation of Localization PerformanceIn this section, we describe a binary test for self-evaluating localization accuracy.It can be used to adjust STSL parameters, such as the localization window W , forrefining the localization procedure, such as data fusion of per-anchor-node basedlocalization (see discussion after (2.20)), or to decide whether an UL node can be usedas a new reference node. For the latter application localization should be extended totracking though, to make sure that location estimates remain accurate when nodesmove. Our self-evaluation test relies on a widely used model that bounds propagationspeed underwater between 1420 m/sec and 1560 m/sec [5]. In particular, given anestimate of the propagation speed, cest, the binary test output ? is computed as? ={1, if 1420 ? cest ? 15600, otherwise}, (2.22)and ? = 1 and ? = 0 indicate accurate and non-accurate localization, respectively.Using (2.12) and since ?l are expected to be close to 1, we obtain the propagation-speed estimate ascest =1NN?i=1?j?i ? pi?2T? pdi, (2.23)where j?i, i = 1, . . . , N ? 1, follow from j?N in (2.20) using relation (2.11). Differentfrom traditional self-evaluation techniques involving the broadcast of a confidenceindex obtained from comparing the measured propagation delay to the estimatedone (e.g. [95, 30]), the advantage of the our method lies in the comparison of cestto a given model of propagation speed, which is independent of the estimation. Ournumerical results (see Section 2.5) show that when localization error is accurate (i.e.,below 10 m), we obtain ? = 1 in more than 99% of the cases, and when localizationerror is non-accurate (i.e., above 10 m), ? = 0 results in 90% of the cases. If morereliable evaluation performance is needed, the proposed test could be combined withother self-evaluation tests.2.3.4 ScalabilityWhen more network nodes are added, often fewer packets are transmitted per node.Hence, performances of both time-synchronization and localization degrade with in-creasing number of UL nodes, NUL. Furthermore, since time-synchronization is per-formed per-anchor node, its performance degrades with increasing number of anchornodes, L. However, since the number of propagation delay measurements, avail-able for localization, increases with L, performance of the localization step by itselfimproves with L. However, since localization depends on the output of the syn-chronization step, the overall performance may improve or degrade with increasingnumber of anchor nodes. In summary, scalability of the STSL algorithm is closelyrelated to the scalability of the underlying communications protocol.31Chapter 2. UWL with Time-Synchronization and Propagation Speed Uncertainties2.3.5 Pseudo-Code for STSLThe operation of the STSL algorithm is summarized in the pseudo-code in Algo-rithm 1. For simplicity, the quantization mechanism introduced in Section 2.3.1 isnot included, and we start when positions are already quantized into locations ul,? andk?. First, equations (2.5) are formed, and an LS estimator is used to estimate clockskew and offset for each anchor node l (lines 2-9). Then, the time-synchronizationstep is concluded by estimating propagation delays for each transmitted or receivedpacket (line 12). The localization step begins with forming equations (2.13) (line 13),followed by an initial LS estimator (line 15). Then, an iterative procedure beginswhere in each step the covariance matrix Q?N and location j?N,k are estimated (lines18-19). The latter is then used to refine the initial estimation by iteratively formingmatrix Q?N till convergence is reached (lines 19-23). The algorithm performs a seriesof LS and WLS estimations with complexity of O (N3 +N) and is executed only onceat the end of the localization window. A software implementation of the algorithmcan be downloaded from [96].2.4 Crame?r-Rao Lower BoundFor the purpose of gauging the performance of the STSL algorithm, in this sectionwe develop analytical expressions to lower bound the performance of any unbiasedUWL estimator, assuming nodes not to be time-synchronized and propagation speedunknown. We start with general expressions for the CRLB, and then apply it to ourspecific localization problem.2.4.1 General Crame?r-Rao Lower BoundConsider a measurement vector y = h(pi,?) + n, where n is a noise vector, andh(pi,?) is some function of a vector of wanted variables, pi, and a vector of nuisancevariables, ?. For an unbiased estimator, the variance of the nth element of pi, pin,can be bounded by the CRLB [97]E[(p?in ? pin)2]? CRB(pin) , (2.24)where CRB(pin) =(I?1)n,n and I is the Fischer information matrix (FIM), whose(n,m)th element isIn,m = ?Ey[?2 lnP (y|pi)?pin?pim], (2.25)and P (y|pi) is the probability density function of y given pi. To calculate P (y|pi) oneneeds to average the nuisance variables, ?, i.e., P (y|pi) = E? [P (y,?|pi)], which makesit hard to calculate (2.25), since often P (y,?|pi) cannot be expressed. Therefore,32Chapter 2. UWL with Time-Synchronization and Propagation Speed UncertaintiesAlgorithm 1 Estimate jN1: {Step 1: Time-synchronization}2: for (l = 1, . . . , L) do3: for (n ? N al , m ? N bl ) do4: if (pn,pm ? Ul,?) ? (j?n, j?m ? K?) then5: Form equations (2.5) using Tn, Rn, Tm and Rm6: end if7: end for8: Estimate Ol and Sl using (2.7)9: end for10: {Step 2: Localization}11: for (i = 1, . . . , N) do12: set T? pdi using (2.9)13: Form equations (2.13)14: end for15: Estimate ??N,1 using (2.16)16: for (k = 1 to Niter) do17: Estimate Q?N using ??N,k18: Estimate j?N,k using (2.20)19: Construct matrix Q???N,k using (2.21)20: if (|Q???N,k| ? |Q???N,k?1| ? ?iter) then21: Return22: end if23: Construct ??N,k+1 using j?N,k24: end forinstead of CRB(pin), the modified Crame?r-Rao bound MCRB(pin) =(I??1)n,nisoften used [98], whereI?n,m = ?Ey,?[?2 lnP (y|pi,?)?pin?pim]. (2.26)In [98] it was shown thatCRB(pin) ? MCRB(pin) . (2.27)Hence, MCRB(pin) may be too loose to compare with.A different approach would be to consider the nuisance variables ? as part of theestimation problem. That is, we consider a new variable vector ? = [piT ,?T ]T andformalize CRB(?n) forIn,m = ?Ey[?2 lnP (y|?)??n??m]. (2.28)33Chapter 2. UWL with Time-Synchronization and Propagation Speed UncertaintiesWe note that CRB(?n) does not bound E[(p?in ? pin)2]but E[(??n ? ?n)2]. Thus,it can only serve as a lower bound for estimators which estimate both pi and ?.2.4.2 Application to STSLSince the STSL algorithm includes a sequence of LS and WLS estimators, it is anunbiased estimator. Thus, we next apply the MCRB(pin) and CRB(?n) bounds forour STSL algorithm. We consider the measurement vector in (2.2) for which y = Ri,pi = [jxN , jyN ], ? = [S1, . . . , SL, O1, . . . , OL, c], and n is as ?i in (2.2). Then, we haveE[(j?xN ? jxN)2+(j?yN ? jyN)2]? CRB(pi1) + CRB(pi2) . (2.29)We note that the variance of y depends on the clock skew, Sli . Thus, although ?i isassumed Gaussian, the often used simplification for the CRB in the Gaussian case(cf. [97]) cannot be used to solve (2.26) and (2.28).Following our discussion in Section 2.4.1, we consider the alternative CRLB for-mulationsCRB = CRB(?1) + CRB(?2) , MCRB = MCRB(pi1) + MCRB(pi2) , (2.30)and compare?MCRB and?CRB to?err =?E[(j?xN ? jxN)2+(j?yN ? jyN)2]. (2.31)2.5 Simulation ResultsIn this section, we present and discuss simulation and sea trial results demonstratingthe performance of the STSL algorithm in different environments. We conducted10, 000 Monte-Carlo simulations of a scenario with two anchor nodes and one ULnode, communicating in a simple TDMA fashion. The three nodes were placeduniformly in a square area of 1 ? 1 km2 and moved between two adjacent packettransmission times at uniformly distributed speed and angle between [?5, 5] knotsand [0, 360] degrees, respectively. We added a zero mean i.i.d. Gaussian noise withvariance ?2 to each of the ToA estimations [see (2.2)]. Furthermore, considering theresults in [99] we added a zero mean i.i.d. Gaussian noise with variance 1 m2 toeach of the distance elements of the motion vectors [see (2.3)] while regarding theirangle components to be accurate. To simulate time-synchronization errors the clockof each of the three nodes had a Gaussian distributed random skew and offset relativeto a common clock with mean values 1 and 0 sec and variances 0.001 and 0.5 sec2,respectively.34Chapter 2. UWL with Time-Synchronization and Propagation Speed Uncertainties30 35 40 45 5010?210?11001/?2 [dB]? err[m] STSLSTSL?mixMultilaterationJLSCRBMCRBFigure 2.1: ?err from (2.31) as a function of 1/?2. Sound speed is known and all nodesare time-synchronized. Vertical bars show 95% confidence intervals of the simulationresults for STSL.We used a quantization threshold ? = 38 meters and a localization window ofW = 20 time-slots. The time-slot duration was selected Tslot = 5 seconds, consideringthe long propagation delay in the UAC (e.g., 4 sec for a range of 6 km). We comparethe performance of the STSL algorithm with those of the multilateration method [27]and the method proposed in [92], which we refer to as the joint localization and syn-chronization (JLS) algorithm. Both benchmark methods use an assumed propagationspeed c?. Furthermore, while the JLS algorithm performs joint time-synchronizationand localization (assuming anchor nodes are time-synchronized but the UL node isnot), the multilateration method assumes all nodes to be time-synchronized. Sinceboth benchmark methods assume static nodes, we used a different simulation envi-ronment for them such that a fair comparison with the STSL algorithm is possible.The simulation environment for the benchmark methods considers fixed nodes andadds virtual anchor nodes according to node movements in the original simulationscenario (i.e., the one used to test the performance of the STSL algorithm). Consider,for example, an anchor node l moving between locations pl1 and pl2 while commu-nicating with a static UL node. To test the benchmark methods, such a scenariowould change into a scenario where two static anchor nodes, l1 and l2, are locatedat pl1 and pl2 , respectively. Allowing a fair comparison between the three testedlocalization methods, the virtual anchor nodes, l1 and l2, have the same local clockas that of the real anchor node l. The implementation code of the STSL algorithmcan be downloaded from [96].First, we consider a scenario where c = c? = 1500 m/sec and all nodes are time-synchronized. Figure 2.1 shows ?err from (2.31) as a function of 1?2 for the three35Chapter 2. UWL with Time-Synchronization and Propagation Speed Uncertainties0 0.1 0.2 0.3 0.4 0.5 0.6 0.710?310?210?1100101102esync [%]? err[m] STSLMultilaterationJLSFigure 2.2: ?err from (2.31) as a function of esync. Sound speed is known and 1?2 =46 dB. Vertical bars show 95% confidence intervals of the simulation results for STSL.methods and the CRB and MCRB from (2.30). For clarity, here and in the followingwe show 95% confidence intervals in error bars only for the STSL algorithm. Theresults show that both benchmark methods achieve better performances than theSTSL algorithm. This is mainly because STSL redundantly estimates c as well asclock offsets and skews, which introduces errors. This is also why the multilaterationmethod achieves slightly better performance than the JLS protocol method. We notethat the MCRB is slightly lower than the CRB and both bounds are quite close tothe STSL error, which implies that although STSL is a heuristic estimator it achievesgood localization results. To show the effect of possible mismatch of our model forthe measurement noise ?i (see (2.2)), Figure 2.1 also includes results for STSL-mix,in which ?i is modeled as a mixture of two distributions. The first distribution (withweight 0.9) is a zero mean Gaussian with variance ?2 and the second (with weight0.1) is a Rayleigh(?) distribution, which accounts for multipath propagation [86].From Figure 2.1, we observe that the performance of STSL-mix decreases comparedto that in the Gaussian-noise case, which is mainly due to the non-zero mean of noise.However, this degradation is fairly moderate demonstrating some robustness of STSLto model mismatch.In Figure 2.2 we compare ?err for the three methods when c = c? = 1500 m/sec, butnodes are not time-synchronized and 1?2 = 46 dB, as a function of esync =S?W+O??WW ,where S? and O? is the average of Sl and Ol, l = 1, . . . , L, respectively. While theperformance of the STSL algorithm is hardly affected by the synchronization error(compared to the results in Figure 2.1), the JLS protocol method, designed for time-synchronized anchor nodes, and the multilateration method suffer from significantestimation errors even for small synchronization errors.36Chapter 2. UWL with Time-Synchronization and Propagation Speed Uncertainties0 10 20 30 40 50 60 70 8010?310?210?1100101102|c - c?| [m/sec]? err[m] STSLMultilaterationJLSFigure 2.3: ?err from (2.31) as a function of |c? c?|. All nodes are time-synchronizedand 1?2 = 46 dB. Vertical bars show 95% confidence intervals of the simulation resultsfor STSL.We now compare performance when c is chosen with uniform distribution betweenthe model boundaries, 1420 m/sec and 1560 m/sec, and the two benchmark methodswere still given the nominal value c? = 1500 m/sec. To understand the effect ofmismatched propagation-speed information on localization accuracy, we compare ?errfrom (2.31) when all nodes are time-synchronized. The results are shown in Figure 2.3as a function of |c ? c?|, again for 1?2 = 46 dB. We observe that for both benchmarkmethods, ?err dramatically increases even for a small difference of |c? c?| = 10 m/sec,which motivates the need to accurately estimate c in UWL. Furthermore, comparedto the results of Figure 2.1 the STSL is almost unaffected by the variations of c.Next, we consider the practical case where all nodes are not time-synchronized(same scenario as for Figure 2.2) and c is unknown (same scenario as for Figure 2.3).We study two of the properties of the STSL algorithm, namely the convergence ofthe refinement iterative process discussed in Section 2.3.3 and the self-evaluationprocess discussed in Section 2.3.3. In Figure 2.4, we demonstrate the convergence ofthe refinement iterative process by showing ?err from (2.31), averaged over all clockoffsets and skews and c instances, as a function of the number of iteration steps andseveral values of 1?2 . The results indicate that a significant performance improvementis achieved after only a few iteration steps. In Figure 2.5 we show the empiricalprobability density function (PDF) of the estimated propagation speed, cest, from(2.23) when c = 1500 m/sec. The results are shown for two cases: 1) when ?err ? 10 mand 2) when ?err ? 10 m. The results show that for small values of ?err, in more than99% of the cases, cest is inside the model boundaries (i.e., 1420 ? cest ? 1560), witha standard deviation of less than 10 m/sec. For large values of ?err, cest seems to be37Chapter 2. UWL with Time-Synchronization and Propagation Speed Uncertainties1 2 3 4 5 6 7 8 9 1010?210?1100101102? err[m]Number of iterations 1/?2 = 30 dB1/?2 = 40 dB1/?2 = 50 dBFigure 2.4: ?err from (2.31) for the STSL algorithm as a function of number ofiteration steps.almost uniformly distributed, with only 10% of the estimations being inside the modelboundaries. However, for some applications (e.g., localization in sparse networks) thismissed-detection probability may be too large. Thus, we conclude that cest can serveas a good indicator to confirm accurate localization, but may be used to complementother self-evaluation techniques to identify non-accurate localization.Finally, in Figure 2.6 we consider the same scenario as for Figure 2.4 and show ?errfrom (2.31) as a function of 1?2 . For clarity, since results for STSL are similar to thoseshown in Figure 2.1, error bars are omitted. We observe that while both benchmarkmethods suffer from a significant error floor, the error for the STSL algorithm de-creases with 1?2 and is the same as in Figure 2.1. Hence, the algorithm compensates forboth synchronization and propagation speed uncertainties. To demonstrate the rela-tion between the number of anchor nodes, L, and the number of UL nodes, NUL (seediscussion in Section 2.3.4), in Figure 2.6 we also include results for L = 3, NUL = 1(STSL, L = 3) and L = 2, NUL = 2 (STSL, 2UL), which had similar standard de-viation to the case of L = 2, NUL = 1. We note that since multilateration doesnot include time-synchronization, and since JLS performs joint time-synchronizationand localization, clearly their performance improves with L and degrades with NUL.Results show that, as expected, also performance of the STSL algorithm slightlyimproves with L and decreases with NUL.38Chapter 2. UWL with Time-Synchronization and Propagation Speed Uncertainties1100 1200 1300 1400 1500 1600 1700 180000.10.20.30.40.50.60.70.80.9c?PDFofc? ?err< 10m?err>10mFigure 2.5: Empirical PDF of cest for c = 1500 m/sec.2.6 Sea Trial ResultsIn this work we assumed 1) node?s clock skew and offset are time-invariant withinthe localization window, 2) propagation speed is time and space invariant for smalldepth differences, 3) propagation delay measurements are affected by a zero-meanGaussian noise, and 4) node movements are relatively slow such that quantization ofnode locations is possible. While the first assumption depends on the system clock,the second and the third depend on the channel. To verify our assumptions andconfirm our results we tested the STSL algorithm in a sea trial along the shores ofHaifa, Israel in August 2010.The sea trial included three drifting vessels, representing three mobile nodes, andlasted for Texp = 300 minutes. In Figure 2.7, we show the UTM coordinates of thenodes during the sea trial. We note that node 3 needed to turn on its engines aroundtime slot 150, which explains the sudden change in its direction and speed. Eachnode was equipped with a transceiver, deployed at 10 meters depth, allowing UWACat 100 bps with a transmission range of 5 km. The nodes communicated in a TDMAnetwork with a time-slot of Tslot = 60 seconds, allowing significant node motionbetween transmission of each packet. Time-slot management was performed at eachnode using an internal clock. These internal clocks were manually time-synchronizedat the beginning of the experiment with an expected clock offset of up to one second.We also note that pre-testing of these clocks showed a clock skew of one second perday.We used GPS receivers as reference for the location of each node as well as itsinertial system to obtain motion samples [see (2.3)]. The localization error of theGPS-based reference locations was reported to be uniformly distributed between 039Chapter 2. UWL with Time-Synchronization and Propagation Speed Uncertainties30 35 40 45 5010?210?11001011021/?2? err[m] STSLSTSL,L=3STSL,2ULMultilaterationJLSCRBMCRBFigure 2.6: ?err from (2.31) for time-synchronization and sound speed uncertainties.and 10 m. To test the effect of this uncertainty in the anchor-node location weconducted simulations similar to the scenario considered in Figure 2.1 with error-freeToA measurement but with anchor-node location uncertainties similar to those ofthe GPS receivers in use. The results showed that using our STSL algorithm suchuncertainty results in an average estimation error of 15 m. Thus, in the sea trial anylocation error below 15 m is considered accurate.2.6.1 Channel and System CharacteristicsAt the beginning and end of the sea trial we measured the propagation speed in waterusing a measuring probe. Both measurements showed that the propagation speed cwas bounded in between 1552 m/sec (for depth of 40 m) and 1548 m/sec (for depthof 1 m) and was on average 1550 m/sec. The small variance of the measurementsof c confirms our assumption that the propagation speed can be considered fixedthroughout the localization window. In the following, for performance evaluation weconsider c = 1550 m/sec, which is within the boundaries of our model (see Section 2.2)but is different from the commonly used value of 1500 m/sec.In Figure 2.8 for each pair of nodes we show T? pddiff , which is the time differencebetween propagation delay estimations at both sides of the communication link in asingle set of receiver-transmitter quantized locations, measured directly from (2.2a)and (2.2b) neglecting the clock skew and offset. For example, for a two-way packettransmission between the quantized locations k? and ul,? with propagation delayestimations T? pd1 and T?pd2 , T?pddiff = |T?pd1 ? T?pd2 |. We note that if nodes are time-synchronized, we would expect T? pddiff to be on the same order of the length of theimpulse response, which was measured as 20 msec on average and did not exceed40Chapter 2. UWL with Time-Synchronization and Propagation Speed Uncertainties20004000600010002000300040000100200300 UTM X (offset by 681000)UTM Y (offset by 3641000) MinutesNode 1Node 2Node 3Figure 2.7: Time-varying location of nodes in the sea trial.30 msec. However, the results show that T? pddiff increases with time and is much greaterthan 30 msec. This implies that nodes suffered from considerable clock skew andoffset. Furthermore, since the values of T? pddiff are different for each pair of nodes,the nodes skew and offset are different, which confirms with our system model (seeSection 2.2).2.6.2 ResultsIn the following, we compare the performance of the STSL algorithm in the sea trialwith that of a method aimed to solve a relaxed sequential time-synchronization andlocalization (R-STSL) problem, in which an a-priori propagation speed c? is given. Inthe R-STSL, time-synchronization is performed similar to the process described inSection 2.3.1, but the localization process is modified as c is known. The results areshown for all three nodes, where each time a different node was considered as the ULnode and the other two nodes were the anchor nodes. We measure the performancein terms of the Euclidean distance between the estimated location and the referenceGPS location, averaged over a sliding localization window of W time-slots, i.e.,??err =1TexpTslot?W + 1TexpTslot?n=W?j?n ? jn?2 , (2.32)where for each UL node location estimation j?n we used ToA and inertial systemmeasurements from time-slot n?W + 1 till n.In Figure 2.9, we demonstrate the effect of mismatched propagation speed, i.e.,c? 6= c, by showing ??err from (2.32) for W = 30 time-slots as a function of |c ? c?|,where c = 1550 m/sec. We note that although such choice of W seems large due41Chapter 2. UWL with Time-Synchronization and Propagation Speed Uncertainties0 50 100 150 200 250 30000.511.522.53Time slotTpd diff [sec] Link (1,2)Link (1,3)Link (2,3)Average increase: 1.2msec / time?slotAverage increase: 13.5msec / time?slotAverage increase: 14.3msec / time?slotFigure 2.8: T? pddiff for all communication links as a function of time slots.to the long time slot duration, the number of transmissions for each node was only10, which is in the same order as considered in our simulations. The results showthat ??err significantly increases with |c ? c?| even for a relatively small difference of10 m/sec.This result, as well as the results in Figure 2.3, validate the need to consider thepropagation speed as an additional variable in UWL. In the following we consider amatched version of R-STSL (MR-STSL), i.e., when c? = c, which in the absence ofbenchmark localization methods that take into account time-synchronization uncer-tainties and availability of inertial system to track short-term node movements, canserve as a lower bound for the STSL.Finally, in Figure 2.10 we show the empirical cumulative density function CDF of??err, averaged over the three nodes, for STSL and MR-STSL andW = 10, 20, 30 time-slots, i.e., in a single localization window an average number of packet transmissionsof 3.3, 6.6, 10 for each node, respectively. We observe that both mean and varianceof ??err improve with W , however, at a cost of delay. We observe that since STSLestimates an additional variable, its performance is worse than that of MR-STSL.However, the difference is not significant. We note that the average ??err for STSLand W = 30 time-slots is 21.5 m, which is close to the expected localization accuracydue to the GPS location uncertainties. Therefore, STSL fully compensates the largeclock skew and offset shown in Figure 2.8, node movements and propagation speeduncertainty.42Chapter 2. UWL with Time-Synchronization and Propagation Speed Uncertainties0 20 40 60 80 100050100150200250300350|c - c?| [m/sec]?? err[m] Node 1Node 2Node 3Figure 2.9: ??err from (2.32) as a function of |c ? c?| for W = 30 time slots. R-STSLmethod.2.7 SummaryIn this chapter, we considered UWL in the practical scenario where nodes are nottime-synchronized and permanently moving, and where the propagation speed isunknown. We introduced a localization algorithm which uses existing self-estimationsof motion vectors of nodes, assumed to be accurate for short periods of time. Thealgorithm utilizes the constant movements of nodes in the channel and relies onpacket exchange to acquire multiple ToA measurements at different locations. We alsopresented a method to self-evaluate the localization accuracy of the node. In addition,we used the applicable Crame?r-Rao lower bounds as references for the performanceof STSL. Considering the problem of establishing a faithful simulation environmentfor the UAC, alongside simulations we tested our algorithm in a designated sea trial.Both simulations and sea trial results demonstrated that our algorithm can cope withtime-synchronization and propagation speed uncertainties in a dynamic environment,and achieves a reasonable localization accuracy using no more than two anchor nodes.43Chapter 2. UWL with Time-Synchronization and Propagation Speed Uncertainties15 20 25 30 35 40 45 50 55 600.30.40.50.60.70.80.91x [m]P(?? err?x) W=10 (STSL)W=20 (STSL)W=30 (STSL)W=30 (MR?STSL)Figure 2.10: Probability that ??err ? x for STSL and MR-STSL.44Chapter 3Spatially Dependent UnderwaterNavigationSince nodes constantly move in the UAC, UWL is only the first step towards under-water navigation. At the presence of motion, underwater navigation must include atracking scheme that uses UWL as an initial estimation and recursively updates thelocation of a TN. A key challenge in UT is motion irregularities, which makes it hardto determine the SSM. In this context, since the ocean current is usually spatiallycorrelated [100], it is reasonable to assume spatial dependencies between the drift mo-tions of network nodes in close proximity. In [101], it was argued that birds in a flockor school of fish obey simple rules for distance and speed relative to others in closeproximity. With this idea in mind, [102] assumes temporal and spatial dependenciesbetween the motion of nodes in an indoor environment and achieves good trackingperformance. For UWL, [103] exploited spatial dependencies for the collaborativelocalization of fleets of vertically sinking drifters assuming equal speed. For track-ing, an acoustic Doppler current profiler is used in [104] to measure ocean currentsat different depths and update the SSM accordingly. Ocean current is treated as aseparate state parameter in [105], and both the self-propelled and drift speed of theAUV are estimated using an EKF. However, while spatial dependencies between thedrifts of the anchor nodes is observed from sea trial results in [105], to the best ofour knowledge it has not been incorporated in a UT scheme yet.In this chapter, we propose a UT scheme that accounts for the effect of oceancurrent, considers sound speed uncertainties, and incorporates Doppler shift mea-surements. To the best of our knowledge, neither of these three components hasbeen considered before for UT. We exploit that the ocean current is spatially corre-lated, and causing correlated drift velocities of nodes participating in the tracking.In particular, by letting anchor nodes report their drift velocities through acousticcommunication, we estimate the drift velocity of the TN as a combination of theformer. We therefore refer to our proposed tracking solution as the drift dependentUT (DD-UT) schemeWe offer two SSM-based tracking solutions, which are based on the EKF andthe unscented Kalman filter (UKF), respectively. The EKF is a modification tothe Kalman filter which linearizes the state-space and measurement model using thecurrent predicted state. While the EKFs are extensively used in both GPS andunderwater navigation [34, 106], they require knowledge of the probability densityfunction of both the model and measurement noise and thus are sensitive to model45Chapter 3. Spatially Dependent Underwater Navigationmismatch [107, 108]. Instead, the UKF approximates the probability density functionby a deterministic sampling of points. If a large amount of data is available, the UKFtends to be more robust than the EKF in its estimation of error and has been provensuperior to the EKF for complex cases such as time series modeling and neuralnetwork training [107]. However, no such comparison has been made for the case ofUT.Our tracking scheme starts from initial estimates of the sound speed, the locationof the TN and its speed. Then, our SSM fuses INS measurements, drift-velocityinformation from anchor nodes, and ranging and Doppler shift estimates to anchornodes, to provide timely estimates of the 2-D location of the TN. Considering thepossibility of weak correlations between drift velocities (for example in the presence ofturbulence), we present two types of confidence indices (CIs). The first one is based onthe distance between TN and anchor node and the corresponding measured Dopplershift, and the second CI is based on the normalized variance of the anchors? velocities.Since these CIs do not depend on the estimated location of the TN, we argue thatthey are unbiased. To evaluate the performance of our DD-UT scheme, we developa hybrid simulator combining the shallow water hydrodynamic finite element model(SHYFEM) for ocean current [109] and the Bellhop ray-tracing numerical model forpower attenuation of sound in water [2]. For a set of bathymetry maps, our simulationprovides time-varying trajectories of drifting nodes along with power attenuation forall communication links. We compare the performance of our DD-UT scheme tobenchmark solutions, as well as to the recursive Crame?r-Rao lower bound applicablefor tracking. To further verify our simulations, we also report results from two seatrials for different bathymetric channel structures, conducted in the MediterraneanSea and in the Indian Ocean.The remainder of this chapter is organized as follows. In Section 3.1, we introduceour system model, which is followed by the discussion of the state-space and measure-ment model in Section 3.2. Our DD-UT scheme is introduced in Section 3.3. Next,results from both simulations and sea trials are presented in Section 3.4. Finally,conclusions are drawn in Section 3.5.3.1 System ModelOur system includes several anchor nodes at known locations and a TN, equippedwith a depth sensor, an INS, and an acoustic modem. The TN is assumed movingwith random acceleration in both the surge and angular direction. For simplicity, weneglect the pitch and roll Euler angles and assume that the vehicle is aligned withthe horizontal plane or its measurements from on-board sensors can be projected tothis plane. Furthermore, instead of the depth-dependent sound speed, c, for eachcommunication link between the anchor node and the TN, we use the average soundspeed c? (see Section 1.1.1). We assume the anchor nodes remain at the same depth,46Chapter 3. Spatially Dependent Underwater Navigationbut the TN can rise or dive in water. Hence, at time instance tk, where k is thesampled time index of the INS system, c?k is time-varying but similar for all anchornodes. To estimate c?k we adopt the widely used model for sound speed [5],ck = 1449+4.6Tk?0.055T 2k +0.0003T 3k +(1.39?0.012Tk)(Sk?35)+0.017zk , (3.1)where Tk is the temperature in degree Celsius, Sk is the salinity in parts-per-thousand,and zk is the depth in meters, and assume that c?k is the mean of ck from (3.1) overthe water column from the transmitting anchor to the TN. We assume a prior UWLprocess (e.g., [110, 14]) that time-synchronizes the TN relative to the anchor nodes,and gives an initial estimation for c?0 and for the TN location and heading.We focus on 2-D location tracking. This is accomplished by the on-board INSgiving timely estimates of the speed of the TN in both the surge and angular direction[34], and by packets transmitted by anchor nodes providing anchor-location and drift-velocity, ToF and Doppler shift information. The INS data rate (usually around100 Hz) is assumed to be much faster than that of the packet exchange. Hence, wedefine intervals ? INS and ? range, representing the time elapsed between two consecutiveacceleration and ranging measurements, respectively, and for simplicity assume ? range? INSis an integer.We define a reference grid system [x, y, z], where x and y are UTM coordinates,and z is the depth in meters. The communication packets carry the 3-D coordinatesof the anchor nodes in the reference grid system, ranck = [xanck , yanck , zanck ]T , and theanchor?s estimated drift vector,vanc,driftk =[vx,anc,driftk , vy,anc,driftk , v?,anc,driftk , ?anck]T,whose elements are the speed in the x, y, and angular directions, and the anchor?sheading direction, respectively. We assume a slowly changing ocean current velocityfield [100], such that vanc,driftk and the drift vector of the TN are correlated. Vec-tor vanc,driftk is measured at the anchor nodes by subtracting the self-propelled mo-tion, vthrustk , from the calculated location-based one. Following [39] (and referencestherein), the thrust velocity can be obtained by measuring the thrust force F force,and solving the differential equationmv?thrustk = F force ? 0.5Cdrag?A(vthrustk )2 , (3.2)where m is the vehicle mass, Cdrag is the vehicle drag coefficient, ? is the density ofthe water, and A is the vehicle cross-section area. Then, vanc,driftk can be estimatedusing a simple IIR filter, whose output is ?vanc,driftk + (1? ?)vanc,driftk?1 with parameter? tuned to the expected rate of change of the drift velocity.3.2 The SSM and Measurement ModelIn this section, we describe our SSM and measurement model. We consider motionat fixed speed and allow for acceleration noise. Model mismatches are considered by47Chapter 3. Spatially Dependent Underwater Navigationincluding drift velocity estimates of anchor nodes.3.2.1 State Space Model (SSM)Let rk = [xk, yk, zk]T and ?k be the 3-D coordinates of the TN in the reference gridsystem and its heading angle, respectively. Furthermore, let vxk , vyk and v?k be thevelocity of the TN in the reference-grid x and y axes and in the angular direction,respectively. Since the TN is assumed to be moving with random acceleration in boththe surge and angular direction, and since the average sound speed, c?k, is assumedtime-varying, we choose the SSV asak =[xk, vxk , yk, vyk , ?k, v?k , c?k]T. (3.3)To set up the SSM we assume a Gaussian distributed acceleration and mismatchfor model (3.1), respectively. Recall that the sound speed in (3.1) depends on watersalinity, temperature, and depth. For small depth changes of a few hundreds ofmeters, we can neglect changes in salinity. Furthermore, up to a water depth of 100 m,Tk ? 3100zk [5]. Thus, neglecting the second and third order terms of the temperaturedependence in (3.1), we can linearly update the average sound speed, c?k, based onthe depth change ?zk = (zk ? zk?1) + (zanck ? zanck?1) at a rate of 0.017 + 4.6 ? 3100 . Theassumed SSM is thereforeak = Bak?1 + u?zk +Nnak , (3.4)where u = [0, 0, 0, 0, 0, 0, 0.155]T , nk =[nxk, nyk, n?k , nc?k]Tis a zero-mean Gaussianvector with covariance matrix Rmodel, and the advance and noise matrices areB =??????????1 ? INS1 0 0 0 0 00 1 0 0 0 0 00 0 1 ? INS1 0 0 00 0 0 1 0 0 00 0 0 0 1 ? INS1 00 0 0 0 0 1 00 0 0 0 0 0 1??????????, N =?????????????(? INS1 )22 0 0 0? INS1 0 0 00 (?INS1 )22 0 00 ? INS1 0 00 0 (?INS1 )22 00 0 ? INS1 00 0 0 1?????????????,respectively.3.2.2 Measurement ModelWhile AUVs can be equipped with a large number of on-board sensors, for longterm missions, strict energy constraints may not permit the use of energy consumingsensors such as DVLs (e.g., [34]). We therefore consider UT using only an INSand occasional packet exchanges with anchor nodes. In the following, we state ourmeasurement model.48Chapter 3. Spatially Dependent Underwater NavigationINSThe on-board INS includes a 3-D accelerometer and angular sensors to measure theTN speed in the surge and angular direction [34],mINSk =??vx?kvy?kv?k??+ nINSk , (3.5)where vx?k and vy?k are the speed of the TN in the x and y directions of the TN localcoordinate system, such that vxk = vx?k cos?k, and vyk = vy?k sin?k. Vector nINSk in (3.5)is a 3-D measurement error vector with zero-mean Gaussian elements and covariancematrix RINS.ToFAt a certain time sample k, the TN receives a packet from an anchor node in itscommunication range and estimates the ToF to the transmitting anchor. We assumethat both the TN and the anchor move slowly relative to the propagation delay inthe channel. Thus, the ToF is modeled bymToFk =?dk?c?k+ nToFk , (3.6)where dk = rk ? ranck and nToFk is the ToF estimation error. While noise nToFk can bebiased due to NLOS false identification, we rely on our method from Chapter 4 toclassify NLOS and LOS ToF measurements, and assume nToFk is a zero-mean Gaussianwith variance ?2ToF.Doppler ShiftAlong with estimating the ToF from received packets, communication signals can beused to evaluate the Doppler shift using, e.g., [111]. The Doppler shift is determinedby the velocity difference of the anchor and the TN. Let vanck = [vx,anck , vy,anck ]T bethe velocity3 of the anchor node in the reference-grid x- and y-axis whose packet isreceived at time instance tk. Velocity vanck can be either reported by the anchor, or,assuming slow changes relative to ? range, be calculated by the TN using the anchor?sprevious reported location. For slow motion, the Doppler shift can be approximatedin terms of the frequency offset?fk =fcc?k?vrelk ? cos ?k , (3.7)3Note that the elements of vanck may be equal to the first two elements of vanc,driftk only if theanchor is not self-propelled.49Chapter 3. Spatially Dependent Underwater Navigationwhere fc is the carrier frequency, vrelk = [vxk , vyk ]T ? vanck , and ?k is the angle betweenvectors vrelk and dk (see (3.6)), such thatcos ?k =dTk ? vrelk?vrelk ??dk?. (3.8)By (3.7) and (3.8),?fkfc= (vxk ? vx,anck ) (xk ? xanck ) + (vyk ? vy,anck ) (yk ? yanck )?dk?c?k. (3.9)We model the Doppler shift measurement bymDopplerk =?fkfc+ nDopplerk , (3.10)where nDopplerk is assumed zero-mean Gaussian noise with variance ?2Doppler.VelocityThe spatial dependencies between the drift motions of the TN and the anchor nodesallow us to estimate the drift velocity of the TN by superimposing the drift velocitiesof the anchor nodes. In particular, by letting anchor nodes report their drift velocitiesand heading directions, vanc,driftk , we obtain an estimate for the TN drift velocity andheading directionv?driftk =?????v?x,driftkv?y,driftkv??,driftk??driftk?????. (3.11)This process will be discussed in detail in Section 3.3.We treat v?driftk as a measurement, modelled asv?driftk = ?dk????vxkvykv?k?k????+ ndriftk , (3.12)where ndriftk is an error vector whose elements are zero-mean Gaussian with covariancematrix Rdrift, and ?dk is a CI whose purpose is to limit the use of anchor velocitieswhen spatial dependencies between the motions of the TN and the anchor nodes aredeemed to be weak (a method to determine ?dk will be presented in Section 3.3). Wealso assume that the TN is drifting or that, if self-propelled, it can compensate forits self-propelled motion using (3.2).50Chapter 3. Spatially Dependent Underwater NavigationMeasurement VectorUsing measurements mINSk from (3.5), mToFk from (3.6), mDopplerk from (3.10), andv?driftk from (3.12), we form the measurement vectoryk =????mINSk?pkmToFk?pkmDopplerk?pk v?driftk????= h(ak) + nmeasurek , (3.13)where ?pk = 1 if a packet is received at time instance tk and 0 otherwise, h(ak) is themeasurement model vector, and nmeasurek is a zero-mean Gaussian measurement noisewith covariance matrix Rmeasure. For the SSV ak in (3.3), the SSM (3.4), and themeasurement vector yk from (3.13), we next derive the CRLB for UT.3.2.3 Crame?r Rao Lower Bound (CRLB)Let A?k = {a?1, . . . , a?k} and Y?k = {y1, . . . ,yk} be the set of SSV estimates (3.3)and measurement vectors (3.13) obtained till time instance tk, respectively. For anyunbiased estimator, the CRLB gives the lower bound on the varianceE[(A?k ?Ak)(A?k ?Ak)T]? J?1(k) , (3.14)whereJ(k) = E[? ?2?2AklogP (Ak,Yk)](3.15)is the inverse of the Fisher information matrix with elements Ji,j(k), P (?) denotesthe probability density function, and E[?] denotes expectation. In [112], it was shownthat if only estimation of ak is of interest, (3.15) can be formulated recursively suchthatJ(k) = J1,k ? JT2,k(J(k ? 1) + J3,k)?1J2,k , (3.16)whereJ1,k = ?E[?2?2aklogP (ak|ak?1)]? E[?2?2aklogP (yk|ak)](3.17a)J2,k = ?E[?2?ak?1?aklogP (ak|ak?1)](3.17b)J3,k = ?E[?2?2ak?1logP (ak|ak?1)]. (3.17c)Since both nak from (3.4) and nmeasurek from (3.13) are modeled to be zero-meanGaussians with corresponding covariance matrices Rmodel and Rmeasure, respectively,51Chapter 3. Spatially Dependent Underwater Navigationand since ?ak?ak?1 = B and introducing Hk =?h(ak)?ak, (3.17) becomesJ1,k =(NRmodelNT)?1 + E[HTk (Rmeasure)?1Hk](3.18a)J2,k = BT(NRmodelNT)?1 (3.18b)J3,k = BT(NRmodelNT)?1B . (3.18c)3.3 The DD-UT SchemeIn this section, we describe our DD-UT scheme. We start with the tracking algorithmand then proceed with the process of estimating the drift velocity of the TN.3.3.1 TrackingOur tracking scheme is based on the SSM (3.4) with measurement vector yk andmeasurement model (3.13). For the purpose of tracking, we consider the use of theKF which is an optimal solution for Gaussian model and measurement noise. Sinceour measurement model is not linear, we adopt the EKF and the UKF. The EKF isformalized byak|k?1 = Bak?1|k?1 + u?zk?1 (3.19a)P k|k?1 = BP k?1|k?1BT +NRmodelNT (3.19b)Kk = P k|k?1HTk(HkP k|k?1HTk +Rmeasure)?1 (3.19c)en = yk ? h(ak|k?1) (3.19d)ak|k = ak|k?1 +Kkek (3.19e)P k|k = (I ?KkHk)P k|k?1 , (3.19f)where I is the identity matrix, P 0|0 is taken as Rmodel, and a0|0 is obtained from aprior UWL process (see Chapter 2). The UKF requires more regression steps, andfor brevity we refer the reader to [108].We initiate Rmodel and Rmeasure as diagonal matrices. For the former, sinceslow acceleration is expected and assuming no bathymetric layer exists (i.e., model(3.1) holds), we use small variance elements relative to the initial SSV (3.3) givenby the UWL process. Similarly, since both ToF and Doppler shift measurementsare expected to be fairly accurate, small values are used for ?2ToF from (3.6) and?2Doppler from (3.10) relative to the initial estimates of range and velocity difference toanchor nodes, respectively. To determine RINS we use the INS specifications, and todetermine ?2drift, we use a self-estimation of the error in the drift velocity estimationintroduced next.52Chapter 3. Spatially Dependent Underwater Navigation3.3.2 Drift Velocity EstimationIn this section we describe the process of obtaining estimates v?driftk as a combina-tion of vanc,driftk measurements, which is then integrated into the measurement vector(3.13). Our drift velocity estimator requires ranging measurements ?dk? = mToFk c?k?1,assumed to be correct. However, we require only a coarse range estimation, mainlyto weight drift velocities from different anchor nodes.Considering the possibility of weakly correlated motion of nodes, we restrict theuse of drift velocities reported from anchor nodes. First, accounting for time changesin the ocean current, temporal restriction is performed by defining a time window,Twin (say 60 seconds), and considering only velocities of anchor nodes whose packetswere received in the last Twin sec. Second, spatial restriction is applied by using aCI, ?rk, indicating whether to use a velocity estimate received from a certain anchornode at time instance tk. This is formalized through the setKk = {k? | tk ? Twin ? tk? ? tk ? ?1,k? = 1} , (3.20)which contains the time instances for which drift velocities reported from anchornodes can be used. Furthermore, we use a second CI, ?dk , which evaluates the accuracyof estimation v?driftk and is incorporated into (3.12).Next, we present three schemes for fusing vanc,driftk measurements. Then, we for-malize the CIs ?rk and ?dk .Nearest NeighborAssuming that the spatial correlation of ocean current decreases with range, thesimplest method to obtain v?driftk is by choosing it to be the velocity of the anchornearest to the TN. Accordingly, we identifyk? = argmink??Kk?dk?? , (3.21)and setv?driftk = vanc,driftk? . (3.22)The main drawback of this method is that not all drift velocity information availablefrom anchor nodes is used. Thus, it is sensitive to noise in estimations vanc,driftk , andto irregularities of the ocean current.Weighted SuperpositionIn the weighted superposition method, velocities of anchor nodes are combined suchthatv?driftk =?k??Kk ?k?vanc,driftk??k??Kk ?k?, (3.23)53Chapter 3. Spatially Dependent Underwater Navigationwhere ?k? is a weight function.Since spatial correlation of ocean current depends on both range and depth (cf.[100]), a pragmatic choice is a normalized weight function?k? =???dk??maxk?Kk(?dk?)????, (3.24)where ? is a pre-determined exponent. However, since ocean current irregularities(e.g., turbulence) occur at small regions [100], anchor nodes located in close prox-imity to each other should be given smaller weights. Thus, the spatial density ofanchor nodes should also have an impact. Following [113], we suggest the alternativeweighting function??k?,k =???dk??maxk?Kk(?dk?)????(1 + hk?,k) , (3.25)which takes into account the locations of the TN and the anchor nodes. In (3.25),hk?,k =1?j?Kk,j 6=k? ?dj??j?Kk,j 6=k?1?dj??(1?(ranck? ? rk)T (rancj ? rk)?dk???dj?). (3.26)The exponent ? in (3.24) and (3.25) is adapted to the expected noise level of vanc,driftk .Our simulation results showed that when noise in estimates vanc,drift is high, goodperformance is obtained for ? = 0, whereas a value ? = 4 should be used for relativelyaccurate vanc,driftk estimations.Least SquaresAssuming that ocean current changes linearly in space, estimating the drift velocityof the TN can also be interpreted as estimating two planes, one each for the x-and y-coordinates of the velocity. Enumerating the elements in Kk from (3.20) ask1, k2, . . . , kLk where Lk = |Kk|, we formulate the problem as the linear combination?????xanck1 yanck1 ?anck1 1xanck2 yanck2 ?anck2 1... ... ...xanckLk yanckLk?anckLk 1?????? ?? ?P????ax ay a?bx by b?cx cy c?dx dy d?????? ?? ?C=?????vx,anc,driftk1 vy,anc,driftk1 v?,anc,driftk1vx,anc,driftk2 vy,anc,driftk2 v?,anc,driftk2... ...vx,anc,driftkLk vy,anc,driftkLkv?,anc,driftkLk?????? ?? ?V,(3.27)where the coefficients of C describe the two planes. Solving for C yieldsC =(P TWP)?1P TWV , (3.28)54Chapter 3. Spatially Dependent Underwater Navigationwhere W is a diagonal weighting matrix whose k th coefficient is ?k from (3.25). Weuse a coarse estimate of the location of the TN to calculate the velocity, such thatv?driftk =(xk?1, yk?1, ?k?1, 1)?C , (3.29)and xk?1, yk?1, and ?k?1 are the previous course estimation of the location of the TN.3.3.3 Confidence Index (CI)Utilizing the spatial dependencies between the motions of the TN and anchor nodeshas a great potential for improving UT. However, since in some cases the spatialcorrelation of ocean current may be weak, there is a need to limit the use of v?driftkthrough a CI. Considering the likely large errors in the TN velocity estimation ifspatial correlation is weak, we use two types of CIs. The first, ?rk, identifies anchornodes with whom spatial dependency of motion is weak, and is incorporated in theestimation methods using the set Kk from (3.20). The second, ?dk , is a measure forthe homogeneity of the current velocity field around the TN.CI for Correlation with Drift Velocity of a Single AnchorIn the above three schemes to evaluate v?driftk , only velocities of anchor nodes whoselast packet was received at time instance tk, k ? Kk, are considered for estimatingv?driftk (see (3.21), (3.23), and (3.27)). Assuming that spatial correlation of oceancurrent degrades with range, relative depth, and relative velocity of the anchor andTN, ?rk is a function of the distance of the TN to the anchor node, the differencebetween the distances of the TN and the anchor to the sea bottom, ?zbottom, and themeasured Doppler shift, ?fk from (3.10). To formalize this, we use three thresholdvalues Thr, Thz, and ThD, and set?rk =[1 if |dk| ? Thr ? |?zbottom| ? Thz ? ?fk c?k?1fc cos ?k?1 ? Thd0 otherwise]. (3.30)We suggest determining thresholds Thd and Thz by training on a set of trajec-tories obtained from real data, or, if not available, simulated by a numerical oceancurrent model (e.g., SHYFEM [109]) for the expected bathymetry of the environ-ment. However, a more educated guess can be made for Thr. In case the bathymtryis known, we can calculate the range for which the velocity of the water current canbe considered dependent. Otherwise, we use the Rossby radius of deformation, Dr,which is the range from a certain location at which the gravitational tide and low-related ocean current considerably change [100]. Since the effect of the gravitationalforce on ocean current is a dominant factor in determining the spatial correlationbetween the motion of the TN and the anchor, we set Thr = Dr. The Rossby radius55Chapter 3. Spatially Dependent Underwater Navigationof deformation is a function of the water depth, p, earth gravity, g, and the Coriolisforce parameter, f0, at the current latitude. For shallow water it is calculated as [100]Dr =?gpf0. (3.31)Since neighter Dr from (3.31) nor Thd and Thz are directly related to the velocity ofthe TN, and so are Thd and Thz, and since only a coarse estimation of the distancedk is needed, we consider ?rk from (3.30) to be an unbiased CI.CI for Detecting Uncorrelated MotionOur second type of CI enables self-evaluation of the estimation v?driftk . Let vk =[vxk , vyk ]T be the true velocity of the TN, and v?driftk (1? 2) be the first two elements ofvector v?driftk . To evaluate the validity of the assumption of spatially correlated oceancurrent for both speed and direction, we use the p-relative distance,ev(k) =?vk ? v?driftk (1? 2)?p?vk?p + ?v?driftk (1? 2)?p(3.32)with p = 2, as the figure of merit. Error (3.32) is a normalized measure combiningerrors in both speed and direction, whose value is 1 if the vectors point in oppositedirections and 0 if the vectors are the same.Our aim is to obtain a CI that detects large errors (3.32). To this end, we useprevious location estimates and anchors? drift velocities and find an error predictionfor a newly drift velocity information vanc,driftLk . Since a large variability of the driftvelocity field around the TN reflects poor spatial dependency of the ocean current,for such error prediction we use the variance of anchor nodes drift velocities. Theweighted variance of velocities is calculated in two steps. First, the weighted meanvelocity vector is calculated by?k =1LkLk?i=1vanc,driftki ?ki , (3.33)where ?k is the weight function from (3.24) or (3.25). Then, the weighted variancevector is computed by?2k =1LkLk?1?i=1(vanc,driftki ?ki ? ?k)2, (3.34)where the square is performed element-wise. The normalized variance is then??2k =??2k???2k?. (3.35)56Chapter 3. Spatially Dependent Underwater NavigationSince spatial dependency of ocean current should be fixed or change slowly overtime [100], we combine previous variances ??2k? , k? ? Kk, in vectors?k? =[??2k? ? ?dk?, ?dk?,??dk?, ??2k? ,???2k? , 1]T, (3.36)where ?dk? is used to match the contribution of the different velocities to the CIwith their contribution to v?driftk . Using vectors ?k? , we form a matrix of reliabilityindicators, ?k =[?k1 , . . . ,?kLk?1]T. Then, introducing ev = [ev(k1), . . . , ev(kLk?1)]Twe calculate?k =(??Tk)?1?kev , (3.37)and form the error predictione?v(Lk) = ?TLk?k . (3.38)Prediction e?v(Lk) is compared with a threshold The to form the CI?dk =[1 e?v(Lk) ? The0 otherwise]. (3.39)Similarly to thresholds Thd and Thz from (3.30), we suggest calculating The bytraining on a set of modeled trajectories.3.4 ResultsIn this section, we show and discuss the performance of our DD-UT scheme in sim-ulations and two sea trials. We compare the following velocity estimation methods:1) nearest neighbor (NN ) from Section 3.3.2, 2) weighted superposition (WSP) fromSection 3.3.2 with weighting function (3.24), 3) weighted superposition with spatialdirection consideration (WSPS ) from Section 3.3.2 with weighting function (3.25), 4)least squares (LS ) from Section 3.3.2 by replacing W from (3.28) with the identitymatrix, 5) and weighted least square (WLS ) from Section 3.3.2. We also compareperformance of the EKF and the UKF when only INS and ToF measurements are inuse and when the sound speed is considered fixed, e.g., [39]. The simulation resultsare also compared with the CRLB (3.18). We consider both evk from (3.32) and thedistance of the estimated location (x?, y?) to the true location,edk =?E[(x?k ? xk)2 + (y?k ? yk)2], (3.40)for tracking performance. Note that edk from (3.40) can readily be compared with thesquare root of J1,1(k) + J3,3(k) from (3.15).57Chapter 3. Spatially Dependent Underwater Navigation(a)0.6 0.70.8 0.91?0.200.20.40.6020040060080010001200 Speed [m/s] (x?axis)Speed [m/s] (y?axis) Time [sec]TNAnchor 1Anchor 2Anchor 3Anchor 4(b)Figure 3.1: Simulated scenario: (a) simulated velocity field at depth 20 m, (b) simu-lated speeds of the anchor nodes and TN.3.4.1 SimulationsOur simulation environment includes a single TN and L anchor nodes. To simulatethe drift motion of nodes as well as the communication link between the anchornodes and the TN, we have developed a MATLAB-based hybrid simulator combiningnumerical models for both ocean current and power attenuation in water. In thiswork we are mostly concerned with tracking in regions of shallow water and in timesteps of a few seconds. For this reason, we have chosen the SHYFEM [109], whichis designed to resolve the hydrodynamic equations for coastal seas and other smallbasins. To model power attenuation in communication links, we use the Bellhopray-tracing numerical model (cf. [2]). The simulator requires a bathymetry map,latitude information, SSP, and initial location of nodes, and produces time-varyingtrajectories of nodes based on a generated current velocity field, as well as time-varying transmission loss values for the channel between the different nodes. Forsimplicity, we assume nodes are not self-propelled4. Example for the current velocityfield and the resulting speed of four network nodes are shown in Figures 3.1a and3.1b, respectively.The structure of the simulated channel is a square area of 4 ? 4 km2 with waterdepth of 100 m, and the bathymetry map includes one or two underwater hills ofuniformly randomly generated width and depth of 3000 m and 2000 m, respectively.Initially, we place nodes uniformly at random in a smaller region of 1 ? 1 km2 onthe water surface and let them drift for 20 min according to the simulated currentvelocity field, where the TN dives at a speed of 0.05 m/sec. We use two sound speed4For a more complex motion pattern, equation (3.2) can be used.58Chapter 3. Spatially Dependent Underwater Navigation(a)0 0.2 0.4 0.6 0.8 100.10.20.30.40.50.60.70.80.91PFP D ?n=0.05 m/sec?n=0.01 m/sec?n=0.2 m/secTh2=0.1Th2=0.2Th2=0.3(b)Figure 3.2: (a) Mean of ev(k) from (3.32) as a function of number of anchor nodes,(b) PD vs. PF for different ?2anc and The values.profiles: one that matches model (3.4) and has a value of 1512 m/sec on the watersurface, and another which has a fixed value of 1512 m/sec. We consider a carrierfrequency of 15 kHz, and the model is set for location 49o16?13.33??N , 126o16?6.4??W(i.e., near Vancouver, BC). We let the anchor nodes transmit in a simple TDMAnetwork, with a time frame of L time slots of duration 10 sec. In its designated timeslot, each anchor transmits a data packet including its 3-D coordinates and its driftvelocity. We assume packets are received error-free as long as the signal-to-noise ratio(SNR) is above 15 dB. The SNR is calculated based on the output transmission lossfrom our simulator for a common power source level of 140 dB//?Pa@1m, bandwidthof 1 kHz, symbol duration of 1 msec, and an ambient noise level of 50 dB//?Pa/Hz.We first investigate the performance of our velocity estimate v?driftk from (3.12) andthe CI ?dk from (3.39). Tracking performance follows next.Velocity Estimate (simulations)In Figure 3.2a, we show the effect of the number of anchor nodes on the velocityestimation. Results are averaged over 81 scenarios of different bathymetry maps and1000 simulation runs for each scenario. For all methods, performance improves withthe number of anchor nodes. As expected, for large L, the performance of the WLSis better than that of the LS method. This improvement confirms our assumptionof range related spatial correlation of the ocean current. However, when L ? 5 noimprovement is observed. One reason for this is the sensitivity of the WLS methodto scenarios where a cluster of anchor nodes is formed and the velocity plane can betilted further due to noises at the measured anchors? locations. The range relatedspatial correlation is also the reason for the better performance of the WSP and WSPSmethods compared to that of the NN. However, the spatial direction considered in59Chapter 3. Spatially Dependent Underwater Navigationthe weighting function of the WSPS method does not seem to have a large impacton performance. We observe that for L < 8, the performances of the NN, WSP, andWSPS methods are better than those of the WLS or LS. This is because of the largeramount of information required for the latter to estimate matrix C from (3.27). Inthe following, we show results only for the WSP method.In Figure 3.2b, we show the detection probability, PD, vs. the false alarm prob-ability, PF, for the CI ?dk from (3.38). To calculate PD and PF, we set a boundaryev(k) = 0.1, such that a correct detection is declared when ev(k) ? 0.1 and ?v ? The,or when ev(k) ? 0.1 and ?v ? The. We compare results for three threshold val-ues, The (see (3.38)). Based on the results of Figure 3.2b, in the following we useThe = 0.2.Tracking (simulations)For tracking purposes, we follow the generated velocity field to calculate the time-varying location of the TN. To test the convergence of the tracking scheme, at thebeginning of each simulation run the TN is given an initial SSV estimate, a?0 (see(3.3)), which includes zero-mean Gaussian noise whose covariance matrix is suchthat the SNR is 10 dB and the signal is regarded as the elements of the true SSVa0. We compare performances of our DD-UT tracking method, tracking without driftvelocity estimation (No Drift), tracking without drift velocity estimation and Dopplershift measurements (No Drift+Doppler), and the latter scheme but without trackingthe sound speed (No Drift+Doppler+Speed). We measure performance using theEuclidian distance edk from (3.40). Since in our system measurement noise componentsin yk from (3.13) have different units, to show effect of measurement error we set theelements of Rmeasure based on a measurement SNR defined asMSNR = E [h(a)2i ](?measurei )2 , (3.41)where h(a)i and ?measurei are the ith element of h(a) from (3.13) and the diagonal ofRmeasure, respectively. We note that the values considered for the MSNR are takenfrom real INS systems and the results obtained in our sea trials (they are on the orderof 50 dB, see [34] and Chapter 4).We start by showing a tracking example for the simulated scenario presented inFigures 3.1a and 3.1b for MSNR = 40 dB. In the simulated scenario, anchor nodes 4and 3 are located roughly 200 m and 700 m away from the TN, respectively, whileanchor nodes 2 and 1 are located more than 1.5 km away. As a result, we observesimilarities between the speeds of the TN and anchor nodes 4 and 3. The scenario in-cludes two sudden changes for the speed of the TN in the x direction at time instances760 sec and 990 sec. To show the effect of link communication errors, we simulatedfailures for all packets received by the TN between time instances 400 sec and 500 sec.In Figure 3.3b, we show tracking performance using the UKF as a function of time,60Chapter 3. Spatially Dependent Underwater Navigation0 200 400 600 800 1000 120015101512151415161518152015221524Time [sec]Propagation Speed [m/sec] True 1True 2UN?DPUN?DP (mismatch)No Drift+Doppler(a)0 200 400 600 800 1000 12000510152025303540Time [sec]ed k [m] DD?UT (WSP)DD?UT (WSP,mismatch)No DriftNo Drift+DopplerNo Drift+Doppler+Speed(b)Figure 3.3: Simulated scenario: (a) speed of the TN, (b) tracking error edk from (3.40)where we used the WSP drift velocity estimation method for our DD-UT scheme.When a change in the velocity of the TN occurs, an increase in edk is observed, fol-lowed by a recovery process. Since the DD-UT scheme uses more information, itconverges roughly 30 sec before the other schemes. Moreover, the convergence is fora slightly lower error edk than that of the reference methods. Comparing the effectof adding Doppler shift information, larger errors and slower recovery is recordedwhen Doppler shift measurements are not in use (No Drift+Doppler) and the veloc-ity changes. In addition, due to the change in the sound speed with depth (roughly9 m/sec over the simulation time), a deterioration of performance of up to 5 m existswhen the sound speed is considered fixed (No Drift+Doppler+Speed). As expected,errors accumulate in time when communication link failure occurs and tracking isperformed solely based on INS measurements. Also here, recovery is faster using theDD-UT scheme. To comment on the sensitivity of our protocol to a mismatchedmodel for the SSP, in Figure 3.3b we also show results (DD-UT (mismatch)) of theDD-UT scheme where the SSM is based on (3.4) but the actual SSP is fixed (whichalso affects the Bellhop-based simulated power attenuation). While slower recoveryfrom velocity changes or communication link failures is observed, performance onlyslightly decreases and results are better than without tracking the sound speed. Toemphasize this, in Figure 3.3a we show tracking performance of c?k, where curvesTrue 1 and True 2 represent the SSP used in (3.4) and the fixed one, respectively.We observe that when the SSP matches (3.4), the estimated c?k closely follows thetrue one. However, inaccuracies exists when Doppler shift measurements are not inuse (No Drift+Doppler) and the additional information extracted from (3.9) is notavailable. When the SSP is fixed and the SSM is mismatched, performance decreasesbut the estimated sound speed is still within 1 m/sec from its true value.Next, we show statistical tracking results. We generate 81 bathymetry maps61Chapter 3. Spatially Dependent Underwater Navigation0 5 10 15 20 25 3000.10.20.30.40.50.60.70.80.91xProb(ed k > x) DD?UT (UKF)DD?UT (EKF)No DriftNo Drift+DopplerNo Drift+Doppler+SpeedCRLB (average)(a)35 40 45 50 55 60 65 700246810121416Measurement SNR [dB]ed k [m] DD?UT (UKF)DD?UT (UKF,mismatch)DD?UT (EKF)No DriftNo Drift+DopplerNo Drift+Doppler+SpeedCRLB(b)Figure 3.4: Tracking performance (simulations): (a) empirical C-CDF forMSNR=50dB, (b) edk as a function of MSNR (thick line shows results for the CRLB).and for each realization we conduct 1000 simulation runs. In each simulation, wedetermine the initial location of the network nodes uniformly at random, and generatenew noise realizations. We note that since the scenarios generated include speed anddirection changes (much like the example shown in Figure 3.1b), we cannot determinea point in time for which tracking converges. Instead, considering convergence timeof up to 100 sec from the start of each simulation, we take the median of all edkresults obtained between time instances 100 sec and 1200 sec. In Figure 3.4b, weshow the effect of the MSNR. The results indicate a relatively moderate effect ofthe measurement noise on performance for the MSNR range considered. Hence,as argued above, the dominant effect in tracking is the ability to accurately modelthe motion of the TN. As expected, performance improves the more information isused. We observe that estimating the drift velocity has more benefit than utilizingDoppler shift measurements. We note that performance of the UKF are notablybetter than that of the EKF, and that this performance gain is more significant forlow MSNR. This is mainly due to the first order approximation of the EKF [108] andthe fast changes in the TN speed and direction relative to the time step interval ofthe filter which are known to cause instabilities [107]. This is emphasized by the goodperformance of the UKF even for a mismatched SSP. The need to track the soundspeed is demonstrated by the significant performance gain of the No Drift+Dopplerscheme compared to that of the No Drift+Doppler+Speed. In Figure 3.4b, we alsoshow the square root of the CRLB J1,1(k)+ J3,3(k) from (3.15). We observe that theCRLB is well approached by our DD-UT scheme. To comment on the variations of theperformance results, in Figure 3.4a we show the empirical complementary cumulativedensity function (C-CDF) of edk for MSNR= 50 dB, as well as the CRLB for reference.We observe small variations of edk for our DD-UT scheme using the UKF. The results62Chapter 3. Spatially Dependent Underwater Navigation?100 0100 2003000200040006000020406080100120 x [m]y [m] time [min]Boat 1Boat 2(a)2000 40006000 8000010002000300002000400060008000 UTM Location (x?axis)UTM Location (y?axis) Time [sec]Anchor 1Anchor 2Anchor 3TN(b)Figure 3.5: Node locations in Cartesian coordinates: (a) Singapore sea trial, (b) Israelsea trial.verify our conclusions from the average performance shown in Figure 3.4b.3.4.2 Sea TrialsDue to the complex sea environment and our reliance on models for the SSM, sensornoise, and on the physical phenomena of spatial correlation of the ocean current, oursimulation performance requires verification in a sea environment. For this reason,we have conducted two sea trials in the different bathymetry channel structures of theMediterranean Sea and the Indian Ocean. The former was performed in August 2010at Haifa, Israel, and included three vessels, representing three mobile nodes, whichdrifted for more than 2 hours. Water depth was around 40 m and the height of thesea waves was around 0.1m. The three nodes deployed transducers at depth of 10 m,and transmitted data packets every 3 minutes using a TDMA protocol includingthe nodes? 2-D UTM coordinates (which was measured at rate of 1 sec using GPSreceivers). The received packets were used to calculate the ToF (applying the methoddescribed in Chapter 4), and the Doppler shift (using the method described in [111]).No acceleration measurements were taken during this sea trial, and instead we use theGPS readings to calculate mINSk (see (3.5)). The time-varying location of the threenodes is shown in Figure 3.5b. The second sea trial was conducted in November2011 at the Singapore strait with water depths of 15 m and ocean wave height ofapproximately 0.5 m. The experiment lasted for more than one hour and includedtwo drifting boats. At each boat, we obtained 3-D acceleration measurements ata rate of f = 4.8 Hz using an Libelium Wasp Mote?s on-board accelerometer, andToF measurements were taken at each node every 20 sec using underwater acousticmodems, manufactured by Evologics GmbH, which were deployed at a depth of 5 m.Throughout the experiment, the locations of the boats were monitored using GPSreceivers at rate of 3 sec. The boats? locations are shown in Figure 3.5a.63Chapter 3. Spatially Dependent Underwater NavigationTable 3.1: Drift velocity estimation results from the Israel and Singapore sea trials.NN WSP LSSea Trial Measure N-I W-I N-I W-I N-I W-IIsrael E{|vk ? v?driftk (1? 2)|} 0.03 0.11 0.03 0.11 0.03 0.24E{|?k ? ??driftk |} 0.16 0.37 0.09 0.25 0.07 0.28Singapore E{|vk ? v?driftk (1? 2)|} 0.09 0.16 - - - -E{|?k ? ??driftk |} 0.17 0.22 - - - -For both sea trials, we use the GPS readings to calculate the drift velocity of theanchor nodes, vanc,driftk , and as the ground truth for the position and velocity of theTN.Velocity Estimate (sea trial)In Table 3.1, we compare the mean of errors |vk ? v?driftk (1 ? 2)| (in m/sec) and|?k ? ??driftk | (in radians) for the NN, WSP, and LS methods described in Section 3.3.To comment on the robustness to motion instabilities, we present results when alldata is considered (W-I ), and when disregarding data with motion irregularities (N-I ). For the Israel sea trial, motion irregularities are considered around time instances2400 sec, 5400 sec, and 7200 sec, and for the Singapore sea trial around time instances100 sec, 1500 sec, and 2900 sec. We note that all the above motion irregularities aredetected by our CIs, ?dk and ?rk, and both ?dk and ?rk are 1 in other time instances.For both sea trials, we observe good estimation results when motion is stable, andmore noisy estimations but still within acceptable range when motion irregularitiesare considered. The results imply that the WSP method outperforms both NN andLS methods.Tracking (sea trial)In this section, we show tracking performance for the two sea trials. Since we usedcommercial GPS receivers whose expected accuracy was ?5 m, we consider a trackingresult for which edk ? 10 m as accurate. In both sea trials, the depth of the transducerswas fixed. Thus, we do not track the sound speed and instead use a predefined soundspeed and modify the SSM (3.4) accordingly. In the Israel sea trial, we have measureda sound speed of c = 1550 m/sec which reduced by 2 m/sec along the water columnof 40 m. In the Singapore sea trial, we rely on yearly measurements showing a fixedsound speed of 1540 m/sec. In addition, in the latter sea trial we did not have accessto Doppler shift measurements, and vector yk from (3.13) is modified to not includemeasurement mDopplerk from (3.10). In the following, tracking starts given an initialestimate a0 (e.g., from a UWL process) such that ed0 = 10 m, and for the Israel sea64Chapter 3. Spatially Dependent Underwater Navigation0 500 1000 1500 2000 2500 3000 3500 400005101520Time [sec]Euclidian Error [m] DD?UT (UKF)DD?UT (EKF)No Drift (UKF)(a)0 1000 2000 3000 4000 5000 6000 70000510152025Time [sec]ed k [m] DD?UT (WSP)DD?UT (NN)DD?UT (LS)No DriftNo Drift+Doppler(b)Figure 3.6: edk from (3.40) vs. time: (a) Singapore sea trial, (b) Israel sea trial.trial c?0 = 1545 m/sec. We next show results for all the acquired data (i.e., withmotion irregularities) and incorporate both types of CIs introduced in Section 3.3.3.Tracking results for the Singapore sea trial are shown in Figure 3.6a. For allmethods, convergence is reached after roughly 200 sec (i.e., using 10 received packetsfrom the anchor). We observe several cases where performance degrades (e.g., aroundtime instances 500 sec, 2000 sec, and after 2800 sec). As the effect is similar to theone observed in Figure 3.3b, these performance fluctuations are explained by bothToF inaccuracies (due to the shallow water environment) and velocity changes (asseen in Figure 3.5a). Similar to the results in Figure 3.4b, better performance isobtained by the UKF compared to the EKF. As seen in the figure between timeinstances 2000 and 2500, the No Drift scheme sometimes achieves better results thanour DD-UT scheme. However, most of the time significant performance gain of upto 5 m is obtained by the DD-UT scheme compared to the No Drift one.In Figure 3.6b, we show tracking performance for the Israel sea trial. For all meth-ods, a sudden spike in tracking error occurs at around 2400 sec. This error is dueto the loss of communication in the network that occurred between time instances2400 sec and 3200 sec. The similar decrease in performance at around 4000 sec isexplained by a sudden change of velocity. Both faster recovery and absolute perfor-mance improvement are observed when drift velocity estimation is used. Moreover,performance significantly improves when Doppler shift measurements are utilized.Comparing the three methods of velocity estimation5, similar to our simulation re-sults, best performance is obtained by the weighted superposition scheme for whichthe tracking error converges to less than 5 m.5We note that the WLS method yielded almost the same results as the LS method65Chapter 3. Spatially Dependent Underwater Navigation3.5 ConclusionsIn this chapter, we considered the problem of UT, where nodes experience irregulardrifting motion. Assuming spatial correlation of the ocean current, we suggested threemethods to estimate the drift velocity of the TN by superimposing drift velocitiesreported by the anchor nodes. This augments the amount of information available forthe TN as an unbiased velocity estimate. To combat errors due to weakly correlatedmotion of nodes, we proposed two types of unbiased confidence indexes. Our schemealso incorporates sound-speed tracking and Doppler shift measurements availablefrom received packets, which, due to the uncertainty of the sound speed and thecontinuous motion of nodes in the channel, is a challenging task. We implementedtracking using the UKF and the EKF, and evaluated performance using numericalmodels for ocean current and power attenuation in the ocean. The results indicatea large gain of our scheme compared to reference approaches which do not considerdrift velocity estimation, Doppler shift measurements, or sound speed uncertainties,and showed that the CRLB is well approached by our scheme. The results fromtwo sea trials for different bathymetry channel structure conducted in Israel and inSingapore support the findings from simulations.66Chapter 4LOS and NLOS Classification forUWLIn the previous chapters, we have presented our approach for UWL and UT. Asexpected and observed in Figures 2.6 and 3.4b, localization accuracy depends onthe measurement noise from internal and external sources. In this chapter, we fo-cus on ranging noise, and specifically on the problem of mistaking NLOS-related ToFmeasurements as LOS measurements. Several works suggested ways to avoid misiden-tification of NLOS and LOS links. In [114], direct sequence spread spectrum (DSSS)signals, which have narrow auto-correlation, are transmitted to allow better separa-tion of paths in the estimated channel response. Averaging ToF measurements fromdifferent signals is suggested in [115], where results show considerable reduction inmeasurement errors. In [15], NLOS-related noise in the UAC is modeled using the Ul-tra Wideband Saleh-Valenzuela (UWB-SV) model, and multipath noise is mitigatedusing the approach introduced in [116]. A way to identify NLOS measurements is pre-sented in [95], where assuming that NLOS-based measurements have larger variancethan LOS-based measurements, measurements which increase the global variance arerejected. In [117], localization accuracy is improved by selecting ToF measurementsbased on minimal statistical mode (i.e., minimal variance and mean). Alternatively,the authors in [118] suggested a method for reducing the effect of NLOS-based noiseby assigning each measurement with a weight inversely proportional to the differencebetween the measured and expected distances from previous localization. However,no systematic way has been suggested to classify NLOS- and LOS-related ToF mea-surements. Instead, we take a more rigorous approach and estimate the distributionof the ToF measurements.The intuition behind our approach is that the continuous motion of nodes in theUAC changes causes shift in the arrival times and energy of signals received throughdifferent propagation paths, causing link variations. This diversity can be exploitedby obtaining multiple propagation delay (PD) measurements and classifying theminto classes of LOS, SNLOS and ONLOS (referring to Figure 1.4, recall SNLOSrefer to sea surface and bottom reflections and ONLOS to reflection off an object). Indoing so, we also take into account the effect of mobility on distance. The proposedclassification can improve the accuracy of UWL by rejecting or correcting NLOS-related PD measurements to obtain a single ToF estimation, or by using them tobound range estimation.To implement our scheme we present a two-step classification algorithm. We67Chapter 4. LOS and NLOS Classification for UWLfirst identify ONLOS-related PD-measurements by comparing PD-based range esti-mations with range estimations obtained from received-signal-strength (RSS) mea-surements. Considering the difficulty in acquiring an accurate attenuation model,our algorithm requires only a lower bound for RSS-based distance estimations. Afterexcluding PD measurements related to ONLOS, we apply a constrained expectation-maximization (EM) algorithm (cf. [91]) to further classify the remaining PD mea-surements into LOS and SNLOS. Through a clustering of PD measurements we mit-igate changes in propagation delay due to mobility of nodes. The EM algorithmalso estimates the statistical parameters of both classes, which can be used to im-prove the accuracy of UWL. Results from extensive simulations and three sea trialexperiments in different areas of the world demonstrate the efficacy of our approachthrough achieving a high detection rate for ONLOS links and good classification ofnon-ONLOS related PD measurements into LOS and SNLOS. To the best of ourknowledge, no prior work considered a machine learning approach for NLOS andLOS classification of multiple PD measurements. Moreover, link classification forranging was not investigated for the special characteristics of the UAC.It should be noted that while our approach can also be adapted to other typesof fading channels, it is particularly suited for UWL for the following reasons. First,our algorithm relies on significant power absorption due to reflection loss in ONLOSlinks, which are typical in the underwater environment. Second, we assume that thedifference in propagation delay between signals traveling through SNLOS and LOSlinks is noticeable, which is acceptable in the UAC due to the low sound speed inwater (approximately 1500 m/sec). Third, our algorithm is particularly beneficial incases where NLOS paths are often mistaken for the LOS path, which occurs in UWL,where the LOS path is frequently either not the strongest or non-existent. Last,we assume that the variance of PD measurements originating from SNLOS links isgreater than that of measurements originating from LOS, which fits channels withlong delay spread such as the UAC.The remainder of this chapter is organized as follows. Our system model and as-sumptions are introduced in Section 4.1. In Section 4.2, we present our approach toidentify ONLOS links. Next, in Section 4.3, we formalize the EM algorithm to classifynon-ONLOS related PD measurements into LOS and SNLOS. Section 4.4 includesperformance results of our two-step algorithm obtained from synthetic UAC envi-ronments (Section 4.4.1) and from three different sea trials (Section 4.4.2). Finally,conclusions are offered in Section 4.5.4.1 System Setup and AssumptionsReferring to Figure 1.4, our system comprises of one or more transmitter-receiverpairs, (u, aj), exchanging a single communication packet of N symbols or impulsesignals, from which a vector X = [x1, . . . , xN ] of PD measurements xi, and corre-68Chapter 4. LOS and NLOS Classification for UWLsponding measured time ti, is obtained using detectors such as in, e.g., [44, 46, 47].We model xi such thatxi = xLOS + ni , (4.1)where xLOS is the PD in the LOS link, and ni is the measurement noise. We note thatdue to the mobility of nodes, xLOS is not strictly constant during the measurementinterval, and this motion bounds the accuracy of ranging.We focus on transmission over short range (on the order of a few km), for whichrefraction of acoustic waves is negligible and propagation delay in the LOS link ison average shorter than in the NLOS links. We assume successively transmittedsignals for PD measurements are separated by guard intervals such that ni in (4.1)can be assumed i.i.d. For each measurement xi a PD-based estimate, dPDi , is obtainedby multiplying xi with an assumed propagation speed, c. In addition, based on anattenuation model for an LOS link, we obtain RSS-based range estimates, dRSSi , i =1, . . . , N , from the received signals.In the following, we introduce our system model for obtaining RSS-based rangemeasurements as well as the assumed PDF for PD measurements.4.1.1 RSS-Based Range MeasurementsLet dLOS denote the distance corresponding to xLOS, i.e., dLOS = xLOSc. For thepurpose of obtaining RSS-based range measurements, we use the popular model [5]TLLOS(dLOS) = PL(dLOS) + AL(dLOS) + \u000F , (4.2)where PL(dLOS) = ? log10(dLOS) is the propagation loss, AL(dLOS) = ?dLOS1000 is theabsorption loss, ? and ? are the propagation and absorption coefficients, respectively,and \u000F is the model noise assumed to be Gaussian distributed with zero mean andvariance ?. Considering the simplicity of the model in (4.2), we do not directlyestimate dRSSi but rather estimate a lower bound dRSS,LBi , for which we apply upperbounds for ? and ? in (4.2) according to the expected underwater environment.For an ONLOS link with distance dONLOS = dONLOS,1 + dONLOS,2, where dONLOS,1and dONLOS,2 are the distance from source to reflector and from reflector to receiver,respectively6, we assume that the power attenuation in logarithmic scale is given by[5]TLONLOS(dONLOS) = TLLOS(dONLOS,1) + TLLOS(dONLOS,2) + RL , (4.3)where RL is the reflection loss of the reflecting object, whose value depends on thematerial and structure of the object and the carrier frequency of the transmittedsignals. Since RL is often large (see examples in [5]), and due to the differencesbetween models (4.2) and (4.3), we further assume thatTLONLOS(dONLOS)\u001D TLLOS(dONLOS) . (4.4)6Referring to the ONLOS link between node pair (u, a2) in Figure 1.4, dONLOS,1 = d21 anddONLOS,2 = d22.69Chapter 4. LOS and NLOS Classification for UWL4.1.2 PDF for PD MeasurementsSince we assume changes in xLOS are bounded by a small transmitter-receiver motionduring the time X is obtained, we can model the PDF of the noisy measurement xias a mixture of M = 3 distributions, corresponding to LOS, SNLOS, and ONLOSlinks, such that (assuming independent measurement noise samples in (4.1))p(X|?) =?xi?XM?m=1kmp(xi|?m) , (4.5)where ? = [?1, k1, . . . ,?M , kM ], ?m are the parameters of the mth distribution, andkm (M?m=1km = 1) is the a-priori probability of the mth distribution. Clearly, p(X|?)depends on multipath and ambient noise in the UAC, as well as on the detector usedto estimate xi. While recent works used the Gaussian distribution for p(xi|?m) (cf.,[86] and [10]), since multipath and ambient noises are hard to model in the UAC,we take a more general approach and model it according to the generalized GaussianPDF [119], such thatp(xi|?m) =?m2?m?(1?m)e?(|xi??m|?m)?m(4.6)with parameters ?m = [?m, ?m, ?m]. We associate the parameter vectors ?1, ?2,and ?3 with distributions corresponding to the LOS, SNLOS, and ONLOS links,respectively. Thus, by (4.1), ?1 = xLOS. Without prior knowledge of the actualdistribution of PD in the LOS and NLOS links, the use of parameter ?m in (4.6) givesour model a desired flexibility, with ?m = 1, ?m = 2, and ?m ? ? corresponding toLaplace, Gaussian, and uniform distribution, respectively. The flexibility and fit ofmodel (4.6) is demonstrated using sea trial results in Section 4.4.2.Following [95] and [117], we assume that PD measurements of NLOS links increasethe variance of the elements of X. Thus, if ?1, ?2, and ?3 are the respective variancesof measurements related to the LOS, SNLOS, and ONLOS links, then we have?1 < ?m, m = 2, 3. (4.7)Since, for the PDF (4.6),?m = (?m)2?(3?m)?(1?m) , (4.8)and by (4.8), ?m does not change much with ?m, constraint (4.7) can be modified to?1 < ?m, m = 2, 3. (4.9)70Chapter 4. LOS and NLOS Classification for UWLFurthermore, let TLIR be the assumed length of the UAC impulse response, whichis an upper bound on the time difference between the arrivals of the last and firstpaths. Then, since ?1, ?2, and ?3 in (4.8) capture the spread of measurements relatedto the LOS, SNLOS, and ONLOS links respectively,??m < TLIR, m = 1, 2, 3 . (4.10)Moreover, the propagation delay through the LOS link is almost always shorter thanthose for any NLOS link7. Hence, we have?1 < ?m < ?1 + TLIR, m = 2, 3. (4.11)Together with assumption (4.7), assumption (4.11) manifests the expected differencesbetween the PD in the LOS and NLOS links. Clearly, the more separable PD mea-surements from LOS and NLOS links are (i.e., the propagation delay difference islarger), the better the classification will be. Since the channel impulse response islonger for deeper channels, classification accuracy is expected to improve with depth.4.1.3 Remark on Algorithm StructureWe offer a two-step approach to classify PD measurements into LOS, SNLOS, andONLOS. First, assuming large attenuation in an ONLOS link, we compare PD-basedand RSS-based range estimates to differentiate between ONLOS and non-ONLOSlinks. Then, assuming PDF (4.6) for PD measurements, we further classify non-ONLOS links into LOS and SNLOS links.The reason for separating classification of ONLOS and SNLOS links is insufficientinformation about the distribution of the two link types. For example, delay inONLOS links may be similar to or different from that of SNLOS links. In the formercase, classification should be made for M = 2 states, while for the latter threestates are required. Since a mismatch in determining the number of states maylead to improper classification, we rely on the expected high transmission loss inONLOS links to first identify these links. Furthermore, a separate identification ofONLOS links can be used as a backup to our LOS/NLOS classifier. That is, wecan still identify the link as ONLOS even when the channel is fixed, and thus PDmeasurements originate from a single link-type. In the following sections, we describeour two-step approach for classifying PD measurements.4.2 Step One: Identifying ONLOS LinksConsidering (4.4), we identify whether measurement xi ?X is ONLOS-related basedon three basic steps as follows:7We note that in some UWACs, a signal can propagate through a soft ocean bottom, in whichcase SNLOS signals may arrive before the LOS signal [5]. However, such scenarios are rare and weassume that on average relation (4.10) holds.71Chapter 4. LOS and NLOS Classification for UWL? Estimation of dPDiWe first obtain the PD-based range estimation as dPDi = c ? xi.? Estimation of dRSS,LBiNext, assuming knowledge of the transmitted power level, we measure the RSSfor the ith received signal/symbol, and estimate dRSS,LBi based on (4.2), replac-ing ? and ? with upper bounds ?UB and ?UB, respectively.? ThresholdingFinally, we compare dPDi with dRSS,LBi . If dRSS,LBi > dPDi , then xi is classified asONLOS. Otherwise, it is determined as non-ONLOS.Note that the use of the above step Thresholding relies on assumption (4.4). Toclarify this, consider an ONLOS link. The RSS-based range estimation, dRSS,LBi , isobtained from an upper bound for an attenuation model (4.2), i.e., from applying anattenuation model for an LOS link to an ONLOS link. Since the latter is expectedto have a much larger power attenuation than the used model, it follows that dRSS,LBiwould be much larger than dPDi . Similarly, consider a non-ONLOS link (i.e., LOS orSNLOS). Here, since we use an upper bound for the attenuation model, we expectdRSS,LBi to be smaller than dPDi .Next, we analyze the expected performance of the above ONLOS link identifica-tion algorithm in terms of (i) detection probability of non-ONLOS links, Prd,non?ONLOS,and (ii) detection probability of ONLOS links, Prd,ONLOS. To this end, since explicitexpression for dLOS cannot be obtained from (4.2), in the following, we use the upperbound d?RSS,LB such thatlog10(d?RSS,LB) =TL?UB. (4.12)We note that (4.12) is a tight bound when the carrier frequency is low or when thetransmission distance is small.4.2.1 Classification of non-ONLOS LinksFor non-ONLOS links, we expect dRSS,LBi ? dLOS. Thus, since by bound (4.12),Pr(dRSS,LBi ? dLOS) ? Pr(d?RSS,LBi ? dLOS), and substituting (4.2) in (4.12), we getPrd,non?ONLOS ? 1?Q((?UB ? ?) log10 (dLOS)? ?dLOS1000?), (4.13)where Q(x) is the Gaussian Q-function.72Chapter 4. LOS and NLOS Classification for UWL4.2.2 Classification of ONLOS LinksWhen the link is ONLOS, we expect dRSS,LBi ? dONLOS. Then, substituting (4.3) in(4.12), and since Pr(dRSS,LBi ? dONLOS) ? Pr(d?RSS,LBi ? dONLOS), we getPrd,ONLOS ? Q(?UB log10 (dONLOS)? ? log10 (dONLOS,1dONLOS,2)? ?dONLOS1000 ? RL?).(4.14)Next, we continue with classifying non-ONLOS links into LOS and SNLOS links.4.3 Step 2: Classifying LOS and SNLOS LinksAfter excluding ONLOS-related PD measurements in Step 1, the remaining elementsof X, organized in the pruned vector Xex, are further classified into LOS (m = 1)and SNLOS (m = 2) links and their statistical distribution parameters, ?m, areestimated. Before getting into the details of our LOS/SNLOS classifier, we firstexplain its basic idea.4.3.1 Basic IdeaThe underlying idea of our approach is to utilize the expected variation in link typeof PD measurements due to mobility of nodes at sea. After identifying ONLOSlinks, this variation means that our set includes PD measurements of different valuesand link types. This allows us to use a machine learning approach to classify themeasurements into two classes, LOS and NLOS. For this purpose, we use the EMalgorithm. While using EM to classify measurement samples into distinct distribu-tions is a common approach, here the distribution parameters in (4.6) should alsosatisfy constraints (4.10), (4.9), and (4.11), where the two latter constraints introducedependencies between the parameters of the LOS and NLOS classes. Furthermore,we incorporate equivalence constraints to group measurements of similar values intoclusters which elements are classed to the same link type, thereby mitigating shiftsin the value of xLOS due to nodes? mobility. As we show later on, this result in anon-convex maximization of the log-likelihood function. For this reason, we presenta heuristic sub-optimal algorithm.In the following, we start by formalizing the equivalence constraints, and formu-lating the log-likelihood function. Next, we formulate a constrained optimizationproblem to estimate the distribution parameters, and present our heuristic approachto efficiently solve it. Given this estimate, we calculate the posterior probability ofeach PD measurement belonging to the LOS and SNLOS class, and classify the ele-ments of Xex accordingly. Finally, we describe how the initial solution, required forthe EM algorithm, is obtained and calculate the hybrid Crame?r-Rao bound to boundthe performance of our classifier.73Chapter 4. LOS and NLOS Classification for UWL4.3.2 Equivalence ConstraintsIn setting equivalence constraints, we assume that the identity and delay of the dom-inant channel path, used for PD detection, is constant within a given coherence time,Tc, and that for a bandwidth B of the transmitted signal, system resolution is lim-ited by ?T = 1B . PD measurements satisfying equivalence constraints are collectedinto vectors ?l, l = 1, . . . , L, where L denotes the number of such equivalence sets.Each PD measurement is assigned to exactly one vector, i.e., ?l have distinct ele-ments. To formalize this, we determine xi (recall that measurement xi correspondsto measurement time ti) and xj being equivalent, denoted as xi ? xj, if|ti ? tj| ? Tc (4.15a)|xi ? xj| ? ?T . (4.15b)To illustrate this, let xi, xj, and xn correspond to the same class (either LOS orNLOS), such that xi ? xj and xj ? xn. Then, although due to motion of nodes,xi and xn may not satisfy condition (4.15b), since ?l, l = 1, . . . , L, are distinctvectors, all three measurements xi, xj, and xn are grouped into same vector ?l andare further classified to the same state. That is, vectors ?p and ?q are merged ifthey have a common element. To form vectors ?l, l = 1, . . . , L, we begin with |Xex|(|?| symbolizes the number of elements in vector ?) initial vectors of single PDmeasurements, and iteratively merge vectors. This process continues until no twovectors can be merged. As a result, we reduce the problem of classifying xi ? Xexinto classifying ?l, which account for resolution limitations and node drifting.4.3.3 Formalizing the Log-Likelihood FunctionLet the random variable ?l be the classifier of ?l, such that if ?l is associated withclass m, m ? {1, 2}, then ?l = m. Also let ? = [?1, . . . , ?L]. Since elements in Xexare assumed independent,Pr(?l = m|?l,?p) =kpmp(?l|?pm)p(?l|?p)=kpm?xi??lp(xi|?pm)?2j=1 kpj?xi??lp(xi|?pj). (4.16)Then, we can write the expectation of the log-likelihood function with respect to theconditional distribution of ? given Xex and the current estimate ?p asL(?|?p) = E [ln (Pr(Xex,?|?)) |Xex,?p] =2?m=1[L?l=1Pr(?l = m|?l,?p)?xi??lln p(xi|?m) +L?l=1Pr(?l = m|?l,?p) ln km],(4.17)where lnx = loge x is the natural logarithmic function.74Chapter 4. LOS and NLOS Classification for UWLAssuming knowledge of ?p, ?p+1 is estimated as the vector of distribution param-eters that maximizes (4.17) while satisfying constraints (4.9), (4.10) and (4.11). Thisprocedure is repeated for Plast iterations, and the convergence of (4.17) to a localmaximum is proven [91]. Then, we calculate Pr(?l = m|?l,?Plast) using (4.16), andassociate vector ?l with the LOS path ifPr(?l = 1|?l,?Plast) > Pr(?l = 2|?l,?Plast) , (4.18)or with an SNLOS path otherwise. Estimation ?Plast and classifications ?l could beused further to improve the accuracy of UWL, e.g., [95, 23]. We observe that the twoterms on the right-hand side of (4.17) can be separately maximized, i.e., given ?p,we can obtain ?p+1m from maximizing the first term, and kp+1m from maximizing thesecond term. Thus (see details in [91]),kp+1m =1LL?l=1Pr(?l = m|?l,?p), m = 1, 2 . (4.19)In the following, we describe the details of our classification procedure for theestimation of ?m, followed by a heuristic approach for the initial estimates ?0.4.3.4 Estimating the Distribution Parameters ?1 and ?2To estimate ?m, we consider only the first term on the right-hand side of (4.17),which for the PDF (4.6) is given byf(?m, ?m, ?m) =L?l=1Pr(?l = m|?l,?p)?xi??lln ?m?ln(2?m)?ln ?(1?m)?(|xi ? ?m|?m)?m.(4.20)Then, considering constraints (4.9), (4.10) and (4.11), we find ?p+1m by solving thefollowing optimization problem:?p+11 ,?p+12 = argmin?1,?2?2?m=1f(?m, ?m, ?m) (4.21a)s.t. : ?1 ? ?2 ? ?1 + TLIR (4.21b)?m??????(3?m)?(1?m) ? TLIR ? 0 , m = 1, 2 (4.21c)?1 ? ?2 ? 0 . (4.21d)We observe that convexity of f(?m, ?m, ?m) depends on ?m. In Appendix A, wepresent an alternating optimization approach (cf. [120]) to efficiently solve (4.21).Next, we present an algorithm to obtain the initial estimation, ?0, whose accuracyaffects the above refinement as well as the convergence rate of the EM algorithm.75Chapter 4. LOS and NLOS Classification for UWL4.3.5 Forming Initial Estimation ?0Our algorithm to estimate ?0 is based on identifying a single group, ?l? , whoseelements belong to the LOS class with high probability, i.e., Pr (?l? = 1) ? 1. Thisgroup is then used as a starting point for the K-means clustering algorithm [91],resulting in an initial classification ?l for ?l, l = 1, . . . , L, to form two classified setsXexm , m = 1, 2. Finally, we evaluate the mean, variance, and kurtosis of the elementsin vector Xexm , denoted as E [Xexm ], Var [Xexm ], and K [Xexm ], respectively, to estimate?0 using the following properties for distribution (4.6):|Xexm ||Xex| = km , E [Xexm ] = ?m ,Var [Xexm ] =?2m?(3?m)?(1?m) ,K [Xexm ] =?(5?m)?(1?m)?(3?m)2 ?3 .(4.22)Since we assume that ?1 < ?2 (see (4.9)), we expect small differences between mea-surements of the LOS link, compared to those of SNLOS links. We use this attributeto identify group ?l? by filtering Xex and calculating the first derivative of the sortedfiltered elements. Group ?l? corresponds to the smallest filtered derivative.4.3.6 DiscussionWe note that the constraints in (4.21) do not set bounds on the values of ?1 and ?2,but rather determine the dependencies between them. This is because, apart fromthe value of TLIR and distribution (4.6), we do not assume a-priori knowledge aboutthe values of km and ?m, m = 1, 2. Without node motion, all elements of Xex belongto one class, but our classifier might still estimate both k1 and k2 to be non-zero,resulting in wrong classification into two classes. In this case, using the average ofthe elements of Xex might give a better estimation of dLOS than ?1.To limit this shortcoming of our classifier, we assume that ?1 and ?2 are distinct ifXex is indeed a mixture of two distributions. To this end, in the last iteration, Plast,we classify Xex as a single class (of unknown type) if the difference |?Plast1 ? ?Plast2 | issmaller than a threshold value, ?v (determined by the system resolution for distinctpaths). Then, if required, we find the distribution parameters of the (single) classby solving a relaxed version of (4.21), setting k1 = 1 and k2 = 0. Nevertheless, wemotivate relevance of our classifier in Section 4.4.2 by showing that scenarios in whichXex is indeed a mixture of two distributions are not rare in real sea environments.4.3.7 Summarizing the Operation of the ClassifierWe now summarize the operation of our classification algorithm, whose pseudo-code is presented in Algorithm 2. First, we evaluate dPDi and dRSS,LBi (lines 1-2). IfdRSS,LBi > dPDi , we classify xi as ONLOS; otherwise, we classify it as non-ONLOS(lines 3-5) and form the vector of non-ONLOS PD measurements, Xex, and groups76Chapter 4. LOS and NLOS Classification for UWLAlgorithm 2 Classifying X1: dPDi := c ? xi2: Calculate dRSS,LBi using RSS measurements, ?UB, ?UB and model (4.2)3: if dRSS,LBi > dPDi then4: Classify xi as ONLOS link5: else6: Exclude ONLOS measurements to form vector Xex and groups ?l satisfying(4.15)7: Estimate ?08: for p := 2 to Plast do9: Calculate kpm, m = 1, 2 using (4.16) and (4.19)10: for i := 1 to Nrepeat do {alternating maximization to solve (4.21) (see Ap-pendix A)}11: Estimate ?p,im , m = 1, 2 and set ?p,i+1m :=?p,im , ?p,i+1m :=?p,im , ?p,i+1m :=?p,im12: end for13: m = 1, 2: ?pm:=?p,Nrepeatm , ?pm:=?p,Nrepeatm , ?pm:=?p,Nrepeatm14: end for15: if |?Plast1 ? ?Plast2 | > ?v then16: Calculate Pr(?l = m|?l,?Plast) and ?l, m = 1, 2, l = 1, . . . , L using (4.16),(4.18)17: else18: Vector Xex consists of a single class. Repeat steps 7-14 for k1 = 1, k2 = 019: end if20: end if?l, l = 1, . . . , L (line 7). Next, we form the initial solution, ?0m (line 7), and run theEM algorithm for Plast iterations (lines 8-14). The procedure starts with estimatingkpm (line 9), followed by an iterative procedure to estimate ?pm for a pre-definednumber of repetitionsNrepeat (lines 10-13). After iteration Plast, we check if vectorXexconsists of two classes (line 15), and determine classifiers ?l, l = 1, . . . , L (line 16);otherwise Xex is classified as a single class (of unknown type), and, if estimating ?mis required, we repeat the above procedure while setting k1 = 1, k2 = 0 (line 18).The software implementation of the above algorithm can be downloaded from [96].Since the data for step 2 in Algorithm 2 is already available through the detectionof each symbol, and since the complexity of the EM algorithm is O (NPlast) (cf.[91]), the complexity of our algorithm is O (NPlastNrepeat). We also note that the EMalgorithm, as well as the alternate maximization process described in Appendix A,provably converge to a local maximum of the log-likelihood function (4.17). In thefollowing, we provide the hybrid Crame?r-Rao lower bound (HCRLB) as a benchmarkfor our classifier.77Chapter 4. LOS and NLOS Classification for UWL4.3.8 Deriving the HCRLBConsider the vector of measurements Xex whose elements are drawn from distribu-tions (4.6) withM = 2 classes (we assume that ONLOS measurements have correctlybeen identified). Our classifier estimates the vector ? = [?1, ?1, ?1, k1, ?2, ?2, ?2] =[?1, . . . , ?7]. We observe that constraints (4.11), (4.9), and (4.10), introduce depen-dencies between pairs (?1, ?5), (?2, ?6), (?r, ?r+1), r = 2, 6, respectively. Thus, wecannot use the conventional Crame?r-Rao Bound to lower bound the variance of anyunbiased estimator of ?. Instead, we apply the HCRLB considering ?1 as a deter-ministic and ?r = [?2, . . . , ?7] a vector of random variables having prior distributions,respectively. The HCRLB is given by [121]EXex,?r|?1[(?Plast ? ?) (?Plast ? ?)T]? H?1(?1) , (4.23)where H(?1) ? <7?7 is the hybrid Fisher information matrix8 (HFIM). Let %i be theclassifier of xi (i.e., %i = ?l if xi ? ?l). Then, the (j, q)th element of the HFIM isH(?1)j,q = E?r|?1 [F (?r, ?1)j,q] + E?r|?1[? ?2??j??qlog p(?r|?1)], (4.24)whereF (?r, ?1)j,q = EXex|?r,?1???|Xex|?i=1?2??j??qlog k%ip(xi|?%i)?? . (4.25)Solving (4.24) requires the calculation ofp(?r|?1) = p(k1)p(?2|?1)p(?2|?1, ?2)p(?1|?1)p(?1)p(?2) . (4.26)Since, as discussed in Section 4.3.6, we do not assume further knowledge aboutthe values of k1 and ?m, m = 1, 2, accounting for constraints (4.7)-(4.11) we as-sume p(?2|?1) is uniform between ?1 and ?1 + TLIR, p(?2|?1, ?2) is uniform be-tween ?1 and TLIR??????(1?2)?(3?2) , p(?1|?1) is uniform between 0 and TLIR??????(1?1)?(3?1) , andp(?m), m = 1, 2, is uniform between 1 and a deterministic parameter, G. Further-more, we assume p(k1) is uniform between 0 and 1. Exact expressions for (4.24) aregiven in Appendix B. For the numerical results presented in the following section weevaluate the HCRLB through Monte-Carlo simulations considering the above uniformdistributions.8Note that while the EM algorithm works on vectors ?l, the actual inputs to our classifier arePD measurements. Thus, in forming the HCRLB, we use xi rather than ?l.78Chapter 4. LOS and NLOS Classification for UWL4.4 Performance EvaluationIn this section, we present results from both computer simulations and sea trials todemonstrate the performance of our classification algorithm. The results are pre-sented in terms of detection probabilities of LOS, SNLOS, and ONLOS links. Inaddition, we measure estimation errors |?pm ? ?m|, |?pm ? ?m|, and |?pm ? ?m|. Wecompare our results to the HCRLB presented in Section 4.3.8, as well as to severalbenchmark methods. The purpose of the simulations is to evaluate the performanceof our classifier in a controlled environment, while results from sea-trial measurementsreflect performance in actual UWACs.4.4.1 SimulationsOur simulation setting includes a Monte-Carlo set of 10000 channel realizations,where two time-synchronized nodes, uniformly randomly placed into a square area of1 km, exchange packets. The setting includes two horizonal and two vertical obsta-cles of length 20 m, also uniformly randomly placed into the square area, such that aLOS always exists between the two nodes. In each simulation, we consider a packetof 200 symbols of duration Ts = 10 msec and bandwidth B = 6 kHz transmitted ata propagation speed of c = 1500 m/sec. To model movement in the channel (dealtwith by forming groups ?l), during packet reception the two nodes move away fromeach other at constant relative speed of 1 m/sec, and dLOS is considered as the LOSdistance between the nodes when the 100th symbol arrives.In our simulations, we use model (4.1) to obtain setX as follows. For each channelrealization and node positions, we find the LOS distance between the two nodes,and determine ?1 = xLOS. Based on the position of nodes and obstacles, we identifyONLOS links as single reflections from obstacles and determine ?3 as the average delayof the found ONLOS links. We use TLIR = 0.1 sec and based on constraint (4.11), werandomize ?2 according to a uniform distribution between ?1 and ?1 + TLIR. For theother distribution parameters ?, we determine ?m, m = 1, 2, 3 as an integer between1 and 6 with equal probability (i.e., G = 6 in (4.26)), and ?m, m = 1, 2, 3 according to(4.8) with ?m uniformly distributed between 0 and (TLIR)2, preserving ?1 < ?2. Basedon model (4.5), we then randomize xi, i = 1, . . . , 200 using distribution (4.6) and auniformly distributed km, m = 1, 2, 3 between 0 and 1, while keeping3?m=1km = 1and setting k3 = 0 if no ONLOS link is identified. Considering the discussion inSection 4.3.6, we use ?v = 1 mc as a detection threshold to check if measurementsin vector Xex correspond to a single link. Additionally, for forming groups ?l (see(4.15)), we use an assumed coherence time T?c = a ? Ts, where a ? {1, 5, 10}, and aquantization threshold ?T = 0.16 msec based on the bandwidth of the transmittedsignals. Note that since the distance between nodes changed by 2 m during receptionof the 200 symbols, if a > 2 m?T ?c ? 8.3, condition (4.15b) is irrelevant, whereas a = 179Chapter 4. LOS and NLOS Classification for UWL10 12 14 16 18 2000.10.20.30.40.50.60.70.80.91Emprirical probability?UB Pd,ONLOSPd,ONLOS (analysis)Pd, non?ONLOSPd, non?ONLOS (analysis)(a) (b)Figure 4.1: Classification results: (a) Prd,non?ONLOS and Prd,ONLOS vs. ?UB, (b) Em-pirical detection probabilities of LOS and SNLOS.results into single-element vectors ?l.To simulate channel attenuation (4.2), we use ? = 15, ? = 1.5 dB/km (consideringa carrier frequency of 15 kHz [5]), and set \u000F to be zero-mean Gaussian with variance5/dB2//?Pa@1m. We use a source power level of 100 dB//?Pa@1m and a zero-mean Gaussian ambient noise with power 20 dB//?Pa@1m, such that the SNR atthe output of the channel is high. Attenuation in LOS and SNLOS links is determinedbased on (4.2), while for ONLOS links we use (4.3) and set RL = 10 dB//?Pa@1m.To obtain the lower bound on RSS-based distance, dRSS,LBi , i = 1, . . . , 200, we use theattenuation model in (4.2) with ?UB = 20 and ?UB = 2 dB/km. An implementationof the simulation environment can be downloaded from [96].First, in Figure 4.1a we show empirical detection probabilities for ONLOS andnon-ONLOS links as a function of ?UB, as well as corresponding results using bounds(4.14) and (4.13). We observe a good match between the empirical results and theanalytical bound for Prd,ONLOS, and that Prd,ONLOS is hardly affected by ?UB. How-ever, Prd,non?ONLOS increases dramatically with ?UB, and the corresponding bound in(4.13) is less tight. This is because choosing ?UB < ? might lead to dRSS,LB > dLOSand neglectance of ? in (4.13) causes analytical inaccuracies. For Prd,ONLOS, however,the large RL is more significant than the effect of ?UB.In Figure 4.1b, we show empirical detection probabilities9 for LOS (LOS (EM))and SNLOS (SNLOS (EM)) links, the total detection probability (ALL (EM)), whichis calculated as the rate of correct classification (of any link). Also shown are classi-fication performance using only the initialization process (init), i.e., before the EMalgorithm is employed, as well as the results for classification without prior identifi-cation of ONLOS links (No ONLOS ID), i.e., without the first step of our algorithm.For the latter, we consider two cases: i) M = 2 and ii) M = 3, where in the9We note that detection probabilities are calculated only when vector X consists of both LOSand SNLOS related PD measurements; classification cannot be made otherwise.80Chapter 4. LOS and NLOS Classification for UWL0 5 10 15 20 25 300.30.40.50.60.70.80.91Emprirical Pr(? err Ts), performance slightlyimproves compared to the case of T?c = Ts. However, a tradeoff is observed as resultsfor T?c = 5Ts are marginally better than for T?c = 10Ts. This is because of erroneous81Chapter 4. LOS and NLOS Classification for UWL0 2 4 6 8 103.23.43.63.84Iteration number, pc|?p 1?? 1|[m]0 2 4 6 8 1023456Iteration number, pc|?p 1?? 1|[m]0 2 4 6 8 1000.511.5Iteration number, p|?p 1??1|[m]0 2 4 6 8 100.050.10.150.20.25Iteration number, p|kp 1?k|[m]Figure 4.3: Estimation error of LOS and SNLOS distribution parameters vs. EMiteration number.assignments to vectors ?l for over-estimated coherence time, T?c. This effect becomesmore significant when in addition to node movements also the channel changes (whichis included in the sea-trial results further below).Convergence of the EM iterative procedure is demonstrated in Figure 4.3, wherewe show average estimation errors of the distribution parameters of the LOS classas a function of the EM iteration step number. We observe that estimates stabilizeafter 10 iteration. While improvement compared to the initialization process (seeSection 4.3.5) is shown for all estimations, the impact of the EM algorithm is mostpronounced for the estimation of k1, which greatly affect classification performance.This improvement is also observed in Figure 4.1b.4.4.2 Sea TrialsWhile our simulations demonstrate good classification performance for our algorithm,the tests relied on the distribution model (4.6), and upper bound on transmission lossmodels (4.2) and (4.3), which might not be faithful representations of realistic UWACsand PD estimators. Thus, we present classification results for UWACs measuredduring three sea trials conducted in Israel and Singapore. One of these experimentswas conducted in a harbor environment to test only ONLOS classification, while theother two were in shallow water to test LOS and SNLOS classification. To acquirePD measurements from recorded sea-trial data, we used the phase-only-correlator(POC) detector as described in [47]. For the ith received signal, xi is estimated asthe first peak at the output of the POC that passes a detection threshold.82Chapter 4. LOS and NLOS Classification for UWL(a) (b)Figure 4.4: Israel harbor experiment: (a) location (picture taken from Google mapson Sep. 09), (b) results for ONLOS link classification.Classifying ONLOS linksIn this section, we show the performance of ONLOS link identification for anexperiment conducted at the Haifa harbor, Israel, in May 2009. The experimentincluded four vessels, each representing an individual node in the network. Here weconsider a subset of the recorded data for which nodes were static. In each vessel,a transceiver was deployed at a fixed depth of 3 m, at a water depth of around6 m. The four nodes were time synchronized using GPS and transmitted with equaltransmission power at a carrier frequency of 15 kHz. Referring to Figure 4.4a, node2 was placed at a fixed location 2A, while nodes 1, 3 and 4 sent packets to node2 while moving between various locations, creating a controlled environment of fivenon-ONLOS and four ONLOS communication links with a maximum transmissiondistance of 1500 m. For each link, (2, j), j ? {1,3,4}, we evaluated (i) dPD as theproduct of an assumed propagation speed of 1550 m/sec and the position of thefirst peak of the POC for the synchronization signal of each received packet, and (ii)dRSS,LB, employing an energy detector over the synchronization signal and using (4.2)for ?UB = 2 db/km and ?UB = 20. We note that results only changed slightly whenalternative methods for obtaining dRSS,LB and dPD were applied.In Figure 4.4b, we present values of dRSS,LB and dPD for each of the 9 commu-nication links. Applying our proposed ONLOS link identification method, all fourONLOS links were correctly classified and there was no false classification of non-ONLOS links. In particular, we observe that for all ONLOS links, dPD is much lowerthan dRSS,LB. Since the latter was obtained from model (4.2) (i.e., a model for LOS)and we assumed spherical propagation (i.e., ?UB = 20 in (4.2)), which is a very looseupper bound on power attenuation in the UAC (e.g., [5], [86]), this difference vali-dates our assumption that the RL level of the reflecting objects (which could havebeen harbor docks, ship hulls, etc.) are sufficiently high to satisfy assumption (4.4).83Chapter 4. LOS and NLOS Classification for UWLClassifying non-ONLOS linksNext, we present results from two separate experiments conducted in open sea:(i) the first along the shores of Haifa, Israel, in August 2010 and (ii) the second in theSingapore straits in November 2011, with water depths of 40 m and 15 m respectively.This is done to demonstrate our classifier?s performance in different sea environments.As communication links were all non-ONLOS links in both experiments, Xex = Xand we only present results for LOS/SNLOS classification.The first sea trial included three vessels, representing three mobile nodes, whichdrifted with the ocean current at a maximum speed of 1 m/sec, and were time-synchronized using the method described in Chapter 2. Throughout the experiment,the node locations were measured using GPS receivers (whose accuracy was up to?5 m), and the sound speed was measured to be c = 1550 m/sec with deviationsof no more than 2 m/sec across the water column. Each node was equipped witha transceiver, deployed at 10 meters depth, and transmitted more than 100 datapackets which were received by the other two nodes. Each packet consisted of 200direct-sequence-spread-sequence (DSSS) symbols of duration Ts = 10 msec and aspreading sequence of 63 chips was used.Table 4.1: Israel sea trial: distribution of estimated values of ?Plastm , m = 1, 2.Measure ? = 1 ? = 2 ? = 3 ? = 4 ? = 5 ? = 6?Plast1 , Tc = Ts 11% 7% 1.8% 1.8% 0% 77%?Plast1 , Tc = 10Ts 16% 11% 1% 1% 0% 67%?Plast2 , Tc = Ts 0% 54% 9% 5% 5% 24%?Plast2 , Tc = 10Ts 0% 57% 6% 5% 5% 25%An estimation parameter of interest is ?Plastm , which determines the type of dis-tribution of the mth class. In Table 4.1, we show the distribution of ?Plast1 and ?Plast2for different values of T?c. As the results are similar for both T?c = Ts and T?c = 10Ts,this implies that clustering PD measurements in vectors ?l does not affect the esti-mated type of distribution. We also observe that the LOS class seems to lean towards?Plast1 = 6, which implies a uniform distribution, while SNLOS measurements clusteraround ?Plast2 = 2, which corresponds to the normal distribution.In Figure 4.5, we show the empirical C-CDF of?err = |cx?? E(di)| , (4.28)where E(di) is the mean of the GPS-based transmitter-receiver distance during thereception of each packet, x? is either ?Plast1 , the average of the PD measurement in X(E(xi)), the minimum of X (min(xi)), or the average of the obtained PD measure-ments after removal of outliers, i.e., the method described in [95] (Outlier).Results84Chapter 4. LOS and NLOS Classification for UWL0 5 10 15 20 25 300.30.40.50.60.70.80.91Emprirical Pr(?err > |an,z|. Otherwise, the10Throughout this chapter upper case bold symbols represent matrices, bold symbols representvectors, ?? represents the measurement of true parameter ?, and ?? represents the estimation of ?.11This assumption is validated using a realistic wave model as well as results from a sea trial.91Chapter 5. A Machine Learning Approach for Underwater DRvessel is moving in almost fixed speed, vi, and projection is not required. However,it is important to note that by (5.1), even when an,z is zero, compensation over both?n and ?n is required.For each time instance tn, we aim to estimate ?n and ?n, and to project thevessel?s local coordinate system about the z-axis followed by its y-axis. The firstprojection converts the set of measurements A? into a projected set A?h with elementsa?hn. The second projection further converts the set of measurements into a projectedset A?h,p with elements a?h,pn . Since after the second projection the local coordinatesystem is expected to be aligned with the horizontal plane, we can readily integrateelements a?h,pn to estimate d?i,j. Furthermore, since the vessel?s orientation with respectto the reference coordinate system is assumed fixed during [tstart, tend], such projectionallows us to estimate the per-time slot change in the vessel?s heading direction withrespect to the initial heading, ?init.By (5.1), the set of measurements A? are affected by the pitch angle of the vessel.Much like the ocean waves, the vessel pitch angle is periodic in nature and likely,?pi4 < ?n