SGP4 model-based moon shadow time calculation method

Abstract

The invention discloses a moon shadow time calculation method based on an SGP4 model, and relates to the technical field of spacecraft orbit dynamics and energy management.The moon shadow time calculation method comprises the steps that a closed feedback loop is formed by extending a Kalman filter state variable and taking the instantaneous noise offset of an atmospheric resistance parameter B * as an independent estimation target; orbital perturbation changes caused by sun activity peaks, sudden geomagnetic disturbance and the like can be tracked in real time, and the stability of penumbra entering and exiting time of the near-earth satellite is remarkably improved; in order to solve the problem of frequent switching of shadow intervals of high-orbit satellites, a time period merging strategy based on an error ellipsoid and shadow area overlapping probability is provided. By quantizing the contribution of the position uncertainty to the shadow state, a pseudo switching event generated by model jitter is filtered out, a continuous shadow time period conforming to engineering practice is output, and invalid operation of an energy system is reduced.

Description

Technical Field

[0001] This invention relates to the field of spacecraft orbital dynamics and energy management technology, and in particular to a method for calculating lunar shadow time based on the SGP4 model. Background Technology

[0002] In the operation of satellites in near-Earth or near-lunar orbits at altitudes of 200-500 km, the ability to accurately predict the timing of entry into and exit from the Earth / Moon's penumbra directly impacts the security of the energy system. Satellites must switch to battery discharge before entering the shadow and promptly restore solar power after being exposed to sunlight again. Current mainstream solutions use TLE data to drive the SGP4 model to calculate the shadow boundary, supplemented by precise ephemeris corrections for orbital errors. The same spacecraft experiences more than ten penumbra entry and exit events every day, and the timing prediction accuracy needs to reach the second level to keep synchronized with power control commands.

[0003] While existing methods can improve orbit accuracy through real-time GNSS positioning or numerical integration, they still have limitations in low-Earth orbit environments with significant non-conservative perturbations (such as atmospheric drag and solar radiation pressure). The physical boundary between the penumbra and illumination itself has a transition zone, and SGP4 is extremely sensitive to empirical coefficients such as B*. When a satellite crosses the penumbra boundary, instantaneous fluctuations in atmospheric density or radiation pressure can induce orbital altitude drift, causing continuous jitter in the calculated entry and exit times. To address this, some schemes use sliding window filtering to smooth the results, but they cannot distinguish between actual orbital changes and noise, and second-level jumps still occur during peak solar activity periods.

[0004] Time jumps directly interfere with energy management: if the prediction of the penumbra start time is delayed, the battery may be over-discharged; if the prediction of the end time is too early, the solar panels may not take over the power supply in time, resulting in a power outage on the bus. Some systems mitigate the risk by increasing time redundancy buffers, but the short period of low Earth orbit means that conservative strategies will compress the payload's operating window. Although machine learning compensation methods have emerged in recent years, they rely on historical samples, and the risk of failure increases when encountering unseen scenarios such as sudden geomagnetic storms. Therefore, there is an urgent need for a lunar shadow time calculation method that still uses SGP4 as its core, can compensate for B* instantaneous deviations online, and is applicable to Earth / Moon shadow scenarios, in order to fundamentally eliminate second-level jumps and ensure the reliable operation of the spacecraft's power system. Summary of the Invention

[0005] In view of the aforementioned existing problems, the present invention is proposed.

[0006] This invention provides a method for calculating lunar shadow time based on the SGP4 model, which solves the problem that the traditional SGP4 model relies on static parameters and cannot handle errors in penumbra time calculation caused by atmospheric disturbances, shadow switching jitter, and lunar topography.

[0007] To solve the above-mentioned technical problems, the present invention provides the following technical solution:

[0008] This invention provides a method for calculating lunar shadow time based on the SGP4 model, which includes:

[0009] Step S1: Obtain the TLE orbital data and current UTC time of the space target;

[0010] Step S2: Calculate the target's orbital position sequence within a preset time window based on the SGP4 model;

[0011] Step S3: Perform dynamic noise compensation for drag parameters on the orbital position sequence.

[0012] Extract the non-conservative perturbation parameter B* from the TLE data to generate the historical error distribution;

