NEW INSIGHTS OF ASTEROID 4179 TOUTATIS USING CHINA CHANG'E-2 CLOSE FLYBY OPTICAL MEASUREMENTS

After departing from the Sun–Earth L2 point, Chang'e-2 had been flying for approximately eight months before flying by Toutatis. The flyby happened at 8:30 (UTC) on 2012 December 13, about one day after Toutatis's periodical flyby of the Earth at 6:40 (UTC) on 2012 December 12. The minimum flyby range from the Earth is 0.046 AU, about 18 times the average distance from the Earth to the Moon. At the moment of flyby, Chang'e-2 was more than 6.9 million kilometers away from the Earth.

To ensure a successful flyby, the Beijing Aerospace Control Center (BACC) had corrected the trajectory of Chang'e-2 several times before the detection, the last of which greatly shortened the theoretical flyby range from hundreds of kilometers to several kilometers, and obtained a relative velocity of 10.73 km s⁻¹. During the flyby, the ground observation systems and control systems were used to track and measure the movement of Chang'e-2 and Toutatis. The spacecraft orbit is tracked through the radiometric data, and the asteroid's orbit is determined through the ground optical telescopes.

A CMOS camera with a focal length of 54 mm, a resolution of 1024-by-1024 pixels, and a field of view (FOV) of 72 by 72 was used in the experiment. After Chang'e-2 departed the Sun–Earth L2 point, imaging tests and parameter calibration were conducted using fixed stars. Chang'e-2 started to take photos of Toutatis after passing by it with a mode of vanishing point gaze (Huang et al. 2013b); the effective imaging time was approximately 100 s, and the imaging frequency was 5 frames per second. Several hours before the flyby, the attitude of Chang'e-2 was adjusted to make sure Toutatis came into the view of the camera and that communication with the Earth was well established. In addition, the constraint on the solar panels–Sun angles was met. The orientation of the Chang'e-2 spacecraft remained fixed inertially throughout the imaging campaign.

When imaging an asteroid with a space-borne camera, the orientation of Toutatis is a function of the relative state/dynamics, which are the position/velocity and the rotation between Chang'e-2 and Toutatis. Strict motion relations of relative displacement and relative rotation based on orbit data provide possibilities to high-precision optical measurement and analysis.

2.2.1. Relative Displacement

2.2.2. Relative Rotation

Three types of coordinate systems are involved herein, namely, an inertial coordinate system, a camera coordinate systems, and a spacecraft body coordinate system. For the convenience of the following analysis, they are uniformly described here. In an inertial space, various tracking data, including the coordinates and velocities of the Sun, the Earth, Toutatis, and Chang'e-2, are unified into the heliocentric inertial system O_i − xyz. The relative motion relation between Toutatis and Chang'e-2 as well as the relative rotational dynamics are described in the original camera coordinate system $O_{c_0 } - xyz$. The optical center of the space-borne camera at the moment of initial imaging is the origin $O_{c_0 }$, and the optical axis of the camera at the moment is the z-axis. Real-time camera coordinate systems $O_{c_k } - xyz$ are built for each imaging position, the origins of which are located at the camera realtime optical centers, respectively. Here, the original camera coordinate system is distinguished from the real-time camera coordinate systems for the former, which is regarded as a fixed frame in the following vision-based calculation. The conversion between the original and real-time camera coordinate systems can be achieved via constraint-based visual calculation. The conversion between the heliocentric inertial system and the camera coordinate systems relies on the body coordinate system O_b − xyz, that between the inertial coordinate system and body coordinate system is based on the attitude transfer matrix, and that between the spacecraft body coordinate system and camera coordinate systems is achieved based on the camera installation matrix. The spatial relationships among the three types of coordinate systems are shown in Figure 1.

Zoom In Zoom Out Reset image size

Assuming that there is no light transmission on Toutatis's surface, when the imaging distance is long, lines of sight and solar rays can, respectively, be considered to be parallel vectors. When Toutatis is treated as a rotating body around the long axis, the visual range of the surface of the equivalent rotating body can be described by the dihedral angle θ formed by the interface of Sun rays and that of lines of sight. The spatial geometric relation of θ is shown in Figure 2, where n₁ are the normal vectors of the boundary plane of Sun rays, n₂ the normal vectors of the boundary plane of lines of sight, and ${\bf l}_{\rm sun}^{c_0 }$, ${\bf l}_{\rm sight}^{c_0 }$, and ${\bf l}_{4179}{^c_0 }$ are the vectors of Sun rays, lines of sight, and Toutatis's long axis in the original camera coordinate system, respectively. The computation of θ is described in Appendix A.1.

Zoom In Zoom Out Reset image size

