Method and system for estimating minimum velocity increment of multi-circle double-pulse orbit transfer under J2 perturbation

Abstract

The invention particularly relates to a method and a system for estimating the minimum velocity increment of multi-circle double-pulse orbit transfer under J2 perturbation. The method comprises the following steps: uniformly representing origin-destination orbit states as six spring equinox elements taking flat longitude as a phase quantity; aiming at a transfer structure of a first pulse, middle drift and a second pulse, determining a limited transfer turn number search window by utilizing a flat longitude drift characteristic and a pulse capability; a closed equality constraint system with the total variable quantity = the first pulse variable quantity + the middle-section drift variable quantity + the second pulse variable quantity is constructed, and it is ensured that the track state is strictly closed; pulse is decomposed into tangential, normal and radial components for cooperative solution, and a phase closed loop is realized by updating a semi-diameter increment; and performing low-dimensional nonlinear optimization in a limited number of turns to quickly obtain minimum speed increment estimation. According to the method, the calculation efficiency is remarkably improved while the precision is guaranteed, and the method is suitable for large-scale sequence optimization in the near-earth orbit multi-target access task.

Description

Technical Field

[0001] This invention relates to the field of aerospace navigation and control technology, specifically to a method and system for estimating the minimum velocity increment of multi-cycle double-pulse orbit transfer under J2 perturbation. Background Technology

[0002] Low Earth orbit (LEO) multi-target rendezvous missions are characterized by "multiple targets, long time domains, and tight time windows": the servicing spacecraft must sequentially visit multiple targets within a given epoch, completing multiple transfers while minimizing the total velocity increment. The sequence order directly determines mission time consumption and fuel consumption, and upper-level sequence optimization often requires tens of thousands of evaluations of the transfer cost between "any two targets." If the transfer between targets is considered as a multi-pulse transfer, each evaluation is equivalent to solving a time-constrained parametric optimization problem; while numerical optimization can provide high-precision velocity increments, it introduces dense propagation and repeated iterations, and the computational burden increases exponentially with the number of targets and candidate time windows, leading to unacceptable time complexity in upper-level sequence optimization. Therefore, developing a method that can quickly and accurately estimate the velocity increment between two targets given origin and destination orbits and transfer durations is a key foundation for supporting large-scale sequence design and real-time mission decision-making.

[0003] Existing fast estimation techniques mainly fall into three categories: First, database-based methods. While online retrieval is rapid, their construction relies on massive amounts of optimized trajectory samples, resulting in extremely high database generation costs. Furthermore, the coverage dimensions are limited, and the ability to adapt to changing scenarios is weak, making it difficult to meet the demands of variable trajectory design. Second, machine learning methods. These methods achieve rapid inference based on neural networks, but they depend on large-scale labeled datasets, requiring training times of tens of hours. The models are also severely limited by the training scenario, with accuracy rapidly degrading when trajectory parameters deviate from the sample distribution, and exhibiting insufficient generalization ability. Third, analytical / semi-analytical methods directly establish a mapping relationship between "trajectory element difference—velocity increment" based on perturbation dynamics approximation. However, most representative works in this category are based on near-circular orbits or small-change assumptions, and their handling of long-term J2 perturbations and orbital plane-phase coupling effects is relatively coarse. This leads to significant phase deviation accumulation and orbital plane misjudgment when dealing with low-Earth orbit long-time domain, multi-turn orbits, and medium-to-high eccentricity orbital transfer conditions. This is due to neglecting the J2 perturbation or excessive near-circular simplification, resulting in a systematic underestimation of velocity increments or rendering the solution infeasible. Consequently, its applicability and robustness are severely limited. Therefore, current technologies still lack a Δ-type algorithm that can uniformly describe phase evolution, possess closed-loop balancing capabilities, and offer efficient search performance while maintaining high computational efficiency when addressing dual-pulse multi-turn orbital transfer problems across the entire eccentricity range under J2 perturbation. V In addition to fast estimation methods, there is an urgent need to develop next-generation estimation strategies to support the practical engineering needs of large-scale on-orbit mission sequence optimization.

[0004] It should be noted that the information disclosed in the background section above is only used to enhance the understanding of the background of the present invention, and therefore may include information that does not constitute prior art known to those skilled in the art. Summary of the Invention

[0005] For existing orbital transfer Δ V The estimation method is effective in handling orbits with large eccentricities and multiple revolutions. J To address the problems of low computational efficiency or limited applicability in J2 perturbation accumulation, this invention provides a multi-cycle double-pulse orbit transfer Δ based on the six roots of the vernal equinox under J2 perturbation. V Rapid estimation methods and systems.

[0006] Other features and advantages of the invention will become apparent from the following detailed description, or may be learned in part by practice of the invention.

