# Hyperpixels: Pixel Filter Arrays of Multivariate Optical Elements for
Optimized Spectral Imaging

Calum Williams 1  Richard Cousins 2  Christopher J. Mellor 3  Sarah E.
Bohndiek  4,5 and George S.D. Gordon  6  1Department of Physics, University of
Exeter, Exeter, EX4 4QL, United Kingdom  2Nano- and Micro-scale Research
Centre, University of Nottingham, UK, University Park, Nottingham NG7 2RD, UK
3School of Physics and Astronomy, University of Nottingham, University Park,
Nottingham NG7 2RD, UK  4Department of Physics, Cavendish Laboratory,
University of Cambridge, JJ Thomson Avenue,Cambridge, CB3 0HE, UK  5, Cancer
Research UK Cambridge Institute, University of Cambridge, Robinson Way,
Cambridge, CB2 0RE, UK  6 Optics & Photonics Research Group, Department of
Electrical and Electronic Engineering, University of Nottingham, Nottingham,
NG7 2RD, UK<EMAIL_ADDRESS>
###### Abstract

We introduce the concept of ‘hyperpixels’ in which each element of a pixel
filter array (suitable for CMOS image sensor integration) has a spectral
transmission tailored to a target spectral component expected in application-
specific scenes. These are analogous to arrays of multivariate optical
elements that could be used for sensing specific analytes. Spectral tailoring
is achieved by engineering the heights of multiple sub-pixel Fabry-Perot
resonators that cover each pixel area. We first present a design approach for
hyperpixels, based on a matched filter concept and, as an exemplar, design a
set of 4 hyperpixels tailored to optimally discriminate between 4 spectral
reflectance targets. Next, we fabricate repeating 2$\times$2 pixel filter
arrays of these designs, alongside repeating 2$\times$2 arrays of an optimal
bandpass filters, perform both spectral and imaging characterization.
Experimentally measured hyperpixel transmission spectra show a 2.4x reduction
in unmixing matrix condition number ($p$=0.031) compared to the optimal band-
pass set. Imaging experiments using the filter arrays with a monochrome sensor
achieve a 3.47x reduction in unmixing matrix condition number ($p$=0.020)
compared to the optimal band-pass set. This demonstrates the utility of the
hyperpixel approach and shows its superiority even over the optimal bandpass
case. We expect that with further improvements in design and fabrication
processes increased performance may be obtained. Because the hyperpixels are
straightforward to customize, fabricate and can be placed atop monochrome
sensors, this approach is highly versatile and could be adapted to a wide
range of real-time imaging applications which are limited by low SNR including
micro-endoscopy, capsule endoscopy, industrial inspection and machine vision.

††journal: oe††articletype: Research Article

##  1 Introduction

Spectral imaging is a powerful technique that captures wavelength-resolved
information (spectra) at every pixel in an image to non-invasively estimate
relative abundances of spectral signatures, representing for example the
chemical composition of spatially heterogeneous samples [1, 2, 3, 4]. The
interaction of light with a given sample, for example through absorption,
reflectance or fluorescence, provides a spectral signature that is determined
by the relative abundance of key chromophores or fluorophores. Applying
‘spectral unmixing’ algorithms to extract chemical abundances can then provide
valuable diagnostic information across a wide range of disciplines, from
biomedical imaging to agricultural inspection [5, 6, 7, 8, 9, 10].

Real-time spectral imaging is commonly realized using multispectral filter
arrays (MSFAs). MSFAs are mosaics of pixelated optical filters integrated atop
an image sensor; images at each spectral band are formed by interpolation
between spatially distributed filters using a demosaicking algorithm [1, 12,
13, 11]. Conventionally, each element in the MSFA is an optical bandpass
filter, usually with a fixed, narrow bandwidth. The center wavelengths of each
spectral band can either be spaced evenly across the spectral range of the
sensor (e.g. 400 – 900 nm for silicon image sensors) for generic applications,
or can be spectrally and spatially optimized to a particular application via
computational processes [4, 13, 12]. One key shortcoming of bandpass filters
is that much of the incoming spectral power—outside of the bandpass region—is
discarded (reflected or absorbed). This can reduce performance, in particular
the signal-to-noise ratio (SNR), which may be problematic in environments with
high noise or interference from other spectral signatures.

If _a priori_ information of the sample is known then an optimal spectral
imaging approach—with regards to SNR—would be to match the spectral
transmission of the filters to the spectrum of each contributing chromophore;
the concept of ‘matched filters’ is widely used for SNR maximization in
telecommunications [14] and other fields. For optical filters, matching has
been achieved using optimized thin-film multilayered elements called
multivariate optical elements [15] in challenging spectroscopy applications
where highly sensitive detectors cannot be easily deployed e.g. in undersea
oil and gas extraction [16, 17]. However, multivariate optical elements have
only been used to make spatially uniform filters and have not been applied as
pixelated filter arrays.

Demonstrating a matched filter concept in an MSFA has been elusive to-date due
to the complexity of creating arbitrary (i.e. not single narrow-bandpass)
spectral characteristics within pixel filter arrays. Previous studies have
achieved some compact, broadband solutions using: nanohole metasurfaces [18,
19]; multilayer structures with spatially varying holes etched into them [20];
metasurfaces designed to perform ‘pixel routing’ where sub-wavelength 3D
structures route light to different parts of the image sensor based on their
wavelength [21]; random or semi-randomized diffusers [22]; and coded apertures
[23]. Nonetheless, it is challenging to tailor spectra with these approaches
due to computationally intensive inverse design processes and designs
typically require high resolution (nanoscale) lithographic techniques with
tight tolerances. Most recently MSFAs with each pixel made up of a Fabry-Perot
etalon with random cavity height (and hence resonant wavelength) have been
demonstrated [24], which also show potential for multiband transmission but
with potential for an easier design and fabrication process.

Building on these previous works, we demonstrate here a new matched filter
concept, ‘hyperpixel’ filter arrays (see Figure 1), that can be deployed as an
MSFA for snapshot multispectral imaging. Hyperpixels comprise multiple
spatially heterogeneous Fabry-Perot resonators, termed _subpixel filters_ ,
located within the bounds of a single image sensor pixel. Each of the subpixel
filters can have a different cavity height, so that when integrated over the
pixel area, the overall transmission is a composite spectral filter based on
the sum of these responses. The hyperpixel approach overcomes the design and
fabrication demands of prior approaches, being easily designed based on cavity
heights and fabricated using grayscale lithography [25]. As a result of the
hyperpixels implementation, spectra can be easily unmixed using only a single
standard matrix multiplication that can be implemented in real-time using
existing cost-effective hardware. We demonstrate here the hyperpixels concept
using 100 subpixel filters to tailor spectra in a 10$\times$10 $\mu$m
hyperpixel targeted to four distinct reflectance spectra and arranged in a
hyperpixel filter array. We demonstrate the potential performance improvement
afforded using a MacBeth ColorChecker chart compared to a standard band-pass
filter array, representing the best possible conventional MSFA.

Figure 1: Hyperpixel concept: a) Conventional color cameras have Bayer filters
applied per-pixel to discern red, green and blue wavelengths. b) Conventional
snapshot multispectral cameras based on multispectral filter arrays (MSFAs)
extend this using bandpass filters centered at many different wavelengths – 16
in this example. c) hyperpixels replace the conventional MSFA filters and
enable different pixels spectral responses to be ‘matched’ to the reflectance
spectrum of a target substance. d) Hyperpixels are implemented by fabricating
several Fabry-Perot resonators within the area of a single pixel on a CMOS
sensor. The area of each resonator dictates the power transmitted in that
particular band, but the total power integrated across the spectrum cannot
exceed that of an equivalent single band-pass filter.

##  2 Methods

###  2.1 Hyperpixel design

Hyperpixels were implemented as an array of thin-film metal-insulator-metal
Fabry-Perot filters bounded within the pixel active area of a CMOS image
sensor (Figure 2a). For alignment simplicity, we assume a 10 $\mu$m pixel size
and divide this up into 100 ‘subpixel filter’ units. The number of subpixel
filter units is adjustable but should be a square number to match the geometry
of the CMOS image sensor pixels. Subpixel filters should also be $>$ 500nm to
avoid performance drop-off experience as the dimension approaches the
wavelength of light [26] and to minimize fabrication artefacts as structures
become smaller. Working within these constraints we choose 100 subpixel
filters as as a compromise between control over spectral shape and quality of
subpixel filter.

Each of these subpixel filters comprises a single cavity filter with Ag
mirrors (thickness 22 nm) and resist cavity (variable thickness $\sim$75-250
nm). A MgF2 encapsulation (capping) layer was used to mitigate oxidization and
surface tarnishing. Silica glass (hundreds of micron thickness) was used as
the substrate, which permits straightforward integration atop commercial
monochrome image sensors or optical alignment (using relay optics) to the
image sensor. To determine the optimal multilayer thicknesses for high
transmission and sufficiently narrow passbands, numerical simulations were
performed using Lumerical FDTD Solutions (Ansys) (Figure 2b).

Following simulations, an initial array with evenly spaced resist thicknesses
from $\sim$50 to 250nm was fabricated as a height calibration target. Actual
heights were measured using a combination of atomic force microscopy (AFM) and
optical micro-spectrophotometery, as described in our previous works [27] and
used to optimize exposure dose by calibrating the relationship between actual
(measured) height and designed height. Height calibration was iterated several
times until heights are matched, because the optimal exposure dose profile is
dependent on numerous factors including the substrate material/s, processing
parameters (i.e. development times, temperatures, proximity effect etc).
Optical micro-spectrophotometery was performed on the sub-pixels to
characterize their optical responses (transmission spectra) and locate
relative positions of center wavelengths, second order resonances, bandwidths
(full width at half maximum, FWHM) and transmission.

To create a hyperpixel from the subpixel filters, we must next determine the
required heights of each of the 100 subpixel filters so that they match the
desired spectral transmission response. To do this, we begin by identifying
the target transmission spectrum required (e.g. Figure 2d). We must then
approximate this as a sum of bandpass filters, which represent the subpixel
filters. The number of independent bandpass filters used need not be the same
as the number of subpixel filters, and can be thought of as a method of
quantizing the heights of the subpixel filters. This quantization process may
be important for future manufacturing. For example if one exposure is required
per grayscale height level (e.g. with a stepper) fewer heights would allow
faster manufacturing. Similarly, there may be tolerances on the grayscale
heights meaning that only finite number of distinct heights are realistically
manufacturable. We therefore choose a quantization level of 32, which enables
up to 32 distinct heights. We find this number enables fast and reliable
convergence of our design optimization algorithm with high fidelity curve
approximation, while still producing good quality fabrication results. The
optimization process may assign certain certain wavelengths weights of zero,
or may assign a certain centre wavelength multiple times so the final
quantization level may even be less than 32.

The optimization algorithm uses experimentally fitted spectral data to map a
resist thickness to a spectral transmission profile and create 32 bandpass
spectra that are then summed together to get the total transmission spectral
profile. The center wavelengths and weights in the final summation are set as
the adjustable parameters that are then optimized using a basin-hopping
algorithm. The error function is a mean-square error with the target spectrum,
which is the expected reflectance spectrum of the target, including the effect
of non-flat spectral profile of the input illumination. The output of this
process is a composite spectrum that closely matches the target (Figure 2d).
The optimization algorithm is implemented using custom Python code (available
at [34]).

Finally, these 32 component spectra must be mapped onto the 100 available
subpixel filters within each hyperpixel. To do this, centre wavelengths are
first rounded to the nearest 3nm reflecting again the achievable height
precision measured using AFM previously [27]. The algorithm produces weights
for each of these components that vary from 0 to 1. In the final design these
weights will be implemented by using the same height value for a greater
number of the 100 available subpixels: greater weight for a particular
component spectrum means more subpixels will have the height that creates that
component. To achieve this, the weights within each hyperpixel are normalized
and rounded to the nearest integer such that they sum to 100. These new
normalized weights are then used to determine how many of the 100 subpixel
filters (i.e. 10$\times$10 array) should have that particular height.