[0013] The drag noise offset is estimated in real time using a Kalman filter.

[0014] The offset is injected into the orbital position coordinate correction;

[0015] Step S4: Calculate the shadow boundary crossing time based on the corrected orbital position and using the Earth or Moon as the occluder;

[0016] Step S5: Output the UTC timestamps for the start / end of the shadow.

[0017] As a preferred embodiment of the lunar shadow time calculation method based on the SGP4 model described in this invention, in step S3, during dynamic noise compensation of the drag parameter, the drag noise offset ΔB* is estimated in real time using a Kalman filter. This offset characterizes the instantaneous noise deviation of the non-conservative perturbation parameter B*. The B* parameter originates from TLE data, and its numerical change reflects the intensity of the non-conservative perturbation's influence on the orbit. The noise offset is estimated based on the residual between real-time GNSS positioning data and SGP4 prediction values, including:

[0018] Initialize state variables and define the Kalman filter state vector. :

[0019] ,

[0020] in, Indicates the discrete time step index, a dimensionless integer. This represents the position deviation components in the ECI coordinate system, in meters. This represents the velocity deviation component in the ECI coordinate system, in meters per second. This represents the drag noise offset, dimensionless, and represents the instantaneous deviation of the B* parameter, set during initialization. covariance matrix Constructed based on historical error distribution: the variance of the position error is set as... The variance of the speed error is set as , Variance set as ,in The historical standard deviation of the B* parameter for targets at the same orbital altitude;

[0021] State prediction based on the SGP4 orbital dynamics model:

[0022] ,

[0023] in, This represents the predicted state value at time k, with a dimension of 7×1. This represents the state transition matrix, with dimensions 7×7, and the structure is as follows:

[0024] ,

[0025] in, Indicates the time step, in seconds. Represents a 3×3 identity matrix. It is a zero matrix with dimensions 3×1. This represents the drag perturbation gradient vector, in units of s⁻², and is calculated as follows:

[0026] ,in represents the nonconservative perturbation acceleration in the SGP4 model, in m / s².

[0027] Covariance prediction:

[0028] ,in The process noise covariance matrix has a dimension of 7×7 and is calibrated based on historical orbital residuals.

[0029] Obtain the residual between GNSS positioning data and SGP4 prediction values:

[0030] ,

[0031] in, The three-axis position residuals in the ECI coordinate system are expressed in meters.

[0032] Observational update equation:

[0033] ,

[0034] ,

[0035] ,

[0036] in, The observation matrix has dimensions of 3×7 and a structure of... , This represents the observation noise covariance, with dimensions of 3×3, determined by the GNSS receiver accuracy parameters. Represents the Kalman gain matrix, dimension: 7×3, from the updated Extracting the seventh element Output as noise offset;

[0037] The estimated Injection trajectory position correction:

[0038] ,

[0039] in, This represents the original position vector of SGP4, in meters. This represents the position error component in the state vector, in meters. express Instantaneous deviation of parameters This represents the sensitivity matrix of position to B*, with dimensions 3×1 and units in meters. The calculation formula is:

[0040] ,in It is a tiny perturbation.

[0041] As a preferred embodiment of the lunar shadow time calculation method based on the SGP4 model described in this invention, step S4 includes:

[0042] Project the corrected orbit position onto the shadow geometry plane;

[0043] Calculate the ratio of the overlap area between the target position error ellipsoid and the shaded area;

[0044] When the overlap ratio of the umbra within a continuous time frame is ≥90%, it is merged into a single continuous shadow period;

[0045] The start / end point of the merged period is used as the shadow entry / exit time.

[0046] As a preferred embodiment of the lunar shadow time calculation method based on the SGP4 model described in this invention, in step S4, when the obstruction is the moon and terrain needs to be considered:

[0047] Load the pre-stored digital elevation model of the moon;

[0048] Perform terrain convolution offset correction on the shadow boundary coordinates: establish a terrain grid with the target-moon center line as the central axis, calculate the elevation standard deviation in the grid as the convolution kernel weight, and scale the terrain offset according to a preset ratio.

[0049] The offset is superimposed onto the coordinates of the base shadow boundary to obtain the corrected trajectory, which includes the error ellipsoid parameters.