[0007] According to a first aspect of the present invention, a method for estimating the minimum velocity increment of multi-cycle double-pulse orbital transfer under J2 perturbation is provided, the method comprising: S1. Initialization Processing: Obtain the initial orbital elements, transfer time, and single-pulse velocity increment upper limit constraint for the service spacecraft and the target spacecraft; convert the initial orbital elements into vernal equinox six-equinox numbers with mean longitude as the phase quantity, obtaining the starting point state vector and the ending point state vector respectively; the vernal equinox six-equinox numbers are represented as... ,in For the track semi-bore, Eccentricity Quantity, Eccentricity Quantity, Inclination angle Quantity, Inclination angle Quantity, Longitude; S2. Derive and calculate using the six orbital elements at the average vernal equinox. and ;in, The rate of change of mean longitude of the initial orbit of the service satellite, The sensitivity of the rate of change of horizontal longitude to the radius of the duct; S3. Determine the finite number of revolutions search window: Based on the single pulse velocity increment upper limit constraint, the starting state vector, and the sensitivity of the longitude drift rate to the half-diameter, calculate the maximum achievable additional number of transfer revolutions. N max and with and N The range of integers is used as a finite cycle search window; S4. Traverse the number of cycles and optimize the solution: For each candidate number of transition cycles in the search window... N : S4.1 Based on the phase closed-loop requirement, utilize the number of loops. N The phase closure equation is used to preliminarily estimate the half-aperture increment required by the first pulse. ; S4.2. Using the first pulse, calculate the four variables (excluding semi-circular diameter and longitude) from the six roots of the vernal equinox. , 、 、 、 Change 、 、 、 As an optimization variable; S4.3. Based on the required half-path increment caused by the first pulse and the current optimized variable values, calculate the tangential, normal, and radial velocity increment components of the first pulse according to the perturbation kinematic equations of the six roots of the vernal equinox. 、 、 And obtain the instantaneous state after the first pulse; S4.4 Using the analytical formula of complex plane rotation mapping, advance the complex quantity of the eccentric vector and the complex quantity of the orbital plane orientation in the instantaneous state after the first pulse, and calculate the end state of the drift segment before the second pulse after the transfer time. S4.5. Based on the closed equation constraint of "total change = change of the first pulse + change of the mid-segment drift + change of the second pulse", calculate the change in the number of the six roots at the vernal equinox required by the second pulse to reach the final state, and then calculate the tangential, normal, and radial velocity increment components of the second pulse. , , ; S4.6. The semi-aperture increment required by the first pulse estimated in step S4.1. The actual pulse components calculated in steps S4.3 and S4.5 are used to update the semi-aperture, resulting in the updated semi-aperture increment. Based on this, steps S4.3 to S4.5 are re-executed for iterative correction; S4.7 Construct an objective function based on the two pulse velocity increments. Within the feasible region of the optimization variables, use a nonlinear optimization algorithm to solve for the optimal value of the variable that minimizes the objective function, and record the corresponding minimum velocity increment. S5. Output the global optimal estimation result: compare all candidate lap numbers. N The minimum velocity increment is used as the final, rapid estimate of the total velocity increment for that segment of the track transfer.

[0008] In some exemplary embodiments, the calculation of the maximum achievable additional number of transfer cycles... N max Specifically:

[0009] In the formula, To serve the spacecraft's rate of change of longitude in its initial orbit, To serve the spacecraft's semi-circular diameter in its initial orbit; It is a mathematical substitution quantity; This represents the maximum tangential velocity increment that a single pulse engine can provide. To the total transfer time from the service spacecraft to the target spacecraft; To serve the overall longitude variation from the spacecraft to the target spacecraft; To serve the spacecraft's velocity in its initial orbit; N max It must be a non-zero integer. If the calculation... N max If the value is less than zero, it proves that the pulse is in the current Δ V t_max There is no solution under the current circumstances; the maximum tangential velocity component Δ needs to be increased. V t_max At the same time, it is necessary to pay attention to Δ V t_max It should not be too large, otherwise it may cause N max Excessive size will affect computational efficiency.

[0010] In some exemplary embodiments, the preliminary estimation of the required half-diameter increment caused by the first pulse... Specifically:

[0011] in, The average angular velocity of the serving star.

[0012] In some exemplary embodiments, the tangential direction is used to match the semi-circular diameter with the average motion; the normal direction is used to correct the orbital surface and calculate the change in horizontal longitude; the radial direction employs closed-loop balancing and is superimposed with adaptive regularization to eliminate in-plane residuals.

[0013] In some exemplary embodiments, the updated semi-pipeline increment Specifically:

[0014] in, , These represent the changes in longitude caused by the first pulse and the second pulse, respectively.

[0015] In some exemplary embodiments, the objective function is expressed as: .

[0016] According to a second aspect of the present invention, a system for estimating the minimum velocity increment of multi-cycle dual-pulse orbital transfer under J2 perturbation is provided, comprising: The initialization module is used to acquire and process input parameters, including the initial orbital elements of the service spacecraft and the target spacecraft, the transfer time, and the upper limit constraint of the single pulse velocity increment, and convert them into the vernal equinox six-element state based on the mean longitude. The preprocessing module is used to calculate the initial phase drift and determine a finite number of transfer cycles search window based on the pulse capability and phase sensitivity; The core estimation module is used to traverse the candidate transition cycles within the search window, and for each cycle, performs the following operations: The initial value of the semi-aperture increment of the first pulse is estimated based on the phase closed-loop requirement; The change in the number of six roots at the vernal equinox caused by the first pulse is used as the optimization variable; Calculate the velocity increment components of the two pulses, where the mid-stage drift is analytically propagated using complex plane rotation mapping; The second pulse increment is solved using closed-form equality constraints. The semi-path increment estimate is corrected by iterative updating; Call the nonlinear optimizer to find the optimal first pulse increment with the goal of minimizing the total velocity increment, and record the minimum estimated value under that number of revolutions; The output module compares the minimum estimate across all candidate lap counts and outputs the globally optimal total velocity increment estimate.

[0017] According to a third aspect of the present invention, a storage medium is provided having a computer program stored thereon, which, when executed by a processor, implements the method for estimating the minimum velocity increment of multi-cycle double-pulse orbital transfer under J2 perturbation as described in the first aspect.

[0018] According to a fourth aspect of the present invention, a computer program product is provided, on which a computer program is stored, wherein when the computer program is executed by a processor, the method for estimating the minimum velocity increment of multi-cycle double-pulse orbital transfer under J2 perturbation as described in the first aspect is implemented.

