CN115080238A - High-precision multi-satellite ephemeris forecasting method based on GPU parallel computation - Google Patents

High-precision multi-satellite ephemeris forecasting method based on GPU parallel computation Download PDF

Info

Publication number
CN115080238A
CN115080238A CN202210744067.0A CN202210744067A CN115080238A CN 115080238 A CN115080238 A CN 115080238A CN 202210744067 A CN202210744067 A CN 202210744067A CN 115080238 A CN115080238 A CN 115080238A
Authority
CN
China
Prior art keywords
satellite
ephemeris
gpu
file
forecasting
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202210744067.0A
Other languages
Chinese (zh)
Inventor
袁春强
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Guosheng Haoxing Suzhou Aerospace Technology Co ltd
Original Assignee
Guosheng Haoxing Suzhou Aerospace Technology Co ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Guosheng Haoxing Suzhou Aerospace Technology Co ltd filed Critical Guosheng Haoxing Suzhou Aerospace Technology Co ltd
Priority to CN202210744067.0A priority Critical patent/CN115080238A/en
Publication of CN115080238A publication Critical patent/CN115080238A/en
Pending legal-status Critical Current

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F9/00—Arrangements for program control, e.g. control units
    • G06F9/06—Arrangements for program control, e.g. control units using stored programs, i.e. using an internal store of processing equipment to receive or retain programs
    • G06F9/46—Multiprogramming arrangements
    • G06F9/50—Allocation of resources, e.g. of the central processing unit [CPU]
    • G06F9/5005—Allocation of resources, e.g. of the central processing unit [CPU] to service a request
    • G06F9/5027—Allocation of resources, e.g. of the central processing unit [CPU] to service a request the resource being a machine, e.g. CPUs, Servers, Terminals
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F9/00—Arrangements for program control, e.g. control units
    • G06F9/06—Arrangements for program control, e.g. control units using stored programs, i.e. using an internal store of processing equipment to receive or retain programs
    • G06F9/46—Multiprogramming arrangements
    • G06F9/50—Allocation of resources, e.g. of the central processing unit [CPU]
    • G06F9/5005—Allocation of resources, e.g. of the central processing unit [CPU] to service a request
    • G06F9/5011—Allocation of resources, e.g. of the central processing unit [CPU] to service a request the resources being hardware resources other than CPUs, Servers and Terminals
    • G06F9/5016—Allocation of resources, e.g. of the central processing unit [CPU] to service a request the resources being hardware resources other than CPUs, Servers and Terminals the resource being the memory
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T1/00—General purpose image data processing
    • G06T1/20—Processor architectures; Processor configuration, e.g. pipelining
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2209/00—Indexing scheme relating to G06F9/00
    • G06F2209/50—Indexing scheme relating to G06F9/50
    • G06F2209/5018—Thread allocation
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

The invention discloses a high-precision multi-satellite ephemeris forecasting method based on GPU parallel computing. Firstly, the high-precision orbit prediction model is provided, and model parameters are divided into four types, namely integral variables, satellite parameters, constants, parameters related to model files and the like. Then, an allocation scheme and a forecast algorithm flow of the GPU video memory are given. The multi-satellite ephemeris forecasting method provided by the invention fully utilizes the parallel computing capability of the GPU, ensures the accuracy of the ephemeris, has higher computing efficiency, can provide input for operations such as multi-satellite sub-satellite point forecasting and coverage forecasting and the like, and can support the application of real-time performance such as collision early warning and the like. Compared with a method based on CPU multi-core parallel or PC cluster parallel, the method has the advantages that the hardware construction cost is low, and the method is suitable for ephemeris forecast of commercial constellations.

Description

