University of Southern California Dissertations and Theses

Ultra-wideband multistatic and MIMO software defined radar sensor networks

(USC Thesis Other)

Ultra-wideband multistatic and MIMO software defined radar sensor networks

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content ULTRA-WIDEBAND MULTISTATIC AND MIMO SOFTWARE DEFINED RADAR SENSOR NETWORKS by Samuel M. Prager A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulllment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (ELECTRICAL ENGINEERING) December 2021 Copyright 2021 Samuel M. Prager Everything takes time. Bees have to move very fast to stay still. David Foster Wallace ii To my family. iii Acknowledgements I am eternally grateful to my faculty advisor M. Moghaddam for encouraging me to pursue my Ph.D. and for being a source of endless support and wisdom while guiding me through this research journey. I thank M. Moghaddam and the rest of my dissertation committee members, K. Chugg, S. Masri, A. Molisch, and J. Stang for their support and guidance. I thank R. Lucas for inspiring me to pursue my Ph.D. and J. Stang for helping to put me on the path towards software dened radio (SDR)-based radars. Many thanks to M. Haynes, D. Hawkins, T. Thrivikraman, M. Lavalle, K. Carpenter and others the NASA Jet Propulsion Laboratory for presenting me with the challenging problems that propelled this research and for helping to inspire solutions. Thank you to J. Fulton, J. Stock, D. O'Leary, and G. Sexstone at the U.S. Geological Survey (USGS) for their collaboration on the UASnow project and assistance in the UAV-SDRadar snow imaging work. Thank you to D. McGrath at Colorado State University for providing GPR and SfM data for the UASnow eld campaign in Cameron Pass, Colorado. Many thanks to R. Akbar and D. Entekhabi at the Massachusetts Institute of Technology and to A. Silva for involving me in the SPCTOR project. A. Silva conceived the architecture to integrate the UAV-SDRadar system into the SoilSCAPE sensor network and R. Akbar developed the path planning algorithms to generate optimal ight paths for the UAV-SDRadar system. Additional thanks to K. Bakian Dogaheh, A. Melebari, and E. Hodges for their assistance in the eld and to the rest of my colleagues in the Microwave Systems, Sensors and Imaging Lab (MiXIL) at the University of Southern California. My Ph.D. work was performed at the University of Southern California funded primarily from the NASA Earth Science Technology Oce through the NASA Earth and Space Science Fellowship 80NSSC18K1421. Additional thanks for support from USGS under award G20AS00019. iv Table of Contents Epigraph ii Dedication iii Acknowledgements iv List of Tables ix List of Figures x Acronyms xiv Abstract xviii 1 Introduction 1 1.1 Signicance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.4 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2 Theoretical Background and Review of Related Work 7 2.1 Software Dened Radar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2 Radar Waveform Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2.1 Linear Frequency Modulated Waveforms . . . . . . . . . . . . . . . . . . . . . 10 2.2.2 Nonlinear Frequency Modulated Waveforms . . . . . . . . . . . . . . . . . . . 11 2.2.3 Stepped Frequency Synthetic Wideband Waveforms . . . . . . . . . . . . . . 12 2.3 Wireless Synchronization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.3.1 NTP, PTP, and Other TWTT-based Methods . . . . . . . . . . . . . . . . . 14 2.3.2 Ultra-wideband Pulse Synchronization . . . . . . . . . . . . . . . . . . . . . . 15 2.3.3 Synchronization using SDR . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.4 Distributed Multistatic and MIMO Radar . . . . . . . . . . . . . . . . . . . . . . . . 17 2.5 Background on Snow Observation Using Radar Systems . . . . . . . . . . . . . . . . 18 2.6 Landmine Detection and UAV-based Radar . . . . . . . . . . . . . . . . . . . . . . . 22 2.7 Half Space Back-Projection Focusing . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.8 Smart Sensor Networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.9 Related Space-borne Missions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3 Software Dened Radar System Development 25 3.0.1 SDRadar Architecture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.0.2 Hardware Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.0.3 MIMO Sensor Network Operation . . . . . . . . . . . . . . . . . . . . . . . . 28 v 4 Ultra-Wideband Synthesis for High-Range-Resolution Software Dened Radar 30 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.2 Synthetic Wideband Waveform Reconstruction . . . . . . . . . . . . . . . . . . . . . 32 4.2.1 Grating Lobe Supression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 4.2.1.1 Non-uniform Frequency Spacing . . . . . . . . . . . . . . . . . . . . 35 4.2.1.2 Grating Lobe Suppression Filter . . . . . . . . . . . . . . . . . . . . 36 4.2.2 Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.2.3 Loopback Verication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.3 Non-Uniform Nonlinear Synthetic Wideband Waveforms . . . . . . . . . . . . . . . . 42 4.3.1 NLFM Waveform Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.3.2 NLFM Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.3.3 Non-Uniform Nonlinear Synthetic Wideband Waveform Reconstruction . . . 45 4.4 UAV Platform Motion Compensation . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.5 Experimental Results and Applications . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.5.1 Re ector Resolution Test I . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.5.2 Re ector Resolution Test II . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 4.5.3 Snow Penetrating Radar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.5.4 SWW Synthetic Aperture Radar . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.6 SWW SDRadar for UAV-based Landmine Detection . . . . . . . . . . . . . . . . . . 60 4.6.1 Sub-Surface SAR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.6.2 Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.6.2.1 Metallic mine . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.6.2.2 Minimum metal mine . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.6.3 Experimental Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 5 Wireless Sub-Nanosecond RF Synchronization for Distributed Ultra-Wideband Software Dened Radar Networks 68 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 5.2 System Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.3 Syntonization and Coarse Synchronization . . . . . . . . . . . . . . . . . . . . . . . . 77 5.4 Fine Synchronization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.4.1 N-Sensor Synchronization Exchange . . . . . . . . . . . . . . . . . . . . . . . 79 5.4.2 Clock Compensation and Transmit Synchronization . . . . . . . . . . . . . . 82 5.4.3 RF Carrier Phase Synchronization . . . . . . . . . . . . . . . . . . . . . . . . 84 5.5 Peak Detection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 5.5.1 Cramer-Rao Lower Bound . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.5.2 Spectral Phase Slope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.5.3 Quadratic Least-Squares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.5.4 Sinc Nonlinear Least-Squares . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 5.5.5 TOF Peak Estimation Algorithm Performance . . . . . . . . . . . . . . . . . 92 5.6 Processing and Exchange of Information Protocol . . . . . . . . . . . . . . . . . . . . 92 5.7 Sensor Localization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 5.8 Experimental Characterization of Synchronization Method . . . . . . . . . . . . . . . 95 5.8.1 3-Array TOF synchronization . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 5.8.2 3-Array TOF Localization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 5.8.3 3-Array TOF Long Range synchronization . . . . . . . . . . . . . . . . . . . . 98 5.8.4 Bistatic Wireless Re ector Test . . . . . . . . . . . . . . . . . . . . . . . . . . 98 5.8.5 3-Sensor Transmit Synchronization . . . . . . . . . . . . . . . . . . . . . . . . 101 vi 5.9 Experimental System Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 5.9.1 Snow Penetrating Radar Test: Monostatic and Bistatic . . . . . . . . . . . . 104 5.9.2 Bistatic Linear Aperture Test . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 5.10 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 5.11 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 6 Characterization of Clock Phase Errors for Distributed Wireless Synchronization Protocol 109 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 6.2 Modeling Oscillator Phase Noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 6.2.1 Simulation of Clock Phase Time Series . . . . . . . . . . . . . . . . . . . . . . 111 6.3 Modeling Synchronized clock PSD . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 6.3.1 Incorporating the SRI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 6.3.2 Prediction of Relative Clock Phase Error at Time of Radar Pulse . . . . . . . 115 6.4 Synchronization Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 6.4.1 TDMA Distributed Sensor Simulation . . . . . . . . . . . . . . . . . . . . . . 118 6.5 Experimental Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 6.6 Distributed MIMO SAR Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 6.6.1 Full MIMO SAR Focusing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 6.6.2 Pair-only MIMO SAR Focusing . . . . . . . . . . . . . . . . . . . . . . . . . . 123 6.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 7 Snow Depth Retrieval with an Autonomous UAV-mounted Software Dened Radar 126 7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 7.2 UAV-SDRadar Payload Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 7.3 SWW Reconstruction and Radar Processing . . . . . . . . . . . . . . . . . . . . . . . 130 7.4 Experimental Design and Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 7.4.1 Snow and Ground Surface Detection Algorithm . . . . . . . . . . . . . . . . . 133 7.4.2 UAV Elevation-Corrected Radargrams Using GNSS . . . . . . . . . . . . . . 134 7.4.3 Comparison with Ground-Based GPR . . . . . . . . . . . . . . . . . . . . . . 137 7.4.4 Snow Pit (Michigan River) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 7.4.5 Preliminary Bistatic Radar Results . . . . . . . . . . . . . . . . . . . . . . . . 143 7.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 7.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 8 UAV-SDRadar Heterogeneous Smart Sensor Network Integration 147 8.1 Autonomous Remote Control of UAV-SDRadar Missions . . . . . . . . . . . . . . . . 147 A Derivation of Frequency Stacking Algorithm 152 B Derivation of Arbitrary NLFM Waveforms 156 C Derivation of TOF Cramer Rao Lower Bound 159 D Derivation of Synchronized Clock Error PSD with Moving Least Squares State Prediction 162 E Snow and Ground Surface Detection Algorithm 169 vii F UAV-SDRadar Snow Measurement Radargrams 173 References 176 viii List of Tables 3.1 USRP E312 Specications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4.1 SWW, 96 % BW eciency, BW t = 5328 MHz . . . . . . . . . . . . . . . . . . . . . . 41 4.2 SWW, 80 % BW Ecency, BW t = 5328 MHz . . . . . . . . . . . . . . . . . . . . . . 41 4.3 NU-SWW, 80 % BW Eciency, BW t = 5328 MHz . . . . . . . . . . . . . . . . . . . 41 4.4 SWW, 50 % BW Ecency, BW t = 5328 MHz . . . . . . . . . . . . . . . . . . . . . . 42 4.5 NU-SWW, 50 % BW Ecency, BW t = 5328 MHz . . . . . . . . . . . . . . . . . . . 42 4.6 Second Re ector Resolution Test Measured Peak Separation . . . . . . . . . . . . . . 55 4.7 Measured sub-surface target depths after focusing . . . . . . . . . . . . . . . . . . . . 67 5.1 3-Array TOF Synchronization Precision. 1000 trials taken over 1000 seconds at 1 GHz 97 5.2 Bistatic Re ector Range Test Target Locations: Actual vs. Measured Range and Echo Characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.3 3-Sensor MIMO TX Sync Statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 7.1 Size and Weight of UAV-SDRadar Peripheral Sensors and Modules . . . . . . . . . . 129 7.2 Statistical comparison of TWT estimated by UAV-SDRadar with ground-based GPR measurements along the Meadow Transect. y Assuming a snow dielectric of 1:41. . . . 140 7.3 Statistical comparison of TWT estimated by UAV-SDRadar with ground-based GPR measurements along the second Meadow Transect. y Assuming a snow dielectric of 1:41.141 7.4 Statistical comparison of TWT estimated by UAV-SDRadar with ground-based GPR measurements along the Forest Transect. y Assuming a snow dielectric of 1:41. . . . . 141 ix List of Figures 1.1 Multi-platform capable dynamic smart sensor network for coherent distributed radar applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.1 Architecture of a prototypical embedded SDR . . . . . . . . . . . . . . . . . . . . . . 8 2.2 General timing diagram for TWTT-based two clock synchronization . . . . . . . . . 13 3.1 SDRadar implemented in USRP E312 hardware . . . . . . . . . . . . . . . . . . . . . 25 3.2 Client-server architecture of SDRadar network . . . . . . . . . . . . . . . . . . . . . 29 4.1 Simulated SWW reconstruction for xed total bandwidth . . . . . . . . . . . . . . . 37 4.2 Fixed total bandwidth SWW autocorrelation as a function of number of sub-pulses . 38 4.3 Fixed sub-pulse bandwidth and frequency spacing SWW autocorrelation . . . . . . . 38 4.4 Loopback Reconstructed SWW Spectrum before and after GLS Filter Correction . . 39 4.5 Loopback SWW as BW t increases . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 4.6 Loopback SWW and NU-SWW shown with and without Hamming GLS lter vs. Theoretical LFM Chirp . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.7 Signal characteristics of an amplitude weighted LFM chirp and the equivalent constant- amplitude NLFM chirp. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.8 Characteristics of NLFM waveforms constructed from equivalent weighted LFM waveforms for a number of common window functions. . . . . . . . . . . . . . . . . . 45 4.9 NU-NLSWW spectrum synthesis for 350 MHz total NU-NLSWW bandwidth, N = 9 sub-pulses, and Hamming window weighting. . . . . . . . . . . . . . . . . . . . . . . 49 4.10 The reconstruction of Nonlinear SWWs using the existing frequency stacking (FS) algorithm and the proposed non-uniform frequency stitching (NUFS) algorithm . . . 49 4.11 NU-NLSWW reconstructed from N = 26 sub-pulses obtained using USRP E312 SDRadar hardware in loopback conguration . . . . . . . . . . . . . . . . . . . . . . 50 x 4.12 Simulated eect of platform motion on SWW . . . . . . . . . . . . . . . . . . . . . . 53 4.13 Corner Re ector Resolution Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 4.14 Second Corner Re ector Resolution Test . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.15 Snow penetrating radar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.16 Test setup and reconstructed SAR imagery of parking structure scene . . . . . . . . 59 4.17 Imaging geometry for an airborne nadir-looking SAR . . . . . . . . . . . . . . . . . . 60 4.18 Simulation scene geometry and B-Scan image . . . . . . . . . . . . . . . . . . . . . . 62 4.19 Raw and focused B-Scan images obtained from simulation of one and two PEC mine-like cylinders buried at 50 cm BGS in dry sand. . . . . . . . . . . . . . . . . . . 63 4.20 BP and HSBP focused B-Scans for simulated VS-50 style minimum-metal AP mine buried at 10 cm and 50 cm BGS in dry sand. . . . . . . . . . . . . . . . . . . . . . . 64 4.21 Metallic mine-like target experimental results . . . . . . . . . . . . . . . . . . . . . . 65 4.22 Unfocused and Focused GPR B-scan imaging results across 1.5 m aperture with metallic-mine like target on ground surface and at 8 cm and 15 cm BGS. . . . . . . 66 5.1 Local time of ight measurements through cables between two sensors . . . . . . . . 74 5.2 Synchronization timing diagram for a single synchronization epoch . . . . . . . . . . 79 5.3 Illustration of synchronization scheme for three sensors . . . . . . . . . . . . . . . . . 81 5.4 Performance comparison of cross-correlation fractional peak estimation methods . . 92 5.5 Synchronization Processing diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 5.6 Architecture and performance of software dened FEC codec architecture implemented in SDR hardware for synchronization message exchange. . . . . . . . . . . . . . . . . 94 5.7 Wireless synchronization test: three array of sensors in equilateral triangle formation 95 5.8 Three sensor array synchronized TOF Matrix measurements for 1000 trials taken over 1000 s demonstrating sub-100 ps precision . . . . . . . . . . . . . . . . . . . . . 96 5.9 Sensor position estimation using three sensor array and TOF obtained using proposed synchronization method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 5.10 3-Sensor synchronization measurements performed at Sante Fe Dam RC Aireld in Irwindale, CA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 5.11 Bistatic radar test using wireless synchronization algorithm to synchronize two sensors100 xi 5.12 3-Sensor MIMO transmit synchronization test results . . . . . . . . . . . . . . . . . . 102 5.13 Monostatic and bistatic snow penetrating radar test . . . . . . . . . . . . . . . . . . 105 5.14 Wireless two-sensor bistatic aperture measurement . . . . . . . . . . . . . . . . . . . 106 6.1 Oscillator output noise PSD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 6.2 Multiple time domain realizations of random clock phase error with lower quality TCVCXO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 6.3 Analytical synchronized output PSDs as a function of T p . . . . . . . . . . . . . . 114 6.4 Simulated clock phase error PSDs after synchronization . . . . . . . . . . . . . . . . 118 6.5 Simulated clock phase error of each sensor before and after synchronization at time of radar transmission . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 6.6 Simulated clock phase error PSDs after synchronization using the least squares (LS) prediction method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 6.7 Simulated clock phase error of each sensor with and without LS prediction . . . . . . 119 6.8 Experimentally obtained PSDs of non-synchronized TOF (input PSD) and synchro- nized radar pulse TOF delay T p for 3 sensors . . . . . . . . . . . . . . . . . . . . . 120 6.9 Simulated 3x3 MIMO SAR scene and unfocused radar data . . . . . . . . . . . . . . 121 6.10 Simulated 3x3 SAR focusing for all pairs of distributed radars with independent oscillators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 6.11 Simulated 3x3 MIMO SAR focusing for distributed radars with independent oscillators123 6.12 Focused 3x3 MIMO SAR target response in azimuth . . . . . . . . . . . . . . . . . . 123 6.13 Simulated 3x2 MIMO SAR focusing for distributed radars with independent oscillators124 6.14 Focused 3x2 MIMO SAR target response in azimuth for bistatic pairs only . . . . . 124 7.1 Illustration of sUAS-based UAV-SDRadar system performing synthetic ultra-wideband radar measurements of snow depth . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 7.2 UAV-SDRadar system and sled GPR team operating at Cameron Pass, CO test site 128 7.3 Diagram of sUAS system used for ight experiments . . . . . . . . . . . . . . . . . . 129 7.4 Radar processing diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 7.5 Photos of Forest and Meadow snow eld transects and GPS positions of sled GPR and UAV-SDRadar for 2020 eld campaign at Cameron Pass, CO . . . . . . . . . . . 132 xii 7.6 Meadow Transect radargrams before and after surface topography removal . . . . . . 135 7.7 Forest Transect radargrams before and after surface topography removal . . . . . . . 136 7.8 Along-prole comparison of GPR and UAV-SDRadar TWT results . . . . . . . . . . 137 7.9 Scatter plot comparison of GPR and UAV-SDRadar TWT results for dierent footprint (FP) sizes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 7.10 Along-prole comparison of UAV-SDRadar TWT measurements with the nearest GPR measurement for a given UAV lat/lon GNSS coordinate . . . . . . . . . . . . . 138 7.11 Scatter plot comparison of GPR TWT measurement closest to the UAV for a given position . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 7.12 Repeated Meadow Transect ights . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 7.13 Snow pit measurement results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 7.14 Bistatic UAV-SDRadar systems at eld site in Winter Park, CO. . . . . . . . . . . . 143 7.15 Preliminary results from bistatic UAV-SDRadar test ight . . . . . . . . . . . . . . . 144 8.1 Hardware diagram of UAV-SDRadar integration into SoilSCAPE/SPECTOR smart heterogeneous sensor network . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 8.2 Process sequence diagram for integration of UAV-SDRadar with SoilSCAPE/SPEC- TOR smart heterogeneous sensor network . . . . . . . . . . . . . . . . . . . . . . . . 148 8.3 Integrated UAV-SDRadar system with onboard lidar altimeter, Nvidia Jetson TX-2 based DJI Manifold SBC, and remote power control module . . . . . . . . . . . . . . 149 8.4 Flight paths and radargram data collected on autonomous ights at SPCTOR eld test site in Walnut Gulch, AZ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150 8.5 Flight 1 repeated over multiple days before and after rainfall at test site in Walnut Gulch, AZ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 F.1 Meadow and Forest Transect radargrams prior to UAV altitude correction . . . . . . 174 F.2 Meadow and Forest Transect radargrams after RTK GNSS/GPS UAV elevation correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175 xiii Acronyms ADC analog-to-digital converter AGL above ground level AWG arbitrary waveform generator AWGN additive white Gaussian noise BP back-projection BRAM block random access memory BW bandwidth COTS commercial-o-the-self CRLB Cramer-Rao Lower Bound CSAC chip scale atomic clock DAC digital-to-analog converter DDC digital down converter DFT discrete Fourier transform DMA direct memory access DOCXO double oven controlled crystal oscillator FEC forward error correction FIFO rst in, rst out FMCW frequency modulated continuous wave FPGA eld programmable gate array GNSS global navigation satellite system GPR ground penetrating radar GPS Global Positioning System HSBP half-space back-projection IFT inverse Fourier transform IMU inertial measurement unit xiv INS inertial navigation system IP Internet Protocol LFM linear frequency modulated LO local oscillator LOS line of sight LS least squares MIMO multiple input multiple output MPSoC Multiprocessor System on a Chip NASA National Aeronautics and Space Administration NL-LS nonlinear least squares NL-SWW nonlinear SWW NLFM nonlinear frequency modulated NTP network time protocol NU-NLSWW non-uniform nonlinear SWW NU-SWW non-uniform SWW NUFS non-uniform frequency stitching OCXO oven controlled crystal oscillator PL programmable logic PLL phase-locked loop PPK post-processing kinematic PPS pulse per second PRI pulse repetition interval PS processing system PSD power spectral density PSL peak sidelobe level PSP Principle of Stationary Phase PTP precision time protocol QPSK quadrature phase-shift keying RF radio frequency xv RFNoC RF Network on Chip RFSoC RF System on a Chip RPC radar pulse controller RTK real-time kinematic RX receive SAR synthetic aperture radar SDR software dened radio SDRadar software dened radar SfM structure from motion SNR signal-to-noise ratio SoC system on a chip SoilSCAPE Soil moisture Sensing Controller and oPtimal Estimator SPCTOR Sensing-Policy Controller and OptimizeR sUAS small unmanned aircraft system SWaP size, weight, and power SWE snow water equivalent SWW synthetic wideband waveform TCP Transmission Control Protocol TDMA time-division multiple access TOA time of arrival TOF time of ight TWT two-way travel time TWTT two-way time transfer TX transmit TxBW time bandwidth product UAV unmanned aerial vehicle UDP user datagram protocol UGV unmanned ground vehicle UHD USRP Hardware Driver xvi USRP Universal Software Radio Peripheral UTC Coordinated Universal Time UWB ultra-wideband UXO unexploded ordnance WSN wireless sensor network xvii Abstract Multistatic and multiple input multiple output (MIMO) radar sensor networks are an area of signicant interest due to spatial diversity gains, improved resolution, and higher dynamic range as well the potential for these systems to reduce overall mission costs and improve survivability. Networks of autonomous small unmanned aerial vehicles (UAVs) and rovers, or unmanned ground vehicles (UGVs), and small CubeSat constellations represent exciting potential as deployment platforms for distributed radar systems. However, the small payload of capacity of such platforms and large number of individual sensors that constitue a distributed network necessitates the development of small, low cost, and ecient radar sensors that are capable of high levels of performance. A compelling application of such sensor networks is in payload directed autonomy and dynamic recongurability for dierent radar sensing modes and objectives, requiring sensors that are also highly exible and intelligent. Furthermore, in order achieve coherent multistatic/MIMO radar operation in a wireless sensor network, precise synchronization to sub-100 picosecond levels between nodes is required and represents a fundamental obstacle to the realization of such systems. This research focuses on the development of software dened radar (SDRadar) algorithms and techniques that address two primary problems in this area. We rst present a synthetic ultra- wideband waveform and reconstruction technique that that allows low-cost software dened radio (SDR)-based SDRadars to eciently achieve cm-level radar resolution from real-time dynamically congurable spectrum allocations. Next, we report a novel O(N) distributed and decentralized wireless synchronization scheme that achieves time synchronization to sub-100 picosecond precision across all nodes and enables time and phase coherent multistatic/MIMO radar operation. The SDRadar system and techniques presented in this work are validated extensively through experimental results. We rst demonstrate the performance of the synthetic wideband waveform reconstruction techniques and the exibility of the implemented synthetic wideband SDRadar in multiple radar imaging modes, including synthetic aperture radar (SAR) and ground penetrating radar (GPR). We then report results from multiple SDRadar synchronization experiments for both two and three sensor networks, demonstrating sub-100 picosecond performance, cm-level localization, xviii and the ability of the method to coherently synchronize sensors across multiple stepped frequency bands for bistatic synthetic wideband radar imaging. We derive an analytic formulation of the synchronization protocol in terms of the sensor input oscillator's phase error power spectral density (PSD) and validate the model with further experimental results. Finally, we describe the integration of the SDRadar as a UAV payload and share experimental results for the UAV-SDRadar system in autonomous UAV-based radar imaging from multiple eld campaigns. xix Chapter 1 Introduction 1.1 Signicance There is, as of yet, no clear answer to the question of what next generation remote sensing instruments and platforms will look like. As scientic observation continues to be dominated by expensive, complex space and air-borne systems, it is vital to gain understanding of the potential for emerging technologies to reduce costs and maximize scientic yield in future missions. The 2014 National Aeronautics and Space Administration (NASA) Science Plan emphasizes the importance of developing and validating such new observation technologies before using them in science missions [1]. Wicks provides the following prediction for next-generation radar sensors in [2], writing: \The sensors will dynamically and automatically change waveform parameters to accomplish [..] goals. Disparate sensors will communicate and share data and instructions in real-time. Intelligent sensor systems will operate within and between sensor platforms such that the integration of multiple sensor data provides information needed to achieve dynamic goals." This research eort seeks to answer this question for a subset of potential next generation instruments: smart radar sensor networks made up of small unmanned aerial vehicles (UAVs), unmanned ground vehicles (UGVs), and in-situ ground-based sensors by developing a low cost, platform-independent smart sensor payload that can act as a node and work dynamically as a non-centralized task-oriented wireless network and demonstrating the ecacy of the network in scientically valuable remote sensing tasks. 1 1.2 Objectives Due to the fundamental geometric spreading that occurs in all radiated electromagnetic waves, as described by the inverse square law, the received power of an active radar is, for a point target, proportional to the inverse fourth power of the target range. This result is described in the well- known radar range equation relating received power P rx to transmit power P tx as a function of antenna gainG, target scattering coecient, wavelength, target rangeR and medium-dependent loss L(R) [3], given in the monostatic form as P rx = P tx G 2 2 (4) 3 R 4 L(R) Thus, low-altitude small UAV and ground-based UGVs and in-situ radar sensors are fundamentally capable of achieving comparable signal-to-noise (SNR) requirements at signicantly lower power than air and space-borne counterparts. If we consider noise to be statistically gaussian with zero mean, coherent pulse averaging can be shown to improve SNR linearly, such that if N avg pulses are coherently averaged, SNR/N avg . For typical space and air-borne radars,N avg is highly constrained by the high platform velocity. However, small UAV, UGV and in-situ ground-based sensors are able to hover or otherwise remain in place, allowing long integration times and resultant SNR gains. Furthermore, such small, low-power sensor platforms enjoy vastly reduced deployment and operational costs. We therefore oer the hypothesis that smart sensor networks made up of UAV, UGV and in-situ radar nodes have broad potential to perform scientic observation of the earth at signicantly reduced cost in comparison to space and air-borne single instrument counterparts. Thus, we arrive at the following objectives in an eort to validate the posed hypothesis: (i) implement a smart radar sensor network capable of multi-mode operation using a single multi- platform compatible instrument and (ii) demonstrate the eectiveness of the proposed smart network in performing useful scientically important remote sensing tasks. To this eect, we propose the following operational modes for the radar sensors: • Sounding: phase-sensitive pulse-compression radar estimates distance to medium boundaries to centimeter-precision. Applications in detection of water-table, sub-glacial lakes, permafrost, etc. 2 (a) Multi-platform (b) UAV swarm (c) Rover eet (d) CubeSat constellation Figure 1.1: Multi-platform capable dynamic smart sensor network for coherent distributed radar applications including synthetic array formation and coherent beam steering • Multistatic: Single transmitter node and multiple receiver nodes • Beam-steering: Multiple phase-coherent MIMO nodes transmit and receive with variable relative delay and phase to enable beam-steering. • Array: Single and/or multiple coherent MIMO nodes transmit and receive, sampling spatial doppler suciently to meet requirements for synthetic aperture formation. • Coherent MIMO: Multiple nodes simultaneously transmit and receive coherently. This is a generalization of both array and beam-steering operational modes. • Statistical MIMO: Multiple nodes transmit and receive orthogonal waveforms simultaneously from disparate locations, resulting in spatial diversity gains. We illustrate the proposed sensor payload operating in a variety of modes on multiple platforms in gure 1.1. There are a wide range of remote sensing based scientic endeavors that could stand to benet from the proposed sensor network. We illustrate the following examples: water-table monitoring with platform-independent network (gure 1.1a), MIMO altimeter radar imaging of sub-surface layers (gure 1.1b), and mapping of sub-glacial lakebeds (gure 1.1c). Software dened radar (SDRadar) based sensors enjoy unparalleled exibility and dynamic- congurability in comparison to traditional designs and have exciting potential as both development and deployment platforms for smart radar sensor networks, which may be deployed across a range of host platforms to perform sensing and imaging tasks optimally. SDRadar allows for dynamic optimization of waveforms, frequency, power and topology. 3 There exist a broad range of challenges that must be overcome in the realization of the proposed eort. Chief among these challenges are • Development of high resolution methods for low-cost radar hardware • Robust and ecient synchronization of nodes (both hardware clock time and Local Oscillator (LO) phase). • Obtaining sucient time and phase synchronization stability necessary for wireless MIMO radar operation. • Enhancing onboard processing capabilities and algorithms for Multistatic and MIMO radar tomographic imaging. Thus, it is the aim of the proposed work to address these challenges and demonstrate the ecacy of their solutions in achieving the stated objectives. 1.3 Contributions This work contains a number of original contributions to the elds of SDRadar, wireless syn- chronization, and small platform distributed radar imaging. We expand upon existing synthetic ultra-wideband techniques with the development of a frequency stacking algorithm for synthetic wide- band waveform (SWW) reconstruction with a grating lobe suppression lter, non-uniform frequency spacing, and platform motion compensation that enables cm-level radar resolution performance using low-cost commercial-o-the-self (COTS) software dened radio (SDR) hardware. We introduce a novel non-uniform nonlinear SWW (NU-NLSWW) that features tunable sidelobe performance and constant waveform amplitude allowing use with highly ecient saturation-mode RF ampliers. We further develop of new non-uniform frequency stitching (NUFS) algorithm for reconstruction of arbitrary SWWs, including non-uniform, nonlinear, and non-uniform nonlinear SWWs with signicantly reduced or eliminated grating lobe contamination. Together these contributions enable the use of COTS SDR-based SDRadar sensors to produce scientically useful data products in high resolution imaging applications. 4 We develop of a novel distributed and decentralized N-sensor wireless synchronization technique that achieves sub-100 picosecond synchronization across wireless sensor networks using only 50 MHz bandwidth signals and achieves coherent-on-transmit operation of all sensors enabling numerous distributed wireless MIMO radar and radio operational modes. We further derive and original analytic formulation of the transfer function applied to input clock phase error power spectral density (PSD) by the decentralized N-sensor wireless synchronization scheme enabling accurate and ecient models and simulations of large-scale distributed radar systems across entire parameter space from error PSDs. We report the rst demonstration of radar retrieval of snow depth from an small unmanned aircraft system (sUAS) UAV platform over forested regions and rst integration of structure from motion (SfM), lidar altimeter, and radar data in snow depth retrieval from an sUAS UAV platform. This information is integrated in an original snow and ground surface detection algorithm to produce two-way time transfer (TWTT) estimates as a proxy for snow depth. We also present preliminary results from the rst reported coherent bistatic radar measurement of snow from two cooperative autonomous sUAS UAV platforms, which use the wireless synchronization scheme described in this work to synchronize in real-time during ight. Finally, we develop an autonomous UAV-SDRadar sub-system and backend for the rst demonstration of autonomous UAV-based radar integration with a wireless sensor network for persistent deployment and measurement as part of the Sensing-Policy Controller and OptimizeR (SPCTOR) project. 1.4 Outline The remainder of this document will be organized as follows: Chapter 2 summarizes existing literature and previous work related to the this research eort. Chapter 3 describes the development of an embedded battery powered SDRadar system based on COTS Universal Software Radio Peripheral (USRP) SDR hardware. Chapter 4 reports novel approached to coherent ultra-wideband stepped- frequency techniques for synthesizing wideband waveforms and obtaining cm-level radar resolution with low-cost SDR hardware as well as experimental demonstrations of the SDRadar system. A novel decentralized wireless synchronization scheme capable of < 100 ps synchronization across a network of SDRadar sensors is presented in Chapter 5, with extensive experimental validation of 5 the technique. Chapter 6 provides an analytic characterization of the distributed synchronization method and presents analytic formulations of synchronized clock error PSD given an arbitrary input oscillator phase error PSD. Eorts to integrate the SDRadar sensor as a UAV payload and results from a eld campaign in Colorado, USA to measure snow depth using the integrated UAV-SDRadar system are presented in Chapter 7. Finally, the integration of the UAV-SDRadar system as an autonomous agent for persistent remote deployment in the eld within a heterogeneous network of sensors as part of the Soil moisture Sensing Controller and oPtimal Estimator (SoilSCAPE) and SPCTOR projects is describd in Chapter 8. 6 Chapter 2 Theoretical Background and Review of Related Work In this chapter, summaries of existing literature related to this eort are provided. Specically, this section will focus on research eorts involving high resolution software dened radar, wireless synchronization, distributed multistatic/MIMO radar imaging, and smart sensor networks. 2.1 Software Dened Radar The summary in this section has been published in [4]. SDRadar has been suggested as a next-generation platform that may address problems in conventional radar systems involving design complexity, cost and recongurability. SDRadar systems are inherently capable of multipurpose operation and real-time reconguration. Common hardware designs reduce cost and system implementation time, making exploration of SDRadar technologies an area of signicant research interest [5{7]. Recently, a number of multi-function radars missions have employed SDRadars and software-dened techniques. For example, weather radars from the U.S Department of Energy and the University of Oklahoma Advanced Radar Research Center, among others, have adopted SDRadar-based receivers for their operational exibility and low cost [8,9]. Due to similar hardware needs and the exibility aorded by software, SDRadars are commonly implemented in commercial SDR systems. The concept of SDR was rst introduced publicly in 1992 by Mitola in [10]. The SDR paradigm calls for digital signal processing to be performed completely in software and eld programmable gate array (FPGA) rmware, with a hardware 7 frontend that performs analog-to-digital (ADC) and digital-to-analog (DAC) conversion as well as radio frequency (RF) mixing, amplication and ltering [11]. A typical SDR architecture, with digital up/down-converters (DUC/DDC) implemented in FPGA fabric, is illustrated in Figure 2.1. With appropriate considerations, commercial SDRs can provide inexpensive and accessible platforms in which SDRadars can be implemented. Figure 2.1: Architecture of a prototypical embedded SDR. SDRadar sensors may be implemented in SDR hardware The common approach to SDRadar design follows the SDR paradigm of generating waveform samples in software on a computer and streaming them over a physical interface to an FPGA for up-sampling and digital-to-analog conversion. In this case, the transmitted signal bandwidth is often limited by the data transfer rate of the interface, requiring a DUC/DDC to reduce bandwidth resulting in relatively poor range resolution performance [12]. In essence, common inexpensive modern processors cannot keep up with FPGA sampling rates. Furthermore, the time delay introduced by data transfer between the FPGA and processor is non-deterministic which poses a problem for time-coherent radar systems that require precise timing to be known [13]. Previous work on SDRadar implemented in SDR hardware has focused on (i) improving resolution by increasing the FPGA-to-processor data transfer rate with more capable computers and faster interfaces [14] and (ii) measuring and correcting non-deterministic time delays that cause a loss of inter-pulse time coherence [13,15]. For repeated waveforms and pulsed operation, this approach consumes excess data bandwidth and unnecessarily introduces a bottleneck in the transmitted signal bandwidth. Furthermore, current SDRadar designs remain constrained by the performance of commercially available SDR hardware. 8 The USRP series of SDR hardware platforms are used widely for SDRadar and have been used across a diverse range of remote sensing applications from ionospheric sounding radar [16] to ground penetrating radar [6]. However, low instantaneous bandwidth and, therefore, low range resolution are cited as critically limiting factors for applications that require high resolution [17]. Despite having low instantaneous bandwidth, USRP SDRs are capable of RF frontend band tuning over a wide range. We explore how to exploit this capability to improve SDRadar resolution. Unlike conventional radars, communications SDR boards are designed without the need for coherence. In our papers [4,18], we demonstrate that these coherency challenged SDRs may be used for high resolution remote sensing applications such as ground penetrating radar (GPR), non destructive testing, environmental monitoring and imaging, and are suitable for deployment on small relatively slow moving platforms such as rotor-based UAVs. As such, this is an enabling technology for swarms of low cost radar sensors that can be dynamically recongurable and have extremely high bandwidth with the potential to use any waveform. 2.2 Radar Waveform Design Fundamentally, a radar system operates by measuring the characteristics, including delay (and therefore range), phase, and/or doppler shift of a transmitted signal or after it has been re ected by some scattering target and received by the radar system. The fundamental relationship between delay and range R in a channel medium with relative dielectric permittivity r (where r = 1 for free-space) is [19] R = p r c 2 (2.1) where c = 2:998e8 m/s is the speed of light in a vacuum. There exists a diverse set of waveforms used commonly in radar applications. Broadly speaking, these waveforms include amplitude, frequency, and phase modulated signals as well as stepped frequency signals and combinations thereof. In general, radar waveforms may be categorized as 9 short-time pulse signals, which have range resolution performance based on the length of the pulse, and longer pulse-compression signals, in which case the range resolution performance is determined by the signal bandwidth [19]. The achievable radar range resolution r is dependent on the signal bandwidth B (which for a short-time pulse signal is related to the pulse-length T p as B = 1=T p ) [20] r = p r c 2B (2.2) The use of pulse-compression waveforms allows the radar system to operate at lower peak-power and time-domain compression of the received signal is commonly realized via application of a matched lter. The matched lter is the optimal lter for maximizing the signal-to-noise ratio (SNR) of a signal in the presence of additive white Gaussian noise (AWGN). The impulse response of the matched lter is dened as the time-reversed complex conjugate of the signal to which it is matched [19]. Therefore, the output of the matched lter for zero delay is represented as the signal autocorrelation function. For a pulse-compression signal having length T p , the application of the matched lter yields a signal processing gain that is equal to the time bandwidth product (TxBW) is T p B n where B n is the receiving system noise bandwidth [20]. For a quadrature complex baseband digital receiver with sampling frequency f s , B n =f s . 2.2.1 Linear Frequency Modulated Waveforms The linear frequency modulated (LFM) chirp waveform is the oldest and most popular pulse compression waveform used in radar applications. Use of the waveform dates back to World War II where it was independently conceived by German, British, and U.S. radar engineers [19]. The complex envelope of an LFM chirp waveform x(t) having bandwidth B and pulse length T p may be represented for some amplitude weighting function A(t) as x(t) =A(t) exp (jKt 2 ) (2.3) 10 where K = B=T p is the chirp rate. In the nominal case, A(t) = rect(t=T p ), which yields a peak sidelobe level (PSL) of 13:5 dB down from the peak signal level. A window function is often chosen for A(t), which apodizes the signal amplitude at the head and tail, to reduce the PSL at the cost of resolution performance [21]. Commmonly selected window function include Hamming, Chebyshev, Tukey, and Kaiser windows among others [19]. 2.2.2 Nonlinear Frequency Modulated Waveforms The summary in this section has been previously published in [22]. Advantages in pulse compression with nonlinear frequency modulated (NLFM) waveforms over LFM are widely known and explored [23{25]. The most common method for reducing LFM waveform side-lobes in radar systems is with the application of an amplitude weighted window function. In modern radar systems, however, highly ecient high power ampliers (HPAs) which operate near saturation, for example class-D ampliers, are desired. In such systems, amplitude weighting of transmit waveforms is undesirable or impossible as the HPA will distort or remove the amplitude weighting when operating eciently in saturation mode. Thus, amplitude weighted window functions are most often applied at the receive matched lter's reference waveform, which tapers the waveforms amplitude near the beginning and end of the pulse. This amplitude apodization, however, is inecient as the power contained around the pulse edges is transmitted and then discarded, leading to a reduction in SNR [25]. NLFM waveforms can be designed such that their autocorrelation functions exhibit signicant improvements in side-lobe level over LFM counterparts, while oering improved eciency and resilience to amplitude distortion. We note, however, that caution must be taken when using NLFM waveforms in the presence of signicant target doppler shifts [23]. We may represent the complex envelope of a general constant amplitude NLFM waveform as x(t) =A expf(j(t))g (2.4) 11 where the phase function (t) of the signal is determined by the time domain frequency function f(t) as (t) = Z t Tp=2 f()d (2.5) Previous work has shown that the Principle of Stationary Phase (PSP) can be used to shape the PSD of a NLFM waveform with constant amplitude to resemble that of a windowed LFM [23]. This work has been extended to obtain solutions for a number of specic common window functions [24]. While these works demonstrate signicant side-lobe level improvements over the non-windowed LFM autocorrelation function, the results do not match the performance of the desired windowed LFM waveforms. This is caused by dierences in the PSD of the NLFM and LFM waveforms, possibly due to the approximations assumed by the PSP [24]. In general, the design strategy for NLFM waveforms is as follows [25]: 1. Consider a windowed LFM waveform with a desired PSD and autocorrelation function. 2. Apply the PSP to obtain an expression for a similar PSD that is a function of phase only. 3. Integrate to obtain a group time delay. 4. Compute the inverse of the group time delay function to obtain the time domain frequency function. 2.2.3 Stepped Frequency Synthetic Wideband Waveforms The summary in this section has been published in [4]. Frequency stepped radar is commonly used as a means of obtaining high range resolution with limited instantaneous bandwidth [19, 20]. However, stepping across a wide band in steps small enough to reconstruct the band from pure tones is inecient in USRP SDR hardware due to local oscillator (LO) tuning times that are long relative to the analog-to-digital converter (ADC) sampling and data transfer rates. It is most time-ecient to use all of the available instantaneous bandwidth at each step and minimize LO re-tunings. However, using wideband sub-pulse waveforms and increasing the frequency step size gives rise to undesirable grating lobes [19,26]. 12 A common method for producing a SWW, often called stepped-chirp or stepped-waveform, aims to reconstruct a wideband LFM chirp from a set of narrowband LFM chirp sub-pulses, which is then pulse compressed to achieve high range resolution [27{29]. The impulse-like shape of the grating-lobes is due to periodicity in the SWW spectrum caused by uniform spacing of sub-pulses in the frequency domain by a xed value f c . It has been shown that by spacing sub-pulses at non-uniform frequencies, the energy in the grating-lobes is spread out in the resulting non-uniform SWW (NU-SWW) [26,30,31]. 2.3 Wireless Synchronization The summary in this section has been published in [32]. The problem of time synchronization in general, and wireless synchronization specically, is relevant to a wide range of distributed communications, computer networking, test and measurement, and remote sensing applications. As such, a signicant amount of research eort has focused on this topic and this problem has been extensively studied in previous literature. However, the precision required for most communications and computer systems is often orders of magnitude less stringent than that required for coherent radar systems [33]. T 1 T 2 T 3 T 4 ∆t 1 ∆t 2 Clock 1 Clock 2 Figure 2.2: General timing diagram for TWTT-based two clock synchronization The problem was rst pondered by Henri Poincar e and Albert Einstein within the context of relativistic event simultaneity [34], [35], wherein the fundamental of TWTT synchronization was introduced. As shown in Fig. 2.2, a signal (or packet) is sent from an initiating source at time T 1 . It is received at time T 2 after a delay of t 1 = T 2 T 1 and responded to (or re ected) after a known delay at time T 3 , arriving back at the source at time T 4 . The time oset of the clocks is then ((T 2 T 1 ) (T 4 T 3 )=2 = (t 1 t 2 )=2 and the propagation delay is 13 ((T 2 T 1 ) + (T 4 T 3 )=2 = (t 1 + t 2 )=2 and the now known delay may be used to correct future time exchanges. The concept of TWTT provides the basis for synchronization of modern systems and networks, including satellites and the internet [36]. 2.3.1 NTP, PTP, and Other TWTT-based Methods In large scale modern computer networks, the network time protocol (NTP) and precision time protocol (PTP) are ubiquitous. In NTP, timestamps are exchanged as user datagram protocol (UDP) packets to estimate round trip delay and clock osets. NTP timestamps are generated in software with non-deterministic latencies and NTP provides synchronization accuracy of 10 s [37]. PTP (IEEE-1588) improves over NTP by generating timestamps in the hardware layer. Clock alignment in PTP is achieved by sending many `sync' packets at the expense of signicant network trac. The implementation is limited by the 125 MHz rate of the hardware timestamp counter to an accuracy of 10 ns [38], [39]. The fundamental process of clock synchronization using two way time transfer is depicted in Fig. 2.2. The vast majority of synchronization protocols in the literature are extensions of the fundamental concepts used in NTP and PTP. The White Rabbit protocol, developed at the European Organization for Nuclear Research (CERN), is an ethernet based protocol that uses hardware-based Digital Dual Mixer Time Dierence (DDMTD) phase detectors and specialized network switches to detect and correct fractional clock phase osets to to achieve 10 ps accuracy [40] across large wired networks. The common wireless sensor network (WSN) synchronizations solutions include reference broad- cast synchronization (RBS) [41], timing synch protocol for sensor networks (TPSN) [42] and ooding time synch protocol (FTSP) [43]. Both RBS and TPSN neglect signal time of ight over the wireless channel. In RBS receiving sensors exchange time of arrival measurements for a synchronization packet sent between two reference nodes and a set of receiving sensors synchronize with one another rather than a transmitter [41]. In pairwise broadcast synchronization (PBS), it is assumed that all nodes can `overhear' pair-wise synchronization messages between other nodes. PBS still requires a hierarchical node of structure [44]. Additionally, PBS synchronization assumes that distances, i.e., signal times of ight between nodes are known a priori. These protocols generally require 14 synchronization to be repeated multiple times in order to obtain high precision. Furthermore, these WSN synchronization protocols are designed to achieve synchronization on the order of microseconds and as such are currently unsuitable for wireless coherent radar applications. In [45] 150 MHz bandwidth frequency modulated continuous wave (FMCW) signals are used to achieve synchronization of 66 ps. Two stations synchronize using the `detect' and `respond' method found commonly in TWTT-based approaches shown in Fig. 2.2. The authors propose a method for network synchronization in which slave nodes wait to detect the signal from the master and then send responses in TDMA fashion, however the authors only demonstrate the performance experimentally for two nodes at far above the expected Cramer-Rao Lower Bound (CRLB). 2.3.2 Ultra-wideband Pulse Synchronization A class of wireless synchronization methods of recent research interest are those based on ultra- wideband (UWB) pulse signals, rather than network packet exchanges. These UWB methods rely on wide bandwidth and high sample rate ADC hardware, often > 1 GHz, to estimate time of arrival (TOA) at nanosecond and sub-nanosecond levels. UWB signaling has been applied to numerous synchronization models, including passive schemes with multiple receivers synchronized to a single transmitter [46], decentralized consensus approaches [47], and distributed sensor localization [48] all using TOA estimated from high sample rate ADCs. Notably, in [49], a chip scale atomic clock (CSAC) and 64 GHz hardware clock packet timestamp counter is used with a propagation-aware time of ight (TOF) protocol to distribute timing and achieve < 5 ns pair-wise sensor synchronization. In [46], a passive method was developed for estimating receiver clock parameters and sensor location using many synchronization messages from a single master transmitter. However, the precision of the protocol is limited to the sampling rate of the receiver ADC and the resolution of a hardware time measurement device. Further a large number of synchronization iterations are required to obtain sub-nanosecond synchronization. In systems without highly stable oscillators, over the time needed to execute the number of synchronizations necessary for sub-nanosecond performance, the clock phase osets would drift. In [47], the Blink algorithm is proposed and higher bandwidth UWB signals and high sampling rate ADCs are used to improve synchronization to the nanosecond level using simple peak sample detection. 15 In [48], ad hoc sensor synchronization is discussed using time of arrival (TOA) measurements. This work adopts a linear model of the local clock and discusses TOA-based distributed localization with unknown internal delays and clock frequency osets in wireless sensor networks. In [50], the problem of clock frequency synchronization is discussed and localization algorithm is proposed. TOA uses a P/N sequence and correlator. Again, this method is limited by the sampling rate of the ADC used in the UWB system. In [49], the Pulsar algorithm uses a CSAC to obtain frequency stability and a propogation aware TOF protocol is developed to distribute timing. The Decawave DW1000 Radio used is able to timestamp packets to 15.6 ps precision \through equivalent time sampling of a pulse stream that is part of the message preamble." which is only possible due to a 64 GHz clock counter running on the hardware. However, a common misconception exists in UWB-based synchronization literature that the time resolution of TOA measurements is bounded by the ADC hardware sampling clock rate. In fact precise TOA estimation to within small fractions of a sample bin are possible, particularly in line of sight (LOS) environments, as dictated by the CRLB [51]. 2.3.3 Synchronization using SDR Research in SDR has seen a recent explosion due to their wide availability, low cost, and ease of development and prototyping for numerous applications. SDR hardware consists of an RF frontend and FPGA or FPGA+Processor system on a chip (SoC) and is capable of agile operation across wide frequency bands. Due to their software basis, SDR platforms oer unparalleled exibility in waveform generation and signal processing operations. Accordingly, wireless synchronization of SDR sensors has been a topic of recent research interest. In [52], wireless synchronization is proposed for SDRs. They address fractional clock phase estimation using a matched lter bank of 16 fractionally delayed Zudo-Chu sequences to estimate residual timing osets of 1=16 sample duration. Propagation delays are not accounted for and slaves are synchronized to a broadcast message sent by a single master. Processing is performed in software, signal bandwidth is limited to 1 MHz, and the achieved residual timing precision is 500 ns, far above the expected CRLB. 16 In [53], USRP E312 embedded SDRs are used to perform timestamp-free synchronization. They obtain< 1s synchronization (.8s with 125 kHz sampling rate; a precision of 1/10 the sample rate). Master waits to detect pulse from slave and schedules a time reversed transmission of the received waveform. Quadratic least-squares tting is used to estimate the peak of the autocorrelation signal to sub-sample resolution. The limited real-time correlation processing capability of the embedded USRP E312 is cited as the reason for choosing such a low 125 kHz sampling rate. In [54], a frequency and phase synchronization architecture for beamforming with commodity SDRs over wireless packet networks is presented. Using multiple transmitters and a single receiver, the authors divide synchronization into a frequency synchronization component and a beamsteering component using one bit feedback from a receiver. An extended Kalmann lter (EKF) is used to estimate frequency/phase oset from CW signals. All of these approaches require continuous real-time signal correlation to detect a pulse trans- mitted by a master sensor and rely on a software implementation of correlators. As such, they are limited by the streaming bandwidth from the FPGA to processor, and must use a digital down converter (DDC) to reduce the signal bandwidth from the full ADC hardware sampling rate so that continuous streaming across the FPGA-processor interface is possible. 2.4 Distributed Multistatic and MIMO Radar Multistatic wireless radar, wherein a single transmitter and multiple receivers operate in a coherent fashion, is most often studied in relation to satellite constellations. According to a 2006 survey and analysis of space-borne bistatic and multistatic radar systems, solving the time and phase synchronization problem is paramount [55]. The most common approaches to multistatic radar synchronization are ultra-high-quality oscillators, such as CSACs, direct path signal based methods, which require accurate positional knowledge a priori, and bidirectional link based methods [56]. The use of a dedicated synchronization link for phase synchronization and residual phase error compensation is explored theoretically in [57] for the bistatic radar case (a subset of multistatic radar with a single transmitter and a single, physically separated receiver). The work in [33] uses the direct path signal to perform time and phase synchronization upon signal reception. However, it relies on Global Positioning System (GPS)/inertial navigation system (INS)/inertial measurement unit (IMU) 17 data for positional knowledge and therefore its precision is limited to that of the peripheral sensors. Furthermore, it requires that the direct path signal is separable from the scattered signal of interest. This is usually impossible for ground-based or low-altitude UAV platforms, in particular when stepped frequency synthetic wideband techniques are necessary to obtain the requisite resolution performance. Finally. this work ignores the problem of spatial synchronization, and assumes the positions of all radar elements are known. multiple input multiple output (MIMO) radar is an extension of multistatic radar (single transmitter, multiple receivers), in which there are multiple transmitters and multiple receivers capable of operating coherently. A detailed overview and analysis of existing approaches to synchronizing MIMO radar elements is given in [58], where broadcast consensus algorithms are said to be scalable and robust for coherent MIMO radar. In [59], the authors present an UWB MIMO radar and propose a synchronization method that enables precise target localization by swapping pulses between the TX/RX elements. However, the method requires a central processor to synchronize the signals in post-processing. Further, the method relies on the ultra-wide bandwidth (500 ps pulse width or 2 GHz bandwidth) custom hardware to obtain high resolution performance, only synchronizes the elements with respect to some target, and requires that the node positions are known and xed. 2.5 Background on Snow Observation Using Radar Systems The summary in this section will be published in [60]. In the context of radar sensors, the most direct measurement is the two-way travel time (TWT) of the radar signal as it propagates through a medium and is re ected back to the receiver, at which point the magnitude and phase of the scattered signal is measured. Snow depth d s is related to the TWT of a radar signal through the snow as d s = TWT 2 c p s (2.6) where s is the relative dielectric permittivity of the snow and c = 3 10 8 m/s is the speed of light in a vacuum. 18 For dry snow, the permittivity is a direct function of snow density, s (gcm 3 ), and the known dielectric of pure ice (typically 3:15) [61{63]. Values for dry snow dielectric permittivity range from 1:0 1:9 as ice volume fraction, an indication of snow density, changes from 0 0:5 [64]. Snow water equivalent (SWE) is dened for a snow layer of height d s as [62] SWE = Z ds 0 s dz [cm] (2.7) In order to capture the high spatial variability of snowpack characteristics on the order of 10 cm, high radar resolution and thus large bandwidth, on the order of 1 GHz, are required [65] [66]. Radar sensors capable of achieving such high resolution, particularly at the lower frequencies required to eectively penetrate snow and ice media as well as vegetation layers, are typically expensive, power-hungry, and heavy. This is a signicant challenge in the widespread deployment of high performance radar sensors on sUAS platforms. The rst report of a ground/snow penetrating radar with high enough bandwidth to resolve snow stratigraphy and low enough weight and power consumption to be own as a small UAV payload being successfully own in a campaign to image snow elds was given in [67] [68]. The authors were able to conduct a single ight \with no in situ measurements to validate the measured snow depth due to harsh weather conditions and dicult accessibility" [68]. They use an UWB radar sensor developed by the German company Ilmsens and ew at a height of 1 m above the snow surface. A limitation of the m-sequence waveform used in the Ilmsens radar its low unambiguous range of 5:75 m, which constrains the maximum UAV height due to the relatively small unambiguous target window [67] [68]. The instrument did not penetrate into the ice with the airborne test. With the in-situ test the instrument was are able to penetrate 14:5 cm and 170 cm at the Artic sea and Kattfjordeidet test sites, respectively [68]. Improvements to the system and results from an extensive measurement campaign in Svalbard, Norway, which were compared with in situ snow depth, are reported in [69]. The authors present results from ve sites consisting of 100 100 m grids and two 100 m transects and report a spatial correlation of :97 with in-situ measurements. The further authors address the range ambiguity of the m-sequence radar and resulting limitations on UAV height by using a laser altimeter and dening an \ambiguity window" of 5:75 m that must contain the totality of the target. However, because of antenna crosstalk, this approach requires the 19 UAV to y in designated altitude zones so that the ambiguous echoed signal does not wrap into the higher power crosstalk signal [69]. The constraint of UAV height to either a maximum distance of 5:75 m above the surface or within designated altitude zones both requires a priori knowledge of the approximate snow depth and is not feasible for highly varied or rough terrain; particularly in the presence of both spare and dense foliage or forests where the laser altimeter is obscured and/or erroneous. However, a lower ight altitude when feasible yields signicant improvements in SNR and the tradeo of unambiguous range limitations for very high spatial sampling pulse repetition frequency (PRF) in the m-sequence radar is justied for suitable environments. Measurements of snow depth over Antartic sea ice from a UAV platform and comparison with manual snow probe depths have been reported recently in [70]. The authors use an array of Vivaldi antennas and vector network analyzer (VNA) to perform ultra-wideband stepped frequency continuous wave (SFCW) measurements of snow depth from a UAV. The Center for Remote Sensing of Ice Sheets (CReSIS) has developed numerous radar instruments for ice and snow sounding. Flexible ultra-wideband FMCW radars operating over selectable bands ranging from 14 MHz to 38 GHz developed by CReSIS have been demonstrated in vertical ice column prole imaging [71] [72] [73]. While the majority of current CReSIS instruments operate from traditional airborne platforms, the next generation CreSIS instruments will be targeted for small multicopters [71]. In [74], a CReSIS radar sensor deployed on a xed wing UAS is used to sound polar ice sheets; the rst successful report of this kind, and [71] describes eorts to integrate the CReSIS UWB Snow Mini Radar onto a small UAS helicopter. As of now, there has been no comparison of multicopter sUAS-based radar snow depth retrievals with ground-based GPR data reported in literature and no demonstration of sUAS radar systems in forested or topographically complex areas. The use of GPR has become widespread in cryospheric research as a means of replacing the laborious and intensive process of digging snow pits and collecting manual snow probe data [75]. Previous literature has proven that GPR surveys can provide accurate and ecient measurements of snow depth as well as high resolution imagery of snow layer stratigraphy [76] [77] [78] [79] [63] [80] [81]. As part of a validation eort for NASA's SnowEx Mission, a large scale eld campaign in Grand Mesa, Colorado was conducted to measure snow depth with GPR [63]. There was good agreement between GPR snow depth retrievals and manual snow depth measurements (correlation coecient 20 =:89), thus reinforcing the perception of GPR as an ecient means of collecting eectively ground- truth snow depth data [75] [76] [80] [78]. Additionally, ground-based GPR results were compared with airborne lidar from the NASA Airborne Snow Observatory (ASO), and stereo photogrammetry from DigitalGlobe WorldView-3 satellite imagery, with reported correlation coecients of 0.9 and 0.7, respectively [63]. Space and airborne [82] [83] [84] [71] [72] [73] as well as tower and xed-line [85] [86] [87] radar systems, particularly higher frequency FMCW radars, have been demonstrated in snow depth measurement, snow stratigraphy imaging, and SWE retrieval. Often times higher frequency FMCW radars are selected to more easily achieve broad bandwidth (and therefore resolution) at the expense of penetration depth. SWE retrieval via inversion of a forward scattering model has been demonstrated for dual-polarized X and Ku band SnowSAR data with correlation coecient of :64 [84]. Snow depth measurements obtained by a 2 18 GHz FMCW radar own on a xed wing aircraft ying at 500 m altitude have been reported with a correlations of :88 with in-situ measurements [86]. While space-borne and high-altitude airborne remote sensing methods are capable of coverage over vast areas, the resulting data products do not have suciently high spatial resolution to capture the ner scale variability and dynamics that are critical to achieving a complete picture of cryospheric processes [82] [65]. Multi-angle optical imagery and lidar-based approaches to snow depth measurements from sUAS platforms have been widely reported [88] [89] [90] [65]. Because of their high operational frequency, optical and lidar sensors can achieve high resolution with very small payload size, making them ideal for deployment on sUAS. Structure from Motion-Multi-View Stereo (SfM-MVS) processing uses a multi-angle sequence of 2D images as well as precise UAV platform knowledge and a series of ground control points (GCPs) to create high resolution 3D surface models of structures and survey area topography. [88] [89]. The fundamental shortcoming of optical and visible spectrum approaches to snow depth measure- ment, however, is that such sensors only measure the snow/ground surface; thus repeat observations are required to detect changes in surface elevation. Furthermore, they are highly sensitive to esti- mated snow densities in order to convert the measured snow depth to SWE. Additionally, in areas 21 with dense foliage or with thick forest canopies, optical signals cannot penetrate to the snow and ground surfaces below to obtain accurate measurements, making sub-canopy snow depth retrievals a persistent challenge [65]. 2.6 Landmine Detection and UAV-based Radar The summary in this section has been published in [18] and will be extended in [91]. The detection and removal of unexploded ordnance (UXO) and landmines is a signicant international issue that poses extensive economic, political and technological challenges on a global scale [92]. The use of ground penetrating radar (GPR) has been explored extensively as a means of detecting landmines for more ecient clearance of mine elds left in civilian areas after cessation of previous con icts. In general, landmines can be divided into three categories: metallic, minimum- metal, and non-metal and may be found on, just beneath, or far below the surface. The vast majority of damage to civilians is caused by small anti-personnel (AP) mines, such as the minimum-metal VS-50, which has a diameter of 9 cm [93]. This poses signicant challenges to detection of land mines with GPR. A high bandwidth is required to obtain the resolution necessary to detect signatures from small targets. Further, target proximity to and rejection of strong surface returns creates additional diculties. As noted by Daniels, "for the case of a weak target adjacent to a strong target [...] there is no accepted denition of resolution" [94]. The majority of current GPR landmine detection systems fall into two categories: vehicle based, such as the CSES HMDS, and handheld, such as the AN/PSS-14 [93]. A handful of airborne SAR-based systems have also been applied to this task. Feasibility studies in previous literature reports that for an airborne SAR system to reliably detect AP mines, a range resolutions of 5 cm is necessary; this would be obtained with 3 GHz of instantaneous bandwidth coverage from 200 MHz to 3 GHz [93]. The Mineseeker System is a large airship-based side-looking UWB SAR that meets these requirements and has demonstrated the capability to detect minimum-metal AP land-mines [93]. The use of small rotor-based unmanned aerial vehicles (UAV)s and inexpensive sensors for AP landmine detection has the potential to signicantly reduce cost, eort, time, and risk associated with mine eld clearance. However, in many cases radar systems capable of such resolution performance, 22 in addition to being expensive, are too heavy and large to be carried as a small rotor-based UAV payload. This has led to interest in the use of small commercial software dened radio (SDR) platforms such as the Universal Software Radio Peripheral (USRP) platforms for such applications. Previously, a USRP SDR mounted on a UAV has been explored as a GPR for landmine detection in [95]. However, this work imposes the USRP sampling bandwidth of 56 MHz as a constraint on the signal bandwidth. This bandwidth is far below the widely accepted requirements for a system to reliably detect small landmines while rejecting false targets and surface clutter. 2.7 Half Space Back-Projection Focusing The summary in this section has been published in [18] and will be extended in [91]. The time-domain back-projection (BP) focusing algorithm is widely used in modern synthetic aperture radar (SAR) systems due to the exibility and precision with which phase and motion errors can be compensated and the high quality of images produced. While more computationally complex than other common SAR focusing algorithms, BP lends itself naturally to parallel processing and pre-computation [96] [97]. Modication of BP that account for refractive and dispersive eects encountered in ground-penetrating SAR have been explored previously for traditional side-looking SAR [98] [99] [100]. In this proposed eort, we formulate a similar correction of BP focusing specic to a low-altitude nadir-looking airborne altimetric ground-penetrating SAR. 2.8 Smart Sensor Networks Initiated in 2010, the SoilSCAPE project is an eort to provide in-situ sensor validation data for air and space-borne remote sensing missions, including the NASA Soil Moisture Active Passive (SMAP) and the Airborne Microwave Observatory of Subcanopy and Subsurface (AirMOSS) missions, using smart sensor network technologies and machine learning [101] [102] [103]. The foundation of the SoilSCAPE wireless network architecture and control is based in part on previous work on sequential decision making in decentralized systems done in [104], which introduces the concept of common information as the key to nding optimal decision strategies in decentralized networks 23 where there is a time-delay in shared information. The next phase of the SoilSCAPE project aims to extend the architecture to heterogeneous sensor networks, including UAV-based SDRadar sensing nodes being developed as part of the proposed eort. 2.9 Related Space-borne Missions MirrorSAR uses optical transponder link referred to as MirrorLink on the receive satellites to relay received echo signal back to the transmitter. The optical transponder analog mixing is done with the same local oscillator that is used for the radar signal, thus preserving signal phase [105]. The TanDEM-X mission uses exchange of synchronization pulses between the transmitting and receiving satellites to compensate for oscillator phase noise and frequency deviations on the ground. [106]. 24 Chapter 3 Software Dened Radar System Development 3.0.1 SDRadar Architecture The SDRadar system described in this work is based on the USRP class of commercialy available SDRs developed by National Instruments (NI) and its subsiduary Ettus Research. USRP SDRs are used throughout academia, government, and industry for developing and testing new and emerging radio-based technologies. USRP platforms sold by Ettus Research are open source, which has enabled us to modify both FPGA rmware and software drivers to implement the SDRadar design. (a) System block diagram (b) Implemented system Figure 3.1: SDRadar implemented in USRP E312 hardware The SDRadar design targets embedded FPGA+ARM Processor SoC hardware, such as Xilinx Zynq series devices, and is comprised of two primary subsystems: the FPGA programmable logic (PL) custom cores and a C++ software server program that runs on the processing system (PS) as shown in Fig. 3.1a. The core components of the FPGA PL subsystem are an arbitrary waveform generator (AWG), a radar pulse controller (RPC), and a direct memory access (DMA)-rst in, rst 25 out (FIFO) buer. The AWG stores arbitrary sequences of complex baseband waveform samples, which are uploaded via software, to onboard block random access memory (BRAM). The RPC controls the timing and ow of transmit (TX)/receive (RX) data to/from the digital-to-analog converter (DAC)/ADC to maintain time-coherence between the TX/RX channels. The internal RPC state machine also implements a custom instruction set that congures and automates SDRadar operation, number of pulses, including pulse repetition interval (PRI), TX/RX window size, and channel selection, as well as demultiplexing data streamed directly from the PS software with radar data from the AWG, enabling simultaneous radio and radar functionality. The RX DMA-FIFO buers RX samples so that they may be transferred from the PL to the PS software [4]. The PL software subsystem is comprised of a driver that controls and communicates with the FPGA and SDRadar hardware, congurable signal processing stages, and an asynchronous Transmission Control Protocol (TCP)/Internet Protocol (IP) server that responds to and executes requests and collection tasks sent by client applications. Heterogeneous client devices with the appropriate SSL/TLS credentials and implementation of the custom request protocol can connect to the SDRadar server via Ethernet, WiFi, or serial radio (including LoRa and Digi XBee, which are reliable at distances > 1 km) and can send commands, initiate radar collections, and receive data. Through the C++ server software, all aspects of the SDRadar operation, including frequency selection, TX waveform, TX power, RX gain, etc can be congured in real-time. The client application has been implemented in C++, Python, and Java. The SDRadar server program is capable of queuing and executing complete complex stepped frequency mission plans autonomously. The SDRadar server is also responsible for streaming data from both internal (GPS, IMU) and optional external (real-time kinematic (RTK) GPS/global navigation satellite system (GNSS), lidar altimeter) sensors, which are stored as timestamped radar pulse metadata. Received complex I/Q radar pulse data are stored onboard in the SDRadar internal SD Card and can optionally be sent directly to the client after the mission completion. Additionally, the SDRadar can broadcast radar data to a UDP port allowing subscribed clients to process and visualize the radar data in real-time [4,60]. 26 The SDRadar design has been implemented and tested on multiple USRP SDRs including the USRP E310, E312, X300, N300, and N310 platforms with minimal modication, and on non-USRP FPGAs including Xilinx Zynq UltraScale+ Multiprocessor System on a Chip (MPSoC) and RF System on a Chip (RFSoC)-based boards with modied backend designs that remain server-side compatible with the existing client protocol. 3.0.2 Hardware Implementation Table 3.1: USRP E312 Specications Parameter Value Sampling Rate, f s 50 MHz Analog Bandwidth, B s 56 MHz Tunable Center Frequency, f c 70 MHz - 6 GHz TX Gain Range 0-89.5 dB RX Gain Range 0-76 dB Size 133 x 68.2 x 31.8 mm Price $ 3199.00 For this work, we use the USRP E312, a battery-powered 2x2 MIMO SDR module shown in Fig. 3.1b. The E312 has a Xilinx Zynq 7020 SoC, which consists of a 7 Series FPGA PL subsystem and PS with an ARM Cortex A9 866 MHz dual-core processor running a linux kernel. The E312 oers a maximum instantaneous bandwidth of 54 MHz over a tunable center frequency range of 70 MHz to 6 GHz and 76 dB and 89.5 dB of programmable frontend gain for the RX and TX channels, respectively. The maximum transfer rate between the PL and PS is limited to 10 MS/s for 32 bit complex samples. We used a sample clock rate of 50 MHz for the FPGA PL. Specications for the USRP E312 are given in Table 3.1 [107]. The AD9361 is a 2x2 RF transceiver with integrated 12-bit ADCs and DACs. It has two independent LOs: one shared by the receive channels (RX-A and RX-B), and one shared by the transmit channels (TX-A and TX-B). The chip's RX subsystem performs DC oset correction, quadrature correction, and digital ltering. We use serial peripheral interface (SPI) control to congure the AD9361 chip appropriately for successful coherent SWW reconstruction. By modifying the register map of the AD9361 so that multi-chip synchronization (MCS) is enabled, the phase relationship between LO sharing channels remains constant while switching between them allowing for coherent calibration. We note that in general, multi-channel USRP devices can be congured such 27 that deterministic relationships exist among channels, allowing for coherent loopback calibration. In single channel SDRs, however, an external switching network, controllable via general purpose input/output (GPIO) or other input/output (IO) bus is necessary [4]. Further software modications to the USRP Hardware Driver (UHD) driver include simplication of the frequency tuning routine for faster local oscillator (LO) phase lock loop (PLL) locking, and fast-switching routines between calibration and data channels. Using the loopback channel, we pre-compute a frequency-to-gain look up table (LUT) in order to atten the transmit power and receiver gain across wide stepped frequency bands. We add wrappers to the custom FPGA PL cores to make them compatible with the RF Network on Chip (RFNoC) infrastructure developed by Ettus Research [108]. When operating the SDRadar in monostatic mode, in order to reduce internal electronic coupling and provide additional gain, we added an external amplier to the TX path. This increases the maximum output power of our SDRadar from 10 dBm, the maximum output power of the E312, to 20 dBm. We used one of the E312's USB ports to power this external amplier. We connected the second RX channel, RX-B, to the TX path through a power splitter and 20 dB attenuator as a phase-coherent signal reference. Fig. 3.1a illustrates the complete SDRadar implementation. In bistatic/multistatic/MIMO modes, the loopback path and external amplier is replaced with two synchronization antennas connected to the TX-B and RX-B channels, which are used to synchronize multiple SDRadar devices. The synchronization technique is described in detail in Chapter 5. We note that this conguration is shown in Fig. 3.1b. We operate the SDRadar in a pseudo-bistatic conguration with two tapered slot Vivaldi TSA600 antennas from RF-Space used for transmit and receive. The antennas have an operational frequency of 600 MHz to 6 GHz and nominal beamwidth of 45 . Extensive antenna testing and pattern measurements were performed by the Antenna Test Lab Company [109]. 3.0.3 MIMO Sensor Network Operation Both the SDRadar and client software are fundamentally designed to support both multi-device network architectures to enable wireless multistatic and MIMO radar operation as well as real-time data broadcast to multiple end-users. The asynchronous server-client architecture implemented in the system is such that a single client can control N SDRadars (or any subset thereof) and N clients can connect to and control a single SDRadar. Furthermore any client may subscribe to the UDP 28 Figure 3.2: Client GUI software controlling multiple SDRadar devices. A single client can control N SDRadars and N clients can control a single SDRadar device. N clients can also subscribe to UDP data stream broadcasts to for real-time data visualization and procesing. data stream of a given SDRadar for real-time data processing and visualization on the client side as illustrated in Fig. 3.2. Therefore, the system as a whole supports can support a exible array of network communication and control topologies, which can be tailored to t the needs of a given operational mode or application. The pulse per second (PPS) signal from each SDRadar's onboard GPS receiver is used to synchronize FPGA hardware clocks to within < 40 ns with respect to Coordinated Universal Time (UTC) allowing each device to operate according to a known time-division multiple access (TDMA) schedule with slot sizes of 100 s when multiple devices are present on the network. The onboard GPS is also used as a local NTP server that, in conjuction with an NTP software daemon, keeps each SDRadar's processor system clock synchronized to within < 1 ms UTC. This coarse level of synchronization, while not precise enough for coherent multistatic or MIMO radar applications, is sucient for large numbers of SDRadar sensors to operate cooperatively within the same network free of interference. 29 Chapter 4 Ultra-Wideband Synthesis for High-Range-Resolution Software Dened Radar In this chapter, we rst describe a frequency stacking algorithm and SWW reconstruction technique for achieving high resolution performance in limited instantaneous bandwidth SDR-based SDRadar systems [4]. We then explore nonlinear FM waveforms and extend the SWW waveform to a new NU-NLSWW waveform [22]. Finally, we propose a novel non-uniform frequency stitching (NUFS) algorithm [60] that is capable of reconstructing both linear SWWs and NU-NLSWWs with minimal grating-lobe contamination, which have previously been regarded as the chief limitation of such waveforms. In this chapter we expand upon work previously published in [4] and work that will be published in [60] which expands upon [22]. 4.1 Introduction Content from this section has been published in [4]. Radar sensors are indispensable in our ability to achieve scientic, commercial and military goals for detection, surveillance, and characterization of discrete or extended targets. Radar systems provide precise power and phase-sensitive distance measurements to scattering targets and media that are hidden from or otherwise inaccessible to optical and other higher-frequency electromagnetic regime sensors. Current and previous generations of radar systems are generally highly specialized, 30 having complex application-specic hardware designs that are nearly impossible to recongure [7]. As such, these radar sensors must be designed, tested and rened on a system-by-system basis, making them exceedingly expensive { and for some applications prohibitively so. Frequency-Stepped radar is commonly used as a means of obtaining high range resolution with limited instantaneous bandwidth [20]. However, stepping across a wide band in steps small enough to reconstruct the band from pure tones is inecient in USRP SDR hardware due to LO tuning times that are long relative to the ADC sampling and data transfer rates. It is most time-ecient to use all of the available instantaneous bandwidth at each step and minimize LO re-tunings. However, using wideband sub-pulse waveforms and increasing the frequency step size gives rise to undesirable grating lobes [110], [19]. We employ a frequency stacking technique based on on the wideband stepped-chirp reconstruction proposed in [27], [28] to reconstruct a SWW. In addition, we explore and employ practical methods for mitigating the widely documented grating-lobes that are characteristic of SWW waveforms. In this work, we demonstrate a tunable SDRadar implementation of SWW reconstruction that is capable of ultra-wideband, high-range resolution performance. Our hardware implementation allows the total bandwidth, pulse length, waveform, frequency bands, and SWW sidelobe levels, to be dynamically tuned so that the algorithm can be optimized to application-specic requirements either remotely or in real-time. With the proposed method and design, an ultra-wideband high-performance coherent SDRadar may be implemented eectively in any low-cost commercial SDR platform. This algorithm is waveform-independent, allowing for any pulse-compression waveform as well as dierent waveforms at each new center frequency. Further, the algorithm can reconstruct SWWs from any set of specied center tuning frequencies, allowing tunable spectrum usage. The market for RF communications hardware is changing rapidly. This method is readily implemented on next-generation SDR hardware that could yield performance improvements. This processing method doesn't require custom developed hardware and is a universal technique. Additionally the algorithm does not require coherent RF hardware, so it can also be used in truly incoherent systems, such as bistatic and multistatic non-colocated receivers. 31 4.2 Synthetic Wideband Waveform Reconstruction The results of this section have been published in [4]. In order to achieve ultra-wideband centimeter-level resolution performance, we have developed an algorithm called frequency stacking, which synthesizes an ultra-wideband radar signal from low instantaneous stepped frequency sub-pulses. The frequency stacking addresses the existing limitations of SWWs, namely the requirement for a large number of frequency steps and signal contamination by inpulse-like grating lobes. This development makes it practical to use commercial SDR boards for high performance radar applications and is the rst demonstration of cm-level resolution performance in a SDRadar implemented in USRP hardware [111] [4]. We obtain high range resolution performance by coherently combining a series of sub-pulses to form a wideband SWW. A common method for producing a SWW, often called stepped-chirp or stepped-waveform, aims to reconstruct a wideband LFM chirp from a set of narrowband LFM chirp sub-pulses, which is then pulse compressed to achieve high range resolution [27], [28], [29]. In Algorithm 4.1, we express the discrete-time SWW reconstruction from N sampled baseband sub-pulse signalsz n [t m ] andv n [t m ], obtained from the RX antenna and coherent reference calibration channels respectively for TX reference waveform w[t m ] each having instantaneous bandwidth B i . The baseband sampling rate and sub-pulse length are denoted as f s and T p respectively. We introduce L =f s T p as the number of sub-pulse samples and f n as the center mixing frequency for sub-pulse n. z n [t m ] =w[t m 2R s c ]e j 2fn 2Rs c e je;n (4.1) v n [t m ] =w[t m ]e je;n (4.2) 32 Due to the generally non-coherent TX and RX LO phase relationship in many commercial radio platforms, including USRP SDRs, we consider terms for the TX phase t;n and RX phase r;n errors present in each sub-pulse. We dene the total contribution of these random phase errors after up/down-conversion as e;n = t;n r;n . This method performs pulse-compression by applying a matched lter to each baseband sub- pulse prior to reconstruction. This has the dual benet of both eliminating the need to perform compression on a massively up-sampled signal, thus freeing computational resources, and allowing parallel processing of the baseband sub-pulses in real time. Under the aforementioned assumptions of a stationary platform and no inter-sub-pulse doppler shifts, the frequency stacking method achieves signicant speed improvements by focusing on compressed SWW reconstruction from the phase coherent matched lter outputs of sub-pulses, making real-time implementation feasible in commercial SDR platforms. Furthermore, this approach is independent of the properties of the LFM chirp and can be applied equivalently to produce a SWW with bandwidth r = c 2Nfc for any selection of arbitrary sub-pulse waveforms. Although the results in this section focuses on LFM chirp waveforms, we explore other waveforms in later sections. The mathematical formulation of the frequency stacking algorithm can be found in Appendix A. 33 Algorithm 4.1 Frequency Stacking 1: Forn = 0;:::;N1, compressz n [t m ] withv n [t m ] at baseband to obtain the series of compressed sub-pulses d n [t m ] =z n [t m ]~v n [t m ]. Take the discrete Fourier transform (DFT) of each and compute the compressed spectra D n [f k ] =Z n [f k ]V n [f k ] 2: Filter the baseband pulse in the frequency domain to bandwidth B s = f c by element-wise multiplication with rect[ f k Bs ]. (Note that in order to reduce spectral notches, we require that B s B i .) 3: Append/pre-append zeros in frequency domain symmetrically about DC to upsample suciently to L up NL. 4: Frequency shift each sub-pulse to f n 5: Sum the n frequency-shifted compressed sub-pulses. D[f k ] = N1 X n=0 D n [f k f n ] rect h f k f n B s i (4.3) 6: Perform an IFFT to obtain the compressed SWW d[t m ] d[t m ] = 1 L up Lup1 X k=0 D[f k ]e j 2km=Lup (4.4) Because SDR hardware cannot in general be assumed to be phase-coherent, care must be taken to ensure that a phase-coherent reference signal v n [t m ] is obtained, either through external loopback or via some other channel. Due to random phase errors caused by ambiguous phase-locked loop (PLL) divider states in many RF transceivers used for communications, the reference v n [t m ] must be accurately measured such that it is phase-coherent with z n [t m ]. Depending on the SNR ofv n [t m ], it may be preferable in step 1 of Algorithm 4.1 to use this reference only to estimate and remove the random phase error of each sub-pulse, via conjugate-phase product, and perform compression directly with the noiseless discrete-time reference waveform w[t m ]. 34 4.2.1 Grating Lobe Supression The Fresnel ripples that occur at the edges of the LFM Chirp's and other waveform's spectra produce well documented grating-lobes that appear at intervals of c 2fc in the pulse-compressed SWW [110] [30]. Several techniques exist for mitigating grating-lobes that result from the stepped- waveform SWW reconstruction including overlapping sub-pulses, use of window functions and spectral weighting techniques [110], [30]. Increasing sub-pulse overlap (or equivalently decreasing the frequency spacing f c and bandpass lter bandwidth B s ) can reduce grating lobes but at the cost of sub-pulse bandwidth eciency (%BW eff = (B s =B i ) 100) as well as a loss of total bandwidth BW t when N and B i remain xed. Further improvements are obtained by windowing each sub-pulse with a tapered cosine or Tukey window, which reduces far out grating lobe levels. The exibility of SDRadar allows for multiple methods of grating lobe suppression to be combined arbitrarily in a manner that maximizes performance. In this work we limit our scope to focus on the following four techniques: (i) overlapping sub-pulses, (ii) sub-pulse window functions, (iii) non-uniform frequency spacing and (iv) grating lobe suppression with an inversion lter. These methods are chosen because they have minimal eects on the total time needed to complete an entire frequency sweep. The characteristics of grating lobes and eectiveness of these dierent techniques in mitigating them are further explored in Section 4.2.2. 4.2.1.1 Non-uniform Frequency Spacing The impulse-like shape of the grating-lobes is due to periodicity in the SWW spectrum caused by uniform spacing of sub-pulses in the frequency domain by a xed value f c . It has been shown that by spacing sub-pulses at non-uniform frequencies, the energy in the grating-lobes is spread out in the resulting non-uniform synthetic wideband waveform (NU-SWW) [110] [30] [31]. We dene a non-uniform frequency spacing for each sub-pulse f c;n (u) for a uniform random variable U n (u) dened over a sample space with points u such that f c;n (u) = f c +U n (u) (4.5) 35 Although dicult to implement in traditional radar systems [112], this technique is readily imple- mented in SDRadar. 4.2.1.2 Grating Lobe Suppression Filter We repeat the frequency stacking algorithm to construct a reference SWW, W [f k ], from either an RX calibration channel data or transmit waveform samples. From this, we obtain a numerically calculated grating-lobe suppression (GLS) lter that can be applied to the SWW to shape the reconstructed spectrum and correct for the hardware transfer function. LetM[f k ] represent the frequency response of some ideal wideband compressed signal (nominally we will use the spectrum of a compressed wideband LFM chirp) with spectral support across the reconstructed total bandwidth BW t =NB s . We dene a lter H z [f k ] as H z [f k ] = 8 > > < > > : M[f k ] W [f k ] rect[ f k BWt ]; W [f k ]6= 0 0; otherwise (4.6) Normalizing for unit lter gain, we obtain our grating-lobe suppression lter H z [f k ] H z [f k ] = H z [f k ] H z [0] (4.7) Apply the inversion lter to D[f k ] to obtain the corrected wideband synthesized pulse spectrum D[f k ] = H z [f k ]D[f k ] (4.8) and take the IFFT to obtain the corrected synthesized compressed pulse d[t m ] 36 4.2.2 Simulation In order to characterize the SWW signal structure, we simulate SWW reconstruction for two cases: (i) xed total SWW bandwidth and (ii) xed sub-pulse bandwidth. In the rst set of simulations, the total bandwidthBW t remains xed at a constant value of 80 MHz and the sub-pulse bandwidths B i and frequency spacing f c are changed as a function of the number of sub-pulses N. This is illustrated in Fig. 4.1 for two values of N. (a) SWW stacked sub-pulse spectra: N=4, Bi = 20 MHz, fc = 20 MHz (b) SWW spectrum: N=4, Bi = 20 MHz, fc = 20 MHz (c) SWW Autocorrelation: N=4, Bi = 20 MHz, fc = 20 MHz (d) SWW stacked sub-pulse spectra: N=20, Bi = 4 MHz, fc = 4 MHz (e) SWW spectrum: N=20, Bi = 4 MHz, fc = 4 MHz (f) SWW Autocorrelation: N=20,Bi = 4 MHz, fc = 4 MHz Figure 4.1: Simulated SWW reconstruction for xed total bandwidth (BW t = 80 MHz). 100 % BW Eciency (B i = B s ). The simulation parameters are chosen to illustrate the spreading of grating-lobes as the number of sub-pulses, N, increases when the total bandwidth, BW t , is constant. The impulse-like grating-lobes that are characteristic of SWWs and spread as a function of of f c are present in the SWW autocorrelation function and are demonstrated in Figures 4.1c and 4.1f. The severity of the SWW signal contamination due to grating-lobes increases with the number of sub-pulses used in the SWW. Hence, increases in the total bandwidth and improvements in range resolution performance gained by frequency stacking come at the cost of high grating-lobes. 37 (a) Non-Uniform SWW (b) Uniform SWW Figure 4.2: Fixed total bandwidth (BW t = 80 MHz) SWW autocorrelation as a function of number of sub-pulses, N = 1; 2;:::; 40 used in stacking. 100 % BW Eciency (B i = B s ). Non-uniform frequency spacing shown to eectively reduce grating-lobe levels. (See Fig. 4.1) (a) Non-Uniform SWW (b) Uniform SWW Figure 4.3: Fixed sub-pulse bandwidth and frequency spacing (B i = 10 MHz, f c = 10 MHz) SWW autocorrelation as a function of number of sub-pulses, N = 1; 2;:::; 40 used in stacking. 100 % BW Eciency (B i =B s ). Non-uniform frequency spacing shown to eectively reduce grating-lobe levels. As discussed in 4.2.1.1 using a non-uniform sub-pulse frequency spacing f c;n (u) = f c +U n (u) is used to mitigate grating lobes by reducing the periodic nature of their spectral occurrence. We use U n (u)2 [2; 2] % bandwidth to constrain the non-uniform frequency oset to within 2% of f n for any sub-pulse. In this case, with f n = 50 MHz, U n (u)2 [1; 1] MHz. We chose this range based on a 48 MHz spacing between 50 MHz bandwidth sub-pulses to reduce notches in the synthetic spectrum. Figures 4.2a and 4.3a demonstrate the eect of introducing a non-uniform frequency oset in simulation for the xed BW t and xed B i cases respectively as a function of N. 38 For comparison, the equivalent uniform SWW reconstructions are shown in Figures 4.2b and 4.3b for these two cases. It is clear that the introduction of non-uniformity in the SWW is eective at mitigating grating-lobe contamination. We note that GLS lter described in Section 4.2.1.2 is not explored in simulation as it guarantees exact reconstruction of the target wideband waveform for simulated data. The eectiveness of the GLS lter is however demonstrated experimentally in Section 4.5. 4.2.3 Loopback Verication In order to validate the performance of the proposed high range resolution SDRadar, we rst perform a test of the SDRadar in a loopback conguration with the RX-A data connected directly to the output of the split TX-A channel through a 20 dB attenuator. Each sub-pulse has an instantaneous bandwidth B i = 50 MHz (sample rate f s = 50 Msps). (a) Before lter (b) After lter Figure 4.4: Loopback Reconstructed SWW Spectrum before and after GLS Filter Correction. Corrected SWW Spectrum shown for spectral-shaping by two GLS lters: rectangular and hamming weighted. We apply a numerically calculated GLS lter as described in section 4.2.1.2 to perform spectrum shaping of the SWW. The lter is derived from the calibration channel SWW and therefore has the property of both spectral discontinuity removal as well as inversion of the digital subsystem transfer function common to both channels. The relationship between the two channels was conrmed to be deterministic through testing and channel swapping. The non-ltered SWW spectrum is shown in Fig. 4.4a. The result obtained by applying GLS lters to the loopback SWW is shown in Fig. 4.4b 39 for two dierent spectral shapes - one with a rectangular window and one with hamming window, demonstrating signicant attening of the SWW spectrum. The y-axis is normalized relative to the ideal full-scale value of an equivalent 16-bit signed integer signal. (a) Compressed SWW (b) 6dB-down main lobe width: theoretical vs. SWW Figure 4.5: Loopback SWW as BW t increases (700 MHz - 6 GHz in 48 MHz steps), N = 1;:::; 111, BW t = 50;:::; 5328 MHz. 96 % BW Eciency We next show the resolution performance across the entire 6 GHz tunable frontend frequency range of the E312. Fig. 4.5 demonstrates main-lobe width convergence as more sub-pulses are added to frequency stacked synthetic pulse. The 6dB-down main-lobe width of the synthetic wideband compressed pulse reconstructed from collected loopback data is shown in Fig. 4.5b. For comparison, the 6dB-down main-lobe width of a theoretical wideband Linear FM Chirp with bandwidthBW t =Nf c is also shown in Fig. 4.5b. For consistency and delity of comparison, the theoretical values for the ideal LFM chirp 6dB-down main-lobe width in Fig. 4.5b are obtained from a simulated wideband LFM chirp using the same algorithm as was used for the SWW. The range resolution capabilities of the presented design remain nearly identical to that of an ideal wideband chirp for a given bandwidth as more sub-pulses are added to the reconstruction. Thus, across the entire frontend tuning capability of the SDRadar, bandwidth may be increased or decreased arbitrarily without any loss in performance as compared to theoretical expectations. As expected, the grating-lobe level increases with the number of sub-pulses used. The proposed lter suppresses far-out grating-lobes to below the noise oor and signicantly reduces grating-lobe levels near the main-lobe, with the rst grating-lobe suppressed by 15 dB. With 96 % BW eciency, 40 (a) 96% BW Eciency (N = 111) SWW (b) 50% BW Eciency (N = 222) NU-SWW Figure 4.6: Loopback SWW and NU-SWW shown with and without Hamming GLS lter vs. Theoretical LFM Chirp with Hamming Window (700 MHz - 6 GHz), BW t = 5:328 GHz the application of this lter alone produces the results shown in Fig. 4.6a. For reference, the compressed pre-lter SWW and theoretical hamming-weighted linear FM chirp with equivalent bandwidth also are plotted in Fig. 4.6a. The peak envelope of each signal is taken to make comparisons more clear. Table 4.1: SWW, 96 % BW eciency, BW t = 5328 MHz Parameter GLS Filter Theoretical None Rect Hamming PGL n/a -21.3 dB -34.7 dB -34.9 dB PSL -13.3 dB -14.6 dB -14.1 dB -20.2 dB -6dB main lobe 3.3 cm 3.4 cm 3.4 cm 5.2 cm Table 4.2: SWW, 80 % BW Ecency, BW t = 5328 MHz Parameter GLS Filter Theoretical None Rect Hamming PGL n/a -25.6 dB -34.7 dB -34.8 dB PSL -13.3 dB -14.6 dB -14.2 dB -20.4 dB -6dB main lobe 3.3 cm 3.4 cm 3.4 cm 5.2 cm Table 4.3: NU-SWW, 80 % BW Eciency, BW t = 5328 MHz Parameter GLS Filter Theoretical None Rect Hamming PGL n/a -22.7 dB -38.0 dB -38.1 dB PSL -13.3 dB -14.6 dB -14.1 dB -20.6 dB -6dB main lobe 3.3 cm 3.4 cm 3.4 cm 5.2 cm 41 Table 4.4: SWW, 50 % BW Ecency, BW t = 5328 MHz Parameter GLS Filter Theoretical None Rect Hamming PGL n/a -30.6 dB -38.9 dB -38.9 dB PSL -13.3 dB -12.9 dB -14.2 dB -20.7 dB -6dB main lobe 3.3 cm 3.4 cm 3.4 cm 5.2 cm Table 4.5: NU-SWW, 50 % BW Ecency, BW t = 5328 MHz Parameter GLS Filter Theoretical None Rect Hamming PGL n/a -30.5 dB -47.6 dB -52.3 dB PSL -13.3 dB -14.7 dB -14.0 dB -20.6 dB -6dB main lobe 3.3 cm 3.4 cm 3.4 cm 5.2 cm Because of the exibility aorded by SDRadar, we implement NU-SWW reconstruction by simply providing a new frequency plan to the server application running on the SDRadar. Fig. 4.6b shows the improvement in grating lobe level obtained by reducing BW Eciency to 50 % and adding a non-uniform frequency step size. Grating lobe levels are signicantly reduced at the cost of bandwidth eciency, which translates to an increased sweep time as more sub-pulses must be collected. Again the envelope of a compressed linear FM chirp with a hamming window and identical time bandwidth product is plotted as a performance benchmark. Note that we have normalized all autocorrelation plots to a 0 dB reference for ease of comparison. Tables 4.1, 4.2,4.5, 4.4 and 4.5 summarize key performance benchmarks, namely, peak grating- lobe (PGL) level, peak side-lobe (PSL) level and -6dB main lobe width for ve dierent SDRadar operating modes and compare the performance of rect and hamming grating-lobe suppression lters. The results are obtained from the average of 100 pulses and are compared to a theoretical wideband LFM chirp having the some total bandwidth as the reconstructed SWW. For applications requiring grating-lobes levels50 dB and able to tolerate BW eciency of 50%, the proposed SDRadar can meet requirements with range resolution equivalent to that of an identical bandwidth theoretical wideband LFM chirp. 4.3 Non-Uniform Nonlinear Synthetic Wideband Waveforms This section is the subject of an upcoming manuscript [60] which expands upon [22]. 42 4.3.1 NLFM Waveform Design Using the design strategy for NLFM waveforms described in Section 2.2.2, a NLFM waveform may be constructed from an LFM waveform with arbitrary amplitude weighting. A derivation of the NLFM waveform generation technique using the PSP can be found in Appendix B. We provide a MATLAB code for generating constant-amplitude NLFM waveform from an LFM waveform that has been amplitude-apodized by an arbitrary weighting in listing 4.1 and available on Github [113]. Listing 4.1: nonlinearfm.m 1 function [x nl , ft ] = nonlinearfm(x, fs ,B) 2 % x : LFM waveform , fs : Sample Frequency (Hz) , B : Bandwidth (Hz) 3 N = numel(x); T=N/fs ; df = fs/N; dt = 1/fs ; 4 f = [−fs /2:df: fs/2−df ]; 5 t = −T/2:dt:(T/2−dt); 6 % Compute Spectrum 7 xfft = fftshift ( fft (x)); 8 xfft = xfft (abs(f)<=(B/2)); 9 % Compute Group Time Delay function 10 Tg = cumsum(abs( xfft ).^2) ; 11 % Solve Boundary Conditions 12 c1 = T/(Tg(end)−Tg(1)); 13 c2 = −T/2−c1 * Tg(1) ; 14 Tg = c1 * Tg+c2; 15 % Invert Tg using interpolation 16 ft = interp1(Tg, f(abs(f)<=(B/2)) , t); 17 % Integrate frequency function to obtain phase 18 phi = 2 * pi * cumsum(ft)/fs ; 19 phi = phi+pi/2−phi(1) ; 20 % Compute nonlinear FM waveform 21 x nl = exp(1i * phi); 4.3.2 NLFM Simulation Results We characterize the waveforms produced by the method described in Section 4.3.1 in simulation and compare the generated NLFM waveforms to the corresponding target LFM waveforms for selected common window functions. The following results are obtained using the MATLAB function made available on Github [113] to construct a NLFM waveform from a sampled version of a signal with a desired PSD with similar PSD and autocorrelation characteristics. We use an LFM chirp having bandwidthB, such that the LFM waveformx(t) for some amplitude weighting functiona(t) is given 43 by x(t) =A(t) exp (jKt 2 ) (4.9) where K =B=T is the chirp rate and T is the pulse length such that the time bandwidth product (TxBW) of the given waveform is TB. Fig. 4.7 illustrates the dierences between amplitude apodized LFM and phase-apodized NLFM waveforms in the time (Fig. 4.7a) and frequency (Fig. 4.7b) domains. We note that the total energy in the NLFM waveform is the same as that in the non-weighted LFM waveform, whereas the total energy of the amplitude weighted LFM waveform is reduced. -5 0 5 Time ( s) -1 -0.5 0 0.5 1 Real Component NLFM Weighted LFM (a) Time-domain waveforms -10 -5 0 5 10 Freq (MHz) 20 25 30 35 40 45 50 55 60 Amplitude (dB) Non-weighted LFM Weighted LFM NLFM (b) Power spectra -1500 -1000 -500 0 500 1000 Range (m) -80 -70 -60 -50 -40 -30 -20 -10 0 Amplitude (dB) -34.619 dB NLFM desired LFM -100 -50 0 50 100 -60 -40 -20 0 (c) Autocorrelation (TxBW = 200) -1500 -1000 -500 0 500 1000 1500 Range (m) -80 -70 -60 -50 -40 -30 -20 -10 0 Amplitude (dB) -42.1015 dB NLFM desired LFM -100 -50 0 50 100 -80 -60 -40 -20 0 (d) Autocorrelation (TxBW = 2000) Figure 4.7: Signal characteristics of an amplitude weighted LFM chirp and the equivalent constant- amplitude NLFM chirp. Increasing the signal TxBW is shown, in general, to produce improved results, wherein the autocorrelation function of the NLFM waveform approaches that of the desired LFM counterpart. However, the degree of this performance improvement varies signicantly for dierent window functions. It is useful to consider the frequency vs. time curves of NLFM waveforms as they provide a complete characterization of the signal. Kaiser-Bessel functions or Kaiser windows are of particular interest because their spectral shape and autocorrelation characteristics are tunable by a single coecient [114]. Many common window functions are similar to a Kaiser window with a particular value of. The NLFM frequency vs. time curves for selected common window functions and Kaiser windows for varying are plotted in Fig. 4.8a and Fig. 4.8b, respectively. 44 -0.5 0 0.5 normalized time -0.5 0 0.5 normalized frequency LFM hamming hann blackman-harris chebyshev (a) Frequency vs. time curve for common window functions -0.5 0 0.5 normalized time -0.5 0 0.5 normalized frequency 0 40 80 120 160 200 (b) Frequency vs. time curve for K Kaiser win- dows 10 2 10 3 10 4 10 5 10 6 TxBW (log scale) -90 -80 -70 -60 -50 -40 -30 PSL (dB) hamming hann blackman-harris chebyshev NLFM desired LFM (c) NLFM PSL vs. TxBW for common window func- tions 10 2 10 3 10 4 10 5 10 6 TxBW (log scale) -70 -60 -50 -40 -30 -20 PSL (dB) 20 30 40 50 60 70 80 90 100 NLFM desired LFM (d) NLFM PSL vs. TxBW K Kaiser windows Figure 4.8: Characteristics of NLFM waveforms constructed from equivalent weighted LFM wave- forms for a number of common window functions. We plot the peak sidelobe level (PSL) as a function of Time-Bandwidth Product (TxBW) for the NLFM phase apodized waveforms in Fig. 4.8c for selected common window functions and in Fig. 4.8d for Kaiser windows K with parameter values 2f20;:::; 100g. For comparison, the PSL of LFM amplitude apodized waveforms are shown for the equivalent window functions as dashed lines. We note that as the PSL characteristic of a given window function decreases, a higher TxBW is required for the equivalent NLFM waveform to achieve the same result. 4.3.3 Non-Uniform Nonlinear Synthetic Wideband Waveform Reconstruction We now extend the frequency stacking (FS) SWW reconstruction technique described in Section 4.2. It has been demonstrated that ultra-wideband waveforms may be coherently synthesized from instantaneous bandwidth-limited stepped-frequency sub-pulses to achieve centimeter-level resolution capability SDRadar using low-cost commercial SDR boards [4,111]. An additional advantage of employing stepped-frequency techniques in SDRadar is that spectrum usage is easily managed by conguring the operational center frequency selection and/or the digital baseband signal used at each frequency to notch out particular bands. The chief limitation of the SWW is the presence of signal-contaminating grating lobes, whose impulse-like shape is due to periodic discontinuities in the SWW spectrum. It has been shown that by spacing sub-pulses at non-uniform frequencies, the energy in the grating-lobes is spread out in the resulting NU-SWW [26] [30] [31]. We note that the NUFS method described in the following 45 section inherently introduces non-uniformity into the SWW spectrum, even when the underlying frequency step sizes are uniform. Other work has been done on active grating lobe cancellation and ltration [4]. In this section, we expand upon work previously reported in [22]. We extend the NU-SWW and NLFM waveform design [22] [23] [25] to a novel NU-NLSWW. The approach to the waveform design is as follows: 1. Consider an ultra-wideband frequency modulated (FM) waveform with a desired power spectral density (PSD) and autocorrelation function. 2. Compute time domain frequency function corresponding to this ultra-wideband FM waveform. 3. Splice the time domain frequency function at non-uniform intervals with osets from a nominal interval f chosen randomly from a uniform distribution to construct a set of limited-bandwidth FM sub-pulses. To reconstruct the NU-NLSWW waveform from narrowband sub-pulses, we propose an extension of the FS algorithm [4] summarized in Section 4.2 that introduces a cost function C n [f] of the form C n [f] = min f 0 :::;f L1 L1 X i=0 (k n;1 jjD n+1 [f i f n+1 ]D n [f i f n ]jj 2 2 (4.10) +k n;2 jj @D n+1 [f i f n+1 ] @f @D n [f i f n ] @f jj 2 2 ) s:t: ff n + f=2f i f n + f=2 +g to be minimized in order to nd an optimal reconstruction. Here is a constraint factor that controls the extent of the search, k n;1 and k n;2 are empirically selected weight factors, D n [f] is the Fourier transform of the n th sub-pulse matched lter output, and f n =f n (f N1 +f 0 )=2 is the baseband frequency oset of each sub-pulse. The cost function incorporates the magnitude and rst-order derivative of the sub-pulse spectra in order to maximize smoothness of adjacent sub-pulse spectra. The sub-pulses are then stitched together at the minimizing frequencies ~ f n . ~ f n = argmin f 0 :::;f L1 C n [f] (4.11) 46 Because the algorithm searches for non-uniform optimal stitching seams in the sub-pulses rather than naively stacking them, we refer to this approach as Non-Uniform frequency stitching (NUFS). The reconstructed NU-NLSWW D[f] then becomes D[f] = N1 X n=0 D n [f f n ]1 ~ fnf< ~ f n+1 (4.12) where 1 af Dev 2 Dev 2 -> Dev 1 TOF Sync 1 s = 7.85 ps Dev 1 Dev 2 TX RX TX RX Figure 5.1: Local time of ight measurements through cables between two devices at 1 GHz carrier frequency. 10 trials each consisting of 100 synchronization pulses taken at a PRI of .01 s are shown. Trials were done at 5 s intervals for 50 s. The synchronized TOF measurement obtained has a standard deviation of 7.8476 ps that i = 1 for all sensors and that i changes slowly over time according to a smooth bounded random process. That is, we assume that i can be treated as a constant for the duration of time required to complete a single synchronization exchange cycle for the entire network; henceforth referred to a synchronization epoch. Assume each sensor is a radio transceiver capable of transmitting and receiving arbitrary complex baseband waveforms having bandwidthB sampled at a sampling ratef s which are up/down-converted to a tunable RF carrier frequency f c using TX/RX local oscillators (LOs) derived from the same oscillator that produces the baseband digital sampling clocks. That is, each sensor i is capable of transmitting a waveform w i (t) from a baseband complex arbitrary waveform s i (t), which may be overheard and received by all sensors in the network as a wireless broadcast. Due to generally non-deterministic LO PLL divider states when tuning to f c , we dene RF carrier phase oset terms tx i and rx i due to signal mixing with the TX LO and RX LO, respectively. The waveform w i ( i ) 74 generated in the clock domain of sensor i w i ( i ) =s i ( i )e j 2fc( i ) e j tx i (5.3) is transmitted in the global clock domain as w i (t) w i (t) =s i (t + i )e j 2fc(t+ i ) e j tx i (5.4) The continuous time waveform s i ( i ) is related to the digital sequence s i [n] by the Whittaker- Shannon interpolation formula for digital-to-analog conversion s i ( i ) = 1 X n=1 s i [n] sinc(f s i n) (5.5) where sinc(x) = sin(x) x is the normalized sinc function. Further, assume that two nodes i and j are separated by line of sight distance R i;j =R j;i which is related to the RF signal TOF by TOF i;j = R i;j c . The signal transmitted by sensor j is then received and down-converted by sensor i, with the RX LO derived from its own clock time i as r i;j (t). r i;j (t) =w j (t R i;j c )e j 2fc(t+ i ) e rx i (5.6) Note that we represent the imaginary unit using the roman character j to distinguish it from the sensor index j. Dene err i;j = tx j rx i as the cumulative phase oset error due to the RX and TX LOs of sensors i and j, respectively. We expand (5.6) to obtain the expression r i;j (t) =s j (t + j R i;j c )e j 2fc( i + j R i;j c )) e j err i;j (5.7) When expressed purely in the clock domain of sensor i, (5.7) becomes r i;j ( i ) =s j i ( i j + R i;j c ) e j 2fc( R i;j c + i j ) e j err i;j (5.8) 75 After digitization with analog-to-digital converters (ADCs) having sample rate f s , the discretized signal received by sensor i is r i;j [n]: r i;j [n] =s j nf s ( i j + R i;j c ) e j 2fc( R i;j c + i j ) e j err i;j (5.9) We consider the continuous time signal d(t) dened as the cross-correlation of r(t) with s(t) d(t) =r(t)s (t), (r?s)(t) (5.10) D(f) =R(f)S (f) (5.11) where the operator indicates convolution. S(f), R(f), and D(f) are the Fourier transforms (FTs) of s(t), r(t), and d(t) respectively, dened for D(f) as D(f) = Z 1 1 d(t)e j 2tf dt (5.12) For d i;j (t), we dene d i;j ( i ) =r i;j ( i )s j ( i ) (5.13) = (s j ?s j ) i ( i j + R i;j c ) (e j 2fc( R i;j c + i j ) e j err i;j ) (5.14) where ? denotes cross-correlation as dened in (5.10). The magnitude of the cross-correlation signal d i;j ( i ) has a global maximum at: t pk:i;j = argmax i jd i;j ( i )j (5.15) = i j +TOF i;j (5.16) at which point the phase is \d i;j (t pk:i;j ) = err i;j 2f c ( i j +TOF i;j ) (5.17) 76 We arrive at a similar expression for the discrete time cross-correlation d i;j [n] of the sampled signals r i;j [n] and s j [n], where the fractional peak index n pk:i;j is related to the peak time as. n pk:i;j =f s t pk:i;j (5.18) There is a common misconception that the time resolution of the cross-correlation peak is limited by the sampling rate [51], and that the sampling clock rate is the limiting factor in how precisely signal arrival time can be measured [47]. However, this is not the case as discrete time signals sampled according to the Nyquist sampling rate contain all information present in their continuous time counterparts. That is to say that estimation of discrete autocorrelation peak location t pk:i;j for a single target is possible at timing resolution that is orders of magnitude beyond the width of a sample clock bin. This is explored Section 5.5. 5.3 Syntonization and Coarse Synchronization Prior to performing synchronization via the proposed method, all device clocks are synchronized in frequency or syntonized. This can be achieved with a phase locked loop (PLL) that is locked to a pulse per second (PPS) signal. Because the device clock will oscillate many times during the period of a single PPS signal, a PLL, digital control loop, and voltage controlled crystal oscillator (VCXO) may be used to create highly stable frequency reference clock from a PPS signal. This type of coarse synchronization of timestamps is achieved to within a few clock cycles by rising edge detection of the same PPS reference. In this work we assume the PPS for coarse synchronization can be obtained from a GPS receiver in a GPS-enabled environment. Generally, the accuracy of GPS PPS signals for commodity GPS receivers is on the order of 10 ns with respect to Coordinated Universal Time (UTC). Thus, coarse synchronization of timestamps to 10 ns is possible by edge detection of a GPS PPS reference. In order for the method described herein to work, the coarse synchronization must be (i) accurate enough such that the master controller, which instructs the sensor network to commence operation, is time synchronized to the sensor network to within 1 s of UTC and (ii) precise enough for the receive windows of all sensors in the network to be suciently aligned such that they are guaranteed 77 to contain only the signal sent by the transmitting sensor(s) as determined by the scheduling and orthogonality schemes of the network (i.e., within 100 s for a TDMA time slot allocation of the same size). Extensive work has been done to explore synchronization to a GPS level of precision in GPS-denied environments [47{50,125]. For GPS-denied environments, syntonization and coarse synchronization must be derived from another source. One solution is a master reference source that distributes a PPS signal as a P/N sequence. Each receiver has a 1-bit correlator that detects this master signal and generates a local PPS from which a timing reference can be derived and used for syntonization and coarse synchronization, similar to [50]. A `propagation-aware' approach, such as [49], is necessary for accurate coarse time synchronization. In general, coarse synchronization is best achieved by an independent system, of which many exist, and as such, further exploration of coarse synchronization in GPS-denied environments is not investigated here. 5.4 Fine Synchronization Based on the previously stated assumptions of syntonization and coarse synchronization, we derive a decentralized method for ne synchronization of wireless sensor nodes in an N-sensor network to nearly 1=1000 th of the sampling clock rate. Assuming that coarse clock synchronization ( 10 ns) and frequency syntonization can be achieved using existing methods such as GPS, this work proceeds based on two key insights: 1. If two radios with random time-varying relative clock osets transmit to one another simulta- neously (according to their own clock), the local delay measured in each is symmetric about the true delay. The average of the two local delay measurements is the minimum variance estimate of the true delay. This approach does not require hardware based signal detection and response in deterministic time to measure TOF delay. This has been shown to apply for integer clock cycles, and holds for fractional clock cycles (clock phase) as well. 2. If sampling clocks are characterized suciently, modifying the baseband digital waveform prior to transmission can result in a transmitted signal that is equivalent to that produced by coherent synchronized clocks. This means that in pulsed applications, e.g., radar, eective coherent synchronization may be achieved without synchronization of hardware clocks. 78 We use a time-division multiple access (TDMA) scheme for this work; however, other orthogo- nality schemes such as code division multiple access (CDMA) or combined TDMA/CDMA can also be used. For a given TDMA TX slot, because the TDMA schedule is known to all sensors in the network, the local integer clock timestamp of any transmitting sensor when a given pulse is sent will be known across the network. 5.4.1 N-Sensor Synchronization Exchange Assume that a wireless network of N sensors have syntonized oscillators and that their timestamps are coarsely synchronized to within a single integer clock cycle of a global reference, and one another. The remaining clock oset for sensor i relative to some global reference is due to a clock phase term i , which represents a fractional clock oset. We dene the synchronization epoch as the period during which a single iteration of the synchronization algorithm is performed as shown in Fig. 5.2 for the two sensor case. We assume that i is constant over the two way synchronization epoch, but varies randomly over larger time scales. Furthermore, we assume that the relationship between the local clock of sensor i, i and some globally `true' time t may be expressed as t = i i . ⌧ proc Radio 1 Radio 2 ⌧ rx 1,2 =TOF+ TDMA + 1 2 ⌧ rx 2,1 =TOF+ 2 1 ⌧ rx 0 1,2 = ⌧ proc +TOF+ TDMA + 1 2 ⌧ rx 0 2,1 = ⌧ proc +TOF+ 2 1 ⌧ 1 = t+ 1 ⌧ 2 = t+ 2 t tx 1 = ⌧ tx 1 1 ⌧ tx 2 = ⌧ tx 1 + TDMA ⌧ proc + ⌧ tx 1 ⌧ proc + ⌧ tx 2 ⌧ 1 ⌧ 2 Encode :{m 1,2 } Encode :{m 2,1 } Decode :{m 2,1 } Decode :{m 1,2 } Both radios have shared knowledge of {m 1,2 } and {m 2,1 } m i,j : message from radio i to radio j TDMA Figure 5.2: Synchronization timing diagram for a single synchronization epoch 79 We denote TOF i;j as the free-space signal time of ight between two sensors i and j due to a physical separation of TOF i;j c where c is the speed of light as described in Section 5.2. Transmissions are performed according to a known TDMA scheme starting at time tx 1 with time slot size TDMA where each sensor transmits at time tx i = tx 1 + (i 1) TDMA according to its local clock. We dene rx i;j as the time that sensor i receives the signal transmitted by sensor j relative to its own clock rx i;j =TOF i;j + (j 1) TDMA + i j (5.19) This is shown for the two sensor case in Fig. 5.2 as the rst signal exchange. For illustrative purposes, TDMA is depicted as being < TOF in this diagram. However, in practice TDMA should be chosen to be >TOF for wireless channels so as to avoid signal interference. We note, however, that from the wireless sensors point of view, these transmissions can be considered to occur simultaneously. We denem i;j as the encoded message containing sensori's measurement of rx i;j after subtraction of the known TDMA time slot oset m i;j =TOF i;j + i j (5.20) A second round of signal transmissions are performed after a known processing time oset proc according to the same TDMA scheme used for the rst exchange. Sensors transmit a second pulse with all messagesfm i;k jk6= ig appended at local time tx 0 i = proc + tx 1 + (i 1) TDMA . We dene rx 0 i;j as the time relative to its own clock that sensor i receives a second signal transmitted by sensor j containing both the synchronization waveform and the set of messagesfm j;k jk6=jg which includes message m j;i . rx 0 i;j = proc +TOF i;j + (j 1) TDMA + i j (5.21) This is shown as the second signal exchange in Fig. 5.2 for the two sensor case. 80 Each of theN sensors now estimates theNN matrix of synchronized time of ight measurements from the shared messages as ~ TOF, which is symmetric with zeroes along the diagonal. ~ TOF i;j = m i;j +m j;i 2 (5.22) Similarly, the relative estimated clock phase osets between sensor i and j are represented as an NN matrix ~ . ~ i;j = m i;j m j;i 2 (5.23) = i j (5.24) We note that equations 5.19-5.22, if reduced to the N = 2 case, are similar to the two sensor expressions given in [124] and exploit the same fundamental concepts of TWTT [36]. Dev 1 Dev 2 Dev 3 m 2,1 m 3,1 (a) Stage 1:TDMA slot 1 Dev 1 Dev 2 Dev 3 m 1,2 m 3,2 (b) Stage 1:TDMA slot 2 Dev 1 Dev 2 Dev 3 m 2,3 m 1,3 (c) Stage 1:TDMA slot 3 Dev 1 Dev 2 Dev 3 Encode: {m 2,1 , m 3,1 } (d) Stage 2:TDMA slot 1 Dev 1 Dev 2 Dev 3 Encode: {m 1,2 , m 3,2 } (e) Stage 2:TDMA slot 2 Dev 1 Dev 2 Dev 3 Encode: {m 1,3 , m 2,3 } (f) Stage 3:TDMA slot 3 Figure 5.3: Illustration of synchronization scheme for three sensors The synchronization process is illustrated for three sensors in Fig. 5.3 for the rst (Figs. 5.3a-5.3c) and second (Figs. 5.3d-5.3f) signal exchange stages. 81 5.4.2 Clock Compensation and Transmit Synchronization Following the completion of the synchronization epoch described in Section 5.4.1, each sensor now transmits a waveform, again using an arbitrary orthogonality scheme. When a signal crosses to/from the clock domain of a given sensor, the time error of the sampling clock and phase error of the mixing LO are imprinted on the signal. From the estimates of the matrices ~ TOF and ~ , which are now known identically across the network, we derive terms for each sensor that estimate and correct its relative clock and RF carrier phase errors prior to transmission, synchronizing all signals in the air; and then again upon reception, synchronizing all signals in each sensor's respective local clock domain. We now choose to align the transmit waveforms to the mean clock phase oset in the network. The fractional delay shift applied to the baseband waveform of sensor i prior to transmission is ~ i ~ i = 1 N X j ~ i;j (5.25) = i 1 N X j j (5.26) We note that transmit waveforms could also be phase aligned to a selected reference sensors rather than the network mean. For example, we may synchronize to sensor 1 by replacing the expression in (5.25) with ~ i = ~ i;1 . For transmit coherence, the waveform is fractionally delayed by ~ i such that the updated local transmission time, denoted as ~ tx i for sensor i is ~ tx i = tx i + ~ i (5.27) Due to up-mixing with the LO generated from the now characterized baseband clock, the sample clock-dependent RF phase error (as described in (5.4)) is pre-compensated by applying a carrier phase correction term of e j 2fc ~ i (5.28) 82 to the now time-delayed waveform prior to transmit. After shifting to the new local transmission times, we can write the true transmission time of each waveform sent by sensor i as t tx i = ~ tx i i (5.29) = tx i 1 N X j j (5.30) Thus all eects of each individual sample clock are removed and signals are transmitted synchronously and coherently in the air. Now the signal sent by sensor j arrives at sensor i at true time t rx i;j t rx i;j =t tx j +TOF i;j (5.31) and is received by sensor i at time rx i;j as measured by its own clock rx i;j =t rx i;j + i (5.32) = tx j 1 N X j 0 j 0 +TOF i;j + i (5.33) Noting that ~ i = i 1 N P j 0 j 0, this becomes rx i;j = tx j +TOF i;j + ~ i (5.34) By adding a fractional delay to the received waveform of ~ i , each sensor will also now coherently receive the coherently transmitted waveforms relative to their own clocks at time ~ rx i;j . ~ rx i;j = rx i;j ~ i (5.35) = tx j +TOF i;j (5.36) 83 Again, due to down-mixing with the LO generated from the receiving baseband clock, the sample clock-dependent RF phase error (as described in (5.6)) is corrected by applying a carrier phase correction term of e j 2fc ~ i (5.37) to all waveforms received by sensor i. The eects of each sensor's random clock phase are removed in both the global time domain upon transmission as well as in each sensor's local clock domain upon reception. Thus both transmitted and received signals are synchronized in time for all sensors in the wireless network. Furthermore, the sample clock-dependent LO RF phases are removed so that the phase relationships between all sensors are stable. Thus, by applying the described clock and carrier phase corrections the signals transmitted (and received) by the network of sensors become indistinguishable from those that would have resulted if all sensor clocks were physically connected and locked to a distributed common reference, and we may therefore consider the network synchronized. 5.4.3 RF Carrier Phase Synchronization Due the generally non-coherent carrier phase relationships between the RF frontends of multiple disparate radios, achieving true coherent operation requires that an additional constant non-clock dependent carrier phase oset of each radio in the network be estimated and corrected. For a signal sent by radio j to radio i, we denote the residual carrier phase oset error as err i;j . We assume that this phase can be modeled as a linear combination of the the transmitting local oscillator (LO) phase, tx j and receiving LO phase tx i , so that err i;j = tx j rx i (5.38) We note that in general for SDR boards, tx i and rx i cannot be assumed to be equal. 84 Using the synchronized time of ights in the network computed from (5.22), the relative clock osets computed from (5.23), and based on the signal model in (5.17), we can estimate the residual uncompensated RF carrier phase error err i;j as err i;j =\d i;j [n pk ] + 2f c ( ~ TOF i;j + ~ i;j ) (5.39) where\() denotes the phase of a complex number. Because there are N 2 measurements of err i;j with only 2N unknowns tx i and rx i , we can nd solutions to the transmit and receive carrier phase errors forN 2. If each sensor is unable to perform a measurement of its own total LO phase oset, i.e., by receiving its own transmitted signal, then the number of equations reduces to N 2 N. Thus carrier phase error estimation is still possible if N 3. In either case, for the resulting linear system matrix to be full column rank, the transmit channel of one sensor must be chosen as the carrier phase zero-reference, i.e., tx 1 = 0. After solving for the residual transmit and receive phase errors, each sensor i may apply a conjugated phase correction term of tx i to its transmit waveform and of rx i to all received waveforms. The regularity with which carrier phase synchronization must be performed depends on both the transceiver characteristics and operating mode. For example, in the AD9361 RFIC used for this research, and in many commercial SDR boards, the LO phase oset changes non-deterministically each time the LO is retuned. Thus for frequency hopping applications [4], carrier phase resyn- chronization is required after each frequency change. However for single frequency applications, performing carrier phase synchronization once as an initial calibration may be sucient. 5.5 Peak Detection The performance of the proposed synchronization algorithm depends on the precision with which time delay can be estimated from a sampled received signal. For a length L discrete complex baseband reference waveforms[n], and the signal transmitted by one radio and received and sampled by another r[n], the time delay is computed from the cross-correlation of s[n] and r[n]. This is the 85 sequence d[n] given by d[n] =r[n]~s [n] DFT ( = =) D[k] =R[k]S [k] (5.40) where the~ operator indicates circular convolution. S[k], R[k], and D[k] are the discrete Fourier transforms (DFTs) of s[n], r[n], and d[n] respectively, dened for D[k] as D[k] = L1 X n=0 d[n]e j 2nk=N (5.41) A number of methods exist for estimation of autocorrelation true peak location to sub-sample precision, including interpolation and slope estimation of the spectral phase [126], [127], [128], [129]. In practice, interpolation is an inecient approach as the accuracy of the estimated peak location is directly dependent on the upsampling factor. In order to obtain picosecond-level precision from a signal sampled at 50 MHz, interpolation by a factor of 20; 000 would be required. For a length 2048 sample sequence, this is not feasible. In this section, we examine three methods for estimating the true sub-sample peak location of the correlation sequence: spectral phase slope estimation, quadratic least-squares (LS) tting, and and a new sinc nonlinear least squares (NL-LS) algorithm. Simulated results are compared with the Cramer-Rao bound for one-way TOF. 5.5.1 Cramer-Rao Lower Bound Here we give the established CRLB for one-way TOF measurements in noisy environments with pulse compression. The CRLB minimum variance one-way TOF measurement 2 TOF for a single pulse non-analytic complex signal (LFM chirp with rect window) is [130], [128] 2 TOF 3 2(B) 2 SNRT p f s (5.42) The CRLB as stated in (5.42) is shown in Fig. 5.6. A derivation of (5.42) is given in Appendix C. We note that in this work, we use critically sampled waveforms such that B =f s . 86 5.5.2 Spectral Phase Slope Shifting of a time domain sequence appears as a linear phase term in its Fourier transform [131]. Thus a fractional delay in a discrete time signal can be retrieved via linear LS estimation of the spectral phase slope. Using the fundamental properties of the Fourier transform, spectral linear phase slope estimation is done as follows: We rst compute the DFT of d[n], D[k], and take the spectral phase term [k] =\D[k] over the sequence of frequencies in the DFT f[k] = fs L (k L 2 ). We compute the linear LS slope estimate of the spectral phase as: = Lhf;ihf; 1ih; 1i Lhf;fihf; 1ihf; 1i (5.43) whereh;i denotes the scalar inner product, and 1 is the unity vector. The sub-sample true peak estimate ~ t pk is then ~ t pk = 2 (5.44) Spectral phase slope estimation is sensitive to noise and performs comparatively poorly in cases where SNR is low as seen in Fig. 5.6. This is because this method does not take advantage of the matched lter gain which is only realized in the time domain. The performance of this method can be improved for low SNR signals by applying a window function centered around the integer peak index n pk of the autocorrelation signal before taking the DFT. For a general window function h[n] of size 2L w + 1, we can dene an amplitude weighting function a[n] a[n] = 8 > > < > > : h[nn pk ]; n2fn pk L w ;n pk +L w g 0; otherwise (5.45) D[k] = L1 X n=0 a[n]d[n]e j 2nk=L (5.46) 87 This is analogous to a bandpass lter and eectively acts to reduces the spectral noise by ltering it in the time domain. In order to preserve the amplitude structure of the autocorrelation signal, a at top window such as a rectangular or tapered cosine (Tukey) should be chosen. This is because the window as dened in (5.45) is in general not centered around the true peak of the autocorrelation signal, and, as a result, a non- at top window top window will introduce bias errors in the true peak estimate. 5.5.3 Quadratic Least-Squares The approach of modeling the autocorrelation peak as a quadratic has been used widely [126], [127], [128]. The peak location is estimated using just the 3 surrounding points, making this estimation algorithm O(1). Because it is imperative that the complete synchronization process execute with as little latency as possible, peak estimation algorithms that execute in constant time are desirable. We nd the maximum peak index n pk of the matched lter output n pk = argmax n jd[n]j (5.47) Using the length-three sample sequence centered about the peak to form the column vectory, where y :y i = log 10 (jd[n pk 1 +i]j) i2f0; 2g (5.48) we perform quadratic least squares (LS) estimation and dierentiate the LS quadratic polynomial to obtain a sub-sample rate estimate of the true peak location ~ n pk . ~ n pk =n pk y 2 y 0 2y 0 4y 1 + 2y 2 (5.49) For sample rate f s the time of the true peak estimate ~ t pk is given by ~ t pk = ~ n pk f s (5.50) 88 At high SNR quadratic LS peak estimation suers from bias that appears as an error that exhibits sinusoidal behavior as a function of fractional waveform delay between integer clock cycles, as seen in Fig. 5.6. Performing upsampling and interpolation by a small integer factor prior to quadratic LS peak estimation is a viable method of reducing this bias error oor to a desired level at the cost of algorithmic eciency. Interpolation may, for example, be performed by zero-padding the DFT, thus increasing the size requirement of the inverse DFT by an integer factor. We note that for high SNR and noiseless signals, spectral slope estimation achieves higher precision than the quadratic LS method as is illustrated in Fig. 5.6. However, the tradeo in algorithmic complexity and performance in realistic signal environments ultimately make quadratic LS peak estimation preferable in practice. 5.5.4 Sinc Nonlinear Least-Squares NL-LS estimation is used widely for peak tting and estimation and has been studied extensively. In [132], NL-LS was used with a hyperbolic and Gaussian functions for peak estimation and tracking for time of arrival (TOA) signals. We note that the peak estimate obtained by NL-LS tting with a Gaussian function is identical to that obtained by quadratic LS estimation of the logarithm of the matched lter output. This result is known as Caruanas algorithm [133], [134] and is due to the quadratic form of the exponential term in the Gaussian function. In this work, we use a sinc function kernel to estimate the autocorrelation peak for pulse compression waveforms using NL-LS. This peak estimation performs well in low-SNR and avoids the biases seen in quadratic least squares estimation. As with the quadratic LS algorithm, only 3 sample points are used, making this algorithm O(1). Next, we derive the sinc-based NL-LS estimation algorithm. 89 Given a known functionf(x;) that depends on input column vectorx as well as parameters in the vector and produces a column vector outputy having the same dimensions asx, we estimate parameter values that minimize the residual error. f(x;) = 0 sinc (x 1 ) 2 (5.51) y :y i =jd[n pk 1 +i]j i2f0; 2g (5.52) x = [1 0 1] T (5.53) = [ 0 1 2 ] T (5.54) We setup a cost function to minimize the residual error S = 2 X i=0 (y i f(x i ;)) 2 (5.55) To solve this we use Gauss-Newton optimization, for which we need the gradients with respect to the model parameters @f(x;) @ 0 = sinc 2 (x 1 ) (5.56) @f(x;) @ 1 = 0 sinc 2 (x 1 ) cos 2 (x 1 ) x 1 (5.57) @f(x;) @ 2 = 0 cos 2 (x 1 ) sinc 2 (x 1 ) 2 (5.58) We initialize 0 =jd[n pk ]j, 1 = 0, and 2 = B fs Now nonlinear least squares tting is performed iteratively. At the m th iteration, the matrix Jacobian J and residual error y are computed J = [ @f(x; m ) @ 0 @f(x; m ) @ 1 @f(x; m ) @ 2 ] (5.59) y =yf(x; m ) (5.60) 90 The matrix equation must then be solved for J = y (5.61) We note that for estimation from 3 sample points, this is a square matrix and this matrix equation may be solved directly. If more sample points are used, the system is over-determined and we instead use the normal equation solution: = (J T J) 1 J T y (5.62) Once is found, the parameters estimates are updated for each iteration m as m+1 = m + (5.63) After the nal iteration, we obtain the estimated parameters ~ and the true peak estimate is given by ~ n pk =n pk + ~ 1 (5.64) ~ t pk = ~ n pk f s (5.65) which includes ~ 1 as the oset term as it corresponds to a shift in the model sinc function stated in (5.51). Because we are able to provide initial peak time values that are close to the true values, the convergence of the algorithm is rapid, usually occurring in < 4 iterations in practice. In general, this method performs well with only three sample points. While the sinc NL-LS is slightly more computationally intensive than quadratic least squares, the performance improvements are signicant as will be shown next. In addition, the three sample sinc NL-LS algorithm is O(1). 91 5.5.5 TOF Peak Estimation Algorithm Performance Performance of each algorithm as a function of SNR is shown in Fig. 5.6 for two types of waveforms: a linear frequency modulted (LFM) chirp in Fig. 5.4a and a pseudorandom noise (P/N) sequence in Fig. 5.4b. In both cases, the proposed sinc NL-LS algorithm performs as well as quadratic LS for low SNR. For high SNR, sinc NL-LS signicantly outperforms quadratic LS and matches the performance of spectral phase slope estimation. Simulated results show that the sinc NL-LS algorithm performance achieves the CRLB for SNR2 [15; 40] dB for both waveforms. -20 0 20 40 SNR (dB) 10 -2 10 0 10 2 log RMSE (ns) phase slope quadratic LS sinc NL-LS CRLB (a) LFM chirp -20 0 20 40 SNR (dB) 10 -2 10 0 10 2 log RMSE (ns) phase slope quadratic LS sinc NL-LS CRLB (b) P/N code Figure 5.4: Performance comparison of cross-correlation fractional peak estimation methods versus SNR. A pulse length N = 1024 and sampling rate f s = 50 MHz are used in both simulations (Time-bandwidth product (TxBW) = 1024). On average, sinc NL-LS estimation converges in < 4 iterations and is limited to a maximum of 5 iterations. The radar time delay TOF CRLB given in (5.42) is shown for comparison. 5.6 Processing and Exchange of Information Protocol In our implementation we encode and exchange time of ight data between the sensors wirelessly using quadrature phase-shift keying (QPSK) modulation. Dene ~ t pk:i;j as the TOF delay as measured by sensori for the synchronization pulse sent by sensorj using the peak detection methods described in Section 5.5. In stage 1 as shown in Fig. 5.5, this value is encoded as a length N m QPSK message m i;j [n] that represents each sensor's combined estimate of channel delay and fractional clock phase 92 Matched Filter Peak Detect True Peak Estimate Matched Filter Peak Detect QPSK Decode True Peak Estimate QPSK Encode Waveform Split Construct TOF Matrix Clock Drift Estimate RF Phase Error Estimate Clock Phase Estimate Matched Filter Symbol Phase Sync Stage 1 Stage 2 Stage 3 TDMA TX TDMA RX TDMA RX TDMA TX Figure 5.5: Synchronization Processing diagram oset. m i;j [n] = QPSK Nsps ( ~ t pk:i;j )[n] (5.66) where QPSK Nsps ()[n] is dened as a non-linear operator that performs quantization of a scalar input and creates a QPSK encoded waveform with the quantized bits encoded as symbols and with N sps samples per symbol. TOF messages obtained from every other sensor are appended to the end of sensor i's synchro- nization waveform s i [n], which has support over the interval [0;L 1] and is zero elsewhere, to produce the stage 2 TX waveform ~ s i [n] as shown in Fig. 5.5. ~ s i [n] =s i [n] + N2 X j=0 m i;j [nLjN m ] (5.67) Now each sensor transmits its stage 2 waveform ~ s i [n] in the allocated TDMA time slot, while every other sensor receives and decodes the messages. Once this operation is complete, every sensor in the network will have complete knowledge of the local pair-wise TOF measurements of all sensors and the synchronized TOF between them as well as the relative RF phase measurements. From this information, each sensor constructs identical TOF (5.22) and clock phase error (5.23) matrices for the entire network as shown in stage 3 of Fig. 5.5. 93 This broadcast messaging scheme is extended to encode additional information including residual carrier phase error measurements from (5.39), which are used estimate and correct TX/RX LO phase osets as described in Section 5.4.3. Because, the proposed synchronization scheme requires error free exchange of messages between sensors, we implement a forward error correction (FEC) coder/decoder (codec) that is software congurable to utilize the lowest complexity and highest rate code possible given the SNR conditions of the operating environment. The architecture of the software dened codec is shown in Fig. 5.6a. The message error rate achieved by two of the implemented encoding schemes (hamming (7,4) and convolutional code with constraint length 3 and polynomial [7,5]) using 8 samples per QPSK symbol are shown as a function of SNR in Fig. 5.6b. We note that convolutional codes with arbitrary constraint lengths and polynomials are supported. (a) Software dened FEC codec architecture for message encoding. (b) Software dened FEC codec decoding error as a function of SNR. Figure 5.6: Architecture and performance of software dened FEC codec architecture implemented in SDR hardware for synchronization message exchange. 5.7 Sensor Localization When all sensors have knowledge of the TOF between all sensors in the network as a TOF matrix, it is possible to use the TOF to estimate the positions of each sensor [48,135{138]. Here we derive a localization algorithm that uses the TOF matrix derived in Eqn. (5.22). Because each sensor has the identical TOF matrix, in principle, each sensor can independently determine the positions of all of the other sensors. 94 Dene X as an N 3 matrix with columns corresponding to the free space coordinates (x;y;z) of each sensor in the network. We denote the i th row of X as x i = (x i ;y i ;z i ). ~ R is the NN matrix of relative RF time of ight range measurements, which is known by each sensor in the network. ~ R i;j =c ~ TOF i;j indicates the line of sight distance between sensor i and j. We note that ~ R is symmetric with zeroes along the diagonal. The relative coordinates of all nodes in the network, ~ X, may be estimated by minimizing the following cost function: ~ X = argmin x;y;z N2 X i=0 N1 X j=i+1 (jjx i x j jj l 2 ~ R i;j ) 2 (5.68) where the operatorjj()jj l 2 indicates the L-2 Euclidean norm. Equation (5.68) may be minimized using any number of Newtonian or gradient based methods, which could be implemented locally on sensor nodes. Once the solution of ~ X is found, each node knows the position of all nodes in the network. 5.8 Experimental Characterization of Synchronization Method In this section, we present results from selected experiments that demonstrate the performance of the synchronization algorithm on the SDRadar platform described in Chapter 3. 5.8.1 3-Array TOF synchronization (a) 1 m separation (b) 5 m separation (c) 10 m separation Figure 5.7: Wireless synchronization test: three array of sensors in equilateral triangle formation 95 In this experiment, we demonstrate the lower-bound performance of the wireless synchronization algorithm for three SDRadar sensors. The sensors are arranged at the vertices of an equilateral triangle. The three experimental setups (1 m, 5 m, and 10 m separation) are shown in Fig. 5.7. The synchronization performed at an RF center frequency of 1 GHz and repeated once per second for 1000 seconds. The performance is dependent on the SNR. A high SNR is maintained during the experiment by increasing the TX gain in the SDRadar as the distance between sensors increases. = 25.8 ps = 35.0 ps = 22.8 ps (a) 1 m sensor separation = 50.5 ps = 39.8 ps = 40.7 ps (b) 5 m sensor separation = 27.3 ps = 31.9 ps = 38.3 ps (c) 10 m sensor separation Figure 5.8: Three sensor array synchronized TOF Matrix measurements for 1000 trials taken over 1000 s ( 17 minutes) demonstrating sub-100 ps precision. Test setup is shown in Fig. 5.7. Statistics are given in Table 5.1. Note that the x-axis is scaled to span 250 ps around the mean . 96 Fig. 5.8a, Fig. 5.8b, and Fig. 5.8c show the synchronization precision as histogram plots of the synchronized TOF measurements for all sensor pairs. Sub-100 picosecond (ps) synchronization precision is obtained in an outdoor LOS environment using the method in Section 5.4, where each reported TOF estimate is obtained from just two transmissions per sensor. In Table 5.1, statistics for the experiments are given. Data is shown for one of the sensors, because the complete TOF sync information is known identically among all sensors in the network. Consequently, ~ TOF i;j = ~ TOF j;i exactly as the proposed synchronization protocol (in the absence of decoding errors) guarantees shared knowledge of all TOF measurements across the network. We use the 1 m separation test to calibrate the sensor TOF mean values reported in Table 5.1. Table 5.1: 3-Array TOF Synchronization Precision. 1000 trials taken over 1000 seconds at 1 GHz Test Standard Deviation Dev 1! Dev 2 Dev 1! Dev 3 Dev 2! Dev 3 1 m Separation .77 cm (25.8 ps) 1.05 cm (35.0 ps) .68 cm (22.8 ps) 5 m Separation 1.52 cm (50.5 ps) 1.19 cm (39.8 ps) 1.22 cm (40.7 ps) 10 m Separation .82 cm (27.3 ps) .96 cm (31.9 ps) 1.15 cm (38.3 ps) Test Mean Dev 1! Dev 2 Dev 1! Dev 3 Dev 2! Dev 3 1 m Separation 1 m (3.33 ns) 1 m (3.33 ns) 1 m (3.33 ns) 5 m Separation 4.5 m (14.8 ns) 5.1 m (16.9 ns) 5.6 m (18.8 ns) 10 m Separation 9.9 m (32.8 ns) 10.1 m (33.7 ns) 10.6 m (35.4 ns) 5.8.2 3-Array TOF Localization (a) Experiment setup 40 35 30 25 20 15 10 5 0 y (m) -30 -20 -10 0 10 20 x (m) Positions (b) Estimated path Figure 5.9: Sensor position estimation using three sensor array and TOF obtained using proposed synchronization method. 97 In this experiment we demonstrate localization using three sensor triangulation. We have xed the position of two sensors with a separation of 10 m. The third SDRadar sensor is moved along a variety of recognizable paths. The position of the moving sensor is estimated from the TOF matrix using constrained solution triangulation given in Section 5.7. The test setup is shown in Fig. 5.9a and the estimated sensor path is shown in Fig. 5.9b, where a 10 point moving average has been applied to the estimated positions. When the third sensor is too close to the axis formed by the two xed sensors, the solution is ill formed, which is shown by the feature in the position paths in the lower left of Fig. 5.9b. 5.8.3 3-Array TOF Long Range synchronization In this experiment, the results of which have been published in [139], we test the synchronization performance over larger distances, with two sensors xed along a 10 m baseline while the third sensor is moved away in 50 m steps out to a maximum distance of 300 m as shown in Fig. 5.10a with GPS coordinates of the sensors shown in Fig. 5.10b. This test was performed at the Sante Fe Dam RC Aireld in Irwindale, CA at a frequency of 2:4 GHz. The aireld is used for recreational ight of model RC aircrafts, of which multiple were present and being actively used throughout the test. Multiple remote controllers with high power 2:4 GHz transmitters operating in close vicinity caused signicant levels of interference and in some cases were observed to saturate the SDRadar receivers. While this degraded synchronization performance, sub-nanosecond synchronization was still achieved throughout the test as shown in Fig. 5.10d. This experiment demonstrates the robustness of the synchronization technique and implemented software dened codec to interference at levels far higher than what would be expected in nominal operational environments. 5.8.4 Bistatic Wireless Re ector Test In order to obtain high resolution, we can coherently combine multiple frequency measurements using stepped-frequency radar techniques [111], [4] to reconstruct a synthetic wideband waveform (SWW). In order to achieve this bistatically, the developed synchronization algorithm is used to wirelessly synchronize two SDRadar sensors. The USRP E312 SDRs used have two RX and two TX 98 (a) Experiment setup (b) GPS coordinates 0 200 400 600 800 1000 0 50 100 150 200 250 300 Distance (m) TOF sync (corrected) dev A -> dev B dev A -> dev C dev B -> dev C 1 2 3 4 5 6 7 Scan (c) Synchronized TOF distance measurements 1 2 3 4 5 6 7 Interval 2000 2200 2400 2600 TOF (ns) Avg TOF Sync (dev A) dev A -> dev B dev A -> dev C 1 2 3 4 5 6 7 Interval 0.2 0.4 0.6 0.8 StdDev (ns) StdDev TOF Sync (dev A) dev A -> dev B dev A -> dev C (d) Interval statistics showing average TOF (top) and TOF standard deviation (bottom) Figure 5.10: 3-Sensor synchronization measurements performed at Sante Fe Dam RC Aireld in Irwindale, CA. TOF Synchronization was performed at 2.4 GHz. Sensor A moved away from two-sensor xed baseline (sensors B and C) in 7 50 m steps. High levels of 2.4 GHz interference due to multiple high power model RC airplane transmitters operating in close vicinity was present throughout the experiment. channels. Using one TX/RX pair for synchronization and the second TX/RX channel for radar measurement, we synthesize ultra-wideband radar pulses by coherently combining a sequence of smaller bandwidth sub-bands. All stepped frequency sub-pulses are synchronized and independently, with the synchronization routine running at each frequency step. The synchronization method is used to (i) correct the clock phase error (time synchronize), (ii) correct the frequency-dependent LO phase error due to 99 sample clock oset (RF carrier phase synchronize), (iii) correct the random phase error due to LO re-tuning which is characteristic of the AD9361 RFIC frontend (RF carrier phase synchronize), and (iv) remove frequency dependent signal path delays due to hardware. Only after all of these steps have been successfully completed is coherent reconstruction of the high resolution SWW possible. (a) Experiment setup 180 175 170 165 160 155 150 145 140 135 Range (m) dB (b) Synchronized bistatic radar measurements Figure 5.11: Bistatic radar test using wireless synchronization algorithm to synchronize two sensors across all 50 MHz bandwidth sub-bands so that coherent stepped frequency radar may be performed (391 MHz total bandwidth synthesized) The two-sensor setup for the experiment is shown in Fig. 5.11a. We synthesize 391 MHz (38:4 cm theoretical radar resolution) of bandwidth, stepping from 1211 to 1552 MHz in 25 MHz steps and collecting 50 MHz bandwidth LFM chirp sub-pulses at each of the 16 frequency steps. The capability of the developed synchronization algorithm and SDRadar to perform coherent stepped frequency radar imaging using a wireless link for synchronization is demonstrated experimentally in Fig. 5.11b. Measurement statistics showing precise localization of re ector targets to within 10 cm are given in Table 5.2. Measured distances provided are from the legs of the camera stand holding the radar to the corner point of the re ector on the ground. The antenna height above the ground surface is 1:1 m. The measured3 dB target echo peak widths are also given in Table 5.2. The mean direct path3 dB peak width for all trials was 32:7 cm and the mean3 dB peak width of the target echo for all re ector positions was 36:6 cm. For reference, the3 dB down peak width for an ideal 391 MHz LFM chirp waveform was estimated as 33:9 cm using the same algorithm, thus 100 demonstrating that proposed synchronization method achieves the time and phase coherence across multiple independent frontend frequency bands necessary to reconstruct a SWW that achieves the theoretical bandwidth resolution performance. Stepped frequency radar requires sub-pulses to be both time and phase coherent in order to realize resolution improvements [4]. Therefore, this experiment demonstrates that coherent wireless bistatic/multistatic radar operation is feasible with the proposed synchronization scheme. Table 5.2: Bistatic Re ector Range Test Target Locations: Actual vs. Measured Range and Echo Characteristics Result Re ector Location 1 2 3 4 5 6 7 True Range (m) 10.77 10.87 10.97 11.07 11.38 11.67 13.42 Measured Range (m) 10.85 10.93 11.09 11.22 11.56 11.89 13.63 Std. Deviation (cm) 1.88 3.92 5.05 3.38 2.79 1.91 3.44 Std. Deviation (ps) 125.6 261.3 336.9 225.1 185.8 127.5 229.2 -3dB Peak Width (cm) 32.2 33.0 35.7 58.2 34.8 30.2 32.4 5.8.5 3-Sensor Transmit Synchronization In this experiment we demonstrate phase coherent transmit synchronization for three sensors operating in MIMO fashion. That is, all sensors transmit mutually orthogonal waveforms (in this case via TDMA scheme) which in turn are received by all sensors. Using the synchronization scheme detailed in Section 5.4, we perform time and phase synchronization of baseband waveforms to the global average clock phase oset prior to transmission for each sensor by applying the time and phase pre-corrections derived in Section 5.4.2. We note that although no MIMO processing is performed using the synchronized signals, the purpose of this test is to demonstrate the precision achieved with the proposed method and to support our claim of its feasibility for wirelessly synchronizing elements of a coherent MIMO or multistatic array. While the demonstration of two sensor synchronization is insucient to prove the validity of a given synchronization scheme forN sensors, demonstration of the three sensor case does indeed prove global synchronization and extends to the N sensor case. This is because in the three sensor case, synchronization of a transmitter with two independent receivers or of two independent transmitters with a third independent receiver can only occur if the transmitted signals are synchronized in the air. 101 20 40 60 80 100 trial -4 -2 0 2 4 phase (rad) dev 1 RX: peak phases dev 1 <- dev 2 (onboard) dev 1 <- dev 2 (offline) dev 1 <- dev 3 (onboard) dev 1 <- dev 3 (offline) 1750 1800 1850 1900 time (ns) 120 130 140 150 160 dB matched filter response dev 1 <- dev 2 dev 1 <- dev 3 Dev 1 ← Dev 2 Dev 1 ← Dev 3 (a) Sensor 1 RX Signals: Dev. 1 Dev. 2 and Dev. 1 Dev. 3 20 40 60 80 trial -4 -2 0 2 4 phase (rad) dev 2 RX: peak phases dev 2 <- dev 1 (onboard) dev 2 <- dev 1 (offline) dev 2 <- dev 3 (onboard) dev 2 <- dev 3 (offline) 1750 1800 1850 1900 time (ns) 120 130 140 150 160 dB matched filter response dev 2 <- dev 1 dev 2 <- dev 3 Dev 2 ← Dev 1 Dev 2 ← Dev 3 (b) Sensor 2 RX Signals: Dev. 2 Dev. 1 and Dev. 2 Dev. 3 20 40 60 80 trial -4 -2 0 2 4 phase (rad) dev 3 RX: peak phases dev 3 <- dev 1 (onboard) dev 3 <- dev 1 (offline) dev 3 <- dev 2 (onboard) dev 3 <- dev 2 (offline) 1750 1800 1850 1900 time (ns) 120 130 140 150 160 dB matched filter response dev 3 <- dev 1 dev 3 <- dev 2 Dev 3 ← Dev 1 Dev 3 ← Dev 2 (c) Sensor 3 RX Signals: Dev. 3 Dev. 1 and Dev. 3 Dev. 2 Figure 5.12: 3-Sensor MIMO transmit synchronization test results. For the signals received by each of the three SDRadar sensors, the waveform phase in radians (left), and the signal time domain matched lter response (center), and the time domain waveforms samples (right) are shown for the signals transmitted by each of the other two sensors. Statistics are given in Table 5.3. Fig. 5.12 shows results from the full 3x3 MIMO transmit synchronization test. The test is performed at a frequency of 1:1 GHz over 100 trials at a PRI of :2 s. The raw time domain signal samples, matched lter response, and the carrier phase are shown for all 6 MIMO signals. The 3 monostatic cases, where a given sensor receives the signal transmitted by itself, are omitted. The reader may observe that the residual RF phase errors are symmetric. That is, the un- compensated RF carrier phase, as shown in Fig. 5.12, i.e., for Dev. 1 Dev. 2 and Dev. 2 Dev. 1 exhibit symmetry. This is due to the relatively large processing latency of the current embedded software implementation over which time, each clock drifts relative to the others. This symmetry indicates that better performance may be achieved by a faster implementation with less time between the synchronization epoch and the radar pulse transmission. This is discussed further in Section 5.10. 102 Fig. 5.12 shows phase coherent transmit to sub-nanosecond time precision and phase precision of 5 (i.e., 10 ). The results presented have not undergone any post-processing and represent the raw data recorded by each sensor. The experimental statistics are given in Table 5.3. The corrections applied to the transmit waveforms assume that local clock osets remain xed from the time the synchronization is performed to when the synchronized radar pulses are transmitted. As previously noted, this is not the case, and the addition of a predictive model of the clock drift, which may be accurately treated as linear in nature over short time spans, would yield improved performance. Implementation of such a state estimation model is the subject of future work. Table 5.3: 3-Sensor MIMO TX Sync Statistics Param. Device exchange 1 2 1 3 2 1 2 3 3 1 3 2 Estimated SNR (dB) 29.5 30.7 35.9 35.7 37.5 31.2 CRLB 2TOF (ps) 4.4 3.8 2.1 2.1 1.7 3.6 Std. Dev. Oine Sync (.02 s eective processing latency) TOF sync (ps) 12.4 7.6 12.4 6.2 7.6 6.2 TX Radar TOA (ps) 29.8 27.4 18.2 31.7 27.5 30.8 TX RF Phase (rad) 0.239 0.192 0.198 0.230 0.243 0.237 Std. Dev. Onboard TX Sync (.16 s processing latency) TOF sync (ps) 14.2 19.0 14.2 11.3 19.0 11.3 TX Radar TOA (ps) 91.7 110.9 79.5 119.0 108.7 114.3 TX RF Phase (rad) 0.525 0.550 0.471 0.619 0.632 0.629 The time between the sync operation and the transmission of the corrected/synchronized TX pulses is:16 s. This delay is due to the synchronization processing being performed in real-time running software on the embedded Zynq-7020 SoC ARM processor. In this same experiment, however, after the synchronization epoch we also transmit and receive two sets of pulses: one through the calibration/synchronization channel (the reference pulse), and one through the radar antenna channel (the data pulse) which are saved directly to le. The pulses on the calibration channel are identical to those sent during the synchronization epoch, however they are saved rather than processed onboard. Because there is no onboard processing performed between TX/RX of the reference pulses and the data pulses, the time between them is:02. This allows us to perform oine processing to `see' what the output of the synchronization would have been if the time between the synchronization exchange and TX synchronized radar pulse transmission was only :02. 103 We provide the results from oine synchronization processing for the :02 second reference pulse and data pulse delay in Table 5.3. Here, the phase precision is 14 (i.e., 28 =:036). We further note that the :02 s delay between the reference and data pulses is primarily due to le write operations for saving the raw data and it is reasonable to expect that a faster implementation of the synchronization processing, which does not require le I/O, could execute with much lower latency. The SNR values reported in Table 5.3 are estimated directly from the raw signal data and are used to calculate the expected CRLB for comparison. The two way TOF CRLB 2TOF reported in Table 5.3 is related the the TOF CRLB TOF given in (5.42) as 2TOF = TOF p 2 (see Appendix C for further discussion). We note that the two-way TOF CRLB is highly sensitive to SNR, which in this case has a considerable degree of uncertainty as it is estimated solely from the raw signal data. 5.9 Experimental System Applications 5.9.1 Snow Penetrating Radar Test: Monostatic and Bistatic In this experiment, two SDRadar sensors are synchronized to perform bistatic imaging of a snowbank in Mammoth Lakes, California. For comparison, a single SDRadar sensor is used in a monostatic conguration to image the same scene. We note that the monostatic result has been previously reported in [4] and is included here for comparison. In both the monostatic and bistatic cases, 2:5 GHz of bandwidth are synthesized from an operational frequency range of 600 3100 MHz using a non-uniform frequency step size as described in [4]. In the monostatic case, a single sensor is suspended from a xed line 1:5 m above the surface of the snowbank and moved across the 25 m scene in :5 m steps. The monostatic test setup and resulting sub-surface radar image, along with a ground truth image of the snow and ice layers present in the snowbank, are shown in Fig. 5.13a. The bistatic snow penetrating radar experiment setup is shown in Fig. 5.13b. One SDRadar sensor (Sensor 1) is suspended from a xed line in a static position at the center of the transect and a second SDRadar sensor (Sensor 2) is moved across the 25 m transect in :5 m steps on a second xed line. The proposed method is used to synchronize the two SDRadar sensors at each frequency tuning step across the 2:5 GHz bandwidth used (72 frequency steps in 36 MHz increments). By correcting clock and carrier phase errors, the entire 2:5 GHz bandwidth is combined coherently to form a SWW 104 (a) Snow layers and depth imaged by a monostatic SDRadar Direction of Sensor Motion Sensor 1 Sensor 2 device B: scaled bistatic 5 10 15 20 25 30 Scan 0 1 2 3 4 5 6 7 8 Depth (m) Ground Layer Direct Path Snow Surface GPS Sync Only With Sync Algorithm 5 10 15 20 25 30 Scan 0 1 2 3 4 5 6 7 8 Depth (m) -5 0 5 10 15 Depth (m) (b) Snow layers and snow depth imaged by bistatic SDRadars. Results obtained with GPS- based synchronization only and the proposed synchronization algorithm are compared. Figure 5.13: Monostatic and bistatic radar test. 2:5 GHz of bandwidth are synthesized from an operational frequency range of 600 3100 MHz using a non-uniform frequency step. In the bistatic case, the SDRadars use the proposed wireless synchronization algorithm to synchronize two sensors across all frequency bands so that coherent stepped frequency radar may be performed. prole at each point along the transect. A comparison of the radar imaging results using GPS-based synchronization only and the method described in Section 5.4 is given in Fig. 5.13b. The bistatic synchronized radar image with annotated features is shown in Fig. 5.13b. Due to the directionality of the antennas, the re ection from the snow surface and near-surface features are prominent only in the center of the image (when the two sensors are close to one another). For near-surface features, the monostatic case shown in Fig. 5.13a demonstrates superior imaging ability. However, the re ection from deeper targets, specically the snow/ground interface as shown in Fig. 5.13b, are signicantly stronger. Further, the physical separation of antennas in the bistatic case improves the dynamic range of the entire system as higher transmit power may be used without saturating the 105 receiver, allowing for better imaging of deeper targets. Thus, in a 2x2 MIMO conguration, wherein both monostatic and synchronized coherent bistatic imaging are performed, it is possible to achieve both high resolution of shallow targets as well as improved imaging of deeper features. 5.9.2 Bistatic Linear Aperture Test (a) Test set up device B: scaled bistatic 10 20 30 40 50 60 70 80 90 100 Scan 0 50 100 150 200 250 300 350 400 450 500 Range (m) 500 m (b) Bistatic radar image Figure 5.14: Wireless two-sensor bistatic aperture measurement In this test, we demonstrate the ability of a wireless two SDRadar bistatic sensor system to image targets at far range. We coherently synthesize 200 MHz total bandwidth (1:5{1:7 GHz in 25 MHz non-uniform steps). At each frequency step, synchronization is performed to enable coherent wideband synthesis. The test setup and target scene are shown in Fig. 5.14a. Sensor 1 is moved in 20 cm increments across the aperture (98 steps total) while Sensor 2 remains xed. We note that after the 68 th step, Sensor 2 is repositioned at the end of the initial 68 step aperture. The resulting bistatic radar image is shown overlain on a satellite image of the test site in Fig. 5.14b. Due to the relatively small size of the aperture relative to the scattering scene, SAR azimuth focusing is not performed. 5.10 Discussion Each sensor's clock drifts over the time between the synchronization operation and the transmission of the corrected waveforms, causing a decoherence of the relative clock states estimated by the synchronization exchange from the relative clock states at the time of signal transmission. The results of the associated reduction of precision are reported in Table 5.3. 106 There are three ways this issue could be addressed, which are independent of the synchronization algorithm itself: 1. higher quality oscillator with greater stability and lower phase noise (relax synchronization repetition frequency requirements) 2. higher performance embedded processor (decrease processing latency for software implementa- tion of synchronization algorithm) 3. implementation of synchronization algorithm in FPGA (remove FPGA-to-processor data transfer bottleneck and perform synchronization processing with deterministic latency) The synchronization processing as shown in Fig. 5.5, is relatively simple and could easily be performed onboard in milliseconds by a more powerful processor or FPGA. In our implementation, the TDMA slot size is 100 s, therefore the lower bound synchronization latency is is 2N 100 s. Furthermore, a limitation of our current hardware demonstration platform is that the two TX/RX channels cannot operate independently in parallel. In a more powerful hardware platform, which is capable of asynchronous multichannel operation, it would be possible to execute the synchronization algorithm on a separate dedicated side-channel asynchronously w.r.t the radar data channel. In such a system, more sophisticated processing of synchronization results such as Kalman ltering could be used to improve clock state estimation and the precision of time and phase synchronization. 5.11 Conclusion In this chapter, we have presented a consensus synchronization algorithm for distributed wireless sensor networks. The algorithm relies on syntonization and course synchronization from GPS signals and then improves synchronization precision by three orders of magnitude. The proposed method requires that each sensor transmits twice in a synchronization epoch, meaning that the entire synchronization process is complete for the entire network of N sensors after 2N transmissions, making it an O(N) algorithm. Using the results from the synchronization procedure, we have formulated time and RF carrier phase corrections to baseband waveforms that may be applied on 107 transmit in order to obtain network wide coherent transmit operation as well as coherent receive operation. Furthermore, the network pair-wise LOS distances are known globally as a result of the proposed scheme enabling decentralized simultaneous positioning of the network. Because the scheme does not rely on real-time signal detection and response or any specialized hardware, it is implementable entirely in software on commercially available SDR platforms and is shown to achieve < 100 ps time synchronization performance with 50 MHz signal bandwidth. We have validated the algorithm and system in numerous eld experiments. We have demonstrated 3-sensor wireless time synchronization to< 50ps (and as low as 10ps in some cases), coherent phase synchronization of =28, cm-level positioning. We have applied the synchronization method across multiple tuning frequencies to perform bistatic stepped-frequency radar coherently to synthesize bandwidths of up to 2:5 GHz. This work has numerous implications for distributed transmit beam forming, wireless sensor localization, and coherent MIMO radar imaging. In particular, this work is an enabling technology for low-cost high performance coherent MIMO radar sensor networks made up of SmallSat/CubeSat sensor constellations and autonomous sensor swarms. 108 Chapter 6 Characterization of Clock Phase Errors for Distributed Wireless Synchronization Protocol In this chapter, we expand upon results previously reported in [140]. 6.1 Introduction The problem of synchronization is critical for the realization of next-generation distributed coherent radio and radar wireless sensor networks. A recent decentralized wireless synchronization protocol was demonstrated to achieve sub-nanosecond coherent synchronization using a wireless network of low-cost commercial-o-the-shelf (COTS) software dened radio (SDR)-based radar sensors [32]. In order to evaluate the method proposed in [32] for use in large-scale distributed radar missions, such as distributed or multiple-input-multiple-output (MIMO) synthetic aperture radar (SAR) from smallSAT/cubeSAT constellations and small unmanned aircraft system (sUAS) swarms, extensive simulation and analysis are critical. To this end, we derive complete analytical expressions for the distributed synchronization scheme that relate an arbitrary input oscillator PSDs to output PSDs following network synchronization. This eort enables ecient modeling of large scale distributed system synchronization and optimization of mission parameters. 109 6.2 Modeling Oscillator Phase Noise The two-sided PSD of an oscillator output phase noiseS TS ;osc (f) can be modeled as [141] [142] [143] [144] S TS ;osc (f) = a 4 f 4 + a 3 f 3 + a 2 f 2 + a 1 f 1 +a 0 (6.1) wherea 4 :::a 0 are determined from the measured oscillator-phase noise characteristics using asymp- totic approximation [141]. These terms are due to 1) random walk frequency noise; 2) frequency icker noise; 3) white frequency noise; 4) icker phase noise; and 5) white phase noise [142]. 10 -2 10 0 10 2 10 4 10 6 Freq (Hz) -150 -100 -50 0 S ? (f) (dBc/Hz) TCVCXO OCXO DOCXO Rubidium + DOCXO Typical USO (Krieger) Figure 6.1: Oscillator output noise PSD S ;osc (f). Power law coecients are estimated from datasheet specications using LS t. Using measurements of the temperature compensated voltage controlled oscillator (TCVCXO) used in the Ettus Research USRP E312 COTS SDR, we empirically determine the parameters a 4 =66 dB a 3 =62 dB a 2 =80 dB (6.2) a 1 =110 dB a 0 =153 dB as being representative of a typical TCVCXO found in COTS SDR systems. We compare this with the values given for a typical ultra-stable oscillator (USO) given in [142]. The changes in a 4 and a 3 with respect to [142] represent a lower quality oscillator and were chosen to t the data that have been measured in the USRP E312 SDR. The corresponding phase noise PSDs are plotted in Fig. 110 6.1. For comparison we also show the PSDs of the Microsemi GPS-1000 oven controlled crystal oscillator (OCXO), GPS-2550 double oven controlled crystal oscillator (DOCXO), and GPS-3500 Rubidium miniature atomic clock with a DOCXO low noise post lter crystal. 6.2.1 Simulation of Clock Phase Time Series As described by Krieger in [142], it is most useful to consider the one-sided PSD S ;osc (f) related to the two-sided PSD as S ;osc (f) = 8 > > < > > : 2S TS ;osc (f); f > 0 0; f < 0 (6.3) The phase noise PSDS ;osc (f) is plotted in Fig. 6.1. We generate time series realizations by dening a length N random spectral sequence z(f) sampled at f s; where z(f) = q f s; NS ;osc (f)e j 2(u;f) (6.4) where (u;f) is the realization at frequency f of a random variable uniformly distributed over the interval [0; 1]. Then the random time series clock phase realization r(t) is r(t) =F 1 fz(f)g [rad] (6.5) whereF 1 fg denotes the inverse Fourier transform (IFT). -0.5 0 0.5 Clk phase error (rad) 0 5 10 15 20 25 t (s) -10 -5 0 5 10 Clk time error (ns) (a) Over 30 seconds -1.5 -1 -0.5 0 0.5 1 1.5 Clk phase error (rad) 0 50 100 150 200 250 t (s) -20 -10 0 10 20 Clk time error (ns) (b) Over 300 seconds Figure 6.2: Multiple time domain realizations of random clock phase error with lower quality TCVCXO PSD from Fig. 6.1. Here we assume an oscillator frequency of 10 MHz Simulated time sequence realizations are shown in Fig. 6.2. 111 6.3 Modeling Synchronized clock PSD Consider that we start with input PSD S (f) of a wide sense stationary (WSS) random process, which has autocorrelation function R() and represents the phase of device clock phase (t) derived from an independent free-running oscillator. We examine the eect of performing the distributed N-sensor synchronization presented in Chapter 5 and in [32]. In this scheme, network synchronization is performed in three steps over a period referred to as the synchronization epoch, which is repeated once every synchronization repetition interval (SRI). In step one, each sensor broadcasts a synchronization waveform in a time-division-multiple-access (TDMA) fashion. Each sensor estimates the sub-sample broadcast signal TOF using the fractional peak estimation methods described in [32], which approach the CRLB for the given waveform and channel SNR. In stage two, each sensor encodes and broadcasts its ToF estimates along with a second synchronization waveform, which is received by all sensors. In stage three, each sensor computes an identical synchronized ToF matrix, a clock phase error matrix, and a RF carrier phase error matrix for the entire network, using the globally shared information. Each sensor is then able to pre-correct its own transmit signal to account for the phase noise of its now characterized local clock [32]. At some time T p after the most recent synchronization epoch, each sensor transmits a radar signal, which is now synchronized in the global domain. Because the clock phase errors are known globally, the network mean clock phase is used as the `master' reference to which all sensors synchronize their signals to, making this approach fully decentralized and hierarchy-free. For a given input oscillator phase error PSD (dened by power-law coecients that depend on the physical characteristics of the oscillator), the following parameters, along with the sync waveform parameters, determine the `synchronization process transfer function' that acts on the input PSD to produce a synchronized output PSD: 1. the pair-wise sensor SNRs 2. the synchronization repetition interval SRI 3. the time delay between the completion of the synchronization epoch and the radar pulse transmission T p 4. the number of sensors in the network N s 112 The eect of performing the distributed N-sensor synchronization with a nite delay T p between the sync epoch and transmission of the next radar pulse can be expressed for a time series realization of the oscillator phase error, (t) as resulting in a `synchronized' oscillator with phase error ~ (t). Knowing that the ToF estimation process can be treated as a mean-zero AWGN process n(t) with variance 2 C as determined by the zero-delay synchronization CRLB [32], we provide an expression for the synchronized oscillator phase error ~ (t) ~ (t) =(t)(t T p ) +n(t) (6.6) The autocorrelation function of this new random process ~ R(t) is then ~ R(t) = [(t)(t T p ) +n(t)] (6.7) [ (t) (T p t) +n (t)] We dene R n;n (t),n(t)n (t), R n; (t),n(t) (t), and R ;n (t),(t)n (t). Assuming that(t) andn(t) are uncorrelated and n(t) is zero-mean, i.e., R n; (t) =R ;n (t) = 0 this is reduced to ~ R(t) = 2R(t)R(t T p )R(t + T p ) +R n;n (t) (6.8) Therefore, using the properties of the Fourier transform, the synchronized PSD ~ S (f) is then ~ S (f) = ~ S 0 (f) +S n (f) (6.9) where ~ S 0 (f) = 2S (f)S (f)e i2fTp S (f)e i2fTp = 2S (f)(1 cos(2fT p )) S n (f) 2 C Simulation results are compared with the analytical expression in (6.9) in Fig 6.3. 113 Figure 6.3: Analytical synchronized output PSDs as a function of T p with a sinc window due to linear interpolation. Here SRI = 1 s and C = 1 ps. 6.3.1 Incorporating the SRI The synchronization process is eectively a downsampling of the clock error PSD by a factor of M =f s; =(1=SRI) where f s; is the sampling rate of the PSD process. If the SRI is insuciently short, the eects of PSD spectrum aliasing cannot be neglected. We express downsampling of the clock phase error by an integer factor M as #M [n] =[nM] (6.10) This can be expressed in the frequency domain via properties of the discrete fourier transform (DFT) as ~ S ;#M (f) = 1 M M1 X i=0 ~ S 0 f if s; M f s; (6.11) 114 Incorporating the eects of aliasing due to an insuciently short SRI, and assuming linear interpo- lation between samples taken every SRI seconds, which has the transfer function sinc(fSRI), we arrive at the following form for the synchronized PSD ~ S sync ;#M (f) = ( ~ S ;#M (f) +S(n))j sinc(fSRI)j 2 (6.12) where sinc(x) = sin(x)=x. Alternatively, in the case of sinc (rather than a linear) interpolation of the clock phase errors, the sinc term in (6.12) is replaced with rect(fSRI=2). 6.3.2 Prediction of Relative Clock Phase Error at Time of Radar Pulse If we assume that the relative clock phase errors are approximately linear over short time scales, as can be seen in Fig. 5.1, we may improve the performance of the synchronization scheme by predicting the clock phase error at some time T p after the last synchronization epoch when the radar pulse is transmitted. We predict the clock phase error at time t by estimating the linear LS slope of the previous N h measurements. We note that the appropriate selection of N h is highly dependent on the SRI and the time scale over which a linear approximation of the phase error is valid for a given oscillator. t = N h P i iSRI i (t) P i iSRI P i i (t) N h P i (iSRI) 2 ( P i iSRI) 2 (6.13) (6.14) where i (t) represents the previously measured clock phase normalized to the oldest measured value i (t) = (t (N h 1i)SRI T p ) (t (N h 1)SRI T p ) (6.15) (t) =(t) +n(t) (6.16) If we normalize the above expression, to the oldest measured values and times, we obtain t = 1 C (N h X i i i (t) X i i X i i (t)) (6.17) 115 where C =SRI N h X i i 2 X i i 2 (6.18) is constant for a given lter size. Noting that N h 1 X i i = N h (N h 1) 2 and N h 1 X i i 2 = N h (N h 1)(2N h 1) 6 (6.19) (6.20) then C =SRI N 2 h (N 2 h 1) 12 (6.21) and t = 12 SRIN 2 h (N 2 h 1) N h X i i i (t) X i i X i i (t) (6.22) Thus, the given the measured value (t T p ) and the linear LS prediction, resulting output clock phase error ~ (t) at the time of radar pulse transmission is ~ (t) =(t) (t T p ) + t T p (6.23) We arrive at the following form for the output synchronized clock error PSD ~ S (f) =S (f) h (2 2 cos(2fT p )) 2T p C N h X i (i N h 1 2 ) cos(2f((N h 1i)SRI + T p )) cos(2f(N h 1i)SRI) + T p C 2 N 2 h X i X j (ij + (N h 1) 2 4 ) exp(j 2f(ij)SRI) N 2 h (N h 1) X i X j i cos(2f(ij)SRI) i +S n (f) h 1 T p C N h (N h 1) + T p C 2 N 3 h (N 2 h 1) 12 i (6.24) 116 A derivation of this expression is given in Appendix D. 6.4 Synchronization Simulation We simulate the synchronization performance for three sensors. The clock phase errors of each sensor is generated from the PSD as described in Section 6.2 as a time-series. The two way time of ight CRLB for a linear frequency modulated (LFM) chirp synchronization waveform, which is derived in Appendix B, is given as [32] 2 2TOF = 3 (2B) 2 SNR (N=f s )f s (6.25) where B is the signal bandwidth, SNR is the signal to noise ration, N is the number of waveform samples, and f s is the sampling rate. Note that N=f s is the pulse duration. We derive the clock phase error CRLB after synchronization by N s sensors to the network mean clock oset. If we include the AWGN measurement noise n i;j which has variance corresponding to the one-way TOF CRLB 2 1TOF = 2 2 2TOF , the relative clock phase error and corresponding correction terms given in Chapter 5 and in [32] have the form ~ i;j = i j + n i;j n j;i 2 (6.26) ~ i = 1 N X j ~ i;j (6.27) = i 1 N X j j + 1 N X j n i;j n j;i 2 (6.28) We compute the variance of this term as Var[ ~ i ] =E[( 1 N X j n i;j n j;i 2 ) 2 ] (6.29) Noting that E[n 2 i;i ] = 0 because ~ is zero along the diagonal, the clock phase error CRLB after synchronization by N s sensors to the network mean clock oset has the form: 2 clk = N s 1 N 2 s 2 2TOF (6.30) 117 6.4.1 TDMA Distributed Sensor Simulation We use the oscillator PSD to generate time-series clock phase errors for each sensor as described in Section 6.2. We note that we atten the PSD for f <:01 Hz in order to control runaway DC bias due to the singularity at f = 0. Using the input clock phases, we simulate a complete distributed radar sensor network, where synchronization is performed as described in Chapter 5 and [32], and used to synchronize each sensors radar signals, which are transmitted at a nite time T p after each synchronization epoch. (a) Tp = 10 ms (b) Tp = 400 ms Figure 6.4: Simulated clock phase error PSDs after synchronization. SNR = 30 dB, SRI =:5 s. (a) Input clock phase error (b) Synchronized clock phase error Figure 6.5: Simulated clock phase error of each sensor before and after synchronization at time of radar transmission (T p = 10 ms). SNR = 30 dB, SRI =:5 s. From simulation, we have observed that when T p is suciently small, the output PSD looks like a band-limited white noise PSD with cuto frequency f c = 1=SRI. On the interval 0<f <f c , the PSD will exhibit suppressed characteristics of the input PSD until the suppression falls beneath the white noise power oor as determined by the SNR, at which point the PSD becomes whitened. 118 This is illustrated for a few values of T p in Fig. 6.4, which shows the synchronized eective output clock PSDs for each sensor. Time series plots of the input clock phase errors for each of the three sensors simulated are shown in Fig. 6.5a, the residual clock phase errors at the radar transmission time T p = 10 ms after the synchronization epoch are shown in Fig. 6.5b. (a) Tp = 10 ms, N h = 2 (b) Tp = 400 ms, N h = 2 Figure 6.6: Simulated clock phase error PSDs after synchronization using the LS prediction method described in Section 6.3.2 from N h = 2 previous clock phase error estimates. SNR = 30 dB, SRI =:5 s. (a) No LS prediction of clock state (b) LS prediction of clock state for N h = 2 Figure 6.7: Simulated clock phase error of each sensor with and without LS prediction for N h = 2 (T p = 400 ms). SNR = 30 dB, SRI =:5 s. We now perform the simulation using the LS prediction method described in Section 6.3.2 for the same values of T p . We use the N h = 2 previous clock phase error measurements made by the synchronization algorithm to estimate the clock phase error states at time T p after the most recent synchronization epoch. As shown in Fig 6.6, this method signicantly reduces the synchronized clock error PSD at lower spectral phase components, improving the performance at the the time of radar pulse transmission. The PSD produced using the expression derived in (6.24) is shown for comparison with the simulated results. 119 A comparison of the synchronized clock phase error time-series for T p with and without the LS clock phase prediction method is shown in Fig. 6.7. 6.5 Experimental Results 10 -2 10 0 10 2 f (Hz) -100 -80 -60 -40 -20 0 20 40 60 dB TOF radar PSDs " T p =520 ms " T p =120 ms " T p =20 ms input analytic Figure 6.8: Experimentally obtained PSDs of non-synchronized TOF (input PSD) and synchronized radar pulse TOF delay T p for 3 sensors over 1000 pulses at SRI=1 s. Here we estimated C = 20 ps according to (6.25). To validate both the analytical formulation and the simulation, we conducted an experiment with three USRP E312 SDRs that implement the distributed wireless synchronization scheme on-board in an embedded Xilinx Zynq System on Chip (SoC) [32] as described in Chapters 3 and 5. We run the 3-sensor synchronization for 1000 pulses at SRI=1 s. We compare the experimental synchronization results with the proposed analytical expressions given in Section 6.3, which were applied to with the measured non-synchronized `input' PSDs in Fig. 6.8 for selected values of T p . As shown, the experimentally measured values agree with the the proposed analytic model. 120 6.6 Distributed MIMO SAR Simulation Using the model described previously in this chapter, we may eciently simulate large scale distributed coherent radar systems that implement the synchronization scheme described in Chapter 5. We generate time-series realizations of synchronized clock phase errors from the derived PSDs for all sensors in the network given the set of synchronization parameters and the input PSDs as dened by characteristics of the oscillators we wish to model. (a) Simulated Scene (b) Unfocused 3x3 Radar Data Figure 6.9: Simulated 3x3 MIMO SAR scene and unfocused radar data following range compression via matched lter. Note that unfocused data is shown here for the ideal (perfectly synchronized) case. We simulate a 3x3 distributed MIMO SAR network and perform MIMO coherent MIMO SAR focusing to demonstrate the performance eects of synchronization or lack thereof. The simulated scene is shown in Fig. 6.9a, which consists of 9 distinct point targets. The sensor network is arranged along a linear aperture at a height of 20 m above ground level (AGL). Each radar transmits to and receives from all other sensors in the distributed array. The raw unfocused radar data is shown in Fig. 6.9b for the perfectly synchronized case. The columns represent the transmitting sensor and the rows represent the receiving sensor. Note that along the diagonal each sensor is receiving its own transmitted signal. In this simulation, synchronization is performed using 50 MHz of bandwidth and the radar operates with 200 MHz of bandwidth at a center frequency of 500 MHz. The aperature is sampled every 5 cm and has a total length 40 m. 121 6.6.1 Full MIMO SAR Focusing Using time domain back-projection each transmit-receive pair SAR image is focused with respect to a global coordinate system that is dened with respect to the network center. The pair-wise focused SAR images are shown in Fig. 6.10a for GPS synchronization only ( 4:4 ns clock error and 14 rad phase error). The pair-wise focused SAR images are shown in Fig. 6.10b synchronization performed over a 0 dB SNR channel with 50 MHz of bandwidth ( 80 ps clock error and:25 rad phase error). (a) GPS sync only (Clk err 4400 ps, phase err 13:9 rad) (b) Sync at 0 dB SNR (Clk err 80 ps, phase err 0:25 rad) Figure 6.10: Simulated 3x3 SAR focusing for all pairs of distributed radars with independent oscillators. Synchronization is performed over a channel with a given SNR and 50 MHz of bandwidth. MIMO SAR focusing is performed using time-domain backprojection. The pair-wise SAR focused images are accumulated into a single MIMO SAR image shown in Fig. 6.11a for the GPS-synchronization case and Fig. 6.11b for the synchronized case. Cross sections of the azimuth target response for both cases are are shown in Fig. 6.12, where the solid black line represents the focused MIMO SAR result. Note that MIMO SAR focusing yields signicant reduction in both SNR and azimuth sidelobe level. With GPS synchronization only, the targets are still resolvable, albeit at signicantly reduced SNR. This is due to the inclusion of along diagonal SAR data wherein each sensor is transmitting to itself and is therefore fully synchronized. 122 3x3 MIMO backprojection 20 40 60 80 100 Ground Range (m) -15 -10 -5 0 5 10 15 Azimuth (m) -30 -20 -10 0 10 20 dB (a) 3x3 MIMO SAR image with GPS sync only (Clk err 4400 ps, phase err 13:9 rad) 3x3 MIMO backprojection 20 40 60 80 100 Ground Range (m) -15 -10 -5 0 5 10 15 Azimuth (m) -20 -10 0 10 20 30 dB (b) 3x3 MIMO SAR image with sync at 0 dB SNR (Clk err 80 ps, phase err 0:25 rad) Figure 6.11: Simulated 3x3 MIMO SAR focusing for distributed radars with independent oscillators. Synchronization is performed over a channel with a given SNR and 50 MHz of bandwidth. MIMO SAR focusing is performed using time-domain backprojection. -20 -15 -10 -5 0 5 10 15 20 Azimuth (m) -100 -50 0 Amplitude (dB) target at 30 m -20 -15 -10 -5 0 5 10 15 20 Azimuth (m) -80 -60 -40 -20 0 Amplitude (dB) target at 56.5685 m -20 -15 -10 -5 0 5 10 15 20 Azimuth (m) -80 -60 -40 -20 0 Amplitude (dB) target at 70 m (a) GPS sync only (Clk err 4400 ps, phase err 13:9 rad) -20 -15 -10 -5 0 5 10 15 20 Azimuth (m) -100 -50 0 Amplitude (dB) target at 30 m -20 -15 -10 -5 0 5 10 15 20 Azimuth (m) -80 -60 -40 -20 0 Amplitude (dB) target at 56.5685 m -20 -15 -10 -5 0 5 10 15 20 Azimuth (m) -80 -60 -40 -20 0 Amplitude (dB) target at 70 m (b) Sync at 0 dB SNR (Clk err 80 ps, phase err 0:25 rad) Figure 6.12: Focused 3x3 MIMO SAR target response in azimuth. Note that the monostatic cases are included in the MIMO packprojection, which results in a degree of focusing for the non-synchronized case. 6.6.2 Pair-only MIMO SAR Focusing We now provide results for Nx(N-1) MIMO SAR focusing where sensors do not receive the signal they have transmitted. That is, only the bistatic pairs of signals are considered. This is an interesting use case for distributed MIMO SAR systems as it allows for the transmitter power of each sensor to be increased, thus improving system-wide SNR, without risking damage to the colocated receiver. 123 3x(3-1) MIMO backprojection 20 40 60 80 100 Ground Range (m) -15 -10 -5 0 5 10 15 Azimuth (m) -40 -30 -20 -10 0 10 dB (a) 3x2 (bistatic pairs only) MIMO SAR image with GPS sync only (Clk err 4400 ps, phase err 13:9 rad) 3x(3-1) MIMO backprojection 20 40 60 80 100 Ground Range (m) -15 -10 -5 0 5 10 15 Azimuth (m) -20 -10 0 10 20 30 dB (b) 3x2 (bistatic pairs only) MIMO SAR image with sync at 0 dB SNR (Clk err 80 ps, phase err 0:25 rad) Figure 6.13: Simulated 3x2 MIMO SAR focusing for distributed radars with independent oscillators. Synchronization is performed over a channel with a given SNR and 50 MHz of bandwidth. MIMO SAR focusing is performed using time-domain backprojection. The focused MIMO SAR results for this case are shown in Fig. 6.13. Now, the GPS- synchronization is shown to be insucient for successful coherent MIMO SAR focusing and the targets are not resolved as seen in Fig. 6.13a. With synchronization at 0 dB SNR, the MIMO SAR system is coherent and can be successfully focused as seen in Fig. 6.13b. -20 -15 -10 -5 0 5 10 15 20 Azimuth (m) -100 -50 0 Amplitude (dB) target at 30 m -20 -15 -10 -5 0 5 10 15 20 Azimuth (m) -80 -60 -40 -20 0 Amplitude (dB) target at 56.5685 m -20 -15 -10 -5 0 5 10 15 20 Azimuth (m) -80 -60 -40 -20 0 Amplitude (dB) target at 70 m (a) GPS sync only (Clk err 4400 ps, phase err 13:9 rad) -20 -15 -10 -5 0 5 10 15 20 Azimuth (m) -100 -50 0 Amplitude (dB) target at 30 m -20 -15 -10 -5 0 5 10 15 20 Azimuth (m) -80 -60 -40 -20 0 Amplitude (dB) target at 56.5685 m -20 -15 -10 -5 0 5 10 15 20 Azimuth (m) -80 -60 -40 -20 0 Amplitude (dB) target at 70 m (b) Sync at 0 dB SNR (Clk err 80 ps, phase err 0:25 rad) Figure 6.14: Focused 3x2 MIMO SAR target response in azimuth when only bistatic pairs are considered (i.e. no sensor receives the signal it transmitted). 124 Cross sections of the azimuth target response for both cases are are shown in Fig. 6.12. Again, that synchronized MIMO SAR focusing yields signicant reduction in both SNR and azimuth sidelobe level, while the targets are not resolvable in azimuth with GPS-only synchronization in this case. 6.7 Conclusion In this work, we have derived an analytical formulation for the PSD transfer function of a previously reported N-sensor distributed wireless synchronization scheme and have validated the model both in simulation and experimentally. The derived expressions relate the input independent sensor clock phase PSDs to output synchronized clock PSDs in terms of the number of sensors, sync repetition interval, time delay between the sync epoch and signal transmission, and the sensor ToF estimation and synchronized clock error CRLBs. This enables accurate modeling and simulation of large distributed radar sensor networks, such as smallSAT/cubeSAT constellations and sUAS swarms, at signicantly improved eciency using the synchronized output PSDs directly to generate time-series realizations of independent sensor clock phase errors following network synchronization via the method described in Chapter 5. Furthermore, the model allows for complete examination of the trade-space between synchronization scheme parameters and oscillator selection enabling optimization of system size, weight, and power (SWaP). 125 Chapter 7 Snow Depth Retrieval with an Autonomous UAV-mounted Software Dened Radar The results of this chapter will be published in [60]. 7.1 Introduction Mapping of the water stored in snowpacks or SWE at high resolution on a global scale is critical to gaining a complete understanding of both surface and subsurface hydrologic processes [82]. Due to the high variability of snow distribution and dynamics across diverse landscapes, there is a critical need for the emergence of technologies that are not only capable of high spatial resolution but also deployable over large spatial extents and across varied and rugged terrain [75]. In mountainous environments of the western United States (U.S.) seasonal snowmelt runo is a major source for stream ow and subsurface recharge across much of the region [145]. Therefore, quantifying and tracking the amount of water stored within seasonal snow is critically important for water resources forecasting and planning. The seasonal snowpack is often highly variable over space and through time in mountainous regions [146] as a result of the driving in uences of meteorology (e.g., precipitation, temperature, vapor pressure, wind, and net radiation) on the snowpack and its interactions with terrain and forest canopy features [147]. These complex processes controlling snow evolution at ne spatial scales make monitoring the spatial variability of the seasonal snowpack a substantial challenge. 126 Satellite and airborne remote sensing platforms oer great potential for monitoring the spatiotem- poral variability of the seasonal snowpack that cannot be done by traditional eld measurements. However, there are still major challenges limiting the direct observation of snow distributions in mountainous and other environments. SWE estimation requires measurement of (i) snow cover extent, (ii) snow depth, and (iii) snow density. Remote sensing of SWE directly at high spatial resolution has not been achieved [87]. The NASA SnowEx Science Plan [82] highlights these challenges associated with remote sensing of snowpack characteristics (including SWE), across the globe, but also outlines various new snow measurement techniques that show promise for improved snow observations. Of these new measurement techniques, radar is listed as a promising technology for snowpack monitoring due its ability to penetrate vegetation and forest cover and the potential for direct SWE retrieval from volume scattering [82]. Furthermore, there has been much recent excitement surrounding the potential of small multicopter aircrafts as a viable path for widespread deployment of snow/ice imaging radar sensors, although the capabilities of such systems are as of yet largely unquantied [71]. Figure 7.1: Illustration of sUAS-based UAV-SDRadar system performing synthetic ultra-wideband radar measurements of snow depth while ying over snow covered meadow at Cameron Pass, Colorado. In this chapter, we demonstrate snow depth retrieval across variable terrain with an autonomous multicopter sUAS consisting of a SWW SDRadar sensor implemented in a COTS SDR, a DJI M600 hexacopter UAV, and ground control station. The UAV-mounted SDRadar (UAV-SDRadar) system 127 is shown in ight over the Meadow Transect in Fig. 7.1. We operate the UAV-SDRadar system, shown on the ground in Fig. 7.2a, over both meadow and forested transects and validate the results with measurements taken on the ground by a sled-based GPR system as seen in Fig. 7.2b. The objectives of this research are to (i) demonstrate snow depth retrieval across varied terrain from a sUAS-based radar system validated with ground-based GPR measurements, (ii) develop radar processing and surface detection algorithms specic to sUAS-based radar platforms to further enable autonomous operation and measurement, (iii) demonstrate the use of widely available open-source COTS SDR hardware for high-resolution radar applications from sUAS platforms using synthetic wideband techniques, and (iv) to quantify performance of sUAS-based radar in varied terrain including both open and forested regions, the latter of which is identied as the rst of seven critical gaps in current remote sensing technologies in the NASA SnowEx Science Plan [82]. (a) UAV-SDRadar sUAS system consisting of a battery powered USRP E312 SDR and two Vivaldi TX/RX antennas mounted on a DJI M600 hexacopter UAV. (b) UAV-SDRadar sUAS system in ight over the Meadow Transect with sled GPR team below. Figure 7.2: UAV-SDRadar system and sled GPR team operating at Cameron Pass, CO test site. The snow depth along the transect was measured by the UAV-SDRadar system as well as by a ground-based sled GPR operated by the team as shown. Note that the landing gear of the UAV retracts in ight. 128 7.2 UAV-SDRadar Payload Integration Table 7.1: Size and Weight of UAV-SDRadar Peripheral Sensors and Modules Module Size Weight DJI Manifold 2 Jetson TX2 SBC 91 x 61 x 35 mm 230 g Emlid M2 RTK GNSS/GPS 56.4 x 45.3 x 14.6 mm 35 g SF11/C lidar altimeter 30 x 56.5 x 50 mm 35 g RFSpace TSA600 Antenna(s) 240 x 330 x 1.5 mm 227 g RF Amp (Minicircuits ZX60-V82-S+) 19.1 x 18.8 x 11.7 mm 23 g USRP E312 (SDRadar) RX-B TX-A RX-A DJI M600 (UAV Aircraft) WLAN SF11/C (Lidar Altimeter) Emlid M2 (GNSS/ GPS) GPS IMU SD Card SDRadar Server Software Flight Control Tablet DroneAmplified Flight Software -Mission control -Camera feed -Radargram feed RF Link DJI Remote Controller SDRadar Control Laptop SDRadar Client C++ Control/Processing Software Ethernet (UDP Data stream for real-time visualization) DJI Radio Link (Control, Telemetry, Video) DJI Manifold 2 (Jetson TX2 Linux SBC) SDRadar Client C++ Control/Processing Software HDMI Video Out (optical imagery) Camera (Situational Awareness) (Only used pre-flight to initiate radar collection and post-flight to download data) HDMI Video Out (real-time radargram) -20 dB TX RX -3 dB Splitter PA Aircraft Ground Station Figure 7.3: Diagram of sUAS system used for ight experiments. The UAV-SDRadar system is shown on the left and the ground control system is shown on the right. Real-time telemetry and (optional) processed radargram live video feed is provided via high power DJI Radio Link. SDRadar conguration, mission control, and data download are done via WiFi Wireless Local Area Network (WLAN). The optional Jetson TX2 SBC (DJI Manifold 2) can process SWW radar data in real-time and relay a live radargram video feed to the ground station. The radar data are also stored locally on the E312 SDRadar SD card for oine processing. 129 The SDRadar system described in Chapter 3 is integrated as a payload on the UAV. A DJI Manifold 2 single board computer (SBC), based on the NVIDIA Jetson TX2 GPU, is also own on the aircraft. The onboard Manifold 2 SBC runs the SDRadar Client C++ software and receives and processes real-time radar data from the UDP stream to form SWW radargrams which are transmitted as an HDMI video stream over the DJI M600 RF Link to the ground station. The complete UAV-SDRadar payload is own on a DJI Matrice M600 Pro hexacopter UAV. The M600 UAV is own autonomously on pre-dened ight paths using control software from Drone Amplied. The UAV runs on six 22.2V 4500 mAh DJI TB47S batteries, with an estimated power consumption of the UAV is 2kW. The maximum payload capability of the UAV is 6 kg. The UAV ight time is 32 minutes with no payload and 16 minutes with a 6 kg payload it at sea level [148]. The elevation at the test site was 9600 ft (2600 m). At this altitude with the 1:25 kg radar payload the UAV ight time was 18 minutes. The radar payload runs o of the internal USRP E312 3.7 V 3200 mAh battery. The radar payload power consumption is 2-6 W, with the battery lasting 2-3 hours in normal operation. Details on the sizes/weights of the peripheral sensors and modules are provided in Table 7.1 and the hardware diagram of the complete system used in this eort is shown in Fig. 7.3. The complete system is shown mounted to the underside of the UAV in Fig. 7.2a. An optical camera is also own for situational awareness while operating the aircraft. 7.3 SWW Reconstruction and Radar Processing The radar processing consists of the following steps: 1. Internal sub-pulse phase calibration 2. Sub-pulse system sky-calibration removal 3. NUFS to form the SWW (described in Section 4.3.3) 4. 2D SWW radargram construction 5. Row-wise complex average subtraction across sweeps to remove direct path signal sidelobes and further reduce grating lobes 130 NUFS SWW Reconstruction Apply NUFS algorithm to form SWW GLS Filter Internal sub-pulse phase error correction Sub-pulse system sky-calibration removal For f c = f 0 ,...,f N-1 For Scan i: i=0,...,N s -1 Radargram Processing Motion compensation using GPS/GNSS and SFM data or surface- detection-based elevation data Row-wise complex mean subtraction Column-wise range magnitude scaling Column-wise 200- point moving average power threshold filter Column-wise mean power subtraction Accumulate Scan i SWW as column in radargram Column-wise median power threshold filter Figure 7.4: Radar processing diagram showing SWW synthesis using the NUFS algorithm (left) and the radargram processing steps used to compensate UAV platform motion and to reduce background noise and speckle. 6. Motion compensation to remove eects of UAV altitude variation 7. Column-wise moving average threshold lter 8. Column-wise mean subtraction for each scan to reduce vertical striping 9. Column-wise Moving median threshold lter to reduce speckle The complete radar signal processing is visualized in Fig. 7.4. We note that the horizontal lines present in the raw radargram at the bottom left of Fig. 7.4 are due to a combination of physical ringing between the TX and RX antennas, the sidelobes of the direct path signal, and minor residual processing artifacts due to the NUFS SWW reconstruction. 131 Forest Transect (Joe Wright Site) Meadow Transect (Michigan River Site) Snow Pit UAV UAV Forest Transect (Joe Wright Site) Meadow Transect (Michigan River Site) Snow Pit Figure 7.5: Photos of Forest (top-left) and Meadow (bottom-left) snow eld transects and GPS positions of sled GPR and UAV-SDRadar along both transects (right) for March 2020 eld campaign at Cameron Pass, CO. The snowpack TWT measurements made by each system are illustrated. Note that the ight path of UAV deviates from the sled-GPR path by varying distances for both transects due to imperfect ight control systems and terrain obstacles. 7.4 Experimental Design and Results We performed a eld campaign at Cameron Pass, Colorado in March 2020 to evaluate the proposed SDRadar system against independent ground-based observations of snow depth. The SDRadar was mounted on a hexacopter UAV as a snow penetrating radar sensor payload capable of high resolution subsurface imaging using the NUFS synthetic wideband and signal processing techniques described in Section 7.3. Two transects were own: a snow covered meadow (\Meadow Transect") in the southern part of the test area and a forested area (\Forest Transect") in the northern part of the test area as shown in Fig. 7.5. A total bandwidth (BW) of 1 GHz is synthesized from the 600 1600 MHz band and a bandwidth eciency of 80% was empirically selected to minimize the number of frequency steps while still having sucient sub-pulse overlap for good NUFS reconstruction characteristics. Additionally 10 sub-pulses are averaged coherently at each step to improve SNR by 10 dB. The total N = 26 sweep time is 0:5 s and is linearly dependent on the number of frequency steps (:018 s/step). The UAV ight speed was 1 m/s, giving a horizontal resolution of 0:5 m and the nominal output power was 10 dBm. The radar sweep PRF is limited by the relatively slow frequency tuning time of the USRP E312 hardware as well as the 10 pulse averaging. In practice, 132 the UAV ight speed must be slow enough such that the scene remains coherent throughout the frequency sweep. Faster ight speeds may be achieved by reducing the total synthesized bandwidth, which is run-time congurable. As a means of validation, we measured the snow depth TWT using a sled equipped with a Sensors and Software commercial pulseEKKO GPR system with 1 GHz antennas and Emlid RS2 GNSS module [149]. The post-processed GNSS coordinates of the ground-based GPR and the UAV SDRadar for the two transects own are shown in Fig. 7.5. Also shown in Fig. 7.5 are the TWT measurements made by each system. The methodology used to obtain these results and their analysis are described in detail later in this Section. 7.4.1 Snow and Ground Surface Detection Algorithm We develop a custom snow and ground surface detection algorithm to track the radar backscatter signatures of both the air/snow interface and the snow/ground interface observed by the UAV- SDRadar. Of the existing retracking algorithms used to detect surface returns in altimetric systems, the cubic spline retracking algorithm is often cited as having the best performance [150] [151]. However, these one-dimensional approaches are designed for high altitude space and airborne systems with wide swaths. They are not sucient for small UAV subsurface imaging radar systems ying at low altitude, as the surface backscatter may vary signicantly in power. This is due to large uctuations in platform elevation relative to the platform height above the ground as dictated by the fundamental radar equation [61]. This is further complicated when ying over forested or densely vegetated areas, due to the sporadic presence of scattering objects very close to the radar sensor. As such, a one-dimensional approach to surface detection will easily lose track of the surface interface when the signal is weak or obstructed, producing erratic results. To this end, we developed a neighbor-aware approach to surface tracking that incorporates independent sensor data from a lidar altimeter and GNSS/GPS receiver as well as the expected continuity and gradient characteristics of the ground and snow surfaces over small spatial scales. The proposed surface detection algorithm is described completely in Appendix E. 133 The results presented herein for the Meadow Transect were obtained entirely using this detection algorithm. For the Forest Transect, which suers from signicant speckle and poor lidar data due to interference from the canopy, 80 manually selected anchor points were added prior to the second iteration to assist the algorithm. While easy for the human brain, the problems of noisy/broken edge detection and curve tracing represent an active area of research in the eld of computer vision and, while highly important, is not the primary focus of this manuscript [152], [153], [154]. Rather, we intend to focus on the performance potential of the described UAV-SDRadar system in measuring physical parameters. 7.4.2 UAV Elevation-Corrected Radargrams Using GNSS Following ultra-wideband synthesis and ltering of the stepped frequency SWW, we correct for the UAV altitude using data collected by an onboard Emlid Reach M2 RTK/ post-processing kinematic (PPK) L1/L2 GNSS/GPS unit. Each column of the SWW radargrams shown in Fig. 7.6 and 7.7 is fractionally shifted via application of a linear phase in the spectral domain in accordance with fundamental Fourier transform theory to compensate for the UAV elevation. With the UAV platform elevation removed, the radargrams shown in Fig. 7.6b and 7.7b are representative of the snow/ground surface topography present at the Meadow and Forested transects, respectively. We note that in the Forested Transect, trees were present throughout and scattering from the forest canopy is responsible for the increased speckle above the snow surface in Fig. 7.7. Additionally, a large tree, which forced the UAV to abruptly gain altitude, is responsible for the feature at 360 m in Fig. 7.7. UAV-based SfM-MVS surveys were conducted at the eld study site with ground control points in the Meadow Transect area as part of an ongoing independent eort. The SfM-MVS surveys were conducted with a DJI Mavic 2 Pro drone own at 60 m AGL with 75% front and side overlap. Drone imagery was processed in Agisoft Metashape using a work ow comparable to [89]. For the meadow transect, we use the 1 cm/pixel resolution SfM snow surface elevation maps generated from the multi-angle SfM imagery of the eld site to correct for surface topography. The SfM-derived surface elevation data as well as the Kalman ltered elevation data obtained from the onboard lidar altimeter are plotted for the Meadow Transect in Fig. 7.6a. The SfM elevation data were retrieved from a geo-referenced map using the RTK GNSS coordinates of the 134 0 20 40 60 80 100 Time delay (ns) 50 100 150 200 250 300 350 Distance along transect (m) 0 5 10 15 Free space range (m) Kalman filtered lidar elevation data RTK+SFM elevation data (a) Radargram after RTK GNSS/GPS UAV elevation correction showing the meadow snow surface topography. Also shown are the Kalman ltered lidar altimeter data and independent 1 cm/pixel SfM-MVS-derived elevation data, both corrected for RTK GNSS/GPS UAV elevation. 0 20 40 60 80 100 Time delay (ns) 50 100 150 200 250 300 350 Distance along transect (m) 0 5 10 15 Free space range (m) (b) Before surface topography removal 0 10 20 30 Time delay (ns) 50 100 150 200 250 300 350 Distance along transect (m) -1 0 1 2 3 4 5 Free space range (m) (c) After surface topography removal (using independent SfM data) Figure 7.6: Meadow Transect radargrams before and after surface topography removal. Detected snow surface and ground layers are colorized in the grayscale radargram images with blue tones indicating weaker re ections and red tones indicating stronger re ections. The ground and surface detections are produced by the algorithm described in Section 7.4.1 UAV platform. For each UAV position, all SfM data within a 0.5 m radius footprint were considered. We also averaged the closest values within a 5 cm radius footprint. The values shown were the average of the minimum 10% values within the 0.5 m footprint. This was done to isolate points corresponding to the snow surface rather than vegetation or other objects. For the Meadow Transect, both the lidar altimeter and SfM map provide good estimates of the snow surface elevation. For the forest transect the lidar altimeter data, even after extensive processing and Kalman ltering, is contaminated signicantly by the dense foliage and forest canopy as seen in Fig. 7.7a. SfM-based 135 0 20 40 60 80 100 Time delay (ns) 50 100 150 200 250 300 350 Distance along transect (m) 0 5 10 15 Free space range (m) Kalman filtered lidar elevation data (a) Radargram after RTK GNSS/GPS UAV elevation correction showing the forest snow surface topography. Also shown is the Kalman ltered lidar altimeter data corrected for RTK GNSS/GPS UAV elevation. 0 20 40 60 80 100 Time delay (ns) 50 100 150 200 250 300 350 Distance along transect (m) 0 5 10 15 Free space range (m) (b) Before surface topography removal 0 10 20 30 Time delay (ns) 50 100 150 200 250 300 350 Distance along transect (m) -1 0 1 2 3 4 5 Free space range (m) (c) After surface topography removal (using surface detec- tion algorithm) Figure 7.7: Forest Transect radargrams before and after surface topography removal. Detected snow surface and ground layers are colorized in the grayscale radargram images with blue tones indicating weaker re ections and red tones indicating stronger re ections. The ground and surface detections are produced by the algorithm described in Section 7.4.1 surface elevation data were only provided for the Meadow Transect, as the primary study area of the SfM survey did not cover the Forest Transect. Furthermore, conducting SfM surveys of heavily forested areas presents signicant challenges due to canopy opaqueness at optical frequencies [65]. However, the RF signal transmitted by the UAV-SDRadar is able to eectively penetrate the forest canopy, with backscatter from both the snow surface and ground below clearly detected. Accordingly, for the Meadow Transect, where the independent elevation sensors provide reliable data, we correct for surface topography using SfM data as shown in Fig. 7.6c. For the forest transect, however, we correct snow surface topography directly using the output of the surface estimation 136 algorithm described in Section 7.4.1 as shown in Fig. 7.7c. In both cases, the TWT results are derived from the ground and snow surface estimated by the proposed surface estimation algorithm. Non-colorized versions of the radargrams can be found in Appendix F. 7.4.3 Comparison with Ground-Based GPR Using the the algorithm described in Section 7.4.1, we obtain the UAV-SDRadar TWT measurements by taking the dierence of the ground and snow surface detections and compare them with the ground-based commercial GPR measurements. As shown in Fig. 7.5, the ight path of UAV deviates from the sled-GPR path by varying amounts for both transects. Furthermore, dierences in spatial sampling rate and radar footprint size necessitate a normalization of the datasets for meaningful comparison. 0 100 200 300 Distance along transect (m) 5 10 15 TWT (ns) UAV GPR (FP=20 m) GPR (FP=10 m) GPR (FP=5 m) GPR (FP=2 m) (a) Meadow Transect 0 100 200 300 Distance along transect (m) 5 10 15 TWT (ns) UAV GPR (FP=20 m) GPR (FP=10 m) GPR (FP=5 m) GPR (FP=2 m) (b) Forest Transect Figure 7.8: Along-prole comparison of GPR and UAV-SDRadar two way travel time (TWT) results. The weighted average is taken over all nearest neighbor GPR measurements falling within a given footprint (FP) size centered around each UAV lat/lon GNSS coordinate as described in (7.1). The FP size refers to the diameter of a circle on the ground surface. (See also Fig. 7.9). Given the set of GPR TWT measurements, T gpr j where the GPS/GNSS-derived position of the GPR isp gpr j = [x gpr j ;y gpr j ; 0] for GPR scan j, we compute a weighted sum of all GPR measurements falling within a given diameter footprint (FP) of the UAV-SDRadar at the GPS/GNSS-derived platform position p uav i = [x uav i ;y uav i ;z uav i ] for each UAV-SDRadar scan i. Based on the radar received power, which, for a point target at range R is proportional to 1=R 4 , we arrive at the following weighting scheme to compute the weighted average of GPR TWT measurements in the i th UAV-SDRadar footprint, T gpr i . 137 ! i;j = z uav i jjp uav i p gpr j jj l2 ! 4 8 j : q (x uav i x gpr j ) 2 + (y uav i y gpr j ) 2 FP=2 (7.1) T gpr i = P j (! i;j T gpr j ) P j ! i;j (7.2) 0 5 10 15 20 GPR TWT (ns) 0 5 10 15 20 UAV TWT (ns) FP = 20 m ; = 0.89 RMSE=0.86ns 0 5 10 15 20 GPR TWT (ns) 0 5 10 15 20 UAV TWT (ns) FP = 10 m ; = 0.90 RMSE=0.93ns 0 5 10 15 20 GPR TWT (ns) 0 5 10 15 20 UAV TWT (ns) FP = 5 m ; = 0.90 RMSE=0.87ns 0 5 10 15 20 GPR TWT (ns) 0 5 10 15 20 UAV TWT (ns) FP = 2 m ; = 0.85 RMSE=0.93ns (a) Meadow Transect 0 5 10 15 20 GPR TWT (ns) 0 5 10 15 20 UAV TWT (ns) FP = 20 m ; = 0.70 RMSE=1.40ns 0 5 10 15 20 GPR TWT (ns) 0 5 10 15 20 UAV TWT (ns) FP = 10 m ; = 0.74 RMSE=1.35ns 0 5 10 15 20 GPR TWT (ns) 0 5 10 15 20 UAV TWT (ns) FP = 5 m ; = 0.68 RMSE=1.43ns 0 5 10 15 20 GPR TWT (ns) 0 5 10 15 20 UAV TWT (ns) FP = 2 m ; = 0.64 RMSE=1.38ns (b) Forest Transect Figure 7.9: Scatter plot comparison of GPR and UAV-SDRadar TWT results for dierent footprint (FP) sizes. The weighted average is taken over all nearest neighbor GPR measurements falling within a given FP size centered around each UAV lat/lon GNSS coordinate as described in (7.1). The FP size refers to the diameter of a circle on the ground surface. (See also Fig. 7.8). 0 100 200 300 Distance along transect (m) 5 10 15 TWT (ns) UAV GPR (dist<15 m): ;=0.88, RMSE=1.00 ns 0 5 10 15 Distance from UAV (m) (a) Meadow Transect 0 100 200 300 Distance along transect (m) 5 10 15 TWT (ns) UAV GPR (dist<15 m): ;=0.61, RMSE=1.53 ns 0 5 10 15 Distance from UAV (m) (b) Forest Transect Figure 7.10: Along-prole comparison of UAV-SDRadar TWT measurements with the nearest GPR measurement for a given UAV lat/lon GNSS coordinate. The distance of the closest GPR data point from the corresponding UAV data point is indicated by line color. GPR data points that are further than 15 m away from the UAV position are not considered. (See also Fig. 7.11). 138 0 5 10 15 20 GPR TWT (ns) 0 5 10 15 20 UAV TWT (ns) ; = 0.88 RMSE=1.00 ns 0 5 10 15 Distance from UAV (m) (a) Meadow Transect 0 5 10 15 20 GPR TWT (ns) 0 5 10 15 20 UAV TWT (ns) ; = 0.61 RMSE=1.53 ns 0 5 10 15 Distance from UAV (m) (b) Forest Transect Figure 7.11: Scatter plot comparison of GPR TWT measurement closest to the UAV for a given position. The distance of the closest GPR data point from the corresponding UAV data point is indicated by color. GPR data points that are further than 15 m away from the UAV position are not considered. (See also Fig. 7.10). Using the weighting scheme described in (7.1), we compare the GPR and UAV-SDRadar measurements for four footprint sizes FP2f20; 10; 5; 2g corresponding to the diameter of a circular footprint on the ground in meters. The along-prole survey results showing the TWT measurements are shown for both transects in Fig. 7.8. Scatter plot representations of the TWT data, with horizontal and vertical axes corresponding to GPR and UAV-SDRadar values, respectively, are presented in Fig. 7.9 for multiple UAV-SDRadar footprint sizes. In Fig. 7.10 we provide along-prole comparisons of the UAV-SDRadar TWT with the closest GPR TWT within 15 m (or equivalently, within a 30 m diameter footprint) of the UAV-SDRadar platform for the corresponding measurement. Scatter plot representations of the TWT data, with horizontal and vertical axes corresponding to GPR and UAV-SDRadar values, respectively, are presented in Fig. 7.11 for the closest GPR measurement. To demonstrate repeatability of measurements, the Meadow Transect was own twice. A comparison of the rst and second meadow ight paths and measured TWT are given in Fig. 7.12a and 7.12b, respectively. The correlation coecient between the two is ights is 0:83 with an RMSE of 1:33 ns. We note that the actual ight paths diverge due to changing weather conditions and imperfect sUAS ight controls. 139 (a) Flight paths 0 100 200 300 Distance along transect (m) 5 10 15 TWT (ns) Meadow Flight 1 Meadow Flight 2 ;=0.83, RMSE=1.33 ns 0 2 4 6 UAV Separation (m) (b) TWT measurements Figure 7.12: Repeated Meadow Transect ights We report the statistical comparison between the UAV-SDRadar and GPR measurements in Table 7.2, Table 7.3, and Table 7.4 for the rst Meadow, second Meadow, and Forest transects, respectively. For the four footprint sizes as well as the closest measurement comparisons, we report the root mean square error (RMSE), correlation coecient (), and un-based Nash-Sutclie eciency (NSE) as well as the number of GPR measurements considered N GPR in each case. All results are statistically signicant with all P-values associated with the correlation coecients lower than 1:4e28 and 1:9e3 achieved for the Meadow and Forest transects respectively. A snow dielectric of 1:41 is used in these statistics based on measurements taken at a snow pit (See Section 7.4.4). The lower correlation of the Forest Transect data is explained by (i) scattering of the signal by forest canopy, (ii) higher ight elevation of UAV AGL to avoid trees resulting in lower SNR, and (iii) divergence of the UAV path from the ground-based GPR path as show in Fig. 7.5. Table 7.2: Statistical comparison of TWT estimated by UAV-SDRadar with ground-based GPR measurements along the Meadow Transect. y Assuming a snow dielectric of 1:41. Parameter Footprint Size (m) 20 10 5 2 closest RMSE (ns) 0.86 0.93 0.87 0.93 1.00 RMSE (cm) y 10.81 11.76 10.97 11.75 12.60 Corr. 0.89 0.9 0.9 0.85 0.88 Bias (ns) 0.44 0.60 0.51 0.53 0.49 unbiased NSE 0.77 0.81 0.81 0.73 0.77 N GPR 533 496 434 314 548 140 Table 7.3: Statistical comparison of TWT estimated by UAV-SDRadar with ground-based GPR measurements along the second Meadow Transect. y Assuming a snow dielectric of 1:41. Parameter Footprint Size (m) 20 10 5 2 closest RMSE (ns) 0.82 0.88 0.84 0.93 1.16 RMSE (cm) y 10.32 11.06 10.66 11.81 14.64 Corr. 0.82 0.8 0.8 0.79 0.75 Bias (ns) -0.02 0.10 -0.12 -0.11 0.02 unbiased NSE 0.67 0.64 0.64 0.61 0.55 N GPR 589 512 439 122 605 Table 7.4: Statistical comparison of TWT estimated by UAV-SDRadar with ground-based GPR measurements along the Forest Transect. y Assuming a snow dielectric of 1:41. Parameter Footprint Size (m) 20 10 5 2 closest RMSE (ns) 1.40 1.35 1.43 1.38 1.53 RMSE (cm) y 17.74 17.05 18.06 17.42 19.28 Corr. 0.7 0.74 0.68 0.64 0.61 Bias (ns) 0.29 0.55 0.62 0.34 0.20 unbiased NSE -0.59 0.24 0.23 0.26 -0.012 N GPR 527 433 359 192 564 7.4.4 Snow Pit (Michigan River) The UAV-SDRadar collected 202 measurements with a total synthesized bandwidth of 1:5 GHz across the 600 2100 MHz frequency bands while hovering over a snow pit located near the Meadow Transect at the Michigan River Site as marked on the map in Fig. 7.5. Two relative dielectric permittivity proles were measured at the snow pit, which had a depth of 115 cm at the prole measurement points. The average snow depth measured at eight points around the perimeter of the pit was 120:9 cm ( = 4:4 cm). Dielectric measurements were made using a WISe Snow liquid water content sensor from A2 Photonic Sensors. In Fig. 7.13, we compare the average radar backscatter with the measured dielectric proles as well as with the gradient of the dielectric proles. We use the average relative snow pit dielectric r = 1:41 to determine the electromagnetic wave propagation velocity v = c= p r and thus to estimate the snow depth from the radar signal TWT measurements as 125:2 cm ( = 12:6 cm). As shown in Fig. 7.13, 141 0.4 0.5 0.6 Depth (m) 15 20 Power (dB) Scan 147 0.4 0.5 0.6 Depth (m) 15 20 Power (dB) Scan 49 20 40 60 80 100 120 140 160 180 200 Scan 0 1 2 Depth (m) 7.2 cm 9.6 cm (a) Radargram collected as UAV-SDRadar hovered over snow pit. 5 10 15 20 25 Power (dB) -0.5 0 0.5 1 1.5 Depth (m) Air Ground 1 1.5 r radar backscatter r (profile A) r (profile B) 5 10 15 20 25 Power (dB) -0.5 0 0.5 1 1.5 Depth (m) Air Ground 0 0.1 0.2 | r | radar backscatter | r | (profile A) | r | (profile B) (b) Snow Pit Prole showing photograph of the snow pit wall (left), the measured dielectric proles r (center), and the dielectric gradient magnitudejrj (right) compared with the measured radar backscatter. Figure 7.13: Snow pit measurement results. The snow pit has a ground truth snow depth of 120:9 cm and average dielectric permittivity r = 1:41. The radar backscatter is the average of 202 radar scans. Dielectric proles A and B were taken at dierent locations in the snow pit. The average depth estimated by the UAV-SDRadar is 125:2 cm. the UAV-SDRadar accurately measures the ground truth depth of the snow at the snow pit to within 4:4 cm. Furthermore, changes in the dielectric gradient, which appear as stratigraphy in the photograph of the snow pit, are present as peaks in the radar backscatter signal. The theoretical radar resolution is r =c=(2BW p r ). Thus for the 1:5 GHz bandwidth SWW in snow with dielectric 1:41, we expect depth resolution performance of 8:4 cm. To assess the resolution achieved by the UAV-SDRadar for imaging subsurface layering, we examine two metrics: (i) the3 dB peak width of a single sub-surface layer re ection and (ii) the separation between two adjacent sub-surface re ections. In Fig. 7.13a, scan 147 of the radargram features a subsurface 142 re ection at 0:55 m. The3 dB peak width of this re ection is 7:2 cm. Scan 49 features two close but distinct re ections separated by a null that is 2:8 dB down from the stronger of the two peaks. The separation of the two peaks is 9:6 cm. This demonstrates that the UAV-SDRadar is capable of measuring not only snow depth but also that the sensor has high enough resolution and sensitivity to measure sub-surface layering in the snow at the expected theoretical performance. 7.4.5 Preliminary Bistatic Radar Results In April 2021, we conducted a second eld campaign at a site near Winter Park, CO with two UAV- SDRadars operating in a bistatic conguration. The two UAV-SDRadars used the synchronization method described in Chapter 5 to wirelessly synchronize during ight and achieve coherent operation across multiple stepped frequencies from 1 1:4 GHz. (a) On ground (b) In ight Figure 7.14: Bistatic UAV-SDRadar systems at eld site in Winter Park, CO. The two UAV-SDRadar systems are shown on the ground along with two other UAV systems equipped with lidar sensors in Fig. 7.14a. The bistatic UAV-SDRadars are shown in ight over the eld site in Fig. 7.14b. We present preliminary results from the bistatic ight in Fig. 7.15. The two UAV-SDRadars ew autonomously across the eld site transect in a straight line parallel to one another. The ight paths of both are shown in Fig. 7.15a. The synchronized bistatic SWW radargram is shown in Fig. 7.15b. Synchronization was performed at each 40 MHz frequency step across 1 1:4 GHz of total bandwidth. The weaker signal at 3:5 m free space range in Fig. 7.15b is the direct path signal between the two UAV-SDRadars. The relatively stronger signal is the bistatic re ection from the snow surface and additional sub-surface features. There is some indication of distinct re ections 143 (a) GPS Coordinates of ight path 0 100 200 300 Time delay (ns) 20 40 60 80 100 120 140 160 180 200 Scan -10 0 10 20 30 40 50 Free space range (m) (b) Synchronized bistatic radargram Figure 7.15: Preliminary results from bistatic UAV-SDRadar test ight. Two UAV-SDRadar systems ew parallel across transect. Wireless synchronization performed in real-time to achieve coherent operation and to form bistatic SWW formed with 400 MHz total bandwidth. from the snow and ground surfaces, however, additional analysis is required to determine if snow depth can be retrieved from this data. Furthermore, we note that the synchronized bistatic signal shown here does not incorporate RTK/PPK GNSS/GPS corrections for UAV position and therefore we expect that coherence and therefore resolution will improve when positional corrections are applied to post-processed data. This will be the subject of future work. 7.5 Discussion The quality of the lidar altimetry data plays a signicant role in determining the ease with which the UAV elevation AGL (and therefore the snow surface) can be estimated. When the UAV is ying slowly over open terrain unobstructed by foliage, the lidar altimetry data were very reliable and well aligned with the snow surface re ection seen in the radargram. When this is not the case, the snow surface must be estimated to a greater degree directly from features present in radargram and possibly with human assistance. It is expected that as this research area matures, better performing surface and feature detection algorithms, particularly those using machine learning (ML)-based approaches will emerge. Because of the dry snowpack present at the studied sites, the re ection from the air-snow interface is signicantly weaker than that from the snow-ground interface at the frequency range used in our observations. This added diculty to the detection of the snow-surface and required the surface detection algorithm to be tuned with this knowledge a priori. 144 In environments with wetter snow, the radar re ection from the snow surface might be stronger than that from the ground below. This may be due to (i) the increase in dielectric contrast between the air and snow (and potential decrease in dielectric contrast between the snow and ground) and (ii) the increased attenuation of the radar signal through the snow due higher water content [155]. For the results presented, the snow surface was more dicult to identify in the radar signal than the ground surface. However in the case of wet snow, the opposite may be true. In such cases, the proposed surface detection algorithm would have to be modied according to the expected relative strength of the returns. A more general broken-edge surface detection algorithm that can identify not only snow and ground surfaces but also identify an arbitrary unknown number of subsurface layers is desirable. ML-based approaches have shown promise for solving such problem [154] [156] [157] but further investigation of their application in low-altitude UAV-based radar imaging problems is necessary. What is shown in both the UAV-SDRadar and GPR data in Section 7.4 is that spatial variability (along a transect) can be high in some areas, with TWT dierences > 5 ns observed over the span of a few meters, and low in others with near constant TWT across over tens of meters. It stands to reason that the spatial variability of snow depth across areas where the paths taken by the UAV and GPR diverge would exhibit similar characteristics as those observed along the measured transects. Therefore, in addition to the naturally expected cases where small osets yield small measurement dierences and vice versa, we expect to see both cases where the UAV-SDRadar and GPR produce dierent measurements for small spatial osets of only a few meters, and cases where they produce similar measurements for spatial osets exceeding 10 m. Therefore, when interpreting the results, it is important to consider both the measurements obtained by each instrument and the relative spatial information. Radar sensors, in theory, provide an obvious and complementary solution to the limitations of optical sensors (i.e., Lidar and SfM-MVS) for measuring snow depth. The lower frequency radar signals are able to penetrate not only foliage and forest canopy cover, but also snow layers and can produce depth prole imagery of a snowpack from a single pass measurement. Furthermore, because it is possible to estimate material properties from backscattered radar signals using electromagnetic scattering models, it is possible in theory to directly retreive SWE from a single-path radar observation [158] [84] [159] [160] [61]. In such a case, information from multiple frequency bands 145 within the synthesized bandwidth would likely be needed to retrieve the depth-dependent dielectric permittivity prole via electromagnetic inverse scattering techniques, which is the subject of our ongoing work. 7.6 Conclusion In this chapter, we have demonstrated successful use of the novel SWW NUFS algorithm in an ultra- wideband SDRadar sensor on a UAV platform to measure distributed snow depth at a eld site in Colorado. Using a low-cost COTS USRP SDR-based SDRadar sensor, we have demonstrated tunable coherent ultra-wideband radar operation from a moving UAV platform with complete removal of grating lobes and other signal processing artifacts. This is the rst successful demonstration of UAV-based radar imaging and retrieval of snow depth in a heavily forested area and the rst integration of UAV-based radar, SfM-MVS, and lidar altimetry data for snowpack mapping. We have validated the UAV-SDRadar results with ground-based GPR sled measurements, making this the rst eld campaign with results from UAV-based radar compared with GPR-derived ground truth. We have also imaged snow layer stratigraphy from a UAV-SDRadar and compared the results with snow pit dielectric measurements. We have shown successful UAV platform motion compensation using RTK/PPK GNSS/GPS sensor and topography compensation using both independent SfM measurements and an original surface detection algorithm that incorporates lidar altimetry data as well as radargram features to identify snow and ground interfaces. Finally, we have reported preliminary results from a coherent stepped-frequency bistatic UAV- based radar measurement over a snow eld. This is the rst report of this kind. Future work will focus on applying post-processing for UAV motion compensation to improve these results and on additional analysis with comparison to ground truth snow depth data. 146 Chapter 8 UAV-SDRadar Heterogeneous Smart Sensor Network Integration 8.1 Autonomous Remote Control of UAV-SDRadar Missions UAV Data Server (UDS) UAVSDRadar Onboard System M600 Aircraft HDMI In UART (API) Telemetry RF Link Battery/Power Controller GPS IMU Ettus E312 RX-B RX-A TX-A TX-B -3 dB Splitter -20 dB TX RX GPS IMU USB WIFI SDRadar Server Data & Command Server Control & Signal Processing Ethernet Lidar Altimeter USB DJI Remote Controller Manual Control (Optional) DJI Manifold (Jetson TX2) HDMI Out Ethernet UART1 USB WIFI SDRadar Client Config/Command Generation SWW Processing DJI Onboard SDK Mission Controller UAV Charging Station (UCS) USB Ethernet USB/Serial UDS-MCU Energy Management System Ethernet Local Coordinator (LC) USB/Serial MCU Energy Management System LC-RPi Comand/ Config Server Ethernet UDS-RPi Command/ Config Server WiFi WLAN Xbee Radio 3/4G Modem Remote Data Server SoilSCAPE Sensor Network End Devices 900 MHz USB XBee Radio 12 V - Solar Panel Figure 8.1: Hardware diagram of UAV-SDRadar integration into SoilSCAPE/SPECTOR smart heterogeneous sensor network 147 UDS-MCU UDS-RPi UCS UAV-SDRadar Power On UDS-RPi Boot Status Check Status Check Data Server Power On UAV [uav-sdradar_cfg valid] [start UAV mission] Power on UCS Power on UAV Verify/Modify Charging State skysense-cli show-charging-state start-charge stop-charge reset Charging State Get Battery Charge and Flight Status [insufficient battery] Copy uav-sdradar _cfg [battery charged] uav-sdradar _cfg Check for valid uav-sdradar_cfg No Yes Execute Waypoint Mission Flight Log and SDRadar Data [confirm] Download SDRadar Data Software Shutdown Software Shutdown Power Off UAV Power Off (Physical switch disconnect) Power Off (WoR-based battery disconnect) [uav-sdradar_cfg not valid] Legend Power Serial Ethernet WIFI/XBee Internal Battery/Flight Status uav-sdradar _status Configure SDRadar (Ethernet) Configure UAV Check TID and confirm handshake uav-sdradar _confirm Download uav-sdradar_cfg message/ file HW SW uav-sdradar _done go.txt get.txt Upload Logs and Data to Data Server Figure 8.2: Process sequence diagram for integration of UAV-SDRadar with SoilSCAPE/SPECTOR smart heterogeneous sensor network The planned new SPCTOR phase of the SoilSCAPE eort, described in Section 2.8, involves the extension of the optimal sensor control architecture to heterogeneous sensing agents and data products. The culmination of the this eort aims to integrate in situ sensor webs with one or more UAV-based SDRadar sensors as agents in eld-scale operations [161]. The proposed works is to demonstrate the integration of one UAV-based SDRadar in the SoilSCAPE architecture. The following further developments are proposed to achieve this goal: (i) Integration of SDRadar mission control with the optimal sensor controller and scheduler, (ii) Development of ight control and mission planning software to execute autonomous UAV ights, and (iii) standardization of UAV SDRadar calibration routines and data products for scientically valuable characterization of environmental parameters. The goal of objective (i) is to implement procedures for the SoilSCAPE Local Coordinator (LC) to initiate and parameterize a UAV ight and SDRadar operational mode. The LC will `wake' the UAV SDRadar sub-system and transmit a standardized mission objective to the UAV ight controller. The UAV ight controller will then generate a series of GPS waypoints and SDRadar data collection scheme and instruct the SDRadar to initiate radar operation. Furthermore, communication protocols and channels will need to be established for the UAV provide updates about resource availability, such as battery levels, to the LC. 148 In objective (ii), software developed for the ight controller will execute the ight plan, precision landing sequence, and SDRadar collection. After the ight plan is executed, the high resolution radar imagery is processed by the GPU-based ight controller and compressed before relay of data to the LC. Time permitting, the UAV ight controller functionality will be extended to include real-time data processing and quality checking for real-time sensor in the loop ight decision making. SDRadar calibration routines, using xed reference targets, will be designed to accomplish objective (iii). Existing work on soil moisture retrieval by the USC MiXIL lab will be leveraged for estimation of environmental parameters. This new phase represents a multi-year sustained eort, and success in the proposed eort will be measured by proof of concept demonstrations of protocols for achieving objectives i-iii. The system diagram of the UAV-SDRadar system integrated with the SoilSCAPE/SPCTOR sensor network control system is shown in Fig. 8.1. The process sequence diagram for initiating an autonomous UAV ight and SDRadar data collection from a remote data server is shown in Fig. 8.2. Figure 8.3: Integrated UAV-SDRadar system with onboard lidar altimeter, Nvidia Jetson TX-2 based DJI Manifold SBC, and remote power control module 149 In August 2021, we performed the rst eld demonstration of the integrated UAV-SDRadar system conducting autonomous ights at an existing sensor network deployment site in Walnut Gulch, Arizona. The autonomous UAV-SDRadar system is shown at the eld test site in Walnut Gulch, Arizona in Fig. 8.3. (a) Flight 1 GPS (b) Flight 1 radargram (c) Flight 2 GPS (d) Flight 2 radargram Figure 8.4: Flight paths and radargram data collected on autonomous ights at SPCTOR eld test site in Walnut Gulch, AZ. Results from the rst eld deployment are shown in Fig. 8.4. The UAV-SDRadar performed multiple autonomous ights and collected radar data using 650 1650 MHz SWW bandwidth as dened by mission plan retrieved from a remote data server. Multiple ights were conducted and repeated across multiple days before and after rain at test site. Flight 1 was conducted multiple times along the same ight path as shown in Fig. 8.4a. A comparison of the radar data collected before and after rainfall is given in Fig. 8.5. This 150 (a) Before rain (b) After rain Figure 8.5: Flight 1 repeated over multiple days before and after rainfall at test site in Walnut Gulch, AZ. The resulting radar data shows sensitivity to the eect of rainfall. demonstrates the the SDRadar sensor is sensitive to the soil moisture change introduced by light rainfall. Complete retrieval of soil moisture proles from the UAV-SDRadar data is an area of future research and will be explored in subsequent works. 151 Appendix A Derivation of Frequency Stacking Algorithm This derivation has been published previously in [4]. Consider a discrete-time baseband complex LFM chirp waveform w[t m ] with instantaneous bandwidth B i and pulse length T p , which is sampled according to the Nyquist theorem at sampling rate f s and has support over the interval t m 2f Tp 2 ; Tp 2 g w[t m ] =A[t m ]e jKt 2 m (A.1) whereK = B i Tp is the chirp rate andA[t m ] is a pulse amplitude weighting function. For a nominal rectangular pulse, we may take A[t m ] = rect[ tm Tp ]. Using the Whittaker-Shannon interpolation formula, the continuous time waveformw(t) produced after digital to analog conversion of w[t m ] can be represented as w(t) = 1 X m=1 w[t m ] sin f s (tt m ) f s (tt m ) (A.2) We dene a series ofN linear FM chirp sub-pulsess n (t) each having bandwidthB i , pulse length T p and center carrier frequency f n =f c0 +nf c for n = 0;:::;N 1 where f c is the sub-pulse frequency spacing. We assume that each sub-pulse is transmitted and received in a stop-and-go 152 fashion suciently separated in time so that sub-pulses can be treated as orthogonal and that target echoes are unambiguous. We note that sub-pulse doppler shifts due to relative target motion are beyond the scope of this work. Due to the generally non-coherent TX and RX LO phase relationship in many commercial radio platforms, including USRP SDRs, we consider terms for the TX phase t;n and RX phase r;n errors present in each sub-pulse. We dene the total contribution of these random phase errors after up/down-conversion as e;n = t;n r;n . Then s n (t) =w(t)e j 2fnt e jt;n (A.3) Letz n (t) be the echo of then th sub-pulse from a scatter located at a distanceR s after demodulation and down-conversion to baseband z n (t) =s n (t 2R s c )e j 2fnt e jr;n (A.4) substituting A.3 and after the signal is digitized with ADCs, we get the discretized signal z n [t m ] z n [t m ] =w[t m 2R s c ]e j 2fn 2Rs c e je;n (A.5) By performing a frequency shift, time shift, phase correction and sinc interpolation with a B s bandwidth lter on the sub-pulse echoes, a synthetic wideband signal z[t m ] that has the characteristics of a wideband LFM chirp may be reconstructed. These operations can be expressed as a set of lters with coecients g n [t m ] given as [27], [28], [29] g n [t m ] = 1 2f c sin[B s (t m T n )] (t m T n ) e jTnfn e j 2(tmTn)fn (A.6) The lters are applied to each frequency shifted baseband sub-pulse and the outputs are summed to obtain the SWW z[t m ] z[t m ] = N1 X n=0 z n [t m ]e j 2fntm e je;n ~g n [t m ] (A.7) 153 where~ denotes the convolution operator. If we make the assumption of a stationary platform throughout the duration of each frequency sweep, the SWW reconstruction is equivalent to frequency-domain stacking of compressed sub-pulse spectra. If we consider pulse compression of the received echo SWW z[t m ] by matched ltering with a SWW v[t m ] similarly constructed from a series of N phase-coherent reference sub-pulses v n [t m ], i.e., v n [t m ] =w[t m ]e je;n (A.8) which is equivalent to the expression in (A.5) with R s = 0. The output of the wideband matched lter d[t m ] is the pulse-compressed SWW and is given by the cross-correlation d[t m ] =z[t m ]~v [t m ] (A.9) We dene D[f k ], Z n [f k ], and V n [f k ] as the discrete Fourier transforms of d[t m ], z n [t m ], and v n [t m ] respectively and G n [f k ] as the frequency response of the g n [t m ] lters given in (A.6). Now d[t m ] may be represented in the frequency domain as D[f k ] = X n Z n [f k f n ]G n [f k ] X n 0 V n 0[f k f n 0]G n 0[f k ] (A.10) Noting that G n [f k ] only has support over f2 [f n B s =2; f n +B s =2] and is zero elsewhere, we can simplify this expression as D[f k ] = X n Z n [f k f n ]V n [f k f n ]G n [f k ]G n [f k ] (A.11) Time-independent complex exponential terms in (A.6) and hence in G n [f k ] will cancel due to multiplication by conjugate phase. This then simplies to D[f k ] = X n Z n [f k f n ]V n [f k f n ] rect h f k f n B s i (A.12) 154 Note, the exponents of the coecients are removed. Because z n [t m ] and v n [t m ] are both baseband complex low-pass signals, if we represent the set of lowpass sub-pulse matched lter outputs d n [t m ] =z n [t m ]~v n [t m ] with their Fourier transform D n [f k ] =Z n [f k ]V n [f k ], then D[f k ] = X n D n [f k f n ] rect h f k f n B s i (A.13) or in the time domain d[t m ] = X n d n [t m ]~ sin[B s t m ] t m e j 2fntm (A.14) 155 Appendix B Derivation of Arbitrary NLFM Waveforms Using the design strategy for NLFM waveforms described in Section 2.2.2, a NLFM waveform may be derived from an LFM waveform with arbitrary amplitude weighting as follows [19,22]. Consider a bandwidth B low-pass signal x(t) x(t) =a(t) exp(j(t)) (B.1) with amplitude a(t) and phase (t). The instantaneous frequency f i at time t i =t is f i = 1 2 @ @t (t i ) (B.2) where B 2 f i B 2 . It has been shown that the PSP may be used to relate the PSD to the amplitude and phase by considering X(f) =A(f) exp(j(f)), the Fourier transform of x(t) [19]. Now: x(t) = Z 1 1 A(f) exp(j((f) + 2ft))df (B.3) X(f) = Z 1 1 a(t) exp(j((t) 2ft))dt (B.4) 156 The PSP tells us that the integral of these oscillating functions is small except for regions where the phase is stationary, that is where: @ @f ((f) + 2ft) = 0 and @ @t ((t) 2ft) = 0 (B.5) Using a Taylor series expansions at the stationary points, the following relationship can be shown for the time domain [19]: jX(f i )j 2 2 a 2 (t i ) j 00 (t i )j (B.6) and in the frequency domain: jx(t)j 2 2 A 2 (f i ) j 00 (f i )j (B.7) We consider some signal z(t) with PSD Z 2 (f)X 2 (f) and constant amplitude A, so that jZ(f i )j 2 2 A j 00 (t i )j and 00 (f) 2 Z 2 (f) A (B.8) whereZ(f) is dened over the signal bandwidthB=2fB=2. The rst derivative of the phase 0 (f) is related to the group time delay function T g (f) as: T g (f) = 1 2 0 (f) (B.9) =c 1 Z f B=2 Z 2 ()d +c 2 (B.10) where c 1 and c 2 are constants determined by the boundary condition requirements that T g (B=2) = T=2 and T g (B=2) =T=2 [24], [25]. The group time delay function is inverted to obtain the instantaneous frequency f(t) of the NLFM f(t) =T 1 g (f) (B.11) 157 and the NLFM signal phase is then (t) = 2 Z t T=2 f()d (B.12) so that z(t) =A exp(j(t)). 158 Appendix C Derivation of TOF Cramer Rao Lower Bound The derivation in this section has been previously published in [32]. The purpose of this appendix is to derive the CRLB. We begin with the general CRLB for a complex autocorrelation signal s(t) in complex white Gaussian noise with variance 2 N . The Fischer information matrix I(t pk ) for a delay parameter t pk is [162] I(t pk ) = 2 2 N Re N1 X n=0 @s[n;t pk ] @t pk 2 (C.1) = 2 2 N Re N1 X n=0 @s(t) @t t=n 2 (C.2) Then the CRLB is 2 TOF 2 N =2 N1 P n=0 @s(t) @t t=n 2 (C.3) where = 1 fs . Approximating the sum with an integral 2 TOF 2 N =2 1 R Tp 0 @s(t) @t 2 dt (C.4) 159 using Fourier theory and Parsevals equation, Z Tp 0 @s(t) @t 2 dt = Z 1 1 (2f) 2 jS(f)j 2 df (C.5) and the signal power P s is P s = 1 T p Z Tp 0 js(t)j 2 dt (C.6) = 1 T p Z 1 1 jS(f)j 2 df (C.7) For a complex linear FM chirp waveform having bandwidth B, we make the approximation jS(f)jjS(0)j rect(f=B) (C.8) and thus Z 1 1 (2f) 2 jS(f)j 2 df (C.9) Z B=2 B=2 (2f) 2 jS(0)j 2 df (C.10) =jS(0)j 2 2 B 3 3 (C.11) similarly Z 1 1 jS(f)j 2 df (C.12) Z B=2 B=2 jS(0)j 2 df (C.13) =jS(0)j 2 B (C.14) therefore 160 2 TOF 3 2 N 2f s jS(0)j 2 2 B 3 (C.15) Noting that jS(0)j 2 = P s T p B (C.16) and that SNR =P s = 2 N we arrive at the approximate CRLB for time delay estimation 2 TOF 3 2 2 B 2 f s T p SNR (C.17) We crucially note that 2 TOF is the CRLB for one way TOF (i.e., TOA or peak) estimation. Because the synchronization method described in this work calculates the synchronized TOF as the mean of two sensors' local TOF measurements (5.22), which represent two independent observations, the CRLB for the two-way TOF is 2 2TOF = 2 TOF =2. 161 Appendix D Derivation of Synchronized Clock Error PSD with Moving Least Squares State Prediction Taking the autocorrelation of t as dened in Section 6.3.2. R ; = t t (D.1) = 1 C 2 N 2 h X i X j ijR i;j + X i i X j j X i X j R i;j N h X i i( X i X j iR i;j + X i X j jR i;j ) (D.2) = 1 C 2 N 2 h X i X j ijR i;j + N 2 h (N h 1) 2 4 X i X j R i;j N 2 h (N h 1) 2 X i X j (i +j)R i;j (D.3) where R i;j (t) = i (t) j (t) (D.4) =R(t) +R n;n (t) +R(t + (ij)SRI)R(t +iSRI)R(tjSRI) +(ij)R n;n (t)(i)R n;n (t)(j)R n;n (t) (D.5) 162 We also dene R ; i =(t) i (t) (D.6) =R(t + (N h 1i)SRI + T p )R(t + (N h 1)SRI + T p ) (D.7) R i ; = i (t) (t) (D.8) =R(t (N h 1i)SRI T p )R(t (N h 1)SRI T p ) (D.9) R ; i = (t T p ) i (t) (D.10) =R(t + (N h 1i)SRI)R(t + (N h 1)SRI) +(i (N h 1))R n;n (t) (D.11) R i ; = i (t) (t + T p ) (D.12) =R(t (N h 1i)SRI)R(t (N h 1)SRI) +(i (N h 1))R n;n (t) (D.13) R ; =(t) t (D.14) = 1 C N h X i iR ; i (t) X i i X i R ; i (t) (D.15) R ; = t (t) (D.16) = 1 C N h X i iR i ; (t) X i i X i R i ; (t) (D.17) R ; = (t T p ) t (D.18) = 1 C N h X i iR ; i (t) X i i X i R ; i (t) (D.19) R ; = t (t + T p ) (D.20) = 1 C N h X i iR i ; (t) X i i X i R i ; (t) (D.21) (D.22) Now the clock phase error is ~ (t) =(t) ( (t T p ) +T p ) (D.23) 163 So that ~ R(t) = ~ (t) ~ (t) (D.24) = 2R(t) +R n;n (t)R(t T p )R(t + T p ) (D.25) T p (R ; +R ; R ; R ; ) + T 2 p R ; = 2R(t) +R n;n (t)R(t T p )R(t + T p ) (D.26) T p C N h X i iR ; i (t) X i i X i R ; i (t) +N h X i iR i ; (t) X i i X i R i ; (t) + T p C N h X i iR ; i (t) X i i X i R ; i (t) +N h X i iR i ; (t) X i i X i R i ; (t) + T p C 2 N 2 h X i X j ijR i;j + X i i X j j X i X j R i;j N h X i i( X i X j iR i;j + X i X j jR i;j ) Expanding terms X i iR ; i (t) = X i i(R(t + (N h 1i)SRI + T p )R(t + (N h 1)SRI + T p )) (D.27) = X i iR(t + (N h 1i)SRI + T p ) (D.28) N h (N h 1) 2 R(t + (N h 1)SRI + T p ) X i iR i ; (t) = X i iR(t (N h 1i)SRI T p ) (D.29) N h (N h 1) 2 R(t (N h 1)SRI T p ) X i i X i R ; i (t) = X i i X i R(t + (N h 1i)SRI + T p )R(t + (N h 1)SRI + T p ) (D.30) = N h (N h 1) 2 X i R(t + (N h 1i)SRI + T p ) (D.31) N h R(t + (N h 1)SRI + T p ) X i i X i R i ; (t) = N h (N h 1) 2 X i R(t (N h 1i)SRI T p ) (D.32) N h R(t (N h 1)SRI T p ) 164 X i iR ; i (t) = X i i(R(t + (N h 1i)SRI)R(t + (N h 1)SRI) +(i (N h 1))R n;n (t)) (D.33) = X i iR(t + (N h 1i)SRI) (D.34) N h (N h 1) 2 R(t + (N h 1)SRI) + (N h 1)R n;n (t) X i iR i ; (t) = X i iR(t (N h 1i)SRI) (D.35) N h (N h 1) 2 R(t (N h 1)SRI) + (N h 1)R n;n (t) X i i X i R ; i (t) = X i i X i (R(t + (N h 1i)SRI)R(t + (N h 1)SRI) (D.36) +(i (N h 1))R n;n (t)) = N h (N h 1) 2 X i R(t + (N h 1i)SRI) (D.37) N h R(t + (N h 1)SRI) +R n;n (t) X i i X i R i ; (t) = N h (N h 1) 2 X i R(t (N h 1i)SRI) (D.38) N h R(t (N h 1)SRI) +R n;n (t) X i X j ijR i;j = X i X j ij(R(t) +R n;n (t) +R(t + (ij)SRI)R(t +iSRI) (D.39) R(tjSRI) +(ij)R n;n (t) ((i) +(j))R n;n (t)) = N 2 h (N h 1) 2 4 R(t) +R n;n (t) (D.40) N h (N h 1) 2 X i iR(t +iSRI) + X i iR(tiSRI) + N h (N h 1)(2N h 1) 6 R n;n (t) + X i X j ijR(t + (ij)SRI) 165 X i X j R i;j = X i X j (R(t) +R n;n (t) +R(t + (ij)SRI)R(t +iSRI) (D.41) R(tjSRI) +(ij)R n;n (t) ((i) +(j))R n;n (t)) =N 2 h R(t) +N h (N h 1)R n;n (t) (D.42) N h X i R(t +iSRI) + X i R(tiSRI) + X i X j R(t + (ij)SRI) X i X j iR i;j = X i X j i(R(t) +R n;n (t) +R(t + (ij)SRI)R(t +iSRI) (D.43) R(tjSRI) +(ij)R n;n (t) ((i) +(j))R n;n (t)) = N h (N h 1) 2 N h R(t) +N h R n;n (t) (D.44) N h X i iR(t +iSRI) N h (N h 1) 2 X i R(tiSRI) + X i X j iR(t + (ij)SRI) X i X j jR i;j = X i X j j(R(t) +R n;n (t) +R(t + (ij)SRI)R(t +iSRI) (D.45) R(tjSRI) +(ij)R n;n (t) ((i) +(j))R n;n (t)) = N h (N h 1) 2 N h R(t) +N h R n;n (t) (D.46) N h X j jR(tjSRI) N h (N h 1) 2 X j R(t +jSRI) + X i X j jR(t + (ij)SRI) Now 166 ~ R(t) = 2R(t) +R n;n (t)R(t T p )R(t + T p ) (D.47) T p C N h X i i(R ; i (t) +R i ; (t)) X i i X i (R ; i (t) +R i ; (t)) + T p C N h X i i(R ; i (t) +R i ; (t)) X i i X i (R ; i (t) +R i ; (t)) + T p C 2 N 2 h X i X j ijR i;j + N 2 h (N h 1) 2 4 X i X j R i;j N 2 h (N h 1) 2 X i X j (i +j)R i;j (D.48) ~ R(t) = 2R(t) +R n;n (t)R(t T p )R(t + T p ) (D.49) T p C N h X i i R(t + (N h 1i)SRI + T p ) +R(t (N h 1i)SRI T p ) N h (N h 1) 2 X i R(t + (N h 1i)SRI + T p ) +R(t (N h 1i)SRI T p ) + T p C N h X i i R(t + (N h 1i)SRI) +R(t (N h 1i)SRI) N h (N h 1) 2 X i R(t + (N h 1i)SRI) +R(t (N h 1i)SRI) +N h (N h 1)R n;n (t) + T p C 2 N 2 h X i X j ijR(t + (ij)SRI) + N 2 h (N h 1) 2 4 X i X j R(t + (ij)SRI) N 2 h (N h 1) 2 X i X j i(R(t + (ij)SRI) +R(t (ij)SRI)) + N 3 h (N 2 h 1) 12 R n;n (t) 167 Using the properties of the Fourier transform, we arrive at the following ~ S (f) =S (f) h (2 2 cos(2fT p )) 2T p C N h X i (i N h 1 2 ) cos(2f((N h 1i)SRI + T p )) cos(2f(N h 1i)SRI) + T p C 2 N 2 h X i X j (ij + (N h 1) 2 4 ) exp(j 2f(ij)SRI) N 2 h (N h 1) X i X j i cos(2f(ij)SRI) i +S n (f) h 1 T p C N h (N h 1) + T p C 2 N 3 h (N 2 h 1) 12 i (D.50) 168 Appendix E Snow and Ground Surface Detection Algorithm This algorithm will be published in [60]. The surface detection algorithm developed for snow depth retrieval from UAV-SDRadar data is summarized as: 1. Apply moving median with a std deviation-based 99:87% condence outlier rejection lter to lidar data 2. Apply Kalmann lter to pre-ltered lidar data to obtain UAV height above ground level (AGL) estimates L i for each scan i. 3. For each scan Z i , compute the mean power Z i and estimate the leading edge of the radar return S 0 i S 0 i = min r j r j jfZ i;j > zc Z i ; L i zc <r j <L i + + zc g (E.1) where zc , zc , and + zc are empirically derived constants representing the mean scale factor, and the upper and lower bounds on the expected deviation of the lidar altimetry data w/r/t the surface return, respectively. 169 4. Compute the 11-point centered moving average across scans of S 0 i to obtain S 0 i S 0 i = 1 2M + 1 l=i+M X l=iM S 0 l (E.2) where M = 5. 5. Estimate the ground surface G 0 i as the highest peak in the received SWW correlation signal G 0 i G 0 i = argmax r j Z i;j jfr j > S 0 i g (E.3) 6. Estimate the average snow depth d 0 = 1 N P i jG 0 i S 0 i j 7. Apply a CA-CFAR threshold [163] to each SWW scan to obtain A i =CFAR CA (Z i ) 8. Compute and store column-wise image gradient of the gaussian ltered radargram. r y A = 2 6 6 6 6 4 1 2 1 0 0 0 1 2 1 3 7 7 7 7 5 (H 5 A) (E.4) whereH 5 is a 5 5 gaussian kernel with = 1. 9. Initialize the ground estimate ~ G 0 i = S 0 i + d 0 and surface estimates as ~ S 0 i = S 0 i . 10. Dene P i;j = 8 > > > > > > > < > > > > > > > : 1 11 ( l=j+5 P l=j5 A i;l ) + max j10<l<j r y A i;l if A i;j >A i;j1 and A i;j >A i;j+1 1 otherwise (E.5) 11. For iteration k := 1 to 2 170 Initialize the cost functions of the rst scan i = 0 ~ S k 0 = argmin r j 1 jr j ~ S k1 0 j 3 P 0;j (E.6) ~ G k 0 = argmin r j 1 jr j ~ G k1 0 j 3 P 0;j (E.7) 12. For each scan i> 1 13. Find all peaks within 30 dB of the maximum and compute the cost functions 14. if k==1 C S i;j = 2 jr j ~ S k i1 j 3 P i;j (E.8) C G i;j = 2 jr j ~ G k i1 j 3 P i;j (E.9) 15. else 16. Estimate three linear least squares ts: LS C using 160 points centered around but not including i,LS L using 40 points to the left of i, andLS R using 40 points to right of i. For each, compute the RMSE of the t. C S i;j = jLS L ~ S k1 i j 2 RMSE L + jLS R ~ S k1 i j 2 RMSE R + jLS C ~ S k1 i j 2 RMSE C 3 P i;j (E.10) C G i;j = jLS L ~ G k1 i j 2 RMSE L + jLS R ~ G k1 i j 2 RMSE R + jLS C ~ G k1 i j 2 RMSE C 3 P i;j (E.11) 17. end if ~ S k i = argmin r j 1 jr j ~ S k1 i j k +C S i;j (E.12) ~ G k i = argmin r j 1 jr j ~ G k1 i j k +C G i;j (E.13) 171 where 1 , 2 , and 3 are empirically derived constant weighting factors. 18. If a manual anchor value is specied for either the ground of surface, use that instead for ~ S i or ~ G i as appropriate. 19. end for 20. Perform cubic interpolation of estimates across scans using the cost function values that are local minima and/or below the median as anchor points. 21. Update snow and ground surface estimates. 22. end for 23. Output surface/ground detections are local peaks nearest to interpolated estimates. 172 Appendix F UAV-SDRadar Snow Measurement Radargrams The data from this section will be published in [60]. In Fig. F.1, we provide radargrams prior to RTK GNSS/GPS UAV elevation correction to show the UAV height AGL across the transects. In Fig. F.2, we provide RTK GNSS/GPS UAV elevation corrected radargrams without colorization of the detected snow/ground surfaces so that the reader may compare the snow/ground surface pixels selected by the proposed algorithm with the raw radargram. 173 0 50 100 Time delay (ns) 50 100 150 200 250 300 350 Distance along transect (m) 0 5 10 15 20 Free space range (m) (a) Meadow Transect 200 250 300 350 400 Time delay (ns) 50 100 150 200 250 300 350 Distance along transect (m) 30 35 40 45 50 55 60 Free space range (m) (b) Forest Transect Figure F.1: Meadow and Forest Transect radargrams prior to UAV altitude correction showing the elevation of the UAV relative to the snow surface. See also Fig. 7.6 and 7.7. 174 100 120 140 160 180 Time delay (ns) 50 100 150 200 250 300 350 Distance along transect (m) 13 18 23 28 Free space range (m) (a) Meadow Transect 240 260 280 300 320 Time delay (ns) 50 100 150 200 250 300 350 Distance along transect (m) 35 40 45 50 Free space range (m) (b) Forest Transect Figure F.2: Meadow and Forest Transect radargrams after RTK GNSS/GPS UAV elevation correction showing the snow surface topography and relative to the UAV elevation AGL. See also Fig. 7.6 and 7.7. 175 References [1] NASA, \2014 Science Plan," National Aeronautics and Space Administration, Tech. Rep., 2014. [2] M. C. Wicks, \Radar the next generation-sensors as robots," in Radar Conference, 2003. Proceedings of the International. IEEE, 2003, pp. 8{14. [3] M. Richards, J. Scheer, and W. Holm, Principles of Modern Radar: Basic Principles., 2010. [4] S. Prager, T. Thrivikraman, M. Haynes, J. Stang, D. Hawkins, and M. Moghaddam, \Ultra- Wideband Synthesis for High-Range-Resolution Software Dened Radar," IEEE Transactions on Instrumentation and Measurement, 2019. [5] W. Wiesbeck, \SDRS: Software-dened radar sensors," in Geoscience and Remote Sensing Symposium, 2001. IGARSS'01. IEEE 2001 International, vol. 7. IEEE, 2001, pp. 3259{3261. [6] J. Ralston and C. Hargrave, \Software dened radar: An open source platform for prototype GPR development," in 2012 14th International Conference on Ground Penetrating Radar (GPR), Jun. 2012, pp. 172{177. [7] T. Debatty, \Software dened RADAR a state of the art," in Cognitive Information Processing (CIP), 2010 2nd International Workshop On. IEEE, 2010, pp. 253{257. [8] J. Meier, R. Kelley, B. M. Isom, M. Yeary, and R. D. Palmer, \Leveraging Software-Dened Radio Techniques in Multichannel Digital Weather Radar Receiver Design," IEEE Transactions on Instrumentation and Measurement, vol. 61, no. 6, pp. 1571{1582, Jun. 2012. [9] B. L. Cheong, R. Kelley, R. D. Palmer, Y. Zhang, M. Yeary, and T.-Y. Yu, \PX-1000: A Solid-State Polarimetric X-Band Weather Radar and Time{Frequency Multiplexed Waveform for Blind Range Mitigation," IEEE Transactions on Instrumentation and Measurement, vol. 62, no. 11, pp. 3064{3072, Nov. 2013. [10] J. Mitola, \Software radios-survey, critical evaluation and future directions," in [Proceedings] NTC-92: National Telesystems Conference, May 1992, pp. 13/15{13/23. [11] F. Harris and W. Lowdermilk, \Software dened radio: Part 22 in a series of tutorials on instrumentation and measurement," IEEE Instrumentation Measurement Magazine, vol. 13, no. 1, pp. 23{32, Feb. 2010. [12] L. K. Patton, \A GNU radio based software-dened radar," 2007. [13] A. E. Prasetiadi and A. B. Suksmono, \A simple delay compensation system in software- dened Frequency Modulated Continuous (FMCW) Radar," in Wireless Conference (European Wireless), 2012 18th European. VDE, 2012, pp. 1{4. [14] S. Costanzo, F. Spadafora, A. Borgia, H. O. Moreno, A. Costanzo, and G. Di Massa, \High Resolution Software Dened Radar System for Target Detection," Journal of Electrical and Computer Engineering, vol. 2013, pp. 1{7, 2013. 176 [15] A. B. Suksmono, \A simple solution to the uncertain delay problem in usrp based sdr-radar systems," arXiv preprint arXiv:1309.4843, 2013. [16] Z. Zhao, M. Yao, X. Deng, K. Yuan, H. Li, and Z. Wang, \A Novel Ionospheric Sounding Radar Based on USRP," IEEE Geoscience and Remote Sensing Letters, vol. 14, no. 10, pp. 1800{1804, Oct. 2017. [17] Y.-K. Kwag, I.-S. Woo, H.-Y. Kwak, and Y.-H. Jung, \Multi-mode SDR radar platform for small air-vehicle Drone detection," in 2016 CIE International Conference on Radar (RADAR). Guangzhou, China: IEEE, Oct. 2016, pp. 1{4. [18] S. Prager and M. Moghaddam, \Application of Ultra-Wideband Synthesis in Software Dened Radar for UAV-based Landmine Detection," in IGARSS 2019 - 2019 IEEE International Geoscience and Remote Sensing Symposium, Yokohama, Japan, pp. 10 115{10 118. [19] N. Levanon and E. Mozeson, Radar Signals. Norwood, MA: John Wiley & Sons, Inc, 2004. [20] D. Wehner, High-Resolution Radar, 2nd ed. Norwood, MA: Artech House, 1995. [21] J. R. Klauder, A. C. Price, S. Darlington, and W. J. Albersheim, \The Theory and Design of Chirp Radars," Bell System Technical Journal, vol. 39, no. 4, pp. 745{808, Jul. 1960. [22] S. Prager, D. Hawkins, and M. Moghaddam, \Arbitrary Nonlinear FM Waveform Construction and Ultra-Wideband Synthesis," in IGARSS 2020 - 2020 IEEE International Geoscience and Remote Sensing Symposium, Waikoloa, HI, 2020. [23] E. Fowle, \The design of FM pulse compression signals," IEEE Transactions on Information Theory, vol. 10, no. 1, pp. 61{67, 1964. [24] R. Ghavamirad, H. Babashah, and M. A. Sebt, \Nonlinear FM waveform design to reduction of sidelobe level in autocorrelation function," in 2017 Iranian Conference on Electrical Engineering (ICEE). Tehran, Iran: IEEE, May 2017, pp. 1973{1977. [25] S. Boukea, Y. Jiang, and T. Jiang, \Sidelobe reduction with nonlinear frequency modulated waveforms," in Signal Processing and Its Applications (CSPA), 2011 IEEE 7th International Colloquium On. IEEE, 2011, pp. 399{403. [26] D. Rabideau, \Nonlinear synthetic wideband waveforms," in Proceedings of the 2002 IEEE Radar Conference (IEEE Cat. No.02CH37322). Long Beach, CA, USA: IEEE, 2002, pp. 212{219. [27] Y. Wang, \BANDWIDTH SYNTHESIS FOR STEPPED CHIRP SIGNAL: A MULTICHAN- NEL SAMPLING PROSPECTIVE," in Radar Conference 2013, IET International, Jan. 2013. [28] W. Chongyu, Y. Yang, and C. Junxian, \Chirp sub-pulse stepped frequency radar signal processing," in Information Processing (ISIP), 2010 Third International Symposium On. IEEE, 2010, pp. 251{254. [29] R. T. Lord and M. R. Inggs, \High range resolution radar using narrowband linear chirps oset in frequency," in Proceedings of the 1997 South African Symposium on Communications and Signal Processing. COMSIG '97, Sep. 1997, pp. 9{12. 177 [30] B. M. Keel, J. A. Saold, M. R. Walbridge, and J. Chadwick, \Nonlinear stepped chirp waveforms with subpulse processing for range side lobe suppression," in Aerospace/Defense Sensing and Controls, R. Trebits and J. L. Kurtz, Eds., Orlando, FL, Aug. 1998, pp. 87{98. [31] D. E. Maron, \Frequency-jumped burst waveforms with stretch processing," in IEEE Interna- tional Conference on Radar, May 1990, pp. 274{279. [32] S. Prager, M. S. Haynes, and M. Moghaddam, \Wireless Subnanosecond RF Synchronization for Distributed Ultrawideband Software-Dened Radar Networks," IEEE Transactions on Microwave Theory and Techniques, vol. 68, no. 11, pp. 4787{4804, Nov. 2020. [33] W. Wang, X. Liang, and C. Ding, \Time and phase synchronisation via direct-path signal for bistatic synthetic aperture radar systems," IET Radar, Sonar & Navigation, vol. 2, no. 1, pp. 1{11, Feb. 2008. [34] A. Einstein, \Zur Elektrodynamik bewegter K orper," Annalen der Physik, vol. 322, no. 10, pp. 891{921, 1905. [35] H. Poincar e, \La Mesure du Temps," Revue de M etaphysique et de Morale, vol. 6, pp. 1{13, 1898. [36] D. Kirchner, \Two-way time transfer via communication satellites," Proceedings of the IEEE, vol. 79, no. 7, pp. 983{990, Jul. 1991. [37] D. L. Mills, Computer Network Time Synchronization: The Network Time Protocol. CRC Press, Mar. 2006. [38] \IEEE standard for a precision clock synchronization protocol for networked measurement and control systems," IEEE Std 1588-2008 (Revision of IEEE Std 1588-2002), pp. 1{300, Jul. 2008. [39] J. C. Eidson, Measurement, Control, and Communication Using IEEE 1588, ser. Advances in Industrial Control. London, UK: Springer, 2006. [40] T. W lostowski, \Precise time and frequency transfer in a White Rabbit network," Master of Science Thesis, Institute of Radioelectronics, Warsaw University of Technology, 2011. [41] J. Elson, L. Girod, and D. Estrin, \Fine-Grained Network Time Synchronization using Reference Broadcasts," in Proceedings of the 5th Symposium on Operating Systems Design and Implementation, Boston, MA, 2002. [42] S. Ganeriwal, R. Kumar, and M. B. Srivastava, \Timing-sync Protocol for Sensor Networks," SenSys'03: Proceedings of the First International Conference on Embedded Networked Sensor Systems, 2004. [43] M. Marti, B. Kusy, G. Simon, and . Ldeczi, \The ooding time synchronization protocol," in Proceedings of the 2nd International Conference on Embedded Networked Sensor Systems - SenSys '04. Baltimore, MD, USA: ACM Press, 2004, p. 39. [44] K.-L. Noh, Q. M. Chaudhari, E. Serpedin, and B. W. Suter, \Novel Clock Phase Oset and Skew Estimation Using Two-Way Timing Message Exchanges for Wireless Sensor Networks," IEEE Transactions on Communications, vol. 55, no. 4, pp. 766{777, Apr. 2007. 178 [45] S. Roehr, P. Gulden, and M. Vossiek, \Method for High Precision Clock Synchronization in Wireless Systems with Application to Radio Navigation," in 2007 IEEE Radio and Wireless Symposium. Long Beach, CA, USA: IEEE, 2007, pp. 551{554. [46] D. Zachariah, S. Dwivedi, P. Handel, and P. Stoica, \Scalable and Passive Wireless Network Clock Synchronization in LOS Environments," IEEE Transactions on Wireless Communica- tions, vol. 16, no. 6, pp. 3536{3546, Jun. 2017. [47] M. Segura, S. Niranjayan, H. Hashemi, and A. F. Molisch, \Experimental demonstration of nanosecond-accuracy wireless network synchronization," in 2015 IEEE International Confer- ence on Communications (ICC). London: IEEE, Jun. 2015, pp. 6205{6210. [48] B. Denis, J.-B. Pierrot, and C. Abou-Rjeily, \Joint distributed synchronization and position- ing in UWB ad hoc networks using TOA," IEEE Transactions on Microwave Theory and Techniques, vol. 54, no. 4, pp. 1896{1911, Jun. 2006. [49] A. Dongare, P. Lazik, N. Rajagopal, and A. Rowe, \Pulsar: A Wireless Propagation-Aware Clock Synchronization Platform," in 2017 IEEE Real-Time and Embedded Technology and Applications Symposium (RTAS). Pittsburg, PA, USA: IEEE, Apr. 2017, pp. 283{292. [50] K. Yu, Y. Guo, and M. Hedley, \TOA-based distributed localisation with unknown internal delays and clock frequency osets in wireless sensor networks," IET Signal Processing, vol. 3, no. 2, p. 106, 2009. [51] R. Exel, \Receiver Design for Time-Based Ranging with IEEE 802.11b Signals," International Journal of Navigation and Observation, vol. 2012, pp. 1{15, 2012. [52] H. Yan, S. Hanna, K. Balke, R. Gupta, and D. Cabric, \Software Dened Radio Implementation of Carrier and Timing Synchronization for Distributed Arrays," arXiv:1811.02098 [eess], Nov. 2018. [53] M. W. Overdick, J. E. Caneld, A. G. Klein, and D. R. Brown, \A software-dened radio implementation of timestamp-free network synchronization," in 2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). New Orleans, LA: IEEE, Mar. 2017, pp. 1193{1197. [54] F. Quitin, M. M. U. Rahman, R. Mudumbai, and U. Madhow, \A Scalable Architecture for Distributed Transmit Beamforming with Commodity Radios: Design and Proof of Concept," IEEE Transactions on Wireless Communications, vol. 12, no. 3, pp. 1418{1428, Mar. 2013. [55] G. Krieger and A. Moreira, \Spaceborne bi- and multistatic SAR: Potential and challenges," IEE Proceedings - Radar, Sonar and Navigation, vol. 153, no. 3, p. 184, 2006. [56] Wen-Qin Wang and Dingde Jiang, \Integrated Wireless Sensor Systems via Near-Space and Satellite Platforms: A Review," IEEE Sensors Journal, vol. 14, no. 11, pp. 3903{3914, Nov. 2014. [57] M. Younis, R. Metzig, and G. Krieger, \Performance Prediction of a Phase Synchronization Link for Bistatic SAR," IEEE Geoscience and Remote Sensing Letters, vol. 3, no. 3, pp. 429{433, Jul. 2006. 179 [58] Yang Yang and R. S. Blum, \Phase Synchronization for Coherent MIMO Radar: Algorithms and Their Analysis," IEEE Transactions on Signal Processing, vol. 59, no. 11, pp. 5538{5557, Nov. 2011. [59] S. Kong, S. Lee, C.-Y. Kim, and S. Hong, \Wireless Cooperative Synchronization of Coherent UWB MIMO Radar," IEEE Transactions on Microwave Theory and Techniques, vol. 62, no. 1, pp. 154{165, Jan. 2014. [60] S. Prager, G. Sexstone, D. McGrath, J. Fulton, and M. Moghaddam, \Snow Depth Re- trieval with an Autonomous UAV-mounted Software Dened Radar," IEEE Transactions on Geoscience and Remote Sensing, vol. Accepted for Publication, 2021. [61] F. T. Ulaby, R. K. Moore, A. K. Fung, and F. T. Ulaby, Microwave Remote Sensing, Fundamentals and Radiometry, ser. Microwave Remote Sensing. Norwood, Mass: ARTECH House, 1981, vol. 1. [62] F. Ulaby, F. F. Moore, Richard K., and A. K. Fung, Microwave Remote Sensing: Active and Passive. Norwood, MA: Artech House, 1986, vol. 2. [63] D. McGrath, R. Webb, D. Shean, R. Bonnell, H.-P. Marshall, T. H. Painter, N. P. Molotch, K. Elder, C. Hiemstra, and L. Brucker, \Spatially Extensive Ground-Penetrating Radar Snow Depth Observations During NASA's 2017 SnowEx Campaign: Comparison With In Situ, Airborne, and Satellite Observations," Water Resources Research, vol. 55, no. 11, pp. 10 026{10 036, Nov. 2019. [64] C. Matzler, \Microwave permittivity of dry snow," IEEE Transactions on Geoscience and Remote Sensing, vol. 34, no. 2, pp. 573{581, Mar. 1996. [65] P. Harder, J. W. Pomeroy, and W. D. Helgason, \Improving sub-canopy snow depth map- ping with unmanned aerial vehicles: Lidar versus structure-from-motion techniques," The Cryosphere, vol. 14, no. 6, pp. 1919{1935, Jun. 2020. [66] A. Tan, K. Eccleston, I. Platt, I. Woodhead, W. Rack, and J. McCulloch, \The design of a UAV mounted snow depth radar: Results of measurements on Antarctic sea ice," in 2017 IEEE Conference on Antenna Measurements & Applications (CAMA). Tsukuba: IEEE, Dec. 2017, pp. 316{319. [67] R. O. R. Jenssen, M. Eckerstorfer, S. K. Jacobsen, and R. Storvold, \DRONE-MOUNTED UWB RADAR SYSTEM FOR MEASURING SNOWPACK PROPERTIES: TECHNICAL IMPLEMENTATION, SPECIFICATIONS AND INITIAL RESULTS," p. 5. [68] R. O. R. Jenssen, M. Eckerstorfer, and S. Jacobsen, \Drone-Mounted Ultrawideband Radar for Retrieval of Snowpack Properties," IEEE Transactions on Instrumentation and Measurement, vol. 69, no. 1, pp. 221{230, Jan. 2020. [69] R. O. R. Jenssen and S. Jacobsen, \Drone-mounted UWB snow radar: Technical improvements and eld results," Journal of Electromagnetic Waves and Applications, vol. 34, no. 14, pp. 1930{1954, Sep. 2020. [70] A. E.-C. Tan, J. McCulloch, W. Rack, I. Platt, and I. Woodhead, \Radar Measurements of Snow Depth Over Sea Ice on an Unmanned Aerial Vehicle," IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, vol. 59, no. 3, p. 8, 2021. 180 [71] E. Arnold, C. Leuschen, F. Rodriguez-Morales, J. Li, J. Paden, R. Hale, and S. Keshmiri, \CReSIS airborne radars and platforms for ice and snow sounding," Annals of Glaciology, vol. 61, no. 81, pp. 58{67, Apr. 2020. [72] M. Stockham, J. Macy, and D. Besson, \Radio frequency ice dielectric permittivity measure- ments using CReSIS data: PERMITTIVITY MEASUREMENTS," Radio Science, vol. 51, no. 3, pp. 194{212, Mar. 2016. [73] B. Panzer, D. Gomez-Garcia, C. Leuschen, J. Paden, F. Rodriguez-Morales, A. Patel, T. Markus, B. Holt, and P. Gogineni, \An ultra-wideband, microwave radar for measur- ing snow thickness on sea ice and mapping near-surface internal layers in polar rn," Journal of Glaciology, vol. 59, no. 214, pp. 244{254, 2013. [74] C. Leuschen, R. Hale, S. Keshmiri, J. Yan, F. Rodriguez-Morales, A. Mahmood, and S. Gogi- neni, \UAS-Based Radar Sounding of the Polar Ice Sheets," IEEE Geoscience and Remote Sensing Magazine, vol. 2, no. 1, pp. 8{17, Mar. 2014. [75] A. Lundberg, N. Granlund, and D. Gustafsson, \Towards automated `Ground truth' snow measurements-a review of operational and new measurement methods for Sweden, Norway, and Finland," Hydrological Processes, pp. n/a{n/a, 2010. [76] W. S. Holbrook, S. N. Miller, and M. A. Provart, \Estimating snow water equivalent over long mountain transects using snowmobile-mounted ground-penetrating radar," GEOPHYSICS, vol. 81, no. 1, pp. WA183{WA193, Jan. 2016. [77] D. McGrath, L. Sass, S. O'Neel, C. McNeil, S. G. Candela, E. H. Baker, and H.-P. Marshall, \Interannual snow accumulation variability on glaciers derived from repeat, spatially extensive ground-penetrating radar surveys," The Cryosphere, vol. 12, no. 11, pp. 3617{3633, Nov. 2018. [78] R. W. Webb, \Using ground penetrating radar to assess the variability of snow water equivalent and melt in a mixed canopy forest, Northern Colorado," Frontiers of Earth Science, vol. 11, no. 3, pp. 482{495, Sep. 2017. [79] R. W. Webb, K. S. Jennings, M. Fend, and N. P. Molotch, \Combining Ground-Penetrating Radar With Terrestrial LiDAR Scanning to Estimate the Spatial Distribution of Liquid Water Content in Seasonal Snowpacks," Water Resources Research, vol. 54, no. 12, Dec. 2018. [80] N. Griessinger, F. Mohr, and T. Jonas, \Measuring snow ablation rates in alpine terrain with a mobile multioset ground-penetrating radar system," Hydrological Processes, Sep. 2018. [81] C. Mitterer, A. Heilig, J. Schweizer, and O. Eisen, \Upward-looking ground-penetrating radar for measuring wet-snow properties," Cold Regions Science and Technology, vol. 69, no. 2-3, pp. 129{138, Dec. 2011. [82] M. Durand, C. Gatebe, E. Kim, N. Molotch, T. H. Painter, M. Raleigh, M. Sandells, and C. Vuyovich, \NASA SnowEx Science Plan: Assessing Approaches for Measuring Water in Earth's Seasonal Snow," National Aeronautics and Space Administration, p. 68, 2017. [83] R. Kwok, B. Panzer, C. Leuschen, S. Pang, T. Markus, B. Holt, and S. Gogineni, \Airborne surveys of snow depth over Arctic sea ice," Journal of Geophysical Research, vol. 116, no. C11, p. C11018, Nov. 2011. 181 [84] J. Zhu, S. Tan, J. King, C. Derksen, J. Lemmetyinen, and L. Tsang, \Forward and Inverse Radar Modeling of Terrestrial Snow Using SnowSAR Data," IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 12, pp. 7122{7132, Dec. 2018. [85] H.-P. Marshall, M. Schneebeli, and G. Koh, \Snow stratigraphy measurements with high- frequency FMCW radar: Comparison with snow micro-penetrometer," Cold Regions Science and Technology, vol. 47, no. 1-2, pp. 108{117, Jan. 2007. [86] J.-B. Yan, D. Gomez-Garcia Alvestegui, J. W. McDaniel, Y. Li, S. Gogineni, F. Rodriguez- Morales, J. Brozena, and C. J. Leuschen, \Ultrawideband FMCW Radar for Airborne Mea- surements of Snow Over Sea Ice and Land," IEEE Transactions on Geoscience and Remote Sensing, vol. 55, no. 2, pp. 834{843, Feb. 2017. [87] N. Yankielun, W. Rosenthal, and R. E. Davis, \Alpine snow depth measurements from aerial FMCW radar," Cold Regions Science and Technology, vol. 40, no. 1-2, pp. 123{134, Nov. 2004. [88] M. Nolan, C. Larsen, and M. Sturm, \Mapping snow depth from manned aircraft on landscape scales at centimeter resolution using structure-from-motion photogrammetry," The Cryosphere, vol. 9, no. 4, pp. 1445{1463, Aug. 2015. [89] J. Goetz and A. Brenning, \Quantifying Uncertainties in Snow Depth Mapping From Structure From Motion Photogrammetry in an Alpine Area," Water Resources Research, vol. 55, no. 9, pp. 7772{7783, Sep. 2019. [90] P. D. Broxton and W. J. D. van Leeuwen, \Structure from Motion of Multi-Angle RPAS Imagery Complements Larger-Scale Airborne Lidar Data for Cost-Eective Snow Monitoring in Mountain Forests," Remote Sensing, vol. 12, no. 14, p. 2311, Jul. 2020. [91] S. Prager, W.-M. Shen, M. Moghaddam, and S. Masri, \Autonomous UAV-based Sensing of Landmines using Synthetic Ultra-Wideband Software Dened Radar," Manuscript in Preparation. [92] G. Tesfamariam and D. Mali, \GPR Technologies for Landmine Detection," TECHNIA - International Journal of Computing Science and Communication Technologies, vol. 5, no. 1, Jul. 2012. [93] D. J. Daniels, \A review of GPR for landmine detection," Sensing and Imaging: An Interna- tional Journal, vol. 7, no. 3, pp. 90{123, Sep. 2006. [94] ||, \Ground penetrating radar for buried landmine and IED detection," in Unexploded Ordnance Detection and Mitigation. Springer, 2009, pp. 89{111. [95] M. R. P. Cerquera, J. D. C. Monta~ no, and I. Mondrag on, \UAV for Landmine Detection Using SDR-Based GPR Technology," in Robots Operating in Hazardous Environments, H. Canbolat, Ed. InTech, Dec. 2017. [96] R. Bamler, \A comparison of range-Doppler and wavenumber domain SAR focusing algorithms," IEEE Transactions on Geoscience and Remote Sensing, vol. 30, no. 4, pp. 706{713, Jul. 1992. [97] L. A. Gorham and L. J. Moore, \SAR image formation toolbox for MATLAB," E. G. Zelnio and F. D. Garber, Eds., Apr. 2010, p. 769906. 182 [98] M. Albahkali and M. Moghaddam, \3D SAR focusing for subsurface point targets," in 2009 IEEE International Geoscience and Remote Sensing Symposium. Cape Town, South Africa: IEEE, 2009, pp. I{60{I{63. [99] T. Jin and Z. Zhou, \Refraction and Dispersion Eects Compensation for UWB SAR Subsur- face Object Imaging," IEEE Transactions on Geoscience and Remote Sensing, vol. 45, no. 12, pp. 4059{4066, Dec. 2007. [100] J. Fortuny-Guasch, \A novel 3-D subsurface radar imaging technique," IEEE Transactions on Geoscience and Remote Sensing, vol. 40, no. 2, pp. 443{452, Feb./2002. [101] M. Moghaddam, D. Entekhabi, Y. Goykhman, K. Li, M. Liu, A. Mahajan, A. Nayyar, D. Shuman, and D. Teneketzis, \A Wireless Soil Moisture Smart Sensor Web Using Physics- Based Optimal Control: Concept and Initial Demonstrations," IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 3, no. 4, pp. 522{535, Dec. 2010. [102] D. Clewley, J. B. Whitcomb, R. Akbar, A. R. Silva, A. Berg, J. R. Adams, T. Caldwell, D. Entekhabi, and M. Moghaddam, \A Method for Upscaling In Situ Soil Moisture Measure- ments to Satellite Footprint Scale Using Random Forests," IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 10, no. 6, pp. 2663{2673, Jun. 2017. [103] R. H. Chen, A. Tabatabaeenejad, and M. Moghaddam, \Analysis of Permafrost Active Layer Soil Heterogeneity in Support of Radar Retrievals," in IGARSS 2018 - 2018 IEEE International Geoscience and Remote Sensing Symposium. Valencia: IEEE, Jul. 2018, pp. 8559{8562. [104] A. Nayyar, \Sequential decision making in decentralized systems," PhD Thesis, University of Michigan, 2011. [105] J. Mittermayer, G. Krieger, and A. Moreira, \Concepts and Applications of Multi-static MirrorSAR Systems," in 2020 IEEE Radar Conference (RadarConf20). Florence, Italy: IEEE, Sep. 2020, pp. 1{6. [106] G. Krieger, A. Moreira, H. Fiedler, I. Hajnsek, M. Werner, M. Younis, and M. Zink, \TanDEM- X: A Satellite Formation for High-Resolution SAR Interferometry," IEEE Transactions on Geoscience and Remote Sensing, vol. 45, no. 11, pp. 3317{3341, Nov. 2007. [107] E. Research, \USRP E312 Data Sheet," https://www.ettus.com/wp- content/uploads/2019/01/USRP E312 Datasheet.pdf, 2019. [108] J. Malsbury and M. Ettus, \Simplifying FPGA design with a novel network-on-chip archi- tecture," in Proceedings of the Second Workshop on Software Radio Implementation Forum - SRIF '13. Hong Kong, China: ACM Press, 2013, p. 45. [109] A. T. Lab, \Antennas by RFSPACE," https://antennatestlab.com/antenna-examples/example- 3-vivaldi-antennas-rfspace, Feb. 2017. [110] D. J. Rabideau, \Nonlinear synthetic wideband waveforms," in Radar Conference, 2002. Proceedings of the IEEE. IEEE, 2002, pp. 212{219. 183 [111] S. Prager, T. Thrivikraman, M. Haynes, J. Stang, D. Hawkins, and M. Moghaddam, \Ultra- wideband synthesis for high-range resolution software dened radar," in 2018 IEEE Radar Conference (RadarConf18). Oklahoma City, OK: IEEE, Apr. 2018, pp. 1089{1094. [112] W.-B. Gao, T. Long, Z.-G. Ding, and Y.-R. Wu, \A Robust Range Grating Lobe Suppression Method Based on Image Contrast for Stepped-Frequency SAR," Sensors, vol. 16, no. 12, p. 2066, Dec. 2016. [113] S. Prager, \Igarss 2020," https://github.com/samprager/igarss2020. [114] J. F. Kaiser and R. W. Schafer, \On the use of the I0-sinh window for spectrum analysis," IEEE Transactions on Acoustics, Speech, and Signal Processing, pp. 105{107, 1980. [115] S. V. Tsynkov, \On the Use of Start-Stop Approximation for Spaceborne SAR Imaging," SIAM Journal on Imaging Sciences, vol. 2, no. 2, pp. 646{669, Jan. 2009. [116] R. Bamler, \A comparison of range-Doppler and wavenumber domain SAR focusing algorithms," IEEE Transactions on Geoscience and Remote Sensing, vol. 30, no. 4, pp. 706{713, Jul. 1992. [117] M. Abramowitz and I. A. Stegun, \Solutions of quartic equations," in Handbook of Mathe- matical Functions with Formulas, Graphs, and Mathematical Tables, 9th ed. New York NY: Dover, 1972, ch. 3.8.3, pp. 17{18. [118] C. Warren, A. Giannopoulos, and I. Giannakis, \An advanced gpr modelling framework: The next generation of gprmax," in Advanced Ground Penetrating Radar (IWAGPR), 2015 8th International Workshop On. IEEE, 2015, pp. 1{4. [119] F. Roth, P. van Genderen, and M. Verhaegen, \RADAR RESPONSE APPROXIMATIONS FOR BURIED PLASTIC LANDMINES," in SPIE, vol. 4758. SPIE, 2002. [120] Mei Leng and Yik-Chung Wu, \Distributed Clock Synchronization for Wireless Sensor Networks Using Belief Propagation," IEEE Transactions on Signal Processing, vol. 59, no. 11, pp. 5404{ 5414, Nov. 2011. [121] J. Du and Y.-C. Wu, \Distributed Clock Skew and Oset Estimation in Wireless Sensor Networks: Asynchronous Algorithm and Convergence Analysis," IEEE Transactions on Wireless Communications, vol. 12, no. 11, pp. 5908{5917, Nov. 2013. [122] B. Ananthasubramaniam and U. Madhow, \On Localization Performance in Imaging Sensor Nets," IEEE Transactions on Signal Processing, vol. 55, no. 10, pp. 5044{5057, Oct. 2007. [123] J. Liang and Q. Liang, \Design and Analysis of Distributed Radar Sensor Networks," IEEE Transactions on Parallel and Distributed Systems, vol. 22, no. 11, pp. 1926{1933, Nov. 2011. [124] S. Lanzisera and K. S. J. Pister, \Burst Mode Two-Way Ranging with Cramer-Rao Bound Noise Performance," in IEEE GLOBECOM 2008 - 2008 IEEE Global Telecommunications Conference. New Orleans, LA, USA: IEEE, 2008, pp. 1{5. [125] M. Gardill, J. Schwendner, and J. Fuchs, \An Approach to Over-the-air Synchronization of Commercial Chirp-Sequence Automotive Radar Sensors," in 2020 IEEE Topical Conference on Wireless Sensors and Sensor Networks (WiSNeT). San Antonio, TX, USA: IEEE, Jan. 2020, pp. 46{49. 184 [126] R. Moddemeijer, \On the determination of the position of extrema of sampled correlators," IEEE Transactions on Signal Processing, vol. 39, no. 1, pp. 216{219, Jan. 1991. [127] G. Jacovitti and G. Scarano, \Discrete time techniques for time delay estimation," IEEE Transactions on Signal Processing, vol. 41, no. 2, pp. 525{533, Feb./1993. [128] I. Jameson, \Time Delay Estimation," Electronic Warfare and Radar Division Defence Science and Technology Organisation, Edinburgh, South Australia, Tech. Rep. DSTO{TR{1705, 2006. [129] Y. Schroder, D. Reimers, and L. Wolf, \Accurate and Precise Distance Estimation from Phase-Based Ranging Data," in 2018 International Conference on Indoor Positioning and Indoor Navigation (IPIN). Nantes: IEEE, Sep. 2018, pp. 1{8. [130] E. Weinstein and A. Weiss, \Fundamental limitations in passive time-delay estimation{Part II: Wide-band systems," IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 32, no. 5, pp. 1064{1078, Oct. 1984. [131] A. Oppenheim and A. Willsky, Signals and Systems. Upper Saddle River, N.J: Prentice Hall, 1997. [132] I. Sharp, K. Yu, and Y. Guo, \Peak and leading edge detection for time-of-arrival estimation in band-limited positioning systems," IET Communications, vol. 3, no. 10, p. 1616, 2009. [133] Richard A. Caruana, Roger B. Searle, Thomas. Heller, and Saul I. Shupack, \Fast algorithm for the resolution of spectra," Analytical Chemistry, vol. 58, no. 6, pp. 1162{1167, 1986. [134] H. Guo, \A Simple Algorithm for Fitting a Gaussian Function [DSP Tips and Tricks]," IEEE Signal Processing Magazine, vol. 28, no. 5, pp. 134{137, Sep. 2011. [135] J. Yan, C. C. J. M. Tiberius, G. J. M. Janssen, P. J. G. Teunissen, and G. Bellusci, \Review of range-based positioning algorithms," IEEE Aerospace and Electronic Systems Magazine, vol. 28, no. 8, pp. 2{27, Aug. 2013. [136] K. Yu and Y. Guo, \Anchor-free localisation algorithm and performance analysis in wireless sensor networks," IET Communications, vol. 3, no. 4, p. 549, 2009. [137] N. Patwari, A. Hero, M. Perkins, N. Correal, and R. O'Dea, \Relative location estimation in wireless sensor networks," IEEE Transactions on Signal Processing, vol. 51, no. 8, pp. 2137{2148, Aug. 2003. [138] A. Cakiroglu and C. Erten, \Fully Decentralized and Collaborative Multilateration Primi- tives for Uniquely Localizing WSNs," EURASIP Journal on Wireless Communications and Networking, vol. 2010, pp. 1{7, 2010. [139] B. Hawkins, M. Anderson, S. Prager, S.-J. Chung, and M. Lavalle, \Experiments with small UAS to support SAR tomographic mission formulation," in IGARSS 2021 - 2021 IEEE International Geoscience and Remote Sensing Symposium, Brussels, Belgium, 2021. [140] S. Prager, M. Moghaddam, and M. Lavalle, \Characterization of Clock Phase Errors for Distributed Wireless Synchronization Protocol," in IGARSS 2021 - 2021 IEEE International Geoscience and Remote Sensing Symposium, Brussels, Belgium, 2021. [141] V. Kroupa, Frequency Stability Introduction and Applications. Wiley-IEEE Press, Sep. 2012. 185 [142] G. Krieger and M. Younis, \Impact of Oscillator Noise in Bistatic and Multistatic SAR," IEEE Geoscience and Remote Sensing Letters, vol. 3, no. 3, pp. 424{428, Jul. 2006. [143] J. A. Barnes, A. R. Chi, L. S. Cutler, D. J. Healey, D. B. Leeson, T. E. McGunigal, J. A. Mullen, W. L. Smith, R. L. Sydnor, R. F. C. Vessot, and G. M. R. Winkler, \Characterization of Frequency Stability," IEEE Transactions on Instrumentation and Measurement, vol. 20, no. 2, pp. 105{120, May 1971. [144] J. Rutman, \Characterization of phase and frequency instabilities in precision frequency sources: Fifteen years of progress," Proceedings of the IEEE, vol. 66, no. 9, pp. 1048{1075, 1978. [145] R. C. Bales, N. P. Molotch, T. H. Painter, M. D. Dettinger, R. Rice, and J. Dozier, \Mountain hydrology of the western United States: MOUNTAIN HYDROLOGY OF THE WESTERN US," Water Resources Research, vol. 42, no. 8, Aug. 2006. [146] J. I. L opez-Moreno, J. Revuelto, S. R. Fassnacht, C. Azor n-Molina, S. M. Vicente-Serrano, E. Mor an-Tejeda, and G. A. Sexstone, \Snowpack variability across various spatio-temporal resolutions: SNOWPACK VARIABILITY ACROSS VARIOUS SPATIO-TEMPORAL RESO- LUTIONS," Hydrological Processes, vol. 29, no. 6, pp. 1213{1224, Mar. 2015. [147] K. Elder, J. Dozier, and J. Michaelsen, \Snow accumulation and distribution in an Alpine Watershed," Water Resources Research, vol. 27, no. 7, pp. 1541{1552, Jul. 1991. [148] DJI, \DJI M600 Pro Specs," https://www.dji.com/matrice600-pro/info#specs, 2021. [149] Sensors and Software, \PulseEKKO GPR," https://www.sensoft.ca/wp- content/uploads/2020/08/pulseEKKO-Brochure-and-Ultra-Receiver.pdf, 2021. [150] E. Ferraro and C. Swift, \Comparison of retracking algorithms using airborne radar and laser altimeter measurements of the Greenland ice sheet," IEEE Transactions on Geoscience and Remote Sensing, vol. 33, no. 3, pp. 700{707, May 1995. [151] S. Hendricks, L. Stenseng, V. Helm, and C. Haas, \Eects of surface roughness on sea ice freeboard retrieval with an Airborne Ku-Band SAR radar altimeter," in 2010 IEEE International Geoscience and Remote Sensing Symposium. Honolulu, HI, USA: IEEE, Jul. 2010, pp. 3126{3129. [152] D. Hutchison, T. Kanade, J. Kittler, J. M. Kleinberg, F. Mattern, J. C. Mitchell, M. Naor, O. Nierstrasz, C. Pandu Rangan, B. Steen, M. Sudan, D. Terzopoulos, D. Tygar, M. Y. Vardi, G. Weikum, S. Alpert, M. Galun, B. Nadler, and R. Basri, \Detecting Faint Curved Edges in Noisy Images," in Computer Vision { ECCV 2010, K. Daniilidis, P. Maragos, and N. Paragios, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2010, vol. 6314, pp. 750{763. [153] N. Or, M. Galun, S. Alpert, A. Brandt, B. Nadler, and R. Basri, \On Detection of Faint Edges in Noisy Images," arXiv:1706.07717 [cs], Jun. 2017. [154] S. Xie and Z. Tu, \Holistically-Nested Edge Detection," arXiv:1504.06375 [cs], Oct. 2015. [155] F. Ulaby, R. Moore, and A. Fung, Microwave Remote Sensing: Active and Passive. Artech House, 1986, vol. 3. 186 [156] L. Jiao, L. Huo, C. Hu, and P. Tang, \Rened UNet: UNet-Based Renement Network for Cloud and Shadow Precise Segmentation," Remote Sensing, vol. 12, no. 12, p. 2001, Jun. 2020. [157] Z. Chu, T. Tian, R. Feng, and L. Wang, \Sea-Land Segmentation With Res-UNet And Fully Connected CRF," in IGARSS 2019 - 2019 IEEE International Geoscience and Remote Sensing Symposium. Yokohama, Japan: IEEE, Jul. 2019, pp. 3840{3843. [158] F. Mazeh, B. Hammoud, H. Ayad, F. Ndagijimana, G. Faour, M. Fadlallah, and J. Jo- maah, \NUMERICAL ANALYSIS OF RADAR RESPONSE TO SNOW USING MULTIPLE BACKSCATTERING MEASUREMENTS FOR SNOW DEPTH RETRIEVAL," Progress In Electromagnetics Research B, vol. 81, pp. 63{80, 2018. [159] C. Grima, I. Koch, J. S. Greenbaum, K. M. Soderlund, D. D. Blankenship, D. A. Young, D. M. Schroeder, and S. Fitzsimons, \Surface and basal boundary conditions at the Southern McMurdo and Ross Ice Shelves, Antarctica," Journal of Glaciology, vol. 65, no. 252, pp. 675{688, Aug. 2019. [160] P. Hoekstra and D. Spanogle, \Backscatter from snow and ice surfaces at near incident angles," IEEE Transactions on Antennas and Propagation, vol. 20, no. 6, pp. 788{790, Nov. 1972. [161] M. Moghaddam, R. Akbar, S. Prager, A. Silva, and D. Entekhabi, \SPCTOR: SENSING POLICY CONTROLLER AND OPTIMIZER," in IGARSS 2020 - 2020 IEEE International Geoscience and Remote Sensing Symposium, 2020, p. 4. [162] S. M. Kay, Fundamentals of Statistical Signal Processing: Estimation Theory. Upper Saddle River, NJ, USA: Prentice-Hall, Inc., 1993, vol. 1. [163] L. L. Scharf, Statistical Signal Processing: Detection, Estimation, and Time Series Analysis. NY: Addison Wesley. 187

Abstract (if available)

Abstract Multistatic and multiple input multiple output (MIMO) radar sensor networks are an area of signiﬁcant interest due to spatial diversity gains, improved resolution, and higher dynamic range as well the potential for these systems to reduce overall mission costs and improve survivability. Networks of autonomous small unmanned aerial vehicles (UAVs) and rovers, or unmanned ground vehicles (UGVs), and small CubeSat constellations represent exciting potential as deployment platforms for distributed radar systems. However, the small payload of capacity of such platforms and large number of individual sensors that constitue a distributed network necessitates the development of small, low cost, and eﬃcient radar sensors that are capable of high levels of performance. A compelling application of such sensor networks is in payload directed autonomy and dynamic reconﬁgurability for different radar sensing modes and objectives, requiring sensors that are also highly ﬂexible and intelligent. Furthermore, in order achieve coherent multistatic/MIMO radar operation in a wireless sensor network, precise synchronization to sub-100 picosecond levels between nodes is required and represents a fundamental obstacle to the realization of such systems. ❧ This research focuses on the development of software deﬁned radar (SDRadar) algorithms and techniques that address two primary problems in this area. We ﬁrst present a synthetic ultra-wideband waveform and reconstruction technique that that allows low-cost software deﬁned radio (SDR)-based SDRadars to eﬃciently achieve cm-level radar resolution from real-time dynamically conﬁgurable spectrum allocations. Next, we report a novel O(N) distributed and decentralized wireless synchronization scheme that achieves time synchronization to sub-100 picosecond precision across all nodes and enables time and phase coherent multistatic/MIMO radar operation. ❧ The SDRadar system and techniques presented in this work are validated extensively through experimental results. We ﬁrst demonstrate the performance of the synthetic wideband waveform reconstruction techniques and the ﬂexibility of the implemented synthetic wideband SDRadar in multiple radar imaging modes, including synthetic aperture radar (SAR) and ground penetrating radar (GPR). We then report results from multiple SDRadar synchronization experiments for both two and three sensor networks, demonstrating sub-100 picosecond performance, cm-level localization, and the ability of the method to coherently synchronize sensors across multiple stepped frequency bands for bistatic synthetic wideband radar imaging. We derive an analytic formulation of the synchronization protocol in terms of the sensor input oscillator’s phase error power spectral density (PSD) and validate the model with further experimental results. Finally, we describe the integration of the SDRadar as a UAV payload and share experimental results for the UAV-SDRadar system in autonomous UAV-based radar imaging from multiple ﬁeld campaigns.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Physics-based bistatic radar scattering model for vegetated terrains in support of soil moisture retrieval from signals of opportunity

PDF

Real-time channel sounder designs for millimeter-wave and ultra-wideband communications

PDF

Improving spectrum efficiency of 802.11ax networks

PDF

Signal processing for channel sounding: parameter estimation and calibration

PDF

Radar remote sensing of permafrost landscapes and active layer soil properties

PDF

Achieving high data rates in distributed MIMO systems

PDF

FPGA based L-band pulse Doppler radar design and implementation

PDF

Design and analysis of reduced complexity transceivers for massive MIMO and UWB systems

PDF

Magnetic induction-based wireless body area network and its application toward human motion tracking

PDF

Physics-based models and microwave sensors for mapping soil carbon in arctic permafrost using long-wavelength radars

PDF

Hybrid beamforming for massive MIMO

PDF

Multidimensional characterization of propagation channels for next-generation wireless and localization systems

PDF

All-optical signal processing toward reconfigurable optical networks

PDF

Solution of inverse scattering problems via hybrid global and local optimization

PDF

Structured codes in network information theory

PDF

Multichannel data collection for throughput maximization in wireless sensor networks

PDF

Propagation channel characterization and interference mitigation strategies for ultrawideband systems

PDF

Design and analysis of large scale antenna systems

PDF

Robust routing and energy management in wireless sensor networks

PDF

Optimal distributed algorithms for scheduling and load balancing in wireless networks

Asset Metadata

Creator Prager, Samuel M. (author)

Core Title Ultra-wideband multistatic and MIMO software defined radar sensor networks

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Electrical Engineering

Degree Conferral Date 2021-12

Publication Date 10/28/2021

Defense Date 08/20/2021

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag distributed radar,MIMO radar,multistatic radar,OAI-PMH Harvest,software defined radar,synthetic aperture radar,wireless synchronization

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Moghaddam, Mahta (committee chair), Chugg, Keith (committee member), Masri, Sami (committee member), Molisch, Andreas (committee member), Stang, John (committee member)

Creator Email prager.sam@gmail.com,sprager@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-oUC16251450

Unique identifier UC16251450

Legacy Identifier etd-PragerSamu-10177

Document Type Dissertation

Format application/pdf (imt)

Rights Prager, Samuel M.

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright. The original signature page accompanying the original submission of the work to the USC Libraries is retained by the USC Libraries and a copy of it may be obtained by authorized requesters contacting the repository e-mail address given.

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA

Repository Email cisadmin@lib.usc.edu