An orbit prediction method for space-based Internet of Things terminals without ephemeris update

Abstract

The present application discloses a method for orbit prediction of a space-based Internet of Things terminal under the condition that the ephemeris is not updated, the method comprising: data initialization, and obtaining TLE data at a preset time; calling a two-row root number fitting prediction model, and calculating two-row root number prediction values ​​based on the prediction model and preset TLE data; and calculating the entry and exit time of the satellite according to the two-row root number prediction values. The technical solution in the present application is helpful to improve the accuracy of orbit prediction of the space-based Internet of Things terminal under the condition that the ephemeris is not updated.

Description

Technical Field

[0001] The present application relates to the technical field of satellite orbit prediction, and in particular to an orbit prediction method for a space-based Internet of Things terminal under the condition that the ephemeris is not updated. Background Art

[0002] Currently, the space-based Internet of Things is under construction, and the number of communication satellites in the sky is limited, so it is not possible to achieve full coverage. However, by using satellite orbit prediction technology to obtain the entry and exit times of satellites in advance, the terminal can be actively controlled to send and receive data after the satellite enters and before it leaves the country, and the terminal can be controlled to sleep during the rest of the time. This can not only greatly improve the communication success rate of the terminal, but also greatly reduce the power consumption of the terminal. Therefore, satellite orbit prediction is an important technology for space-based Internet of Things terminals to improve the communication success rate and reduce power consumption.

[0003] The existing orbit prediction algorithm based on the SGP4 model relies on the two-row root parameters released by the North American Aerospace Defense Command. The ground measurement and control network of the North American Aerospace Defense Command tracks and determines the target in space, and then measures the flight orbit of the target, and then publishes it on the website in the format of two-row root parameters. The two-row root data consists of only two lines of data, which include parameters such as ephemeris number, TLE duration, BSTAR drag modulation coefficient, orbit inclination, ascending node right ascension, orbit eccentricity, perigee angular distance, mean anomaly angle, and mean motion.

[0004] The existing orbit prediction algorithm for space-based IoT terminals is to fix the pitch angle threshold and the latest TLE data downloaded from the Internet into the flash memory of the terminal device. Each time the orbit is predicted, the fixed TLE data is read from the flash memory and substituted into the SGP4 model to calculate the satellite position vector. Then, based on the latitude and longitude information of the IoT terminal at the current moment, the terminal position vector is calculated, and the pitch angle of the ground IoT terminal and the satellite is calculated based on the position vectors of the satellite and the terminal. The terminal then determines whether the satellite has entered the country by judging whether the calculated pitch angle value exceeds the set threshold, thereby obtaining the satellite's entry and exit time in advance.

[0005] Since the orbit prediction algorithm of the above-mentioned space-based Internet of Things terminal relies on two rows of roots, in order to obtain more accurate prediction results, it is necessary to frequently update the TLE data fixed in the terminal.

[0006] In the prior art, for a space-based IoT terminal device, there are mainly three methods for updating the TLE data fixed in the terminal: (1) manually writing flash updates on-site through relevant instructions; (2) relying on some means of traditional networks for wireless remote updates; (3) relying on satellite communication links for downlink updates.

[0007] However, since most space-based IoT terminals are deployed in wild areas such as forests, oceans, and mines, these areas are basically in harsh environments, not covered by traditional networks, and inaccessible to humans. As a result, it becomes very difficult for space-based IoT terminals with no or little downlink to frequently update TLE data. Over time, the error in orbit prediction based on fixed TLE parameters will become larger and larger.

[0008] Therefore, a method for predicting the orbit of a space-based Internet of Things terminal is needed, which can improve the orbit prediction accuracy of the space-based Internet of Things terminal within a certain period without updating the ephemeris. Summary of the invention

[0009] The purpose of this application is to improve the accuracy of orbit prediction of space-based Internet of Things terminals without updating the ephemeris.

[0010] The technical solution of the present application is: to provide an orbit prediction method for a space-based Internet of Things terminal under the condition that the ephemeris is not updated, the method comprising: data initialization, and obtaining TLE data at a preset time; calling a two-row root number fitting prediction model, and calculating two-row root number prediction values ​​based on the prediction model and preset TLE data; and calculating the entry and exit times of the satellite according to the two-row root number prediction values.

[0011] In any of the above technical solutions, further, data initialization specifically includes: initialization forecast start time, forecast days, terminal longitude and latitude information, satellite entry elevation angle threshold.

[0012] In any of the above technical solutions, further, according to the two-row root prediction values, the entry and exit times of the satellite are calculated, specifically including: substituting the two-row root prediction values ​​into the SGP4 satellite orbit prediction model to calculate the satellite position vector of the satellite in the geocentric earth-fixed coordinate system from the forecast start time to the end time; calculating the terminal position vector of the terminal in the geocentric earth-fixed coordinate system from the forecast start time to the end time according to the terminal longitude and latitude information; calculating the elevation angle from the terminal to the satellite from the forecast start time to the end time according to the satellite position vector and the terminal position vector, and when it is determined that the elevation angle is greater than or equal to the satellite entry elevation angle threshold, calculating the entry and exit times of the satellite.