High-precision multi-satellite ephemeris forecasting method based on GPU parallel computation
Technical Field
The invention relates to the fields of astronomical mechanics and aerospace measurement and control, in particular to a high-precision multi-satellite ephemeris forecasting method based on GPU parallel computing.
Background
With the development of human space technology, more and more satellites or constellations are transmitted to a space orbit, high-precision multi-satellite ephemeris forecast is a core technology for monitoring the in-orbit running state of the satellites or the constellations and is a necessary function of a satellite measurement and control center, and typical operations comprise satellite down-pointing forecast, satellite-to-ground coverage forecast, collision forecast and the like which all depend on ephemeris forecast. The existing high-precision multi-satellite ephemeris forecast is mainly based on the multi-core concurrency capability of a server CPU or a concurrent computing mode of a multi-PC cluster, the hardware cost is high, and the operation and maintenance management of commercial constellations is not utilized.
With the development of the computational capability and the storage capacity of the GPU, particularly after the CUDA parallel framework introduced by NVIDIA corporation, the GPU is gradually applied to other fields from graphics development, such as the fields of hydrodynamics simulation, image processing, aerospace, and the like. At present, the latest GPU-RTX3090Ti of NVIDIA corporation, about 1.5 ten thousand RMB, already has 10752 core processors, an acceleration frequency of 1.86GHz and a video memory capacity of 24 GB. A large number of core processors and high-capacity video memories can enable programs which need to be operated on a server or a cluster to be operated on a single GPU even if the computing power of a single core is slightly weaker than that of a PC, and on the premise that the computing efficiency meets the actual requirement, the hardware cost is greatly reduced. The advent of high-level GPUs undoubtedly provides a low-cost implementation environment for the parallel forecasting of satellite ephemeris. In recent years, researchers at home and abroad have realized an analytical model such as a two-body model and an SGP4 on a GPU, but a GPU forecast version based on a high-precision orbit model has not been reported yet, and besides the reason that the GPU display memory capacity is insufficient in the past years, the main reason is that the inquiry of a JPL ephemeris (solar ephemeris and lunar ephemeris) and the calculation of an atmospheric density model depend on a software package provided at home and are extremely difficult to realize on the GPU again.
Disclosure of Invention
In order to overcome the defects, the invention provides a high-precision multi-satellite ephemeris forecasting method based on GPU parallel computing, which makes full use of the parallel computing capability of the GPU, ensures the accuracy of ephemeris, has higher computing efficiency, can provide input for operations such as multi-satellite down-satellite point forecasting and coverage forecasting, and can support the application of real-time performance such as collision early warning.
The technical scheme adopted by the invention for solving the technical problem is as follows:
a high-precision multi-satellite ephemeris forecasting method based on GPU parallel computing is characterized in that an ephemeris forecasting process of each satellite is used as an independent thread and deployed to a GPU to operate, and high-precision parallel forecasting of multi-satellite ephemeris is achieved, and the method specifically comprises the following steps:
step S1, a high-precision orbit prediction model is given, namely, the gravity F of the center of mass of the earth is considered in the inertial system of the earth 0 Global non-spherical perturbation force f g Solar attraction f s Moon attraction force f m Atmospheric resistance f d And solar light pressure radiation pressure f p The motion mechanics model of (2):
Figure BDA0003716390670000021
wherein the content of the first and second substances,
Figure BDA0003716390670000022
and
Figure BDA0003716390670000023
respectively representing the position, the speed and the acceleration vector of the satellite at the moment t, wherein m is the mass of the satellite;
step S2, dividing the parameters of the high-precision orbit prediction model into four types, namely an integral variable, a satellite parameter, an astronomical constant and a model file, wherein:
the integral variable comprising the position of the satellite
Figure BDA0003716390670000031
And velocity
Figure BDA0003716390670000032
The satellite parameters comprise satellite mass m and aspheric order N g Satellite surface to mass ratio
Figure BDA0003716390670000033
Coefficient of atmospheric resistance C d Satellite reflection factor C r ；
The astronomy constants include an earth gravity constant mu and an earth reference radius R e "Shitai" medicineConstant of positive attraction mu s Moon gravitational constant μ m Solar radiation pressure P s ；
Step S3, pre-fetching an earth gravitational field parameter file, generating a solar ephemeris file, a lunar ephemeris file and an atmospheric density file, and loading the files into a GPU memory;
the solar ephemeris file and the lunar ephemeris file are generated based on a JPL ephemeris, and the atmospheric density file is generated for each satellite in advance by utilizing the existing atmospheric model software;
step S4, allocating video memory for the integral variable, the satellite parameter and the parameter related to the model file, and copying the initial value of the integral variable, the satellite parameter, the astronomical constant and the parameter related to the model file from the GPU memory to the GPU video memory;
step S5, taking the forecast of each satellite as a thread, solving the motion mechanics model on the GPU based on the KSG integrator value, and storing the calculation result into the position and speed array of each satellite;
and step S6, copying the position and speed array of each satellite from the GPU video memory to the memory, and storing the ephemeris forecast result.
As a further improvement of the present invention, in step S5, when solving the motion mechanics model, a pre-calculation method is used to generate a solar ephemeris file and a lunar ephemeris file, and after loading the files into the GPU video memory, the positions of the sun and the moon are obtained by using an inquiry method.
As a further improvement of the present invention, in step S5, when the motion mechanics model is solved and the atmospheric resistance is calculated, an atmospheric density file of each satellite is generated in a pre-calculation manner, and is loaded into the GPU for display, and then the atmospheric density at the position of the satellite is obtained in an interpolation manner.
As a further improvement of the invention, in the steps S1-S6, the ephemeris forecasting process runs completely on the GPU and is calculated in parallel by using a CUDA framework.
As a further improvement of the invention, the satellite comprises an artificial satellite and an in-orbit object.
As a further improvement of the invention, the in-orbit object comprises a large space station and space debris of unequal size.
The invention has the beneficial effects that: the method adopts the GPU to realize high-precision multi-satellite ephemeris forecast, fully utilizes the parallel computing capability of the GPU, ensures the ephemeris precision, has higher computing efficiency, can provide input for operations such as multi-satellite undersatellite point forecast, coverage forecast and the like, can support the application of real-time performance such as collision early warning and the like, has lower hardware construction cost compared with a CPU multi-core parallel or PC cluster parallel based method, and is suitable for ephemeris forecast of commercial constellations.
Detailed Description
The following describes the specific embodiments in detail by taking a single GPU as an example. The embodiments described herein are only used for intuitively explaining the present invention, and are not used for limiting the present invention, and the present invention may also be implemented by using a plurality of GPUs, especially in the case of a large number of satellites or a short integration step.
A high-precision multi-satellite ephemeris forecast method based on GPU parallel computing is disclosed, wherein the multi-satellite ephemeris forecast refers to the positions and the speeds of a plurality of known satellites at initial time, and the satellites are predicted in a future time period [ T ] according to a kinetic model of the satellites 1 ,T 2 ]The position and speed of the vehicle. The multi-satellite ephemeris forecasting method takes the ephemeris forecasting process of each satellite as a single thread to be deployed to a GPU for operation, so that the multi-satellite ephemeris high-precision parallel forecasting is realized, and the method specifically comprises the following steps:
firstly, a high-precision orbit prediction model is given, namely, in the earth inertial system, the gravity F of the earth mass center is considered 0 Global non-spherical perturbation force f g Solar attraction f s Moon attraction force f m Atmospheric resistance f d And solar light pressure radiation pressure f p The motion mechanics model of (2):
Figure BDA0003716390670000051
wherein the content of the first and second substances,
Figure BDA0003716390670000052
and
Figure BDA0003716390670000053
respectively, the position, velocity and acceleration vectors of the satellite at time t, and m is the mass of the satellite.
Wherein, the gravity of the center of mass F of the earth 0 Is the main term, other forces are small amounts compared to it, but not negligible for high accuracy predictions. The initial velocity and the initial acceleration of the satellite are respectively
Figure BDA0003716390670000054
And
Figure BDA0003716390670000055
formally, acceleration vector
Figure BDA0003716390670000056
Which can be abbreviated as t, can be abbreviated as,
Figure BDA0003716390670000057
and
Figure BDA0003716390670000058
but in practice it is also related to many other parameters.
From the two body problem, gravity of center of mass F of the earth 0 Comprises the following steps:
Figure BDA0003716390670000059
where μ is the earth constant, and F 0 The main parameters of interest include satellite position
Figure BDA00037163906700000519
Satellite mass m and earth gravity constant μ.
Since the earth is not a regular sphere, or even an ellipsoid, its gravitational effect on spacecraft cannot be easily handled asThe central gravitational force. Irregular latitude of earth to geocentric
Figure BDA00037163906700000510
Longitude of center of earth is lambda and distance of center of earth is
Figure BDA00037163906700000518
The perturbation function of the gravitational field generated by the particle is U:
Figure BDA00037163906700000511
wherein R is e Is the earth's reference radius; distance between the earth and the center
Figure BDA00037163906700000512
X, Y and Z are position components of the satellite in an inertial coordinate system;
Figure BDA00037163906700000513
and
Figure BDA00037163906700000514
is a normalized gravity coefficient; n is a radical of g Is the order of the asphericity;
Figure BDA00037163906700000515
to relate to
Figure BDA00037163906700000516
The associated legendre polynomial of (a). The offset derivative is calculated by the formula (3) to calculate the aspheric perturbation force f of the earth g ：
Figure BDA00037163906700000517
From the above formula, with f g The main parameters of interest include satellite position
Figure BDA0003716390670000061
Radius of reference of the earth R e Earth's gravityConstant mu, aspheric order N g And gravity coefficient set
Figure BDA0003716390670000062
The earth gravitational field model adopted in the implementation is EGM2008, and a gravitational coefficient set S can be loaded from an earth gravitational field parameter file.
Gravitational forces of the sun and moon (f) s And f m ) Relative position of sun or moon to earth, and non-spherical perturbation force f of earth g Similarly, it can be abbreviated as relating to satellite position
Figure BDA0003716390670000063
Position of the sun
Figure BDA0003716390670000064
Position of moon
Figure BDA0003716390670000065
Constant of solar attraction mu s And the moon gravitational constant mu m Function of (c):
Figure BDA0003716390670000066
in the above formula, the position of the sun
Figure BDA0003716390670000067
And the position of the moon
Figure BDA0003716390670000068
Defined under the earth's inertial system.
Figure BDA0003716390670000069
And
Figure BDA00037163906700000610
respectively coming from a solar ephemeris file and a lunar ephemeris file, wherein the two ephemeris files record the position of the sun or the moon line by line, and the time step is 1 minute. The first number of each row is UTC time, and the last three are solar or lunar inertial systemsThe position component of (a). The implementation pre-manufactures the solar and lunar ephemeris files based on the JPL ephemeris, and loads the solar and lunar positions to the video memory at one time according to the forecast starting and stopping time during operation.
For low orbit satellites, atmospheric drag is one of the main perturbations. Atmospheric resistance f d Is the coefficient of atmospheric resistance C d Atmospheric density ρ d (location with satellite)
Figure BDA00037163906700000611
Related), the areal quality ratio of the satellite
Figure BDA00037163906700000612
(A is the effective cross-sectional area of the satellite) and the satellite velocity
Figure BDA00037163906700000613
Function of (c):
Figure BDA00037163906700000614
in the above formula, the atmospheric density ρ d The change mechanism is very complex and not only changes with altitude, but also relates to changes of solar activity, time, season, latitude and geomagnetic activity. In order to simplify the problem, the implementation does not implement the calculation process of a typical atmosphere model (such as MSIS) again on the GPU, but uses the existing atmosphere model software to generate an atmosphere density file for each satellite in advance, loads the atmosphere density file into a GPU video memory as required during running, and calculates the atmosphere density file by using nearest neighbor interpolation based on the position and time.
Given the initial position and velocity of the satellite, the satellite trajectory is an elliptic curve defined under the inertial frame, assuming that the satellite motion belongs to a two-body model. Over a longer period of time in the future, such as a month, the actual positions of the satellites are substantially distributed along this elliptical curve without excessive bias. Uniformly sampling on an elliptic curve to obtain a discrete satellite position set
Figure BDA0003716390670000071
N r For the number of sampling positions, the time sequence is obtained by uniformly sampling a long time period in the future
Figure BDA0003716390670000072
Taking the time step length as 1 minute, N t The number of time samples, the atmospheric density ρ of the satellite d Can approximate the discrete function phi s Represents:
Figure BDA0003716390670000073
in the above formula, discrete function
Figure BDA0003716390670000074
Equivalent to one N r ×N t And storing the matrix into a file to obtain the atmospheric density file of the satellite.
The satellite is subjected to direct radiation pressure of the sun in the movement process, the sunlight pressure is a corrosion factor v (when the satellite is in a terrestrial shadow, v is 0, otherwise v is 1), and the solar radiation pressure P at the position of the earth orbit radius from the sun s Satellite reflection factor C r Satellite surface-to-mass ratio
Figure BDA0003716390670000075
Satellite position
Figure BDA0003716390670000076
Position of the sun
Figure BDA0003716390670000077
Function of (c):
Figure BDA0003716390670000078
satellites of different heights may employ different kinetic models
Figure BDA0003716390670000079
Of low-earth-orbit satellitesThe main perturbation forces are non-spherical perturbation forces, atmospheric resistance and sun-moon attraction, and a higher non-spherical order, such as N, is required g 60; the medium orbit satellite does not need to consider atmospheric resistance, sun-moon attraction and aspheric perturbation force, and the aspheric order is set to be N g 40; the high-orbit satellite does not need to consider atmosphere, sunlight pressure and sun-moon gravitation, and the aspheric order is set to be N g ＝20。
Given the initial position and velocity of the satellite, incorporating an accurate kinetic model
Figure BDA0003716390670000081
And its parameters (see table 1), ephemeris forecast uses a numerical integrator to perform the prediction calculations of satellite position and velocity. The implementation employs a KSG integrator, which is a multi-step integrator of fixed order and fixed step size. Ephemeris predictions are a sequence of positions and velocities defined in the earth's inertial frame.
As shown in Table 1, the present embodiment is a mechanical model
Figure BDA0003716390670000082
The parameters of (2) are divided into 4 types, such as integral variable (predicted value), satellite parameter, constant and parameter related to model file.
The first is the integration variable, i.e., the position and velocity of the satellite, and since the KSG integrator is a multi-step integrator, i.e., the predictor is a linear combination of historical data, the present implementation allocates data storage for the entire forecast time period for each satellite at once. According to the forecast of 1 minute step length and 7 days (other time step lengths and time period forecasts are similar), the video memory required by N satellites is as follows: n (60 × 24 × 7 × 8)/(1024 × 1024) ≈ 0.54 × N (mb), the first number 7 indicates a 7-day forecast, the second number 7 indicates 7 double-precision floating-point numbers including time, position, and velocity, and the number 8 indicates an 8-byte double-precision floating-point number.
The second category is five satellite parameters, i.e.
Figure BDA0003716390670000083
Requiring the advance allocation of these parameters to each satelliteAnd storing and copying the corresponding numerical value from the memory to the video memory in advance. If N is present g The integer is 4 bytes, and the other are double-precision floating point numbers, so the video memory required by the N satellites is as follows: n × (8+4+8+8+8)/(1024 × 1024) ≈ 0.000036 × N (mb).
The third category is five astronomical constants, i.e., μ, R, associated with the earth, sun or moon e ,μ s ,μ m ,P s And macro definition is adopted in the program, so that video memory is not consumed.
The fourth class is 4 parameters associated with the model file. Wherein the first is the gravity coefficient set S, if N g 60, the consumed video memory does not exceed N g *(N g +1)/2 × 8/(1024 × 1024) ≈ 0.02 (MB). The second and third parameters are the sun and moon positions, and the required video memory is: 2 (60 × 24 × 7 × 4 × 8)/(1024 × 1024) ≈ 0.62(MB), numeral 2 refers to both sun and moon, and numeral 4 refers to 4 double-precision floating-point numbers including time and position components. Fourth is the atmospheric density ρ d Rho, as well as the integral variables and satellite parameters d And is also a satellite-specific parameter, and video memory needs to be allocated to each satellite in advance. If the position sampling number of the elliptical orbit is N r When 360, the required display memory is: n (N) r 60 × 24 × 7 × 1 × 8)/(1024 × 1024) ≈ 27.69 × n (mb). The fourth type requires video memory in common: 27.69 × N +0.64 (MB).
Mechanical model
Figure BDA0003716390670000091
The total required display memory for the 4 types of parameters is about 28.23 × N (mb), where N is the number of satellites. The apparent atmospheric density was about 27.69 × N (time step 1 minute), accounting for 98%. The 1 minute step size 7 days for 500 satellites forecasts requires display 14115 (MB). It can be seen that if the time step of the atmospheric density sampling is properly enlarged to be k times of the original time step, for example, 5 times, i.e., from 1 minute to 5 minutes, the memory consumption will be correspondingly reduced by k times, and the number of the satellites that can be accommodated is also approximately increased to be k times of the original time step.
TABLE 1 acting force and its parameter table
Figure BDA0003716390670000092
Figure BDA0003716390670000101
Given the initial position and velocity of the satellite, incorporating an accurate kinetic model
Figure BDA0003716390670000102
The table is table 1, and ephemeris forecast utilizes a numerical integrator based prediction calculation of satellite position and velocity. The implementation employs a KSG integrator, which is a multi-step integrator of fixed order and fixed step size. Ephemeris predictions are a sequence of positions and velocities defined in the earth's inertial frame.
In the implementation, ephemeris forecast of N satellites is realized by adopting a CUDA parallel framework, and the method is totally divided into 4 steps: the method comprises the steps of firstly, loading an earth gravitational field parameter file, a solar ephemeris file, a lunar ephemeris file and an atmospheric density file to a memory. And secondly, distributing a video memory for the integral variable, the satellite parameter and the parameter related to the model file, and copying the initial value of the integral variable, the satellite parameter and the parameter related to the model file from the memory to the video memory. And thirdly, taking the forecast of each satellite as a thread, solving a formula (1) on the GPU based on the KSG integrator value, and storing the calculation result into the position and speed array of each satellite. And fourthly, copying the position and speed array of each satellite from a video memory to a memory, and storing the ephemeris forecast result.
In the previous description, numerous specific details were set forth in order to provide a thorough understanding of the present invention. The foregoing description is only a preferred embodiment of the invention, which can be embodied in many different forms than described herein, and therefore the invention is not limited to the specific embodiments disclosed above. And that those skilled in the art may, using the methods and techniques disclosed above, make numerous possible variations and modifications to the disclosed embodiments, or modify equivalents thereof, without departing from the scope of the claimed embodiments. Any simple modification, equivalent change and modification of the above embodiments according to the technical essence of the present invention are within the scope of the technical solution of the present invention.

