---
author:
- Max Mahlke
- Enrique Solano
- Hervé Bouy
- Benoit Carry
- 'Gijs A. Verdoes Kleijn'
- Emmanuel Bertin
bibliography:
- 'bib.bib'
title: 'The ssos Pipeline: Identification of Solar System Objects in Astronomical Images'
---

Introduction
============

The identification and characterisation of Solar System objects (SSOs) are key steps in the understanding of the formation of the Solar System. Early evolutionary processes like planetary migration left their imprints in the distributions of the orbital parameters, chemical compositions, and size distributions of the bodies of the Solar System. The dynamical properties of the different populations have given rise to numerous models of planetary formation and evolution [@tsiganis2005origin; @morbidelli2005chaotic; @gomes2005origin; @walsh2011low; @2017SciA....3E1138R]. Multi-wavelength data is required to determine the spectral properties of the SSO populations and thereby constrain and refine these scenarios [@demeo2014solar].

Besides the science case, there are more practical implications of SSO detection and identification. Near-Earth Objects (NEOs) pose a threat to Earth, as reminded by impacts such as over the Russian city of Chelyabinsk in 2013 [@brown2013500]. The cited event raised awareness for the problem that sub-kilometre NEOs can have devastating local impacts. While dedicated efforts like the Catalina Sky Survey likely completely surveyed the population of NEOs larger than 1km in diameter, smaller NEOs remain a constant risk [@1998BAAS...30.1037L; @2008Natur.453.1178H]. Early detection and a detailed compositional study is required to choose among different response strategies.

To derive model constraints and assess risks, observational data is necessary. Since the discovery of the first asteroid (1) Ceres in 1801 at the Palermo Observatory, numerous efforts haven been taken to survey the SSO populations (see @2015aste.book..795J for a review). Over 790,000 objects have been detected and their orbital parameters were determined, though 9% carry ephemeris uncertainties equal to or larger than 1deg in January 2019, underlining the need for further observations.[^1] The physical and compositional characterisation of these objects is lacking far behind, due to the large increase of both required effort and possible parameter space between the identification and physical characterisation of SSOs.

Meanwhile, astronomical observatories are producing more imaging data each night than can be analysed in real-time. The Large Synoptic Survey Telescope (LSST), the next-generation wide-field telescope currently being built in Chile, will produce on average 15TB of raw data each observation night [@juric2015lsst]. Serendipitous observations of SSOs are a well-known by-product in these images and their astrometry and photometry is sometimes recovered. Automatic SSO detection pipelines have been developed for surveys where the required time and funding are available.

The Pan-STARRS survey demonstrates the full use of SSO discovery potential through pipelines. Their Moving Object Processing System (MOPS) automatically identifies SSOs in the survey images and reports them to the Minor Planet Centre (MPC) [@panstarrs2013]. In October 2017, Pan-STARRS observed the first interstellar object, 1I/2017 U1, named ‘Oumuamua [@meech2017]. The photometric pipeline of the Sloan Digital Sky Survey (SDSS) automatically flags moving objects based on a $\chi^2$-fit of their coordinates if the object has been observed at least three times [@2000AJ....120.1579Y; @2001ASPC..238..269L]. The fourth data release of the Moving Object Catalogue (MOC) contains 471,569 moving objects, with an estimated completeness of 95% and a contamination of $\sim$6% [@ivezic2002asteroids]. Many of these moving objects could only later be associated to asteroids with certainty thanks to the increased number of SSO observations, highlighting the potential of legacy archive science [@CARRY2016340]. A citizen-science projected led by @solano2014precovery recovered 2,351 new measurements of 551 NEOs in the images of SDSS. 73% of these measurements were recovered from images taken before May 2007, meaning the objects were not identified by the automated SDSS pipeline.

In space, the NASA *WISE* mission [@wise] observed about 150,000 small bodies while surveying the entire celestial sphere in the mid-infrared, and reported their albedo [@masiero2011]. The ESA *Gaia* mission [@gaia], repeatedly surveying the sky down to magnitude 20 in the visible, has already reported over 2,000,000 observations of 14,099 small bodies [@spoto2018]. The final catalog will contain multi-epoch astrometry, visible photometry and spectra for about 300,000 small bodies [@2007-EMP-101-Mignard]. The ESA *Euclid* space survey telescope [@2011arXiv1110.3193L] with an expected launch in 2022 will provide both infrared and optical observations of around 150,000 SSOs, allowing for detailed compositional studies [@Carry2018]. An automated identification of serendipitous SSO observations is investigated, e.g. by means of a convolutional neural network [@2018arXiv180710912L].

Yet, for smaller missions or PI observations on large facilities, the resources for the development of an automated detection pipeline may not be available. A template solution may lower this threshold sufficiently to make the identification of SSOs feasible. As the telescope and observation conditions differ vastly from observatory to observatory, a mission-agnostic pipeline can only have one basic requirement: at least three images need to cover an overlapping area of the sky within a reasonable amount of time to identify moving sources. This work introduces the [`ssos`]{} pipeline for the identification of SSOs and the recovery of their astrometric and photometric properties from astronomical images. As versatile tool, it is intended to be a template pipeline for all kinds of observatories, with minimal set-up and hardware requirements.

The [`ssos`]{} pipeline was developed from the Kilo-Degree Survey (KiDS) SSO identification pipeline described in @Mahlke2017. This proof-of-concept, survey-specific software recovered more than 20,000 SSO candidates with a degree of contamination below 0.05%. However, it was tailored specifically to the KiDS observation strategy, instrument characteristics, and data product file specifications. By applying the alpha-version to images from other observatories, these dependencies were identified and removed or weakened. Among other changes, the [`ssos`]{} pipeline no longer relies on regular image cadence or specific FITS keywords in the header. Common observation cases like SSOs blended with other sources, missed SSO detections, long temporal baselines of observations, and single- or multi-band observations are treated better, leading to more SSOs being identified. The [`python`]{} implementation runs faster, has improved error-handling, and the distribution via the Python Package Index[^2] offers a simple installation process.