[0013] In any of the above technical solutions, further, the two rows of root number prediction values ​​include: mean anomaly angle value, and calculating the two rows of root number prediction values ​​specifically includes:

[0014] According to the established mean anomaly fitting prediction model, the mean anomaly value at the forecast start time is calculated. The calculation formula of the mean anomaly value is:

[0015] m1=(y0+((tx-t0)*n1*360 / (24*3600)))%360.0

[0016] Where m1 is the mean anomaly, y0 is the mean anomaly at the reference point, tx-t0 is the number of seconds from the forecast start time to the reference point, n1 is the average number of circles around the earth per day, and % is the modulo operation.

[0017] In any of the above technical solutions, further, the average number of orbits around the earth per day n1 is the average value obtained by calculating the average value of the predicted number of orbits around the earth per day n2 of the satellite from the reference point time to the forecast start time at intervals of days, wherein the calculation formula for the predicted number of orbits around the earth per day n2 of the satellite is:

[0018] n2=p11*t2^3+p12*t2^2+p13*t2+p14

[0019] Where n2 is the number of times the satellite orbits the earth every day, t2 is the number of seconds from the current time to the earliest sample historical TLE, and p11, p12, p13, and p14 are the first satellite fitting parameters.

[0020] In any of the above technical solutions, further, the two rows of root number prediction values ​​include: perigee argument prediction value, and calculating the two rows of root number prediction values ​​specifically includes:

[0021] According to the established perigee argument fitting prediction model, the perigee argument prediction value at the forecast start time is calculated. The calculation formula of the perigee argument prediction value is:

[0022] ω1=(p21*t+p22)%360.0

[0023]

[0024] Where ω1 is the calculated process value, ω2 is the predicted value of the perigee argument, t is the number of seconds from the forecast start time to the earliest sample historical TLE, and p21 and p22 are the fitting parameters of the second satellite.

[0025] In any of the above technical solutions, further, the method also includes: correcting the predicted value of the perigee argument, and calculating the entry and exit time of the satellite according to the corrected predicted value of the perigee argument, wherein the calculation formula of the corrected predicted value of the perigee argument is:

[0026]

[0027] Wherein, ω3 is the predicted value of the corrected perigee argument angle, % is the modulo operation, x1, x2, x3 are the corrected values, ω2_1, ω2_2, ω2_3, ω2_4 are the corrected angle interval values.

[0028] In any of the above technical solutions, further, the two-row root number prediction value includes: orbital eccentricity value, and calculating the two-row root number prediction value specifically includes:

[0029] According to the established orbital eccentricity fitting prediction model, the orbital eccentricity value at the forecast start time is calculated. The calculation formula of the orbital eccentricity value is:

[0030] ecc=a0+a1*cos(w*t)+b1*sin(w*t)-c1*t

[0031] Where ecc is the orbit eccentricity, t is the number of seconds from the forecast start time to the earliest sample historical TLE, and a0, a1, b1, c1, and w are the eccentricity fitting parameters.

[0032] In any of the above technical solutions, further, the two rows of predicted values ​​of the root number include: the right ascension value of the ascending node, and the two rows of predicted values ​​of the root number are calculated, specifically including:

[0033] According to the established ascending node right ascension fitting prediction model, the right ascension value of the ascending node at the start of the forecast is calculated. The calculation formula of the right ascension value of the ascending node is:

[0034] Ω=(p31*t+p32)%360.0

[0035] Where Ω is the right ascension of the ascending node, t is the number of seconds from the forecast start time to the earliest sample historical TLE, and p31 and p32 are the fitting parameters of the third satellite.

[0036] In any of the above technical solutions, further, the two rows of element number prediction values ​​include: orbital inclination value, and calculating the two rows of element number prediction values ​​specifically includes:

[0037] According to the established orbital inclination fitting prediction model, the orbital inclination value at the forecast start time is calculated. The calculation formula of the orbital inclination value is:

[0038] x=t%31536000

[0039] y0=i0-Δi*(t / 31536000)

[0040] i=p41*x^9+p42*x^8+p43*x^7+p44*x^6+p45*x^5+p46*x^4+p47*x^3+p48*x^2+p49*x+y0

[0041] Where x is the number of seconds used for calculation, i is the orbital inclination, t is the number of seconds from the forecast start time to the earliest sample historical TLE, y0 is the initial angle of the orbital inclination with a periodicity of one year, % is the modulo operation, i0 is the orbital inclination value corresponding to the earliest sample TLE data, Δi is the total change value of the orbital inclination corresponding to the sample TLE data in one year, and p41, p42, p43, p44, p45, p46, p47, p48, and p49 are the fitting parameters of the fourth satellite.

[0042] The beneficial effects of this application are:

[0043] Compared with the existing method of simply fixing the TLE data of a certain moment in the space-based Internet of Things terminal and performing orbit prediction based on the fixed TLE parameters, the technical solution in the present application effectively improves the accuracy of orbit prediction of the space-based Internet of Things terminal under the condition that the ephemeris is not updated. The technical solution in the present application enables the space-based Internet of Things terminal to not need to frequently update the latest TLE parameter data, and realizes a high-precision orbit prediction under the condition that the ephemeris is not updated, thereby improving the communication success rate of satellite communication relying on orbit prediction;

[0044] The technical solution in this application can be extended to the application of space-based Internet of Things terminal equipment that adopts the SGP4 orbit prediction algorithm and has no ephemeris update method, such as space-based Internet of Things data collection equipment with little or no downlink and inaccessible to human power, which will help popularize the space-based Internet of Things. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] The advantages of the above and / or additional aspects of the present application will become apparent and easily understood in the description of the embodiments in conjunction with the following drawings, in which:

[0046] Figure 1 is a schematic flow chart of a method for orbit prediction of a space-based Internet of Things terminal under the condition that the ephemeris is not updated according to an embodiment of the present application;

[0047] Figure 2 is a detailed flow chart of a trajectory prediction method according to an embodiment of the present application;

[0048] Figure 3 is a schematic diagram of a BSTAR drag modulation coefficient fitting and a sample curve according to an embodiment of the present application;

[0049] Figure 4 is a schematic diagram of orbital inclination fitting and sample comparison curves according to an embodiment of the present application;

[0050] Figure 5 is a schematic diagram of a curve comparing the right ascension of the ascending node with a sample according to an embodiment of the present application;

[0051] Figure 6is a schematic diagram of orbital eccentricity fitting and sample comparison curves according to an embodiment of the present application;

[0052] Figure 7 is a schematic diagram of perigee argument fitting and sample comparison curves according to an embodiment of the present application;

[0053] Figure 8 is a schematic diagram of a mean anomaly fitting and sample comparison curve according to an embodiment of the present application;

[0054] Fig. 9 is a schematic diagram of a fitting and sample comparison curve of the number of circles around the earth / day according to an embodiment of the present application;

[0055] Fig.10 is a schematic diagram of a fitting and sampling curve before perigee argument correction according to an embodiment of the present application;

[0056] Fig.11 It is a schematic diagram of the fitting and sampling curves after the perigee argument correction according to an embodiment of the present application.

[0057] Among them, the dotted line in the figure is the fitting data curve, and the solid line in the figure is the sampling data curve. DETAILED DESCRIPTION

[0058] In order to more clearly understand the above-mentioned purposes, features and advantages of the present application, the present application is further described in detail below in conjunction with the accompanying drawings and specific implementation methods. It should be noted that the embodiments of the present application and the features in the embodiments can be combined with each other without conflict.

[0059] In the following description, many specific details are elaborated to facilitate a full understanding of the present application. However, the present application may also be implemented in other ways different from those described herein. Therefore, the protection scope of the present application is not limited to the specific embodiments disclosed below.

[0060] It is understandable to those skilled in the art that the errors caused by not updating the two-row root parameters for a long time have the most direct impact on the space-based Internet of Things terminals that rely on orbit forecast results to control data communication, which is a greatly reduced communication success rate. Therefore, in applications where the two-row root parameters cannot be updated frequently, how to improve the orbit forecast accuracy within the longest possible period is crucial to the communication success rate and power consumption control of the space-based Internet of Things terminals.

[0061] like Figure 1As shown, this embodiment provides an orbit prediction method for a space-based Internet of Things terminal under the condition that the ephemeris is not updated. The method analyzes the historical TLE data of the predicted satellite (i.e., satellite ephemeris, also known as two-line orbital data (Two-LineOrbital Element, TLE)), focusing on 7 parameters related to orbit prediction: BSTAR drag modulation coefficient, orbital inclination, ascending node right ascension, orbital eccentricity, perigee angular distance, mean anomaly, and average motion. The historical data of these parameters are fitted to establish a prediction model, so that the latest TLE parameter value of the satellite at any time can be calculated. The fitted predicted parameters are substituted for the original solidified TLE data and substituted into the SGP4 model for calculation, so that the accuracy of orbit prediction can be greatly improved within a certain period. The method includes:

[0062] Step 1: Initialize the data and obtain the TLE data at the preset time;

[0063] Step 2, calling the two-row root number fitting prediction model, and calculating the two-row root number prediction value based on the prediction model and the preset TLE data;

[0064] Step 3, calculate the entry and exit times of the satellite based on the two rows of root number prediction values.

[0065] In this embodiment, the method is based on the orbit prediction algorithm of the SGP4 / SDP4 model, and the space-based Internet of Things terminal with the STM32 series MCU as the control chip is used as the prototype. The Internet of Things terminal includes satellite communication, GNSS data acquisition, sensor data acquisition and processing, data storage, orbit prediction and other functional modules. In this method, a large amount of historical data is analyzed in the early stage to establish a fitting prediction model for each parameter of the two-row root, and in the orbit prediction algorithm, the parameter values ​​predicted by the two-row root fitting prediction model are used instead of the TLE data of a certain moment of the solid setting for calculation, thereby realizing the improvement of the orbit prediction accuracy of the space-based Internet of Things terminal within a certain period under the condition that the ephemeris is not updated.

[0066] like Figure 2 As shown, the above method in this embodiment specifically includes:

[0067] Step 1: The orbit prediction module initializes the forecast start time, forecast days, terminal longitude and latitude information, and satellite entry elevation angle threshold, and reads the TLE data of the preset time set in the flash memory, and parses the TLE data of the preset time, where the preset time can be the latest time updated before the terminal is put into field use.

[0068] Among them, the start time of the forecast cannot be earlier than the earliest epoch time in the TLE fitting prediction model and the epoch time of the two rows of roots that are fixed and set; the longitude and latitude information of the terminal comes from the acquisition and analysis of the GNSS module data; the satellite entry elevation angle threshold can be set, and is generally set to 10 degrees for space-based Internet of Things terminals. The mean anomaly angle corresponding to the TLE data fixed in the terminal flash memory is the initial value of the mean anomaly angle at the reference point time for calculating the mean anomaly angle in step 2.

[0069] It should be noted that the data required for the calculation process of this embodiment can be derived and analyzed based on preset TLE data, such as the mean anomaly angle at the reference point time, and the specific process will not be repeated here.

[0070] Step 2: Call the two-row root number fitting prediction model, and calculate the two-row root number prediction value based on the prediction model and the preset TLE data. The two-row root number prediction value includes multiple parameters. The specific steps are as follows:

[0071] For a more specific description, taking the NORAD satellite No. 44486 target star as an example, a two-row root number fitting prediction model was established based on the historical two-row root numbers from January 1, 2020 to December 31, 2020.

[0072] (1) According to the established BSTAR drag modulation coefficient fitting prediction model, the BSTAR drag modulation coefficient at the forecast start time is calculated;

[0073] Specifically, according to the BSTAR drag modulation coefficient fitting prediction model of the target star of NORAD satellite No. 44486, the BSTAR drag modulation coefficient is calculated, and the corresponding calculation formula is:

[0074] bstar=(a11*sin(w11*t+w12)+a12*sin(w13*t+w14)

[0075] +a13*sin(w15*t+w16)+a14)*(1.0e-09)

[0076] Where t is the number of seconds from the forecast start time to the earliest sample historical TLE duration, e is the symbol of scientific notation, bstar is the calculated BSTAR drag modulation coefficient, and a11, a12, a13, a14, w11, w12, w13, w14, w15, and w16 are BSTAR drag modulation fitting coefficients.

[0077] It should be noted that for different satellites, the values ​​of the BSTAR drag modulation fitting coefficients are different. In this embodiment, the above formula can be obtained by fitting:

[0078] bstar=(5.215e+03*sin(5.9e-08*t-1.258)+4.342e+04*sin(4.607e-06*t+1.706)+3.968e+04+1525*sin(1.122e-06*t-2.416))*(1.0e-09).

[0079] In this embodiment, the BSTAR drag modulation coefficient fitting and sample comparison curve are as follows: Figure 3 shown.

[0080] (2) Calculate the orbital inclination value at the forecast start time based on the established orbital inclination fitting prediction model;

[0081] Specifically, the calculation formula corresponding to the orbital inclination fitting prediction model of the NORAD satellite number 44486 target star is:

[0082] y0=i0-Δi*(t / 31536000)

[0083] i=p41*x^9+p42*x^8+p43*x^7+p44*x^6+p45*x^5+p46*x^4+p47*x^3+p48*x^2+p49*x+y0

[0084] Where x is the number of seconds used for calculation, i is the orbital inclination, t is the number of seconds from the forecast start time to the earliest sample historical TLE, y0 is the initial angle of the orbital inclination with a periodicity of one year, % is the modulo operation, i0 is the orbital inclination value corresponding to the earliest sample TLE data, Δi is the total change value of the orbital inclination corresponding to the sample TLE data in one year, and p41, p42, p43, p44, p45, p46, p47, p48, and p49 are the fitting parameters of the fourth satellite.

[0085] In this embodiment, the above formula can be converted to:

[0086] x=t%31536000

[0087] y0=97.592918046467157-0.035302310512776*(t / 31536000)