Claims (6)

1. A high-precision multi-satellite ephemeris forecasting method based on GPU parallel computation is characterized by comprising the following steps: the method comprises the following steps of taking the ephemeris forecasting process of each satellite as a separate thread, deploying the ephemeris forecasting process to a GPU for operation, and realizing high-precision parallel forecasting of the multi-satellite ephemeris, wherein the method specifically comprises the following steps:
step S1, a high-precision orbit prediction model is given, namely, the gravity F of the center of mass of the earth is considered in the inertial system of the earth 0 Global non-spherical perturbation force f g Solar attraction f s Moon attraction force f m Atmospheric resistance f d And solar light pressure radiation pressure f p The motion mechanics model of (2):
Figure FDA0003716390660000011
wherein the content of the first and second substances,
Figure FDA0003716390660000012
and
Figure FDA0003716390660000013
respectively representing the position, the speed and the acceleration vector of the satellite at the moment t, wherein m is the mass of the satellite;
step S2, dividing the parameters of the high-precision orbit prediction model into four types, namely an integral variable, a satellite parameter, an astronomical constant and a model file, wherein:
the integral variable comprising the position of the satellite
Figure FDA0003716390660000014
And velocity
Figure FDA0003716390660000015
The satellite parameters comprise satellite mass m and aspheric order N g Satellite surface to mass ratio
Figure FDA0003716390660000016
Coefficient of atmospheric resistance C d Satellite reflection factor C r ；
The astronomy constants include an earth gravity constant mu and an earth reference radius R e Sun attraction constant mu s Moon gravitational constant μ m Solar radiation pressure P s ；
Step S3, pre-fetching an earth gravitational field parameter file, generating a solar ephemeris file, a lunar ephemeris file and an atmospheric density file, and loading the files into a GPU memory;
the solar ephemeris file and the lunar ephemeris file are generated based on a JPL ephemeris, and the atmospheric density file is generated for each satellite in advance by utilizing the existing atmospheric model software;
step S4, allocating video memory for the integral variable, the satellite parameter and the parameter related to the model file, and copying the initial value of the integral variable, the satellite parameter, the astronomical constant and the parameter related to the model file from the GPU memory to the GPU video memory;
step S5, taking the forecast of each satellite as a thread, solving the motion mechanics model on the GPU based on the KSG integrator value, and storing the calculation result into the position and speed array of each satellite;
and step S6, copying the position and speed array of each satellite from the GPU video memory to the memory, and storing the ephemeris forecast result.
2. The GPU parallel computation-based high-precision multi-satellite ephemeris forecasting method of claim 1, comprising: in step S5, when the kinematic mechanical model is solved, a pre-calculation method is used to generate a solar ephemeris file and a lunar ephemeris file, and the solar ephemeris file and the lunar ephemeris file are loaded into the GPU video memory and then the positions of the sun and the moon are obtained by using an inquiry method.
3. The method for forecasting the high-precision multi-satellite ephemeris based on GPU parallel computing as claimed in claim 1, wherein: in step S5, when the motion mechanics model is solved and the atmospheric resistance is calculated, an atmospheric density file of each satellite is generated in a pre-calculation manner, and is loaded into the GPU for display, and then the atmospheric density at the position of the satellite is obtained in an interpolation manner.
4. The method for forecasting the high-precision multi-satellite ephemeris based on GPU parallel computing as claimed in claim 1, wherein: in the steps S1 to S6, the ephemeris forecasting process runs entirely on the GPU and is calculated in parallel by using the CUDA framework.
5. The method for forecasting the high-precision multi-satellite ephemeris based on GPU parallel computing as claimed in claim 1, wherein: the satellite comprises a man-made satellite and an on-orbit object.
6. A high-precision multi-satellite ephemeris forecasting method based on GPU parallel computing according to claim 5, characterized in that: the in-orbit objects include large space stations and unequal-size space debris.
CN202210744067.0A 2022-06-27 2022-06-27 High-precision multi-satellite ephemeris forecasting method based on GPU parallel computation Pending CN115080238A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210744067.0A CN115080238A (en) 2022-06-27 2022-06-27 High-precision multi-satellite ephemeris forecasting method based on GPU parallel computation

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210744067.0A CN115080238A (en) 2022-06-27 2022-06-27 High-precision multi-satellite ephemeris forecasting method based on GPU parallel computation

