AUTONOMOUS ONLINE NAVIGATION AND PARAMETER ESTIMATION IN THE VICINITY OF SMALL CELESTIAL BODIES

Abstract

An exemplary system and method are disclosed for a standalone vision-based system for autonomous online navigation, e.g., around an unknown target small body. The exemplary system (also referred to as AstroSLAM) is predicated on the formulation of the SLAM problem as an incrementally growing factor graph. Using the GTSAM library and the iSAM2 engine, the exemplary system and method combine sensor fusion with the prior orbital motion to provide navigation and parameter estimation that improves the performance over a baseline SLAM implementation. The exemplary system and method can incorporate orbital motion constraints into the factor graph by devising a relative dynamics factor that can link the relative pose of the spacecraft to the problem of predicting trajectories stemming from the motion of the spacecraft in the vicinity of the small body.

Description

BACKGROUND

Precise relative navigation techniques, incorporating increased levels of autonomy, will be a key enabling element of future small celestial body orbiter missions. Firstly, good navigation can inform safe and efficient path planning, control execution, and maneuvering. In near small celestial body missions, achieving fuel-efficiency during non-critical maneuvers and guaranteeing execution of safety-critical maneuvers requires precise knowledge of the relative position and orientation of the spacecraft with respect to the small celestial body. Secondly, precise navigation situates the acquired science data. Indeed, scientists and mission planners design science acquisition phases based on the expected scientific value of instrument data acquired at predetermined times, on specific orbits and with specific spacecraft orientations. Thirdly, precise navigation facilitates the detailed mapping and shape reconstruction of the target small celestial body, since good knowledge of the spacecraft relative position and orientation with respect to the target, as well as a good knowledge of the Sun light direction, are crucial in commonly used shape reconstruction solutions. Finally, good estimates of the spacecraft state enable precise characterization of the target small celestial body's spin state, mass moment values and gravitational model.

