text new page (beta)
English (pdf)
Article in xml format
Article references
Automatic translation
Send this article by e-mail
Cited by SciELO 
Similars in
SciELO 
Permalink
1. Introduction
The geostationary satellites play a crucial role in science and telecommunications. The geostationary environment (GEO) is now populated not only by active satellites but also by those satellites whose fuel has been consumed, as well as fragments resulting from breakup of satellites. These bodies called Space Debris (SD) together with active satellites are grouped in the GEO forming several populations distributed along a torus in the Laplace stable plane which is inclined about 7.3° with respect to the equatorial plane, with the nodes near the equinoxes called the GEO ring (McKnight & Di Pentino 2013; Flohrer 2017). A precise inventory of these populations, together with accurate determinations of their spatial distribution and orbital dynamics is relevant to plan for positioning new geostationary satellites.
Observations are fundamental inputs for accurate orbit determinations, allowing to correct for the uncertainties of theoretical models based on simplified assumptions such as the Kepler motion.
For this reason, several surveys have been conducted during the last three decades by telescope and radar networks from different geographical positions. Radar observations give accurate distance estimations, while telescope observations allow for an accurate determination of the remaining two spatial coordinates. Particularly, the radar surveys are focused on low Earth orbit (LEO) objects while the optical surveys cover other kinds of objects, including GEO objects (Pelton 2016).
Some previous telescope surveys are those of Molotov et al. (2008); Dicky et al. (1993), intended to warn against a SD collision; Alby et al. (2004); Laas-Bourez et al. (2009, 2011), focused on orbit determination; Sabol et al. (2004) dealing with the detection of SD smaller than 1 cm; Schildknecht et al. (2004); and Silha et al. (2017) about disentangling of different populations of SD. Radar observations were carried out by Mehrholz (2001) to characterize objects in LEO. Currently, none of these surveys are public since they are specifically designed for particular needs of governments, with the exception of the International Scientific Optical Network (ISON) (Molotov et al. 2008).
Nowadays the knowledge of SD in the GEO ring has been extended by developing different theoretical models and from new observations. It is known that the main sources of SD in this region were two catastrophic events: the breakup of the Ekran-2 in 1978 and the breakup of the Titan rocket in 1992 (Rykhlova et al. 1997; Schildknecht et al. 2004; Bordovitsyna & Aleksandrova 2010). Nevertheless, it is necessary to understand the dynamics of SD more accurately in order to improve aspects such as alert against possible collisions, repositioning GEO satellites or placing new satellites that come into service. Reliable observations to support the development and validation of dynamical models are important in order to build catalogues with the most accurate information supporting mission analysis and design.
The main goal of this paper is to present the astrometric method applied to a preliminary set of observations as a proof-of-concept, to show the astrometric precision achieved, the highest reported to date. In order to reach this high precision, first we perform an astrometric reduction process adequate to the observations obtained of the CIDA Survey of the GEO Ring, which is a public program carried out with the 1-m reflector telescope at the Venezuela National Observatory (VNO) near the terrestrial Equator. This is a favorable condition for an efficient and accurate astrometric observation of orbiters in the GEO ring. Second, we compute and correct for the systematic errors introduced by the optics and electronics of the observational instrument represented by the distortion pattern (Stock 1981; Abad 1993).
The paper is organized as follows: In § 2 we describe the survey strategy. In § 3 we present and describe the astrometric method applied for coordinate determination, which allows us to obtain a low mean relative error. In § 4 we summarize our analysis and conclusions.
2. The CIDA survey of the geo ring
From its canonical definition, a perfect geostationary orbiter has an orbit with null eccentricity (e) and inclination (i), and furthermore a semimajor axis (a) of 42, 164 km, making the orbital period equivalent to a sidereal day (Spalding 1995). However, most of active satellites are geosynchronous with significant nominal inclinations and semimajor axis intentionally designed by control space agencies. Additionally, the position of orbiters in the GEO are perturbed by the joint action of conservative and non-conservative forces, the most important of which are the strongly non-spherical gravitational potential of the Earth, the solar radiation pressure (SRP), and the gravitational attraction of the Sun and the Moon (Battin 1999; Chovotov 2002). As a result the GEO ring is composed by orbiters having e ≤ 0.2, 39,664.0 ≤ a/km ≤ 45,314.0 and i ≤ ± 70° as supported by observational data and theoretical modeling (e.g. McKnight & Di Pentino 2013; Anderson et al. 2016). Nevertheless, most of these objects have inclinations between i ± 20° (e.g. Flohrer 2017). In consequence, the projection of the GEO ring on the sky results in a torus of 40° width centered at the stable plane of Laplace.
The goal of the CIDA Survey of the GEO Ring presented here is to perform an accurate, low-cost, long-term and large-scale astrometric survey of 1/3 of the GEO ring. The observations are performed with the 1-meter Carl Zeiss reflector telescope with geodetic coordinates lat = 8◦47′21.9516′′, lon = 70◦52′12.7776′′ and alt = 3584m, considering the model of the WGS84 ellipsoid. The telescope is provided with a hybrid anastigmatic aspherical primefocus corrector (Della Prugna & Schenner 2009) and a ProLine 4260 − E2v56 CCD FLI camera of 2k × 2k pixels, resulting in a field of view (FOV) of 0.33 square degrees with a plate scale of 0.546 ′′/pix.
The survey follows a simple observational strategy: the observations are performed at zenith angles z < 60◦ in order to minimize the effects of the atmospheric distortions on the astrometric accuracy. Considering the geographical latitude of the telescope and the constraint on the zenith angle, the CIDA survey of the GEO ring covers 1/3 of the complete GEO ring. It covers the range 10.90 W < λ/◦ < 130.8 W including the unstable λu2 = 11.5◦W and the stable λs2 = 105.3◦W points. Then, the 120◦ in longitude of the GEO ring which are observable from the VNO are divided into 120 strips along right ascension (α) with a length of 120◦ in α and a width of 20′ in declination (δ). We denote those as the α-strips. Each α-strip is divided into 360 observing windows with a size of 20′ × 20′, corresponding to the FOV of the instrument. Therefore, the complete survey is composed a total of 43,200 observing windows.
During a typical night of observation, the survey focuses on a single α-strip, and the telescope is fixed at certain hour angle (HA) in the range −4 < HA/h < 4 ( without sidereal tracking). The result is a set of 40 images with 10 s of integration time each one. Then, the telescope is moved to the next observing window in the same α-strip and a new set of 40 images is obtained. This observational strategy was designed in order to ensure the overlaping between consecutive observations needed to use the block-adjustement method (Stock 1981) for the astrometric reduction.
With the proposed time sampling i.e., 40 consecutive images of the same observing window in a period of approximately 11.03 minutes the relative motion of the different orbiters in the GEO ring can be obtained. These observed objects may correspond to SD or active satellites whose orbits are being corrected. From February, 2015 to March, 2017 a set of 20, 000 images was obtained in the α-strip centered at δ = 0◦. Figure 1 shows the number of images obtained in each observing window. The time sampling chosen in this work allows a preliminary characterization of the relative motion of the orbiters in the GEO ring.
Fig. 1 Number of observations for each observing window corresponding to objects in the interval 8.6W < λ/°< 117.5W and −0.165 < δ/°< 0.165. The available number of observations allows the determination of the relative motion of orbiters in the GEO ring over different time scales.
Due to the fact that objects placed inside the GEO ring are quasi-fixed with respect to an Earthrotating reference frame, they appear in the 10 s integration images as point sources while, due to the sidereal motion, the stars cross through the FOV appearing as streak sources whose total length depends on the integration time (Warner 2006). The width of the streak sources as well as the FWHM of the point sources depend on the seeing conditions (Warner 2006). An example of the resulting image is shown in Figure 2. The integration time was defined in order to improve the sensitivity limit allowing the detection of orbiters of small size, and increasing the number of detected stars that are useful for the definition of an accurate astrometric reference frame. Additionally, the detector is aligned according to the coordinates (α, δ) in order to obtain the stellar images aligned along the direction of α, which helps to improve the accuracy of the astrometric solutions.
Fig. 2 Example of an astrometric observations performed in the window λ=77◦W and δ=0◦ during the night of January, 14th 2016 with an exposure time of 10s. Black rectangles enclose some stars used as part of the astrometric reference frame. Black circles indicates the satellites EchoStar1 and QuetzSat1 placed inside the GEO ring. Black diamonds indicate other satellites or SD crossing the FOV producing trails in any direction.
The combination of instrumental sensitivity and the selected integration time allows the detection of point sources with visual magnitudes mV < 20, which corresponds to SD belonging to the GEO ring with mean sizes as small as 0.3 m, according to the estimated limiting magnitudes and sizes from the USSPACECOM catalogue4 (assuming an albedo of 0.08).
3. Astrometric processing pipeline
The large volume of data acquired during the survey is reduced in order to detect and to distinguish different kinds of objects present in each image. We use an adequate astrometric treatment to compute their equatorial coordinates in an accurate and efficient way. In this section we describe the procedure, techniques, and general scheme of the astrometric processing pipeline used in this work.
3.1. Reduction
The images are reduced by subtracting the readout level of the detectors obtained from bias frames, and by correcting the pixel-to-pixel quantum efficiency differences by dividing the images by normalized dome flat frames. Both bias and flat frames are obtained at the beginning and the end of each observing night. No dark-current correction is applied because the instrument produces a negligible dark current during the selected integration time. Other effects such as hot pixels and cosmic rays are detected and substituted by the mean value of the surrounding pixels. All these reductions are performed through IRAF based sub-routines included in the IMRED package.
3.2. Detections
The detection of stars as well as objects from the GEO ring was performed as follows. Each image is collapsed in the direction of x; this means that for a fixed y (δ coordinate) all counts in the pixels of all rows along the x axis are added. The resulting count addition is plotted against the y axis, as shown in Figure 3. In this representation, both stars and objects from the GEO ring look like point sources and the corresponding y coordinates of their centroids are obtained by fitting the distribution to a Gaussian point spread function (PSF) (Sinachopoulos & Seggewiss 1990).
Fig. 3 Distribution of the signal along the δ (y) coordinate for a typical observation. Orbiters in the GEO ring and stars are seen as PSF. The convolution process is applied where there is a significant signal.
In order to identify streak and/or point sources in every row where a significant signal is detected, we apply the convolution method. For each yi, where i is the i−th row, a box with the length of a streak and the width as function of the FWHM of the peak (shown in Figure 3) is defined. Then, the box is located at yi and moved along the x axis pixel to pixel as shown in Figure 4. Panels 1, 2 and 3 of this figure show a sequence of this displacement over a yi row. The counts inside the box along the x coordinate are collapsed where the corresponding value x of the sum is the beginning of the box (left side).
Fig. 4 Convolution process over a row with significant signal. Panels 1, 2 and 3 represent a sequence of a displacement of a defined convolution box with a length of 280 pixels (length of the streaks proportional to exposure time of 10 s) and a width related to the FWHM of the streak. Bottom panel shows the resulting convolution. The pyramid represents the convolution for a complete streak, being maximum when the box covers the entire streak, labeled as (a); the respective coordinate x is represented by the left side of the box. If the convolution box is passing through a portion of the CCD with no significan signal, the convolution has a minimum value (b). When the box finds a point source, the signal is constant over the path through it, and the resulting convolution shows a plateau (c).
Consequently, when the convolution finds a streak, the counts added inside the box are not fixed but vary along the path through the streak, resulting in a pyramid pattern as is illustrated in Figure 4. The maximum occurs when the box encloses the complete streak, labeled as (a) in Panel 1 and bottom of Figure 4. On the other hand, if the convolution finds a point source the number of counts is fixed while the box is moved through it and the result of the convolution is a flat pattern, labeled as (c) in Panel 3 and bottom of Figure 4. In Panel 2 of this figure there is no source with significant signal, and the convolution shows the minimum value, labeled as (b). Therefore, using the convolution method over an image it is possible to distinguish what kind and how many object exist in each y and x coordinate of the CCD. The (x, y) coordinates obtained in the convolution process act as initial conditions in the determination of centroids computed with the respective fitting functions.
Once the convolution is made the fitting process, for each object detected begins. To generate a profile to be adjusted, an area centered at the (x, y) coordinates of each object is defined as (FWHM×FWHM) for a point source and (5/3 box length×FWHM) for a streak. Then, this area is collapsed in the y direction and the count number for each column is computed, adding the counts in every pixel along the column, that is, for a fixed x and over all pixels with y inside the defined area. In this representation, point sources in the image appears as a PSF profile, while a star generates a flattened peak, which we named a tepui distribution (Abad et al. 2004). In Figure 5 we show an example of a PSF (left panel) and a tepui (right panel) profiles.
Fig. 5 Distribution of the signal along the α (x) coordinate. Left panel shows a point source PSF. Right panel shows a streak represented as a tepui.
The detection of such patterns and the computation of their corresponding centroids in the x coordinate were performed, respectively, by fitting a Gaussian PSF and a tepui function. The latter is defined by Abad et al. (2004) as:
where A is the amplitude after applying a normalization factor, b1 and b2 represent the tilt of the edges of the image profile and c is related to the semiwidth of the profile of the image. The x centroid of the PSF is computed as the x coordinate where the maximum of the PSF occurs, and the corresponding x centroid of the tepui function as the x coordinate of the mid-point of its length.
Figure 6 shows the complete tepui profile of a star in the x coordinate and the corresponding fitting to a tepui function. Following this procedure, we have obtained mean accuracies of 1.3 × 10-2 with a standard deviation (STD) of 3.3 × 10-3 for the faint limit, and 6.7 × 10-3 with a STD of 3.4 × 10-3 for the bright limit measured in pixels in the x direction. Analogously, we have obtained 4.6 × 10-3 with a standard deviation (STD) of 3.2 × 10-3 for the faint limit and 4.2×10-3 with a STD of 3.1×10-3 for the bright limit, measured in pixels in the y direction.
3.3. The Astrometric Distortion Pattern
It is well known that mainly the optic and electronic elements of astronomical instrumentation introduce a characteristic astrometric distortion pattern (DP). It can be calculated from the overlapping bewteen guided observations or with the data of the observations obtainded in this survey. In order to characterize and correct for such distortions, we use the block-adjustment method developed by Stock (1981). This procedure uses the spatial overlapping between consecutive observations to find the solution. It compares the observed position of a set of reference stars with their expected position obtained through the projection of their equatorial coordinates onto the focal plane of the telescope, allowing to find the linear solution and the associated residuals for the overlapping and reference equations.
This method allows the computation of the DP. The linear terms, such as the rotation of the detector (related to the equatorial coordinates), the displacement of the focal distance, and the changes of the focal scale, are taken into account in the linear solution. Abad (1993) showed that there are two ways to calculate the DP with the application of the blockadjustment method: (1) introducing coefficients as common unknowns in the equations, and (2), analyzing the final residuals as function of the measured coordinates. We use here the second approach, since second and greater order terms are included, representing the deformations of the telescope, which means that a previous knowledge of them is not required.
The residuals obtained are translated as systematic errors and represent non-linear terms of the solution. There is a systematic trend that depends on the position of the measured coordinates (x,y), calculated through an interpolation of the sliding weighted polynomial described in Abad (1993). This polynomial calculates numerically the value that corresponds to each coordinate. In an iterative process, new residuals are calculated for each step until convergence is reached. The resulting final mean errors are equivalent to 0.18′′ and 0.10′′ for right ascension and declination respectively, considering for the solution of the DP stars brighter than 13 mag from the UCAC-4 catalogue (Zacharias et al. 2013). Comparing this with the mean errors of 0.29′′ in α and 0.28′′ in δ obtained before the DP correction, a considerable decrease is noted. Is important to remark that if the optic and/or the electronic changes, it is necessary to determine a new DP.
The DP calculated is shown in Figure 7. The vector field represents the individual residuals in arcseconds after the application of a linear solution, as well as the coordinates (x, y) in pixels. Notice the strong distortion at the corners of the CCD due to field distortion, the biggest component of the DP, which is radial and highest at the borders of the optical field.
Fig. 7 Distortion pattern of the telescope and camera system used during the CIDA Survey of the GEO ring. The modulus and direction of each arrow are given by the composition of x and y residues, in pixels, respectively. The plate scale is of 0.5´´/pix, which corresponds to 100 m in the GEO ring.
We remark that the uncertainty of 200 m in angular position at the mean distance of the GEO ring translates into an uncertainty of 1´´, so a DP correction is mandatory to obtain a good astrometric reduction. Also if offers the advantage of allowing to estimate an orbit propagated with a powerful numerical and/or analytic integrator, since a high precision in the determination of coordinates as initial condition is obtained. In our case, the DP produces an uncertainty of between 30 and 25 m according to the coordinates. Once the topocentric equatorial coordinates are calculated, they are corrected by the DP in order to eliminate the systematic errors produced mainly by the optical and electronic systems.
3.4. The Topocentric Equatorial Coordinates
In addition to the reduction process and the accurate computation of the centroids, we considered the following issues in order to achieve a high precision in the computation of topocentric coordinates referred to ICRS : (1) the definition and accurate determination of the instant of time to which the measurements from an image are referred; (2) the choice of an accurate reference frame system and, (3) the corrections of the systematic errors due to astrometric distortions introduced by the instruments.
Another issue in the reduction process is the irregularity at the borders of the streaks caused by the inhomogeneity in the opening and closing of the CCD shutter (Montojo et al. 2011). As a consequence, the beginning and the end of the tepui function may have different slopes and the coordinates of the streak may be not well determined. In order to avoid this issue we chose the time instant equal to half of the integration time to measure the coordinates x and y.
The equatorial coordinates were referred to an external and independent system of the Earthrotating reference frame. We selected the UCAC-4 catalogue, for which the internal uncertainty is between 15 − 100 mas per coordinate, depending on magnitude, and a mean number between 7 < N < 60 stars per image were detected during the survey, depending on galactic latitude.
In order to determine the equatorial coordinates, it is necessary to measure at least 3 stars, that is, to detect a minimum of three complete streak sources in each image. Therefore, the exposure time of 10 s was defined according to the CCD size in order to guarantee the required number of comparison stars.
In Figure 8 we show the distributions of the final errors of the equatorial coordinates α (top panel) and δ (bottom panel) obtained from a sample of 40 observations in a window, taking into account all the stars included in the solution matrix. It is noteworthy that the mean error in right ascension is α=0.16´´ and δ=0.11´´ for declination. We want to remark that our results show the smallest error reported to date. For example, Montojo et al. (2011) also use a tepui function for the detection of stars, but the DP is not included. They report mean errors of α=0.5´´ and δ= 0.2´´. Sun & Zhao (2013) obtain errors of α=1.5´´ and δ=0.6´´ using the mathematic morphology method. Núñez et al. (2015) reach errors of 0.8´´ for both equatorial coordinates α and δ using the image deconvolution method.
4. Conclusions
In this paper we have presented a method that allows to reduce the mean astrometric errors in positions of objects in the GEO ring. To test this method we used the CIDA survey of the GEO ring, an accurate, low cost and long-term astrometric survey covering one third of the GEO ring with a telescope of medium aperture, and without system tracking. The observations used were located in the interval 8.6W ≤ λ/◦ ≤ 117.5W and −0.165 ≤ δ/◦ ≤ +0.165.
The astrometric processing had a low computational cost. The reduction of one field of 40 observations was carried out in approximately 10 min, automatically, on a personal Intel Core i7 laptop computer. We obtained a precision better than 0.1 pix in the detection of stellar sources and objects in the GEO ring.
With the computation of the DP and the elimination of systematic errors we obtained a high precision for the topocentric equatorial coordinates. The errors in angular position are as low as (0.16′′, 0.11′′) in α and δ respectively, showing an improvement by a large factor compared to the (0.5′′, 0.6′′), (1.5′′, 0.6′′) and (0.8′′,0.8′′) values reported previously by Montojo et al. (2011), Sun & Zhao (2013) and Núñez et al. (2015), respectively. Our errors are equivalent to 25 m of displacement at the mean distance of the GEO ring. It is not necessary to recalculate the DP unless the optical system changes.
The contribution of our work is to use optical sensors coupled with a rigorous astrometric processing in order to obtain a high precision for the topocentric equatorial coordinates for objects in the GEO ring. The represents a considerable improvement over previous researches.
