and multiple bandpasses. An example scenario is presented where the models are proxy geostationary Earth orbit (GEO) satellites with different bus con...

1 downloads 0 Views 480KB Size

SHAPE, SURFACE PARAMETER, AND ATTITUDE PROFILE ESTIMATION USING A MULTIPLE HYPOTHESIS UNSCENTED KALMAN FILTER Charles J. Wetterer, * Bobby Hunt, † Kris Hamada, ‡ John L. Crassidis § and Paul Kervin ** Multiple hypothesis testing using an underlying Unscented Kalman Filter (UKF) to estimate state parameters has been previously demonstrated where the state includes the attitude, angular rate, position, velocity, and surface parameters of the space object. This algorithm is extended to include multiple components, and multiple bandpasses. An example scenario is presented where the models are proxy geostationary Earth orbit (GEO) satellites with different bus configurations (e.g. cylinder vs. rectangular prism) and in different controlled states (e.g. nadir-tracking vs. inertial).

INTRODUCTION The utility of using brightness (photometric flux intensity) measurements to determine a space object (SO)’s attitude and shape has been long established [1]. Ref. [2] first demonstrated how these brightness measurements can also be used in an estimation filter, such as the unscented Kalman filter (UKF), to estimate the attitude and angular rates of a SO. An estimation filter that fuses both angles (line-of-sight) and brightness measurements for the purpose of orbit, attitude, and shape determination was recently demonstrated [3] and expanded to include simultaneous surface parameter estimation as well [4]. Estimation of the SO’s mass has also been demonstrated [5-6]. Whereas angles measurements are direct observations of the SO’s orbital position, the brightness measurement is dependent on orbital position and the SO’s shape, surface parameters and attitude, and thus provides an indirect observation of these other attributes. The physical correlation between the SO’s shape/attitude and its orbital position are caused by the various nongravitational forces and torques, such as radiation pressures that produce a linear and angular acceleration on the SO. These radiation pressures must be consistent with the surface bidirectional reflectance distribution functions (BRDFs) [7]. The work presented here expands the algorithm of Refs. [3-6] to include multiple bandpass observations and multiple SO components with possibly different attitude profiles. First, the *

Principal Scientist, Integrity Applications Incorporated - Pacific Defense Solutions, 1300 N. Holopono St. #116, Kihei, HI 96753, Email: [email protected] † Integrity Applications Incorporated - Pacific Defense Solutions, Email: [email protected] ‡ Integrity Applications Incorporated - Pacific Defense Solutions, Email: [email protected] § CUBRC Professor in Space Situational Awareness, Department of Mechanical & Aerospace Engineering, University at Buffalo, Amherst, NY 14260-4400, Email: [email protected] ** Air Force Research Laboratory, 535 Lipoa Parkway, Suite 200, Kihei HI 96753.

1

complete list of parameters included in the SO’s “augmented” state corresponding to all parameters that could be included in the filter estimation process is presented. It is important to note that through radiation pressure, the BRDF parameters and other potential parameters in the augmented state (e.g. mass) influence the dynamics that affect both the attitude and orbit and thus influence the evolution of the system from point to point. As with the “classical” state parameters, these new parameters are intrinsic and unique properties of the system. Next, the multiple hypothesis testing using the UKF as the engine is reviewed. Generally, this is referred to as Multiple Model Adaptive Estimation (MMAE) [8], but alternatively is simply an implementation of an Unscented Particle Filter (UPF) [9] where each particle is a different hypothesis. Finally, an example scenario is examined in detail. BUILDING THE AUGMENTED STATE An Earth orbiting SO’s attitude, position, size, shape, mass and surface characteristics are all needed for high fidelity orbit propagation and in calculating associated measurements that are remotely observable, such as angles and brightness. All of these parameters make up the SO’s “augmented” state. In this paper, the quaternion, which is based on the Euler angle/axis parameterization and contains four values denoted by q = [q0 , q1 , q2 , q3 ] , is used to specify the SO’s attitude. For a SO with multiple components (e.g. a 3-axis stabilized satellite with a bus and independently tracking solar panels), a simplistic method where each component’s attitude is specified with a separate quaternion is used. In addition, each component’s angular rate is denoted by ω = ω x , ω y , ω z as defined in the SO’s body frame.