[0088] i=3.82614660121e-68*x^9-3.476294247930169e-60*x^8+1.150674884118187e-54*x^7+9.649466411807340e-45*x^6-4.016130435245 305e-37*x^5+6.379324007326139e-30*x^4-3.380564016008760e-2 3*x^3-2.197045289343381e-17*x^2-1.363845987310940e-09*x+y0.

[0089] In this embodiment, the orbital inclination fitting and sample comparison curve are as follows: Figure 4 shown.

[0090] It should be noted that technical researchers have analyzed the orbital inclination data curves obtained from TLE data of multiple low-orbit satellites for more than two years and found that the curve change pattern of the orbital inclination repeats with a period of about one year. Therefore, selecting one year's data as a sample and using one year as a segment to perform segmented fitting prediction of the orbital inclination will be more accurate than directly fitting the data.

[0091] In the above orbital inclination fitting prediction model, x is the number of seconds from the current moment in the current segment interval (1 year) to the starting point of the segment, that is, the number of seconds used for calculation, 0.035302310512776 is the angle reduction value of the orbital inclination in one year, and y0 is the initial angle of the orbital inclination in each segment interval (i.e. 1 year) of the prediction model.

[0092] In this way, a fitting prediction model is established only for the orbital inclination data of the segmented interval (1 year). By combining the calculation of x, y0 and the fitting prediction model, the orbital inclination i at any time greater than the earliest sample historical TLE duration can be calculated.

[0093] (3) Calculate the right ascension of the ascending node at the start of the forecast based on the established right ascension fitting prediction model;

[0094] Specifically, the calculation formula corresponding to the fitting prediction model of the ascending node right ascension of the NORAD satellite number 44486 target star is:

[0095] Ω=(p31*t+p32)%360.0

[0096] Where Ω is the right ascension of the ascending node, t is the number of seconds from the forecast start time to the earliest sample historical TLE, and p31 and p32 are the fitting parameters of the third satellite.

[0097] In this embodiment, the above calculation formula for the target star of NORAD satellite number 44486 can be converted to:

[0098] Ω=(1.139345573666438e-05*t+77.972487600318090)%360.0.

[0099] In this embodiment, the ascending node right ascension fitting and sample comparison curve are as follows: Figure 5 shown.

[0100] It should be noted that through analyzing the sampling data of the right ascension of the ascending node, technical researchers found that the right ascension of the ascending node changes periodically between 0-360 degrees, and its change pattern tends to linear change.

[0101] Therefore, when establishing a fitting prediction model for the right ascension of the ascending node, the sampling data of the right ascension of the ascending node are first integrated, and the segmented linear curves are integrated into a continuous linear curve, that is, the sampling data are added with the corresponding angle value 360*(n-1) in sequence according to the segments, where n represents the number of segments. For example, the second segment of the curve adds 360 degrees on the basis of the sampling value, and the third segment of the curve adds 720 degrees on the basis of the sampling value, and so on. A set of integrated sampling values ​​is obtained, and finally a continuous linear curve is formed. The integrated continuous linear curve is fitted with data and modulo 360 degrees to obtain a 0-360 degree segmented fitting prediction curve.

[0102] (4) Calculate the orbit eccentricity value at the forecast start time based on the established orbit eccentricity fitting prediction model;

[0103] Specifically, the calculation formula corresponding to the orbital eccentricity fitting prediction model of the NORAD satellite number 44486 target star is:

[0104] ecc=a0+a1*cos(w*t)+b1*sin(w*t)-c1*t

[0105] Where ecc is the orbit eccentricity, t is the number of seconds from the forecast start time to the earliest sample historical TLE, and a0, a1, b1, c1, and w are the eccentricity fitting parameters.

[0106] It should be noted that for different satellites, the values ​​of the eccentricity fitting parameters a0, a1, b1, c1, and w are different.

[0107] In this embodiment, the calculation formula of the orbit eccentricity fitting prediction model can be converted to:

[0108] ecc=0.002119+1.17414927e-04*cos(6.896e-07*t)+0.00015943473*sin(6.896e-07*t)-(6.651882e-13)*t.

[0109] In this embodiment, the orbital eccentricity fitting and sample comparison curve are as follows: Figure 6 shown.

[0110] It should be noted that, through analyzing the sampling data of orbit eccentricity, technical researchers found that the change curve of orbit eccentricity is a sine curve with a gradually downward offset. Therefore, when establishing a fitting prediction model for orbit eccentricity, Fourier fitting is performed on the orbit eccentricity data, and the change value of the offset is added to the fitting model, that is, the offset parameter 0.002119-(6.651882e-13)*t in the above model.

[0111] (5) Calculate and correct the predicted value of the perigee argument at the forecast start time based on the established perigee argument fitting prediction model;

[0112] Specifically, the calculation formula of the perigee argument fitting prediction model of the NORAD satellite No. 44486 target star is:

[0113] ω1=(p21*t+p22)%360.0

[0114]

[0115]

[0116] Wherein, ω1 is the calculated process value, ω2 is the predicted value of the perigee argument, t is the number of seconds from the forecast start time to the earliest sample historical TLE, p21 and p22 are the fitting parameters of the second satellite, ω3 is the corrected predicted value of the perigee argument, % is the modulo operation, x1, x2, x3 are the corrected values, ω2_1, ω2_2, ω2_3, ω2_4 are the corrected angle interval values.

[0117] It should be noted that for different satellites, the correction values ​​x1, x2, x3 and the correction angle interval values ​​ω2_1, ω2_2, ω2_3, ω2_4 have different values.

[0118] In this embodiment, the above formula can be converted into:

[0119] ω1=(-3.951714645710906e-05*t+1.586901637646532e+03)%360.0

[0120]

[0121]

[0122] Where ω1 is the calculated process value, ω2 is the calculated perigee argument prediction value, ω3 is the corrected perigee argument prediction value, t is the number of seconds from the forecast start time to the earliest sample historical TLE, and % is the modulo operation.

[0123] In this embodiment, the perigee argument fitting and sample comparison curve are as follows: Figure 7 shown.

[0124] It should be noted that through analyzing the sampling data of the perigee argument, technical researchers found that the perigee argument changes periodically between 360-0 degrees, and its changing pattern tends to linear change.

[0125] Therefore, when establishing a fitting prediction model for the perigee argument, the sampling data of the perigee argument are first integrated, and the segmented linear curve is integrated into a continuous linear curve, that is, the sampling data is sequentially added with the corresponding angle value 360*(nm) according to the segment, where n represents the total number of segments, and m represents the mth segment. For example, the total number of segments n=5, the first segment curve adds 360*(5-1) degrees on the basis of the sampling value, the second segment curve adds 360*(5-2) degrees on the basis of the sampling value, and so on, to form a set of integrated sampling values, and finally form a continuous linear curve. Then the integrated continuous linear curve is fitted with data and modulo 360 degrees. At the same time, if the calculated initial perigee argument value ω1<0, then ω2=ω1+360.0, the fitting prediction curve of the 360-0 degree segment can be obtained.

[0126] According to the fitting curve, it can be seen that the curve obtained by linear fitting has a certain error with the sampling curve. Therefore, the fitted piecewise linear curve is corrected in each segment, which can improve the accuracy of fitting prediction. The comparison curves before and after correction are as follows: Fig.10 and Fig.11 shown.

[0127] (6) Calculate the mean anomaly value at the forecast start time according to the established mean anomaly fitting prediction model;

[0128] Specifically, the calculation formula of the mean anomaly fitting prediction model of the NORAD satellite No. 44486 target star is:

[0129] m1=(y0+((tx-t0)*n1*360 / (24*3600)))%360.0

[0130] Wherein, m1 is the mean anomaly, which is the predicted value, y0 is the mean anomaly at the reference point, that is, the initial value of the mean anomaly. In actual application, its value is the mean anomaly corresponding to the TLE data fixed in the terminal flash memory, unit: degree, tx-t0 is the number of seconds from the forecast start time to the reference point time, n1 is the average number of circles around the earth per day, where tx is the number of seconds from the forecast start time to 0:0:00 on January 1, 1970, and t0 is the number of seconds from the reference TLE time to 0:0:00 on January 1, 1970.

[0131] In this embodiment, the mean anomaly angle fitting and sample comparison curve are as follows: Figure 8 shown.

[0132] In order to establish the fitting prediction model of the mean anomaly angle and calculate the mean anomaly angle value at the forecast start time, it is also necessary to predict the number of times the satellite orbits the earth every day. In addition, the average value of the predicted number of times the satellite orbits the earth every day, n2, is calculated with the interval of days from the reference point to the forecast start time, and the average value of the number of times the satellite orbits the earth every day, n1, is obtained, and then substituted into the calculation formula of the mean anomaly angle for calculation, and the satellite mean anomaly angle m1 at the forecast start time is obtained.

[0133] (7) Based on the established daily orbital number prediction model, the number of orbits the satellite will make every day at the start of the forecast is calculated. The corresponding calculation formula is:

[0134] n2=p11*t2^3+p12*t2^2+p13*t2+p14

[0135] Where n2 is the number of times the satellite orbits the earth every day, t2 is the number of seconds from the current time to the earliest sample historical TLE, and p11, p12, p13, and p14 are the first satellite fitting parameters.

[0136] It should be noted that, for different satellites, the values ​​of the first satellite fitting parameters p11, p12, p13, and p14 are different.

[0137] In this embodiment, the above formula can be converted into:

[0138] n2=2.814e-26*t2^3+2.003e-19*t2^2+1.094e-10*t2+15.07105

[0139] In the formula, n2 is the number of times the satellite orbits the earth every day, and t2 is the number of seconds from the forecast start time to the earliest sample historical TLE duration time, that is, the number of seconds from the forecast start time to the TLE duration time corresponding to the TLE parameter with the earliest TLE duration in the historical TLE data sample used for data fitting.

[0140] In this embodiment, the fitting and sample comparison curves of the number of circles around the earth every day are as follows: Fig. 9 shown.

[0141] Step 3: Substitute the two rows of root prediction values ​​predicted in step 2 into the SGP4 satellite orbit prediction model to calculate the satellite position vector in the Earth-centered Earth-fixed coordinate system from the forecast start time to the end time, which specifically includes calculating the satellite's average motion and semi-major axis, calculating atmospheric drag and gravitational perturbations, calculating long-period terms and short-period terms, calculating the instantaneous orbital elements, and calculating the satellite's position in the Earth-centered Earth-fixed coordinate system.

[0142] Step 4: Based on the terminal longitude and latitude information obtained from the GNSS module in step 1, calculate the terminal position vector of the terminal from the forecast start time to the end time in the Earth-centered Earth-fixed coordinate system.

[0143] Step 5: Calculate the elevation angle from the terminal to the satellite from the forecast start time to the forecast end time based on the position vectors of the terminal and the satellite in the Earth-centered Earth-fixed coordinate system from the forecast start time to the forecast end time.

[0144] The entry and exit times of the satellite are calculated by determining whether the elevation angle is greater than or equal to the satellite entry elevation angle threshold in step 1.

[0145] After the orbit prediction is completed, the space-based Internet of Things terminal will automatically wake up from sleep mode at the time of satellite entry according to the orbit prediction results and communicate with the satellite; instead of the satellite passing by, the terminal is in a low-power sleep state (it will also wake up automatically when data is collected).

[0146] The orbit prediction results based on the fitted prediction model are compared with the orbit prediction results of the TLE data at a fixed time (04:50:30 on December 29, 2020), as shown in Table 1:

[0147] Table 1

[0148]

[0149] The statistical forecast results show that in the absence of ephemeris updates, orbit prediction based on the fitting prediction model has a much higher accuracy within a certain period than orbit prediction based on simply solidifying the TLE data at a certain moment.

[0150] The above describes the technical solution of the present application in detail in combination with the accompanying drawings. The present application proposes a method for orbit prediction of a space-based Internet of Things terminal under the condition that the ephemeris is not updated, including: data initialization, and obtaining TLE data at a preset time; calling a two-row root number fitting prediction model, and calculating two-row root number prediction values ​​based on the prediction model and preset TLE data; and calculating the entry and exit time of the satellite according to the two-row root number prediction values. The technical solution in the present application helps to improve the accuracy of orbit prediction of space-based Internet of Things terminals under the condition that the ephemeris is not updated.

[0151] The steps in this application can be adjusted in order, combined, and deleted according to actual needs.

[0152] The units in the device of the present application can be combined, divided and deleted according to actual needs.

[0153] Although the present application is disclosed in detail with reference to the accompanying drawings, it should be understood that these descriptions are merely exemplary and are not intended to limit the application of the present application. The scope of protection of the present application is defined by the appended claims and may include various modifications, alterations and equivalents made to the invention without departing from the scope and spirit of the present application.

Claims

1. A method for orbit prediction of a space-based Internet of Things terminal under the condition that the ephemeris is not updated, characterized in that: The method includes Data initialization, and obtaining TLE data at the preset time; Calling a two-row root number fitting prediction model, and calculating two-row root number prediction values ​​based on the prediction model and the preset TLE data; Calculate the satellite's entry and exit times based on the two rows of predicted root values; The two rows of predicted values ​​of the root numbers include: the mean anomaly value and the predicted value of the argument of perigee; Calculating the mean anomaly angle value specifically includes: calculating the mean anomaly angle value at the forecast start time according to the established mean anomaly angle fitting prediction model, and the calculation formula of the mean anomaly angle value is: m1=(y0+((tx-t0)*n1*360 / (24*3600)))%360.0 Where m1 is the mean anomaly, y0 is the mean anomaly at the reference point, tx-t0 is the number of seconds from the forecast start time to the reference point, n1 is the average number of circles around the earth per day, and % is the modulo operation; Calculating the predicted value of the perigee argument specifically includes: calculating the predicted value of the perigee argument at the forecast start time according to the established perigee argument fitting prediction model, and the calculation formula of the predicted value of the perigee argument is: ω1=(p21*t+p22)%360.0 Wherein, ω1 is the calculated process value, ω2 is the predicted value of the perigee argument, t is the number of seconds from the forecast start time to the earliest sample historical TLE, and p21 and p22 are the fitting parameters of the second satellite.

2. The orbit prediction method of a space-based Internet of Things terminal under the condition that the ephemeris is not updated as claimed in claim 1, characterized in that: The data initialization specifically includes: Initialize the forecast start time, forecast days, terminal longitude and latitude information, and satellite entry elevation angle threshold.

3. The orbit prediction method of a space-based Internet of Things terminal under the condition that the ephemeris is not updated as claimed in claim 2, characterized in that: The step of calculating the entry and exit time of the satellite according to the two rows of predicted root values ​​specifically includes: Substituting the two rows of root number prediction values ​​into the SGP4 satellite orbit prediction model, the satellite position vector of the satellite in the Earth-centered Earth-fixed coordinate system from the prediction start time to the end time is calculated; Calculate the terminal position vector from the forecast start time to the end time in the Earth-centered Earth-fixed coordinate system according to the terminal latitude and longitude information; According to the satellite position vector and the terminal position vector, the elevation angle from the terminal to the satellite from the forecast start time to the end time is calculated, and when it is determined that the elevation angle is greater than or equal to the satellite entry elevation angle threshold, the entry and exit times of the satellite are calculated.

4. The orbit prediction method of a space-based Internet of Things terminal under the condition that the ephemeris is not updated as claimed in claim 1, characterized in that: The average value n1 of the number of orbits around the earth per day is the average value obtained by calculating the average value of the predicted number of orbits around the earth per day n2 of the satellite from the reference point time to the forecast start time at intervals of days, wherein the calculation formula for the predicted number of orbits around the earth per day n2 of the satellite is: n2=p11*t2^3+p12*t2^2+p13*t2+p14 Where n2 is the number of times the satellite orbits the earth every day, t2 is the number of seconds from the current time to the earliest sample historical TLE, and p11, p12, p13, and p14 are the first satellite fitting parameters.

5. The orbit prediction method of a space-based Internet of Things terminal under the condition that the ephemeris is not updated as claimed in claim 1, characterized in that: The method further includes: correcting the predicted value of the argument of perigee, and calculating the entry and exit time of the satellite according to the corrected predicted value of the argument of perigee, wherein the calculation formula for correcting the predicted value of the argument of perigee is: Wherein, ω3 is the predicted value of the corrected perigee argument angle, % is the modulo operation, x1, x2, x3 are the corrected values, ω2_1, ω2_2, ω2_3, ω2_4 are the corrected angle interval values.

6. The orbit prediction method of a space-based Internet of Things terminal under the condition that the ephemeris is not updated as claimed in any one of claims 1 to 3, characterized in that: The two rows of predicted values ​​of the elements include: orbital eccentricity value, and the two rows of predicted values ​​of the elements are calculated, specifically including: According to the established orbital eccentricity fitting prediction model, the orbital eccentricity value at the forecast start time is calculated. The calculation formula of the orbital eccentricity value is: ecc=a0+a1*cos(w*t)+b1*sin(w*t)-c1*t Where ecc is the orbit eccentricity, t is the number of seconds from the forecast start time to the earliest sample historical TLE, and a0, a1, b1, c1, and w are the eccentricity fitting parameters.

7. The orbit prediction method of a space-based Internet of Things terminal under the condition that the ephemeris is not updated as claimed in any one of claims 1 to 3, characterized in that: The two rows of predicted values ​​of the elements include: the right ascension value of the ascending node, and the calculation of the two rows of predicted values ​​of the elements specifically includes: According to the established ascending node right ascension fitting prediction model, the right ascension value of the ascending node at the forecast starting time is calculated. The calculation formula of the right ascension value of the ascending node is: Ω=(p31*t+p32)%360.0 Where Ω is the right ascension of the ascending node, t is the number of seconds from the forecast start time to the earliest sample historical TLE, and p31 and p32 are the fitting parameters of the third satellite.

8. The orbit prediction method of a space-based Internet of Things terminal under the condition that the ephemeris is not updated as claimed in any one of claims 1 to 3, characterized in that: The two rows of predicted values ​​of the number of elements include: orbital inclination value, and the calculation of the two rows of predicted values ​​of the number of elements specifically includes: According to the established orbital inclination fitting prediction model, the orbital inclination value at the forecast start time is calculated, and the calculation formula of the orbital inclination value is: x=t%31536000 y0=i0-Δi*(t / 31536000) i=p41*x^9+p42*x^8+p43*x^7+p44*x^6+p45*x^5+p46*x^4 +p47*x^3+p48*x^2+p49*x+y0 Where x is the number of seconds used for calculation, i is the orbital inclination, t is the number of seconds from the forecast start time to the earliest sample historical TLE, y0 is the initial angle of the orbital inclination with a periodicity of one year, % is the modulo operation, i0 is the orbital inclination value corresponding to the earliest sample TLE data, Δi is the total change value of the orbital inclination corresponding to the sample TLE data in one year, and p41, p42, p43, p44, p45, p46, p47, p48, and p49 are the fitting parameters of the fourth satellite.

Images