Keywords

Commercial Camera; Optical Observation; Estimation Algorithms; Genetic Algorithms; Simulated Annealing.

Introduction

Commercial camera utilization for space applications has attracted the attention of vast number of researches. Accurate spacecraft attitude determination requires a good estimation of the camera optical parameters. For star sensors, it is common to estimate focal length, optical distortion, and principal point. These parameters are considered to be critical for star identification algorithms which are utilized to identify certain stars located at the captured image. This process is usually characterized by a huge computational effort which is very sensitive to these parameters. Samaan M (2012) utilizes a commercial off-the-shelf camera as a low cost star tracker. Camera optical parameters and lens distortion are estimated in Zhou F (2015), Dzamba T (2009), discusses the problem of characterization of field curvature and lens astigmatism aberrations. Pal M, Bhat M (2009), solves the problem of spacecraft attitude determination independent of the problem of camera calibration except for the distortion of the camera lens. The star spot location is estimated in Liu HB (2011) based on Kalman filter only. Samaan MA (2001) uses only two methods to estimate the camera focal length. Samaan MA (2003) used the least squares to optimally estimate the focal length.

Spacecraft in-orbit failure could be divided into, mechanical failure, software failure, electrical failure, and unknown [8]. When a spacecraft is totally lost due to any of these reasons, it is necessary to determine whether it is still in-orbit or not. If such a case is encountered, the optical observation of the spacecraft is an optimum choice to assure that the spacecraft is still in-orbit, see [9] for more details. In addition, optical observations could help to increase the accuracy of orbit predictions. Nowadays, commercial cameras could be used instead of complex, heavy, and expensive telescopes for optical observations. Spacecraft orbit estimation based on optical observations captured by a commercial camera is considered to be a challenge that has never been posed before. The nature of the process of spacecraft optical observation is somehow different than that of star identification. Thus, the critical parameters of the two processes are not the same. For, the problem of spacecraft optical observation the camera focal length is considered to be the main parameter of utmost importance.

The contribution of this research is to use approximately costless algorithms of commercial camera calibration for the purpose of spacecraft orbit determination. We should also take into account that the problem of spacecraft orbit estimation, based on optical observations captured by a commercial camera, is considered to be a challenge. Thirteen different algorithms are utilized such as Kalman filter, unscented Kalman filter, derivative free implementation of the extended Kalman filter, genetic algorithms,and simulated annealing are utilized to accurately determine the camera focal length based on actual measurements provided by the camera. The results of all of these algorithms are compared to each other in terms of an error cost function that needs to be minimized.

Process and Measurement Model

The camera focal length is commercially provided by its manufacturer. The utilized camera model in this research is supposed to have a minimum focal length of 4mm. The process model is given by

Where the state vector, X, consists from the camera focal length. Figure.1 shows the observation geometry associated with the calibration process. The calibration process starts by taking an image for at least two reference stars in the sky using free software such as. The star azimuth angle is named, Φc, and the star elevation angle is named, λc [10] expressed in the sensor reference frame of axes. Camera attitude angles are defined with respect to North West Zenith (NWZ) system of axes, shown in Figure. 2. In the case of zero attitude angles, the camera body axes are coinciding with the NWS coordinate system. Thus, the star unit vector expressed in the NWZ reference frame is given by

---

---

---

The first and the second star unit vector expressed in the North West Zenith (NWZ) system of axes is given as

Recall that, Φ1, λ1, Φ2, and λ2 are given by any free software such as Stellarium. Therefore, V1, and V2 are totally known. Similarly, these vectors are represented in the camera body axes (denoted by the subscript c), as

Measurements of the focal length could be calculated from either of the two equations

and Where,
U : Is the measured horizontal distance of a reference object imaged by the camera at the focal plane.
V : Is the vertical distance of a reference object imaged by the camera at the focal plane.
Φc : Object azimuth angle with respect to camera body axes.
λc : Object elevation angle with respect to camera body axes.

Thus, we have two sets of measurements.

The Kalman Filter

Traditionally, the Kalman filter was used to estimate system states and reduce the effect of noise and disturbances. The best estimate is chosen such that, the expected value of the error squares’ sum is minimum. So,

Now, suppose that we have a continuous process model described by

Where u the system is input, and B is the input matrix. Equation (8) could be transformed to its equivalent discrete form with sampling interval ΔT as

Reordering

Therefore

with

And
wk-1 : Is the process noise characterized by.

Qk, is the non-negative definite covariance matrix of the process noise. The measurement process is modeled as

And vk is the measurement noise characterized by.

Rk the non-negative definite covariance matrix of the measurement noise. The time and measurement update equations of the Kalman filter are given as

With zˆk : Estimated measurements.
A k: State transition matrix.
Xˆk : Priori estimate of the state.
Pk : A posteriori covariance of the estimation error.
Xˆk : A posteriori state estimate.
H k : Measurement matrix.

Unscented Kalman filter (USKF)

The unscented Kalman filter algorithm follows the fundamental steps of the Extended Kalman Filter (EKF). The difficulties associated with the traditional EKF are alleviated by the USKF. The filter basic structure given in [11], and [12] is briefly reviewed in this section. The USKF prediction stage starts by forming (2N+1) sigma points as follows

With ξ defined as an adjustment coefficient calculated from

And The filter time update steps for the sigma points are

With the weights calculated from

Thus, the state estimate Xˆk , is given by the relation

Accordingly, the estimated observation corresponding to χj,kis given as

The equations which represent the filter measurement update are given by

Roh K (2007) also describes the basic structure of the USKF with some modifications given below.

With, λ, defined as a scaling parameter that could be computed using the relation

Given that ( α = 1), ξ is a secondary scaling parameter given as

With

The filter time update steps for the sigma points are still the same as equation (27). With the weights calculated from

(β = 2 For Gaussian distribution). Therefore, the state estimate Xˆk , is given by the relation

The estimated observation corresponding to χj,k is given throughout the relation The equations which represent the filter measurement update are now given by

Derivative free Implementation of the Extended Kalman filter (DFEKF)

Derivative free implementation of the extended Kalman filter was firstly developed in [14]. If the number of filter states are given as N, then N vectors, ΔX1 , ΔX2, …. ΔXN, could be formed as follows,

The constant αD is a scaling parameter. The estimated observations at each YI,k could now be calculated from Thus, the filter measurement update stage is given by

Filter Covariance

The process noise covariance, Q k, in its discrete form is related to its continuous form, Q , by the following equation [15]

Similarly, the measurement noise covariance, R k , in its discrete form is related to its continuous form, R(t) , by the following equation [16]

We should also from ref. [16] note that, for a multi-rate Kalman filter (i.e, different measurement rates coming from different sensors), R-1k →0, and hence Kk →0 .

Observability Analysis

The filter state transition matrix based on Floquet theory for a discrete time system is given by [17]

For a stable filter, the A0 matrix Eigen values must all have a magnitude less than unity. The rate of convergence is indicated by the smallness of the Eigen values.

Finding the Global Optimum Solutions and Formulation of the Cost Function

It is well known that there exists several algorithms to find global optimum solutions (i.e. finding the minimum or maximum) of a certain cost function. Genetic algorithms, and simulated annealing could be used efficiently for such purpose. These algorithms will not be reviewed herein. The reader should refer to standard text books such as [18] to find more relevant information. Camera attitude with respect to NWZ coordinate system is represented by three rotations φc , λc , θc around three independent axes, Z, Y, and X respectively. Thus, for the first star, we could write down the following transformation matrices

The same transformations are applied also to the second star. So, if a unit vector pointing to the north, Vn , it follows that

And similarly

Based on the above treatment, the cost function to be minimized could be expressed as

At this point, we have two solutions. The first solution is to search for the minimum of the cost function in the 4 dimensional space corresponding to f, Φc, λc, and θc. The second solution is to search only for the focal length that minimize this cost function and obtain Φc, λc, and θc using the TRIAD method described in [19]. Both solutions are investigated herein. The minimum of the cost function given by equation (77), could be obtained using numerical algorithms such as genetic algorithms (GA), or simulated annealing (SA).

Experimental validation

Estimation algorithms