In recent years, with ever improving navigation solutions, space missions have successfully performed daring firsts in navigation around small celestial bodies. Orbiter Near-Earth Asteroid Rendezvous (NEAR) Shoemaker's controlled asteroid touchdown, Hayabusa I & II's touchdown and successful sample return, Dawn's orbiting of two celestial bodies in a single mission and the recent Origins, Spectral Interpretation, Resource Identification, Security Regolith Explorer's (OSIRIS-REx's) Touch-and-Go (TAG) operation leveraging Natural Feature Tracking (NFT) relying on high navigation solution accuracy during descent, are only few of the most notable feats accomplished thanks to autonomous navigation.

OSIRIS-REx proximity operations at the near-Earth asteroid Bennu, in particular, pushed the boundaries of what can be accomplished using primarily ground-in-the-loop navigation techniques. Like similar small celestial body missions, OSIRIS-REx's proximity to Bennu, as well as the asteroid's small size and low gravitational attraction relative to perturbing forces, drove the need for frequent and timely navigation updates in order to achieve mission objectives. These updates drove operations complexity and cadence, challenging the flight team, and heavily utilizing Deep Space Network (DSN) assets. OSIRIS-REx proximity operations navigation and TAG also relied on detailed local and global topographic maps constructed from image and LiDAR data, with ground sample distances ranging from 75 cm down to 8 mm. Building these maps required dedicated, months-long observation and data collection campaigns and a substantial amount of effort by the Altimetry Working Group on the ground, as well as multiple iterations with the navigation team. It also required downlinking tens of thousands of images and hundreds of gigabytes of LiDAR data from the spacecraft through the DSN.

It is recognized that the high-risk nature of missions in proximity of small bodies, along with a lack of autonomy in current mission procedures, severely limits the possibilities in mission design. Indeed, ground-segment operators are intimately involved in all in-situ tasks, which ultimately rely on extensive human-in-the-loop verification, as well as ground-based computations for estimation, guidance, and control. In addition, long round-trip light times and severely limited bitrate in communications render ground-in-the-loop processes extremely tedious. In tandem, the incorporation of autonomous capabilities has the potential to improve navigation performance and reduce operational complexity for future missions.

SUMMARY

An exemplary system and method are disclosed for a standalone vision-based system for autonomous online navigation, e.g., around an unknown target small body. The exemplary system (also referred to as AstroSLAM) is predicated on the formulation of the SLAM problem as an incrementally growing factor graph. Using the GTSAM library and the iSAM2 engine, the exemplary system and method combine sensor fusion with the prior orbital motion to provide navigation and parameter estimation that improves the performance over a baseline SLAM implementation. The exemplary system and method can incorporate orbital motion constraints into the factor graph by devising a relative dynamics factor that can link the relative pose of the spacecraft to the problem of predicting trajectories stemming from the motion of the spacecraft in the vicinity of the small body.

In an aspect, a non-transitory computer-readable medium is disclosed having instructions stored thereon, wherein execution of the instructions by a processor causes the processor to determine the navigation trajectory for a vehicle around a small celestial body (e.g., asteroid) by using an on-board camera of the vehicle, wherein the determination employs an inertial equations of motion of the vehicle and of the small celestial body in conjunction with solar gravitational force and the solar radiation pressure forces in a unified factor-graph data model.

In some embodiments, the instructions are configured to implement a front-end image processing system that provides the motion of the vehicle in part using tracked surface features of the image processing system.

In some embodiments, the factor-graph data model includes, at least: a center of mass and principal axis frame of the small body, wherein the principal axis frame of the small body is based on an a-priori unknown mass moment; a center of mass and body-fixed frame of the vehicle; and an inertial frame and point of the solar gravitational force and the solar radiation pressure forces associated with a barycenter of a solar system.

In some embodiments, the factor-graph data model is solved using a Bayesian estimator.

In some embodiments, the unified factor-graph data includes an encoding (e.g., RelDyn) for constraints derived from relative and dynamic kinematics per:

In another aspect, a method is disclosed of operating a spacecraft using the above-discussed computer-readable medium.

In another aspect, a system comprising a controller is disclosed configured to direct navigation of a spacecraft, the controller using the above-discussed computer-readable medium.

In some aspects, the techniques described herein relate to a computer-implemented including: receiving a plurality of images of a small celestial body by a spacecraft, wherein each image of the plurality of images is associated with a time; initializing a plurality of spacecraft states and small celestial body states by the spacecraft; based on the received plurality of images, generating a 3D map of the small celestial body including a plurality of landmarks by the spacecraft, wherein each landmark is associated with a location in the 3D map; generating a factor graph corresponding to the plurality of landmarks and the plurality of spacecraft states and small celestial body states by the spacecraft; receiving a new image by the spacecraft; based on the new image and some or all of the plurality of images, identifying one or more new landmarks in the 3D map and a changed location for at least one of the plurality of landmarks by the spacecraft; updating the factor graph based on the one or more new landmarks and the changed location by the spacecraft; and solving a navigation problem based on the factor graph by the spacecraft to generate relative orientations of the spacecraft and the small celestial body, relative positions of the spacecraft and the small celestial body, and relative velocities of the relative positions of the spacecraft and the small celestial body; performing path planning for the spacecraft by the spacecraft using some or all of the relative orientations of the spacecraft and the small celestial body, the relative positions of the spacecraft and the small celestial body, the relative velocities of the relative positions of the spacecraft and the small celestial bod, and the 3D map; and causing the spacecraft to travel based on the path planning by the spacecraft.

In some aspects, the techniques described herein relate to a computer-implemented method, further including: receiving visual odometry constraints and relative kinematic and dynamical constraints; and solving the navigation problem based on the factor graph and the received visual odometric constraints and the relative kinematic and dynamical constraints.

In some aspects, the techniques described herein relate to a computer-implemented method, wherein the dynamical constraints include the RelDyn factors.

In some aspects, the techniques described herein relate to a computer-implemented, further including receiving prior information about the small celestial object and initializing the plurality of small celestial body states based on the prior information.

In some aspects, the techniques described herein relate to a computer-implemented method, wherein the prior information includes a spin state estimate.

In some aspects, the techniques described herein relate to a computer-implemented method, further including solving the navigation problem based on the factor graph by the spacecraft to generate an estimate of a center-of-mass position of the small celestial body.

In some aspects, the techniques described herein relate to a computer-implemented method, wherein the small celestial body includes an asteroid or a comet.

In some aspects, the techniques described herein relate to a spacecraft including: at least one computing device; and a computer-readable medium with computer-executable instructions stored thereon that when executed by the at least one computing device cause the at least one computing device to: receive a plurality of images of a small celestial body, wherein each image of the plurality of images is associated with a time; initialize a plurality of spacecraft states and small celestial body states; based on the received plurality of images, generate a 3D map of the small celestial body including a plurality of landmarks, wherein each landmark is associated with a location in the 3D map; generate a factor graph corresponding to the plurality of landmarks and the plurality of spacecraft states and small celestial body states; receive a new image; based on the new image and some or all of the plurality of images, identify one or more new landmarks in the 3D map and a changed location for at least one of the plurality of landmarks; update the factor graph based on the one or more new landmarks and the changed location; and solve a navigation problem based on the factor graph to generate relative orientations of the spacecraft and the small celestial body, relative positions of the spacecraft and the small celestial body, and relative velocities of the relative positions of the spacecraft and the small celestial body; perform path planning for the spacecraft using some or all of the relative orientations of the spacecraft and the small celestial body, the relative positions of the spacecraft and the small celestial body, the relative velocities of the relative positions of the spacecraft and the small celestial bod, and the 3D map; and cause the spacecraft to travel based on the path planning.

In some aspects, the techniques described herein relate to a spacecraft, further including computer-executable instructions that when executed by the at least one computing device cause the at least one computing device to: receive visual odometry constraints and relative kinematic and dynamical constraints; and solve the navigation problem based on the factor graph and the received visual odometric constraints and the relative kinematic and dynamical constraints.

In some aspects, the techniques described herein relate to a spacecraft, wherein the dynamical constraints include the RelDyn factors.

In some aspects, the techniques described herein relate to a spacecraft, further including receiving prior information about the small celestial object and initializing the plurality of small celestial body states based on the prior information.

In some aspects, the techniques described herein relate to a spacecraft, wherein the prior information includes a spin state estimate.

In some aspects, the techniques described herein relate to a spacecraft, further including solving the navigation problem based on the factor graph by the spacecraft to generate an estimate of a center-of-mass position of the small celestial body.

In some aspects, the techniques described herein relate to a spacecraft, wherein the small celestial body includes an asteroid or a comet.

In some aspects, the techniques described herein relate to a non-transitory computer-readable medium with computer-executable instructions stored thereon that when executed by the at least one computing device cause the at least one computing device to: receive a plurality of images of a small celestial body, wherein each image of the plurality of images is associated with a time; initialize a plurality of spacecraft states of a spacecraft and small celestial body states; based on the received plurality of images, generate a 3D map of the small celestial body including a plurality of landmarks, wherein each landmark is associated with a location in the 3D map; generate a factor graph corresponding to the plurality of landmarks and the plurality of spacecraft states and small celestial body states; receive a new image; based on the new image and some or all of the plurality of images, identify one or more new landmarks in the 3D map and a changed location for at least one of the plurality of landmarks; update the factor graph based on the one or more new landmarks and the changed location; and solve a navigation problem based on the factor graph to generate relative orientations of the spacecraft and the small celestial body, relative positions of the spacecraft and the small celestial body, and relative velocities of the relative positions of the spacecraft and the small celestial body; perform path planning for the spacecraft using some or all of the relative orientations of the spacecraft and the small celestial body, the relative positions of the spacecraft and the small celestial body, the relative velocities of the relative positions of the spacecraft and the small celestial bod, and the 3D map; and cause the spacecraft to travel based on the path planning.

In some aspects, the techniques described herein relate to a computer-readable medium, further including computer-executable instructions that when executed by the at least one computing device cause the at least one computing device to: receive visual odometry constraints and relative kinematic and dynamical constraints; and solve the navigation problem based on the factor graph and the received visual odometric constraints and the relative kinematic and dynamical constraints.

In some aspects, the techniques described herein relate to a computer-readable medium, wherein the dynamical constraints include the RelDyn factors.

In some aspects, the techniques described herein relate to a computer-readable medium, further including receiving prior information about the small celestial object and initializing the plurality of small celestial body states based on the prior information.

In some aspects, the techniques described herein relate to a computer-readable medium, wherein the prior information includes a spin state estimate.

In some aspects, the techniques described herein relate to a computer-readable medium, further including solving the navigation problem based on the factor graph by the spacecraft to generate an estimate of a center-of-mass position of the small celestial body.

BRIEF DESCRIPTION OF THE DRAWINGS

The foregoing summary, as well as the following detailed description of illustrative embodiments, is better understood when read in conjunction with the appended drawings. For the purpose of illustrating the embodiments, there are shown in the drawings example constructions of the embodiments; however, the embodiments are not limited to the specific methods and instrumentalities disclosed. In the drawings:

FIG. 1 is an illustration of an exemplary environment for navigating a spacecraft in the presence of a small celestial body;

FIG. 2 is an illustration of relative navigation problems frame definitions and vector quantities;

FIG. 3 is an illustration of a factor graph encoding the SLAM problems and the RelDyn factors;

FIG. 4 is an illustration of a full AtsroSLAM factor graph;

FIG. 5 is an illustration of step of the SLAM problem;

FIG. 6 is an illustration of a front-end pipeline factor graph following the initialization of the SLAM problem;

FIG. 7 is an illustration of an operational flow of an implementation of a method for performing path planning for a spacecraft in the presence of a small celestial body; and

FIG. 8 shows an exemplary computing environment in which example embodiments and aspects may be implemented.

DETAILED DESCRIPTION

Some references, which may include various patents, patent applications, and publications, are cited in a reference list and discussed in the disclosure provided herein. The citation and/or discussion of such references is provided merely to clarify the description of the present disclosure and is not an admission that any such reference is “prior art” to any aspects of the present disclosure described herein. In terms of notation, “[n]” corresponds to the n^th reference in the list. All references cited and discussed in this specification are incorporated herein by reference in their entireties and to the same extent as if each reference was individually incorporated by reference.

FIG. 1 is an illustration of an exemplary environment 100 for navigating a spacecraft 102 in the presence of a small celestial body 101. As described above, good navigation or path planning is extremely important for a spacecraft 102 especially when the spacecraft is near a small celestial body 102. Example small celestial bodies 102 (as referred to herein as objects or asteroids) may include asteroids and comets, for example.

To facilitate the path planning, the spacecraft may include a navigation engine 120 and one or more sensors 110. The navigation engine 120 may be implemented using one or more general purpose computing devices such as the computing device 800 illustrated with respect to FIG. 8. Example sensors 110 may include one or more cameras or image generating sensors. These sensors 110 may generate images of the small celestial object 101 and may include monocular visible spectrum images. Other types of images may be supported. Each image may be associated with time or timestamp that reflects when the image was captured by the sensor 110.

To facility path planning, a front-end component of the navigation engine 120 may process images of the small celestial body 101 that are received from the sensors 110. The front-end component may process images to detect features in the images. These features are referred to herein as landmarks. These landmarks on the small celestial body 101 may then be matched between subsequent images as new images are received.

The front-end component of the navigation engine 120 may further initialize one or more states for the small celestial object 101 and the spacecraft 102. These states may include spin-state estimates, velocity estimates, and location estimates. Initially, the front-end component of the navigation engine 120 may initialize the states using prior information known about the small celestial body 101. This prior information may be measurements received from one or more ground-based systems, for example.

The front-end component of the navigation engine 120 may generate a 3D map of the landmarks detected in the images. The 3D map may be of the small celestial object 101 and may be a point cloud generated by triangulating the landmarks.

A back-end component of the navigation engine 120 may generate and maintain a factor graph. The factor graph may be initially generated by the back-end component based on the landmarks in the 3D map, the initialized states of the spacecraft 102, and orbital motion constraints. As will described further below, these orbital motion constraints can be incorporated into the factor graph by devising a relative dynamics factor (“RelDyn”) that can link the relative pose of the spacecraft 102 to the problem of predicting trajectories stemming from the motion of the spacecraft 102 in the vicinity of the small celestial body 101.

The back-end component of the navigation engine 120 receives a new image of the small celestial body 101 and based on the locations of the landmarks in the received image and the factor graph, performs an optimization to generate a solution that includes estimates of past relative positions of the positions of the small celestial body 101 and the spacecraft 102, past relative orientations of the small celestial body 101 and the spacecraft 102, and past relative velocities of the small celestial body 101 and the spacecraft 102. This data, and the updated 3D map, may be used to provide path planning for the spacecraft 101 with respect to the small celestial body 101. The details of the proposed system and methods are described further below.

1. Problem Statement

In this section, are discussed the relevant theory and the problem statement pertaining to the asteroid relative navigation problem incorporating monocular SLAM, an appropriate motion prior and sensor fusion.

Firstly, the notations and conventions used throughout the disclosure are summarized in Section 1.1. Secondly, the defined notation is contextualized within the problem of a spacecraft navigating around an asteroid in Section 1.2.

1.1. Notation

Given the affine space (³, ³), the translation vector between any two points X, Y∈³ is denoted r_YX(Y−X)∈³, read “from X to Y”. Any frame is a tuple (X, {{right arrow over (x)}_i}_i=1³) where the point X∈³ denotes the origin of the frame and the set of unit directions {{right arrow over (x)}_i}_i=1³, where {right arrow over (x)}_i∈², i=1, . . . , 3, constitutes the right-handed orthonormal basis of the frame. The expression of any vector v∈³ in a given frame =(*, {{right arrow over (x)}_i}_i=1³) is denoted by v[v·{right arrow over (x)}₁ v·{right arrow over (x)}₂ v·{right arrow over (x)}₃]^T∈³. For any two frames and , we denote the rotation [{right arrow over (y)}₁^x {right arrow over (y)}₂^x {right arrow over (y)}₃^x]^T∈SO(3) such that v= for any v∈³. This notation is consistent with the composition rule = for any three frames . Given a coordinate vector v∈³ expressed in frame , we denote the corresponding homogeneous coordinates v[(v)^T 1]^T∈³. For any two frames =(X, {{right arrow over (x)}_i}_i=1³) and y=(Y, {{right arrow over (y)}_i}_i=1³), we denote the homogeneous transformation ∈SE(3) by:

$\begin{array}{cc}\stackrel{\Delta }{=}\left[\begin{array}{cc}& r_{\mathrm{YX}}^{\varkappa }\\ 0_{1\times 3}& 1\end{array}\right],& \left(1\right)\end{array}$

such that r_PX= for any point P∈. This notation is consistent with the composition rule T= for any three frames

Given a Lie algebra of dimension n associated (at the identity) to a matrix Lie group , and given a set of basis vectors {E_i}_i=1ⁿ of the matrix Lie algebra (n), we denote the hat operator [⋅]{circumflex over ( )}:ⁿ→ which maps any n-vector x=[x₁ . . . x_n]^T∈ⁿ to an element [x]{circumflex over ( )}∈, by [x]{circumflex over ( )}=Σ_i=1ⁿ x_iE_i. We denote its inverse vee operator [⋅]^∨:→ⁿ, which extracts the coordinates x=[x₁ . . . x_n]^T from [x]{circumflex over ( )}∈(n) in terms of E_i. We denote the exponential map exp: →, mapping an element [x]{circumflex over ( )}∈ g to the element exp([x]{circumflex over ( )})∈ in the neighborhood of the identity element. Specifically, for the SO(3) group of rotations, the exponential map at the identity exp: so(3)→SO(3) which associates any tangent vector (skew-symmetric matrix) [x]{circumflex over ( )}∈so(3) to a 3D rotation as the matrix exponential, is given by

$\begin{array}{cc}\exp \left(\left[x\right]^\right)=I_3+\left[x\right]^+\frac{1}{2}{\left(\left[x\right]^\right)}^2+\frac{1}{6}{\left(\left[x\right]^\right)}^3+\dots =I_3+\frac{\sin x}{x}\left[x\right]^+\frac{1-\cos x}{{x}^2}{\left(\left[x\right]^\right)}^2.& \left(2\right)\end{array}$

Finally, we denote the logarithm map (at the identity) by log: →, which maps an element in the neighborhood of the identity element of group to an element in the associated Lie algebra . In the case when =SO(3) and =so(3), the logarithm map log: SO(3)→so(3) is the inverse of the operation given in (2), and is a bijective mapping as long as ∥x∥<π[Chirikjian (2011)].

The following facts are useful in deriving first-order approximations for the effect of noise on stochastic differential equations, as well as for deriving partial derivatives of functions of elements of SO(3).

Fact 1 Given the hat operator [⋅]{circumflex over ( )}:³→so(3), for any R∈SO(3), and coordinate vector x∈³,

R[x]{circumflex over ( )}R ^T =[Rx]{circumflex over ( )},  1.

R[x]{circumflex over ( )}=[Rx]{circumflex over ( )}R,  2.

[x]{circumflex over ( )}R=R[R^Tx]{circumflex over ( )}.  3.

Note that Fact 1 is a consequence of the rotational invariance of the cross product in ³:

Fact 2 Given the exponential map exp: so(3)→SO(3), for any R∈SO(3) and any coordinate vector x∈³,

R exp([x]{circumflex over ( )})R^T=exp([Rx]{circumflex over ( )}),  1.

R exp([x]{circumflex over ( )})=exp([Rx]{circumflex over ( )})R,  2.

exp([x]{circumflex over ( )})R=R exp([R^Tx]{circumflex over ( )}).  3.

Note that Fact 2 is a direct consequence of Fact 1 applied to the matrix exponential expansion in (2).

Fact 3 Given the exponential map exp: so(3)→SO(3), for any vector x∈³ having sufficiently small length, it follows that:

exp([x]{circumflex over ( )})≃I₃+[x]{circumflex over ( )}.  (3)

Note that Fact 3 is a consequence of dropping higher order terms in the matrix exponential expansion of equation (2).

Fact 4 Let T: SO(3)׳→SE(3) be function that produces the pose corresponding to a rotation matrix and a translation vector pair, such that:

$T\left(R,r\right)\stackrel{\Delta }{=}\left[\begin{array}{cc}R& r\\ 0& 1\end{array}\right],R\in \mathrm{SO}\left(3\right),r\in {ℝ}^3.$

Define the perturbation vectors δϕ, δr∈³. We may write the perturbed pose

$T\left(R\exp \left[\delta \varphi \right]^,r+\delta r\right)=\left[\begin{array}{cc}R\left(\exp \left[\delta \varphi \right]^\right)& r+\delta r\\ 0& 1\end{array}\right].$

Multiplying the inverse of the pose T(R, r) on the left, we obtain the perturbation pose:

$\begin{array}{c}\delta T={\left(T\left(R,r\right)\right)}^{-1}T\left(R\exp \left[\delta \varphi \right]^,r+\delta r\right)\\ =\left[\begin{array}{cc}{R}^{\top }& -R^{\top }r\\ 0& 1\end{array}\right]\left[\begin{array}{cc}R\exp \left(\left[\delta \varphi \right]^\right)& r+\delta r\\ 0& 1\end{array}\right]\\ =\left[\begin{array}{cc}\exp \left(\left[\delta \varphi \right]^\right)& R^{\top }\delta r\\ 0& 1\end{array}\right]\end{array}.$

Taking the logarithm map of the SE(3) followed by the vee operator, we obtain the corresponding pose perturbation vector:

$\delta P={\log \left(\delta T\right)}^{\bigvee }=\left[\begin{array}{c}\delta \varphi \\ R^{\top }\delta r\end{array}\right],$

thus yielding the Jacobian:

$J_R=\left[\begin{array}{c}{I}_3\\ 0_{3\times 3}\end{array}\right],\mathrm{and}{J}_r=\left[\begin{array}{c}{0}_{3\times 3}\\ R^{\top }\end{array}\right].$

It then follows that:

$\left[\begin{array}{c}\delta \varphi \\ \delta r\end{array}\right]=\left[\begin{array}{cc}{I}_3& 0\\ 0& R\end{array}\right]\delta p.$

1.2 Problem Definitions

Let A, S, O, W∈³ and assume that the point A corresponds to the center of mass of the small celestial body, S corresponds to the center of mass of the spacecraft, O corresponds to the center-of-mass of the Sun, and W corresponds to an inertial point in space, e.g., non-accelerating, constant velocity. Initially, we distinguish three frames of interest: the inertial frame (W, {{right arrow over (n)}_i}_i=1³), and the spacecraft body-fixed frame (S, {{right arrow over (s)}_i}_i=1³), and an arbitrarily chosen small body-fixed frame (G, {{right arrow over (g)}_i}_i=1³). Note that since frame is fixed with respect to small celestial body, and assuming that the small body is a rigid body, the position of the small celestial body's center of mass point A with respect to the point G in frame coordinates, denoted , is fixed.

Given an inertial point W, the inertial absolute position vector of the observing spacecraft is denoted r_SW and that of the small celestial body r_AW. The relative position vector of the observing spacecraft is then given as r_SA=r_SW−r_AW. The inertial relative velocity vector of the observing spacecraft is denoted by

$v_{\mathrm{SA}}\stackrel{△}{=}{ }^{}\frac{d}{\mathrm{dt}}\left(r_{\mathrm{SA}}\right).$

In the same fashion, the inertial relative acceleration is denoted

$a_{\mathrm{SA}}\stackrel{△}{=}{ }^{}\frac{d}{\mathrm{dt}}\left(v_{\mathrm{SA}}\right).$

To later exploit knowledge of the relative orbital pose, we decompose the relative position vector r_SA using the intermediary point G fixed in the g-frame coordinates, such that r_SA=r_SG+r_GA=r_SG−r_AG. Expressing in the body-fixed g-frame coordinates and writing as a function of the relative position vector, we get the relative pose translation vector:

++  (4)

Any reference hereafter to the relative pose of the spacecraft 102, for a given spacecraft frame =(S, {{right arrow over (s)}}_i=1³), designates the transformation:

$\begin{array}{cc}=\left[\begin{array}{cc}& \\ 0_{1\times 3}& 1\end{array}\right],& \left(5\right)\end{array}$

which encodes the relative rotation []^T∈SO(3) of the spacecraft 102 with respect to the frame and the coordinates =[r_sg·{right arrow over (g)}_i r_sg·{right arrow over (g)}₂ r_sg·{right arrow over (g)}₃]∈³ of the spacecraft position vector relative to the point G as expressed in the frame. FIG. 2 provides an illustration 200 of the problem's relevant points, vectors, and frames, as defined above.

1.3. Scene Mapping Considerations

Let t₀ be the initial time, let t≥t₀, let (t_k)_k=0ⁿ⊂[t₀, t] be the sequence of sensor acquisition times, and let (t_k) describe the pose of the spacecraft 102 as expressed in the frame at each time index 0≤k≤n. Then, the sequence ()_k=0ⁿ∈Π_k=1ⁿSE(3) describes the discrete trajectory of the relative pose of the spacecraft 102. We define the camera sensor frame =(S, {{right arrow over (c)}_i}_i=1³), with fixed pose with respect to the spacecraft frame , and obtain the sequence ()_k=1ⁿ of all camera poses, also known as frames, where the relative camera pose is ==, (k=0, . . . , n).

A landmark L∈³ is defined as a notable 3D point in the scene, which is potentially triangulated during SLAM using camera observations. We denote by Ψ_k={L_i∈³, i=1, . . . , m_k} the set of all landmarks accumulated up until time index k=0, . . . , n, also called the map at time index k. To each landmark L∈Ψ_k corresponds a position vector r_LG∈³, whose coordinates in the frame, denoted , are fixed since the small body is presumed to be a rigid body. Let r_k be the total number of image feature points detected in the camera image captured at time t_k, k=0, . . . , n. We collect all detected feature points in the set γ_k={P_i∈², i=1, . . . , r_k} where, to each P_i∈γ_k are associated the ideal 2D image coordinates y_ik∈².

Assume a point P_i∈γ_k corresponds to the image projection of a scene landmark L_j∈Ψ_k, as captured at time index k. The 2D image coordinates y_ik relate to the 3D position coordinates through the pinhole camera model relationship, given by:

$\begin{array}{cc}\left[\begin{array}{c}\lambda y_{\mathrm{ik}}\\ \lambda \end{array}\right]=K,& \left(6\right)\end{array}$

where λ>0 is a scaling factor, and where:

$K\stackrel{△}{=}\left[\begin{array}{cccc}{f}_x& 0& c_x& 0\\ 0& f_y& c_y& 0\\ 0& 0& 1& 0\end{array}\right]$

is the camera intrinsic matrix, with f_x, f_y, c_x, c_y scalars corresponding to the known camera focal lengths and optical center offsets along the two image dimensions. The real measured feature point coordinates y_ik^m are then defined such that:

y _ik ^m =y _ik+ν_y,  (7)

where ν_y˜(0, Σ_y^m) corresponds to the feature point measurement noise, with associated 2×2 covariance matrix Σ_y^m. We collect all real camera frame-landmark observations up until time index k=0, . . . , n in the set {y_ij^m∈²: i=1, . . . , r_j=0, . . . , k}. We also define }∈SE(3): i=0, . . . , k} as the set of all possible spacecraft relative poses discretized at times (t_k)_k=0ⁿ, and define _k{∈³: L∈Ψ_k} as the set of all landmark coordinates mapped up until time index k=0, . . . , n.

