ABSTRACT
Modern air and missile defense (AMD) systems are required to encounter targets at ranges beyond the horizon using multiple sensors and assets for prosecuting the engagement� Search, detection, transition to track, fire control, midcourse guidance, and terminal homing, all parts of the air defense interceptor engagement, may be performed from different assets that will most likely be separated by hundreds of kilometers� These engagements are made possible by accomplishing an accurate location of each asset that is part of the engagement and passing data� The precise relative instantaneous location of assets can be accomplished using satellite-aided (e�g�, GPS) navigation and targeting and accurate models of the earth geoid� The data passed onto the data links will include target positions, velocities, launch positions, target and interceptor relative positions, and velocity� Accurate flight and state representation parameters are paramount to developing firing doctrines and fire-control solutions and computing midcourse and terminal guidance commands� This section provides the fundamental mathematics of locating and translating solutions from one reference frame to another and accurately computing flight performance in a target environment�
9.1.1 A WGS-84 Oblate, Rotating Earth Model
The WGS-84 oblate, rotating earth (ORE) model [1] is likely the most widely used by the air and missile defense community and is provided in some detail here� The WGS-84 model accounts for rotation effects such as tangential and centripetal acceleration and oblateness effects like geodetic versus geocentric positioning and nonuniform gravity� Moreover, the ORE model includes rotationally induced forces and effects encountered by a missile in flight� Although the WGS-84 ORE is not a precise representation of the actual earth, because it does not take into account surface features such as mountains and higher-order shape anomalies, its ellipsoidal representation
is more accurate than a spheroid with constant radius that is typically insufficient for serious AMD design and analysis�
9.1.1.1 Transformation Matrices: Coordinate Frames and Position
From Etkin [1,2], there are many frames of reference encountered in flight dynamics� For example, the earth-centered inertial and launch-centered inertial Cartesian (ECIC and LCIC, respectively) frames, the locally level vehicle (LLV) carried reference frame that is oriented to the earth surface fixed Cartesian (ESFC) frame and sometimes with the LCIC, and other various coordinate systems are used in the design and development of air and missile defense systems� The body frame is centered at the missile center of gravity, with the positive x-axis out the nose of the missile� The wind frame has the same origin as the body frame but with x-axis along the direction of the velocity vector (wind)� The locally level vertical frame has the same missile center of gravity origin, but with the x-axis directed north toward the pole of the earth, and the y-axis directed east and the z-axis directed down, toward the center of the earth�
Transformation between any of these systems requires a multiplication by a transformation matrix that relates the orientation of each frame to one another and an addition of the distance between the origins of each reference frame being transformed� For example, to transform from LCIC to ECIC frames, a multiplication by the transformation matrix is required, along with the addition of distance between the two frame origins in Cartesian coordinates�
Figure 9�1 [1-3] provides the rotation sequence for the derivation of the transformation matrix between ECIC and LCIC given in the following equations:
¢ ¢ ¢
é
ë
ê ê ê
ù
û
ú ú ú
= é
ë
ê ê ê
ù
û
ú ú ú
= ( ) ( )
- ( ) X
Y
Z
B
X
Y
Z
B
; cos sin sin co
0 s ( )
é
ë
ê ê ê
ù
û
ú ú ú
0 0 0 1
(9�1)
N
E
D
B
X
Y
Z
B é
ë
ê ê ê
ù
û
ú ú ú
= ¢ ¢
é
ë
ê ê ê
ù
û
ú ú ú
= - ( ) ( )
m m
m
’ ;
sin cos
cos
0 0 1 0
( ) ( )
é
ë
ê ê ê
ù
û
ú ú ú0 sin m
(9�2)
X
Y
Z
B
N
E
D
B az
é
ë
ê ê ê
ù
û
ú ú ú
= é
ë
ê ê ê
ù
û
ú ú ú
= ( ) ( )
- (l l l l l;
cos sin sin
) ( ) é
ë
ê ê ê
ù
û
ú ú ú
cos l 0 0 0 1
(9�3)
The transformation shown in Figure 9�1 only involves rotations to orient the axes in angle along the LCIC frame launch azimuth relative to the ECIC
frame and does not involve any position translation� Thus, with this matrix, a point in space defined in the ECIC coordinate system can be defined in the LCIC frame� The first rotation () shown in Figure 9�1 is a rotation along the earth longitudinal angle about the ECIC vertical z-axis to locate the launch frame in the easterly plane� The second rotation (μ) is a rotation about the ECIC X-Y plane in earth latitude to locate the launch frame in the northerly plane� After carrying out the two previous rotations, the resulting axis is a north-east-down (NED) system, with each axis pointing in one of those three directions� The launch frame, however, may not be aligned with the NED frame, and thus, the third rotation (λ) aligns the launch frame along the firing azimuth� The LCIC z-axis always aligns toward the center of earth� The alignment demonstrated in Figure 9�1 assumes a spherical earth model noting that there is no difference between geodetic and geocentric latitude�
Combining all three rotation matrices by multiplication, the overall ECIC to launch the transformation matrix is achieved, shown in Figure 9�2 [3], and the resulting transformation matrix is given in the following equations:
X
Y
Z
B B B
X
Y
Z
B B B B B
b b baz
é
ë
ê ê ê
ù
û
ú ú ú
= é
ë
ê ê ê
ù
û
ú ú ú
b b b
b b b
é
ë
ê ê ê
ù
û
ú ú ú
(9�4)
b
b
= - ( ) ( ) ( ) - ( ) ( )
= - ( ) ( ) (
cos sin cos sin sin
cos sin sin
l m l
l m
) + ( ) ( )
= ( ) ( )
= ( ) ( ) ( ) -
sin cos
cos cos
sin sin cos cos
l
l m
l m l
b
b
21 ( ) ( )
= ( ) ( ) ( ) + ( ) ( )
= - ( )
sin
sin sin sin cos cos
sin cos
b
b
l m l
l m
m
m
m
( )
= - ( ) ( )
= - ( ) ( )
= - ( )
b
b
b
cos cos
cos sin
sin
(9�5)
With a nonrotating earth, the position of the launch frame with respect to ECIC coordinates is fixed, assuming that the launch frame does not translate during the flight of the missile� There are no rotation effects in the
transformation from one frame to the other� Modeling a rotating earth, however, requires the position of a stationary launch frame on the surface of the earth to change with respect to the ECIC frame� The earth surface rotates with a fixed angular velocity about the axis of rotation (normally assumed to be straight through the pole)� This change in position is dependent on time relative to an initial time-t0 seconds�
At t0, the launch frame is fixed relative to the Greenwich meridian assumed as the plane through which the ECIC x-axis passes� At time t0 + τ, the launch frame’s angular position has changed as a function of the earth rotation rate vector, and the position vector change is also related to the location of the frame relative to the equator (μ)�
The [B] transformation matrix will need to be modified by substituting + ω(τ − t0) for to account for the change in position of the launch frame with time due to the rotation of the earth’s surface about its center, where τ is the time from initial conditions� The effect is illustrated in Figure 9�3 [3]�
The new transformation matrix will be labeled BRE to denote the rotating earth transformation from ECIC� Notable is the fact that this timedependent term only needs to be added to the ECIC-to-Launch Frame transformation matrix or similar earth surface fixed frames and any other transformations dependent on earth rotation effects� The new general model is depicted in Figure 9�4 [3]:
X
Y
Z
B
X
Y
Z
B
b b b
b az
ë
ê ê ê
ù
û
ú ú ú
= é
ë
ê ê ê
ù
û
ú ú ú
b b
b b b
é
ë
ê ê ê
ù
û
ú ú ú
(9�6)
b t t
b
= - ( ) ( ) + ×( ) - ( ) + ×( )
= -
cos sin cos sin sin
c
l m w l w D D
os sin sin sin cos
cos cos
l m w l w
l
( ) ( ) + ×( ) + ( ) + ×( )
= ( )
t t
b
D D
13 m
l m w l w
( )
= ( ) ( ) + ×( ) - ( ) + ×( )b t t
b
sin sin cos cos sin D D
= ( ) ( ) + ×( ) + ( ) + ×( )
= - ( )
sin sin sin cos cos
sin
l m w l w
l
t t
b
D D
23 cos
cos cos
cos sin
m
m w
m w
( )
= - ( ) + ×( )
= - ( ) + ×( )
b t
b t
b
D
D
33 RE = - ( )sin m
(9�7)
9.1.1.2 Transformation Matrices: Velocity and Acceleration
Modeling a rotating earth requires the introduction of a time-dependent coordinate frame component when developing transformations having to do with either velocity or acceleration [1,2]� Applying the chain rule to the transformation of a position vector pa in one frame of reference to another pb given in Equation 9�8 to obtain the velocity transformation given in Equation 9�9 reveals that it is necessary to determine if the transformation matrix is time dependent:
p L pb b a
a= × (9�8)
p L p L pb ba a ba a= + (9�9)
Notice that when the derivative of the transformation matrix exists, it is necessary to transform either velocity or acceleration vectors�
Another effect of a rotating earth that will require transformation is the additional centripetal and tangential accelerations [1]� Tangential acceleration exists since there is a change in speed of the interceptor due to the rotation of the earth, and centripetal acceleration is created from the change in direction of the vehicle velocity due to rotation� Acceleration terms are given in the following equations:
Centripetal acceleration = ´ ´( )w we e ECICp (9�10)
Tangential acceleration = ´we bv (9�11)
Rotation also introduces an initial velocity relative to the earth’s center for the vehicle before launch� Since the earth’s surface rotates with respect to ECIC coordinates, a missile at rest before launch still has velocity with respect to the ECIC frame� The equation used to calculate this initial missile velocity is given in the following:
v pinit e ECIC= ´w (9�12)
9.1.1.3 Oblateness Effects, Nonuniform Gravity
The latitude angles extending from the earth’s center are different for an oblate earth model compared to a spherical model [4], and thus, various equations involving those angles must be defined� As shown in Figure 9�5 [3], there are two latitudes for an oblate model, the geocentric latitude from the earth’s center, identified as μc, and the geodetic latitude, μd, on the semimajor axis directly below the earth’s surface at the specified point� For a spherical model, these two were equivalent since a direction perpendicular to the earth’s surface at any given point always intersected the earth’s center� The term latitude on an oblate body, such as the earth, refers to the geodetic value�
The coordinate pair (vs, zs) refers to the earth surface geodetic subvehicle point� The coordinate pair (vp, zp) refers to the vehicle point at a geodetic altitude, H, above the surface� The oblate earth model also introduces a flattening factor, f, which describes the deviation of the geoid from a perfect sphere or the ellipticity of the earth and can significantly influence calculations and equations having to do with the flight equations of motion� One effect of an
oblate earth is a nonuniform gravitational field [4]� In a spherical model, the gravity field vector is a constant toward the center of the earth regardless of latitude and longitude� With an oblate model, gravity is highly dependent on latitude yet still independent of longitude, assuming the model is flattened at the poles� Kaplan [4] gives the earth gravitational potential function V
as 9�8:
V GM
r a r
P C m S m n
= + æ
è ç
ö ø ÷ ( ) +( )
ë ê ê
ù
û ú ú
(9�13)
The term GM is the universal gravitational constant, r represents the distance from the earth center, and the subscripts (m, n) represent the order and degree, respectively, of the associated Legendre function (P) and the gravitational coefficients (C, S)� The overbar represents normalization according to the following:
n m
n m n k
+( ) -( ) × +( )
é
ë ê ê
ù
û ú ú
× ( )! ! 2 1
·
(9�14)
The associated Legendre function is defined in the following equation:
P
d
sinm m m
m( ) = ( ) × ( )
× ( )éë ùû (9�15)
The Legendre polynomial is defined in the following:
P
n d
n sin
! sin sinm
m m( ) = ×
( ) -( )1
2 12 (9�16)
A perfectly oblate earth model that is symmetric about the z-axis has no zonal harmonics beyond a degree, n = 2� The larger the value, n, the less impact the magnitude of each coefficient has on the overall equation for gravity potential� For sensible atmospheric flight and ballistic flights less than 2000 km in range, a value of n = 2 should be reasonable� The so-called, J2, Jeffry constant typically implies that an oblate spheroid earth model is being assumed� With the J2 assumption, Equation 9�13 can be simplified and is rewritten in the following:
V
GM r
a r
P Cc
= + æ è ç
ö ø ÷ ( )( )
é
ë ê ê
ù
û ú ú
1 2
The -( )C2 0, term in this equation is in normalized form and can be related to the J2 term and is defined [5] as follows for an order value m = 0, when the constant Cn m, is assumed as shown in the following equation:
J C nn n= - × +,0 2 1 (9�18)
Therefore, we can calculate the J term with n = 2 and m = 0 in the following equation:
J2 0 0010826269= . (9�19)
Using the definition for P c2 0, sinm( ), Equation 9�17 reduces to the following equation, with m mc df= -( ) ( )éë
ù û
1 [1]:
V
GM r
J a r
= - × æ
è ç
ö ø ÷ × ( ) - }{
é
ë ê ê
ù
û ú ú
1 0 5 3 12 2
2. sin m (9�20)
The gradients of the potential function are evaluated in ECIC coordinates to give the following x, y, and z components of gravity where XE, YE, and ZE are the vehicles’ respective ECIC coordinates:
G GM r
J a r
XE r
G GM r
J a
= + æ è ç
ö ø ÷ -( )
é
ë ê ê
ù
û ú ú
= +
1 3 2
1 5
1 3 2
sin m
r ZE
r
G GM r
J a r
æ è ç
ö ø ÷ -( )
é
ë ê ê
ù
û ú ú
= + æ è ç
ö ø ÷ -
1 5
1 3 2
1 5
sin
si
m
n2 mc ZEr( ) é
ë ê ê
ù
û ú ú
(9�21)
9.1.1.4 Geodetic and Geocentric Latitude Relationship
The mathematical relationship of the geocentric and geodetic latitude is obtained� Refer to Figure 9�5� The two points, one on the earth’s surface (vs, zs) and one directly above that point (vp, zp), can be described in two separate but related methods� The relation between the two geodetic and geocentric coordinates is the flattening factor that gives the relation between the two angles� In the following equations, we will derive the relationship between geodetic and geocentric latitude and produce another equation to calculate vehicle position in ECIC coordinates given its longitude and geodetic latitude and altitude from Sooy [3]� From Figure 9�5, the following equation is defined:
tan( )md
z v
= (9�22)
The equation of an ellipse from Figure 9�5 is provided in the following equation:
v a
z b
2 1+ = (9�23)
The earth equatorial radius is a, and b is the polar radius� If we set b = (1 − f)a, we can rearrange Equation 9�23 and obtain the following equation:
z f a v2
2 2 21= -( ) -( ) (9�24)
The slope of the line tangential to the ellipse for any given v and z is given in the following:
m
dz dv
v z
f= = - -( )1 2 (9�25)
The slope of the line perpendicular to Equation 9�25 is the negative reciprocal� The perpendicular slope passing through (vs, zs) is then given by the following equation:
m z
v f s
^ = -( )1 2
(9�26)
The equation of the slope of Equation 9�26 is also the equation for tan(μc)� Thus, with substitution, we have a relation between tan(μc) and tan(μd) given as follows:
tan( ) tan( )m md c s
m f
z
v f = =
-( ) =
-( ) ^
1 1 2 2 (9�27)
9.1.1.5 Geodetic to Geocentric Latitude: Vehicle Position
Most often, vehicle location is provided in longitude, geodetic latitude, and altitude� Flight performance computation is most often accomplished in geocentric or ECIC coordinates� To calculate vehicle position in ECIC coordinates given geodetic coordinates, we will solve Equation 9�27 for zs and substitute this into the equation for the ellipse to give us the following equations:
v a
v f
1 1+
-( ) ( ) -( )
= tan m
(9�28)
Solving for vs (magnitude) gives
v a
f s
= + -( ) ( )1 1 2 2tan m
(9�29)
Rearranging vs (sign and magnitude)
v a
f s
= ( )
( ) + -( ) ( ) cos
cos sin
m
m m2 2 21
(9�30)
Equation 9�30 can be simplified to the form given in Equation 9�31� Substituting vs and f into Equation 9�27 and solving zs provide the subvehicle point on the oblate earth:
v a
f f s
= ( )
+ -( ) ( ) cos
sin
m
m1 2 2 (9�31)
Assuming the vehicle is at some altitude (H) above the earth surface, the geodetic subvehicle point (vs, zs) will need to be adjusted to locate the vehicle in ECIC coordinates� From [4], the altitude-adjusted differential values for (vs, zs) are defined in the following equations:
Dv H d= ( )cos m (9�32)
Dz H d= ( )sin m (9�33)
The coordinates of the elevated point (vp, zp) are from [4] and given in the following equations:
v v v a
f f Hp s
d= + = ( )
+ -( ) ( ) + ( )D cos
sin cos
m
m m
1 2 2 (9�34)
z z z v f Hp s s d d= + = -( ) ( ) + ( )D 1 2
tan sinm m (9�35)
The ECIC coordinates of a vehicle can then be found by using the following equations when given only longitude (), geodetic latitude (μd), and geodetic altitude (H):
x vECI p= ( )cos (9�36)
y vECI p= ( )sin (9�37)
z z v f HECI p s= = -( ) ( ) + ( )1 2
tan sinm m (9�38)
To transform the launch point (launch frame coordinates) and multiple vehicle positions to ECIC coordinates, an additional transformation will be required due to the oblate earth representation� This new transformation results from the difference in the angle between geodetic and geocentric latitude� For the spherical earth case, only three rotations were involved� Notice in Figure 9�5 that the oblate earth model now requires a fourth rotational transformation to align the vehicle’s downward z-axis vector with the earth’s center� The angle between the line perpendicular to the surface at a given point and the line from that point to the center of the earth is equal to the difference between the geodetic and geocentric latitude angles� This additional rotation is accomplished with another transformation matrix and accomplished in the same manner as previous transformation matrices�
9.1.1.6 Latitude, Longitude, and Altitude Calculation
To compute latitude, longitude, and altitude during flight is often required when designing midcourse guidance or when simply examining time-space correlated flight performance� Assume the engagement computational tool begins with latitude, longitude, and altitude of the launch frame with respect to the ECIC system� The initial target and interceptor position are acted on computationally in geocentric launch coordinates, but it is desired to reference those time-space correlated interceptor positions to geodetic assets� Thus, it is necessary to convert those changes in launch coordinates to ECIC positions using transformations and then calculate the new geodetic latitude, longitude, and altitude of each asset� These new values can then be used by the guidance algorithm(s) to direct the interceptor to the target and enable the user to correlate the trajectory of each vehicle in geodetic space� From [4], the following equations are given, where XE, YE, and ZE are the missile position in ECIC coordinates, e2 is the eccentricity squared, and RΘ is the equatorial radius of the earth:
md ZE
e XE YE =
-( ) + é
ë
ê ê
ù
û
ú ú
(9�39)
= é
ëê ù ûú
arctan YE XE
(9�40)
H XE YE ZE e e
R
e = + +
+ ( ) -( ) - - ( ) 2 2 2
2 2 2 2 21 2 1sin sinm m Q (9�41)
The assumed values for the WGS-84 ORE model [5] are defined as the following parameters and values� The first four parameters are the fundamental ones that can be used to recalculate all other coefficients as updates are required:
where RQ is the equatorial radius of the earth = 6�378137 * 106 m wQ is the earth’s angular rotational rate = 7�29211585530 * 10-5 rad/s C2,0 = −484�1654663 * 10-6 μ = 3�896004418 * 105 km3/s2 e2 is the eccentricity squared = 0�006694385000 f is the flattening factor = 0�003352813000
The ORE provides a realistic ellipsoidal model to more accurately evaluate the flight performance of the AMD system and interceptor and specifically
may be required when developing midcourse guidance designs� The earth’s rotation introduces several new calculations including centripetal and tangential acceleration and time-dependent transformation matrices� Oblateness introduces a nonuniform gravity field and various new equations to calculate latitude, longitude, and altitude� A comprehensive treatment of earth models, coordinate systems, and transformations can be found in Valado [5]�
9.1.1.7 Forces and Moments and Equations of Motion
We define aerodynamic forces and moments consistent with AIAA, R-0041992 [6]� The aerodynamic reference system is defined in Figure 9�6� Normal, side, and axial forces (NF, YF, AF) and rolling, pitching, and yawing moments (Lm, Mm, Nm) are functions of independent variables angle of attack (α), angle of sideslip (β), aerodynamic roll angle (ϕ), Mach number and control surface deflection angles (δ), and the control surface aerodynamic incidence angle in the A-plane (i′) and B-plane (i)� Subsequently, it is necessary to define the relationships between the angles and the body axis velocity components to ensure proper interpretation of aerodynamic requirements�
The mathematic relationships in the following equation hold for Figure 9�6:
a b f
a
= = =
= + =
- - -tan ( ); tan ( ); tan ( )
tan ( )
tan
w u v u w u
v w u
/ / /
a b+
= + +
tan2
2 2 2V u v w
(9�42)
The total velocity vector is defined as V, and the total angle of attack is defined as αT�
Two reference frames are of interest before developing the complete sixdegree-of-freedom (6DOF) equations of motion (EOM)� We should assume that there is a locally level vehicle (LLV) carried reference frame that is oriented to the earth surface fixed Cartesian (ESFC) reference system, and we should assume that there is a missile body fixed reference system (B)� Euler angles (ψM, ΘM, ΦM) represent the yaw, pitch, and roll orientation angles of B relative to LLV�
The reader should refer to Etkin [2, pp� 114-117] for details concerning this section� Vehicle body relative Euler angles (yaw, ψb; pitch, Θb; roll, ϕb) are normally employed to measure the angular displacement of the vehicle body attitude relative to the ESFC or LLV frame and will be followed here� Once the body Euler angles are found, they are used in a direction cosine matrix to compute the remaining required flight relationships�
The rate of change of the Euler angles between the ESFC/LLV frame axes and the vehicle body axes is defined by the Euler rate equations in terms of the body rotational rates (pb, roll rate; qb, pitch rate; rb, yaw rate)� The body rates are typically measured by the inertial reference unit (IRU) with rate gyroscopes� Homing interceptors are typically roll attitude stabilized with the exception of some shorter-range point defense systems� Moreover, as long as a relatively tight roll attitude control is maintained, the roll rate can be assumed zero for roll attitude-stabilized systems and the Euler rate equations can be simplified by setting Fb = 0�
The LLV to vehicle body axis direction cosine matrix [A] is defined in terms of the Euler angles and given in Equation 9�43:
[ ]A a a a
a a a
a a a LLV B =
é
ë
ê ê ê
ù
û
ú ú ú
(9�43)
a
a
a
a
=
=
= -
= -
cos cos
cos sin
sin
sin sin cos
Q Y
Q Y
Q
F Q Y cos sin
sin sin sin cos cos
sin cos
F Y
F Q Y F Y
F Q
a
a
a
= +
=
= +
= -
cos sin cos sin sin
cos sin sin sin cos
F Q Y F Y
F Q Y F Y
a b b33 = cos cosF Q
(9�44)
The [ ]A LLVB matrix is an orthogonal transformation, and therefore, the direction cosine matrix to reverse the direction and go from B to LLV frame is its transpose and given in the following equation:
[ ]A a a a
a a a
a a a B LLV =
é
ë
ê ê ê
ù
û
ú ú ú
(9�45)
The LLV velocity is found by multiplying the [ ]A BLLV matrix and the vehicle body frame velocity (u, v, w) components shown in the following equations:
R
R
R
A
u
v
w
é
ë
ê ê ê ê
ù
û
ú ú ú ú
= é
ë
ê ê ê
ù
û
ú ú ú
[ ] (9�46)
R ua va wa
R ua va wa
R ua va wa
= + +
= + +
= + +
(9�47)
It is the purpose of flight analysis to solve the earlier set of equations to understand the flight behavior of the interceptor under boost, midcourse guidance, or terminal homing to assess design alternatives� This is typically formulated in an engagement simulation where a set of parallel target equations exist, and they are solved relative to one another� It is instructive here to first establish the mechanization to solve the primary vehicle set of equations, and then generality applies to characterize a target in a similar fashion if this detail is warranted�
We first establish the effect that the rotating earth has on the vehicle body� Assuming an ORE model, we establish the following equation as the set of rotational rates imparted on the vehicle body axes from the earth:
p a a
q a a
r
= × × - ×( )
= × × - ×( )
w m m
w m m
cos sin
cos sin
E e d da a= × × - ×( )w m m31 33cos sin
(9�48)
The resulting translational force equations can now be developed along the body axes� They are given in the following equation� The applied forces along the body axes include aerodynamic, propulsion, and environmental (wind)
effects� The ORE gravity vector components gX, Y, Z are along the LLV frame and rotated using the body Euler angles to align with the vehicle body:
u X m
g r r v q q w
v Y m
g
= - × + +( ) × - +( ) ×
= - × ×
sin
cos sin
q
q f + +( ) × - +( ) ×
= - × × + +( ) × -
p p w r r u
w Z m
g q q u p
E b b cos cosq f E bp v+( ) ×
(9�49)
Rotational equations of motion are developed assuming rigid-body dynamics and that a set of coordinate axes where the principal axes exist and the inertia matrix is a diagonal� Therefore, there are no products of inertia and the inverse of the inertia matrix is also a diagonal and Euler’s equations of motion hold� Moments or torques are the result of force vectors applied about (at a distance from) the center of gravity and the body responds to in the form of angular accelerations� These vehicle body rotational equations of motion hold regardless of the earth model implemented and are given in the following equation:
p L I
I I q r I
q M I
I I r p I
r
= + -( ) ×
= + -( ) ×
= N I
I I p q I
ZZ +
-( ) ×
(9�50)
The Euler angle rotation sequence is not commutative, and it is important to retain whatever sequence throughout the flight dynamics derivation and analysis� The rotation sequence is given here as (yaw, pitch, roll):
P
Q
R
p
q
r
A b
é
ë
ê ê ê
ù
û
ú ú ú
= é
ë
ê ê ê
ù
û
ú ú ú
-[ ] +( )
-
- +
w m
m
w
cos
( )
é
ë
ê ê ê ê
ù
û
ú ú ú úsinmd
(9�51)
Y F F Q
Q F F
F F
Q R
Q R
P Q R
= +
= -
= + +
( sin cos ) cos
cos sin
( sin cos )tanF Qb b
(9�52)
The Euler angle approach used to define vehicle orientation in space contains a singularity by inspection of the previous equation� As the second, or
pitch rotation (Θ), approaches ±90°, the yaw (Ψ) and roll (Φ) Euler angle rates become undefined� To eliminate this problem, flight dynamics has evolved to using quaternions [1] to compute body attitude�
Quaternions are mathematically defined but have no physical meaning and in themselves have no usefulness in the study of flight dynamics� Quaternions are therefore only useful from a computational perspective� They must be converted back to Euler angles before any subsequent analysis can be accomplished� Specifically, quaternions are mathematical functions of the direction cosines of the vehicle orientation�
The initial vehicle orientation needs to be specified in Euler angles and then are converted to four quaternion values (e0, e1, e2, e3) that are used to initialize the quaternion rates when combined with the body rates� These quaternion rates are then integrated and converted back to Euler angles for interpretation and to begin the next computational iteration� The sequence for implementing quaternions is given in the following equations:
e
e
1 2
= × × + × ×[ ]
=
cos( ) cos( ) cos( ) sin( ) sin( ) sin( )y q f y q f
1 2
cos( ) cos( ) cos( ) sin( ) sin( ) sin( )
co
y q f y q fb b b b b b
e
× × - × ×[ ]
= s( ) sin( ) cos( ) sin( ) cos( ) sin( )
cos(
y q f y q fb b b b b b
e
× × + × ×[ ]
= -3 1 2
y q f y q fb b b b b b) sin( ) sin( ) sin( ) cos( ) cos( )× × + × ×[ ]
(9�53)
e e P e Q e R
e e P e R e Q
e e Q e P e R
1 2
1 2
1 2
= + +[ ]
= + +[ ]
= + -[ ]
e e R e Q e P3 0 1 2 1 2
= + +[ ]
(9�54)
y
q
e e e e e e e e
e e
= × -( ) + - -
é
ë ê
ù
û ú
= - × -
tan
sin
2 e e
e e e e e e e e
( )éë ùû
= × -( ) + - -
é
ë ê
ù
û ú
-j tan
(9�55)
Clutter is a source of noise that interferes with target detection� By definition, clutter is any kind of unwanted echoes that originates from raindrops, birds, water, and ground surfaces� The occurrence of clutter may be a serious problem since it could result in echoes stronger than the intended target� It is often possible to discard much of the clutter by simply ignoring echoes from slow-moving objects using a Doppler gate, that is, by removing echoes not having the expected frequency� Using an appropriate frequency for the radar will also help decrease clutter since the reflectivity is frequency dependent� It is not always realistic to remove all clutter, and possible clutter must be taken into account in some real-life situations� Fortunately though, clutter can be removed with good results when using Doppler techniques and slow or stationary surface returns can be neglected�
9.2.1 Radar Returns
The radar range equation is used to determine both the target and clutter powers returned from a given range� The monostatic power returned from a target or clutter is a function of the radar parameters, range, and radar cross section as follows:
P r
PG A r L
[ ] = s
p4 2 4 (9�56)
where Pt is the peak transmitter power Gt is the transmit gain of the radar antenna Ae is the effective receiving area of the radar antenna σ is the radar cross section of the target or clutter r is the slant range of the target or clutter L is total transmit, receive, and processing losses
A complete and rigorous derivation of the monostatic radar range equation is presented by Blake [7]� For distributed clutter, the clutter radar cross section is also a function of range� The following sections describe how the clutter radar cross section can be estimated for both surface and volume clutter�
9.2.2 Surface Clutter Returns
The geometry for a simple surface clutter model is shown in Figure 9�7 for a sea/land environment where mountains rise from the seashore�
The model is general in the sense that it accounts for rising terrain overland when determining the surface area intercepted by the main beam of the radar antenna at height hr. The range of the surface clutter patch from the radar is rc� The radar antenna and clutter heights with respect to a mean sea-level reference are hr and hc, respectively� The perpendicular distance from the radar to a line tangent to the clutter surface at the point where the antenna main beam axis intercepts the illuminated patch is the adjusted height, ha, of the radar for an equivalent flat land scenario� In the model, the main beam axis is always assumed to intersect the illuminated clutter patch at each range�
Using the model geometry, the surface length of the clutter patch in both the azimuth and elevation extents of the antenna main beam can be determined� The clutter patch surface area is computed assuming the intersection of the main beam and clutter surface is an elliptical surface� The clutter surface length intercepted in the azimuth direction is a function of the antenna azimuth beamwidth, while the clutter surface length intercepted in the elevation direction is a function of either the pulse width or the antenna elevation beamwidth depending on the specific geometry� When the elevation surface length is a function of the pulse width, this is referred to as a pulse-limited clutter patch� Beam-limited clutter patch refers to the case of functional dependency on the antenna beamwidth� Beam limiting usually occurs for closerange clutter, with a transition to pulse limiting with increasing clutter range�
To calculate the elevation surface length intercepted by the antenna beam for the pulse-limited case, the angle φ must be determined� The elevation surface length intercepted can be determined from
z t jep c= 1 2/ sec (9�57)
where ζep is the elevation direction surface length intercepted c is the speed of light τ is the pulse width
The angle φ is also necessary to compute the surface length intercepted for the beam-limited case, that is,
z q jeb c Er= csc (9�58)
where θE is the antenna 3 dB elevation beamwidth� When the pulse-limited length, ζep, is greater than the beam-limited length, the surface clutter is beam limited� The clutter patch surface area is given for both the following cases:
S rp ep c E= z q pulse limited (9�59)
S rb eb c A= p z q/ 4 beam limited (9�60)
The surface clutter patch area is a truncated elliptical region for the pulselimited case and elliptical for the beam-limited case�
The radar cross section of the distributed clutter illuminated can be estimated:
s sc oS= (9�61)
where σo is the mean value of the areal backscatter coefficient for the clutter� Typical values for the areal backscatter coefficient for different clutter types can be found in the literature [9]�
Once the clutter patch radar cross section is determined, the received clutter power can be determined by substitution of the clutter patch radar cross section, σc, into the radar range equation� The clutter power returned for the pulse-limited case decays as the inverse cube of range� At the transition range, the beam-limited and pulse-limited clutter patch areas are equal, providing a smooth transition�
9.2.3 Volume Clutter Returns
For volume clutter such as rain, it is necessary to determine the illuminated clutter volume at each range� The illuminated clutter volume can be treated as two separate cases: (1) the beam (radar range-angle cell) is completely filled with clutter and (2) the beam is partially filled with clutter� When the beam is filled with clutter, the beam shape in angle and the pulse width in slant range bound the clutter volume� Assuming an elliptically shaped antenna beam, the clutter volume illuminated for the beam-filled case can be expressed as
V r r cc E c A= p q q t/ /4 2( )( )( ) (9�62)
Since the gain of the transmit and receive beams is not uniform over the areal beamwidth of the antenna, the aforementioned expression is usually reduced by the factor 2(ln 2) assuming a Gaussian two-way antenna pattern [9], that is,
V r r cc E c A= p q q t/( (ln ))( )( )( )16 2 (9�63)
The radar cross section of the volume clutter illuminated is estimated as
s sc vV= (9�64)
where σv is the mean value of the volumetric backscatter coefficient� Typical values for the volumetric backscatter coefficient for various clutter types can be found in the literature [8]� Once the illuminated clutter volume radar cross section is determined, the received clutter power can be determined from the radar range equation� Beam-filled volume clutter returns decay as the inverse square of range�
When the antenna beam (range-angle cell) is not filled with clutter, the clutter volume is no longer bounded in elevation by the beam shape� For nonbeam-filled clutter, two phenomena can occur simultaneously or individually: (1) the beam shape exceeds the clutter maximum height in elevation and/or (2) part of the beam intercepts the surface of the earth� Figures 9�8 and 9�9 illustrate the two cases in a cross-sectional view perpendicular to the main beam axis assuming an elliptically shaped two-way antenna areal beamwidth�
To determine the cross-sectional area illuminated, it is necessary to know the areas of the triangular region and elliptic sector shown in Figure 9�9� The illuminated area of the truncated semiellipse of Figure 9�9 can be determined straightforwardly when the height, h, from the main beam axis to the truncation point is known� The base length of the triangular sector, g, can be determined from geometry as
g a h a b= -[ ( ) ] /2 2 2 1 2/ (9�65)
a rc A= q /2 (9�66)
b rc E= q /2 (9�67)
where a and b are the ellipse semiaxis lengths� At this point, the angle φ can be determined as
j p= - -/ /2 1tan ( )g h (9�68)
Next, the area of the elliptic sector is determined by computing the incremental area of the ellipse as a function of the angle φ� This elliptic sector is determined by integrating in cylindrical coordinates� The result of this integration results in the area of the elliptic sector as
A ab a bE = -( ) tan [( ) tan ]/ /2 1 j (9�69)
The semi-ellipse area illuminated is just twice the sum of the areas of the triangular and elliptic sector regions:
A gh ab a bSE = + -( ) tan [( ) tan ]1 / j (9�70)
The illuminated volume of the partially filled beam can now be approximated as
V A A cSEUPPER SELOWER= +( ) t/( ln )4 2 (9�71)
where the length of the truncated elliptic cylinder is cτ/2� This result assumes that the antenna beam axis is parallel to the earth’s surface at each clutter range�
9.2.3.1 Clutter Processing Considerations
Ambiguous clutter returns occur in an MTI or PD waveform when clutter exists beyond the unambiguous range of the waveform� The unambiguous range is the range associated with the duration of the pulse repetition interval (PRI) or time between successive pulses� For a PRI of 0�1 μs, the unambiguous range is 15 km� Ambiguous range clutter returns compete with targets in the unambiguous range interval for detection� All targets of interest are assumed to lie within the unambiguous range� In addition, the ambiguous or folded clutter return may not be present in each of the pulse returns processed� When this occurs, the clutter filter is no longer matched to reject the clutter return, resulting in decreased clutter rejection� This effect can be quite significant for clutter filters that only process several pulses�
The following sections discuss topics that directly relate to the ability of the clutter processing to reject ambiguous clutter effectively� The topics discussed include the clutter processing interval, maximum clutter range, fill pulses, and clutter rejection degradation�
9.2.4 Coherent Processing Interval
The coherent processing interval (CPI) can be best described by referring to a simple example� In this example, the scenario is limited to three range frames, which contain either a target or a clutter that is processed using three-pulse MTI� A range frame is simply a range interval equal to the unambiguous range, rua, of the waveform� The target and clutter scenario is shown in Figure 9�10 with the target denoted by the caret symbol and clutter denoted by the asterisks� For this example, point clutter is used; however, the results apply to distributed clutter as well� A target is located at a range of 0�5 times, the unambiguous range, and clutter is located at 1�5 and 2�5 times, the unambiguous range�
The superscripts refer to the range frame where the clutter or target physically resides and subscripts refer to the pulse number within the burst from which the radar return is received� For a five-pulse MTI waveform (see Figure 9�10), the time coincident returns in the first three processing intervals are shown in Table 9�1, where the coherent processing interval covers one unambiguous range interval for three contiguous range frames� The number of the CPI corresponds to the pulse number where clutter processing begins, that is, CPI #1 begins at pulse 1, CPI #2 begins at pulse 2, and so on� The data are processed out to the unambiguous range of the last pulse in the burst, and typically, the first available coherent processing interval is used� It is observed that the mixture of target and clutter returns is vastly different depending on which CPI is selected for processing�
For example, in CPI #1, the target return is present in each pulse return; however, the clutter return from range frame #2 is in two returns and the clutter from range frame #3 is only in one return� Close examination of Table 9�1 reveals that the target appears in each range frame processed in CPI #1, that is, a target within the unambiguous range interval will always appear in each frame of CPI #1� Further, the number of missing clutter returns in the first available CPI is equal to the displacement in range frames of the CPI from the actual clutter range� Since missing clutter returns result in degraded clutter rejection, it is desirable to fill the coherent processing interval with clutter� The task of the clutter processing waveform designer is to provide an adequate number of fill pulses for the specified clutter scenario as illustrated
in Figure 9�11� Since real-world clutter environments can be quite complex, a set of waveforms may be required that can be adapted to the clutter environment to optimize clutter rejection and target detection�
9.2.5 Fill Pulses
Fill pulses are used to fill in the missing clutter returns� They are extra pulses in addition to the number required for a particular clutter filter� In this example case, three pulses are required with additional pulses used to fill in clutter returns� Assuming four pulses are transmitted and processing does not begin until after the second pulse is transmitted (CPI #2 processing), the time coincident returns are given in Table 9�2�
Note that CPI #2 (the first available CPI) is now filled except for a single missing clutter return in frame 1 at long range, C3, which will cause some degradation in clutter rejection against the C3 clutter component� If only the second time around clutter is present, then CPI #2 is the preferred processing interval; however, the third time around clutter can be effectively canceled by processing CPI #3 with a two-pulse MTI�
Now, consider a five-pulse burst where processing begins after pulse 3 is transmitted� This corresponds to the case of using two fill pulses� The time coincident radar returns are summarized in Table 9�3�
Now, the clutter returns are completely filled in CPI #3� Clearly, for nth time around clutter, n − 1 fill pulses are needed to completely fill in the clutter in the first available processing interval after the fill pulses are transmitted�
9.2.6 Maximum Clutter Range
The maximum clutter range is defined as the maximum range at which a clutter return can be received and processed� This also corresponds to the minimum range at which a jammer can intercept and repeat a single pulse of the waveform that is received and processed by the radar� The maximum clutter range is proportional to the total number of pulses in the burst and is given as
r c M FP PRI r M FPmc ua= + = +[ ( ) ] ( )/2 (9�72)
where FP is the number of fill pulses M is the number of pulses coherently processed
This is the range beyond which clutter returns and responsive jamming are not received and processed by the radar�
9.2.7 Clutter Rejection Degradation
Range ambiguous clutter returns degrade clutter rejection when the ambiguous clutter returns are not completely filled� Specifically, one or more range frames may not contain the return(s) from a range ambiguous clutter source� When this occurs, one or more of the pulses do not contain the same clutter returns as the other pulses, resulting in reduced clutter rejection (the clutter notch depth is reduced)� For example, in Figure 9�12, the frequency response of a five-pulse MTI filter with Chebyshev weighting is shown for 0 missing returns (1st time clutter), 1 missing return (2nd time clutter), and 2 missing returns (3rd time clutter)�
Multiple time around clutter reduces the ability of the MTI filter to reject the low-velocity radar returns associated with the clutter that is slow moving or stationary� The low-velocity region of the clutter filter is sometimes referred to as the clutter notch� In general, the clutter may be distributed
in velocity and overlap the transition and/or passband of the clutter filter� When this is the case, the clutter cannot be completely rejected unless an MTI filter with a wider clutter notch is used�
9.2.8 Range-Velocity Visibility
A target is said to be visible in a given range-velocity cell if the signal-toclutter pulse noise ratio at the detector meets or exceeds the required level to achieve a specified probability of detection� Range eclipsing and blind velocities result from implementation practices and ultimately determine the range-velocity visibility of a waveform, while ineffective clutter rejection can result in greatly degraded subclutter visibility�
Eclipsed ranges result from turning off or blanking the receiver while transmitting each pulse in the burst since it is difficult to obtain enough isolation to receive radar returns without degraded sensitivity during highpower transmission� Blind velocities occur because the clutter filter frequency response rejects targets that lie within the clutter notch� The effects of eclipsing and blind velocities can be reduced by the use of multiple bursts with different pulse repetition frequencies�
Loss in target visibility due to insufficient clutter rejection is very much dependent on the specific clutter environment� Percent target visibility is a figure of merit that is used in designing clutter rejection waveforms� Some
range-velocity regions may be given a higher priority depending on the application and need to be weighted accordingly� The following paragraphs discuss eclipsing, blind velocities, and visibility in clutter in more detail�
9.2.9 Eclipsing
Eclipsing, in general, may result from at least two sources: (1) receiver blanking during high-power transmission and (2) when the radar antenna has nonreciprocal phase shifters such as in some phased array designs, the receiver is additionally blanked during the time it takes to switch the phase shifters from receive to transmit and then from transmit to receive� The total blanking time determines the eclipsed range interval:
r c te pc= +( )t /2 (9�73)
where tpc is the phase shifter cycling time� Eclipsing occurs each time a pulse in the burst is transmitted�
9.2.10 Blind Velocities
Blind velocities result because the clutter filter response repeats at multiples of the pulse repetition frequency (PRF)� Since the unambiguous velocity is equal to the PRF, the blind velocity interval will repeat at multiples of the unambiguous velocity that is given as
v c PRIfua c= /( )2 (9�74)
where fc is the radar microwave frequency� The blind velocity interval is determined by the width of the notch region of the MTI filter�
9.2.11 Visibility in Clutter
The percent of the instrumented range-velocity space visible for a given waveform is a figure of merit often used to gauge performance� A target is visible if it is not located at an eclipsed range or blind velocity and the signalto-clutter plus noise ratio is sufficient to achieve the required probability of detection� When the MTI waveform consists of more than one burst PRF in order to resolve velocity ambiguities and to effectively reduce or eliminate blind velocities, a particular target range-velocity combination is said to be visible if the requirements are met for one or more of the multiple PRFs� In other words, if the requirement is met for at least one PRF for a given target range-velocity combination, the target is visible�
Since some target range-velocity combinations may have a higher priority for detection than others, such as a high-velocity target at close range as
opposed to a low-velocity target at long range, it may be desirable to weight the range-velocity space� For the general case of nonuniform weighting, the percent coverage can be written as
%
( , ) ( , )
C W r v W r v
v= S S
(9�75)
where Wv are the weights for the visible range-velocity space and the denominator is the sum of all weights over the entire space under consideration�
Various combinations of range and velocity weighting should provide a reasonable amount of flexibility� Three range weighting alternatives are considered� The range weightings are (1) linearly decreasing with range, (2) exponential decreasing, and (3) a combination of uniform and monotonically decreasing weights as shown in Figure 9�13�
The slope for the linearly decreasing range weighting is specifiable such that the range axis intercept point always occurs at the maximum range� For the case of exponential decay weighting, the rate of decay is specifiable� For combinational weighting, a uniform weight is applied out to specifiable range rt and an exponentially decreasing weight is applied from range rt to the maximum range�
The two velocity weightings considered are (1) uniform and (2) linearly increasing with velocity with the velocity intercept occurring at zero velocity� The range and velocity weights are assumed to be separable, so the combined weight is just their product:
W W Wrv r v= (9�76)
This approach allows greater emphasis to be given to high-velocity, shortrange targets that stress the capabilities of the radar and weapon system�
9.3.1 Receiver Noise
Receiver noise is a range-dependent noise source� The effect of receiver noise and the guidance command generation needs to be considered and possibly filtered before seeker angle measurements can be used in a guidance law� Receiver noise comes from a number of different sources, both natural and man-made, and according to Nesline [10], receiver noise is the thermal noise produced at the radar receiver with range dependence governed by the radar range equation as shown in the following equation:
s srn RN
TMR R
= ×0 0
(9�77)
where σRN0 is the receiver noise standard deviation at the reference range, R0, given a receiver noise correlation time constant, TRN0� The white reference receiver noise angular error standard deviation is a function of its power spectral density (PSD) and is found using the following equation:
sRN
02 =
× F
(9�78)
RTM is the time-dependent range between the target and the interceptor missile� Internal noise is primarily introduced in the early stages of the receiver,
in that it passes through all amplifying stages of the receiver and will therefore dominate noise introduced later in the receiver� Other sources of radio signals are the ground, space, etc�, but in airborne applications these are negligible compared to the internal noise [9]�
9.3.1.1 Glint
Much has been written [10-15] on Glint and its associated effects on missile terminal homing accuracy and miss distance� Glint can be defined as the RCS fluctuation-induced apparent LOS angle and angle rate change� Glint is not dependent on the radar or specific tracking techniques (conical scan, monopulse, etc�) and is principally the result of coherent scattering from a complex target (e�g�, body, wings, tail)� When considering the reflections from a complex target, they are not simply from a point but are a collection of point scatterers all contributing to the reflections collected by the tracking radar�
Target motion in the form of maneuvering or stabilization motions (short or long period motions for example) contributes glint by changing the relative phase and magnitude relationship between each independent scatterer� The resultant effect of this motion is to cause an apparent random wander of the radar center of the target being tracked� The wander occurs in both angle and range and will have rate components� Glint is particularly a problem for seekers as it is a low-frequency phenomenon and usually falls in the passband of the servos driving the antenna and cannot be easily filtered� Compounding the problem is that as the intercept range closes to zero, the angle error increases infinitely per lateral fluctuation deviation resulting in the saturation of the servo system�
Glint, for complex targets with many scatters, is statistically characterized by James [16], Barton [17], and Garnell [18] as a Gaussian probability density function (PDF) with the mean and standard deviation of the spectrum tied to the target dimensions� According to Garnell, glint is largely a phenomenon in the target yaw plane and the RMS value can be characterized as 1/5 its wingspan, while James claims values between 1/3 and 1/6 are appropriate� Gordon [15] characterizes glint standard deviation as follows:
sgl
k L R
int = ×
(9�79)
where k (0�1-0�2) is the proportionality constant from empirical study L is the target wingspan R is the range to target
Although Gordon predicts that glint follows a time-correlated, student-t distribution, he concludes that using a Gaussian distribution in an extended
nonlinear Kalman filter is appropriate� Gordon should be consulted for a sophisticated approach to modeling glint�
Garnell proposes to model the glint spectrum as white noise passed through a first-order lag with correlation time constant (Tg) and shown in the following:
F( )w
w =
+ K
T g
2 21 m /rad/s2 (9�80)
where K k Tg g g2 2 2= × ×p / is the mean square value of the glint Tg is typically in the numeric ranges between 0�1 and 0�25
Gordon, referencing Barton, expresses Tg as follows:
T
L g
a @
× × l w3 4.