Once the subpixel filters heights are determined, they are placed in
increasing height order in a snaking pattern to minimize height steps between
adjacent subpixel filters, which should enable more reliable fabrication
(Figure 2c).

Figure 2: Design process for hyperpixels. a) Structure of constituent Fabry-
Perot cavities making up each hyperpixel. $t_{m1}$ = 22nm, $t_{m2}$ = 22nm and
$t_{ins}$=75–200nm. b) Simulation data showing the theoretical relationship
between cavity thickness and centre transmission wavelength. c) Example
spatial design of hyperpixel showing ordering by increasing height of sub-
pixels for more reliable fabrication. d) Example transmission spectrum
produced by hyperpixel design showing a close match to target spectrum.

###  2.2 Filter array design

To design a test set of filter arrays to evaluate performance of hyperpixels
against conventional bandpass filters, we identified a subset of target
spectra from a MacBeth ColorChecker chart (Figure 3a) and designed hyperpixels
optimally matched to them. We combined manufacturer provided spectral
reflectance data with the spectrum of the illumination source to be used (OSL2
Halogen source, ThorLabs) to determine the total expected reflected spectral
profile. These spectra were used as the target design spectra for the
hyperpixels.

Figure 3: Designing hyperpixels for optimal spectral unmixing. a) Selected
spectral reflectance targets on a Macbeth color checker chart, showing
spectral reflectance profiles. b) Unit cell of designed hyperpixel filter
array, including 8 10$\times$10 $\mu$m pixel designs: 4 hyperpixel designs
(comprising 10$\times$10 1$\mu$ subpixel filters) and 4 bandpass filters with
optimized centre wavelengths. c) Layout of the full filter array as
fabricated, including height calibration targets around the edge. f) Simulated
spectra of bandpass filter pixels from filter array. e) Simulated spectra of
hyperpixel filters from filter array.

Four hyperpixel designs were produced as shown in Table 1.

Table 1: Specifications of designed set of four hyperpixel filters. | HP filter 1 | HP filter 2 | HP filter 3 | HP filter 4
---|---|---|---|---
Height (nm) | No. subpixel filters
87 | 2 | 3 | 1 | 2
93 | 3 | 10 | 1 | 8
96 | 4 | 8 | 0 | 26
99 | 0 | 0 | 2 | 0
102 | 11 | 6 | 12 | 16
108 | 8 | 13 | 23 | 8
111 | 24 | 3 | 0 | 9
114 | 0 | 0 | 27 | 0
117 | 24 | 27 | 25 | 10
123 | 22 | 29 | 9 | 14
126 | 1 | 1 | 0 | 5
132 | 1 | 0 | 0 | 2

As a benchmark we also designed 4 bandpass filters optimized for spectral
unmixing of the same 4 selected targets. Specifically, for the set of 4
bandpass filters, parameterized by 4 centre wavelengths, we multiply the
filter spectral transmission profile with the expected spectral reflectance of
each of the target areas. The output is the $4\times 4$ matrix which, when
inverted, produces an ‘unmixing’ matrix that recovers the correct abundances
of the target spectra. The condition number of this matrix represents how
‘well’ the unmixing can be performed in the presence of noise and so we
minimize this quantity as our optimization objective to produce optimal
performance of the bandpass filters [28]. Similarly, the condition number of
the hyperpixel filter set can also be computed. The ratio of these condition
numbers is then computed to determine whether the hyperpixels should perform
better or worse than the optimal bandpass filters for a particular set of
spectral targets:

$r=\mathrm{cond}(\mathbf{H}_{bandpass})/\mathrm{cond}(\mathbf{H}_{hyperpixel})$
(1)

where $\mathbf{H}_{bandpass}$ is the $4\times 4$ unmixing matrix of the
optimized bandpass filters and $\mathbf{H}_{hyperpixel}$ is the unmixing
matrix for the hyperpixel filters. For the set of spectral reflectance targets
in the MacBeth color chart, we find that $r$ can vary from 0.1 (bandpass
filters perform better) up to to 1.6 (hyperpixels perform better).

Finally, we create a layout for fabrication of the filter arrays with the 4
hyperpixels adjacent to the bandpass filters (Figure 3b). By placing these
regions together in a repeating array, we ensure that any localized
fabrication variations (e.g. proximity effect from e-beam lithography, resist
thickness etc) apply equally to all the different designs. The unit cells have
dimension 20$\mu$m$\times$20$\mu$m and are repeated into an array of 1mm
$\times$ 2mm (Figure 3c). Around the edge of the array are height calibration
targets with a linear ramp of 0 to 250nm so that fabricated heights be
validated with AFM or stylus profiler.

Using our previously described calibration method, which combines AFM and with
optical micro-spectrophotometery and simulations (Figure 2b), we map the
target heights to expected spectral profiles for the bandpass filters, as
shown in Figure 3d. Similarly, we can also produce expected spectra for the
hyperpixels by including second-order resonances using our previously
described characterization approach [27].

###  2.3 Fabrication

22 mm high precision glass coverslips (ThorLabs) were used as the substrates.
Before processing, these were successively cleaned in ethyl lactate, acetone,
methanol, and isopropanol (IPA), then dehydrated on a contact hotplate at
150oC for 3 minutes. For the first mirror layer, a thermal evaporator was used
to deposit: a 1 nm Ti adhesion layer followed by 22 nm Ag and MgF2 capping
layer. The Ag was evaporated at a high deposition rate to improve the film
quality.

As soon as possible after evaporation, a film of AR-7720 (Allresist GmBh) was
spin-coated at 4000 rpm for 60 s. The resist was diluted using PGMEA
(1-methoxy-2-propanol acetate), so an initial film thickness of 165 nm was
produced. The sample was then placed in an electron beam lithography tool,
NanoBeam nB5 (NanoBeam Ltd.). To vary the height of each pixel, the dose was
varied using an earlier dose test which gave the contrast as 0.55. The exposed
sample then had a post-exposure bake on a hotplate at 1100C for 60 s followed
by a bake in an oven at 90oC for 2 hours; this latter step reduces the surface
roughness of the film. The resist was developed in a TMAH (Tetramethylammonium
hydroxide) based developer for 90 s and then rinsed with DI (deionized) water.
This resulted in a resist pattern with varying heights. The second mirror
layer was thermally evaporated, consisting of 22 nm of Ag and MgF2
encapsulation layer to prevent silver tarnishing, and is kept to $<$ 30 nm.

###  2.4 Spectroscopic characterization

For spectral characterization, a modified Olympus BX-51 optical microscope was
used to analyse the hyperpixel filter arrays. A UV-visible spectrometer (Ocean
Optics HR2000+) was used to collect the spectral transmission data, while a
digital camera (Lumenera Infinity-2 2MP CCD) was used to image the surface.
The spectra collected were normalized using OceanView (Ocean Optics) software
by collecting the bright and dark reference spectra. The bright spectra is
obtained by using a non-metal coated area of the filter array i.e. just the
light source transmitting through the glass.

The spectra of each sub-pixel was collected using a series 100x objective and
centering the collection area over the center of the pixel. 4 spectra (at
different spatial positions in the adjacent array areas) were collected and
averaged to reduce noise.

In order to examine the variations in fabrication quality across the sample,
14 different regions across the fabricated filter array were sampled and
spectra recorded from each bandpass and hyperpixel in the unit cell. The
spectra from the bandpass filters were further analysed by fitting a sum of
two Gaussian curves to identify the first- and second-order Fabry-Perot
resonance wavelengths and the FWHM. Direct measurement of the hyperpixel
spectra is challenging because the circular aperture of the spectrometer
collection fibre is mismatched with the square aperture of the pixel.
Therefore, the bandpass filter spectra are used as a measure of deviation from
designed height in each of the 14 sampled areas.

First, the centre wavelength of each of the four bandpass areas is measured
and plotted against the target height from the original design. A line is
fitted to this. Then the target heights of the each of the subpixel filters in
the hyperpixel filters are input into this fitted line and an adjusted actual
centre wavelength for each of the subpixel filters is found. This adjusted
centre wavelength, along with other parameters from the BP areas such as the
FWHM and position of second-order resonance, are used to construct a
transmission for each subpixel area. These are then summed to infer the
overall hyperpixel filter response in this area of the filter array.

Based on these measured bandpass and inferred hyperpixel spectra, an unmixing
matrix is formed by multiplying these transmission spectra by the ideal
reflected spectra expected from the target areas on the color chart. The
condition numbers of these matrices for the hyperpixel and bandpass areas are
then computed and compared using a paired $t$-test.

###  2.5 Spectral imaging characterization

Our validation experiment was to image a real Macbeth color chart target using
the fabricated filter arrays placed on top of a CMOS sensor image sensor. An
optical imaging relay system was constructed for spectral imaging. The system
contains a first imaging lens: fixed-focal length camera lens (f=35mm,
Navitar), and a set of relay lenses (Thorlabs, AC254-075-A-ML). The target
scene, illuminated by a halogen light source coupled to a ring light
illuminator (Thorlabs OSL2IR), is imaged onto a monochrome image sensor
(Basler, a2A3840-45umPRO) via a conjugate plane located between two 4x planar
achromat objectives (Thorlabs, RMS4X). The conjugate plane contains the
hyperpixels mounted to a XYZ-translation stage (Thorlabs, MBT616D). Imaging
the scene through the filter pixels in this way mimics a conventional snapshot
multispectral camera whereby MSFAs are integrated directly atop the pixels.

###  2.6 Image processing

Raw monochrome images from the camera were first straightened and cropped such
that only the color chart area was included. Images were then demosaicked to
produce multispectral image data cubes by manually identifying the same pixel
type in 5-10 separate locations across the image then fitting a grid to these
locations to determine the periodicity, location of each pixel type and offset
(custom interactive MATLAB code, available from [34]). Once this grid is
determined, the raw pixels within each filter array pixel are summed to form a
single point in the hypercube. The demosaicking and interpolation process
represents a reduction in resolution from the raw images: the raw images are
$3840\times 2160$ but each pixel in the pixel filter array covers around 11.05
sensor pixels as the pixel sizes were not matched to the sensors. The raw
demosaicked images were therefore $174\times 49$ pixels. Because the image
height is scaled down twice as much due to the rectangular shape of the unit
cell, the images in the hypercube are resampled to $1920\times 1080$ to
preserve aspect ratio.

A second processing step was then performed on the demosaicked images to
identify each of the different spectral reflectance areas within the color
chart. First, the flat spectral reflectance areas (white and black) were
identified and used in a similar fashion to white and dark references in
spectroscopy; pixels within these areas are averaged to create white and dark
reference, $r_{white}$ and $r_{dark}$ respectively, for each of the channels
in the hypercube. The dark reference was subtracted from each of the channels
in the hypercube and then the white reference was used as a normalization
value to divide values in each channel. Specifically:

$\hat{p}(x,y)=\frac{p(x,y)-r_{dark}}{r_{white}-r_{dark}}$ (2)

where $\hat{p}(x,y)$ is the corrected pixel value at position $(x,y)$,
$p(x,y)$ is the raw pixel value and $r_{white}$ and $r_{dark}$ are the white
and dark reference values respectively.

Because of the fixed power constraint of the hyperpixel design, we expect that
all pixels in the white reference areas should transmit approximately the same
total power and so this normalization is applied independently to each channel
in the BP and hyperpixel filter arrays.

Finally, the four target spectral areas are identified. The pixel values
within each of these 4 areas are averaged and this is repeated for the 4
channels representing the hyperpixels and bandpass filters. The result is 16
values for each of the bandpass and hyperpixel cases, which are arranged into
the $4\times 4$ unmixing matrices, similarly to the spectral analysis case.
The condition numbers of these matrices represent quality of unmixing. For
visual comparison, the 4 target regions in the hypercube images are unmixed,
using the computed spectral unmixing matrix, and are plotted in a $2\times 2$
array for each of the hyperpixel and bandpass cases. To quantify the effect of
fabrication variation across the sample on performance, the unmixing matrices
are computed for 10 different areas across the pixel filter array and a
$t$-test is used to establish whether a statistically significant reduction in
condition number is achieved.