By exploiting the multi-view geometry constraints derived by capturing observations _k of the landmarks Ψ_k at poses _k, as well as the constraints derived from the intrinsic motion of the spacecraft 102 around the small body 101 and other sensor measurements, further detailed in Section 3, we wish to find a solution to the trajectory _k along with the set of mapped landmark coordinates _k, on-the-fly, for k=0, . . . , n in a sequential and incremental manner. The method employed to solve this problem is described in Section 4. First, we provide a brief primer on the Bayesian estimation, the SLAM problem formulated in the Bayesian framework, and the relevance of representing the structure of the problem using a factor graph.

2. Estimating Relative Motion Using IMU Measurements

We now further explain the improvements of our work compared to prior art methods, specifically: (a) why the IMU preintegration scheme used in prior works is inadequate and why relative dynamics need to be considered; and (b) why the kinematic rotation factor used in prior art methods may be modified to account for non-inertial and large motion of the target center-of-mass point.

2.1. Using IMU Accelerometer Measurements to Perform Dead-Reckoning in the Presence of Unknown Gravity Vector

In the scenario described herein: (1) in terms of gravitational forces, the observing spacecraft 102 is affected by the combined pull of the target small celestial body 101 and of the Sun, while the small celestial body 101 is only affected by the gravitational acceleration of the Sun. This results in a net non-negligible relative gravitational interaction between the “chief” (the small celestial body 101) and the “deputy” (the observer S/C), which, in turn, affects the relative motion of the observer spacecraft. However, an IMU accelerometer is insensitive to gravitational forces. The double integrator dynamics do not hold over extended arcs of time, and a full non-linear treatment provides higher precision; (2) there are persistent external non-gravitational forces acting on the observer spacecraft 102 which are small in magnitude, and consequently may not be detectable given the signal-to-noise ratio (SNR) of an IMU measurement signal, but which may still affect the long-term trajectory of the observing spacecraft over an inspection arc. The IMU measurements required for the IMU factors incorporated in prior art methods may not capture the accumulated effect of such small-scale non-gravitational forces; (3) unlike spacecraft inspection arcs, the duration of small celestial body inspection arcs is typically extended—that is, their length is a non-negligible portion of the period of revolution of the spacecraft around the small celestial body. Assuming that the center-of-mass of the target body only moves negligibly between time instances, or that its motion is inertial throughout the inspection arc, may induce considerable relative navigation errors during estimation.

When we analyze the SPHERES-VERTIGO platform-inspired problem formulation, we note that the inspector spacecraft (deputy) and the ISS habitat module in which the experiment is run (chief) are both subjected to the same Earth gravitational acceleration which maintains them in orbit, without any other gravitational interaction between the two. The ISS frame, being nadir-pointing, may be assimilated to the Hill frame, and consequently the relative motion of the observer spacecraft with respect to the Hill frame is fundamentally described by the Clohessey-Wiltshire (CW) equations of motion. We note that the target spacecraft (another deputy) has a similar relationship to the ISS frame (chief), but its coordinate is coincident with that of the Hill frame and assumed to be inertial. In this specific context, given the very small difference in Earth-centered inertial position between the chief (ISS module) and the deputy (observer) spacecraft, and the short time-frame of the considered experiment, the Hill frame may be considered as quasi-inertial and the motion of the observer may indeed be approximated by the solution to double-integrator dynamics. Only in this limited scenario, the simple accelerometer model where gravity is set to zero and acceleration measurements are simulated as thruster-specific forces, seems valid.

In the latter case, the accelerometer measurements may be used to perform “dead-reckoning,” that is, to compute the change in the inertial velocity and the inertial position of the observer spacecraft with respect to the quasi-inertial world frame due to the influence of external accelerations, such as the action of thruster firing. Such dead-reckoning may be performed using a pre-integration factor, which attempts to encode the residual between the current state estimate and the state obtained by preintegration of these accelerometer measurements starting from the previous navigational state.

However, as is generally known that an IMU's accelerometer does not measure pure accelerations, but rather externally applied non-gravitational specific forces, that is, inertial accelerations sans gravity-induced accelerations.

The IMU measurements do not provide any information about the direction or the magnitude of the gravity vector in freefall (unforced orbital) motion. We can intuit this fact by considering a simple model for the measurement of an accelerometer centered at point M and with frame , given as:

=−)+ν,  (8)

where O is an inertial point, is an inertial frame, c is a measurement sensitivity constant, is the rotation between frames and , and ν is a Gaussian noise term. It follows that, in an unforced (freefall) motion, while assuming the IMU frame is coincident with the spacecraft body frame , only gravity affects the IMU, i.e. =, and thus the accelerometer measurement over time is null on average. In other words, the IMU accelerometer is unusable for the purpose of generating estimates of inertial velocity and change in inertial position by dead-reckoning in the presence of unknown gravitational accelerations. To obtain an accurate estimate of inertial motion, needs to be known.

2.2. Using IMU Accelerometer Measurements to Perform Dead-Reckoning of Motion Owing to Small-Scale Persistent Perturbing Forces

Admittedly, accelerometer measurements may be used in our work to obtain estimates of motion owing only to non-gravitational external forces being applied on the S/C, such as actuated thruster profiles. However, such forces need to be strong enough to generate a large enough signal-to-noise ratio to be useful. One common external force which does not meet the latter requirement, but whose effect on the trajectory is non-negligible over the length of an inspection arc, is solar radiation pressure (SRP). Its magnitude is typically orders of magnitude smaller than the gravitational effect, and it does not meet the criterion of having enough SNR to be detectable by the accelerometer.

In light of these caveats, incorporating IMU factors to process such IMU measurements, given that we need to augment the state with the additional IMU bias states and incur additional computation time, may not be justifiable. Therefore, an admissible “odometric” factor for our scenario would then have to rely on explicitly modeling and tracking the gravity vector, to ascertain the true relative motion between the observer spacecraft and the target small celestial body. Moreover, since the center-of-mass of the small celestial body may not be considered as inertial owing to the gravitational pull by the Sun over the inspection arc, and its inertial velocity may not be known with good precision during proximity operations, it is preferable to model the relative dynamics of the spacecraft-small celestial body pair, by avoiding the modeling of the inertial dynamics of both spacecraft and target small celestial body separately.

3. Relative Kinematic and Orbital Odometric Factors

In this section, we detail the theoretical approach for devising an adaptation of Setterfield's kinematic rotation factor as well as our new RelDyn factor, which enforces a strong odometric constraint, namely, the equations governing the motion of the spacecraft relative to the small celestial body, and which incorporates fusion of other sensor measurements. We also discuss the relevant sources of noise, and the propagated system of equations.

3.1. Factor Graph Encoding and ISAM2 Algorithm

In modern renditions of SLAM, the problem stated in Section 1.2 is formulated using a probabilistic inference framework, predicated on Bayesian estimation, which we evoke for our solution in this section. The use of the factor graph formulation for modeling the small body navigation problem is motivated by the fact that the typical Bayesian estimation method used to solve the SLAM estimation problem is amenable to a graph representation due to the sparsity in the structure of the cost. Factor graphs are un-directed bi-partite probabilistic graphical models constituted of factor nodes and variable nodes, with edges connecting variable-factor pairs. The structure captured by the edges and nodes of the factor graph encodes the structure of the estimation problem's posterior probability density function, by exploiting the fact that the latter can be factorized as a product of many functions, each depending on a subset of the variables of the problem. By exploiting sparsity in the structure of the joint density function, the factor graph formulation can render very large estimation problems tractable in terms of computation. Factor functions can be derived and emplaced in the factor graph based on the problem-specific constraints which we wish to include.

Given a constructed factor graph and a specific variable ordering, variable elimination is performed, transforming the factor graph into a chordal Bayes net. We note that, in theory, the chordal Bayes net may then be written as an equivalent square root information matrix, which may then be solved using algebra (back-substitution) when the factors are linear, and by Gauss-Newton iteration when the factors are non-linear. The procedure creates the resulting Gauss-Newton iteration directly from the factor graph and then attempts to solve for updated values of the variable set.

However, in more recent and efficient renditions of factor graph solvers, Bayes trees are used to store internally the structure of the problem, while the factor graph is used principally as an intermediary for modeling. A Bayes tree is a type of directed junction tree in which the nodes store the cliques of the chordal Bayes net it represents and the edges are the separator set of variables which separate the cliques.

In iSAM2, specifically, after a first step variable elimination step is performed on the initial factor graph, a Bayes tree is constructed from the chordal Bayes net and stored for the next time step, contrary to the classical approach of storing the factor graph itself for the next stage. When a new subgraph relating to a set of new variables and measurements is available for insertion, iSAM2 updates the Bayes tree by (1) converting the part of the Bayes tree relating to this new subgraph into a factor graph, (2) appending the new factor subgraph—usually at the top of the tree, i.e. affecting the most recent cliques—then (3) eliminating the factor subgraph to produce a new Bayes net, and (4) reassesses the affected cliques, (5) updates the appropriate variables. iSAM2 updates variables involved in non-linear factors by performing a non-linear optimization step at a current linearization point. To select which variables undergo optimization update, iSAM2 first marks all variables in a subset of the graph variables—that is, the ones being added to the graph—which clear the “wildfire” threshold to update, after which a new linearization point is computed and subsequent all cliques which involve the updated variable get marked for update as well.

These points of the procedure are important in that they dictate the effect on the iSAM2 runtime and complexity when we insert variables which may render the graph fully connected, greatly affecting the size of the cliques marked for update and optimization. We further discuss these in Section 4.1. Next, we devise the factor which serves as a motion constraint in the small body circumnavigation problem.

3.2. The Relative Kinematic Factor Formulation

A kinematic rotation factor is incorporated it in the problem factor graph. Unlike the kinematic rotation factors, however, we do not assume that the target's (i.e., the body 101) center-of-mass point A is inertial or that the motion of the target's center-of-mass between two successive time indices is negligible. Therefore, for any two time indices i, j, i≠j, r_AjAi≠0, and we can write:

r _{S j S i} =r _{S j G j} +r _{G j A j} +r _{A j A i} +r _{A i G i} +r _{G i S i} ,  (9)

We re-write the left hand side of equation (9) by transiting through points A_j and A_i in the vector chain, such that:

r _{S j S i} =r _{S j A j} +r _{A j A i} +r _{A i S i} ,  (10)

We replace equation (10) into the equation (9), cancel out terms, and manipulate to obtain the kinematic relationship, stated as:

r _{S i A j} −r _{S i A i} =r _{S j G j} +r _{G j A j} +r _{A i G i} +r _{G i S i} .  (11)

We note that in lieu of a inertial position vectors r_SiW with respect to an inertial world frame origin W used is other works, we have the relative vector r_SiAi. While acknowledging that the center-of-mass vector r_AG is fixed when expressed in the -frame coordinates, such that

$r_{A_i{G}_i}^{{𝒢}_i}=r_{A_j{G}_j}^{{𝒢}_j}=r_{\mathrm{AG}}^{𝒢},$

and expressing all vectors in the inertial frame, we get the kinematic relationship:

$\begin{array}{cc}-r_{S_i{A}_i}^{𝒥}=\left(-\right)+R_{{\mathrm{𝒥𝒮}}_i{R}_{{𝒮}_i{𝒢}_i}}\left(r_{\mathrm{AG}}^{𝒢}-\right).& \left(12\right)\end{array}$

Further assuming that we have access to a filtered estimate of the inertial attitude =exp([ν_R]{circumflex over ( )}), ν_R,i˜(0, τ_R) from the spacecraft's inertial star tracker at any desired time, we may re-write equation (12) as:

$\begin{array}{cc}-r_{S_i{A}_i}^{𝒥}=\exp \left(\left[-v_{R,j}\right]^\right)\left(-\right)+\exp \left(\left[-v_{R,i}\right]^\right)\left(r_{\mathrm{AG}}^{𝒢}-r_{S_i{G}_i}^{{𝒢}_i}\right).& \left(13\right)\end{array}$

Using Facts 1 and 3, and manipulating equation (13), we obtain the translation residual:

$\begin{array}{cc}{\epsilon }_{\mathrm{ij}}^{\mathrm{RelKin},r}\stackrel{\Delta }{=}--\left(-\right)-\left(-\right)=\left[{\hat{R}}_{{\mathrm{𝒥𝒮}}_j}{R}_{{𝒮}_j{𝒢}_j}\left(r_{S_j{G}_j}^{{𝒢}_j}-r_{\mathrm{AG}}^{𝒢}\right)\right]^{\hat{R}}_{{\mathrm{𝒥𝒮}}_j}{v}_{R,j}+\left[{\hat{R}}_{{\mathrm{𝒥𝒮}}_i}{R}_{{𝒮}_i{𝒢}_i}\left(r_{\mathrm{AG}}^{𝒢}-r_{S_i{G}_i}^{{𝒢}_i}\right)\right]^{\hat{R}}_{{\mathrm{𝒥𝒮}}_i}{v}_{R,i}.& \left(14\right)\end{array}$

While assuming that [ν_R,iν_R,j^T]=0_3×3, i≠j, and using the shorthand

$B_i\stackrel{\Delta }{=}\left[\left(-\right)\right]^,$

an expression for the covariance of ϵ_ij^RelKin,r may be found by computing:

[ϵ_ij^RelKin,r(ϵ_ij^RelKin,r)^T]=B_iΣ_RB_i^T+B_jΣ_RB_j^T.  (15)

In parallel, we can write the rotation kinematic relationship using the rotation composition rule, given as:

=.  (16)

We rewrite the left-hand side of the equation using the intermediary frame, such that:

=.  (17)

Using an Euler on-manifold integration scheme with the simplifying assumption of piecewise constant (t)=, t∈[t_i, t_j), a kinematic relationship between and holds, such that:

=exp([(t_j−t_i)]{circumflex over ( )}).  (18)

Substituting equation (18) into equation (17), and then substituting that result into equation (16), we obtain the kinematic relationship:

exp([(t_j−t_i)]{circumflex over ( )})=.  (19)

Further assuming that we have access to a filtered estimate of the inertial attitude =exp([ν_R,i]{circumflex over ( )}), ν_R,i˜(0, Σ_R) from the spacecraft's inertial star tracker, we may re-write equation (19) as:

exp([(t_j−t_i)]{circumflex over ( )})=exp([ν_R,i]{circumflex over ( )})exp([−ν_R,j]{circumflex over ( )}).

Using Fact 2, and manipulating the right-hand side of the latter equation, we obtain:

exp([(t_j−t_i)]{circumflex over ( )})=exp([νR,i]{circumflex over ( )})exp([−ν_R,j]{circumflex over ( )}).

Taking the logarithm map followed by the vee operator on both sides of the equation, and then applying the Baker-Campbell-Hausdorff formula while omitting any second-order noise term, we obtain the residual:

ϵ_ij^RelKin,Rlog(exp([(t_j−t_i)]{circumflex over ( )}))^∨=ν_R,i−ν_R,j.  (20)

Since (a) rotation matrices are full rank and norm preserving, (b) the star tracker measurement noise is typically isotropic, e.g. Σ_R=σ_R²I₃, and (c) that the star tracker measurement noises ν_R,j,ν_R,i are independent for i≠j, that is [ν_R,jν_R,i^T]=0_3×3, it follows that ϵ_ij^RelKin,R˜(0,2τ_R).

We combine the two residuals to obtain the full relative kinematic factor:

$\begin{array}{cc}{ϵ}_{\mathrm{ij}}^{\mathrm{RelKin}}=\left[\begin{array}{c}{ϵ}_{\mathrm{ij}}^{\mathrm{RelKin},R}\\ {ϵ}_{\mathrm{ij}}^{\mathrm{RelKin},r}\end{array}\right]~\left(0,\left[\begin{array}{cc}2{\sum }_R& 0\\ 0& B_i{\sum }_R{B}_i^{\top }+B_j{\sum }_R{B}_j^{\top }\end{array}\right]\right).& \left(21\right)\end{array}$

We incorporated this adapted factor into our factor graph formulation.

3.3. The RelDyn Factor Formulation

We derive a relative dynamics factor based on the moments of the distribution describing the dispersion of the solution realizations of a system of nonlinear stochastic differential equations, evaluated at discrete times, for the purpose of improving the performance of baseline SLAM for small celestial body imagery.

Given a sequence of spacecraft states Ξ_n=(x_k)_k=0ⁿ at discrete times (t_k)_k=0ⁿ, for k=1, . . . , n, we want to devise factors which relate to the odometric probability distribution p(x_k,x_k+1). Note that, using Bayes' rule, the joint distribution p(x_k,x_k+1) can be rewritten as either p(x_k,x_k+1)=p(x_k+1|x_k)p(x_k) or p(x_k,x_k+1)=p(x_k|x_k+1)p(x_k+1), depending on which prior, p(x_k) or p(x_k+1), is readily available at the time of computation. In our case, we typically know the state of the spacecraft and the distribution of the uncertainty of that state at the beginning of a segment, and therefore we pick p(x₀)˜({circumflex over (x)}₀,Σ₀) as a known prior, where {circumflex over (x)}₀ is the known state mean and Σ₀ is the known state covariance. We associate the factor ϕ₀^prior(x₀) to the prior p(x₀) accordingly. Now, the RelDyn factor ϕ_k^RelDyn can be formulated as:

ϕ_k^RelDyn(x_k,x_k+1)∝p(x_k+1|x_k)=(0,P_k),  (22)

where P_k is a covariance matrix derived from the propagation of the moments of the distribution describing the realizations of the stochastic differential equations.

Along with projection factors, denoted ϕ_i^proj, i=1, . . . , r_k, a factor graph of the SLAM problem with motion priors can be constructed, as conceptually illustrated in the graph 300 of FIG. 3.

In light of this goal, a brief discussion of the RelDyn factor is provided throughout this section. To compute the value of the RelDyn, we need to model the equations of the relative motion induced by the relative navigation problem in the presence of gravitational attraction between the observer and target bodies.

3.4. Equations of Relative Motion

Recalling the vectors and frames od FIG. 2, we assume that the frame is inertially fixed, the frames and are rotating with associated inertial angular velocity vectors , and , and that frame is a target small celestial body-fixed frame. The relative angular velocity between the spacecraft frame and frame is =−. The orientation of the frame relative to the frame is encoded in the rotation matrix ∈SO(3), which satisfies the kinematic relationship:

=[]{circumflex over ( )},  (23)

where ∈³ is the relative angular velocity, expressed in the g frame. Rewriting (23) in terms of the frame and frame angular velocities, and using Fact 1, we obtain:

=[−]{circumflex over ( )}=[]{circumflex over ( )}−[]{circumflex over ( )}.  (24)

We can derive the position vector of the spacecraft at point S relative to the arbitrary point G by writing:

r _SG =r _SA +r _AG.

Taking the time derivative in the g frame gives

${ }^{}\frac{d}{\mathrm{dt}}\left(r_{\mathrm{SG}}\right)={ }^{}\frac{d}{\mathrm{dt}}\left(r_{\mathrm{SA}}\right)+{ }^{}\frac{d}{\mathrm{dt}}\left(r_{\mathrm{AG}}\right)={ }^{}\frac{d}{\mathrm{dt}}\left(r_{\mathrm{SA}}\right)-{\omega }_{}\times r_{\mathrm{SA}}.$

Expressing in the coordinates, we get:

=−[]{circumflex over ( )}.  (25)

Finally, we derive the translational dynamics of the spacecraft relative to the small body center of mass. Our equations of motion for the relative position and relative velocity are then given as:

=  (26)

=.  (27)

It thus suffices to find an expression for .

We assume herein that the spacecraft 102 is subjected to the gravitation force of the small body 101, denoted by F_sa, the gravitation of the Sun, denoted by F_s⊙, the solar radiation pressure (SRP), denoted by F_SRP, as well as the spacecraft actuation thrust force, denoted by F_s. Assuming that the spacecraft's known mass, denoted by m_s, is fixed, the linear acceleration of the spacecraft with respect to the Sun's origin O, is given by a_SO=1/m_s(F_sa+F_s⊙+F_SRP+F_s). The small body's gravitation force is obtained by computing:

$\begin{array}{cc}{F_{\mathrm{sa}}=\frac{\partial U\left(r\right)}{\partial r}❘}_{r=r_{\mathrm{SA}}}& \left(28\right)\end{array}$

For an appropriate gravity field potential function U(r). Assume that U(r) is parameterized using spherical harmonics. Then, when the probe 102 is relatively distant from the small body 101, the spherical term of the potential dominates, in which case the attractive force is given by F_sa=−μ_am_s/∥r_SA∥³r_SA. Consider the Sun as a point-mass central body. Then, F_s⊙=−μ_⊙m_s/∥r_SO∥³r_SO, where r_SO is the position vector of the spacecraft 102 center of mass S with respect to the Sun origin O. Assume that the solar radiation pressure is a function of its position vector with respect to the Sun, F_SRP=F_SRP(r_SO).

In turn, we assume that the mass of the small body 101, denoted m_a is fixed, and that the only force acting on the small body 101 is the Sun's gravitational force. Then, the linear acceleration of the small body center of mass with respect to the Sun's origin O, is given by a_AO=1/m_aF_a⊙, where F_a⊙=−μ_⊙m_a/∥r_AO∥³r_AO. The relative dynamics of the spacecraft-small celestial body system are obtained computing a_SA=a_SO−a_AO, leading to the relationship:

$\begin{array}{cc}{a}_{\mathrm{SA}}=-\frac{{\mu }_a}{{r_{\mathrm{SA}}}^3}{r}_{\mathrm{SA}}-\frac{{\mu }_{\odot }}{{r_{\mathrm{SO}}}^3}{r}_{\mathrm{SO}}+\frac{1}{m_s}{F}_{\mathrm{SRP}}\left(r_{\mathrm{SO}}\right)+\frac{1}{m_s}{F}_s+\frac{{\mu }_{\odot }}{{r_{\mathrm{AO}}}^3}{r}_{\mathrm{AO}}.& \left(29\right)\end{array}$

Assume that r_SO≅r_AO, given the very large distance between the spacecraft-small celestial body system and the Sun, and rewrite equation (29) by making explicit in terms of the state variables and input variables of interest, yielding:

$\begin{array}{cc}{a}_{\mathrm{SA}}^{𝒥}=-\left(\frac{{\mu }_a}{{r_{\mathrm{SA}}}^3}+\frac{{\mu }_{\odot }}{{r_{\mathrm{AO}}}^3}\right)r_{\mathrm{SA}}^{𝒥}+\frac{1}{m_s}{F}_{\mathrm{SRP}}^{𝒥}\left(r_{\mathrm{AO}}^{𝒥}\right)+\frac{1}{m_s}{R}_{\mathrm{𝒥𝒮}}{F}_S^{𝒮}.& \left(30\right)\end{array}$

Additionally, for the sake of readability, we hereafter define the shorthand notation Q q=, w r, v, R, s, d, c, g(d)1/m_s(d), f1/m_sF_s, which allow us to restate the equations of motion as:

$\begin{array}{cc}\stackrel{.}{Q}=Q\left[s\right]^-\left[w\right]^Q,& \left(31\right)\end{array}$ $\begin{array}{cc}\stackrel{.}{q}={\mathrm{QR}}^{\top }v-\left[w\right]^{\mathrm{QR}}^{\top }r& \left(32\right)\end{array}$ $\begin{array}{cc}\stackrel{.}{v}=-\left(\frac{{\mu }_a}{{r}^3}+\frac{{\mu }_{\odot }}{{d}^3}\right)r+g\left(d\right)+\mathrm{Rf}& \left(33\right)\end{array}$ $\begin{array}{cc}\stackrel{.}{r}=v.& \left(34\right)\end{array}$

3.5. Modeling Assumptions

We note that we may view R and s, along with f, as inputs to the dynamical system described in (31)-(34). Here we consider R={circumflex over (R)}exp([ν_R]{circumflex over ( )}) as the already-filtered estimate value {circumflex over (R)} of the spacecraft's orientation provided by the inertial star tracker instrument, perturbed on the right by its associated uncertainty ν_R, and s=ŝ+ν_s is the filtered and unbiased estimate value blacks of its angular velocity, provided by the rate gyro instrument, perturbed by the associated uncertainty ν_s. This inclusion allows us to fuse the known filtered measurements of a star tracker system and a rate gyro into the overall navigation solution. Moreover, attitude and angular rate filtered measurements may be obtained at high rate and may therefore be associated in time so to closely coincide with the moment an image is taken and processed.