[0050] As a preferred embodiment of the lunar shadow time calculation method based on the SGP4 model described in this invention, in step S4, when calculating the lunar shadow boundary, a terrain convolution offset correction is performed. This step is only enabled when the target's distance from the lunar center angle is <0.5°. Based on the pre-stored lunar digital elevation model (DEM), a terrain grid is established with the line connecting the target and the lunar center as the central axis, and the convolution kernel weights are calculated through the following process:

[0051] Establishing the terrain grid and loading elevation data:

[0052] A local coordinate system is established with the line connecting the target lunar centers as the z-axis, and a rectangular grid region is defined on a plane perpendicular to the z-axis; the grid size is fixed to the side length. Square cells; load DEM elevation data to grid nodes, and associate each grid cell with an elevation value. ,in Represents grid row and column indices; grid elevation matrix Represented as:

[0053] ,

[0054] in, Represents the number of rows and columns of the grid, a dimensionless integer, determined by the region size and cell side length. Sure, Represents grid cells The elevation value, in meters. This represents the side length of the grid cell, in meters. It is set as a constant by default and is determined based on the resolution of the lunar DEM.

[0055] Calculate the standard deviation of elevation within the grid:

[0056] For each grid cell Extract all elevation values ​​within its neighborhood window, window size Calculate the standard deviation of elevation. As a measure of local topographic relief; the standard deviation formula is:

[0057] ,

[0058] in, Represents grid cells Standard deviation of elevation, in meters. Represents the number of elevation points within the neighborhood window, a dimensionless integer. Indicates the position within the neighborhood window The elevation value, in meters. This represents the average elevation within the neighborhood window, in meters. The formula is:

[0059] ,

[0060] Indicates the size of the neighborhood window, a dimensionless odd number, defaulted to 3, covering 3×3 cells. and The row and column indexes for the grid neighborhood;

[0061] Define the convolution kernel weights:

[0062] Elevation standard deviation Normalized to convolution kernel weights This is used to weight the impact of terrain on the offset of the shadow boundary; the weighting formula is:

[0063] ,

[0064] in, Represents grid cells Convolution kernel weights, dimensionless scalars This represents the maximum standard deviation over the entire grid area, in meters. This represents the safety threshold constant, in meters, and is preset to [value missing]. rice;

[0065] weight matrix The construction is as follows:

[0066] ,

[0067] Correcting weight-driven convolution offsets:

[0068] Based on the weight matrix Perform convolution to calculate the shadow boundary offset caused by the terrain. The convolution formula is:

[0069] ,

[0070] in, This represents the terrain convolution offset, in meters. Represents grid cells The influence of geometric distance on the base shadow edge, in meters, is calculated from the DEM elevation gradient.

[0071] As a preferred embodiment of the lunar shadow time calculation method based on the SGP4 model described in this invention, the terrain convolution offset correction is only enabled when the target's distance from the lunar center angle is <0.5°.

[0072] As a preferred embodiment of the lunar shadow time calculation method based on the SGP4 model described in this invention, wherein: in the noise compensation step S3:

[0073] The historical error distribution was constructed using the standard deviation of the B* parameter of targets at the same orbital altitude for a preset number of days in the past.

[0074] The observations from the Kalman filter are taken from the residuals between real-time GNSS positioning data and SGP4 predictions.

[0075] As a preferred embodiment of the lunar shadow time calculation method based on the SGP4 model described in this invention, the calculation of the overlap area ratio includes:

[0076] Establish an error ellipsoid centered on the target position;

[0077] The ellipsoid is discretized into multiple sampling points, and the number of sampling points is preset to a fixed value;

[0078] The proportion of sampling points falling into the umbra or penumbra is statistically analyzed.

[0079] As a preferred embodiment of the lunar shadow time calculation method based on the SGP4 model described in this invention, the output timestamp in step S5 is truncated to UTC integer seconds.

[0080] As a preferred embodiment of the lunar shadow time calculation method based on the SGP4 model described in this invention, the Kalman filter in step S3 includes a state variable definition process, wherein the state variables include position error and velocity error components.