[

]

A single position and velocity corresponding to the SO’s center of mass, denoted by I r = [x, y , z ] and v I = [v x , v y , v z ]I respectively, are used whether the SO has multiple components or not. This is the classic 6DOF representation of the SO’s orbit. I

Figure 1. The Geometry of Reflection. In addition to the mass of each component, the parameters associated with the SO’s surface materials are also contained in the augmented state. Each material that makes up the SO could reflect light differently. Thus, a SO might have many parameters that specify its surface properties. The function that defines how light is reflected from an opaque surface with a given surface normal direction ( Nˆ ), illumination direction ( Lˆ with angles θi and ϕi from Nˆ ), and observer direction ( Vˆ with angles θr and ϕr from Nˆ ) as shown in Figure 1 is called the bidirectional reflectance distribution function (BRDF) and is given by

2

f r (θ i , fi ;θ r , f r ; λ ) =

dLr (θ r , f r ) dE i (θ i , fi )

(1)

where dLr is the reflected radiance in Wm-2sr-1 and dEi is the irradiance in Wm-2. The bisector ˆ = (L ˆ +V ˆ) L ˆ +V ˆ with angles α and β vector between the illumination source and the observer is H from Nˆ . There are many different reflectance models that could be used for the BRDF, such as that of Ashikhmin-Shirley [10] or Cook-Torrance [11]. Shape models were built as the sum of facets where each facet has a position in the body frame and is specified by a particular area and normal vector. The SO’s brightness is also calculated by summing the contribution to the brightness by each facet using the BRDF in the equation

mobject = mSun

N facets ˆ k ⋅ Lˆ N ˆ k ⋅V ˆ ∑ Ak ( f r )k N k =1 − 2.5 log10 r2

(

)(

)

(2)

where Ak is the area of the kth facet, r is the distance between the SO and observer, and mSun is the apparent magnitude of the illumination source (in this case the Sun) as listed in Table 1 for various photometric bandpasses. Table 1. Solar magnitudes bandpass mSun Flux (W/m2) U -25.96 159 B -26.09 171 V -26.74 175 R -27.15 211 I -27.49 270 J -27.3 261 K -25.1 121 In this paper, a simplistic method for accommodating multiple bandpasses is employed. Certain parameters associated with the BRDF model are geometric in nature (e.g. the exponential factors in the Ashikhmin-Shirley model), and so are independent of bandpass. Other parameters require a separate, although related, value for each bandpass (e.g. for the Ashikhmin-Shirley model, this includes the diffuse reflectance and the specular reflectance at normal incidence). ALGORITHM DESCRIPTION MMAE as used in Refs. [3-6] (or equivalently the UPF without resampling) is used in this paper. At the end of each step, a hypothesis’s weight is calculated using a likelihood function and the sum of these weights are normalized to unity. The UKF employed for each hypothesis is that of the quaternion-based UKF of Crassidis and Markley [12-13]. The augmented state of each component of the SO’s orientation and rotation rate, and the SO’s position, velocity, mass, and other parameters associated with the BRDF is given by

[

x = q1T

q TN c

ω1T

ωTN c

(r ) (v ) I T

3

I T

m1 m N c

pT

]

T

(3)

where Nc is the number of components that make up the SO. The Ashikhmin-Shirley BRDF is used in this paper with the exponential factors equal (n = nu = nv) and in two bandpasses (B and V), and so

p = [n ρ B

dB

ρV

dV ]

T

(4)