##  3 Results

###  3.1 Pixel-wise spectroscopic characterization

Experimentally fabricated filter arrays (Figure 4a) were used to obtain
spectral response (transmission) results for 14 different spatial regions
(Figure 4) from the centre of the array to the outer edges. The measured
transmission spectra from the bandpass filters (Figure 4b) indicate a
relatively large variation in central peak position across the sample, though
some areas do closely match the designed structures. The inferred spectra for
the hyperpixels (Figure 4c. were reconstructed from the corresponding
wavelength shift of adjacent bandpass filters. By aggregating the observed
shifts in centre wavelength of the bandpass filters (Figure 4d) a normal
distribution fit can be made, which gives a standard deviation of 47.7 nm in
centre wavelength across our sample.

Due to a fabrication specification error, the hyperpixel spectral are blue-
shifted by 50-100nm (dashed lines in Figure 4c) compared to the target designs
(solid lines in Figure 4c). This may reduce performance, however a close
inspection of the spectra measured from hyperpixels in different filter array
regions (dotted lines in Figure 4c) shows that the variation in spectral shift
across the filter array is comparable in magnitude to this blue shift and in
many cases actually compensates. At a minimum we therefore expect that at
least some areas of the sample will exhibit hyperpixel spectra matching the
design specification. Further, we note that since this fabrication blueshift
did not occur with the bandpass filters these ought to closely match the
optimal designs (save for variation across the sample). Any observed
performance advantage of hyperpixels is therefore still being measured against
an optimal benchmark. Such variations can be minimized with tighter control
over fabrication tolerances and processes in future.

Despite the limitations in the fabrication process, which led to inaccuracies
in the filter properties compared to design, the measured condition number for
the unmixing matrices of bandpass and hyperpixel filters shows that
hyperpixels produce statistically significantly lower condition numbers than
optimized bandpass filters (ratio = 2.41x, 95% CI = 1.10 – 5.30, $p$=0.031)
for this set of reflectance targets (Figure 4e).

Figure 4: Spectral characterization of fabricated hyperpixel filter array. a)
Microscope image showing repeating $2\times 4$ array of hyperpixels and
bandpass filters. b) Measured spectra from bandpass filters. Solid lines:
designed spectra, dotted lines: measured spectra from 14 different areas on
filter array. c) Reconstructed spectra from hyperpixel filters. Dashed lines:
designed spectra, solid lines: blue-shifted spectra following fabrication,
dotted lines: inferred spectra from 14 different areas on filter array. d)
Histogram showing error between designed centre wavelength and measured centre
wavelength for bandpass filters across 14 locations on filter array (normal
distribution fit: $\mu$=0, $\sigma$=47.4nm. e) Comparison of condition number
of unmixing matrix shows a statistically significant reduction (ratio = 2.41x,
95% CI = 1.10 – 5.30) of the hyperpixels compared to optimal bandpass filters.

###  3.2 Spectral imaging characterization

Spectral imaging data were next obtained using the MacBeth Color Checker chart
scene (Figure 5a, obtained using a smartphone camera; target reflectance
squares are circled). The raw unprocessed image through the hyperpixel filter
array (Figure 5b; mosaic visible on zoomed inset) was demosaiced to produce 2
spectral hypercubes: bandpass (Figure 5c) and hyperpixel (Figure 5d). Each
hypercube comprises 4 planes representing the different spectral channels.
From these hypercubes, relevant sections of the image were cropped to
highlight the target spectral reflectance areas (Figure 5d and e).

Figure 5: Imaging experiments with hyperpixel filter array. a) Conventional
RGB image of sample scene. b) Raw image recorded with filter array placed atop
CMOS sensor, inset shows the mosaic pattern. c) Example demosaicked hypercube
of images for bandpass filters. d) Example demosaicked hypercube of images for
hyperpixel filters. e) The appropriate spectral reflectance target areas are
cropped from the full image frames, and then are spectrally unmixed using the
appropriate 4$\times$4 unmixing matrix to produce 4 images. Under ideal
conditions of perfect unmixing, only one of the 4 spectral reflectance targets
should contain power in any recovered frame. f) Unmixed images for a section
of the filter array where the hyperpixel outperforms the optimal band-pass
filter in terms of condition number, and also a section where hyperpixel
performance is lower than optimal bandpass. g) Over 10 different regions of
the filter array we find that the condition number of the unmixing matrix is
significantly reduced for hyperpixels compared to bandpass (ratio = 3.47x, 95%
CI = 1.28 – 9.40, $p$ = 0.020).

Due to variations in fabrication quality across the filter array, we find that
relative benefit of the hyperpixels compared to bandpass filters varies. We
present an example with greater reduction in condition number (Figure 5f).
Though reconstruction is still far from ideal, we observe that for
hyperpixels, more power is directed to the appropriate pixel than for the
bandpass filters, in which the power is largely directly to the first 2 pixels
only (particularly in the 3rd and 4th columns of Figure 5f). This visible
benefit reduces somewhat in a ‘bad’ example (Figure 5g), where the hyperpixels
exhibits lower condition number than the bandpass. However, we still note that
the hyperpixels perform visibly better than the bandpass filters, e.g. in the
4th column of Figure 5 where more power appears in the bottom right square.

To robustly quantify any change in condition number in spite of this
variation, we perform a $t$-test across 10 different measured areas of the
filter array (Figure 5h). We find a significant reduction in condition number
for hyperpixels compared to bandpass filters (ratio = 3.47x, 95% CI = 1.28 –
9.40, $p$ = 0.020).

##  4 Discussion

We have presented a novel approach to matched filter design for pixel filter
arrays that utilizes numerous small ‘subpixel filter’ Fabry-Perot resonators
within an image sensor pixel area to create ‘hyperpixels’, whose spectral
response matches that of target spectra in a sample. Though our proposed
design operates under a fixed-power constraint—since any given region can
still only transmit one spectral band—there should still theoretically be a
discernible advantage of this approach when performing spectral unmixing as it
better exploits the broad spectral shape of many realistic targets. Indeed,
our experimental tests with fabricated hyperpixel filter arrays show, using
both measured spectra and imaging, that hyperpixels reduce the condition
number of the spectral unmixing matrix at least 2-3 times compared to optimal
bandpass filters. We have thus shown that these designs provide a measurable
practical benefit. Considering that our experimental filter arrays are non-
optimal (i.e. mismatch with designed structures), our results suggest the
potential impact under optimized conditions could be far greater.

To progress the hyperpixel approach towards real-world application, there are
some challenges that must be addressed: The first challenge is achieving more
uniform performance across the fabricated sample. This could be improved by
tighter control of the fabrication process, including exposure and development
steps for the grayscale lithography aspect. Using a parallelized exposure
scheme such as mask-based photolithography with either a single grayscale mask
or a stepper (with dose modulation) may improve uniformity and is also
scalable to larger-volume manufacturing.

The second challenge is to modify the design process to better optimize for
the the fixed power constraint. Due to the division of pixel area into
spectral bands, there is inherently a fixed power constraint: the integrated
area under the transmission spectrum of a hyperpixel is equal to that of a
bandpass filter. However, previous simulation work has demonstrated that even
under this fixed power constraint, a $>$ 5-fold improvement in signal to noise
ratio (SNR) under realistic noise conditions can be achieved [27]. If there
were no power constraint, we would expect hyperpixels to always perform equal
or better than bandpass, with the limiting case being that the two are
equivalent (e.g. for bandpass spectral reflectance targets). However, because
of the fixed power constraint of our design there are scenarios in which
bandpass filters perform better [27]. Conversely, there are also many
scenarios where the hyperpixels, even under the total power constraint,
outperform optimal bandpass filters and these will vary depending on the
target.

Our design process for hyperpixels relies on the same assumptions used in
designing matched filters, which a fixed total power constraint. Though this
offers simplicity and clearly provides a measurable benefit (as demonstrated),
alternative approaches such as water-pouring algorithm or a generic full-
spectrum optimization may further maximize condition number [29].

Linked to this is the third challenge: improving power transmittance. This
could be achieved by using high-index all-dielectric metasurfaces to provide
the tailored spectral response, though scaling such approaches to
manufacturability scale may be challenging. However, approaches such as
nanoimprint lithography make this option increasingly feasible for production.
Spectral responses could also be achieved through multilayer reflectors, as
used for multivariate optical elements, although this may require numerous
additional fabrication steps to create arrays [30, 31, 32, 33].

In terms of applications, we note that hyperpixels do not provide benefit in
all scenarios. For example if the expected spectral response of the targets
resemble bandpass filters, e.g. detecting fluorescence from multiple
fluorophores, the optimal hyperpixel design will be similar to the bandpass
filter design. However, the hyperpixel approach offers advantages in specific
applications where the target spectra span large spectral bandwidths. Because
of the high customisability, we anticipate hyperpixels could find use for low-
volume applications where real-time performance is needed and unmixing is
typically noisy. By using techniques such as photolithography or nanoimprint
lithography, one could make standalone filter arrays in moderate quantities to
be bonded to pre-packaged monochrome CMOS sensors rather than fabricating
directly on sensors during manufacturing. Potential applications may include
machine vision where high volumes of items with _a priori_ known spectra are
encountered, or biomedical spectral imaging where common analyte spectra such
as oxygenated and deoxygenated hemoglobin are imaged across a wide range of
applications.

##  5 Conclusion

Overall, we conclude that hyperpixel filter arrays can provide a measurable
performance benefit over optimal band-pass filter arrays despite the
aforementioned fabrication limitations. They are therefore a highly promising
technology for bespoke spectral imaging systems.

Funding The authors acknowledge support from a UKRI Future Leaders Fellowship
(MR/T041951/1, GSDG), Cancer Research UK (A9545/A29580, SEB) and EPSRC
(EP/R003599/1, SEB), and a NanoPrime grant from the Nottingham Nano- and
Micro-scale Research Centre (EP/R025282/1).

Acknowledgments The authors thanks the Nottingham Nano- and Micro-scale
Research Centre (nmRC) for fabricating all samples.

Disclosures The authors declare no conflicts of interest.

Data Availability Statement Data in this study available at [34].