[0081] The beneficial effects of this invention are as follows: By expanding the state variables of the Kalman filter, this invention takes the instantaneous noise offset of the atmospheric drag parameter B* as an independent estimation target, forming a closed feedback loop; compared with traditional static filtering or interpolation methods, this mechanism can track orbital perturbation changes caused by solar activity peaks, sudden geomagnetic disturbances, etc. in real time, significantly improving the stability of near-Earth satellite penumbra entry and exit times.

[0082] This invention addresses the problem of frequent switching of shadow regions for high-orbit satellites by proposing a time-segment merging strategy based on the overlap probability between the error ellipsoid and the shadow region. By quantifying the contribution of positional uncertainty to the shadow state, it filters out spurious switching events caused by model jitter, outputting continuous shadow periods that conform to engineering realities, and reducing ineffective operations in the energy system.

[0083] This invention pioneers a terrain convolution kernel weighting mechanism in lunar shadow calculation. By dynamically adjusting the offset weight of terrain on the shadow boundary through the local elevation standard deviation, the correction amount adapts to lunar surface terrain features (such as steep changes at crater edges or flat areas in lunar maria). Without significantly increasing computational load, this invention effectively compensates for the optical path difference between the ideal spherical model and the real terrain, making it particularly suitable for scenarios requiring sub-second precision, such as lunar obstacle avoidance by landers.

[0084] This invention is compatible with existing space data infrastructure (TLE, GNSS, DEM) and requires no additional hardware. By deploying lightweight models (such as Kalman filtering and convolution kernel calculation) at the edge computing layer, it reduces reliance on cloud computing power and meets the real-time requirements of onboard processors. Attached Figure Description

[0085] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the following description of the embodiments will be briefly introduced. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0086] Figure 1 This is a flowchart illustrating the lunar shadow time calculation method based on the SGP4 model in Example 1. Detailed Implementation

[0087] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

[0088] Many specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways different from those described herein, and those skilled in the art can make similar extensions without departing from the spirit of the invention. Therefore, the invention is not limited to the specific embodiments disclosed below.

[0089] Secondly, the term "one embodiment" or "embodiment" as used herein refers to a specific feature, structure, or characteristic that may be included in at least one implementation of the present invention. The phrase "in one embodiment" appearing in different places in this specification does not necessarily refer to the same embodiment, nor is it a single or selective embodiment that is mutually exclusive with other embodiments.

[0090] Example 1, referring to Figure 1 This embodiment provides a method for calculating lunar shadow time based on the SGP4 model, extending TLE for near-lunar orbit, including the following steps:

[0091] Step S1: Obtain the TLE orbital data and current UTC time of the space target;

[0092] Step S2: Calculate the target's orbital position sequence within a preset time window based on the SGP4 model;

[0093] Step S3: Perform dynamic noise compensation for drag parameters on the orbital position sequence.

[0094] Extract the non-conservative perturbation parameter B* from the TLE data to generate the historical error distribution; the B* parameter is a scalar in the TLE data that characterizes the intensity of non-conservative perturbations (including solar radiation pressure and atmospheric drag). An increase in its value indicates that the drag has a stronger impact on the orbit, and it is related to the target mass, cross-sectional area and atmospheric density; (in near-Earth orbit, it mainly characterizes atmospheric drag, and in near-lunar orbit, it can be extended to the empirical coefficient of non-conservative perturbations such as solar radiation pressure).

[0095] The drag noise offset is estimated in real time using a Kalman filter.

[0096] The offset is injected into the orbital position coordinate correction;

[0097] In step S3, noise compensation:

[0098] The historical error distribution is constructed using the standard deviation of the B* parameter of targets at the same orbital altitude for a preset number of days (corresponding to targets with the same gravitational field of the obstructing body);

[0099] The observations from the Kalman filter are taken from the residuals between real-time GNSS positioning data and SGP4 predictions;

[0100] The Kalman filter in step S3 includes a state variable definition process, wherein the state variables include position error and velocity error components;

[0101] In step S3, during dynamic noise compensation of the drag parameter, the drag noise offset ΔB* is estimated in real time using a Kalman filter. This offset characterizes the instantaneous noise deviation of the non-conservative perturbation parameter B*. The B* parameter originates from TLE data, and its numerical change reflects the intensity of the non-conservative perturbation's influence on the orbit. The noise offset is estimated based on the residual between real-time GNSS positioning data and SGP4 prediction values, including:

[0102] Initialize state variables and define the Kalman filter state vector. :

[0103] ,

[0104] in, Indicates the discrete time step index, a dimensionless integer. This represents the position deviation components in the ECI coordinate system, in meters. This represents the velocity deviation component in the ECI coordinate system, in meters per second. This represents the drag noise offset, dimensionless, and represents the instantaneous deviation of the B* parameter, set during initialization. covariance matrix Constructed based on historical error distribution: the variance of the position error is set as... The variance of the speed error is set as , Variance set as ,in The historical standard deviation of the B* parameter for targets at the same orbital altitude;

[0105] State prediction based on the SGP4 orbital dynamics model:

[0106] ,

[0107] in, This represents the predicted state value at time k, with a dimension of 7×1. This represents the state transition matrix, with dimensions 7×7, and the structure is as follows:

[0108] ,

[0109] in, Indicates the time step, in seconds. Represents a 3×3 identity matrix. It is a zero matrix with dimensions 3×1. This represents the drag perturbation gradient vector, in units of s⁻², and is calculated as follows:

[0110] ,in represents the nonconservative perturbation acceleration in the SGP4 model, in m / s².

[0111] Covariance prediction:

[0112] ,in The process noise covariance matrix has a dimension of 7×7 and is calculated based on the statistical characteristics of orbital residuals of targets at the same orbital altitude over the past 30 days, calibrated according to historical orbital residuals.

[0113] Obtain the residual between GNSS positioning data and SGP4 prediction values:

[0114] ,

[0115] in, The three-axis position residuals in the ECI coordinate system are expressed in meters.

[0116] Observational update equation:

[0117] ,

[0118] ,

[0119] ,

[0120] in, The observation matrix has dimensions of 3×7 and a structure of... , This represents the observation noise covariance, with dimensions of 3×3, determined by the GNSS receiver accuracy parameters. Represents the Kalman gain matrix, dimension: 7×3, from the updated Extracting the seventh element Output as noise offset;

[0121] The estimated Injection trajectory position correction:

[0122] ,

[0123] in, This represents the original position vector of SGP4, in meters. This represents the position error component in the state vector, in meters. express Instantaneous deviation of parameters This represents the sensitivity matrix of position to B*, with dimensions 3×1 and units in meters. The calculation formula is:

[0124] ,in For minute perturbations;

[0125] Specifically, this Kalman filter design expands the state variables to directly estimate drag noise offset, establishing a quantitative relationship between non-conservative perturbations and orbital position errors. The drag gradient term in the state transition matrix couples the physical model with noise characteristics, while observation updates utilize real-time GNSS data to correct system drift. In low Earth orbit scenarios, this method effectively tracks instantaneous changes in B* parameters caused by solar activity or sudden atmospheric disturbances, suppressing second-level jumps in penumbra entry and exit time calculations. The filter output is directly used for position correction, forming a closed-loop compensation mechanism, which is more adaptable to complex space environments compared to traditional open-loop correction.

[0126] Step S4: Calculate the shadow boundary crossing time based on the corrected orbital position and using the Earth or Moon as the occluder;

[0127] Step S4 includes:

[0128] Project the corrected orbit position onto the shadow geometry plane;

[0129] Calculate the ratio of the overlap area between the target position error ellipsoid and the shaded area;

[0130] When the overlap ratio of the umbra within a continuous time frame is ≥90%, it is merged into a single continuous shadow period;

[0131] The start / end point of the merged time period is used as the shadow entry / exit time;

[0132] The calculation of the overlap area ratio includes:

[0133] Establish an error ellipsoid centered on the target position;

[0134] The ellipsoid is discretized into multiple sampling points, and the number of sampling points is preset to a fixed value;

[0135] Statistical analysis of the proportion of sampling points falling within the umbra or penumbra region;

[0136] In step S4, when the obstruction is the moon and terrain needs to be considered:

[0137] Load the pre-stored digital elevation model of the moon;

[0138] Perform terrain convolution offset correction on the shadow boundary coordinates: establish a terrain grid with the target-moon center line as the central axis, calculate the elevation standard deviation in the grid as the convolution kernel weight, and scale the terrain offset according to a preset ratio.

[0139] The offset is superimposed on the coordinates of the base shaded boundary to obtain the corrected trajectory that includes the error ellipsoid parameters;

[0140] In step S4, when calculating the lunar shadow boundary, a terrain convolution offset correction is performed. This step is only enabled when the target's distance from the lunar center angle is <0.5°. Based on the pre-stored lunar digital elevation model (DEM), a terrain mesh is established with the line connecting the target and the lunar center as the central axis. The convolution kernel weights are calculated through the following process:

[0141] Establishing the terrain grid and loading elevation data:

[0142] A local coordinate system is established with the line connecting the target lunar centers as the z-axis, and a rectangular grid region is defined on a plane perpendicular to the z-axis; the grid size is fixed to the side length. Square cells are used, with the grid covering a 5km × 5km area centered on the lunar projection point (based on lunar DEM data resolution optimization); DEM elevation data is loaded into the grid nodes, and each grid cell is associated with an elevation value. ,in Represents grid row and column indices; grid elevation matrix Represented as:

[0143] ,

[0144] in, Represents the number of rows and columns of the grid, a dimensionless integer, determined by the region size and cell side length. Sure, Represents grid cells The elevation value, in meters. This indicates the side length of the grid cell in meters. It is set as a constant by default and is based on the resolution of the lunar DEM. This step ensures that the grid covers the shadow boundary influence area (the orthophoto area centered on the lunar center projection point), and the central axis alignment reduces coordinate transformation errors.

[0145] Calculate the standard deviation of elevation within the grid:

[0146] For each grid cell Extract all elevation values ​​within its neighborhood window, window size Calculate the standard deviation of elevation. As a measure of local topographic relief; the standard deviation formula is:

[0147] ,

[0148] in, Represents grid cells Standard deviation of elevation, in meters. Represents the number of elevation points within the neighborhood window, a dimensionless integer. Indicates the position within the neighborhood window The elevation value, in meters. This represents the average elevation within the neighborhood window, in meters. The formula is:

[0149] ,

[0150] Indicates the size of the neighborhood window, a dimensionless odd number, defaulted to 3, covering 3×3 cells. and The grid neighborhood is indexed by rows and columns, and the standard deviation is calculated to reflect the terrain roughness. The value corresponds to significant topographical changes (such as mountains or impact craters).

[0151] Define the convolution kernel weights:

[0152] Elevation standard deviation Normalized to convolution kernel weights This is used to weight the impact of terrain on the offset of the shadow boundary; the weighting formula is:

[0153] ,

[0154] in, Represents grid cells Convolution kernel weights, dimensionless scalars This represents the maximum standard deviation over the entire grid area, in meters. This represents the safety threshold constant, in meters, and is preset to [value missing]. rice;

[0155] weight matrix The construction is as follows:

[0156] ,

[0157] This weight is directly used in the convolution operation to scale the terrain offset and enhance the contribution of high-undulation areas to the shadow boundary.

[0158] Correcting weight-driven convolution offsets:

[0159] Based on the weight matrix Perform convolution to calculate the shadow boundary offset caused by the terrain. The convolution formula is:

[0160] ,

[0161] in, This represents the terrain convolution offset, in meters. Represents grid cells The geometric distance effect on the edge of the base shadow (equivalent distance of optical path difference caused by terrain undulation (along the target-sun line)) is calculated in meters by the DEM elevation gradient and finally superimposed on the coordinates of the base shadow boundary to complete the correction.