where ρ is the diffuse reflectance treated independently (assumed equal to the specular reflectance at normal incidence) and d is the diffuse fraction (where s = 1 – d is the associated specular fraction). In practice, the actual state that is estimated is a subset of this augmented state, and unconstrained “proxy values” are used in place of parameters with constraints (e.g. reflectance is restricted to a value between 0 and 1) as done in Refs [4] and [6]. EXAMPLE The particular example is that of discriminating the surface materials, shape and attitude profile of a geosynchronous Earth orbit (GEO) satellite. The orbit was set to exactly geostationary (a = 42364.16932 km, e = 0, i = 0, M0 = 0, ω = 0, Ω = 0). Observations were simulated starting at 2010 Dec 21 at 5:00:00 UT, 1500 observations every 24 s (for a total of 10 hours), with an observation site corresponding to the top of Haleakala on Maui (latitude = 20.71 deg, longitude = 156.26 deg, altitude = 3.0586 km). The GEO is assumed to be one of four shape configurations, as shown in Figure 2.

Figure 2. GEO shape configurations GEOsat1 (upper left) has a bus comprised of a cylinder with 23 facets on its cylindrical surface with a height of z = 0.5 meters and a radius of r = 1 π meters. GEOsat2 (upper right) has a bus that is a rectangular prism of dimensions 1 × 1 × 0.5 meters such that its z-facing surface area is identical to the z-facing surface area of GEOsat1. GEOsat3 (lower left) has a bus comprised of

4

an ellipsoid with 25 facets, and the 10 facets with a +z-axis component have the same projected area as GEOsat1 and GEOsat2. GEOsat4 (lower right) has a bus identical to GEOsat2, except for the addition of a small “antenna” that is tilted 15o towards the +x-axis. Each shape has identical solar panel arrays tilted ±10o along the x-axis. For all shape configurations, the satellite bus mass is 900 kg while the solar panels mass is 100 kg. The satellite bus (and antenna for GEOsat4) is assumed to be comprised of a uniform material, and so all facets share common BRDF parameters. Likewise, the solar panels are assumed to be a different, uniform material, and so have a separate set of BRDF parameters. The Ashikhmin-Shirley BRDF model is used in the B and V bandpasses with “truth” parameters listed in Table 2. Table 2. BRDF truth parameters bus Solar panel nu = nv 100 500 ρB = (F0)B 0.40 0.10 dB 0.70 0.10 ρV = (F0)V 0.50 0.15 dV 0.60 0.05 Each GEOsat is placed in one of three different attitude/rate profiles. Profile #1 corresponds to the satellite having a single component that is inertial pointing with the z-axis aligned towards the Sun and the y-axis aligned close to the north celestial pole (NCP) such that the solar panels and bus glints are observed during the observation period. Profile #2 corresponds to the satellite having two components with the bus tracking nadir (rotation about the y-axis at a rate of 2π/86400 rad/s) and the solar panels inertial pointing as in Profile #1. The bus’ attitude in Profile #2 exactly matches the bus’ attitude in Profile #1 at the time of zero phase angle when the bus glints. Profile #3 corresponds to the satellite having one component that is nadir-tracking, where again the attitude exactly matches the attitude in Profile #1 at the time of zero phase angle.

Figure 3. Light curves for Model #1 (upper left) to Model #6 (lower right)

5

Figure 4. Light curves for Model #7 (upper left) to Model #12 (lower right) There are thus twelve different shape/profile combinations possible (e.g. Model #1 = GEOsat1/inertial, Model #2 = GEOsat1/nadir-inertial, Model #3 = GEOsat1/nadir, Model #4 = GEOsat2/inertial, etc…). The corresponding light curves for each model using the truth parameters in Table 1 are shown in Figs 3 and 4. A measurement noise of 0.1 magnitudes has been added to each light curve. At first glance, many of these light curves look identical. The most pronounced differences are between those models with the GEOsat3 shape and all the others. This is because GEOsat3 does not have an exactly +z-axis facing side that produces a glint at zero phase angle. For the other light curves, there are, however, still small differences corresponding to the different shapes and the different profiles. For example, Figure 5 plots the difference between Model #2 (GEOsat1/nadir-inertial) and Model #5 (GEOsat2/nadir-inertial). The only difference occurs at the beginning of the light curve where the sides of the bus are illuminated and visible and thus contribute to the brightness measurements.

Figure 5. Difference in light curves for Model #2 and Model #5