Publications (1)

Publication Number Publication Date
CN115080238A true CN115080238A (en) 2022-09-20

Family

ID=83255998

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210744067.0A Pending CN115080238A (en) 2022-06-27 2022-06-27 High-precision multi-satellite ephemeris forecasting method based on GPU parallel computation

Country Status (1)

Country Link
CN (1) CN115080238A (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116015401A (en) * 2022-12-06 2023-04-25 中科南京移动通信与计算创新研究院 Satellite constellation model construction method, system, equipment and storage medium

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116015401A (en) * 2022-12-06 2023-04-25 中科南京移动通信与计算创新研究院 Satellite constellation model construction method, system, equipment and storage medium

Similar Documents

Publication Publication Date Title
CN111680354B (en) Method for calculating self-intersection point of orbit of near-earth regression orbit satellite subsatellite point and photographing point
CN108820260B (en) Medium-term orbit forecasting method and device for low-orbit spacecraft and storage medium
Doornbos Thermospheric density and wind determination from satellite dynamics
CN107066641B (en) The numerical computation method and system of extensive space junk Distribution evolution
Dachwald et al. Impact of optical degradation on solar sail mission performance
Marcos Accuracy of atmospheric drag models at low satellite altitudes
CN110595485A (en) Low-orbit satellite long-term orbit forecasting method based on two-line root number
Wang et al. Propagation errors analysis of TLE data
CN114327919B (en) Space target collision early warning method and system
Lee et al. Maximizing photovoltaic power generation of a space-dart configured satellite
CN115080238A (en) High-precision multi-satellite ephemeris forecasting method based on GPU parallel computation
CN115258197A (en) Spacecraft orbit terminal point prediction method and device, processor and electronic equipment
CN111127295B (en) SGP4 orbit model integrated parallel method based on GPU
Smith et al. Ionospheric drag for accelerated deorbit from upper low earth orbit
Srivastava et al. Satellite ephemeris prediction for the Earth orbiting satellites
Golikov THEONA—a numerical-analytical theory of motion of artificial satellites of celestial bodies
CN114580181B (en) Huge constellation coverage performance parallel computing method based on CUDA
Aleksandrova et al. Numerical modeling of motion of near-Earth objects in a parallel computing environment
CN116125503A (en) High-precision satellite orbit determination and prediction algorithm
Ravanbakhsh et al. Multidisciplinary design optimization approach to conceptual design of a LEO earth observation microsatellite
Allasio et al. GOCE mission: design phases and in-flight experiences
Saghari et al. “Wasted performance” minimization of the multi-purpose mini-satellite platform for an EO mission using a dynamic simulation-based model
Podgorny et al. Investigation of the mechanism of a solar flare by means of MHD simulations above the active region in real scale of time: The choice of parameters and the appearance of a flare situation
Wei et al. Redesign of high-precision reference orbit for interferometric SAR satellite with injection error
Schäff Low-thrust multi-revolution orbit transfers

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination