Imagery Synthesis for Drone Celestial Navigation Simulation

Abstract

Simulation plays a critical role in the development of UAV navigation systems. In the context of celestial navigation, the ability to simulate celestial imagery is particularly important, due to the logistical and legal constraints of conducting UAV flight trials after dusk. We present a method for simulating night-sky star field imagery captured from a rigidly mounted ‘strapdown’ UAV camera system, with reference to a single static reference image captured on the ground. Using fast attitude updates and spherical linear interpolation, images are superimposed to produce a finite-exposure image that accurately captures motion blur due to aircraft actuation and aerodynamic turbulence. The simulation images are validated against a real data set, showing similarity in both star trail path and magnitude. The outcomes of this work provide a simulation test environment for the development of celestial navigation algorithms.

1. Introduction

Celestial navigation in uncrewed aerial vehicles (UAV) has been a topic of interest for over half a century (see, for example, [1]). The significance of this mode of navigation has been overshadowed, however, by the ubiquity of global navigation satellite systems and the integration of compact micro-electromechanical attitude sensors into aviation platforms. Nonetheless, celestial navigation has unique advantages due to its independence from critical infrastructure and robustness to external interference. We see recent works, such as [2,3] integrating celestial imaging into their navigation solutions. Modern UAVs must typically conform to size, weight and power constraints and, to this end, benefit from a strapdown celestial implementation, as opposed to an actively stabilized alternative. In a strapdown configuration, the imaging sensor has no control authority over the vehicle, and therefore requires a larger field of view, and longer exposure intervals, to track stars during motion. We propose here a method for simulating the imagery captured from such a strapdown celestial system.

Celestial imagery is commonly used in spacecraft to obtain a highly accurate attitude reference. This technique is less commonly used, however, in low altitude aircraft navigation. Aircraft are subjected to many sources of noise that spacecraft are not, such as light pollution, cloud cover, atmospheric diffraction, aerodynamic turbulence, engine vibration and control/actuation, which all impact the signal strength of a celestial image obtained from an aircraft. These effects are exacerbated by the need for long-exposure imagery when operating at low-altitude (less than 500 m). The standard approach to this problem is to use a stabilized viewing platform with a telescopic lens [1], which limits the aforementioned attenuation. Such an approach is costly and adds significant weight to an airframe. The design of this simulation has arisen from the desire for a low cost, low weight, “strapdown” [4] celestial navigation solution.

As with all avionic navigation solutions, simulation plays an important role in the system design and precedes the implementation. The intent of this work is to provide a means of simulating imagery from a camera with a wide-angle lens, rigidly mounted to a fixed wing airframe with no active stabilization. Preliminary testing indicated that despite the increased motion blur, longer exposure images are consistently able to capture stars at lower levels of illumination. Consequently, the need arises for a simulation that can replicate the effects of motion blur due to these longer exposure images. In addition, there is benefit to tuning the simulation based on a reference image captured from the ground. This provides a quick solution for users to encapsulate their camera sensor and lens characteristics in the simulation environment, without the logistical constraints of night flying.

The use of star field simulation is most commonly seen in star identification research and development [5]. In this field, simulation is used to obtain baseline performance metrics for newly designed algorithms. An example of star field simulation for this purpose can be seen in [6], in which, rather than rendering stars, their position and magnitude are generated directly, with the addition of Gaussian noise. Other works tend to follow a similar design, seen in [7,8]. These simulations are intended to replicate imagery captured from spacecraft, which is not typically affected by rotational motion, nor atmospheric attenuation. Recent work considers the effects of star smearing [9] and the effect this has on the observability of stars. This work leverages from the simulation concepts presented in [10], later followed by [11]. These studies assume the angular velocity of the camera to be constant; however, we can see in Section 3 that for aerial vehicles, this assumption is invalid. Advancements were made in [12], highlighting the importance of modeling in star sensor design and calibration. In each of these cases, testing was conducted using a turntable, and as indicated in [10], this approach is not capable of running in real-time due to the large number of integral calculations involved.

This work offers two significant contributions to this field of research. We present a simple and effective framework for the real-time simulation of long-exposure images from non-stabilized UAV-mounted hardware using spherical linear interpolation, and a method for calibrating the simulation based on a ground reference image. Concepts from this work may be extended to aid in simulation design for spacecraft in highly dynamic situations.

2. Simulation Architecture

The position of stars and planets in the sky are represented in the International Celestial Reference Frame (ICRF). The location of a celestial body is expressed in terms of right ascension, $\alpha$, and declination, $\delta$, as seen in Figure 1. Stars are assumed to be infinitely far away, lying on the celestial sphere. Consequently, translational motion has no effect on the apparent position of the stars.

The aircraft position is represented in the Earth-centred, Earth-fixed WGS84 system. This is the standard reference frame for GPS positional data. Aircraft attitude is represented in the local North-East-Down (NED) coordinate system. The camera is assumed to be mounted to the aircraft, sharing a location with the vehicle and differing in orientation by a single rotational transformation. The camera coordinate system is oriented with z positive in the optical axis, x positive towards the right of the image plane and y positive towards the bottom of the image plane.

A star catalogue must be used as a reference for the location of stars. While there are many star catalogues available, we selected the Yale Bright Star Catalogue (BSC) due to its minimal size. The BSC contains records of stars down to magnitude 6.5, totalling 9110 stars. This magnitude threshold is sufficient for most aircraft camera systems (including stabilized systems) [1]. For ease of implementation, the ASCII-format catalogue was converted into an SQLite database. Indices were created for magnitude, right ascension and declination, including additional composite indices to allow fast querying of the database.

This simulation initially corrects for celestial phenomena which perturb the right ascension and declination in the star catalogue, before entering a simulation loop. The simulation loop performs the following steps:

1.

Update position, time and attitude of simulation;

2.

Calculate star homogeneous coordinates in camera frame of reference;

3.

Project each star onto the image plane;

4.

Render the star.

2.1. Initial Corrections

On initialization, adjustments are made for the right ascension and declination of stars due to annual proper motion (the apparent motion of stars), precession (changes in the Earth’s rotational axis over time), nutation (axial changes due to the Moon’s gravitational pull) and aberration (due to the velocity of the Earth’s orbit).

Proper motion is provided in star catalogues and is compensated by computing the change in right ascension and declination since the given epoch, such that:

$$\begin{array}{c}\widehat{\alpha }={\alpha }_0+\dot{\alpha }T\\ \widehat{\delta }={\delta }_0+\dot{\delta }T\end{array}$$

(1)

where ${\alpha }_0$ and ${\delta }_0$ are the right ascension and declination at the epoch. $\dot{\alpha }$ and $\dot{\delta }$ are the annual proper motion, typically expressed in arcseconds per year, and T is the time (in years) since the epoch (J2000 in our case). Subsequently, right ascension and declination are corrected due to precession, following the method outlined in [13]. That is, we use the polynomial approximation for the precession angles $\zeta$, z and $\gamma$:

$$\begin{array}{c}\zeta =2306.2181t+0.30188t^2+0.017998t^3\\ z=2306.2181t+1.09468t^2+0.018203t^3\\ \gamma =2004.3109t-0.42665t^2-0.041833t^3\end{array}$$

(2)

where t is the number of centuries since the J2000 epoch. Right ascension $\alpha$ and declination $\delta$ are then found given:

$$\begin{array}{c}tan(\alpha -z)=\frac{sin(\widehat{\alpha }+\zeta )}{cos\gamma cos(\widehat{\alpha }+\zeta )-sin\widehat{\gamma }tan\widehat{\delta }}\\ sin\left(\delta \right)=sin\gamma cos\widehat{\delta }cos(\widehat{\alpha }+\zeta )+cos\gamma sin\widehat{\delta }\end{array}$$

(3)

Following [13], we correct for the effects of nutation. Nutation is comprised of nutation in longitude, $\Delta {\lambda }_n$, and nutation in obliquity, $\Delta ϵ$. These quantities can be approximated to within 0.5 arcseconds by the following equations (expressed in arcseconds):

$$\begin{array}{cc}\hfill \Delta {\lambda }_n& =-17.20sin\mathsf{\Omega }+1.32sin2L-0.23sin2L^{\prime }+0.21sin2\mathsf{\Omega }\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill \Delta ϵ& =9.2cos\mathsf{\Omega }+0.57cos2L+0.10cos2L^{\prime }-0.09cos2\mathsf{\Omega }\hfill \end{array}$$

(5)

where $\mathsf{\Omega }$, the longitude of the ascending node of the Moon’s mean orbit on the ecliptic, expressed in degrees, is approximated as:

$$\mathsf{\Omega }=125.04452-1934.136261t$$

(6)

where t is expressed in Julian centuries, as above. L and $L^{\prime }$, the mean longitudes of the Sun and Moon, respectively, expressed in degrees, are given by:

$$\begin{array}{c}L=280.4665+36000.7698t\\ L^{\prime }=218.3165+481267.8813t\end{array}$$

(7)

The mean obliquity of the ecliptic can be found given the following equation, expressed in degrees:

$${ϵ}_0=23.439291-13.004166t\times {10}^{-3}-0.1639t^2\times {10}^{-6}+0.50361t^3\times {10}^{-6}$$

(8)

Subsequently, the true obliquity of the ecliptic $ϵ$ is given by:

$$ϵ={ϵ}_0+\Delta ϵ$$

(9)

The resulting corrections due to nutation, $\Delta {\alpha }_n,\Delta {\delta }_n$ for the star’s right ascension and declination, respectively, are given by:

$$\begin{array}{cc}\hfill \Delta {\alpha }_n=& (cosϵ+sinϵsin\alpha tan\delta )\Delta {\lambda }_n-(cos\alpha tan\delta )\Delta ϵ\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill \Delta {\delta }_n=& (sinϵcos\alpha )\Delta {\lambda }_n+(sin\alpha )\Delta ϵ\hfill \end{array}$$