6

Figure 6 plots the difference between Model #5 and Model #11 (GEOsat4/nadir-inertial). Here, the only difference is due to the contribution of the small antenna that produces a glint midway through the pass.

Figure 6. Difference in light curves for Model #5 and Model #11 Figure 7 plots the difference between Model #11 and Model #12 (GEOsat4/nadir). Here, the glints associated with the nadir-tracking solar panels occur slightly earlier (for the first glint) and slightly later (for the second glint) than the inertial solar panels, giving the magnitude difference a distinctive signature.

Figure 7. Difference in light curves for Model #5 and Model #11 In this example, it was assumed the solar panel surface parameters are known but the bus surface parameters are unknown. As such, three different hypotheses for each of the five unknown surface parameter were used corresponding to 35 = 243 different hypotheses for each base model. So, hypotheses 1-243 use Model #1, hypotheses 244-486 use Model #2, and so forth. This yields a total of 2916 hypotheses in our estimation filter. The attitude and angular rates for each component and the BRDF surface parameters of the bus were included in the estimation state. The orbital position and velocity are not in the state, and so the orbit is assumed to be perfectly known. Additionally, the radiation pressure forces and torques and any internal and intercomponent torques are assumed to be negligible. The measurement function of Eq. (26) includes only the two brightness measurements for the B and V bandpasses. The initial state and initial uncertainty used in creating the covariance matrix for each of the hypotheses are shown in Table 3. For the attitude and angular rate, it is assumed the values are known to a high degree of precision, and the truth of that particular hypothesis (depending on the

7

profile) is used as the initial state. For the BRDF parameters in the state, multiple values are listed corresponding to the specific values used in the various hypotheses. Table 3. BRDF truth parameters state uncertainty truth 0.067 deg q truth 6.7x10-5 deg/s ω nu = nv [100 200 300] 5 ρB = (F0)B, ρV = (F0)V [0.1 0.5 0.9] 0.1 dB, dV [0.1 0.5 0.9] 0.1 Two tests were run. In the first test, Model #6 (GEOsat2/nadir) was used to generate the observations. The closest light curve match is between Model #6 and Model #3 where there is a slight difference caused by the illumination of the side of the cylindrical bus versus the rectangular bus similar to that shown in Figure 5. Figure 8 plots the weights of the 2916 hypotheses as a function of time over the entire pass. Red lines correspond to hypotheses of the correct Model (hypotheses 1216-1458). By the end of the pass, a single hypothesis (#1265) of the correct shape and attitude profile has all the weight. Further, the BRDF parameters of this hypothesis at the final time step, shown in Table 4, have converged to the correct values with a reasonable uncertainty.

Figure 8. Weights as function of time for test #1 (Model #6 truth) Table 4. BRDF parameter test #1 results truth

Estimated #1265 nu = nv 100 101.8 ± 2.1 ρB = (F0)B 0.40 0.398 ± 0.003 dB 0.70 0.705 ± 0.004 ρV = (F0)V 0.50 0.502 ± 0.005 dV 0.60 0.610 ± 0.003 In the second test, Model #11 (GEOsat4/nadir-inertial) was used to generate the observations. The closest light curve match is between Model #11 and Model #5 where there is a slight difference caused by the antenna as shown in Figure 6. Figure 9 plots the weights of the 2916 hypothe-

8