When the orientation of the long axis is limited in a rough degree range, the dihedral angle θ of the visible surface of the rotating body was calculated to be from 122° to 180°. From the spatial geometric relationship, the conclusion can be thereby made that the long axis of Toutatis can fortunately be completely reflected on the Chang'e-2 optical images, whereas the edge of the upper right part perpendicular to the long axis is partly shadowed.

Based on data from several radar observations since 1992, Takahashi et al. (2013) studied Toutatis's inertia moment and verified the high-resolution shape model (Hudson et al. 2003) using the observation in 2012. The shape model provides a good reference for marking Toutatis's long axis in optical images captured by Chang'e-2. First, based on the solar line-of-sight angle, the illuminated surface of the shape model was simulated. Then, the attitudes of the shape model's three axes were rotated to achieve the most correspondence with the observed optical images of Toutatis. Then, the first order estimate of the direction of the axis of Toutatis was identified in the optical images. It is emphasized here that the long axis cannot be extremely precisely identified due to uncertainties in the shape model and the resolution of the optical images. In this paper, the two landmarks A and B along the long axis are considered the endpoints of the long axis. Figure 3 shows the attitude and landmarks, respectively, in the shape model and Chang'e-2 image.

Zoom In Zoom Out Reset image size

The relative motion between Chang'e-2 and Toutatis leads to a special imaging relationship. That is, all lines between homonymy points at different spatial coordinates intersect at an infinite point P_∞ that indicates the translational direction in space. Furthermore, on the image planes, the lines through homonymy landmarks at different positions intersect at one point, epipole e, which is the projection of P_∞. Based on the relation between the image sequence and the spatial motion, the constraints and the computational algorithm can be developed.

During the imaging period, the principal point coordinate and the focal length of the space-borne camera are deemed as known and stable with even exposure at each time of imaging. According to the observing geometry, three strict constraints can be proposed for accurate modeling.

The scale constraint. During any two times of continuous imaging, the relative motion of the spacecraft and Toutatis is unchanged. Therefore, Δa_{m, n}, the relative displacement vector of the camera between the picture ${Img}_m$ and the picture ${Img}_n$, can be expressed as

$$\begin{equation} \Delta {\boldsymbol a}_{m,n} = \left({n - m} \right) \cdot {\bf v} \cdot \delta t, \end{equation} \tag{ 1 }$$

where δt is the time interval between two camera exposures that could be accurately controlled by the high-precision clock equipment on the spacecraft, and v is the relative velocity vector between Chang'e-2 and Toutatis. When tracking Chang'e-2 using ground communication, errors, including systematic error and random error, are inevitable. Nevertheless, the measurement accuracy of the relative velocity in a short time could reach the mm s⁻¹ level, which provides the baseline accuracy for the whole measurement system.

The attitude constraint. According to the relative rotation relation between Chang'e-2 and Toutatis, if a unit matrix I_{[3 × 3]} of the first image Img₀ is established as the reference, the camera rotation matrix R_k of any other image Img_k could also be uniformly expressed as

$$\begin{equation} R_k = I_{\left[ {3 \times 3} \right]}. \end{equation} \tag{ 2 }$$

The epipole constraint. Theoretical values of e_{m, n}, the epipole coordinate, between any two optical images Img_m and Img_n in the sequence are identical and can be uniformly denoted by

$$\begin{equation} {\bf e}_{m,n} = {\bf e} = [e_u,e_v ]^T. \end{equation} \tag{ 3 }$$

Incorporating the scale constraint, the attitude constraint, and the epipole constraint into the solution of multiple view geometry, the coordinate of a spatial point P in the original camera coordinate systemcan be obtained. Appendix A.2 gives a detailed description of the solution. The camera matrices P_m and P_n can be obtained by combining with a pure translational imaging model. The specific process refers to the multi-view geometry principle (Hartley & Zisserman 2004). After getting the camera matrix of each optical frame and coordinates of homonymy image points of the spatial point P, the exact spatial position of P in the original camera coordinate system can be obtained.

According to the above geometric relationship, the computation accuracy is directly determined by the accuracy of the coordinates of key feature points, including the homonymy landmarks and the epipole. Although the measurement accuracy of each individual point is low, the specific imaging model and constraint relations during Chang'e-2's flyby of Toutatis could help to improve the estimation accuracy of the feature points. The epipole can be estimated using the two-point algorithm with the least squares criterion (Li et al. 2000). Accurate registration of homonymy landmarks can be realized, according to a motion constraint based filter algorithm. The following first theoretically gives a motion constraint within the sequence of homonymy image landmarks. Then, based on the theoretical relation, the pixel coordinates of the homonymy landmarks are estimated accurately using a polynomial-based filter.

