- Email: [email protected]
Space collision probability computation based on on-board optical cues
Accepted Manuscript Space collision probability computation based on on-board optical cues Meng Yu, Shuang Li, Shu Leng PII:
S0094-5765(18)30336-9
DOI:
https://doi.org/10.1016/j.actaastro.2018.06.053
Reference:
AA 6972
To appear in:
Acta Astronautica
Received Date: 14 February 2018 Revised Date:
8 June 2018
Accepted Date: 27 June 2018
Please cite this article as: M. Yu, S. Li, S. Leng, Space collision probability computation based on onboard optical cues, Acta Astronautica (2018), doi: https://doi.org/10.1016/j.actaastro.2018.06.053. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Space Collision Probability Computation Based on On-board Optical Cues Meng Yu1, 2* Shuang Li1, 2† and Shu Leng1 1. College of Astronautics, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China
210016, China
RI PT
2. Advanced Space Technology Laboratory, Nanjing University of Aeronautics and Astronautics, Nanjing
Abstract: This paper presents a novel on-board space collision analysis method for space situational
SC
awareness. The framework is developed under the following assumptions: 1) A satellite can be equipped
M AN U
with on-board sensors for space object recognition. 2) No a–priori knowledge of the space objects is provided. A space object size and relative state estimation method is firstly proposed, wherein optical cues acquired from onboard sensors are utilized to achieve the estimation. Then, the unscented transform approach is employed to calculate the probability density function (PDF) of collision probability based on
TE D
the estimate information. Monte Carlo simulations and an experimental test demonstrate that the proposed approach can achieve high-precision on-board collision probability estimation with an error less than 3%.
1. Introduction
EP
Keywords: Optical cue; Collision analysis; Probability density function; Gaussian uncertainty.
AC C
The problem of vast growing population of resident space objects (RSO) poses a great challenge to long-term safe space operation [1][2]. As the near-earth space environment is becoming more crowded than ever, fast and accurate collision analysis has attracted much attention in the field of space situational awareness. The collision between space objects depends on many attributes, for instance, size, shape, orbital movement, self-rotation of space objects [3]. Considering the lack of a-priori precise knowledge *
Lecture, Department of Astronautics Engineering, Nanjing University of Aeronautics and Astronautics.
†
Professor, Department of Astronautics Engineering, Corresponding Author, Email: [email protected], Tel:
+86(25)84896039. -1-
ACCEPTED MANUSCRIPT about the motion state of space objects, these aforementioned attributes are often associated with uncertainties, which leads to the current probabilistic modeling treatment that generally existed in space collision analysis methods [4].
RI PT
The probabilistic modeling manner aims at evaluating the probability of collision between two space objects. Such a method commonly relies on many assumptions, for instance, a RSO is usually assumed to be influenced by planetary gravitations and J2 perturbation effects. At the same time, its shape is also
SC
approximated as a sphere or an ellipsoid with a determined size in order to simplify the derivation [8]. In
M AN U
addition, many works that evolve computing the collision probability assumed a limited encounter time between two objects, which is also known as short-term encounters [10], wherein the relative positions and velocities are assumed to admit linear forms with Gaussian-distributed uncertainties. For complicated scenarios such as a long-term encounter or RSOs with possible maneuvers, the computation of collision
TE D
probability might be more legitimate when taking into account the non-Gaussian and nonlinear distributed uncertainties in propagating the relative position and velocity [8]. In general, the collision analysis of space objects can be classified into two categories. The first approach
EP
relies on Monte Carlo (MC) sampling [14]. Such a MC based method samples from the probability density
AC C
functions (PDFs) of two RSOs, the relative motion between two RSOs is estimated from these samples and propagated over a certain epoch of interest to record the collision events. The MC sampling methods can be applied to both simple and complicated encounter scenarios with a high accuracy, However, such a technique requires significant computational resources, thus it is currently computationally intractable for a computer system onboard the space vehicle. For lessening the computational burden, Graphics Processing Unit (GPU) can be used to parallelize the computation and then improve real-time performance [14]. The second method focuses on statistically inferring the propagation of RSOs’ relative states over a certain time
-2-
ACCEPTED MANUSCRIPT range, wherein the uncertainty associated with RSOs’ states are usually described by Gaussian or non-Gaussian probability density functions (PDF) [9][16]. The collision probability can then be calculated by integrating the probability density function of relative states over a specific time range and examining
RI PT
the intersection between the relative states and hardball regions ( a hardball region is usually defined as a sphere, whose radius is the sum of actual radii of two RSOs [8]). Specifically, the collision probability can be computed via either determining the volume swept out by the hardball [18] or solving the amount of
SC
influx flow of the PDF of the relative states into the hardball [1][21]. These methods share a similar
M AN U
structure, such that they all capture the relative state uncertainty in a probabilistic manner and use efficient integration methods to calculate the collision probability [9].
To date, most of the researches in the field of space collision probability analysis are based on ground detection and tracking. Some recent works in the field of space situational awareness suggest that efforts
TE D
should be made to enhance the quick response to potential space collision events [23], and space operations might benefit from employing on-board sensors to detect RSOs and conduct collision analysis. Motivated by such demand, this paper proposes a novel approach for on-board space collision analysis using on-board
EP
optical sensors. The proposed approach employs a monocular camera for space object detection and a laser
AC C
detection and ranging device (e.g. LADAR) for measuring the relative distance. Then, the acquired optical data is incorporated to derive the PDF of the relative state between two RSOs. After that, an unscented transform based Gaussian mixture model (GMM) is adopted to deduce the PDF of the collision probability. We follow the common assumption that the RSO can be modeled by an ellipsoid [10][13][14]. Unlike the priori work in the related field of research, this work does not assume a-priori knowledge of the RSO’s size information. Alternatively, we propose to utilize the optical cues acquired from on-board sensors to estimate such information.
-3-
ACCEPTED MANUSCRIPT The rest of this paper is organized as follows. Section 2 gives the brief review, coordinates and problem formulation. Section 3 details the incorporation of camera and LADAR ranging in hardball radius determination and relative state estimation. Section 4 introduces a framework of optical sensor based space
RI PT
collision analysis. Simulation results and discussion in terms of practical concerns of the proposed approach are presented in Section 5. Finally, Section 6 discusses the conclusions of this paper.
2. Coordinate and problem formulation
SC
We first introduce the coordinate frames used in this study. A space-encounter scenario is presented in
M AN U
Fig. 1 as an illustration, denote by S a satellite that is well-tracked by ground stations is operating in the low-earth orbit (LEO), and assume a space object that might have a short or long term encounter with S is detected by optical sensors onboard the satellite S . The sensors consist of a monocular camera and a LADAR, which are attached to one side of the satellite with calibrated installation matrices. The inertial
TE D
frame is denoted by {I} , the satellite body frame is denoted by {S} , the camera and LADAR body frame are denoted by {C},{L} , respectively, and the 2D image frame of the camera is denoted by {im} ,
zS
zI
xS
AC C
EP
which is perpendicular to the y − axis of the satellite body frame.
yS
OC OL
zd
yd
Od
O im
xd
yI
xI
Fig. 1 Illustration of inertial {I}, satellite {S}, object {d}, camera {C}, LADAR {L} frame
-4-
ACCEPTED MANUSCRIPT The position and velocity of space object and satellite in the inertial frame are denoted by
( rdI , r&dI ) , ( rSI , r&SI ) ,
respectively. Denote by x RI = [ rRI , r&RI ] the relative state between them, it can be T
defined as follows:
rRI = rdI − rSI r&RI = r&dI − r&SI
CS and LADAR LC S , the relative state can be expressed
as: T
x = [ r , r& S
L R
]
T
= ( C S )( C I ) x L
S
(2)
I R
M AN U
where
L R
SC
x RC = [ rRC , r&RC ] = ( C C S )( S C I ) x RI L R
RI PT
C
Given the installation matrix of camera
(1)
C I is the matrix that denotes the rotation from the inertial frame {I} to the satellite’s body frame
{S}, as is shown in Fig. 1, the space object in the 2D image frame {im} is represented by a shaded ellipse, whose center is located at pdim = [ µd , vd ]T in {im} , which can be obtained by the pin-hole
TE D
transformation [24] as follows:
µd =
fim pRC ( x) + µ0 pRC ( z )
(3)
EP
f pC ( y) vd = im C R + v0 pR ( z )
where fim represents the focal length of the camera,
[ µ0 , v0 ]
T
AC C
intersects the image frame, rRC = [ pRC ( x ), pRC ( y ), pRC ( z ) ]
T
is the coordinate at which the optical axis
denotes the relative position expressed in the
camera frame as state in Eq. (2). The on-board optical measurement at time step t can be expressed as:
Z d (t) = {[ µd (t ), vd (t ) ] , Ld (t ) }
(4)
where Ld (t ) is the LADAR range measurement, which is a scalar. Given the definition of coordinate frames and measurement model, the problem of collision between the satellite and the space object can be formulated as follows. A space object d is determined as a collision threat to satellite if the relative position rR (t ) falls -5-
ACCEPTED MANUSCRIPT within the hardball region Λ H at time t > t0 , where the hardball radius H Sd is defined as the sum of the radii of satellite and space object. The collision probability Pc (t ) is a function of time t , hardball region
Λ H and relative state xRI that is defined in Eq. (2), Pc (t ) can be calculated by integrating the PDF
t
Pc (t ) = ∫ pc (λ )d λ
(5)
−∞
where pc (t ) is a probability density function that is defined as:
∫
A vR ≤ 0
v R p ( x RI ( t ; t0 ) ) dvt dA
(6)
SC
pc ( t ) = ∫
RI PT
pc (t ) over time t [9]:
where vR = vt ⋅ rˆ , A denotes the surface area of the hardball H dS , rˆ denotes the normal vector with
M AN U
respect to the area element dA , vt = r&dI − r&SI denotes the relative velocity, and vR is a scalar that denotes the amount of relative velocity projected onto the normal vector rˆ . The integration over vR ≤ 0 indicates the influx of probability distribution for the relative state into the hardball region, which can be
TE D
interpreted as the boundary of the hardball that is crossed by the relative velocity. vt . As illustrated in Fig. 2, the collision probability can be calculated by integrating over the crossed surface of hardball.
rˆ
dA
rR
EP
vt
θ
AC C
φ
hardball
Fig. 2 Illustration of influx to the hardball region
By employing the spherical integration rule, Eq. (5) can be expressed as [19]:
Pc ( t ) =
t0 + T 2 π
π 2
∫ ∫ ∫ ( H Sd )
t0
0 − π2
2
cos θ
∫
vt ⋅rˆ ≤ 0
vt ⋅ rˆ
∞
∫ρ
S
( xS , t ; t0 ) ρ d ( xd , t ; t0 ) dxS dv R dθ dφ dt (7)
−∞
where
θ , φ represent the spherical angles, H Sd is the radius of the hardball, T is the total integration
time,
ρ S , ρd are the PDF for the state of the satellite and the space object, respectively. The process of
-6-
ACCEPTED MANUSCRIPT spherical integration is illustrated in Fig. 3:
dA dθ
rdθ r cos θd φ
φ
RI PT
θ
r = H Sd
dA = r2 cos θ dθ dφ
dφ
Fig. 3 Illustration of the spherical integration used in Eq. (7)
As shown in Fig. 3, wherein the area integral can be denoted by dA = (H S ) cos θ dθ dφ . Note that in Eq. (7), the differential term for
SC
d 2
ρS and ρ d are dx S and dxd , which can be reorganized as:
M AN U
dxS dxd = dx S d ( xS + x R ) = dxS dx R = dxS drR dv R = ( H Sd ) cos θ dxS dv R drdθ dφ = dAdr 2
(8)
dr dAdt = vt ⋅ rˆ dAdt dt
=
By such manner, Eq. (8) can be readily related to Eq. (6). The computation of Pc requires the
TE D
knowledge of H Sd , ρ S , ρ d , the determination of which will be detailed in Section 3.
3. Optical Cue based metric determination
AC C
collision analysis.
EP
In this section, we propose an approach to determine the key metrics of RSO as illustrated in Fig. 4 for
de br is
z c (L )
Ladar ranging
α max
{C }({L})
φL
ϕL
tr av
e ers
z on
e θd
hardball
αL
a ccess ible zone
yc (L )
xc ( L ) Fig. 4 Space object size determination based on optical cues
As shown in Fig. 4. We assume that a space object can be precisely detected by object detection
-7-
ACCEPTED MANUSCRIPT methods (related works on object detection can be found in [25]), and the ranging direction of on-board LADAR can be controlled via sending the azimuth/elevation command:
ϕcom = arccos
(9)
f im
RI PT
φcom =
f im2 + µ d2 lˆ
f + µ d2 2 im
where pim = [ µd , vd ]T denotes the central coordinate of the detected object in the image frame;
µd
vd ]T represents a three-dimensional vector from the satellite to the space object.
M AN U
and lˆ = [ f im
SC
ϕcom , φcom are the elevation and azimuth angle control command of the on-board LADAR, respectively;
Denote by Lˆ d the LADAR bearing angle, note that sometimes the continuous bearing adjustment of LADAR might be too computationally intensive, thus some deviations between the LADAR bearing angle
Lˆ d and lˆ should be allowed. We first define the deviation angle between Lˆ d and lˆ as α% L , and
TE D
define an accessible zone (see Fig. 4), such that only the LADAR measurement that falls within the accessible zone is considered valid for distance ranging. The angle of field of view (FoV) for the accessible zone can be determined by
rim | lˆ |
EP
α max = arctan
(10)
AC C
As illustrated in Error! Reference source not found., denote by lin the location at which the laser
ˆ in = cin − lin can be beam hits the space object; and cin the center of the space object, a vector p generated to estimate the radius of the space object Rd . Denote by
θ d the angle between pˆ in and Lˆ d ,
it can be related to α% L by noting that:
θd ≈
π α% L ⋅ 2 α max
(11)
Defining the radius of the space object as Rd . We can estimate Rd based upon the following
-8-
ACCEPTED MANUSCRIPT geometrical relationship:
rim = Rd
| lˆ |
(12)
Rd2 + Ld + 2 Rd Ld cos θ d 2
which yields 2
Solving for the root of Eq. (13), and noting that | lˆ |>> rim yields:
rim | Ld | | lˆ |2 −rim2 sin 2 θ d + rim2 | Ld | cos θ d | lˆ |2 −r 2 im
(14)
SC
Rd =
(13)
RI PT
( rim2 − | lˆ |2 ) Rd2 + 2rim2 | Ld | cos θ Rd + rim2 Ld = 0
Then, the hardball radius can be obtained by H Sd = Rd + RS , where RS is the radius of the satellite.
Rd2 + Ld + 2 Rd Ld cos θ d as the distance from the space object to the satellite, the
M AN U
Defining rSD =
2
relative position in the camera frame can be expressed as
rRc =
lˆ ⋅ rSD | lˆ |
(15)
TE D
Given rRc , the relative velocity in the camera frame can be approximated by solving the differential of position:
EP
r&Rc (t ) = lim
rRc ( t + ∆t ) − rRc ( t )
∆t →0
∆t
(16)
AC C
4. Computation of collision probability The computation of collision probability is carried out given the estimation of the relative state
x Rc = [ rRc , r&Rc ] and the hardball radius H Sd . Such an estimate would inevitably contain measurement T
noises, which results in a non-linear uncertainty. We employ the unscented transform (UT) [28] to cope with this problem due to its robustness against non-linear noise. First, define the joint PDF of the states of the satellite and the space object as
p ( xS , xd ) = p ( x S , x S + x R ) -9-
(17)
ACCEPTED MANUSCRIPT The PDF of the relative state x R can be expressed by marginalizing over xS , yielding:
p ( x R ) = ∫ p ( x S , x S + x R ) d xS
(18)
To facilitate the computation of collision probability, the UT-transform [28] is employed to
variable, (e.g. x R , x S or xd ), such that:
mi = x +
(
mi + n = x −
(n + κ ) Pxx
(
)
(n + κ ) Pxx
W0 = κ / (n + κ )
Wi = 1 / 2(n + κ ), i = 1,..., n
i
)
i
Wi + n = 1 / 2(n + κ ), i = 1,..., n
SC
m0 = x
RI PT
characterize the PDF. Specifically, 2n + 1 sigma points are generated based on a n − dimensional
(19)
M AN U
Where x and Pxx denote the mean and covariance of a variable, respectively; W0~2 n denotes
the weights for the generated sigma points m0~2 n . Using the UT-transform, PDFs for the states of the satellite and the space object can be expressed as: US
p ( xS ) = ∑ Wi S pg ( xS ; miS , PxxS ) j =1
Ud
(20)
TE D
p ( xd ) = ∑ Wi pg ( xd ; m , P ) b
d j
d xx
i =1
where U S = U d = 2n + 1 . pg ( x, m, P ) denotes a standard Gaussian PDF. Assuming the satellite state
EP
and object state are independent, we can break their joint PDF into a product of two PDFs, such that:
p( x S , xd ) = p( x S ) p( xd ) U s Ud
AC C
= ∑∑ Wi SWi d pg ( x S ; miS , PxSx ) pg ( xd ; m dj , Pxxd )
(21)
i =1 j =1
The authors in [21] proved that
∫ p (x g
S
; m S , PxxS ) pg ( xS + x R ; md , Pxxd ) dx S = pg ( x R ; uE , Σ E )
(22)
where uE = mS − md , Σ E = PxxS + Pxxd . Substituting Eq.(21), (22) into Eq. (18), the PDF for the relative state can be rewritten as U s Ud
p ( x R ) = ∑∑ Wi SW jd pg ( x R ; uE(ij ) , PE(ij ) ) i =1 j =1
- 10 -
(23)
ACCEPTED MANUSCRIPT Following the method proposed in [21], the PDF for relative state is decomposed into the product of PDFs for the relative position and velocity to maintain the computational stability. Specifically, uE , PE are decomposed into:
u
Σ(rvij ) Σ(vijv )
(24)
RI PT
Σ(rrij ) ur(ij ) ( ij ) = (ij ) , PE = T ( Σ(rijv ) ) uv
( ij ) E
Substituting Eq. (24) into Eq. (23), and following the decomposing rule proposed in [15], it shows that
SC
p( xR ) can be expressed as U s Ud
p ( x R ) = ∑∑ Wi SW jd pg (rR ; ur(ij ) , Σr(ij ) ) pg ( v '; ( uv(ij ) ) ', ( Σv(ij ) ) ' ) i =1 j =1
M AN U
where x R = [rR , v R ]T , and
(25)
v ' = v R − ( Σ(rvij ) ) ( Σ(rrij ) ) rR −1
T
( uv(ij ) ) ' = Σ(vvij ) ( Σ(rrij ) )
−1
ur(ij )
( Σ(vij ) ) ' = Σ(vvij ) − Σ(vrij ) ( Σ(rrij ) )
−1
(26)
Σ(rijv )
TE D
Substituting Eq. (25) into Eq. (6) and reallocate the differential elements, such that the integration is taken over the hardball region, the collision probability density function can be expressed as US Ud
pc ( x R ) = ∑∑ Wi SW jd ∫ pg (rR ; ur(ij ) , Σ (rij ) )Θ ( rR , v R ) dA
(27)
A
EP
i =1 j =1
AC C
where Θ ( rR , v R ) is the integration over the relative velocity that is expressed as:
Θ ( rR , v R ) = ∫ ( ij )
vn ≤ 0
vn( ij ) pg ( v '; ( uv( ij ) ) ', ( Σ (vij ) ) ') dv '
(
where vn( ij ) is a scalar and admits the form vn( ij ) = v '⋅ rˆ + rˆT Σ (rvij )
(28)
) ( Σ( ) ) T
ij rr
−1
rR . Given the hardball
radius H Sd estimated by the approach proposed in Section 3, and using the spherical transformation (see Eq. (7)), pc ( x R ) can be expressed as
pc ( xR ) = ( H Sd )
2
US Ud
∑∑W
i
i =1 j =1
S
W jd ∫
2π
0
∫
π 2
− π2
pg ( H Sd rˆ; ur(ij ) , Σ(rij ) )Ξ ( rˆ, uE , PE ) cos θ dθ dφ (29)
As is defined in [21], Ξ ( rˆ, uE , PE ) is expressed as - 11 -
ACCEPTED MANUSCRIPT Ξ ( rˆ, uE , PE ) =
Ξ 02 ( rˆ ) Ξ 0 ( rˆ ) σ ( rˆ ) Ξ 0 ( rˆ ) exp − 2 1 − erf − 2 2π 2σ ( rˆ ) 2σ ( rˆ )
(30)
where
σ
2
( rˆ ) = rˆ
T
(
Σ( ij ) vv
−(
Σ (ij ) rv
)( T
Σ (ij ) rr
)
−1
Σ( ij ) rv
) rˆ
(31)
RI PT
T −1 Ξ 0 ( rˆ ) = rˆ T uv( ij ) + ( Σ(rvij ) ) ( Σ(rrij ) ) ( H Sd rˆ − ur(ij ) )
Thus the calculation of Θ ( rR , v R ) can be conducted analytically. The rest of the integrals over the
θ and φ in Eq. (29) corresponds to integrating the relative position over the hardball
SC
spherical element
sphere. The spherical integration over the relative position cannot be performed analytically. Alternatively,
M AN U
it can be given by numerical approximation, for example, using a numerical quadrature integration approach [9]. In this study, we choose the Lebedev and Laikov’s method [29] to approximate the integration over the relative position. Lebedev and Laikov’s method is a kind of quadrature approach that can approximate the surface integral over a sphere by a summation of weighted N quadrature points nˆ k .
pc (t ) = ( H dS )
2
TE D
In this manner, the PDF for collision probability at time t can be expressed as Us Ud
N
i =1 j =1
k =1
∑∑ Wi,StW jd,t ∑ wk pg ( H dS nˆ k ; ur(ij,t ) , Σ(rij,t) ) Ξ ( nˆ k , uE ,t , PE ,t )
(32)
EP
The PDF for collision probability can be solved by substituting Eq. (32) into Eq. (5) assuming that
AC C
Pc (t0 ) = 0 at t = t0 , the probability of collision can be expressed as
Pc (t ) = ∫ pc ( λ )d λ t
0
(33)
Such integration over time can be efficiently approximated using quadrature rules, such as the trapezoidal integral method.
5. Simulation results and discussion 5.1 Numerical simulation To demonstrate the effectiveness of the proposed approach for the on-board computation of the
- 12 -
ACCEPTED MANUSCRIPT probability of collision, three space encounter scenarios are considered. Denote by S a near-earth satellite, and B a space object orbiting nearby the satellite S . The space object B is assumed to possess a regular ellipsoid shape, its size is set to small/medium/large in the 1st/2nd/3rd scenario, respectively, which
velocity setup for satellite in the inertial frame is given as:
SC
−6074.6621 −2438.9919 m / s −258.4168 0.9829 0.1294 0.1305 I QS = 0.1125 −0.9851 0.1294 0.1453 −0.1125 −0.9829
(34)
M AN U
[ pS 0
−4017.3356 v S 0 ] = 6435.5823 km 6207.7969
RI PT
might correspond to a low/medium/high collision probability with the satellite S . The initial position and
where I QS is a matrix that defines the rotation from the satellite frame to the inertial frame. We assumed that the on-board camera has just observed a space object at t = t0 . The satellite is assumed to be
TE D
well-tracked by ground stations, its observation uncertainty is characterized by a single Gaussian model. The covariance P0 = diag ( PS ,r0 , PS ,v0 , PS ,q0 ) of this Gaussian distribution at the initial observation time
t = t0 in the inertial frame is given as
EP
PS ,r0 = diag ([ 0.08,0.08, 0.08]) km2 , PS ,v0 = diag ([1.5,1.5,1.5]) m 2 / s 2 , PS ,q0 = diag ([ 0.1, 0.1, 0.1]) deg 2
(35)
AC C
The covariance is chosen according to the results acquired from the ground tracking of near-earth satellites [30]. The space object is assumed to be observed onboard for the first time at
t = t0 with no
a-priori knowledge of its shape and size, indicating an unknown magnitude of uncertainty in its initial state
estimation. To cope with this issue, the uncertainty of its state is characterized by a Gaussian mixture model acquired from the UT-transform. The initial position and velocity for the space object in the camera frame were given by:
- 13 -
ACCEPTED MANUSCRIPT −1430.9141 rBC0 = 4104.0239 m, −5154.4521
0.0651 r&BC0 = −0.1677 m / s 0.2153
(36)
We set the simulation time to 60 mins and assumed that the trajectories for satellite and space object in
RI PT
all three scenarios were identical. In this manner, the only difference in each scenario is the size and shape of the space object. Consider the standard two-body problem, the time history of the relative distance acquired by propagating the initial means is presented in Fig. 5. 7000
5000
400 350 300 250 200 150
M AN U
relative distance (m)
6000
SC
relative re lati vedistance distance (m) (m)
450
100 50
4000
2440
2460
2480
2500
2520
2540
2560
2580
Simulation time (s)
Simulation time (s)
3000 2000
0 0
TE D
1000
500
1000
1500
2000
2500
3000
3500
4000
Simulation time (s)
Fig. 5 Ground truth of the time history for the relative distance
EP
It is shown from Fig. 5 that the true minimum distance between satellite and the space object is
AC C
44.54m, indicating a close fly-by of the space object with respect to the satellite. We additionally developed a simple space object simulator in Matlab software for simulating the on-board optical cue and set the hardball radius to 1m, 5m, 25m in 1st/2nd/3rd scenarios, respectively. An example of a space object fly-by is given in Fig. 6:
- 14 -
RI PT
ACCEPTED MANUSCRIPT
detection
SC
Fig. 6 Example of a space object close fly-by in synthesized images, Hough transform is used for object
In this example, a space object was assumed to sweep over the target satellite. Fig. 6 presents four
M AN U
synthesized images, which are captured by the on-board camera during the fly-by. [ µ d , vd ] represents the center of the space object, and rim represents the image radius of the space object. The space object is assumed to appear as an ellipse in a synthesized image, classic contour detection algorithms such as Hough
TE D
transform [31] or Canny edge fitting [32] can be utilized to detect the space object. The focal length of camera is set at 0.5 cm, the pixel size is set to 4 um. Denote by
Z B 0 = ( µ0 , v0 )T , L0
T
the initial measurement of the space object, where pim ,0 = [ µ0 , v0 ]
EP
center of the object in the image frame, and
AC C
and the relative state in the camera frame
T
is the
L0 is the LADAR ranging distance. The actual object radius
[ rRc , r&Rc ]
T
can be estimated by following Eqs. (10) to (16),
the
corresponding measurement noise covariance R is taken to be diagonal with image location with a value of
2
0.3 pixel 2 and LADAR ranging measurement with a value of 0.2% L , such that
R0 = diag ([0.3 pix 2 , 0.3 pix 2 , 0.2% L ]) . 2
The time history of estimation error for the object radius and relative state in this example is depicted in Fig. 7:
- 15 -
RI PT
ACCEPTED MANUSCRIPT
SC
Fig. 7 Time history of radius and relative state estimation error
M AN U
It can be seen from Fig. 7 that the estimation error for object radius and relative state remains within 1% during the entire simulation, and the estimation for radius error exhibits an unstable jitter around 0.4%, which is partially due to the non-Gaussian noises added to the simulated image. The estimation for relative position and velocity is relatively stable with an error below 0.2%. Simulations of other examples also
satellite and the space object.
TE D
indicate that the radius estimation error is approximately proportional to the relative position between the
The initial measurement covariance is used to determine the weights, mean, and covariance in
EP
UT-transform. After that, the Lebedev and Laikov’s method [29] is applied to generate the Gaussian
AC C
mixture models in each time step, and then apply Eq. (32) to compute the PDF of collision probability, note that using more integral points from Lebedev and Laikov’s method can improve the overall accuracy, however, at the expense of more computational effort. In this simulation, we choose N k = 194 in Lebedev and Laikov’s method, which is a compromise between accuracy and computational cost. The time history of collision probability and its PDF in three scenarios is presented in Fig. 8, a)~f):
- 16 -
PDF of collision probability in 1st scenario.
c)
PDF of collision probability in 2nd scenario.
e)
PDF of collision probability in 3rd scenario.
b) collision probability in 1st scenario
M AN U
SC
Collision Probability, 2nd scenario
a)
RI PT
ACCEPTED MANUSCRIPT
EP
TE D
d) collision probability in 2nd scenario
f)
collision probability in 3rd scenario
AC C
Fig. 8 Collision probability and its PDF in three scenarios
We adopted the trapezoidal integration method to compute the collision probability. As is shown in Fig. 8 a)~f), the final collision probability for the first/second/third scenario is 0.1207% , 3.5142% and 18.6421% according to Fig. 8 b), d) and e). The PDF for collision probability peaks at around 2500s throughout three scenarios and returns to zero after roughly 3000s, which is consistent with the time history of relative distance (see Fig. 6), indicating an accurate reflection of the trend of the relative state variation. In addition, the PDF of collision probability exhibits another peak at around 3000s in Fig. 8 e), which is
- 17 -
ACCEPTED MANUSCRIPT mostly due to the uncertainty in object detection. As we set the radius of hardball to 25m in this scenario (the radius of satellite is 3m, and the radius of space object is 22m), the space object possessed a large part of the synthesized image, which brings ambiguity to the estimate of the object’s center. The maximum error
RI PT
of hardball radius estimation reaches 3% in this scenario, which is nearly 3 times bigger than that in other two scenarios. The improvement over the radius estimation of a space object with irregular shape or a large size is left to future work.
SC
To verify the collision probability computation, an extensive Monte Carlo simulation is carried out.
M AN U
The space object’s orbit samples are drawn from the initial uniform distributions p ( x S0 ) , p ( x B0 ) in each scenario and propagated over the entire simulation time. For each MC run in each scenario, the orbital elements of satellite and object B are propagated for 60 mins at a time step of 0.5 s. The MC simulation stops if the position between object B and satellite S is smaller than the hardball radius, which is
TE D
regarded as a space collision event. The collision probability is calculated by dividing the total number of Monte Carlo runs by the total number of collision events. For each scenario, 1 million MC runs are carried out, the variation of collision probability along with the number of MC runs for three scenarios are shown
Collison probability, P c
AC C
MC runs.
EP
in Fig. 9 a), b), c), respectively, demonstrating that the collision probability converges after roughly 4 × 105
a) MC based Collision probability for 1st scenario
b) MC based Collision probability for 2nd scenario - 18 -
ACCEPTED MANUSCRIPT 0.25
0.2
0.15
0.1
0 0
2
4
6
8
10
Number of MC Runs for the 3rd scenario
RI PT
0.05
12 10 5
c) MC based Collision probability for 3rd scenario
SC
Fig. 9 Collision probability based on Monte Carlo simulation for three scenarios (1 million trails)
In the first scenario, 1187 collision events were triggered for 1 million MC runs, which yields a
M AN U
collision probability of Pc ,1st = 0.1187% , the computation of collision probability for the 2nd and 3rd scenario is the same as in the first scenario, which are Pc ,2 nd = 3.4233% , Pc ,3rd = 19.1232% . The proposed algorithm differs by 1.65%, 2.65% and 2.51% from Monte Carlo results for the three scenarios,
TE D
respectively, demonstrating a plausible accuracy in collision probability calculation. In addition, since the proposed approach relies on the assumption that space object B can be approximated by a regular ellipsoid, its application in practice might benefit from introducing a more accurate object recognition
EP
algorithm.
AC C
5.2 Experimental Test
To further validate the applicability of the proposed approach, a semi-physical experiment was additionally conducted. The experimental set-up is exhibited in Fig. 10
- 19 -
SC
RI PT
ACCEPTED MANUSCRIPT
M AN U
Fig. 10 Illustration of the complete experimental set-up
As shown in Fig. 10, a Laser Tracker serves as the on-board optical sensor suite, a spherically mounted retro-reflector (SMR) ball is mounted on the last joint of a 6D manipulator to simulate the space object, and the 6D manipulator is employed to approximate the orbital movement of the space object.
TE D
Specifically, we first simulate an orbit and downsize such orbit based upon the parameters of the 6D manipulator (for example, its total length is 1.76m), such that the last joint of this manipulator is
EP
programmed to follow the scaled orbital movement. The laser tracker is set to remain stationary to track the target SMR ball, from which the size information and relative state is computed and then employed in
AC C
collision probability estimation framework. The experiment lasts 20 mins, the data acquisition frequency of the optical sensor is set to 1Hz. Some snapshots of the target SMR ball movement during the experiment are shown in Fig. 11:
Fig. 11 Moments of the target ball in motion during the experiment - 20 -
ACCEPTED MANUSCRIPT Note that the target SMR ball is an accessory of the Laser Tracker, therefore the Laser Tracker can automatically recognize the SMR ball and range the relative distance. The simulated orbit and the location
RI PT
of the optical sensor is shown in Fig. 12. The time history of the relative distance is presented in Fig. 13
simulated orbit Optical Sensor Location
t0 1.5
tend
1 0.5
1
4
SC
Z axis (m)
2
0
3
-1
2
-2
1 -3
0
X axis (m)
M AN U
Y axis (m)
-1
-4
EP
TE D
Fig. 12 Illustration of the simulated orbit and the sensor location
AC C
Fig. 13 Time history of the relative distance (measured by the Laser Tracker)
As shown in Fig. 12 and Fig. 13, the minimum relative distance is 0.1682m, which occurs at 265s. The relative distance data is captured by the Laser Tracker with an accuracy of 0.1um, which serves as the ground truth for validating the accuracy the PDF estimation. We set the radius of the hardball to 0.20m, and follow the procedure established in section 5.1 to calculate the collision probability rate. Specifically, we used the optical cues acquired during the beginning to infer the collision probability and its PDF, results of which are shown in Fig. 14:
- 21 -
ACCEPTED MANUSCRIPT 0.4 0.35 0.3 0.25 0.2 0.15
0.05 0 0
240
RI PT
0.1
480
720
960
1200
Experiment Time (s)
Fig. 14 Experiment test result: PDF of collision probability (left) and collision probability (right)
SC
As depicted in Fig. 14, the time history of PDF of collision probability generally agrees with that of
M AN U
the relative distance in Fig. 13. Note that the PDF of collision probability would in theory peak at the location where the relative distance between two RSOs reaches its minimum. In this experiment, the PDF peaks at 267s, which only differs from the ground truth by 2s, thus showing a plausible estimation accuracy. Unfortunately, since it is impractical to actually let the target ball hit the satellite in this experiment. We
TE D
cannot validate the accuracy of collision probability, a more credible experimental test for validating the accuracy of the collision probability estimate is left for future work. 5.3 Practical Concerns
EP
The aforementioned on-board collision probability computation framework aims at enabling the
AC C
quick response to space threatens, especially the newly spawned objects. To implement this framework in practice, the capability of on-board optical sensors must be primarily considered. Choice of on-board Sensors The space object detection prefers a camera with a long focal length. When dealing with the detection of small and irregular space objects, a camera with long focal length will yield a higher physical resolution than the one with short focal length. For example, suppose the relative distance between a satellite and a space debris varies from 10km to 4km, the actual radius of the space debris and the pixel block is 5~20m
- 22 -
ACCEPTED MANUSCRIPT and 3e-6m，respectively. The relationship between the camera focal length and the image area of the object
M AN U
SC
RI PT
is shown in Fig. 15:
Fig. 15 Relationship between the focal length (m) and object’s image area (Number of pixels)
TE D
In Fig. 15, the actual radius of the space debris is set to 5m, 10m, 15m, 20m, respectively. It shows that the image area of the space debris increases as the focal length is enlarged from 5mm to 0.3m. Taking the squared case (actual debris radius is 5m) highlighted in Fig. 15 for example, and assuming the relative
EP
distance is 10km, an accurate estimate of the debris radius would need at least 300 pixels, which requires the focal length to be at least 0.15m. However, using a camera with long focal length also brings difficulty
AC C
to sensor installment and mass distribution. Moreover, a narrow FOV camera also limits the detection rate of space objects. Thus in practice, a balance between accuracy and feasibility must be properly stroked. Besides camera, the detection range of the onboard LIDAR is also an important factor. To save enough time for making a collision avoidance maneuver decision, a spacecraft needs to detect the space objects as early as possible. Therefore for safety concerns, an onboard LADAR with a long detection range is favored. Recent developments of onboard LIDAR design show that the LIDAR system will be able to detect a 10cm size space object given a relative distance of roughly 100km [33]. A major practical concern of an onboard LIDAR system is its high power consumption, an onboard system might not allow a long range LIDAR sensor to operate in full-time. To cope with this issue, sequential images acquired from the onboard camera
- 23 -
ACCEPTED MANUSCRIPT can be employed to estimate the relative state while the LIDAR is offline. This alternative method has the advantage of much lower power consumption, but it is also less accurate than the LIDAR system. An autonomous sensor switching mechanism can be designed to remedy this issue. For instance, we can wake up the LIDAR ranging system when the estimated relative distance is below a threshold, and actively guide
RI PT
the LIDAR to measure the relative distance to improve the calculating of the collision probability. Warning Time for Collision Avoidance
The collision avoiding maneuver needs to be carried out if the calculated collision probability is high. Such a maneuver requires an enough advanced warning time. To estimate the warning time, one can think
SC
of the epoch at which the collision probability rate peaks as the time where a collision event might happen, thus estimating the warning time can be achieved by finding the peaks of the calculate PDF for collision
M AN U
probability. Fig. 16 shows the peaks of the collision probability rate that corresponds to the case in Fig. 8, e):
EP
TE D
peaks of the PDF
Fig. 16 Example of determining the warning time given the calculated collision probability rate
AC C
Since the collision probability in this case is high (19.1312%), the spacecraft should make a maneuver to avoid collision. It shows in Fig. 16 that the PDF first peaks at 2498s, and then 2730s, and 2969s, thus the advanced warning time for a spacecraft making necessary maneuvers is 2498s. A collision avoiding maneuver can be planned by minimizing the predicted collision probability. For the scenario where the fuel is limited to make the required maneuvers, the maneuver can focus on at least postponing the collision moment to save time for other solutions.
6. Conclusion In this study, we proposed an optical cue-based on-board collision probability analysis method for
- 24 -
ACCEPTED MANUSCRIPT space situational awareness. Optical cues are employed to estimate the hardball radius and relative state between two RSOs, results of which are then incorporated into a probability collision analysis structure that utilizes the UT-transform method to yield the PDF of collision probability. Finally, the trapezoidal
RI PT
integration method is employed to compute the collision probability. The effectiveness of proposed approach was verified by a numerical simulation and an experimental test. Future work will include recognizing more complicated space objects, developing collision maneuver strategies and their
SC
incorporation in a more credible experimental test.
M AN U
Acknowledgement
The work described in this study was supported by the National Natural Science Foundation of China (Grant No. 61701225 and 11672126), the Opening Grant from the Key Laboratory of Space Utilization, Chinese Academy of Sciences (Grant No. LSU-KJTS-2017-01 and LSU-2016-04-01), China Postdoctoral
TE D
Science Foundation (Grant No. 2017M620210). The authors fully appreciate their financial support. We would also like to thank the reviewers for their instructive suggestions.
EP
References
[1] Akella, Maruthi R., and Kyle T. Alfriend. Probability of collision between space objects. Journal of
AC C
Guidance, Control, and Dynamics, 23 (No.5) (2000) pp: 769-772. DOI: 10.2514/2.4611 [2] Kessler, Donald J., and Burton G. Cour Palais. "Collision frequency of artificial satellites: The creation of a debris belt." Journal of Geophysical Research: Space Physics 83.A6 (1978): 2637-2646. DOI: 10.1029/JA083iA06p02637 [3] Nishida S I, Kawamoto S, Okawa Y, et al. Space debris removal system using a small satellite. Acta Astronautica, 65 (No.2) (2009), pp. 95-102. DOI: 10.1016/j.actaastro.2009.01.041 [4] Pardini C, Hanada T, Krisko P H. Benefits and risks of using electrodynamic tethers to de-orbit - 25 -
ACCEPTED MANUSCRIPT spacecraft. Acta Astronautica, 64 (No.5) (2009), pp. 571-588. DOI: 10.1016/j.actaastro.2008.10 [5] McKnight, Darren S., and Frank R. Di Pentino. New insights on the orbital debris collision hazard at GEO. Acta Astronautica 85 (2013), pp. 73-82. DOI: 10.1016/j.actaastro.2013.04.005
RI PT
[6] Braun V, Flohrer T, Krag H, et al. Operational support to collision avoidance activities by ESA’s space debris office. CEAS Space Journal, 8(No.3) (2016), pp: 177-189. DOI: 10.1007/s12567-016-0119-3
SC
[7] Smirnov, N. N., Kiselev, A. B., Smirnova, M. N., Nikitin, V. F. Space traffic hazards from orbital debris mitigation strategies. Acta Astronautica, 109 (2015), pp. 144-152. DOI:10.1016/j.actaastro.2014.09.014 [8] Carpenter, J. R., et al. Sequential probability ratio test for collision avoidance maneuver decisions
M AN U
based on a bank of norm-inequality-constrained epoch-state filters. AAS/AIAA Astrodynamics Specialist Conference, Girdwood, AK. 2011.
[9] DeMars, Kyle J., Yang Cheng, and Moriba K. Jah. Collision probability with Gaussian mixture orbit
10.2514/1.62308
TE D
uncertainty. Journal of Guidance, Control, and Dynamics 37 (No. 3) (2014): 979-985. DOI:
[10] Patera, R. P., Space Vehicle Conflict-Avoidance Analysis, Journal of Guidance, Control, and
EP
Dynamics, 30 (No. 2), (2007), pp. 492–498 DOI: 10.2514/1.24067 [11] Alfriend, K. T., Akella, M. R., Frisbee, J., Foster, J. L., Lee, D.-J., and Wilkins, M., Probability of
AC C
Collision Error Analysis, Space Debris, 1 (No.1), (1999), pp. 21–35. DOI: 10.1023/A:1010056509803 [12] Khutorovsky, Z. N., Boikov, V. F., and Kamensky, S. Y., Direct Method for the Analysis of Collision Probability of Artificial Space Objects in LEO: Techniques, Results, and Applications, Proceedings of the First European Conference on Space Debris, European Space Agency, Paris, France, April 1993, pp. 491–508 [13] McKinley, D. P., Development of a Nonlinear Probability of Collision Tool for the Earth Observing System, AIAA/AAS Astrodynamics Specialist Conference and Exhibit, AIAA Paper 2006-6295, 2006 - 26 -
ACCEPTED MANUSCRIPT [14] Vittaldev V, Russell R P. Space object collision probability via Monte Carlo on the graphics processing unit. The Journal of the Astronautical Sciences, 64(No. 3) (2017) pp. 285-309. DOI: 10.1007/s40295-017-0113-9
RI PT
[15] Morselli, Alessandro, et al. A high order method for orbital conjunctions analysis: Monte Carlo collision probability computation. Advances in Space Research 55 (No.1) (2015) pp. 311-333. DOI: 10.1016/j.asr.2014.09.003
SC
[16] Serra, R., Arzelier, D., Joldes, M., Lasserre, J. B., Rondepierre, A., & Salvy, B. Fast and accurate
M AN U
computation of orbital collision probability for short-term encounters. Journal of Guidance, Control, and Dynamics, 39 (No.5) (2016) pp. 1009-1021. DOI: 10.2514/1.G001353 [17] Dolado-Perez, J. C., Carmen Pardini, and Luciano Anselmo. Review of uncertainty sources affecting the long-term predictions of space debris evolutionary models. Acta Astronautica 113 (2015) pp. 51-65.
TE D
DOI: 10.1016/j.actaastro.2015.03.033
[18] Bérend, N., Estimation of the Probability of Collision Between Two Catalogued Orbiting Objects, Advances in Space Research, 23 (No. 1) (1999), pp. 243–247. DOI: 10.1016/S0273-1177(99)00009-5
EP
[19] Coppola, V. T., Evaluating the Short Encounter Assumption of the Probability of Collision Formula,
AC C
Proceedings of the 22nd AAS/AIAA Space Flight Mechanics Meeting, Vol. 143, Advances in the Astronautical Sciences, Univelt, San Diego, CA, 2012, pp. 2179–2190. [20] Patera, R. P., General Method for Calculating Satellite Collision Probability, Journal of Guidance, Control, and Dynamics, 24( No. 4), 2001, pp. 716–722. DOI: 10.2514/2.4771 [21] Coppola, V. T., Including Velocity Uncertainty in the Probability of Collision Between Space Objects, Proceedings of the 22nd AAS/AIAA Space Flight Mechanics Meeting, Vol. 143, Advances in the Astronautical Sciences, Univelt, San Diego, CA, 2012, pp. 2159–2178.
- 27 -
ACCEPTED MANUSCRIPT [22] Patera, R. P., Method for Calculating Collision Probability Between a Satellite and a Space Tether, Journal of Guidance, Control, and Dynamics, 25 (No. 5) (2002) pp. 940–945 DOI: 10.2514/2.4967 [23] Lidtke A A, Lewis H G, Armellin R, et al. Considering the collision probability of Active Debris
RI PT
Removal missions. Acta Astronautica, 131 (2017) pp. 10-17. DOI:10.1016/j.actaastro.2016.11.012 [24] Segal, Shai, Pini Gurfil, and Kamran Shahid. In-orbit tracking of resident space objects: A comparison of monocular and stereoscopic vision. IEEE Transactions on Aerospace and Electronic Systems 50
SC
(No.1) (2014) pp. 676-688. DOI: 10.1109/TAES.2013.120006
M AN U
[25] Maccone, Claudio. Advantages of Karhunen–Loève transform over fast Fourier transform for planetary radar and space debris detection. Acta Astronautica 60 (No.8-9) (2007) pp.775-779. DOI: 10.1016/j.actaastro.2006.08.015
[26] Thorn, E., Korkiakoski, V., Grosse, D., Bennet, F., Rigaut, F., d'Orgeville, & Smith, C..
TE D
Stereo-SCIDAR system for improvement of adaptive optics space debris-tracking. Advanced Maui Optical and Space Surveillance Technologies Conference, 2017 [27] Gamage S, Cheng Y. CubeSat Detection Using Convolutional Neural Networks. 2018 Space Flight
EP
Mechanics Meeting. 2018: 2234.
AC C
[28] Julier, Simon J. The scaled unscented transformation. Proceedings of the American Control Conference, (Vol. 6) (2002), pp. 4555-4559 [29] Lebedev, V., and Laikov, D., Quadrature Formula for the Sphere of the 131st Algebraic Order of Accuracy, Doklady Mathematics, 59 (No. 3) (1999) pp. 477–481 [30] Mehrholz D, Leushacke L, Flury W, et al. Detecting, tracking and imaging space debris. ESA Bulletin 109 (0376-4265) (2002) pp. 128-134. [31] Ballard, Dana H. Generalizing the Hough transform to detect arbitrary shapes. Readings in computer
- 28 -
ACCEPTED MANUSCRIPT vision. 1987. pp. 714-725. DOI: 10.1016/0031-3203(81)90009-1 [32] Canny, John. A computational approach to edge detection. Readings in Computer Vision. 1987. pp. 184-203. DOI: 10.1109/TPAMI.1986.4767851
RI PT
[33] Quinn, M. N., Jukna, V., Ebisuzaki, T., Dicaire, I., Soulard, R., Summerer, L., Mourou, G. (2015). Space-based application of the CAN laser to LIDAR and orbital debris remediation. The European
AC C
EP
TE D
M AN U
SC
Physical Journal Special Topics, 224(13), 2645-2655. DOI: 10.1140/epjst/e2015-02577-5
- 29 -
ACCEPTED MANUSCRIPT 1. Optical-cues is efficiently utilized to estimate the object information. 2. The GMM models are employed to form the analytical form of collision analysis.
AC C
EP
TE D
M AN U
SC
RI PT
3. MC simulation and an experimental test verify the effectiveness of this approach.
S0094-5765(18)30336-9
DOI:
https://doi.org/10.1016/j.actaastro.2018.06.053
Reference:
AA 6972
To appear in:
Acta Astronautica
Received Date: 14 February 2018 Revised Date:
8 June 2018
Accepted Date: 27 June 2018
Please cite this article as: M. Yu, S. Li, S. Leng, Space collision probability computation based on onboard optical cues, Acta Astronautica (2018), doi: https://doi.org/10.1016/j.actaastro.2018.06.053. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Space Collision Probability Computation Based on On-board Optical Cues Meng Yu1, 2* Shuang Li1, 2† and Shu Leng1 1. College of Astronautics, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China
210016, China
RI PT
2. Advanced Space Technology Laboratory, Nanjing University of Aeronautics and Astronautics, Nanjing
Abstract: This paper presents a novel on-board space collision analysis method for space situational
SC
awareness. The framework is developed under the following assumptions: 1) A satellite can be equipped
M AN U
with on-board sensors for space object recognition. 2) No a–priori knowledge of the space objects is provided. A space object size and relative state estimation method is firstly proposed, wherein optical cues acquired from onboard sensors are utilized to achieve the estimation. Then, the unscented transform approach is employed to calculate the probability density function (PDF) of collision probability based on
TE D
the estimate information. Monte Carlo simulations and an experimental test demonstrate that the proposed approach can achieve high-precision on-board collision probability estimation with an error less than 3%.
1. Introduction
EP
Keywords: Optical cue; Collision analysis; Probability density function; Gaussian uncertainty.
AC C
The problem of vast growing population of resident space objects (RSO) poses a great challenge to long-term safe space operation [1][2]. As the near-earth space environment is becoming more crowded than ever, fast and accurate collision analysis has attracted much attention in the field of space situational awareness. The collision between space objects depends on many attributes, for instance, size, shape, orbital movement, self-rotation of space objects [3]. Considering the lack of a-priori precise knowledge *
Lecture, Department of Astronautics Engineering, Nanjing University of Aeronautics and Astronautics.
†
Professor, Department of Astronautics Engineering, Corresponding Author, Email: [email protected], Tel:
+86(25)84896039. -1-
ACCEPTED MANUSCRIPT about the motion state of space objects, these aforementioned attributes are often associated with uncertainties, which leads to the current probabilistic modeling treatment that generally existed in space collision analysis methods [4].
RI PT
The probabilistic modeling manner aims at evaluating the probability of collision between two space objects. Such a method commonly relies on many assumptions, for instance, a RSO is usually assumed to be influenced by planetary gravitations and J2 perturbation effects. At the same time, its shape is also
SC
approximated as a sphere or an ellipsoid with a determined size in order to simplify the derivation [8]. In
M AN U
addition, many works that evolve computing the collision probability assumed a limited encounter time between two objects, which is also known as short-term encounters [10], wherein the relative positions and velocities are assumed to admit linear forms with Gaussian-distributed uncertainties. For complicated scenarios such as a long-term encounter or RSOs with possible maneuvers, the computation of collision
TE D
probability might be more legitimate when taking into account the non-Gaussian and nonlinear distributed uncertainties in propagating the relative position and velocity [8]. In general, the collision analysis of space objects can be classified into two categories. The first approach
EP
relies on Monte Carlo (MC) sampling [14]. Such a MC based method samples from the probability density
AC C
functions (PDFs) of two RSOs, the relative motion between two RSOs is estimated from these samples and propagated over a certain epoch of interest to record the collision events. The MC sampling methods can be applied to both simple and complicated encounter scenarios with a high accuracy, However, such a technique requires significant computational resources, thus it is currently computationally intractable for a computer system onboard the space vehicle. For lessening the computational burden, Graphics Processing Unit (GPU) can be used to parallelize the computation and then improve real-time performance [14]. The second method focuses on statistically inferring the propagation of RSOs’ relative states over a certain time
-2-
ACCEPTED MANUSCRIPT range, wherein the uncertainty associated with RSOs’ states are usually described by Gaussian or non-Gaussian probability density functions (PDF) [9][16]. The collision probability can then be calculated by integrating the probability density function of relative states over a specific time range and examining
RI PT
the intersection between the relative states and hardball regions ( a hardball region is usually defined as a sphere, whose radius is the sum of actual radii of two RSOs [8]). Specifically, the collision probability can be computed via either determining the volume swept out by the hardball [18] or solving the amount of
SC
influx flow of the PDF of the relative states into the hardball [1][21]. These methods share a similar
M AN U
structure, such that they all capture the relative state uncertainty in a probabilistic manner and use efficient integration methods to calculate the collision probability [9].
To date, most of the researches in the field of space collision probability analysis are based on ground detection and tracking. Some recent works in the field of space situational awareness suggest that efforts
TE D
should be made to enhance the quick response to potential space collision events [23], and space operations might benefit from employing on-board sensors to detect RSOs and conduct collision analysis. Motivated by such demand, this paper proposes a novel approach for on-board space collision analysis using on-board
EP
optical sensors. The proposed approach employs a monocular camera for space object detection and a laser
AC C
detection and ranging device (e.g. LADAR) for measuring the relative distance. Then, the acquired optical data is incorporated to derive the PDF of the relative state between two RSOs. After that, an unscented transform based Gaussian mixture model (GMM) is adopted to deduce the PDF of the collision probability. We follow the common assumption that the RSO can be modeled by an ellipsoid [10][13][14]. Unlike the priori work in the related field of research, this work does not assume a-priori knowledge of the RSO’s size information. Alternatively, we propose to utilize the optical cues acquired from on-board sensors to estimate such information.
-3-
ACCEPTED MANUSCRIPT The rest of this paper is organized as follows. Section 2 gives the brief review, coordinates and problem formulation. Section 3 details the incorporation of camera and LADAR ranging in hardball radius determination and relative state estimation. Section 4 introduces a framework of optical sensor based space
RI PT
collision analysis. Simulation results and discussion in terms of practical concerns of the proposed approach are presented in Section 5. Finally, Section 6 discusses the conclusions of this paper.
2. Coordinate and problem formulation
SC
We first introduce the coordinate frames used in this study. A space-encounter scenario is presented in
M AN U
Fig. 1 as an illustration, denote by S a satellite that is well-tracked by ground stations is operating in the low-earth orbit (LEO), and assume a space object that might have a short or long term encounter with S is detected by optical sensors onboard the satellite S . The sensors consist of a monocular camera and a LADAR, which are attached to one side of the satellite with calibrated installation matrices. The inertial
TE D
frame is denoted by {I} , the satellite body frame is denoted by {S} , the camera and LADAR body frame are denoted by {C},{L} , respectively, and the 2D image frame of the camera is denoted by {im} ,
zS
zI
xS
AC C
EP
which is perpendicular to the y − axis of the satellite body frame.
yS
OC OL
zd
yd
Od
O im
xd
yI
xI
Fig. 1 Illustration of inertial {I}, satellite {S}, object {d}, camera {C}, LADAR {L} frame
-4-
ACCEPTED MANUSCRIPT The position and velocity of space object and satellite in the inertial frame are denoted by
( rdI , r&dI ) , ( rSI , r&SI ) ,
respectively. Denote by x RI = [ rRI , r&RI ] the relative state between them, it can be T
defined as follows:
rRI = rdI − rSI r&RI = r&dI − r&SI
CS and LADAR LC S , the relative state can be expressed
as: T
x = [ r , r& S
L R
]
T
= ( C S )( C I ) x L
S
(2)
I R
M AN U
where
L R
SC
x RC = [ rRC , r&RC ] = ( C C S )( S C I ) x RI L R
RI PT
C
Given the installation matrix of camera
(1)
C I is the matrix that denotes the rotation from the inertial frame {I} to the satellite’s body frame
{S}, as is shown in Fig. 1, the space object in the 2D image frame {im} is represented by a shaded ellipse, whose center is located at pdim = [ µd , vd ]T in {im} , which can be obtained by the pin-hole
TE D
transformation [24] as follows:
µd =
fim pRC ( x) + µ0 pRC ( z )
(3)
EP
f pC ( y) vd = im C R + v0 pR ( z )
where fim represents the focal length of the camera,
[ µ0 , v0 ]
T
AC C
intersects the image frame, rRC = [ pRC ( x ), pRC ( y ), pRC ( z ) ]
T
is the coordinate at which the optical axis
denotes the relative position expressed in the
camera frame as state in Eq. (2). The on-board optical measurement at time step t can be expressed as:
Z d (t) = {[ µd (t ), vd (t ) ] , Ld (t ) }
(4)
where Ld (t ) is the LADAR range measurement, which is a scalar. Given the definition of coordinate frames and measurement model, the problem of collision between the satellite and the space object can be formulated as follows. A space object d is determined as a collision threat to satellite if the relative position rR (t ) falls -5-
ACCEPTED MANUSCRIPT within the hardball region Λ H at time t > t0 , where the hardball radius H Sd is defined as the sum of the radii of satellite and space object. The collision probability Pc (t ) is a function of time t , hardball region
Λ H and relative state xRI that is defined in Eq. (2), Pc (t ) can be calculated by integrating the PDF
t
Pc (t ) = ∫ pc (λ )d λ
(5)
−∞
where pc (t ) is a probability density function that is defined as:
∫
A vR ≤ 0
v R p ( x RI ( t ; t0 ) ) dvt dA
(6)
SC
pc ( t ) = ∫
RI PT
pc (t ) over time t [9]:
where vR = vt ⋅ rˆ , A denotes the surface area of the hardball H dS , rˆ denotes the normal vector with
M AN U
respect to the area element dA , vt = r&dI − r&SI denotes the relative velocity, and vR is a scalar that denotes the amount of relative velocity projected onto the normal vector rˆ . The integration over vR ≤ 0 indicates the influx of probability distribution for the relative state into the hardball region, which can be
TE D
interpreted as the boundary of the hardball that is crossed by the relative velocity. vt . As illustrated in Fig. 2, the collision probability can be calculated by integrating over the crossed surface of hardball.
rˆ
dA
rR
EP
vt
θ
AC C
φ
hardball
Fig. 2 Illustration of influx to the hardball region
By employing the spherical integration rule, Eq. (5) can be expressed as [19]:
Pc ( t ) =
t0 + T 2 π
π 2
∫ ∫ ∫ ( H Sd )
t0
0 − π2
2
cos θ
∫
vt ⋅rˆ ≤ 0
vt ⋅ rˆ
∞
∫ρ
S
( xS , t ; t0 ) ρ d ( xd , t ; t0 ) dxS dv R dθ dφ dt (7)
−∞
where
θ , φ represent the spherical angles, H Sd is the radius of the hardball, T is the total integration
time,
ρ S , ρd are the PDF for the state of the satellite and the space object, respectively. The process of
-6-
ACCEPTED MANUSCRIPT spherical integration is illustrated in Fig. 3:
dA dθ
rdθ r cos θd φ
φ
RI PT
θ
r = H Sd
dA = r2 cos θ dθ dφ
dφ
Fig. 3 Illustration of the spherical integration used in Eq. (7)
As shown in Fig. 3, wherein the area integral can be denoted by dA = (H S ) cos θ dθ dφ . Note that in Eq. (7), the differential term for
SC
d 2
ρS and ρ d are dx S and dxd , which can be reorganized as:
M AN U
dxS dxd = dx S d ( xS + x R ) = dxS dx R = dxS drR dv R = ( H Sd ) cos θ dxS dv R drdθ dφ = dAdr 2
(8)
dr dAdt = vt ⋅ rˆ dAdt dt
=
By such manner, Eq. (8) can be readily related to Eq. (6). The computation of Pc requires the
TE D
knowledge of H Sd , ρ S , ρ d , the determination of which will be detailed in Section 3.
3. Optical Cue based metric determination
AC C
collision analysis.
EP
In this section, we propose an approach to determine the key metrics of RSO as illustrated in Fig. 4 for
de br is
z c (L )
Ladar ranging
α max
{C }({L})
φL
ϕL
tr av
e ers
z on
e θd
hardball
αL
a ccess ible zone
yc (L )
xc ( L ) Fig. 4 Space object size determination based on optical cues
As shown in Fig. 4. We assume that a space object can be precisely detected by object detection
-7-
ACCEPTED MANUSCRIPT methods (related works on object detection can be found in [25]), and the ranging direction of on-board LADAR can be controlled via sending the azimuth/elevation command:
ϕcom = arccos
(9)
f im
RI PT
φcom =
f im2 + µ d2 lˆ
f + µ d2 2 im
where pim = [ µd , vd ]T denotes the central coordinate of the detected object in the image frame;
µd
vd ]T represents a three-dimensional vector from the satellite to the space object.
M AN U
and lˆ = [ f im
SC
ϕcom , φcom are the elevation and azimuth angle control command of the on-board LADAR, respectively;
Denote by Lˆ d the LADAR bearing angle, note that sometimes the continuous bearing adjustment of LADAR might be too computationally intensive, thus some deviations between the LADAR bearing angle
Lˆ d and lˆ should be allowed. We first define the deviation angle between Lˆ d and lˆ as α% L , and
TE D
define an accessible zone (see Fig. 4), such that only the LADAR measurement that falls within the accessible zone is considered valid for distance ranging. The angle of field of view (FoV) for the accessible zone can be determined by
rim | lˆ |
EP
α max = arctan
(10)
AC C
As illustrated in Error! Reference source not found., denote by lin the location at which the laser
ˆ in = cin − lin can be beam hits the space object; and cin the center of the space object, a vector p generated to estimate the radius of the space object Rd . Denote by
θ d the angle between pˆ in and Lˆ d ,
it can be related to α% L by noting that:
θd ≈
π α% L ⋅ 2 α max
(11)
Defining the radius of the space object as Rd . We can estimate Rd based upon the following
-8-
ACCEPTED MANUSCRIPT geometrical relationship:
rim = Rd
| lˆ |
(12)
Rd2 + Ld + 2 Rd Ld cos θ d 2
which yields 2
Solving for the root of Eq. (13), and noting that | lˆ |>> rim yields:
rim | Ld | | lˆ |2 −rim2 sin 2 θ d + rim2 | Ld | cos θ d | lˆ |2 −r 2 im
(14)
SC
Rd =
(13)
RI PT
( rim2 − | lˆ |2 ) Rd2 + 2rim2 | Ld | cos θ Rd + rim2 Ld = 0
Then, the hardball radius can be obtained by H Sd = Rd + RS , where RS is the radius of the satellite.
Rd2 + Ld + 2 Rd Ld cos θ d as the distance from the space object to the satellite, the
M AN U
Defining rSD =
2
relative position in the camera frame can be expressed as
rRc =
lˆ ⋅ rSD | lˆ |
(15)
TE D
Given rRc , the relative velocity in the camera frame can be approximated by solving the differential of position:
EP
r&Rc (t ) = lim
rRc ( t + ∆t ) − rRc ( t )
∆t →0
∆t
(16)
AC C
4. Computation of collision probability The computation of collision probability is carried out given the estimation of the relative state
x Rc = [ rRc , r&Rc ] and the hardball radius H Sd . Such an estimate would inevitably contain measurement T
noises, which results in a non-linear uncertainty. We employ the unscented transform (UT) [28] to cope with this problem due to its robustness against non-linear noise. First, define the joint PDF of the states of the satellite and the space object as
p ( xS , xd ) = p ( x S , x S + x R ) -9-
(17)
ACCEPTED MANUSCRIPT The PDF of the relative state x R can be expressed by marginalizing over xS , yielding:
p ( x R ) = ∫ p ( x S , x S + x R ) d xS
(18)
To facilitate the computation of collision probability, the UT-transform [28] is employed to
variable, (e.g. x R , x S or xd ), such that:
mi = x +
(
mi + n = x −
(n + κ ) Pxx
(
)
(n + κ ) Pxx
W0 = κ / (n + κ )
Wi = 1 / 2(n + κ ), i = 1,..., n
i
)
i
Wi + n = 1 / 2(n + κ ), i = 1,..., n
SC
m0 = x
RI PT
characterize the PDF. Specifically, 2n + 1 sigma points are generated based on a n − dimensional
(19)
M AN U
Where x and Pxx denote the mean and covariance of a variable, respectively; W0~2 n denotes
the weights for the generated sigma points m0~2 n . Using the UT-transform, PDFs for the states of the satellite and the space object can be expressed as: US
p ( xS ) = ∑ Wi S pg ( xS ; miS , PxxS ) j =1
Ud
(20)
TE D
p ( xd ) = ∑ Wi pg ( xd ; m , P ) b
d j
d xx
i =1
where U S = U d = 2n + 1 . pg ( x, m, P ) denotes a standard Gaussian PDF. Assuming the satellite state
EP
and object state are independent, we can break their joint PDF into a product of two PDFs, such that:
p( x S , xd ) = p( x S ) p( xd ) U s Ud
AC C
= ∑∑ Wi SWi d pg ( x S ; miS , PxSx ) pg ( xd ; m dj , Pxxd )
(21)
i =1 j =1
The authors in [21] proved that
∫ p (x g
S
; m S , PxxS ) pg ( xS + x R ; md , Pxxd ) dx S = pg ( x R ; uE , Σ E )
(22)
where uE = mS − md , Σ E = PxxS + Pxxd . Substituting Eq.(21), (22) into Eq. (18), the PDF for the relative state can be rewritten as U s Ud
p ( x R ) = ∑∑ Wi SW jd pg ( x R ; uE(ij ) , PE(ij ) ) i =1 j =1
- 10 -
(23)
ACCEPTED MANUSCRIPT Following the method proposed in [21], the PDF for relative state is decomposed into the product of PDFs for the relative position and velocity to maintain the computational stability. Specifically, uE , PE are decomposed into:
u
Σ(rvij ) Σ(vijv )
(24)
RI PT
Σ(rrij ) ur(ij ) ( ij ) = (ij ) , PE = T ( Σ(rijv ) ) uv
( ij ) E
Substituting Eq. (24) into Eq. (23), and following the decomposing rule proposed in [15], it shows that
SC
p( xR ) can be expressed as U s Ud
p ( x R ) = ∑∑ Wi SW jd pg (rR ; ur(ij ) , Σr(ij ) ) pg ( v '; ( uv(ij ) ) ', ( Σv(ij ) ) ' ) i =1 j =1
M AN U
where x R = [rR , v R ]T , and
(25)
v ' = v R − ( Σ(rvij ) ) ( Σ(rrij ) ) rR −1
T
( uv(ij ) ) ' = Σ(vvij ) ( Σ(rrij ) )
−1
ur(ij )
( Σ(vij ) ) ' = Σ(vvij ) − Σ(vrij ) ( Σ(rrij ) )
−1
(26)
Σ(rijv )
TE D
Substituting Eq. (25) into Eq. (6) and reallocate the differential elements, such that the integration is taken over the hardball region, the collision probability density function can be expressed as US Ud
pc ( x R ) = ∑∑ Wi SW jd ∫ pg (rR ; ur(ij ) , Σ (rij ) )Θ ( rR , v R ) dA
(27)
A
EP
i =1 j =1
AC C
where Θ ( rR , v R ) is the integration over the relative velocity that is expressed as:
Θ ( rR , v R ) = ∫ ( ij )
vn ≤ 0
vn( ij ) pg ( v '; ( uv( ij ) ) ', ( Σ (vij ) ) ') dv '
(
where vn( ij ) is a scalar and admits the form vn( ij ) = v '⋅ rˆ + rˆT Σ (rvij )
(28)
) ( Σ( ) ) T
ij rr
−1
rR . Given the hardball
radius H Sd estimated by the approach proposed in Section 3, and using the spherical transformation (see Eq. (7)), pc ( x R ) can be expressed as
pc ( xR ) = ( H Sd )
2
US Ud
∑∑W
i
i =1 j =1
S
W jd ∫
2π
0
∫
π 2
− π2
pg ( H Sd rˆ; ur(ij ) , Σ(rij ) )Ξ ( rˆ, uE , PE ) cos θ dθ dφ (29)
As is defined in [21], Ξ ( rˆ, uE , PE ) is expressed as - 11 -
ACCEPTED MANUSCRIPT Ξ ( rˆ, uE , PE ) =
Ξ 02 ( rˆ ) Ξ 0 ( rˆ ) σ ( rˆ ) Ξ 0 ( rˆ ) exp − 2 1 − erf − 2 2π 2σ ( rˆ ) 2σ ( rˆ )
(30)
where
σ
2
( rˆ ) = rˆ
T
(
Σ( ij ) vv
−(
Σ (ij ) rv
)( T
Σ (ij ) rr
)
−1
Σ( ij ) rv
) rˆ
(31)
RI PT
T −1 Ξ 0 ( rˆ ) = rˆ T uv( ij ) + ( Σ(rvij ) ) ( Σ(rrij ) ) ( H Sd rˆ − ur(ij ) )
Thus the calculation of Θ ( rR , v R ) can be conducted analytically. The rest of the integrals over the
θ and φ in Eq. (29) corresponds to integrating the relative position over the hardball
SC
spherical element
sphere. The spherical integration over the relative position cannot be performed analytically. Alternatively,
M AN U
it can be given by numerical approximation, for example, using a numerical quadrature integration approach [9]. In this study, we choose the Lebedev and Laikov’s method [29] to approximate the integration over the relative position. Lebedev and Laikov’s method is a kind of quadrature approach that can approximate the surface integral over a sphere by a summation of weighted N quadrature points nˆ k .
pc (t ) = ( H dS )
2
TE D
In this manner, the PDF for collision probability at time t can be expressed as Us Ud
N
i =1 j =1
k =1
∑∑ Wi,StW jd,t ∑ wk pg ( H dS nˆ k ; ur(ij,t ) , Σ(rij,t) ) Ξ ( nˆ k , uE ,t , PE ,t )
(32)
EP
The PDF for collision probability can be solved by substituting Eq. (32) into Eq. (5) assuming that
AC C
Pc (t0 ) = 0 at t = t0 , the probability of collision can be expressed as
Pc (t ) = ∫ pc ( λ )d λ t
0
(33)
Such integration over time can be efficiently approximated using quadrature rules, such as the trapezoidal integral method.
5. Simulation results and discussion 5.1 Numerical simulation To demonstrate the effectiveness of the proposed approach for the on-board computation of the
- 12 -
ACCEPTED MANUSCRIPT probability of collision, three space encounter scenarios are considered. Denote by S a near-earth satellite, and B a space object orbiting nearby the satellite S . The space object B is assumed to possess a regular ellipsoid shape, its size is set to small/medium/large in the 1st/2nd/3rd scenario, respectively, which
velocity setup for satellite in the inertial frame is given as:
SC
−6074.6621 −2438.9919 m / s −258.4168 0.9829 0.1294 0.1305 I QS = 0.1125 −0.9851 0.1294 0.1453 −0.1125 −0.9829
(34)
M AN U
[ pS 0
−4017.3356 v S 0 ] = 6435.5823 km 6207.7969
RI PT
might correspond to a low/medium/high collision probability with the satellite S . The initial position and
where I QS is a matrix that defines the rotation from the satellite frame to the inertial frame. We assumed that the on-board camera has just observed a space object at t = t0 . The satellite is assumed to be
TE D
well-tracked by ground stations, its observation uncertainty is characterized by a single Gaussian model. The covariance P0 = diag ( PS ,r0 , PS ,v0 , PS ,q0 ) of this Gaussian distribution at the initial observation time
t = t0 in the inertial frame is given as
EP
PS ,r0 = diag ([ 0.08,0.08, 0.08]) km2 , PS ,v0 = diag ([1.5,1.5,1.5]) m 2 / s 2 , PS ,q0 = diag ([ 0.1, 0.1, 0.1]) deg 2
(35)
AC C
The covariance is chosen according to the results acquired from the ground tracking of near-earth satellites [30]. The space object is assumed to be observed onboard for the first time at
t = t0 with no
a-priori knowledge of its shape and size, indicating an unknown magnitude of uncertainty in its initial state
estimation. To cope with this issue, the uncertainty of its state is characterized by a Gaussian mixture model acquired from the UT-transform. The initial position and velocity for the space object in the camera frame were given by:
- 13 -
ACCEPTED MANUSCRIPT −1430.9141 rBC0 = 4104.0239 m, −5154.4521
0.0651 r&BC0 = −0.1677 m / s 0.2153
(36)
We set the simulation time to 60 mins and assumed that the trajectories for satellite and space object in
RI PT
all three scenarios were identical. In this manner, the only difference in each scenario is the size and shape of the space object. Consider the standard two-body problem, the time history of the relative distance acquired by propagating the initial means is presented in Fig. 5. 7000
5000
400 350 300 250 200 150
M AN U
relative distance (m)
6000
SC
relative re lati vedistance distance (m) (m)
450
100 50
4000
2440
2460
2480
2500
2520
2540
2560
2580
Simulation time (s)
Simulation time (s)
3000 2000
0 0
TE D
1000
500
1000
1500
2000
2500
3000
3500
4000
Simulation time (s)
Fig. 5 Ground truth of the time history for the relative distance
EP
It is shown from Fig. 5 that the true minimum distance between satellite and the space object is
AC C
44.54m, indicating a close fly-by of the space object with respect to the satellite. We additionally developed a simple space object simulator in Matlab software for simulating the on-board optical cue and set the hardball radius to 1m, 5m, 25m in 1st/2nd/3rd scenarios, respectively. An example of a space object fly-by is given in Fig. 6:
- 14 -
RI PT
ACCEPTED MANUSCRIPT
detection
SC
Fig. 6 Example of a space object close fly-by in synthesized images, Hough transform is used for object
In this example, a space object was assumed to sweep over the target satellite. Fig. 6 presents four
M AN U
synthesized images, which are captured by the on-board camera during the fly-by. [ µ d , vd ] represents the center of the space object, and rim represents the image radius of the space object. The space object is assumed to appear as an ellipse in a synthesized image, classic contour detection algorithms such as Hough
TE D
transform [31] or Canny edge fitting [32] can be utilized to detect the space object. The focal length of camera is set at 0.5 cm, the pixel size is set to 4 um. Denote by
Z B 0 = ( µ0 , v0 )T , L0
T
the initial measurement of the space object, where pim ,0 = [ µ0 , v0 ]
EP
center of the object in the image frame, and
AC C
and the relative state in the camera frame
T
is the
L0 is the LADAR ranging distance. The actual object radius
[ rRc , r&Rc ]
T
can be estimated by following Eqs. (10) to (16),
the
corresponding measurement noise covariance R is taken to be diagonal with image location with a value of
2
0.3 pixel 2 and LADAR ranging measurement with a value of 0.2% L , such that
R0 = diag ([0.3 pix 2 , 0.3 pix 2 , 0.2% L ]) . 2
The time history of estimation error for the object radius and relative state in this example is depicted in Fig. 7:
- 15 -
RI PT
ACCEPTED MANUSCRIPT
SC
Fig. 7 Time history of radius and relative state estimation error
M AN U
It can be seen from Fig. 7 that the estimation error for object radius and relative state remains within 1% during the entire simulation, and the estimation for radius error exhibits an unstable jitter around 0.4%, which is partially due to the non-Gaussian noises added to the simulated image. The estimation for relative position and velocity is relatively stable with an error below 0.2%. Simulations of other examples also
satellite and the space object.
TE D
indicate that the radius estimation error is approximately proportional to the relative position between the
The initial measurement covariance is used to determine the weights, mean, and covariance in
EP
UT-transform. After that, the Lebedev and Laikov’s method [29] is applied to generate the Gaussian
AC C
mixture models in each time step, and then apply Eq. (32) to compute the PDF of collision probability, note that using more integral points from Lebedev and Laikov’s method can improve the overall accuracy, however, at the expense of more computational effort. In this simulation, we choose N k = 194 in Lebedev and Laikov’s method, which is a compromise between accuracy and computational cost. The time history of collision probability and its PDF in three scenarios is presented in Fig. 8, a)~f):
- 16 -
PDF of collision probability in 1st scenario.
c)
PDF of collision probability in 2nd scenario.
e)
PDF of collision probability in 3rd scenario.
b) collision probability in 1st scenario
M AN U
SC
Collision Probability, 2nd scenario
a)
RI PT
ACCEPTED MANUSCRIPT
EP
TE D
d) collision probability in 2nd scenario
f)
collision probability in 3rd scenario
AC C
Fig. 8 Collision probability and its PDF in three scenarios
We adopted the trapezoidal integration method to compute the collision probability. As is shown in Fig. 8 a)~f), the final collision probability for the first/second/third scenario is 0.1207% , 3.5142% and 18.6421% according to Fig. 8 b), d) and e). The PDF for collision probability peaks at around 2500s throughout three scenarios and returns to zero after roughly 3000s, which is consistent with the time history of relative distance (see Fig. 6), indicating an accurate reflection of the trend of the relative state variation. In addition, the PDF of collision probability exhibits another peak at around 3000s in Fig. 8 e), which is
- 17 -
ACCEPTED MANUSCRIPT mostly due to the uncertainty in object detection. As we set the radius of hardball to 25m in this scenario (the radius of satellite is 3m, and the radius of space object is 22m), the space object possessed a large part of the synthesized image, which brings ambiguity to the estimate of the object’s center. The maximum error
RI PT
of hardball radius estimation reaches 3% in this scenario, which is nearly 3 times bigger than that in other two scenarios. The improvement over the radius estimation of a space object with irregular shape or a large size is left to future work.
SC
To verify the collision probability computation, an extensive Monte Carlo simulation is carried out.
M AN U
The space object’s orbit samples are drawn from the initial uniform distributions p ( x S0 ) , p ( x B0 ) in each scenario and propagated over the entire simulation time. For each MC run in each scenario, the orbital elements of satellite and object B are propagated for 60 mins at a time step of 0.5 s. The MC simulation stops if the position between object B and satellite S is smaller than the hardball radius, which is
TE D
regarded as a space collision event. The collision probability is calculated by dividing the total number of Monte Carlo runs by the total number of collision events. For each scenario, 1 million MC runs are carried out, the variation of collision probability along with the number of MC runs for three scenarios are shown
Collison probability, P c
AC C
MC runs.
EP
in Fig. 9 a), b), c), respectively, demonstrating that the collision probability converges after roughly 4 × 105
a) MC based Collision probability for 1st scenario
b) MC based Collision probability for 2nd scenario - 18 -
ACCEPTED MANUSCRIPT 0.25
0.2
0.15
0.1
0 0
2
4
6
8
10
Number of MC Runs for the 3rd scenario
RI PT
0.05
12 10 5
c) MC based Collision probability for 3rd scenario
SC
Fig. 9 Collision probability based on Monte Carlo simulation for three scenarios (1 million trails)
In the first scenario, 1187 collision events were triggered for 1 million MC runs, which yields a
M AN U
collision probability of Pc ,1st = 0.1187% , the computation of collision probability for the 2nd and 3rd scenario is the same as in the first scenario, which are Pc ,2 nd = 3.4233% , Pc ,3rd = 19.1232% . The proposed algorithm differs by 1.65%, 2.65% and 2.51% from Monte Carlo results for the three scenarios,
TE D
respectively, demonstrating a plausible accuracy in collision probability calculation. In addition, since the proposed approach relies on the assumption that space object B can be approximated by a regular ellipsoid, its application in practice might benefit from introducing a more accurate object recognition
EP
algorithm.
AC C
5.2 Experimental Test
To further validate the applicability of the proposed approach, a semi-physical experiment was additionally conducted. The experimental set-up is exhibited in Fig. 10
- 19 -
SC
RI PT
ACCEPTED MANUSCRIPT
M AN U
Fig. 10 Illustration of the complete experimental set-up
As shown in Fig. 10, a Laser Tracker serves as the on-board optical sensor suite, a spherically mounted retro-reflector (SMR) ball is mounted on the last joint of a 6D manipulator to simulate the space object, and the 6D manipulator is employed to approximate the orbital movement of the space object.
TE D
Specifically, we first simulate an orbit and downsize such orbit based upon the parameters of the 6D manipulator (for example, its total length is 1.76m), such that the last joint of this manipulator is
EP
programmed to follow the scaled orbital movement. The laser tracker is set to remain stationary to track the target SMR ball, from which the size information and relative state is computed and then employed in
AC C
collision probability estimation framework. The experiment lasts 20 mins, the data acquisition frequency of the optical sensor is set to 1Hz. Some snapshots of the target SMR ball movement during the experiment are shown in Fig. 11:
Fig. 11 Moments of the target ball in motion during the experiment - 20 -
ACCEPTED MANUSCRIPT Note that the target SMR ball is an accessory of the Laser Tracker, therefore the Laser Tracker can automatically recognize the SMR ball and range the relative distance. The simulated orbit and the location
RI PT
of the optical sensor is shown in Fig. 12. The time history of the relative distance is presented in Fig. 13
simulated orbit Optical Sensor Location
t0 1.5
tend
1 0.5
1
4
SC
Z axis (m)
2
0
3
-1
2
-2
1 -3
0
X axis (m)
M AN U
Y axis (m)
-1
-4
EP
TE D
Fig. 12 Illustration of the simulated orbit and the sensor location
AC C
Fig. 13 Time history of the relative distance (measured by the Laser Tracker)
As shown in Fig. 12 and Fig. 13, the minimum relative distance is 0.1682m, which occurs at 265s. The relative distance data is captured by the Laser Tracker with an accuracy of 0.1um, which serves as the ground truth for validating the accuracy the PDF estimation. We set the radius of the hardball to 0.20m, and follow the procedure established in section 5.1 to calculate the collision probability rate. Specifically, we used the optical cues acquired during the beginning to infer the collision probability and its PDF, results of which are shown in Fig. 14:
- 21 -
ACCEPTED MANUSCRIPT 0.4 0.35 0.3 0.25 0.2 0.15
0.05 0 0
240
RI PT
0.1
480
720
960
1200
Experiment Time (s)
Fig. 14 Experiment test result: PDF of collision probability (left) and collision probability (right)
SC
As depicted in Fig. 14, the time history of PDF of collision probability generally agrees with that of
M AN U
the relative distance in Fig. 13. Note that the PDF of collision probability would in theory peak at the location where the relative distance between two RSOs reaches its minimum. In this experiment, the PDF peaks at 267s, which only differs from the ground truth by 2s, thus showing a plausible estimation accuracy. Unfortunately, since it is impractical to actually let the target ball hit the satellite in this experiment. We
TE D
cannot validate the accuracy of collision probability, a more credible experimental test for validating the accuracy of the collision probability estimate is left for future work. 5.3 Practical Concerns
EP
The aforementioned on-board collision probability computation framework aims at enabling the
AC C
quick response to space threatens, especially the newly spawned objects. To implement this framework in practice, the capability of on-board optical sensors must be primarily considered. Choice of on-board Sensors The space object detection prefers a camera with a long focal length. When dealing with the detection of small and irregular space objects, a camera with long focal length will yield a higher physical resolution than the one with short focal length. For example, suppose the relative distance between a satellite and a space debris varies from 10km to 4km, the actual radius of the space debris and the pixel block is 5~20m
- 22 -
ACCEPTED MANUSCRIPT and 3e-6m，respectively. The relationship between the camera focal length and the image area of the object
M AN U
SC
RI PT
is shown in Fig. 15:
Fig. 15 Relationship between the focal length (m) and object’s image area (Number of pixels)
TE D
In Fig. 15, the actual radius of the space debris is set to 5m, 10m, 15m, 20m, respectively. It shows that the image area of the space debris increases as the focal length is enlarged from 5mm to 0.3m. Taking the squared case (actual debris radius is 5m) highlighted in Fig. 15 for example, and assuming the relative
EP
distance is 10km, an accurate estimate of the debris radius would need at least 300 pixels, which requires the focal length to be at least 0.15m. However, using a camera with long focal length also brings difficulty
AC C
to sensor installment and mass distribution. Moreover, a narrow FOV camera also limits the detection rate of space objects. Thus in practice, a balance between accuracy and feasibility must be properly stroked. Besides camera, the detection range of the onboard LIDAR is also an important factor. To save enough time for making a collision avoidance maneuver decision, a spacecraft needs to detect the space objects as early as possible. Therefore for safety concerns, an onboard LADAR with a long detection range is favored. Recent developments of onboard LIDAR design show that the LIDAR system will be able to detect a 10cm size space object given a relative distance of roughly 100km [33]. A major practical concern of an onboard LIDAR system is its high power consumption, an onboard system might not allow a long range LIDAR sensor to operate in full-time. To cope with this issue, sequential images acquired from the onboard camera
- 23 -
ACCEPTED MANUSCRIPT can be employed to estimate the relative state while the LIDAR is offline. This alternative method has the advantage of much lower power consumption, but it is also less accurate than the LIDAR system. An autonomous sensor switching mechanism can be designed to remedy this issue. For instance, we can wake up the LIDAR ranging system when the estimated relative distance is below a threshold, and actively guide
RI PT
the LIDAR to measure the relative distance to improve the calculating of the collision probability. Warning Time for Collision Avoidance
The collision avoiding maneuver needs to be carried out if the calculated collision probability is high. Such a maneuver requires an enough advanced warning time. To estimate the warning time, one can think
SC
of the epoch at which the collision probability rate peaks as the time where a collision event might happen, thus estimating the warning time can be achieved by finding the peaks of the calculate PDF for collision
M AN U
probability. Fig. 16 shows the peaks of the collision probability rate that corresponds to the case in Fig. 8, e):
EP
TE D
peaks of the PDF
Fig. 16 Example of determining the warning time given the calculated collision probability rate
AC C
Since the collision probability in this case is high (19.1312%), the spacecraft should make a maneuver to avoid collision. It shows in Fig. 16 that the PDF first peaks at 2498s, and then 2730s, and 2969s, thus the advanced warning time for a spacecraft making necessary maneuvers is 2498s. A collision avoiding maneuver can be planned by minimizing the predicted collision probability. For the scenario where the fuel is limited to make the required maneuvers, the maneuver can focus on at least postponing the collision moment to save time for other solutions.
6. Conclusion In this study, we proposed an optical cue-based on-board collision probability analysis method for
- 24 -
ACCEPTED MANUSCRIPT space situational awareness. Optical cues are employed to estimate the hardball radius and relative state between two RSOs, results of which are then incorporated into a probability collision analysis structure that utilizes the UT-transform method to yield the PDF of collision probability. Finally, the trapezoidal
RI PT
integration method is employed to compute the collision probability. The effectiveness of proposed approach was verified by a numerical simulation and an experimental test. Future work will include recognizing more complicated space objects, developing collision maneuver strategies and their
SC
incorporation in a more credible experimental test.
M AN U
Acknowledgement
The work described in this study was supported by the National Natural Science Foundation of China (Grant No. 61701225 and 11672126), the Opening Grant from the Key Laboratory of Space Utilization, Chinese Academy of Sciences (Grant No. LSU-KJTS-2017-01 and LSU-2016-04-01), China Postdoctoral
TE D
Science Foundation (Grant No. 2017M620210). The authors fully appreciate their financial support. We would also like to thank the reviewers for their instructive suggestions.
EP
References
[1] Akella, Maruthi R., and Kyle T. Alfriend. Probability of collision between space objects. Journal of
AC C
Guidance, Control, and Dynamics, 23 (No.5) (2000) pp: 769-772. DOI: 10.2514/2.4611 [2] Kessler, Donald J., and Burton G. Cour Palais. "Collision frequency of artificial satellites: The creation of a debris belt." Journal of Geophysical Research: Space Physics 83.A6 (1978): 2637-2646. DOI: 10.1029/JA083iA06p02637 [3] Nishida S I, Kawamoto S, Okawa Y, et al. Space debris removal system using a small satellite. Acta Astronautica, 65 (No.2) (2009), pp. 95-102. DOI: 10.1016/j.actaastro.2009.01.041 [4] Pardini C, Hanada T, Krisko P H. Benefits and risks of using electrodynamic tethers to de-orbit - 25 -
ACCEPTED MANUSCRIPT spacecraft. Acta Astronautica, 64 (No.5) (2009), pp. 571-588. DOI: 10.1016/j.actaastro.2008.10 [5] McKnight, Darren S., and Frank R. Di Pentino. New insights on the orbital debris collision hazard at GEO. Acta Astronautica 85 (2013), pp. 73-82. DOI: 10.1016/j.actaastro.2013.04.005
RI PT
[6] Braun V, Flohrer T, Krag H, et al. Operational support to collision avoidance activities by ESA’s space debris office. CEAS Space Journal, 8(No.3) (2016), pp: 177-189. DOI: 10.1007/s12567-016-0119-3
SC
[7] Smirnov, N. N., Kiselev, A. B., Smirnova, M. N., Nikitin, V. F. Space traffic hazards from orbital debris mitigation strategies. Acta Astronautica, 109 (2015), pp. 144-152. DOI:10.1016/j.actaastro.2014.09.014 [8] Carpenter, J. R., et al. Sequential probability ratio test for collision avoidance maneuver decisions
M AN U
based on a bank of norm-inequality-constrained epoch-state filters. AAS/AIAA Astrodynamics Specialist Conference, Girdwood, AK. 2011.
[9] DeMars, Kyle J., Yang Cheng, and Moriba K. Jah. Collision probability with Gaussian mixture orbit
10.2514/1.62308
TE D
uncertainty. Journal of Guidance, Control, and Dynamics 37 (No. 3) (2014): 979-985. DOI:
[10] Patera, R. P., Space Vehicle Conflict-Avoidance Analysis, Journal of Guidance, Control, and
EP
Dynamics, 30 (No. 2), (2007), pp. 492–498 DOI: 10.2514/1.24067 [11] Alfriend, K. T., Akella, M. R., Frisbee, J., Foster, J. L., Lee, D.-J., and Wilkins, M., Probability of
AC C
Collision Error Analysis, Space Debris, 1 (No.1), (1999), pp. 21–35. DOI: 10.1023/A:1010056509803 [12] Khutorovsky, Z. N., Boikov, V. F., and Kamensky, S. Y., Direct Method for the Analysis of Collision Probability of Artificial Space Objects in LEO: Techniques, Results, and Applications, Proceedings of the First European Conference on Space Debris, European Space Agency, Paris, France, April 1993, pp. 491–508 [13] McKinley, D. P., Development of a Nonlinear Probability of Collision Tool for the Earth Observing System, AIAA/AAS Astrodynamics Specialist Conference and Exhibit, AIAA Paper 2006-6295, 2006 - 26 -
ACCEPTED MANUSCRIPT [14] Vittaldev V, Russell R P. Space object collision probability via Monte Carlo on the graphics processing unit. The Journal of the Astronautical Sciences, 64(No. 3) (2017) pp. 285-309. DOI: 10.1007/s40295-017-0113-9
RI PT
[15] Morselli, Alessandro, et al. A high order method for orbital conjunctions analysis: Monte Carlo collision probability computation. Advances in Space Research 55 (No.1) (2015) pp. 311-333. DOI: 10.1016/j.asr.2014.09.003
SC
[16] Serra, R., Arzelier, D., Joldes, M., Lasserre, J. B., Rondepierre, A., & Salvy, B. Fast and accurate
M AN U
computation of orbital collision probability for short-term encounters. Journal of Guidance, Control, and Dynamics, 39 (No.5) (2016) pp. 1009-1021. DOI: 10.2514/1.G001353 [17] Dolado-Perez, J. C., Carmen Pardini, and Luciano Anselmo. Review of uncertainty sources affecting the long-term predictions of space debris evolutionary models. Acta Astronautica 113 (2015) pp. 51-65.
TE D
DOI: 10.1016/j.actaastro.2015.03.033
[18] Bérend, N., Estimation of the Probability of Collision Between Two Catalogued Orbiting Objects, Advances in Space Research, 23 (No. 1) (1999), pp. 243–247. DOI: 10.1016/S0273-1177(99)00009-5
EP
[19] Coppola, V. T., Evaluating the Short Encounter Assumption of the Probability of Collision Formula,
AC C
Proceedings of the 22nd AAS/AIAA Space Flight Mechanics Meeting, Vol. 143, Advances in the Astronautical Sciences, Univelt, San Diego, CA, 2012, pp. 2179–2190. [20] Patera, R. P., General Method for Calculating Satellite Collision Probability, Journal of Guidance, Control, and Dynamics, 24( No. 4), 2001, pp. 716–722. DOI: 10.2514/2.4771 [21] Coppola, V. T., Including Velocity Uncertainty in the Probability of Collision Between Space Objects, Proceedings of the 22nd AAS/AIAA Space Flight Mechanics Meeting, Vol. 143, Advances in the Astronautical Sciences, Univelt, San Diego, CA, 2012, pp. 2159–2178.
- 27 -
ACCEPTED MANUSCRIPT [22] Patera, R. P., Method for Calculating Collision Probability Between a Satellite and a Space Tether, Journal of Guidance, Control, and Dynamics, 25 (No. 5) (2002) pp. 940–945 DOI: 10.2514/2.4967 [23] Lidtke A A, Lewis H G, Armellin R, et al. Considering the collision probability of Active Debris
RI PT
Removal missions. Acta Astronautica, 131 (2017) pp. 10-17. DOI:10.1016/j.actaastro.2016.11.012 [24] Segal, Shai, Pini Gurfil, and Kamran Shahid. In-orbit tracking of resident space objects: A comparison of monocular and stereoscopic vision. IEEE Transactions on Aerospace and Electronic Systems 50
SC
(No.1) (2014) pp. 676-688. DOI: 10.1109/TAES.2013.120006
M AN U
[25] Maccone, Claudio. Advantages of Karhunen–Loève transform over fast Fourier transform for planetary radar and space debris detection. Acta Astronautica 60 (No.8-9) (2007) pp.775-779. DOI: 10.1016/j.actaastro.2006.08.015
[26] Thorn, E., Korkiakoski, V., Grosse, D., Bennet, F., Rigaut, F., d'Orgeville, & Smith, C..
TE D
Stereo-SCIDAR system for improvement of adaptive optics space debris-tracking. Advanced Maui Optical and Space Surveillance Technologies Conference, 2017 [27] Gamage S, Cheng Y. CubeSat Detection Using Convolutional Neural Networks. 2018 Space Flight
EP
Mechanics Meeting. 2018: 2234.
AC C
[28] Julier, Simon J. The scaled unscented transformation. Proceedings of the American Control Conference, (Vol. 6) (2002), pp. 4555-4559 [29] Lebedev, V., and Laikov, D., Quadrature Formula for the Sphere of the 131st Algebraic Order of Accuracy, Doklady Mathematics, 59 (No. 3) (1999) pp. 477–481 [30] Mehrholz D, Leushacke L, Flury W, et al. Detecting, tracking and imaging space debris. ESA Bulletin 109 (0376-4265) (2002) pp. 128-134. [31] Ballard, Dana H. Generalizing the Hough transform to detect arbitrary shapes. Readings in computer
- 28 -
ACCEPTED MANUSCRIPT vision. 1987. pp. 714-725. DOI: 10.1016/0031-3203(81)90009-1 [32] Canny, John. A computational approach to edge detection. Readings in Computer Vision. 1987. pp. 184-203. DOI: 10.1109/TPAMI.1986.4767851
RI PT
[33] Quinn, M. N., Jukna, V., Ebisuzaki, T., Dicaire, I., Soulard, R., Summerer, L., Mourou, G. (2015). Space-based application of the CAN laser to LIDAR and orbital debris remediation. The European
AC C
EP
TE D
M AN U
SC
Physical Journal Special Topics, 224(13), 2645-2655. DOI: 10.1140/epjst/e2015-02577-5
- 29 -
ACCEPTED MANUSCRIPT 1. Optical-cues is efficiently utilized to estimate the object information. 2. The GMM models are employed to form the analytical form of collision analysis.
AC C
EP
TE D
M AN U
SC
RI PT
3. MC simulation and an experimental test verify the effectiveness of this approach.