Interference mitigation for GPS based attitude determination

MISSING IMAGE

Material Information

Title:
Interference mitigation for GPS based attitude determination
Physical Description:
169 leaves : ill. ; 29 cm.
Language:
English
Creator:
Markel, Matthew David, 1968-
Publication Date:

Subjects

Subjects / Keywords:
Global Positioning System   ( lcsh )
Artificial satellites -- Attitude control systems   ( lcsh )
Adaptive antennas   ( lcsh )
Electrical and Computer Engineering thesis, Ph. D   ( lcsh )
Dissertations, Academic -- Electrical and Computer Engineering -- UF   ( lcsh )
Genre:
bibliography   ( marcgt )
theses   ( marcgt )
non-fiction   ( marcgt )

Notes

Thesis:
Thesis (Ph. D.)--University of Florida, 2002.
Bibliography:
Includes bibliographical references (leaves 164-167).
General Note:
Printout.
General Note:
Vita.
Statement of Responsibility:
by Matthew David Markel.

Record Information

Source Institution:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 028821579
oclc - 50390937
System ID:
AA00022669:00001

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 anti-jam 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. Cramer-rao 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 co-advisors. 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 off-campus 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 ANTI-JAM ATTITUDE DETERMINATION RECEIVER ....... 18

3.1 Introduction ...... ... .. ....... . ... ........ 18
3.2 A GPS Receiver for Anti-Jam 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 Dual-Frequency 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 CRAMtR-RAO 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 Sine-space 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 (+), DF-W (diamond with solid line),
DF-U (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 (+), DF-W (diamond with solid line),
DF-U (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 (+), DF-W (diamond with solid line),
DF-U (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 DF-W 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 DF-U 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 (+), DF-W (diamond with solid line),
DF-U (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 (+), DF-W (diamond with solid line),
DF-U (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 (+), DF-W (diamond with solid line),
DF-U (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 DF-W 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 DF-U 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 (+), DF-W (diamond with solid
line), DF-U (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 (+), DF-W (diamond with solid line),
DF-U (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 (+), DF-W (diamond with solid line),
DF-U (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 DF-W 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 DF-U 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 (+), DF-W (diamond with
solid line), DF-U (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 (+), DF-W (diamond with solid line),
DF-U (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 (+), DF-W (diamond with
solid line), DF-U (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 DF-W 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 DF-U 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

anti-jam capabilities of GPS position location, but until this work no extension of anti-jam

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, jam-resistant 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

Cram6r-Rao 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 anti-jam 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 anti-jam 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 anti-jam 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 anti-jam 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 anti-jamming 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 multi-sensors for anti-jam 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 anti-jam 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 Cram6r-Rao 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 anti-jam 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 maximum-likelihood 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, hand-held 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 Sequence-Spread Spectrum (DS-SS) 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 local-level 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 "North-East-Down" (NED) or "East-North-Up" (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 "Front-Right-Below" (FRB) convention. The attitude is then

defined as the relation between the local-level 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 "3-1-2" 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 = Q-1 (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 non-collinear 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 multi-user 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, zero-mean

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 so-called 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 ANTI-JAM 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 Sequence-Spread Spectrum (DS-SS) 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 single-antenna 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 anti-jam 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 Anti-Jam 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 sensor-to-sensor base-
line vectors, provide the data necessary for attitude determination in a low interference

environment.
Historically, anti-jam antenna arrays have operated as pre-processors to a GPS

receiver [3,30-32]. The function of the anti-jam 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 anti-jam antenna is often designed

as an applique, the post-antenna 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 pre-processor type anti-jam 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 array-processing levels of anti-jam 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 In-phase (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 anti-aliasing 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 multi-sensor 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 anti-jam

array pre-processor) 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 top-level architecture for a multi-sensor GPS

receiver that combines the anti-jam benefit of adaptive beamforming with development
of data required for anti-jam 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 pre-processor

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 re-derives 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 multi-source

direction finding problem is reviewed, to introduce the multi-source 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[7-1,7-2,...,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 so-called space-time 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
space-time 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)]HR-t [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 space-time 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 space-time 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 de-spread and integrated, increasing the signal to interference ratio.

This is accomplished when the estimated temporal steering vector is used to de-spread the

signal. De-modulation 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)Rt-1u1 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 local-level frame is included in

the estimation process.

An Alternate Array Response Vector
We begin by re-examining 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 non-colocated
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 re-parameterization 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 space-satellite 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 multi-user 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
W-2
W 2= (4.42)


WL
such that

U = Z + W (4.43)

Finally, let Rs. denote the ML x ML space-satellite 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 space-satellite 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)jHC-1[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 ufC-luj j(*a(v, q)HC lu -
q=


jui-CT 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)jHC--lu 12
q= argmax (vq)C-la(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 HC-1/2p C-I1/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[PIC-1] (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 space-satellite 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 non-jammed 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 antenna-frame 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 u-iC 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) = uC-1/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)
NV--oo /

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 H-1/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 Cauchy-Schwarz 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 uHC-1/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 { uHC-I/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 { >oC-P Idi(l)-7y (5.24)

where wl, defined in (4.37), may be written as

w, = u, 7-a(() (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'HC-1/2 [9-1/5.2
q = q-2 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 psuedo-inverse, a&, is found by direct evaluation to be

((1= (wH-1 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:

-- [P-i( 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 ujCT-l2[ r t (4H'H -1~~) p-t
i =2 Re H 12 it) Hd ()Pd,()at + P ] C-l2u/ } (5.42)
8qOqiO-'k
Now recall that asymptotically,
ul - 7^a, (5.43)
and therefore
ufCt-1/2 I (5.44)
resulting in T2 being pre-multiplied by yaHP0 which is of course zero. Therefore, (5.42)
asymptotically approaches
82G(q) 2Re{- 'yiI2'2' [1 (it)H-dk(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 h|21(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~q-f (5.61)

The Cram6r-Rao Bound, PCR is a lower limit on the covariance matrix of the estimation
error for unbiased estimators, and an estimator that achieves the Cram6r-Rao 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 Cramer-Rao 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)(4-qf [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 Cramer-Rao 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 non-linear in the attitude q, and therefore the attitude cannot be found

directly through matrix inversion or psuedo-inversion. 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 unit-norm 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 no-noise 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 = q-4 (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 psuedo-inverse 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)C-a(0) (6.16)

where 7-y 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 = C-a (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 multi-user 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 local-level 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 attitude-space. 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 attitude-space, 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 attitude-space, 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 6-10

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 = 7LIA-L2 (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 Cram6r-Rao 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
DS-SS 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 space-satellite 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 space-satellite 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, Ras-a, 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
Rss-a = (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

dual-frequency 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" dual-frequency 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)
Pa-j()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 dual-frequency MLAE.
Several observations may be made about the dual-frequency 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 dual-frequency 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 dual-frequency 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 dual-frequency 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

dual-frequency 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


iE-20

-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
1-4 3681 5052 683
1-5 3389 4690 153
1-6 5028 3313 55
2-3 6169 1409 28
2-4 733 4670 523
2-5 3875 1743 89
2-6 4358 1826 25
3-4 12095 3667 159
3-5 1772 1654 9
3-6 1 463 1788 01








7.4 Dual-Frequency MLAE Performance Increase
In this section we analytically examine the performance of the dual-frequency MLAE.
Specifically, we will show that the dual-frequency 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 dual-frequency MLAE has the same desirable properties as the single
frequency MLAE.
We begin be examining the dual-frequency 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(v-L,A2, q) L+1 C< < 2L

1 <= ) (7.23)
U(-L,2) L+1

c= (7.24)
(C-L,2) L+1 (, | C(i) l<_<_L
j ( 1)< ( -(7.25)
[C(C-L,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" dual-frequency
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 l-1/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 Cramer-Rao 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 dual-frequency MLAE performed as well
or better than the single frequency MLAE using data from either frequency, consider that
the Fisher information matrix for the dual-frequency, FIMDF is the sum of FIM for each
frequency. That is,
f2L1
FIMDF = 2Re { L I2I5H(C)P (q)6(0)

= 2Re IY|c21H(()P )(q) l(


+ 2 Re { |_7(c2|H(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 dual-frequency 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
Cram6r-Rao Bound for the dual-frequency 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 dual-frequency 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 "widelane-like"

quality, and allows for better quality attitude estimates without introducing ambiguity.