By modeling the sensor dynamics we avoid having to represent the Newton-Euler dynamics of the spacecraft 102. Indeed, attitude control inputs, such as the actual actuated torques applied on the spacecraft, are often not known with precision at the instant when they are applied on the spacecraft 102, but may be estimated a-posteriori. Yet, these would be required at the moment of integrating the vehicle's Newton-Euler equations of motion if implemented on-the-fly.

Furthermore, following our assumption that the small celestial body 101 is in single-axis rotation, we set the spin state vector w₀w(t). Along with the unknown gravity parameter μ_a, we may view w₀ as a parameter of the dynamical system.

For the purposes of the RelDyn factor modeling, we may now define the state tuple x(Q, q, v, r)∈SO(3)׳׳׳X, the input tuple u(R, s, f)∈SO(3)׳׳׳U, and the parameter tuple p(w₀,μ_a)∈³×_>0P. To be consistent with earlier representations, we note that by virtue of Fact 4, given any Q, q, we may construct the corresponding pose T=T(Q, q).

Recall that the RelDyn factor function encodes the odometric constraint by means of a likelihood or prior distribution function predicated on the state and the known accumulated observations at discrete times. Therefore, we proceed as follows: (1) first, we identify and quantify the probability distribution which relates the dispersion of the state evolving over time. We wish to be able to evaluate this distribution at a sequence of a-priori unknown discrete times along the trajectory. We conserve two moments of this distribution, as is classically done in Gaussian process procedures; (2) we propagate the first moment (the mean state) by using the non-linear equations of motion, through an on-manifold integration technique; This provides us a means to evaluate the residual of the factor by comparing the propagated state against the next state; (3) we propagate the process covariance by linearizing the model dynamics around specific points, using a piece-wise constant assumption, and then integrate the linear model in the Lie algebra of the manifold. This provides us with a discrete process covariance matrix, which we use as a weighting in the factor formulation.

3.6. Stochastic Differential Equations of Motion

We assume that the dispersion of the state at any given time is due to the accumulated effect of exogenous perturbations on the trajectory across time. Specifically, we wish to exploit the stochastic form of the equations of motion to derive the distribution of trajectories as a function of time. Furthermore, to be able to evaluate the distribution at any desired time, we require the stochastic differential equations to be in continuous form. Finally, we desire to obtain the relevant stochastic differential equations directly from the equations of motion, as further detailed below, by admitting realistic and physical perturbations through the input channels.

To this end, let the 3-dimensional white noise Gaussian process ν_*(t) (where *=R, s, f), such that [ν_*(t)]=0 and [ν_*(t)ν_*^T(τ)]=W_*δ(t−τ), where W_* is the spectral density matrix of the continuous-time noise ν_*. Additionally, we assume that ν_* for all *=R, s, f are mutually uncorrelated.

Consider the equations of motion given in (31)-(34) and substitute in u, where for û({circumflex over (R)},ŝ,{circumflex over (f)})∈U, we have

u=({circumflex over (R)} exp([ν_R]{circumflex over ( )}),ŝ+ν_s,{circumflex over (f)}+ν_f).

To establish the governing stochastic differential equations of the system, we first define the 3-dimensional Wiener process ε_*(t), (*=R, s, f), such that

ε_*(t)=∫₀^tν_*(τ)dτ,ε_*(0)=0,  (35)

with the increment dε_*(t)=ε_*(t+dt)−ε_*(t)=ν_*(t) dt, satisfying [ε_*(t)]=0, [(ε_*(t)−ε_*(τ))(ε_*(t)−ε_*(τ))^T]=W_*|t−τ|, which when τ→t yields the relationship [dε_*dε_*^T]=W_*dt. Note also that dε_*dt=0.

Given any two states x₂=(Q₂, q₂, v₂, r₂)∈X and x₁=(Q₁, q₁, v₁, r₁)∈X, we define the error Δ_X(x₁,x₂)∈T_x1X between x₁ and x₂ centered at x₁ such that:

$\begin{array}{cc}{\Delta }_X\left(x_1,x_2\right)\stackrel{△}{=}\left[\begin{array}{c}{\log \left(Q_1^{\top }{Q}_2\right)}^{\bigvee }\\ q_2-q_1\\ v_2-v_1\\ r_2-r_1\end{array}\right].& \left(36\right)\end{array}$

To obtain the stochastic differential equations resulting from this noise injection in the input, we evaluate the stochastic increment dx(t)=Δ_X(x(t), x(t+dt)), we use Facts 2, 3 and 4 and we separate the equations, while dropping the time dependence for readability, yielding

$\begin{array}{cc}d\kappa ={\left(Q^{\top }\mathrm{dQ}\right)}^{\bigvee }=\left(\hat{s}-Q^{\top }{w}_0\right)\mathrm{dt}+d{\epsilon }_s,& \left(37\right)\end{array}$ $\begin{array}{cc}\mathrm{dr}=\left(Q{\hat{R}}^{\top }v-\left[w_0\right]^Q{\hat{R}}^{\top }r\right)\mathrm{dt}+\left(Q\left[{\hat{R}}^{\top }v\right]^-\left[w_0\right]^Q\left[{\hat{R}}^{\top }r\right]^\right)d{\epsilon }_R& \left(38\right)\end{array}$ $\begin{array}{cc}\mathrm{dv}=\left(-\left(\frac{{\mu }_a}{{r}^3}+\frac{{\mu }_{\odot }}{{d}^3}\right)r+g\left(d\right)+R\hat{f}\right)\mathrm{dt}-\hat{R}\left[\hat{f}\right]^d{\epsilon }_R+\mathrm{Rd}{\epsilon }_f,& \left(39\right)\end{array}$ $\begin{array}{cc}\mathrm{dr}=\mathrm{vdt}.& \left(40\right)\end{array}$

By decomposition of the covariances W_*=L_*L_*⁻¹, (*=R, s, f), and by defining the standard Wiener process ε, where [dε(t)]=0 and [dε(t)dε(t)^T]=I₉dt, obtain the system of stochastic differential equations (41):

dκ=f _Q(x,û,p)dt+[0_3×3L_Qs0_3×3]dε,

dq=f _q(x,û,p)dt+[L_qR0_3×30_3×3]dε,

dv=f _v(x,û,p)dt+[L_vR0_3×3L_vf]dε,

dr=f _r(x)dt,  (41)

• • where x=(Q, q, r, v)∈X, û=(R,ŝ,{circumflex over (f)})∈U, and p=(w₀,μ_a)∈P, and
  • where:
  • L_Qs=I₃, L_vR=−{circumflex over (R)}[{circumflex over (f)}]{circumflex over ( )},L_vf={circumflex over (R)},
  • L_qR=Q[{circumflex over (R)}^Tv]{circumflex over ( )}−[w₀]{circumflex over ( )}Q[{circumflex over (R)}^Tr]{circumflex over ( )}.Collecting into dx=[dκ^T dq^T dv^T dr^T]^T, and into f(x, û, p)=[f_Q^T(x,û,p) f_q^T(x,û,p) f_v^T(x,û,p) f_r^T(x)]^T, we can rewrite (41) in succinct form as:

$\begin{array}{cc}\mathrm{dx}=f\left(x,\hat{u},p\right)\mathrm{dt}+L\left(x,\hat{u},p\right)d\epsilon .& \left(42\right)\end{array}$ $\mathrm{where}$ $\begin{array}{cc}L=\left[\begin{array}{ccc}{0}_{3\times 3}& L_{\mathrm{Qs}}& 0_{3\times 3}\\ L_{\mathrm{qR}}& 0_{33}& 0_{3\times 3}\\ L_{\mathrm{vR}}& 0_{3\times 3}& L_{\mathrm{vf}}\\ 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}\end{array}\right].& \left(43\right)\end{array}$

Due to the non-linearity of the considered dynamical system in equation (42), the distribution of the Itô process x(t) may require more than two moments to fully describe at any time. However, as a stepping stone to obtain a practical representation of the distribution p(x(t)|x(t₀)), we are interested in computing the first two moments {circumflex over (x)}(t)=({circumflex over (Q)}(t),{circumflex over (q)}(t),{circumflex over (v)}(t),{circumflex over (r)}(t))[x(t)] and Σ(t)[Δ_X({circumflex over (x)}(t),x)Δ_X^T({circumflex over (x)}(t),x)] at each time t.

Naturally, we want the relative dynamics factor to be proportional to the distribution p(x_k+1|x_k), as shown in the beginning of the Section 3.3. Noting that the time step between t k and t_k+1 is not known a-priori during the navigation segment, it is crucial to compute a process noise covariance P_k=P(t_k), which appropriately scales with the time step length, using a discretization scheme. Such a treatment is similar to the prediction step of an Extended Kalman Filter with first-order approximation assumptions.

First, for the continuous case, we define the RelDyn factor residual ϵ^RelDyn such that:

ϵ^RelDyn(x(t),x(t+dt))dx(t)−f(x(t),û(t),p)dt.

Then, the factor ϕ^RelDyn, as a function proportional to (0, dP(t)), where:

dP(t)=[(dx(t)−f(x(t),ū(t),p))(⋅)^T],

may be simply defined as:

ϕ^RelDyn(x(t),x(t+dt))exp((ϵ^RelDyn(x(t),x(t+dt)))^TdP⁻¹(t)ϵ^RelDyn(*)).

Let the partition {t_k}_k=0^N of the time interval [t₀, t_f] be such that t₀<t₁< . . . <t_k< . . . <t_N=t_f. Given some û(t), t∈[t₀, t_f], we generate the sequence of predictions {{circumflex over (x)}_k}_k=0^N such that:

$\begin{array}{cc}{\Delta }_X\left({\hat{x}}_{k+1},{\hat{x}}_k\right)={\int }_{t_k}^{t_{k+1}}\hat{f}\left(\hat{x}\left(\tau \right),\hat{u}\left(\tau \right),p\right)d\tau ,\hat{x}\left(t_k\right)={\hat{x}}_k,& \left(44\right)\end{array}$

is a valid discretization of the first-order approximation {circumflex over ({dot over (x)})}(t)≈{circumflex over (f)}({circumflex over (x)}(t),û(t),p) for all k=1, . . . , N−1, and such that:

$\begin{array}{cc}{P}_k={\int }_{t_k}^{t_{k+1}}L\left(\hat{x}\left(\tau \right),\hat{u}\left(\tau \right),p\right)L^{\top }\left(\hat{x}\left(\tau \right),\hat{u}\left(\tau \right),p\right)\mathrm{dt},x\left(t_k\right)=x_k,& \left(45\right)\end{array}$

characterizes the discrete process noise over the time-span [t_k, t_k+1]. Then, using the definition in equation (36) and the equation (44), we define the RelDyn residual ϵ_k^RelDyn for the discrete case, as:

$\begin{array}{cc}{ϵ}_k^{\mathrm{RelDyn}}\left(x_k,x_{k+1}\right)\stackrel{△}{=}\left[\begin{array}{c}{\log \left({\left(Q_k^{\top }{Q}_{k+1}\right)}^{\top }{\Omega }_k^{\mathrm{prop}}\right)}^{\bigvee }\\ q_{k+1}-q_k-{\int }_{t_k}^{t_{k+1}}{f}_q\left(x\left(\tau \right),\hat{u}\left(\tau \right),p\right)d\tau \\ v_{k+1}-v_k-{\int }_{t_k}^{t_{k+1}}{f}_v\left(x\left(\tau \right),\hat{u}\left(\tau \right),p\right)d\tau \\ r_{k+1}-r_k-{\int }_{t_k}^{t_{k+1}}{f}_r\left(x\left(\tau \right),p\right)d\tau \end{array}\right],& \left(46\right)\end{array}$

where x(t_k)=(Q_k, q_k, v_k, r_k) and where Ω_k^prop is defined such that, for a partition

$t_k=s_0<s_1<\dots <s_n=t_{k+1},{\Omega }_k^{\mathrm{prop}}\stackrel{△}{=}\underset{n\to \infty }{\lim }{\prod }_{i=0}^{n-1}\exp \left({\int }_{s_i}^{s_{i+1}}\left[f_Q\left(x\left(\tau \right),\hat{u}\left(\tau \right),p\right)\right]^d\tau \right),$

x(t_k)=x_k, and is dubbed the McKean-Gangolli injection method, allowing for the stochastic process on the Lie group SO(3) to be written as a product integral. With appropriate correction factors, the limit can be truncated to numerically perform on-manifold integration with minimal error accumulation stemming from the approximation, as discussed in Section 4.6 below.

The factor ϕ_k^RelDyn for the discrete case may now be written as:

ϕ_k^RelDyn(x_k,x_k+1)((ϵ_k^RelDyn(x_k,x_k+1))^TP_k⁻¹ϵ_k^RelDyn(x_k,x_k+1)).

3.7. Smoothability of Chosen State Variables

It has been shown in that batch-style maximum a-posteriori solution to the factorization of a joint distribution over the variables to be estimated, which was discussed in Section 3, is equivalent to the solution obtained from optimal fixed-interval smoothing. Typically, the optimal smoothed state is defined as a linear combination of the state of the forward filter and the state of the backward filter at each time, with optimally chosen weights. A state is said to be smoothable if an optimal smoother provides a state estimate superior to that obtained when the final optimal filter estimate is extrapolated backwards in time. Furthermore, it has been shown that, given linear forward and backward filters, only those states which are controllable by the noise driving the system state vector are smoothable.

