Artificial celestial body orbit estimating device, orbit estimating method, and orbit estimating program

Abstract

Provided are an artificial celestial body orbit estimating device, orbit estimating method and orbit estimating program capable of estimating orbital elements with high precision even from observations of short arcs obtained by means of astronomical surveys. Provided is a method for estimating the orbits of a plurality of artificial celestial bodies by means of Bayesian estimation employing MCMC, on the basis of astronomical surveys, the method comprising: a low-precision posterior distribution calculation step 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 means of astronomical surveys; an artificial celestial body identification step for identifying observations of the same artificial celestial body by comparing the posterior distributions for different observations with one another, among the observations obtained from the astronomical surveys; and a high-precision posterior distribution calculation step for obtaining a high-precision posterior distribution of the orbital elements using the observations of three or more short arcs for the same artificial celestial body.

Description

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

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

[0002] Currently, a large amount of space debris is distributed around the Earth. As shown in Figure 1, when operating an artificial satellite 1, to avoid these pieces of 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 result in economic losses due to a shortened lifespan caused by 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 shorter the distance and duration of the temporary orbital change, and therefore the smaller the economic loss. Accurately estimating the orbit of an artificial celestial object can reduce the distance of the temporary orbital change.

[0003] Traditionally, when estimating the orbits of artificial celestial bodies such as artificial satellites and space debris, Keplerian elements (orbital elements based on Kepler's laws) have been estimated from optical observations or radio radar observations. Optical observations have advantages over radio radar, such as higher positioning accuracy and easier observation of artificial celestial bodies in high orbits, such as those beyond geostationary orbit. However, they do not provide distance information, and only two angles, right ascension (RA) and declination (Dec), are obtained. Orbit determination based on observations that only provide the angle of motion of an artificial celestial body typically uses methods based on celestial mechanics, such as Gauss's method or Goddings' method. In recent years, methods that can be applied even when orbit determination using these methods is difficult have been developed (see, for example, Patent Document 1). Patent Document 1 discloses a method for determining initial launch using angle-only observations.

[0004] Meanwhile, the use of astronomical survey data, such as wide-field video surveys, is gaining attention as a means of discovering and observing artificial celestial bodies. Wide-field video surveys have traditionally focused on transient phenomena such as supernova explosions, the detection of optically identifiable gravitational wave objects, and moving celestial bodies such as planets, comets, and near-Earth objects. Wide-field video surveys observe a wide area of ​​the sky at night from a designated observation location on the ground using short-duration video recordings to detect numerous artificial celestial bodies. By observing a wide area every night using short-duration video recordings, these observations detect thousands of moving celestial bodies per night, and the majority of these detected moving celestial bodies are believed to be artificial. If such observational data can be used to discover artificial celestial bodies and determine their orbits, it could be useful for the safe operation of artificial celestial bodies.

[0005] Patent No. 6106341 specification

[0006] However, conventional methods for estimating Keplerian elements from optical or radio radar observations have difficulty in determining orbits from short arcs 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 duration of each video observation is short, and therefore the arcs traced by artificial celestial bodies are far smaller than one degree, making it difficult to determine orbits using conventional methods.

[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 during a sky sweep.

[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 that may include a memory unit that stores data for each of the multiple artificial celestial bodies obtained by the sky surveying observations, a calculation unit that includes: 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 in 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 time and position 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 orbital elements, the low-precision posterior distribution calculation unit may calculate the average position and average velocity at the average time as tangential elements from the time and position of each of multiple 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 brightness of each of the plurality of artificial celestial objects.

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

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

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

[0015] In the first aspect of the present invention, when sampling equinoctial elements using MultiNest, the high-accuracy posterior distribution calculation unit may obtain the Prior by obtaining the observation points of each 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 at each of the obtained observation points for a 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 intersection of equinoctial elements obtained in matching of multiple fields of view using a coarse posterior distribution.

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

[0019] A second aspect of the present invention is an orbit estimation method for multiple artificial celestial bodies, which performs orbit estimation using Bayesian estimation with MCMC based on sky surveying observations, and comprises: a low-precision posterior distribution calculation step of obtaining low-precision posterior distributions of orbital elements and celestial body sizes from short arcs of the orbits of each of the multiple artificial celestial bodies obtained by sky surveying observations; a man-made celestial body identification step of identifying 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 step of obtaining high-precision posterior distributions of orbital elements using observations of three or more short arcs of the same artificial celestial body.

[0020] A third aspect of the present invention is an orbit estimation program for an orbit estimation device for multiple artificial celestial bodies, comprising a memory unit for storing data on each of multiple artificial celestial bodies obtained by sky sweeping observations, and a calculation unit, and performing orbit estimation using Bayesian estimation with MCMC based on the sky sweeping observations. The orbit estimation device is configured to realize a low-precision posterior distribution calculation function for obtaining a low-precision posterior distribution of orbital elements and celestial body size from short arcs of the orbits of each of the multiple artificial celestial bodies obtained by sky sweeping observations, a man-made celestial body identification function for identifying observations of the same artificial celestial body by comparing posterior distributions for different observations from 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.

[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 can estimate orbital elements with high accuracy even from observations of a short arc in a sky sweep.

[0022] [Correction based on Rule 91 18.12.2024] A schematic diagram explaining how an artificial satellite's orbit changes due to space debris. A schematic diagram explaining the orbit of an artificial celestial body described by Keplerian elements, where FIG. 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, FIG. 2(b) shows the relationship between the position of the elliptical orbit in three dimensions and ω, I, and Ω, and FIG. 2(c) shows the relationship between the position of the artificial celestial body on the elliptical orbital plane at time t and M, true anomaly ν, and eccentric anomaly E. A schematic diagram explaining the observation position of an artificial celestial body. A block diagram showing an example of an orbit estimation device according to an embodiment of the present invention. A flowchart explaining the overall flow of an example of an orbit estimation method according to an embodiment of the present invention. A flowchart explaining a method for calculating a low-precision posterior distribution in the orbit estimation method according to an embodiment of the present invention. α ｉ , δ ｉ and the average position α ａｖｇ , δ ａｖｇ , average speed v α，ａｖｇ , v δ，ａｖｇ9(a) and 9(b) are graphs showing an example of the relationship between the log likelihood and the tangential element. 9(c) are graphs showing an example of an actual observation, where FIG. 9(a) is a trace in burn-in, FIG. 9(b) is the likelihood, and FIG. 9(c) is the acceptance probability. 9(c) are graphs showing an example of the posterior distribution of sampled tangential elements. 9(a) is a graph of the posterior distribution of Keplerian elements. 9(b) is a graph of the posterior distribution of equinoctial elements. 13(a) is a trace in burn-in, FIG. 13(b) is the time derivative of the trace shown in FIG. 13(a), FIG. 13(c) is the likelihood, and FIG. 13(d) is the acceptance probability, in sampling with two parameters. 13(c) is a graph showing an example of the posterior distribution of sampled tangential elements. Graph of the posterior distribution of Keplerian elements Graph of the posterior distribution of equinoctial elements In sampling with four parameters, Figure 17(a) shows the trace in burn-in, Figure 17(b) shows the time derivative of the trace shown in Figure 17(a), and Figure 17(c) shows the average velocity v α，ａｖｇ , and Fig. 17(d) shows the average velocity v δ，ａｖｇ 17(e) is the likelihood, and FIG. 17(f) is the acceptance probability. A graph showing an example of the posterior distribution of sampled tangential elements. A graph of the posterior distribution of Keplerian elements. A graph of the posterior distribution of equinoctial elements. A schematic diagram explaining the relationship between an artificial celestial object, an observation device, and the angle of sunlight. A graph showing the prior distribution of albedo. A graph showing a calculated posterior distribution. A graph showing the posterior distribution when observations of three fields of view are all of the same celestial object. A graph showing the posterior distribution when one field of view is of a different celestial object. A graph showing an example of creating a prior for equinoctial elements using three observation points. A graph showing an example of the posterior distribution of equinoctial elements sampled by MultiNest. 10 is a graph showing an example of the posterior distribution of Keplerian elements sampled by MultiNest.

[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 limit 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) An orbit estimation device according to this embodiment can accurately estimate orbital elements from data on short arcs of the orbits of multiple artificial celestial bodies obtained by sky surveying.

[0026] The trajectory estimation method according to this embodiment can be realized by causing the trajectory estimation device to execute a trajectory estimation program that describes the processing content for realizing the trajectory estimation method according to this embodiment.

[0027] First, the Keplerian elements and equinoctial elements (celestial equatorial elements) that describe the orbit of an artificial celestial body are explained below. Keplerian elements are the parameters at epoch t: a: semimajor axis, e: eccentricity, ω: argument of periapsis, I: orbital inclination, Ω: right ascension of the ascending node, and M: mean anomaly. Figure 2 shows the orbit of an artificial celestial body 20 described by these parameters. Figure 2(a) is a diagram showing the relationship between the positions of 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 Ω, along with perigee 22, apogee 23, orbital plane 24, equatorial plane 25, and orbital direction 26 of the artificial celestial body 20. FIG. 2(c) shows the relationship between the position of the artificial celestial body at time t on the elliptical orbital plane and M, true anomaly ν, and eccentric anomaly E.

[0028] Keplerian elements determine an elliptical orbit that is based on the two-body problem of 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 ｙ =esin(ω+Ω), h ｘ = tan (i / 2) cos (Ω), h ｙ =tan(i / 2)sin(Ω),I ν = M + ω + Ω.

[0029] In reality, the orbit of an artificial celestial body changes 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 the orbital elements from optical or radio wave observations, Keplerian elements that take these perturbations into account can be determined by using orbital propagation algorithms such as SDP4 and SGP4 that handle these perturbations.

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

[0031] One example of a sky survey is wide-field video survey. Wide-field video survey is an observation method in which a predetermined range of field of view is captured in monochromatic visible light using multiple image sensors mounted on a telescope. In a specific example, a field of view of approximately 20 square degrees is covered by 84 CMOS sensors. Two frames per second are captured per field of view, with 12 frames of 0.5-second exposure. In other words, a 6-second video is captured per field of view. By repeatedly capturing images by slightly moving the telescope around the field of view, the gaps between the image sensors are filled. Furthermore, by changing the angle of the telescope and repeatedly capturing video, 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 video data acquired by wide-field video survey observations. 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 probability of being truly moving objects above a certain value. Thousands of moving objects are detected per night, and the majority of these are thought to be known or unknown artificial celestial bodies.

[0033] The typical apparent speed of a Near Earth Object (NEO) is about 3 arcsec [s －１ ]. Apparent speed ~ 50 arcsec [s －１ ] 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 NEO contamination, such as space debris or artificial satellites.

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

[0035] The data of artificial celestial objects obtained by wide-field video surveys are the observation position, i.e., angle and brightness, of a given frame per field of view. The observation position of an artificial celestial object is determined by the imaging time t ｉ RA (α ｉ ), and Dec(δ ｉ ), brightness is visible light magnitude m ｉ 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 ０ , t １ , and t ２ The artificial celestial body 30 is moving along an orbit 34. ０ The observation position of the artificial celestial body 30 in RA (α ０ ), and Dec(δ０ ) is expressed as the image capture time t １ , and t ２ The same is true for .

[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 comprises 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 (computational resources) such as a personal computer (PC), a mainframe, a workstation, a cloud computing system, etc. The I / O 43 is an interface, 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 performed by the CPU 41. The input unit 45 is an input device such as a keyboard or a mouse, and is realized by a device that receives input from a user and transmits it to the CPU 41.

[0040] The block diagram of Figure 4 shows the functional units within the CPU 41. When each functional unit of the CPU 41 is realized by software, the CPU 41 realizes it by executing instructions of a program, which is 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 for orbital elements and celestial body sizes from observations of one field of view of the sky sweeping observation, based on observation data from the sky sweeping observation acquired by the I / O 43 .

[0042] The artificial celestial body 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 body by comparing the posterior distributions for different observations, i.e., performs 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 now be described. The orbit estimation method according to this embodiment is a method for identifying artificial celestial objects and determining their orbits when optical observations of a short arc obtained over multiple nights in sky surveys, including wide-field video surveys, are obtained. Figure 5 shows a flowchart illustrating 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 the observation of one field of view during 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 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 the details of step S502 shown in FIG.

[0051] (Calculating oscillating parameters) In step S601, oscillating parameters are calculated. First, the time t obtained from the observation of one field of view is calculated. ｉ , and position αｉ , δ ｉ Based on the following equation, the average time t ａｖｇ The average position α in ａｖｇ , δ ａｖｇ and the average speed v α，ａｖｇ , v δ，ａｖｇ Calculate the average time t ａｖｇ , and the average position α ａｖｇ , δ ａｖｇ is given by the following formula:

[0052]

[0053]

[0054]

[0055] Also, α ｉ , δ ｉ Therefore, the velocity at each time v α，ｉ , v δ，i but

[0056]

[0057]

[0058] is given by the average velocity v α，ａｖｇ , v δ，ａｖｇ teeth,

[0059]

[0060]

[0061] is given by

[0062] The tangent element is the average time t ａｖｇ The geocentric distance r, the time derivative of the geocentric distance r dr / dt, and the average position α ａｖｇ , δ ａｖｇ , average speed v α，ａｖｇ , v δ，ａｖｇ It should be noted 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, α ｉ , δ ｉand the average position α ａｖｇ , δ ａｖｇ , average speed v α，ａｖｇ , v δ，ａｖｇ In the graph shown in FIG. 7, data 71 shown in gray scale in the graph is α ｉ , δ ｉ The grayscale intensity is ｉ The horizontal and vertical axes correspond to α ｉ , δ ｉ The grayscale indicates the difference from the average time. The data 72 in the center of the graph is the average position α ａｖｇ , δ ａｖｇ and the vector from data 72 to data 73 is (v α，ａｖｇ , v δ，ａｖｇ ) multiplied by 0.25 seconds.

[0064] (Procedure for Obtaining Orbital Elements from Tangential Elements) In step S602, since the orbital elements (Keplerian Elements) are composed of six parameters, just like the tangential elements, and each parameter is in a variable transformation relationship with another, 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. Geocentric distance r, RAα at the observer center Ｇｅｏ , Decδ Ｇｅｏ Then the vector from the center of the Earth to the observer is

[0066]

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

[0068]

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

[0070]

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

[0072] Tangent elements r, dr / dt, α ａｖｇ , δ ａｖｇ , v α，ａｖｇ , v δ，ａｖｇ Given the above, we can calculate the position of the celestial body in the geocentric coordinate system, and further calculate the velocity of the celestial body from these small displacements. That is, the position vector in the geocentric coordinate system

[0073]

[0074] and the velocity vector

[0075]

[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 MCMC (Markov Chain Monte Carlo) sampling of the tangential trajectory elements. Note that Log Likelihood is maximum likelihood estimation, a method in statistics for point-estimating the parameters of the probability distribution to which given data follows. Furthermore, MCMC is a general term for various algorithms that sample a probability distribution by creating a Markov chain that has the desired probability distribution as its balanced distribution.

[0078] First, for each sample, the observation time t ｉ RA (α ｍ，ｉ ), Decδ ｍ，ｉ Calculate the observation α ｉ , δ ｉ are the α of the model, respectively. ｍ，ｉ , δ ｍ，ｉ These differences are assumed to be normally distributed around the range of the positioning accuracy of the observations.

[0079]

[0080]

[0081] It should be noted that the difference in the RA direction needs to be multiplied by cos(δ).

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

[0083]

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

[0085] The orbital elements obtained from the tangential orbital elements are free of perturbations, while from these orbital elements α ｍ，ｉ , δ ｍ，ｉ When calculating LL, perturbations can be taken into account and an algorithm such as SGP4 can be used. The difference between taking perturbations into account and not taking them into account appears as an offset in the global position within the timescale range 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, Δα ｉ , Δδ ｉ is divided into global and local components. Global component Δα ｇｌｏｂ , Δδ ｇｌｏｂ is averaged over i, and the following formula is used:

[0087]

[0088]

[0089] Calculation is based on the local component Δα ｌｏｃ，ｉ , Δδ ｌｏｃ，ｉ is the displacement from the global component, and is expressed by the following formula:

[0090]

[0091]

[0092] The global component is given by the typical difference σ ｇｌｏｂ Assuming that the local component follows a normal distribution with the standard deviation of the positioning accuracy of the observation, Log Likelihood is

[0093]

[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 considered.

[0095] In an example of actual observation, the position determination accuracy of the observation is 0.3 arcsec. ｇｌｏｂ The value of is estimated to be about 20 arcmin as a typical value. In the following, when calculating the position of the artificial celestial body from the orbital elements, perturbations are taken into account using SGP4.

[0096] (Setting the initial value of 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 α ａｖｇ , δ ａｖｇ , average speed v α，ａｖｇ , v δ，ａｖｇ For the tangent elements r and dr / dt, the 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 values ​​that can be taken is limited. For example, the tangent element r ranges from an altitude of about 100 km to an altitude of about 380,000 km above the moon, and dr / dt ranges from ±5 km / s. －１ 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 the other tangential elements. A point that gives a closed orbit or a point with a large LL can be used as the initial value for MCMC.

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

[0099] (Setting Prior Distribution) In step S605, a prior distribution can be set when performing MCMC sampling, which will be described later. Bayesian estimation begins 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 sampling the tangent elements, the two-dimensional distribution of the tangent elements r and dr / dt and the average position α ａｖｇ , δ ａｖｇ , average speed v α，ａｖｇ , v δ，ａｖｇ 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 orbital elements of known artificial celestial bodies, etc. The mean position α ａｖｇ , δ ａｖｇ , average speed v α，ａｖｇ , v δ，ａｖｇ For example, a normal distribution can be used with a standard deviation obtained from 12 frame observations.

[0101] (Sampling of Tangent Elements) In step S606, when sampling tangent elements using MCMC, sampling is performed using the following three patterns (1) to (3): (1) Sampling with one parameter of tangent element r (2) Sampling with two parameters of tangent element r and dr / dt (3) Sampling with tangent element r, dr / dt, and average velocity vα，ａｖｇ , v δ，ａｖｇ In (1), we assume that dr / dt = 0 and the average position α ａｖｇ , δ ａｖｇ , average speed v α，ａｖｇ , v δ，ａｖｇ is fixed as the observed value. In (2), the unsampled mean position α ａｖｇ , δ ａｖｇ , average speed v α，ａｖｇ , v δ，ａｖｇ is fixed as the observed value. In (3), the unsampled mean position α ａｖｇ , δ ａｖｇ is fixed to the observed value. While (2) and (3) are suitable for sampling various cases, including orbits with large eccentricities, they tend to stochastically overestimate the eccentricity when the orbit of the target celestial body is close to a circular orbit. (1) cannot completely eliminate the tendency to stochastically overestimate the eccentricity when the orbit of the target celestial body 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 tangent element r MCMC sampling is performed on one parameter of tangent element r. As described above, Metropolis-Hastings (MH method) sampling or Gibbs sampling is used from the obtained initial value, and sampling is performed from the Log Likelihood equation (Equation 19). The number of steps until convergence to the posterior distribution is less than 100, and 200 steps are sufficient for burn-in. Furthermore, 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 acceptance 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, stable convergence was confirmed by setting the width of the proposal distribution so that the acceptance probability was approximately 20 to 40%. Figures 9(a), 9(b), and 9(c) show examples of the trace, likelihood, and acceptance probability in burn-in. The likelihood shown in Figure 9(b) is normalized so that its maximum value is 1. The grayscale of the data in Figures 9(a), 9(b), and 9(c) indicates the parallel sampling, i.e., each chain started from a different initial value. Around the 30th step, the influence of the initial value disappears, and the decrease in the acceptance probability stops, as shown in Figure 9(c).

[0104] FIG. 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 initial values. Since only one parameter is sampled, all parameters except r are fixed. The accuracy of determining r is tens of thousands of km.

[0105] Once the tangent elements are sampled, the corresponding orbital elements can be obtained. Figure 11 shows a graph of the posterior distribution of Keplerian elements, and Figure 12 shows a graph of the posterior distribution of equinoctial elements. As in Figure 10, the graph enclosed by dashed line A shows the one-dimensional posterior distribution of each parameter, the graph enclosed by dashed line B shows the two-dimensional posterior distribution of each parameter pair, and the graph enclosed by dashed line C shows the two-dimensional Log Likelihood. The dotted lines in the figure indicate the values ​​of known artificial celestial bodies matched by position and velocity. In tangent element space, all parameters except r have a value of 1, but in orbital element space, the orbital elements take on values ​​corresponding to each tangent element set, so all parameters have a broad 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 zero, e = 0.00433. Because the orbit is nearly circular, the argument of perigee ω is indeterminate, and therefore the mean anomaly M is also indeterminate. In Figure 11, the distributions of these parameters do 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 nearly circular. In fact, as shown in Figure 12, in the posterior distribution of 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. Examples of the results of sampling with two parameters, r and dr / dt, are shown in Figures 13 to 16. The details of Figures 13 to 16 are the same as the details of sampling with one parameter shown in Figures 9 to 12.

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

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

[0110]

[0111] Here, m is the apparent magnitude of the sun, which is -26.47 in visible light (SDSS g band). Regarding the directivity of the reflection of the artificial celestial body, if we assume that it scatters isotropically only in the direction in which sunlight hits it in a spherical model, then F(ψ) = 1 / 2π. Furthermore, let s be the typical size of the artificial celestial body, and A = πs ２ As shown in Figure 21, when the artificial celestial body 210 is spherical, if the angle between the line of sight 212 to the artificial celestial body as seen from the observation device 211 and the 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]

[0113] A in the above formula is A ｓｈｉｎｅ Substituting this,

[0114]

[0115]

[0116] Therefore, the size of the artificial celestial body is

[0117]

[0118] It is expressed as:

[0119] From the above equation, once the distance d and albedo ρ are determined, the size of the celestial body can be calculated from the observed magnitude m. Once 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 albedo ρ is obtained using the brightness of observed values ​​that match known celestial bodies in actual observations. The brightness fluctuations of the same celestial body, except for the distance and phase angle, are due to the rotation and attitude of the celestial body, etc. These are all pushed into the albedo term, and the above equation is used to obtain the range of variation in ρ. By setting the mean 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 directly as the prior distribution, or a normal distribution fitted to this ρ distribution can be used as the prior distribution. The normal distribution fitted to this ρ distribution is shown as a dotted line in the graph in Figure 22. Hereinafter, the normal distribution fitted to this ρ distribution will be referred to as the prior distribution.

[0120] The observed magnitude, the sampled tangent element, and the prior distribution of albedo give a posterior distribution of object sizes. The observed magnitude is the average m ａｖｇ as

[0121]

[0122] is used and calculated using equation (24).

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

[0124] (Matching of Multiple Visual Fields Using Rough Posterior Distributions) In step S503 shown in FIG. 5 , multiple visual fields are matched using a posterior distribution obtained by MCMC from the observation of one visual field. First, any combination of three visual field observations 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 visual fields, i.e., the common portion of the posterior distributions, is calculated, and the area or volume of the obtained common portion is referred to as 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. Furthermore, 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 seven-dimensional space. Either the area or 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 body, and small otherwise. If the degree of agreement exceeds an appropriate threshold, we proceed to obtain a high-precision posterior distribution. Among the seven parameters, the h parameter, which determines the orbital plane, is ｘ , h ｙ In this graph, 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 intersections when all three observations are of the same celestial object, and Figure 25 shows examples when one field of view is different. Figures 24(a) to 24(f) and 25(a) to 25(f) show examples of equinoctial elements, while 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 on the left vertical axis, and the intensity of the grayscale color indicates different fields of view. The intersections are indicated by dots and are represented on the right vertical axis. The vertical dotted lines indicate the values ​​of the matched known celestial objects.

[0126] (Acquisition of High-Precision Posterior Distribution) In step S504, 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 benefit to sampling tangent elements. Therefore, here, the equinoctial elements themselves are sampled.

[0127] MultiNest has the advantage that convergence to the posterior distribution is relatively robust even when the log likelihood is non-convex or when there are many peaks. Taking advantage of this advantage, equinoctial element sampling can be performed relatively robustly. Addressing multi-peaks is essential for robust trajectory estimation. This is because, for example, when there is a reasonable fit in all fields of view and when there is a good fit only in a certain field of view (12 frames), the log likelihood becomes large, i.e., multi-peaks occur. Furthermore, when the observation times between fields of view are separated by, for example, about one day, a small change in the parameters results in a large change in log likelihood, i.e., the gap between multi-peaks is large, making it difficult for random walk sampling such as MCMC to move between multi-peaks.

[0128] Furthermore, 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, which allows the effects of perturbations to be reflected in the posterior distribution of the orbital elements that are finally obtained.

[0129] (Setting Priors) When using MultiNest, priors (prior distributions) of the parameters to be sampled are calculated. 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 may be obtained from observations of three fields of view.

[0130] When obtaining the prior of equinoctial elements from observations of three fields of view, the observation α of field of view j ｊ，ｉ , δ ｊ，ｉ From representative point α ｊ，ｏ , δ ｊ，ｏ This can be the first or last observation of 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 manner as sampling of the tangent elements using MCMC in the following three patterns (1) to (3): (1) Sampling with one parameter of the tangent element r; (2) Sampling with two parameters of the tangent element r and dr / dt; (3) Sampling with the tangent element r, dr / dt, and the average velocity v α，ａｖｇ , v δ，ａｖｇ Here, the reason why tangential elements can be used is that the observation points in each field of view are narrowed down to one, thereby reducing the multi-peak nature.

[0132] After sampling the tangential elements, the posterior distribution of the equinoctial elements obtained from sampling one, two, and four parameters may be added together. This posterior distribution (probability density function) is used as the prior of the equinoctial elements provided to MultiNest. In practice, the inverse function (percentile point function, ppf) of the cumulative probability density function is provided.

[0133] Figure 26 shows an example of creating a prior for equinoctial elements using three observation points. The graphs on the left side of Figure 26 are probability density functions, with histogram 261 obtained from sampling of tangential elements and step 262 being a numerical correction. The graphs on the right side of Figure 26 are cumulative density functions calculated from the probability density functions.

[0134] Three-field observation α ｊ，ｉ , δ ｊ，ｉ Using all, equinoctial elements are sampled. Here, since SDP4 is also used for orbital propagation, an equation that does not consider offsets due to perturbations is adopted as the log likelihood. By sampling with MultiNest, the posterior distribution of equinoctial elements is obtained. FIG. 27 shows an example of the posterior distribution of equinoctial elements sampled with MultiNest. Also, FIG. 28 shows an example of the posterior distribution of Keplerian elements sampled with MultiNest. For example, for a, the width of the distribution is on the order of 10 m, and it can be seen that the posterior distribution is obtained with high accuracy.

[0135] As mentioned above, the present invention naturally includes various embodiments not described herein. Therefore, the technical scope of the present invention is defined only by the invention-specifying matters according to the scope of the claims that are appropriate from the above description.

[0136] REFERENCE SIGNS LIST 1 Artificial satellite 2 Space debris 3, 4, 5 Orbit 6 Forecast circle 20, 30 Artificial celestial body 21 Earth 22 Perigee 23 Apogee 24 Orbital plane 25 Equatorial plane 26 Orbital direction 33 Telescope 34 Orbit 40 Orbit estimation device 41 CPU 42 Memory unit 43 I / O (input / output interface) 44 Display unit 45 Input unit 411 Low-precision posterior distribution calculation unit 412 Matching unit 413 High-precision posterior distribution calculation unit 71, 72, 73 Data 210 Artificial celestial body 211 Observation device 212, 213 Line of sight 214 Sunlight

Claims

1. An orbit estimation device for multiple artificial celestial bodies that performs orbit estimation using Bayesian estimation with MCMC based on sky sweeping observations, comprising: a memory unit that stores data for each of the multiple artificial celestial bodies obtained by the sky sweeping observations; 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 sweeping observations; a man-made celestial body identification unit that identifies observations of the same artificial celestial body by comparing the posterior distributions for different observations among the sky sweeping 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.

2. The orbit estimation device according to claim 1, characterized in that the low-precision posterior distribution calculation unit calculates the orbital elements from the time and position of each of the multiple artificial celestial bodies obtained by the sky sweeping observation, and samples the calculated orbital elements.

3. The orbit estimation device described in claim 2, characterized in that, when calculating the orbital elements, the low-precision posterior distribution calculation unit calculates an average position and average velocity at the average time as tangential elements from the time and position of each of the multiple artificial celestial bodies, and calculates the orbital elements from the calculated average position and average velocity.

4. The orbit estimation device according to claim 1, characterized in that the low-precision posterior distribution calculation unit obtains a low-precision posterior distribution of the celestial body sizes from the brightness of each of the plurality of artificial celestial bodies.

5. The orbit estimation device described in claim 1, characterized in that the artificial celestial body identification unit: selects any combination of observations from three fields of view; compares the low-precision posterior distribution for the selected combination; calculates the degree of agreement between the combinations of observations from the three fields of view; and, if the calculated degree of agreement exceeds a predetermined threshold, considers the selected combinations to be those of the same artificial celestial body.

6. The orbit 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. The orbit estimation device according to claim 6, characterized in that, when sampling the equinoctial elements using MultiNest, the high-precision posterior distribution calculation unit obtains the Prior by obtaining each observation point of the combination of observations of the three fields of view.

9. The orbit estimation device described in claim 8, characterized in that the high-precision posterior distribution calculation unit, after obtaining the observation points for each combination of observations of the three fields of view, samples tangential elements of each of the obtained observation points, and adds up the posterior distributions of the equinoctial elements.

10. The trajectory estimation device described in 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 the common part of the equinoctial elements obtained in matching of multiple views using a coarse posterior distribution.

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

12. A method for estimating the orbits of multiple artificial celestial bodies, which estimates the orbits of the multiple artificial celestial bodies by Bayesian estimation using MCMC based on sky surveying observations, comprising: a low-precision posterior distribution calculation step for obtaining 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 step for identifying observations of the same artificial celestial body by comparing the posterior distributions for different observations among the sky surveying observations; and a high-precision posterior distribution calculation step for obtaining a high-precision posterior distribution of orbital elements using observations of three or more short arcs of the same artificial celestial body.

13. An orbit estimation device for multiple artificial celestial bodies, comprising a memory unit for storing data on each of multiple artificial celestial bodies obtained by sky sweeping observations, and a calculation unit, and performing orbit estimation by Bayesian estimation using MCMC based on the sky sweeping observations, for realizing in the orbit estimation device: a low-precision posterior distribution calculation function for obtaining a low-precision posterior distribution of orbital elements and celestial body size from short arcs of the orbits of each of the multiple artificial celestial bodies obtained by the sky sweeping observations; an artificial celestial body identification function for identifying observations of the same artificial celestial body by comparing the posterior distributions for different observations among 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.

Images