A fast visibility prediction algorithm for multi-beam communication satellites

Abstract

The present invention proposes a fast visibility prediction algorithm for multi-beam communication satellites. The method uses a precise orbit predictor suitable for on-orbit orbit prediction of low-orbit satellites, and improves the orbit prediction accuracy within one day to the sub-beam visibility accuracy of 0.1s. At the same time, a terminal critical communication condition solution method based on the ground inclination angle of the sub-satellite point trajectory is adopted. This method meets the precise search range considering the J2 average perturbation, and avoids the invalid solution process caused by the use of excessively large empirical proportional coefficients in engineering to expand the terminal critical communication conditions. Finally, a multi-beam visibility judgment function based on pitch angle and azimuth angle is proposed for communication satellites, and a polynomial fitting method is used to avoid iterative solution, thereby speeding up the calculation speed without affecting the accuracy.

Description

Technical Field

[0001] The present invention provides a high-precision fast prediction algorithm for the visibility of a multi-beam communication satellite, which relates to a fast orbit prediction model suitable for ground terminals, a method for determining the visible range based on the ground inclination angle of the sub-satellite point trajectory, and a method for determining the multi-beam visibility based on piecewise linear interpolation. The invention belongs to the field of satellite communication technology. Background Art

[0002] Low-orbit satellites generally refer to satellites with orbital altitudes of 200km-2000km. Low-orbit satellites are mainly used for ground target reconnaissance and ground terminal communications. For communication satellites, the lower orbital altitude gives them the advantages of short transmission delay and small path loss. Low-orbit satellite constellations generally refer to large satellite systems composed of multiple low-orbit satellites that can cooperate to complete specific tasks. Low-orbit satellite constellations are currently the most promising satellite mobile communication systems. High-precision visibility forecasting is one of the key technologies for establishing communication between satellites and ground terminals. Ground terminals that integrate high-precision visibility forecasting algorithms establish long-term, stable data transmission and communication cooperation tasks with satellites in unmanned conditions, and achieve the advantages of low power consumption and long endurance through standby during non-communication time.

[0003] There have been many studies on the visibility forecast of ground imaging, but these studies mainly focus on the modeling of the imaging area of ​​the space camera, and the time accuracy of the forecast reaches the second level. There is no relevant research on the visibility forecast of the sub-beam communication of communication satellites. In addition, for low-orbit satellites, the communication time of a single beam is only 10s-200s. In order to make full use of communication resources, higher requirements are put forward for the forecast accuracy of multi-beam communication. Since the current visibility forecast model does not establish a high-precision rapid orbit forecast model for low-orbit satellites, its visibility forecast accuracy is limited.

[0004] In summary, based on the huge research and application value of low-orbit satellite communication systems, the present invention proposes a high-precision and fast visibility prediction algorithm for multi-beam communication satellites. Summary of the invention

[0005] (I) Purpose of the invention

[0006] The present invention first derives the orbit prediction algorithm relative to the Earth-centered Earth-fixed system, then introduces the visible range determination method based on the ground inclination angle of the sub-satellite point trajectory, and further proposes a multi-beam visibility determination method based on piecewise linear interpolation. The purpose is to propose a high-precision and fast visibility prediction algorithm for multi-beam communication satellites, and provide technical support for the construction of intelligent satellite communication terminals.

[0007] (II) Technical solution

[0008] The high-precision and fast visibility prediction algorithm for multi-beam communication satellites of the present invention is implemented in the following steps:

[0009] Step 1: Preparation

[0010] First, the high-precision orbit prediction algorithm of the satellite relative to the Earth-centered Earth-fixed system is derived:

[0011] High-precision orbit prediction usually integrates the space position velocity in an inertial system, and requires the use of very complex mechanical models, including high-order gravity field models, three-body gravity models, tidal gravity models, atmospheric drag, solar light pressure, etc. Due to the existence of model errors, the accuracy of orbit prediction will become worse and worse as the prediction time increases, which will affect the accuracy of the over-the-top time and communication beam prediction. If a complex dynamic model is used for orbit prediction, although the orbit prediction error can be reduced, it will take up a large amount of computing resources. At present, many autonomous orbit predictors on satellites use the simplified conventional perturbation model (SGP4) developed by NORAD. The position accuracy of the commonly used SGP4 model prediction is 10km, and the visibility prediction error for low-orbit satellites reaches 1-2s. For higher prediction time accuracy requirements, the present invention uses the orbit perturbation Cowell method as a predictor, and provides the orbit dynamics equation under the earth-centered earth-fixed system.

[0012] Differential equations of motion for low-orbit communication satellites:

[0013]

[0014] In formula (1), r is the satellite position, and f is the central gravity and various perturbations acting on the satellite.

[0015] The low-orbit communication satellite is in the Earth-centered Earth-fixed coordinate system S e Differential equations of motion:

[0016]

[0017] In formula (2), the superscript e represents the vector in the Earth-centered Earth-fixed coordinate system S e The component array below. e , are the satellite positions, relative to S e Speed, relative to S e acceleration.

[0018] Analyze the impact of different types of perturbation acceleration on orbit prediction, and determine the most suitable rapid orbit prediction model for low-orbit satellites: for low-orbit satellites with orbit altitudes greater than 500km, a 9*9 order gravity field model can be used, and the orbit prediction position error within 24 hours is less than 1km.

[0019] The 9×9 order Earth gravitational potential described by the normalized spherical harmonic coefficients is:

[0020]

[0021] In formula (3), r, φ, λ represent the distance from the center of the earth, the latitude of the center of the earth, and the longitude of the center of the earth, respectively; μ is the gravitational constant of the earth, and R E is the Earth's equatorial radius, and is the normalized nth-order and mth-order spherical harmonic coefficient, is the normalized nth order and mth degree associated Legendre polynomial, the normalized coefficient

[0022] According to the recursive relationship of Legendre polynomials and the addition theorem of trigonometric functions, the recursive terms in the above earth gravitational potential are combined (Cunningham), and the earth gravitational potential model can be written as:

[0023]

[0024] In formula (4),

[0025] Gravitational acceleration of low-orbit satellites:

[0026]

[0027] Since the reference coordinate system of the Earth gravitational potential spherical harmonic model is the Earth-centered Earth-fixed system, the above Earth gravitational acceleration is actually the projection of the absolute acceleration of the satellite in the Earth-centered Earth-fixed coordinate system. Therefore, the relative acceleration relative to ECEF used for the numerical integration of orbital dynamics can be expressed as:

[0028]

[0029] In formula (6), the superscript e represents the component array in the Earth-centered Earth-fixed coordinate system, is the absolute acceleration, is the implicate acceleration, is the Coriolis acceleration.

[0030] The simplified Earth rotation model does not consider polar motion and nutation. It is believed that the Earth-centered Earth-fixed coordinate system is obtained by rotating the Earth-centered equatorial inertial system around the Z axis by a certain angle (Greenwich right ascension). Therefore, the calculation formulas for the entangled acceleration and the Coriolis acceleration can be written as:

[0031]

[0032] a c =2ω E ×v r =[-2ωE v y 2ω E v x z] (8)

[0033] In formula (7) and (8), ω E is the angular velocity of the Earth's rotation.

[0034] Then, a method for determining the critical communication conditions of the terminal based on the ground inclination of the sub-satellite point trajectory is introduced:

[0035] Assuming that the earth is a sphere and the low-orbit satellite is in a circular orbit, the communication range of the satellite is a cone, and its projection on the earth is a circle. The receiving range refers to the spherical radius of the projected circle, which can be expressed by the center angle of the sphere in a spherical triangle.

[0036] As the instruction manual Figure 1 As shown: S represents the location of the measuring station, M represents the location of the satellite, O represents the center of the earth, and N is the sub-satellite point of satellite M.

[0037] In the OSM triangle, by the law of sines:

[0038]

[0039] In formula (9), l OM =|r m |, l MS =|r s -r m |, l OS =R E +h, ∠OMS=η.

[0040] Radius of receiving range:

[0041]

[0042] In formula (10), η is the maximum pitch angle that meets the communication conditions in this system.

[0043] The best observation (communication) time refers to the point where the projection of the satellite on the local horizontal plane of the terminal is closest to the terminal. This point satisfies that the velocity direction of the sub-satellite point is perpendicular to the baseline direction.

[0044] As the instruction manual Figure 2 As shown: In the Earth-centered Earth-fixed coordinate system OXYZ, N represents the North Pole of the Earth, S represents the location of the terminal, N1 represents a point on the subsatellite point trajectory of the orbital plane P1, O1 represents the intersection of the latitude plane where N1 is located and the Z axis, and satisfies the spherical great circle SN1 = γ; the ground inclination angle at N1 is i E , and satisfies That is, the dihedral angle between plane O1N1N2 and the equatorial plane is iE ; Plane O1SN1 is perpendicular to plane O1N1N2. Therefore, N1 represents the critical communication position of the terminal shown in the figure.

[0045] According to the above geometric relationship, the angle between plane O1SN1 and the equatorial plane is π / 2-i E , we get the following system of equations:

[0046]

[0047] In formula (11), represents the normal vector of plane O1SN1, and j represents the unit vector of the Z axis.

[0048] Assume the longitude and latitude of S is (λ s ,φ s ), the longitude and latitude of N1 are (λ1, φ1), and according to the cosine theorem of spherical triangles, SN1 can be expressed as:

[0049] cos SN1=cos NS cos NN1+sin NS sin NN1 cos∠SNN1

[0050] Normal vector of plane O1SN1 in:

[0051]

[0052] S stands for sin and C stands for cos.

[0053] Substituting into formula (11) we can get:

[0054]

[0055] Solve the quartic equation about sinφ1 in equation (12) simultaneously:

[0056]

[0057] After solving, the latitude and longitude (λ1, φ1) of the sub-satellite point of the critical communication position are obtained.

[0058] Next, the calculation method of the ground inclination of the sub-satellite point trajectory is derived.

[0059] As the instruction manual Figure 3 As shown: In the spherical triangle BMT, B is the ascending node, M is the satellite position, and T is the intersection of the meridian plane through M with the equator. BM is a circular arc on the satellite orbit plane P{Ω}, angle B is the orbit inclination i, side BM is the latitude argument u of M, and side MT is the declination δ of M M , the edge BT is equal to α M -Ω M , where the right ascension of M is αM , the right ascension of the ascending node of the orbit is Ω M .

[0060] From the sine and cosine theorems of spherical right triangles:

[0061]

[0062] When satellite M is in the critical communication position, that is, corresponding to the attached Figure 2 When the subsatellite point N1 is in the middle, the right ascension and declination (α M , δ M ) has the following relationship with the longitude and latitude (λ1, φ1) of the sub-satellite point N1:

[0063]

[0064] In formula (14): α G (T) is the Greenwich right ascension at the corresponding moment of the current satellite position M.

[0065] When only the long-term influence of the Earth's oblateness is considered, only the right ascension of the ascending node of the satellite's orbital elements will change. For this type of orbit family P{Ω} (satellite orbits with the same semi-major axis, orbital inclination, eccentricity, and perigee argument), considering a simplified Earth model, it has the following properties: (1) the inclination of the ground track at the same declination of each orbit in the family is equal; (2) the time taken to fly from the same declination through the same latitude argument in each orbit in the family is equal.

[0066] Based on the orbit family property (2) above, we can use the orbit corresponding to the design value of the ascending node right ascension to calculate the ground inclination of the sub-satellite point track corresponding to the critical communication position N1 of the actual orbit.

[0067] Known δ M =φ1, let α′ M -Ω0=α M -Ω M , substituting into formula (14), we get:

[0068]

[0069] Considering the average perturbation effect of J2, we get:

[0070]

[0071] In formula (16), W Ω The average change rate of the ascending node of the circular orbit considering the J2 perturbation is: w u The average rate of change of the circular orbit latitude argument considering the J2 perturbation is:

[0072] Furthermore, the ground inclination angle of the sub-satellite point trajectory at the critical communication position N1 is obtained:

[0073]

[0074] By combining the above equations (12) and (17), the longitude and latitude (λ1, φ1) of the terminal critical communication position N1 and the corresponding ground track inclination angle i can be calculated: E .

[0075] Furthermore, the terminal critical communication condition based on the ground inclination of the sub-satellite point trajectory is derived from the property (2) of the above orbit family. The condition is expressed by the longitude range on the terminal latitude line that the satellite crosses.

[0076] As the instruction manual Figure 4 As shown in the figure: S represents the location of the terminal, M(N1) represents the critical communication location, M2 is the intersection of the orbital plane P1 and the terminal latitude, and N2 is the sub-satellite point corresponding to M2. The dotted line MM2 corresponds to the orbital plane of the satellite, and the solid line NN1 corresponds to the sub-satellite point trajectory of the orbital plane. The right ascension and declination of M2 and M are respectively (α M , δ M ), the longitude and latitude of N2 and N1 are (λ2, φ s ), (λ1, φ1).

[0077] Temporarily ignoring the non-spherical gravitational perturbation, we can obtain from equation (12):

[0078]

[0079] Formula (18) is solved to obtain:

[0080]

[0081] Considering the influence of the average perturbation of J2, the ascending node precesses in one cycle. Considering its average rate of change, λ2 is corrected to obtain:

[0082]

[0083] From the solution of formula (12), we can know that the terminal critical communication position N1 corresponds to two points, so N2 also has two corresponding points Therefore, the terminal critical communication condition based on the ground inclination of the sub-satellite point trajectory is the longitude range of the terminal latitude line crossed by the satellite Indicates that the visibility search continues for the satellite position when the condition is met.

[0084] Finally, the sub-beam visibility judgment function is constructed.

[0085] For communication satellites, it is generally believed that the satellite system is oriented relative to the satellite orbit system, that is, the X-axis points to the direction of satellite flight, the Z-axis points to the center of the earth, and the Y-axis points to the opposite direction of the orbital angular momentum. The beam definition is determined by the position of the transmitting antenna in the system, and the elevation angle and azimuth angle of the satellite-terminal baseline in the system can be used to determine the communication range of the terminal in a certain beam. At the same time, in a communication window, the satellite's position information has very weak nonlinearity, and the continuous function of the satellite position can be expressed by the piecewise interpolation fitting method, and the accurate visibility range can be obtained by further solving the judgment function of the sub-beam visibility.

[0086] It is known that satellite M satisfies the terminal critical communication condition based on the ground inclination of the sub-satellite point trajectory at time T, and its position in the Earth-centered Earth-fixed system is The orbit prediction model given above is used to predict the orbit for a period of time, and the pitch angle judgment condition is used to determine whether to terminate the orbit prediction. The prediction time interval is 10s, and a series of satellite positions are obtained.

[0087] For the above Using the quadratic polynomial fitting method, we can obtain the continuous function of the satellite M position with respect to time: Construct the full beam visibility judgment function based on the pitch angle information:

[0088]

[0089] In formula (20), r S It is the Earth-centered Earth-fixed coordinate of the terminal's location.

[0090] When F(t)>cosη, η is the maximum elevation angle that satisfies the communication conditions in this system. Let F(t)=cosη be rooted at (t1, t2), then the satellite full beam visibility range is T∈(t1, t2).

[0091] Furthermore, a sub-beam visibility judgment function is constructed according to the elevation angle and azimuth angle information corresponding to different beams.

[0092] The satellite-terminal baseline is a component array in the Earth-centered Earth-fixed system. The attitude transformation of the communication satellite's local system relative to the Earth-centered Earth-fixed system is a function of the position. Therefore, the component array of the satellite-terminal baseline in the local system can be written as:

[0093]

[0094] In formula (21), the superscript b represents the satellite system, and the superscript e represents the Earth-centered Earth-fixed system.

[0095] Furthermore, the azimuth information and elevation information of the satellite-terminal baseline can be expressed as:

[0096]

[0097] In formula (22), p(t) and a(t) represent the elevation angle and azimuth angle, respectively.

[0098] Sub-beam visibility judgment function:

[0099]

[0100] In formula (23), beam i (η) and beam i (ψ) correspond to the elevation and azimuth boundaries of the i-th sub-beam, respectively.

[0101] The sub-beam visibility range can be obtained by solving the root of equation (23), but this method involves solving the attitude transformation matrix and solving the root of the inverse trigonometric function. To further reduce the amount of calculation, the above sampling points can be Substituting into formula (22), we get a series of p i and a i , fit the new pitch angle and azimuth sampling points. Since the pitch angle and azimuth show high nonlinearity, the piecewise linear interpolation method is used to obtain its continuous function and Formula (23) can be rewritten as:

[0102]

[0103] Solving the root of equation (24) yields the visibility range of the communication satellite sub-beams.

[0104] At this point, the preparation work for the high-precision and fast visibility prediction algorithm for multi-beam communication satellites has been completed.

[0105] Step 2: Initialize orbital parameters of low-orbit communication satellite

[0106] At the current communication moment, the navigation receiver carried by the low-orbit satellite calculates the position and velocity of the satellite in the Earth-centered Earth-fixed system and broadcasts it to the ground terminal. The terminal completes the clock synchronization with the low-orbit communication satellite.

[0107] Step 3: Communication beam, visibility forecast duration, and terminal position initialization

[0108] According to the design parameters of the communication satellite, the elevation and azimuth angles of the multi-beam are initialized, and the total duration of the multi-beam visibility forecast is determined, which is set to 24 hours by default. In addition, the longitude and latitude of the terminal are obtained from the prior information.

[0109] Step 4: Calculate the satellite's osculating orbit elements

[0110] According to the position of the low-orbit communication satellite in the Earth-centered Earth-fixed system received by the terminal, the satellite's osculating orbit elements are calculated, and the orbital period and the angular velocity of the ascending node's right ascension precession considering the influence of J2 perturbation are further obtained, and the westward longitude of the satellite's sub-satellite point track in one orbital period is obtained. The terminal critical communication conditions based on the ground inclination of the sub-satellite point track are also calculated.

[0111] Step 5: Predict the satellite position in the first orbital period starting from the current time

[0112] Using the orbit prediction method that takes into account the influence of J2 perturbation, a series of satellite positions are obtained with a calculation step of 50s to 200s. The satellite coordinates in the Earth-centered Earth-fixed system are converted into the corresponding longitude and latitude of the sub-satellite point in the Earth coordinate system.

[0113] Step 6: Find the time corresponding to the satellite position closest to the terminal latitude in the first cycle

[0114] By obtaining a series of latitude information in step 5, find out the satellite position closest to the latitude of the terminal, and obtain the time corresponding to the position, which is recorded as the terminal latitude crossing time t0.

[0115] Step 7: Calculate the terminal latitude crossing time that meets the critical communication conditions based on the westward longitude

[0116] Based on the terminal latitude crossing time and westward longitude in steps 4 and 6, the terminal latitude crossing time t that meets the critical communication conditions within the visibility forecast duration is obtained.

[0117]

[0118] In formula (25), t0 represents the time when the terminal latitude crosses in the first cycle, T represents an orbital period, Δλ represents the difference between the current satellite longitude and the longitude corresponding to the critical communication condition, and λ W Indicates westward longitude, [x] Z Indicates rounding, whether it is rounding up or down depends on the east-west boundary of the critical communication condition.

[0119] Step 8: High-precision orbit prediction

[0120] According to the low-orbit satellite orbit parameters in step 2, use the high-precision orbit predictor derived in step 1, use an integration step of 100s, and solve the satellite position velocity at the time of terminal latitude crossing that meets the critical communication conditions.

[0121] Step 9: Short step sampling and judgment function fitting

[0122] Based on the satellite position velocity in step eight, a high-precision orbit predictor is used with a small step size (5s-20s) to integrate and calculate a series of satellite position velocities, and calculate the pitch angle and azimuth angle of the satellite-terminal baseline in the satellite system. The termination condition of the orbit prediction is set to the full-beam visibility judgment function F(t) in step one.

[0123] Furthermore, a series of satellite-terminal baselines are piecewise linearly fitted at the elevation angle and azimuth angle in the satellite system to obtain the continuous functions of the elevation angle and azimuth angle in step 1. and

[0124] Step 10: Beam visibility solution

[0125] Substituting the elevation and azimuth judgment conditions of the multiple communication beams in step 3 into the equation (23) described in step 1, the sub-beam visibility interval can be solved.

[0126] Through the above steps, a fast prediction algorithm for the visibility of multi-beam communication satellites is proposed. This method uses a precise orbit predictor suitable for on-orbit orbit prediction of low-orbit satellites, and improves the orbit prediction accuracy within one day to the sub-beam visibility accuracy of 0.1s. At the same time, a method for solving the terminal critical communication conditions based on the ground inclination of the sub-satellite point trajectory is adopted. This method meets the precise search range considering the J2 average perturbation, so that steps eight to ten can solve the terminal critical communication conditions more accurately, avoiding the invalid solution process caused by using too large empirical proportional coefficients in engineering to expand the terminal critical communication conditions. In addition, a multi-beam visibility judgment function based on pitch angle and azimuth angle is proposed for communication satellites, and a polynomial fitting method is used to avoid iterative solution, which speeds up the calculation speed without affecting the accuracy.

[0127] (III) Advantages

[0128] The advantages of the fast visibility prediction algorithm of a multi-beam communication satellite provided by the present invention are:

[0129] ① The algorithm proposed in the present invention involves a multi-beam visibility judgment function based on the pitch angle and the azimuth angle, which can give the visibility prediction results of each beam according to different communication states.

[0130] ② The algorithm proposed in the present invention involves an Earth-centered Earth-fixed orbit predictor, whose visibility forecast in the short term reaches 0.1s, which is 1-2 orders of magnitude higher than the visibility forecast result of the commonly used SGP4 orbit predictor. At the same time, the optimized orbit predictor can also adapt to the computing power of the terminal in terms of solution speed.

[0131] ③ The algorithm proposed in the present invention relates to a method for solving the critical communication conditions of the terminal based on the ground inclination of the sub-satellite point trajectory, which avoids the invalid solution process caused by using an excessively large empirical proportional coefficient to expand the critical communication conditions of the terminal in engineering, and reduces the amount of calculation. BRIEF DESCRIPTION OF THE DRAWINGS

[0132] Figure 1 It is a schematic diagram of solving the radius of the receiving range in the present invention.

[0133] Figure 2 , Figure 3 , Figure 4 It is a schematic diagram of solving the critical communication conditions of the terminal based on the ground inclination angle of the sub-satellite point trajectory.

[0134] Figure 5 It is a flow chart of the implementation steps of the present invention. DETAILED DESCRIPTION

[0135] The following will be combined with the attached Figure 5 The specific implementation process of the present invention is further described in detail with reference to the technical solution.

[0136] Step 1: Input parameter initialization

[0137] At the current communication moment, the navigation receiver carried by the low-orbit satellite calculates the position and velocity of the satellite in the Earth-centered Earth-fixed system and broadcasts it to the ground terminal. The terminal completes the clock synchronization with the low-orbit communication satellite.

[0138] According to the design parameters of the communication satellite, the elevation and azimuth angles of the multi-beam are initialized, and the total duration of the multi-beam visibility forecast is determined, which is set to 24 hours by default. In addition, the longitude and latitude of the terminal are obtained from the prior information.

[0139] This step corresponds to the attached Figure 1 The first box in .

[0140] Step 2: Calculate the satellite's orbital elements and related orbital parameters

[0141] According to the position of the low-orbit communication satellite in the Earth-centered Earth-fixed system received by the terminal, the satellite's osculating orbit elements are calculated, and the orbital period and the angular velocity of the ascending node's right ascension precession considering the influence of J2 perturbation are further obtained, and the westward longitude of the satellite's sub-satellite point track in one orbital period is obtained. The terminal critical communication conditions based on the ground inclination of the sub-satellite point track are also calculated.

[0142] This step corresponds to the attached Figure 1 The second box in .

[0143] Step 3: Predict the first cycle terminal latitude crossing time t0

[0144] Using the orbit prediction method that takes into account the influence of J2 perturbation, a series of satellite positions are obtained with a calculation step of 50s to 200s. The satellite coordinates in the Earth-centered Earth-fixed system are converted into the corresponding longitude and latitude of the sub-satellite point in the Earth coordinate system.

[0145] By obtaining a series of latitude information, the satellite position closest to the latitude of the terminal is found, and the time corresponding to the position is obtained, which is recorded as the terminal latitude crossing time t0.

[0146] This step corresponds to the attached Figure 1 The third box in .

[0147] Step 4: Calculate the terminal latitude crossing time that meets the critical communication conditions

[0148] The terminal latitude crossing time t that meets the critical communication condition is obtained from the terminal latitude crossing time and the westward longitude in steps 2 and 3.

[0149] This step corresponds to the attached Figure 1 The fourth box in .

[0150] Step 5: High-precision orbit prediction

[0151] According to the low-orbit satellite orbit parameters in step 2, the high-precision orbit predictor in the above technical solution is used, and the integration step is 100s to solve the satellite position speed at the time of terminal latitude crossing that meets the critical communication conditions.

[0152] This step corresponds to the attached Figure 1 The fifth box in .

[0153] Step 6: Baseline elevation and azimuth short-step sampling

[0154] Based on the satellite position velocity in step 5, a high-precision orbit predictor is used with a small step size (5s-20s) to integrate and calculate the satellite position velocity at the next moment, and calculate the pitch angle and azimuth angle of the satellite-terminal baseline in the satellite system.

[0155] This step corresponds to the attached Figure 1 The sixth box in .

[0156] Step 7: Termination conditions for precise orbit prediction

[0157] The termination condition of the orbit prediction is set to the full-beam visibility judgment function F(t)>cosη in the above technical solution.

[0158] This step corresponds to the attached Figure 1 The seventh box in .

[0159] Step 8: Determine function fit

[0160] Perform piecewise linear fitting on the elevation angle and azimuth angle of the satellite-terminal baseline generated in steps 6 and 7 in the satellite system to obtain the continuous function of the elevation angle and azimuth angle in the above technical solution. and

[0161] This step corresponds to the attached Figure 1 The eighth box in .

[0162] Step 9: Beam visibility solution

[0163] Substituting the elevation and azimuth boundary conditions of the multiple communication beams in step 1 into the equation described in formula (23) in the above technical solution, the sub-beam visibility interval can be solved.

[0164] This step corresponds to the attached Figure 1 The ninth box in .

[0165] In summary, the algorithm flow of the present invention is as shown in the attached Figure 5 As shown. A fast visibility prediction algorithm for multi-beam communication satellites is proposed. This method uses a precise orbit predictor suitable for on-orbit orbit prediction of low-orbit satellites, and improves the orbit prediction accuracy within a day to the sub-beam visibility accuracy of 0.1s. At the same time, a method for solving the terminal critical communication conditions based on the ground inclination of the sub-satellite point trajectory is adopted. This method meets the precise search range considering the J2 average perturbation and avoids the invalid solution process caused by the use of excessively large empirical proportional coefficients in engineering to expand the terminal critical communication conditions. Finally, a multi-beam visibility judgment function based on pitch angle and azimuth angle is proposed for communication satellites, and a polynomial fitting method is used to avoid iterative solution, which speeds up the calculation without affecting the accuracy.