## References

  * [1] Q. Li, X. He, Y. Wang, H. Liu, D. Xu, and F. Guo, “Review of spectral imaging technology in biomedical engineering: Achievements and challenges,” 18, 100901.
  * [2] G. Lu and B. Fei, “Medical hyperspectral imaging: A review,” 19, 010901.
  * [3] J. Spigulis, D. Jakovels, and L. Elste, “Towards single snapshot multispectral skin assessment,” p. 82160L.
  * [4] J. R. Bauer, A. A. Bruins, J. Y. Hardeberg, and R. M. Verdaasdonk, “A Spectral Filter Array Camera for Clinical Monitoring and Diagnosis: Proof of Concept for Skin Oxygenation Imaging,” 5, 66.
  * [5] D. J. Waterhouse, C. R. Fitzpatrick, B. W. Pogue, J. P. O’Connor, and S. E. Bohndiek, “A roadmap for the clinical implementation of optical-imaging biomarkers,” 3, 339–353.
  * [6] N. T. Clancy, G. Jones, L. Maier-Hein, D. S. Elson, and D. Stoyanov, “Surgical spectral imaging,” 63, 101699.
  * [7] A. Siri Luthman, D. Waterhouse, L. Ansel-Bollepalli, J. Yoon, G. Gordon, J. Joseph, M. Di Pietro, W. Januszewicz, and S. Bohndiek, “Bimodal reflectance and fluorescence multispectral endoscopy based on spectrally resolving detector arrays,” 24.
  * [8] Y. Liu, H. Pu, and D.-W. Sun, “Hyperspectral imaging technique for evaluating food quality and safety during various processes: A review of recent applications,” 69, 25–35.
  * [9] K. Clarke, A. Hennessy, A. McGrath, R. Daly, S. Gaylard, A. Turner, J. Cameron, M. Lewis, and M. B. Fernandes, “Using hyperspectral imagery to investigate large-scale seagrass cover and genus distribution in a temperate coast,” 11, 4182.
  * [10] T. Adão, J. Hruška, L. Pádua, J. Bessa, E. Peres, R. Morais, and J. Sousa, “Hyperspectral Imaging: A Review on UAV-Based Sensors, Data Processing and Applications for Agriculture and Forestry,” 9, 1110.
  * [11] C. Williams, G. Gordon, S. Gruber, T. Wilkinson, and S. Bohndiek, “Grayscale-to-color: Single-step fabrication of bespoke multispectral filter arrays,” .
  * [12] T. W. Sawyer, M. Taylor-Williams, R. Tao, R. Xia, C. Williams, and S. E. Bohndiek, “Opti-MSFA: A toolbox for generalized design and optimization of multispectral filter arrays,” 30, 7591.
  * [13] R. Wu, Y. Li, X. Xie, and Z. Lin, “Optimized Multi-Spectral Filter Arrays for Spectral Reconstruction,” 19, 2905.
  * [14] G. Turin, “An introduction to matched filters,” 6, 311–329.
  * [15] S. Cianci, J. Bland-Hawthorn, and J. W. O’Byrne, “Multiband interference filters for differential imaging,” 4008, 1368.
  * [16] M. P. Nelson, J. F. Aust, J. A. Dobrowolski, P. G. Verly, and M. L. Myrick, “Multivariate Optical Computation for Predictive Spectroscopy,” 70, 73–82.
  * [17] S. Kozhukh, S. Zhukovsky, I. Yurkevich, S. Stepkova, S. Gaponenko, A. Lunevich, and Y. Glukhov, “Design of high-transmission multiband multilayer filters for Raman spectroscopy,” 6, 061704.
  * [18] J. Xiong, X. Cai, K. Cui, Y. Huang, J. Yang, H. Zhu, W. Li, B. Hong, S. Rao, Z. Zheng, S. Xu, Y. He, F. Liu, X. Feng, and W. Zhang, “Dynamic brain spectrum acquired by a real-time ultraspectral imaging chip with reconfigurable metasurfaces,” 9, 461.
  * [19] A. Ahamed, C. Bartolo-perez, A. S. Mayet, S. Ghandiparsi, S.-y. Wang, and M. S. Islam, “Reconstruction-based spectroscopy using CMOS image sensors with random photon-trapping nanostructure per sensor,” .
  * [20] M. S. Mirotznik, “Photonic crystal filters for multi-band optical filtering on a monolithic substrate,” 3, 031506.
  * [21] X. Zou, Y. Zhang, R. Lin, G. Gong, S. Wang, S. Zhu, and Z. Wang, “Pixel-level Bayer-type colour router based on metasurfaces,” 13, 1–7.
  * [22] K. Monakhova, K. Yanny, N. Aggarwal, and L. Waller, “Spectral DiffuserCam: Lensless snapshot hyperspectral imaging with a spectral filter array,” 7.
  * [23] A. Wagadarikar, R. John, R. Willett, and D. Brady, “Single disperser design for coded aperture snapshot spectral imaging,” 47, B44.
  * [24] M. Yako, Y. Yamaoka, T. Kiyohara, C. Hosokawa, A. Noda, K. Tack, N. Spooren, T. Hirasawa, and A. Ishikawa, “Video-rate hyperspectral camera based on a CMOS-compatible random array of Fabry–Pérot filters,” 17, 218–223.
  * [25] C. Williams, G. S. D. Gordon, T. Wilkinson, and S. E. Bohndiek, “Snapshot multispectral imaging using metal-insulator-metal filter arrays fabricated using grayscale lithography,” .
  * [26] T. Goossens, “Pixel-integrated thin-film filter simulation and scaling trade-offs,” .
  * [27] M. Taylor-Williams, R. B. Cousins, C. Williams, S. E. Bohndiek, C. J. Mellor, and G. S. D. Gordon, “Spectrally tailored ’hyperpixel’ filter arrays for imaging of chemical compositions,” in _Optical Biopsy XX: Toward Real-Time Spectroscopic Imaging and Diagnosis,_ R. R. Alfano, S. G. Demos, and A. B. Seddon, eds. (SPIE), p. 23.
  * [28] L. Gao and R. T. Smith, “Optical hyperspectral imaging in microscopy and spectroscopy - a review of data acquisition,” 8, 441–456.
  * [29] A. Goldsmith and P. Varaiya, “Capacity of fading channels with channel side information,” 43, 1986–1992 (Nov./1997).
  * [30] D. Yi, L. Kong, J. Wang, and F. Zhao, “Fabrication of multispectral imaging technology driven MEMS-based micro-arrayed multichannel optical filter mosaic,” p. 792711.
  * [31] C. M. Hooper, S. E. Bohndiek, and C. Williams, “Optimizing achromaticity in metalenses, and development of a layered thin-film metalens,” in _High Contrast Metastructures XI,_ J. A. Fan, C. J. Chang-Hasnain, and W. Zhou, eds. (SPIE), p. 31.
  * [32] S.-W. Wang, X. Chen, W. Lu, L. Wang, Y. Wu, and Z. Wang, “Integrated optical filter arrays fabricated by using the combinatorial etching technique,” 31, 332.
  * [33] O. Soyemi, D. Eastwood, L. Zhang, H. Li, J. Karunamuni, P. Gemperline, R. A. Synowicki, and M. L. Myrick, “Design and Testing of a Multivariate Optical Element: The First Demonstration of Multivariate Optical Computing for Predictive Spectroscopy,” 73, 1069–1079.
  * [34] https://doi.org/10.17639/nott.7397


# A Case Study Analysis for Designing a Lunar Navigation Satellite System with
Time-Transfer from Earth-GPS

Sriramya Bhamidipati  Tara Mina and Grace Gao  Stanford University

## biography

Sriramya Bhamidipatiis a postdoctoral scholar in the Aeronautics and
Astronautics Department at Stanford University. She received her Ph.D. in
Aerospace Engineering at the University of Illinois, Urbana-Champaign in 2021,
where she also received her M.S in 2017. She obtained her B.Tech. in Aerospace
from the Indian Institute of Technology, Bombay in 2015. Her research
interests include GPS, power and space systems, artificial intelligence,
computer vision, and unmanned aerial vehicles.

Tara Minais a Ph.D. candidate in the Electrical Engineering Department at
Stanford University. She received her B.S. degree in Electrical Engineering
from Iowa State University in 2017 and completed her M.S. degree in Electrical
and Computer Engineering from the University of Illinois at Urbana-Champaign
in 2019.

Grace Gaois an assistant professor in the Department of Aeronautics and
Astronautics at Stanford University. Before joining Stanford University, she
was an assistant professor at University of Illinois at Urbana-Champaign. She
obtained her Ph.D. degree at Stanford University. Her research is on robust
and secure positioning, navigation and timing with applications to manned and
unmanned aerial vehicles, robotics, and power systems.

## Abstract

Recently, there has been a growing interest in the use of a SmallSat platform
for the future Lunar Navigation Satellite System (LNSS) to allow for cost-
effectiveness and rapid deployment. However, many design choices are yet to be
finalized for the SmallSat-based LNSS, including the onboard clock and the
orbit type. As compared to the legacy Earth-GPS, designing an LNSS poses
unique challenges: (a) restricted Size, Weight, and Power (SWaP) of the
onboard clock, which limits the timing stability; (b) limited lunar ground
monitoring stations, which engenders a greater preference toward stable LNSS
satellite orbits.

In this current work, we analyze the trade-off between different design
considerations related to the onboard clock and the lunar orbit type for
designing an LNSS with time-transfer from Earth-GPS. Our proposed time-
transfer architecture combines the intermittently available Earth-GPS signals
in a timing filter to alleviate the cost and SWaP requirements of the onboard
clocks. Specifically, we conduct multiple case studies with different grades
of low-SWaP clocks and various previously studied lunar orbit types. We
estimate the lunar User Equivalent Range Error (UERE) metric to characterize
the ranging accuracy of signals transmitted from an LNSS satellite. Using the
Systems Tool Kit (STK)-based simulation setup from Analytical Graphics, Inc.
(AGI), we evaluate the lunar UERE across various case studies of the LNSS
design to demonstrate comparable performance as that of the legacy Earth-GPS,
even while using a low-SWaP onboard clock. We further perform sensitivity
analysis to investigate the variation in the lunar UERE metric across
different case studies as the Earth-GPS measurement update rates are varied.

##  1 Introduction

Over a half century since the end of the Apollo program, NASA will return
humans to the Moon. In the coming decade, through the Artemis missions,
National Aeronautics and Space Administration (NASA) will land the first woman
and person of color on the lunar south pole [1]. Additionally, many space
agencies are participating in an international effort to establish a
sustainable human presence on the Moon, which serves as a crucial platform to
support future deep space exploration [2]. Indeed, we are entering a second
space race, with international space agencies planning over forty lunar
missions in the next decade and with crucial involvement from the commercial
space industry companies, including SpaceX and Blue Origin. To support the
increasing plans for crewed and robotic activities, future lunar missions will
require access to reliable and precise Position, Navigation, and Timing (PNT)
services everywhere on the Moon.

Recently, the NASA Goddard Space Flight Center (GSFC) and the European Space
Agency (ESA) have conceptualized GPS-like satellite constellations around the
Moon, named LunaNet [3] and Moonlight [4], respectively. Furthermore, there
has been an emerging interest in the use of a SmallSat platform for these PNT
constellations to allow for cost-effectiveness and rapid deployment [3]. These
lunar PNT constellations by NASA and the ESA will assist in the overarching
effort of establishing a sustainable human presence on the Moon, by providing
global PNT and communication services to lunar users. In particular, in the
next decade, these initiatives seek to satisfy needs expressed by the Global
Exploration community, with a targeted position accuracy of less than $50$
$\mathrm{m}$ for lunar users [4].

Given that the lunar PNT constellation initiatives are in the preliminary
design phases, the design a Lunar Navigation Satellite System (LNSS) involves
finalizing many key design considerations, including:

  * •

Lunar satellite orbit. Several types of lunar orbits that have previously been
investigated include Low Lunar Orbit (LLO), Prograde Circular Orbit (PCO),
Near-Rectilinear Halo Orbit (NRHO), and Elliptical Lunar Frozen Orbit (ELFO).
In particular, ELFOs refer to a specific category of frozen orbits providing a
greater coverage of the lunar poles, wherein frozen orbits represent the
orbits that maintain nearly constant orbital parameters for extensive periods
of time, without requiring station keeping [5, 6]. Although PCOs are not
frozen, they maintain multi-year stability with orbital parameters exhibiting
predictable, repeatable behavior [6]. NRHOs are highly elliptical orbits with
nearly constant visibility of Earth and the lunar poles [7]. Due to the low
orbiting altitude, LLOs have shorter orbital periods around the Moon, and
there exist a few inclinations in which LLOs are also considered to be frozen
or quasi-frozen [5].

  * •

Onboard clock. The choice of onboard clock is critical for designing a
navigation system as its grade (which depends on timing stability) directly
affects the ranging precision offered to lunar users. Among various clock
choices is the commercial Chip Scale Atomic Clock (CSAC) with radiation-
tolerance and low Size, Weight and Power (SWaP), which has been specifically
developed for space applications [8]. Another potential clock choice is the
Deep Space Atomic Clock (DSAC), which has been recently designed by NASA to
provide greater long-term timing stability and to assist in spacecraft
radionavigation [9].