[0019] According to a fifth aspect of the present invention, an electronic device is provided, comprising: Processor; and Memory for storing the executable instructions of the processor; The processor is configured to implement the J2 perturbation-based multi-cycle dual-pulse orbit transfer minimum velocity increment estimation method described in the first aspect by executing the executable instructions.

[0020] The method and system for estimating the minimum velocity increment of multi-cycle double-pulse orbit transfer under J2 perturbation provided in the embodiments of the present invention have the following advantages compared with the prior art: (1) The proposed method uses the six roots of the vernal equinox as the basis for describing the orbital state and uses the mean longitude as the only phase quantity, thus unifying the description of the transition process from near-circular orbits to medium-high eccentricity orbits. This description method effectively avoids the singularity problem of classical Kepler roots in near-circular orbits or near low inclination angles. Since the method does not rely on the near-circular assumption throughout the modeling and solution process, it is effective in cases where the eccentricity is 0 ≤ e It remains continuously solvable and numerically robust across the entire range of values ​​less than 1.

[0021] (2) The proposed method uses the average advancement of complex plane rotation mapping to transform variables f, g, h, k The intermediate evolution is replaced by analytical-level mapping instead of long-time numerical integration. The proposed method directly generates a consistent intermediate state before the second pulse between two pulses, significantly reducing the computation time and error accumulation during drift propagation.

[0022] (3) The proposed method uses closed-form equations for constraint and unification. p, f, g, h, k, λ The constraint formula for the change is "total change = change of the first pulse + change of the mid-segment drift + change of the second pulse", achieving strict closure at the algorithm level. The proposed method completes the constraint loop on the variables by updating the half-path change caused by the first pulse twice in the phase constraint problem. The proposed method uses closed-loop balancing for the radial component and superimposes adaptive regularization, and limits the search range to a finite number of cycles with a phase sensitivity window, thus maintaining robustness in ill-conditioned geometry and multi-cycle cases.

[0023] (4) The proposed method decomposes each pulse into three components: tangential, normal, and radial, which respectively undertake the physical functions of scale matching, plane transformation and phase jump, and fine in-plane registration. The proposed method uses the change generated by the first pulse... The second pulse increment is solved using a nonlinear optimizer, and the closed-form equation is automatically closed, significantly reducing the optimization dimensionality and iterative complexity. The proposed method constructs the objective function from the total velocity increment of the orbital transfer synthesized from the two pulses, and selects the minimum Δ value within the candidate orbital number. V This results in a robust and widely applicable estimation algorithm, suitable for frequent use in upper-level multi-objective sequence optimization.

[0024] It should be understood that the above general description and the following detailed description are exemplary and explanatory only, and are not intended to limit the invention. Attached Figure Description

[0025] The accompanying drawings, which are incorporated in and constitute a part of this specification, illustrate embodiments consistent with the invention and, together with the description, serve to explain the principles of the invention. It is obvious that the drawings described below are merely some embodiments of the invention, and those skilled in the art can obtain other drawings based on these drawings without any inventive effort.

[0026] Figure 1 A schematic diagram illustrating the method flow of an exemplary embodiment of the present invention is provided. Detailed Implementation

[0027] Exemplary embodiments will now be described more fully with reference to the accompanying drawings. However, these exemplary embodiments can be implemented in many forms and should not be construed as limited to the examples set forth herein; rather, they are provided so that the invention will be more comprehensive and complete, and will fully convey the concept of the exemplary embodiments to those skilled in the art. The described features, structures, or characteristics may be combined in any suitable manner in one or more embodiments.

[0028] Furthermore, the accompanying drawings are merely illustrative of the invention and are not necessarily drawn to scale. The same reference numerals in the drawings denote the same or similar parts, and therefore repeated descriptions of them will be omitted. Some block diagrams shown in the drawings are functional entities and do not necessarily correspond to physically or logically independent entities. These functional entities can be implemented in software, in one or more hardware modules or integrated circuits, or in different network and / or processor devices and / or microcontroller devices.

[0029] This invention proposes a method based on the six roots of the vernal equinox. J 2. Perturbation-induced multi-cycle double-pulse orbital transfer Δ V A rapid estimation method and system are proposed for near-Earth orbit pulse transfer estimation scenarios, balancing applicability and online computational efficiency. First, a unified model of the origin and destination orbits is constructed using the vernal equinox six-root number, and mean longitude, rather than true longitude, is used as the unique phase quantity, thus ensuring that the phase is within the range of 0 ≤...e The continuity and differentiability of the phase description are maintained within the full eccentricity range of <1. Under the structure of "first pulse, mid-segment drift, second pulse", the mid-segment is analytically advanced using a complex plane rotation mapping: the evolution of variables is equivalent to the phase rotation of a complex plane, and the natural drift rate of longitude and its sensitivity to the semi-circle are calculated. Based on this, an initial guess of the semi-circle increment of the first pulse and a finite multi-cycle search window are given. Closed mapping is used to replace long-interval numerical integration, thereby significantly reducing the calculation time of the drift segment. Subsequently, each pulse is decomposed into three components: tangential, normal, and radial. The tangential component is used to match the semi-circle and the average motion; the normal component is used to correct the orbital plane and calculate the change in longitude; the radial component uses closed-loop balancing and is superimposed with adaptive regularization to eliminate in-plane residuals. Subsequently, a closed equation was constructed around the six roots at the vernal equinox: "Total change = Change of the first pulse + Change of the mid-segment drift + Change of the second pulse." By updating the semi-path increment of the initially guessed first pulse, the state and phase closure were ensured simultaneously without introducing additional penalty terms, thus enhancing numerical stability. Finally, the objective function was constructed from the total velocity increment of the orbital transfer from the two pulses, and a nonlinear optimizer was used to solve it for four variables in the initial six roots. The candidate values ​​were traversed within the revolution window, and the minimum Δ was selected. V This is used as the estimated value for this segment, thereby achieving high-precision Δ for three-component coordination. V Fast calculation. The above process ensures 0 ≤ e While being applicable to the entire domain, it also has the combined advantages of high computational efficiency, strict closure, and simple implementation.