In , the [`ssos`]{} pipeline is described from a conceptual point of view. The pipeline implementation and computational performance are discussed in . In , the application and results of two test studies are presented, highlighting the versatility of the [`ssos`]{} pipeline. Conclusions and an outlook on future development goals are given in .

(data) [Images, (Weights, Masks)]{}; (sextractor) \[below =of data\] [*SExtractor*]{}; (sextractor\_settings) \[right =of sextractor\] [SExtractor Configuration]{}; (scamp) \[below =of sextractor, yshift=-0.5cm\] [*SCAMP*]{}; (ref\_cat) \[left =of scamp\] [Reference Catalogue\
from VizieR Service]{}; (scamp\_settings) \[right =of scamp\] [SCAMP Configuration]{};

(filter) \[below =of scamp, yshift=-1cm\] [[`python`]{} *Filter Routine*]{}; (filter\_settings) \[right =of filter, yshift=0.4cm\] [[`ssos`]{} Configuration]{}; (hyg) \[below =of filter\_settings, yshift=1.0cm\] [Bright-Sources Catalogue]{};

(candidates) \[below =of filter, yshift=-1\] [SSO Candidates]{}; (cutouts) \[right =of candidates, yshift=-1\] [Cutouts]{}; (skybot) \[left =of candidates\] [Known-SSOs Catalogue\
from IMCCE]{};

(data) edge \[bend left =0\] (sextractor); (sextractor\_settings) edge \[bend right =0\](sextractor); (sextractor) edge \[bend left =0\] node\[right = .3 cm\] [Source Catalogues]{}(scamp);

(ref\_cat) edge \[bend left =0\] (scamp); (scamp\_settings) edge \[bend right =0\](scamp); (scamp) edge \[bend left =0\] node\[right = .3 cm, align=left\] [Source Catalogues\
Full and Merged Catalogue]{}(filter);

(filter) edge \[bend left =0\](candidates); (candidates) edge \[bend left =0\](cutouts); (skybot) edge \[bend left =0\](candidates);

(filter\_settings) – (filter)coordinate\[pos=0.5\](mm); (filter\_settings) –(filter\_settings.west)-| (mm) |- (filter.east); (hyg) – (filter)coordinate\[pos=0.5\](mm); (hyg) –(hyg.west)-| (mm) |- (filter.east);

(c4) ; (data.east) to node \[xshift=2em\] (c4); (data.east) – (c4) |- (cutouts.east);

Pipeline design {#sec:pipeline}
===============

The focus of the pipeline design is placed on versatility. As a template solution for distinct telescope and observation set-ups, the requirements for the imaging data should be minimal. This versatility is achieved by relying on the widely used astr*O*matic[^3] software packages SExtractor and SCAMP for source detection and association . Both routines are robust, fast, and can handle many kinds of imaging data thanks to a large degree of configurability. After source detection with SExtractor, SCAMP is used to compute the global astrometric solution of the multi-epoch images. The resulting catalogues contain all stars, galaxies, imaging artefacts, and SSOs present in the images. A set of filter algorithms written in [`python`]{} then identifies SSOs in this population primarily based on their apparent motion and trail appearance. The use of [`python`]{} allows for simple automation, extensive error handling, and informative output to guide the user. The [`ssos`]{} pipeline itself serves as a wrapper for these three steps and is written in [`python`]{}.

All SSO selection filters are optional and configurable. Their parameters allow for different degrees of strictness in SSO candidate rejection. In general, the parameter selection steers the final set of SSO candidates either towards purity, i.e. no false-positive identifications, or completeness, i.e. all SSOs in the images are recovered. While in theory not mutually exclusive properties, in practice getting a complete and pure set of candidates is hardly feasible, and the user has to place the emphasis considering e.g. the final size of the output sample and the possibility of visual confirmation of the nature of all candidates. On the other hand, experience shows that artefacts are readily identified as outliers e.g. in colour-colour diagrams in the subsequent analysis of the SSO data.

This section steps through the pipeline processes and explains them on a physical, non-computational level. An overview of all steps and dependencies in the [`ssos`]{} pipeline is given in . Before executing the pipeline, the telescope metadata such as pixel scale, detector gain, and the FITS file properties such as the designation of the keywords of right ascension and declination coordinates in the headers must be entered in the configuration files of SExtractor, SCAMP, and the [`ssos`]{} pipeline itself.[^4]

Source detection: SExtractor
----------------------------

SExtractor offers source detection in large astronomical images requiring only marginal setup. The source detection is described here briefly, the reader is referred to the official documentation[^5] and for details.

In agreement with the design principle of the [`ssos`]{} pipeline, the minimal input for source detection with SExtractor are the FITS images. Corresponding weight images and pixel masks may be supplied, but the routine offers background estimation and several ways to reject artificial sources introduced by bad pixel apart from these auxiliary files.

A low-resolution, gridded map of the sky background is first computed by SExtractor by estimating the mode of the histogram of pixel values in each grid mesh. A full resolution background model is then derived through bi-cubic spline interpolation of this low-resolution map. After subtracting the background model, sources are detected by thresholding. By default, a source requires at least 5 contiguous pixel with more than 1.5 standard deviations of the local background to be extracted. Neighbouring sources that have merged in the image are deblended by re-thresholding the flux distribution in the source pixels in exponentially-spaced levels and evaluating the area and intensity of peaks at each level. Photometric measurements and morphometric parameters are extracted by calculating an adaptive elliptical aperture from the second order moments of the source’s profile.

Each step of the source detection can be adapted by the user via a configuration file. For a reliable detection of SSOs, special focus should be placed on the source deblending parameters. As SSOs move through the field-of-view (FoV), they may get close in angular distance to other sources. If these mergers are not correctly deblended, the SSO is not detected in the image and the subsequent filter routine may reject all other detections of this SSO on the grounds of non-linear motion.

The source detection by SExtractor effectively converts the imaging pixel data to catalogue data, creating one catalogue of detected sources per input image. The next step aims to associate the detections of individual sources at distinct epochs by cross-matching these catalogues.

Source association: SCAMP {#sec:scamp}
-------------------------

Standard procedure in astronomical observations suggests a dithering pattern when observing a part of the sky to bridge gaps in the camera CCDs and avoid biases e.g. introduced by hot pixel. The FoVs of the images are therefore shifted and possibly rotated and scaled with respect to each other. An astrometric calibration is required to trace source observations from one image to another, which is computed with SCAMP. Again, this computation is described in brevity, see @2006ASPC..351..112B and @Bouy2013 for details.

SCAMP uses the centroid source positions retrieved by SExtractor to compute the astrometric solution. A reference catalogue is retrieved from the Vizier database and detected and reference sources are cross-correlated for each input image. Relative shifts, scaling, and rotations between the FoVs are thereby corrected for and the reprojected, detected sources are cross-matched in overlapping images. Finally, the astrometric solution is derived by a $\chi^2$-minimisation of the quadratic sum of differences in the positions of overlapping detected sources in pairs of images.

SCAMP computes the proper motions of all sources by computing a linear, weighted $\chi^2$-fit of the source coordinates over observation epoch. The reported uncertainty on the proper motions is computed from the covariance matrix of the best-fit parameters. Sources which display non-linear motion therefore have larger errors on their proper motion. Outlier detections in right ascension or declination space are identified by means of $\chi^2$-reduction. The source observation which increases the $\chi^2$ of the fit the most is seen as outlier and removed from the fit. This procedure is repeated until the reduced $\chi^2$ is below 6 or 20% of the source detections were rejected.

As SExtractor, SCAMP offers a large degree of configuration. For reliable SSO detection, proper setting of the [`CROSSID_RADIUS`]{} is paramount. It defines the maximum distance in arcseconds between two source detections that SCAMP will consider for cross-matching.[^6] Therefore, the cross-match radius sets an intrinsic upper limit on the proper motion of detected sources. For each pair of subsequent exposures, source detections which are displaced by more than this radius will not be regarded as belonging to the same source. Setting the [`CROSSID_RADIUS`]{} paramter too large, however, leads to an increase in mis-matched sources. A brief investigation of the optimal value is presented in .

SSO identification: [`python`]{} filter routine
-----------------------------------------------

The difficult tasks of source detection and association across observation epochs are handled by the widely-used SExtractor and SCAMP packages, with the [`ssos`]{} pipeline serving as a wrapper to handle the catalogue and image dependencies. The following filter routine is novel and has been developed specifically for the identification of SSO in source catalogues recovered from a series of images [@Mahlke2017].

The main characteristic of SSOs to differentiate them from other sources, astrophysical or artificial, is their apparent motion across the sky. Even over short observation times, they describe in first-order linear trails in images. This feature particularly sets them apart from stars and galaxies. Meanwhile, artefacts like hot pixel or diffraction spikes can mimic the linear apparent motion if a linear dithering pattern is used. The differentiating feature here is the consistent appearance of the SSO trail in images taken under similar observation conditions. The following six filter algorithms are based on these premises.

#### Number of Detections

The number of detections for each source is between 1 and *N*, where *N* is the total number of input images. Typically, stars and galaxies have close to $N$ detections as they do not leave the FoV unless the dithering pattern excludes them. Imaging artefacts like cosmic rays (CRs) make up the majority of sources with 1-3 detections as they only appear in singular images and may be associated by chance to another artefact by SCAMP. An example distribution of the number of detections per source in a series of 12 images is shown in . The peak of sources with 5 detections is specific for this survey set-up (J-PLUS, see ), which contains 5 broad-band in the set of 12 filters, increasing the chance of observing a source five times.

The number of detections for SSOs is largely dependent on the exposure time and observed bands of the survey. Filtering sources with few detections ($\leq$ 3) is still advisable to remove a large fraction of imaging artefacts. Furthermore, the ability to correctly evaluate a source’s motion and appearance over time increases with each detection, strongly reducing the probability of a false-positive detection.

#### Bad Pixel

Hot pixels in CCDs mimic linearly moving sources of constant appearance in case of a linear dithering pattern. While there are several ways to eliminate these artificial sources from the sample prior to this step,[^7] a simple yet effective method is provided at this point.

Bad-pixel sources are characterised by their fixed position in the CCD-space. The centroid position recovered by SExtractor lies on or adjacent to the bad pixel itself, while non-artificial sources should not fall onto the same pixel in every exposure. The [`DELTA_PIXEL`]{} filter parameter defines the minimum spread in pixel-space that sources have to exhibit to pass this filter. shows a rejected source, a hot pixel, in the upper left corner, while an accepted SSO projected onto CCD-space is shown below.

#### Proper Motion

Motion is the primary characteristic differentiating an SSO from extrasolar sources. A strong separation of non-SSOs and SSOs is the signal-to-noise ratio (SNR) of the source candidate’s proper motion. As explained in Sec. \[sec:scamp\], SCAMP computes the proper motion with a weighted, linear $\chi^2$-fit. A large uncertainty derived from the best fit parameters’ covariance matrix indicates non-linear motion. A lower limit on the SNR therefore rejects sources exhibiting non-linear proper motion. The SNR is defined as

$$\mathrm{SNR}_\mu = \mu / \sigma_\mu,$$

where $\mu$ is the absolute value of the proper motion and $\sigma_\mu$ the combined error of the proper motion components as derived by SCAMP.

shows that setting the [`PM_SNR`]{} filter parameter is equivalent to introducing a lower limit on the proper motion. In general, sources with larger proper motion have a larger SNR. Extrasolar sources, however, will exhibit random coordinate fluctuations on the order of the seeing, leading to non-linear apparent motion and hence a small SNR in the proper motion calculation.

The absolute value of the apparent motion should be restricted with care. As shown e.g. in @Carry2018, the proper motion space of SSOs encompasses several orders of magnitudes, from a few hundredths of [$^{\prime\prime}$]{}/h for distant Kuiper-Belt objects (KBOs) to over 1000[$^{\prime\prime}$]{}/h for the close population of NEOs. Instead, the possibility to filter the proper motion values outside a user-specified range is intended for focusing either on fast- or slow-moving sources like NEOs or KBOs, or to overcome persistent spurious detections in highly contaminated images. shows a typical distribution of proper motion values of sources for an input set of 12 images, split into the right-ascension and declination components. The cluster towards zero proper motion is made up of background stars and galaxies.

#### Linear Motion

The filter on the linearity of the apparent motion is the most effective and reliable step in the source selection. While the proper motion computation by SCAMP is already exploited in the previous filter, it excludes outlier detections in coordinate-space from the fit, weakening the rejection. Furthermore, SSO observations with large gaps between observations may fool the linear motion filter due to the long baseline of the fit. These cases need to be identified and handled differently and make this dedicated filter step necessary.

For each SSO candidate, the right ascension and declination coordinates are fit over epoch using a weighted, linear $\chi^2$-fit. The coefficient of determination parameter $R^2$ is then defined by regarding the total sum of the squares over the sum of the squared residuals,

$$R^2 = 1-\frac{\sum_i(y_i-f_i)^2}{\sum_i (y_i-\bar{y})^2},$$

where $f_i$ describes the linear model function evaluated at the observation epochs $i$, $y_i$ are the coordinate data points and $\bar{y}$ the mean coordinate value. The fraction compares the variance of data and model to the variance of the data. If the fit describes the data well, i.e. the source exhibits linear apparent motion in this coordinate over time, $R^2$ will be close to unity. By setting the [`R_SQU_M`]{} filter parameter, sources with $R^2$ values smaller than the limit will be rejected. This is analogous to the proper motion computation by SCAMP, however, no source detections are dismissed during the fitting procedure.

Linear fits can be fooled by long temporal baselines. Any source candidate observed several times in the first hour of the observation night and one more during morning hours exhibits an $R^2$ value close to unity. To detect and treat these occurrences, source detections are checked for outliers in epoch-space. Outlier detections are defined using the median absolute deviation (MAD) of the observation epochs $E_i$,

$$\mathrm{MAD}=\mathrm{median}(|E_{i} - \mathrm{median}(E)|),$$

describing the median time step between the single observation epochs and the median observation epoch. Using the [`OUTLIER_THRESHOLD`]{} filter parameter, the user sets the upper limit of observation time between two observations in units of MADs. If any gap in observation time is larger than the set threshold, the source detections are divided into two subgroups at this gap. This may occur at multiple epochs for a single source, specifically if the observation strategy includes frequent changes between narrow- and wide-band filter observations.

Each detection subgroup is then checked for linear motion in both coordinate dimensions by the fitting described above. If all linear fits pass the $R^2$ limit, the subgroup is accepted. If a subgroup does not pass or if it contains fewer than 3 detections, it is flagged as outlier detection. If all subgroups of a source are flagged as outlier detections, the source is rejected.

illustrates this for an SSO with 5 detections, which erroneously was associated with a CR detection by the SCAMP cross-match. After identifying the CR outlier in epoch-space, the motion of the SSO is correctly evaluated as linear and the source is accepted, with the CR detection flagged as outlier.

Sources with flagged outliers can then be accepted, visually confirmed, or, as a more conservative approach, only kept if matched with a known SSO. Alternatively, all outlier detections can be rejected, while the remaining source detections are kept.

#### Trail consistency

SSOs imaged under similar observation conditions like seeing, binning, exposure times, and observation band should leave consistent trails in the images. Meanwhile, artefacts do not have any bounding conditions on their appearances, so they will likely fluctuate in size. The morpohmetric trail parameters can therefore be examined for consistency to differentiate these sources.

The consistency of the trail parameters with the hypothesis of a constant SSO imprint in the image is quantified using the semi-major and -minor axis lengths of the ellipse fit to the source by SExtractor, denoted by *a* and *b* respectively. If the observation conditions are constant, the standard deviation should be in the order of the seeing.

For each axis, the average weighted uncertainty $\bar{\sigma}$ is computed and compared to the standard deviation std of the data. The [`RATIO`]{} filter parameter sets an upper limit on the observed standard deviation with respect to the uncertainties of the morphometric parameters,

$$\begin{aligned}
    \begin{split}
        \bar{\sigma}_{a} &= \frac{1}{\sqrt{\sum_i \frac{1}{\sigma_{a,i}^2}}}\\
        {\texttt{RATIO}}_{a} &= \frac{\bar{\sigma}_{a}}{\mathrm{std}(a)},
    \end{split}\end{aligned}$$

and analogously for the semi-minor axis $b$. A smaller [`RATIO`]{} allows for larger dispersion in the ellipse parameters. depicts the relation between the weighted average uncertainty and the standard deviation of the semi-major axes of sources in a random set of single-band observations by the Gran Telescopio Canarias. The majority of sources is of astrophysical nature and exhibits small standard deviations of their axes parameters even as the average uncertainties increase. Artificial sources are located towards the bottom of the figure, indicating large size fluctuations with respect to their morphometric uncertainties.

In practice, the [`RATIO`]{} value is set by examining the distribution of the values plotted in or by trial-and-error analysis of a set of known SSOs and artefacts in the images.

The trail filter should be applied to multi-band data with caution. SSOs display exhibit wavelength-dependent fluxes depending on their spectral energy distribution (SED), leading to a change in their trail appearance under constant exposure time if the observation band is changed.

#### Bright-Sources Catalogue

Diffraction spikes and ghosts introduced by bright sources are tough challenges for the filter pipeline. They can mimic linear motion given a linear dithering pattern. The trail consistency filter can mitigate this contamination, but in the case of multi-band data or when aiming for a more conservative approach, filtering the SSO candidates in close vicinity to bright sources is the most reliable way to reject these artefacts. The [`DISTANCE`]{} parameter sets the minimum distance in arcseconds to sources that the mean candidate coordinates must have in order to pass the filter, as depicted in . The reference source catalogue can be either the same catalogue as retrieved by SCAMP, with optional lower- and upper magnitude limits applied on the reference sources via the [`MAG_LIMITS`]{} parameter, or a user-provided, local catalogue.

Optional analyses
-----------------

The [`ssos`]{} pipeline produces a catalogue of SSO candidates detections. Depending on the filter parameters and quality of the images, the output sample will be contaminated by false-positive detections to a certain percentage. To investigate the amount of artefacts surviving the pipeline, cutout images of every SSO detection can optionally be extracted using the SWARP software [@2002ASPC..281..228B]. These images allow for quick visual confirmation of the results and further aid in the identification of common artefacts and necessary filter parameter tuning.

Another standard analysis procedure is the cross-match of all SSO candidate detections with the SkyBoT service [@2006ASPC..351..367B]. The [`ssos`]{} pipeline automatically queries SkyBoT providing the image FoV and observation epoch. The service returns the metadata of SSOs that are predicted to cross the FoV at the given epoch by comparing the observation metadata to pre-computed SSO ephemerides. If a known SSO is within a user-defined cross-match radius, it designation and other properties are added to the source candidate’s metadata. In case there is more than one known SSO within the cross-match radius, the one exhibiting the closest co-linear motion to the SSO candidate is chosen by minimizing the difference in the proper motion angles between the candidate source and the known SSO.

Finally, an additional SExtractor run can be performed for all source detections to compute fixed-aperture magnitudes for each SSO, required e.g. for colour indices. SExtractor is run in dual-image mode, with a user-specified filter as reference filter. The fixed-aperture flux and magnitude are added to the output data.

Implementation and performance {#sec:implementation}
==============================

The [`ssos`]{} pipeline is written in [`python3`]{}. The source code is open-source[^8] and extensive documentation is available online.[^9] The code is further available from the Python Package Index, simplifying the install process and the handling of the [`python`]{} dependencies. The astr*O*matic software SExtractor, SCAMP, and SWARP need to be installed separately. As such, the pipeline can be installed on personal computers for the configuration setup and testing, as well as on computer grids for the batch-processing of images.

The development of the source code so far focused on simplicity of use. All pipeline steps are initiated using the [`ssos`]{} executable from the command line. The pipeline settings are read from a configuration file. The user interface by design resembles closely that of SExtractor and SCAMP, to facilitate the application for users already familiar with these software.

Upon instantiation, the [`ssos`]{} pipeline creates and manages all necessary directory structures, image and catalogue dependencies for each step. In the usual application mode, the user should not have to specify more than the directory of the FITS images and the path to the configuration file.

Some development effort has been placed on the execution speed. Improvements were achieved by profiling the python code and reducing the number of value lookups in catalogues as far as possible. Nevertheless, the main workload is handled by the efficient SExtractor and SCAMP runs, which are the analysis’ bottlenecks in terms of speed and can be multi-threaded.

Execution time depends on a variety of factors, such as the number of images searched, the pixel dimensions of each image, the number of sources present, the chosen filter steps and parameters, and finally the computer hardware. As a reference, the complete analysis of 6,132 J-PLUS images with dimensions of 2deg$^2$ in 9,500$\times$9,500 pixel took $\sim$30 hours on an 2 GHz Intel Core i5 computer architecture with 8GB RAM (MacBookPro Late 2016), around 3.5 minutes for each of the 511 input fields made up of 12 images each.

Test studies {#sec:test}
============

Applying the alpha-version of the pipeline described in @Mahlke2017, 20,221 SSO candidates with an estimated false-positive content of less than 0.05% were recovered in 346deg$^2$ of the KiDS survey. Since then, the pipeline has been re-structured to facilitate the application not only to surveys but all kinds of imaging data products. While the well-defined observation strategy of survey telescopes simplifies the application of the [`ssos`]{} pipeline due to the set dithering patterns and exposure times, it is not required for successful identification of SSOs. The versatility of the [`ssos`]{} pipeline is now illustrated using two case studies: images from the GTC OSIRIS Broad Band DR1 Catalogue, made up of PI-lead observations, and survey data from the J-PLUS DR1.

GTC OSIRIS Broad Band DR1 Catalogue {#sec:gtc}
-----------------------------------

The [`ssos`]{} pipeline is applied to 8,096 images acquired with the Gran Telescopio Canarias (GTC) OSIRIS, which were used to derive the GTC OSIRIS Broad Band DR1 catalogue (Cortés-Contreras, in preparation). OSIRIS (Optical System for Imaging and low-Intermediate-Resolution Integrated Spectroscopy) is an optical imaging and spectroscopy instrument part of GTC with a vignetted FoV of 7.8${\ensuremath{^{\prime}}}\times$7.5${\ensuremath{^{\prime}}}$. Observations in the DR1 catalogue are taken in the Sloan system broadband filters $g'$, $r'$, $i'$, and $z'$. The SSO identification is carried out on the astrometrically and photometrically corrected images provided by the Spanish Virtual Observatory. The observation footprints are depicted in .

Given the small FoV and in general large ecliptic latitude of the footprints, the output sample of SSO candidates is small, allowing for illustrative studies of the optimal pipeline filter configuration and the influence of the SCAMP cross-match radius parameter [`CROSSID_RADIUS`]{}, as visual confirmation of the results is feasible.

The pointed observations apply varying observation strategies. A single observation run can include different time-spans, singular or a multitude of exposures, different exposure times, and filter configurations. To apply the [`ssos`]{} pipeline, the 8,096 exposures are therefore divided into logical groups. Each logical group is made up of all exposures with overlapping FoVs taken in a single night. If a group is made up of 3 or fewer exposures, it is discarded, as the linear motion cannot be detected reliably. 712 exposures are rejected on the basis of this criterion. An additional 402 exposures had to be discarded due to a large contamination by imaging artefacts, introduced during the detector read-out or caused by telescope tracking failures.

In total, 6,982 footprints are searched, divided into 420 logical groups. The observation epochs span the time from April 2009 to July 2014. The smallest number of exposures per group is 4, the largest group consists of 218 images, while 316 groups have 20 or fewer exposures. The average exposure time per group ranges from 17s to 1h.

For each exposure, a weight image is computed from the flat field images available online and passed to SExtractor for source detection.

The small FoV presents a challenge for SCAMP, as the reference catalogues might contain too few sources to compute the astrometric solution. The standard reference catalogue is set to SDSS-R9. For 129 groups, the reference had to be switched to the GSC-2.3 catalogue due to coverage reasons [@2008AJ....136..735L; @2012ApJS..203...21A].

A further challenge arises due to the absence of an underlying observation strategy. The cross-match radius has to be chosen on a per-group basis as the cadence between single exposures for each observation night varies. The radius is computed on-the-fly so that the average cadence between exposures allows for the observation of SSOs with proper motions up to 50[$^{\prime\prime}$]{}/h, 100[$^{\prime\prime}$]{}/h, and 150[$^{\prime\prime}$]{}/h. These three values are chosen to test the hypotheses that choosing a small radius misses fast moving SSOs if there are gaps in the observation epochs, while slower SSOs are recovered more reliably. An incomplete but pure output sample is expected. On the other hand, larger [`CROSSID_RADIUS`]{} values should yield a higher chance of recovering a complete SSO set, however, the amount of contamination due to spurious detections should increase as well.

For each of the three runs with varying cross-match radii, three different filter parameter settings are set, again aiming to investigate the effect on the completeness versus purity relation. The *default* pipeline settings are listed in \[app:settings\]. Only sources with more than 3 detections are accepted and the trail consistency filter is off. The *HYG* database is chosen as reference catalogue for bright sources in the FoV. It contains all stars in the Hipparcos, Yale Bright Star, and Gliese catalogues.[^10]

In a second step, denoted *trail*, the trail consistency filter is turned on and [`RATIO`]{} is set to 0.25, to reject artefacts more conservatively. A final run, *strict*, increases the [`RATIO`]{} to 0.5, [`R_SQU_M`]{} from 0.95 to 0.97, and rejects sources with 6 or fewer detections.

With the three chosen cross-match radii and three filter parameter settings, a total of 9 analysis runs is analysed and the output samples are visually inspected. Between the SSO candidates resulting from the *default* runs, there is a large overlap as the input images remain the same, but the different cross-match radii include some sources and exclude others.

After visually confirming the nature of all SSO candidates, 204 distinct SSOs are identified. One SSO was recovered in 53 exposures. No [`CROSSID_RADIUS`]{} setting successfully recovered all 204 SSOs in the *default* setting. displays the results of the 9 analysis runs. The most complete output set is achieved with the [`CROSSID_RADIUS`]{} set to 100[$^{\prime\prime}$]{}/h, however, the difference between the sets is within a few percent. On the other hand, the number of recovered false-positives is proportional to the radius. The contamination increases from 19% to 30% between the 50[$^{\prime\prime}$]{}/h and 150[$^{\prime\prime}$]{}/h runs. For the two stricter parameter settings, the differences in the output samples disappear. Ultimately, a [`CROSSID_RADIUS`]{} of $\sim$100[$^{\prime\prime}$]{}/h is recommended for most use cases, as the majority of SSOs should not exhibit proper motions faster than this, and the degree of contamination is limited.

  Comet   Inner MB   Middle MB   Outer MB   Trojan
  ------- ---------- ----------- ---------- --------
  4       21         17          20         1

  : The distribution of cross-matched SSOs over SSO classes. MB describes Main-Belt asteroids.[]{data-label="table:skybot_gtc"}

Finally, the number of known SSOs in the output sample is determined. GTC images from early observation epochs may carry an error of up to 15 minutes on their [`DATE_OBS`]{} header keyword. While this presents no problem for the identification pipeline, the SkyBoT cross-matching may be inaccurate. Therefore, the [`CROSSMATCH_RADIUS`]{} parameter in the SkyBoT queries was increased from 10[$^{\prime\prime}$]{} to 40[$^{\prime\prime}$]{}.

63 of the 204 SSOs are matched to known SSOs using the SkyBoT service. Their classes are given in . The 4 observed comets are recovered from targeted observations, as well as one observation of (300163) 2006 VW139. All targeted observations of SSOs are successfully recovered, though the targets of some observation runs were left undefined by the observer. The vast majority of serendipitously observed SSOs are from the Main-Belt, with one Jupiter Trojan observed as well. compares the absolute proper motion values of the cross-matched SSOs as computed by SkyBoT and as recovered by the [`ssos`]{} pipeline. A successful cross-match is indicated by a ratio close to unity between these two values, which is the case for most matches. Ratios far from unity point to co-incidental cross-matches.

All 2,828 detections of the 204 detected SSOs were reported to the MPC. Further scientific exploitation of this dataset will be presented in Cortés-Contreras (in preparation).

J-PLUS DR1 {#sec:jplus}
----------

The J-PLUS (Javalambre Photometric Local Universe Survey) is an on-going optical survey executed at the Observatorio Astrofísico de Javalambre (OAJ) [@cenarro2018j]. The observation strategy entails 36 observations in 12 bands taken within a few hours, using the JAST/T80 telescope with an FoV of 2deg$^2$. The whole survey aims to image 8,500deg$^2$. depicts the footprints of the J-PLUS first data release from July 2018 containing 511 fields imaged between November 2015 and January 2018.[^11] The observations imaged 1022deg$^2$ of the sky with limiting magnitudes ranging from 20 to 23.

For each of the 511 J-PLUS fields, the 3 dithers per band are combined to 12 images, which are retrieved from the online repository, including their weights and image masks. When combining the 3-step dither pattern, a frame of zero value pixels is padded around a reference exposure. This frame is later removed with the supplied image masks. The weight images are supplied to SExtractor during source detection.

At a first inspection of the images, it is found that the reduction pipeline successfully removes most imaging artefacts. The [`CROSSID_RADIUS`]{} is therefore set to 125[$^{\prime\prime}$]{}/h, as the amount of artificial sources is expected to remain low. Applying the [`ssos`]{} pipeline on a test set of 6 tiles showed that the contamination of the final dataset is almost exclusively made up of stars. There is a probability that stars will mimic linear motion when in reality their centres are stochastically displaced due to seeing, which in combination with the large FoV of the JAST/T80 camera leads to this contamination.

The reference catalogue for the astrometric matching by SCAMP is the 2MASS All-Sky Catalogue of Point Sources [@2003yCat.2246....0C].

The [`ssos`]{} pipeline is applied with two filter settings. The *default* setting is analogous to the one used in the GTC study, \[app:settings\], except that the bad pixel filter is turned off due to the low level of contamination. The *strict* setting increases the required number of detections per source to 5 and the [`R_SQU_M`]{} parameter to 0.97, to reject the stellar contaminants. The HYG database is used again as reference catalogue for bright sources in the FoV.

The results of the two applications are given in . The *default* sample contains 4,606 SSO candidates, with an estimated contamination of (2.0$\pm$0.2)%. The *strict* sample is made up of 3,696 candidates and has a false-positive ratio (FPR) below 0.68%. The number of false-positive detections for each sample is derived from visual inspection of a random sub-sample of 500 (441) SSOs for the *default* and the *strict* sample respectively. 10 artefacts are found in the *default* sub-sample, none of which is present in the *strict* sub-sample.

Of the 4,606 SSO candidates in the *default* sample, 3,932 are matched within 10[$^{\prime\prime}$]{}with a known SSO by SkyBoT. The distribution of the cross-matched SSOs over the SSO populations is given in . The majority of known SSOs is identified in the Main-Belt, while also members of the close NEAs and far KBOs are recovered. displays the relation between the computed and recovered proper motion values of the cross-matched SSOs, showing a good agreement in general.

              SSO Candidates   FPR
  ----------- ---------------- ----------------
  *default*   4,606            (2.0$\pm$0.2)%
  *strict*    3,696            $\leq$0.68%

  : The number of SSO candidates and the derived false-positive ratio (FPR) for the two pipeline runs on images of the J-PLUS survey. []{data-label="table:jplus_results"}

  NEA   MC   MB      Trojan   KBO
  ----- ---- ------- -------- -----
  3     31   3,887   10       1

  : The distribution of cross-matched J-PLUS SSOs over the SSO classes. NEA describes near-Earth asteroids, MCs are Mars-Crosser, MB are asteroids in the Main-Belt, while KBO are Kuiper-Belt objects.[]{data-label="table:skybot_jplus"}

After photometric calibration of the recovered magnitudes, 34,014 detections of 4,596 SSOs were reported to the MPC.

Surveys like J-PLUS are prime candidates for the application of the [`ssos`]{} pipeline, as the strict observation strategy definition offers precise settings of the filter parameters. Especially imaging campaigns covering areas close to the ecliptic will often yield a large number of SSO candidates, rendering visual confirmation of the nature of all candidates impossible. Therefore, a detailed discussion of the statistical derivation of the false-positive ratio is given in \[app:fpr\].

Upper limits for the completeness of the SSOs sample could be derived by regarding the fraction of recovered, known SSOs using SkyBoT. The large ephemeris uncertainties of SSOs with few recorded observations and the uncertainty on the predicted visual magnitude have to be accounted for.

Conclusion and Outlook {#sec:conclusion}
======================

Presented is an easy-to-implement, light-weight SSO recovery pipeline which can be applied to increase the science-yield of any kind of imaging data. The [`ssos`]{} pipeline offers a high degree of versatility thanks to the underlying SExtractor, SCAMP, and [`python`]{} software. A flexible set of filter algorithms developed in various use cases offers a wide range of SSO identification methods.

The two presented test cases illustrate the potential of the pipeline. The observation strategies and telescope set-ups differ substantially, yet the [`ssos`]{} pipeline is successfully applied in both cases with only minor changes in the set-up. The wide range of classes of recovered SSOs in the J-PLUS survey, from close NEAs to distant KBOs, displays that the method does not discriminate between the populations, as the only prerequisite are subsequent observations of an overlapping area in the sky.

The [`ssos`]{} pipeline becomes more versatile with each new application case, as different observation strategies or data product formats are accounted for. Future development will focus on easier filter parameter testing as well as increasing the science yield. The former may be achieved by automatically evaluating the number of source candidates in each step for different filter setups, and comparing them to a sample of previously identified artefacts and SSOs. The most promising improvement in terms of SSO science would be the successful cross-matching of SSOs recovered in subsequent observation nights, greatly increasing the recovered arc of the orbit. Finally, the prior knowledge of the ephemerides of known SSOs provided by SkyBoT can be exploited more by cross-matching sources before applying the filter pipeline, and extracting matches with co-linear movement right away.

Acknowledgements {#acknowledgements .unnumbered}
================

MM is funded by the European Space Agency under the research contract C4000122918.

MM would like to thank the EDAS Helpdesk at ESAC for their advice on the derivation of the false-positive sample ratios.

Based on data from the GTC Public Archive at CAB (INTA-CSIC).

Based on observations made with the JAST/T80 telescope for the J-PLUS project at the Observatorio Astrofísico de Javalambre, in Teruel, owned, managed and operated by the Centro de Estudios de Física del Cosmos de Aragón.

This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France.

H. Bouy acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 682903, P.I. H. Bouy), and from the French State in the framework of the “Investments for the future” Program, IdEx Bordeaux, reference ANR-10-IDEX-03-02 .

Estimating the false-positive ratio {#app:fpr}
===================================

Applying the [`ssos`]{} pipeline to a large set of imaging data often results in an amount of SSO candidates where visual confirmation of their true nature is no longer feasible. The number of false-positive detections has to be estimated by limiting the inspection to a smaller subset. In the following, the statistical foundation of this estimation is described in brevity. Further information can be found e.g. in @jaynes2003probability [@rice2003mathematical].

The recovered SSO candidates form a dichotomous population of size $N$. They can be divided into $K$ artefacts and $(N-K)$ SSOs.[^12] A random sample of size $n$ drawn without replacement from this population forms a vector of Bernoulli random variables $\textbf{X}=(X_{1}, X_{2}, \dots, X_{n})$, where

$$X_i =
  \begin{cases}
                                   0 & \text{if SSO} \\
                                   1 & \text{if artefact} \qquad i\in\{1, \dots ,n\}
  \end{cases}$$

The number of artefacts $k$ in the sample $n$ is then $$k = \sum_{i=1}^n X_i$$

The probability density function of $k$ is described by the hypergeometric distribution $$h(k\,|\,N, K, n) = \frac{\binom{K}{k}\binom{N-K}{n-k}}{\binom{N}{n}}.$$ The notation here emphasizes that the probability distribution $h(k\,|\,N, K, n)$ of the sample statistic $k$ depends, among others, on the unknown population parameter $K$. In the following, this notation is dropped and the dependency is implicit.

Each population member has a probability of $n/N$ of being included in the sample. The linearity of the expectation then yields the expected number of artefacts $E(k)$ in the sample, given by the contribution of the expectation of all $K$ artefacts; $$E(k) = K\frac{n}{N}.$$

The derivation of the sample variance $s^2(k)$ is more challenging as the population members are dependent variables. While the result is given here, the reader is referred to the literature for detailed steps: $$\label{equ:vark}
    s^2(k) = \frac{nK(1-\frac{K}{N})(N-n)}{N(N-1)}.$$

The sample statistic $k$ is now used to estimate the population parameter $K$. An intuitive estimator of $K$ is defined by using the method of moments; the ratio of artefacts in the sample is assumed to be equal to that in the population: $$\frac{k}{n} \approx \frac{K}{N} \implies K \approx k\frac{N}{n}$$ The expectation of $K$ is then given by $$\mu(k\frac{N}{n}) = K.$$ As the estimator is unbiased, its variance is equal to the mean square error of $K$, where is used: $$\label{equ:sigma_K}
    \sigma^2(k\frac{N}{n}) = s^2(k)\cdot\Big(\frac{N}{n}\Big)^2= (N-K)\frac{K}{N}\frac{N-n}{N-1}$$ It is apparent that, as the sample size $n$ increases towards $N$, the variance of $K$ decreases. The intuitive conclusion is that a larger test sample leads to a more accurate estimation of the population contamination.

The description above is valid for all population parameters $N, K$ and sample parameters $n$. However, special care should be taken if the sample does not contain any artefact; i.e. $k=0$. In this case, yields an uncertainty of 0, which is underestimated unless $n=N$.

Instead, the upper limit of the the 95% confidence interval on the false-positive ratio of the population can be estimated by approximating the hypergeometric distribution with the binomial distribution. This is valid for large $N/n$ values. The upper limit is then given by

$$\mathrm{Pr}(k=0) = (1-p)^n \overset{!}{\geq}0.05$$

If all inspected sources in a sample of size $n>30$ are SSOs, the false-positive ratio is therefore well approximated by $3/n$, which is known as the rule of three [@eypasch1995probability].

To demonstrate the application, the derivation of the false-positive ratio of the J-PLUS survey pipeline runs presented in is given in .

            N       n     k    K          FPR
  --------- ------- ----- ---- ---------- ----------------
  Default   4,606   500   10   92$\pm$9   (2.0$\pm$0.2)%
  Strict    3,696   441   0    0$\pm$0    $\leq$0.68%

  : Calculation of sample contamination for the two J-PLUS pipeline runs presented in . FPR refers to the false-positive ratio. Note that for the Default run, the values refer to the expectation and standard deviation, while for the Strict run, the rule of three was used to derive the upper limit of the 95% confidence interval.[]{data-label="table:jplus_fpr"}

[`ssos`]{} pipeline settings file {#app:settings}
=================================

[^1]: <https://minorplanetcenter.net/>

[^2]: <https://pypi.org/project/ssos/>

[^3]: <https://www.astromatic.net/>

[^4]: Until FITS header keywords are truly standardized, this tedious task cannot be reliably handled automatically.

[^5]: <https://sextractor.readthedocs.io/en/latest/>

[^6]: While the [`CROSSID_RADIUS`]{} is defined in arcsecond in SCAMP, in the context of SSO identification it is more convenient to convert it to an equivalent proper motion limit by dividing by the time between the beginning of the two subsequent exposures. From here on, the [`CROSSID_RADIUS`]{} is therefore given in [$^{\prime\prime}$]{}/h.

[^7]: Bad pixel masks and increasing the minimum number of pixel per source in the SExtractor configuration are effective ways to reject artificial sources caused by broken CCD pixel, but the former might not be available and the latter might impede the detection of faint SSOs.

[^8]: <https://github.com/maxmahlke/ssos>

[^9]: <https://ssos.readthedocs.io>

[^10]: <http://www.astronexus.com/hyg>

[^11]: <http://archive.cefca.es/catalogues/jplus-dr1/>

[^12]: *Artefacts* includes any kind of false-positive result, not necessarily only source detections caused by imaging artefacts.