Furthermore, designing a SmallSat-based LNSS involves unique challenges as
compared to the legacy Earth-GPS, which lead to additional design limitations
listed as follows: a) Limited size of LNSS satellites. A SmallSat platform
limits its payload capacity, including the SWaP of the onboard clock. Given
that lower SWaP clocks tend to have worse timing stabilities [8], the SWaP
limitation on clock directly affects the timing stability; b) Limited ability
to monitor LNSS satellites. Given that a limited number of ground monitoring
stations can be established on the Moon and that resources on Earth for
monitoring the lunar constellation are limited, it is desirable for the LNSS
satellites to require less maintenance, including fewer station-keeping
maneuvers and clock correction maintenance; and c) Increased orbital
perturbations in lunar environment. Because the Moon has a highly non-uniform
distribution of mass [10], its gravitation field is more anisotropic as
compared to that of Earth’s. In addition, Earth’s gravity can significantly
impact satellites in high-altitude orbits around the Moon, thereby limiting
the set of feasible and stable lunar orbits.

Given these challenges for designing a PNT constellation in the lunar
environment, one may consider the potential in leveraging the existing Earth’s
legacy GPS, which is equipped with higher grade atomic clocks and an extensive
ground monitoring network. At lunar distances of about $385000$
$\mathrm{k}\mathrm{m}$, the Earth-GPS signal is significantly attenuated, and
the Earth-GPS satellites directed toward Earth are largely occluded by Earth
and often the Moon. This limits the Earth-GPS signal availability at lunar
distances to come from only the Earth-GPS transmit antenna’s side lobes and
the small, unoccluded parts of the main lobe. NASA’s magnetospheric multiscale
mission, or MMS, has used these largely attenuated and intermittently
available Earth-GPS signals to successfully compute position estimate in space
[11]. In fact, the MMS broke the Guinness World Record for the highest
altitude of achieving an Earth-GPS fix in 2016 for traveling at distances of
about one-fifth of the way to the Moon [12, 11] and then again in 2019 for
distances of about halfway to the Moon [13]. Several simulated works also have
demonstrated the feasibility of using Earth-GPS at lunar distances [7, 14,
15]. Through GPS Antenna Characterization Experiment (ACE) study, NASA has
characterized the GPS antenna gain patterns at high elevation angles from
boresight for space users [16]. NASA has also developed a spaceborne Earth-GPS
receiver [17], which will be tested on the lunar surface for the first time in
2023 [18, 19]. Additionally, in 2023, the ESA will launch the Lunar Pathfinder
communication satellite to the Moon, which will utilize a spaceborne, high-
sensitivity Earth-GPS receiver to provide a position fix for the first time in
lunar orbit [20].

In our prior work [21], we designed an LNSS architecture that harnessed the
legacy Earth-GPS system to provide precise timing corrections to the onboard
clock, as depicted in Fig. 1. Our proposed time-transfer technique leverages
intermittently available Earth-GPS signals to alleviate the SWaP requirements
of the onboard clocks and to mitigate the need for an extensive ground
monitoring infrastructure on the Moon. We also designed a lunar User
Equivalent Range Error (UERE) metric, which is proportional to the Root-Mean-
Square (RMS) timing error, to analyze the ranging accuracy of an LNSS
satellite. We further demonstrated that our proposed method achieves a low
UERE of less than $10$ $\mathrm{m}$ while using a low-SWaP Chip Scale Atomic
Clock (CSAC) for an LNSS satellite in an Elliptical Lunar Frozen Orbit (ELFO).

Given that many design choices, including the grade of the onboard clock and
the orbit type, still need to be finalized for the future SmallSat-based LNSS,
in this work, we extend our prior work [21] to analyze the LNSS performance
using our proposed time-transfer architecture (in Fig. 1) from Earth-GPS under
various case studies. Specifically, we simulate an LNSS satellite in various
lunar orbit types, including the ELFO, LLO, PCO, and NRHO, and equipped with
different grades of onboard clocks. For a given LNSS satellite orbit, we
examine an Earth-GPS Continual Outage Period (ECOP) metric to analyze the
visibility effects of Earth-GPS on the timing stability of the onboard clock.
We further estimate the lunar UERE metric to perform comparison across
different case studies. Through our analysis, we illustrate the trade-off
between the different design considerations of the onboard clock and the orbit
type for an LNSS design that leverages Earth-GPS time-transfer. We
additionally examine how the lunar UERE is affected by different rates of
collecting available Earth-GPS measurements. This analysis provides insights
into the extent of Earth-GPS signal tracking and processing required to
provide sufficient ranging precision. In particular, less frequent use of
Earth-GPS measurements would allow the LNSS satellite to continually switch
off the onboard Earth-GPS receiver for longer periods of time in order to save
power. Across the various case studies, we observe that with our time-transfer
architecture, the LNSS can achieve comparable performance as that of the
legacy Earth-GPS, even while using a low-SWaP onboard clock.

Figure 1: Architecture of our proposed time-transfer from Earth-GPS [21],
which utilizes intermittently available Earth-GPS signals to correct the lower
grade clocks onboard the LNSS satellite.

###  1.1 Key Contributions

The key contributions of our paper are listed as follows:

  1. 1.

We design various case studies related to the grade of onboard clocks and the
orbit types for analyzing the trade-offs in designing an LNSS with time-
transfer from Earth-GPS. In particular, we investigate five clock types with
diverse SWaP characteristics that range from a low-SWaP CSAC to a high-SWaP
Deep Space Atomic Clock (DSAC) [9] developed by NASA. We also investigate four
previously studied lunar orbit types, namely ELFO, NRHO, LLO and PCO.

  2. 2.

We perform a comparison analysis across various case studies by investigating
the associated RMS timing errors, which are in-turn governed by the duration
for which no Earth-GPS satellites are visible (ECOP metric) and the geometric
configuration between the Earth-GPS constellation and the LNSS satellite
(occultations due to Earth and the Moon).

  3. 3.

We evaluate the lunar UERE metric for each case study and demonstrate
comparable measurement ranging accuracy as that of the legacy Earth-GPS, even
while using a low-SWaP onboard clock.

  4. 4.

We additionally examine how the lunar UERE is affected by different rates of
collecting Earth-GPS measurements. This analysis provides insight into the
extent of Earth-GPS signal processing required, and, correspondingly, the
amount of power required to operate the onboard Earth-GPS receiver, in order
to provide sufficient ranging precision

The rest of the paper is organized as follows. Section 2 summarizes our
previously proposed time-transfer architecture from Earth-GPS to LNSS [21] and
later describes the modifications incorporated in our time-transfer to conduct
further analysis. Section 3 provides a high-level overview of various case
studies and also describes our high-fidelity lunar simulation setup that
involves modeling the onboard clock and the orbit for each case study. Section
4 discusses the implications of our case study analysis in designing an LNSS.
Section 5 provides concluding remarks.

##  2 Our Time-Transfer from Earth-GPS to LNSS

In this section, we summarize our prior work [21] on time-transfer from Earth-
GPS, wherein, we consider an LNSS satellite equipped with a Earth-GPS receiver
and an onboard clock that can provide short-term timing stability. We designed
a timing Kalman filter [22, 23] to update the LNSS satellite clock with the
intermittently available Earth-GPS signals and formulated the lunar UERE
metric to characterize the ranging accuracy of transmitted navigation signals.
We additionally provide an overview on the aspects of further analysis
conducted in this work. In particular, we analyze the sensitivity of the lunar
UERE metric in different simulated case studies, including modifications of
the measurement update rate for our proposed timing Kalman filter.

At any LNSS satellite, our proposed filter maintains the LNSS clock estimate
at each time epoch $t$ by propagating the following timing state vector:
$x_{t}\coloneqq\begin{bmatrix}b_{t}&\dot{b}_{t}\end{bmatrix}^{\top}$, where
$b_{t}$ is the clock bias state in $\mathrm{m}$ and $\dot{b}_{t}$ is the clock
drift in $\mathrm{m}\mathrm{s}^{-1}$, with units converted from the timing
domain through multiplication by the speed of light $c=299792458$ m/s.

To maintain the LNSS clock estimate, our timing Kalman filter performs a time
update every $T_{\text{pred}}$ seconds, based on the clock error propagation
model. For this, we define the associated process noise covariance $Q$ in
terms of the power spectral density coefficients $h_{0}$, $h_{-1}$, $h_{-2}$
from the clock Allan deviation plot [22]. These Power Spectral Density (PSD)
coefficients reflect the short-term and long-term stability of the onboard
clock [24]. To perform time-transfer from Earth-GPS, we first determine if any
Earth-GPS signals are visible by examining the received carrier-to-noise-
density ratio $C/N_{0}$. Then, our timing Kalman filter conducts a measurement
update for the available Earth-GPS measurements with sufficiently high
$C/N_{0}$.

During the measurement update step, we determine the expected pseudorange and
pseudorange rates from the visible Earth-GPS satellites to form a measurement
vector of residuals. In particular, we formulate the measurement residual
vector by leveraging the position aiding from available ephemeris information
for the LNSS satellite. We additionally model the measurement covariance
matrix as a time-dependent diagonal matrix [21], based on the tracking errors
of the receiver Delay Lock Loop (DLL) [25, 26, 27] and Phase Lock Loop (PLL)
[28, 27] as well as the Earth-GPS UERE [25] and the expected error in
available LNSS satellite ephemeris. With our measurement vector and modeled
measurement covariance, our filter applies corrections to the predicted timing
state via the standard Kalman filter expressions to obtain the updated state
$\hat{x}_{t}$ and covariance.

Based on the RMS error in our filter estimate, we formulate a lunar UERE
metric that characterizes the accuracy of the LNSS ranging signals for lunar
users. On the Moon, any atmospheric delays are minimal. Moreover, we consider
multipath effects experienced by users on the lunar surface to be negligible
due to the lack of building infrastructure and foliage. As a result, the final
lunar UERE can be computed in terms of four most significant error components
as follows:

$\displaystyle{\sigma_{\text{UERE, LNSS}}}=\sqrt{\sigma_{\text{clk,
LNSS}}^{2}+\sigma_{\text{gd, LNSS}}^{2}+\sigma_{\text{eph,
LNSS}}^{2}+\sigma_{\text{rec, LNSS}}^{2}},$ (1)

where the errors due to the differential group delay $\sigma_{\text{gd,
LNSS}}$ and the receiver noise $\sigma_{\text{rec, LNSS}}$ will depend on the
final LNSS signal structure and lunar user receiver. Note that, because the
timing filter uses position aiding from the LNSS satellite, the lunar
ephemeris error component $\sigma_{\text{eph, LNSS}}$ directly impacts the
pseudorange residual measurement received at the LNSS, which will thus also
affect the LNSS clock error $\sigma_{\text{clk, LNSS}}$.

In our current work, we investigate how the LNSS clock error component
$\sigma_{\text{clk, LNSS}}$ differs for various grades of onboard clocks and
for various types of lunar orbits, while also analyzing the corresponding
impact on the overall lunar UERE. We additionally investigate the impact to
the lunar UERE for a potentially reduced measurement update rate when Earth-
GPS signals are available, with a sampling period of
$T_{\text{meas}}=mT_{\text{pred}}$ seconds, where $m$ is a positive integer.
Indeed, a larger choice of $m$ corresponds to less frequent GPS measurement
updates, which allows the spaceborne Earth-GPS receiver onboard an LNSS
satellite to be switched off for longer durations of time to save power.

##  3 Overview of Case Studies on Clocks and Orbits

We perform an extensive case study analysis to examine the trade-off between
different choices of the onboard clocks and the orbit types that can be
considered for designing an LNSS with time-transfer from Earth-GPS. In
particular, we design a high-fidelity simulation of an LNSS satellite for each
orbit type using the Systems Tool Kit (STK) software by the Analytical
Graphics, Inc. (AGI) [29]. For each modeled orbit type, we design case studies
in MATLAB by simulating various grades of onboard clocks. For each case study,
we consider the start time epoch to be $9$ Nov $2025$ $00:00:00.000$ UTC and
the experiment time duration to be $2$ months (which equals $61$ days).

We first provide an overview of case studies related to the onboard clocks and
the lunar orbit types investigated in this work. Thereafter, we summarize the
simulation steps executed in the STK software and MATLAB for modeling the
transmission and reception of Earth-GPS signals in each case study, which is
based on the validation framework designed for our earlier work [21].

As mentioned in Section 1, many prior studies [30, 31] have investigated
various types of lunar orbits, primarily based on their stability and the
duration for which their stability can be ensured. Our choice of lunar orbit
types for this work include ELFO, NRHO, LLO and PCO. We create realistic
simulations of an LNSS satellite in different lunar orbits by leveraging the
High Precision Orbit Propagator (HPOP) in the STK software [29]. The HPOP
generates and propagates accurate position and velocity solutions of the LNSS
satellite by accounting for precise force models of Earth, the Sun and the
Moon [32].

Table 1: Keplarian parameters of the three (ELFO, LLO and PCO) among the four lunar orbit types considered in our case study analysis. | Orbit Type
---
| Altitude
---
($\mathrm{k}\mathrm{m}$)
Eccentricity |  | Inclination
---
(∘)
| Argument of
---
Perigee (∘)
| RAAN
---
(∘)
| Mean Anomaly
---
(∘)
ELFO | 9750.5 | 0.7 | 63.5 | 90 | 0 | 0
LLO | 100 | 0 | 28.5 | 0 | 0 | 0
PCO | 3000 | 0 | 75 | 0 | 90 | 0

We model the three orbit types, namely ELFO, LLO and PCO, in the STK software
using classical orbit mechanics [33], according to which objects orbiting in
space requires six elements (in our case six Keplarian parameters) to fully
characterize their position and velocity at any point in time. Specifically,
we refer to prior literature [34, 35, 6] on ELFO, LLO and PCO to define their
corresponding six Keplarian parameters at the start time epoch, which includes
semi-major axis, eccentricity, inclination, argument of perigee, Right
Ascension of the Ascending Node (RAAN), and mean anomaly. Table 1 lists the
associated Keplarian parameters of the three lunar orbit types, while Fig. 2
shows their associated illustration in the Moon-inertial frame.

Figure 2: Illustration of the three (ELFO, LLO and PCO) among the four lunar
orbit types considered in our case study analysis. All these orbits are
visualized in the Moon-inertial frame, where ELFO indicated in orange is at an
altitude of $9750.5$ $\mathrm{k}\mathrm{m}$, LLO represented in magenta is at
$100$ $\mathrm{k}\mathrm{m}$, and PCO represented in blue is at $3000$
$\mathrm{k}\mathrm{m}$.

(a) NRHO orbit type in Moon-inertial frame

(b) NRHO orbit type in Moon-centered Earth-Moon rotating frame

Figure 3: The NRHO (L2, radius of perigee=$4500$ $\mathrm{k}\mathrm{m}$, South
family) visualized in (a) the Moon-Inertial frame; and (b) the Moon-Centered
Earth-Moon rotating frame. The x-axis of the Earth-Moon rotating frame is
along the instantaneous Earth-Moon position vector, the z-axis is along the
instantaneous angular momentum vector of the Moon’s orbit around the Earth,
and the y-axis completes the orthogonal system.

In contrast, for modeling the fourth lunar orbit type in our case study
analysis, namely L2 South NRHO, we specify the initial conditions of position
$[r_{x},r_{y},r_{z}]$, and velocity $[v_{x},v_{y},v_{z}]$ in the Moon-centered
Earth-Moon rotating frame. Defining the initial conditions in an Earth-Moon
rotating frame is particularly useful when discussing halo orbits in the
Earth-Moon system [36]. In the Earth-Moon rotating frame, the x-axis is along
the instantaneous Earth-Moon position vector, the z-axis is along the
instantaneous angular momentum vector of the Moon’s orbit around the Earth,
and the y-axis completes the orthogonal system. Specifically, we refer to
prior literature [36] that solves for an NRHO in the ephemeris model (L2,
radius of perigee to be $4500$ $\mathrm{k}\mathrm{m}$, South family) using a
forward/backward shooting process to provide the following initial state
vector at an initial time epoch of $8$ Nov $2025$ $23:22:07.10353~{}$TDB:
$r_{x}=-125.952$ $\mathrm{k}\mathrm{m}$, $r_{y}=120.961$
$\mathrm{k}\mathrm{m}$, $r_{z}=4357.681~{}$km, $v_{x}=-0.042$
$\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}$, $v_{y}=1.468$
$\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}$ and $v_{z}=-0.003$
$\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}$. The illustrations of the designed
NRHO both in the Moon-inertial frame (similar frame as that in Fig. 2) and in
the Moon-Centered Earth-Moon rotating frame are shown in Figs. 3(a) and 3(b),
respectively.

Given the interest in SmallSat platform for the future LNSS, we choose various
case studies of onboard clock types while keeping in mind that the limited
payload capacity restricts the SWaP of the onboard clock. Based on prior
literature [8], our choice of clock types for this work include Microchip Chip
Scale Atomic Clock (CSAC), Microchip Micro Atomic Clock (MAC), Stanford
Research Systems (SRS) PRS-10, Excelitas Rubidium Atomic Frequency Standard
(RAFS) and NASA’s Deep Space Atomic Clock (DSAC) [9]. The specifications of
these clock types are listed in Table 2, which have been arranged in the
increasing order of their SWaP for convenience. For each clock type, we
simulate the true clock error model in MATLAB to have a constant drift in the
clock bias and thereafter, propagate the clock states forward in time using a
first-order state transition matrix. We assign the true clock drift (constant
value) based on the known specifications of Time Deviation (TDEV) observed at
the end of a day, which are reported in Table 2. Note that, for any clock
type, TDEV refers to the expected error in reported time after a certain
holdover time that essentially depends on their respective values of Allan
deviation and frequency drift. For time update in our timing filter described
previously in Section 2, we utilize the PSD coefficients listed in Table 2 to
model the corresponding process noise covariance matrix $Q$.

Table 2: SWaP characteristics, time deviation and PSD coefficients of the five clock types considered in our case study analysis. SWaP and TDEV values are taken from [8] while PSD coefficients were computed from their respective Allan deviation plots. Clock type |  Size ($\mathrm{c}\mathrm{m}^{3}$) |  Weight ($\mathrm{k}\mathrm{g}$) |  Power ($\mathrm{W}$) |  TDEV per day ($\mu\mathrm{s}$) | PSD coefficients
---|---|---|---|---|---
$h_{0}$ | $h_{-1}$ | $h_{-2}$
Microchip CSAC | 17 | 0.035 | 0.1 | 1.5$\times 10^{-6}$ | 1.3$\times 10^{-20}$ | 1.0$\times 10^{-24}$ | 3.7$\times 10^{-29}$
Microchip MAC | 50 | 0.084 | 5 | 1.7$\times 10^{-7}$ | 4.7$\times 10^{-22}$ | 1.2$\times 10^{-25}$ | 1.7$\times 10^{-30}$
SRS PRS 10 | 155 | 0.6 | 14.4 | 7.0$\times 10^{-8}$ | 1.3$\times 10^{-22}$ | 2.3$\times 10^{-26}$ | 3.3$\times 10^{-31}$
Excelitas RAFS | 1645 | 6.35 | 39 | 4.8$\times 10^{-9}$ | 8.0$\times 10^{-24}$ | 0 | 0
NASA DSAC | 17000 | 16 | 47 | 4.0$\times 10^{-11}$ | 1.8$\times 10^{-27}$ | 0 | 0

For a given LNSS satellite clock and orbit type in each case study, we model
the simulation scenario in the STK software and MATLAB to compute the
following statistics that are later given as input to our timing filter: the
$C/N_{0}$ and measurement residual vector for visible Earth-GPS satellites. We
summarize the key modeling aspects below, while a more detailed explanation
regarding the same has been discussed in our earlier work [21].

First, we design our simulated Earth-GPS constellation to be made up of $31$
satellites with $8$ satellites from Block IIR, $7$ from IIRM, $12$ from IIF
and $4$ from Block III. We model the transmit antenna of Earth-GPS satellites
by utilizing the transmit power and antenna gain patterns of the L1 C/A
signals, which are made available from the NASA GPS ACE study [16]. Next, we
simulate a spaceborne Earth-GPS receiver with a steering antenna pointed
towards the Earth so as to maximize the visibility of Earth-GPS signals at the
LNSS satellite. Based on prior literature [35, 27], we consider a high-gain
antenna with $14$ $\mathrm{d}\mathrm{B}\mathrm{i}$ at $0^{\circ}$ off-
boresight angle, and a $3~{}$dB beamwidth of $12.2^{\circ}$. We consider an
Earth-GPS satellite to be visible when the received $C/N_{0}$ value is greater
than $15$ $\mathrm{d}\mathrm{B}\mathrm{-}\mathrm{H}\mathrm{z}$ for a
continuous time duration of at least $40$ $\mathrm{s}$.

Finally, we simulate received measurements at the LNSS satellite by
incorporating the true clock bias and drift in the true range and range rate
between the visible Earth-GPS and LNSS satellites, respectively. Note that the
$C/N_{0}$, true range and range rate values are extracted from the STK
simulation while the true clock bias and drift are obtained from the simulated
clock error model in MATLAB. To formulate the residual vector given to
measurement update explained previously in Section 2, we induce stochastic
errors based on the simulated uncertainties from the receiver tracking loops.

##  4 Our Case Study Analysis: Results and Discussion

For designing an LNSS with time-transfer from Earth-GPS, we analyze the trade-
off across different case studies related to the onboard clock and the orbit
type, which were listed earlier in Section 3.

To characterize the lunar UERE discussed in Section 2, we consider the same
group delay and receiver noise error magnitudes as with the Earth-GPS system,
i.e. $\sigma_{\text{gd, LNSS}}=0.15$ $\mathrm{m}$ and $\sigma_{\text{rec,
LNSS}}=0.1$ $\mathrm{m}$. Given that the future LNSS will have greater
limitations on ground monitoring infrastructure as compared to Earth-GPS, we
perform case study analysis by scaling the error component due to the
broadcast ephemeris in the lunar UERE as $\sigma_{\text{eph, LNSS}}=3$
$\mathrm{m}$. This value is essentially one order of magnitude higher in
comparison to that of Earth-GPS. We design our timing filter with the time
update executed every $T_{\text{pred}}=60$ $\mathrm{s}$. For the measurement
update, details regarding the reduced update rate with a sampling period of
$T_{\text{meas}}=mT_{\text{pred}}$ for measurement update will be discussed
below.

###  4.1 Validation Metrics

To perform comparison analysis across different case studies, we define the
following four validation metrics:

  1. 1.

Satellite visibility, which indicates the percentage time in the entire
experiment duration for which the number of visible Earth-GPS satellites are
greater than a pre-specified threshold. We examine two satellite visibility
conditions described as follows: a) the percent of time for which at least $1$
Earth-GPS satellite is visible, since this is the minimum number required to
estimate the clock bias and drift; b) the percent time for which at least $4$
Earth-GPS satellites are visible, since this is the minimum number to estimate
the full state vector that includes position, clock bias, velocity and clock
drift;

  2. 2.

Maximum ECOP to identify the region of maximum continuous time when no Earth-
GPS satellites are visible;

  3. 3.

Root-Mean-Square (RMS) errors in clock estimates to analyze the performance of
our time-transfer from Earth-GPS for reduced measurement update rate with a
sampling period of $T_{\text{meas}}=mT_{\text{pred}}$, where $m=5$; and

  4. 4.

Lunar UERE metric that characterize the ranging measurement accuracy of
signals transmitted by an LNSS satellite. Recollect from an earlier
explanation in Section 2 that the lunar UERE metric depends on the RMS error
in clock bias.

We consider a case study to be desirable if it exhibits either some or all of
the following: greatest satellite visibility, least maximum ECOP, and least
lunar UERE. To perform sensitivity analysis of our case studies, we compute
the lunar UERE metric for different reduced measurement update rates with
$m=1,5,30,60$. Note that, we perform $m=1$ for baseline comparison (not
reduced rate) since it depicts the case where the measurement update step is
executed whenever Earth-GPS satellites are visible, and thus depicts the
original framework of timing filter proposed in our prior work [21].

###  4.2 Across Orbit Types: Variation in Satellite Visibility and Maximum
ECOP

For the four orbit types considered in our case study, we showcase the number
of visible Earth-GPS satellites by blue in Figs. 4(a)-4(h) while the
highlighted red vertical bars indicate the regions of ECOP. Based on Table 3,
we observe that NRHO achieves the greatest at least one satellite visibility,
which is $99.9\%$ of the total time, and also the greatest at least four
satellite visibility, which is $92.9\%$ of the total time. Furthermore, NRHO
also exhibits the least maximum ECOP of only $420$ $\mathrm{s}$, while the
other orbit types experience ECOP of at least $2880$ $\mathrm{s}$. These
observations related to NRHO seem reasonable since an LNSS satellite in NRHO
operates at high altitudes ranging between $4500$
$\mathrm{k}\mathrm{m}$-$700000$ $\mathrm{k}\mathrm{m}$ above the Moon’s
surface, and thus experiences fewer occultations from Earth and the Moon.
Comparing Figs. 4(b), 4(d), 4(f) and 4(h), we observe that LLO, which has a
low altitude of $100$ $\mathrm{k}\mathrm{m}$, exhibits the smallest time
spacing between consecutive regions of ECOP.

Table 3: Comparison analysis across different orbit types, where NRHO exhibits the least maximum ECOP as well as the greatest visibility of both at least one and at least four Earth-GPS satellites. Orbit Type |  Max ECOP ($\mathrm{s}$) | Satellite Visibility (%)
---|---|---
$\geq 1$ | $\geq 4$
ELFO | 3360 | 99.2 | 92.1
NRHO | 420 | 99.9 | 92.9
LLO | 2880 | 61.5 | 55.6
PCO | 3900 | 98.1 | 90.8

(a) ELFO: Full satellite visibility

(b) ELFO: Zoomed-in satellite visibility

(c) NRHO: Full satellite visibility

(d) NRHO: Zoomed-in satellite visibility

(e) LLO: Full satellite visibility

(f) LLO: Zoomed-in satellite visibility

(g) PCO: Full satellite visibility

(h) PCO: Zoomed-in satellite visibility

Figure 4: Earth-GPS satellite visibility and maximum ECOP across different
orbit types, where blue dotted lines indicate the number of visible Earth-GPS
satellites and the red vertical bars indicate regions of ECOP; a), c), e) and
g) show the satellite visibility for the entire time duration in ELFO, NRHO,
LLO and PCO, respectively; b), d), f) and h) show the zoomed-in satellite
visibility for a shorter time segment of $2$ days. We observe that NRHO not
only exhibits a least maximum ECOP of $420$ $\mathrm{s}$ but also the greatest
at least one satellite visibility of $99.9\%$.

###  4.3 Across Clock Types: Variation in RMS Timing Errors

Table 4 compares the RMS estimation error in clock bias and drift across
different case studies. Intuitively, RMS timing error provides insights
regarding the component of lunar UERE metric, whose results will be discussed
in the next two subsections.

Table 4: Comparison analysis across different orbit types, where NRHO exhibits the least maximum ECOP as well as the greatest visibility of both at least one and at least four Earth-GPS satellites. | RMS Timing Error
---
in Clock Type
ELFO | NRHO | LLO | PCO
Microchip CSAC | Bias ($\mu\mathrm{s}$) | 0.0405 | 0.0383 | 0.0679 | 0.0408
Drift ($\mathrm{n}\mathrm{s}\mathrm{/}\mathrm{s}$) | 0.0407 | 0.0401 | 0.0419 | 0.0404
Microchip MAC | Bias ($\mu\mathrm{s}$) | 0.0213 | 0.0215 | 0.0278 | 0.0217
Drift ($\mathrm{n}\mathrm{s}\mathrm{/}\mathrm{s}$) | 0.0081 | 0.0081 | 0.0081 | 0.0083
SRS PRS-10 | Bias ($\mu\mathrm{s}$) | 0.0176 | 0.0177 | 0.0213 | 0.0178
Drift ($\mathrm{n}\mathrm{s}\mathrm{/}\mathrm{s}$) | 0.0045 | 0.0046 | 0.0045 | 0.0046
Excelitas RAFS | Bias ($\mu\mathrm{s}$) | 0.0123 | 0.0121 | 0.0143 | 0.0121
Drift ($\mathrm{n}\mathrm{s}\mathrm{/}\mathrm{s}$) | 0.0015 | 0.0015 | 0.0015 | 0.0015
NASA’s DSAC | Bias ($\mu\mathrm{s}$) | 0.0044 | 0.0044 | 0.0047 | 0.0037
Drift ($\mathrm{n}\mathrm{s}\mathrm{/}\mathrm{s}$) | 6.4347$\times 10^{-5}$ | 6.5224$\times 10^{-5}$ | 6.2655$\times 10^{-5}$ | 5.6112$\times 10^{-5}$

(a) Clock bias error for entire experiment duration

(b) Clock bias error for a smaller time segment

(c) Clock drift error for entire duration

(d) Clock drift error for a smaller time segment

Figure 5: Comparison of estimation errors in clock bias and drift across
different orbit types with an onboard Microchip CSAC, where ELFO is indicated
in orange, NRHO in green, LLO in magenta, and PCO in blue; a) and c)
demonstrate the errors in clock bias and clock drift for the entire experiment
duration, respectively, while b) and d) demonstrate the zoomed-in errors in
clock bias and clock drift for a smaller time segment of $2$ days. We observe
that the three orbit types, namely ELFO, NRHO and PCO, demonstrate a
comparable RMS timing error of $<0.0408$ $\mu\mathrm{s}$ in clock bias and
$<0.0407$ $\mathrm{n}\mathrm{s}\mathrm{/}\mathrm{s}$ in clock drift, while the
LLO exhibits slightly higher RMS error in clock bias of $0.0679$
$\mu\mathrm{s}$.

For a given Microchip CSAC $(17$ $\mathrm{c}\mathrm{m}^{3}$$\cdot 0.0035$
$\mathrm{k}\mathrm{g}$$\cdot 0.12$ $\mathrm{W}$$)$, whose low SWaP
characteristics were reported earlier in Table 2, Fig. 5 provides an
illustration of the variation in timing errors across the four orbit types,
namely ELFO, NRHO, LLO and PCO. In Fig. 5, we follow the same color coding as
that in Figs. 2 and 3 to denote different orbit types, with ELFO indicated in
orange, NRHO in green, LLO in magenta and PCO in blue. As described previously
in the description of validation metrics, we reduce the measurement update
rate by setting the sampling period to $T_{\text{meas}}=mT_{\text{pred}}$,
where $m=5$. We observe that the three orbit types, namely ELFO, NRHO and PCO,
demonstrate a comparable RMS timing error of $<0.0408$ $\mu\mathrm{s}$ in
clock bias and $<0.0407$ $\mathrm{n}\mathrm{s}\mathrm{/}\mathrm{s}$ in clock
drift, while the case study based on LLO shows a higher RMS error in clock
bias of $0.0679$ $\mu\mathrm{s}$. This observation implies that the RMS timing
error not only depends on the least maximum ECOP and the greatest satellite
visibility but also on orbital parameters, namely eccentricity, inclination
and altitude, that govern the geometric configuration between the Earth-GPS
and LNSS. From Table 4, we also observe that, as the SWaP and the timing
stability of the onboard clock increases, the variation in RMS error across
orbit types becomes less significant, wherein the RMS error in clock bias and
drift for DSAC are $<0.0044$ $\mu\mathrm{s}$ and $<6.5224\times 10^{-5}$
$\mathrm{n}\mathrm{s}\mathrm{/}\mathrm{s}$, respectively.

###  4.4 Across Case Studies and Earth-GPS Measurement Update Rates:
Sensitivity Analysis of Lunar UERE

Fig. 6(a) demonstrates the variation in lunar UERE metric across case studies
conducted while considering no reduction in the measurement update rate, i.e.,
we set $T_{\text{meas}}=T_{\text{pred}}=60$ $\mathrm{s}$ or $m=1$. We validate
that our proposed time-transfer achieves a low lunar UERE of $<10$
$\mathrm{m}$ for most case studies except the one involving the LLO and the
Microchip CSAC, which exhibits a value of $18.4$ $\mathrm{m}$. These
observations imply that, if the desired lunar UERE to be maintained by an LNSS
satellite is, for instance, $<10$ $\mathrm{m}$ at all times, we can easily opt
for an onboard clock that falls in the lower end of the SWaP spectrum, such as
Microchip CSAC $(17$ $\mathrm{c}\mathrm{m}^{3}$$\cdot 0.0035$
$\mathrm{k}\mathrm{g}$$\cdot 0.12$ $\mathrm{W}$$)$ or Microchip MAC $(50$
$\mathrm{c}\mathrm{m}^{3}$$\cdot 0.0084$ $\mathrm{k}\mathrm{g}$$\cdot 0.5$
$\mathrm{W}$$)$, instead of a high-SWaP clock, such as Excelitas RAFS or
NASA’s DSAC. Given a desired lunar UERE to maintain, we can also wisely choose
the orbit type that is easy to maintain and has a longer lifespan, such as
LLO, PCO or ELFO, over the more complex NRHO that requires constant station-
keeping maneuvers to maintain stability. Furthermore, for future
investigations, we observe the potential for a heterogenous, SmallSat-based
lunar PNT constellation, wherein we can choose the grade of the onboard clock
based on the orbit of each LNSS satellite to satisfy a desired lunar UERE. For
instance, to maintain a desired lunar UERE of $<5$ $\mathrm{m}$, we observe
from Fig. 6(a) that we can design a heterogenous LNSS constellation based in
ELFO and LLO, wherein the satellites in LLO can be equipped with a PRS-10
clock while the ones in ELFO can be equipped with a lower-SWaP Microchip MAC.

(a) $m=1$ or $T_{\text{meas}}=1$ $\mathrm{m}\mathrm{i}\mathrm{n}$

(b) $m=5$ or $T_{\text{meas}}=5$ $\mathrm{m}\mathrm{i}\mathrm{n}$

(c) $m=30$ or $T_{\text{meas}}=30$ $\mathrm{m}\mathrm{i}\mathrm{n}$

(d) $m=60$ or $T_{\text{meas}}=60$ $\mathrm{m}\mathrm{i}\mathrm{n}$

Figure 6: Sensitivity analysis of lunar UERE metric across different case
studies for reduced measurement update rates where (a) $m=1$ is the baseline
case with no reduction of the measurement update rate, i.e.,
$T_{\text{meas}}=T_{\text{pred}}$; (b) $m=5$; (c) $m=30$; and (d) $m=60$.

Figures 6(b)-6(d) show the variation in lunar UERE metric for three cases of
reduced update rates, where $T_{\text{meas}}=mT_{\text{pred}}$ with $m=5$,
$m=30$ and $m=60$ respectively. As the Earth-GPS measurement update rate
increases for low-SWaP clocks with sampling periods from $T_{\text{meas}}=60$
$\mathrm{m}\mathrm{i}\mathrm{n}$ ($m=60$) to $T_{\text{meas}}=1$
$\mathrm{m}\mathrm{i}\mathrm{n}$ ($m=1$), we observe an increased sensitivity
of the lunar UERE metric across lunar orbit types, i.e., difference in value
between LLO and others increases. We also demonstrate that the estimated lunar
UERE is $<30$ $\mathrm{m}$ for a reduced Earth-GPS measurement update rate
with a sampling period of up to $T_{\text{meas}}=30$
$\mathrm{m}\mathrm{i}\mathrm{n}$, which is comparable in order of magnitude to
that of the baseline case with $m=1$ (in Fig. 6(a)) as well as the legacy
Earth-GPS. Thus, we validate that our LNSS design that utilizes time-transfer
from Earth-GPS at a reduced measurement update rate to lower the SWaP
requirements of the onboard clock. Furthermore, we can re-observe the
potential for heterogeneous, SmallSat-based PNT constellation explained above,
wherein we can intelligently choose not only the grade of the onboard clock
but also the reduced update rate of the timing filter based on the LNSS
satellite orbit so as to satisfy a desired lunar UERE.

##  5 Conclusion

We performed an exhaustive case study analysis for designing a SmallSat-based
Lunar Navigation Satellite System (LNSS) with time-transfer from Earth-GPS,
wherein we investigate the trade-off between different design considerations
related to the onboard clock and the lunar orbit type. In our time-transfer,
we alleviated the Size, Weight and Power (SWaP) requirements of the onboard
clocks by leveraging the intermittently available Earth-GPS signals to provide
timing corrections, and also designed a lunar User Equivalent Range Error
(UERE) metric to characterize the ranging accuracy of the LNSS satellite.

Using high-fidelity simulations of an LNSS satellite in the System Tool Kit
(STK) software by Analytical Graphics, Inc. (AGI), we designed multiple case
studies comprising five onboard clocks and four lunar orbit types. We
validated that the least maximum Earth-GPS Continual Outage Period (ECOP) of
only $420$ $\mathrm{s}$ is observed for a Near-Rectilinear Halo Orbit (NRHO)
since it experiences fewer occultations from Earth and the Moon given its high
altitude. We demonstrated that a low lunar UERE of $<30$ $\mathrm{m}$ can be
achieved using low-SWaP onboard clocks (e.g., Microchip chip scale atomic
clock, Stanford Research Systems PRS-10) even for reduced Earth-GPS
measurement update rates with sampling period of up to $30$
$\mathrm{m}\mathrm{i}\mathrm{n}$. Through a case study analysis of our time-
transfer from Earth-GPS, we demonstrated that lower-SWaP onboard clocks and
easier-to-maintain lunar orbit types can be chosen over the others to still
achieve a desired lunar UERE across the entire LNSS constellation.

## ACKNOWLEDGEMENTS

We would like to thank the Analytical Graphics, Inc. (AGI) Educational
Alliance Program (EAP) for providing us with the System Tool Kit (STK)
software license to perform this research work. We would also like to thank
Keidai Ilyama for the insightful discussions related to this work and Daniel
Neamati for reviewing this paper.

## References

  * [1] Smith, M., Craig, D., Herrmann, N., Mahoney, E., Krezel, J., McIntyre, N., and Goodliff, K., “The Artemis program: An overview of NASA’s activities to return humans to the Moon,” in _2020 IEEE Aerospace Conference_. IEEE, 2020, pp. 1–10.
  * [2] Laurini, K. C. and Gerstenmaier, W. H., “The global exploration roadmap and its significance for NASA,” _Space Policy_ , vol. 30, no. 3, pp. 149–155, 2014.
  * [3] Israel, D. J., Mauldin, K. D., Roberts, C. J., Mitchell, J. W., Pulkkinen, A. A., La Vida, D. C., Johnson, M. A., Christe, S. D., and Gramling, C. J., “LunaNet: A flexible and extensible lunar exploration communications and navigation infrastructure,” in _2020 IEEE Aerospace Conference_. IEEE, 2020, pp. 1–14.
  * [4] Cozzens, T., “Galileo will help lunar pathfinder navigate around Moon,” _GPS World_ , 2021. [Online]. Available: https://www.gpsworld.com/galileo-will-help-lunar-pathfinder-navigate-around-moon/?utm_source=Navigate!+Weekly+GNSS+News&utm_medium=Newsletter&utm_campaign=NCMCD210331002&oly_enc_id=6244C4859912F3Z
  * [5] Folta, D. and Quinn, D., “Lunar frozen orbits,” in _AIAA/AAS Astrodynamics Specialist Conference and Exhibit_ , 2006, p. 6749.
  * [6] Whitley, R. and Martinez, R., “Options for staging orbits in cislunar space,” in _2016 IEEE Aerospace Conference_. IEEE, 2016, pp. 1–9.
  * [7] Schonfeldt, M., Grenier, A., Delépaut, A., Giordano, P., Swinden, R., Ventura-Traveset, J., Blonski, D., and Hahn, J., “A system study about a lunar navigation satellite transmitter system,” in _2020 European Navigation Conference (ENC)_. IEEE, 2020, pp. 1–10.
  * [8] Schmittberger, B. L. and Scherer, D. R., “A review of contemporary atomic frequency standards,” _arXiv preprint arXiv:2004.09987_ , 2020.
  * [9] Ely, T. A., Burt, E. A., Prestage, J. D., Seubert, J. M., and Tjoelker, R. L., “Using the deep space atomic clock for navigation and science,” _IEEE transactions on ultrasonics, ferroelectrics, and frequency control_ , vol. 65, no. 6, pp. 950–961, 2018.
  * [10] Melosh, H., Freed, A. M., Johnson, B. C., Blair, D. M., Andrews-Hanna, J. C., Neumann, G. A., Phillips, R. J., Smith, D. E., Solomon, S. C., Wieczorek, M. A. _et al._ , “The origin of lunar mascon basins,” _Science_ , vol. 340, no. 6140, pp. 1552–1555, 2013.
  * [11] Winternitz, L. B., Bamford, W. A., Price, S. R., Carpenter, J. R., Long, A. C., and Farahmand, M., “Global positioning system navigation above 76,000 km for NASA’s magnetospheric multiscale mission,” _NAVIGATION, Journal of the Institute of Navigation_ , vol. 64, no. 2, pp. 289–300, 2017.
  * [12] Johnson-Groh, M., “NASA’s MMS breaks Guinness world record,” _NASA_ , Nov 2016. [Online]. Available: https://www.nasa.gov/feature/goddard/2016/nasa-s-mms-breaks-guinness-world-record
  * [13] Baird, D., “Record-breaking satellite advances NASA’s exploration of high-altitude GPS,” _NASA_ , Apr 2019. [Online]. Available: https://www.nasa.gov/feature/goddard/2019/record-breaking-satellite-advances-nasa-s-exploration-of-high-altitude-gps
  * [14] Winternitz, L., Bamford, B., Long, A., and Hassouneh, M., “GPS based autonomous navigation study for the Lunar Gateway,” 2019.
  * [15] Cheung, K.-M., Lee, C., and Heckman, D., “Feasibility of ‘weak GPS’ real-time positioning and timing at lunar distance,” in _2020 IEEE Aerospace Conference_. IEEE, 2020, pp. 1–7.
  * [16] Donaldson, J. E., Parker, J. J., Moreau, M. C., Highsmith, D. E., and Martzen, P. D., “Characterization of on-orbit GPS transmit antenna patterns for space users,” _NAVIGATION, Journal of the Institute of Navigation_ , vol. 67, no. 2, pp. 411–438, 2020.
  * [17] Winternitz, L., Moreau, M., Sirotzky, S. _et al._ , “Navigator GPS receiver for fast acquisition and weak signal space applications,” in _Proceedings of the 17th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS 2004)_ , 2004, pp. 1013–1026.
  * [18] “ESA, NASA race to the Moon for first lunar GNSS fix,” _Inside GNSS_ , Mar 2021. [Online]. Available: https://insidegnss.com/esa-nasa-race-to-the-moon-for-first-lunar-gnss-fix/
  * [19] Kraft, R., “NASA enlists commercial partners to fly payloads to Moon,” _NASA_ , Sept 2020. [Online]. Available: https://www.nasa.gov/feature/nasa-enlists-commercial-partners-to-fly-payloads-to-moon
  * [20] “Galileo will help Lunar Pathfinder navigate around Moon,” Mar 2021. [Online]. Available: https://www.esa.int/Applications/Navigation/Galileo_will_help_Lunar_Pathfinder_navigate_around_Moon
  * [21] Bhamidipati, S., Mina, T., and Gao, G., “Design considerations of a Lunar Navigation Satellite System with time-transfer from Earth-GPS,” in _Proceedings of the 34th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS+ 2021)_ , 2021, pp. 950–965.
  * [22] Krawinkel, T. and Schön, S., “Benefits of receiver clock modeling in code-based GNSS navigation,” _GPS Solutions_ , vol. 20, no. 4, pp. 687–701, 2016.
  * [23] Martoccia, D., Bernstein, H., Chan, Y., Frueholz, R., and Wu, A., “GPS satellite timing performance using the autonomous navigation (Autonav),” in _Proceedings of the 11th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GPS 1998)_ , 1998, pp. 1705–1712.
  * [24] Van Dierendonck, A. J. and McGraw, J., “Relationship between Allan variances and Kalman filter parameters,” in _Proceedings of the 16th Annual Precise Time and Time Interval Systems and Applications Meeting_ , 1984, pp. 273–293.
  * [25] Kaplan, E. D. and Hegarty, C., _Understanding GPS/GNSS: Principles and applications_. Artech House, 2017.
  * [26] Capuano, V., Botteron, C., Leclère, J., Tian, J., Wang, Y., and Farine, P.-A., “Feasibility study of GNSS as navigation system to reach the Moon,” _Acta Astronautica_ , vol. 116, pp. 186–201, 2015.
  * [27] Capuano, V., Basile, F., Botteron, C., and Farine, P.-A., “GNSS-based orbital filter for Earth Moon Transfer Orbits,” _The Journal of Navigation_ , vol. 69, no. 4, pp. 745–764, 2016.
  * [28] Borio, D., Sokolova, N., and Lachapelle, G., “Doppler measurement accuracy in standard and high-sensitivity global navigation satellite system receivers,” _IET Radar, Sonar & Navigation_, vol. 5, no. 6, pp. 657–665, 2011.
  * [29] “Systems Tool Kit (STK),” _AGI_ , 2021. [Online]. Available: https://www.agi.com/products/stk
  * [30] Delépaut, A., Giordano, P., Ventura-Traveset, J., Blonski, D., Schönfeldt, M., Schoonejans, P., Aziz, S., and Walker, R., “Use of GNSS for lunar missions and plans for lunar in-orbit development,” _Advances in Space Research_ , vol. 66, no. 12, pp. 2739–2756, 2020.
  * [31] Schönfeldt, M., Grenier, A., Delépaut, A., Giordano, P., Swinden, R., and Ventura-Traveset, J., “Across the lunar landscape: Towards a dedicated lunar PNT system,” _Inside GNSS_ , 2020.
  * [32] Chao, C., Cook, J. W., Cox, J., Starchville, T. F., Thompson, R. C., and Wagner, L., “Independent verification and validation for Analytical Graphics, Inc. of three astrodynamic functions of the Satellite Tool Kit: Version 4.1. 0,” _The Aerospace Corporation ATR-2000 (7605)-1_ , 2000.
  * [33] Montenbruck, O. and Gill, E., _Satellite orbits: models, methods and applications_. Springer Science & Business Media, 2012.
  * [34] Ely, T. A. and Lieb, E., “Constellations of elliptical inclined lunar orbits providing polar and global coverage,” _Journal of the Astronautical Sciences_ , vol. 54, no. 1, pp. 53–67, 2006.
  * [35] Delépaut, A., Schönfeldt, M., Giordano, P., Blonski, D., Sarnadas, R., Ries, L., and Ventura-Traveset, J., “A system study for cislunar radio navigation leveraging the use of realistic galileo and gps signals,” in _Proceedings of the 32nd International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS+ 2019)_ , 2019, pp. 1199–1219.
  * [36] Williams, J., Lee, D. E., Whitley, R. J., Bokelmann, K. A., Davis, D. C., and Berry, C. F., “Targeting cislunar near rectilinear halo orbits for human space exploration,” 2017.