Orbit estimation apparatus, orbit estimation method and orbit estimation program for artificial celestial body

Abstract

To provide an orbit estimation apparatus, an orbit estimation method and an orbit estimation program for an artificial celestial body, which can accurately estimate an orbit element, even by estimating a short circular arc through an astronomical survey.SOLUTION: An orbit estimation method for a plurality of artificial celestial bodies, which performs orbit estimation by Bayes' estimation using MCMC, on the basis of an astronomical survey, comprises: a low-accuracy posterior distribution calculating step of obtaining a low-accuracy posterior distribution of orbit elements and sizes of celestial bodies from short circular arcs of respective orbits of a plurality of artificial celestial bodies obtained through an astronomical survey; an artificial celestial body identifying step of identifying a survey of the same artificial celestial body by comparing the posterior distributions with each other, with respect to different surveys, of surveys through the astronomical survey; and a high-accuracy posterior distribution calculating step of obtaining a high-accuracy posterior distribution of the orbit elements using surveys of three or more short circuit arcs of the same artificial celestial body.SELECTED DRAWING: Figure 5

Description

[Technical field]

[0001] The present invention relates to an orbit estimation device, an orbit estimation method, and an orbit estimation program for an artificial celestial body. [Background technology]

[0002] Currently, a large amount of space debris is distributed around the Earth. As shown in Figure 1, when operating an artificial satellite 1, in order to avoid these space debris 2, it is necessary to estimate the orbit of the space debris 2 and temporarily change the artificial satellite 1 from its original orbit 3 to, for example, orbits 4 and 5 so as to avoid the forecast circle 6. Temporary orbital changes cause economic losses due to shortened lifespan due to excess fuel consumption and business suspension during the orbital change period. As shown in Figures 1(a) and 1(b), the smaller the forecast circle 6, the smaller the distance of the temporary orbital change and the shorter the change period, and therefore the smaller the economic loss. By accurately estimating the orbit of an artificial celestial object, the distance of the temporary orbital change can be reduced.