The performance of the prescribed estimation algorithms is evaluated based on actual measurements. The estimation error is plotted in Figure.1 for the three algorithms which are namely: the Kalman filter (KF), the derivative free Kalman filter (DFEKF), and the unscented Kalman filter (USKF). The error is measured with respect to a simulated camera that has a focal length equal to that was given at the obtained image prosperities written by the camera for each image. Three camera models are used in this research, Benq GH600, Canon PowerShot SX150 IS, and Samsung DV100 Digital camera. Figure.3. Shows the estimated focal length using KF, DFEKF, USKF, in addition to the measured focal length based on λ measurements (LM) and Φ measurements (PM) respectively of the camera model Benq GH600. All of the calibration images have a theoretical focal length of 4mm as indicated by the prosperities of the images.

As clear in these figures the filters have succeeded in filtering the high noise associated with measurements. In addition, the estimator is able to converge quickly despite of large initial estimation error, which is thrown in purpose to prove the high performance of the estimator. The maximum Eigen values of A0 the matrix is 2.9×10-152 which is a very small value that indicates a high rate of convergence. If the measurements were smoothed, the resulting estimated focal length is shown in Figure 4.

---

---

---

---

Iteration based methods

The required parameters to be determined from the image are camera attitude angles and focal length. The proposed algorithm scans a range of focal lengths, and for each foal length solves for the attitude angles using the triad method, and finally computes an error. At this point there are several algorithms are adopted as follows:

Algorithm A: The range of focal lengths is selected and scanned by the simply evaluating the error corresponding to each focal length, and the focal length with minimum error is chosen based on simple gridding technique.

Algorithm B: The error is plotted versus the focal length, and the focal length is selected graphically from the graph, as shown in Figure.5. As clear in this figure the error function is even.

Algorithm C: The focal length is determined based on genetic algorithms.

Algorithm D: The focal length is determined based on simulated annealing.

Algorithm E: The focal length and camera attitude are determined based on simple gridding.

Algorithm F: The focal length and camera attitude are determined based on genetic algorithms.

Algorithm G: The focal length and camera attitude are determined based on simulated annealing.

A comparison among all of the utilized algorithms is given in Table 1.

As shown in Table.1. The best performance (in terms of the selected cost function) is achieved by algorithms G, and F respectively. Algorithms A, B, C, D, E, and F show medium performance. Data smoothing has enabled estimation algorithms to enhance their performance. There are also some important otes those must be mentioned regarding estimation algorithms. The first note is that the cost function that is minimized by these estimation algorithms is given by equation (7) not equation (77). So, their performance is not the best in terms of the formed cost function defined in equation (7).

Therefore, these algorithms obtain lower performance in terms of the cost function defined by equation (77) than the other algorithms. The second note that should be mentioned is that these algorithms converge quickly, nearly in the second time step despite of large initial estimation error because of the strong observability indicated by the small Eigen values of equation (71). Consequently, it is considered to be a good design approach to initialize the process of focal length estimation based on estimation algorithms and then proceed to a fine tuning process by using iteration based algorithms such as algorithm F. On the other hand, estimation algorithms are complex and require lots of mathematical treatment. The simplicity of medium performance algorithms such as Algorithm B, may represent an advantage even though they are characterized by medium performance. The algorithms presented in Table 1 are applied for two camera models which are Canon PowerShot SX150 IS (Table 2), and Samsung DV100 (Table 3) Digital camera respectively. As clear in both tables, algorithms B, and F, usually achieve a very good performance. This result is identical to what is obtained in Table 1.

---

---

---

---

Conclusion

The problem of camera focal length determination has turned out to be of much importance for the purpose of spacecraft orbit observation based on commercial camera instead of large, heavy, complex, and expensive telescopes. Thirteen different algorithms are examined extensively to solve this problem. A comparison among these algorithms showed that estimation algorithms are able to converge despite of large initial estimation error. The performances of these algorithms are enhanced if they are preceded by the process of measurement smoothing. On the other hand estimation algorithms are characterized by complexity and low performance in terms of the cost function defined by equation (77). Some medium performance algorithms such as Algorithm B are characterized by their simplicity compared to estimation algorithms. Algorithms G, and F, exhibit high performance but they are based on iteration. Accordingly, a good design approach could be achieved if estimation algorithms are used first, and then a fine tuning process is performed by using an iteration algorithm such as algorithm F.