The motion constraint for the filter. The 2D image coordinate [u_P(t), v_P(t)]^T of any spatial point P on Toutatis's surface and its spatial 3D coordinate ${\bf X}_P^{c_0 } \left(t \right) = \left[ {x_P \left(t \right),y_P \left(t \right),z_P \left(t \right)} \right]^T$ in the original camera coordinate system meet the following continuous time-varying relationship, the detailed derivation of which is given in Appendix A.3.

$$\begin{equation} u_P \left(t \right) = \frac{{f_u x_P \left(t \right) + u_0 z_P \left(t \right) - \frac{{Vt}}{K}f_u f_v e_u }}{{z_P \left(t \right) - \frac{{Vt}}{K}f_u f_v }} \end{equation} \tag{ 4 }$$

$$\begin{equation} v_P \left(t \right) = \frac{{f_v y_P \left(t \right) + v_0 z_P \left(t \right) - \frac{{Vt}}{K}f_u f_v e_v }}{{z_P \left(t \right) - \frac{{Vt}}{K}f_u f_v }} \end{equation} \tag{ 5 }$$

$$\begin{equation} {K} = \sqrt {({f_v ({e_u - u_0 })})^2 + ({f_u ({e_v - v_0 })})^2 + ({f_u f_v })^2 }, \end{equation} \tag{ 6 }$$

where u_P(t) and v_P(t) are the pixel coordinates of the row and column direction at time t, f_u and f_v are the camera equivalent focal lengths in the two directions, u₀ and v₀ are coordinates of the camera principle point, and V = ||v|| represents the relative speed.

The relationship between the homonymy landmarks [u_i, v_i]^T, i = 1, 2, ..., n of each image can be accurately established through Equations (4)–(6), which lay the theoretical foundation for the accuracy of the homonymy landmarks. However, since ${\bf X}_P^{c_0 } (t)$ is unknown, the exact coordinates of the landmarks cannot be calculated via Equations (4)–(6) directly. Fortunately, there was a favorable condition that ${\bf X}_P^{c_0 } (t)$ was constant in the original camera coordinate system. Therefore, to analyze the time-varying rule of the image coordinates of homonymy landmarks, different constant values can be assigned to ${\bf X}_P^{c_0 } (t)$. Using the scale constraint, the image points x_k = [u_P[k], v_P[k]]^T, k = 1, 2, ..., n change continuously and stably with time sample k, which can be approximated with a polynomial.

The filter for the homonymy landmarks was designed according to the rule of approximate polynomials. First, ${\bf X}_P^{c_0 } (t)$ was set as a constant, and the theoretical formula in Equations (4)–(6) was used to precisely calculate pixel coordinates of the image point sequence x_k; then, a high-order polynomial was used to fit x_k, and the values of the polynomial functions were used as an estimate of x_k. Figure 4 shows a comparison of mean square errors of image coordinates of different sequence lengths via polynomials of different orders. Each statistical point was based on 200 random sample points of ${\bf X}_P^{c_0 } (t)$. The figure shows that for sequences of the same length, the error converges as the order of the polynomial is increased; for sequences of the same fitting order, the error diverges with an increase in the sequence scale. When the order of the polynomial is above 7, in the sequence scale of 60 points, the mean square error is less than 0.004 pixel. Therefore, fitting polynomials with orders of more than seven can be used as filter functions for an accurate correction of the homonymy landmarks.

Zoom In Zoom Out Reset image size

The aggregate uncertainties of the optical imaging during flyby need to be addressed. The source errors include uncertainties in the relative dynamics of Toutatis and spacecraft and the measurement models. The uncertainties in the relative dynamics are determined by that of the motion relation in Section 2. Their effects on the measurement system have proved to be negligible. Therefore, the measurement uncertainty is primarily determined by the uncertainty of the measurement model.

According to the measurement model, the main errors include the matching error of homonymy landmarks, the estimation error of the epipole, and transmission errors. To precisely analyze the degree of error, a ground testing platform is developed based on the relative imaging relationship and actual parameters including relative velocity, exposure time, and camera parameters.

The ground testing platform, on one hand, could generate theoretical spatial coordinates of simulated landmarks and corresponding homonymy image points at each exposure time. On the other hand, after adding noise to the simulated data, it can realize a series of processes including filtering and spatial position calculating, as well as a comprehensive analysis of transfer relations of various error sources.

The spatial coordinates of Toutatis's endpoints are basic landmarks for calculating the dimension and the orientation of the long axis. Based on the exact measurement model and Monte Carlo method, the uncertainty of the endpoints can be obtained statistically via the ground testing platform after estimation of the epipole, filtration of the homonymy image points, and error transmission.

The cumulative variance of a parameter can be computed as shown in the following equation. The uncertainty of the dimension d_AB and the orientation α_AB can be obtained from the coordinates of endpoints and their uncertainties

