UFDC Home  UF Institutional Repository  UF Theses & Dissertations  Internet Archive   Help 
Material Information
Subjects
Notes
Record Information

Table of Contents 
Title Page
Page i Page ii Dedication Page iii Acknowledgement Page iv Page v Table of Contents Page vi Page vii List of Tables Page viii List of Figures Page ix Page x Page xi Page xii Page xiii Page xiv Page xv Abstract Page xvi Page xvii Chapter 1. Introduction Page 1 Page 2 Page 3 Page 4 Chapter 2. Background Page 5 Page 6 Page 7 Page 8 Page 9 Page 10 Page 11 Page 12 Page 13 Page 14 Page 15 Page 16 Page 17 Chapter 3. An antijam attitude determination receiver Page 18 Page 19 Page 20 Page 21 Chapter 4. Maximum likelihood attitude estimation Page 22 Page 23 Page 24 Page 25 Page 26 Page 27 Page 28 Page 29 Page 30 Page 31 Page 32 Page 33 Page 34 Page 35 Page 36 Page 37 Page 38 Page 39 Page 40 Page 41 Chapter 5. Properties of the maximum likelihood attitude estimator Page 42 Page 43 Page 44 Page 45 Page 46 Page 47 Page 48 Page 49 Page 50 Page 51 Page 52 Page 53 Chapter 6. Attitude from direction finding Page 54 Page 55 Page 56 Page 57 Page 58 Page 59 Page 60 Chapter 7. Wider baselines and dual frequency use Page 61 Page 62 Page 63 Page 64 Page 65 Page 66 Page 67 Page 68 Page 69 Page 70 Page 71 Page 72 Page 73 Page 74 Page 75 Page 76 Page 77 Page 78 Page 79 Page 80 Page 81 Page 82 Page 83 Chapter 8. Simulations and results Page 84 Page 85 Page 86 Page 87 Page 88 Page 89 Page 90 Page 91 Page 92 Page 93 Page 94 Page 95 Page 96 Page 97 Page 98 Page 99 Page 100 Page 101 Page 102 Page 103 Page 104 Page 105 Page 106 Page 107 Page 108 Page 109 Page 110 Page 111 Page 112 Page 113 Page 114 Page 115 Page 116 Page 117 Page 118 Page 119 Page 120 Page 121 Page 122 Page 123 Page 124 Page 125 Page 126 Page 127 Chapter 9. Conclusions Page 128 Page 129 Page 130 Page 131 Page 132 Appendix A. Cramerrao bound Page 133 Page 134 Page 135 Page 136 Page 137 Page 138 Page 139 Page 140 Appendix B. The interference covariance matrix Page 141 Page 142 Page 143 Page 144 Page 145 Page 146 Page 147 Page 148 Page 149 Page 150 Page 151 Page 152 Page 153 Page 154 Page 155 Page 156 Page 157 Page 158 Page 159 Appendix C. Resolution of the attitude ambiguity Page 160 Page 161 Page 162 Page 163 References Page 164 Page 165 Page 166 Page 167 Biographical sketch Page 168 Page 169 Page 170 Page 171 
Full Text 
INTERFERENCE MITIGATION FOR GPS BASED ATTITUDE DETERMINATION B y .. ..... ....... MATTHEW DAVID MARKEL A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2002 Copyright 2002 by Matthew David Markel This work is dedicated to my wife and best friend, Colleen. Her unending support through every facet of this program was truly an inspiration. On an uncountable number of occasions she has been a source of both encouragement and enlightenment. ACKNOWLEDGMENTS I express my great appreciation to my parents, David and Linda Markel, for their encouragement and support. For as long as I can remember, they have encouraged me to make the most of the gifts I have been given. My success is a direct consequence of their encouragement. It goes without saying that I am in the debt of my wife, Colleen, to whom this work is dedicated. She has graciously given of herself to allow me the means to pursue this program. On more days than I can count she, without complaint, has sacrificed so that I could complete this program. She was understanding when I was not present because of school and when I was physically present although elsewhere in thought. I will spend the rest of my life repaying her for what she has without question provided me these last three years. I have been extremely fortunate to have had Professors Henry Zmuda and Eric Sutton as my coadvisors. Their insight, encouragement, and talent have been a genuine blessing. The attitude and dedication they have shown throughout this work have made this research an exciting and fulfilling endeavor for me. I will truly miss our weekly research meetings. I am grateful as well for the other members of my research committee, Professors Pasquale Sforza and Tan Wong. In addition to serving on the my research committee, I was fortunate to take a class taught by Professor Wong, where several key concepts used in this dissertation were presented. Professor Sforza, as director of the Graduate Engineering Research Center (GERC), has succeeded in the difficult task of providing a center that provides a university setting beneficial to the unique challenges of a military base environment. I am also thankful for the generosity of my dissertation reviewers: Dr. Timothy J. (TJ) Klausutis, Mr. James Ciccarelli, and Mr. Jerry Weed. Each provided a wealth of constructive comments and critiques, and their contributions without doubt have made this a stronger dissertation. I have been fortunate to take classes taught by the talented and dedicated Dr. Klausutis. More than the academic material, however, Dr. Klausutis has been a source of much needed information and inspiration concerning the Ph.D. process and experience. Mr. Ciccarelli, in addition to reviewing the dissertation, has contributed to our research team since June 2001, providing both a sounding board for the technical details and an independent review of my simulation code. Mr. Weed was the first to encourage me to pursue this Ph.D. For years he has been the instigator of stimulating technical discussion between us concerning nearly every aspect of radar and signal processing, forcing us each to a higher level of understanding. The logistics of taking classes, receiving and presenting assignments and tests, registration, and general communication with the Department and Graduate School is often taken for granted by students on campus. However, for the offcampus students at the GERC, it would be a significant challenge were it not for the professionalism of the GERC staff, especially Ms. Judi Shivers. I am extremely grateful for the help each has provided me during both my master's and doctoral work. Sverdrup Technology, through its Edwin "Bud" George Fellowship, has financed my academic pursuit. Without their support and flexibility I would not have been successful. This fellowship is a benefit to the company, its employees, and the symbiotic relationship among Sverdrup, Eglin Air Force Base, and the University of Florida. It is my sincere hope that it endures and continues to attract students. Finally, I would like to thank God for the talents and blessings He has provided. I truly believe that with Him all things are possible. TABLE OF CONTENTS page ACKNOWLEDGMENTS ................................. iv LIST OF TABLES ..................................... viii LIST OF FIGURES .................................... ix ABSTRACT ........................................ xvi CHAPTERS 1 INTRODUCTION ............................... 1 1.1 M otivation . . . . . . . . . . . . . . . . 1 1.2 Contributions of This Dissertation ...................... 3 1.3 Dissertation Overview .............................. 4 2 BACKGROUND ................................ 5 2.1 Introduction . . . . .. .. . . . . . . . .. .. 5 2.2 GPS Overview .................................. 5 2.3 Conventional Attitude Determination Using GPS ............. 7 2.3.1 Attitude and Coordinate Frames ................... 7 2.3.2 Parameterization of Attitude ..................... 8 2.3.3 Conventional Phase Difference Model .............. 10 2.4 Direction Finding .............................. 13 3 AN ANTIJAM ATTITUDE DETERMINATION RECEIVER ....... 18 3.1 Introduction ...... ... .. ....... . ... ........ 18 3.2 A GPS Receiver for AntiJam Position Location and Attitude . .. 18 4 MAXIMUM LIKELIHOOD ATTITUDE ESTIMATION .......... 22 4.1 Introduction .. .... ... ... .. .. .. .. ... ...... 22 4.2 Notation, Preliminaries, and Signal Model................ 23 4.3 Maximum Likelihood Direction Estimation ............. 27 4.4 New Concepts .............................. 30 4.5 Maximum Likelihood Attitude Estimation ............. 32 4.6 Discussion . . . . . .. .. . . . . .. . .. . 37 4.7 Comments on Searches .......................... 40 5 PROPERTIES OF THE MAXIMUM LIKELIHOOD ATTITUDE ESTI M ATO R .. . . .. . . .. . . . .. . . .. .. .. 42 5.1 Introduction . . . . . . . . . . . . . . ... .. 42 5.2 Consistency . . . . . . . . . . . . . . ... .. 43 5.3 Lemmas on the MLAE . . . . . . . . . . . ... .. 44 5.4 B ias . . . . . . . . . . . . . . . . . . 51 5.5 Efficiency . . . . . . . . . . . . . . . ... .. 52 5.6 Conclusions . . . . . . . . . . . . . . ... .. 53 6 ATTITUDE FROM DIRECTION FINDING . . . . . . . ... .. 54 6.1 Introduction . . . . . . . . . . . . . . ... .. 54 6.2 Conceptual Approach . . . . . . . . . . . ... .. 54 6.3 Equal Satellite Weighting . . . . . . . . . . ... .. 56 6.4 Satellite Weighting via Adapted SINR . . . . . . . ... .. 58 6.5 Conclusions . . . . . . . . . . . . . . ... .. 59 7 WIDER BASELINES AND DUAL FREQUENCY USE ........... 61 7.1 Introduction . . . . . . . . . . . . . . ... .. 61 7.2 Dual Frequency Maximum Likelihood Attitude Determination . 63 7.3 Reduction in False Attitudes . . . . . . . . . ... .. 68 7.4 DualFrequency MLAE Performance Increase . . . . . ... .. 80 7.5 Sum m ary . . . . . . . . . . . . . . . ... .. 82 8 SIMULATIONS AND RESULTS . . . . . . . . . . ... .. 84 8.1 Introduction . . . . . . . . . . . . . . ... .. 84 8.2 Simulation Methodology . . . . . . . . . . ... .. 84 8.3 Mean Total Error . . . . . . . . . . . . ... .. 87 8.4 Study 1: Single Jammer with "Random" Location . . . ... .. 88 8.5 Study 2: Varying Attitude Update Rate . . . . . . ... .. 107 8.6 Study 3: Wider Baselines and Dual Frequency . . . . ... .. 121 8.7 Conclusions . . . . . . . . . . . . . . ... .. 126 9 CONCLUSIONS . . . . . . . . . . . . . . ... .. 128 9.1 Summ ary . . . . . . . . . . . . . . . ... .. 128 9.2 Future Work . . . . . . . . . . . . . . ... .. 130 APPENDIX . . . . . . . . . . . . . . . . . . . 133 A CRAMtRRAO BOUND . . . . . . . . . . . . ... .. 133 B THE INTERFERENCE COVARIANCE MATRIX . . . . . ... .. 141 B.1 Introduction . . . . . . . . . . . . . . ... .. 141 B.2 Derivation of Interference Statistics ................... 144 B.3 Estimation of Interference Statistics ................... 156 C RESOLUTION OF THE ATTITUDE AMBIGUITY .. ......... 160 REFERENCES . . . . . . . . . . . . . . . . . . . 164 BIOGRAPHICAL SKETCH . . . . . . . . . . . . . . ... .. 168 LIST OF TABLES Table page 7.1 Number of false solutions, i.e. points below the "possible minimum" thresh old that do not converge to the true solution. Resolution in each Euler an gle is two degrees. Scenario is the same as that involving the first jammer in the random jammer study of Chapter 8 . . . . . . . .... .. 79 8.1 Mean total error performance comparison of the four estimators in an un jammed environment. Update rate is 50 Hz . . . . . . . .... .. 106 8.2 Mean total error performance comparison of the four estimators in an un jammed environment. Update rate is 12.5 Hz . . . . . . ... .. 106 LIST OF FIGURES Figure page 2.1 A two sensor interferometric system. The baseline vector vector dcan be determined by the known LOS vector 6 and the phase difference between the signals measured at the two sensors x0 and x . . . . . .... .. 11 2.2 Mean Total Error Performance of the conventional attitude estimator vs. signal to jammer ratio (SJR), shown as the solid line. The dotted line rep resents the performance of the attitude estimator in an unjammed sce nario. Performance of the phase difference approach to attitude estima tion is degraded even when the SJR is high . . . . . . . .... .. 13 3.1 Receiver block diagram. Demodulated data are provided to the attitude es timator (shown in the circles) as well as to the beamformers. The beam former output is used for position location and updating the code and Doppler tracking loops . . . . . . . . . . . . . . . ... ..20 4.1 Visualization of the possible (ambiguous) attitudes corresponding to a single source. The possible array attitudes can be found by "twirling" the LOS vector, which is assumed affixed to the array, while keeping it pointed at the satellite source . . . . . . . . . . . . . . ... .. 33 4.2 Possible roll and pitch Euler angles arising from the ambiguity in the new array response vector (yaw angles not shown). The true roll is 25 degrees, and the true pitch is 15 degrees, as indicated by the diamond. The locus of possible attitudes is determined by the LOS to the source and the true attitude. A different source LOS would have a different ambiguity. . 34 4.3 Contribution to the value of the maximum likelihood attitude estimator of equation (4.49) from the first of seven satellites in view. This satellite pro vides little information about yaw, and the most in the positive pitch  negative yaw to negative pitch positive yaw dimension . . . ..... .. 38 4.4 Contribution to the value of the maximum likelihood attitude estimator from the fifth of seven satellites in view. This satellite provides little informa tion in the positive pitch negative yaw to negative pitch positive yaw dimension, incidentally the dimension satellite in which satellite one pro vided the most information . . . . . . . . . . . .... .. 39 4.5 Contribution to the value of the maximum likelihood attitude estimator from the seventh of seven satellites in view . . . . . . . . ... .. 40 4.6 Total value of the metric including contributions from all satellites in view. Notice that the areas of weak information from any particular satellite have been filled in by the contributions of other satellites . . . ... .. 41 7.1 Normalized likelihood value vs. direction of arrival for a 15 element ULA re ceiving a narrowband signal at a frequency corresponding to a 2A spacing. True direction is on boresight, i.e. 0 degrees . . . . . . ... .. 69 7.2 Normalized likelihood value vs. direction of arrival for a 15 element ULA re ceiving a narrowband signal at a frequency corresponding to a spac ing. True direction is on boresight, i.e. 0 degrees . . . . .... .. 69 7.3 Likelihood value vs. direction of arrival for a 15 element ULA incorporat ing data collected at two frequencies corresponding to 2A and 5A spac ing. True direction is on boresight, i.e. 0 degrees. The grating lobes are reduced in amplitude from the true direction of arrival, producing a clear indication which direction is correct . . . . . . . . . ... .. 70 7.4 Three dimensional plot of all attitudes whose metric value is below the 10dB threshold, using Satellites 1 and 4, and frequency LI . . . . ... .. 72 7.5 Two dimensional view of all attitudes whose metric value is below the 10dB threshold, using Satellites 1 and 4, and frequency LI . . . . ... .. 73 7.6 Two dimensional view of those attitudes whose metric value is below the 10dB threshold that are not contiguous to the true attitude, (i.e. false so lutions) using Satellites 1 and 4, and frequency LI . . . . . .... .. 73 7.7 False attitude solutions using Satellites 1 and 4, and frequency L2 . ... .. 74 7.8 False attitude solutions using Satellites 1 and 4, and frequency the dual fre quency M LAE . . . . . . . . . . . . . . . ... .. 74 7.9 False attitude solutions using Satellites 1 and 5, and frequency LI . ... .. 75 7.10 False attitude solutions using Satellites 1 and 5, and frequency L2 . ... .. 75 7.11 False attitude solutions using Satellites 1 and 5, and the dual frequency MLAE. 76 7.12 False attitude solutions using Satellites 1 and 6, and frequency Ll . ... .. 76 7.13 False attitude solutions using Satellites 1 and 5, and frequency L2 . ... .. 77 7.14 False attitude solutions using Satellites 1 and 5, and the dual frequency MLAE. 77 7.15 False attitude solutions using Satellites 2 and 3, and frequency Ll . ... .. 78 7.16 False attitude solutions using Satellites 2 and 3, and frequency L2 . ... .. 78 7.17 False attitude solutions using Satellites 2 and 3, and the dual frequency MLAE. 79 8.1 Sensor locations for the three antenna topologies. Quad (asterisk) ,Y (dia mond), and Hex (x). The Y antenna is actually a subset of the Hex an tenna . . . . . . . . . . . . . . . . . . .. 86 8.2 Sensor gain pattern used in the attitude simulation . . . . . .... .. 87 8.3 Sinespace plot of the satellite (asterisk) and the 12 jammer locations (X) for the random jammer study. Only one jammer is simulated at a time; in the first scenario the jammer location is at Ky = .7, and the remaining scenarios use the jammer positions shown in order clockwise. The outer circle represents the horizon, and the area inside the inner circle repre sents the portion of the sky visible to the sensors . . . . . ... .. 89 8.4 Standard deviation of antenna roll estimates and CRB for the Quad antenna (star), Y antenna (square), and Hex antenna (diamond). The 12 points along the abscissa correspond to the 12 jammer locations identified in Fig ure 8.3. The update rate is 50 Hz . . . . . . . . . ... .. 90 8.5 Standard deviation of antenna pitch estimates and CRB for the Quad an tenna (star), Y antenna (square), and Hex antenna (diamond). The 12 points along the abscissa correspond to the 12 jammer locations identified in Figure 8.3. The update rate is 50 Hz . . . . . . . . ... .. 90 8.6 Standard deviation of antenna yaw estimates and CRB for the Quad an tenna (star), Y antenna (square), and Hex antenna (diamond). The 12 points along the abscissa correspond to the 12 jammer locations identified in Figure 8.3. The update rate is 50 Hz . . . . . . . . ... .. 91 8.7 Mean total angle error using the Quad antenna. The 12 points along the abscissa correspond to the 12 jammer locations identified in Figure 8.3. Performance is shown for the MLAE (+), DFW (diamond with solid line), DFU (asterisk), and conventional attitude estimation algorithms (dia mond with dashed line). The update rate is 50 Hz . . . . . ... .. 92 8.8 Mean total angle error using the Y antenna. The 12 points along the ab scissa correspond to the 12 jammer locations identified in Figure 8.3. Per formance is shown for the MLAE (+), DFW (diamond with solid line), DFU (asterisk), and conventional attitude estimation algorithms (dia mond with dashed line). The update rate is 50 Hz . . . . . ... .. 93 8.9 Mean total angle error using the Hex antenna. The 12 points along the ab scissa correspond to the 12 jammer locations identified in Figure 8.3. Per formance is shown for the MLAE (+), DFW (diamond with solid line), DFU (asterisk), and conventional attitude estimation algorithms (dia mond with dashed line). The update rate is 50 Hz . . . . . ... .. 94 8.10 Mean total angle error using the MLAE algorithm, comparing the single jam mer at "random" locations to the unjammed scenario. The solid lines are with the jammer, and the dashed lines unjammed. Performance is shown for the Hex antenna (+), Y antenna (asterisk), and the Quad antenna (di amond). The update rate is 50 Hz . . . . . . . . . ... .. 95 8.11 Mean total angle error using the DFW algorithm, comparing the single jam mer at "random" locations to the unjammed scenario. The solid lines are with the jammer, and the dashed lines unjammed. Performance is shown for the Hex antenna (+), Y antenna (asterisk), and the Quad antenna (di amond). The update rate is 50 Hz . . . . . . . . . .... .. 96 8.12 Mean total angle error using the DFU algorithm, comparing the single jam mer at "random" locations to the unjammed scenario. The solid lines are with the jammer, and the dashed lines unjammed. Performance is shown for the Hex antenna (+), Y antenna (asterisk), and the Quad antenna (di amond). The update rate is 50 Hz . . . . . . . . . .... .. 96 8.13 Mean total angle error using the conventional algorithm, comparing the sin gle jammer at "random" locations to the unjammed scenario. The solid lines are with the jammer, and the dashed lines unjammed. Performance is shown for the Hex antenna (+), Y antenna (asterisk), and the Quad antenna (diamond). The update rate is 50 Hz . . . . . . ... .. 97 8.14 Standard deviation of antenna roll estimates and CRB for the Quad antenna (star), Y antenna (square), and Hex antenna (diamond). The update rate is 12.5 Hz . . . . . . . . . . . . . . . . . . ..98 8.15 Standard deviation of antenna pitch estimates and CRB for the Quad an tenna (star), Y antenna (square), and Hex antenna (diamond). The up date rate is 12.5 Hz . . . . . . . . . . . . . ... .. 99 8.16 Standard deviation of antenna yaw estimates and CRB for the Quad an tenna (star), Y antenna (square), and Hex antenna (diamond). The up date rate is 12.5 Hz . . . . . . . . . . . . . ... .. 99 8.17 Mean total angle error using the Quad antenna. The 12 points along the abscissa correspond to the 12 jammer locations identified in Figure 8.3. Performance is shown for the MLAE (+), DFW (diamond with solid line), DFU (asterisk), and conventional attitude estimation algorithms (dia mond with dashed line). The update rate is 12.5 Hz . . . . ... .. 100 8.18 Mean total angle error using the Y antenna. The 12 points along the ab scissa correspond to the 12 jammer locations identified in Figure 8.3. Per formance is shown for the MLAE (+), DFW (diamond with solid line), DFU (asterisk), and conventional attitude estimation algorithms (dia mond with dashed line). The update rate is 12.5 Hz . . . . ... .. 101 8.19 Mean total angle error using the Hex antenna. The 12 points along the ab scissa correspond to the 12 jammer locations identified in Figure 8.3. Per formance is shown for the MLAE (+), DFW (diamond with solid line), DFU (asterisk), and conventional attitude estimation algorithms (dia mond with dashed line). The update rate is 12.5 Hz . . . . ... .. 102 8.20 Mean total angle error using the MLAE algorithm, comparing the single jam mer at "random" locations to the unjammed scenario. The solid lines are with the jammer, and the dashed lines unjammed. Performance is shown for the Hex antenna (+), Y antenna (asterisk), and the Quad antenna (di amond). The update rate is 12.5 Hz . . . . . . . . . ... .. 103 8.21 Mean total angle error using the DFW algorithm, comparing the single jam mer at "random" locations to the unjammed scenario. The solid lines are with the jammer, and the dashed lines unjammed. Performance is shown for the Hex antenna (+), Y antenna (asterisk), and the Quad antenna (di amond). The update rate is 12.5 Hz . . . . . . . . . ... .. 103 8.22 Mean total angle error using the DFU algorithm, comparing the single jam mer at "random" locations to the unjammed scenario. The solid lines are with the jammer, and the dashed lines unjammed. Performance is shown for the Hex antenna (+), Y antenna (asterisk), and the Quad antenna (di amond). The update rate is 12.5 Hz . . . . . . . . . ... .. 104 8.23 Mean total angle error using the conventional algorithm, comparing the sin gle jammer at "random" locations to the unjammed scenario. The solid lines are with the jammer, and the dashed lines unjammed. Performance is shown for the Hex antenna (+), Y antenna (asterisk), and the Quad antenna (diamond). The update rate is 12.5 Hz . . . . . . ... .. 104 8.24 Satellite (asterisk) and Jammer Locations (x) for the varying update rate study. The jammer indicated by the arrow is the only jammer used in case 1. In case 2, all three jammers appear . . . . . . . . ... .. 108 8.25 Standard deviation of antenna roll estimates and CRB for the Quad antenna (star), Y antenna (square), and Hex antenna (diamond) vs. update rate. One jammer in view . . . . . . . . . . . . . ... .. 109 8.26 Standard deviation of antenna pitch estimates and CRB for the Quad an tenna (star), Y antenna (square), and Hex antenna (diamond) vs. update rate. One jammer in view . . . . . . . . . . . . ... .. 109 8.27 Standard deviation of antenna yaw estimates and CRB for the Quad an tenna (star), Y antenna (square), and Hex antenna (diamond) vs. update rate. One jammer in view . . . . . . . . . . . . ... .. 110 8.28 Mean total angle error using the Quad antenna vs. update rate for 1 jam mer. Performance is shown for the MLAE (+), DFW (diamond with solid line), DFU (asterisk), and conventional attitude estimation algorithms (diamond with dashed line) . . . . . . . . . . . ... .. 111 8.29 Mean total angle error using the Y antenna vs. update rate for 1 jammer. Performance is shown for the MLAE (+), DFW (diamond with solid line), DFU (asterisk), and conventional attitude estimation algorithms (dia mond with dashed line) . . . . . . . . . . . . .... 112 8.30 Mean total angle error using the Hex antenna vs. update rate for 1 jammer. Performance is shown for the MLAE (+), DFW (diamond with solid line), DFU (asterisk), and conventional attitude estimation algorithms (dia mond with dashed line) . . . . . . . . . . . . ... .. 112 8.31 Mean total angle error using the MLAE algorithm, comparing the single jam mer at "random" locations to the unjammed scenario, vs. update rate. The solid lines are with the jammer, and the dashed lines unjammed. Per formance is shown for the Hex antenna (+), Y antenna (asterisk), and the Quad antenna (diamond) . . . . . . . . . . . . .... .. 113 8.32 Mean total angle error using the DFW algorithm, comparing the single jam mer scenario to the unjammed scenario, vs. update rate. The solid lines are with the jammer, and the dashed lines unjammed. Performance is shown for the Hex antenna (+), Y antenna (asterisk), and the Quad antenna (di am ond). . . . . . . . . . . . . . . . . . . 114 8.33 Mean total angle error using the DFU algorithm, comparing the single jam mer scenario to the unjammed scenario, vs. update rate. The solid lines are with the jammer, and the dashed lines unjammed. Performance is shown for the Hex antenna (+), Y antenna (asterisk), and the Quad antenna (di am ond). . . . . . . . . . . . . . . . . . . 114 8.34 Mean total angle error using the conventional algorithm, comparing the sin gle jammer scenario to the unjammed scenario, vs. update rate. The solid lines are with the jammer, and the dashed lines unjammed. Performance is shown for the Hex antenna (+), Y antenna (asterisk), and the Quad antenna (diamond) . . . . . . . . . . . . . . ... .. 115 8.35 Mean total angle error using the Quad antenna vs. update rate for 3 jam mers. Performance is shown for the MLAE (+), DFW (diamond with solid line), DFU (asterisk), and conventional attitude estimation algo rithms (diamond with dashed line) . . . . . . . . . .... .. 116 8.36 Mean total angle error using the Y antenna vs. update rate for 3 jammers. Performance is shown for the MLAE (+), DFW (diamond with solid line), DFU (asterisk), and conventional attitude estimation algorithms (dia mond with dashed line) . . . . . . . . . . . . ... .. 117 8.37 Mean total angle error using the Hex antenna vs. update rate for 3 jam mers. Performance is shown for the MLAE (+), DFW (diamond with solid line), DFU (asterisk), and conventional attitude estimation algo rithms (diamond with dashed line) . . . . . . . . . .... .. 117 8.38 Mean total angle error using the MLAE algorithm, comparing the three jam mer scenario to the unjammed scenario, vs. update rate. The solid lines are with the jammer, and the dashed lines unjammed. Performance is shown for the Hex antenna (+), Y antenna (asterisk), and the Quad antenna (di am ond) . . . . . . . . . . . . . . . . . . . 118 8.39 Mean total angle error using the DFW algorithm, comparing the three jam mer scenario to the unjammed scenario, vs. update rate. The solid lines are with the jammer, and the dashed lines unjammed. Performance is shown for the Hex antenna (+), Y antenna (asterisk), and the Quad antenna (di am ond). . . . . . . . . . . . . . . . . . . 119 8.40 Mean total angle error using the DFU algorithm, comparing the three jam mer scenario to the unjammed scenario, vs. update rate. The solid lines are with the jammer, and the dashed lines unjammed. Performance is shown for the Hex antenna (+), Y antenna (asterisk), and the Quad antenna (di am ond). . . . . . . . . . . . . . . . . . . 119 8.41 Mean total angle error using the conventional algorithm, comparing the three jammer scenario to the unjammed scenario, vs. update rate. The solid lines are with the jammer, and the dashed lines unjammed. Performance is shown for the Hex antenna (+), Y antenna (asterisk), and the Quad antenna (diamond) . . . . . . . . . . . . . . ... .. 120 8.42 Sensor Locations used in the wide baseline study. The three antenna topolo gies are the Quad (asterisk),Y (diamond), and Hex (x) . . . . ... .. 122 8.43 Mean total angle error using the Quad antenna vs. update rate for 1 jam mer. Performance is shown for the single frequency MLAE at LI (dia mond) and L2 (+), and for the dual frequency MLAE (asterisk) . . . 122 8.44 Mean total angle error using the Y antenna vs. update rate for 1 jammer. Performance is shown for the single frequency MLAE at LI (diamond) and L2 (+), and for the dual frequency MLAE (asterisk) . . . ... .. 123 8.45 Mean total angle error using the Hex antenna vs. update rate for 1 jammer. Performance is shown for the single frequency MLAE at LI (diamond) and L2 (+), and for the dual frequency MLAE (asterisk) . . . ... .. 123 8.46 Mean total angle error using the Quad antenna vs. update rate for 3 jam mers. Performance is shown for the single frequency MLAE at LI (dia mond) and L2 (+), and for the dual frequency MLAE (asterisk) . . . 124 8.47 Mean total angle error using the Y antenna vs. update rate for 3 jammers. Performance is shown for the single frequency MLAE at LI (diamond) and L2 (+), and for the dual frequency MLAE (asterisk) . . . ... .. 124 8.48 Mean total angle error using the Hex antenna vs. update rate for 3 jam mers. Performance is shown for the single frequency MLAE at LI (dia mond) and L2 (+), and for the dual frequency MLAE (asterisk) . . . 125 B.1 The receiver model used in this dissertation. The output of the chip matched filter is oversampled by a factor of N,, (N, is 4 in this figure) and these samples are decimated into sequences composed of 1 sample per chip. The decimated sequences are then despread using the sequency A(i) i = 1,2,... N to form early, punctual, late, and other temporal gates . . . ... .. 143 B.2 One sided jammer PSD . . . . . . . . . . . . ...... 146 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy INTERFERENCE MITIGATION FOR GPS BASED ATTITUDE DETERMINATION By Matthew David Markel May 2002 Chair: Henry Zmuda Major Department: Electrical and Computer Engineering With its low transmit power and considerable orbit distance from receivers, the global positioning system (GPS) waveforms are known to be susceptible to jamming. Efforts to date have investigated various receiver designs and adaptive antenna arrays to improve the antijam capabilities of GPS position location, but until this work no extension of antijam capabilities to attitude determination (determination of the orientation of a body in space) existed. This dissertation provides a comprehensive investigation into the new field of robust, jamresistant attitude determination using the Global Positioning System. We present a generic adaptive antenna array and receiver design that provide the necessary data for both position location and attitude determination in a jammed envi ronment. Using the data this receiver provides we develop a new maximum likelihood attitude estimator (the MLAE) that provides attitude determination capability even in jammed environments. In order to accomplish this, new ways of viewing the standard array processing and direction finding tenants are presented. This leads to a new in terpretation and parameterization of the array response vector (spatial steering vector) and secondary data covariance matrix. The MLAE optimally includes information from all satellites in view, and its performance is shown herein to asymptotically achieve the Cram6rRao Bound (CRB), i.e. the MLAE asymptotically achieves the performance limit for unbiased estimators. As an approximation to this estimator, a second estimator is developed that, at the expense of performance, may provide increased computational capability. Two versions of this second estimator are presented that combine direction finding with attitude estimation. A method of optimally incorporating data from both GPS frequencies into the estimation is developed. This approach maintains the jammer mitigation capabilities of the previous estimators and allows for a larger array spacing, therefore increasing the accuracy of the attitude estimates. Simulation based performance results of the various estimators are presented for various antenna topologies and interference scenarios. The simulation results indicate that the new algorithms provide significant improvement over conventional attitude determination methods. In addition, the MLAE is shown to provide better performance than conventional methods of attitude determination even in unjammed environments. CHAPTER 1 INTRODUCTION 1.1 Motivation It is well known that adaptive antenna arrays (often referred to as "Smart Antennas" for communications systems) provide significant resistance to unintentional interference and intentional jamming for both signal extraction and direction finding. Global Position ing System (GPS) based attitude determination (AD) systems utilize multiple sensors as well to extract attitude through carrier phase differences between the sensors [1]. However, present AD algorithms typically only offer the jamming resistance inherent in the receiver (i.e. gain from the spread spectrum waveform and perhaps coupling with some form of INS), which for even low powered jammers may not be enough to prevent corruption of the attitude estimates. Therefore, it seems appealing to exploit the similarities between the two fields of adaptive array processing and GPS based attitude determination to in crease the antijam capabilities of GPS attitude systems. Indeed, this work centers around a substantially different method of approaching the attitude determination task than is typically employed. The scenario that originally motivated this research involved a medium to long range platform that must fly for extended periods in a severe jamming environment. The platform incorporates an adaptive antenna array to provide GPS jamming resistance for position location and navigation. In response to this scenario, a method has been developed that makes use of an antijam antenna array to provide GPS based attitude information, as well as the typical position location data, even during periods of jamming. Three distinct factors motivated this research into antijam attitude determination. 1. Threat from GPS jammers. The susceptibility to jamming and inadvertent interference is high due to the GPS satellites' low transmit power and considerable orbital altitude. Even a low power jammer several miles from the GPS user can degrade accuracy or affect acquisition. With GPS jammers being marketed as commercial items [2], one could control access to GPS over an entire region by deploying several of these low cost/low power jammers throughout an area and activating them only when desired to deny the use of GPS. 2. Related research. Previous research in direction finding algorithms using antenna arrays for radar and communications systems provides a substantial mathematical basis and a likelihood of tractability for the approaches developed in this dissertation. At the macro level, attitude determination and direction finding are similar. For example, they both incorporate multiple sensors and use carrier phase interferometry to develop their estimates. Since some direction finding algorithms provide significant performance in external interference environments, it is natural to to look to this field with the desire to gain similar antijam capabilities for attitude determination. The primary algorithm developed in this dissertation does not use direction finding directly. However, it is not unlike many direction finding algorithms in the sense that it defines steering vectors that map unknown parameters into a model for observation data from the antenna array, and the unknown parameters are estimated by maximizing (or minimizing) a function of the steering vectors and received data. Suboptimal algo rithms developed later in the dissertation directly incorporate direction finding, or more specifically assume that the directions are measured by some method. 3. Hardware trends. State of the art systems are incorporating adaptive antenna arrays with the GPS receivers to provide additional antijamming resistance for position location, providing much of the hardware and design philosophy for attitude location as well (see, for example, Falcone et al. [3]). A possible direct application of this work would be to incorporate the algorithms developed here into a GPS system that already incorporates multisensors for antijam position location and navigation in order to additionally provide a GPS based attitude estimate. Responding to these three motivations (i.e. need, related research, and hardware technology trends), it is natural to investigate a system that makes maximum use of the presence of the antijam sensor array by providing attitude as well as position location. 1.2 Contributions of This Dissertation This work draws on several separate fields of research. The following background chapter provides an overview of the Global Positioning System, attitude and associated "conventional" (i.e. not robust to interference) methods of attitude determination and direction finding because all of these fields were used to develop new interference resistant attitude estimation methods. These results are novel because external interference mitigation for GPS based attitude estimation is a completely undeveloped field [4], and the results derived here provide significant improvement over conventional attitude determination algorithms in both jammed and unjammed environments. The constraints and challenges of this task are different from the "typical" antenna array direction finding and communications applications. 1 To address these challenges the direction finding task is recast as an "attitude finding" problem by redefining such standard and familiar terms as the antenna array response vector (often called the spatial steering vector) and data matrix. This provides increased insight into methods to address the problem and facilitates development of a new maximum likelihood attitude estimation method. This estimator is shown to be asymptotically unbiased, consistent, and efficient (i.e. asymptotically achieves the Cram6rRao Bound). Another contribution is the extension of maximum likelihood direction/attitude estimation to incorporate both GPS carrier frequencies, as discussed in Chapter 7. Using the method developed in this dissertation, incorporation of the additional GPS frequency provides a signal to noise ratio (SNR) gain, as opposed to methods in use today in which the SNR is reduced when the additional frequency is used. Finally, the simulation based performance results provide a quantification of the performance of this work. 1 Of course, caution is used with the generalization of "typical," since every application has some unique attributes. 1.3 Dissertation Overview The remainder of this dissertation is organized as follows. Chapter 2 contains an overview of the multiple research fields from which this dissertation draws. Chapter 3 describes a generic receiver that provides the data necessary to estimate both position and attitude, and the differences between this receiver and those typically employed for either antijam position location or attitude determination. Chapter 4 derives a novel algorithm for attitude estimation in an external interference environment. This estimator takes the form of a maximumlikelihood estimator. Chapter 5 presents several statistical properties of this estimator. In Chapter 6, two variants of suboptimal algorithms are developed for attitude estimation. These algorithms have attractive implementation and computational features, but their performance is not as good as the estimator of Chapter 4. Chapter 7 derives a method of improving the attitude estimator of Chapter 4 by using both GPS frequencies. Chapter 8 presents quantification of the performance of these methods through simulation results. Finally, Chapter 9 contains a summary of this work and a discussion of areas of research related to this work, but beyond the scope of this dissertation. CHAPTER 2 BACKGROUND 2.1 Introduction This chapter presents an overview of the major fields which provide the foundation for the contributions of this dissertation. Interference mitigation for GPS based attitude determination involves the synthesis of several different fields, so a brief introduction and summary of them is provided here. This chapter is organized as follows. Section 2.2 provides an overview of the Global Positioning System. Section 2.3 discusses methods and signal models used in attitude determination. These two sections provide the framework for the term "conventional" GPS based attitude determination, which is used in this work to describe the methods of attitude estimation used today. Finally, a survey of direction finding is provided in section 2.4, with an emphasis on maximum likelihood methods. 2.2 GPS Overview From personal, handheld civilian units to integrated military systems, the Global Positioning System (GPS) is an important link in the worldwide navigation infrastructure. The Global Positioning System is designed to provide a user with the ability to calculate his position in three dimensional space (i.e. latitude, longitude, and altitude) through reception of signals from a constellation of 24 dedicated satellites orbiting the earth. The Global Positioning System conceptually consists of three segments. The space segment involves the satellites in orbit about the earth. These satellites provide a low power output signal that may be received on or above the earth. The control segment tracks and communicates with the satellites in order to update their navigation informa tion. Finally, the user segment represents the receiver and associated equipment required to decode and interpret the signals sent by the satellites and convert this into the user's position and velocity. With an impressive market penetration, the number of "users" occupying the user segment is steadily increasing. The ability to calculate position arises from measuring the range to several satellites. If one knows the range to a fixed object, the locus of possible user positions is on the surface of a sphere with its center at the object position and its radius equal to the range to the object. Therefore, if the user knows his range to three satellites, the positions of the satellites (calculated using the satellite's ephemeris parameters), and has a precise measurement of time, then he can solve for his position as the intersection of the three spheres. In practice, less precise clocks are used in the user segment, requiring a minimum of four satellites to obtain a position. Obviously, errors and uncertainty in the ability to measure the range to the satellite directly affect the accuracy of the position measurement. The GPS satellites incorporate Direct SequenceSpread Spectrum (DSSS) modu lation. This allows all satellites to occupy the same spectrum and, by monitoring the code tracking loop, enables users to accurately calculate range to the satellite. Two levels of accuracy are currently available through the use of two different spreading sequence architectures. The Standard Positioning Service, which is available to all users, incorpo rates a 1023 length sequence with a chip rate of 1.023 MHz called the Coarse Acquisition Code (i.e. C/A Code). The Precise Positioning Service (PPS) uses the P(Y) code that operates at a chip rate of 10.23 MHz. The P(Y) codes have very long natural periods (over 6 months); however, they are truncated and restarted every week. While the C/A Code waveform is publicly known and easily obtainable, the P(Y) code sequence is only available to authorized users. Since knowledge of the spreading sequence is required for demodulation and data detection, the PPS is essentially available only to military applications. The C/A waveform (code and data sequence) modulates a carrier frequency of 1.57542 GHz (often referred to as the LI frequency). Since several satellites are simultane ously transmitting at the same frequency, Code Division Multiple Access (CDMA) is used to allow unambiguous satellite reception. Each satellite uses a different sequence from a set of 10th order Gold codes. Of the 1025 possible Gold codes available in this family, 37 have been identified for use on satellites. The P(Y) code is broadcast on two frequencies, the LI frequency defined above (in phase quadrature with the C/A code) and the L2 frequency of 1.2276 GHz. At the time of this work, an effort is underway to enhance GPS with additional frequencies for civilian use and a different waveform, denoted as the M code, for military and other restricted use. The navigation message data bits (which are modulated by the spreading sequences described above) are transmitted at a rate of 50 Hz. The GPS navigation data consist of synchronization information, satellite ephemeris data, and other satellite parameters. These data are designed to be packaged into 30 second (i.e. 1500 bit) time slices. An important observation is that unlike many other communication systems where the only information is the demodulated data stream, GPS uses the code timing and tracking information to develop an estimate of range to the satellites. This range, coupled with the other information obtained from the GPS data stream, allows a formulation of user position. 2.3 Conventional Attitude Determination Using GPS 2.3.1 Attitude and Coordinate Frames Attitude determination involves two coordinate systems and the transformation that relates them. The locallevel coordinate system represents an initial, at rest, or other nominal orientation of the platform of interest. This frame is often referenced in terms of a "NorthEastDown" (NED) or "EastNorthUp" (ENU) convention. The line of sight (LOS) vectors to the satellites are assumed known in this coordinate system. The second coordinate system is the antenna frame of reference, and represents the orientation of the GPS antenna at a particular point in time. The antenna frame and the body frame' are related by a known transform; if one is known, the other can be easily found. For this reason the antenna frame is sometimes referred to in this work as the body frame and is 1 The body frame is defined as the frame of reference along principal platform axes, such as along the fuselage and wings. Typically the body frame is of higher interest from a guidance and navigation perspective. usually referenced using a "FrontRightBelow" (FRB) convention. The attitude is then defined as the relation between the locallevel and the antenna frames of reference at any point in time. 2.3.2 Parameterization of Attitude Some common parameterizations of the attitude include Euler angles, quaternions, and direction cosine matrices. Throughout this work, this attitude is referred to generi cally as q, to remain independent of any particular parameterization.2 However, at times a specific parameterization is required either to clearly demonstrate a point or prevent ambiguities. When specific parameterizations of this attitude are employed they follow the notation and description below. q = the attitude of the antenna S= the antenna roll Euler angle e = the antenna pitch Euler angle S= the antenna yaw Euler angle q = [qi, q2, q3, q4J, the attitude quaternion Q = the 3x3 direction cosine matrix that transforms local level to antenna frame The Euler angles represent successive rotations about the three coordinate axes. The roll angle, 4,, represents rotation about the local level x axis, the pitch angle 6e about the y axis, and the yaw angle 6 the z axis. The order of rotation determines the final orientation. In this work we employ a "312" rotation, meaning the rotations are performed first in yaw, then roll, and then pitch. A drawback with the Euler angle 2 Notice that the unparameterized attitude q in this work is not represented in bold type, even though at least three parameters are required to specify the attitude. This is to emphasize that when this notation is used, the attitude may be specified in any manner. representation is the possibility of encountering singularities; however, in general they provide an intuitive, clear interpretation of the rotations. Another method of defining attitude is with the direction cosine matrix Q. Each term in the direction cosine matrix represents the cosine of the angle between an axis of the local level frame and the antenna frame. For example, Q(2,2) is the cosine of the angle between the y axis of the local level frame and the y axis of the antenna frame [5]. The direction cosine matrix is orthonormal. QT = Q1 (2.1) IQI = 1 (2.2) This orthogonal nature of the direction cosine matrix is very useful, since an inverse rotation (say, transforming from the body frame back to the local level frame) can be performed without inverting the matrix. The final method of specifying attitude is with a quaternion. In order to define the quaternion's ability to represent a rotation, first consider that between any two coordinate systems sharing a common origin, a single axis exists in which measurements along that axis are the same in both coordinate systems. The relation between the two coordinate systems (i.e. the attitude) is then defined by this axis and an angle of rotation about it. If we take the axis as V and the angle of rotation as w, then V. V = V (2.3) Where where IVI = 1 (2.4) and using these the quaternion q may be written as V. q= sin(f) V 25 q 2 Y(2.5) V. cos( ) and q = 1 (2.6) The identity quaternion, which performs no rotation, is 0 0 qIdentity (2.7) 0 1 For an excellent discussion and history of attitude and its parameterizations, the reader is referred to Schleppe [5]. 2.3.3 Conventional Phase Difference Model In order to describe the operation of a conventional attitude determination system, we first describe the signal model employed and the assumptions it incorporates. Consider the two GPS sensor scenario of Figure 2.1. After downconversion, despreading, and integration the signal received at sensor x0 consists of a contribution from the satellite and a noise term. XO = ge3 + no (2.8) The response for the second sensor may be written in a similar manner. xi = ge +60) + ni (2.9) Notice that the second sensor also receives the same satellite signal only shifted in time (i.e. phase). This phase difference is related to the projection of the baseline vector donto the LOS and the wavelength: S27r  &g=  vid A (2.10) The noise terms no and nl are assumed uncorrelated and moreover are assumed to be small compared to the received gain from the satellite signal g. The attitude determi ~\f d= (V/27) 86/ . V Sensor Xo d Sensor Xl Figure 2.1: A two sensor interferometric system. The baseline vector vector dcan be determined by the known LOS vector 6 and the phase difference between the signals mea sured at the two sensors x0 and xi. nation problem is to estimate the baseline vector d. It is clear to see that by taking the difference A0 between the phase measured at x0 and xl, then the following relationship between the orientation of the vector d, the known LOS vector to the satellite V, and the phase difference may be obtained as [6] (2.11) (2.12) (2.13) where 4)(.) implies taking the phase angle of (.), i.e. O= (xl) (x0o) JO + k(27r) (D() = j ln(.) (2.14) and ; implies "is approximately equal to." The reason that equation (2.13) holds is due to the previous assumption that the amplitude of the interference terms is much smaller than the amplitude of the satellite terms, and therefore the satellite components dominate the phase of the sum of satellite signal and noise. The k(27r) term represents the integer portion of the phase difference be tween the two sensors and is called the "integer ambiguity," since it cannot be measured. By using at least three sensors (in a noncollinear arrangement) and two satellites [7], the complete attitude of the antenna may be found, once the integer ambiguity values have been found. Several methods (see for example Sutton [8] and its references) exist for resolving the integer ambiguities. Since the attitude estimation process makes use of knowledge of the LOS vectors to the satellites in the local level frame, a question to consider is the sensitivity to error in the antenna position estimate. It turns out that small errors (within typical GPS accuracies) in position estimation will not have a substantial impact on the LOS angles for essentially every ground and air based application. This is due to the large distance to the satellites of approximately 23,000km from the earth's center. Even for very high altitudes of 10Okft, the angular error is negligible. The extension of equation (2.13) to multiple sensors and satellites provides a straight forward and computationally simple method to estimate the antenna attitude when the assumptions of the signal model are met. However, when an interference source such as a jammer is present, the accuracy quickly degrades. Figure 2.2 indicates the level of degradation a jammer can cause to an attitude determination system. The parameter "mean total error" is the average of the total errors from 500 realizations of the attitude estimation, where the total error is defined as the angular rotation between the estimated attitude (using the attitude estimator defined above) and the true attitude.3 Notice that even with relatively high signal to jammer ratios, the performance is degraded. The 3 Chapter 8 provides a mathematical definition of the mean total error. fundamental thrust of this dissertation is the development of algorithms that increase the resistance of the attitude estimator to external interference. Mean Total Error 251 i \ Jammer Present [ ammr Asnt 20 0 ...........2 15 ........ ...... ........ ................. ............... .... ........ 1 0 . .. . . .. . .. .... . . .. . .. . . .. . . . . . . . . 10 ............ ... .... . .. ...... 5\ 0 I I !I 5 0 5 10 15 20 25 SJR at Phase Measurement (dB) Figure 2.2: Mean Total Error Performance of the conventional attitude estimator vs. signal to jammer ratio (SJR), shown as the solid line. The dotted line represents the per formance of the attitude estimator in an unjammed scenario. Performance of the phase difference approach to attitude estimation is degraded even when the SJR is high. 2.4 Direction Finding The direction finding problem is the determination of the directions of arrival (DOA) of multiple sources impinging on an array of sensors. Since the problem may be viewed as estimating how radiated energy is distributed over space, it is often referred to as "spatial spectral estimation" [9]. Depending upon the application, DOA may be considered as a single parameter (e.g. azimuth angle) or two parameters (azimuth and elevation angle). When the application only requires a single direction parameter, simple antenna topologies like the uniform linear array (ULA) may be employed. Several algorithms make use of this topology to provide computationally efficient estimates of direction. However, since attitude estima tion is intrinsically a multidimensional problem requiring full angular information to the satellite source, little insight can be obtained from these single angle parameter algorithms for the attitude determination application. Therefore, we have limited the discussion of this section to those algorithms that can provide both parameters of direction to the source. This is not to say that the single DOA parameter algorithms are without merit, just that they are not applicable for attitude estimation. Through the last several decades advances in direction finding theory have been ex tensive, so much so that even a survey of the major contributions is difficult. One method of gaining insight into these advances is to partition direction finding approaches into the two classes of maximum likelihood and suboptimal. As is often the case, suboptimal algorithms offer computational savings, but at the cost of reduced performance. The first algorithms we consider are the maximum likelihood estimators. These optimal estimators may be further partitioned by considering the parameters estimated. The earliest works in this field were improvements to radar and sonar systems that considered the angle and complex amplitude estimation problem for a single source in a jammed environment [10]. These works focused on maximum likelihood estimators approximated by adaptive monopulse ratios (i.e. adapted delta to adapted sum ratio) [11, 12]. Although four unknowns exist in the scenario addressed by these algorithms (two parameters of the complex gain and two of the angle), only a two dimensional search (over the two angle parameters) of the maximum likelihood metric is required, since the complex gain at the solution point can be explicitly expressed in terms of the angle. More recently, adaptive array technology is beginning to be applied to communication systems to aid in multiuser wireless systems [13]. In this extension of the estimation problem, the unknown parameters estimated are not only the directions and a complex gain, but also an entire signal waveform for multiple signals impinging simultaneously on the array. For N snapshots of data and L sources, there are (2N + 2) x L real parameters to estimate (the complex gain is incorporated into the unknown waveform). As with the single source problem, the waveform parameters may be expressed in terms of the unknown angles, resulting in a search of the likelihood function over the 2L angular parameters. Two distinct signal models have been used in the multiple source plus waveform estimation problem, with differing results. The deterministicc" or "conditional" maximum likelihood approach considers the source waveforms to be unknown deterministic signals in an additive Gaussian noise environment. The noise is often assumed temporally uncorrelated, complex circularly symmetric, and spatially white. The mechanization of the maximum likelihood metric for estimating the directions is to fit the subspace spanned by the matrix of array response vectors (also called steering vectors) to the sources to the measurements in a least squares sense [14,15]. Several direction finding algorithms are cast as variations of the weighted subspace fitting algorithm in Viberg and Ottersten [15]. Some important results have been found concerning the conditional maximum likelihood signal model assumption. First, the direction estimates obtained from this estimator are consistent [16]. Second, because the number of waveform parameters required to estimate increases with each new snapshot of data, the waveform estimates are not consistent, and the estimator is not statistically efficient, i.e. does not achieve the Cram6r Rao Bound (CRB), for a finite number of sensors [16]. However, as the number of sensors increases, or the signal to noise ratio increases, this estimator asymptotically approaches the CRB [17]. The "stochastic" or "unconditional" maximum likelihood approach is derived from a model where both the signals and the noise are stationary, temporally white, zeromean complex Gaussian random processes [18]. It has been shown that the unconditional maximum likelihood method is statistically more efficient than the conditional ap proach [17,19]. A more computationally efficient implementation that uses a Cholesky decomposition instead of an eigenvalue decomposition is presented in Swindlehurst [18]. The next class to consider is estimation of the directions when the signal waveforms are known a prior to the estimation process. This situation occurs in communications systems when a known special preamble is added to a packet of data. In addition, the GPS waveforms fall into this category once code lock is obtained. As expected, significant performance gains may be obtained when the signal waveforms are (to within a complex gain) known as compared to when they must additionally be estimated. In Li [20] and Li and Compton [21] the CRBs for several cases where the signal waveform is known to within a complex gain (and the gain is either known or partially known) are derived and compared to the CRB when the signal waveform is completely unknown. Finally, Li et al. show that the estimation process of the directions to the multiple sources may be performed independently when each of the source waveforms are uncor related with each other [22,23]. This is the socalled decoupled maximum likelihood, or "DEML" estimator. The DEML is shown to be consistent and asymptotically efficient for uncorrelated waveforms. In addition, the performance of this estimator does not degrade as either the number of signals increases or the angular separation between the sources decreases. The concepts presented in these two works provide the most direct applicability of direction finding to the attitude estimation process derived later in this work, since the GPS P(Y) waveforms are indeed uncorrelated and known, and the number of source (i.e. satellite) waveforms impinging upon the antenna array is typically greater than the number of sensors. Several suboptimal direction finding algorithms have been introduced in addition to the maximum likelihood methods discussed above, and a few of the more popular algo rithms are now briefly discussed. The Multiple Signal Classification (MUSIC) algorithm is based on an eigenvalue decomposition of the sample covariance matrix [24]. The algo rithm partitions this eigenvalue decomposition into "signal" and "noise" subspaces and, by observing that steering vectors to the sources present will lie in the signal subspace, estimates directions by minimizing projections of steering vectors onto the noise subspace. Typically, this will involve a two dimensional search across angle. However, for some array topologies a significantly simpler rooting algorithm named PRIME may be employed [25]. An algorithm that relaxes the requirement for a calibrated array (one where the steering vector is known for all possible direction of arrival angles) is the Estimation of Signal Pa rameters via Rotational Invariance Techniques (ESPRIT) [26]. Although this is typically an azimuth only algorithm and would therefore not meet our discussion criteria, it has been extended to two parameter angle estimation by Swindlehurst and Kailath [27]. 17 Finally, for two excellent survey articles on a variety of direction finding algorithms, the reader is referred to Godara [28] and Krim and Viberg [29]. In addition to a summary of the algorithms, Godara [28] presents several implementation issues common in array processing such as array calibration (knowledge of the exact array response vector), failed sensors, and errors in beamformer weights. CHAPTER 3 AN ANTIJAM GPS RECEIVER FOR ATTITUDE DETERMINATION 3.1 Introduction In this chapter we develop a functional GPS receiver architecture that provides data for both position location and attitude determination in a jammed environment. We apply the following notation when describing the receiver. A hardware channel is the portion of the receiver from the antenna through the analog to digital converter, containing the mixers for down conversion, filters, amplifiers, etc. A satellite channel performs the Direct SequenceSpread Spectrum (DSSS) demodulation, Doppler demodulation, and integration for one satellite. Outputs from the satellite channel are used to keep the satellite waveform in code and Doppler track, as well as to provide data for position location. As an example, a simple singleantenna GPS receiver that can simultaneously track 5 satellites would have (without incorporating any multiplexing) one hardware channel and 5 satellite channels. Using this notation, we review "conventional" attitude determination and antijam receivers, and then contrast them with a proposed receiver designed to implement the algorithms of Chapters 4 and 6. 3.2 A GPS Receiver for AntiJam Position Location and Attitude Attitude determination from GPS is based on measuring the satellite time (phase) differences of arrival between multiple sensors. Using the terminology above, an attitude determination system with M antennas that tracks L satellites would have M hardware channels, and (again ignoring multiplexing) L satellite channels per hardware channel. Since phase differences, not the absolute phase of the carrier sinusoid, are used to deter mine attitude, the M hardware channels may be physically separate and independent receivers using different clocks, local oscillators, etc. The output phases from each of the L satellite channels for all M sensors, with knowledge of the sensortosensor base line vectors, provide the data necessary for attitude determination in a low interference environment. Historically, antijam antenna arrays have operated as preprocessors to a GPS receiver [3,3032]. The function of the antijam antenna is to weight and combine the voltages received at the various sensors into one time signal, which is then passed to the single hardware channel. With the interference attenuated by the adaptive array, the tracking loops following the satellite channels are able to keep the satellites in track and provide data for position location. Note that since the antijam antenna is often designed as an applique, the postantenna portion of the GPS receiver operates as if a single sensor had provided the time signal, even though it is actually the linear combination of data from several sensors. It is clear that a preprocessor type antijam antenna, although it uses multiple sensors, does not provide the data required for attitude determination. The data signal presented to the hardware chanel is the weighted combination of the data collected by the multiple sensors. Even if the weights are known, determining the satellite phase differences across sensors from the weighted combination of these sensors is an underdetermined problem. However, the receiver can collect the data required for both position location and attitude determination by increasing the number of hardware channels and incorporating the weight application into the receiver. Figure 3.1 shows the top level architecture for one possible GPS system that provides adaptive arrayprocessing levels of antijam performance for both position location and attitude determination. This system begins with an array of M sensors. The output of each sensor is downconverted from the GPS carrier frequency to an intermediate frequency (IF), using Inphase (I) and Quadrature (Q) downconversion. The same mixing signals (from the same local oscillator) are used in each of the M hardware channels. The signals in each channel pass through an antialiasing matched filter and are sampled in I and Q. For ease of notation and discussion, we define a Vector Demodulation block as a col lection of M satellite channels, one per hardware channel, with each block demodulating Figure 3.1: Receiver block diagram. Demodulated data are provided to the attitude es timator (shown in the circles) as well as to the beamformers. The beamformer output is used for position location and updating the code and Doppler tracking loops. the same satellite spreading sequence. It is useful to visualize this process as a block operation on all M signals simultaneously, since the same satellite demodulation sequence is applied to all M signals from the hardware channels. The outputs of this vector demod ulator are M x 1 data vectors for the various code delays required to track the satellite signal in time, e.g. an early, punctual, and late delay. More than three code delays may be implemented to increase satellite acquisition time or to provide a quicker reacquisition time if the satellite signal is temporarily lost. For the multisensor receiver to process L satellites, it needs L vector demodulators. The outputs of the L vector demodulators, ul, are used separately for two functions, updating the tracking loops to provide position location information, and as input to the attitude determination algorithm. Although the demodulator output vectors have received the spread spectrum process ing gain, in a jammed environment the signal to interference plus noise ratio (SINR) for each vector element may still be too small for satellite tracking until beamformer weights are applied. By projecting the M x 1 data vectors onto the adaptive beamformer weights, the SINR is increased to an acceptable level to calculate the discriminants used for code and Doppler tracking. Ignoring the effects of finite precision arithmetic, the beamformer weights may be applied before the satellite channel demodulation (as with the antijam array preprocessor) or after satellite channel demodulation (as is shown here) with equiv alent results. In addition, the method proposed here allows for different weights for each satellite, improving tracking performance over a single set of weights for all satellites. The second use of the data vectors is for attitude determination. As discussed above, attitude determination requires the data from multiple sensors to determine the time (phase) differences. Therefore it must operate on data obtained prior to beamforming, even though this may contain significant interference power, since as mentioned earlier the beamforming operation destroys the individual sensor phases. The algorithms developed in the next section estimate attitude even with the interference present. In summary, this section developed the toplevel architecture for a multisensor GPS receiver that combines the antijam benefit of adaptive beamforming with development of data required for antijam attitude determination algorithms of the next section. This system is like conventional attitude determination systems in that it uses data from multiple sensors and hardware channels. However, it differs in that adaptive weights are applied after the satellite channel demodulation, and that the data passed to attitude determination algorithms are the complex data vectors for each satellite and sensor, not just the phase information. This system is similar to an adaptive array preprocessor architecture, but differs in that multiple hardware channels (instead of a single channel) are employed, and the array weights are applied after demodulation. CHAPTER 4 MAXIMUM LIKELIHOOD ATTITUDE ESTIMATION 4.1 Introduction In this chapter we derive a new method of for GPS based attitude estimation. This method uses the receiver construct developed in Chapter 3. Motivated by Figure 2.2, the desire is to develop an attitude estimator that is jammer resistant. This development begins by evaluating the mathematics of jammer resistant direction finding algorithms, and using this background rederives the approaches in terms of the attitude estimation field. This chapter is organized as follows. Section 4.2 provides a discussion of the signal models and notation used throughout this chapter. In section 4.3 the multisource direction finding problem is reviewed, to introduce the multisource maximum likelihood method. Specifically in this section the task at hand is to develop direction estimates to each satellite in view. Section 4.4 presents some new and modified concepts that are required for the new attitude estimator, which is developed in section 4.5. Finally, section 4.6 provides a discussion of this estimator and some graphical illustrations of its performance. 4.2 Notation, Preliminaries, and Signal Model We begin this section by reviewing the matrix and vector notation used throughout this work. AT = the transpose of the matrix A AH = the conjugate transpose of A A B = the Kronecker (Tensor) product of A and B Ij = the j x j identity matrix 0j = the j x j matrix of zeros A> 0 = A is positive definite Z E Cxb = Z is a complex matrix of size a x b 71 0 0 diag[7,y,..., ,L] 0 "'. 0 0 0 7LJ ql q2= quaternion multiplication TR(A) = trace of A In addition, the following variables are consistently used in the following roles. L = the number of satellites being tracked M = the number of sensors in the array N = the number of snapshots of data available (4.1) Consider the following popular baseband model for the signals received by the M element array of antennas. The total data, x(t), received by the M elements are composed of the contributions from each of the L satellites, pl(t), and the interference n(t): L x(tn)= pi(t.n)+n(tn), n = 1,...,N (4.2) /=1 where x(t), p1(t), and n(t) E CMx1. We now build up pl(t) by its components. Let 7yi represent an unknown complex gain and 01 the direction corresponding to the Ith satellite source. The direction 09, which is always referenced in the antenna frame, is completely specified by two parameters. To illustrate this, consider an antenna located at the origin of the antenna coordinate system. The vector from the antenna to any source can be described in spherical coordinates by a distance and two angular parameters, say r, WI, and W2. Then the two angles W, and W2 comprise 0, and the LOS unit vector v = [vz' Vy v]T in antenna coordinates is found from: cos(WI) cos(W2) v = cos(Wi) sin(p2) (4.3) S sin(Wi) The array response vector (spatial steering vector), a(01), is defined as the response of the array to a signal impinging from direction 01, including both the gain and phase response of the array. Using this definition the array response vector can be written in terms of the components of v and a series of vectors bi, i = 1,..., M from the array reference point to the sensor locations. 2w bTV ej7 1.U e2 bT, e 2 a(o) = (4.4) .2wbT The signal from the Ith satellite, s1(t), contains the sampled spreading waveform di(t), and the sampled Doppler ei'dit referenced at the array reference point, defined as the origin of the antenna frame. si(tn) = di(tn)ejwdt" n = 1,... ,N (4.5) Lastly, 0i1(t) is the low rate (50 Hz for GPS) data sequence that modulates the fast rate spreading waveform. Using these definitions, the vector representation of the signal received from source I at the M sensors is p =(tn) = 7ri(tn)a(Oi)s((tn), n = 1,..., N (4.6) For this analysis, we assume that /31(t) is constant over the period of interest, and include this constant in the unknown gain 7j. If the data bit is not constant but its values are known, they can be included in the waveform sequence si(tn). Equation (4.2) may now be written in a more compact form as L x(tn) = +7 a(0i)sL(tn)+n(tn) 1=1 = A(8)rs(tn) + n(tn) n = 1,...,N (4.7) where E = [01,02,...,OL] A(E) = [a(01),a(02),...,a(L)] r = diag[71,72,...,7L] s(tn) = [s1(tn),s2(tn),...,sL(t)]T (4.8) and the time varying components of the satellite signals at time tn, are contained in s(tn). The signal portion of the received data from each source may be interpreted as the Kronecker product of the Mxl spatial steering vector a(0z) and the lxN temporal steering vector yj for each source (we use yj for the temporal steering vector to minimize confusion with the column vector s(t), which contains information from all sources at time t). yt= [si(ti),si(t2),.. ,sI(tN)] (4.9) This product is the socalled spacetime steering vector v(01, YL) [33], v(01, Y) = [YL(tl)aT(OL), y(t2)aT(Ol),..., yl(tN)aT(OIT)]T (4.10) which may be written concisely as v(01,yj) = a(0) y (4.11) Using this notation allows us to write the MN x 1 received data vector in terms of the spacetime steering vectors, complex gains, and additive interference. L x= 7,v(0i,y) + n (4.12) 1=1 where x and n E CMN x I as below. x = [xT(tl),xT(t2),...,xT(tN)]T (4.13) n = [nT(t),nT(t2), n...,T(tN)]T (4.14) The model for the interference is left relatively unstructured, allowing for thermal noise, multiple jammers, and other interfering signals. We assume that the interference noise waveform vector, n(t), n E CMX1, is zero mean, wide sense stationary, and circularly symmetric complex Gaussian with covariance matrix R,, where the subscript "s" implies the spatial, i.e. sensor to sensor covariance. The interference is uncorrelated from all satellite signals. Moreover, we assume the interference is temporally uncorrelated. Since thermal noise will exist in the sensors' data, it is safe to assume that R.(0) is positive definite. These conditions on the interference may be concisely stated as E[n(t)] = 0 (4.15) E[n(t,)nT(tj)] = OM (4.16) E[n(ti)nH(tj)] R R.(t, t) = R,(0)6t.t, (4.17) R,(0) > 0 (4.18) These assumptions are common for analysis involving jamming or other interfering signals [21,22] and are less restrictive than those of many works on direction finding for communication systems [16,34] which require R, to be a (scaled) identity matrix, i.e. R. = o2I. A final assumption is that the interference covariance is known. In practice, it is not known, but can be estimated. A method of estimating the statistics of the interference is presented in Appendix B. 4.3 Maximum Likelihood Direction Estimation In general, the task is to estimate 6L parameters from the N samples of the data vector x(t), these being two parameters of the direction 0, the real and imaginary com ponents of the complex gain 71, and the appropriate delay and Doppler of the temporal steering vector for each source. However, the problem may be simplified by assuming that accurate delay and Doppler estimates are provided from the receiver code and phase tracking loops, providing Yk as the estimate of the temporal steering vector to source I. This leaves 4L real parameters to estimate from the MN data samples. This estimation begins by examining the likelihood function. In general, the space time likelihood function conditioned on the unknown parameters 8 and r is f(x r) = )RIexp [[x T(,r)]HRt [x T(,r)]] (4.19) where L T(e,r) = yiv(,yi) (4.20) 1=1 The parameters and I that maximize (4.19) are defined to be the maximum likelihood (ML) estimates of 8 and r [35]. The spacetime interference covariance matrix R8t contains the interference covariances between all sensors m = 1,..., M at all times t = 1,..., N, such that the ijth M by M block is E[n(tj)nH(t3)] [36]. By employing the assumption made in (4.17) that the interference is temporally white, the spacetime covariance matrix takes on a block diagonal structure, which serves to decompose the argument of the exponential into the sum of the individual time components. R, 0 ... 0 0 RA 0 0 R.t = (4.21) 0 ... '. 0 0 ... 0 R, Taking the negative of the natural logarithm of (4.19), and dropping the constant terms produces the familiar log likelihood ratio (LLR). N LLR = E[x(t) A(e)ry(t,,)]JHR;l[x(t,) A(e)ry(tn)] (4.22) n=l Since maximizing (4.19) is equivalent to minimizing (4.22), the estimates 1 and t are found as N E,t = arg rmin [x(t) A(e)ry(t)]HR;1l[x(tn) A(e)ry(tn)] (4.23) n=l A significant factor in the estimation of the directions to the satellites (and later, the attitude) is the knowledge of the satellite waveform. Since the waveforms are known, the received data may be despread and integrated, increasing the signal to interference ratio. This is accomplished when the estimated temporal steering vector is used to despread the signal. Demodulation of the Doppler and spreading sequence is achieved by projecting the data from each sensor onto the estimated temporal steering (row) vector ky for each source. This process occurs in the "vector demodulators" of Figure 3.1. As shown in figure 3.1, let ul represent this normalized projection onto the M x N data matrix X, which is formed by reorganizing the MN x 1 received data vector x such that the columns correspond to the time snapshots as shown below. ( 2 u1 1 (4.24) where X = [x(t1),x(t2),..., x(t)] (4.25) "H E, = kikH N = (tn)(tn) (4.26) n=1 Each of these projections are normalized by the energy E, in the lth source. This is a common normalization and is chosen so that the power in the integrated signal snapshots remains constant, while the signal to interference plus noise ratio (SINR) increases linearly with the number of samples integrated. That is, if SINR1(N) is the SINR in the lth satellite component of the data vector ul for an integration of N samples, then SINR1(N + 1) _N+1 (4.27) SINR, (N) N Unless otherwise stated, the remainder of this work assumes that the receiver maintains track on all satellites used for direction or attitude estimation in both time (code) and frequency (Doppler). This allows the substitution k, = y, for 1 = 1,..., L. This assumption is justified because the adaptive beamformer, which follows the satellite channel, mitigates the effects of jammers. We conclude this section by developing concise expressions for the unknown parame ters. To review, the problem is to determine, for each of the L satellites, the 2 components of the direction vector and the 2 components (real and imaginary) of the complex gain. It can be shown [22] that if the signal waveforms are uncorrelated, as is the case with the waveforms employed by the satellites, that this multidimensional estimation problem decomposes to a series of decoupled estimations of individual source parameters. At the stationary point corresponding to the (as yet unknown) estimate of direction 01, the complex gain is found to be aH(i)R,1u( a"(,)R.a(,) (4.28) and the estimator for 01 is lH(O)Rt1u1 2 0= argmax H ) l= ,..,L (4.29) 0, aH (00)R. la(01) " Equation (4.29) provides the means to independently estimate the directions of arrival of each of the satellites in view, and concludes the review of direction finding. 4.4 New Concepts Direction finding, in loose terms, states that the information provided by each source (satellite), through the differences in time or phase of the source signal as received across the sensors, is a direction. However, for the GPS attitude determination task consider the new concept that each satellite directly provides an (ambiguous) estimate of the antenna attitude when the known direction to the satellite in the locallevel frame is included in the estimation process. An Alternate Array Response Vector We begin by reexamining the array response vector a. Previously this was defined as a function of the two parameter valued 0, which was derived from the line of sight vector v (in the antenna frame) to the source, as shown in equation (4.4). Indeed, for the last 30 years or so of array processing the argument of the array response vector has been the angular direction to the source. However, another equally valid method of parameterizing this expression is in terms of the line of sight vector in the local level frame and the antenna attitude. These two LOS vectors (using the subscripts B and LL to denote body frame and local level frame, respectively) are related by the direction cosine matrix that transforms vectors in local level frame to body frame. The direction cosines that comprise the direction cosine matrix are, of course, defined by the body attitude, q. VB = Q(q) ViLL (4.30) The new array response vector parameterized on antenna attitude is then written as ej 2 bTQ(q)VLL ewbTQ(q)vLL a(VLL, q)  (4.31) ej bTQ(q)qLL For the GPS attitude determination application, this method of parameterization is preferred since the line of sight vectors are known in inertial frame, and the body attitude is the desired quantity. It is important to note that when the antenna is actually at the orientation incorporated in the second method, the two array response vectors are identical: a(01) =a(v, q) (4.32) This redefinition of the array response vector is a critical step in the development of the estimator of the next section. Notice that the directions to the sources in the antenna frame no longer appear in the array response vector, and supports the earlier claim that this algorithm is not strictly a direction finding algorithm. However, this new parameterization introduces a complication not found when parameterized by direction: attitude ambiguity. Attitude Ambiguity The attitude ambiguity resulting from the new definition of the array response vector is the manifestation of the fact that attitude cannot be uniquely resolved when using information from only a single satellite source. A useful visualization of this ambiguity is obtained by imagining the LOS vector from the array reference to the satellite source as a "stick," and this stick is "glued" to the array of sensors at the array reference point. With this physical model the family of possible attitudes is obtained by "twirling the stick" while keeping it pointed at the satellite source, as shown in Figure 4.1. Mathematically, this ambiguity can be obtained from the quaternion representation of attitude. Let q represent the true attitude of the antenna. The locus of possible attitudes can be determined by the true attitude, the LOS to the satellite source, and a scalar variable w where 7r < w < 7r. Using these factors, the ambiguity in attitude from the new array response vector definition may be represented as a(q, v)=a((4(w)*q), v) (4.33) where, from equations (2.5) and (4.3), sin(f) E(s)2= (4.34) cos( ) and represents quaternion multiplication. Fortunately, this ambiguity can be completely resolved by using two noncolocated sources. This is again consistent with standard GPS attitude determination theory, which states that the minimum number of satellites required for attitude determination is two [7]. Appendix A proves this fact. Since another common phenomenon in GPS based attitude determination is integer ambiguity, it is worthwhile to take a moment to clearly delineate the difference between the integer and attitude ambiguities. The integer ambiguity arises when the sensors are spaced farther apart than the Nyquist sampling limit. When this occurs, the spatial spectrum is undersampled, and possibly several discrete attitudes would produce the same phase differences across the antenna array. This is the same phenomena as the antenna design term grating lobes. As shown above and in Figure 4.2, the attitude ambiguity in the array response vector is actually a continuum of possible attitudes, not a finite discrete set as in the case of the integer ambiguity. 4.5 Maximum Likelihood Attitude Estimation With the reparameterization of the array response vector, we now proceed with the development of the maximum likelihood attitude estimator. The observables for this algorithm are the vector outputs of the satellite channels, i.e. the demodulated and integrated data vectors ul. There are two reasons this algorithm only considers the demodulated and integrated data vectors. First, the spreading sequence and Doppler Array Reference *" Figure 4.1: Visualization of the possible (ambiguous) attitudes corresponding to a single source. The possible array attitudes can be found by "twirling" the LOS vector, which is assumed affixed to the array, while keeping it pointed at the satellite source. demodulation is typically implemented in hardware [37]. Second, the GPS signal is too weak to evaluate before significant integration. Recall from equations (4.12) and (4.24) that ul, I = 1,2,... L contains contributions from all satellites and interference. We denote the signal contribution from the lth satellite as z, and the interference (including jammers, thermal noise, and the multiple access interference from other satellites) as wi, so that ut = z1 + w1 (4.35) where 1 o H z, = ;a( v, q )ylk" 1 = 71a(vl,q) (4.36) and wi [n(ti),n(t2),...,n(tN)]+ 1 a(vp ,)y] Sq 1 (4.37) El p=i, pi Possible Ambiguous Attitudes 1 50 0 50 100 Roll Angle (deg) Figure 4.2: Possible roll and pitch Euler angles arising from the ambiguity in the new ar ray response vector (yaw angles not shown). The true roll is 25 degrees, and the true pitch is 15 degrees, as indicated by the diamond. The locus of possible attitudes is determined by the LOS to the source and the true attitude. A different source LOS would have a different ambiguity. and q represents the actual (true) attitude of the antenna array. As was assumed pre viously, the signals are assumed in track such that y is essentially yj. It is clear that uL is a consistent estimator of zi. Just as with the thermal noise and jammer signals we consider w, to be a zero mean, circularly symmetric, complex vector of Gaussian random variables.1 E[wl] = 0 [wIwfT] = 0 (4.38) (4.39) Using the these, we define the spacesatellite data vector, U, of the received data as 1 Using this model for the multiple access interference is similar to approaches used for bit error analysis of multiuser communication systems. Ul U2 U U2= (4.40) UL Z and W, the components of U, are similarly defined as Zl Z2 Z = (4.41) ZL Wi W2 W 2= (4.42) WL such that U = Z + W (4.43) Finally, let Rs. denote the ML x ML spacesatellite covariance formed from W as R,, = E [WWH] (4.44) As with the direction finding example earlier, the attitude estimation problem is cast in terms of the log likelihood function. However, now the LLR is a parameterized by the antenna attitude and the complex gains, and as with any likelihood function the true parameter (i.e. the true attitude q) is replaced with a variable (q), the parameter to be estimated: LLR = [U Z(r, q)]H R, [U Z(r, q)] (4.45) Using the results of Appendix B, the ML x ML spacesatellite interference covariance matrix asymptotically reduces to a block diagonal structure of M x M blocks: C1 0 ... 0 0 C2 0 0 R,, = (4.46) 0 ... **. 0 0 ... 0 CL Therefore as N becomes large, the LLR asymptotically approaches the following sum of terms: L LLR ZE[ui ya(vi, q)]gCl [ult 7ya(vi, q)] (4.47) 1=1 Notice that we have reparameterized the LLR in terms of the known satellite LOS directions in the inertial frame, the desired attitude q, and now the unknown complex gains 7y. As with the direction finding application, a closed form expression for 7Y at a station ary point may be obtained by setting the partial derivative of (4.47) with respect to 71 equal to zero and solving for 7y in terms of q. a7 = a(q) '(, q)C, q) (4.48) Substituting this into equation (4.47) produces in equation (4.49) the Maximum Likeli hood Attitude Estimator (MLAE): L S= argmin [ui jja(vi, q)jHC1[ul 7Ia(vi, q)] (4.49) 1=1 One way to view this expression is that it produces a "metric" value at every possible attitude, and the attitude with the lowest metric value is chosen as the estimate. We now present three alternate forms of the MLAE. The first alternate method of expressing the MLAE is found by expanding the summation above: L = min ufCluj j(*a(v, q)HC lu  q= juiCT la(v1, q) + 7171*a(vl, q)HCtla(vi, q) (4.50) and substituting in the expression for the estimate of the complex gain, , from (4.48) to produce =arg max L a(vj'q)jHClu 12 q= argmax (vq)Cla(v, q)4.51 A simpler form of the expression may be written using the shortened form of the array response vector: at = a(v1,q) (4.52) A "whitened" array response vector may be written in shortened form as well: C C l/2a(vz, q) (4.53) where 1/ C21/2 C(4.54) C12 (C 1/2) (4.55) Using these definitions, (4.51) may be rewritten as L 1 c,1/2 )2 . = argmax lu CIl/2(q)2 (4.56) q 1= k(q)12 L axgmax HC1/2p CI1/2U = argmax uCl PI(q)C1/2uL (4.57) /=1 where P&I(q) = (q)A(q) (4.58) Finally, a very compact expression for this estimator for q may be written as L = arg min E TR[PIC1] (4.59) q=1 where PI = [u1 1La(vi, q)] [ut Ita(vi, q)]H (4.60) 4.6 Discussion The most salient observation of this estimator is that the likelihood metric is defined as a function of the antenna (body) attitude. Directions to sources in the antenna frame are not a part of this estimator. Regardless of how parameterized (Euler angles, quaternion, rotation matrix, etc.), the attitude that minimizes (4.49) is the maximum likelihood estimate of the antenna attitude, and since the spacesatellite covariance matrix is (for large sample numbers) block diagonal, the estimator decomposes to a sum of terms, each of which contains contributions from a single satellite. Satellite 1 600 . 500 ..... ...... .... 4 0 0 . . ...... . 300 200 100  * 5 5 Yaw degrees 5 5 Pitch degrees Figure 4.3: Contribution to the value of the maximum likelihood attitude estimator of equation (4.49) from the first of seven satellites in view. This satellite provides little in formation about yaw, and the most in the positive pitch negative yaw to negative pitch  positive yaw dimension. To provide insight into the operation of the estimator, examine these contributions in a neighborhood near the true attitude. When the LOS to an interference source is close to that of the satellite, intuitively the contribution of this satellite to the final solution should be small. Indeed, this is the case with this estimator. In the neighborhood of the true attitude, the contribution to the metric of satellites close to an interference source varies much less than the contribution from those with a large spatial separation from interfering signals. Visually, this causes the metric to appear flatter across attitude, and therefore contribute less to the final attitude estimate. The shape of the likelihood metric contribution from a particular satellite is not uniformly steep or flat across all attitude Satellite 5 600 ...... ... 500.. ..... ' 400 300 100~ ..... 0 0 Yaw degrees 5 5 Pitch degrees Figure 4.4: Contribution to the value of the maximum likelihood attitude estimator from the fifth of seven satellites in view. This satellite provides little information in the positive pitch negative yaw to negative pitch positive yaw dimension, incidentally the dimension satellite in which satellite one provided the most information. dimensions, but varies with the amount of attitude information that satellite provides in each dimension. For example, in a nonjammed scenario a satellite directly above a level antenna array provides little information about the antenna yaw angle. Similarly, the location of jammers may decrease the amount and type of attitude information a satellite may provide. In effect, the estimator is weighting each satellite's contributions to the attitude estimate by the amount and type of information that they are able to provide. This phenomena may be observed in Figures 4.3, 4.4, 4.5, and 4.6. Figures 4.3, 4.4, and 4.5, depict the contribution of three individual satellites to the likelihood "metric" as a function of yaw and pitch. Only the two attitude dimensions yaw and pitch are shown for graphical reasons; for the plots the metric is evaluated at the true roll angle. The true attitude in Euler angles is [0 0 0], i.e. the local level frame. Note that the individual satellite contributions contain regions where they provide high quality information (i.e. a steep slope across a dimension) and regions which contribute little information (little slope or fiat). When the individual contributions are combined, one satellite's weakness may 40 Satellite 7 600 ...... ., ... 5 0 0 ..... . .. .. 400 300 \ 200 100 1 5 05 Yaw degrees 5 5 Pitch degrees Figure 4.5: Contribution to the value of the maximum likelihood attitude estimator from the seventh of seven satellites in view. be compensated by another's strength, producing the steep slope in all dimensions easily observable from Figure 4.6. The nature of the likelihood metric has implementation benefits as well. This separable by satellite structure provides a framework for parallel computations, where for a given attitude under test, the terms of the likelihood for each satellite could be computed simultaneously on multiple processors. In addition, should the receiver lose lock on a particular satellite signal, its contribution to the metric can easily be removed. 4.7 Comments on Searches There are many ways the two dimensional search of (4.29) or the three dimensional search of (4.49) can be implemented. In general, the likelihood surface does not monotoni cally converge to a single global minimum, but may have several local minima. However, if a reasonably accurate initial estimate of the antenna attitude is available, it can be used to provide a starting point for a limited search over an uncertainty region that contains only one minimum. One method of searching this surface is to evaluate (4.29) or (4.49) at several points in a coarse grid spanning the uncertainty region. The point with the largest All Satellites 4000 3000 ........... : 2000. 1000. ...... . 55 00 O, ......... i Yaw degrees 5 5 Pitch degrees Figure 4.6: Total value of the metric including contributions from all satellites in view. Notice that the areas of weak information from any particular satellite have been filled in by the contributions of other satellites. likelihood value is then chosen as the center of a finer grid search, where the fine grid size is chosen from the required angular resolution for the particular application. Another method is to use a variation of the alternating maximization proposed by Ziskind and Wax [14] (for further reference see Press et al. [38]). For a two dimensional direction finding search the process is as follows: First 0. is fixed at some initial estimate and the LLR is maximized with respect to 6y across the uncertainty area. Then Oy is fixed at the value that maximized the likelihood function and the LLR is maximized with respect to 0_. This iterative method alternates movement throughout the two dimensional likelihood space (parameterized by 0, and Oy) along the O0 and Oy axes, and will converge to a local minimum. The extension to a three dimensional attitude search is straightforward. Again, if the uncertainty area is reasonably small and contains the true global maximum, then this method will converge to it. If no a priori information exists to provide an initial estimate of attitude (and therefore no initial estimate of 0 for the search of equation (4.29)), then the problem is significantly more tedious. In this case use of a more exotic method of searching, such as genetic algorithms may be employed [39,40]. CHAPTER 5 PROPERTIES OF THE MAXIMUM LIKELIHOOD ATTITUDE ESTIMATOR 5.1 Introduction Chapter 4 derived a new method of estimating antenna attitude using the demodu lated (despread) data vectors for the satellites in view, the maximum likelihood attitude estimator (MLAE). It showed that by reparamterizing the array response vector in terms of attitude and not antennaframe direction, the estimator decomposed into a series of components, one per satellite source, parameterized over the antenna attitude. This attitude estimator is fundamentally different from the direction finding works [22,23], where each source direction estimation is decoupled from another, even though in both the attitude determination and direction finding cases the source wave forms are known and uncorrelated. However, in this chapter we show that many of the desirable statistical properties of the direction finding estimators apply to the new attitude estimator as well. Specifically, in this chapter we show that the MLAE is asymptotically consistent, unbiased, and efficient. Recall from Chapter 4 that the MLAE attitude estimate (from equation (4.56), which is repeated here for convenience) is the attitude that maximizes a sum of terms obtained by despreading the sensor data matrix with each satellites spreading waveform: L = max uiC PA (q) C T12 ut (5.1) 1=1 where Pj(q) = 5(q)a,(q) (5.2) Now consider the right side of equation (5.1) as the function G(q) of the unknown attitude: L G(q) = uC1/2P(q)C1/2U (5.3) i=1 G(q) is simply the value of the estimator "metric," parameterized on the various hypothe sized attitudes of the GPS antenna array, as discussed in Chapter 4. This chapter is organized as follows. First, using the function G(q) defined above, this chapter proves that the MLAE is a consistent estimator. Consistency is a requirement for the development of the remainder of the chapter. To facilitate proof that the MLAE is unbiased and efficient, a series of lemmas on the asymptotic properties of various aspects of the MLAE are presented. Using these lemmas, this chapter then presents two theorems that show that the MLAE is asymptotically unbiased and statistically efficient. The chapter concludes with a summary and a brief discussion of the validity of asymptotic properties. 5.2 Consistency A consistent estimator has the desirable property that the estimation error decreases as the number of data samples increases. For the MLAE, this corresponds to the number of samples used for integration, N. Theorem 1 The MLAE is a consistent estimator for antenna attitude. That is, 0 lim q = q (5.4) NVoo / Proof: Consider the contribution of only the first satellite to G(q) in equation (5.3), and define this contribution as Gi (q): Gi(q) = uC H1/2P (a)C,1/2U (5.5) 1l "1 "1() u(5) Now from Chapter 4 ul is a consistent estimator of 7ya(v1, q) where q represents the true attitude. So asymptotically, (5.5) will converge to the maximum of aH (v1, q)C1 la(vi, q) or using the whitened array response vector notation, the maximum of ( ( q)152 I1711 2 a(v, q)'dv, q (5.7) where d(vi, q) is defined in (4.53) One could look at equation (5.7) and ascertain, by the CauchySchwarz inequality, that the attitude that maximizes G1 is q. However, this would not be entirely correct. Indeed, q would be one attitude to maximize G1. But since the array response vector is ambiguous when parameterized in attitude, an entire family of attitudes will maximize G1. Let ql(wl) represent the family of attitudes that asymptotically maximize (5.7). Then using (4.33) and (4.34), ql(w1) can be written as q1(w1) = qi() *q (5.8) where represents quaternion multiplication. Now consider the remaining satellites' contributions. In a similar fashion as with satellite one, the family of attitudes q2(w2) that asymptotically minimize the second through Lth satellite's contribution to (4.49) can be found as q.(w.) =4t,(wi)*q i= 2,...,L (5.9) For the L satellites, there now exist L families of attitudes that separately minimize each satellites contribution to (4.49). In general, the parameter value that minimizes a sum of functions is not necessarily the value that minimizes any particular function. However, for the attitude determination case, Appendix C states that each of these families of minima intersect in a single point, the true attitude q. So asymptotically, the estimated attitude is the true attitude, and therefore the MLAE is consistent. U 5.3 Lemmas on the MLAE We begin by considering the function G(q) of equation (5.3) near the estimated atti tude q. A popular method for analyzing G(q), and therefore the estimator performance, is through a Taylor series expansion of the gradient vector g'(q) near 4 (see, for exam ple [16]). The gradient vector is a 3 x 1 vector, since there are three independent attitude parameters (for example, the three Euler angles ,, e, and p). From (5.1), 4 maximizes G(q), and therefore the gradient g'(4) at this stationary point is zero. If the difference between the estimated and true attitude is small,1 then the higher order terms in the Taylor series expansion may be ignored, leaving: 0 = g'() = g'() + g"()( ) + (5.10) such that )= [g()] g'(q) (5.11) where g'(q) aG8q Lq o(q) (5.12) [ ql 90 Gg and ( 82G(g) 92G(q) 82(9) Oq18ql OqiOq92 9OqlOq31 g"(q) A 82G( ) 82G() 12G(q) (5.13) S%2991 8q289q2 89q289q3 82G() 82G(q) 82G(q) 8q389q2 8q38q2 Oq3893 / and qj, q2, and q3 are the three attitude parameters. We now develop four lemmas describing the terms of equation (5.11). Lemma 1 The gradient vector g'(q) of equation (5.12), evaluated at q = q, asymptotically may be written as g'(q)  2 Re ZYf)H(l)P } (5.14) U=l M Proof: This lemma is proven by straightforward evaluation of the terms in (5.12). We begin by evaluating the partial derivative of G(q) with respect to the ith attitude parameter. G() = 2 Re UHcl1/2Pj 1(q)/2() (q)Cl/2l (5.15) U/=l I 9=U where a(q) = (a'(q)Afi(q))1ii'(q) (5.16) 1 The error is small if the integration time (i.e. N) is sufficiently large since the MLAE is a consistent estimator, as previously shown. and we define d(l) = O (q) (5.17) o9qi and da(l) = C1'2di(l) (5.18) and have used the fact that the derivative of the projection matrix is [41] aPp(q) = P ~(q)adi( ) +t (P (q)dai()aa(q)) (5.19) PTq,()a (q 84+ In (5.12) and the remainder of this chapter we use the shortened array response vector notation of (4.52). Using the method of [22,23], the second ul in (5.15) may be replaced by its asymptotic value 71a1(q). G 2 Re{Z uHC1/2P ai(1)t(O)C/2Ya() (5.20)  R e Pu di a (q) C 7t q) (5.20) it q= (OJ q 1=1I = 2Re { uHCI/2P(O)a()(q)7,1(/) } (5.21) = 2Re { uc/2P (0)d,() (5.22) and observing that p() C1/2al(q) is zero, equation (5.22) may be equivalently written as aG((q')) R 2Re) Hu_a/q) R Yi/)) p1Ho di(/)y} (5.23) ^tqa q ilq) = 2 Re { >oCP Idi(l)7y (5.24) where wl, defined in (4.37), may be written as w, = u, 7a(() (5.25) and W C = C/2wi (5.26) Therefore, the ith term of the gradient vector may be asymptotically written as 2 Re Z7'd,'(l)P ) Wj (5.27) Now define the m x 3 matrices D(l) = [di() d2(l) d3(l)] (5.28) and D() = [a(l) d2() ad3(l)] (5.29) Using (5.27) and (5.29), the entire gradient vector g'(g) is found as ,( q) Gq OG(q) O T q[ l 8q2 9q3 J 2Re { ,ZD()P()pwl (5.30) Ufl q) which concludes the proof. Lemma 2 The Hessian matrix g" of equation (5.13) asymptotically may be written as (L g"(q) = 2 Re {I 12DH(l)PIf (l) (5.31) Proof: This lemma is proven by evaluation of the terms in (5.13). We begin by considering the partial derivative of (5.15) with respect to the kth attitude parameter. a2 G(q) LV'HC1/2 [91/5.2 q = q2 Re uL [Pd,(1)ai] C,1/u} (5.32) q=aqk iq aaqk where a] = (i&ij) fH (5.33) Now examine [P:da(l)Al], the term to be differentiated. From the definition of a Pg, this may be rewritten as L [P ()sl] = a [ia,(); P,a()j [ Pa I] (5.34) where we define = di( (5.35) Using the chain rule to evaluate the derivatives of the second term in (5.34), we get = (a.(1)), + a,() (al)( (5.36) aqk where ()' implies differentiation with respect to the kth attitude parameter. Using (5.19) and (5.36), the second term of (5.34) is found as (9P. i 4 +P= a P!, ) 9qk 9qk 9 (P) + (')  t~ t () H ~p L (t = Pdk(l)ajtdi(I)a + dk H &Id a +Pkl (aim), a + P5d.(0 (at), (5.37) k k/ The derivative of the psuedoinverse, a&, is found by direct evaluation to be ((1= (wH1 H = (iiHf)1 d'( I (f)2 (a'(H)4 + 1 Adak()) AI (g'&L)1 (da(j)P if'd(a1)&) (5.38) Substituting (5.38), (5.37), and (5.36) into (5.34) produces the desired partial derivative:  [Pi( j P (d,(a),Al Pdk(l)aidi(l)a(l aqk k ( dt)H d(/)Pidi(/)at + ((Hp)Pd(l)d(l)Hp, +(aat)' Pgdi(l)a'dak(I)at (5.39) which may be consolidated and rewritten more simply as [Pa()it] (i)H a()p ( i( padP(l)a + P t (5.40) Oqk = d l)Pd() + aPl by defining 02 to be 412 (d<(I) dI(lald,()aI + (afa,)' d.(Od(1)(Pa + (a d(l)4dk(l)4 (5.41) Substituting (5.40) into (5.32) produces 092G(q) Z_._ Ret ujCTl2[ r t (4H'H 1~~) pt i =2 Re H 12 it) Hd ()Pd,()at + P ] Cl2u/ } (5.42) 8qOqiO'k Now recall that asymptotically, ul  7^a, (5.43) and therefore ufCt1/2 I (5.44) resulting in T2 being premultiplied by yaHP0 which is of course zero. Therefore, (5.42) asymptotically approaches 82G(q) 2Re{ 'yiI2'2' [1 (it)Hdk(1)P~di(I)a + p1"2] )i (5.45) 2 Re [y,2?[ 1 d()P d,(1i }(5.46) 2 Re i{I Y2af()pad(l)},a (5.47) since = 1. To complete the proof, notice that the Hessian matrix of g"(q) may then qbe concisely written as SL H1 2g"(q) =Re 2 ReI ) H(l)Pid(l) (5.48) where D() is defined in (5.29). Pi (5.47) since fitk = 1. To complete the proof, notice that the Hessian matrix of g"(q) may then be concisely written as fL1 911(q) = 2 Re ^ I^ Pfi, )(1) (5.48) where f)(1) is defined in (5.29). 0 Lemma 3 The expected value E[g'(q)] is asymptotically 0. Proof: By using the result of Lemma 1, taking the expected value of (5.14) asymptotically produces E[g'(g)] 2 Re { Df(l)P(q)E[wL] = 0 (5.49) since, from (4.38), w; is zero mean. This may also be seen by substituting the asymptotic value of ut, yjaj, into (5.15) and observing that 7[*aP is zero. U 1(0 TH Lemma 4 The expected value E g'(q) g'(q)) may be asymptotically written as E [g'() (g'())T] 2 Re {I h21(l)Pf)(l) (5.50) 1 ] U=lJ Proof: Consider the i,j term of the matrix g"(q). Using (5.27) and taking the expected value, this produces E LG(_)9G() Oqi Oq3 aq, aq_,J E 4E Z Re { Pe} R {wIP,(q)d(p)} (5.51) 1=1 P=l l ) 90 Using the identity Re(A)Re(B) = Re(AB + AB*) [23] E [oG(q) OG(/)] L L [Y* ydaH(P (0**pn ( ) E 0 q j =2 EEReE[J^p() pw , +7psd'j(p)P ,(w)wPwI P(o)daj(1)] (5.52) and employing the conditions on the interference w from (4.39), (5.52) may be written as [Oa(q) Oa(q)l L2 E fe a(p (aS3) E aqi ]qj 2 Re l(5.53) 51 Equation (5.53) provides the means to populating the entire 3 x 3 matrix E [g'(0) (g'(q))T]. By using (5.29), this may be written asymptotically as E [g'() (g'(q))T] 2 Re f IfI2 )(l)P (l)} (5.54) which concludes the proof. U 5.4 Bias The bias of an estimator is the difference between the expected value of the estimate and the true value of the parameter being estimated. Clearly, having an unknown bias in an estimate may have deleterious impacts. Theorem 2 The MLAE attitude estimator is asymptotically unbiased, i.e. E [4] q = 0 (5.55) Proof: The bias may be evaluated by taking the expected value of equation (5.11), since E [] q = E q ]. E [(4 q)]  E [ g()] ,g,)] (5.56) Substituting in the asymptotic values from Lemma (2) provides E [( q)] E 2 Re {I 1 2H(l)PLf)(}1) g'(q) (5.57) Re Y 1f H L 11 = =1 )Rel PD()} 1 E [g'(g)] (5.58) Now using the result from Lemma (3) provides the following expression for the asymptotic bias [ fL \*I J [0] = 0 (5.60) which shows that the MLAE, asymptotically, is unbiased. 5.5 Efficiency A measure of the quality of estimator performance is provided by its mean squared error (MSE). The MSE for the attitude estimate provided by the MLAE (i.e. the covari ance matrix of the MLAE estimation error), Pq, is found from P=E q q~qf (5.61) The Cram6rRao Bound, PCR is a lower limit on the covariance matrix of the estimation error for unbiased estimators, and an estimator that achieves the Cram6rRao Bound is said to be statistically efficient [9]. That is, for the parameter vector 0, PS > PCR(O) (5.62) where the equality holds if 9 is an efficient estimate of 9. Theorem 3 The MLAE asymptotically achieves the CramerRao Bound on estimation error, and is therefore statistically efficient. That is, P, PcR(q) (5.63) Proof: The estimation error covariance matrix may be evaluated using (5.11): E[( )Q_ )T] = E [g,,"()l 1 g,() (g,())T [g,,()] 1] (5.64) since g"(q) is symmetric. Using lemma (2), we may replace the terms involving g"(q) by their deterministicc) asymptotic values, represented as g,(g): [g,(q)] (5.65) Ew[(4 0)(4qf [ga(0)]" E [g'(0) g 2)(. where (from lemma 2) L } e'(q) = 2 Re jYj HW^~p: ) (5.66) Lemma (4) has previously addressed the quantity involved in the expectation, allowing us to write E [(a  )(2 ] i [g:(0)]1 [2Re{2DH(1)PDb(1)} [g:(g)] (5.67) Noting that the first two bracketed terms in (5.67) are inverses (with opposite signs), we get the following expression for the asymptotic covariance of the estimation error: Pq E = F L [ 2 Re jx: I^2fHlPiiby) \ (5.68) rom (A.61), the CramerRao Bound on the attitude estimation error is IF ( L 1'1 PcR(q) = 2 Re 7LI2DIH(1)P )D(I) (5.69) Observing that P. and PcR(q) are asymptotically equal completes the proof that the MLAE is asymptotically efficient. U 5.6 Conclusions In this chapter we have proven that the Maximum Likelihood Attitude Estimator (MLAE) is statically consistent, and asymptotically unbiased and efficient (i.e. achieves the CRB). These are three useful properties that serve to ensure applicability of the algorithm for realistic applications. However, it is worthwhile to inquire as to whether an asymptotic (i.e. large sample) property has value for this application. The simple answer for GPS attitude determination is "yes." The GPS signal strength is extremely low, so significant amounts of despreading and integration are required to achieve a useful SNR. That is, GPS is simply not used in a small sample manner. Therefore the asymptotic properties of the MLAE are useful perfor mance indicators, since the expected region of employment is a large sample application. This is not to say that in all situations the MLAE will achieve the CRB, especially in low signal to interference plus noise (SINR) conditions, but that the difference between the MLAE performance and the CRB will decrease as integration time increases. CHAPTER 6 ATTITUDE FROM DIRECTION FINDING 6.1 Introduction In Chapter 4 we derived a jammer resistant maximum likelihood estimator for the antenna attitude, the MLAE, and in Chapter 5 this estimator was shown to be asymptotically unbiased, consistent, and efficient. One drawback to the MLAE is the computational burden, since its implementation involves a three dimensional search across attitude. In this chapter we examine methods that, due to their computational efficiency, may prove useful in a resource constrained tactical environment. As is the case with suboptimal methods, the computational savings comes at the cost of decreased estimator performance. In Chapter 2 attitude determination and direction finding are presented as related fields. Using this relation, this chapter develops attitude estimators that operate on the estimated directions to the satellites. This chapter is organized as follows. Section 6.2 introduces the general approach and form these estimators take. In sections 6.3 and 6.4 two formulations of the estimator are presented. Finally, section 6.5 contains a discussion of the attributes and costs of these algorithms. 6.2 Conceptual Approach The approaches developed in this chapter employ two steps. First, the directions of arrival of each satellite are found. Several methods are presented in section 2.4 that may be used for this purpose. Once the satellite directions in the antenna frame are estimated, the body orientation is then found by minimizing the error between the satellite directions (known in the local level frame) rotated into the antenna frame and the direction estimates. The rotation that minimizes the error is chosen as the antenna attitude. To set up these direction finding based approaches, we define the following matrices of directions where the subscripts B, and LL denote the body frame and local level frame, respectively: A is the matrix of measured LOS directions to the satellites, B is the known matrix of LOS vectors in the local level frame, and W is a weighting matrix. A = [' ...L;L]B (6.1) B = [... VLILL (6.2) W = diag[wi,w2,...W. ] (6.3) where vi, i = 1,2,..., L is a 3 x 1 LOS direction vector. Then using the measured and known directions to the satellites, the attitude found from the following minimization. 4 = argmin lAW Q(q)BWIIF (6.4) q This attitude determination algorithm may be summarized as follows: Step 1 For each of the L satellites, estimate the directions 01, which correspond to the LOS unit vector ;LB. Use an algorithm from Section 2.4 or 4.3. Step 2 Use these LOS vectors to estimate attitude from equation (6.4). Two issues exist in implementing (6.4), once the directions are known. First, the equation is nonlinear in the attitude q, and therefore the attitude cannot be found directly through matrix inversion or psuedoinversion. Second, the choice for the weighting matrix is undefined. We address the first issue through an iterative method based on a linearized version of (6.4). To address the second issue, two options of the weighting matrix are presented. For the remainder of this chapter we parameterize attitude only as a unitnorm quaternion to avoid the convergence issues and singularities that are possible when using Euler angles. 6.3 Equal Satellite Weighting As in Markel et al. [4], one choice for the weighting matrix is the identity matrix. With this choice every satellite contributes equally to the solution, and (6.4) may be written as = argmin 1A Q(q)BIIf (6.5) q Conceptually, the true LOS vectors in body frame may be considered to have been rotated by the true attitude qt. Using this concept we define Z(q) as the matrix of LOS vectors in the local level frame rotated by the true attitude, and Z(q) as the matrix of LOS vectors in the local level frame rotated by the attitude represented by the arbitrary quaternion q. Z(q) = Q(q)B (6.6) Z(q) = Q(q)B (6.7) Now A, the matrix of measured LOS vectors in body frame, may be considered to be an estimate of Z(q), say Z(q), since the directions are measured in the body frame and known in the local level frame. In order to extract the attitude from the measured directions in A, i.e. to estimate q" from the estimate of Z(q), substitute Z(q) (i.e. A) and Z(q) into (6.5). 4 arg min j [Z(q) Z(q) [F (6.8) In the nonoise case it is easy to see that the quaternion that minimizes (6.8) is q = q, since Z(q) will equal Z(q). This quaternion 4 of (6.8) may be found through an iterative process. This procedure follows the general form of Cohen [7] (which was solving for attitude in Euler angles by using phase differences), but has been rewritten for this application. We begin by expressing Z(q) through its Taylor series expansion. 4 Z(t) = Z(4) + Z Z(/l (q 4) + {higher order terms} (6.9) ,= qi Ignoring the higher order terms, (6.9) may be written as E = X q (6.10) where E = vec(Z(q) Z(4) (6.11) X = [vec(O ), vec(OZ()), vec(0 ),q3 vec( )]K  (6.12) 6q = q4 (6.13) Jql Sq2 (6.14) Sq3 6q4 and vec(.) represents stacking all columns of a matrix into a single column vector. The quaternion update vector is found from the psuedoinverse of X. 6q = (XHX) XHE (6.15) The iterative algorithm for solving (6.8) may be summarized as follows: Step 1 Choose an initial value for the estimated attitude quaternion 4. This may be the last attitude estimated from the previous data set, or if no a priori information exists, the unit quaternion. Step 2 Calculate the quaternion update Sq from (6.15). Step 3 Calculate the new quaternion by adding 4 to 6q. Step 4 Normalize the new quaternion to unit norm. Step 5 Return to step 2 until convergence occurs, i.e. when 6q is arbitrarily small. It is important to include Step 4, since the above method does not intrinsically ensure that the quaternion remain unit norm. This algorithm will be referred to in this work as DFU, since it uses direction finding to obtain the attitude estimate, and does not use an weighting matrix that is the identity matrix (i.e. is unweighted). 6.4 Satellite Weighting via Adapted SINR The algorithm developed in this section uses the same iterative approach as the DFU algorithm, however now the weighting matrix is not the identity matrix. The choice employed here is to weight each satellite's contribution to the attitude estimate by its adapted SINR. The adapted SINR is found from the array response vector to the source, the interference covariance matrix, and the complex gain [42].1 SINRAdapted = hj2aH (0)Ca(0) (6.16) where 7y is the complex gain, a(0) is the array response vector from the source located at direction 0, and C is the interference covariance matrix. In practice, these values may be replaced with their estimates. Throughout this chapter the method of determining the directions to the satellites in the antenna frame has been left unspecified. However, if the maximum likelihood algorithms of (4.28) and (4.29) of Section 4.3 are employed to estimate the directions to the satellites, then calculation of the adaptive SINR is essentially provided by the algorithm already. Notice that the adaptive SINR can be viewed in terms of the whitened array response vector d. ii = Ca (6.17) and therefore SINRAdapted = H[Y2jij2 (6.18) The greater the attenuation from the whitening, the lower the adapted SINR. We will refer to this algorithm as DFW, since it uses direction finding and a weighting matrix. 1 The term adapted SINR is actually not used in Applebaum [42], but is a common term today used to represent the signal to interference plus noise ratio after application of the optimal beamformer weights derived in Applebaum [42]. 6.5 Conclusions Three facets of this method of attitude estimation warrant further discussion. The first is the relationship among the satellites. In a multiuser wireless communications system where the individual sources are essentially unconstrained in their possible location, determining the directions to each source separately is acceptable, and in most cases desirable as it reduces the dimensionality of the search. However, the GPS source (satellite) locations are constrained in location by their known flight paths. The constellation, in the locallevel frame, is known by the user receiver, as is the receiver position, providing additional information that could be used in the direction estimation process. Indeed, the direction to one satellite is not completely independent from the directions to all others. Therefore, to estimate the directions independently, as in (4.29), is a computationally appealing but suboptimal method of attitude estimation. In Chapter 8, simulation results quantify the degradation of performance of these suboptimal algorithms from the optimal MLAE. The second area of discussion centers around the contribution of each satellite. In section 4.6 a graphical example was used to show that the contribution of each satellite is not the same across attitudespace. Inherent in the metric of equation (4.59) is the (optimal) mechanization of how each satellite contributes across attitude space to the attitude estimate, even in the presence of interference and jammers. With the LOS mapping approach, this optimality is lost. Using equal weighting, each satellite not only contributes equally across attitudespace, they all contribute equally relative to each other. With the adapted SINR weighting satellites near the jammer sources will not be allowed to contribute as much to the attitude estimate as those with larger angular separations. However, the contribution from any satellite is still equal across attitudespace, since only a scalar weighting is applied to the satellite. One can clearly see how this is a suboptimal approach by trying to represent the information of Figures 4.3, 4.4, and 4.5 with a single scalar per figure. The third discussion topic is the relative performance of these two new attitude estimators. In general, weighting the satellite contributions by the adapted SINR in the attitude minimization process outperforms equal satellite weighting. As shown in Chapter 8, this is often the case for the scenarios simulated. However, in some of the cases examined the weighted estimator has a slightly larger mean total error than the unweighted. This seems to occur when the adapted SINR's are below approximately 610 dB, which occurs at the fastest update rates (i.e. smallest integration times), with the fewest sensors, and smallest baseline lengths. It seems that in this case, the weighting puts too great a confidence in the noisy measurements of a few satellite directions. Although not investigated in this work, one possible improvement to this process would be to use a switchable weighting scheme based on the estimated adapted SINR. As discussed in Chapter 9, this scheme would use the weights only when the adapted SINR's crossed an empirically derived threshold. CHAPTER 7 WIDER BASELINES AND DUAL FREQUENCY USE 7.1 Introduction It is well known that as the separation between sensors increases the variance of the direction estimation error decreases! However, this increase in sensor spacing comes at a cost. If the sensors are spaced farther than the Nyquist criterion, grating lobes appear, resulting in an ambiguity of direction. A similar phenomena arises for the attitude determination application. With the conventional AD approach of using phase differences, only the fractional portion of the total phase may be measured. However, the total phase difference between any two sensors receiving a satellite signal is composed of the fractional part (i.e. the measured portion), and the integer portion. Since the integer portion is not measured, yet must be accounted for to properly determine the antenna attitude, it is referred to as the integer ambiguity by the attitude determination community. Depending on the sensor spacing and the number of satellites used to determine attitude, the integer ambiguity may be resolved uniquely. However, this typically requires an intensive search of all the possible ambiguities. One method to eliminate many of the ambiguities is through a process called "widelaning." The widelane phase measurement is formed from the difference of the phase measurements from the two GPS frequencies, i.e. the Ll and L2 frequencies: 4widelane = 4Ll (L2 (7.1) 1 This fact is easily demonstrated by the CRB of the direction estimate v, for a two sensor array: ao2, > 1/[SNR ir2d2], where d is the distance between the sensors. This gives an effective wavelength of [37] Awidelane = 7LIAL2 (7.2) ALl AL2 which for the two GPS frequencies of 1.228 and 1.575 GHz corresponds to .86 meters, or approximately four times the wavelength corresponding to the highest GPS frequency. This significantly reduces the size of the integer ambiguity search. The price paid for the increased effective wavelength is a doubling of the receiver noise, since the receiver noise in each frequency channel is uncorrelated. The Maximum Likelihood Attitude Estimator (MLAE) of Chapter 4 is, in many cases, unaffected by the integer ambiguity phenomenon. Since the MLAE searches attitude with a fixed array configuration, it can often perform unambiguously even when the sensors are farther apart than the Nyquist criterion. Provided that the uncertainty in attitude over which to search is not too large, the MLAE will converge to a single local minima within the uncertainty when several satellites are being tracked, even if the sensor spacing is greater than the Nyquist criterion. This occurs because inclusion of the additional satellites (over the minimum two required) causes the MLAE metric at the bogus attitudes to indicate that these attitudes are less likely than the true attitude. If the MLAE is cast as a minimization as in equation (4.49), the metric value at the bogus attitudes will not typically be as small as the value at the true attitude. However, a situation similar to the integer ambiguity problem does exist for the MLAE when the number of satellites being tracked is small, the sensor spacing is greater than the Nyquist criterion, and the attitude uncertainty (i.e. the region to search) is large. In this case the MLAE metric at several regions of attitude space may be close to the value at the true attitude, resulting in local minima of (4.49) that could result in false solutions to a minimization search. Therefore, it would be desirable to develop a dual frequency version of the MLAE developed in chapter 4 that reduces the number of false attitude solutions like the conventional attitude determination process of widelaning. Unfortunately, a simple "phase difference" approach like widelaning is not applicable to the MLAE. Since the MLAE does not operate on phases alone but on the complex values of the data vectors, a more elaborate estimation process is required. In this chapter we will develop a method of incorporating the second GPS frequency into the MLAE. We will show that this dual frequency extension provides two positive attributes: a decrease in estimation error and a reduction in possible false attitudes. The following section derives the dual frequency MLAE. Similar to the initial MLAE derivation, this begins with a new look at the array response vector. Following the derivation we present a section that examines the reduction in false attitudes by using the second GPS frequency. Finally, we show that this dual frequency attitude estimator asymptotically achieves the Cram6rRao Bound, and produces higher quality estimates than the single frequency version. 7.2 Dual Frequency Maximum Likelihood Attitude Determination Before this section develops the dual frequency attitude estimator, first review the single frequency attitude estimator of Chapter 4. This estimator used the L demodulated data vectors (each M x 1), the L interference covariance matrices, and an array response vector parameterized on the antenna attitude (and using the known directions to the satellites in the local level frame) to define the log likelihood surface. The minima or maxima (depending on whether equation (4.49) or (4.51) are used) of this function is chosen to be the attitude estimate. Although the development of this estimator is performed using baseband signals, the carrier frequency of the GPS signals appears in the array response vector of equation (4.31) through the division by the wavelength A. Assuming that the conditions on the interference from Chapter 4 are met for both GPS frequencies, the derivation of this estimator will apply equally well to either frequency, since the statistical properties of the DSSS satellite waveforms are the same for both the Ll and L2 GPS frequency. Without introducing any additional unknowns into the estimator, the array response vector may be written as a function of the unknown attitude q, the known satellite direction in the local level frame VLL, and now the known wavelength A: eJ 'bTQ(q)VnLL e' bTQ(q)LLL ei r 2 LL a(VLL,A,q)= (7.3) 'bT' Q(q)vLL Now consider a situation that makes use of this modified array response vector. A receiver capable of demodulating both frequencies of the GPS P(Y) code can generate 2L data vectors, L per frequency for each integration period. To distinguish the data vectors gathered from the different frequencies, the notation representing the data vector obtained from demodulating with the spreading sequence from satellite I is expanded from u, to u(lf), where f denotes which of the two possible carrier frequencies. Using this expanded notation, we define the augmented spacesatellite data vector, Ua, of the received data similar to equation (4.40): u(i,1) U(2,1) U(L,1) U.a = (7.4) U(1,2) U(2,2) U(L,2) Letting U, represent the spacesatellite data vector for frequency i, i = 1, 2, equation (7.4) may be written in block form as Ua= u (7.5) U2 Now consider the statistics of the interference. Similar to equation (7.5), the aug mented interference vector (containing thermal noise, jammers, and the multiple access interference from other GPS satellites, as in Chapter 4) may be written in block form as Wi Wa = (7.6) W2 where W, represents the ML x linterference vector for frequency i, i = 1,2. Assuming that the jammer waveforms are uncorrelated frequency to frequency, Rasa, the covari ance of Wa, is block diagonal just as R,,, the covariance of W for a single frequency in equation (4.46) is block diagonal. This implies that the entire interference at one frequency/satellite channel is uncorrelated from the interference in all other frequen cies/satellite channels. Letting C1,1 represent the covariance of the interference in the Ith satellite channel, I = 1,2,..., L and the fth frequency, i = 1,2, the augmented interference covariance matrix may be written as R,,i 0 Rssa = (7.7) 0 R..2 C1,1 0 *.. 0 0 C2,1 0 0 0 ... ** 0 o ... 0 CL,( (7.8) Ci,2 0 ... 0 0 C2,2 0 0 0 ... '.. 0 0 *..**. 0 CL,2 Using the above definitions, the asymptotic log likelihood ratio of (4.47) may be written in dualfrequency form as L LLR 7 Z[u(,1) (;,i)a('I, AI, q)]HC')[u(,1)  (t,1)a(/, A,, q)] + i=1 L [U(1,2) 7(1,2)a(vI, A2, q)]HC2)[u(2) 7,2a(vi, A2, q)] (7.9) (=1 As with the single frequency attitude estimator, a closed form expression for 7(1,f) at a stationary point may be obtained by setting the partial derivative of (7.9) with respect to 7(,f) equal to zero and solving for 7(Ij) in terms of q. a H(vi, A f, q)C'l)u~j TW ) = 1/7 f, i , (7.10) SaH (v, Af q)Cf)a(vL, Af q) (7.10) Substituting in the complex gain expression of (7.10) provides the following minimiza tion for the attitude q: L S= argmin Z[u(L,1) (L,1)a(m, A,, q)]HC(1ll)[u(,1) y(,1)a(Lt,A1, q)] + 1=1 L E[U(,2) (1,2)a(v(, A2, q)]Hc12)[u(,2) 7,2a(v, A2, q)J (7.11) 1=1 By expanding and removing terms that do not vary with attitude, and substituting in the expression for 7(yj) from equation (7.10), equation (7.11) may be written as a maximization: x" laaA^ Cvi), A2 *' q= argmax f [a2 2q_] + a I (1,^2)U(1,2)1 1 (7.12) aH ,A2,q)C.)a(v, A2, q) J 7 As with the single frequency MLAE, a simpler form of the expression may be written using the shortened form of the array response vector: a(jf)(q) = a(v,A1/,q) (7.13) And as in (4.53) we may define the "whitened" dualfrequency attitude array response vector: a(Q,f) = C(t1) a(vt, A1, q) (7.14) where COX C1/2 0 1/2 (7.15) CV (lf) (lf) (1f)5) C1/2 1C2 H (7.16) (1,1) = \ "' (,f)) Using these definitions, (7.12) may be rewritten as L 2 lU(H C(/2 SL 211f U (11) f2&(Ilf) (q) j' (7.17) argmax (717) q 1=1 f= (,,)(q L 2 = max ZE u )C (1f) .(q) ( 1f)2 u(.) (7.18) 1=1 f=1 where () a(,)(q) a(,f)(q) Paj()q) R ( 1) (q) (7.19) The attitude estimators of equations (7.11), (7.12), and (7.31) will be referred to in this work as the dualfrequency MLAE. Several observations may be made about the dualfrequency MLAE. First, as with the single frequency version, an expression for the complex gain may be found using data from a single satellite and single frequency. However, the attitude estimate includes information from all satellites and both frequencies. Although coupled, this could be implemented in a parallel structure, having perhaps one computer calculating the metric for a set of attitudes using the data from the first frequency while another processes the data from the second frequency. The metric values could then be combined for the set of attitudes investigated, and the attitude with the minimum combined metric value would be chosen as the estimate. Although the metric was derived using the same satellites on both bands, this in not a requirement. In fact, none of the satellites need to appear at both frequencies for the algorithm to operate. Similarly, this approach performs if the jammer scenario is different on the two frequencies. For example, if one band was jammed and the other unjammed, this algorithm would perform better than a single frequency MLAE operating at either the jammed or even the unjammed frequency. 7.3 Reduction in False Attitudes Consider the situation where only a few satellites are visible, a large uncertainty in attitude exists, and to achieve higher resolution attitude estimates the sensors are spaced farther apart than a half wavelength. In practice, this scenario may occur either on startup or if INS aiding is either not available or not operating. In this case, there may be several areas of attitude space that lead to local minima of equation (7.11) (or maxima of equations (7.12) and (7.31)). To find the attitude estimate, the extrema of each area would have to be found and sorted, with the attitude estimate chosen as the attitude corresponding to the global extrema. It is clear that a reduction in the number of areas that may contain an extrema (perhaps crossing a threshold and therefore requiring that the extrema be found and its metric value obtained for comparison to others) would be desirable. To visualize how incorporation of the second frequency reduces the number of false attitude positions, consider the simpler direction finding application. Figure 7.1 shows the (normalized) likelihood ratio for a received signal impinging on the array from broadside for a 15 element uniform linear array (ULA) with a 2A spacing between elements. In this example no jammers exist, i.e. the likelihood value is proportional to the quiescent array pattern. Notice that several directions produce likelihood values equal to the maximum at the correct direction. For a radar these grating lobes cause two deleterious effects: they produce an ambiguity in the direction estimation process, and they allow undesired signals from these ambiguous directions to pass through the antenna with the same gain as the desired signal. Now consider the likelihood value for the same array receiving signals at a carrier frequency equal to 1.25 times greater. From Figure 7.2 it can be seen that the locations of the grating lobes have changed. The extension of this example to the dual frequency MLAE is obtained by summing the likelihood values obtained from these two frequencies, as shown in Figure 7.3. In this case the grating lobes are attenuated, since only the true direction appears at the same angular location as frequency changes. Array Response for Frequency 1 I IM I' 40 20 0 20 40 Degrees from Antenna Boresight Figure 7.1: Normalized likelihood value vs. direction of arrival for a 15 element ULA re ceiving a narrowband signal at a frequency corresponding to a 2A spacing. True direction is on boresight, i.e. 0 degrees. Array Response for Frequency 2 80 60 40 20 0 20 40 60 80 Degrees from Antenna Boresight Figure 7.2: Normalized likelihood value vs. direction of arrival for a 15 element ULA receiving a narrowband signal at a frequency corresponding to a 2 A spacing. True direc 1.25tion is on foresight, i.e. 0 degrees. tion is on boresight, i.e. 0 degrees. Combined Array Response for Both Frequencies '21 0 0 Z 80 60 40 20 0 20 40 60 80 Degrees from Antenna Boresight Figure 7.3: Likelihood value vs. direction of arrival for a 15 element ULA incorporating data collected at two frequencies corresponding to 2A and 2 A spacing. True direction is on boresight, i.e. 0 degrees. The grating lobes are reduced in amplitude from the true direction of arrival, producing a clear indication which direction is correct. Now consider how this applies to the attitude determination case. In the following example, the true attitude (in Euler angles) is [0,0,0], and the attitude uncertainty is assumed to be 88 degrees in each coordinate. Three likelihood surfaces are generated at a two degree resolution: the MLAE (equation (4.49)) using frequency Ll the MLAE (equation (4.49)) using frequency L2, and the dual frequency MLAE (7.11). Each of these estimators chooses the attitude that minimizes the likelihood expression, and each likelihood surface is calculated using two satellites. A threshold of 10dB below the mean value of each surface was chosen, and all points in attitude space below this threshold were considered possible candidates for the attitude estimate. In order to show the reduction in false values of attitude, any point that had a contiguous connection to the true attitude was removed in the following plots. Any point with such a connection would lead, by a simple minimization scheme, to the true attitude. Therefore, for comparison purposes we wish to examine only those points that would lead to false values of attitude, i.e. a local minimum that is not at the true attitude. The intent of this is to evaluate the potential decrease in false attitudes when using the dualfrequency MLAE over using a single frequency MLAE. Figure 7.4 shows a three dimensional representation (over the three Euler coordinates of attitude) of the attitudes where the likelihood value was below the threshold using the Ll frequency. To provide an easier visualization of these points, Figure 7.5 shows a two dimensional version of the same plot, obtained by projecting all points onto the O = 0 plane. As discussed above, the attitude values that connect to the true attitude are removed, leaving only points that would lead to a a false attitude. To illustrate this, figure 7.5 shows all attitudes with metric values below the threshold, and figure 7.6 shows only those that do not connect to the true attitude. The relatively few points leading to the true attitude have been "thinned out." The jammer scenario used is the same as that involving the first jammer in the random jammer study of Chapter 8. The four sensors used are spaced in a "quad" formation with spacing of approximately .6 meters. Figure 7.7 presents the false attitudes for the MLAE using the L2 frequency in the same manner. In Figure 7.8, which contains the likelihood points below the threshold for the dualfrequency MLAE, it is apparent that a significant reduction in the number of points leading to false attitude solutions exists. Similar to the grating lobe reduction in the direction finding application above, the dualfrequency MLAE "smoothes" the likelihood value away from the true value of attitude, reducing the number of points with likelihood values that may be considered candidates for the global minimum. To illustrate the reduction in false attitudes, three additional data sets are presented. Figures 7.9, 7.10, and 7.11 present the single frequency MLAE for LI and L2, and the dualfrequency MLAE using satellites 1 and 5. Similarly, figures 7.12, 7.13, and 7.14 present data using satellites 1 and 6, and figures 7.15, 7.16, and 7.17 present data using satellites 2 and 3. Table 7.1 contains the number of false attitude possibilities for the three estimators for several more satellite combinations. Raw Points for Satellites 1 and 4. Frequency 1 50 ""5 0. 100 ^<^ 550^^ ^ 50 100 0 50 50 50 0 Pitch (Degrees) 100 100 Roll (Degrees) Figure 7.4: Three dimensional plot of all attitudes whose metric value is below the 10dB threshold, using Satellites 1 and 4, and frequency Li. Raw Points for Satellites 1 and 4. Frequency 1 .1 i. . 20 0 20 Roll (Degrees) Figure 7.5: Two dimensional view of all attitudes whose metric value is below the 10dB threshold, using Satellites 1 and 4, and frequency L1. Thinned Points for Satellites 1 and 4. Frequency 1 lv o * : H It t ." :'1 . ................ 80 60 40 20 0 20 40 Roll (Degrees) 60 80 Figure 7.6: Two dimensional view of those attitudes whose metric value is below the 10dB threshold that are not contiguous to the true attitude, (i.e. false solutions) using Satellites 1 and 4, and frequency L1. Satellites 1 and 4. Frequency 2 80 60 40 20 0 20 40 60 80 Roll (Degrees) Figure 7.7: False attitude solutions using Satellites 1 and 4, and frequency L2. Satellites 1 and 4. Both Frequencies 8 0 ... .. ... .. . .. .. . i . .. .. . .. .... .. 4 0 .. . . . * .. .. . i . . . . . . :. . I i .. i . . . .i . . . . . . . . . . i *'lj il J ,2 ,.. ... .. .. ... .... ... ... . ............ + .. .. ... ...... . , .. 40 : : .. " a . . . . i . .. . . . .. . . . .. . . . . . . . . .. + . . . . . . i . . .P. S 26 0 .. .. .. , .. ....... .. .. .. .... .. . ... ... .... .. .........I.. " 0 ... ... .^ ".4 6 . . .. .. ... ; .. . . ; . . .. ; ,. ... .; '" . . . . .. .: .. .. . .. .. .. . 80 60 40 20 0 20 40 60 80 Roll (Degrees) Figure 7.8: False attitude solutions using Satellites 1 and 4, and frequency the dual fre quency MLAE 75 Satellites 1 and 5. Frequency 1 80 ....... ! I ..... : : I... in : I 2 0 ^^ i J ,^ 40 I I : iiijpii: 11 60 .... .I1, .A t l 80 IIjAh "1, "'i ,iJl *II u Ih: h 11 0 ... H .... HH i 1 ____ . 80 60 40 20 0 20 40 60 80 Roll (Degrees) Figure 7.9: False attitude solutions using Satellites 1 and 5, and frequency Li. Satellites 1 and 5. Frequency 2 40 N^ ilirf 80 jii,, 'Iiir "* *i ii *'ll jI!'jII*"ii'i 11 B~I 2. ..... 2 0 .... Il ..... o i l ,.,irl,,,...iSn ,i 0, .... .... ..... f\ .2 It ,x 7{***t t" 2 A .... .. z0 ..il,..:t ]!l : ::"::: t .. .............. i ........ ......... .. 80 60 40 20 0 20 40 60 80 Roll (Degrees) Figure 7.10: False attitude solutions using Satellites 1 and 5, and frequency L2. Satellites 1 and 5. Both Frequencies .1 . . .tl:" if I *Ijf ill :1: 80 60 40 20 0 20 40 60 80 Roll (Degrees) Figure 7.11: False attitude solutions using Satellites 1 and 5, and the dual frequency MLAE. Satellites 1 and 6. Frequency 1 80 60 40 20 0 20 40 60 80 Roll (Degrees) Figure 7.12: False attitude solutions using Satellites 1 and 6, and frequency Li. Satellites 1 and 6. Frequency 2 'HI tt"1,".t i I' .i .. I .. ,h ....siI 1 Is.b .. It l2I l l .j :ia .... ..S..H~f lIi ..".. 20 .. .... i .0 ....1 _2 I l ;:,,, ... .... . ~j ,, ;,^His 'O ... . ....... ..! ... 80 ... .I...... 80 60 40 20 0 20 Roll (Degrees) Figure 7.13: False attitude solutions using Satellites 40 60 80 1 and 5, and frequency L2. 80 60 40 r 20 0a iE20 40 60 80 Satellites 1 and 6. Both Frequencies i.. . . 11~ 80 60 40 20 0 20 40 60 80 Roll (Degrees) Figure 7.14: False attitude solutions using Satellites 1 and 5, and the dual frequency MLAE. ! ...... .... .. . .i. . i i I I I I Satellites 2 and 3. Frequency 1 80 60 40 20 0 20 40 60 80 Roll (Degrees) Figure 7.15: False attitude solutions using Satellites 2 and 3, and frequency Li. Satellites 2 and 3. Frequency 2 8 0 :: . ... . . i t i .. . ... ... . . . t i ,,.., , 60 j(IIin t I, 60 li ..... .. .. .. 4 0 ..... In k;. .. .... ...... ji.. ... ... ... .. ....... S. I h,:ot S 2 0 . i It 1 ... ..! i ...... .. . ..... ..... ... .. . .. . 20t JJJ *1: W i .' .' 3i ii. .i 40ga o 1 j 20I > ^ t . ,I i ... ., i . 0 2 .. .. !. .I : ...... .. ... . .. t .. ..... .. .. .. .' . ; ... .. .. ..! .... i .. . . 60 .. ;,, i itl: I 8 0 ..... .... ... .......l .. .. ..... ..... l i l i ," _! "^ ii 80 60 40 20 0 20 40 60 80 Roll (Degrees) Figure 7.16: False attitude solutions using Satellites 2 and 3, and frequency L2. Satellites 2 and 3. Both Frequencies 20 : 40 . I I I 80 60 40 20 0 20 Roll (Degrees) 40 60 80 Figure 7.17: MLAE. False attitude solutions using Satellites 2 and 3, and the dual frequency Table 7.1: Number of false solutions, i.e. points below the "possible minimum" threshold that do not converge to the true solution. Resolution in each Euler angle is two degrees. Scenario is the same as that involving the first jammer in the random jammer study of Chapter 8. Satellite Combination MLAE LI MLAE L2 Dual Freq. MLAE 14 3681 5052 683 15 3389 4690 153 16 5028 3313 55 23 6169 1409 28 24 733 4670 523 25 3875 1743 89 26 4358 1826 25 34 12095 3667 159 35 1772 1654 9 36 1 463 1788 01 7.4 DualFrequency MLAE Performance Increase In this section we analytically examine the performance of the dualfrequency MLAE. Specifically, we will show that the dualfrequency estimator may be written in the form similar to that used for the single frequency estimator. Then, using the results of Chapter 5, we show that the dualfrequency MLAE has the same desirable properties as the single frequency MLAE. We begin be examining the dualfrequency MLAE from equation (7.11), repeated here for convenience. L = argmin [u((,I) (lj)a(v, A,, q)]HC)[u(L,j) i(j,1)a(vj, A1,,q)] + q 1=1 L Z[U(L,2) "(1,2)a(vL, A2, q)](C')u(,2) 1,2a(v, A2,q)] (7.20) 1=1 This equation may be rewritten in simpler notation as 2L q = argmin [ri( (q)]H t [f 7Ck(q)] (7.21) q where we have defined ( a(v,Ai,,q) 1 < C< L (7.22) a(vL,A2, q) L+1 C< < 2L 1 <= ) (7.23) U(L,2) L+1 c= (7.24) (CL,2) L+1 j ( 1)< ( (7.25) [C(CL,I) L + 1 < < 2L With this substitution, the frequency is no longer an argument of the array response vector, the data vector, or the covariance matrices. The terms relating to the second frequency are simply indexed as the last L C's. With this notation, we can rewrite the maximization equation (7.12) as 2L [~qH1j2 x=argmax y A "I(q)H6 C l 12 (7.26) q l (.(q)H26)(q) As in Chapter 4, we can simplify further by defining the whitened" dualfrequency attitude array response vector: k (q) = 12k (q) (7.27) where Cc = 1/2" 1/2 (7.28) 12 = (C12H (7.29) Using these definitions, (7.12) may be rewritten as S2L If l1/2()J2 S= argmax 1 (7.30) q = jk(q)2 2L max i(C (P q)e"2 (7.31) < =1 where Pk(q) = (q) 4(H/(q) (7.32) P (H) = ( a(q) Notice that equation (7.31) is in exactly the same form as the single frequency version of equation (4.56). With this form we can use the same asymptotic formulas for efficiency as in Chapter 5 and Appendix A. The only difference is that caution must be used in calculating the derivative matrices, to ensure that the appropriate effective array spacing is used depending on whether ( is less than or greater than L. Using the results of Chapter 5 and Appendix A, the CramerRao Bound on the attitude estimation error using two frequencies is P 2L 1 ]1 PcR(q) = 2 Re I k (7.33) I =l1J where, as with the single frequency, ai o6q) (7.34) dqi d() = (12i() (7.35) and 50 [d = )dOa2( ) 3(] (7.36) To justify the observation made earlier that the dualfrequency MLAE performed as well or better than the single frequency MLAE using data from either frequency, consider that the Fisher information matrix for the dualfrequency, FIMDF is the sum of FIM for each frequency. That is, f2L1 FIMDF = 2Re { L I2I5H(C)P (q)6(0) = 2Re IYc21H(()P )(q) l( + 2 Re { _7(c2H(C)P&q)I5(C)} (7.37) 1 =L+I = FIML1 +FIML2 (7.38) where FIMLI and FIML2 are the Fisher information matrices for the two single frequency estimations. Since FIML1, and FIML2 are each always positive definite, FIMDF must be positive definite and greater2 than FIMLI and FIML2. 7.5 Summary In this chapter we have derived the dualfrequency MLAE. This estimator was shown to optimally incorporate information from the two frequencies that carry the GPS P(Y) wavefrom. This estimator has the same desirable properties of consistency and asymptotic efficiency as were demonstrated for the single frequency MLAE. By examining the Cram6rRao Bound for the dualfrequency MLAE, it is clear that it will perform as well as 2 That is, FIMDF FIMLI > 0 and FIMDF FIML2 > 0. 83 or better than the single frequency MLAE using data from either of the GPS frequencies. In addition, the dualfrequency MLAE significantly reduces the number of false attitudes (attitudes with metric values numerically close to the value at the true attitude) when the sensors are spaced farther apart than the Nyquist criterion. This gives it a "widelanelike" quality, and allows for better quality attitude estimates without introducing ambiguity. 