[0162] Specifically, this step dynamically generates convolution kernel weights using the elevation standard deviation to address the non-ideal deviation of shadow boundaries caused by lunar surface topographic undulations (such as mountains or canyons). The standard deviation calculation captures local terrain roughness, and the normalized weights ensure that high-undulation areas dominate the correction, avoiding noise introduced by smooth terrain. The convolution operation couples the physical characteristics of the terrain with the shadow geometry model, effectively correcting second-level time errors caused by light diffraction effects in scenarios with a small distance from the lunar center (such as the lander's approach phase). The weighting mechanism adapts to different lunar regions (such as highlands and lunar maria), reducing reliance on fixed parameters. Only standard DEM data is required during implementation, improving the overall consistency of shadow entry and exit times and supporting high-precision mission planning.

[0163] Terrain convolution offset correction is only enabled when the target's distance from the lunar center angle is less than 0.5°;

[0164] Step S5: Output the UTC timestamps for the start / end of the shadow; the output timestamps in step S5 are truncated to UTC integer seconds.

[0165] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A method for calculating lunar shadow time based on the SGP4 model, characterized in that, include, Step S1: Obtain the TLE orbital data and current UTC time of the space target; Step S2: Calculate the target's orbital position sequence within a preset time window based on the SGP4 model; Step S3: Perform dynamic noise compensation for drag parameters on the orbital position sequence. Extract the non-conservative perturbation parameter B* from the TLE data to generate the historical error distribution; The drag noise offset is estimated in real time using a Kalman filter. The offset is injected into the orbital position coordinate correction; Step S4: Calculate the shadow boundary crossing time based on the corrected orbital position and using the Earth or Moon as the occluder; Step S5: Output the UTC timestamps for the start / end of the shadow.

2. The method for calculating lunar shadow time based on the SGP4 model as described in claim 1, characterized in that, In step S3, during dynamic noise compensation of the drag parameter, the drag noise offset ΔB* is estimated in real time using a Kalman filter. This offset characterizes the instantaneous noise deviation of the non-conservative perturbation parameter B*. The B* parameter originates from TLE data, and its numerical change reflects the intensity of the non-conservative perturbation's influence on the orbit. The noise offset is estimated based on the residual between real-time GNSS positioning data and SGP4 prediction values, including: Initialize state variables and define the Kalman filter state vector. : ， in, Indicates the discrete time step index, a dimensionless integer. This represents the position deviation components in the ECI coordinate system, in meters. This represents the velocity deviation component in the ECI coordinate system, in meters per second. This represents the drag noise offset, dimensionless, and represents the instantaneous deviation of the B* parameter, set during initialization. covariance matrix Constructed based on historical error distribution: The variance of the position error is set as... The variance of the speed error is set as , Variance set as ,in The historical standard deviation of the B* parameter for targets at the same orbital altitude; State prediction based on the SGP4 orbital dynamics model: ， in, This represents the predicted state value at time k, with a dimension of 7×1. This represents the state transition matrix, with dimensions 7×7, and the structure is as follows: ， in, Indicates the time step, in seconds. Represents a 3×3 identity matrix. It is a zero matrix with dimensions 3×1. This represents the drag perturbation gradient vector, in units of s⁻², and is calculated as follows: ,in represents the nonconservative perturbation acceleration in the SGP4 model, in m / s². Covariance prediction: ,in The process noise covariance matrix has a dimension of 7×7 and is calibrated based on historical orbital residuals. Obtain the residual between GNSS positioning data and SGP4 prediction values: ， in, The three-axis position residuals in the ECI coordinate system are expressed in meters. Observational update equation: ， ， ， in, The observation matrix has dimensions of 3×7 and a structure of... , This represents the observation noise covariance, with dimensions of 3×3, determined by the GNSS receiver accuracy parameters. Represents the Kalman gain matrix, dimension: 7×3, from the updated Extracting the seventh element Output as noise offset; The estimated Injection trajectory position correction: ， in, This represents the original position vector of SGP4, in meters. This represents the position error component in the state vector, in meters. express Instantaneous deviation of parameters This represents the sensitivity matrix of position to B*, with dimensions 3×1 and units in meters. The calculation formula is: ,in It is a tiny perturbation.

3. The method for calculating lunar shadow time based on the SGP4 model as described in claim 1, characterized in that, Step S4 includes: Project the corrected orbit position onto the shadow geometry plane; Calculate the ratio of the overlap area between the target position error ellipsoid and the shaded area; When the overlap ratio of the umbra within a continuous time frame is ≥90%, it is merged into a single continuous shadow period; The start / end point of the merged period is used as the shadow entry / exit time.

4. The method for calculating lunar shadow time based on the SGP4 model as described in claim 1, characterized in that, In step S4, when the obstruction is the moon and terrain needs to be considered: Load the pre-stored digital elevation model of the moon; Perform terrain convolution offset correction on the shadow boundary coordinates: establish a terrain grid with the target-moon center line as the central axis, calculate the elevation standard deviation in the grid as the convolution kernel weight, and scale the terrain offset according to a preset ratio. The offset is superimposed onto the coordinates of the base shadow boundary to obtain the corrected trajectory, which includes the error ellipsoid parameters.

5. The method for calculating lunar shadow time based on the SGP4 model as described in claim 4, characterized in that, In step S4, when calculating the lunar shadow boundary, a terrain convolution offset correction is performed. This step is only enabled when the target's distance from the lunar center angle is <0.5°. Based on the pre-stored lunar digital elevation model (DEM), a terrain mesh is established with the line connecting the target and the lunar center as the central axis. The convolution kernel weights are calculated through the following process: Establishing the terrain grid and loading elevation data: A local coordinate system is established with the line connecting the target lunar centers as the z-axis, and a rectangular grid region is defined on a plane perpendicular to the z-axis; the grid size is fixed to the side length. Square cells; load DEM elevation data to grid nodes, and associate each grid cell with an elevation value. ,in Represents grid row and column indices; grid elevation matrix Represented as: ， in, Represents the number of rows and columns of the grid, a dimensionless integer, determined by the region size and cell side length. Sure, Represents grid cells The elevation value, in meters. This represents the side length of the grid cell, in meters. It is set as a constant by default and is determined based on the resolution of the lunar DEM. Calculate the standard deviation of elevation within the grid: For each grid cell Extract all elevation values ​​within its neighborhood window, window size Calculate the standard deviation of elevation. As a measure of local topographic relief; the standard deviation formula is: ， in, Represents grid cells Standard deviation of elevation, in meters. Represents the number of elevation points within the neighborhood window, a dimensionless integer. Indicates the position within the neighborhood window The elevation value, in meters. This represents the average elevation within the neighborhood window, in meters. The formula is: ， Indicates the size of the neighborhood window, a dimensionless odd number, defaulted to 3, covering 3×3 cells. and The row and column indexes for the grid neighborhood; Define the convolution kernel weights: Elevation standard deviation Normalized to convolution kernel weights This is used to weight the impact of terrain on the offset of the shadow boundary; the weighting formula is: ， in, Represents grid cells Convolution kernel weights, dimensionless scalars This represents the maximum standard deviation over the entire grid area, in meters. This represents the safety threshold constant, in meters, and is preset to [value missing]. rice; weight matrix The construction is as follows: ， Correcting weight-driven convolution offsets: Based on the weight matrix Perform convolution to calculate the shadow boundary offset caused by the terrain. The convolution formula is: ， in, This represents the terrain convolution offset, in meters. Represents grid cells The influence of geometric distance on the edge of the base shadow, in meters, is calculated from the DEM elevation gradient.

6. The method for calculating lunar shadow time based on the SGP4 model as described in claim 4, characterized in that, The terrain convolution offset correction is only enabled when the target's angle from the lunar center is less than 0.5°.

7. The method for calculating lunar shadow time based on the SGP4 model as described in claim 1, characterized in that, In step S3, noise compensation: The historical error distribution was constructed using the standard deviation of the B* parameter of targets at the same orbital altitude for a preset number of days in the past. The observations from the Kalman filter are taken from the residuals between real-time GNSS positioning data and SGP4 predictions.

8. The method for calculating lunar shadow time based on the SGP4 model as described in claim 3, characterized in that, The calculation of the overlap area ratio includes: Establish an error ellipsoid centered on the target position; The ellipsoid is discretized into multiple sampling points, and the number of sampling points is preset to a fixed value; The proportion of sampling points falling into the umbra or penumbra is statistically analyzed.

9. The method for calculating lunar shadow time based on the SGP4 model as described in claim 1, characterized in that, The output timestamp of step S5 is truncated to an integer UTC second.

10. The method for calculating lunar shadow time based on the SGP4 model as described in claim 1, characterized in that, The Kalman filter in step S3 includes a state variable definition process, where the state variables include position error and velocity error components.

Images