$$\begin{equation} \sigma _{\tau _{AB} } ^2 = \sum\limits_{i = 1}^6 {\left({\frac{{\partial \tau _{AB} }}{{\partial p_i }} \cdot \sigma _{pi} } \right)^2 }, \end{equation} \tag{ 7 }$$

where τ_AB represents the function of Euclidean distance d_AB or the orientation angle α_AB between endpoint A and B, p_i denotes the coordinates of A and B in the original camera coordinate system, and σ_pi denotes the variance of each coordinate.

Based on the above mathematical formulation, the optical measurement of Toutatis's long axis is achieved as follows. First, the endpoints of the long axis are identified. Then, the extracted endpoints are matched on the sequence images. In order to obtain enough precision under a limited image resolution, a filter is invoked. Then, the spatial coordinates of the endpoints are computed under the constraints. Finally, the dimensions and the orientation of the long axis on the optical images are calculated.

A total of 30 continuous images with relatively high resolution were used for final processing. The initial matching of homonymy image points is achieved on the biquadratic interpolation images by virtue of corner-based matching algorithms, ensuring a maximum error not greater than 1 pixel and a standard deviation not greater than 0.35 pixel. After filtering of homonymy points' coordinates based on the motion constraint, the mean error declines from 0.39 pixel to 0.06 pixel, and the standard deviation declines from 0.4 pixel to 0.05 pixel. The spatial position error of the endpoints caused by matching errors could be calculated statistically on the basis of the measurement transfer relation. After filtering, the order of magnitude of the calculated standard deviation of the spatial position decreases from hundreds of meters to ∼10 m. When errors are described by the rms, errors on the x, y directions in the original camera coordinate system are ∼2 m and that on the z direction is ∼60 m; thus, the error is mainly reflected in the z direction, namely, the depth direction. Table 1 shows estimated coordinates of endpoints A and B before and after filtering.

The epipole is obtained from several pairs of random images. In order to ensure the reliability of the result, the image pairs with a small imaging distance and fewer homonymy points are filtered out and matching operators with scale invariant characteristics are adopted. Table 2 is the epipole estimation result in three color channels of the RGB format based optical images. Based on the transfer relation, the epipole estimation error of one times the STD leads to less than 30 m of calculation errors for the endpoints.

Figure 5 shows the measurement results of identified endpoints A and B of six sequential frames arbitrarily selected from the above 30 continuous images under UTC time coordinates. Table 3 shows the measurement results of typical dimension features and the orientation of the long axis. Typical dimensions include the maximum dimension and the length of the two lobes by taking the long axis as the reference distance. The orientation of the long axis is described through the angles of direction cosine between the long axis and three coordinate axes of the original camera coordinate system. Based on the analysis in the analytical conclusion of Section 3, due to the shadowing of the images in the direction perpendicular to the long axis, the relevant dimension in this direction could be measured on the optical images but does not represent the real scale of Toutatis.

Zoom In Zoom Out Reset image size

The dimension of the long axis measured using the optical method is 4354 m with an uncertainty of 56 m, a difference of 246 m compared to 4.6 ± 0.1 km (Hudson et al. 2003) and 106 m compared to 4.46 ± 0.1 km (Takahashi et al. 2013; Huang et al. 2013a). The longest dimension of Toutatis reflected by the image is 4391 m, and the error bar of 56 m is shorter than 100 m. The reasons for the difference include: (1) there exist deviations between identified long axis endpoints and the actual positions, whether they are based on the shape model or extracted from the optical image; (2) there may be a little shadow shelter at the bottom right of the optical images, and point B is slightly offset from the actual end point; and (3) there do exist certain deviations as to the measurement results of these two methods. The dimension of 4391 ± 56 m agrees with the length range of 4.75 km ± 10% estimated by Huang et al. (2013a), while the error bar has been greatly improved. The optical measurements presented here fully integrate much useful information and constraints of Chang'e-2 flying by at a close distance. Therefore, the measurement accuracy of Toutatis's dimensions is improved by comparing it with ground-based remote observations. Dimensions of Toutatis's three axes can be estimated through a proportional relation between the results of the optical measurement and that of the refined shape model, which is 4354 × 1835 × 2216 ± ∼ 56 m. This gives us more accurate estimate of the size of the asteroid.

The orientation of Toutatis's long axis in the original camera coordinate system has been estimated through the canny-based matching between Chang'e-2 optical images and images from the shape model by Zou et al. (2014), which is (− 338, 330, 471) ± 1° in 3-1-2 Euler angles, and (12791, 12300, 12481) ± 1° in the expression of direction cosine angles in the original camera coordinate system. Comparing it with the orientation of (12613 ± 029, 12298 ± 021, 12663 ± 046) calculated through optical approaches in Table 3, the deviations of the three direction cosine angles are (− 178, − 002, 182), respectively. This, on one hand, shows that Toutatis's attitudes measured through the two methods are basically consistent. On the other hand, through further comparison, some obvious differences could still be found between the optical images and the radar-based images, such as the tail of the larger lobe. Thus, some errors are introduced if considering the shape model based image as the reference for the canny operator based matching, especially in the depth direction. Therefore, the discrepancies in orientation between the two methods are insignificant overall, but the deviation of α_{l − z} is larger comparatively. The spatial orientation contains rotation information, which is discussed in Section 6.2.

Based on the accurate coordinates of endpoints A and B, the shadow in the optical images could be precisely calculated. The dihedral angle of the visible part between the boundary plane of Sun rays and that of lines of sight calculated through the long axial vector and Sun ray vector is 12998, as shown in Figure 6(a). That is to say, if we consider a cylinder of equivalent size, about 27.78% of the body is hidden due to shadowing. Figure 6(b) is a prediction of the shadow in the optical images combined with a simulation image based on the shape model. Comparing Figure 6 and the previous Figure 3(a), the calculated shadow is consistent with the simulation result based on the shape model.

Zoom In Zoom Out Reset image size

Toutatis's dimension can be obtained through the optical method, but due to incomplete coverage of Chang'e-2 images, the volume of the asteroid cannot be estimated directly. Fortunately, Ostro and Hudson's radar based shape model provides more accurate relative information on the shape. According to the cube law, we can scale Toutatis's volume based on the shape model by adjusting axes of Toutatis from 4.46 × 1.88 × 2.27 km to 4354 × 1835 × 2216 m. The new estimated volume is 7.5158 km³, 0.1830 km³ less than 7.6988 km³ by the refined radar based shape model. In addition, the newest radar observation and Chang'e-2 images show that Ostro and Hudson's shape model corresponds to Toutatis's actual shape over ∼97% of its volume (Takahashi et al. 2013), so if a more accurate shape model can be established based on the radar observation of 2012, the estimation of the volume will be more accurate.

Furthermore, based on the estimated volume, the mass of Toutatis can be recalculated. According to analysis results from the radar echo, Toutatis's density is between 2100 kg m⁻³ and 2500 kg m⁻³ (Scheeres et al. 1998; Ostro et al. 1999; Birlan 2002), so the mass is between 1.578 × 10¹³ kg and 1.879 × 10¹³ kg. Naturally, if the exact mass of the asteroid can be estimated according to other approaches, such as the Yarkovsky effect (Chesley et al. 2003), an accurate density can be calculated based on the volume, and the bulk composition and internal makeup of Toutatis can be further analyzed.

The equivalent rotating body at the time of the Chang'e-2 flyby can be simulated in the earth plane of sky as shown in Figure 7(a). Comparing this with the ground radar observations on 2012 December 12 as shown in Figure 7(b)⁹, it can be seen that endpoint A at the head of the smaller lobe, point C at the neck, and the sheltered endpoint B on the optical image correspond well to those on the radar image, indicating that the two images are identical in spatial orientation. Nevertheless, there exists a certain deviation as to the long axis orientations in the two figures. One reason for this deviation is the extracting error of the long axis endpoints. Another is that the imaging time of the optical image was about half a day later than that of the radar observation; therefore, after the radar observation in Figure 7(b), Toutatis had spun further, as described by Hudson & Ostro (1995).

Zoom In Zoom Out Reset image size

Furthermore, such disparity provides an opportunity for an accurate analysis of the spin period of Toutatis. The spin period can be calculated by taking the average rotation rate over the observation period. The change in attitude can be obtained by directly comparing landmarks such as A, B, and C in optical and radar images. If the images from the radar observation around 2012 December 13 can be obtained, then corresponding feature points of Toutatis can be found. Such an analysis would help improve the existing rotational dynamical models of Toutatis.

In the original camera coordinate system, let the line-of-sight vector from the spacecraft to the body be ${\bf l}_{\rm sight}^{c_0 }$, the vector of Sun rays be ${\bf l}_{\rm sun}^{c_0 }$, and the vector of Toutatis's long axis be ${\bf l}_{4179}^{c_0 }$; then, the visible surface of the rotating body can be calculated using the following three steps.

First, in the original camera coordinate system, ${\bf l}_{\rm sun}^{c_0 }$ is determined by ${\bf X}_{\rm sun}^{c_0 }$, the Sun centroid coordinate, and ${\bf X}_M^{c_0 }$, the equivalent rotating body center, as

$$\begin{equation} {\bf l}_{\rm sun}^{c_0 } = {\bf X}_{\rm sun}^{c_0 } - {\bf X}_M^{c_0 }. \end{equation} \tag{ A.1.1 }$$