In this context, we evaluate the relevance of our choice of state variables and the noise input and derived stochastic equations of motion in Section 3.6 in terms of smoothability. However, since the concept of smoothability pertains to linear smoothing, direct application of the criterion to our problem is not possible. Instead, by demonstrating small-time local controllability through the noise by analysing a local linearization of the equations (31) to (33) at some ({circumflex over (x)}_k,û_k), we can assess the smoothability of our chosen state variables at that point.

Specifically, we can compute the local linearization matrices:

$\begin{array}{cc}{F_k\stackrel{△}{=}\frac{\partial f}{\partial x}❘}_{\begin{array}{c}x={\hat{x}}_k\\ u={\hat{u}}_k\\ p=\hat{p}\end{array}}={\left[\begin{array}{cccc}\frac{\partial f_Q}{\partial \kappa }& 0& 0& 0\\ \frac{\partial f_q}{\partial \kappa }& 0& \frac{\partial f_q}{\partial v}& \frac{\partial f_q}{\partial r}\\ 0& 0& 0& \frac{\partial f_v}{\partial r}\\ 0& 0& \frac{\partial f_r}{\partial v}& 0\end{array}\right]}_{\begin{array}{c}x={\hat{x}}_k\\ u={\hat{u}}_k\\ p=\hat{p}\end{array}},\mathrm{where}\frac{\partial f_Q}{\partial \kappa }=\left[Q^{\top }{w}_0\right]^,\left[\kappa \right]^\in T_Q\mathrm{SO}\left(3\right),\frac{\partial f_q}{\partial \kappa }=-Q\left[R^{\top }v\right]^+\left[w_0\right]^Q\left[R^{\top }r\right]^\frac{\partial f_q}{\partial v}={\mathrm{QR}}^{\top },\frac{\partial f_q}{\partial r}=-\left[w_0\right]^{\mathrm{QR}}^{\top },\frac{\partial f_r}{\partial v}=I_3,\frac{\partial f_v}{\partial r}=\frac{{\mu }_a}{{r}^5}\left(3{\mathrm{rr}}^{\top }-{r}^2{I}_3\right)-\frac{{\mu }_{\odot }}{{d}^3}{I}_3& \left(47\right)\end{array}$

and L_k=L({circumflex over (x)}_k,û_k,p), where L is as defined in (43).

We now apply the PBH controllability test to show that our system does have small-time local smoothability at all points of interest in the state space X.

For simplifying purposes, since both {circumflex over (Q)} and {circumflex over (R)} are full rank, we assume {circumflex over (Q)}={circumflex over (R)}=I₃ without loss of generality. Specifically, we see that the matrix

$\frac{\partial f_Q}{\partial \kappa }=\left[w_0\right]^$

is skew-symmetric and is therefore of rank 2. It follows that λ=0 is an eigenvalue of F_k. We compute the matrix [F_k−λI₃,L_k], which for λ=0 takes the form

$\left[F_k,L_k\right]=\left[\begin{array}{ccccccc}\left[{\hat{w}}_0\right]^& 0& 0& 0& 0& I_3& 0\\ \left[\hat{v}\right]^+\left[\hat{w}\right]^\left[\hat{r}\right]^& 0& I_3& \left[\hat{w}\right]^& \left[\hat{v}\right]^-\left[\hat{w}\right]^\left[\hat{r}\right]^& 0& 0\\ 0& 0& 0& I_3& -\left[\hat{f}\right]^& 0& I_3\\ 0& 0& I_3& 0& 0& 0& 0\end{array}\right].$

We observe by inspection that, as long as {circumflex over (v)} and ŵ₀ are not null at the same time, i.e., that the spacecraft 102 has some non-null relative velocity in the small celestial body fixed frame, then the PBH controllability matrix has rank 12, and thus the system is small-time locally controllable through the noise.

Therefore, it is crucial to include appropriate and physically justified noise terms in the considered dynamics to support smoothability. Furthermore, the inclusion of such noise avoids the degeneracy of the discrete process noise Gaussian distribution, with covariance P_k, as computed in Section 3.3. This argument supports the injection of noise through the input channel in the derivation of the stochastic equations of motion in Section 3.3.

4. Algorithm and Implementation Details

We are now ready to finally state the full solution to the problem and provide the key implementation details. Equations (44) and (45) constitute constraints on the evolution of the mean system state and its associated system noise covariance between time instances t_k and t_k+1. We wish to enforce these constraints. Practically, at each time instance t_k, k=1 . . . N, we let {tilde over (x)}_k be the current best guess of {circumflex over (x)}_k and {tilde over (P)}_k be the current best guess of P_k. Given {tilde over (x)}_k, {tilde over (x)}_k+1 and û(t), t∈[t_k, t_k+1], we compute the residual due to equation (44), then the guess matrix {tilde over (P)}_k, and finally the value of the factor ϕ_k^RelDyn({tilde over (x)}_k,{tilde over (x)}_k+1). Note that the integrals in equations (46) and (45) are computed numerically using an on-manifold integration scheme, detailed later in this section. We note that the inclusion of the RelDyn factor introduces the spin state w₀ and the gravity parameter μ_a as variables into the factor graph.

As both the prior and projection factors are developed in detail in other works (see, for example), we simply restate them here using our notation, such that:

ϕ^prior(x₀)exp(Δ_X({circumflex over (x)}₀,x₀)^TΣ₀⁻¹Δ_X({circumflex over (x)}₀,x₀)),

ϕ_*^proj(x_k,_i)exp((y_ik^m−y_ik)^T(Σ_y^m)⁻¹(y_ik^m−y_ik)),

where each y_ik^m∈_k is obtained as per the relationship in equation (6). We incorporate these prior factors on the initial state tuple. We finally incorporate the relative kinematics factors introducing the center-of-mass position parameter vector c.

The structure of the problem may now be organically schematized and analyzed as a factor graph. We have illustrated two initial steps of the AstroSLAM factor graph 400 in FIG. 4, where the RelKin factors and their related dependencies are illustrated by the nodes 401 and 403 and their edges, and the RelDyn factors and their related dependencies are illustrated by the nodes 405 and 407.

Next, we provide the implementation details of the algorithm, while underlining the specificities of the asteroid navigation problem. We first discuss the issue of complexity growth in the graph induced by the inclusion of global parameters. Then we discuss the state and map initialization steps of the algorithm, in which we insert prior factors for initial poses and kinematic variables of both the spacecraft and small celestial body, as well as generate an initial estimate for the map Ψ₁.

Finally, we discuss the overall architecture of AstroSLAM. Our processing pipeline consists of a front-end system and a back-end system predicated on the iSAM2 engine and GTSAM library, along with an initialization step and a loop closure detection step.

5.1. Considerations Regarding Complexity Growth

Since we want to estimate the navigational states and simultaneously update the center of mass position c, the spin vector w₀, and the gravity parameter μ_a estimates on the fly, we have introduced these as variables in the factor graph, by means of the relative kinematic factor and RelDyn factor formulations. It is clear, however, that the inserted edges render the graph fully connected, since each state tuple x_i, i=1, . . . , n is linked to the single parameter μ_a by means of the RelDyn factor ϕ_i−1^RelDyn and the factor ϕ_i+1^RelDyn. Due to this, the clique sizes in the Bayes tree determined at each solution step of iSAM2 would get larger as the graph grows incrementally, and the updates in the Bayes tree would involve cliques that are ever more deep in the tree. If not addressed, this issue essentially forces iSAM2 to perform expensive batch optimization at every time-step, with the size of the batch update growing at every time-step. One intuitive workaround, is to solve for the variables related to the inertial factors in an incremental fashion at every step, but only sporadically optimize the kinematic rotation factors using a batch solve step, so to avoid the repeated and expensive batch optimization. However, we cannot implement such a method for our procedure, as our relative odometry RelDyn factor also depends on the global parameter μ_a.

In our approach, we somewhat emulate the procedure of a fixed-lag smoother, by defining epochs i=1, . . . , m to which we associate specific global variables c⁽ⁱ⁾, w₀⁽ⁱ⁾ and μ_a⁽ⁱ⁾. Over the epoch, the iSAM2 optimizer performs a batch update of the variable only related to that epoch at each timestep, by virtue of the fully connected property of the graph induced by the relative kinematic factor and the RelDyn factors. At every new epoch i, we marginalize the global parameter variables c^(i-1), w₀^(i-1), and μ_a^(i-1) from the previous epoch, and insert a new prior factor for each of the current global parameter variables c⁽ⁱ⁾, w₀⁽ⁱ⁾, and μ_a⁽ⁱ⁾. When inserting relative kinematic factors or RelDyn factors during epoch i, we associate these factors to their corresponding epoch's global parameter variables c_i and μ_a⁽ⁱ⁾. A new epoch is induced at a fixed timestep interval (e.g., every 10 steps).

4.2. State Initialization

We leverage estimates from pre-encounter Earth-based measurements and approach phase sensor measurements to perform the initialization of both the spacecraft and small celestial body states. Inertial position measurements of the spacecraft 102, modelled as +ν_r, where ν_r˜(0_3×1,Σ_r^m), are based on Earth-relative radiometric ranging and bearing measurements, a method of localization widely practiced in deep space mission spacecraft tracking using communication station networks, such as the DSN. Accurate ground-based navigation estimates using Earth-relative range and bearing measurements, such as uplink-downlink pulse ranging and Delta-Differenced One-Way Ranging (DDOR) are available to be used for initialization. It is important to note that our algorithm requires DSN-type measurements only during the initialization phase, to anchor the initial pose as described below.

Inertial position estimates of the small celestial body, denoted by , are predicated on orbit determination (OD) performed using Earth-based telescopic measurements. An initial relative position vector {circumflex over (r)}_S0G estimate may be obtained by differencing spacecraft and small celestial body inertial position estimates. Alternatively, relative optical navigation (OpNAV) performed using pre-encounter sensor measurements, accurate in the order of several hundreds of meters for large small-bodies and tens of meters for small small-bodies, may be used to provide such an initial prior on the relative position.

To further simplify the initialization step, we assume that the small celestial body 101 is a stable single axis rotator. Thus, we align the arbitrary frame such that =w₀{right arrow over (g)}₃, where w₀ is the magnitude of the angular velocity vector of the small celestial body 101. The unit vector is classically parameterized by the spin pole tilt angles relative to the J2000 ecliptic plane. If sufficient pre-encounter observations are available, a prior for these angles can be estimated and their associated uncertainty computed.

Typically, a prominent and salient feature on the surface is hand-picked as the prime meridian direction, thus fixing the {right arrow over (g)}₁ axis. In our case, given the initial prior discussed earlier, the prime meridian {right arrow over (g)}₁ may be initialized by first computing {right arrow over (g)}₂ following =[]{circumflex over ( )}/∥∥, and then computing =[]{circumflex over ( )}. Consequently, R₀=[] is determined. Note that, for this method to work, the spacecraft relative position at initialization time cannot be coincident with the small celestial body spin axis. In practice, an on-board star tracker system is leveraged to obtain an orientation measurement , modelled as exp([ν_R]{circumflex over ( )}),ν_R˜(0_3×1,Σ_R^m). Orientation measurements are usually known with very good accuracy and little uncertainty.

We may now establish the prior factor for relative pose by combining the information we have available at initialization. We first define the relative measurement at time k=0, using (4), Oven by

$\begin{array}{c}=\left[\begin{array}{cc}& \\ 0_{1\times 3}& 1\end{array}\right]=\left[\begin{array}{cc}\exp \left(\left[v_R\right]^\right)& +v_r\\ 0_{1\times 3}& 1\end{array}\right]\\ =\left[\begin{array}{cc}& r_{S_0{G}_0}^{{𝒢}_0}\\ 0_{1\times 3}& 1\end{array}\right]\left[\begin{array}{cc}\exp \left(\left[v_R\right]^\right)& v_r\\ 0_{1\times 3}& 1\end{array}\right]\end{array},$

where is the relative position vector at time k=0 obtained by differencing inertial position estimates of the spacecraft 102 and the small celestial body 101, or by means of an OpNav solution. It then follows that:

$\begin{array}{cc}=\left[\begin{array}{cc}\exp \left(\left[{\epsilon }_1\right]^\right)& {\epsilon }_2\\ 0_{1\times 3}& 1\end{array}\right],& \left(48\right)\end{array}$

where ε₁=ν_R and ε₂=ν_r. Using the logarithm map of the SE(3) group of homogeneous transformations at the identity and the vee operator, as detailed in Section 1.1, we can write:

$\begin{array}{cc}{\left[\log \left(\right)\right]}^{\bigvee }=\left[\begin{array}{c}{\epsilon }_1\\ {\epsilon }_2\end{array}\right]~\left(0_{6\times 1},{\sum }_{T,0}^m\right),& \left(49\right)\end{array}$

where, by first-order linear approximation, we have Σ_T,0^m=J_RΣ_R^mJ_R^T+J_r()Σ_r^mJ_r^T(), computed using Fact 4 where R∈SO(3). A prior factor ϕ₀^prior() is emplaced in the graph. This factor encodes the residual between the pose ∈SE(3) and the measurement , with covariance Σ_T,0^m.

4.3. Map Initialization

Initialization of the map is delayed until time index k=1, at which point at least two images of the target with sufficient parallax are captured. Local image features γ₀ and γ₁ are extracted and undergo data association, with outlier rejection, producing a set of 2D-2D correspondences. A strict outlier rejection criterion is used to obtain a subset of reliable correspondences. It is now possible to apply a typical 8-point algorithm using the inlier 2D-2D correspondences to find a guess pose transformation {circumflex over (T)}°∈SE(3) such that:

$\hat{T}°\stackrel{△}{=}\left[\begin{array}{cc}\hat{R}°& \hat{t}°\\ 0& 1\end{array}\right],$

where {circumflex over (R)}°∈SO(3) and {circumflex over (t)}°∈³, describe the estimated change in relative orientation and relative position between the poses of the camera frame at time indices k=0 and k=1. From the frame kinematics and composition rules, and knowing that {circumflex over (R)}°=, it follows that =. We note that ≈exp([−Δt₀]{circumflex over ( )}) It follows that

$\approx \exp (\Delta t_0]^)\hat{R}°.$

Additionally, we know that:

$\begin{array}{cc}t°==\lambda \left(+-\left(+\right)\right)=\lambda \left(\exp \left(\left[\Delta t_0\right]^\right)-\right),& \left(50\right)\end{array}$

with λ>0, an unknown scaling factor.

Thus, in order to establish reliable relative pose estimates {circumflex over (T)}₀ and {circumflex over (T)}₁, a good prior knowledge of the small celestial body 101 inertial angular velocity is required. Additionally, the ambiguity in the scale parameter λ still remains. Typically, scale can be established using off-nadir altimeter measurements over a range of orbital configurations in a multi-arc solution, as perform by the OSIRIS-REx mission around asteroid Bennu. However, the latter solution, in its current format, is not amenable to on-the-fly autonomy, as the batch estimation process requires a compilation of data from multiple mission arcs. Furthermore, altimeter measurements have to be used at close range to the target small celestial body 101, again restricting autonomy since judiciously pre-designed maneuvers have to be executed to first obtain the necessary altimeter measurements at close range. If a good initial estimate of the angular velocity w₀ of the small celestial body is available, the map scale can be determined by combining DSN-like inertial position measurements with the orientation matrix and position vectors extracted from the relative poses T₀ and T₁, as shown in Equation (50). The availability of such prior knowledge currently depends heavily on Earth-based and pre-encounter measurements of the target small celestial body 101. Regarding this matter, the ambiguity in scale due to the parameter λ is determining in the estimation of an appropriate spherical gravity parameter term μ_a. An error in μ_a will then affect the propagation step of the smoothing process. For the purpose of this work, we have assumed that the scale parameter λ is well known.

Having an estimate of poses {circumflex over (T)}₀ and {circumflex over (T)}₁, we compute the triangulation of the landmarks Ψ₁ using measurements γ₀ and γ₁. We then generate guess values for the estimated landmark positions {}_L∈Ψ1. At this step, all appropriate factors are inserted based on the 2D-to-3D correspondences, resulting in the factor graphs 500 and 501 illustrated in FIG. 5, where the shorthands _i, L_i∈Ψ_k and T_k, k=0, . . . , n are used for brevity.

4.4. Front-End

The front-end system includes feature detection and matching, and encodes the structure of a typical monocular bundle adjustment pipeline in a factor graph based on index mappings defined in Section 1.3.

At every new image acquisition, first, we store the corresponding wall clock time, as the difference between two consecutive image acquisition times is necessary for propagating the RelDyn factors inserted later in the factor graph. We then detect ORB features in the new image. ORB features perform well in practice and they are fast to compute. These are also good placeholders for more robust automatic features to be implemented in the future. We use brute force nearest neighbor search based on the Hamming distance for binary feature descriptors to perform feature matching against previous image frames (e.g., n−1, n−2) and populate an initial 2D-to-2D correspondence map. We also enforce an essential matrix constraint as a geometric check for matched features per a RANSAC procedure, thus eliminating outlier associations in the 2D-to-2D correspondence map. We append the corresponding matches to a list of feature tracks maintained since the beginning of the algorithm procedure.

All tracked feature points in the current frame n may be separated into three subsets: a set of sightings for new landmark points to be triangulated, a set of sightings of known landmark points to be tracked, and a set of outlier points to be rejected. The outlier set of feature points is discarded.

The first subset of the tracked features corresponds to measurements of surface point landmarks that have not been previously sighted and which need to be triangulated using all accumulated sightings in the feature track. Leveraging the latest guess pose , a guess value for the 3D position of the landmarks is generated using triangulation. We delay triangulation and insertion of new landmarks into the graph based on the number of sightings of the associated landmark. If that number is above a predetermined threshold, say 3, then we attempt triangulation of the point using all of the measurements in the sighting of that surface point. If the landmark point has been successfully triangulated using the feature point sightings that were used for the triangulation, appropriate 2D-to-3D correspondences capturing the new data association are accumulated. Projection factors between the new landmark variable and all camera poses where the landmark was sighted are added to the factor graph.

The second subset of tracked features corresponds to new sightings of surface point landmarks previously triangulated. This subset is populated by verifying if the currently tracked feature has an associated landmark already in the 2D-to-3D correspondences maps. The feature-landmark matched pairs undergo a reprojection error test, after which the surviving pairs establish new GTSAM projection factors inserted into the graph between known landmarks and relating to the most recent frame n. If there are enough tracked features from time index k=n−1, then visual tracking is successful and the PnP algorithm is used to guess the camera pose value from matched correspondences. Alternatively, when RelDyn factors are included and if the PnP solution is of poor quality, we use the prediction from the motion model at time t_n to guess a new camera pose .

For every new frame inserted, the graph is incrementally augmented with the new variables and factors by the front-end, as is illustrated in the graphs 600 and 601 of FIG. 6.

4.5. Back-End

The underlying structure of the navigation problem is thus captured by encoding visual SLAM measurement constraints, relative kinematic constraints, Earth-relative inertial position measurement constraints as factors, and RelDyn factors in a single factor graph, as we explained in Section 3.1.

In practice, the optimization associated with the small celestial body relative navigation problem is performed incrementally in the factor graph framework using the iSAM2 algorithm. The back-end implements iSAM2 algorithm's procedure, which was discussed earlier in Section 3.1. iSAM2 evaluates the losses at each factor marked for update, computes the associated Jacobians at the guess values, and performs the minimization for inference, yielding a new estimate solution. Note that since the process noise covariance related to the RelDyn factor, as detailed in Section 3.6, is not fixed, but is rather the solution to a dynamical equation, the out-of-the-box GTSAM factor template had to be modified to allow for a variable process noise covariance to be incorporated accordingly. This procedure deviates from the pre-integrated IMU accelerometer factor, since our case, the integrated quantities are non-linear in terms of the state. Indeed, at each evaluation of the factor error, the discrete process noise covariance, also a non-linear function of the state, is propagated accordingly. This covariance is associated with a Gaussian process distribution, and once computed, we replace GTSAM's internal representation of the covariance of the factor error with new process noise covariance.

4.6. Computation of the RelDyn Factor Residual

We wish to compute numerically integrated quantities between successive states, so to evaluate the RelDyn factor error established in Section 3.3. To this end, we perform on-manifold discrete Crouch-Grossman geometric integration to obtain accurate predictions which also respect the constraints of the relative rotation matrix Q(t)∈SO(3). Given {t_k}_k=0^N, for every k=0, . . . , N−1, let Δt_kt_k+1−t_k, we use an N_s-stage Crouch-Grossman geometric integration scheme and compute the propagated quantities:

$\begin{array}{cc}{\overline{\Omega }}_k^{\mathrm{prop}}\stackrel{\Delta }{=}{\prod }_{i=1}^{N_s}\exp \left({\left[b_i\Delta t_k{\tilde{f}}_Q^{\left(i\right)}\right]}^{\wedge }\right),& \left(51\right)\end{array}$ $\begin{array}{cc}\left[\begin{array}{c}\Delta {\tilde{q}}_k^{prop}\\ \Delta {\tilde{v}}_k^{prop}\\ \Delta {\tilde{r}}_k^{prop}\end{array}\right]\stackrel{\Delta }{=}{\sum }_{i=1}^{N_s}{b}_i\Delta t_k\left[\begin{array}{c}{\tilde{f}}_{q,k}^{\left(i\right)}\\ {\tilde{f}}_{v,k}^{\left(i\right)}\\ {\tilde{f}}_{r,k}^{\left(i\right)}\end{array}\right],& \left(52\right)\end{array}$

whereby the shorthands:

$\begin{array}{cc}\widetilde{f_{*,k}^{\left(i\right)}}\stackrel{\Delta }{=}f*\left({\tilde{x}}_k^{\left(i\right)},u_k^{\left(i\right)},p\right),*=Q,q,v,r& \left(53\right)\end{array}$ $\begin{array}{cc}{\tilde{x}}_k^{\left(i\right)}\stackrel{\Delta }{=}\left({\tilde{Q}}_k^{\left(i\right)},{\tilde{q}}_k^{\left(i\right)},{\tilde{v}}_k^{\left(i\right)},{\tilde{r}}_k^{\left(i\right)}\right),& \left(54\right)\end{array}$ $\begin{array}{cc}{û}_k^{\left(i\right)}\stackrel{\Delta }{=}\hat{u}\left(t_k+c_i\Delta t_k\right),& \left(55\right)\end{array}$ $\begin{array}{cc}{\tilde{Q}}_k^{\left(i\right)}\stackrel{\Delta }{=}{\tilde{Q}}_k{\prod }_{j=1}^{i-1}\exp \left({\left[a_{ij}\Delta t_k{\tilde{f}}_{Q,k}^{\left(j\right)}\right]}^{\wedge }\right),& \left(56\right)\end{array}$ $\begin{array}{cc}\left[\begin{array}{c}{\tilde{q}}_k^{\left(i\right)}\\ {\tilde{v}}_k^{\left(i\right)}\\ {\tilde{r}}_k^{\left(i\right)}\end{array}\right]\stackrel{\Delta }{=}\left[\begin{array}{c}{\tilde{q}}_k\\ {\tilde{v}}_k\\ {\tilde{r}}_k\end{array}\right]+{\sum }_{j=1}^{i-1}{a}_{ij}\Delta t_k\left[\begin{array}{c}{\tilde{f}}_{q,k}^{\left(j\right)}\\ {\tilde{f}}_{v,k}^{\left(j\right)}\\ {\tilde{f}}_{r,k}^{\left(j\right)}\end{array}\right]& \left(57\right)\end{array}$

are defined and computed for every i=1, . . . , N_s. The coefficients a_ij, b_i, c_i, j<i=1, . . . , N_s are obtained from an appropriate Butcher table. The number of rows in the table represents the number of time stages, with the {c_i} column representing the coefficients of steps in time t_i=t_k+c_iΔt_k where are functions are to be evaluated, the number of columns is representing the number of half-steps, with {c_i} representing the weight of each half-step, and the coefficients {a_ij} representing the weight of the intermediary stage i when computing the half-step j.

The RelDyn residual ϵ_k^RelDyn({circumflex over (x)}_k,{circumflex over (x)}_k+1) is now computed such that:

$\begin{array}{cc}{ϵ}_k^{\mathrm{RelDyn}}\left({\tilde{x}}_k,{\tilde{x}}_{k+1}\right)\stackrel{\Delta }{=}\left[\begin{array}{c}\log {\left({\left({\tilde{Q}}_k^T{\tilde{Q}}_{k+1}\right)}^{\top }{\overline{\Omega }}_k^{prop}\right)}^{\bigvee }\\ {\tilde{q}}_{k+1}-{\tilde{q}}_k-\Delta {\tilde{q}}_k^{prop}\\ {\tilde{v}}_{k+1}-{\tilde{v}}_k-\Delta {\tilde{v}}_k^{prop}\\ {\tilde{r}}_{k+1}-{\tilde{r}}_k-\Delta {\tilde{r}}_k^{prop}\end{array}\right].& \left(58\right)\end{array}$

4.7 Process Noise Propagation

To compute P_k we introduce the matrix κ_k such that, given any pair ({circumflex over (x)}_k,û_k).

${\Delta }_X\left(x\left(t_{k+1}\right),{\hat{x}}_k\right)=\Phi \left(t_{k+1},t_k\right){\Delta }_X\left(x\left(t_k\right),{\hat{x}}_k\right)+{\int }_{t_k}^{t_{k+1}}\Phi \left(t_{k+1},\tau \right)L\left(\hat{x}\left(\tau \right),p\right)d\epsilon \left(\tau \right)=\Phi \left(t_{k+1},t_k\right){\Delta }_X\left(x\left(t_k\right),{\hat{x}}_k\right)+{\Lambda }_k{\epsilon }_k,$

with initial conditions {circumflex over (x)}(t_k)={circumflex over (x)}_k, and û(t_k)=û_k, and where ε_k˜(0, I₉) and {dot over (Φ)}(t, t_k)=F({circumflex over (x)}(t),û(t), p)Φ(t, t_k), Φ(t_k, t_k)=I₉. We see that P_k=Λ_kΛ_k^T. It follows that:

$P_k={\int }_{t_k}^{t_{k+1}}\Phi \left(t_{k+1},\tau \right)L\left(\hat{x}\left(\tau \right),p\right)L^{\top }\left(\hat{x}\left(\tau \right),p\right){\Phi }^{\top }\left(t_{k+1},\tau \right)d\tau ,\hat{x}\left(t_k\right)={\hat{x}}_k,\hat{u}\left(t_k\right)={\hat{u}}_k.$

We further approximate the linearized system by assuming that the matrices F, L shown in Section 4.7 are piecewise constant over the interval [t_k, t_k+1], i.e., F_kF({circumflex over (x)}(t_k),û(t_k), p), L_kL({circumflex over (x)}(t_k),û(t_k), p). Then, we construct the following matrix:

$=\left[\begin{array}{cc}{F}_k& L_k{L}_k^{\top }\\ 0_{12\times 12}& F_k^{\top }\end{array}\right],$

take the matrix exponential, yielding:

${{=}_{}}^{}\left[\begin{array}{cc}{\Phi }_k& {\Gamma }_k\\ 0_{12\times 12}& {\Phi }_k^{\top }\end{array}\right],$

and we extract from the right-hand side the desired submatrices, i.e., the discrete state transition matrix Φ_k=(t_k+1, t_k) of the piecewise-constant linearization, as well as the covariance of the discretized process noise Λ_kε_k, given as P_k=Γ_kΦ_k^T.

5.8 Loop Closure

To perform loop closure, we leverage a bag-of-words (BoW) representation. We convert each image i=1, . . . , N to a bag-of-words vector v_i and compute the similarity metric

$\begin{array}{cc}s\left(v_i,v_j\right)\stackrel{\Delta }{=}1-\frac{1}{2}❘\frac{v_i}{❘v_i❘}-\frac{v_j}{❘v_j❘}❘.& \left(59\right)\end{array}$

We compare all prior images that are at least 10 frames away from the current frame. When the similarity score for two images is greater than a threshold we perform an additional geometric check and then add a factor between the poses corresponding to the detected loop.

FIG. 7 is an illustration of an operational flow of an implementation of a method 700 for performing path planning for a spacecraft in the presence of a small celestial body. The method 700 may be implemented by the navigation engine 120 of the spacecraft 102.

At 710, a plurality of images are received. The plurality of images may be received by the navigation engine 120 of the spacecraft 102 from a sensor 110 such as a camera. Each image of the plurality of images may be associated with a time. Each image may be of a small celestial body 101 near the spacecraft 102. The small celestial body 101 may a comet or an asteroid.

At 720, a plurality of spacecraft and small celestial body states are initialized. The states may be initialized by a front-end component of the navigation engine 120. The states may be initialized using prior information known about the small celestial body 101 and/or the spacecraft 102.

At 730, a 3D map is generated from the plurality of images. The 3D map be generated by the front-end component of the navigation engine 120 from the plurality of images. The 3D map may be a 3D map of the small celestial body 101 and may be include a plurality of landmarks determined from the images of the plurality of images.

At 740, a factor graph is generated. The factor graph may be generated by a back-end component of the navigation engine 120. The may be generated based on the landmarks in the 3D map, the initialized states of the spacecraft 102, and one or more RelDyn factors that link the relative pose of the spacecraft 102 to the problem of predicting trajectories stemming from the motion of the spacecraft 102 in the vicinity of the small celestial body 101. Example factor graphs are shown in the graphs 500 and 501 of FIG. 5.

At 750, a new image is received. A new image may be received by the back-end component of the navigation engine 120. The new image may be received from the sensor 110 and may be an image that was not part of the plurality of images.

At 760, new landmarks or changed location are determined based on the new image. The new landmarks or changed location may be determined by the back-end component of the navigation engine 120.

At 770, the factor graph is updated. The factor graph may be updated by the back-end of the navigation engine 120 to reflect the new landmarks or changed locations in the new image.

At 780, the navigation problem is solved. The navigation problem may be solved by the back-end of the navigation engine 120 performing an optimization on the factor graph. The results of the optimization may include estimates of past relative positions of the positions of the small celestial body 101 and the spacecraft 102, past relative orientations of the small celestial body 101 and the spacecraft 102, and past relative velocities of the small celestial body 101 and the spacecraft 102.

At 790, path planning is performed based on the output of the navigation problem. The path planning may be performed by the navigation engine 120 for the spacecraft 102. The path planning may include control of the spacecraft 102 in the vicinity of the small celestial body 101.

FIG. 8 shows an exemplary computing environment in which example embodiments and aspects may be implemented. The computing device environment is only one example of a suitable computing environment and is not intended to suggest any limitation as to the scope of use or functionality.

Numerous other general purpose or special purpose computing devices environments or configurations may be used. Examples of well-known computing devices, environments, and/or configurations that may be suitable for use include, but are not limited to, personal computers, server computers, handheld or laptop devices, multiprocessor systems, microprocessor-based systems, network personal computers (PCs), minicomputers, mainframe computers, embedded systems, distributed computing environments that include any of the above systems or devices, and the like.

Computer-executable instructions, such as program modules, being executed by a computer may be used. Generally, program modules include routines, programs, objects, components, data structures, etc. that perform particular tasks or implement particular abstract data types. Distributed computing environments may be used where tasks are performed by remote processing devices that are linked through a communications network or other data transmission medium. In a distributed computing environment, program modules and other data may be located in both local and remote computer storage media including memory storage devices.

With reference to FIG. 8, an exemplary system for implementing aspects described herein includes a computing device, such as computing device 800. In its most basic configuration, computing device 800 typically includes at least one processing unit 802 and memory 804. Depending on the exact configuration and type of computing device, memory 804 may be volatile (such as random access memory (RAM)), non-volatile (such as read-only memory (ROM), flash memory, etc.), or some combination of the two. This most basic configuration is illustrated in FIG. 8 by dashed line 806.

Computing device 800 may have additional features/functionality. For example, computing device 800 may include additional storage (removable and/or non-removable) including, but not limited to, magnetic or optical disks or tape. Such additional storage is illustrated in FIG. 8 by removable storage 808 and non-removable storage 810.

Computing device 800 typically includes a variety of computer readable media. Computer readable media can be any available media that can be accessed by the device 400 and includes both volatile and non-volatile media, removable and non-removable media.

Computer storage media include volatile and non-volatile, and removable and non-removable media implemented in any method or technology for storage of information such as computer readable instructions, data structures, program modules or other data. Memory 804, removable storage 808, and non-removable storage 810 are all examples of computer storage media. Computer storage media include, but are not limited to, RAM, ROM, electrically erasable program read-only memory (EEPROM), flash memory or other memory technology, CD-ROM, digital versatile disks (DVD) or other optical storage, magnetic cassettes, magnetic tape, magnetic disk storage or other magnetic storage devices, or any other medium which can be used to store the desired information and which can be accessed by computing device 800. Any such computer storage media may be part of computing device 800.

Computing device 800 may contain communication connection(s) 812 that allow the device to communicate with other devices. Computing device 800 may also have input device(s) 814 such as a keyboard, mouse, pen, voice input device, touch input device, etc. Output device(s) 816 such as a display, speakers, printer, etc. may also be included. All these devices are well known in the art and need not be discussed at length here.

The term “about,” as used herein, means approximately, in the region of, roughly, or around. When the term “about” is used in conjunction with a numerical range, it modifies that range by extending the boundaries above and below the numerical values set forth. In general, the term “about” is used herein to modify a numerical value above and below the stated value by a variance of 10%. In one aspect, the term “about” means plus or minus 10% of the numerical value of the number with which it is being used.

Similarly, numerical ranges recited herein by endpoints include subranges subsumed within that range (e.g., 1 to 5 includes 1-1.5, 1.5-2, 2-2.75, 2.75-3, 3-3.90, 3.90-4, 4-4.24, 4.24-5, 2-5, 3-5, 1-4, and 2-4). It is also to be understood that all numbers and fractions thereof are presumed to be modified by the term “about.”

The following patents, applications, and publications, as listed below and throughout this document, describes various application and systems that could be used in combination the exemplary system and are hereby incorporated by reference in their entirety herein.

Claims

1. A computer-implemented comprising: receiving a plurality of images of a small celestial body by a spacecraft, wherein each image of the plurality of images is associated with a time;initializing a plurality of spacecraft states and small celestial body states by the spacecraft;based on the received plurality of images, generating a 3D map of the small celestial body including a plurality of landmarks by the spacecraft, wherein each landmark is associated with a location in the 3D map;generating a factor graph corresponding to the plurality of landmarks and the plurality of spacecraft states and small celestial body states by the spacecraft;receiving a new image by the spacecraft;based on the new image and some or all of the plurality of images, identifying one or more new landmarks in the 3D map and a changed location for at least one of the plurality of landmarks by the spacecraft;updating the factor graph based on the one or more new landmarks and the changed location by the spacecraft; andsolving a navigation problem based on the factor graph by the spacecraft to generate relative orientations of the spacecraft and the small celestial body, relative positions of the spacecraft and the small celestial body, and relative velocities of the relative positions of the spacecraft and the small celestial body;performing path planning for the spacecraft by the spacecraft using some or all of the relative orientations of the spacecraft and the small celestial body, the relative positions of the spacecraft and the small celestial body, the relative velocities of the relative positions of the spacecraft and the small celestial bod, and the 3D map; andcausing the spacecraft to travel based on the path planning by the spacecraft.

2. The computer-implemented method of claim 1, further comprising: receiving visual odometry constraints and relative kinematic and dynamical constraints; andsolving the navigation problem based on the factor graph and the received visual odometric constraints and the relative kinematic and dynamical constraints.

3. The computer-implemented method of claim 2, wherein the dynamical constraints comprise the RelDyn factors.

4. The computer-implemented of claim 1, further comprising receiving prior information about the small celestial object and initializing the plurality of small celestial body states based on the prior information.

5. The computer-implemented method of claim 4, wherein the prior information comprises a spin state estimate.

6. The computer-implemented method of claim 1, further comprising solving the navigation problem based on the factor graph by the spacecraft to generate an estimate of a center-of-mass position of the small celestial body.

7. The computer-implemented method of claim 1, wherein the small celestial body comprises an asteroid or a comet.

8. A spacecraft comprising: at least one computing device; anda computer-readable medium with computer-executable instructions stored thereon that when executed by the at least one computing device cause the at least one computing device to:receive a plurality of images of a small celestial body, wherein each image of the plurality of images is associated with a time;initialize a plurality of spacecraft states and small celestial body states;based on the received plurality of images, generate a 3D map of the small celestial body including a plurality of landmarks, wherein each landmark is associated with a location in the 3D map;generate a factor graph corresponding to the plurality of landmarks and the plurality of spacecraft states and small celestial body states;receive a new image;based on the new image and some or all of the plurality of images, identify one or more new landmarks in the 3D map and a changed location for at least one of the plurality of landmarks;update the factor graph based on the one or more new landmarks and the changed location; andsolve a navigation problem based on the factor graph to generate relative orientations of the spacecraft and the small celestial body, relative positions of the spacecraft and the small celestial body, and relative velocities of the relative positions of the spacecraft and the small celestial body;perform path planning for the spacecraft using some or all of the relative orientations of the spacecraft and the small celestial body, the relative positions of the spacecraft and the small celestial body, the relative velocities of the relative positions of the spacecraft and the small celestial bod, and the 3D map; andcause the spacecraft to travel based on the path planning.

9. The spacecraft of claim 8, further comprising computer-executable instructions that when executed by the at least one computing device cause the at least one computing device to: receive visual odometry constraints and relative kinematic and dynamical constraints; andsolve the navigation problem based on the factor graph and the received visual odometric constraints and the relative kinematic and dynamical constraints.

10. The spacecraft of claim 9, wherein the dynamical constraints comprise the RelDyn factors.

11. The spacecraft of claim 8, further comprising receiving prior information about the small celestial object and initializing the plurality of small celestial body states based on the prior information.

12. The spacecraft of claim 11, wherein the prior information comprises a spin state estimate.

13. The spacecraft of claim 8, further comprising solving the navigation problem based on the factor graph by the spacecraft to generate an estimate of a center-of-mass position of the small celestial body.

14. The spacecraft of claim 8, wherein the small celestial body comprises an asteroid or a comet.

15. A non-transitory computer-readable medium with computer-executable instructions stored thereon that when executed by the at least one computing device cause the at least one computing device to: receive a plurality of images of a small celestial body, wherein each image of the plurality of images is associated with a time;initialize a plurality of spacecraft states of a spacecraft and small celestial body states;based on the received plurality of images, generate a 3D map of the small celestial body including a plurality of landmarks, wherein each landmark is associated with a location in the 3D map;generate a factor graph corresponding to the plurality of landmarks and the plurality of spacecraft states and small celestial body states;receive a new image;based on the new image and some or all of the plurality of images, identify one or more new landmarks in the 3D map and a changed location for at least one of the plurality of landmarks;update the factor graph based on the one or more new landmarks and the changed location; andsolve a navigation problem based on the factor graph to generate relative orientations of the spacecraft and the small celestial body, relative positions of the spacecraft and the small celestial body, and relative velocities of the relative positions of the spacecraft and the small celestial body;perform path planning for the spacecraft using some or all of the relative orientations of the spacecraft and the small celestial body, the relative positions of the spacecraft and the small celestial body, the relative velocities of the relative positions of the spacecraft and the small celestial bod, and the 3D map; andcause the spacecraft to travel based on the path planning.

16. The computer-readable medium of claim 15, further comprising computer-executable instructions that when executed by the at least one computing device cause the at least one computing device to: receive visual odometry constraints and relative kinematic and dynamical constraints; andsolve the navigation problem based on the factor graph and the received visual odometric constraints and the relative kinematic and dynamical constraints.

17. The computer-readable medium of claim 16, wherein the dynamical constraints comprise the RelDyn factors.

18. The computer-readable medium of claim 15, further comprising receiving prior information about the small celestial object and initializing the plurality of small celestial body states based on the prior information.

19. The computer-readable medium of claim 18, wherein the prior information comprises a spin state estimate.

20. The computer-readable medium of claim 15, further comprising solving the navigation problem based on the factor graph by the spacecraft to generate an estimate of a center-of-mass position of the small celestial body.