Claims

1. A fast visibility prediction algorithm for a multi-beam communication satellite, characterized in that: The steps are as follows: Step 1: Preparation: A multi-beam visibility judgment function based on piecewise linear interpolation is given, a high-precision orbit prediction algorithm relative to the Earth-centered Earth-fixed system is derived, and the terminal critical communication condition based on the ground inclination angle of the sub-satellite point trajectory is derived; Step 2: Initialize the orbital parameters of the low-orbit communication satellite. At the current communication moment, calculate the position velocity of the low-orbit communication satellite in the Earth-centered Earth-fixed system according to the navigation receiver carried by the low-orbit communication satellite, further obtain the orbital parameters of the low-orbit communication satellite, and broadcast them to the terminal, so that the terminal completes the clock synchronization with the low-orbit communication satellite. Step 3: Initialize the communication beam, visibility forecast duration, and terminal position. According to the design parameters of the low-orbit communication satellite, initialize the elevation angle and azimuth angle boundaries of the multi-beam, and determine the total duration of the multi-beam visibility forecast, which is set to 24 hours by default. In addition, the longitude and latitude of the terminal are obtained from the prior information. Step 4: Calculate the osculating orbit elements of the low-orbit communication satellite. According to the position of the low-orbit communication satellite in the Earth-centered Earth-fixed system received by the terminal, calculate the osculating orbit elements of the low-orbit communication satellite, further obtain the orbital period and the angular velocity of the right ascension of the ascending node considering the influence of the second-order band harmonic J2 perturbation term, obtain the westward longitude of the sub-satellite point track of the low-orbit communication satellite in one orbital period, and calculate the terminal critical communication condition based on the ground inclination of the sub-satellite point track; Step 5: Predict the position of the low-orbit communication satellite in the first orbital period starting from the current time, use a high-precision orbit prediction algorithm that takes into account the influence of J2 perturbation, adopt an appropriate calculation step size, obtain a series of low-orbit communication satellite positions, and convert the coordinates of the low-orbit communication satellite in the Earth-centered Earth-fixed system into the corresponding longitude and latitude of the sub-satellite point in the Earth coordinate system; Step 6: Find the time corresponding to the position of the low-orbit communication satellite closest to the terminal latitude in the first orbital period, find the position of the low-orbit communication satellite closest to the terminal latitude through the longitude and latitude corresponding to the sub-satellite point obtained in step 5, and obtain the time corresponding to the position, which is recorded as the terminal latitude crossing time; Step 7: Calculate the terminal latitude crossing time that meets the critical communication condition according to the westward longitude. The terminal latitude crossing time that meets the critical communication condition within the visibility forecast duration is obtained by the terminal latitude crossing time and the westward longitude in steps 4 and 6. Step 8: high-precision orbit prediction, based on the low-orbit communication satellite orbit parameters in step 2, using the high-precision orbit prediction algorithm derived in step 1, using a suitable integration step size, to solve the low-orbit communication satellite position velocity at the time of terminal latitude crossing that meets the critical communication condition; Step nine, short-step sampling and judgment function fitting, using the low-orbit communication satellite position velocity in step eight, using a high-precision orbit prediction algorithm, using a small step size, integrating and calculating a series of low-orbit communication satellite position velocities, calculating the elevation angle and azimuth angle of the low-orbit communication satellite and the terminal baseline in the satellite system, further, performing piecewise linear fitting on the elevation angle and azimuth angle of a series of low-orbit communication satellites and the terminal baseline in the satellite system, and obtaining the continuous functions of the elevation angle and azimuth angle in the sub-beam visibility judgment function. and Step 10, solving the sub-beam visibility, substituting the elevation angle and azimuth angle boundaries of the multi-beam in step 3 into the sub-beam visibility judgment function, the sub-beam visibility interval can be solved.

Images