[0003] Conventionally, when estimating the orbit of an artificial celestial body such as an artificial satellite or space debris, Keplerian elements (orbital elements based on Kepler's law) have been estimated from optical observations or observations by radio radar. Optical observation has advantages such as higher positioning accuracy compared to radio radar, and easier observation of artificial celestial bodies in high orbits such as beyond geostationary orbit. On the other hand, distance information cannot be obtained, and only two angles, right ascension (Right Ascension; RA) and declination (Declination; Dec), can be obtained. In orbit determination based on observations that only obtain the angle of the motion of an artificial celestial body, methods based on celestial mechanics such as the Gauss method and Goddings method are generally used. In recent years, methods that can be applied even when orbit determination is difficult using these methods have been developed (for example, see Patent Document 1). Patent Document 1 discloses a method for initial startup determination by observing only angles.

[0004] On the other hand, the use of astronomy observation data such as wide-field video survey observations has been attracting attention as a means of discovering and observing artificial celestial bodies. Wide-field video survey observations have traditionally focused on sudden phenomena such as supernova explosions, detection of optical counterparts of gravitational waves, and moving celestial bodies such as planets, comets, and near-Earth objects. Wide-field video survey observations observe a wide range of the sky at night from a designated observation position on the ground by taking short videos, and detect a large number of artificial celestial bodies. In this observation, by observing a wide range with short videos every night, several thousand moving celestial bodies are detected in one night, and the majority of these detected moving celestial bodies are thought to be artificial celestial bodies. If such observation data can be used to discover artificial celestial bodies and determine their orbits, it can be useful for the safe operation of artificial celestial bodies. [Prior art documents] [Patent documents]

[0005] [Patent Document 1] Patent No. 6106341 specification Summary of the Invention [Problem to be solved by the invention]

[0006] However, with conventional methods for estimating Keplerian elements from optical observations or radio radar observations, it is difficult to perform robust orbit determination from a short arc or when observation points are far apart, and it is also difficult to obtain orbital elements that reflect perturbations. In sweeping observations such as wide-field video surveys, the time for each video observation is short, and therefore the arc of an artificial celestial body is much smaller than one degree, making it difficult to apply conventional methods to orbit determination.

[0007] In view of the above problems, the present invention aims to provide an orbit estimation device, an orbit estimation method, and an orbit estimation program for an artificial celestial body that can estimate orbital elements with high accuracy even from observations of a short arc in a sky survey. [Means for solving the problem]

[0008] A first aspect of the present invention is an orbit estimation device for multiple artificial celestial bodies that performs orbit estimation using Bayesian estimation with MCMC based on sky surveying observations, and includes a memory unit that stores data for each of the multiple artificial celestial bodies obtained by the sky surveying observations, and a calculation unit including: a low-precision posterior distribution calculation unit that obtains a low-precision posterior distribution of orbital elements and celestial body sizes from short arcs of the orbits of each of the multiple artificial celestial bodies obtained by the sky surveying observations; a man-made celestial body identification unit that identifies observations of the same artificial celestial body by comparing posterior distributions for different observations among the sky surveying observations; and a high-precision posterior distribution calculation unit that obtains a high-precision posterior distribution of orbital elements using observations of three or more short arcs of the same artificial celestial body.

[0009] In the first aspect of the present invention, the low-precision posterior distribution calculation unit may calculate orbital elements from the times and positions of each of the multiple artificial celestial objects obtained by sky surveying, and sample the calculated orbital elements.

[0010] In a first aspect of the present invention, when calculating the orbital elements, the low-precision posterior distribution calculation unit may calculate an average position and an average velocity at the mean time as tangential elements from the time and position of each of a plurality of artificial celestial bodies, and calculate the orbital elements from the calculated average position and average velocity.

[0011] In the first aspect of the present invention, the low-precision posterior distribution calculation unit may obtain a low-precision posterior distribution of celestial object sizes from the brightnesses of the plurality of artificial celestial objects.

[0012] In a first aspect of the present invention, the artificial celestial object identification unit may select any combination of observations from three fields of view, compare low-precision posterior distributions for the selected combinations, calculate the area of ​​the common part of the combinations of observations from the three fields of view, and if the calculated area exceeds a predetermined threshold, consider the selected combinations to be those of the same artificial celestial object.

[0013] In the first aspect of the present invention, the high-accuracy posterior distribution calculation unit may sample the equinoctial elements.

[0014] In the first aspect of the present invention, the high-precision posterior distribution calculation unit may use MultiNest when sampling the equinoctial elements.

[0015] In the first aspect of the present invention, when the high-precision posterior distribution calculation unit executes sampling of the equinoctial elements using MultiNest, the high-precision posterior distribution calculation unit may obtain the prior by obtaining each observation point of a combination of observations of three fields of view.

[0016] In the first aspect of the present invention, the high-precision posterior distribution calculation unit may obtain the prior by sampling tangential elements of each of the obtained observation points for each combination of observations of three fields of view, and then summing up the posterior distributions of the equinoctial elements.

[0017] In the first aspect of the present invention, when the high-precision posterior distribution calculation unit performs sampling of equinoctial elements using MultiNest, the prior may be obtained from the distribution of each parameter obtained from the area or volume of the common part of the equinoctial elements obtained in matching of multiple fields of view using a coarse posterior distribution.

[0018] In the first aspect of the present invention, the high-precision posterior distribution calculation unit may repeat sampling for all combinations of observations of three fields of view when sampling the equinoctial elements.

[0019] A second aspect of the present invention is an orbit estimation method for multiple artificial celestial objects, which performs orbit estimation by Bayesian estimation using MCMC based on sky sweeping observations, and comprises a low-precision posterior distribution calculation step of obtaining a low-precision posterior distribution of orbital elements and celestial object sizes from short arcs of the orbits of each of the multiple artificial celestial objects obtained by the sky sweeping observations, a man-made celestial object identification step of identifying observations of the same artificial celestial object by comparing posterior distributions for different observations among the sky sweeping observations, and a high-precision posterior distribution calculation step of obtaining a high-precision posterior distribution of orbital elements using observations of three or more short arcs of the same artificial celestial object.

[0020] A third aspect of the present invention is an orbit estimation program comprising a memory unit for storing data relating to each of a plurality of artificial celestial bodies obtained by sky sweeping observations, and a calculation unit, and in an orbit estimation device for a plurality of artificial celestial bodies that performs orbit estimation by Bayesian estimation using MCMC based on the sky sweeping observations, the orbit estimation device is provided with a low-precision posterior distribution calculation function for obtaining a low-precision posterior distribution of orbital elements and celestial body sizes from short arcs of the orbits of each of the plurality of artificial celestial bodies obtained by sky sweeping observations, an artificial celestial body identification function for identifying observations of the same artificial celestial body by comparing posterior distributions for different observations among the observations made by the sky sweeping observations, and a high-precision posterior distribution calculation function for obtaining a high-precision posterior distribution of orbital elements using observations of three or more short arcs of the same artificial celestial body. Effect of the Invention

[0021] According to the present invention, it is possible to provide an orbit estimation device, an orbit estimation method, and an orbit estimation program for an artificial celestial body that are capable of estimating orbital elements with high accuracy even from observations of a short arc in a sky sweep. [Brief description of the drawings]

[0022] [Figure 1] FIG. 1 is a schematic diagram illustrating how an artificial satellite's orbit is changed due to space debris. [Diagram 2]These are schematic diagrams explaining the orbit of an artificial celestial body described by Keplerian elements. Figure 2(a) shows the relationship between the positions of the Earth and the artificial celestial body on the elliptical orbital plane and a and e. Figure 2(b) shows the relationship between the position of the elliptical orbit in three dimensions and ω, I, and Ω. Figure 2(c) shows the relationship between the position of the artificial celestial body on the elliptical orbital plane at t and M, true anomaly ν, and eccentric anomaly E. [Diagram 3] FIG. 2 is a schematic diagram illustrating the observation position of the artificial celestial body. [Figure 4] 1 is a block diagram showing an example of a trajectory estimation device according to an embodiment of the present invention. [Diagram 5] 1 is a flowchart illustrating an overall flow of an example of a trajectory estimation method according to the present embodiment. [Figure 6] 10 is a flowchart illustrating a method for calculating a low-precision posterior distribution in the trajectory estimation method according to the present embodiment. [Figure 7] 11 is a graph showing an example of the relationship between αi, δi and average positions αavg, δavg, and average velocities vα,avg, vδ,avg, calculated in the trajectory estimation device according to the embodiment of the present invention. [Figure 8] 4 shows an example of a calculation result of Log Likelihood in the trajectory estimation device according to the embodiment of the present invention. [Figure 9] 9 is a graph showing an example of actual observations, where FIG. 9(a) shows the trace in burn-in, FIG. 9(b) shows the likelihood, and FIG. 9(c) shows the adoption probability. [Figure 10] 13 is a graph illustrating an example of the posterior distribution of sampled tangent elements. [Figure 11] This is a graph of the posterior distribution of Keplerian elements. [Figure 12] This is a graph of the posterior distribution of equinoctial elements. [Figure 13] For sampling with two parameters, Figure 13(a) shows the trace in burn-in, Figure 13(b) shows the likelihood, and Figure 13(c) shows the acceptance probability. [Figure 14]13 is a graph illustrating an example of the posterior distribution of sampled tangent elements. [Figure 15] Graph of the posterior distribution of Keplerian elements [Figure 16] Graph of the posterior distribution of equinoctial elements [Figure 17] For sampling with four parameters, Figure 17(a) shows the trace in burn-in, Figure 17(b) shows the likelihood, and Figure 17(c) shows the acceptance probability. [Figure 18] 13 is a graph illustrating an example of the posterior distribution of sampled tangent elements. [Figure 19] Graph of the posterior distribution of Keplerian elements [Figure 20] Graph of the posterior distribution of equinoctial elements [Figure 21] 1 is a schematic diagram illustrating the relationship between an artificial celestial body, an observation device, and the angle of sunlight. [Figure 22] 1 is a graph showing a prior distribution of albedo. [Figure 23] 13 is a graph showing the calculated posterior distribution. [Figure 24] This is a graph showing the posterior distribution when all three observations are of the same celestial object. [Diagram 25] This is a graph showing the posterior distribution when one field of view has different celestial objects. [Figure 26] 13 is a graph showing an example of creating the prior of equinoctial elements using three observation points. [Figure 27] 13 is a graph showing an example of the posterior distribution of equinoctial elements sampled by MultiNest. [Figure 28] 1 is a graph showing an example of the posterior distribution of Keplerian elements sampled by MultiNest. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0023] Next, an embodiment of the present invention will be described with reference to the drawings. In the description of the drawings relating to the embodiment, the same or similar parts are denoted by the same or similar reference numerals. However, the drawings are schematic.

[0024] Furthermore, the embodiments are merely examples of devices and methods for embodying the technical idea of ​​the present invention, and the technical idea of ​​the present invention does not specify the configuration, arrangement, layout, etc. of each component to those described below. The technical idea of ​​the present invention can be modified in various ways within the technical scope defined by the claims.

[0025] (Embodiment) The orbit estimation device according to this embodiment can estimate orbital elements with high accuracy from data on short circular arcs of the orbits of multiple artificial celestial bodies obtained through sky surveying.

[0026] The trajectory estimation method according to this embodiment can be realized by causing an trajectory estimation device to execute an trajectory estimation program in which processing contents for realizing the trajectory estimation method according to this embodiment are described.

[0027] First, the Keplerian elements and equinoctial elements (celestial equatorial elements) that describe the orbit of the artificial celestial body will be described below. The Keplerian elements are parameters at epoch t, ​​a: semimajor axis, e: eccentricity, ω: argument of anomaly, I: orbital inclination, Ω: right ascension of ascending node, and M: mean anomaly. Figure 2 shows the orbit of the artificial celestial body 20 described by these parameters. Figure 2(a) is a diagram showing the relationship between the positions of the Earth 21 and the artificial celestial body 20 on the elliptical orbital plane and a and e. Figure 2(b) is a diagram showing the relationship between the position of the elliptical orbit in three dimensions and ω, I, and Ω, together with the perigee 22, apogee 23, orbital plane 24, equatorial plane 25, and orbital direction 26 of the artificial celestial body 20. Figure 2(c) is a diagram showing the relationship between the position of the artificial celestial body at t on the elliptical orbital plane and M, true anomaly ν, and eccentric anomaly E.

[0028] Keplerian elements determine the elliptical orbit assuming a two-body problem involving the Earth and an artificial celestial body. On the other hand, it is difficult to express an orbit that is close to a perfect circle or an orbit that is close to parallel to the equatorial plane using Keplerian elements. Therefore, equinoctial elements are also used, which can express the orbit of any artificial celestial body orbiting the Earth, including these orbits. Equinoctial elements are determined by the parameters a: semi-major axis radius, ex=ecos(ω+Ω), e y = esin(ω+Ω), h x = tan(i / 2)cos(Ω), h y = tan(i / 2)sin(Ω), I ν =M+ω+Ω.

[0029] In reality, the orbits of artificial celestial bodies change due to perturbations caused by the unevenness of the Earth's gravity, the gravity of the Moon and the Sun, the radiation pressure of the Sun, the air resistance of the upper atmosphere, etc. When determining orbital elements from optical or radio observations, Keplerian elements that take perturbations into account can be obtained by using orbit propagation algorithms such as SDP4 and SGP4, which handle these perturbations.

[0030] In this embodiment, the orbital elements are estimated from data from multiple nights of sky sweeping observation. Sky sweeping observation will be described below.

[0031] One example of a sky sweep is wide-field video survey observation. Wide-field video survey observation is an observation method in which a given range of field of view is captured in monochromatic visible light using multiple image sensors installed on a telescope. In a specific example, a field of view of about 20 square degrees is covered by 84 CMOS sensors. Two frames per second are taken per field of view, with 12 frames of 0.5-second exposure. In other words, a six-second video is obtained per field of view. By moving the telescope slightly around one field of view and repeatedly taking images, the gaps between the image sensors are filled. Furthermore, by changing the angle of the telescope and repeatedly taking videos, the entire sky is swept at least once per night.

[0032] Moving objects are detected within each image sensor in one field of view from the video data acquired by the wide-field video survey observation. During this detection, since the video data contains a lot of noise and "seeing dance" (fluctuations in star images due to atmospheric turbulence), a random forest classifier is used to extract objects with a certain probability or higher that they are truly moving objects. Thousands of moving objects are detected in one night, and the majority of them are thought to be known or unknown man-made objects.

[0033] The typical apparent speed of a Near Earth Object (NEO) is ~3 arcsec [s -1 ]. Apparent speed is ~50 arcsec [s -1 ] corresponds to an altitude of about 10,000 km. The video data from the wide-field video survey captures objects moving faster than this speed, but artificial objects moving slower than this speed are detected as contamination of NEOs, such as space debris and artificial satellites.

[0034] The arc of a detected artificial object is much smaller than 1 degree in a 6 second observation over one field of view; for example, at an altitude of 10,000 km, the 6 second arc is 1 / 12 degree.

[0035] The data of artificial objects obtained by the wide-field video survey are the observation position, i.e., angle and brightness, of a given frame per field of view. The observation position of an artificial object is determined by the imaging time t i RA(α i ), and Dec(δ i ), brightness is visible light class m i Here, i is an integer between 1 and the number of frames, for example, an integer between 1 and a maximum of 12 in the case of 12 frame observation.

[0036] FIG. 3 shows the telescope 33 and the image capture time t 0 , t 1 , and t 2 The artificial celestial body 30 moves along an orbit 34. 0The observation position of artificial celestial body 30 in RA(α 0 ), and Dec(δ 0 ) is expressed as the image capture time t 1 , and t 2 The same applies to.

[0037] The orbit estimation device according to this embodiment will be described with reference to Fig. 4. As shown in Fig. 4, the orbit estimation device 40 according to this embodiment is composed of a CPU 41 for executing various calculations, a storage unit 42 for storing data including processing programs and data on artificial celestial bodies, an I / O (input-output interface) 43, a display unit 44, and an input unit 45.

[0038] The trajectory estimation device 40 may be various electronic computers (computing resources) such as a personal computer (PC), a mainframe, a workstation, a cloud computing system, etc. The I / O 43 is an interface, a buffer, etc. for communication (transmission and reception).

[0039] The display unit 44 is a display device such as a display, and displays the results of calculations by the CPU 41. The input unit 45 is an input device such as a keyboard, a mouse, etc., and is realized by a device that accepts input from a user and transmits it to the CPU 41.

[0040] The block diagram of Fig. 4 shows the functional units in the CPU 41. When each functional unit of the CPU 41 is realized by software, the CPU 41 realizes it by executing the instructions of a program, which is the software that realizes each function. In detail, the CPU 41 includes a low-precision posterior distribution calculation unit 411, an artificial celestial object identification unit 412, a high-precision posterior distribution calculation unit 413, etc.

[0041] The low-precision posterior distribution calculation unit 411 calculates low-precision posterior distributions of orbital elements and celestial object sizes from one field of view of the sky survey, based on observation data from the sky survey acquired by the I / O 43 .

[0042] The artificial celestial object identification unit 412 uses a rough posterior distribution obtained by MCMC sampling from observations of one field of view to identify observations of the same artificial celestial object by comparing the posterior distributions for different observations, i.e., performing matching.

[0043] The high-precision posterior distribution calculation unit 413 calculates high-precision posterior distributions of orbital elements from observations of multiple fields of view in the sky survey.

[0044] An orbit estimation method according to an embodiment of the present invention will be described below. The orbit estimation method according to this embodiment is a method for identifying an artificial celestial object and determining its orbit when optical observations of a short arc obtained over multiple nights in a sky survey including a wide-field video survey are obtained. Figure 5 shows a flowchart explaining the overall flow of an example of the orbit estimation method according to this embodiment. As shown in Figure 5, the orbit estimation method according to this embodiment is carried out in the following four steps.

[0045] In step S501, observation data from a sky survey is acquired.

[0046] In step S502, a low-precision posterior distribution of the orbital elements and celestial object size is calculated from one field of view of the sky survey. In this embodiment, the posterior distribution of the orbital elements is obtained by MCMC sampling of the tangential orbital elements defined below. Furthermore, the posterior distribution of the celestial object size is calculated from the posterior distribution of the orbital elements and the observed brightness.

[0047] In step S503, a rough posterior distribution obtained by MCMC sampling from observations of one field of view is used to identify observations of the same artificial celestial body by comparing the posterior distributions for different observations, i.e., matching is performed.

[0048] In step S504, a highly accurate posterior distribution of the orbital elements is calculated from the observations of multiple fields of view in the sky survey.

[0049] The steps of the trajectory estimation method according to the present embodiment shown in FIG. 5 will now be described in detail with reference to FIGS.

[0050] FIG. 6 shows details of step S502 shown in FIG.

[0051] (Calculate tangent elements) In step S601, osculating parameters are calculated. First, the time t i , and position α i , δ i Based on the following equation, the average time t avg The average position α avg , δ avg and the average speed v α,avg , v δ,avg Calculate the average time t avg , and the average position α avg , δ avg is given by the following formula:

[0052]

number

[0053]

number

[0054]

number

[0055] In addition, α i , δ i Therefore, the velocity at each time v α,i , v δ,i but

[0056]

number

[0057]

number

[0058] and the average velocity v α,avg , v δ,avg teeth,

[0059]

number

[0060]

number

[0061] is given by:

[0062] The tangent element is the average time t avg The distance r from the center of the earth, the time derivative dr / dt of the distance r from the center of the earth, and the average position α avg , δ avg , average speed v α,avg , v δ,avg Note that r and dr / dt are defined in the geocentric frame, and the other four parameters are defined in the topocentric frame.

[0063] In Fig. 7, α i , δ i and the average position α avg , δ avg , average speed v α,avg , v δ,avg In the graph shown in FIG. 7, data 71 shown in gray scale in the graph is α i , δ i The grayscale intensity changes at time t i The horizontal and vertical axes correspond to α i , δ i The grayscale indicates the difference from the average time. The data 72 in the center of the graph is the average position α avg , δ avgAnd the vector from data 72 to data 73 is (v α,avg ,v δ,avg ) multiplied by 0.25 seconds.

[0064] (Procedure for obtaining orbital elements from tangential elements) In step S602, the orbital elements (Keplerian Elements) are composed of six parameters like the tangential elements, and since each parameter has a variable transformation relationship with each other, the orbital elements are calculated based on the obtained tangential elements. The procedure is explained below.

[0065] Calculate the position of a celestial object relative to the center of the Earth. Distance from the Earth's center r, RAα at the observer's center Geo , Decδ Geo Then, the vector from the center of the Earth to the observer is

[0066]

number

[0067] and the vector from the observer to the object is

[0068]

number

[0069] and the vector from the center of the Earth to the celestial body is

[0070]

number

[0071] Now we can solve the vector S for d to determine the vector s, which gives us S, the position of the celestial body in the geocentric coordinate system.

[0072] Tangent elements r, dr / dt, α avg , δ avg , vα,avg , v δ,avg If we give the position of a celestial body in the geocentric coordinate system, we can calculate the velocity of the celestial body from these small displacements. That is, the position vector in geocentric coordinates

[0073]

number

[0074] and the velocity vector

[0075]

number

[0076] Once the position vector pGeo and the velocity vector vGeo are obtained, the angular momentum vectors can be used to find the orbital elements in the Earth-celestial body two-body problem.

[0077] (Definition of Log Likelihood in MCMC) In step S603, the Log Likelihood is defined when sampling the tangential orbital elements using MCMC (Markov Chain Monte Carlo). Note that Log Likelihood is maximum likelihood estimation, which is a method in statistics for point estimating the parameters of the probability distribution that given data follows. Also, MCMC is a general term for various algorithms that sample probability distributions by creating a Markov chain that has the desired probability distribution as an equilibrium distribution.

[0078] First, for each sample, the observation time t i RA(α m,i ), Decδ m,i Calculate the observation α i , δ i are the α of the model, respectively. m,i , δ m,iWe assume that the differences are normally distributed around the range of the observation's positioning accuracy.

[0079]

number

[0080]

number

[0081] Note that the difference in the RA direction must be multiplied by COS(δ).

[0082] Furthermore, if the observation error is σ, the Log Likelihood is

[0083]

number

[0084] where f(x,σ) has mean 0 and variance σ 2 gives the natural logarithm of the probability density function at x of the normal distribution.

[0085] The orbital elements obtained from the tangential orbital elements are free of perturbations, while from these orbital elements α m,i , δ m,i When calculating, perturbations can be taken into account and algorithms such as SGP4 can be used. The difference between taking perturbations into account and not taking them into account appears as a global position offset in the range of the time scale considered here. Therefore, if the above formula is applied as is, LL becomes very small due to the global offset. Therefore, a term that absorbs the global offset is added to common sense.

[0086] First, Δα i , Δδ i is divided into global and local components. glob , Δδglob is averaged over i, and the following formula is used:

[0087]

number

[0088]

number

[0089] Calculate based on the local component Δα loc,i , Δδ loc,i is the displacement from the global component, and is expressed as follows:

[0090]

number

[0091]

number

[0092] The global component is the typical difference σ glob Assuming that the distribution follows a normal distribution with standard deviation of , and the local component follows a normal distribution with standard deviation of the observation positioning accuracy, the Log Likelihood is

[0093]

number

[0094] The third and fourth terms on the right-hand side can be interpreted as terms that absorb the global offset caused by whether or not perturbations are taken into account.

[0095] As an example of actual observation, the position determination accuracy of the observation is ~0.3 arcsec. Also, σ globThe typical value of is estimated to be about 20 arcmin. In the following, when calculating the position of an artificial celestial body from the orbital elements, perturbations are taken into account using SGP4.

[0096] (Setting the initial values ​​for MCMC) In step S604, the initial value of MCMC is set. The initial value of MCMC in sampling of tangent elements is the average position α avg , δ avg , average speed v α,avg , v δ,avg For , the values ​​obtained from observations can be used as is. On the other hand, for the tangent elements r and dr / dt, initial values ​​must be given externally.

[0097] The initial values ​​of the tangent elements r and dr / dt are calculated at several points on these grids, and are chosen from those that give a closed orbit. For artificial celestial bodies, the range of these elements is limited. For example, the tangent element r ranges from an altitude of about 100 km to an altitude of about 380,000 km for the lunar orbit, and dr / dt ranges from ±5 kms. -1 This is the realistic range. After cutting an appropriate number of grids (for example, about 20) within this range, the orbital elements are calculated using values ​​obtained from observations for other tangent elements. The points that give a closed orbit or points with a large LL can be used as the initial values ​​for MCMC.

[0098] The advantage of realizing the posterior distribution of orbital elements by sampling tangential elements is that, in observations, LL is a convex function with respect to the tangential elements. Since LL is a convex function, robust sampling is possible using MCMC. Figure 8 shows an example of LL calculation at points on a grid of r and dr / dt. In the graph of Figure 8, the horizontal axis is r, the vertical axis is dr / dt, the dots are the points where calculations were performed, the grayscale density indicates the value of LL, and the area 80 indicated by dark gray in the graph is the area with a large LL.

[0099] (Setting the prior distribution) In step S605, a prior distribution can be set when performing MCMC sampling, which will be described later. Bayesian estimation starts with determining a prior distribution. When an approximate prior distribution has been obtained from past literature, preliminary data, experience, or the like, such a distribution can be used as the current prior distribution. The prior distribution refers to information that has been held in advance, and the likelihood is acquired information.

[0100] First, in the sampling of the tangent elements, the two-dimensional distribution of the tangent elements r and dr / dt, and the average position α avg , δ avg , average speed v α,avg , v δ,avg Each one-dimensional distribution can be used. However, the influence of the prior distribution on the posterior distribution is small. The prior distributions of the tangent elements r and dr / dt can be calculated from the distribution of the orbital elements of known artificial celestial bodies. The mean position α avg , δ avg , average speed v α,avg , v δ,avg For example, a normal distribution can be used using the standard deviation obtained from observations of 12 frames.

[0101] (Tangential element sampling) In step S606, when sampling tangent elements using MCMC, sampling is performed in each of the following three patterns (1) to (3). (1) Sampling with one parameter of the tangent element r (2) Sampling with two parameters of tangent element r and dr / dt (3) Tangent element r, dr / dt, average speed v α,avg , v δ,avg Sampling with four parameters In (1), we assume that dr / dt=0 and the average position α avg , δ avg , average speed v α,avg , v δ,avg is fixed as the observed value. In (2), the unsampled mean position α avg , δ avg , average speed vα,avg , v δ,avg is fixed as the observed value. In (3), the unsampled mean position α avg , δ avg is fixed to the observed value. (2) and (3) are suitable for sampling various cases, including orbits with large eccentricities, but they tend to stochastically overestimate the eccentricity when the orbit of the target object is close to a circular orbit. (1) also cannot completely eliminate the tendency to stochastically overestimate the eccentricity when the orbit of the target object is close to a circular orbit, but since dr / dt = 0, it can give values ​​close to the true orbital elements when the tangent element r is close to the true value.

[0102] (1) Sampling with one parameter of the tangent element r MCMC sampling is performed on one parameter, r, of the tangent elements. As mentioned above, Metropolis-Hastings (MH) sampling or Gibbs sampling is used from the obtained initial value, and sampling is performed from the Log Likelihood formula (Equation 19). The number of steps required to converge to the posterior distribution is less than 100, and 200 steps are sufficient for burn-in. In addition, calculation time can be reduced by setting multiple initial values ​​and sampling in parallel starting from each initial value.

[0103] The width of the proposal distribution can be set so that the adoption probability in burn-in falls within an appropriate range. As an example of actual observation, in the case of 6-second 12-frame observation data, it was confirmed that stable convergence was achieved by setting the width of the proposal distribution so that the adoption probability was about 20 to 40%. Figures 9(a), 9(b), and 9(c) show examples of trace, likelihood, and adoption probability in burn-in. The likelihood shown in Figure 9(b) is normalized so that the maximum value is 1. The grayscale of each data in Figures 9(a), 9(b), and 9(c) indicates the chains in parallel sampling, i.e., chains that started from different initial values. Around the 30th step, the influence of the initial values ​​disappears, and the decrease in the adoption probability stops, as shown in Figure 9(c).

[0104] Figure 10 shows an example of the posterior distribution of sampled tangent elements. The graph surrounded by dashed line A is the one-dimensional posterior distribution of each parameter, and the graph surrounded by dashed line B is the two-dimensional posterior distribution of each parameter pair. The graph surrounded by dashed line C shows the Log Likelihood in two dimensions. The points in the graphs surrounded by dashed lines B and C indicate the initial values. Since this is a sampling of one parameter, all parameters except r are fixed. The accuracy of determining r is tens of thousands of km.

[0105] Once the tangential elements are sampled, the corresponding orbital elements can be obtained. Figure 11 shows the graph of the posterior distribution of the Keplerian elements, and Figure 12 shows the graph of the posterior distribution of the equinoctial elements. As in Figure 10, the graph surrounded by dashed line A shows the one-dimensional posterior distribution of each parameter, the graph surrounded by dashed line B shows the two-dimensional posterior distribution of each parameter pair, and the graph surrounded by dashed line C shows the two-dimensional Log Likelihood. The dotted lines in the figure show the values ​​of known artificial celestial bodies matched by position and velocity. In the tangential element space, all parameters except r have a value of 1, but in the orbital element space, the orbital elements corresponding to each tangential element set have a spread distribution.

[0106] The celestial objects shown in Figures 10 to 12 match the positions and velocities of known artificial celestial objects, and their eccentricities are close to 0, at e = 0.00433. Because the orbit is close to a circular one, the argument of perigee ω is indefinite, and therefore the mean anomaly M is also indefinite. In Figure 11, the distribution of these parameters does not cover the values ​​of the matched known celestial object, i.e., values ​​close to the true values. This is because the orbit of the celestial object is close to a circular orbit. In fact, as shown in Figure 12, in the posterior distribution of the equinoctial elements, the posterior distribution of all parameters covers the values ​​of the known celestial object.

[0107] (2) Sampling with two parameters: r and dr / dt Sampling with two parameters, r and dr / dt, is performed in the same manner as sampling with one parameter. An example of the results of sampling with two parameters, r and dr / dt, is shown in Figures 13 to 16. The details of Figures 13 to 16 are the same as the details of sampling with one parameter in Figures 9 to 12.

[0108] (3) r, dr / dt, average speed v α,avg , v δ,avg Sampling with four parameters r, dr / dt, average speed v α,avg , v δ,avg The sampling for four parameters is performed in the same way as for one parameter. r, dr / dt, average velocity v α,avg , v δ,avg An example of the results of sampling with these four parameters is shown in Figures 17 to 20. The details of Figures 17 to 20 are the same as the details of Figures 9 to 12 in the case of sampling with one parameter.

[0109] (Calculation of the posterior distribution of object sizes) In step S607, after sampling of the tangential elements, the posterior distribution of the size of the object is calculated from the brightness of the object. The brightness of the artificial object depends on conditions such as the distance to the artificial object, the size, albedo, shape, and rotation of the artificial object. If the albedo of the artificial object is ρ, the cross-sectional area is A, the distance from the observer to the artificial object is d, and the function representing the directivity of reflection is F(ψ), then the observed magnitude of the artificial object is:

[0110]

number

[0111] Here, m is the apparent magnitude of the sun, which is -26.47 in visible light (SDSSg band). Regarding the directivity of reflection from an artificial celestial body, if we assume that the object scatters isotropically only in the direction in which the sunlight hits it in a spherical model, then F(ψ) = 1 / 2π. Furthermore, if the typical size of an artificial celestial body is s, then A = πs 2 As shown in Fig. 21, when an artificial celestial body 210 is spherical, if the angle between a line of sight 212 to the artificial celestial body as seen from an observation device 211 and a line of sight 213 to the sun is β, then the angle (phase angle) γ of the part of the artificial celestial body 210 illuminated by sunlight 214 is equal to β. Therefore, the area of ​​the artificial celestial body illuminated by the sun is

[0112]

number

[0113] In the above formula, A is shine Substituting:

[0114]

number

[0115]

number

[0116] Therefore, the size of the artificial object is

[0117]

number

[0118] It is expressed as:

[0119] From the above formula, if the distance d and albedo ρ are determined, the size of the celestial body can be calculated from the observed magnitude m. If the tangential elements, i.e., orbital elements, are determined, the distance d can be calculated by performing orbital calculations. The albedo ρ is assumed to be a prior distribution. The prior distribution of the albedo ρ is obtained using the brightness of objects observed in actual observations that match known celestial bodies. The brightness fluctuations of the same celestial body, other than the distance and phase angle, are due to the rotation and attitude changes of the celestial body, etc., but these are all pushed onto the albedo term and the above formula is used to obtain the fluctuation range of ρ. By setting the average value of ρ to a typical value of 0.1, the distribution of ρ is obtained as shown in Figure 22. This distribution of ρ can be used as the prior distribution as it is, or the distribution of ρ can be applied to a normal distribution as the prior distribution. The distribution of ρ applied to a normal distribution is shown as a dotted line in the graph in Figure 22. In the following, the distribution of ρ applied to a normal distribution is referred to as the prior distribution.

[0120] The posterior distribution of object sizes is obtained from the observed magnitude, the sampled tangent elements, and the prior distribution of albedo. The observed magnitude is the average m avg As

[0121]

number

[0122] and calculate using equation (24).

[0123] The obtained posterior distributions are shown in Figure 23. The posterior distributions of object sizes are obtained for the one-parameter, two-parameter, and four-parameter cases, as well as the tangent element.

[0124] (Multiple-view matching using coarse posterior distribution) In step S503 shown in FIG. 5, a posterior distribution obtained by MCMC from the observation of one field of view is used to match multiple fields of view. First, a combination of observations of any three fields of view is selected. For this combination, the posterior distributions obtained by MCMC for one, two, and four parameters are compared. For the posterior distribution of seven parameters, which is the equinoctial elements plus size, the area or volume of the overlapping portion of the observations of the three fields of view, that is, the common portion of the posterior distributions, is calculated, and the area or volume of the common portion obtained is called the degree of agreement. The area of ​​the common portion of the posterior distributions is the area of ​​the common portion of the one-dimensional posterior distributions of the seven parameters. Also, the volume of the common portion of the posterior distributions is the volume of the common portion of the posterior distributions of the seven parameters in a seven-dimensional space. Either the area or the volume of the common portion may be used as the degree of agreement.

[0125] The degree of agreement of the posterior distribution is large if the observations in the three fields of view are all of the same celestial object, and small otherwise. If the degree of agreement exceeds a certain threshold, we proceed to obtain a high-precision posterior distribution. Among the seven parameters, the h x , h y In the case of the three fields of view, there tends to be a large difference between when all three observations are of the same celestial object and when they are not. Figure 24 shows examples of the posterior distribution and common parts when all three observations are of the same celestial object, and Figure 25 shows examples of the posterior distribution and common parts when one field of view is different. Figures 24(a) to 24(f) and 25(a) to 25(f) show examples of the posterior distribution of equinoctial elements, and Figures 24(g) and 25(g) show examples of the posterior distribution of seven parameters including size. The histograms in each graph are probability density distributions represented by the vertical axis on the left, and the intensity of the grayscale color indicates different fields of view. The common parts are shown by dots and are represented by the vertical axis on the right. The vertical dotted lines show the values ​​of the matched known objects.

[0126] (Obtaining highly accurate posterior distributions) In step S504, the equinoctial elements are sampled using the observed values ​​of RA and Dec. This sampling is performed using MultiNest, a type of nested sampling. Since the log likelihood is not a convex function, there is no advantage to sampling the tangent elements. Therefore, here, the equinoctial elements themselves are sampled.

[0127] MultiNest has the advantage that the convergence to the posterior distribution is relatively robust even when the log likelihood is non-convex or has many peaks. Taking advantage of this advantage, equinoctial elements sampling can be performed relatively robustly. Dealing with multi-peaks is essential for robust orbit estimation. This is because, for example, when the log likelihood is reasonably accurate in all views and when it is good only in a certain view (12 frames), the log likelihood becomes large, i.e., multi-peaks. Furthermore, when the observation times of the views are separated by, for example, about one day, the log likelihood changes significantly when the parameters are changed slightly, i.e., the gap between the multi-peaks is large, and it is difficult to move between the multi-peaks with random walk sampling such as MCMC.

[0128] In addition, by sampling the equinoctial elements themselves, it is possible to use an algorithm that handles perturbations such as SDP4 to calculate the RA and Dec of the model from the equinoctial elements. This allows the effect of perturbations to be reflected in the posterior distribution of the orbital elements that are finally obtained.

[0129] (Prior setting) To use MultiNest, we calculate the priors for the parameters to be sampled. The priors may be obtained from the distribution of each parameter obtained from the area or volume of the intersection of equinoctial elements obtained in matching multiple fields of view using a rough posterior distribution, or from observations of three fields of view.

[0130] When obtaining the priors of equinoctial elements from three observations, the observation α of the jth field of view is j,i , δ j,i From representative point α j,o , δ j,o This can be the first or last observation in each field of view, or the average. At this stage, three observation points (j=1,2,3) are obtained.

[0131] Furthermore, using these three observation points, sampling of the tangent elements may be performed in the same procedure as that in which sampling of the tangent elements is performed by MCMC in each of the following three patterns (1) to (3). (1) Sampling with one parameter of the tangent element r (2) Sampling with two parameters of tangent element r and dr / dt (3) Tangent element r, dr / dt, average speed v α,avg , v δ,avg Sampling with four parameters Here, the reason why the tangential elements can be used is that the observation points in each field of view are narrowed down to one, thereby reducing multi-peak characteristics.

[0132] After sampling the tangent elements, we may add the posterior distribution of the equinoctial elements obtained by sampling 1, 2, and 4 parameters. This posterior distribution (probability density function) is the prior of the equinoctial elements given to MultiNest. In practice, we give the inverse function of the cumulative probability density function (percentile point function, ppf).

[0133] Figure 26 shows an example of creating the priors of equinoctial elements using three observation points. All the graphs on the left side of Figure 26 are probability density functions, and histograms 261 in the graphs are obtained by sampling tangential elements, and step 262 is a numerical correction. All the graphs on the right side of Figure 26 are cumulative density functions calculated from the probability density functions.

[0134] 3-field observation α j,i , δ j,i The equinoctial elements are sampled using all of the above. Here, because SDP4 is also used for orbital propagation, an equation that does not consider offsets due to perturbations is adopted as the log likelihood. The posterior distribution of the equinoctial elements is obtained by sampling with MultiNest. Figure 27 shows an example of the posterior distribution of equinoctial elements sampled with MultiNest. Figure 28 shows an example of the posterior distribution of Keplerian elements sampled with MultiNest. For example, for a, the distribution width is on the order of 10 m, indicating that the posterior distribution is obtained with high accuracy.

[0135] As described above, the present invention naturally includes various embodiments not described here. Therefore, the technical scope of the present invention is defined only by the invention-specifying matters related to the scope of the claims appropriate from the above description. [Explanation of symbols]

[0136] 1 satellite 2. Space debris 3, 4, 5 orbits 6 Forecast Circle 20, 30 Artificial celestial bodies 21 Earth 22 Perigee 23 Apogee 24 Raceway surface 25 Equatorial plane 26 Circulation Direction 33 Telescope 34 orbit 40 Orbit estimation device 41 CPU 42 Storage section 43 I / O (Input / Output Interface) 44 Display section 45 Input section 411 Low precision posterior distribution calculation unit 412 Matching Department 413 High precision posterior distribution calculation unit 71, 72, 73 Data 210 Artificial celestial bodies 211 Observation Equipment 212, 213 line of sight 214 Sunlight

Claims

1. An orbit estimation device for multiple artificial celestial bodies that estimates orbits by Bayesian estimation using MCMC based on sky surveying observations, comprising: a storage unit for storing data relating to each of the plurality of artificial celestial bodies obtained by the sky survey; A calculation unit, a low-precision posterior distribution calculation unit that obtains a low-precision posterior distribution of orbital elements and celestial object sizes from short arcs of the orbits of each of the plurality of artificial celestial objects obtained by the sky sweep; an artificial celestial object identification unit that identifies observations of the same artificial celestial object by comparing the posterior distributions of different observations among the sky surveys; a high-precision posterior distribution calculation unit that obtains a high-precision posterior distribution of orbital elements using observations of three or more short arcs of the same artificial celestial object; A calculation unit including: A trajectory estimation device comprising:

2. The low-precision posterior distribution calculation unit Calculating the orbital elements from the time and position of each of the plurality of artificial celestial bodies obtained by the sky sweeping observation; Sampling the calculated orbital elements 2. The trajectory estimation device according to claim 1 .

3. When calculating the orbital elements, the low-accuracy posterior distribution calculation unit Calculating an average position and an average velocity at an average time as tangential elements from the time and positions of each of the plurality of artificial celestial bodies; Calculating the orbital elements from the calculated average position and average velocity 3. The trajectory estimation device according to claim 2.

4. The low-precision posterior distribution calculation unit A low-precision posterior distribution of the size of the plurality of artificial celestial objects is obtained from the brightness of each of the plurality of artificial celestial objects.

2. The trajectory estimation device according to claim 1 .

5. The artificial celestial body identifying unit includes: Select any combination of three observation fields, For the selected combination, the low-precision posterior distribution is compared, and a degree of agreement between the combinations of observations of the three fields of view is calculated; If the calculated degree of coincidence exceeds a predetermined threshold, the selected combination is regarded as belonging to the same artificial celestial body.

2. The trajectory estimation device according to claim 1 .

6. 6. The trajectory estimation device according to claim 5, wherein the high-precision posterior distribution calculation unit samples equinoctial elements.

7. The trajectory estimation device according to claim 6 , wherein the high-precision posterior distribution calculation unit executes sampling of the equinoctial elements by using MultiNest.

8. 7. The orbit estimation device according to claim 6, wherein when the sampling of the equinoctial elements is performed using MultiNest, the high-precision posterior distribution calculation unit obtains a Prior by obtaining each observation point of the combination of observations of the three fields of view.

9. The high-precision posterior distribution calculation unit obtains the Prior by calculating, after obtaining each observation point of the combination of observations of the three fields of view, Sampling the tangential elements of each of the obtained observation points; Add up the posterior distributions of the equinoctial elements 9. The trajectory estimation device according to claim 8.

10. The trajectory estimation device according to claim 6, characterized in that, when the sampling of the equinoctial elements is performed using MultiNest, the high-precision posterior distribution calculation unit obtains the Prior from the distribution of each parameter obtained from the area or volume of a common part of the equinoctial elements obtained in matching of multiple views using a coarse posterior distribution.

11. 6. The orbit estimation device according to claim 5, wherein the high-precision posterior distribution calculation unit repeats the sampling for all combinations of observations of the three fields of view when sampling equinoctial elements.

12. A method for estimating the orbits of multiple artificial celestial objects, which estimates the orbits of the artificial celestial objects by Bayesian estimation using MCMC based on a sky survey, comprising the steps of: a low-precision posterior distribution calculation step of obtaining a low-precision posterior distribution of orbital elements and celestial object sizes from short arcs of the orbits of each of the plurality of artificial celestial objects obtained by the sky sweeping observation; an artificial celestial object identification step of identifying observations of the same artificial celestial object by comparing the posterior distributions of different observations among the sky surveys; a high-precision posterior distribution calculation step for obtaining a high-precision posterior distribution of orbital elements using three or more short arc observations of the same artificial celestial body; 23. A method comprising:

13. A trajectory estimation device for a plurality of artificial celestial bodies, comprising a memory unit for storing data on each of a plurality of artificial celestial bodies obtained by a celestial survey, and a calculation unit, and performing trajectory estimation by Bayesian estimation using MCMC based on the celestial survey, comprising: a low-precision posterior distribution calculation function for obtaining low-precision posterior distributions of orbital elements and celestial object sizes from short arcs of the orbits of each of the plurality of artificial celestial objects obtained by the sky sweep; an artificial celestial object identification function for identifying observations of the same artificial celestial object by comparing the posterior distributions of different observations among the sky surveys; A high-precision posterior distribution calculation function that obtains a high-precision posterior distribution of orbital elements using three or more short arc observations of the same artificial celestial body; An orbit estimation program to achieve this.

Images