[0030] The method of the present invention specifically includes the following steps: Step 1: Construct a near-Earth orbit pulse transfer estimation scenario For time-constrained and propulsion-constrained missions estimating single-pair single-pulse transfers in low Earth orbit, the following scenario is constructed: the initial orbital elements and transfer times of the service and target spacecraft are known, and the space environment considers Earth. J The orbital precession effect of the perturbation term is uniformly expressed using the six roots of the vernal equinox, and the true longitude is replaced with mean longitude for easier linearization. The engineering boundaries include the upper limit of the single pulse amplitude and the number of transfer cycles required for multi-cycle transfers. N With a given transition duration Δ t The transfer motion from the origin to the destination adopts a three-stage structure of "first pulse, mid-stage drift, and second pulse," where the pulse is applied only at the origin and destination ends, and the drift segment propels the spacecraft in the complex plane according to J2-mean dynamics. The scenario inputs are: the initial orbital elements of the service spacecraft and the target spacecraft, and the transfer time. Δt There is an upper limit constraint on the size of a single pulse. The scene output is: the total velocity increment Δ of the two pulses. VEstimation. To ensure that the costs of different origin-end combinations are both comparable and compatible, the scenario layer only defines the state representation, time base, and constraint boundaries; closed-loop equations and incremental analysis are implemented uniformly in the algorithm section, thereby providing standardized single-segment transition costs and feasible region criteria for the upper-level one-to-many sequence transition planning.

[0031] Step 2: Problem Description 2.1 Motion Model First, a motion model is constructed. To simplify the subsequent state propulsion evolution, considering only the long-term effect of the J2 perturbation, an improved vernal equinox orbital element is adopted: the true longitude in the traditional element is replaced with a mean longitude that is linearly related to time to describe the motion of the spacecraft and the target, specifically expressed as follows: . X The six components represent the orbital semi-circular diameter, eccentricity, etc. x Components, eccentricity y Quantity, tilt angle x Components, tilt angle y Components, mean longitude. The dynamic equations at this point are as follows: (1) In the formula, Re For the Earth's radius, μ is the Earth's gravitational constant. The mean longitude is chosen here instead of the true longitude of the traditional six roots of the vernal equinox because mean longitude, like the mean anomalous angle, exhibits linear variation characteristics, facilitating subsequent calculations.

[0032] To calculate the change in orbital elements caused by the pulse, it is also necessary to establish the instantaneous perturbation kinematic equations for the six elements at the vernal equinox. When establishing the perturbation kinematic equations, the perturbation acceleration is... u Decomposed into radial components u r tangential component u t and the normal component of the orbital plane u n Then the instantaneous perturbation kinematic equation of the six roots at the vernal equinox is: (2) In the formula, L For true longitude, n These are the average motion parameters. i For the track inclination angle, E The angles are for the nearest point; the others are mathematical substitutions for ease of calculation and will not be elaborated upon in the subsequent derivation. Their specific definitions are as follows: (3) For both sides of equation (2) in a very short time dtIntegrating within the inner quadrant, we can obtain the perturbation motion equation based on the pulse velocity increment: (4) In the formula, Δ V r Δ V t Δ V n These are the radial, tangential, and normal components of the velocity increment, which are the basic elements for constructing the velocity increment estimation expression.

[0033] Because the kinematic expression of the six roots of the vernal equinox does not depend on e With explicit participation, this method is applicable in 0 ≤ e <1 is absent in the entire range e It is unique and can uniformly handle transition estimation with different eccentricities.

[0034] 2.2 Constraint Model To address the problem of estimating the orbital transfer velocity increment, this method deconstructs the transfer process into three successive stages: "first pulse, mid-stage drift, and second pulse." Based on this, to accurately satisfy the terminal constraints, the total change in the six roots at the vernal equinox is modeled as a linear superposition of the dynamic contributions from the above three stages, i.e., constructing a closed-form equation constraint system of "total change = change in the first pulse + change in the mid-stage drift + change in the second pulse."

[0035] Specifically, p The change is composed of contributions from the tangential components of the two pulses, while the mid-range drift remains unchanged. p ; f, g, h, k The change is obtained by adding the three components of the first pulse at the starting phase, the complex plane rotation drift mapping based on J2 average in the middle section, and the three components of the second pulse after one phase closure, and is strictly closed to the target change in the form of an equation. λ The changes are composed of the natural drift in the middle section, the instantaneous jumps caused by the two normal pulses, and several terms of the entire cycle. Therefore, the problem constraints can ultimately be modeled as: (5) In the formula, subscripts 1 and 2 represent the first pulse and the second pulse, respectively; subscripts drift This represents the drift segment.

[0036] In practice, the six root variables at the vernal equinox are all implemented by "first calculating the intermediate value caused by the first pulse and drift, and then solving the necessary increment of the second pulse from the target change in a closed loop"; the radial component is directly eliminated in the plane residual by closed-loop balancing within each segment. The above equation constraints constitute the unique closed relationship of the entire estimation, which can ensure the consistent closure of the six roots without additional penalty terms or artificial upper bounds.

[0037] Step 3: Design Δ V Fast estimation algorithm Taking into account the requirements of multi-objective on-orbit missions for high-frequency evaluation and online decision-making, this step focuses on "determining the search window based on longitude and mid-course drift propulsion" and "constructing the constraint model algorithm and Δ..." V The two main lines of estimation and analytical solution construct Δ V Fast estimation algorithm: First, using longitude as the sole phase quantity and combining it with long-term perturbation averaging theory, a preliminary guess of the semi-aperture increment of the first pulse is directly given. Simultaneously, the sensitivity of the phase to orbital scales and the single-pulse capability jointly limit the search window to a finite number of revolutions. Between two pulses, an analytical mapping of complex plane rotation is used to advance the eccentricity vector and orbital surface parameters, replacing long-term numerical integration and ensuring computational efficiency across the entire eccentricity range. Second, constraints are constructed around the closed equation "total change = change from the first pulse + change from mid-range drift + change from the second pulse," explicitly accounting for the instantaneous phase jump caused by the normal pulse. A secondary update compensates for the semi-aperture change caused by the first pulse, thereby controlling... λ The parameter drift rate is ultimately constrained using mid-range drift. In pulse modeling, each pulse is decomposed into three components: tangential, normal, and radial. The radial component employs closed-loop balancing and is superimposed with an adaptive regularization term. Finally, within the candidate cycle number set, only the four increments of the first segment are optimized using a nonlinear optimizer. The second segment is automatically closed by a closed-loop relationship. Costs are weighted according to segment velocity, and the optimal solution is selected, achieving fast, robust, and reproducible ΔE in large time domains and multi-cycle scenarios. V Estimate.

[0038] 3.1 Search Window Determination Based on Longitude and Mid-course Drift Propulsion To maintain estimation efficiency and robustness in the large time domain and multi-cycle scenarios, the approach focuses on two aspects: "determining the search window" and "mid-course drift propulsion." First, using longitude as the sole phase quantity and combining long-term perturbation averaging theory, an initial estimate of the first pulse's semi-path increment is calculated. Based on the phase's sensitivity to orbital scales and single-pulse capability, a finite cycle search window is generated, and the geometric basis for tangential, normal, and radial decompositions is uniformly solved at the beginning of the segment. Subsequently, an analytical mapping of complex-plane rotation propels the eccentricity vector and orbital parameters between two pulses. This closed mapping replaces long-time numerical integration, significantly reducing computational overhead and being applicable across the entire eccentricity range. The output of this step is a constrained candidate set of cycles and a consistent intermediate state before the second pulse, providing a reliable starting point for subsequent constraint closure and three-component balancing.

[0039] 3.1.1 Define the search window Mean longitude represents the phase angle of a spacecraft or target in its orbit. Without maneuvering, the phase angle of the target and spacecraft at the end of their orbits is...t f phase difference Δ λ It can be represented as: (6) In the formula, , The subscript 0 represents the state variable of the initial trajectory, that is, the state before the first pulse; the subscript... f This represents the state variable of the final trajectory after the mission is completed, i.e., the state after the second pulse; The natural longitude drift rate affected by J2 perturbation is given in equation (1); , The argument of latitude, The perigee argument, The right ascension of the ascending node.

[0040] From the sixth equation of equation (2), we can see that the variables in the six roots of the vernal equinox are... p, f, g, h, k All of these factors affect the rate of change of mean longitude. Therefore, to achieve zero phase difference between the spacecraft and the target orbit at the end of the orbit, the orbital elements need to be precisely adjusted using the first pulse. Essentially, this process involves controlling the rate of change of mean longitude by altering the elements. This, in turn, affects the mid-range drift trend of the phase difference, ultimately achieving closed-loop constraint. The change in the spacecraft's phase change rate is calculated as follows: (7) In the formula, N This represents the additional number of orbits the servicing spacecraft performs relative to a scenario without a pulse. Because... The principal quantity is the average motion , Half-bore p The order of magnitude of the partial derivatives is usually significantly larger than the partial derivatives with respect to other variables, therefore it can be simplified to considering only the adjustment of the semi-path. p The effect on the rate of phase change, among which right p The partial derivatives are: (8) In the formula, n These are the average motion parameters. Therefore, the change in half-pipe diameter caused by the first pulse can be obtained: (9) in, The average motion parameter is the average angular velocity of the serving satellite.

[0041] At the same time, according to the first equation of formula (4), the tangential component of the first pulse can also be obtained. The relationship between it and the resulting change in the semi-major axis: (10) in, As a mathematical substitution quantity, it is calculated using the six-root state quantity of the servicing spacecraft at the vernal equinox in its initial orbit. The specific formula is as follows: ; In the formula, V 0 represents the spacecraft's velocity before the first pulse. If it is stipulated that the magnitude of the tangential component of the first pulse does not exceed Δ... V t_max By solving equations (9) and (10) simultaneously, the maximum number of laps can be derived. N max expression: (11) In the formula, N max It must be a non-zero integer. If the calculation... N max If the value is less than zero, it proves that the pulse is in the current Δ V t_max There is no solution under the current circumstances; the maximum tangential velocity component Δ needs to be increased. V t_max At the same time, it is necessary to pay attention to Δ V t_max It should not be too large, otherwise it may cause N max Excessive size negatively impacts computational efficiency. Ultimately N The range of values ​​for is: (12) After determining the transfer semi-path and search window of the first pulse, further analytical derivation of the orbit transfer estimation can be performed.

[0042] 3.1.2 Mid-stage drift propulsion In the transfer scenario of "first pulse, mid-course drift, second pulse," the mid-course drift phase typically lasts for multiple orbital cycles. If a stepwise numerical integration method is used to advance the average perturbation equations of the six roots at the vernal equinox, full-time integration must be performed for every candidate transfer cycle. Due to numerical stability requirements, the integration step size must be significantly smaller than the orbital cycle, resulting in an extremely large total number of integration steps and a computational load that increases proportionally with the number of transfer cycles. More critically, long-time-domain numerical integration introduces truncation and rounding errors. These errors accumulate over multiple orbital cycles, causing deviations in the final orbital state and severely affecting the accuracy of phase and orbital plane parameter closure, thereby reducing the reliability of velocity increment estimation. Therefore, it is necessary to use a computationally more efficient analytical mapping to replace long-time-domain stepwise numerical integration.

[0043] In the J2 average sense, semi-pipe diameter pThe long term is zero ( ),and f, g, h, k It mainly exhibits uniform precession characteristics, and its precession rate is only related to the initial state and the Earth's constant. Based on this, the initial state can be regarded as a constant, thus during the drift period Δ t The internally coupled nonlinear dynamic equations are transformed into a linear system with constant coefficients in the complex plane. This simplification makes... f, g, h, k The evolution can be analytically expressed through a single complex plane rotation, forming a closed-form solution. This method preserves... J While mitigating the main long-term effects of perturbation, it significantly reduces computational complexity.

[0044] Here, we introduce a complex plane for modeling, let: (13) in, Let be an "eccentric vector" on the complex plane, with real part f virtual part g These correspond to the two orthogonal components of the eccentricity vector in the orbital plane in the vernal equinox reference frame, and their magnitudes represent the eccentricity. e ; Let be the "orbital plane orientation vector" on the complex plane, with real part . h virtual part k The equivalent parameters corresponding to the inclination angle and the right ascension of the ascending node have moduli related to the orbital inclination angle. Embedding the two pairs of real variables into the complex plane, according to the J2 average perturbation theory, the evolution of the complex variables under the above definition can be directly given by the following rotational mapping: (14) Its analytical solution is: (15) In the formula, , These represent the rotational angular rates: (16) in, e Let be the eccentricity, and its expression is: .

[0045] The drift segment begins after the first pulse ends and ends before the second pulse begins; therefore, let the starting time variable be: (17) The rotation angle at this time is: (18) Finally, after rotation in the complex plane, the variable at the end of the drift segment is: (19) This closed-form progression is a purely algebraic operation; for anyN The computational cost is extremely low, significantly reducing the overall evaluation time.

[0046] 3.2 Algorithm Construction of Constraint Model and Δ V Estimate analytical solution 3.2.1 Algorithm Construction of Constraint Model To accurately satisfy the terminal constraints, the total change of the six roots at the vernal equinox is modeled as a dynamic contribution consisting of three stages: the first pulse, the mid-stage drift, and the second pulse. The instantaneous effect of the pulse is first-order linear, and the transmission of the "small deviation" caused by the mid-stage drift to the first pulse is also first-order linear. The series connection of linear components remains linear. Therefore, the contributions of the three stages can be linearly added and directly solved, thus constructing a closed equation constraint system of "total change = change of the first pulse + change of the mid-stage drift + change of the second pulse".

[0047] Before the first pulse, the six roots at the starting vernal equinox are: (20) The target vernal equinox point has six roots: (twenty one) Since the pulse is an instantaneous action, a state variable will change instantaneously after the first pulse, i.e., at the start of the drift segment. Let the instantaneous variable after the pulse be: (twenty two) And there are: (twenty three) In the formula, The change generated by the first pulse can be calculated using equation (4).

[0048] Similarly, after the drift segment and before the second pulse, the variables are: (twenty four) And there are: (25) In the formula, This represents the change caused by mid-range drift.

[0049] Finally, the second pulse will correct the six roots to the target end state, completing the closed loop for the entire constraint. At this point, we have: (26) In the formula, The change generated by the first pulse can be calculated using equation (4).

[0050] Therefore, there are constraints on the change in the six roots at the vernal equinox for the entire pulse transfer process: (27) Within the framework of the six roots of the vernal equinox J 2. Average half-pipe diameter p No drift occurs, but f, g, h, k The mid-range evolution can be quickly given by a closed mapping of complex plane rotation, thus these five quantities can achieve equation closure within the structure of "first pulse, mid-range drift, second pulse". In contrast, the mean longitude λ is not an independent state; its expression is: (28) Pulse-induced phase jump and It is not an independent variable; its magnitude and sign are determined by the normal component. The only determination is the geometric coupling: (29) in ,and Again received h, k Used to implement Constrained by the balancing constraints, therefore the Δ caused by it λ The jump is predetermined, and there is no means to decouple it and adjust it independently; similarly, the tangential and radial components are... λ The immediate impact is only through f, g The geometric terms indirectly reflect this and have been used for balancing. It can no longer stand alone as λ Leave room for degrees of freedom. Therefore, it is evident that if we maintain the variables... p, f, g, h, k Forced modification under the premise of closure and This will directly damage the h, k or f, g The equations cannot be balanced by adjusting them at this point. λ The pulse change thus satisfies the closed constraint.

[0051] Therefore, the only continuous, adjustable channel that satisfies phase closure and does not violate the five root closed constraints is the natural drift term. Given a transition time Δ t First, we calculate... and By adjusting the tangential action of the first pulse Change The average value is used to make up for the remaining phase difference, and the updated semi-pipeline is updated. The expression is: (30) The closed constraint of the longitude is completed by the constraint of equation (30).

[0052] 3.2.2 Δ V Analytical Construction and Optimization of Objective Construction Given the number of laps N The change in tangential half-pipe diameter obtained from the phase closed equation Under the premise of only the four vernal equinoxes of the first pulse, the increment Δ f 1 、 Δ g 1 、 Δ h 1 、 Δ k 1. To optimize the variables, construct the three-component velocities of the two pulse segments and establish the L2 norm objective function.

[0053] Based on the corresponding transformation of equation (4), the expressions for each component can be obtained. The tangential component is given by the phase closed equation. By unique definition and mathematical processing of the first line of equation (4), we can obtain: (31) In the formula, subscripts I Indicates the first I Sub-pulse I =1,2.

[0054] The amplitude of the normal component is determined by By squaring both sides of the equations in the fourth and fifth lines of formula (4), adding them together, and then taking the square root, we can obtain: (32) To determine the magnitude of the radial velocity increment, it is necessary to... and right The closed-form balancing is performed, and the residual between the residual and the target value is calculated. This residual is then solved using the closed-form balancing method. To improve numerical robustness, adaptive regularization is introduced during the solution process.

[0055] (33) Therefore, the in-plane residual that the radial component needs to compensate for can be expressed as: (34) Simultaneously radial term Also correct A contribution is expressed as: (35) but , Must meet This can be solved by constructing a least squares problem. The variables are then represented in matrix form: (36) Taking Tikhonov regularization, the regularized least squares solution is: (37) From this point on All three components have been calculated, combined with the known number of transfer cycles. N And the four optimization variables Δ that were hypothesized. f 1 、 Δ g 1 、 Δ h 1 、 Δ k 1. Construct a four-variable objective function: (38) Since all four optimization variables represent the change in orbital elements caused by the pulse, and their numerical magnitudes are relatively small, a compact feasible region can be defined based on their physical meaning.

[0056] (39) For this smooth, box-only constrained low-dimensional objective, a nonlinear optimizer algorithm is used to find a local minimum on the current optimization variable, and the minimum cost is recorded as the estimate of the combination. Finally, the minimum among all combinations is taken as the estimate of the transfer velocity increment.

[0057] Example This invention proposes a method for modeling total eccentricity based on the six roots of the vernal equinox. J 2-Perturbation Dual-Pulse Multi-Cycle Near-Earth Orbit Transfer Δ V A rapid estimation method for the kinematics of the six roots at the vernal equinox. J Under the average dynamics model, a three-stage balancing process ("first pulse, mid-stage drift, second pulse"), explicit phase closure, and multi-cycle finite search are achieved, balancing robustness and online efficiency. For example... Figure 1 As shown, the specific process is as follows: Step 1: Initialization process, calculate the six root numbers of the starting and ending vernal equinoxes, and calculate their difference to obtain the transition time Δ. t The change in the six roots within; Step 2: Derive and calculate the orbital elements using the six roots of the average vernal equinox. and ;in, The rate of change of mean longitude in the initial orbit of the service satellite (before the first pulse). The effect of the rate of change of longitude on the radius of the duct (rate); Step 3, based on single tangential capability and Calculate the maximum number of transfer cycles Nmax And obtain the traversed intervals When performing a traversal loop, N Values ,like N max If the value is less than 0, the segment is marked as infeasible and the process exits; where, the single tangential capability is the maximum tangential velocity increment that the engine can provide. ; Step 4: Use equation (9) to preliminarily predict and estimate the semi-aperture increment generated by the first pulse. Let the change in the six roots of the vernal equinox be determined by the first pulse. To optimize the variables, for each candidate N Take the value given by the current optimization solver x Calculate the three-directional components of the first pulse and the change in longitude generated by the first pulse. The value of is used to obtain the size of the first pulse; Step 5: Use complex plane rotation mapping to quickly calculate the mid-segment drift change of the six roots; Step 6: Solve for the required balancing increment for the second pulse based on the closed-form variation constraint, and similarly use... p , f, g, h, k The calculation of the changes in velocity increments in the tangential, normal, and radial directions of the second pulse, and the change in longitude generated by the second pulse. The value of is used to obtain the size of the second pulse; Step 7, obtain and Substituting back into step 4, equation (9) is updated to equation (30), and the updated semi-pipeline increment is recalculated. Use the updated version Repeat steps 5 and 6 to obtain the solutions for the first and second pulses; Step 8: Minimize the objective function using nonlinear optimization. J Record the number of laps. N Minimum orbital transfer Δ V The estimated value is used; the minimum value within the traversed interval is taken as the global output.

[0058] Other embodiments of the invention will readily occur to those skilled in the art upon consideration of the specification and practice of the invention herein. This application is intended to cover any variations, uses, or adaptations of the invention that follow the general principles of the invention and include common knowledge or customary techniques in the art not disclosed herein. The specification and embodiments are to be considered exemplary only, and the true scope and spirit of the invention are indicated by the claims.

[0059] It should be understood that the present invention is not limited to the precise structure described above and shown in the accompanying drawings, and various modifications and changes can be made without departing from its scope. The scope of the invention is defined only by the appended claims.

Claims

1. A method for estimating the minimum velocity increment of multi-cycle double-pulse orbit transfer under J2 perturbation, characterized in that, The method includes: S1. Initialization Processing: Obtain the initial orbital elements, transfer time, and single-pulse velocity increment upper limit constraint for the service spacecraft and the target spacecraft; convert the initial orbital elements into vernal equinox six-equinox numbers with mean longitude as the phase quantity, obtaining the starting point state vector and the ending point state vector respectively; the vernal equinox six-equinox numbers are represented as... ,in For the track semi-bore, Eccentricity Quantity, Eccentricity Quantity, Inclination angle Quantity, Inclination angle Quantity, Longitude; S2. Derive and calculate using the six orbital elements at the average vernal equinox. and ;in, The rate of change of mean longitude of the initial orbit of the service satellite, The sensitivity of the rate of change of horizontal longitude to the radius of the duct; S3. Determine the finite number of revolutions search window: Based on the single pulse velocity increment upper limit constraint, the starting state vector, and the sensitivity of the longitude drift rate to the half-diameter, calculate the maximum achievable additional number of transfer revolutions. N max and with and N The range of integers is used as a finite cycle search window; S4. Traverse the number of cycles and optimize the solution: For each candidate number of transition cycles in the search window... N : S4.1 Based on the phase closed-loop requirement, utilize the number of loops. N The phase closure equation is used to preliminarily estimate the half-aperture increment required by the first pulse. ; S4.

2. Using the first pulse, calculate the four variables (excluding semi-circular diameter and longitude) from the six roots of the vernal equinox. , 、 、 、 Change 、 、 、 As an optimization variable; S4.

3. Based on the required half-path increment caused by the first pulse and the current optimized variable values, calculate the tangential, normal, and radial velocity increment components of the first pulse according to the perturbation kinematic equations of the six roots of the vernal equinox. 、 、 And obtain the instantaneous state after the first pulse; S4.4 Using the analytical formula of complex plane rotation mapping, advance the complex quantity of the eccentric vector and the complex quantity of the orbital plane orientation in the instantaneous state after the first pulse, and calculate the end state of the drift segment before the second pulse after the transfer time. S4.

5. Based on the closed equation constraint of "total change = change of the first pulse + change of the mid-segment drift + change of the second pulse", calculate the change in the number of the six roots at the vernal equinox required by the second pulse to reach the final state, and then calculate the tangential, normal, and radial velocity increment components of the second pulse. , , ; S4.

6. The semi-aperture increment required by the first pulse estimated in step S4.

1. The actual pulse components calculated in steps S4.3 and S4.5 are used to update the semi-aperture, resulting in the updated semi-aperture increment. Based on this, steps S4.3 to S4.5 are re-executed for iterative correction; S4.7 Construct an objective function based on the two pulse velocity increments. Within the feasible region of the optimization variables, use a nonlinear optimization algorithm to solve for the optimal value of the variable that minimizes the objective function, and record the corresponding minimum velocity increment. S5. Output the global optimal estimation result: compare all candidate lap numbers. N The minimum velocity increment is used as the final, rapid estimate of the total velocity increment for that segment of the track transfer.

2. The method according to claim 1, characterized in that, The calculation of the maximum achievable additional transfer cycles... N max Specifically: In the formula, To serve the spacecraft's rate of change of longitude in its initial orbit, To serve the spacecraft's semi-circular diameter in its initial orbit; It is a mathematical substitution quantity; This represents the maximum tangential velocity increment that a single pulse engine can provide. To the total transfer time from the service spacecraft to the target spacecraft; To serve the overall longitude variation from the spacecraft to the target spacecraft; To serve the spacecraft's velocity in its initial orbit; N max It must be a non-zero integer. If the calculation... N max If the value is less than zero, it proves that the pulse is in the current Δ V t_max There is no solution under the current circumstances; the maximum tangential velocity component Δ needs to be increased. V t_max At the same time, it is necessary to pay attention to Δ V t_max It should not be too large, otherwise it may cause N max Excessive size will affect computational efficiency.

3. The method according to claim 1, characterized in that, The preliminary estimate of the required half-diameter increment due to the first pulse Specifically: in, The average angular velocity of the serving star.

4. The method according to claim 1, characterized in that, The tangential direction is used to match the semi-circular diameter with the average motion; the normal direction is used to correct the orbital surface and calculate the change in longitude; the radial direction uses closed-loop balancing and is superimposed with adaptive regularization to eliminate in-plane residuals.

5. The method according to claim 1, characterized in that, The updated semi-aperture increment Specifically: in, , These represent the changes in longitude caused by the first pulse and the second pulse, respectively.

6. The method according to claim 1, characterized in that, The objective function is expressed as: .

7. A system for estimating the minimum velocity increment of multi-cycle double-pulse orbit transfer under J2 perturbation, characterized in that, include: The initialization module is used to acquire and process input parameters, including the initial orbital elements of the service spacecraft and the target spacecraft, the transfer time, and the upper limit constraint of the single pulse velocity increment, and convert them into the vernal equinox six-element state based on the mean longitude. The preprocessing module is used to calculate the initial phase drift and determine a finite number of transfer cycles search window based on the pulse capability and phase sensitivity; The core estimation module is used to traverse the candidate transition cycles within the search window, and for each cycle, performs the following operations: The initial value of the semi-aperture increment of the first pulse is estimated based on the phase closed-loop requirement; The change in the number of six roots at the vernal equinox caused by the first pulse is used as the optimization variable; Calculate the velocity increment components of the two pulses, where the mid-stage drift is analytically propagated using complex plane rotation mapping; The second pulse increment is solved using closed-form equality constraints. The semi-path increment estimate is corrected by iterative updating; Call the nonlinear optimizer to find the optimal first pulse increment with the goal of minimizing the total velocity increment, and record the minimum estimated value under that number of revolutions; The output module compares the minimum estimate across all candidate lap counts and outputs the globally optimal total velocity increment estimate.

8. A storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, it implements the method for estimating the minimum velocity increment of multi-cycle double-pulse orbital transfer under J2 perturbation as described in any one of claims 1 to 6.

9. A computer program product, comprising a computer program, characterized in that, When the computer program is executed by the processor, it implements the method for estimating the minimum velocity increment of multi-cycle double-pulse orbital transfer under J2 perturbation as described in any one of claims 1 to 6.

10. An electronic device, comprising: processor; as well as Memory for storing the executable instructions of the processor; The processor is configured to implement the method for estimating the minimum velocity increment of multi-cycle dual-pulse orbit transfer under J2 perturbation as described in any one of claims 1 to 6 by executing the executable instructions.

Images