# Robust Real-time RGB-D Visual Odometry in Dynamic Environments
via Rigid Motion Model

Sangil Lee, Clark Youngdong Son, and H. Jin Kim Sangil Lee, Clark Youngdong
Son, and H. Jin Kim are with the Department of Mechanical and Aerospace
Engineering, Seoul National University, Seoul, 08826, Korea, Republic of
{sangil07, clark.y.d.son<EMAIL_ADDRESS>
###### Abstract

In the paper, we propose a robust real-time visual odometry in dynamic
environments via rigid-motion model updated by scene flow. The proposed
algorithm consists of spatial motion segmentation and temporal motion
tracking. The spatial segmentation first generates several motion hypotheses
by using a grid-based scene flow and clusters the extracted motion hypotheses,
separating objects that move independently of one another. Further, we use a
dual-mode motion model to consistently distinguish between the static and
dynamic parts in the temporal motion tracking stage. Finally, the proposed
algorithm estimates the pose of a camera by taking advantage of the region
classified as static parts. In order to evaluate the performance of visual
odometry under the existence of dynamic rigid objects, we use self-collected
dataset containing RGB-D images and motion capture data for ground-truth. We
compare our algorithm with state-of-the-art visual odometry algorithms. The
validation results suggest that the proposed algorithm can estimate the pose
of a camera robustly and accurately in dynamic environments.

##  I INTRODUCTION

Visual odometry is a fundamental process of recognizing the pose of the camera
itself using video input[1, 2]. Although various visual odometry algorithms
already show satisfactory performances in well-conditioned environments and
well-defined datasets such as TUM[3] and KITTI[4], most of them assume that
the world the camera is looking at is stationary, thus making it possible to
estimate the pose of the camera by virtue of the motion of the images taken.
However, most real environments involve dynamic situations such as residential
roads, crowded places, other robots for cooperation. Although some visual
odometry algorithms [5] which utilize the principle of RANSAC can regard
pixels whose motion is disparate as outliers, they have a limit that non-
stationary objects should occupy small areas in the image plane. Thus, they
cannot be employed to estimate the motion of a camera in dynamic environments
including large non-stationary objects.

In this paper, we aim for a real-time robust visual odometry by separating
stationary parts from the image via motion model update. However, the existing
motion segmentation algorithms have expensive computation loads[7] or
constraints on the shape[8] or the number[9] of objects. Here, we design a
fast motion segmentation algorithm with an adequate performance applicable for
real-time visual odometry algorithm.

Figure 1: The 3D trajectory on the uav-flight-circular sequence.

In order to differentiate between non-stationary parts and stationary
background, we utilize scene flow vectors which are distributed uniformly in
the image. Particularly, we choose grid-based scene flow to take advantage of
both dense and sparse methods; a dense flow that calculates temporary motions
for all pixels provides a high resolution, meanwhile, a sparse flow that
calculates temporary motions for distinctive features has a lighter
computational load. Then, motion segmentation is performed to differentiate
between non-stationary parts and stationary background with grid-based
temporary motions, and the pose of a camera is estimated using static parts.

Overall, our algorithm can estimate the ego-motion robustly and accurately in
highly dynamic environments while separating non-stationary parts from an
image with no prior information such as the shape[8], the number of dynamic
objects[9], or the movement of objects[10]. Fig. 1 shows the performance of
the proposed algorithm on the dataset collected from a multirotor flight.
Moreover, the proposed algorithm shows significantly low runtime of average 19
ms at VGA resolution, thus it can be applied to real-time tasks.

Figure 2: A schematic diagram of the proposed algorithm.

###  I-A Related Work

Most existing visual odometry algorithms [11, 12, 13, mur2017orb, 15] assume
stationary environments. However, in real applications, there are often a
number of non-stationary objects such as vehicles. To deal with such
environments, some research has attempted to improve robustness against
dynamic situations. Such efforts can be categorized into two types:
dichotomous and model-based.

The first type of research has tried to exclude a region which has a different
motion from the major movement in a similar way to RANSAC as mentioned above.
A. Dib et al.[5] utilizes RANSAC for direct visual odometry in dynamic
environments. They minimize photometric shift over six random patches in an
image, unlike the naïve direct visual odometry whose optimization process
covers the whole image. In [6], they first categorize patches by depth value
and subtract a category which has a quite different motion from a majority
movement of the camera by comparing standard deviation of motion vectors,
including or excluding that category. They still assume that the stationary
background occupies a large area on the image so that the camera movement is
estimated via the background.

The second type of research has exploited statistical models to classify
motions belonging to independent objects. BaMVO [10] estimates background by
choosing pixels whose depth does not change unusually. It is effective in
situations where a dynamic object such as a pedestrian moves in parallel with
the principal axis of a camera. However, the performance can degenerate when
the dynamic object moves perpendicular to the principal axis while reducing
the depth transition. Joint-VO-SF[17] formulates a minimization problem with
3D points from RGB-D images. They estimate camera pose and scene flow
accurately, but the average runtime is 80 ms on an i7 multi-core CPU at QVGA
resolution, which is slightly slow for real-time implementation with current
on-board computers. StaticFusion[scona2018staticfusion] is the enhanced
version of Joint-VO-SF. It adds dense 3D modeling of only the static parts of
the environment and reduces the overall drift though frame-to-model alignment.

Among the motion segmentation methods, H. Jung et al.[9] propose randomized
voting to extract independent motions using epi-polar constraints with the
average computational time of 300ms per frame. However, the number of moving
objects should be known. In [21], MCD5.8ms has an advantage in that it detects
a moving object using a dual-mode model with a low computational load while
showing the execution time per frame of 5.8 ms. Since they calculate the
homography using naïve visual odometry, however, the performance can
deteriorate when a moving object occupies more than half of an image.

###  I-B Contributions

In this paper, we focus on the environment where there are rigid moving
objects. Furthermore, we propose a dichotomy-type algorithm to prevent the
pose estimates from being polluted due to pixels with vague states between
static and dynamic.

Our main contributions can be summarized as follows:

  1. 1.

We propose a real-time robust visual odometry algorithm in dynamic
environments. It estimates the motion of both stationary and non-stationary
parts robustly.

  2. 2.

We design a dual-mode motion model for rigid motion segmentation. It
distinguishes between the background and moving objects with no prior
information.

##  II BACKGROUND

###  II-A Rigid Transformation

The proposed algorithm is based upon the property that 3D motions of 3D points
belonging to the same rigid object spatially have the same rigid motion
temporally. We calculate the motion, $\mathbf{H}$, by the least-squares rigid
transformation[23]. Moreover, for the evaluation of the motion, we can define
the rigid transformation error as below:

$\displaystyle\mathbf{E}(\mathbf{H},\tilde{\mathbf{X}}^{(j)},\tilde{\mathbf{X}}^{(k)})\coloneqq\textrm{diag}((\mathbf{H}\tilde{\mathbf{X}}^{(j)}-\tilde{\mathbf{X}}^{(k)})^{T}(\mathbf{H}\tilde{\mathbf{X}}^{(j)}-\tilde{\mathbf{X}}^{(k)})),$
(1)

where
$\tilde{\mathbf{X}}^{(l)}=[\tilde{\mathbf{x}}_{1}^{(l)},\tilde{\mathbf{x}}_{2}^{(l)},\ldots,\tilde{\mathbf{x}}_{n}^{(l)}]\in\mathbb{R}^{4\times
n}$ is a 3D point set, and
$\tilde{\mathbf{x}}_{*}^{(l)}=[x,y,z,1]^{T}\in\mathbb{R}^{4\times 1}$ is a 3D
point in the $l$-th frame represented by homogeneous coordinates.

###  II-B Advanced Grid-based Scene Flow

In order to extract a 3D motion of pixels, we utilize Lucas-Kanade optical
flow. For expanding the dimension of optical flow, we use depth changes from
the adjacent depth images to obtain depth-directional flow. However, since
there can be invalid depth pixels, we improve the quality of depth values
through preprocessing. We fill pixels whose depth value is invalid with
similar value in the vicinity in case of narrow holes and eliminate pixels
whose depth value is abnormally high or low. Afterward, we interpolate the
depth of tracked points.

##  III MOTION SPATIAL SEGMENTATION

In the motion spatial segmentation procedure, we first fetch a total of $n$
grid-based scene flow vectors, $\mathbf{V}_{i}$, with grid cell size,
$w_{\rm{grid}}$, and derive segments through their motions. In the procedure,
the number of points $m$ and heuristic threshold $th_{\rm{inlier}}$ are used
for generating motion hypothesis. To be specific, $m$ number of points are
used to generate motion hypothesis, and the threshold value is used to decide
whether motions of the $m$ points belong to the same movement. The motion
spatial segmentation procedure is divided into three steps: motion hypothesis
search, refinement, and clustering. Especially, the searching and the refining
processes are executed iteratively until a total of $h_{\textrm{max}}$
hypotheses are found or the mean of entropies $\mathbf{S}$ is decreased enough
and saturated. Fig. 3 describes the motion spatial segmentation process.

Figure 3: The description of the motion segmentation on the place-items
sequence.

###  III-A Motion Hypothesis Search

In the search process, a point is randomly selected by weighting with entropy,
$\mathbf{S}$, whose $i$-th value is defined as follows:

$S_{i}=1-\min_{1\leq h\leq h_{\rm{max}}}{\exp(-\lambda\hskip
1.42262pt\mathcal{E}_{i}-\delta)},~{}i=1\ldots n,$ (2)

where

$\displaystyle\mathcal{E}_{i}=\left[E_{i}(\mathbf{H}_{1},\tilde{\mathbf{X}}^{(j)},\tilde{\mathbf{X}}^{(k)}),\ldots,E_{i}(\mathbf{H}_{h_{\rm{max}}},\tilde{\mathbf{X}}^{(j)},\tilde{\mathbf{X}}^{(k)})\right].$
(3)

In the above equations, $E_{i}\left(\cdot\right)$ is the $i$-th element in the
output of the $\mathbf{E}\left(\cdot\right)$ function defined in Eq. 1. Also,
$\lambda$ and $\delta$ are parameters for defining the characteristic of the
entropy.

As will be discussed in more detail in the following section, the entropy
measures how well the hypothesis is generated in the corresponding pixels. A
pixel with the high entropy tends to belong to an inappropriate motion
hypothesis, which results in a high probability to be chosen, consequently
improving the convergence rate. Then we select $(m-1)$-points randomly near
the chosen 1-point. Finally, the rigid motion hypothesis $\mathbf{H}$ can be
estimated from the selected $m$-points[23].

###  III-B Motion Hypothesis Refinement and Clustering

Motion hypothesis refinement aims to calculate precise rigid transformation
and find multiple regions with the same rigid motion. We first calculate rigid
transformation error using Eq. 1 for all of the $n$ grid-based scene flow
vectors. This process is designed to consider the cases where a part of a
rigid object could appear in multiple regions on an image at the same time, in
which case, we regard the multiple regions as belonging to the same object.
Then we re-estimate a refined motion hypothesis using the increased
$N$-points. For more details, please refer to [lee2017real].

In order to find distinct motions from the motion hypothesis set
$\mathcal{H}$, we use the existing clustering algorithm, in particular,
density-based spatial clustering of applications with noise (DBSCAN), since
DBSCAN does not require the number of clusters and is robust to noise. In the
process, hypotheses are refined as the algorithm reorganizes
$h_{\textrm{max}}$ hypotheses into $g$ distinct motions, which is represented
by vector $\mathbf{G}\in\mathbb{N}^{n\times 1}$ whose elements are distributed
from 1 to $g$ in an natural number. Thus, a $\mathbf{G}^{(k)}$ is the
segmentation results of the $k$-th frame at the motion spatial segmentation
stage.

##  IV MOTION SEGMENT TRACKING

In addition to distinguishing objects which move independently from one
another in a frame, we propose a motion segment tracking algorithm to update
the region of dynamic objects persistently. Since the above motion spatial
segmentation algorithm runs in a similar way as the naïve segmentation
technique based on two frames, it has no idea which is which at different
times. Thus, we propose the motion segment tracking algorithm following the
motion spatial segmentation. The following subsection provides a segment
matching algorithm and a probabilistic approach based on the dual-mode simple
Gaussian model [21].

###  IV-A Segment Matching

The segment matching algorithm calculates a kind of correlation coefficient
and finds corresponding segment pairs $\hat{\mathcal{P}}$ between different
frames. The correlation coefficient and segment pairs, $\mathcal{P}^{(k,l)}$,
between the current $k$-th and the previous $l$-th frames are defined as
follows:

$\displaystyle C_{ij}^{(k,l)}$
$\displaystyle=\frac{N_{i,j}^{(k,l)}}{\sqrt{N_{i}^{(k)}N_{j}^{(l)}}},$ (4)
$\displaystyle\mathcal{P}^{(k,l)}$
$\displaystyle=\left\\{\mathcal{P}^{(k,l)}_{i}\right\\}=\left\\{\left(\mathcal{P}^{(k,l)}_{i,\textrm{left}},\mathcal{P}^{(k,l)}_{i,\textrm{right}}\right)\right\\}$
(5)
$\displaystyle=\left\\{\left(i,~{}\underset{j}{\textrm{argmax}}C_{ij}^{(k,l)}\right)\right\\},$

where $N_{i}^{(k)}$ is the number of grid-based points belonging to the $i$-th
segment in the $k$-th frame, and $N_{i,j}^{(k,l)}$ is the number of grid-based
points belonging to both the $i$-th and $j$-th segments in the $k$-th and
$l$-th frames, respectively. Also, the corresponding segment pairs,
$\hat{\mathcal{P}}^{(k)}$, whose correlation score is the maximum, is defined
as follows:

$\hat{\mathcal{P}}^{(k)}=\mathcal{P}^{(k,l_{m})},$ (6)

where $l_{m}$ is

$l_{m}=\underset{l}{\textrm{argmax}}(score^{(k,l)}),$ (7)

a score between the $(k,l)$-th frames is

$score^{(k,l)}=\frac{\sum_{i=1}^{g}{\left(N_{i}^{(k)}\times\underset{j}{\max}{C}_{ij}^{(k,l)}\right)}}{\sum_{i=1}^{g}{N_{i}^{(k)}}},$
(8)

and $g$ is the number of distinct motions. And then, we rearrange $\mathbf{G}$
to $\hat{\mathbf{G}}$ through $\hat{\mathcal{P}}^{(k)}$ as follows:

$\hat{G}^{(k)}_{i}\equiv\hat{G}^{(k)}_{\hat{\mathcal{P}}^{(k)}_{i,\textrm{left}}}=\hat{G}^{(l_{m})}_{\hat{\mathcal{P}}^{(k)}_{i,\textrm{right}}},$
(9)

so that we track the object that appeared previously or add a new segment on
the unmatched pairs of the current frame. From now, we will call the value of
$\hat{\mathbf{G}}$ as a label in the following paragraphs, while we have
called the value of $\mathbf{G}$ as a segment.

###  IV-B Dual-mode Motion Model

Since the RGB-D camera has a limited depth range, it may fail to calculate
scene flow vectors in some grid cells. Thus, we use a discrete statistical
model which has a probability vector and an age as properties in order to
track the static and dynamic parts in the image sequences persistently under
the assumption that the static and dynamic elements do not appear or disappear
abruptly on a frame by frame. The $i$-th element of the probability vector
means the likelihood of corresponding pixel belonging to the $i$-th label. As
shown in the right dashed box of Fig. 2, the algorithm first compensates the
model through the previously estimated ego-motion to update it with the
measurement corresponding to the identical 3D point. Then, we update the
probability vector and age of the model based on certain criteria. Finally,
labels can be selected as indices which indicate the maximum value in the
probability vector. In this paper, each of $n$ grid cells has two models, i.e.
apparent and candidate models. The apparent model indicates the estimated
label currently, and the candidate model implies a hidden label which can
appear later in the apparent model. The candidate model is designed to
recognize again an object as static parts when the moving object stops.

####  IV-B1 Model Update

We denote the probability vector of a grid cell $i$ as
$P_{i}^{(k)}\in\mathbb{R}^{n_{\textrm{obj}}\times 1}$ in the $k$-th frame, and
the age of the grid cell in the $k$-th frame as $\alpha_{i}^{(k)}$, where
$n_{\rm{obj}}=15$ is the designated maximum number of identified objects in a
frame. Then, the probability vector and age are updated as follows:

$\displaystyle P_{i}^{(k+1)}$
$\displaystyle=\frac{\tilde{\alpha}_{i}^{(k)}}{\tilde{\alpha}_{i}^{(k)}+1}\tilde{P}_{i}^{(k)}+\frac{1}{\tilde{\alpha}_{i}^{(k)}+1}\check{G}_{i}^{(k+1)}$
(10) $\displaystyle\alpha_{i}^{(k+1)}$
$\displaystyle=\tilde{\alpha}_{i}^{(k)}+1,$ (11)

where $\check{G}^{(k)}\in\mathbb{R}^{n_{\textrm{obj}}\times 1}$ is a binary-
valued vector whose $\hat{G}_{i}^{(k)}$-th element only is one and the others
are zeros, and $\hat{G}_{i}^{(k)}$ is the result of segment matching algorithm
in the $k$-th frame. $\tilde{P}_{i}^{(k)}$, $\tilde{\alpha}_{i}^{(k)}$ are the
compensated parameters of the dual-mode motion model which will be discussed
in Section IV-B2.

Contrary to [21] which measures pixel intensity to update models, we take a
temporarily matched label, which is a result of segment matching, as
measurements. However, our measurements might be incorrect, if scene flow may
too fast to find an appropriate match and the properly-matched pair could not
be found. Thus, we make some modification for the update of both models. The
probability vector, ${}^{A}P$, and age, ${}^{A}\alpha$, of the apparent model
are updated as follows:

$\displaystyle{}^{A}P_{i}^{(k+1)}$
$\displaystyle=\begin{cases}{}^{A}\tilde{P}_{i}^{(k)},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\text{if
$\check{G}_{i}^{(k+1)}\neq\check{G}_{i}^{(k)}$}\vspace{3mm}\\\
\frac{1}{{}^{A}\tilde{\alpha}_{i}^{(k)}+1}({}^{A}\tilde{\alpha}_{i}^{(k)}{}^{A}\tilde{P}_{i}^{(k)}+\check{G}_{i}^{(k+1)}),~{}~{}\text{otherwise}\end{cases}\vspace{1mm}$
$\displaystyle{}^{A}\alpha_{i}^{(k+1)}$
$\displaystyle=\begin{cases}{}^{A}\tilde{\alpha}_{i}^{(k)},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\text{if
$\check{G}_{i}^{(k+1)}\neq\check{G}_{i}^{(k)}$}\vspace{2mm}\\\
bm\min{({}^{A}\tilde{\alpha}_{i}^{(k)}+1,\alpha_{\rm{max}})},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\text{otherwise}\end{cases}$

and the probability vector, ${}^{C}P$, and age ,${}^{C}\alpha$, of the
candidate model are updated as follows:

$\displaystyle{}^{C}P_{i}^{(k+1)}$
$\displaystyle=\begin{cases}{}^{C}\tilde{P}_{i}^{(k)},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\text{if
$\check{G}_{i}^{(k+1)}=\check{G}_{i}^{(k)}$}\vspace{3mm}\\\
\frac{1}{{}^{C}\tilde{\alpha}_{i}^{(k)}+1}({}^{C}\tilde{\alpha}_{i}^{(k)}{}^{C}\tilde{P}_{i}^{(k)}+\check{G}_{i}^{(k+1)}),~{}~{}\text{otherwise}\end{cases}\vspace{1mm}$
$\displaystyle{}^{C}\alpha_{i}^{(k+1)}$
$\displaystyle=\begin{cases}{}^{C}\tilde{\alpha}_{i}^{(k)},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\text{if
$\check{G}_{i}^{(k+1)}=\check{G}_{i}^{(k)}$}\vspace{2mm}\\\
\min{({}^{C}\tilde{\alpha}_{i}^{(k)}+1,\alpha_{\rm{max}})},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\text{otherwise}\\\
\end{cases}$

Also, the candidate and apparent models are swapped if the age of the
candidate model reaches the maximum age, $\alpha_{\rm{max}}$, or is larger
than that of the corresponding apparent model. Because we treat a foreground
object as a static element when the object stops. Thus, after the previously
moving object stops, the corresponding apparent model is not updated whereas
the candidate model is updated. Consequently, the age of only the candidate
model increases and both models will be swapped with each other when the age
of the candidate model becomes saturated or larger than that of the apparent
model.

A label which is the output of the motion segment tracking is obtained from
the probability vector of the apparent model with several criteria. The label
is updated when the corresponding apparent model is initialized or updated. By
doing so, we prevent the algorithm from prejudging the label of an unobserved
or unmeasured grid cell while maintaining the previous label of the grid cell.
Finally, the label is obtained from the indices that represent the maximum
value in the probability vector of the apparent model.

####  IV-B2 Compensate Model

In the previous section, we use the probability vector for updating and
determining the label of a grid cell. These processes for the motion model
assume that each grid cell represents a fixed point in the world coordinate
consistently. However, in the case of a non-stationary camera, the result of
the segment matching following the spatial segmentation cannot be used
directly for updating the apparent and candidate model. In order to update the
motion model, therefore, we compensate the model by warping.

|  |  |
---|---|---|---
Figure 4: Variations of the parameters. The relative pose error are denoted as
magenta boxplot with 3-$\sigma$ whiskers. Average computational time is
represented as a black solid line with respect to the right y-axis.

Since each model has simple parameters such as probability vector and age, we
use area-weighted interpolation. The current model is proportionally
compensated with scene flow vectors in Section II-B. For grid cells that have
valid scene flow vectors, these motion models are compensated individually on
the two-dimensional image plane. We denote the set of grid cells overlapping
with the current grid cell $i$ in the $(k+1)$-th frame as
$\mathcal{A}_{i}^{(k+1)}$, weights for interpolation as $\omega_{j}$, and the
region of overlap between the current grid cell $i$ and the previous
corresponding grid cell as $R_{j}$ where $j$ is the index of element in
$\mathcal{A}_{i}^{(k+1)}$. Then, the compensated probability vector and age
are obtained as follows:

$\small{\tilde{P}_{i}^{(k)}=\sum_{j\in\mathcal{A}_{i}^{(k+1)}}{\omega_{j}P_{j}^{(k)}},~{}\tilde{\alpha}_{i}^{(k)}=\sum_{j\in\mathcal{A}_{i}^{(k+1)}}{\omega_{j}\alpha_{j}^{(k)}},}$
(12)

where $\omega_{j}$ are set to be proportional to $R_{j}$ and normalized to
$\sum_{j}{\omega_{j}}=1$.

##  V EVALUATION RESULTS

This section is divided into three parts: analysis of parameters, description
of the dataset, and quantitative evaluation of the odometry. In the
validation, we used RGB-D images of 640$\times$480 pixels. Our algorithm is
run on an Intel i7-7500U CPU at 2.7 GHz and Ubuntu 16.04 LTS.

###  V-A Parameters Analysis

Here, we discuss the value of parameters introduced in the paper. These
parameters are obtained from indoor environments including our dataset.
Therefore, they need to be tuned for a new environment including outdoor
environments.

Fig. 4 shows variations of the performance and computational loads according
to parameters (from left): grid spacing, i.e. size of grid cell
$w_{\rm{grid}}$, the maximum number of motion hypothesis $h_{\rm{max}}$, rigid
transformation error threshold in motion hypothesis generation
$th_{\rm{inlier}}$, the saturated age of motion model $\alpha_{\rm{max}}$. In
each plot, one parameter varies while others remain the designated values
given in the last paragraph of this subsection. As $w_{\rm{grid}}$ increases,
the spatial density of the scene flow reduces, thus the computation time
decreases reciprocally and Relative Pose Error (RPE) increases slightly.
Therefore, this is advantageous in the presence of a small dynamic object.
Also, if we choose a small value of $h_{\textrm{max}}$, it is difficult to
extract all meaningful motions due to few tries in motion hypothesis
searching, whereas a large value of $h_{\textrm{max}}$ can produce an
erroneous label with a high computational load. Thus, if there are a lot of
dynamic objects, it is recommended that the value of $h_{\textrm{max}}$ be not
small. $th_{\rm{inlier}}$ is related to hypothesis searching; a low threshold
makes the algorithm to extract motion hypothesis strictly. A small value of
$\alpha_{\rm{max}}$ is vulnerable to erroneous label, whereas a high value
makes the algorithm too insensitive not to notice that the previously moving
object has stopped.

From the above boxplot analysis, the parameters were set to:
$w_{\rm{grid}}=16$, $h_{\rm{max}}=20$, $th_{\rm{inlier}}=3\times 10^{-5}$,
$\alpha_{\rm{max}}=5$, which are well-tuned values for our indoor dataset. The
parameters are kept the same over the dataset for consistency.

Rest of the parameters are set as follows:

  * •

The number of points and the searching radius in Section III-A are set to
$m=7$ and $r_{\rm{search}}=2$ pixels.

  * •

The characteristic parameters and the minimum threshold of search weighting in
Section III-B are set to $\lambda=10^{3}$ and $\delta=10^{-2}$, heuristically.

  * •

DBSCAN parameters are set to $p_{\rm{min}}=1$, $\epsilon=0.005$.

Figure 5: The validation environment for the uav-flight-circular sequence.

###  V-B Dataset

To the best of our knowledge, there is no well-known RGB-D dataset captured in
a situation where a rigid moving object appears. TUM dataset[3] contains
moving people but human body is not a rigid object and KITTI[4] does not
provide depth images. Therefore, in order to evaluate the proposed algorithm,
we collected a dataset using ASUS Xtion RGB-D camera and Vicon motion capture
system for ground truth (Fig. 5). Each sequence is composed of 16-bit depth
images and 8-bit color images with a size of 640$\times$480\. There are two
kinds of sequences depending on whether they are captured from the static or
non-stationary camera as shown in Table I. In the case of sequences using the
static camera, we regard the origin of the world coordinate as the true
position of the camera. On the other hand, for the non-stationary camera, we
validate the performance of our visual odometry algorithm with Vicon
measurements. The detailed description and download link of the dataset are
available online at:

http://sangillee.com/_pages/icsl-de-dataset

TABLE I: Evaluation of visual odometry algorithms on our dataset.

Environment Sequence Relative Pose Error [m/s] Proposed DVO ORB-VO BaMVO
Joint-VO-SF StaticFusion Static camera & Dynamic environment one-object-static
0.0053 0.6002 0.0008 0.0017 0.0053 0.2042 two-object-static 0.0221 0.1673
0.0725 0.0063 0.0172 0.1603 place-items 0.0039 0.0550 0.0017 0.0462 0.0344
0.0085 Dynamic camera & Dynamic environment fast-object 0.0249 0.4240 0.3428
0.2022 0.0405 0.1726 slow-object 0.0600 0.1962 0.1772 0.1248 0.0724 0.1311
close-approach 0.0469 0.2101 0.0931 0.0992 0.0707 0.1360 leading-pioneer
0.0996 0.1449 0.0679 0.2790 0.1444 0.5503 uav-flight-static 0.0231 0.3880
0.2389 0.2342 0.0279 0.3812 uav-flight-circular 0.0290 0.2586 0.1512 0.4034
0.0737 0.3109

TABLE II: Evaluation of visual odometry algorithms on TUM dataset.

Environment Sequence Relative Pose Error [m/s] Proposed DVO ORB-VO BaMVO
Joint-VO-SF StaticFusion Dynamic camera & Static environment fr1/xyz 0.0266
0.1379 0.0139 0.1763 0.0174 0.0549 fr1/rpy 0.0420 0.0406 0.0303 0.1858 0.0384
0.0889 fr1/desk 0.0535 0.0675 0.0409 0.2653 0.0291 0.1718 fr1/floor 0.0306
0.1172 0.0129 0.1247 0.0266 0.4150 Dynamic camera & Dynamic environment
fr3/walking_static 0.0374 0.2022 0.1669 0.0939 0.0709 0.0146 fr3/walking_xyz
0.2358 0.2980 0.2434 0.1887 0.2064 0.0913 fr3/sitting_xyz 0.0909 0.0367 0.0110
0.0442 0.0444 0.0325

###  V-C Visual Odometry

For the quantitative comparison between the proposed algorithm and the current
state-of-the-art visual odometry algorithms, we use RPE with a one second
drift as proposed in [3] since RPE is well-suited for evaluating the drift of
visual odometry. Open-source algorithms were executed with default settings.

Table I and Table II show the evaluation results for each sequence in our
dataset and TUM dataset, respectively. The algorithm with the best result in
each sequence is shown in bold, and the algorithm with an error greater than
0.1 m/s is in red. The proposed algorithm is compared with a well-known or
state-of-the-art visual odometry algorithms such as DVO[12], ORB-
VO[mur2017orb], BaMVO[10], Joint-VO-SF[17], and
StaticFusion[scona2018staticfusion]. Particularly, BaMVO, Joint-VO-SF, and
StaticFusion were designed to be robust in dynamic environments likewise ours.
In order to verify the performance as visual odometry, we evaluate a modified
ORB-SLAM2, which is unable to detect and correct loop closure. We refer to
this modified ORB-SLAM2 as ORB-VO.

On the one-object-static sequence, algorithms show outstanding performance
except for DVO. Because DVO is a direct dense method, it performs optimization
across all pixels, thereby its performance is seriously influenced by dynamic
elements. Next, when there appear two moving objects, the performance of our
algorithm, BaMVO, and Joint-VO-SF is superior to the other algorithms. The
performance of BaMVO is degraded when the object moves perpendicular to the
principal axis in the uav-flight-static, uav-flight-circular sequences. ORB-VO
tends to track features of a dynamic object when the object is observed for a
long time in the fast-object sequence even if there are a lot of features on
static backgrounds. Besides, Joint-VO-SF and StaticFusion fail optimization
sometimes and does not perform well on some sequences. On the other hands, our
algorithm shows a balanced and sufficient performance over all tested
sequences. In the static environment sequences of TUM dataset, ORB-VO shows
reliable performance superior to the other algorithms. Since the dynamic
environment sequences of TUM dataset include the non-rigid human body, the
proposed algorithm does not have an advantage over the existing algorithms,
but it does not have a bad performance either. Please see the distribution of
red texts for validating the reliability of the proposed algorithm.

TABLE III: Enhanced visual odometry algorithms for robustness.

Environment Relative Pose Error [m/s] Proposed $\times$ ORB-VO ORB-VO one-
object-static 0.0009 0.0008 two-object-static 0.0242 0.0725 place-items 0.0053
0.0017 fast-object 0.0232 0.3428 slow-object 0.1553 0.1772 close-approach
0.0728 0.0931 leading-pioneer $\times$ 0.0679 uav-flight-static 0.0078 0.2389
uav-flight-circular 0.0714 0.1512

---
Figure 6: Some of the sequence images and the segmentation results that the
proposed algorithm provides. Filled circles mean grid cell that has accurate
scene flow and valid label. (Recommended to print out in color)

Moreover, since the motion segmentation and the estimation parts of our
algorithm are not strongly coupled with each other, it is possible to
incorporate the proposed motion segmentation into the existing visual odometry
algorithms in order to improve their robustness in dynamic environments. As
shown in Table III, we validate the combination of ORB-VO and our motion
segmentation. We can see that the proposed motion segmentation enhances the
odometry performance compared with the original version. The leading-pioneer
sequence contains only a small amount of valid depth for the entire time, so
the feature points are not extracted enough causing failure. One way to
improve the performance in this case is to convert from the grid-base to the
dense method, and the accurate dense optical flow must be preceded first.

###  V-D Runtime Comparison

The median runtimes of the compared algorithms are:

  *   * •

Proposed: 53 ms1,
267 ms2

  * •

DVO2: 1.032 sec

  * •

ORB-VO1: 33 ms

  * •

BaMVO1: 278 ms

  * •

Joint-VO-SF1: 92 ms

  * •

StaticFusion1: 1.221 sec

They are evaluated with 1 C++ implementation on a laptop computer as described
in the first paragraph of Section V or 2 Matlab on a desktop computer (Intel
i5-3770 at 3.4 GHz). Also, we make algorithms fetch RGB-D images of the same
VGA size for fair comparison.

##  VI CONCLUSIONS

In this paper, we proposed a real-time robust visual odometry algorithm via
rigid motion segmentation using grid-based scene flow. The proposed algorithm
is considerably more robust and accurate than the state-of-the-art visual
odometry algorithms. For robustness, the proposed spatial motion segmentation
uses scene flow to generate and search distinct motions with no prior
information such as the shape or number of objects. Besides, temporal
segmentation initializes and updates a dual-mode motion model of the grid cell
so that our algorithm differentiates stationary background and dynamic objects
robustly. Finally, the ego-motion is estimated by the use of scene flow vector
fields belonging to the stationary background. An additional benefit of the
proposed algorithm is that it can be combined with the existing visual
odometry algorithms to improve their robustness in dynamic environments.
Furthermore, the proposed approach can estimate the motion of moving objects
unlike the other existing algorithms, so it can be employed as a part of an
efficient dynamic obstacle avoidance algorithm for an autonomous robot by
using the kinematic information of moving objects.

## References

  * [1] D. Nistér, O. Naroditsky, and J. Bergen, “Visual odometry,” in _Computer Vision and Pattern Recognition, 2004. CVPR 2004. Proceedings of the 2004 IEEE Computer Society Conference on_ , vol. 1. IEEE, 2004, pp. I–652.
  * [2] D. Scaramuzza and F. Fraundorfer, “Visual odometry [tutorial],” _IEEE Robotics & Automation Magazine_, vol. 18, no. 4, pp. 80–92, 2011.
  * [3] J. Sturm, N. Engelhard, F. Endres, W. Burgard, and D. Cremers, “A benchmark for the evaluation of rgb-d slam systems,” in _Intelligent Robots and Systems (IROS), 2012 IEEE/RSJ International Conference on_. IEEE, 2012, pp. 573–580.
  * [4] A. Geiger, P. Lenz, C. Stiller, and R. Urtasun, “Vision meets robotics: The kitti dataset,” _The International Journal of Robotics Research_ , vol. 32, no. 11, pp. 1231–1237, 2013.
  * [5] A. Dib and F. Charpillet, “Robust dense visual odometry for rgb-d cameras in a dynamic environment,” in _Advanced Robotics (ICAR), 2015 International Conference on_. IEEE, 2015, pp. 1–7.
  * [6] S.-J. Jung, J.-B. Song, and S.-C. Kang, “Stereo vision-based visual odometry using robust visual feature in dynamic environment,” _The Journal of Korea Robotics Society_ , vol. 3, no. 4, pp. 263–269, 2008.
  * [7] R. Sabzevari and D. Scaramuzza, “Monocular simultaneous multi-body motion segmentation and reconstruction from perspective views,” in _Robotics and Automation (ICRA), 2014 IEEE International Conference on_. IEEE, 2014, pp. 23–30.
  * [8] Y.-H. Tsai, M.-H. Yang, and M. J. Black, “Video segmentation via object flow,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2016, pp. 3899–3908.
  * [9] H. Jung, J. Ju, and J. Kim, “Rigid motion segmentation using randomized voting,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2014, pp. 1210–1217.
  * [10] D.-H. Kim and J.-H. Kim, “Effective background model-based rgb-d dense visual odometry in a dynamic environment,” _IEEE Transactions on Robotics_ , vol. 32, no. 6, pp. 1565–1573, 2016.
  * [11] F. Steinbrücker, J. Sturm, and D. Cremers, “Real-time visual odometry from dense rgb-d images,” in _Computer Vision Workshops (ICCV Workshops), 2011 IEEE International Conference on_. IEEE, 2011, pp. 719–722.
  * [12] C. Kerl, J. Sturm, and D. Cremers, “Dense visual slam for rgb-d cameras,” in _Intelligent Robots and Systems (IROS), 2013 IEEE/RSJ International Conference on_. IEEE, 2013, pp. 2100–2106.
  * [13] C. Forster, Z. Zhang, M. Gassner, M. Werlberger, and D. Scaramuzza, “Svo: Semidirect visual odometry for monocular and multicamera systems,” _IEEE Transactions on Robotics_ , vol. 33, no. 2, pp. 249–265, 2017.
  * [14] R. Mur-Artal and J. D. Tardos, “Orb-slam2: an open-source slam system for monocular, stereo and rgb-d cameras,” _arXiv preprint arXiv:1610.06475_ , 2016.
  * [15] J. Engel, V. Koltun, and D. Cremers, “Direct sparse odometry,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 2017.
  * [16] B. Kitt, F. Moosmann, and C. Stiller, “Moving on to dynamic environments: Visual odometry using feature classification,” in _Intelligent Robots and Systems (IROS), 2010 IEEE/RSJ International Conference on_. IEEE, 2010, pp. 5551–5556.
  * [17] M. Jaimez, C. Kerl, J. Gonzalez-Jimenez, and D. Cremers, “Fast odometry and scene flow from rgb-d cameras based on geometric clustering,” in _Proc. International Conference on Robotics and Automation (ICRA)_ , 2017.
  * [18] R. Sabzevari and D. Scaramuzza, “Multi-body motion estimation from monocular vehicle-mounted cameras,” _IEEE Transactions on Robotics_ , vol. 32, no. 3, pp. 638–651, 2016.
  * [19] L. Zappella, E. Provenzi, X. Lladó, and J. Salvi, “Adaptive motion segmentation algorithm based on the principal angles configuration,” _Computer Vision–ACCV 2010_ , pp. 15–26, 2011.
  * [20] E. Elhamifar and R. Vidal, “Sparse subspace clustering,” in _Computer Vision and Pattern Recognition, 2009. CVPR 2009. IEEE Conference on_. IEEE, 2009, pp. 2790–2797.
  * [21] K. Moo Yi, K. Yun, S. Wan Kim, H. Jin Chang, and J. Young Choi, “Detection of moving objects with non-stationary cameras in 5.8 ms: Bringing motion detection to your mobile device,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops_ , 2013, pp. 27–34.
  * [22] B. D. Lucas, T. Kanade, _et al._ , “An iterative image registration technique with an application to stereo vision,” 1981.
  * [23] D. W. Eggert, A. Lorusso, and R. B. Fisher, “Estimating 3-d rigid body transformations: a comparison of four major algorithms,” _Machine vision and applications_ , vol. 9, no. 5-6, pp. 272–290, 1997.
  * [24] K. Yamaguchi, “Mexopencv,” _Collection and a development kit of matlab mex functions for OpenCV library, available at http://www. cs. stonybrook. edu/~ kyamagu/mexopencv_ , 2013.
  * [25] M. Ester, H.-P. Kriegel, J. Sander, X. Xu, _et al._ , “A density-based algorithm for discovering clusters in large spatial databases with noise.” in _Kdd_ , vol. 96, no. 34, 1996, pp. 226–231.


Graduate Student Member, IEEE

Senior Member, IEEE

Baris Donmez and Gunes Karabulut Kurt are with the Department of Electrical
Engineering and Poly-Grames Research Center, Polytechnique Montréal, Montréal,
QC H3T 1J4, Canada, (e-mail: baris.donmez@polymtl.ca).

# Continuous Power Beaming to Lunar Far Side from EMLP-2 Halo Orbit

BARIS DONMEZ  Polytechnique Montréal, Montréal, QC, Canada  GUNES KARABULUT
KURT  Polytechnique Montréal, Montréal, QC, Canada

(2024)

###### Abstract

This paper focuses on FSO-based wireless power transmission (WPT) from Earth-
Moon Lagrangian Point-2 (EMLP-2) to a receiver optical antenna equipped with
solar cells that can be located anywhere on the lunar far side (LFS).
Different solar-powered satellite (SPS) configurations which are EMLP-2
located single stable satellite and EMLP-2 halo orbit revolving single,
double, and triple satellites are evaluated in terms of 100% LFS surface
coverage percentage (SCP) and continuous Earth visibility. It is found that an
equidistant triple satellite scheme on EMLP-2 halo orbit with a semi-major
axis length of 15,000 km provides full SCP for LFS and it is essential for the
continuous LFS wireless power transmission. In our proposed dynamic cislunar
space model, geometric and temporal parameters of the Earth-Moon systems are
used in affine transformations. Our dynamic model enables us to determine the
full coverage time rate of a specific region such as the LFS southern pole.
The outcomes show that the equidistant double satellite scheme provides
SCP=100% during 88.60% time of these satellites’ single revolution around the
EMLP-2 halo orbit. Finally, the probability density function (PDF) of the
random harvested power $P_{H}$ is determined and it validates the simulation
data extracted from the stable EMLP-2 satellite and revolving satellite around
EMLP-2 halo orbit for minimum and maximum LoS distances. Although the pointing
devices to mitigate random misalignment errors are considered for the stable
and revolving SPSs, better pointing accuracy is considered for the stable
satellite. Our simulations show that the probability of $P_{H}\leq$41.6 W is
around 0.5 for the stable satellite whereas the CDF=0.99 for the revolving
satellite case for a transmit power of 1 kW.

Acquisition, tracking, and pointing (ATP), coverage analysis, Earth-Moon
Lagrangian Point (EMLP), energy harvesting, L2 halo orbit, lunar far side
(LFS), misalignment error, satellite, wireless power transfer (WPT).

##  1 INTRODUCTION

The Moon is the closest satellite to the Earth and hence it can be utilized as
a base station for deep space exploration missions. For instance, a spacecraft
can satisfy its needs (i.e., maintenance and battery recharge) before
launching interplanetary missions. Moreover, asteroid and meteor impacts on
the moon lead to the existence of platinum-group minerals on the lunar soil
whereas solar winds enable Helium-3, which is very scarce on the Earth, to be
found in the lunar regolith. Due to these rare and valuable resources on the
Moon, lunar mining is considered vital since there is a finite available
resource on Earth [1, 2, 3]. However, the tidal locking between the Moon and
the Earth leads to a permanent unobservability of one side of the Moon (i.e.,
the lunar far side (LFS)) from the Earth. Since the Moon orbits around the
world, there are 14 Earth days of total darkness, followed by the same amount
of time for complete sunlight. This causes a mission interruption due to the
run-out of batteries. Therefore, a solution for recharging the lunar equipment
and vehicles while they are working on the LFS in complete darkness is
necessary. Wireless power transfer (WPT) can be an option to overcome this
problem. Since the considered distances are significantly large in point-to-
point space link, and the losses that occur due to atmospheric attenuation,
scintillation, and fiber-coupling are negligible [4], the free-space optics
(FSO) communication technology can be considered as a proper candidate.

In our previous work, an FSO-based inter-satellite energy harvesting system is
modeled by considering realistic laser power and solar cell conversion
efficiencies. A solar-powered satellite (SPS) beams power to a small satellite
such as 1U (0.1$\times$0.1$\times$0.1 m) or 12U (0.2$\times$0.2$\times$0.3 m)
with adaptive beam divergence to maintain spot diameters of 0.1 and 0.2 m,
respectively. Independent random misalignment errors on both the transmitter
and receiver sides are mitigated by proposing an acquisition, tracking, and
pointing (ATP) system module [5]. On the other hand, existing lunar orbits in
the literature need to be assessed to find adequate cislunar orbits. The
following characteristics of lunar orbits are critical in the decision-making
process: the vicinity to the LFS, surface coverage percentage (SCP), and
stability (i.e., station-keeping necessity). There are smaller and larger moon
orbit options and they offer advantages and disadvantages so a trade-off must
be made while making a decision.

Smaller moon orbit options are low lunar orbit (LLO), elliptical lunar orbit
(ELO), prograde circular orbit (PCO), and frozen lunar orbit (FLO). The period
of LLO is around 2 hours due to the altitude of 100 km thus LLO is favorable
for lunar surface access (i.e., Apollo lunar exploration program). On the
other hand, ELO, FLO, and PCO are elliptic orbits having substantial
differences between perilune and apolune. The amplitude range of ELO, FLO, and
PCO are 100–10,000, 880–8,800 km, and 3,000–5,000 km, respectively. Their
orbital periods are similar and change between 11 to 14 hours. Moreover, these
lunar orbits have different inclinations such that LLO can have any whereas
ELO has only an equatorial inclination. Although these lunar orbits have
different inclinations, they offer limited SCP since these are smaller lunar
orbits. In general, smaller cislunar orbits compel satellites to maneuver
frequently for station-keeping and these orbits cannot provide continuous LoS
communication with the Earth due to the regular blackouts[6].

Larger moon orbit options are near rectilinear orbit (NRO), distant retrograde
orbit (DRO), and Earth-Moon Lagrangian point (EMLP) halo orbits. Since these
orbits are significantly large, the altitudes can reach even beyond 70,000 km
and their period can be up to 14 days. NROs can be located over either the
north or south pole of the moon and act as a bridge between EMLP L1 and L2
thus, they provide polar coverage. Unlike NROs, the DRO has a zero out-of-
plane amplitude (i.e., Earth-Moon plane) thus, it provides lunar equator
coverage only. It should be noted that a couple of hours of occultation occurs
in every period of DRO. On the other hand, DRO and PCO provide foreseeable
behaviors, and hence negligible corrective maneuvers are required [6].

Five equilibrium, or libration, points are established due to the
gravitational forces of moving two celestial bodies acting on each other and
that is shown by the French mathematician Joseph Lagrange in 1772. Three of
these points which are L1, L2, and L3 are collinear whereas L4 and L5 create
equilateral triangles with these two rotating planets. In the case of Earth
and Moon bodies, EMLPs L4 and L5 provide stability whereas L1, L2, and L3
require station-keeping. Since our focus is on the LFS, the halo orbit design
around the EMLP L2 point which is around 65,000 km far from the lunar
centroid, offers the best coverage among all and has a period of 8–14 days
[6].

In the existing literature, FSO-based wireless power transfer in space is
reviewed and future challenges are discussed [7]. The recent progress in the
power conversion efficiencies (PCEs) of transmitting laser diodes and laser
power converters (LPCs) for various wavelengths is summarized. Moreover,
Farquhar focuses on finding a solution to the interrupted communication link
between the Earth and the Lunar Far Side (LFS) by proposing a halo orbit on
the EMLP-2 which is one of the libration points [8]. The ratio between the
semi-minor and semi-major axes to minimize the stationkeeping is computed as
0.343. The minimum semi-minor axis to pass over the lunar occultation and
establish a continuous link with the Earth is computed as 3671 km [9].

Furthermore, [10] focuses on the feasibility of WPT from EMLP-1 and -2 to a
manned rover which requires 30 kW of harvested power. In addition, the
required transmitter aperture diameter is computed as 50 m since the distance
between these L1 and L2 points to the moon is beyond 60,000 km, and hence, a
very narrow beam divergence angle is required. The laser wavelength of 800 nm
and EHCE of 22.5 percent are preferred. In [11], the feasibility of WPT to the
lunar habitat and a rover is analyzed. Three satellites located 2,000 km away
from the lunar surface form a triangle with equal angular separation. The
laser wavelength of 850 nm and EHCE of 32 percent are considered in [11].
Besides, the required powers for the habitat and rover are 1 MW and 75 kW,
respectively. In [12], laser power versus transmitter aperture diameter (
i.e., up to 3 m) and then versus solar cell receiver diameter (i.e., up to 5
m) are analyzed for 10 kW transferred power and 1 km lunar surface to surface
WPT case. In [13], a polar LLO, an FLO, and a DRO in which three and ten
satellite constellations are evaluated to provide 100 kW to the lunar base.
The maximum laser power of 30 kW is considered for each solar-powered
satellite.

Figure 1: System model of lunar far side optical power transmission and energy
harvesting by a receiver equipped with solar cells.

There is a research gap on continuous energy harvesting on the lunar far side
that enables a receiver (Rx) optical antenna equipped with solar cells to
convert optical power to electrical power for the utilization of the research
equipment during the 14 Earth days of eclipse. To fill this gap we select
EMLP-2 and design a halo orbit with appropriate major and minor axis lengths
to minimize the stationkeeping cost while providing continuous Earth
visibility (i.e., deep space network stations (DSNSs)). We consider a dynamic
cislunar space model with a minimum number of solar-powered satellites (SPS)
that provide full LFS surface coverage so a receiver anywhere on LFS can be
power-beamed continuously. In addition, the losses induced by the misalignment
fading which are mitigated by using the beaconless pointing module (i.e.,
celestial references, internal sensors) [14] due to the propagation time of
0.2 s are considered.

Our paper investigates different EMLP-2 halo orbits and various SPS
configurations that offer full LFS coverage and Earth visibility in our
dynamic cislunar space model. The impact of the random misalignment error is
considered in energy harvesting made at the receiver equipped with solar cells
that can be located anywhere on LFS. The key contributions of this study can
be listed as follows:

  * •

We construct an EMLP-2 halo orbit by considering the semi-minor axis
$A_{y}=0.343\times A_{z}$ [8] and $A_{y}\geq 3671$ km [9] so there will be
less stationkeeping cost, and continuous Earth visibility, respectively. We
analyze the surface coverages of the LFS for a single EMLP-2 stable satellite,
a single EMLP-2 halo orbiting satellite, an EMLP-2 halo orbiting two
satellites, and an EMLP-2 halo orbiting three satellites for various semi-
major axis lengths $A_{z}$. We conclude that equidistant three satellites
revolving in EMLP-2 halo orbit with $A_{z}=15\times 10^{3}$ km achieve
continuous power beaming to a receiver optical antenna anywhere on the lunar
far side due to SCP = 100% and visibility to the Earth as well.

  * •

We establish a time-based Earth-Moon-EMLP-2 halo orbit system presented in
Fig. 3 and it enables us to compute the full SCP duration of a specific region
such as the promising south pole [15] of the LFS (Table 4). Aside triple
satellite scheme that provides continuous LFS full coverage, even a two-
satellite configuration provides full coverage during 88.60% of a full cycle
around the EMLP-2 halo orbit. We show the contribution of each satellite to
the overall SCP in detailed 3D figures at a sampled time for each scenario
(Fig. 7).

  * •

We determine the probability density function (PDF) of the random harvested
power and then compute its cumulative distribution function (CDF). It
validates our simulation data extracted from the stable EMLP-2 satellite and a
revolving halo orbit satellite whose pointing accuracies are different due to
the different amounts of mechanical vibrations that occur on their
transmitters (Figs. 8). The probability of $P_{H}\leq$41.6 W is around 0.5 for
the stable satellite whereas the CDF=0.99 for the revolving satellite case for
1 kW transmit power.

The remainder of this paper is organized as follows. In Section 2, we explain
our dynamic cislunar space model then, the PDF of random harvested power is
computed. In Section 3, the performances of various EMLP-2 satellite schemes
are evaluated based on full SCP and continuous Earth visibility first. The 3D
plots showing the contributions of each satellite to the overall SCP of LFS at
a sample time are presented. The full SCP rates of various SPS configurations
are presented by using our dynamic system model. Then, the CDF of stochastic
harvested powers which are extracted from the data are validated by the
theoretical CDF obtained in Section 2. Finally, we conclude our paper in
Section 4.

##  2 SYSTEM MODEL

In our time-based system model, different solar-powered satellite (SPS)
configurations which are stable single-satellite on EMLP-2, and revolving
single, double, and triple-satellite around the EMLP-2 halo orbit are
evaluated based on different semi-major axis length $A_{z}$ and SCP of the
LFS. The distances between the satellites are equal in each multi-satellite
scheme. The SPS-1 works actively to transfer laser power to the receiver
equipped with solar cells which can be located anywhere on LFS to provide
uninterrupted lunar missions in proximity. However, if continuous coverage
cannot be satisfied by the first SPS then the second, and even the third
satellite can start transferring laser power since a three-satellite scheme
with a specific semi-major axis length can provide continuous SCP=100% for LFS
as proved in Section 3.

The harvested power of a free space LoS optical link is expressed by [5, 16,
4] as follows

${{P}_{H}}\\!=\\!{{P}_{T}}{{\left(\frac{\lambda}{4\pi
R}\right)}^{2}}{{\eta}_{T}}{{\eta}_{H}}{{G}_{T}}{{L}_{T}}({\gamma_{T}}){{G}_{R}}{{L}_{E}}{{L}_{S}}{{L}_{C}},$
(1)

where ${{P}_{H}}$ is the electrical harvested power, ${{P}_{T}}$ is the
electrical (input) transmit power, ${{\eta}_{T}}$ and ${{\eta}_{H}}$ are the
wavelength-dependant FSO transmitter PCE and solar cell EHCE, respectively.
Moreover, ${{L}_{T}}({\gamma_{T}})$ is the transmitter misalignment loss
factor, ${{G}_{T}}$ and ${{G}_{R}}$ are the transmitter and receiver gains,
${{L}_{A}}$ is the atmospheric attenuation loss, ${{L}_{S}}$ is the
scintillation loss, and ${{L}_{C}}$ denotes the fiber coupling loss. We
consider ${{L}_{E}}={{L}_{S}}={{L}_{C}}=1$ as in [4].

The gains at the transmitting and receiving optical antennas are approximated
by, respectively [16].

$\displaystyle{{G}_{T}}\approx{{\left(\frac{\pi{{d}_{T}}}{\lambda}\right)}^{2}},$
(2)
$\displaystyle{{G}_{R}}\approx{{\left(\frac{\pi{{d}_{R}}}{\lambda}\right)}^{2}},$
(3)

where ${d}_{T}$ and ${d}_{R}$ are the transmitter and receiver optical antenna
diameters, respectively.

The required divergence beam angle of a laser diode can be computed when the
spot diameter (i.e., $d_{R}$) and the LoS distance $d$ are known [17] as
follows

$\phi={{{d}_{R}}}/{d}\;,$ (4)

The corresponding aperture diameter of the transmitter can be found as [18]

${{d}_{T}}\cong\lambda/\phi\;,$ (5)

The transmitter misalignment loss can be calculated as [16]

${{L}_{T}}({\gamma_{T}})=\exp\left(-{{G}_{T}}\gamma_{T}^{2}\right),$ (6)

where ${\gamma_{T}}$ is the radial misalignment error at transmitter.

Figure 2: Grid point method for lunar coverage analysis.

###  2.1 Misalignment Error Angle Model

In modeling the misalignment error for the stable satellite on the EMLP-2, and
a satellite revolving around the EMLP-2 halo orbit, the random misalignment
induced by the amount of mechanical vibrations is considered in an ascending
order. In other words, the maximum possible deviation from the focus point is
considered larger for a revolving satellite than the stable one. In addition,
the misalignment loss factor ${L}({\gamma})$ is inversely proportional to the
LoS distance $R$ between a satellite and the Rx optical antenna. Thus, as the
distances are in the order of thousands of kilometers, utilization of the ATP
is inevitable to mitigate the pointing errors.

ATP mechanisms that are used for the EMLP-2 satellite and any of the EMLP-2
orbiting identical satellites mitigate the pointing errors in our proposed
models and provide high accuracy in the order of nanoscale [19]. These errors
are considered independent and have different maximum deviation values as
presented in Table 1. Monte Carlo simulations are conducted to determine many
different misalignment loss factors in (6) and hence the harvested power
values can be modeled statistically.

Table 1: Simulation Parameters Parameter | Value
---|---
Transmit Power ($P_{T}$) | 1 kW [5]
Laser wavelength ($\lambda$) | 1064 nm [5]
Laser diode PCE (${{\eta}_{T}}$) | 51% [5]
Beam divergence angle (${{\eta}_{T}}$) | Adaptive [5]
EHCE (${{\eta}_{H}}$) | 50.8% [20]
| EMLP-2 stable satellite pointing
---
accuracy with ATP (${{\sigma}_{\gamma}}$)
5 nrad [14, 19]
| EMLP-2 halo orbit revolving satellite
---
pointing accuracy with ATP (${{\sigma}_{\gamma}}$)
50 nrad [14, 19]
| Optical receiver aperture diameter ($d_{r}$)
---
1 m [12]

The azimuth pointing error angle ${\gamma}_{a}$ and elevation pointing error
angle ${\gamma}_{e}$ can be modeled with a Gaussian distribution (i.e., zero-
mean) as given [21]

$\displaystyle{{f}_{X}}(x)=\frac{1}{\sqrt{2\pi{{\sigma}^{2}}}}\exp\left(-\frac{{{\left(x-\mu\right)}^{2}}}{2{{\sigma}^{2}}}\right),$
$\displaystyle{{f}_{{{X}_{1}}}}(\left.{{x}_{1}}={{\gamma}_{a}}\right|{{\mu}_{1}}=0,{{\sigma}_{1}}={{\sigma}_{a}})=\frac{1}{\sqrt{2\pi\sigma_{a}^{2}}}\exp\left(-\frac{\gamma_{a}^{2}}{2\sigma_{a}^{2}}\right),$
(7)
$\displaystyle{{f}_{{{X}_{2}}}}(\left.{{x}_{2}}={{\gamma}_{e}}\right|{{\mu}_{2}}=0,{{\sigma}_{2}}={{\sigma}_{e}})=\frac{1}{\sqrt{2\pi\sigma_{e}^{2}}}\exp\left(-\frac{\gamma_{e}^{2}}{2\sigma_{e}^{2}}\right),$

where $\sigma_{a}^{2}$ and $\sigma_{e}^{2}$ are the variances of azimuth and
elevation misalignment angles, respectively. These random pointing error
angles are independent and identically distributed.

Therefore, the transmitter radial pointing error angle $\gamma$ can be modeled
with the Rayleigh distribution by assuming
${{\sigma}_{\gamma}}={{\sigma}_{a}}={{\sigma}_{e}}$ due to the symmetry as
expressed [21]

$\displaystyle Y=\sqrt{{{X}_{1}}+{{X}_{2}}},$
$\displaystyle{{f}_{Y}}(y)=\frac{y}{{{\sigma}^{2}}}\exp\left(-\frac{{{y}^{2}}}{2{{\sigma}^{2}}}\right),$
(8)
$\displaystyle{{f}_{Y}}(\left.y=\gamma\right|\sigma={{\sigma}_{\gamma}})=\frac{\gamma}{\sigma_{\gamma}^{2}}\exp\left(-\frac{{{\gamma}^{2}}}{2\sigma_{\gamma}^{2}}\right),$

Transmitter misalignment loss can be modeled statistically by considering the
functions of the random variables approach [22, 23]. The PDF of the harvested
power random variable is determined as follows

${{f}_{H}}(h)=\frac{{{c}_{2}}}{\left|{{c}_{1}}{{c}_{2}}\right|h}\psi\exp\left(-\psi\left(\frac{\ln\left({h}/{{{c}_{2}}}\;\right)}{{{c}_{1}}}\right)\right),$
(9)

where ${{c}_{1}}=-{{G}_{T}}$, ${{c}_{2}}={{P}_{T}}{{\left(\frac{\lambda}{4\pi
R}\right)}^{2}}{{\eta}_{T}}{{\eta}_{H}}{{G}_{T}}{{G}_{R}}{{L}_{E}}{{L}_{S}}{{L}_{C}}$,
and $\psi={1}/{\left(2{\sigma_{\gamma}^{2}}\right)}\;$.

###  2.2 Surface Coverage Analysis

The selection of a lunar orbit is vital since none of the aforementioned
orbits offers perfect characteristics. Smaller orbits are unfavorable due to
the lack of coverage area and instability in general, despite the shorter
satellite-to-LFS distances. On the other hand, larger orbits provide
sufficient LFS and Earth coverage and stability, however, the maximum LoS
distances between a satellite and a receiver optical antenna equipped with
solar cells can even exceed 70,000 km for different EMLP-2 halo orbits [24].
Thus, a higher order of laser power [25] can be transmitted by a satellite to
an Rx optical antenna on LFS.

Table 2: Geometric Earth-Moon System Model Parameters Geometric Parameters | Values
---|---
Earth radius |  $6,371$ km [26]
Moon radius |  $1,737.4$ km [26]
| Lunar orbit radius
---
$385,000$ km [27]
| Earth-Moon Lagrange Point 2
---
(EMLP-2) distance to Moon
$64,500$ km [8]
Earth obliquity |  $23.44^{\circ}$ [28]
Lunar obliquity |  $6.68^{\circ}$ [28]
Lunar obliquity to ecliptic plane |  $1.54^{\circ}$ [28]
| Lunar orbital plane inclination
---
to Ecliptic plane
$5.14^{\circ}$ [28]

Our goal is to provide continuous energy harvesting from the SPS to an easily
transportable circular receiver (i.e., limited $d_{R}$), which can even be
located on the curve where LFS and lunar near side (LNS) coincide. Therefore,
different satellite schemes and orbit sizes [8, 6] are evaluated to determine
a benchmark model that offers full coverage (i.e., SCP=100%) for the LFS. The
semi-minor axis length $A_{y}$ is chosen larger than 3671 km to provide
continuous Earth visibility and has a relation with $A_{y}=0.343\times A_{z}$
to reduce the energy requirement to remain in the halo orbit [9, 8]. This
scheme is considered as a benchmark and the LoS distances of links are used
for modeling the harvested power statistically.

The uniform sampling method is used to determine the grid points on the lunar
surface by considering the sampling interval for latitude angle, $\theta_{M}$,
and longitude angle, $\varphi_{M}$, as
$\Delta{{\theta}_{M}}=\Delta{{\varphi}_{M}}={{1}^{\circ}}$. The grid point
approach [29, 30, 31] is applied to find the surface coverage percentage (SCP)
of the lunar far side which has a latitude and longitude ranges of (-90∘,90∘).
However, the latitude range for the LFS southern pole is considered as (-90∘,
-80∘) [15].

The largest possible lunar surface coverage can be achieved when the elevation
angle $\epsilon$, between the satellite and the horizon line, becomes zero.
Since we assume that there is no constraint on the moveable laser transmitter
of the satellites, the largest possible surface coverages can be achieved in
each instant of our dynamic system model.

When the largest possible coverage is attained, the nadir angle $\phi_{S}$
takes its maximum value. The nadir angle $\phi_{S}$ can be determined by the
elevation angle $\epsilon$ as follows [29]

$\sin\phi_{S}=\frac{{{r}_{M}}}{{{r}_{M}}+h}\cos\epsilon,$ (10)

where $r_{M}$ denotes the moon radius which is 1737.4 km, $h$ is the altitude.

The elevation angle is considered as $\epsilon=0$ and hence the nadir angle
$\phi_{S}$ can be rewritten as follows

$\phi_{S}=\arcsin\left(\frac{{{r}_{M}}}{{{r}_{M}}+h}\right),$ (11)

Since the elevation angle $\epsilon$, nadir angle $\phi_{S}$, and central
angle $\beta$ equals to $90^{\circ}$ as shown in Fig. 2, and $\epsilon$ is
also considered zero therefore $\beta=90^{\circ}-\phi_{S}$.

The coverage condition of a grid point, $P$, on the lunar far side surface
from a satellite located at point $S$ in the EMLP-2 orbit is given as [30]

$\arccos\left(\frac{\overrightarrow{MP}\,\cdot\,\overrightarrow{MS}}{\left\|\overrightarrow{MP}\right\|\left\|\overrightarrow{MS}\right\|}\right)<\beta,$
(12)

where $M$ is the lunar centroid thus, $\overrightarrow{MP}$ and
$\overrightarrow{MS}$ are the vectors that represent the displacements between
$M$ and $P$, and between $M$ and $S$, respectively. As per Fig. 2, $P_{1}$ and
$P_{2}$ are illustrated as the covered and uncovered points, respectively.

The surface coverage percentage can be computed precisely by considering each
grid area as
${r_{M}^{2}}\cos{{\theta}_{M}}\Delta{{\theta}_{M}}\Delta{{\varphi}_{M}}$, and
SCP can be computed by [31]

$SCP=\frac{\sum\nolimits_{i\in
S}{{{\zeta}_{i}}\,r_{M}^{2}\cos{{\theta}_{M}}\Delta{{\theta}_{M}}\Delta{{\varphi}_{M}}}}{\sum\nolimits_{i\in
S}{r_{M}^{2}\cos{{\theta}_{M}}\Delta{{\theta}_{M}}\Delta{{\varphi}_{M}}}},$
(13)

where ${{\zeta}_{i}}$ is a function that takes the value of 1 when $i^{th}$
grid point is covered, elsewhere takes zero.

Table 3: Temporal Earth-Moon System Model Parameters Temporal Parameters | Values
---|---
Earth’s rotation | 1 day [32]
Moon’s rotation | 27 days [27]
Moon’s revolution | 27 days [27]
| Earth-Moon Lagrange Point 2
---
(EMLP-2) halo orbit period
8 days [6]
Total simulation time |  $648$ hours
Sampling time |  $1$ hour
Angular velocity of Earth |  $15^{\circ}$/hour [32]
Angular velocity of Moon |  $0.56^{\circ}$/hour [27]
Laser propagation time |  $\approx 0.2$ sec

##  3 PERFORMANCE EVALUATION

The SCP performances of different satellite configurations and orbit sizes are
evaluated to determine a benchmark that offers continuous LFS SCP=100% and
Earth visibility. Then we compare the full coverage duration per period
performances of these satellite configurations for a specific region such as
the lunar south pole where ongoing researches are conducted. Since the
misalignment errors are random, we model the harvested power by the Rx optical
antenna equipped with solar cells statistically, based on the varying
distances of the orbiting satellites in our benchmark model.

###  3.1 Simulation Settings

The simulation parameters are presented in Table 1. The laser diode with 1064
nm wavelength is selected as in [5] since it is suitable for energy harvesting
and deep-space missions. The average PCE is considered as 51% for the 1064 nm
laser source [33]. As EHCE needs to be consistent with the laser source
wavelength, InGaAs-made solar cells which offer 50.8$\%$ are selected in our
system model [20]. The circular Rx optical antenna diameter is considered as 1
m which facilitates its installation and transportation [12]. It is a
challenging task to achieve a very narrow beam divergence angle to obtain a
smaller spot diameter at far distances in the order of thousands of kilometers
since larger transmitter aperture diameters are required such that $d_{T}=$ 50
m is computed in [10].

The dynamic system model of the benchmark scheme which consists of the
equidistant three satellites is presented in Fig. 3. Both geometric and
temporal Earth-Moon system model parameters shown in Tables 2 and 3 are
considered in our simulations. In addition, Earth visibility can be attained
with the closest deep space network station (DSNS) [34] rotating with the
Earth-centered Earth-fixed (ECEF) reference frame [35] as well. The affine
transformations in which intrinsic rotation and translation transforms [36,
37, 38] are utilized to determine the time-varying locations of the system
elements at every new hour until 648 hours, or one revolution of the Moon
around the Earth. This dynamic system allows us to determine a specific
region’s full coverage (i.e., SCP=100%) duration within a period.

###  3.2 Results and Discussions

Figure 3: Time-based triple satellite Earth-Moon system model

(a)

(b)

(c)

Figure 4: EMLP-2 satellite schemes (a) Single satellite, (b) Double satellite,
(c) Triple satellite

Various satellite schemes and halo orbit sizes are evaluated first to define a
benchmark that offers continuous LFS SCP=100% and Earth visibility (i.e.,
DSNS) in our dynamic cislunar space model presented in Fig. 3. The overall SCP
values for different satellite scenarios whose initial positions are presented
in Fig. 4 are obtained at each sampling time during the orbital period of 192
hours (Table 3). Single and double satellite scenarios cannot achieve SCP=100%
for different $A_{z}$ values [6] but the three-satellite configuration does.
However, not all three satellite schemes can guarantee continuous full LFS SCP
as is presented in Fig. 5 thus, different $A_{z}$ values are evaluated as
well. For instance, halo orbit with $A_{z}=5\times 10^{3}$ is not only
unsuccessful in providing SCP=100% but also fails in enabling continuous Earth
visibility as demonstrated in Fig. 6. Therefore, a three-satellite scheme
moving around the halo orbit with $A_{z}=15\times 10^{3}$ km is accepted as
our benchmark to compare with others since it becomes successful in terms of
continuous full SCP of LFS and Earth visibility as proved in both Figs. 5 and
6.

To elaborate on the SCP impact of each satellite which shows an effort to
minimize the uncovered LFS regions that remain from the preceding satellite,
the 3D LFS figures are presented in Fig. 7. It can be inferred that the
satellite on EMLP-2 is not successful in power beaming to the contours of LFS
although it provides SCP=96.864% whereas in orbiting single satellite
scenario, each and every point of the LFS can be covered but for a while only.
In multi-satellite scenarios, the second satellite shows an effort to attain
overall SCP=100% whereas our benchmark model achieves an overall full LFS SCP
throughout a period. To compare the SCP performances of these configurations
based on the time, the ratio of the duration for full coverage of a specific
region such as the promising southern pole [15], over the satellite orbital
period of EMLP-2 is considered. It is found that an equidistant 2-satellite
scheme provides SCP=100% at 88.60% of the orbital period of 8 days.

Figure 5: LFS SCP values for various $A_{z}$ in the triple satellite scheme
Figure 6: Earth visibility from the EMLP-2 halo orbit with different semi-
major axis lengths $A_{z}$

(a)

(b)

(c)

(d)

Figure 7: Minimum surface coverage (SCP) scenarios based on lunar reference
frame for different satellite schemes (a) stable single satellite, (b)
revolving single satellite, (c) revolving two satellites, (d) revolving three
satellites

The pointing accuracy of a stable EMLP-2 satellite is better than that of a
satellite revolving around EMLP-2 as considered in Table 1. In other words,
the statistical model of the random misalignment radial angle spreads more for
a revolving satellite. The minimum and maximum covered distances are computed
along with the surface coverage analyses first for the former and the latter
schemes. The CDFs of harvested powers extracted from the simulations are
validated with the theoretical results as presented in Fig. 8. The probability
of $P_{H}\leq$1.6 W is around 0.1 for the stable satellite whereas the
CDF=0.98 for the revolving satellite case for 1 kW transmit power. It can be
inferred that the harvested power is expected to be around $P_{H}\leq$1.6W in
revolving satellite scenarios whereas much higher power can be generated from
a stable satellite but not on the contour of LFS.

Table 4: Time rate comparison for LFS south pole SCP=$100\%$ LFS SP Full SCP
Duration
---
1-Satellite | 2-Satellite | 3-Satellite
44.56% | 88.60% | 100%

(a)

(b)

Figure 8: Statistical models of harvested power $P_{H}$ for (a) stable EMLP-2
satellite, and (b) revolving EMLP-2 halo satellite

##  4 CONCLUSIONS

In this paper, we evaluated different satellite configurations and EMLP-2
orbit sizes that provide continuous full lunar far-side surface coverage and
Earth visibility while adhering to the semi-major and -minor axis length ratio
that minimizes the stationkeeping cost. Single and equidistant multiple
satellite schemes in Fig. 4 were investigated with different semi-different
major axis lengths to define a benchmark with a minimum number of satellites.
We found that a three-satellite scenario with a semi-major axis length of
15,000 km attains continuous LFS full coverage and Earth visibility as proved
in Fig. 5 and 6, respectively.

Our dynamic cislunar space model presented in Fig. 3 was generated by
considering the geometric and temporal features of the Earth-Moon system. This
model enables us to investigate the time-based performances of different
scenarios having the same EMLP-2 orbit as the benchmark scheme. It was found
that single and double satellite schemes provide full coverage of the LFS
southern pole at 44.56% and 88.60% of the orbital period of 192 hours (Table
4).

The impact of each satellite on overall SCP is elaborated in Fig. 7 and it is
found that a stable single satellite cannot cover the edges where LFS and LNS
coincide however a revolving satellite on the halo orbit can power beam a
receiver optical antenna equipped with solar cells anywhere on LFS if not all
the time. The two-satellite scenario shows that equidistant satellites almost
achieve a full coverage whereas our benchmark scenario with three satellites
becomes successful with a little contribution of the third one.

Finally, the stochastic harvested powers extracted from the stable EMLP-2
satellite and a revolving halo orbit satellite scenarios in our simulations
are validated by the theoretical cumulative distribution function computed
from the previously determined probability density function (Figs. 8). The
probability of $P_{H}\leq$1.6 W is around 0.1 for the stable satellite whereas
the CDF=0.98 for the revolving satellite case for 1 kW transmit power.

## References

  * [1] M. Y. Abdelsadek _et al._ Future space networks: Toward the next giant leap for humankind _IEEE Transactions on Communications_ , vol. 71, no. 2, pp. 949–1007, 2023.
  * [2] M. Duggan, X. Simon, and T. Moseman Lander and cislunar gateway architecture concepts for lunar exploration In _IEEE Aerospace Conference_ , 2019, pp. 1–9.
  * [3] D. A. Peacock Mining on the moon _Min Eng_ , vol. 2017, pp. 23–31, 2017.
  * [4] H. Kotake, Y. Abe, M. Sekiguchi, T. Fuse, H. Tsuji, and M. Toyoshima Link budget design of adaptive optical satellite network for integrated non-terrestrial network In _IEEE International Conference on Space Optical Systems and Applications (ICSOS)_ , 2022, pp. 240–247.
  * [5] B. Donmez, I. Azam, and G. Karabulut Kurt Mitigation of misalignment error over Inter-Satellite FSO energy harvesting In _Workshop: Mega Constellation Satellite Network for 6G (PIMRC’23 WS MCSN-6G)_ , Toronto, Canada, Sep. 2023, p. 5.
  * [6] R. Whitley and R. Martinez Options for staging orbits in cislunar space In _IEEE Aerospace Conference_ , 2016, pp. 1–9.
  * [7] Y. Zheng _et al._ Wireless laser power transmission: Recent progress and future challenges _Space Solar Power and Wireless Transmission_ , 2024.
  * [8] R. W. Farquhar _The utilization of halo orbits in advanced lunar operations_. National Aeronautics and Space Administration, 1971, vol. 6365.
  * [9] R. W. Farquhar Lunar communications with libration-point satellites. _Journal of Spacecraft and Rockets_ , vol. 4, no. 10, pp. 1383–1384, Oct. 1967.
  * [10] M. D. Williams _et al._ Power transmission by laser beam from lunar-synchronous satellites to a lunar rover In _Intersociety Energy Conversion Engineering Conference_ , Aug. 1992.
  * [11] M. D. Williams, J. H. Kwon, G. H. Walker, and D. H. Humes Diode laser satellite systems for beamed power transmission 1990.
  * [12] T. W. Kerslake Lunar surface-to-surface power transfer In _AIP Conference Proceedings_ , vol. 969, no. 1. American Institute of Physics, 2008, pp. 466–473.
  * [13] F. Lopez, A. Mauro, S. Mauro, G. Monteleone, D. E. Sfasciamuro, and A. Villa A lunar-orbiting satellite constellation for wireless energy supply _Aerospace_ , vol. 10, no. 11, p. 919, 2023.
  * [14] H. Kaushal and G. Kaddoum Optical communication in space: Challenges and mitigation techniques _IEEE Communications Surveys & Tutorials_, vol. 19, no. 1, pp. 57–96, 2017.
  * [15] J. Stopar and H. M. Meyer Topography and Permanently Shaded Regions (PSRs) of the Moon’s South Pole (80°S to Pole) May 2019. [Online]. Available: https://repository.hou.usra.edu/handle/20.500.11753/1255
  * [16] A. Polishuk and S. Arnon Optimization of a laser satellite communication system with an optical preamplifier _J. Opt. Soc. Am. A_ , vol. 21, no. 7, pp. 1307–1315, Jul. 2004.
  * [17] A. G. Alkholidi and K. S. Altowij Free Space Optical Communications — Theory and Practices In _Contemporary Issues in Wireless Communications_ , M. Khatib, Ed., Rijeka, 2014.
  * [18] Z. Ghassemlooy, W. Popoola, and S. Rajbhandari _Optical Wireless Communications: System and Channel Modelling with MATLAB®, Second Edition_. CRC press, 2019.
  * [19] M. M. Freitas, M. A. Fernandes, P. P. Monteiro, F. P. Guiomar, and G. M. Fernandes Requirements and Solutions for Robust Beam Alignment in Fiber-Coupled Free-Space Optical Systems _Photonics_ , vol. 10, no. 4, 2023.
  * [20] N. A. Kalyuzhnyy _et al._ Optimization of photoelectric parameters of ingaas metamorphic laser ($\lambda$= 1064 nm) power converters with over 50% efficiency _Solar Energy Materials and Solar Cells_ , vol. 217, p. 110710, 2020.
  * [21] S. Arnon, S. Rotman, and N. S. Kopeika Beam width and transmitter power adaptive to tracking system performance for free-space optical communication _Appl. Opt._ , vol. 36, no. 24, pp. 6095–6101, Aug. 1997.
  * [22] A. Leon-Garcia _Probability and Random Processes for Electrical Engineering_. Pearson Education India, 1994.
  * [23] S. Chan _Introduction to Probability for Data Science_. Michigan Publishing, 2021.
  * [24] A. Yuan, Y. Xu, Q. Zhang, R. Wang, and Z. Yang Optimal energy efficiency for relay selection via power allocation in l2 halo-orbital cislunar communication network _IEEE Transactions on Vehicular Technology_ , vol. 71, no. 3, pp. 3155–3170, 2022.
  * [25] Ursi white paper on solar power satellite (sps) systems and report of the ursi inter-commission working group on sps 2007.
  * [26] D. R. Williams Moon Fact Sheet . [Online]. Available: https://nssdc.gsfc.nasa.gov/planetary/factsheet/moonfact.html
  * [27] S. C. Editors Moon Facts - NASA Science . [Online]. Available: https://science.nasa.gov/moon/facts/
  * [28] C. Olthoff, D. Kaschubek, and M. Killian Dynamic thermal interactions between spacesuits and lunar regolith in permanently shaded regions on the moon _Acta Astronautica_ , vol. 203, pp. 351–369, 2023.
  * [29] S. Cakaj _Ground Station Design and Analysis for LEO Satellites: Analytical, Experimental and Simulation Approach_. John Wiley & Sons, 2022.
  * [30] Z. Song, G. Dai, M. Wang, and X. Chen A novel grid point approach for efficiently solving the constellation-to-ground regional coverage problem _IEEE Access_ , vol. 6, pp. 44 445–44 458, 2018.
  * [31] Z.-Y. Gao and X.-Y. Hou Coverage analysis of lunar communication/navigation constellations based on halo orbits and distant retrograde orbits _The Journal of Navigation_ , vol. 73, no. 4, pp. 932–952, 2020.
  * [32] S. C. Editors Chapter 2: Reference Systems - NASA Science . [Online]. Available: https://science.nasa.gov/learn/basics-of-space-flight/chapter2-1/
  * [33] H.-Y. Lin, M.-Y. Wu, J. Tang, M.-J. Zhao, and D. Sun Compact Nd:YVO4 laser 1087.5 nm pumped by an edge-emitting LD _Optik_ , vol. 185, pp. 414–417, 2019.
  * [34] S. D. Slobin Coverage and geometry _DSN No. 810_ , vol. 5, pp. 810–005, 2014.
  * [35] J. Farrell _Aided Navigation: GPS with High Rate Sensors_. McGraw-Hill, Inc., 2008.
  * [36] O. Rodrigues Des lois géométriques qui régissent les déplacements d’un système solide dans l’espace, et de la variation des coordonnées provenant de ces déplacements considérés indépendamment des causes qui peuvent les produire _Journal de mathématiques pures et appliquées_ , vol. 5, pp. 380–440, 1840.
  * [37] R. P. Paul _Robot Manipulators: Mathematics, Programming, and Control: The Computer Control of Robot Manipulators_. Richard Paul, 1981.
  * [38] D. House and J. C. Keyser _Foundations of Physically Based Modeling and Animation_. CRC Press, 2016.

[] Baris Donmez (Graduate Student Member, IEEE) received the B.Sc. and M.Sc.
degrees (with high honors) in electrical and electronics engineering, in 2009
and 2022, respectively, from FMV Işık University, Istanbul, Turkiye. Now, he
is currently working toward a Ph.D. degree in electrical engineering at
Polytechnique Montréal, Montréal, QC, Canada

His research interests include communication and energy harvesting systems in
deep-space networks.

[] Gunes Karabulut Kurt (Senior Member, IEEE) is a Canada Research Chair (Tier
1) in New Frontiers in Space Communications and Associate Professor at
Polytechnique Montréal, Montréal, QC, Canada. She is also an adjunct research
professor at Carleton University, ON, Canada. Gunes is a Marie Curie Fellow
and has received the Turkish Academy of Sciences Outstanding Young Scientist
(TÜBA-GEBIP) Award in 2019. She received her Ph.D. degree in electrical
engineering from the University of Ottawa, ON, Canada.