ses as a function of time over the entire pass. Red lines correspond to hypotheses of the correct Model (hypotheses 2431-2673). By the end of the pass, a single hypothesis (#2480) of the correct shape and attitude profile has 98.2% of the weight. Again, the BRDF parameters of this hypothesis at the final time step, shown in Table 5, have converged to the correct values with a reasonable uncertainty.

Figure 9. Weights as function of time for test #2 (Model #11 truth) Table 5. BRDF parameter test #2 results truth nu = nv ρB = (F0)B dB ρV = (F0)V dV

100 0.40 0.70 0.50 0.60

Estimated #2480 102.6 ± 2.3 0.396 ± 0.003 0.707 ± 0.005 0.498 ± 0.006 0.611 ± 0.005

CONCLUSION Multiple hypothesis testing with the unscented Kalman Filter (UKF) is a powerful estimation tool. The algorithm as used in previous work was extended to include multiple bandpasses and multiple space object (SO) components, and the ability to pick and choose what parameters, including those associated with the SO’s surface properties, are included in the estimation state. An example demonstrating the discrimination between various shape models and various attitude profiles was also presented. ACKNOWLEDGEMENTS This work was conducted under the Integrated Supervised and Unsupervised Processing (ISUP) Phase I Small Business Innovative Research (SBIR) contract FA9451-12-M-0311 sponsored by the Air Force Research Laboratory (AFRL). Some of the material was also produced along with the SBIR Phase II contract FA9453-11-C-0154, also sponsored by AFRL.

9

REFERENCES [1] Calef, B., Africano, ,J. Birge, B., Hall, D. and Kervin, P., "Photometric Signature Inversion," Proc. SPIE 6307, 63070E, 2006 [2] Wetterer, C.J. & Jah, M., “Attitude estimation from light curves,” Journal of Guidance, Control, and Dynamics 32: 1648-1651, 2009 [3] Linares, R., Crassidis, J. L., Jah, M. K., and Kim, H., “Astrometric and Photometric Data Fusion for Resident Space Object Orbit, Attitude, and Shape Determination via MultipleModel Adaptive Estimation,” AAS/AIAA Astrodynamics Specialists Conference, Vol. AIAA8341, 2010. [4] Wetterer, C. J., Crassidis, J. L., Linares, R., Chow, C. C., and Jah, M. K., “Simultaneous Position, Velocity, Attitude, Angular Rates, and Surface Parameter Estimation Using Astrometric and Photometric Observations,” Proceedings of the FUSION Conference, Istanbul, Turkey, July 2013 [5] Linares, R., Jah, M. K., Leve, F. A., Crassidis, J. L., and Kelecy, T., “Astrometric and Photometric Data Fusion For Inactive Space Object Feature Estimation,” Proceedings of the International Astronautical Federation, Cape Town, South Africa, Sept. 2011, Paper ID: 11340. [6] Linares, R., Crassidis, J. L., Wetterer, C. J., Hill, K., and Jah, M. K., “Astrometric and Photometric Data Fusion for Mass and Surface Material Estimation using Refined Bidirectional Reflectance Distribution Functions-Solar Radiation Pressure Model,” Advanced Maui Optical and Space Surveillance (AMOS) conference, Wailea, Maui, HI, September 2013 [7] Wetterer, C.J., Linares, R., Crassidis, J. L., Kelecy, T.M., Ziebart, M.K., Jah, M.K., and Cefola, P.J., “Refining Space Object Radiation Pressure Modeling with Bidirectional Reflectance Distribution Functions,” Journal of Guidance, Control, and Dynamics, submitted, 2012 [8] Magill, D. T. Optimal adaptive estimation of sampled stochastic processes. IEEE Trans. Auto. Contr., 10:434–439, 1965. [9] van der Merwe, R., Doucet, A., de Freitas, J.F.G, and Wan, E., “The Unscented Particle Filter,” Technical Report CUED/F-INFENG/TR 380, Cambridge University Engineering Department, 2000. [10] Ashikhmin, M., and Shirley, P., "An Anisotropic Phong BRDF Model," Journal of Graphics Tools, Vol. 5, No. 2, 2000, pp.25-32. [11] Cook, R. L., and Torrance, K. E., “A reflectance model for computer graphics,” ACM Transactions on Graphics, 1(1):7–24, Jan. 1982. [12] Julier, S. J., Uhlmann, J. K., and Durrant-Whyte, H. F., “A New Approach for Filtering Nonlinear Systems,” Proceedings of the American Control Conference, Seattle, WA, 1628– 1632, June 1995 [13] Crassidis, J. L. and Markley, F. L., “Unscented Filtering for Spacecraft Attitude Estimation,” Journal of Guidance, Control and Dynamics 26, 536–542, 2003

10