(11)

We also consider the effects of aberration, as presented in [13]. Given the constant of aberration, $\kappa =20.29552$ arcseconds, the true longitude of the Sun ${\lambda }_s$, eccentricity of the Earth’s orbit e, and the longitude of the perihelion, $\rho$, we can compute the corrections $\Delta {\alpha }_a,\Delta {\delta }_a$ for the star’s right ascension and declination with the following equations:

$$\begin{array}{cc}\hfill \Delta {\alpha }_a=-\kappa & \left(\frac{cos\alpha cos{\lambda }_{s}cosϵ+sin\alpha sin{\lambda }_s}{cos\delta }\right)+e\kappa \left(\frac{cos\alpha cos\rho cosϵ+sin\alpha sin\rho }{cos\delta }\right)\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill \Delta {\delta }_a=-\kappa & \left. ['cos{\lambda }_{s}cosϵ(tanϵcos\delta -sin\alpha sin\delta )+cos\alpha sin\delta sin{\lambda }_s\right]\hfill \\ \hfill & +e\kappa \left. ['cos\rho cosϵ(tanϵcos\delta -sin\alpha sin\delta )+cos\alpha sin\delta sin\rho \right]\hfill \end{array}$$

(13)

where:

$$\begin{array}{cc}\hfill e& =16.708634\times {10}^{-3}-42.037t\times {10}^{-6}-0.1267t^2\times {10}^{-6}\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill \rho & =102.9375+1.71946t+0.46t^2\times {10}^{-3}\hfill \end{array}$$

(15)

and the true longitude of the Sun is calculated as:

$${\lambda }_s=L_0+C$$

(16)

where the mean longitude of the Sun, $L_0$, is given by:

$$L_0=280.46646+36000.76983t+0.3032t^2\times {10}^{-3}$$

(17)

and the Sun’s equation of the center, C, is given by:

$$\begin{array}{cc}C=\left[\right(191& 4.602-4.817t-0.014t^2)sinM\hfill \\ \hfill +(& 19.993-0.101t)sin2M\hfill \\ \hfill & +0.289sin3M]\times {10}^{-3}\hfill \end{array}$$

(18)

and the mean anomaly of the Sun, M, is given by:

$$M=357.52911+35999.05029t-0.1537t^2\times {10}^{-3}$$

(19)

The correction terms, $\Delta {\alpha }_a$, $\Delta {\delta }_a$, $\Delta {\alpha }_n$ and $\Delta {\delta }_n$ are added to $\alpha$ and $\delta$, yielding the corrected celestial coordinates of the star.

2.2. Updating Position and Time

Conversion from terrestrial time (as measured with the Gregorian calendar) to celestial time (typically expressed in Julian days) is a necessary step in simulating celestial bodies. The BSC represents the location of stars with respect to the J2000 epoch. The conversion from Gregorian date to Julian day is detailed in [13]. The basic steps are shown in Algorithm 1.

Algorithm 1 Conversion from Gregorian date to Julian day

let Y be the current year                          ▹ Integer, Gregorian let M be the current month                     ▹ Integer [1..12], Gregorian let D be the current day                         ▹ Decimal, Gregorian if$M\le 2$then     $Y=Y-1$     $M=M+12$ end if let A = $int\left(\frac{Y}{100}\right)$ let B =$2-A+int\left(\frac{A}{4}\right)$ J = $int\left(365.25\right(Y+4716\left)\right)+int\left(30.6001\right(M+1\left)\right)+D+B-1524.5$    ▹ Julian Days

This conversion expresses the current time with respect to the number of Julian days that have elapsed since 4716 BC. Furthermore, the J2000 epoch is expressed in relation to the Gregorian date 1 January 2000 and can be found by subtracting the Julian day at that time from the current Julian date:

$$J2000=J_D-2451545.0$$

(20)

Sidereal local time, analogous to celestial longitude, can be calculated given WGS84 longitude and the current time of day expressed in decimal hours. Given longitude $\lambda$, Julian date $J_D$ expressed in the J2000 epoch, and decimal hours $H_D$, the local sidereal time $LST$ is given by [13]:

$$LST=100.46+\left(0.985647J_D\right)+\lambda +15H_D$$

(21)

$LST$ is typically limited to the range $[0,360]$ degrees. The hour-angle $\omega$ of a celestial body can then be simply calculated from $LST$ and right ascension $\alpha$ as:

$$\omega =LST-\alpha$$

(22)

The hour-angle of a body can be used to compute its azimuth and elevation, expressed in local NED coordinates. The local elevation $\theta$ of a celestial body given hour angle $\omega$, declination $\delta$ and latitude $\varphi$ is given by:

$$\theta =\mathrm{asin}(sin\delta sin\varphi +cos\delta cos\varphi cos\omega )$$

(23)

Subsequently, using elevation from Equation (23), the local azimuth $\psi$ is given by the following equation:

$$\psi =\mathrm{atan}2(sin\omega ,cos\omega sin\varphi -tan\delta cos\varphi )+\pi$$

(24)

where the addition of $\pi$ radians converts azimuth to a representation that is positive East of North [13].

The observed elevation of a celestial body is altered due to the effects of atmospheric refraction. Consequently, objects in the sky appear at a greater elevation than they would without the atmospheric effects. This effect is exaggerated at lower elevations (closer to the horizon), which leads to an angular displacement of up to 0.5${}^{\circ }$. For cameras in the visible light spectrum, the refractive distance R (expressed in arcminutes) can be approximated to within 4 arcseconds [14], given the following formula:

$$R=\frac{1.02}{tan\left(\theta +\frac{10.3}{\theta +5.11}\right)}$$

(25)

If a greater level of accuracy is required, alternative methods such as that seen in [15] can be used. The apparent elevation ${\theta }^{\prime }$ is then given by:

$${\theta }^{\prime }=\theta +R$$

(26)

The inverse of this formula, Equation (27) [16], allows for the correction of observations:

$$R=\frac{1}{tan\left({\theta }^{\prime }+\frac{7.31}{{\theta }^{\prime }+4.4}\right)}$$

(27)

Formulas (25) and (27) assume an atmospheric pressure of 1010 millibars, and an air temperature of 10 °C. According to [13], an approximate scale factor may be applied given pressure P at the Earth’s surface, and air temperature T °C, given by the following formula:

$$\frac{P}{1010}\times \frac{283}{273+T}$$

(28)

Finally, we map the celestial sphere onto a unit sphere for the purposes of image plane projection. Given the azimuth and elevation of a star, the corresponding unit vector in local NED coordinates is given by:

$$\begin{array}{cc}\hfill x& =cos\psi cos{\theta }^{\prime }\hfill \\ \hfill y& =sin\psi cos{\theta }^{\prime }\hfill \\ \hfill z& =-sin{\theta }^{\prime }\hfill \end{array}$$

(29)

For resource constrained systems, the rate at which these unit vectors are updated should be chosen relative to the precision required by the simulation. For reference, a geostationary camera with an update rate of 1 Hz will be accurate to within ±15 arcseconds (4.17 $\times {10}^{-3}$ degrees). For an aircraft at a latitude of $\pm 45$${}^{\circ }$, travelling East/West at mach 1, an update rate of 1Hz will be accurate to within $\pm 30$ arcseconds ( $8.{37}^{\circ }$$\times {10}^{-3}$ ).

2.3. Updating Attitude

Representing celestial bodies in the local NED frame of reference simplifies the transformation from aircraft attitude to camera attitude. Aircraft attitude should come either directly from the onboard attitude reference or from a simulator. The work presented here uses the open source Ardupilot toolchain. Attitude log data collected from real flights were used for this research, so as to correlate real images with their simulated counterparts. The roll $p_a$, pitch $q_a$ and yaw $r_a$ Euler angles of the aircraft in local NED coordinates were logged at 30 Hz. These Euler angles can be represented as a rotation matrix, through a yaw–pitch–roll rotation sequence. The rotation matrix, ${\mathit{C}}_{\mathit{a}/\mathit{l}}$, transforms objects in the local NED frame to the aircraft body frame, where ${\mathit{C}}_{\mathit{a}/\mathit{l}}$ is given by:

$$\left. [\begin{array}{ccc}\mathrm{c}\left(q_a\right)\mathrm{c}\left(r_a\right)& \mathrm{c}\left(q_a\right)\mathrm{s}\left(r_a\right)& -\mathrm{s}\left(q_a\right)\\ -\mathrm{c}\left(p_a\right)\mathrm{s}\left(r_a\right)+\mathrm{s}\left(p_a\right)\mathrm{s}\left(q_a\right)\mathrm{c}\left(r_a\right)& \mathrm{c}\left(p_a\right)\mathrm{c}\left(r_a\right)+\mathrm{s}\left(p_a\right)\mathrm{s}\left(q_a\right)\mathrm{s}\left(r_a\right)& \mathrm{s}\left(p_a\right)\mathrm{c}\left(q_a\right)\\ \mathrm{s}\left(p_a\right)\mathrm{s}\left(r_a\right)+\mathrm{c}\left(p_a\right)\mathrm{s}\left(q_a\right)\mathrm{c}\left(r_a\right)& -\mathrm{s}\left(p_a\right)\mathrm{c}\left(r_a\right)+\mathrm{c}\left(p_a\right)\mathrm{s}\left(q_a\right)\mathrm{s}\left(r_a\right)& \mathrm{c}\left(p_a\right)\mathrm{c}\left(q_a\right)\end{array}\right]$$

(30)

where c $\left(x\right)$ and s $\left(x\right)$ represent $cosx$ and $sinx$, respectively. The camera is mounted to the aircraft, with roll $p_c$, pitch $q_c$ and yaw $r_c$ expressed in the aircraft body frame. The z axis of the camera is parallel to the optical axis, positive in the direction of the image plane. The y and x axes are orthogonal and directed towards the bottom and the right of the image plane, respectively. The rotation matrix, ${\mathit{C}}_{\mathit{c}/\mathit{a}}$, transforms objects in the aircraft body frame to the camera frame, where ${\mathit{C}}_{\mathit{c}/\mathit{a}}$ is given by:

$$\left. [\begin{array}{ccc}\mathrm{c}\left(q_c\right)\mathrm{c}\left(r_c\right)& \mathrm{c}\left(q_c\right)\mathrm{s}\left(r_c\right)& -\mathrm{s}\left(q_c\right)\\ -\mathrm{c}\left(p_c\right)\mathrm{s}\left(r_c\right)+\mathrm{s}\left(p_c\right)\mathrm{s}\left(q_c\right)\mathrm{c}\left(r_c\right)& \mathrm{c}\left(p_c\right)\mathrm{c}\left(r_c\right)+\mathrm{s}\left(p_c\right)\mathrm{s}\left(q_c\right)\mathrm{s}\left(r_c\right)& \mathrm{s}\left(p_c\right)\mathrm{c}\left(q_c\right)\\ \mathrm{s}\left(p_c\right)\mathrm{s}\left(r_c\right)+\mathrm{c}\left(p_c\right)\mathrm{s}\left(q_c\right)\mathrm{c}\left(r_c\right)& -\mathrm{s}\left(p_c\right)\mathrm{c}\left(r_c\right)+\mathrm{c}\left(p_c\right)\mathrm{s}\left(q_c\right)\mathrm{s}\left(r_c\right)& \mathrm{c}\left(p_c\right)\mathrm{c}\left(q_c\right)\end{array}\right]$$

(31)

Given rotation matrices ${\mathit{C}}_{\mathit{a}/\mathit{l}}$ and ${\mathit{C}}_{\mathit{c}/\mathit{a}}$, the transformation of a unit vector, $\mathit{u}$, from the local NED frame to the camera frame, can be computed as:

$${\mathit{u}}^{\mathit{c}}={\mathit{C}}_{\mathit{c}/\mathit{a}}{\mathit{C}}_{\mathit{a}/\mathit{l}}\mathit{u}$$

(32)

Equation (32) allows the unit vectors of the celestial bodies computed in Section 2.2 to be represented in the camera frame of reference.

2.4. Projection

We assume that the camera intrinsic matrix, $\mathit{K}$, is known. The intrinsic properties of a camera can be found through a calibration method such as that described in [17], yielding matrix:

$$\mathit{K}=\left[\begin{array}{ccc}{f}_x& s& x_0\\ 0& f_y& y_0\\ 0& 0& 1\end{array}\right]$$

(33)

where $f_x$ and $f_y$ are the x and y focal lengths in pixel units, $x_0$ and $y_0$ are the x and y pixel locations of the principal point, respectively, and s describes the sensor skewness (typically 0 for digital sensors).

We first convert each celestial body unit vector into homogeneous coordinates. For components x, y and z of unit vector ${\mathit{u}}^{\mathit{c}}$, the homogeneous point $\mathit{P}$ in the camera frame of reference is calculated as:

$$\mathit{P}=\left[\begin{array}{c}\frac{x}{z}\\ \frac{y}{z}\\ 1\end{array}\right]$$

(34)

Objects whose depth value, z, is less than 0 are ignored due to being positioned behind the camera object. Finally, given no translational component to our camera, the pixel location, $\mathit{x}$, of the object is found as:

$$\mathit{x}=\mathit{K}\mathit{P}$$

(35)

Lens distortion may also be included in the model. Lens distortion models are typically non-linear, and expressed as a function of displacement from the principal point. Various models exist for lens distortion, and should be chosen according to the level of precision required [18]. If using a lens distortion model, this should be applied after the star is rendered, so as to capture the resultant eccentricities from the distortion. We do not model lens distortion in our simulation; instead, we rectify all images prior to analysis, such that any residual distortion is negligible.

2.5. Calibration

The apparent pixel intensity of a star is determined by the apparent star magnitude, atmospheric conditions, lens properties and sensor properties. We present here a method which precludes the need for detailed modeling, through a single-image calibration process. We use a reference image captured from a stationary aircraft on the ground and fit an exponential curve to define the relationship between star magnitude and pixel intensity.

The relationship between observed brightness and apparent star magnitude is given by the equation:

$$m_x-m_r=-2.5{log}_{10}\left(\frac{B_x}{B_r}\right)$$

(36)

where $m_x$ is the observed star magnitude, $m_r$ is the reference star magnitude, $B_x$ is the observed star brightness and $B_r$ is the reference star brightness. This magnitude scale is designed such that a magnitude difference of $-5$ correlates with a brightness factor of 100. That is to say, a magnitude 1 star is 100 times brighter than a magnitude 6 star. For this work, we take the brightness of a star to be its maximum pixel value. Images are converted from the blue–green–red (BGR) colour space to the hue–saturation–value (HSV) colour space, and the value channel is used as a greyscale image.

By rearranging Equation (36) we can compute the brightness of an observed star, given that we have a reference brightness $B_r$, reference magnitude $m_r$ and observed magnitude $m_x$:

$$B_x=B_r{10}^{\frac{m_x-m_r}{2.5}}$$

(37)

The choice of reference star is an important factor, as the magnitude of stars are typically considered to be unreliable [19]. Many factors can cause the apparent magnitude of a star to differ from the catalogue, including spectral attenuation caused by the atmosphere, camera characteristics, atmospheric refraction, as well as the luminescent characteristics of the star itself. We assume that the magnitude error follows a zero-mean Gaussian distribution. Following this assumption, we select the star with magnitude and brightness that minimizes the Frobenius norm of the difference between observed and calculated brightness:

$$min\left({\left[\sum _{i}abs{(B_i^{\prime }-B_i)}^2\right]}^{1/2}\right)$$

(38)

where $B_i^{\prime }$ is a vector containing all calculated star brightnesses, and $B_i$ is a vector containing all measured star brightnesses. The vector $B_i^{\prime }$ is computed for each star in the reference image by choosing $B_r$ and $m_r$ from the reference star, and recomputing $B_x$ for all stars in the image using Equation (37). Stars whose brightness is saturated (i.e., have a maximum pixel value of 255) are excluded from this process. Figure 2 shows the output from this calibration process, conducted on an image captured from a grounded aircraft, using a PiCamHQ with 500 ms exposure interval. We note that this procedure is most effective with a larger number of visible stars, such that the sample better represents the population. In this example, stars were automatically detected and matched to the database using a star identification algorithm. As an alternative to this automated process, one could manually label each star in the reference image.

Stars are effectively point light sources. The camera lens will tend to defocus a star, such that it appears larger than one pixel. The apparent size of a star is affected by magnitude only to the extent that otherwise-undetectable pixels exceed the noise floor. Consequently, we can treat each star as having constant size, and changes to the star’s brightness will yield a larger or smaller effective area. We assume that the pixel intensity of a star follows a two-dimensional Gaussian point spread function, with diagonal covariance matrix whose elements are equal. That is, a star image is circularly symmetric about its center. In practice, this may not be the case, and we note that such methods might be used for calibrating lens and sensor distortion (this is, however, out of scope for this work). Following this assumption, we measure the standard deviation across the x and y axes of each normalized star detection in the reference image and use the median value across all stars as the reference size for rendering in the simulation. We normalize by scaling the peak pixel value to 1. The standard deviation ${\sigma }_B$ is calculated as:

$${\sigma }_B=\frac{1}{n}\sum _n\mathrm{std}\left({\mathbf{B}}_n\odot \frac{1}{\max \left({\mathbf{B}}_n\right)}\right)$$

(39)

where ${\mathbf{B}}_n$ is the histogram of intensities with respect to pixel position, std$\left(x\right)$ is the standard deviation of sample set x, the symbol ⊙ represents element-wise multiplication, and max $\left({\mathbf{B}}_n\right)$ is a scalar value equal to the maximum-valued element of ${\mathbf{B}}_n$.

A graph showing the standard deviation in pixel intensity vs. star magnitude can be seen in Figure 3, demonstrating the approximate uniformity of apparent projected star size across various magnitudes. Figure 4 shows an example of a Gaussian star render, reconstructed from the standard deviation and star brightness.

Simulation noise levels are calibrated from the reference image. Sources such as moonlight and atmospheric light pollution contribute to Gaussian noise observed in an image and consequently reduce the signal-to-noise ratio and observability of stars. The mean and standard deviation is measured in the value channel of the image. We ignore sections of the image in which stars have been detected so as to avoid bias introduced by the stars themselves. Thus, we have the noise function:

$$X\sim N(\mu ,{\sigma }^2)$$

(40)

2.6. Rendering

A long exposure image can be generated by superimposing multiple short-exposure images. Doing so requires a fast attitude update so as to reproduce the motion of the aircraft throughout the exposure period. The required temporal resolution is affected by the aircraft dynamics; however we found that a 5 ms (200 Hz) attitude update results in the contiguous rendering of stars for wide-angled lenses, even during aggressive maneuvers. A 200 Hz attitude measurement is not typically available on low-cost hardware, so we interpolate attitude measurements using the spherical linear interpolation algorithm described in [20] to achieve the desired rate. This method requires orientations to be represented as quaternions. The details are omitted here for brevity; however, most computer graphics libraries support this operation, transforming the direction cosine matrix in Equation (32) to quaternion format. We can retrieve any arbitrary orientation between two quaternions with the following equation:

$$\mathrm{Slerp}(q_1,q_2,u)=q_1{\left(q_1^{-1}{q}_2\right)}^u,\{u\in \mathbb{R}:0\le u\le 1\}$$

(41)

where Slerp is the spherical linear interpolation function; a value of $u=0$ returns $q_1$, a value of $u=1$ returns $q_2$ and intermediate values of u provide interpolation along the shortest path between $q_1$ and $q_2$ on the unit sphere. Interpolation is performed such that attitudes are captured at 5 ms intervals from the beginning of the camera exposure interval, to the end, yielding a total of $n=200\Delta t$ attitude references, where $\Delta t$ is the camera exposure time in seconds.

The long-exposure image is constructed by superimposing n floating point images. For each attitude $a_i$, the pixel location $\mathit{x}$ of each star $S_j$ in the database is found from Equation (35). A discrete Gaussian kernel, $\mathit{G}$ is constructed using the standard deviation found in Equation (39). The kernel is programmatically generated (see [21]) such that the value $G_i$ of element i is given by:

$$G_i=\alpha exp\left(\frac{-{(i-\frac{(k-1)}{2})}^2}{\left(2{\sigma }^2\right)}\right),\{i\in \mathbb{Z}:0\le i<k\}$$

(42)

where the kernel size $k=\lceil 6{\sigma }^2\rceil$ is odd and is selected to contain a minimum of $99.7\%$ of the total star energy, and $\alpha$ is selected such that ${\sum }_i{G}_i=1$. The kernel is subsequently scaled, such that the maximum element at index $i=\frac{k-1}{2}$ is equal to the peak pixel value, $B_x$, calculated from Equation (37) using the magnitude $m_x$ of star $S_j$, as well as the reference magnitude and intensities, $m_r$ and $B_r$ respectively. Finally, the kernel is scaled down by a factor of n. Assuming the photon flux density is constant across the exposure window, a single short-exposure window contains a fraction $\frac{1}{n}$ of the total star energy. This scaling of kernel G generated in Equation (42) is given by:

$$G_i=G_i\left(\frac{B_x}{n\ast G_{max}}\right),\{i\in \mathbb{Z}:0\le i<k\}$$

(43)

Note that $\mathit{G}$ is stored as a floating point array. The short-exposure star is drawn to the image canvas by centering the kernel $\mathit{G}$ on pixel $\mathit{x}$ and rotating 180${}^{\circ }$ about the centre, adding the values of $\mathit{G}$ to the canvas, so as to render the star symmetrically. This process is repeated for each star in the database $S_j$ at the current attitude.

Once each of the short-exposure images have been rendered, the Gaussian noise defined in Equation (40) is added to the canvas. Finally, the image is converted from a single-channel floating point image to an 8-bit single-channel image. The resulting image contains stars rendered with the motion blur caused by camera movement throughout the exposure interval.

3. Results

We compare here the simulation output with images captured during a flight to evaluate the performance of the simulator. Attitude logs from the flight test were recovered, and used to generate these simulation images. We followed the procedure outlined in Section 2 for image generation. We used a Phantom FX-61 airframe (Figure 5), with a PixHawk v2 autopilot, a Raspberry Pi 4 companion computer for image capture and storage and a Raspberry Pi High Quality Camera sensor, mounted with the official 6mm wide-angle lens. The camera was rigidly mounted to the autopilot, so as to mechanically couple sources of vibration and deflection. An approximate transform from aircraft body frame to camera frame is used for this test, at a yaw angle of $-90$${}^{\circ }$, pitch angle of 90${}^{\circ }$ and roll angle of 0${}^{\circ }$ (given a yaw-pitch-roll Euler sequence).

A flight was conducted capturing images every 3 s, with an exposure interval of 500 ms and ISO set to 800. Ground images were captured prior to launching the aircraft. The aircraft climbed to an altitude of 150 m above ground level, and loitered for several minutes before landing. A total of 130 images were captured for analysis. Attitudes throughout each exposure window were recorded as image metadata at a rate of 30 Hz, and retrieved from the flash logs after the flight. A ground image was selected for performing the simulation calibration outlined in Section 2.5. Simulation images were subsequently generated using these calibration parameters.

A visual comparison of simulation output against in-flight images can be seen in Figure 6 and Figure 7. By observation, we can see that the shape and trajectory of the simulated light trail closely matches reality. Figure 6 shows a cluster of stars, captured while the aircraft was yawing at a rate of 8${}^{\circ }$ per s. The difference in star trail direction between stars is due to alignment of the yaw axis, approximately directed toward the centre of the image. Figure 7 pictures the brightest visible stars in the image, highlighting the effectiveness of the spherical linear interpolation, and its ability to reproduce trajectories with visual magnitude similar to what is observed in reality.

We evaluate the performance of the simulator based on its ability to replicate in-flight star intensities, given a ground calibration image. Figure 8 shows individual stars which were identified across multiple images in-flight. The intensity of the simulated star was plotted against the intensity of the observed star, such that points lying on the line $y=x$ represent a perfect match between reality and simulation. Three stars are shown here; these stars were identified frequently in both reality and simulation. We can see that, in practice, there tends to be error in the pixel intensity. The mean error is near-zero with a value of 3.38 pixels. The mean absolute percentage error in pixel intensity is measured at 47.4% (this is an average of per-star absolute percentage errors), which is similar to the mean absolute percentage error of 51% in the brightness ground calibration, as expected. Using Equation (36), a 47.4% error in star magnitude correlates to an absolute magnitude error of 0.42. This is within comparable range to the noise level simulated in other works, such as [7,8], which artificially add magnitude noise with standard deviation in the range of 0.3–0.9 to their simulation.

Furthermore, we consider the difference between the centroids of stars detected in both simulation and reality. A region of interest (ROI), $\mathit{R}$, is chosen for each star such that $\mathit{R}$ is the smallest grayscale image that contains the star. We compute the weighted centre $D_x,D_y$ for both real and simulated images, expressed with respect to the centre of the ROI:

$$\left[\begin{array}{c}{D}_x\\ D_y\end{array}\right]=\left[\begin{array}{c}\frac{w}{2}\\ \frac{h}{2}\end{array}\right]-\frac{1}{M}\sum _{x,y}{R}_{x,y}\left[\begin{array}{c}x\\ y\end{array}\right]$$

(44)

$$M=\sum _{x,y}{R}_{x,y}$$

(45)

where w and h are the width and height of the region of interest, respectively, x and y are the column and row indices of image R and scalar $R_{x,y}$ is the pixel intensity of $\mathit{R}$ at pixel $(x,y)$.

We compute the L2 norm of the difference between simulation centroid and real centroid to find the distance L:

$$L=\sqrt{{(D_x^s-D_x^r)}^2+{(D_y^s-D_y^r)}^2}$$

(46)

for simulation centroid $D_x^s,D_y^s$ and real centroid $D_x^r,D_y^r$. A histogram containing the computed centroid errors can be seen in Figure 9. For reference, we also compute a baseline estimate, which assumes the centroid is located at the centre of the ROI (analagous to simulating stars as a straight line with uniform intensity). The mean simulation error is measured to be 0.92 pixels, with a median error of 0.68 pixels. By comparison, the mean baseline error is measured at 1.23 pixels, with a median error of 0.93 pixels. This corresponds to a 25.2% reduction in mean centroid error.

4. Discussion

The temporal correlation between the camera and attitude sensor, as well as the attitude sensor’s accuracy and resolution, and individual differences between the database and observed star magnitudes, pose limitations to the accuracy with which simulation images can be generated. The tendency towards a low mean error, however, is an indication that there is little systematic error propagating from the simulation architecture, and that large sample sets will provide a statistically accurate representation of star intensity.

We can see in Figure 8 that despite a low mean error across the sample set, individual stars tend to be subjected to biases in intensity. Future work could make use of multiple ground images to map the intensity of individual stars, so as to reduce bias within individual stars. Furthermore, it is apparent that brighter stars will tend to be detected more frequently than dimmer stars. One may be able to determine an appropriate magnitude threshold for the observability of stars in-flight and bias the calibration towards the brighter stars, which are more likely to be detected.

We note that the ground calibration process is effective only for low-altitude applications. The simulation does not account for changes in atmospheric attenuation due to altitude. Higher-altitude flight will also result in changes to atmospheric refraction. Furthermore, the ground calibration process is subjected to light pollution, which again is a function of altitude. The simulation of higher altitude flight should make use of atmospheric models to account for these disparities between ground and high-altitude observations. The level of image noise due to moonlight is assumed to be constant here; however, in practice there is some degree of variation as the viewing angle changes with respect to the position of the moon. This is most noticeable within a 5${}^{\circ }$ viewing angle [22], but less significant at greater angles.

Validation of this simulation was conducted with a fixed wing aircraft; however, the simulation architecture is applicable for any airframe which is capable of reporting its attitude. This might also be used for simulating motion artefacts from two-axis gimbals. While the Raspberry Pi Camera HQ is fit with a rolling shutter, long exposure images are achieved by a series of shorter exposures, similar to the process followed in this simulation. The intra-frame motion is not captured by this simulation; however, this effect appears to be negligible. If the characteristics of a rolling shutter are known, one could replicate this effect by interpolating at a faster rate and selecting an appropriate attitude given the time at which the shutter exposes the star. It is common, however, for charged-couple device (CCD) cameras to be used for celestial imaging. These cameras utilize a global shutter, which exposes all pixels simultaneously and hence are not subjected to the rolling shutter effect.

5. Conclusions

The intent of this work was to design and validate a simulation architecture to support the development of strapdown celestial navigation solutions in lightweight, low-altitude aircraft. An architecture for replicating the effects of long-exposure imagery was designed and implemented by superimposing multiple short-exposure images from aircraft attitude data. Additionally, a method for augmenting low-rate attitude data was proposed and validated. Simulation calibration was achieved through a single ground reference image, producing results which match reality within reasonable tolerance. The simulation architecture provides a tool for baseline testing star detection and identification algorithms. Future work will extend the capability of this simulator from low-altitude to high altitude through atmospheric modeling.