The Application of a Particle Filter to Urban Ground Target Localization, Tracking, and Intercept by Emily A. Doucette A dissertation submitted to the Graduate Faculty of Auburn University in partial ful llment of the requirements for the Degree of Doctor of Philosophy Auburn, Alabama May 7, 2012 Keywords: Particle Filter, Urban Environment, Path Planning Copyright 2012 by Emily A. Doucette Approved by Andrew J. Sinclair, Chair, Associate Professor of Aerospace Engineering John E. Cochran, Jr., Professor and Department Head of Aerospace Engineering David A. Cicci, Professor of Aerospace Engineering Gilbert L. Crouse, Associate Professor of Aerospace Engineering George T. Flowers, Dean of Graduate School Abstract Ground vehicle localization is a problem of signi cance in an urban setting given the recent global con icts, security interests, and rapid growth of sensor networks. This task is often di cult due to the lack of information regarding a target vehicle?s position, velocity, and destination. Field operatives can provide binary measurements of a target?s presence in an area, and these measurements can be processed to obtain estimates of the target?s location. A particle lter is more suitable for this application than a Kalman lter due to its ability to handle non-Gaussian distributions and non-di erentiable measurement models, however it is computationally expensive. Suppose there is a mobile ground vehicle of known description but unknown position, velocity, or destination that is to be found, tracked, and intercepted by an unmanned aerial vehicle. The vehicle is known to be in an urban environment, and full knowledge of that environment (roads, obstacles, intersection constraints, and speed limits) is available. There are numerous issues of interest within this problem. A particle lter in an urban environment was developed to locate, track, and intercept a ground vehicle given soft binary measurements (measurements from human sources). Two particular issues are studied in this work: the e ect of a sophisticated particle dynamic model on target localization and tracking, and the development of a real time path planning routine in the particle lter framework to enable target interception. The contributions of this work are threefold. First, the importance and impact of an accurate particle time update on target localization and tracking is validated. Secondly, a thorough investigation into the e ect of particle spatial resolution in the presence of imperfect measurements is made that will prove valuable for future particle lter applications. Finally, the path planning routine o ers reduced computational expense when compared to existing ii systems and lends itself to unmanned aerial vehicle implementation. Proper exploitation and implementation of the particle lter framework prove vital in the complete characterization of the urban tracking problem. iii Acknowledgments The author would like to acknowledge and sincerely thank Dr. Andrew J. Sinclair for his guidance, support, patience, and encouragement throughout her time in graduate school. The author would also like to thank the SMART Scholar Program, Dr. J. Willard Curtis, and the Munitions Directorate of Air Force Research Laboratory for the opportunity to investigate this work. Also, the author sincerely appreciates the faculty of the Aerospace Engineering Department at Auburn University, namely Dr. R. Steven Gross, for many years of invaluable advisement and instruction and for all teaching opportunities granted. Personal gratitude is bestowed upon the sisters and alumnae of the Zeta Xi Chapter of Alpha Xi Delta for the opportunity to realize lifelong friendships, leadership skills, and potential. The author expresses heartfelt thanks to her dear family and loyal friends for their constant encouragement and unwavering support. The powerful prayers and outpouring of love provided by late grandmother Imelda S. Gares and the strength and conviction to remain true to one?s self provided by grandmother Molly M. Waller have always enabled the author to overcome any obstacle. This work is dedicated to the author?s beloved parents, Judith C. Doucette and the late Roland Doucette Jr., whose strong faith, self sacri ce, rm belief in hard work, dedica- tion to one?s best, encouragement, kindness, humor, and love are ever present, unmatched, inspirational, and forever held dear. iv Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2 The Nonlinear Dynamic System Filter Problem . . . . . . . . . . . . . . . . . . 4 2.1 The Kalman Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1.1 The Extended Kalman Filter . . . . . . . . . . . . . . . . . . . . . . 8 2.2 The Particle Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.3 Entropy of a Probability Density Function . . . . . . . . . . . . . . . . . . . 14 2.4 Applications to Target Localization and Tracking . . . . . . . . . . . . . . . 18 3 The E ect of the Particle Motion Model . . . . . . . . . . . . . . . . . . . . . . 21 3.1 Problem Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.1.1 Target Motion Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.1.2 Measurement Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.2 Dispersion Motion Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.3 Tra c Motion Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.4 Measurement Update . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.4.1 Particle Weight Update . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.4.2 Resample Step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.5 Motion Model Comparison Results . . . . . . . . . . . . . . . . . . . . . . . 42 3.6 Importance of Spatial Resolution . . . . . . . . . . . . . . . . . . . . . . . . 46 3.6.1 Particle Redistribution . . . . . . . . . . . . . . . . . . . . . . . . . . 47 v 3.6.2 False Rate Inclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.6.3 Reduced Number of E ective Particles . . . . . . . . . . . . . . . . . 57 3.7 Summary of Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4 Accounting for an Imperfect Sensor . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.1 Sensor Belief Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.2 Number of E ective Particles Study . . . . . . . . . . . . . . . . . . . . . . . 84 4.3 Bayesian Risk . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 4.3.1 Risk Assessment Results . . . . . . . . . . . . . . . . . . . . . . . . . 102 5 Target Intercept . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 5.1 UAV Path Following . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 5.2 UAV Measurement Update . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 5.3 Candidate Path Formation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 5.4 Path Selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 5.4.1 Maximum Likelihood . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 5.4.2 Information Utility Function . . . . . . . . . . . . . . . . . . . . . . . 137 5.5 Method Comparison Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150 A Expected Travel Time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 vi List of Figures 2.1 The Resampling Process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2 E ect of Resampling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.3 (a) Low density pdf and (b) High Density pdf . . . . . . . . . . . . . . . . . . . 15 2.4 Piecewise linear approximation of a PDF . . . . . . . . . . . . . . . . . . . . . . 16 3.1 Map of Urban Environment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.2 Geometry of a Simple Car . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.3 (a) Straight Path and (b) Winding Path . . . . . . . . . . . . . . . . . . . . . . 25 3.4 Determine Region of Highest Particle Weight . . . . . . . . . . . . . . . . . . . 26 3.5 Measure Region of Highest Particle Weight . . . . . . . . . . . . . . . . . . . . 27 3.6 Geometry of Unicycle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.7 Intersection of Two Roads . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.8 U-turn course adjustment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.9 Geometric characterization of an intersection . . . . . . . . . . . . . . . . . . . 32 3.10 Example of a Required Lane Change . . . . . . . . . . . . . . . . . . . . . . . . 33 3.11 De nition of q Variable . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 vii 3.12 Geometry During a Sample Turn . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.13 Geometric Turning Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.14 Determine Region of Highest Particle Weight . . . . . . . . . . . . . . . . . . . 38 3.15 Measurement Update Example Scenario . . . . . . . . . . . . . . . . . . . . . . 39 3.16 Post-measurement Update for Example Scenario . . . . . . . . . . . . . . . . . 41 3.17 Resampled Distribution for Example Scenario . . . . . . . . . . . . . . . . . . . 42 3.18 Average Mean Square Error After Target Found, Straight Path . . . . . . . . . 45 3.19 Average Mean Square Error After Target Found, Winding Path . . . . . . . . . 45 3.20 (a) Time = 0 sec. (b) Time = 200 sec. . . . . . . . . . . . . . . . . . . . . . . . 46 3.21 Redistribution Region . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.22 Average Mean Square Error After Target Found, Straight Path, Tra c Motion Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.23 Average Mean Square Error After Target Found, Straight Path, Tra c Motion Model, Final 30 seconds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.24 Average Mean Square Error After Target Found, Winding Path, Tra c Motion Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.25 Average Mean Square Error After Target Found, Winding Path, Tra c Motion Model, Final 30 seconds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.26 E ect of false report rate belief on particle spatial resolution and weight distribution 54 3.27 E ect of false report rate belief on particle spatial resolution and weight distribution 55 viii 3.28 Average Mean Square Error After Target Found, Straight Path, 90% Con dence 56 3.29 Average Mean Square Error After Target Found, Winding Path, 90% Con dence 57 3.30 E ect of the number of e ective particles on particle spatial resolution, Ne > 90%N . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.31 E ect of the number of e ective particles on particle spatial resolution, Ne > 90%N . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.32 E ect of the number of e ective particles on particle spatial resolution, Ne > 80%N . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.33 Dispersion Model, Average Mean Square Error After Target Found, Straight Path, 80%N and 90%N E ective Particle Thresholds . . . . . . . . . . . . . . . 62 3.34 Tra c Motion Model, Average Mean Square Error After Target Found, Straight Path, 80%N and 90%N E ective Particle Thresholds . . . . . . . . . . . . . . . 62 3.35 Dispersion Model, Average Mean Square Error After Target Found, Winding Path, 80%N and 90%N E ective Particle Thresholds . . . . . . . . . . . . . . . 63 3.36 Tra c Motion Model, Average Mean Square Error After Target Found, Winding Path, 80%N and 90%N E ective Particle Thresholds . . . . . . . . . . . . . . . 63 4.1 Impact of a Sensor with 10% False Report Rate, Tra c Model, Winding Path, N = 1000, False Report Belief 10% . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.2 Dispersion Model, Straight Path, False Report Belief Study, Perfect Sensor . . . 70 4.3 Tra c Model, Straight Path, False Report Belief Study, Perfect Sensor . . . . . 71 4.4 Tra c Model, Straight Path, False Report Belief Study, Perfect Sensor, nal 20 seconds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 ix 4.5 Dispersion Model, Winding Path, False Report Belief Study, Perfect Sensor . . . 73 4.6 Tra c Model, Winding Path, False Report Belief Study, Perfect Sensor . . . . . 73 4.7 Tra c Model, Winding Path, False Report Belief Study, Perfect Sensor, Final 20 seconds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.8 Dispersion Model, Straight Path, False Report Belief Study, 10% False Report Sensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.9 Tra c Model, Straight Path, False Report Belief Study, 10% False Report Sensor 77 4.10 Dispersion Model, Winding Path, False Report Belief Study, 10% False Report Sensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 4.11 Tra c Model, Winding Path, False Report Belief Study, 10% False Report Sensor 79 4.12 Dispersion Model, Straight Path, False Report Belief Study, 20% False Report Sensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.13 Tra c Model, Straight Path, False Report Belief Study, 20% False Report Sensor 81 4.14 Dispersion Model, Winding Path, False Report Belief Study, 20% False Report Sensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 4.15 Tra c Model, Winding Path, False Report Belief Study, 20% False Report Sensor 83 4.16 Dispersion Model, Straight Path, Number of E ective Particles Study, Perfect Sensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 4.17 Tra c Model, Straight Path, Number of E ective Particles Study, Perfect Sensor 85 4.18 Dispersion Model, Winding Path, Number of E ective Particles Study, Perfect Sensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 x 4.19 Tra c Model, Winding Path, Number of E ective Particles Study, Perfect Sensor 87 4.20 Dispersion Model, Straight Path, Number of E ective Particles Study, 10% False Report Sensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 4.21 Tra c Model, Straight Path, Number of E ective Particles Study, 10% False Report Sensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 4.22 Dispersion Model, Winding Path, Number of E ective Particles Study, 10% False Report Sensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 4.23 Tra c Model, Winding Path, Number of E ective Particles Study, 10% False Report Sensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 4.24 Dispersion Model, Straight Path, Number of E ective Particles Study, 20% False Report Sensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 4.25 Tra c Model, Straight Path, Number of E ective Particles Study, 20% False Report Sensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 4.26 Dispersion Model, Winding Path, Number of E ective Particles Study, 20% False Report Sensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 4.27 Tra c Model, Winding Path, Number of E ective Particles Study, 20% False Report Sensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 4.28 Measurement Update with Risk Assessment Example Scenario . . . . . . . . . . 101 4.29 E ect of Risk Assessment, Straight Path, 10% False Report Sensor . . . . . . . 104 4.30 E ect of Risk Assessment, Winding Path, 10% False Report Sensor . . . . . . . 105 4.31 E ect of Risk Assessment, Straight Path, 20% False Report Sensor . . . . . . . 106 xi 4.32 E ect of Risk Assessment, Winding Path, 20% False Report Sensor . . . . . . . 107 4.33 Dispersion Model Tracking Performance over 5 minutes, Straight Path . . . . . 109 4.34 Tra c Motion Model Tracking Performance over 5 minutes, Straight Path . . . 110 4.35 Dispersion Model Tracking Performance over 5 minutes, Winding Path . . . . . 111 4.36 Tra c Motion Model Tracking Performance over 5 minutes, Winding Path . . . 112 5.1 Vector eld for straight-line following . . . . . . . . . . . . . . . . . . . . . . . . 116 5.2 Vector elds for various values of k . . . . . . . . . . . . . . . . . . . . . . . . . 117 5.3 UAV Path Following via Vector Fields Example . . . . . . . . . . . . . . . . . . 119 5.4 UAV Sensor Requirements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 5.5 UAV sensor span while traveling down-road . . . . . . . . . . . . . . . . . . . . 121 5.6 UAV sensor tilt while traveling down-road . . . . . . . . . . . . . . . . . . . . . 122 5.7 Node Network . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 5.8 UAV planning horizon . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 5.9 Node Map . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 5.10 Example: Nodes and Roads for Candidate Paths . . . . . . . . . . . . . . . . . 127 5.11 Edge of Directed Graph from Root Node to Node 2 . . . . . . . . . . . . . . . . 128 5.12 Edges of Directed Graph from Root Node to Node 5 . . . . . . . . . . . . . . . 129 5.13 Edges of Directed Graph from Root Node to Node 6 . . . . . . . . . . . . . . . 130 xii 5.14 Directed Graph after nding Node 3 . . . . . . . . . . . . . . . . . . . . . . . . 131 5.15 Directed Graph after nding a second route to Node 6 . . . . . . . . . . . . . . 132 5.16 Directed Graph after nding Node 4 . . . . . . . . . . . . . . . . . . . . . . . . 133 5.17 Directed Graph of Example Path Planning Problem . . . . . . . . . . . . . . . . 134 5.18 Two Part Road Weight . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 xiii List of Tables 3.1 Time to Detect Results for Straight Path, Perfect Sensor . . . . . . . . . . . . . 43 3.2 Time to Detect Results for Winding Path, Perfect Sensor . . . . . . . . . . . . . 43 3.3 Computational Expense Comparison . . . . . . . . . . . . . . . . . . . . . . . . 44 3.4 Time to Detect Results for Straight Path, N = 1000 . . . . . . . . . . . . . . . 53 3.5 Time to Detect Results for Winding Path, N = 1000 . . . . . . . . . . . . . . . 56 3.6 Time to Detect Results for Straight Path, N = 1000 . . . . . . . . . . . . . . . 64 3.7 Time to Detect Results for Winding Path, N = 1000 . . . . . . . . . . . . . . . 64 4.1 Time to Detect Results for Winding Path, N = 1000, False Report Belief 10% . 69 4.2 Time to Detect Results, Straight Path, Perfect Sensor, False Report Belief 0%-50% 72 4.3 Time to Detect Results, Winding Path, Perfect Sensor, False Report Belief 0%-50% 75 4.4 Time to Detect Results, Straight Path, 10% False Report Sensor, False Report Belief 0%-50% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 4.5 Time to Detect Results, Winding Path, 10% False Report Sensor, False Report Belief 0%-50% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.6 Time to Detect Results, Straight Path, 20% False Report Sensor, False Report Belief 0%-50% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 4.7 Time to Detect Results, Winding Path, 20% False Report Sensor, False Report Belief 0%-50% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 4.8 Time to Detect Results, Straight Path, Perfect Sensor, Nthr = (100% 50%)N . 86 4.9 Time to Detect Results, Winding Path, Perfect Sensor, Nthr = (100% 50%)N, False Report Belief 0% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 4.10 Time to Detect Results, Straight Path, Perfect Sensor, Nthr = (100% 50%)N, False Report Belief 10% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 xiv 4.11 Time to Detect Results, Winding Path, 10% False Report Sensor, Nthr = (100% 50%)N, False Report Belief 10% . . . . . . . . . . . . . . . . . . . . . . . . . . 92 4.12 Time to Detect Results, Straight Path, 20% False Report Sensor, Nthr = (100% 50%)N, False Report Belief 20% . . . . . . . . . . . . . . . . . . . . . . . . . . 94 4.13 Time to Detect Results, Winding Path, 20% False Report Sensor, Nthr = (100% 50%)N, False Report Belief 20% . . . . . . . . . . . . . . . . . . . . . . . . . . 96 4.14 Tuned Values for Sensor Belief and Number of E ective Particles . . . . . . . . 103 4.15 Time to Detect and False Start Results, Straight Path, 10% False Report Sensor 105 4.16 Time to Detect and False Start Results, Winding Path, 10% False Report Sensor 106 4.17 Time to Detect and False Start Results, Straight Path, 20% False Report Sensor 107 4.18 Time to Detect and False Start Results, Winding Path, 20% False Report Sensor 108 4.19 Time to Actually Detect, Straight Path . . . . . . . . . . . . . . . . . . . . . . 110 4.20 Time to Actually Detect, Winding Path . . . . . . . . . . . . . . . . . . . . . . 112 5.1 Path Following Constants . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 5.2 Path to Node 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 5.3 Path to Node 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 5.4 Path to Node 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 5.5 Path to Node 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 5.6 Paths to Node 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 5.7 Path List . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 5.8 Final Path List . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 5.9 Time to Intercept and Run Time Results, without UAV measurements . . . . . 140 5.10 Time to Intercept and Run Time Results, with UAV measurements . . . . . . . 141 5.11 Comparison of ML and IUF routines using UAV measurements without regional sensor data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 5.12 Time to simulate one second . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 5.13 Comparison of ML and IUF routine using only large regional sensor data . . . . 142 xv Chapter 1 Introduction Modern statistical and engineering practice requires the ability to estimate and predict the behavior of nonlinear systems with accuracy, precision, and e ciency. This need, coupled with increased computational capabilities, has lead to a focused interest in nonlinear ltering over the past 30 years. Most problems of interest require a sequential estimate of a set of variables that may or may not change over time, called the state of the dynamic system. The state completely describes the dynamic system under investigation. Estimates of the state are made by analysis of measurements taken on that system. Measurements are assumed to be imperfect or noisy. One scenario of signi cance to which nonlinear ltering has been applied is that of a ground target vehicle in an urban environment. The problem of locating, tracking, and intercepting a ground target in an urban environment presents numerous challenges. Tradi- tional measurements such as range and range rate are not often readily available at regular intervals in such a setting. Information may come in the form of binary measurement reports from sources such as human observers, xed cameras, or intermittent satellite images. This nontraditional measurement format can provide valuable information, although it is coupled with increased uncertainty. Therefore, it must be processed in an e cient and e ective manner in order to account for its irregular frequency and non-di erentiable measurement model. Additionally, prior knowledge of terrain and road networks provide vital insight into the target?s state. Such information cannot be properly accounted for by estimation routines based on Gaussian assumptions. This dissertation is centered on the numerous bene ts the particle lter o ers in solving certain di cult nonlinear ltering problems that exist, namely in urban warfare. Chapter 2 1 describes the most common approaches to the nonlinear ltering problem and the di erences in their probabilistic representations of the state. The particle lter is compared to its more commonly used counterpart, the Kalman Filter, and its variations. The ease of adaptability of the particle lter to any motion model and the ability to represent knowledge of the state with a probability density function of arbitrary shape decided its application to this problem. Section 2.4 details previous work done in the eld of particle lters, path planning, and risk assessment, as it relates to this dissertation. Chapters 3 and 4 describe the bene ts of improving and tuning each step in the particle ltered estimation process of prediction, correction, and resampling. Chapter 3 compares two methods for modeling the prediction step: a simplistic model and a heuristic model based o of conventional tra c rules. Both dynamic models are compared in regard to the time to detect the target, their ability to track the target, and their computational expense. Chapter 4 introduces a sensor model that provides false reports. Throughout the dynamic model study, the e ect of additional variables are also investigated. The measurement model in question involves a human component, which presents the possibility of a false report. In the presence of a false report, the spatial resolution of the particle cloud is compromised based on bad information. False reports diminish tracking performance, but improvements may be made by adjusting built in parameters in the particle lter framework. The measurement update step updates particle weights in accordance with the amount of con dence the lter has in the measurements it receives. The amount the particle lter trusts a measurement may be tuned to reduce error when a sensor of known false alarm rate is in play. Section 4.1 details a study to determine the e ect of the amount of con dence the particle lter has in a sensor?s report. This con dence parameter is present in the likelihood function in the measurement update, or correction, step. The study involves Monte Carlo runs done with three sensor models and eleven possible con dence values for each sensor. 2 The frequency of resampling is also a particle lter design point. Resampling is done when the number of e ective particles, or particles with signi cant weight, falls below a de- sired threshold. When imperfect measurements are known to be possible, resampling after each measurement is bad practice as the spatial resolution is compromised. An investigation of the e ect of the number of e ective particles was done to tune the lter to an optimal value based on the sensor?s estimated performance. Section 4.2 explains how the frequency of resampling was treated as a control variable through tuning the number of e ective parti- cles. Following the determination of the most bene cial values of con dence and number of e ective particles for each sensor model, a Bayesian risk assessment was included as outlined in Section 4.3 as an additional means of dealing with possibly false reports. The discussion of the second phase of this dissertation begins in Chapter 5. An un- manned aerial vehicle (UAV) is introduced into the urban scenario with the mission of intercepting the ground target. A path planning routine was developed by the use of re- ceding horizon control. A novel path-selection metric is developed to provide for target interception in the particle lter framework. A new cost function is created to avoid high computational expense of previously used minimum entropy formulations. The path planner selects paths with high particle weight, taking anticipated travel time into account. In addi- tion, a minimum distance between a path and the location of the heaviest particle is desired. This method selects paths with maximum likelihood, while reducing entropy, with far fewer calculations than previous work. This will prove most e ective in smaller platforms where computational power is limited, as in small UAV?s to be utilized in urban warfare. The result of this work is a robust simulation with real world applications and bene ts to the modern war ghter. 3 Chapter 2 The Nonlinear Dynamic System Filter Problem The objective of nonlinear ltering is to recursively estimate the state of a dynamic system based on measurements. The estimation of the state is typically broken into two steps: prediction and correction. The prediction step models the propagation of the state through time between measurements. The equations of motion of the state variables are propagated through time until a measurement is taken. Process noise is included in the state equations to account for inaccurate modeling and unforseen disturbances. This noise causes the accuracy of the state estimate to diverge from the true state until a measurement is taken. In the ground-target scenario, this divergence can be attributed to uncertainty in a target?s position, velocity, or destination. The correction step incorporates measurement data into the state estimate. Each measurement is taken into account in accordance with its con dence to improve the state estimate. These models can be represented probabilistically and therefore lend themselves to the Bayesian framework. In this rigorous general approach, a probability density function (pdf) is used to store all knowledge and belief about the true state. During the prediction step, the state pdf is generally translated, deformed, and broadened, whereas following the corrector step, it is typically tightened. The state pdf is updated with the new measurement by applying Bayes? theorem. Suppose the target state vector is xk 2Rn, where the k is the time index, R is the set of real numbers, and n is the dimension of the state vector. The state evolves according to a discrete-time stochastic model in a rst-order Markov process. xk = fk 1 (xk 1;vk 1) (2.1) 4 In Eq. (2.1), the function fk 1 describes the behavior of the state at xk 1, and vk 1 is the process noise. Because the state evolves according to a Markov process, the future state is only dependent upon the current state. zk = hk (xk;wk) (2.2) Measurements are related to the state by the measurement vector of dimension m. In Eq. (2.2), the function hk is a known and possibly nonlinear function of the state and it includes some measurement noise, wk. Both the measurement noise and process noise are assumed to be white, independent, and to have known probability density functions. The sequence of all measurements, Zk fzi;i = 1;:::;kg, is used to obtain ltered estimates of the state. For implementation in the Bayesian framework, a posterior pdf must be constructed to quantify the belief in the state xk given Zk. This pdf, p(xk jZk), is calculated recursively using the prediction correction process. It is assumed that p(x0 jZ0) is known where z0 is the set of no measurements and therefore independent of all noise. To begin recursive estimation, it is assumed that p(xk 1 jZk 1) is available. The pre- diction step is carried out by using the system equation (2.1) to create the predicted or prior density of the state at time k via the Chapman-Kolmogorov equation: p(xk jZk 1) = Z p(xk jxk 1)p(xk 1 jZk 1)dxk 1 (2.3) The probabilistic model of the state evolution, or the transitional density, p(xk jxk 1), is de ned by the system equation (2.1) and the known statistical parameters of vk 1. When a measurement, zk, becomes available, the correction stage begins. By use of measurement zk, the prior density p(xk jZk 1) is modi ed and the posterior density of the current state 5 p(xk jZk) is obtained. p(xk jZk) = p(xk jzk;Zk 1) = p(zk jxk;Zk 1)p(xk jZk 1)p(z k jZk 1) (2.4) = p(zk jxk)p(xk jZk 1)p(z k jZk 1) p(zk jZk 1) = Z p(zk jxk)p(xk jZk 1)dxk (2.5) From this posterior density, an optimal state estimate may be computed. Such forms for the estimate include the the minimum mean-square error (2.6) and the maximum a posteriori estimate (2.7), or the conditional mean of xk and the maximum of p(xk jZk), respectively. ^xMMSEkjk Efxk jZkg= Z xkp(xk jZk)dxk (2.6) ^xMAPkjk arg maxx k p(xk jZk) (2.7) Recursive propagation of the posterior density through Eqs. (2.3) and (2.5) provide the basis for the optimal Bayesian solution, however it is only conceptual in that it cannot be generally determined analytically. The storage and representation of the entire pdf requires a nite dimensional representation of an in nite dimensional function. Consequently, approx- imations or suboptimal Bayesian algorithms are conventionally used in practice. Namely, the Kalman lter, the extended Kalman lter, the Unscented Kalman lter, and the particle lter are methods of note in tracking applications [1]. 2.1 The Kalman Filter In 1960, R. E. Kalman published a paper that provided a rigorous mathematical ap- proach for sequentially processing observations of a linear dynamic system [6]. This approach was conveniently introduced at a time when computational power and technology were on 6 the rise, allowing for its successful implementation in numerous applications and consequen- tial popularity. The set of recursive equations outlined in the 1960 paper are known as the Kalman lter [2]. The Kalman lter requires certain restrictive assumptions. It is assumed that the pos- terior density is Gaussian at every time step. This provides for the posterior density to be characterized by its mean and covariance, but it also requires that the system equation xk be a linear function of xk 1 and vk 1, the measurement equation zk is a linear function of xk and wk, and the noise terms vk 1 and wk are drawn from Gaussian densities of known mean and variance. If these restrictions hold, the Kalman lter yields the optimal solution [1]. In the Kalman lter framework, Eqs. (2.1) and (2.2) become: xk = Fk 1xk 1 + vk 1 (2.8) zk = Hkxk + wk (2.9) where Fk 1 (n x n) and Hk (m x n) are known. Also, the sequences vk 1 and wk are mutually independent zero-mean white Gaussian noise with covariances Qk and Rk, respectively. Due to the assumption that all posterior pdf0s are Gaussian, the pdf update equations become: p(xk 1 jZk 1) = N xk 1; ^xk 1jk 1;Pk 1jk 1 (2.10) p(xk jZk 1) = N xk; ^xkjk 1;Pkjk 1 (2.11) p(xk jZk) = N xk; ^xkjk;Pkjk (2.12) whereN(x; ;P) is a Gaussian density with argument x, mean , and covariance P, de ned by: N(x; ;P) j2 Pj 1=2 exp 12 (x )T P 1 (x ) (2.13) 7 The Kalman lter measurement update equations for the mean and covariance are given by the equations below. ^xkjk 1 = Fk 1^xk 1jk 1 (2.14) Pkjk 1 = Qk 1 + Fk 1Pk 1jk 1FTk 1 (2.15) Kk = Pkjk 1HTk HkPkjk 1HTk + Rk 1 (2.16) ^xkjk = ^xkjk 1 + Kk zk Hk^xkjk 1 (2.17) Pkjk = [I KkHk] Pkjk 1 (2.18) The Kalman Filter provides optimal solutions for a certain class of problems, but it was expanded upon in order to be applied to nonlinear systems and/or measurement models. This lead to the development of the Extended Kalman Filter and the Unscented Kalman Filter [1]. 2.1.1 The Extended Kalman Filter In practical situations, the strict restrictions on the Kalman Filter are not satis ed. The Extended Kalman lter (EKF) operates on the basic premise that the true state is su ciently close to the estimated state. This provides for the use of a linearized rst order Taylor series expansion to represent the error dynamics. The posterior pdf p(xk jZk) is assumed Gaussian and therefore Equations (2.10) through (2.12) are valid. Although the EKF is not considered optimal like the Kalman Filter, it is a well known technique for nonlinear system estimation [2]. In the Extended Kalman Filter framework, Equations (2.1) and (2.2) become: xk = fk 1 (xk 1) + vk 1 (2.19) zk = hk (xk 1) + wk (2.20) 8 The random noise sequences vk 1 and wk are assumed zero mean, mutually indepen- dent, white Gaussian with covariances Qk 1 and Rk, respectively. The following recursive equations are used to compute the mean and covariance of the state: ^xkjk 1 = fk 1 ^xk 1jk 1 (2.21) Pkjk 1 = Qk 1 + ^Fk 1Pk 1jk 1^FTk 1 (2.22) Kk = Pkjk 1 ^HTk h^ HkPkjk 1 ^HTk + Rk i 1 (2.23) ^xkjk = ^xkjk 1 + Kk zk hk ^xkjk 1 (2.24) Pkjk = Pkjk 1 Kk ^ HkPkjk 1 ^HTk + Rk KTk (2.25) The local linearizations of fk 1 and hk are ^Fk 1 and ^Hk. They are the Jacobian of the corresponding nonlinear equations evaluated at ^xk 1jk 1 and ^xkjk 1, respectively. The EKF is known as an analytic approximation because the aforementioned Jacobians must be derived analytically. Clearly, if the nonlinear functions fk 1 and hk are discontinuous or nondi er- entiable, the EKF may not be applied. In addition, the EKF is restricted by its assumption that the posterior pdf p(xk jZk) is Gaussian. Models with severe nonlinearities will result in increased errors due to this assumption. As a result, discrete sampling approaches, such as the Unscented Kalman Filter and the Particle Filter, were developed. These methods are more suited for urban tracking problems. 2.2 The Particle Filter The Particle Filter is a suboptimal lter that performs sequential Monte Carlo estima- tion utilizing independent random samples of probability densities. These samples, called particles, are point mass representations of the posterior pdf of the state. Like the Kalman lter, the Particle lter (PF) gained heightened popularity with an increase of computational capabilities. Although the framework for sequential Monte Carlo estimation originated in statistics in the 1950?s, the method was also limited by the use of sequential importance 9 sampling, which degenerates over time. The introduction of the resampling step in the early 1990?s [3], coupled with increased computational capabilities, lead to the drastic increase in PF implementation. Monte Carlo methods for the numerical evaluation of multidimensional integrals form the basis for the PF approach. Suppose the solution to the following multidimensional integral is required, where x2Rn: I = Z g(x)dx (2.26) Monte Carlo (MC) integration methods factorize g(x) = f(x) (x) such that (x) may be interpreted as the posterior pdf in the Bayesian framework. It is assumed possible to draw N independent samplesfxi;i = 1;:::;Ng, where N 1, distributed according to (x). This yields the unbiased MC estimate of I, IN, that converges to I with error of the order O(N 1=2). I IN = 1N NX i=1 f(xi) (2.27) The rate of convergence is independent of the dimension of the integrand because the samples, xi, are chosen from regions of high importance relative to the state space. It is, however, often di cult to sample from the posterior pdf as it is nonstandard, multivariate, and not fully known. Consequently, the MC integration method of importance sampling is used sequentially to draw samples that represent the posterior pdf. The importance- sampling method samples and weighs the points from a known importance density function to simulate the samples from an unknown distribution. This method represents the posterior pdf with a set of random samples and associated weights, then estimates the state based on the samples and weights. When importance sampling is applied to the nonlinear ltering problem described in Section 2, the result is the Sequential Importance Sampling (SIS) algorithm, a sequential MC lter also known as the Particle Filter. 10 An initial estimate of the state pdf, p(x0), is assumed known, may be arbitrarily chosen, and is described by the set of particles xi0 with associated weights wi0, such that PNi=1wik = 1 at any time k. The PF follows the same prediction-correction pattern as the aforementioned lters. The state pdf at time k 1 is approximated by a sum of delta functions: p(xk 1 jZk 1) NX i=1 wik (xk 1 xik 1) (2.28) The prediction step consists of propagating the particles forward in time according to the dynamic model, xk = fk 1(xk 1;vk 1). This yields the predicted state pdf, p(xk jZk 1). When a measurement, zk, is received, the state pdf is updated in the correction step. p(xk jzk) = p(zk jxk)p(xk jZk 1)p(z k jZk 1) (2.29) The measurement likelihood function, p(zk jxk), is used to update the particles weights according to Bayes? Rule. wik = w i k 1p(zk jx i k)P N j=1 w j k 1p z k xj k (2.30) It has been shown that the variance of the importance weights can only increase over time in the presence of measurements [4]. This degeneracy phenomenon results in one particle holding all of the weight after a certain number of recursive steps. It is desirable to keep track of the number of e ective particles that hold a non-negligible weight, Ne < N. This is done through the following calculation: Ne = 1PN i=1(w i k)2 (2.31) The process of resampling is used when Ne falls below the desired threshold, Nthr. Although resampling is computationally expensive because it is not parrallelizable, it is necessary to avoid particle degeneracy. Resampling eliminates particles with relatively low weight and creates duplicates of particles with relatively high weight. The new set of particles is obtained 11 by resampling with replacement N times from the approximation of p(xk jzk) given by the set of particle states and weights,fxik;wig. The resulting particle set is an independent and identically distributed sample from the posterior pdf. There are multiple resampling meth- ods, but one of the most computationally e cient,O(N), and simple methods is systematic resampling. The pseudocode for systematic reampling is given below [1], where xik is the state of the ith particle, wik is the weight of the ith particle, and CSW is the cumulative sum of the particle weights. Initialize the CSW : c1 = w1k for i = 2 : N do Construct CSW : ci = ci 1 +wik end for Start at the bottom of the CSW : j = 1 Draw a starting point : u1 U[0;N 1] for i = 1 : N do Move along the CSW : ui = u1 + (i 1)=N while ui >ci do j = j + 1 end while parent(i) = j end for Assign sample : xk = xk(:;parent) Assign weight : wk = 1=N Figure 2.1 pictorially describes the process for determining which particles should be eliminated and which should be duplicated. The cumulative sum of particle weights (CSW) is compared to the random variable ui;i = 1;:::;N that is uniformly distributed on the interval [0;1]. If a particle?s weight does not signi cantly increase the cumulative sum when compared to the next element in the uniform distribution, that particle will be eliminated. 12 Conversely, if a particle?s weight signi cantly increases the cumulative sum when compared to the next element in the uniform distribution, that particle will be duplicated at least once. Given su cient computing power and trustworthy sensor input, resampling should be Figure 2.1: The Resampling Process done following each measurement to maintain su cient particle resolution. Figure 2.2 [1] illustrates the introduction of a measurement and the e ect of resampling. More particles are introduced in areas of higher probability and particles are removed from areas of low probability. Figure 2.2: E ect of Resampling It is evident that the PF framework does not make assumptions regarding the shape of the posterior pdf. Also, it puts no restriction on the forms of the dynamic or measurement 13 models. Therefore, it is utilized in this work to estimate the location of a moving ground vehicle. 2.3 Entropy of a Probability Density Function Entropy is often used as a measure of uncertainty of an estimate in the eld of target tracking, [24, 26, 28]. It is a positive scalar quantity that represents the amount of un- certainty in an estimate. Entropy reduction has been used in cost functions in problems regarding sensor assignment to determine the best con guration to reduce uncertainty in target location. For a general ltering density, p(xkjZk), the entropy is given by: H (p (xkjZk)) = Z X log2 (p(xkjZk))p(xkjZk)dxk (2.32) For an n-dimensional Gaussian variable with covariance matrix P, as in the Kalman lter framework, the entropy H is given by: H = log p (2 e)njPj (2.33) When put into the particle lter framework, the resulting expression for entropy is: H (p (xkjZk)) NX i=1 wik log2(wik) (2.34) This result is incorrect as it gives no mention to the location and local density of the particle distribution. For example, the two sample distributions illustrated by Fig. 2.3 each contain four particles of equal weight (0.25) that are used to represent the location of a target [5]. 14 (a) (b) Figure 2.3: (a) Low density pdf and (b) High Density pdf The distribution in Fig. 2.3(a) is more spread out than that in (b), corresponding to more uncertainty in target location given the distribution. However, both distributions would have an entropy value of 2 if Eq. (2.34) is used. The development of an approximate expression for the entropy of a non-Gaussian pdf as is often present in particle lter implementations has been the subject of recent research. Driessen et al. [5] used a 2-D bimodal distribution in a target tracking application to validate the following entropy approximation: H (p(xkjZk)) = Z X log2 (p(zkjxk)p(xkjZk 1))p(xkjZk)dxk+ log2 (p(zkjZk 1)) H (p(xkjZk)) log NX i=1 p zkjxik wi ! NX i=1 log p zkjxik NX j=1 p xikjxik 1 wjk 1 !! wik 15 Skoglar et. al [26] implemented and tested the following expression on a 1-D Gaussian sum density: H (p(xkjZk)) NX i=1 wik logp zkjxik NX i=1 wik log NX i=1 wik 1p xjkjxik 1 + log NX i=1 wik 1p zikjxik 1jk In addition, Ryan [28] derived a piecewise linear approximation of a pdf represented by a particle set by using particle locations as vertices of triangular elements, and the height of the elements being p(xikjZk). Figure 2.4 illustrates the evolution of the particle cloud after a measurement update and the piecewise linear approximation of the pdf implemented by Ryan [28] to calculate the entropy of the distribution. Figure 2.4: Piecewise linear approximation of a PDF As there is yet to be a standard expression for the entropy of a distribution represented by particles, an indirect method for reducing entropy is desired. Ho mann and Tomlin [24] 16 present an information utility function, I(zkjxk), that when maximized, will reduce the ex- pected entropy of the target state, H (p(xkjZk)). H (p(xkjZk)) = H (xk) I (zk; xk) where H (xk) = Z p(xk) log2p(xk)dx H (p(xkjZk)) = Z p(xk;zk) log2p(xkjzk)dxdz I (zk; xk) = Z p(xk;zk) log2 p(xk;zk)p(x k)p(zk) Additionally, Ho mann and Tomlin express I (zk; xk) as: I (zk; xk) = H (zk) H (p(zkjxk)) which can be approximated by: H (zk) Z ( NX i=1 wi kp z kjxik log 2 " NX i=1 wi kp z kjxik #) dz H (zkjxk) Z NX i=1 wi kp z kjxik log 2p z kjxik dz The Information Utility Function is based on the idea that minimizing the posterior un- certainty is equivalent to maximizing the di erence between the uncertainty that any par- ticular observation will be made H (zk), and the uncertainty of the measurement model, H (p(zkjxk)). The terms above are readily available in a particle lter and may be used to select sensor placement to indirectly reduce the entropy of a particle distribution. 17 2.4 Applications to Target Localization and Tracking The most computationally e cient algorithms used for ground target tracking are based on Kalman ltering [5, 6, 7]. However, the aforementioned strict assumptions in Kalman ltering do not hold in urban target tracking, as the system is nonlinear, non-Gaussian, and possibly multi-modal. A particle lter requires no assumption about the distribution of the uncertainties and linearity of the system dynamics, consequently a non-Gaussian posterior distribution is admissible [1, 3, 4, 8]. The urban terrain is characterized by a road network, tra c signs and signals, speed limits, crossings, and tra c participants such as vehicles, bicycles, and pedestrians. Solutions that do not utilize these features may lead to suboptimal performance. Exploitation of prior knowledge of the terrain attributes, for instance, road maps and tra c information, are therefore highly desirable to signi cantly improve the target tracking performance [9]. The earliest usage of particle lters for target tracking with road network information can be found in [10], where Arulampalam et al. presented an algorithm termed Variable Structure Multiple Model Particle Filter (VS-MMPF) with ground moving target indicator (GMTI) radar using non-standard information such as road maps and terrain-related visibil- ity conditions. Subsequently, Agate and Sullivan [11] demonstrated the feasibility of jointly tracking and identifying targets using a particle lter. Yang et al. [12] explored several ways to include terrain database and road maps to assist the tracking of ground targets and dis- cussed motion models for brake to stop, road intersections, and target interactions. Further, Ulmke and Koch [13] developed a GMTI based target tracking approach using the discrete Gauss-Markov target dynamics and road-map information. Ekman and Sviestins [14] intro- duced the Multiple Likelihood Model Particle Filter (MLM-PF) for target tracking in the presence of hard constraints such as speed or acceleration demonstrating excellent tracking performance. Concurrently, Kamrani and Ayani [15] presented a method for path planning of an unmanned aerial vehicle with the task of tracking a moving target on a road network. 18 Orguner et al. [16] considered the problem of tracking targets, which can move both on- road and o -road, while the results suggested the preference of Interacting Multiple Model Particle Filter over Bootstrap Multiple Model Particle Filter. Kreucher et al. [23] utilized a particle lter and the expected information gain as the metric to assign a sensor to track multiple targets. Each second, a 100 m x 100 m section of the simulation space was inspected and a binary response of target presence was given. With the intention of reducing computational time, Hong et al. [17] presented a Multi- rate Interacting Multiple Model Particle Filter that reduced by half the computational cost in comparison to the Multiple Model Particle Filter for terrain based ground target track- ing. Thereafter, Kravaritis and Mulgrew [18] introduced the Variable-mass Particle Filter (VMPF) for terrain-aided tracking problem by allocating particles with variable masses to the modes as against VS-MMPF. The simulation results in [18] showed the improved e - ciency of the VMPF. Skoglar et al. [19] proposed a Rao-Blackwellized Particle Filter to reduce the dimension of a state space in road target tracking application, and showed the use of prior information about the probability of detection can be used to improve the estimation performance. The goals of the search problem are well framed in the context of information theo- retic costs. Information theoretic cost metrics have been used to manage sensors [20], and have led to algorithms to control sensor networks for information gathering over an area by parameterizing the motion of collectives of vehicles [21]. The optimal probing control law to minimize Shannon entropy for the dual control problem has been shown to be the input that maximizes mutual information [22]. The expected alpha-divergence of a particle lter distribution was used for sensor management, and specialized to select modes for binary sen- sors, though scalability in sensor network size was not addressed, and the Shannon entropy was only found in the limit of the presented equations [23]. An earlier version of entropy approximation techniques was presented in [24]. The continuation of that work, presented in [25], develops an entropy-based metric to maximize the mutual information in a mobile 19 sensor network. Signi cant work has been done in the eld of sensor assignment. Spletzer and Taylor select the di erence in the expected value of the pdf and each particle?s location as the metric for path selection [31]. Wang et al. study a linearized system with a Kalman Filter in the Bayesian framework to optimize when to accept an estimate and the risk associated with it based on Renyi information divergence. [32]. O?Reilly applied Bayesian likelihood ratio tests to the detection of faults within a sensor by observing potential changes in the variance of the estimate [33]. 20 Chapter 3 The E ect of the Particle Motion Model A ground vehicle is to be located and tracked in an urban environment. This chapter will investigate the e ect of a high delity particle motion model on localization and tracking performance given unconventional measurements. Two motion models will be compared to explore the bene ts and development of a sophisticated dynamic model. The performance of each motion model will be studied in three ways: computational expense, the amount of time to detect the target, and the tracking capabilities post-detection. The particle lter is not restricted in the form of the dynamic models or probability distributions. This chapter will investigate the multi-resolution problem of maintaining su cient spatial resolution while distributing particle weight according to the likelihood of target presence. A thorough description of the simulation environment will be discussed, as well as the necessary assumptions that were made. The motion model utilized to simulate the target ground vehicle will be introduced, followed by the measurement model and the two particle motion models. Section 3.5 presents comparison results for the two motion models. Reduced performance by the high- delity model in the amount of time it takes to detect the target motivated further investigation into the role of particle spatial resolution by means of particle redistribution and allowing for the belief of the presence of false measurement reports. 3.1 Problem Description Suppose there is a mobile ground vehicle of known description but unknown position, velocity, and destination that is to be found, tracked, and intercepted by an unmanned aerial vehicle. The vehicle is known to be in an urban environment, and full knowledge of that environment (roads, obstacles, intersection constraints, and speed limits) is available. Urban 21 warfare often relies on information from passersby, agents in the eld, and tra c cameras. The inclusion of such non-traditional and non-di erentiable measurements in addition to the shape of a complex road network necessitates the use of a particle lter to estimate the location of the target vehicle. In this problem, measurements are received from eld operatives in the form of binary responses of target presence within a sector of the city. There are numerous issues of interest within this problem. The e ect of a sophisticated particle motion model on target localization time and tracking capabilities in the presence of unconventional measurements will be studied in this chapter. Figure 3.1: Map of Urban Environment The urban environment utilized in this work is illustrated by Fig. 3.1, where obstacles are orange and roads are white. Complete knowledge of the map is assumed, including the location of obstacles and roads. Several assumptions are made throughout the simulation. Roads have known speed limits related to their width (as consistent with practice), are assumed to have constant heading, and permit two way travel. Large avenues have a speed limit of 35 mph and smaller roads have a speed limit of 25 mph. The target is restricted to 10 mph of the 22 speci ed road speed limit to simulate a realistic tra c scenario. Although the same map is used throughout this work, the particle lter framework is suitable for any road network [29]. 3.1.1 Target Motion Model It is known that the target vehicle is located within the urban environment, with no bias given to any particular road initially. Therefore, the initial particle distribution, p(x0), consists of N particles of equal weight scattered throughout the roadways. The x and y locations of each particle is drawn from the distribution 1000U[0;1], under the restriction that the point [x;y] lies within a road. A uniform distribution was used to create the initial particle distribution because there is equal probability that the target lies in any location throughout the map. for i = 1 : N do rand U[0;1] xik = 1000rand rand U[0;1] yik = 1000rand point = [xik;yik] while point2Obstacle do rand U[0;1] xik = 1000rand rand U[0;1] yik = 1000rand point = [xik;yik] end while wik = 1=N end for 23 It is assumed that the target is a mobile ground vehicle that is modeled as a simple car illustrated by Fig. 3.2 [30], with position [x,y], wheelbase L, turn angle , and heading measured clockwise from the positive y-axis. Figure 3.2: Geometry of a Simple Car min = Ltan max (3.1) The wheelbase, L, is set to 2.8 m and a minimum turning radius of 11.5 m is utilized throughout the simulation, as consistent with a standard sedan. Equation (3.1) is used to compute a maximum turning angle, max, of 13:7 . The target is required to remain on the right side of the road and stay within 10 mph of the road speci c speed limit. The speed at each time step is updated with noise to account for unknown accelerations. The value for the noise, k, is a random draw from the distribution U[0;1]. Vk+1 = Vk + Vk( k 0:5)2 (3.2) The resulting range of possible values for Vk+1 is the set [0:75Vk;1:25Vk]. If the speed escapes the bounds of 10 mph of the speed limit, it is adjusted by the following control sequence: Vmin = Vroad 10 Vmax = Vroad + 10 24 while Vk+1 Vmax do Vk+1 = Vk Vkj( k 0:5)=2j end while end while In addition, the target may not reverse and must remain within the simulated map, and therefore is required to make a U-turn at the map?s edge. Although the target?s path is unknown to the estimation routine, two candidate paths were chosen for simulation, one predominantly straight and one winding, and are illustrated by Fig. 3.3. (a) (b) Figure 3.3: (a) Straight Path and (b) Winding Path 3.1.2 Measurement Model It is assumed that a soft measurement source (a measurement from a human) is able to provide information on each 200 m 200 m sector depicted by the grid squares in Fig. 3.4, and that there is no cost associated with utilizing any sensor. Every three seconds, a measurement is taken at the sector containing the highest particle weight sum. The highlighted square region in Fig. 3.5 represents the sector being measured. The structure of the sensor grid was arbitrarily chosen. The lter framework is also applicable to overlapping 25 sensor regions, sparse sensor coverage, non-stationary sensors, and sensor regions of irregular shape. In addition, the measurement report rate of one measurement per 3 seconds was chosen to observe the e ect of receiving reports at a reduced frequency. Any frequency could be used in this framework. Finally, only one sensor region was polled per measurement. However, if in practice more regions may be polled at once, the particle lter framework presented herein is suitable for such an application. Figure 3.4: Determine Region of Highest Particle Weight The polled sensor returns a binary response, where 0 means the target is not present and 1 means the target is present somewhere in the associated region. No speci c location within the region is given. This information is fused into the measurement likelihood function to update particle weights. A sensor that always reported the truth was initially used. However, the assumption of a perfect sensor is unrealistic, especially as the sensor is modeled after a human operator. Because of this, false positive and false negative rates of 10 % were later introduced into the sensor model. Methods used to account for an imperfect soft binary sensor will be discussed in Chapter 4. 26 Figure 3.5: Measure Region of Highest Particle Weight 3.2 Dispersion Motion Model The e ect of the particle motion model was the rst major investigation in this work. The initial particle dynamic model was that of a constant-velocity point-mass unicycle, as illustrated by Figure 3.6 where r = 0. This model caused the particles to randomly disperse through the road network until they reached an obstacle, at which point their heading was rotated 180 deg. The dispersion model was chosen for comparison because of its ease of implementation and its tendency to ll the roadway. A similar model was used by Kreucher et al. in [23] to track a target with the same dynamic model on a random walk, with no obstacles obstructing its motion. The following equations describe the particle motion in the dispersion model: xik+1 = xik +Vik sin ikdt (3.3) yik+1 = yik +Vik cos ikdt (3.4) Vik+1 = Vik (3.5) ik+1 = ik + ( k 0:5) 10 (3.6) 27 where the ith particle?s position at time k is denoted [xik;yik], its heading ik is measured clockwise from the positive y-axis, the time step dt is 1 second, and the noise k is a random draw from the distribution U(0;1), thus limiting the change in heading per second to the range [ 9 ;9 ]. Figure 3.6: Geometry of Unicycle 3.3 Tra c Motion Model In addition to the dispersion model, an alternative rst-order Markov model was devel- oped to exploit the knowledge of vehicle tendencies in an urban environment. In an urban setting, state variables may not be solely time dependent, but also location dependent. The tra c model accounts for road conditions such as in-lane travel, speed limits, two-way tra c, changing lanes to execute a future turn, and pauses at intersections. A complete characteri- zation of the road network is assumed available. This requires categorizing regions into four categories: road, intersection, obstacle, or o -map. Each road has two possible headings that re ect two-way tra c. The required heading is based on a particle?s position on the road, as to mimic two way tra c ow and driving on the right side of the road. Figure 3.7 illustrates the required headings of two roadways and their intersection at the grey region. Each particle stores a variable in its state vector that is related to the required heading of the road it is on, or the required heading of the road it is turning on to. This variable, q, becomes important during a turning maneuver. 28 Figure 3.7: Intersection of Two Roads Because state variables in an urban setting are both time and location dependent, each particle stores a variable in its state that corresponds to the motion model it is to use. The value of this variable mik is based on whether the particle is on a road, going straight in an intersection, turning, left, turning right, or making a U-turn. The motion models for each of these modes are described below. While on the rth road with speed limit Vr and mik = 2, the ith particle is propagated forward in time using the following discrete-time dynamic model: xik+1 = xik +Vik sin ikdt (3.7) yik+1 = yik +Vik cos ikdt (3.8) Vik+1 = Vik + V i k( i k 0:5) 2 (3.9) Subject to Vr 10mph Vik+1 Vr + 10mph (3.10) ik+1 = ik (3.11) where xik and yik de ne the position relative to the origin, which is located in the bottom left corner of the map in Fig. 3.1, ik is the heading angle measured clockwise from the positive y-axis, the time step dt is 1 second, and k is a random draw from the distribution N(0;1). Because the road network used in this work is modeled after a grid-like street pattern, the roadways are straight and a constant heading is desirable along them. However, a directed 29 graph could be used to determine a road?s heading in order to apply this model to any road network [29]. As mentioned in Section 3.1.1, the target is required to make a U-turn at the map?s edge. All particles are subject to the same restriction. If when propagated forward to the next time step, a particle would be located outside of the map?s boundaries, a U-turn is performed. During this maneuver, it is necessary to know which road the particle is on so that the proper heading and position may be obtained following the U-turn. Figure 3.8: U-turn course adjustment Figure 3.8 illustrates how a particle?s position and heading are adjusted when it is to cross the map?s boundary, drawn in red. In this example, the center of the road is the dashed line and the particle is initially on the lower half of the road, heading to the right. During the U-turn, it moves to the other side of the road. The distance between the particle and the edge of map is d1, and the distance between the particle and the center of the road is d2. The particle?s proximity to the center of the road is preserved following a U-turn as that distance is used to determine what turns a particle is capable of performing. The heading variable is changed by radians and the x and y position variables are adjusted based on heading and the values of d1 and d2, respectively. The equations for the U-turn motion model depend on the particle?s initial heading. The equations below are used during the U-turn motion model for the example in Fig. 3.8 where the particle?s initial 30 heading is =2. The three additional sets of equations for the U-turn motion model are derived in a similar manner. xik+1 = xik (3.12) yik+1 = yik + 2d2 (3.13) Vik+1 = Vik + V i k( i k 0:5) 2 (3.14) Subject to Vr 10mph Vik+1 Vr + 10mph (3.15) ik+1 = ik + (3.16) If a particle will be located in an intersection at the next time step, it must make a decision on which way to go based on what is permitted by the road network. At the start of the simulation, a list of possible turns from every direction at each intersection is generated based on the map. This turn list provides for obstacle avoidance. When a particle enters an intersection, its location is used to determine which intersection it is in and its heading is used to determine the turning possibilities. The list of possible turn choices is compiled and a random selection is made to indicate the chosen turn direction, either left, right, or straight. Next, the particle?s proximity to the obstacles on its right and left, velocity, span-wise location on its current road, and the width of all possible future roads are taken into account to determine what adjustments must be made to complete the chosen turn. Figure 3.9 illustrates the various parameters that are taken into account to prepare for a turn, where the particle?s initial position is in the bottom right with a heading of = 0 in the positive y direction. First, the width of the possible future road, denoted as r2 in Fig. 3.9, is calculated. In the map used in this simulation, the width of the road is determined by the size of the intersection and is assumed constant along the length of the road. If either rright or rleft are greater than r2, a lane change would be necessary to successfully make the turn in their respective directions. In the example illustrated by Fig. 3.10, the ith particle at position 31 Figure 3.9: Geometric characterization of an intersection [xik;yik] has chosen to turn right, but its distance from the nearest obstacle to the right, rright, is larger than half of the width of its future road, r2. If it were to turn in its current position, it would be on the wrong side of the road when it reached the new road. Therefore, it must merge to the right. The particle will merge to its new distance from the obstacle, rnew, which is equal to min + 1. xik+1 = x1;obs ( min + 1) (3.17) The lane change maneuver depends on the chosen turn direction and the particle?s heading. Similar update equations are derived for the remaining con gurations. Lane changes are permitted when approaching an intersection to provide for improved tracking. This was motivated by the resampling step in the particle- lter framework. As a result of resampling, particles of signi cant weight are duplicated into multiple particles of smaller weight. Because the noise in the velocity term will only separate duplicated particles 32 Figure 3.10: Example of a Required Lane Change along the length of a road and not along the width, lane changes allow the particle cloud to more fully describe the possible motion of a target vehicle entering an intersection, without violating the simulated vehicle?s physical constraints, namely min. The appropriate lane is determined by the particle?s velocity and chosen turn direction. Once a turn direction is randomly chosen and a lane adjustment is made, if necessary, the particle?s turning parameters are computed. While turning, a particle?s heading will change a total of =2 and the particle?s path during a turn will follow an arc of constant radius, i, which is the distance between the ith particle an the nearest obstacle in its chosen turning direction. The two choices for i are rright and rleft, as illustrated by Fig. 3.9. Because the dynamic system is a rst order Markov process, the future state depends only on the current state. Therefore, only the state at time k is stored. In addition, the time step used in the simulation is 1 second. This presented a challenge during turning maneuvers, as the amount of time it takes to complete a turn is not always an integer or 33 constant. This required the use of the aforementioned variable q to determine when the appropriate amount of turning had occurred. The value of a particle?s q2[1;2;3;4] and it is related to a particle?s heading by Eq. (3.18) and is illustrated by Fig. 3.11. ik = (qik 1) 2 (3.18) Figure 3.11: De nition of q Variable If a particle chooses to turn right, its mode variable mik becomes 1. Similarly, if a particle chooses to turn left, mik is set to -1. When a particle makes a decision to turn, qik is updated by use of the mode variable mik, as indicated in Eq. (3.20), and is always restricted to the set [1;2;3;4]. The value of qik is used to update a particle?s heading during a turn. While on a road, qik remains constant as the heading is constant. qik = qik +mik (3.19) qik = mod(qik;4) (3.20) Given the selected turn direction and turning radius, the number of time steps to complete a turn, Ns, is then computed based on the arc angle to be traversed per time step, . Because the number of time steps to complete a turn depends on the velocity, no noise is introduced 34 into the velocity during a turn, thus keeping it constant. tan ik = L i k (3.21) i = mV i k L tan i k (3.22) Nis = =2 i + 1 (3.23) Figure 3.12 illustrates the parameter that is utilized during a sample turn where Ns = 5, m = 1, and the value of q changes from 1 to 2. The position of the ith particle along the arc during a turn at each time step k is illustrated by the black dot. Figure 3.12: Geometry During a Sample Turn The angle i is always less than or equal to i as it is the remainder of the arc to be traversed at the nal time step of a turn. This simulates steering to stay in the same lane. i = 2 m i(Ns 1) (3.24) 35 Figure 3.13: Geometric Turning Parameters Once the above parameters are determined, the time update parameters illustrated by Fig. 3.13 are determined. The time dependant variable k is the angle between the heading at which the particle entered the intersection, k0, and the next way point along the arc, [xk+1;yk+1]. The constant h is the distance between [xk;yk] and [xk+1;yk+1]. These param- eters are computed by the following set of equations, where k0 is obtained by determining the previous value of q: k0 = [(q m) 1] 2 (3.25) k = k0 + [(k 1) + 1=2] (3.26) h = m sin( )=cos( 2 ) (3.27) Finally, the state variables for the ith particle are updated by the algorithm below for the sample turn in Fig. 3.12. The variable k is stored in each particle?s state vector in order to keep track of how far into a turn theith particle is. if ki 90%N Ne , falls below a desired threshold,Nthr, the particle cloud is resampled. Equation 2.31 is repeated below for convenience. Ne = 1PN i=1(w i k)2 (3.41) If the required number of e ective particles is set less than the total number of particles (N), the distribution may not be resampled after each measurement. This could provide for reduced detection time by preserving some particles This could be important in the presence of a sensor with a false report rate. By not resampling after each measurement, the spatial resolution of the particle cloud is not immediately lost after a measurement, reducing the e ect of a false report. For example, consider the initial particle distribution is illustrated by Fig. 3.30 (a) where the required number of e ective particles is set to 90%N and a non- detection is reported in the region where 400 x 600 and 200 y 400. Figure 3.30 (b) illustrates that the particles in the region of non-detection are not automatically eliminated because the distribution was not resampled. An investigation was done into the e ect of reducing the required number of e ective particles on the spatial resolution of the particle distribution. Two thresholds for the required number of e ective particles were used: 90%N and 80% N. The sensor always reported the 58 truth and the false report rate belief was set to 0%. Six measurements are taken, one every three seconds, and the particle distribution is shown after each measurement. The target is found in the rst measurement, and the initial contraction and then gradual expansion of the particle cloud over time is shown. Particles of signi cant weight are blue and particles of lower weight are red. The rectangular region where 400 x 600 and 0 y 200 is measured at each measurement. Both thresholds cause the particle distribution to contract to within the measured region after the rst measurement. In addition, both thresholds cause the distribution to not be resampled after the second measurement. After the fourth measurement, the distribution is not resampled in the case of the 80% threshold, while it is resampled under a 90%N threshold. In both cases, when the target leaves the measured area and is not found, particles are not preserved in the region of non detection. Because the target is near the edge of the measured region but not detected after the fth measurement, it would be bene cial to preserve some particles in the region of non-detection, as this could improve tracking performance. In addition to spatial resolution, tracking performance and time to detect data was collected. The sensor always reported the truth, and the particle distribution was only re- sampled when the number of e ective particles fell below the speci ed threshold. Thresholds for the number of e ective particles of 80%N and 90%N were used, and the performance was measured in the case of both a winding and a straight target path. 59 (a) Measurement 1: target found, resample (b) Measurement 2: target found, no resample (c) Measurement 3: target found, resample (d) Measurement 4: target found, resample (e) Measurement 5: target not found, resample (f) Measurement 6: target not found, resample Figure 3.31: E ect of the number of e ective particles on particle spatial resolution, Ne > 90%N 60 (a) Measurement 1: target found, resample (b) Measurement 2: target found, no resample (c) Measurement 3: target found, resample (d) Measurement 4: target found, no resample (e) Measurement 1: target not found, resample (f) Measurement 6: target not found, resample Figure 3.32: E ect of the number of e ective particles on particle spatial resolution, Ne > 80%N 61 0 10 20 30 40 50 60 0 0.5 1 1.5 2 2.5 3 3.5 4 x 10 5 Time (sec) Mean Square Error (m 2 ) Dispersion 100% N Dispersion 90% N Dispersion 80% N Figure 3.33: Dispersion Model, Average Mean Square Error After Target Found, Straight Path, 80%N and 90%N E ective Particle Thresholds 0 10 20 30 40 50 60 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 x 10 4 Time (sec) Mean Square Error (m 2 ) Traffic 100% N Traffic 90% N Traffic 80% N Figure 3.34: Tra c Motion Model, Average Mean Square Error After Target Found, Straight Path, 80%N and 90%N E ective Particle Thresholds 62 0 10 20 30 40 50 60 05 1015 x 10 4 Time (sec) Mean Square Error (m 2 ) Dispersion 100% N Dispersion 90% N Dispersion 80% N Figure 3.35: Dispersion Model, Average Mean Square Error After Target Found, Winding Path, 80%N and 90%N E ective Particle Thresholds 0 10 20 30 40 50 60 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 x 10 4 Time (sec) Mean Square Error (m 2 ) Traffic 100% N Traffic 90% N Traffic 80% N Figure 3.36: Tra c Motion Model, Average Mean Square Error After Target Found, Winding Path, 80%N and 90%N E ective Particle Thresholds 63 Table 3.6: Time to Detect Results for Straight Path, N = 1000 Particle Ne Time to Model Threshold Detect Dispersion 100% N 53:59 sec Tra c 100% N 29:31 sec Dispersion 90% N 56:61 sec Tra c 90% N 35:38 sec Dispersion 80% N 60:93 sec Tra c 80% N 29:89 sec Table 3.7: Time to Detect Results for Winding Path, N = 1000 Particle Ne Time to Model Threshold Detect Dispersion 100% N 66:13 sec Tra c 100% N 69:54 sec Dispersion 90% N 73:02 sec Tra c 90% N 44:83 sec Dispersion 80% N 65:64 sec Tra c 80% N 50:95 sec 3.7 Summary of Results A high delity tra c model was developed to propagate a particle lter in time in an urban environment to track a ground target vehicle. The sensor model used to locate and track the target vehicle was non-di erentiable and provided only a binary response of target presence within a large region. In all cases, it is observed that the Tra c Motion Model pro- vides for superior target tracking, and computational time may be appropriately decreased without signi cant loss of delity by adjusting the number of particles used, keeping spatial resolution in mind. These results are encouraging, given the ease of introducing a false report belief into the particle lter framework and when paired with the previously obtained tracking performance 64 results. By decreasing the particle lter?s con dence in the sensor, a softer pdf is obtained. This forces a broader spatial resolution of particle cloud and decreases the probability of losing the target after it is detected. It is nonintuitive to detune a lter to improve per- formance, but because the sensor model in this problem requires an added focus on spatial resolution, it is necessary. This result also brings to light the coupling between the performance due to the dy- namics model, sensor model, and particle management. Chapter 4 will explore the e ect of an imperfect sensor on the performance of each model. The false report rate belief may be tuned to maintain the spatial resolution of the particle distribution after a false detection is reported. In addition, decreasing the required number of e ective particles could provide for su cient spatial resolution in the presence of a missed-detection. Chapter 4 discusses this in a structured and systematic approach to tuning the spatial resolution to achieve improved performance in the presence of various sensor models. The importance of the dynamic model has been solidi ed in the presence of a perfect sensor. 65 Chapter 4 Accounting for an Imperfect Sensor Following the validation of the simulation?s stability and the e ectiveness of the Tra c Motion Model, a more realistic sensor model was introduced. Because the regional sensor is modeled after a human operator, the possibility of a false report must be addressed. The idea of adjusting the particle lter?s belief in the sensor?s false report rate was introduced in Section 3.6.2. By increasing the particle lter?s belief in the false report rate without actually simulating false measurements, the spatial resolution of the particle cloud was preserved such that the target was not always lost when it was just outside of a measured region. This decreased the time to detect the target. The bene t provided by adjusting the believed false report rate provided insight into how to improve performance in the presence of actual false reports. The previous chapter indicates the vital role of the motion model in tracking perfor- mance and the exibility of the particle lter framework when sensor data is true. The results given in this chapter indicate the importance of the measurement update in the presence of untrustworthy measurements. A thorough investigation into the tuning of the existing parameters in the measurement update was done in order to establish the e ect of each parameter. The detrimental impact of accepting a false report to be true is presented, in addition to the idea of establishing a risk metric that could be used to determine if a measurement is likely false. The sensor model was modi ed to provide false reports at some measurements. For example, if a sensor is said to have a 10% false report rate (Rfalse), each measurement has a 10% chance of reporting the opposite of the truth. This results in identical values for missed detection and false detection rates. False measurements are achieved through the 66 following steps: measure a region, determine the truth, take a random draw from a uniform distribution between 0 and 1, if the random value is greater than 1 - Rfalse, the opposite of the truth is reported. Figure 4.1 illustrates the diminished tracking performance that results from the intro- duction of a sensor with a 10% false report rate. It should be noted that in each case, the particle lter believed the sensor had a false report rate of 10%. The particle lter?s belief in the false report rate a ects how much a particle?s weight is changed based on a measurement. It does not a ect if the particle lter accepts or rejects a measurement. Each measurement is accepted to be true. As mentioned in Section 3.5, the tracking performance is measured for 60 seconds following the rst time a detection is reported. By introducing a 10% false report rate into the sensor model, the tracking performance could begin being measured before the target is actually found. In the 10% false report example in Fig. 4.1, 60% of the initial \detection" reports were false reports. This resulted in seemingly diminished tracking performance when compared to a perfect sensor. The results in Table 4.1 indicate that the average values of detection time decreased, as would be expected with the possibility of a reported detection actually being a false positive. For example, if measurements are taken every 3 seconds by a sensor with a false positive rate of 10%, one would expect a maximum time to detect around 30 seconds. This measure is less than the time to detect for the Dispersion and Tra c Motion Models presented in Section 3.5, which indicates a high probability of a false positive occurring before the target is actually found. Trusting a false positive is a signi cant cause of decreased tracking performance. Sec- tions 4.1 and 4.2 examine the e ect of tuning parameters that exist within the standard particle lter framework in an e ort to improve tracking performance in the presence of false reports. Section 4.1 investigates the measurement update step and the e ect of the particle lter?s belief of the sensor?s false report rate. Section 4.2 studies the e ect of adjusting the resampling step by tuning the e ective number of particles. Both parameters play into the 67 0 10 20 30 40 50 60 0 0.5 1 1.5 2 2.5 3 3.5 4 x 10 5 Time (sec) Mean Square Error (m 2 ) 10% False Report Perfect Sensor Figure 4.1: Impact of a Sensor with 10% False Report Rate, Tra c Model, Winding Path, N = 1000, False Report Belief 10% importance of maintaining proper spatial resolution of the particle cloud based on a sensor?s perceived accuracy. After the most bene cial values for the sensor belief and number of e ective particles is determined for each sensor model, a risk assessment is introduced. The practice of assessing the cost of resampling after a measurement and accepting a positive measurement through a Bayesian risk analysis is discussed in Section 4.3. 4.1 Sensor Belief Study As previously discussed in Section 3.6.2, the inclusion of a believed false report rate into the particle lter?s measurement update caused particle weights in regions of non-detection to be signi cantly reduced, but not always eliminated in the resampling process. Conversely, some particles outside of a region of detection were no longer immediately eliminated until after a repeated detection. Recall, the particle update process directly involves the particle 68 Table 4.1: Time to Detect Results for Winding Path, N = 1000, False Report Belief 10% Particle Model Target Path Sensor Model Time to Detect Dispersion Straight Perfect 57:9 sec Dispersion Winding Perfect 67:7 sec Tra c Straight Perfect 28:8 sec Tra c Winding Perfect 32:2 sec Dispersion Straight 10% False Report 17:43 sec Dispersion Winding 10% False Report 20:9 sec Tra c Straight 10% False Report 16:5 sec Tra c Winding 10% False Report 22:35 sec lter?s belief in the sensor?s false report rate, Rfalse. p(zk jxk) = 1 Rfalse (4.1) p(xk jZk 1) = NinX i=1 wiin (4.2) p(zk jZk 1) = (1 Rfalse) NinX i=1 wiin + (4.3) (Rfalse)(1 NinX i=1 wiin) p(xk jzk) = p(zk jxk)p(xk jZk 1)p(z k jZk 1) (4.4) wik = w i kp(xk jzk)P win (4.5) An investigation into the e ects of the sensor accuracy on the amount of time required to detect the target and the tracking performance was included in this study. Three sensor models were used for comparison: a perfect sensor, a sensor with 10% false report rate, and a sensor with 20% false report rate. Figures are shown for the converged average error over 200 Monte Carlo runs. The target traversed both a straight and a winding path, and was initialized the same way at the start of each run. Particles were not redistributed in this study. 69 The following gures illustrate that regardless of the sensor model, the Tra c Motion Model tracks the target better than the Dispersion Motion Model. Some additional consid- eration will be given to the \optimal" value for the particle lter?s belief of the false alarm rate as the target?s path is never truly known. 0 10 20 30 40 50 60 0 0.5 1 1.5 2 2.5 3 3.5 x 10 5 Time (sec) Mean Square Error (m 2 ) 0% False Report 5% False Report 10% False Report 15% False Report 20% False Report 25% False Report 30% False Report 35% False Report 40% False Report 45% False Report 50% False Report Figure 4.2: Dispersion Model, Straight Path, False Report Belief Study, Perfect Sensor 70 0 10 20 30 40 50 60 0 0.5 1 1.5 2 2.5 3 3.5 x 10 5 Time (sec) Mean Square Error (m 2 ) 0% False Report 5% False Report 10% False Report 15% False Report 20% False Report 25% False Report 30% False Report 35% False Report 40% False Report 45% False Report 50% False Report Figure 4.3: Tra c Model, Straight Path, False Report Belief Study, Perfect Sensor 40 45 50 55 60 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 x 10 4 Time (sec) Mean Square Error (m 2 ) 0% False Report 5% False Report 10% False Report 15% False Report 20% False Report 25% False Report 30% False Report 35% False Report 40% False Report 45% False Report 50% False Report Figure 4.4: Tra c Model, Straight Path, False Report Belief Study, Perfect Sensor, nal 20 seconds 71 The gures above indicate that in the presence of a perfect sensor, improved tracking performance is obtained when the particle lter assumes a false report rate exists, regardless of motion model. The Dispersion Model yields slightly improved performance with a false report rate belief of 5% and 10%. As expected, the Tra c Motion Model provides improved tracking performance when compared to the Dispersion Model, and with the inclusion of a false report rate belief of 15%. All of the results follow the same trend, except for the 50% false report belief cases. In those cases, no bene t is given to region of detection, therefore additional measurements do not improve tracking performance. Table 4.2: Time to Detect Results, Straight Path, Perfect Sensor, False Report Belief 0%- 50% False Report Time to Detect Time to Detect Belief Dispersion Tra c 0% 57:93 sec 34:41 sec 5% 51:42 sec 30:91 sec 10% 49:23 sec 32:64 sec 15% 53:41 sec 36:19 sec 20% 51:91 sec 31:41 sec 25% 44:86 sec 33:21 sec 30% 50:86 sec 33:55 sec 35% 51:33 sec 33:45 sec 40% 47:53 sec 33:65 sec 45% 52:87 sec 34:66 sec 50% 284:58 sec 208:36 sec The results in Table 4.2 indicate that some bene t is obtained by including a false report belief in the measurement updated in the case of both motion models. The improved tracking capabilities yielded by the Tra c Motion Model over the Dispersion Model consistently aids in reducing the time to detect the target. In addition, the 50% false report belief results indicate a signi cantly increased time to detect, as the particle weights are not altered signi cantly based on any measurement. 72 0 10 20 30 40 50 60 0 0.5 1 1.5 2 2.5 3 3.5 4 x 10 5 Time (sec) Mean Square Error (m 2 ) 0% False Report 5% False Report 10% False Report 15% False Report 20% False Report 25% False Report 30% False Report 35% False Report 40% False Report 45% False Report 50% False Report Figure 4.5: Dispersion Model, Winding Path, False Report Belief Study, Perfect Sensor 0 10 20 30 40 50 60 0 0.5 1 1.5 2 2.5 3 3.5 4 x 10 5 Time (sec) Mean Square Error (m 2 ) 0% False Report 5% False Report 10% False Report 15% False Report 20% False Report 25% False Report 30% False Report 35% False Report 40% False Report 45% False Report 50% False Report Figure 4.6: Tra c Model, Winding Path, False Report Belief Study, Perfect Sensor 73 40 45 50 55 60 0.5 1 1.5 2 2.5 3 3.5 4 4.5 x 10 4 Time (sec) Mean Square Error (m 2 ) 0% False Report 5% False Report 10% False Report 15% False Report 20% False Report 25% False Report 30% False Report 35% False Report 40% False Report 45% False Report 50% False Report Figure 4.7: Tra c Model, Winding Path, False Report Belief Study, Perfect Sensor, Final 20 seconds In agreement with the results in the case of a straight path, the gures above indicate that the Tra c Model outperforms the Dispersion Model?s tracking performance. In ad- dition, tuning the sensor?s belief to 0% false report rate to match the sensor?s actual false report rate provides improved tracking performance for both motion models. Tuning the false report rate belief to 10% also provides improved tracking in the case of the Tra c Motion Model, as consistent with the results in Chapter 3. The time to detect results are consistent with that of a straight path, in that the Dispersion Model provides for faster localization. However, both models bene t from the inclusion of a false report rate in the measurement update. This occurs because the winding path traverses through the same quadrant of the map for a signi cant amount of time. Therefore, if the target is located just outside of a measured region, it would bene t from some of the particles in that measured region being preserved. This would provide for improved particle spatial resolution in the proximity of the target and therefore decrease the time to detect the target. 74 Table 4.3: Time to Detect Results, Winding Path, Perfect Sensor, False Report Belief 0%- 50% False Report Time to Detect Time to Detect Belief Dispersion Tra c 0% 72:13 sec 38:85 sec 5% 73:43 sec 42:53 sec 10% 70:50 sec 39:61 sec 15% 63:06 sec 38:03 sec 20% 67:44 sec 40:96 sec 25% 62:20 sec 39:03 sec 30% 73:71 sec 38:40 sec 35% 62:28 sec 37:43 sec 40% 66:85 sec 40:96 sec 45% 74:14 sec 42:53 sec 50% 349:67 sec 265:15 sec The following gures and tables illustrate the e ect of including a sensor with a false report rate of 10% and 20%, respectively. Tracking performance is greatly diminished, and the mean square error is greater than the area of a regional sensor in all cases. False start rates are included for comparison. The false-start rate is de ned as the percentage of times the rst reported detection of the target is a false report, thus causing the tracking performance to begin being measured when the target was not actually found. 75 0 10 20 30 40 50 60 1 1.5 2 2.5 3 3.5 x 10 5 Time (sec) Mean Square Error (m 2 ) 0% False Report 5% False Report 10% False Report 15% False Report 20% False Report 25% False Report 30% False Report 35% False Report 40% False Report 45% False Report 50% False Report Figure 4.8: Dispersion Model, Straight Path, False Report Belief Study, 10% False Report Sensor Figure 4.8 illustrates that in the case of a target on a straight path and a 10% false report rate sensor, the inclusion of a 5% false report belief results in improved performance for the Dispersion Model. However, the best performers converge to an error that is close to an order of magnitude larger than that of the perfect sensor, and well outside the bounds of a sensor region. 76 0 10 20 30 40 50 60 1 1.21.41.61.8 2 2.22.42.62.8 3 x 10 5 Time (sec) Mean Square Error (m 2 ) 0% False Report 5% False Report 10% False Report 15% False Report 20% False Report 25% False Report 30% False Report 35% False Report 40% False Report 45% False Report 50% False Report Figure 4.9: Tra c Model, Straight Path, False Report Belief Study, 10% False Report Sensor Tuning the false report belief to that of the sensor?s actual false report rate provided for optimal tracking performance when compared to the other false report belief values for the Tra c Motion Model. Although the performance is much worse than that with results from having perfect sensor data, the Tra c Motion Model provides for better target tracking than the Dispersion Model, despite a consistently higher false start rate, as given by Table 4.4. The higher false start rate was expected for the Tra c Motion Model in the presence of a sensor with a false report rate, as its time to detect with a perfect sensor was higher than that of the Dispersion Model for a perfect sensor. Therefore, because it would take longer to localize the target given a perfect sensor, the probability of receiving a false report before the target is actually found is higher than that of the Dispersion Model. 77 Table 4.4: Time to Detect Results, Straight Path, 10% False Report Sensor, False Report Belief 0%-50% False Report Time to Detect % False Starts Time to Detect %False Starts Belief Dispersion Dispersion Tra c Tra c 0% 18:58 sec 61:50 19:84 sec 65:00 5% 18:57 sec 55:00 20:26 sec 64:50 10% 20:04 sec 55:00 19:56 sec 57:50 15% 18:78 sec 63:50 21:25 sec 63:50 20% 20:01 sec 56:00 17:56 sec 67:50 25% 18:61 sec 62:00 21:21 sec 62:50 30% 20:07 sec 60:00 19:08 sec 62:50 35% 19:05 sec 58:50 19:77 sec 68:50 40% 17:51 sec 69:50 19:47 sec 58:50 45% 20:35 sec 65:00 21:77 sec 64:50 50% 28:02 sec 85:00 27:45 sec 90:00 0 10 20 30 40 50 60 1 1.5 2 2.5 3 3.5 4 4.5 x 10 5 Time (sec) Mean Square Error (m 2 ) 0% False Report 5% False Report 10% False Report 15% False Report 20% False Report 25% False Report 30% False Report 35% False Report 40% False Report 45% False Report 50% False Report Figure 4.10: Dispersion Model, Winding Path, False Report Belief Study, 10% False Report Sensor 78 0 10 20 30 40 50 60 1 1.5 2 2.5 3 3.5 4 x 10 5 Time (sec) Mean Square Error (m 2 ) 0% False Report 5% False Report 10% False Report 15% False Report 20% False Report 25% False Report 30% False Report 35% False Report 40% False Report 45% False Report 50% False Report Figure 4.11: Tra c Model, Winding Path, False Report Belief Study, 10% False Report Sensor The results for a winding path are similar to those of a straight path. The Tra c Motion Models shows improved performance when tuned to have a 25% false report belief, while the Dispersion Model bene ts from a 5% false report belief. Again, the tracking performance is substantially worse than when the truth is always reported. Despite the higher false start rate, the Tra c Motion Model produces better tracking results than the Dispersion Model in the case of a winding path with a 10% false report rate sensor model. 79 Table 4.5: Time to Detect Results, Winding Path, 10% False Report Sensor, False Report Belief 0%-50% False Report Time to Detect % False Starts Time to Detect %False Starts Belief Dispersion Dispersion Tra c Tra c 0% 20:78 sec 55:50 20:79 sec 62:50 5% 20:73 sec 57:00 21:37 sec 70:00 10% 17:76 sec 65:00 20:91 sec 73:50 15% 18:28 sec 61:00 19:80 sec 61:00 20% 20:91 sec 59:50 21:06 sec 71:50 25% 18:36 sec 65:00 20:61 sec 64:50 30% 19:96 sec 69:00 20:02 sec 64:50 35% 20:37 sec 65:00 22:00 sec 65:50 40% 20:52 sec 62:00 21:19 sec 65:00 45% 21:03 sec 58:00 18:49 sec 69:00 50% 29:05 sec 94:50 31:25 sec 92:50 0 10 20 30 40 50 60 1 1.5 2 2.5 3 3.5 x 10 5 Time (sec) Mean Square Error (m 2 ) 0% False Report 5% False Report 10% False Report 15% False Report 20% False Report 25% False Report 30% False Report 35% False Report 40% False Report 45% False Report 50% False Report Figure 4.12: Dispersion Model, Straight Path, False Report Belief Study, 20% False Report Sensor 80 0 10 20 30 40 50 60 2.5 3 3.5 4 x 10 5 Time (sec) Mean Square Error (m 2 ) 0% False Report 5% False Report 10% False Report 15% False Report 20% False Report 25% False Report 30% False Report 35% False Report 40% False Report 45% False Report 50% False Report Figure 4.13: Tra c Model, Straight Path, False Report Belief Study, 20% False Report Sensor The increase in the sensor?s false report rate from 10% to 20% nearly doubled the tracking error for both motion models. The gures above indicate that the Dispersion Model achieves increased performance when tuned to believe the sensor?s false report rate is 5% or 35%, but does not yield convergent tracking performance. Whereas, the Tra c Motion Model?s performance is heavily diminished in the presence of an average false start of 83%, it is still able to produce tracking error that decreases over time, and the best performer is that of a false report belief of 0% and 40%. 81 Table 4.6: Time to Detect Results, Straight Path, 20% False Report Sensor, False Report Belief 0%-50% False Report Time to Detect % False Starts Time to Detect %False Starts Belief Dispersion Dispersion Tra c Tra c 0% 13:06 sec 77:00 13:35 sec 82:00 5% 12:45 sec 80:00 14:79 sec 82:00 10% 12:36 sec 74:00 13:53 sec 84:00 15% 12:72 sec 77:50 12:15 sec 81:50 20% 13:26 sec 78:50 12:42 sec 83:50 25% 13:89 sec 76:00 14:01 sec 83:00 30% 11:24 sec 78:00 12:81 sec 87:50 35% 14:41 sec 76:00 13:72 sec 82:00 40% 11:53 sec 78:50 13:65 sec 84:50 45% 12:82 sec 76:50 13:68 sec 88:00 50% 13:72 sec 93:50 14:91 sec 90:50 0 10 20 30 40 50 60 2 2.5 3 3.5 4 4.5 5 x 10 5 Time (sec) Mean Square Error (m 2 ) 0% False Report 5% False Report 10% False Report 15% False Report 20% False Report 25% False Report 30% False Report 35% False Report 40% False Report 45% False Report 50% False Report Figure 4.14: Dispersion Model, Winding Path, False Report Belief Study, 20% False Report Sensor 82 0 10 20 30 40 50 60 2.42.62.8 3 3.23.43.63.8 4 x 10 5 Time (sec) Mean Square Error (m 2 ) 0% False Report 5% False Report 10% False Report 15% False Report 20% False Report 25% False Report 30% False Report 35% False Report 40% False Report 45% False Report 50% False Report Figure 4.15: Tra c Model, Winding Path, False Report Belief Study, 20% False Report Sensor Similar to results in the case of a straight path, the results obtained for a target on a winding path indicate that the Dispersion Model should be tuned to a 30% false report belief and the Tra c Model to 0% or 15% false report belief. Both motion models su ered from the elevated false start rate. The Tra c Motion model yielded convergent tracking error, while tracking error diverged in the case of the the Dispersion Model. Detuning the particle lter?s belief in the sensor?s false report rate provided increased spatial resolution of the particle cloud. This improved the time to detect metric in the winding target scenario, as it allowed for all particles in a region of non-detection to not be eliminated immediately, providing for increased particle presence in the proximity of a target just outside of a measured region. In the presence of a sensor with a false report rate, the tracking performance is signi cantly diminished, regardless of the motion model that is used. But, the Tra c Motion Model out performed the Dispersion Model?s tracking performance in the presence of each sensor model. Section 4.2 will explore the e ect of decreasing the 83 Table 4.7: Time to Detect Results, Winding Path, 20% False Report Sensor, False Report Belief 0%-50% False Report Time to Detect % False Starts Time to Detect %False Starts Belief Dispersion Dispersion Tra c Tra c 0% 13:05 sec 78:50 13:35 sec 87:50 5% 12:57 sec 74:50 14:40 sec 82:50 10% 13:62 sec 77:50 13:39 sec 83:00 15% 13:23 sec 70:50 12:75 sec 84:50 20% 12:34 sec 80:00 12:81 sec 86:50 25% 12:73 sec 75:50 13:57 sec 86:00 30% 11:92 sec 74:50 13:36 sec 83:00 35% 13:13 sec 73:00 12:36 sec 84:50 40% 12:34 sec 78:50 12:87 sec 86:50 45% 12:37 sec 75:50 13:29 sec 84:00 50% 14:13 sec 97:50 13:77 sec 95:50 required number of e ective particles to further study the e ect of spatial resolution on the tracking performance and false alarm rate in the presence of a sensor with a false report rate. 4.2 Number of E ective Particles Study For example, if a detection is reported, weight is shifted to the particles within the measured region and the region will be measured again, as it continues to maintain the highest particle weight sum. However, if the number of e ective particles remains below a desired threshold, the particle cloud is not resampled. If the particles are not resampled and the measurement was a false detection, the following measurement has a high probability of reporting the truth and the particle cloud will maintain its previous coverage of the simulation space. This will improve tracking performance, but will not provide for a means of rejecting false reports, as all measurements are still believed to be true. All simulations were done with 1000 particles, and all curves re ect the average value of the mean square error for a case. 84 0 10 20 30 40 50 60 0 0.5 1 1.5 2 2.5 3 3.5 4 x 10 5 Time (sec) Mean Square Error (m 2 ) 50% N 55% N 60% N 65% N 70% N 75% N 80% N 85% N 90% N 95% N 100% N Figure 4.16: Dispersion Model, Straight Path, Number of E ective Particles Study, Perfect Sensor 0 10 20 30 40 50 60 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 x 10 4 Time (sec) Mean Square Error (m 2 ) 50% N 55% N 60% N 65% N 70% N 75% N 80% N 85% N 90% N 95% N 100% N Figure 4.17: Tra c Model, Straight Path, Number of E ective Particles Study, Perfect Sensor 85 In the presence of a perfect sensor and a target on a straight path, improved tracking results are achieved for the Dispersion Model when the particle distribution is resampled after each measurement, as consistent with the Nthr = N. The Tra c Model displays its best performance when the threshold is set at 85%N. The detection times indicate that faster target localization results from resampling following each measurement and settingNthr = N. This result could be used in future work to adapt the required number of e ective particles to the phase of the application: localization or tracking. Table 4.8: Time to Detect Results, Straight Path, Perfect Sensor, Nthr = (100% 50%)N Resample Time to Detect Time to Detect Threshold Dispersion Tra c 100% N 51:71 sec 33:33 sec 95% N 58:87 sec 32:86 sec 90% N 56:61 sec 35:38 sec 85% N 57:54 sec 35:97 sec 80% N 60:93 sec 29:89 sec 75% N 50:64 sec 36:54 sec 70% N 60:41 sec 34:14 sec 65% N 48:47 sec 36:91 sec 60% N 52:09 sec 38:07 sec 55% N 53:19 sec 30:39 sec 50% N 48:34 sec 27:73 sec 86 0 10 20 30 40 50 60 05 1015 x 10 4 Time (sec) Mean Square Error (m 2 ) 50% N 55% N 60% N 65% N 70% N 75% N 80% N 85% N 90% N 95% N 100% N Figure 4.18: Dispersion Model, Winding Path, Number of E ective Particles Study, Perfect Sensor 0 10 20 30 40 50 60 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 x 10 4 Time (sec) Mean Square Error (m 2 ) 50% N 55% N 60% N 65% N 70% N 75% N 80% N 85% N 90% N 95% N 100% N Figure 4.19: Tra c Model, Winding Path, Number of E ective Particles Study, Perfect Sensor 87 In the case of a winding path, the Dispersion Model bene ts from tuning its required number of e ective particles to 95%N while the Tra c Model bene ts from a 90%N threshold. The times to detect at the tuned value of required number of e ective particles are similar to the minimal value for each motion model in the case of a winding target path. Table 4.9: Time to Detect Results, Winding Path, Perfect Sensor, Nthr = (100% 50%)N, False Report Belief 0% Resample Time to Detect Time to Detect Threshold Dispersion Tra c 100% N 59:37 sec 42:11 sec 95% N 60:99 sec 44:55 sec 90% N 73:02 sec 44:83 sec 85% N 53:10 sec 44:46 sec 80% N 65:64 sec 50:95 sec 75% N 54:11 sec 46:13 sec 70% N 64:29 sec 45:18 sec 65% N 61:65 sec 42:63 sec 60% N 61:88 sec 45:71 sec 55% N 68:51 sec 50:26 sec 50% N 58:06 sec 49:74 sec The inclusion of a sensor with false report rate is expected to diminish tracking per- formance and reduce the detection time, as false positives are believed true. The e ect of reducing the required number of e ective particles is illustrated by the gures below. A sensor with a false report rate of 10% and 20% were studied, and in each case, the particle lter?s belief of the sensor?s false report rate was set to its actual false report rate. 88 0 10 20 30 40 50 60 1 1.5 2 2.5 3 3.5 x 10 5 Time (sec) Mean Square Error (m 2 ) 50% N 55% N 60% N 65% N 70% N 75% N 80% N 85% N 90% N 95% N 100% N Figure 4.20: Dispersion Model, Straight Path, Number of E ective Particles Study, 10% False Report Sensor 0 10 20 30 40 50 60 1 1.21.41.61.8 2 2.22.42.62.8 3 x 10 5 Time (sec) Mean Square Error (m 2 ) 50% N 55% N 60% N 65% N 70% N 75% N 80% N 85% N 90% N 95% N 100% N Figure 4.21: Tra c Model, Straight Path, Number of E ective Particles Study, 10% False Report Sensor 89 Although the tracking performance is much weaker than that resulting from the use of a perfect sensor, the number of e ective particles play a role in the spatial resolution of the particle cloud and can increase tracking performance. The Dispersion Model does not bene t from the reduction in the amount of required e ective particles, but the Tra c Motion Model achieves improved tracking performance with a required number of e ective particles of 50%N. Table 4.10: Time to Detect Results, Straight Path, Perfect Sensor, Nthr = (100% 50%)N, False Report Belief 10% Resample Time to Detect % False Time to Detect % False Starts Threshold Starts Dispersion Tra c Tra c 100% N 15:60 sec 55.50 21:37 sec 56.50 95% N 21:03 sec 57.00 20:45 sec 58.00 90% N 16:95 sec 56.50 18:06 sec 57.00 85% N 17:53 sec 63.50 20:23 sec 60.50 80% N 19:21 sec 58.00 23:34 sec 62.50 75% N 19:06 sec 50.00 20:33 sec 63.00 70% N 19:91 sec 58.00 20:59 sec 61.00 65% N 20:52 sec 59.00 19:56 sec 60.00 60% N 17:12 sec 53.50 19:59 sec 66.00 55% N 16:98 sec 54.50 18:73 sec 59.50 50% N 16:78 sec 55.50 19:98 sec 57.50 90 0 10 20 30 40 50 60 1.61.8 2 2.22.42.62.8 3 3.23.4 x 10 5 Time (sec) Mean Square Error (m 2 ) 50% N 55% N 60% N 65% N 70% N 75% N 80% N 85% N 90% N 95% N 100% N Figure 4.22: Dispersion Model, Winding Path, Number of E ective Particles Study, 10% False Report Sensor 0 10 20 30 40 50 60 1.41.61.8 2 2.22.42.62.8 3 3.2 x 10 5 Time (sec) Mean Square Error (m 2 ) 50% N 55% N 60% N 65% N 70% N 75% N 80% N 85% N 90% N 95% N 100% N Figure 4.23: Tra c Model, Winding Path, Number of E ective Particles Study, 10% False Report Sensor 91 In the case of a target on a winding path, the Dispersion Model exhibits its best tracking results when the required number of particles is set to 70%N, whereas the Tra c Model does not bene t from a decreased number of e ective particles. Table 4.11: Time to Detect Results, Winding Path, 10% False Report Sensor, Nthr = (100% 50%)N, False Report Belief 10% Resample Time to Detect % False Time to Detect % False Starts Threshold Starts Dispersion Tra c Tra c 100% N 18:19 sec 59.00 18:06 sec 60.50 95% N 19:45 sec 55.00 21:49 sec 69.50 90% N 18:52 sec 59.00 22:58 sec 66.00 85% N 20:85 sec 55.50 21:99 sec 72.00 80% N 17:59 sec 54.50 18:93 sec 68.50 75% N 20:96 sec 61.50 18:47 sec 62.50 70% N 15:42 sec 57.00 19:85 sec 66.50 65% N 18:05 sec 56.50 22:83 sec 59.50 60% N 17:97 sec 55.50 22:75 sec 68.50 55% N 19:80 sec 52.50 22:38 sec 64.00 50% N 19:43 sec 56.50 19:14 sec 69.50 92 0 10 20 30 40 50 60 1.21.41.61.8 2 2.22.42.62.8 3 3.2 x 10 5 Time (sec) Mean Square Error (m 2 ) 50% N 55% N 60% N 65% N 70% N 75% N 80% N 85% N 90% N 95% N 100% N Figure 4.24: Dispersion Model, Straight Path, Number of E ective Particles Study, 20% False Report Sensor 0 10 20 30 40 50 60 1 1.5 2 2.5 3 3.5 x 10 5 Time (sec) Mean Square Error (m 2 ) 50% N 55% N 60% N 65% N 70% N 75% N 80% N 85% N 90% N 95% N 100% N Figure 4.25: Tra c Model, Straight Path, Number of E ective Particles Study, 20% False Report Sensor 93 The gures above indicate that the Dispersion Model achieves improved tracking per- formance with a 70%N e ective particle threshold, while the Tra c Motion Model performs best with a 50% e ective particle threshold. The false alarm rates and diminished tracking performance are consistent with the results in Section 4.1. It remains true that believing a false report to be true is damaging to the tracking performance. Table 4.12: Time to Detect Results, Straight Path, 20% False Report Sensor, Nthr = (100% 50%)N, False Report Belief 20% Resample Time to Detect % False Time to Detect % False Starts Threshold Starts Dispersion Tra c Tra c 100% N 12:36 sec 78.50 12:49 sec 77.50 95% N 12:55 sec 76.00 13:45 sec 78.50 90% N 12:07 sec 82.00 11:52 sec 83.00 85% N 12:16 sec 76.00 11:34 sec 85.50 80% N 12:60 sec 81.00 13:11 sec 79.00 75% N 12:54 sec 80.50 12:84 sec 78.00 70% N 12:81 sec 75.00 13:12 sec 76.00 65% N 11:43 sec 80.50 12:46 sec 77.50 60% N 11:73 sec 73.50 13:11 sec 82.50 55% N 11:05 sec 76.50 11:83 sec 82.50 50% N 12:30 sec 80.00 14:48 sec 79.50 94 0 10 20 30 40 50 60 2 2.5 3 3.5 4 4.5 x 10 5 Time (sec) Mean Square Error (m 2 ) 50% N 55% N 60% N 65% N 70% N 75% N 80% N 85% N 90% N 95% N 100% N Figure 4.26: Dispersion Model, Winding Path, Number of E ective Particles Study, 20% False Report Sensor 0 10 20 30 40 50 60 2 2.22.42.62.8 3 3.23.4 x 10 5 Time (sec) Mean Square Error (m 2 ) 50% N 55% N 60% N 65% N 70% N 75% N 80% N 85% N 90% N 95% N 100% N Figure 4.27: Tra c Model, Winding Path, Number of E ective Particles Study, 20% False Report Sensor 95 Similar to the results for a target on a straight path, the Dispersion Model bene ts from a 75%N e ective particle threshold and the Tra c Model bene ts from a 50%N e ective particle threshold. Table 4.13: Time to Detect Results, Winding Path, 20% False Report Sensor, Nthr = (100% 50%)N, False Report Belief 20% Resample Time to Detect % False Time to Detect % False Starts Threshold Starts Dispersion Tra c Tra c 100% N 11:89 sec 79.50 12:79 sec 78.00 95% N 13:75 sec 79.00 13:90 sec 81.00 90% N 12:81 sec 80.00 12:90 sec 82.00 85% N 13:20 sec 75.50 12:63 sec 85.00 80% N 14:75 sec 76.00 11:20 sec 81.50 75% N 12:87 sec 79.50 11:70 sec 83.00 70% N 13:58 sec 78.00 13:87 sec 81.50 65% N 14:05 sec 77.00 13:74 sec 83.50 60% N 13:38 sec 81.00 12:15 sec 84.00 55% N 12:27 sec 83.50 12:70 sec 81.50 50% N 11:86 sec 78.50 14:22 sec 78.50 The number of e ective particles may be tuned to improve performance, however false sensor reports continue to diminish tracking performance. Decreasing the required number of e ective particles reduces the frequency of resampling and therefore forces a broader probability density function and widens the spatial resolution of the distribution. When a false detection is reported, signi cant weight is shifted to the measured region. That region will be measured at the next measurement as it contains the highest sum of particle weight. If a false detection is reported and the required number of e ective particles is low (50- 75%N), signi cant particle coverage outside of the measured region is maintained. Upon the next measurement of the detection region, it is likely that another false detection will not be reported, causing particle weight to be shifted from that region back to the particles outside the measured region that were preserved by not resampling. Then, the regional sensor polling problem continues and other areas of the map are searched for the target. In this situation, however, the tracking performance would have begun to be measured at the 96 report of the false detection. This results in perceived poor tracking performance, as the target was not actually found at the start of the 60 second observation period. After a false detection is reported, the 60 second observation period is often spent searching the rest of the map for the target, therefore preventing tracking performance similar to that of a perfect sensor. The frequency of this increases as the false report rate of the sensor increases. Because of the signi cant e ect of false reports on the tracking performance, an addi- tional means of improving performance was investigated. A measure of risk was developed to determine the potential impact of accepting a measurement and resampling the particle distribution. By the inclusion of a means to reject a possibly false measurement, the im- portance of the number of e ective particles and sensor belief will come to light. Without reducing the number of e ective particles and sensor belief in the presence of a sensor with a false report rate, the rst false detection would result in completely loss of particle resolu- tion outside of the measured region, and therefore would result in a severely increased true detection time. 4.3 Bayesian Risk A Bayesian Risk Assessment was used to evaluate the risk of accepting a measurement as true, based on an example where the expected change in the estimate variance was used to determine the cost of a accepting a measurement [33]. In this application, the expected ux of particle weight into or out of the measured region will be used as the cost of resampling after a measurement. In the formulation presented by O?Reilly in [33], two hypotheses are considered: the measurement is true (H0) and the measurement has been altered (H1). A ratio of the likelihoods of each hypothesis given the data, d, is computed and denoted as LR, the likelihood ratio. LR = p(djH1)p(djH 0) (4.6) 97 In order to determine which hypothesis to accept, the likelihood ratio is compared to a threshold, denoted . if p(djH1)p(djH 0) < ; then H1 is accepted (4.7) if p(djH1)p(djH 0) ; then H0 is accepted The threshold may be formulated in the Bayesian framework based on the cost of accepting a certain measurement and the probability associated with each hypothesis. = (C10 C00)PH0(C 01 C11)PH1 (4.8) The cost variable CXY is the cost of accepting hypothesis HX if the hypothesis HY is actually true. There is usually no cost associated with accepting a true hypothesis, in which case the terms C00 and C11 are zero. The framework outlined by the previous equations provide a means for assigning a value of risk to each measurement. Risk is based on the likelihood of a measurement being true, taking into account the amount of particle weight that would shift if a measurement is accepted. If a measurement is deemed to risky to accept, the particle lter accepts the measurement but does not resample the distribution. This preserves all particles outside of the measured region. The risk assessment was used to control when to resample instead of whether or not to accept a measurement. This was done because if a measurement is correct, a repeated measurement will con rm it and when the risk decreases to an acceptable level, the distribution is resampled. If a measurement was false, a repeated measurement will often bring this to light and the spatial resolution of the particle cloud will not have been compromised. To examine the risk, the likelihood ratio must be established in terms associated with the particle lter. The likelihood ratio is computed based on the expressions below for the likelihoods that the target is actually present within the measured region (p(present)) or not 98 (p(not present)). p(present) = (1 Rfalse) NinX i=1 wiin + (Rfalse)(1 NinX i=1 wiin) (4.9) p(not present) = (Rfalse) NinX i=1 wiin + (1 Rfalse)(1 NinX i=1 wiin) When a detection is reported, the likelihood ratio is given by: LR = p(not present)p(present) (4.10) When a non-detection is reported, the likelihood ratio is given by: LR = p(present)p(not present) (4.11) Expressions for the cost of accepting a measurement are derived below. The terms C00 and C11 in Eq. (4.8) are set to zero because there is no cost associated with accepting a true measurement. The terms PH0 and PH1 depend on the received measurement. If a detection is reported, PH0 is the sum of the weights of the particles inside the measured region, and PH1 is the sum of the weights of the particles outside of the measured region. If a non-detection is reported, PH0 and PH1 are the sum of the weights of the particles outside the measured region and the sum of the particles outside the measured region, respectively. The cost associated with accepting a false measurement is quanti ed by the amount of particle weight that would shift if a measurement was accepted. In the case of a detection 99 being reported, the cost of accepting the measurement if it is false is computed by: in = (PNini=1 wiin) p(+jpresent) Nin (4.12) out = (1 PNini=1 wiin p( jpresent) Nout (4.13) = out Pw in in(1 PNini=1 wiin) (4.14) In the case of a non-detection being reported, the cost of accepting the measurement if it is false is determined by the following expressions: in = (1 PNini=1 wiin) p( jnot present) Nin (4.15) out = PNini=1 wiin p(+jnot present) Nout (4.16) = out(1 PNin i=1 w i in) inPNini=1 wiin (4.17) In a situation where all particle weight is located in a region and a non-detection is reported, it is obvious that the measurement is most likely false. Also, if a region of very low particle weight is polled and a detection is reported, it would be risky to believe such a measurement. However, in such instances where there is signi cant particle weight in each region, it is not immediately clear if a measurement should be accepted or rejected. An example of such a case is outlined for completion. The risk associated with a reported detection and a reported non-detection will be examined if a measurement is taken of the region in the bottom right corner in Fig. 4.28. In the example illustrated by Fig. 4.28, it is believed that the sensor has a false report rate of 10% and all particles have equal weight of 1/8. 100 Figure 4.28: Measurement Update with Risk Assessment Example Scenario If a detection is reported in the bottom right region of Fig. 4.28, the likelihood ratio is computed as follows: p(present) = (0:9)(3=8) + (0:1)(5=8) (4.18) p(not present) = (0:1)(3=8) + (0:9)(5=8) LR = p(not present)p(present) = 1:5 The cost ratio, , is then calculated. p(+jpresent) = (1 Rfalse) Pw in p(present) = 0:8438 (4.19) in = j( Pw in) p(+jpresent)j Nin = 0:1563 p( jpresent) = (1 Rfalse)(1 Pw in) p(not present) = 0:9375 (4.20) out = j(1 Pw in p( jnot present)j Nout = 0:0625 = out Pw in in(1 Pwin) = 0:2400 Because is less than the likelihood ratio, the particle cloud would not be resampled as the risk is too high. If the sensor reports that the target is not detected in the bottom right corner region, that is a false report. However, there is signi cant particle presence in other sensor regions, 101 so it is not immediately clear to the particle lter that the measurement is not likely true. The likelihood ratio is the inverse of that of a detection, and is therefore 2/3. The cost ratio, , is then calculated. p( jnot present) = Rfalse Pw in p(not present) = 0:0625 (4.21) in = (1 PNini=1 wiin) p( jnot present) Nout = 0:1125 p(+jpresent) = Rfalse(1 PNin i=1 w i in) p(not present) = 0:0625 (4.22) out = (1 PNini=1 wiin p( jpresent) Nin = 0:0729 = out(1 PNin i=1 w i in) inPNini=1 wiin = 1:0802 Because the likelihood ratio is less than , the particle distribution would be resampled if the number of e ective particles is below the desired threshold. 4.3.1 Risk Assessment Results The Bayesian Risk Assessment discussed above was implemented in the cases of the two sensors with false report rates. The following values were used for the false report belief and the number of e ective particles: The values in Table 4.14 were determined in Chapter 4. It is evident that the number of e ective particles should be reduced in the presence of an imperfect sensor and the false report belief should be set to a value greater than or equal to the true false report rate. Although the results in Section 4.1 indicate that the false report belief for the Tra c Motion Model in the presence of a sensor with a 20% false repot rate should be set to 0%, the false report belief for the Tra c Motion Model was set to 20% in this section. This was done because in Section 4.2, the false report belief was set to the true false report rate of each sensor, and superior tracking performance occurred in the case of the 20% false report sensor. 102 Table 4.14: Tuned Values for Sensor Belief and Number of E ective Particles Motion Target Sensor False Report Number of Model Path Model Belief E ective Particles Dispersion Straight 10% False Report 5% 60%N Dispersion Winding 10% False Report 5% 70%N Tra c Straight 10% False Report 10% 50%N Tra c Winding 10% False Report 25% 65%N Dispersion Straight 20% False Report 35% 55%N Dispersion Winding 20% False Report 30% 60%N Tra c Straight 20% False Report 20% 50%N Tra c Winding 20% False Report 20% 50%N In an e ort to reduce the false alarm percentage, the time to detect was not established until two consecutive detection measurements were reported. Any time a false alarm is recorded, the second detection in the sequence is a false detection. Data regarding the time to detect the target, tracking performance, and the decrease in false start rates is presented to show the importance of the Bayesian Risk Assessment when a sensor with a false report rate is supplying information. In each gure, the results are presented for the optimal values found for the sensor belief and number of e ective particles found in Chapter 4 for each motion model. Results are given without the risk assessment, with the risk assessment, and the results for a perfect sensor are also plotted for comparison. 103 0 10 20 30 40 50 60 0 0.5 1 1.5 2 2.5 3 3.5 4 x 10 5 Time (sec) Mean Square Error (m 2 ) Traffic Tuned No Risk Traffic Tuned Risk Traffic Perfect Sensor Dispersion Tuned No Risk Dispersion Tuned with Risk Dispersion Perfect Sensor Figure 4.29: E ect of Risk Assessment, Straight Path, 10% False Report Sensor Figure 4.29 illustrates the bene t to the tracking performance provided by the risk assessment in the case of the Tra c Motion Model. By not resampling after a possibly false measurement, the spatial resolution of the particle cloud is preserved, and the tracking performance greatly improves when compared to the error values when the risk assessment is not included. The Dispersion Model yields poor tracking results in the case of a straight path, regardless of the particle lter?s measurement update and inclusion of risk assessment. The frequency of believing a false detection greatly decreased for both motion models. As expected, this increased the time to detect when compared to the case of a perfect sensor and the case of when all measurements are immediately accepted. 104 Table 4.15: Time to Detect and False Start Results, Straight Path, 10% False Report Sensor Motion Risk Sensor Time to Percent False Model Assessment Model Detect Starts Tra c No 10% False Report 22:53 sec 70:50 Tra c Yes 10% False Report 70:92sec 19:50 Tra c No Perfect 28:80 sec 0 Dispersion No 10% False Report 16:47 sec 45.00 Dispersion Yes 10% False Report 33:69 sec 10:00 Dispersion No Perfect 57:90 sec 0 0 10 20 30 40 50 60 0 0.5 1 1.5 2 2.5 3 3.5 4 x 10 5 Time (sec) Mean Square Error (m 2 ) Traffic Tuned no Risk Traffic Tuned with Risk Traffic Perfect Sensor Dispersion Tuned no Risk Dispersion Tuned with Risk Dispersion Perfect Sensor Figure 4.30: E ect of Risk Assessment, Winding Path, 10% False Report Sensor Figure 4.32 illustrates the improved tracking performance that results from the imple- mentation of the risk assessment in both motion models. The Tra c Motion Model continues to provide improved tracking performance. 105 Table 4.16: Time to Detect and False Start Results, Winding Path, 10% False Report Sensor Motion Risk Sensor Time to Percent False Model Assessment Model Detect Starts Tra c No 10% False Report 19:68 sec 61:00 Tra c Yes 10% False Report 71:62 sec 21:00 Tra c No Perfect 32:20 sec 0 Dispersion No 10% False Report 17:29 sec 55.00 Dispersion Yes 10% False Report 39:19 sec 14:00 Dispersion No Perfect 67:70 sec 0 0 10 20 30 40 50 60 0 0.5 1 1.5 2 2.5 3 3.5 4 x 10 5 Time (sec) Mean Square Error (m 2 ) Traffic Tuned no Risk Traffic Tuned Risk Traffic Perfect Sensor Dispersion Tuned no Risk Dispersion Tuned Risk Dispersion Perfect Senso r Figure 4.31: E ect of Risk Assessment, Straight Path, 20% False Report Sensor The tracking performance in the case of a straight path improved with the addition of the risk assessment, and the requirement to only start observing tracking behavior after two successive reported detections. The Dispersion Model did not experience as great a decrease in the number of false starts as the Tra c Motion Model, but the tracking performance was improved. The Tra c Motion Model had a higher percentage of false starts, but the presence of the risk assessment aided in tracking performance for both motion models. 106 Table 4.17: Time to Detect and False Start Results, Straight Path, 20% False Report Sensor Motion Risk Sensor Time to Percent False Model Assessment Model Detect Starts Tra c No 10% False Report 22:53 sec 70:50 Tra c Yes 10% False Report 56:02sec 58:00 Tra c No Perfect 28:80 sec 0 Dispersion No 10% False Report 16:47 sec 45.00 Dispersion Yes 10% False Report 40:78 sec 42:00 Dispersion No Perfect 57:90 sec 0 0 10 20 30 40 50 60 0 0.5 1 1.5 2 2.5 3 3.5 4 x 10 5 Time (sec) Mean Square Error (m 2 ) Traffic Tuned no Risk Traffic Tuned Risk Traffic Perfect Sensor Dispersion Tuned no Risk Dispersion Tuned Risk Dispersion Perfect Senso r Figure 4.32: E ect of Risk Assessment, Winding Path, 20% False Report Sensor In the case of a winding path, both motion models bene t from the inclusion of the risk bene t. The Tra c Motion Model experienced a higher percentage of false starts, but was able to produce better tracking results in all cases than the Dispersion Model. In the presence of a sensor with a false report rate, the inclusion of the risk assess- ment presented above helped to improve tracking performance by preventing resampling after a possibly false measurement is received. As the false report rate increased, the risk 107 Table 4.18: Time to Detect and False Start Results, Winding Path, 20% False Report Sensor Motion Risk Sensor Time to Percent False Model Assessment Model Detect Starts Tra c No 10% False Report 19:68 sec 67:50 Tra c Yes 10% False Report 65:08 sec 61:00 Tra c No Perfect 32:20 sec 0 Dispersion No 10% False Report 17:29 sec 55.00 Dispersion Yes 10% False Report 35:19 sec 40:50 Dispersion No Perfect 67:70 sec 0 assessment?s bene t to tracking performance decreased. A false report is the result of three possibilities: two false detections reported consecutively, a false detection followed by a true detection, or a true detection followed by a false detection. The rst case is most detrimental, but is statistically less probable (1% chance in the case of a 10% false report sensor). The second and third cases are undesirable, but because a true detection is involved, the target is within the vicinity of the measured region and was likely located just on the edge of the measured region. In that instance, the required number of e ective particles and the sensor belief play critical roles in the preservation of the particle spatial resolution outside of the measured region. Further study of ltering out false detections is necessary. With the inclusion of the risk assessment, much improved tracking performance was obtained in the case of a sensor with a false report rate of 10%, and some improvement was obtained in the case of a 20% false report sensor. Only the tuned values of the number of e ective particles and the sensor belief were used in the simulations in this section because in the presence of false measurements, they provided for superior tracking performance. A longer observation period would give further insight into the convergence rate of the risk assessment method to the perfect sensor results. An additional investigation into the long term tracking performance of each particle motion model. This was done to observe the tracking performance throughout the entire simulation, not just after detection. This 108 was motivated by the order of magnitude increase in the tracking error resulting from the increasing the sensor?s false report rate by 10% increments. 0 50 100 150 200 250 300 0123456 x 10 5 Time (sec) Mean Square Error (m 2 ) 0% False Report 10% False Report 20% False Report Figure 4.33: Dispersion Model Tracking Performance over 5 minutes, Straight Path 109 0 50 100 150 200 250 300 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 x 10 5 Time (sec) Mean Square Error (m 2 ) 0% False Report 10% False Report 20% False Report Figure 4.34: Tra c Motion Model Tracking Performance over 5 minutes, Straight Path Table 4.19: Time to Actually Detect, Straight Path Motion Sensor Time to Model Model Detect Tra c Perfect 50:10 sec Tra c 10% False Report 79:65 sec Tra c 20% False Report 101:97 sec Dispersion Perfect 48:59 sec Dispersion 10% False Report 64:99 sec Dispersion 20% False Report 75:78 sec It is observed that regardless of the sensor model, the Dispersion Model yields a cyclical uctuation in tracking performance for a target on a straight path, with peaks of decreasing amplitude as time progresses. The Tra c Motion Model?s tracking performance is, as ex- pected, diminished by increased sensor inaccuracy. The 10% false report rate sensor model yields tracking performance that decreases in error as time progresses toward the error values from the perfect sensor. The 20% false report rate sensor error values decrease much more 110 slowly than the results from more accurate sensor models. The peaks and valleys in the error values occur approximately every 20 to 30 seconds. That time period is approximately how long it would take the target vehicle to traverse through a sensor region. This behavior of the error value could relate to the e ect of entering a new sensor region and potentially losing track on the target or regaining track on the target. For both motion models, the time to truly detect a target increases as the frequency of false measurements increases. 0 50 100 150 200 250 300 0123456 x 10 5 Time (sec) Mean Square Error (m 2 ) 0% False Report 10% False Report 20% False Report Figure 4.35: Dispersion Model Tracking Performance over 5 minutes, Winding Path 111 0 50 100 150 200 250 300 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 x 10 5 Time (sec) Mean Square Error (m 2 ) 0% False Report 10% False Report 20% False Report Figure 4.36: Tra c Motion Model Tracking Performance over 5 minutes, Winding Path Table 4.20: Time to Actually Detect, Winding Path Motion Sensor Time to Model Model Detect Tra c Perfect 36:57 sec Tra c 10% False Report 66:78 sec Tra c 20% False Report 116:70 sec Dispersion Perfect 61:50 sec Dispersion 10% False Report 82:29 sec Dispersion 20% False Report 80:88 sec Similar to the case of a straight path, tracking performance is diminished for both motion models as sensor accuracy decreases. The tracking performance for the Dispersion Model is less cyclical in the case of a winding path, as the Dispersion Model is more suited to track a target on a winding path. However, the tracking performance diverges for each sensor model at approximately 225 seconds after the start of the simulation. The Tra c Motion Model exhibits similar behavior in the case of both target paths, and outperforms 112 the Dispersion Model in all cases. Again, the time to truly detect the target in the presence of false measurements increases as the sensor?s accuracy decreases. The true times to detect give insight into why the false alarm rates in Sections 4.1 and 4.2 were so high. The steady decrease in tracking error yielded by the Tra c Motion Model in the case of each sensor model motivated its use in the next phase of the urban tracking problem. A UAV is to intercept the ground vehicle that was being tracked in the aforementioned studies. The UAV will utilize the information provided by the particle lter to plan its path to intercept the target. The same sensor models will be used to study the e ectiveness of the Tra c Motion Model?s tracking performance in a target intercept problem. 113 Chapter 5 Target Intercept A constant velocity UAV is to be introduced into the urban scenario. The mission of the UAV is to plan its path and intercept the target based only on the information provided to it by the particle lter. In this work, the mission will be considered complete when the UAV is within 100 m of and approaching the target, regardless of what road the UAV and target are on. The dynamic motion model for the UAV is based on the vector eld following method found in [34], and is discussed in Section 5.1. Initially, the UAV?s path will depend on only measurements provided by the soft binary regional sensor. This reinforces the importance of the particle lter?s tracking performance. The e ect of adding a measurement source to the UAV will also be investigated. The measurement model that will be added to the UAV is discussed in Section 5.2. The UAV is to plan its path to traverse above the roadways and must remain between the obstacles. The receding horizon control method for calculating the UAV?s path is discussed in Section 5.4. The method for path selection developed in this work is compared to a previously implemented minimum entropy method to show similar performance and computational cost reduction. 5.1 UAV Path Following A constant speed dynamic model for the UAV was desired for this work. The vector eld following method described herein was chosen because it was developed for constant speed MAV?s in the presence of unknown wind [34]. It has been experimentally validated in a similar urban environment. In addition, it is assumed that the UAV is aware of its 114 ground track position. Although the time step for the particle lter is 1.0 sec, a more re ned integration routine was necessary for the UAV. A fourth order Runge-Kutta integration with a time step of 0.1 sec was used to integrate the UAV equations of motion. The state vector below is used to describe the UAV, where V is the speed and is the heading angle. The speed is held constant while on a road and during a turn. _x = V sin (5.1) _y = V cos (5.2) It is assumed that the UAV is equipped with an autopilot that implements a course-hold loop and that the resulting dynamics are represented by the rst order system: _ = ( c ) (5.3) where c is the commanded heading angle and is a known positive constant that charac- terizes the response speed of the autopilot loop. Figure 5.1 illustrates the vector eld leading to a straight path with path heading angle of path = 0, as angles are measured clockwise from the positive y-axis. The time update equations to follow this path angle will be discussed below. Quadrant adjustments to the given equations are necessary for the three remaining path heading angles involved in the map in this study. 115 Figure 5.1: Vector eld for straight-line following The lateral di erence between the UAV and the desired path is denoted d. The resulting vector eld approach will direct the UAV to approach the path at angle 1 when d is very large. As d approaches zero, the heading angle will approach the path angle. The desired heading angle, d is computed as a function of the distance d by the following expression: d(d) = 12 tan 1(kd) + path (5.4) Figure 5.2 illustrates the e ect of the value of k on the di erence in the sharpness of the approach to the path heading. Large values of k result in short, abrupt transitions, whereas small values of k result in long, smooth transitions in the desired heading angle. 116 Figure 5.2: Vector elds for various values of k By restricting 1 to the range 12 (0; =2), it can be shown by using the Lyapunov function W1(d) = (1=2)d2 that y!0 asymptotically if = d(d) [34]. A sliding mode approach is used in [34] to ensure the set S = (d; ) : = d(d) is positively invariant and the system trajectory reachesSin nite time. De ning ~ d d(d) and di erentiating yields: _~ d = _ _ d(d) (5.5) = ( c ) + 1 k1 + (ky)2)V sin (5.6) (5.7) By de ning a second Lyapunov function, W2 = (1=2)~ 2 and di erentiating: _W2 = ~ _~ (5.8) = ~ ( _ _ d(d)) (5.9) = ~ ( c ) + 1 k1 + (ky)2)V sin (5.10) 117 If the control signal is chosen as: c = 1 12 k1 + (ky)2)V sin sign ( d) (5.11) where sign(x) = 8 >>> >< >>> >: 1; if x> 0 0; if x = 0 1; if x< 0 (5.12) and the gain parameter > 0 and controls the shape of the trajectories onto the sliding surface. Large values of drive ~ to zero quickly. This results in _W2 j~ j, which yields that ~ ! 0 in nite time. To avoid the chattering associated with the sign function, the control signal is chosen as: c = 1 12 k1 + (ky)2)V sin sat d (5.13) where de nes the width of the boundary region around the sliding surface in radians and: sat(x) = 8> < >: x; if jxj 1 sign(x); otherwise (5.14) The values for the constants used throughout the path following routine listed in Table 5.1 were used in experimental tests in [34]. 118 Table 5.1: Path Following Constants Constant Value V 20 m/s 1:65 rad/sec2 k 0:02/m 1 =2 rad/sec =2 rad2/sec 1:0 rad Figure 5.3 illustrates a sample path that was implemented into the urban scenario used throughout this work. The UAV?s starting point is the green square and stopping point is the red square. The UAV successfully traversed down the center of roadways, as it was commanded, and was able to complete smooth turns, and then re-center itself on a new roadway. The minimum turning radius encountered on this path was 20 m. Figure 5.3: UAV Path Following via Vector Fields Example 119 5.2 UAV Measurement Update The sensor on the UAV is modeled after the perfect binary sensor discussed in previous chapters. It is assumed that the UAV will take one measurement per second and that the eld of regard spans the width of the road the UAV is traversing for one second, creating the sensor footprint illustrated below. Figure 5.4: UAV Sensor Requirements The UAV is commanded to travel down the center of the road at a constant altitude of 30 m. The maximum width of a road in the road network is 100 m, resulting in a maximum required pan angle of approximately 60 to the right and to the left of the center of the UAV. This is illustrated by Fig. 5.5, where the dark arrow in the center of the road indicates the direction of travel. 120 Figure 5.5: UAV sensor span while traveling down-road It is assumed that the UAV travels faster than the target vehicle, therefore a constant speed of 20 m/s was chosen. In addition, to reinforce the importance of path planning, the UAV was constrained to the road network by a forced constant altitude of 30 m, which is less than the height of the obstacles. A constant commanded speed of 20 m/s, altitude of 30 m, and a measurement frequency of 1 Hz results in a maximum tilt angle of approximately 20 . Figure 5.6 illustrates the 20 angle required to cover the required ground. 121 Figure 5.6: UAV sensor tilt while traveling down-road A constant commanded altitude of 30 m results in a maximum sight distance requirement of 62 m for the UAV sensor. The UAV sensor will be assumed to have a 0% false report rate. The particle cloud will not be resampled after each UAV measurement to conserve computational expense as the UAV sensor footprint observes a relatively small subset of particles each second. This small number of particles will contribute to but not greatly a ect the number of e ective particles. This results in a sensor footprint that covers the width of the road and the length of the road traversed by the UAV each second. 5.3 Candidate Path Formation As the UAV is required to remain within the road network of the urban environment used in this study, it must make a decision when it reaches an intersection as to which way to proceed. The UAV?s path may contain up to three segments that must consist of one road between two intersections; the UAV may not traverse over an obstacle. This problem is similar to a directed graph where each intersection is a node and each road is an edge. The 122 direction of each edge is based on the UAV?s heading at the time the path is planned. The UAV may only travel forwards and is forced to stay on the map, requiring a U-turn at the map edge. The node network used in this work is illustrated by Fig. 5.7. Figure 5.7: Node Network Only a subset of the total set of nodes is used to plan the UAV?s path as it approaches an intersection. The subset of nodes are chosen based on their proximity to and relative direction from the UAV at the time of path planning. This is done to limit the complexity of the path planning problem, as the UAV replans its path at each new intersection it encounters. Planning paths to nodes far beyond UAV?s position would be a waste of computational expense, therefore a planning horizon is de ned. The path planning routine collects all permissable paths within a speci ed semi-circular region, known as the planning horizon, in the direction of the UAV?s current heading. The radius for the planning horizon used in this work is 350 m, and it is measured from the root node. The nearest intersection ahead of the UAV?s current path at the time of planning is de ned as the root node. The aforementioned path planning parameters are illustrated by Fig. 5.8. The UAV, illustrated by the green circle, is located at the point [550,400], with a heading angle of 123 3 =2 as indicated by the arrow, and is approaching an intersection located at [450,400]. In Fig. 5.8, the red circle is the target, the black circle is the root node. The blue semi-circle illustrates the UAV?s planning horizon with a radius of 350 m. The seven nodes within the planning horizon are denoted with an X. These nodes would be used to create a list of candidate paths, from which the best path for the UAV will be selected. The UAV then follows the selected path until it approaches another intersection. Figure 5.8: UAV planning horizon The remaining nodes outside of the planning horizon, denoted by blue circles, are in- cluded in Fig. 5.8 for completeness. The nodes outside of the planning horizon are not included in the formation of the candidate paths. The collection of intersections and roads within the planning horizon are used to form a directed graph, where the intersections are nodes and the roads are the edges. It is assumed that the UAV has complete knowledge of the map within its planning horizon, and is therefore aware of the location of each node and dimension and location 124 of each road within the horizon. Paths must begin at the root node and may consist of a maximum of three node points (intersections) beyond the root node to nodes within the planning horizon. The roads eligible for path creation, given the planning horizon and root node location in Fig. 5.8, are highlighted in yellow in Fig. 5.9. In the scenario given in Fig. 5.9, there are 10 possible paths to consider, as paths may consist of up to three roads and multiple routes to a node are permitted. An e ective means to collect paths between nodes in a directed grid-like map is outlined below. Figure 5.9: Node Map An algorithm was developed to build a directed graph to represent all candidate paths from the root consisting of a maximum of three roads. This framework was motivated by the use of a directed graph to represent map features such as intersections and curved roadways [29]. In general, the the path formation process begins at the root node and searches the four cardinal directions ( 2[0; =2; ;3 =2]) in order of increasing magnitude for nodes within the planning horizon. This is done by rst ranking the set of nodes in a 125 given direction within the planning horizon in order of increasing distance from the root. When the nearest node in one direction is found, the path to that node is stored and the four cardinal directions are searched from that node for additional nodes. This process continues until all paths within the planning horizon consisting of a maximum of three roads have been collected. First, the direction = 0 is searched from the root for a node within the planning horizon. If a node is found in that direction, the following quantities are stored to the characterize the path from the root to the node: which road is between the root and the node, the distance between the root and the node, the direction of travel from the root to the node, and which node was found. From the node that was found, the direction = 0 is searched for a node within the planning horizon. If a node is found in the search direction, the path to that node is stored, and a third node beyond the root is searched for in the direction = 0. If at any point during this process a node is not found in the search direction, the search direction proceeds in order of increasing values of until a node is found. If no node is found in any direction, the routine returns to the node it searched to get to the current node, and searches the remaining cardinal directions. When each direction has been searched from the root node, and consequently all nodes within the planning horizon have been searched for connections in each direction, all of the candidate paths have been collected. Each path is characterized by which roads and intersections it involves, the distance between nodes, and the direction each road would be traversed. The example below will be used to describe this process in more detail. 126 Figure 5.10: Example: Nodes and Roads for Candidate Paths Consider the road network in Fig. 5.10, where roads and intersections are drawn in white, obstacles in orange, node ID numbers are circled and road ID numbers are not circled. The search angle is measured clockwise from the positive y-axis. In this example, the node network in Fig. 5.10 is made up of only the nodes within the planning horizon. Node 4 is the only node located at the edge of the map, and therefore a U-turn would be required to depart from Node 4. The root node is denoted with an X. The intersections and roads within the planning horizon have been collected and will be used to determined the candidate paths. In this example, it is assumed that the UAV velocity is constant and each road is of equal dimension, therefore the distance between nodes connected by only one road is equal. This results in an equal travel time (T) for the UAV between any two nodes connected by only one road. Beginning at the root node, the nearest node is searched for in in the direction of = 0. Node 2 is found, and the path from the root to Node 2 is stored as detailed by 127 Table 5.2, and the rst edge of the directed graph is formed. Figure 5.11 illustrates the root node as Node 1, the road ID number as the number 1 next to the edge connected the two nodes, the direction of travel indicated by the arrow, and the end node ID as Node 2. Table 5.2: Path to Node 2 Path Route Edge Road Travel Direction End to Node ID Number ID Time of Travel Node ID 2 1 1 1 T 0 2 Figure 5.11: Edge of Directed Graph from Root Node to Node 2 After storing the rst path to Node 2, a node in the planning horizon is searched for from Node 2 in the direction of = 0. Node 5 is found and a path to Node 5 is stored, consisting of the path from the root to Node 2 and the recently found road between Nodes 2 and 5, as in Table 5.3. Another edge is added to the directed graph, as illustrated by Fig. 5.12. 128 Table 5.3: Path to Node 5 Path Route Edge Road Travel Direction End to Node ID Number ID Time of Travel Node ID 5 1 1 1 T 0 2 2 4 T 0 5 Figure 5.12: Edges of Directed Graph from Root Node to Node 5 From Node 5, the direction = 0 is searched for a node. Because there is no node in that direction, the next direction ( =2) is searched, and Node 6 is found. The path to Node 6 is stored consisting of the three nodes from the root: nodes 2, 5, and 6. Because the path to Node 6 consists of the maximum number of roads (3), no further nodes will be searched for from Node 6. The directed graph is updated to include the recently found path to Node 6. 129 Table 5.4: Path to Node 6 Path Route Edge Road Travel Direction End to Node ID Number ID Time of Travel Node ID 6 1 1 1 T 0 2 2 4 T 0 5 3 5 T =2 6 Figure 5.13: Edges of Directed Graph from Root Node to Node 6 As no more nodes may be searched for from Node 6, the routine will return to the node it found immediately prior to Node 6, which is Node 5. The next direction to search from Node 5 is = . That would require a U-turn, which is not permitted when not at the map?s edge. The routine then searches the nal direction from Node 5, = 3 =2, and there are no nodes in that direction. Therefore, the routine returns to the node it found prior to Node 5, which is Node 2. As the direction = 0 has already been searched from Node 2, the routine searches the next direction for Node 2, = =2. Node 3 is found, and the path to Node 3 is stored in the path list and the directed graph is updated. Next, the direction 130 of = 0 is searched from Node 3 and Node 6 is found, resulting in a second route to Node 6. The path list is updated accordingly. Table 5.5: Path to Node 3 Path Route Edge Road Travel Direction End to Node ID Number ID Time of Travel Node ID 3 1 1 1 T 0 2 2 2 T =2 3 Figure 5.14: Directed Graph after nding Node 3 131 Table 5.6: Paths to Node 6 Path Route Edge Road Travel Direction End to Node ID Number ID Time of Travel Node ID 6 1 1 1 T 0 2 2 4 T 0 5 3 5 T =2 6 6 2 1 1 T 0 2 2 2 T =2 3 3 6 T 0 6 Figure 5.15: Directed Graph after nding a second route to Node 6 Again, no additional nodes are searched for from Node 6 because the paths to Node 6 consist of the maximum number of edges. Therefore, the routine returns to the previous node, Node 3. As there are no nodes in the directions =2 or from Node 3, and Node 2 would require an illegal U-turn, the routine retraces its path to Node 3 to Node 2. Searching for a node in the direction of from Node 2 would require an illegal U-turn, therefore the 132 next direction is searched. The direction 3 =2 is searched from Node 2 and Node 4 is found and the path to Node 4 is stored and the directed graph is updated. Table 5.7: Path List Path Route Edge Road Travel Direction End to Node ID Number ID Time of Travel Node ID 4 1 1 1 T 0 2 2 3 T 3 =2 4 Figure 5.16: Directed Graph after nding Node 4 Node 4 is at the edge of the map, therefore a U-turn will be permitted, resulting in another path to Node 2. As all nodes within the planning horizon have been searched in each direction, the candidate path list is complete and is given in Table 5.8. There are a total of seven candidate paths in the nal directed graph illustrated by Fig 5.17. 133 Table 5.8: Final Path List Path Route Edge Road Travel Direction End to Node ID Number ID Time of Travel Node ID 2 1 1 1 T 0 2 2 2 1 1 T 0 2 2 3 T 3 =2 5 3 3 T =2 6 3 1 1 1 T 0 2 2 2 T =2 3 4 1 1 1 T 0 2 2 3 T 3 =2 4 5 1 1 1 T 0 2 2 4 T 0 5 6 1 1 1 T 0 2 2 4 T 0 5 3 5 T =2 6 6 2 1 1 T 0 2 2 2 T =2 3 3 6 T 0 6 Figure 5.17: Directed Graph of Example Path Planning Problem Following the collection of the candidate paths, a means of selecting the path with the most bene t was developed. Section 5.4 outlines the maximum likelihood method for 134 selecting the best candidate path that was developed in this work, and compares that method to an entropy reduction method presented in previous work. 5.4 Path Selection Once the candidate paths have been collected, the best path for the UAV must be selected. It is desirable to choose a path that will enable the UAV to intercept the target as quickly as possible based only on the information provided by the particle lter. Numerous metrics exist for path selection such as minimal mean square error, maximal Information gain, minimal Shannon entropy, and maximal likelihood of detection. In Chapters 3 and 4, the large regional sensor measured the region of highest particle weight. The same metric will be used to select the best candidate path. In this work, it will be shown that selecting the path that holds the highest particle weight at the time the UAV will traverse it provides for e ective target intercept. In addition, comparable performance to a entropy reduction method is achieved with a signi cant reduction in computational time. The time to intercept the target and the amount of time to simulate one second will be used as metrics to compare the two methods presented herein. In the path planning comparison below, a path may consist of a maximum of three roads. Both methods will make their path selections based on similar cost functions in an e ort to maximize the amount of information gained, minimize the distance between the UAV and the expected value of the particle cloud, and minimize the travel time of the UAV. Jpath = NedgesX i=1 1 Path Valuei +di=Nedges +tpath;i (5.15) The second two terms in the summation above will be the same for both methods. In the maximum likelihood method, the rst term will involve the sum of the particle weights to be encountered on a path. In the minimum entropy method, the rst term will involve a mutual information utility function. 135 Computational time is taken into account when the UAV approaches an intersection to ensure enough time is available to select the next path and have a decision made when the UAV arrives at an intersection. Although path are planned to have a length of up to three roads, corresponding to upwards of 20 seconds in the future, paths are be replanned at each intersection to account for additional measurement data from outside sources, as consistent with receding horizon control. Therefore, four seconds before reaching the next intersection, a new path is planned. 5.4.1 Maximum Likelihood In previous work, path selection was based on maximizing the reduction in the entropy of the pdf [24], [26], [28]. In this work, path bene t is based on the weight of the particles along a path at the appropriate time, the distance between the nal node on a path and the expected value of the particle cloud, and the amount of travel time required by each path. Jpath = NedgesX i=1 R Nwedgei + d i rplan Q + tpath2 (5.16) The cost function to be minimized is de ned by Eq. (5.16), where Nedges is the number of edges, or roads, in a path, N is total number of particles, di is the distance between the nal node in the path and the expected value of the particle cloud, rplan is the planning radius for the UAV, and tpath is the amount of time it would take the UAV to complete the path. The parameters R and Q may be used to adjust the impact of the weight term and the distance term on the cost of a path. In this work, R is set to unity and Q is 2. A path is chosen through a simple search for a minimum cost using the min function in MATLAB. The particle weight of a road was calculated by propagating the particle cloud forward in time without measurements from outside sources (human operatives). The weight of a particles on a road is measured at two times: at the midpoint of the road and at the end node. At each of those times, the region of the road the UAV has just passed is measured for particle presence. The time of interest to which the particle cloud is propagated is determined by 136 the speed of the UAV and the length of the road, which is stored in the candidate path list. Figure 5.18 illustrates the two regions of two roads on a candidate path. The blue regions would be measured at the midway point of the roads, times 5 and 15, respectively. The orange regions highlight the second half of the road and the end node, which are measured for particle weight at times 10 and 20. The road the a particle is on is stored in its state vector. This reduces computational time by avoiding searching each road for N particles at each time of interest. Figure 5.18: Two Part Road Weight If there are no particles on a road, the road weight is set to 1=N to avoid dividing by zero in Eq. (5.16). The particle weight on the road of interest is then calculated to determine target likelihood at the time when the mobile sensor will be on that road. This is a computationally expensive process that is also necessary in minimum entropy methods. 5.4.2 Information Utility Function A minimum entropy path selection routine as described in [25] was implemented to provide for comparison to the maximum likelihood approach described above. As mentioned 137 in Section 2.3, the Information Utility Function (IUF) is computed by: I (zk; xk) = H (zk) H (p(zkjxk)) The IUF is based on the idea that minimizing the posterior uncertainty is equivalent to maximizing the di erence between the uncertainty that any particular observation will be made H (zk), and the uncertainty of the measurement model, H (p(zkjxk)). The entropy terms are approximated by the following expressions: H (zk) Z ( NX i=1 wi kp z kjxik log 2 NX i=1 wi kp z kjxik )dz H (zkjxk) Z NX i=1 wi kp z kjxik log 2p z kjxik dz Because in this work, only one measurement is received per time step, the integral falls away, reducing the IUF to: I (zk; xk) NX i=1 wi kp z kjxik log 2 NX i=1 wi kp z kjxik + NX i=1 wi kp z kjxik log 2p z kjxik The weights of the particles wik;i = 1;:::;N, are readily available. The measurement likelihood p(zkjxik) is computed based on the UAV sensor model and the ith particle?s po- sition, as discussed in Section 3.4. A false report belief of 1% was used to avoid taking the logarithm of zero. 5.5 Method Comparison Results Results for the maximum likelihood (ML) and Information Utility Function (IUF) sensor routing methods are presented below. The target was initialized at a random location and 138 the UAV was initialized at a constant location to begin each run, and the target traversed a random path throughout each run. For a metric of comparison, an evaluation was done of the expected amount of time it would take the UAV to travel between two random locations in the map (See Appendix A). The expected value of the UAV?s travel time from a random point in the environment to a random stationary point is 33.33 seconds. An average time to detect for the perfect regional sensor for a target on a random path starting at a random location was found to be 47.31 seconds. This time was added to the expected value of travel time between two points to be used as a metric of comparison for target interception time. Therefore, if the target were to be stationary, the expected intercept time should be approximately 80 seconds. The Tra c Motion Model was used to propagate the particles in time and 1000 particles were used. The risk assessment was included and the sensor belief and number of e ective particles were tuned to the optimized results obtained in Chapter 4. The results in Table 5.9 indicate that the ML and IUF routines yield similar intercept times, while the IUF routine outperforms the ML routine slightly in each case. After observing the similarity in the intercept times, it can be concluded that there is not a signi cant di erence in the paths chosen by the maximum likelihood routine and the IUF routine. This also implies that the ML routine indirectly reduces the entropy of the particle cloud as the IUF routine does, but at a much reduced computational expense. This is a desirable result as the path planning is to be done on board a small UAV, where computational power is limited. The results in Table 5.9 indicated that in the case of a perfect sensor, the UAV inter- cepted the target at approximately 81 seconds in both cases, which slightly greater than the expected value of the time to intercept a stationary target of 80 seconds. This shows the importance of path planning and the e ectiveness of the chosen methods. Additionally, the detection time in the case of a perfect sensor is less than the average value of 47.31 seconds that was obtained without the presence of a UAV. This is because although the UAV did not 139 take measurements, it could intercept the target on its own before the target was located by the large regional sensor. Table 5.9: Time to Intercept and Run Time Results, without UAV measurements Sensor Path Selection Time to % False Average Model Routine Intercept Alarms Detection Time Perfect ML 81:28 sec 0 37:62 sec Perfect IUF 80:44 sec 0 32:45 sec 10% False Report ML 111:35 sec 12:50 45:16 sec 10% False Report IUF 123:13 sec 15:50 43:44 sec 20% False Report ML 150:15 sec 40:50 36:67 sec 20% False Report IUF 147:75 sec 44:50 39:91 sec When a sensor with a false report rate is utilized, the time to intercept increased for both the ML and IUF routines beyond the 80 second mark to intercept a stationary target. This delayed interception time is due to false reports incorrectly shifting particle weight, despite the e orts of the risk assessment, throughout the path planning process. The ML routine outperforms the IUF method in the case of a 10% false report rate sensor in regard to intercept time. This reinforces the conclusion that the ML routine is an e ective means of path planning that requires reduced computational expense when compared to the IUF routine. In addition, the false alarm rates for both imperfect sensor models are much lower than those presented in Section 4.3. This occurs because the UAV is capable of intercepting the target before the large regional sensor. Table 5.10 lists the results given the use of the UAV as a measurement source. As expected, the inclusion of the measurements from the UAV decreased the time to intercept. The most signi cant decrease in interception time occurred in the case of the 10% false report rate sensor. In addition, when a perfect regional sensor is used, the average detection time decreased in the case of both ML and IUF routines when compared to the average value of 47.31 seconds achieved without the UAV. The reduction in the time to intercept the target provided by UAV measurements was less pronounced in the case of a perfect sensor and a 140 Table 5.10: Time to Intercept and Run Time Results, with UAV measurements Sensor Path Selection Time to % False Average Model Routine Intercept Alarms Detection Time Perfect ML 77:65 sec 0 34:21 sec Perfect IUF 81:53 sec 0 35:01 sec 10% False Report ML 105:98 sec 17:26 43:30 sec 10% False Report IUF 102:52 sec 14:00 39:88 sec 20% False Report ML 149:90 sec 43:50 36:08 sec 20% False Report IUF 146:42 sec 43:00 31:50 sec sensor with a 20% false report rate. This illustrates the strong e ect of the large regional sensor on the particle cloud. The regional sensor is able to observe large areas (40,000 m2) every 3 seconds, and there is complete coverage of the simulation space at all times. The UAV is restricted to observe the area within its vicinity, with maximum coverage of 6,000 m2 over 3 seconds. To further illustrate this, the simulation was run with only measurements from the UAV and no regional sensor. The particle distribution was resampled every 3 seconds and a measurement was taken by the UAV every second. Table 5.11: Comparison of ML and IUF routines using UAV measurements without regional sensor data Time to Routine Intercept Maximum Likelihood 225:70 sec Information 201:30 sec The results in Table 5.11 illustrate the importance of the large regional sensor, and there- fore the inclusion of non-traditional measurement sources into the target intercept problem. In addition, the presence of the UAV is vital to the successful localization and tracking in the presence of a regional sensor with a false report rate. The results in Chapter 4 indicate that given the regional sensor alone, poor tracking results are obtained in the presence of false reports. This is often the result of accepting a false report to be true, and therefore possibly 141 losing the target completely. The inclusion of the UAV and an intelligent path planning rou- tine provides for consistent target localization. Further bene t to intercept time reduction would be obtained with the inclusion of additional UAV?s. The computational expense for each scenario is listed in Table 5.12. The run time is decreased by use of the ML method both with and without UAV measurements. In addition, the use of UAV measurements does not add signi cant computational expense to either routine. Table 5.12: Time to simulate one second Path Selection Without UAV With UAV Routine Measurements Measurements ML 0.3671 sec 0.3839 sec IUF 0.4358 sec 0.4518 sec Finally, the ML and IUF methods were compared by just using the large regional sensor to locate the target. No UAV was simulated and no paths were planned. A comparison was done to determine which region to poll based o of either the sum of the particle weights, which is comparable to the maximum likelihood method, or which region would provide the highest value of the Information Utility function. The work presented in previous chapters polled the sensor in the region with the highest particle weight sum. The Information Utility Function approach was implemented in the regional sensor framework for comparison. The times to detect the target were compared. A perfect sensor was used and only the Tra c Motion Model with 1000 particles was used in this study. The results give in Table 5.13 indi- Table 5.13: Comparison of ML and IUF routine using only large regional sensor data Time to Routine Detect Maximum Likelihood 47:31 sec Information 49:96 sec 142 cate that the maximum likelihood routine yields similar results to the maximum information routine, with a signi cantly decreased computational expense. The results presented above illustrate that the particle weight to be encountered along a path is a su cient and computationally e cient metric for path ranking. The particle lter framework provides for this intuitive metric without the need to perform additional calculations. Proper tuning of particle lter parameters such as sensor con dence and the number of e ective particles provide for improved target tracking, and is therefore necessary when particle weight is to be used as the path selection metric. 143 Chapter 6 Conclusion The development of a particle lter based method of localizing, tracking, and inter- cepting a mobile ground target in an urban scenario given nontraditional and inconsistent measurement data was presented. The large regional sensor model was a signi cant chal- lenge in this work. The target?s location could only be measured to be within the region of the detection, and the time span between measurement reports reinforced the importance of a sophisticated dynamic model. In the presence of a perfect sensor, the dynamic model used to propagate the particle cloud in time plays a vital role in tracking performance. In addition, the importance of particle spatial resolution must be taken in to consideration to improve tracking performance and reduce detection time in the presence of the aforemen- tioned measurement model. In the presence of a sensor with a false report rate, the measurement update becomes more important than the dynamic model update. The particle lter?s belief in the sensor?s false report rate, coupled with the frequency at which the distribution is resampled must be taken in to account in order to maintain proper spatial resolution of the particle cloud in the presence of a sensor with a false report rate. Additionally, means must be taken to assess when a measurement is possibly false. The ux of particle weight into or out of a region may be used to assign a cost to accepting a measurement through the Bayesian risk assessment framework. A risk assessment is necessary to provide proper tracking in the presence of false measurements. A risk assessment and proper tuning of the number of e ective particles and sensor belief proved bene cial in maintaining proper spatial resolution in the presence of a sensor with a false report rate. 144 These conclusions proved vital in the navigation of a UAV to the target by providing motivation for a new path selection metric. A metric for path selection in the particle lter framework was developed based on the cumulative sum of particle weight to be encountered on a path. This metric provides insight to future likelihood of target intercept while taking advantage of the particle lter framework. A comparison was done to an existing method where the new metric matched performance and required reduced computational expense. This work will bene t the modern war ghter, in addition to researchers considering a particle lter due to a non-di erentiable measurement model. Additionally, it is observed that given set computational capabilities, a tradeo may be made to give priority to the dynamic model or the path selection metric. The amount of time to simulate one second with the Dispersion model and the Information Utility Function is approximately equal to that of the use of the Tra c Motion Model and the Maximum Likelihood method. Signi cantly improved tracking performance results from using the Tra c Motion Model, while target intercept time is not as signi cantly delayed by use of the Maximum Likelihood method. By developing a sophisticated dynamic model, a less computationally intensive path selection metric is applicable and successful. Further improvement to the work discussed herein would be achieved by a more robust risk assessment routine to reject false measurements. Currently, if the risk is deemed unac- ceptable, the measurement is still accepted but the particle distribution is not resampled. Investigation into the practice of rejecting measurements given the current framework could provide for improved tracking performance and false alarm rate reduction. In addition, further study into the e ect of various values of the coe cients R and Q in the cost function could provide decreased intercept times. In the presented results, the distance term in the cost function is the distance between the UAV and the expected value of the particle cloud. This is an e ective metric after the target has been found by the regional binary sensor, but ine ective before that time of detection because the expected value tends to be in the center of the city until the target is detected. An improvement to this would be 145 to adjust Q to include the distance term only after the target is found, or to use the distance between the UAV and the center of the region of highest weight until the target is found. Future expansion of this problem could include the idea of particle interaction, which may provide for improved tracking performance. Particles are statistically required to be uncorrelated, and therefore may not necessarily be aware of what other particles are choosing to do. However, real time tra c updates are often readily available, and would t nicely into the presented framework to allow for adaptive particle motion models, without violating the statistical properties of the particle lter. Finally, mobile sensor units could be introduced into the existing framework with min- imal e ort. This could provide more a more realistic sensor network structure. In addition, the present framework permits the acceptance of multiple measurement reports at once or at varying intervals. This would provide for the inclusion of measurements from additional sources such as security cameras and satellite images. 146 Bibliography [1] Ristic, B., Arulampalam, S., and Gordon, N., Beyond the Kalman Filter: Particle Filters for Tracking Applications, Artech House, Boston, 2004. [2] Crassidis, J.L. and Junkins, J.L., Optimal Estimation of Dynamic Systems, Chapman & Hall/CRC, Boca Raton, Florida, 2004. [3] Gordon, N.J., Salmond, D.J., and Smith, A.M.F, \Novel approach to nonlinear/non- Gaussian Bayesian state estimation," IEE Proceedings-F, vol. 140, no.2, pp. 107-113, 1993. [4] Doucet, A., Godsill, S., and Gordon, N., \An introduction to sequential Monte Carlo Methods," in Sequential Monte Carlo Methods in Practice, Springer, New York, 2001. [5] Boers, Y., Driessen, J.N., \Particle lter based detection for tracking," Proceedings of the American Control Conference, vol. 6, pp. 43934397, 2001. [6] Kalman, R.E., \A new approach to linear ltering and prediction problems," Transac- tions of the ASME Journal of Basic Engineering, (82 (Series D)), 3545, 1960. [7] Kalman, R.E., Bucy, R.S., \New results in linear ltering prediction theory," Transac- tions of the ASME Journal of Basic Engineering, 83(1), 95108, 1961. [8] Liu, J.S., Monte Carlo Strategies in Scienti c Computing, Springer, Verlag, 2001. [9] Agate, C.S., Wilkerson, R.M., Sullivan, K.J., \Utilizing negative information to track ground vehicles through move-stop-move cycles," Proceedings of SPIE on Signal Pro- cessing, Sensor Fusion, and Target Recognition, vol. 5429, pp. 273283, 2004. [10] Arulampalam, M.S., Gordon, N., Orton, M., Ristic, B., \A variable structure multiple model particle lter for GMTI tracking," Proceedings of the International Conference on Information Fusion, vol. 2, pp. 927934, 2002. [11] Agate, C.S., Sullivan, K.J., \Road-constrained target tracking and identi cation using a particle lter," Proceedings of the SPIE Conference on Signal and Data Processing of Small Targets, vol. 5204, pp. 532543, 2003. [12] Yang, C., Bakich, M., Blasch, E., \Nonlinear constrained tracking of targets on roads," Proceedings of the International Conference on Information Fusion, pp. 235242, 2005. [13] Ulmke, M., Koch, W., \Road-map assisted ground moving target tracking," IEEE Transactions on Aerospace and Electronic Systems, 42(4), 12641274, 2006. 147 [14] Ekman, M., Sviestins, E., \Multiple model algorithm based on particle lters for ground target tracking" In: Proceedings of the International Conference on Information Fusion, pp. 18, 2007. [15] Kamrani, F., Ayani, R. \Simulation-aided path planning of UAV" In: Proceedings of the 2007 Winter Simulation Conference, pp. 13061314 (2007) [16] Orguner, U., Schon, T.B., Gustafsson, F., \Improved target tracking with road network information" In: Proceedings of the IEEE Aerospace Conference, pp. 111 (2009) [17] Hong, L., Cui, N., Bakich, M., Layne, J.R.: Multirate interacting multiple model parti- cle lter for terrain-based ground target tracking. IEE Proceedings on Control Theory and Applications 153(6), 721731 (2006) [18] Kravaritis, G., Mulgrew, B.: Variable-mass particle lter for road-constrained vehicle tracking. EURASIP Journal on Advances in Signal Processing 2008, 113 (2008) [19] Skoglar, P., Orguner, U., Tornqvist, D., Gustafsson, F.: Road target tracking with an approximative Rao-Blackwellized particle lter. In: Proceedings of the International Conference on Information Fusion, pp. 1724 (2009) [20] F. Zhao, J. Shin, and J. Reich, \Information-driven dynamic sensor collaboration, IEEE Signal Processing Mag., pp. 6172, March 2002. [21] N. E. Leonard, D. A. Paley, F. Lekien, R. Sepulchre, D. M. Fratantoni, and R. E. David, \Collective motion, sensor networks, and ocean sampling, Proc. IEEE, vol. 95, no. 1, pp. 4874, January 2007. [22] X. Feng, K. A. Loparo, and Y. Fang, \Optimal state estimation for stochastic systems: An information theoretic approach, IEEE Trans. Automat. Contr., vol. 42, no. 6, June 1997. [23] C. Kreucher, K. Kastella, and A. O. Hero III, \Information based sensor management for multitarget tracking, in Proc. SPIE, vol. 5204, Bellingham, WA, August 2003, pp. 480489. [24] G. M. Ho mann, S. L. Waslander, and C. J. Tomlin, \Mutual information methods with particle lters for mobile sensor network control, in Proc. 45th IEEE Conf. Decision and Control, San Diego, CA, December 2006, pp. 10191024. [25] Ho mann, G.M., Tomlin, C.J., \Mobile Sensor Network Control Using Mutual Infor- mation Methods and Particle Filters", IEEE Transactions on Automatic Control, vol. 5, Issue 1, p. 32-47, 2010. [26] Skoglar, P.; Orguner, U.; Gustafsson, F. , \On information measures based on particle mixture for optimal bearings-only tracking," Aerospace conference, 2009 IEEE , vol., no., pp.1-14, 7-14 March 2009. 148 [27] Boers) Boers, Y., Driessen, H., Bagchi, A., and Mandal, P., \Particle Filter Based Entropy", 13th Conference on Information Fusion (FUSION), 2010, vol., no., pp.1-8, 26-29 July 2010. [28] Ryan, A., \Information-theoretic tracking control based on particle lter estimate," in Proceedings AIAA Guidance, Navigation, and Control Conf., Honolulu, HI, August, 2008. [29] Berg-Yuen, P.E.K, Mehta, S.S., Waden, D.L., Curtis, J.W., \Particle Filter-Based Ground Target Tracking in a Constrained Road Network", 3rd International Confer- ence on the Dynamics of Information Systems, 2011. [30] LaValle, S., Planning Algorithms, Cambridge University Press, 2006. [31] Spletzer, J.R. and Taylor, C.J., \Dynamic Sensor Planning and Control for Optimally Tracking Targets", The International Journal of Robotics Research, No. 1, Jan. 2003, p. 7-20. [32] Wang, Y., Hussein, I.I., and Erwin, R.S., \Risk Based Sensor Management for Integrated Detection and Estimation", AIAA Journal of Guidance, Control, and Dynamics,2011, (in press). [33] O?Reilly, P.G., \Local Sensor Fault Detection Using Bayesian Decision Theory", Con- ference Proceedings: UKACC International Conference on Control, Sept. 1998, Vol. 1, p. 247-251. [34] Nelson, D.R.; Barber, D.B.; McLain, T.W.; Beard, R.W.; , \Vector Field Path Following for Miniature Air Vehicles," Robotics, IEEE Transactions on , vol.23, no.3, pp.519-529, June 2007. [35] Grinstead, C.M. and Snell, J.L., Introduction to Probability, American Mathematical Society, 1998. 149 Appendices 150 Appendix A Expected Travel Time A performance metric is desired for the estimated travel time between two random points in a 1000 m x 1000 m urban environment. A UAV approaching a stationary target at a constant speed is considered. Assuming travel is required in a grid-like fashion, travel is only permitted in one direction at at time. This metric could be determined by dividing the expected value of distance traveled between two random points in the environment by the constant speed of travel. The following derivation of the expected value of the distance between two random points is used in the comparison of the path planning metrics. The expected value of a random variable is given by: E(X) = Z 1 1 xf(x)dx (A.1) Equation (A.1) is to be used to determined the expected value of the distance between to random points in a uniform distribution. Therefore, an expression for the probability density function for the distance between two random points is required. Consider two independent random numbers from the interval [0;d], X and Y, where X represents the x-coordinate of the UAV starting position and Y represents the target?s x-coordinate. The probability density functions of X and Y are given by: fX(x) = fY (x) = 8 >>< >>: 1 if 0 x d 0 if otherwise (A.2) The di erence between X and Y is de ned as Z = X Y. Because X and Y are independent, the density of their sum is the convolution of their densities [35]. Therefore, the density 151 function of Z = X Y is: fZ(z) = Z 1 1 fX(z +y)fY (y)dy (A.3) Because fY (y) = 1 if 0 y d and 0 otherwise, fZ(z) reduces to fZ(z) = Z d 0 fX(z +y)dy (A.4) Now, the integrand is 0 unless 0 z+y d (i.e., unless z y d z), where it is 1. So, if 0 z d, fZ(z) = Z 1 z 0 dy = d z (A.5) and if d z 0, fZ(z) = Z d z dy = d+z (A.6) and fZ(z) = 0 otherwise. The previous example was done in terms of a uniform distribution between 0 and an unknown constant d that corresponds to the dimension of the space. Because distance is always positive, the expected value of the absolute value of z is taken. E(jzj) = Z 1 1 jzjfZ(z)dz = Z 0 d z(d+z)dz + Z d 0 z(d z)dz = d=3 The value of d in this study is 1000 m, therefore the expected value of the distance in the x-direction between two points is x = 1000=3 = 333:33 m. The same would result in the y-direction. The grid-like structure of the road network in this work requires travel in one direction and then another. This results in a total expected value for the traveled distance 152 between two points in the given 1000 m x 1000 m x-y plane ( ) as: = x + y = 333:33 m + 333:33 m = 666:66 m Given a constant speed of 20 m/s and assuming the target remained stationary, this would result in an expected travel time of 33:33 sec. 153