Then, define plane S₁ as the plane that passes the main axis ${\bf l}_{4179}^{c_0 }$ and is co-planed with the vector of Sun rays ${\bf l}_{\rm sun}^{c_0 }$ and calculate the normal vector ${\bf n}_{_{sun - l_r } }^{c_0 }$ of the plane S₁ with the following equation:

$$\begin{equation} {\bf n}_{_{sun - l_r } }^{c_0 } = {\bf l}_{\rm sun}^{c_0 } \times {\bf l}_{4179}^{c_0 }.\end{equation} \tag{ A.1.2 }$$

Then, define plane $S_1 ^ \bot$ as the plane that passes the main axis ${\bf l}_{4179}^{c_0 }$ and is perpendicular to the plane S₁, and calculate the normal vector ${\bf n}_{_{\left({sun - l_r } \right) \bot } }^{c_0 }$ of the plane $S_1 ^ \bot$ using the equation

$$\begin{equation} {\bf n}_{_{\left({sun - l_r } \right) \bot } }^{c_0 } = {\bf n}_{_{\left({sun - l_r } \right)} }^{c_0 } \times {\bf l}_{4179}^{c_0 }. \end{equation} \tag{ A.1.3 }$$

Finally, calculate the plane $S_1 ^ \bot$, which passes the center point ${\bf X}_M^{c_0 }$, and with ${\bf n}_{_{\left({sun - l_r } \right) \bot } }^{c_0 }$ as its normal vector, using the equation

$$\begin{equation} \left({{\bf n}_{_{({sun - l_r }) \bot } }^{c_0 } } \right)^T \left({{\bf X}^{c_0 } - {\bf X}_M^{c_0 } } \right) = 0. \end{equation} \tag{ A.1.4 }$$

Figure 8 is the sketch map of the boundary plane of Sun rays.

Zoom In Zoom Out Reset image size

Figure 9 is the sketch map of the boundary plane of the line of sight. The boundary plane of lines of sight is calculated in the original camera coordinate system. Strictly speaking, any position on the Toutatis surface can correspond to a vector of the line of sight, yet its analysis and calculation are extremely complex. Because the distance from the spacecraft to Toutatis at the imaging time is much larger than the dimension of Toutatis itself, the vector of lines of sight ${\bf l}_{\rm sight}^{c_0 }$ can be represented by the parallel beams passing through the origin and the center ${\bf X}_M^{c_0 }$ as

$$\begin{equation} {\bf l}_{{\rm sight}}^{c_0 } = {\bf X}_M^{c_0 }. \end{equation} \tag{ A.1.5 }$$

Zoom In Zoom Out Reset image size

Define S₂ as the plane that passes the cylinder main axis ${\bf l}_{4179}^{c_0 }$ and co-planes with the line of sight (directory vector ${\bf l}_{{\rm sight}}^{c_0 }$), and calculate the normal vector ${\bf n}_{_{sight - l_r } }^{c_0 }$ of the plane S₂ according to the equation

$$\begin{equation} {\bf n}_{_{sight - l_r } }^{c_0 } = {\bf l}_{\rm sight}^{c_0 } \times {\bf l}_{4179}^{c_0 }.\end{equation} \tag{ A.1.6 }$$

Then, calculate the normal vector ${\bf n}_{_{\left({sight - l_r } \right) \bot } }^{c_0 }$ of $S_2 ^ \bot$, which is defined as the plane passes the cylinder main axis ${\bf l}_{4179}^{c_0 }$ and is perpendicular to S₂:

$$\begin{equation} {\bf n}_{_{\left({sight - l_r } \right) \bot } }^{c_0 } = {\bf n}_{_{\left({sight - l_r } \right)} }^{c_0 } \times {\bf l}_{4179}^{c_0 }. \end{equation} \tag{ A.1.7 }$$

Finally, according to the below equation, calculate the plane $S_2 ^ \bot$, which passes the point ${\bf X}_M^{c_0 }$, and with ${\bf n}_{_{\left({sight - l_r } \right) \bot } }^{c_0 }$ as its normal vector:

$$\begin{equation} \left({{\bf n}_{_{({sight - l_r }) \bot } }^{c_0 } } \right)^T \left({{\bf X} - {\bf X}_M^{c_0 } } \right) = 0. \end{equation} \tag{ A.1.8 }$$

The visible surface of the rotating body is directly determined by the dihedral angle θ between the boundary plane of the Sun rays and that of the line of sight. The angle θ can be calculated through the angle between the two normal vectors ${\bf n}_{_{\left({sun - l_r } \right) \bot } }^{c_0 }$ and ${\bf n}_{_{\left({sight - l_r } \right) \bot } }^{c_0 }$ with the following relationship:

$$\begin{equation} \cos ({\pi - \theta }) = \frac{{{\bf n}_{_{({sun - l_r }) \bot } }^{c_0 } {\bf n}_{_{({sight - l_r }) \bot } }^{c_0 } }}{{\left\| {{\bf n}_{_{({sun - l_r }) \bot } }^{c_0 } } \right\| \cdot \left\| {{\bf n}_{_{({sight - l_r }) \bot } }^{c_0 } } \right\|}}. \end{equation} \tag{ A.1.9 }$$

For convenience, let n₁ represent ${\bf n}_{_{\left({sun - l_r } \right) \bot } }^{c_0 }$, and n₂ represent ${\bf n}_{_{\left({sight - l_r } \right) \bot } }^{c_0 }$. Thus, we yield the following equations:

$$\begin{equation} \theta = \pi - \cos ^{ - 1} \left({{\bf n}_1 {\bf n}_2 } \right) \end{equation} \tag{ A.1.10 }$$

$$\begin{equation} {\bf n}_1 = \left({{\bf l}_{\rm sun}^{c_0 } \times {\bf l}_{4179}^{c_0 } } \right) \times {\bf l}_{4179}^{c_0 } \end{equation} \tag{ A.1.11 }$$

$$\begin{equation} {\bf n}_2 = \left({{\bf l}_{\rm sight}^{c_0 } \times {\bf l}_{4179}^{c_0 } } \right) \times {\bf l}_{4179}^{c_0 }. \end{equation} \tag{ A.1.12 }$$

The part within θ could have been illuminated by the Sun and shot by the camera at the same time, which makes it visible in the optical images. Other parts are regarded as shadowed or regions that are not detected by the spacecraft camera.

Supposing P is a point on Toutatis's surface, ${\bf X}_P^{c_0 }$ is the spatial coordinate of P in the original camera coordinate system, m_P[3 × 1] is the homogeneous image coordinate of P in Img_m, and n_P[3 × 1] the homogeneous image coordinate of P in Img_n. There exists the following relationship between the image coordinates and ${\bf X}_P^{c_0 }$:

$$\begin{equation} {\bf X}_P^{c_0 } = P_{\bf m} ^{\bf + } {\bf m}_P{\left[ {3 \times 1} \right]} \end{equation} \tag{ A.2.1 }$$

$$\begin{equation} {\bf X}_P^{c_0 } = P_{\bf n} ^{\bf + } {\bf n}_P{\left[ {3 \times 1} \right]}, \end{equation} \tag{ A.2.2 }$$

where P_m is the camera matrices of Img_m, P_n is the camera matrices of Img_n,$P_m^ +$ is the pseudo-inverse of P_m, and $P_n ^ +$ is the pseudo-inverse of P_n.

If the camera position of Img_m is taken as the origin and the relative world coordinate system is built based on the reference of the optical axis direction at this moment, the camera matrices of P_m and P_n are expressed as follows, respectively:

$$\begin{equation} \left\{ \begin{array}{@{}l} P_m = K[ {I_{[ {3 \times 3} ]} |{\bf 0}} ] \\ \ms P_n = K\left[ {R_n |{\bf a}_n } \right] \\ \end{array} \right., \end{equation} \tag{ A.2.3 }$$

whereR_n is the camera rotation matrix of Img_n, a_n is the translation vector at Img_n's imaging position, and K is the inner orientation matrix of the camera, which is expressed as

$$\begin{equation} K = \left[ {\begin{array}{*{20}c} {f_u } & 0 & {u_0 } \\ \ms 0 & {f_v } & {v_0 } \\ \ms 0 & 0 & 1 \\ \end{array}} \right]. \end{equation} \tag{ A.2.4 }$$

Use the homogeneous vectors [u_P(t), v_P(t), 1]^T to formulate the image coordinates m_P(t) of the motion point P that is in line with the relation model of movement and imaging herein. Use the homogeneous vectors [x_P(t), y_P(t), z_P(t), 1]^T to formulate the coordinate ${\bf X}_P^{c_0 } \left(t \right)$ in the original camera coordinate system. Based on the multiple-view geometrical relationship, there is a mapping relation between m_P(t) and ${\bf X}_P^{c_0 } \left(t \right)$ as follows:

$$\begin{equation} {\bf m}_P \left(t \right) = K\left({I,0} \right)\left[ {\begin{array}{*{20}l} R & { - R\tilde C} \\ {0^T } & 1 \\ \end{array}} \right]{\bf X}_P^{c_0 } \left(t \right). \end{equation} \tag{ A.3.1 }$$

According to the attitude constraint, whenever t > 0, the above relation can be simplified as follows:

$$\begin{equation} {\bf m}_P \left(t \right) = K(I, - {\bf \tilde C}){\bf X}_P^{c_0 } \left(t \right) \end{equation} \tag{ A.3.2 }$$

$$\begin{equation} {\bf \tilde C} = \left[ {V_x \left(t \right)t,V_y \left(t \right)t,V_z \left(t \right)t} \right]^T, \end{equation} \tag{ A.3.3 }$$

where V_x(t)，V_y(t)，V_z(t), respectively, refer to the velocity components of the motion point P in the original camera coordinate system in the direction of x, y, and z.

According to the scale constraint, ${\bf \tilde C}$ in Equation (A.3.3) can be expressed as

$$\begin{equation} {\bf \tilde C} = Vt{\bf T}_{{\rm norm}}, \end{equation} \tag{ A.3.4 }$$

where T_norm is the normalized velocity vector of the motion point P in the original camera coordinate system, namely, the normalized coordinate of the Chang'e-2 spacecraft in the original camera coordinate system.

Introduce the essential matrix E according to the multiple view geometry, and establish the following relationship between T_norm and E.

$$\begin{equation} {\bf T}_{{\rm norm}} = \left[ {T_X,T_Y,T_Z } \right]^T \left/\sqrt {T_X ^2 + T_Y ^2 + T_Z ^2 }\right. \end{equation} \tag{ A.3.5 }$$

$$\begin{equation} \left[ {\bf T} \right]_ \times {\rm = }\left[ {\begin{array}{*{20}c} 0 & { - T_Z } & {T_Y } \\ \ms {T_Z } & 0 & { - T_X } \\ \ms { - T_Y } & {T_X } & 0 \\ \ms \end{array}} \right] = E.\end{equation} \tag{ A.3.6 }$$

Based on the essential matrix's definition, combined with the epipole constraint, the formulation of E can be expressed as

$$\begin{equation} E = K{^\prime}{^T} FK = K{^\prime}{^T} \left[ {\bf e} \right]_ \times K,\end{equation} \tag{ A.3.7 }$$

where K and K' are camera inner orientation matrixes. As the same camera was used during this flyby, the formulation of K and K' becomes

$$\begin{equation} K = K^{\prime}.\end{equation} \tag{ A.3.8 }$$

Substitute Equation (A.3.8) into Equation (A.3.7); E can be derived thereby as

$$\begin{equation} E = \left[ {\begin{array}{*{20}ccc} 0 & { - f_v f_u } & {e_v f_u - v_0 f_u } \\ \ms {f_u f_v } & 0 & {u_0 f_v - e_u f_v } \\ \ms { - f_u \left({e_v - v_0 } \right)} & {f_v \left({e_u - u_0 } \right)} & 0 \\ \end{array}} \right].\end{equation} \tag{ A.3.9 }$$

Then, the vector [T_X, T_Y, T_Z]^T is yielded, and each component is as follows:

$$\begin{equation} \left[ {\begin{array}{*{20}c} {T_X } \\ \ms {T_Y } \\ \ms {T_Z } \\ \end{array}} \right] = \left[ \begin{array}{l} f_v \left({e_u - u_0 } \right) \\ \ms f_u \left({e_v - v_0 } \right) \\ \ms f_u f_v \\ \end{array} \right].\end{equation} \tag{ A.3.10 }$$

Subsequently, according to Equations (A.3.4) and (A.3.5), ${\bf \tilde C}$ is derived:

$$\begin{equation} {\bf \tilde C} = \frac{{Vt}}{{K}}\left[ \begin{array}{l} f_v \left({e_u - u_0 } \right) \\ \ms f_u \left({e_v - v_0 } \right) \\ \ms f_u f_v \\ \end{array} \right], \end{equation} \tag{ A.3.11 }$$

where ${K} = \sqrt {\left({f_v \left({e_u - u_0 } \right)} \right)^2 + \left({f_u \left({e_v - v_0 } \right)} \right)^2 + \left({f_u f_v } \right)^2 }.$

Substituting Equation (A.3.11) into Equation (A.3.2), the conclusion can be derived:

$$\begin{equation} \left[ {\begin{array}{*{20}c} {u_P \left(t \right)} \\ \ms {v_P \left(t \right)} \\ \ms \end{array}} \right] = \left[ \begin{array}{l} \dsty\frac{{f_u x_P \left(t \right) + u_0 z_P \left(t \right) - \frac{{Vt}}{{K}}f_u f_v e_u }}{{z_P \left(t \right) - \frac{{Vt}}{{K}}f_u f_v }} \\ \ms \ms \dsty\frac{{f_v y_P \left(t \right) + v_0 z_P \left(t \right) - \frac{{Vt}}{{K}}f_u f_v e_v }}{{z_P \left(t \right) - \frac{{Vt}}{{K}}f_u f_v }} \\ \end{array} \right]. \end{equation} \tag{ A.3.12 }$$