Double-satellite formation satellite orbit determination method and device based on inter-satellite relative optical imaging and microwave measurement

Abstract

The invention provides a double-satellite formation satellite orbit determination method and device based on inter-satellite relative optical imaging and microwave measurement, and the method comprises the steps: building a satellite motion equation based on orbit dynamics, and obtaining a state equation; establishing an observation equation based on inter-satellite relative optical imaging and microwave measurement; measuring the accurate orbit of the reference satellite A; the reference star A performs optical imaging on the target star B to obtain azimuth information of the star B relative to the star A; carrying out microwave distance measurement on the reference satellite A and the target satellite B to obtain relative distance information; calculating orbit parameters of the satellite B, and obtaining position coordinates and a satellite orbit of the satellite B according to the observation coordinate value of the satellite B relative to the satellite A and the state equation; and obtaining an orbit determination result of the target satellite at any moment. The economic benefit of satellite in-orbit operation can be improved.

Description

Technical Field

[0001] This invention belongs to the field of formation spacecraft orbit determination, specifically relating to a method and apparatus for determining the orbit of a two-satellite formation satellite based on inter-satellite relative optical imaging and microwave measurement. Background Technology

[0002] Conventional methods for determining the orbits of spacecraft flying in formation involve individually determining the orbits of each satellite through radio ranging, velocity measurement, and VLBI angle measurement. This method requires long-term coordination and tracking by ground control stations, and the required control resources increase as the number of satellites increases. When ground control resources are insufficient, it becomes difficult to determine the orbits of satellites simultaneously. For low Earth orbit satellites, GNSS navigation constellations can also provide some orbit determination and positioning. However, when satellites are in deep space orbits such as lunar orbits, GNSS signals are extremely weak, making it difficult for satellites to effectively receive and use GNSS signals for orbit determination and positioning.

[0003] Chinese patent application CN102305630A discloses an autonomous orbit determination method for SAR satellites based on extended Kalman filtering. This method establishes observation equations using the distance between the SAR and ground markers and the Doppler frequency shift between the SAR and ground markers as observation quantities, and then establishes a recursive equation for extended Kalman filtering to obtain the satellite's position and velocity vectors. This method is designed for SAR satellites, requires the establishment of ground markers, and is not applicable to satellites flying in formation.

[0004] Chinese patent application CN118483779A discloses a rapid satellite positioning method based on the unscented Kalman algorithm. This method includes a positioning system that obtains the direction of the target system based on the azimuth and elevation angles of the satellite, captures the desired low-Earth orbit satellite in a satellite coordinate system relative to the satellite, an algorithm for calculating the intersection of a straight line with the Earth's sphere and providing the latitude and longitude coordinates of the target satellite's radiation source, and an unscented Kalman-based filtering algorithm for rapid modeling and estimation. However, it only proposes the filtering algorithm and does not address scenarios applicable to deep-space formation satellite flight, nor does it cover methods for obtaining the relative position between two satellites.

[0005] Existing research on autonomous orbit determination for navigation constellations based on inter-satellite ranging proposes an autonomous orbit determination algorithm that utilizes an improved Kalman filter to fuse dynamic information and inter-satellite ranging information, and corrects the initial integral value in real time. It employs IGS precise ephemeris simulation of inter-satellite measurements. However, this method does not use optical imaging methods to determine inter-satellite azimuth information.

[0006] In summary, although existing orbit positioning technologies for deep-space dual-satellite formation flights still have certain limitations, they can only determine the orbit of a single satellite by relying on high-power deep-space ground stations, which is highly dependent on ground resources. Summary of the Invention

[0007] To address the aforementioned technical problems, this invention provides a method and apparatus for determining the orbits of two-satellite formations based on inter-satellite relative optical imaging and microwave measurement. This method reduces reliance on ground-based telemetry and control resources without relying on external navigation signals (such as GNSS). Given a reference satellite orbit, the relative position and velocity between the two satellites are autonomously measured through optical imaging and microwave ranging between the reference and the target satellite. Combined with precise orbit determination of the target satellite obtained through ground-based processing, the orbit is transferred from the reference satellite to the target satellite. This invention effectively improves the efficiency of constellation satellite formation flight, reduces reliance on external navigation systems, enhances system independence and anti-interference capabilities, reduces system costs and energy consumption, and improves the economic benefits of satellite on-orbit operation.

[0008] To achieve the above objectives, the present invention adopts the following technical solution:

[0009] A method for determining the orbits of binary satellite formations based on inter-satellite relative optical imaging and microwave measurements includes the following steps:

[0010] Step 1: Establish the satellite motion equations based on orbital dynamics, and then obtain the state equations;

[0011] Step 2: Establish observation equations based on inter-satellite relative optical imaging and microwave measurements;

[0012] Step 3: Determine the precise orbit of reference star A;

[0013] Step 4: Use reference star A to perform optical imaging of target star B to obtain the azimuth information of star B relative to star A;

[0014] Step 5: Perform microwave ranging between reference satellite A and target satellite B to obtain relative distance information;

[0015] Step 6: Calculate the orbital parameters of satellite B. Based on the observed coordinates of satellite B relative to satellite A, obtain the position coordinates and satellite orbit of satellite B according to the state equation.

[0016] Step 7: Process the observation information using the extended Kalman filter algorithm, and perform time and measurement updates. Repeat steps 3 to 6, and perform step 7 after each loop to perform time and measurement updates, so as to obtain the orbit determination result of the target satellite at any time.

[0017] The present invention also provides a device for determining the orbit of a binary satellite formation based on inter-satellite relative optical imaging and microwave measurement, comprising the following modules:

[0018] The state equation establishment module establishes satellite motion equations based on orbital dynamics, thereby obtaining the state equations;

[0019] The observation equation establishment module establishes observation equations based on inter-satellite relative optical imaging and microwave measurements.

[0020] The measurement module determines the precise orbit of the reference satellite A.

[0021] The imaging module uses reference star A to perform optical imaging of target star B, obtaining the orientation information of star B relative to star A.

[0022] The ranging module performs microwave ranging between the reference satellite A and the target satellite B to obtain relative distance information.

[0023] The calculation module calculates the orbital parameters of satellite B. Based on the observed coordinates of satellite B relative to satellite A, it obtains the position coordinates and satellite orbit of satellite B according to the state equation.

[0024] The orbit determination module uses the extended Kalman filter algorithm to process observation information and performs time and measurement updates to obtain the orbit determination result of the target satellite at any time.

[0025] This invention also provides an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, it implements the steps of the above-described method for determining the orbit of a binary satellite formation based on inter-satellite relative optical imaging and microwave measurement.

[0026] This invention also provides a non-transitory computer-readable storage medium storing a computer program thereon, which, when executed by a processor, implements the steps of the above-described method for determining the orbit of a dual-satellite formation based on inter-satellite relative optical imaging and microwave measurement.

[0027] Beneficial effects:

[0028] This invention is primarily used to determine the relative positions of two satellites in a constellation formation when telemetry and control resources are insufficient. For example, after precise orbit determination of one satellite, optical imaging of satellite B is performed using satellite A, and the relative positions of satellites A and B are calculated from the images. Simultaneously, combined with the attitude information of satellite A itself, the relative orbital position of satellite B is calculated through continuous iteration using the Kalman rate wave algorithm. Its advantages are mainly reflected in the following aspects:

[0029] 1. High Autonomy: This invention reduces reliance on ground-based orbit determination systems or other external systems, instead utilizing inter-satellite relative optical imaging and microwave measurements for autonomous positioning between the reference satellite and the target satellite. This enhances the autonomy and independence of the binary satellite formation and reduces dependence on external systems to some extent.

[0030] 2. High-precision relative positioning: This invention combines optical imaging and microwave measurement technologies. Optical imaging provides high-precision angle measurement information, while microwave measurement provides high-precision distance measurement information. The combination of these two technologies significantly improves the relative positioning accuracy of binary star formations.

[0031] 3. High adaptability: This invention is applicable to binary satellite formations with different orbital altitudes, speeds, and attitudes. Whether in low-to-medium orbit, high orbit, or deep space orbit, this method is applicable.

[0032] 4. Enhance satellite formation coordination capabilities: The accurate positioning information of this invention helps to enhance the coordination capabilities of satellite formations, enabling satellites in the formation to work more closely together and improve overall performance.

[0033] 5. Reduced operating costs: Because this invention does not rely on ground equipment or other external systems, it can reduce the operating costs of satellite formations. At the same time, the high-precision positioning information also helps reduce the waste of satellite resources caused by positioning errors.

[0034] 6. Promoting the development of space orbit determination: This invention combines optical imaging and microwave measurement technologies, providing new ideas and methods for the development of space orbit determination. This contributes to the continuous innovation and development of space orbit determination technology. Attached Figure Description

[0035] Figure 1 This is a flowchart of a method for determining the orbit of a binary satellite formation based on inter-satellite relative optical imaging and microwave measurement, according to an embodiment of the present invention. Detailed Implementation

[0036] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention. Furthermore, the technical features involved in the various embodiments of this invention described below can be combined with each other as long as they do not conflict with each other.

[0037] This invention combines optical imaging technology and microwave measurement technology to achieve accurate transmission of orbital data for satellite formations. Optical imaging technology captures images of the relative positions between satellites using onboard cameras, while microwave measurement technology measures the distance and velocity between two satellites by transmitting and receiving microwave signals. The combination of these two technologies provides precise information on the relative position, distance, and velocity between two satellites, thus enabling accurate positioning of a two-satellite formation. For ease of description, the reference satellite in the two-satellite formation is referred to as satellite A, and the target satellite being measured is referred to as satellite B.

[0038] In terms of hardware, satellite A in the constellation is equipped with microwave ranging equipment and optical imaging equipment, while satellite B is equipped with Ka-band communication payload and optical beacon. Satellites A and B fly in formation in lunar orbit. Satellite A conducts continuous microwave ranging with satellite B through its microwave ranging equipment and continuously photographs and images satellite B's optical beacon using its imaging equipment.

[0039] Inter-satellite microwave ranging results were obtained using the microwave ranging equipment on satellite A. Simultaneously, images of satellite B were captured using an imaging device, and the position and orientation of satellite B relative to satellite A were determined by analyzing the images. The microwave ranging data and the position and orientation data obtained from image analysis were used as inputs, and a Kalman filter algorithm was employed for data fusion processing.

[0040] like Figure 1 As shown, the method for determining the orbit of a binary satellite formation based on inter-satellite relative optical imaging and microwave measurement according to an embodiment of the present invention includes the following steps:

[0041] Step 1: Establish the satellite motion equations based on orbital dynamics, and then obtain the state equations:

[0042] Let the state vector be... ,in These are the satellite's position coordinates in a certain coordinate system. The corresponding velocity component is t, which represents time, and the superscript T indicates the transpose of the matrix.

[0043] Assuming the satellite's motion conforms to a dynamic model, the state equation can be expressed as: ,in It is a state transition function. It is process noise, usually assumed to be Gaussian white noise, and its covariance matrix is... Δt is the time interval between two consecutive state updates. For simple two-body problems, It can be derived from Newton's laws of motion and the law of universal gravitation.

[0044] Step 2: Establish observation equations based on inter-satellite relative optical imaging and microwave measurements:

[0045] The equation for observing satellite A is:

[0046] Observation vector ;

[0047] The observation equation is ;

[0048] in It is the observation matrix. It is observation noise, also assumed to be Gaussian white noise, with a covariance matrix of... ;

[0049] Distance equation Represented as:

[0050] ;

[0051] in, These are the location coordinates of the ground station. This is ranging noise. The superscript T indicates the transpose of the matrix.

[0052] The equation for observing satellite B is:

[0053] Observation vector ;

[0054] The observation equation is ;

[0055] in It is the observation matrix. The noise from the B-satellite observation is also assumed to be Gaussian white noise, and its covariance matrix is... ;

[0056] The distance equation can also be expressed as:

[0057] ;

[0058] in These are the position coordinates of star A. It is ranging noise.

[0059] Step 3: Determine the precise orbit of reference star A:

[0060] The ground station receives observation data from satellite A and calculates the distance dA(t) between satellite A and the station, as well as other orbital information, based on the signal propagation time, where t represents time. Satellite A's orbit has six elements, including its semi-major axis. eccentricity ,inclination Right ascension of ascending node Angular distance from perigee And true near point angle .

[0061] Step 4: Reference satellite A performs optical imaging of target satellite B to obtain the azimuth information of B relative to A. Simultaneously, the imaging equipment on satellite A takes pictures of the optical beacon of B at equal time intervals. The images are processed, and by identifying the pixel position and attitude of the optical beacon in the image, combined with A's own attitude information and the intrinsic and extrinsic parameters of the imaging equipment, the direction vector of B relative to A is calculated. The attitude of satellite A can be determined using onboard gyroscopes, solar sensors, Earth sensors, and other measurement equipment. Specifically, this includes the following sub-steps:

[0062] 1) Image Acquisition and Preprocessing: The imaging equipment of satellite A takes pictures of the optical beacon of satellite B at the same time intervals to obtain image data. The captured images are preprocessed, including noise reduction and enhancement, to improve the accuracy of subsequent processing.

[0063] 2) Pixel position recognition: (e.g., Pix2Pose algorithm, which translates the target through two-dimensional target detection to align the bounding box with the target center and eliminate background interference) to identify the pixel position of the optical beacon in the image.

[0064] 3) Attitude Estimation: The attitude of the optical beacon on satellite B is estimated by combining the attitude information of satellite A itself. This involves converting the pixel coordinates in the image into three-dimensional spatial coordinates of satellite B relative to satellite A. Using the homography transformation principle, the mapping relationship between the target point in the local coordinate system and the corresponding target point in the pixel coordinate system plane is established, resulting in the homography matrix H.

[0065] 4) Imaging equipment parameter calibration: This requires the imaging equipment's intrinsic parameters (such as focal length, principal point coordinates, etc.) and extrinsic parameters (such as the imaging equipment's position and attitude relative to the A-satellite coordinate system). These parameters can be obtained through the calibration process.

[0066] 5) Direction Vector Calculation: Using the intrinsic and extrinsic parameters of the imaging device, the pixel coordinates are converted to coordinates in the camera coordinate system. Converting the coordinates in the camera coordinate system to those in satellite A's coordinate system requires using satellite A's attitude information. Finally, the direction vector of satellite B relative to satellite A is calculated. .

[0067] 6) Adjust the attitude of satellite A based on the direction vector results.

[0068] Step 5: Microwave ranging is performed between reference satellite A and target satellite B to obtain relative distance information. The microwave ranging equipment on satellite A transmits microwave signals to satellite B at regular time intervals and receives the signals transmitted by satellite B. Using a two-way comparative ranging method, the distance is determined based on the propagation time of the microwave signals. Calculate the distance between satellite A and satellite B. (where c is the speed of light). Specifically, it includes the following steps:

[0069] 1) The dual-satellite ranging module is powered on and initialized, and ranging frames are sent unidirectionally to complete time synchronization.

[0070] 2) The two satellite ranging modules use radio frequency signals to perform clock correction and tracking processes.

[0071] 3) During the distance maintenance process, the two-way comparison method is also used to measure the time difference, and the synchronization error gradually decreases during the convergence process.

[0072] Step Six: Calculate the orbital parameters of satellite B. Based on the observed coordinates of satellite B relative to satellite A, the position coordinates and orbit of satellite B can be obtained according to the state equation. Specifically, this includes the following steps:

[0073] 1) Determine the position vector of satellite B relative to satellite A ( ) and velocity vector ( ).

[0074] 2) Position vector of star B relative to star A It can be represented as:

[0075] ;

[0076] Where d is the relative distance between stars A and B obtained in step five. α is the azimuth angle of star B relative to star A, and β is the elevation angle of star B relative to star A, both obtained from optical imaging analysis in step four.

[0077] The velocity vector of star B relative to star A It can be represented as:

[0078] ;

[0079] in , , These are the components of the velocity of star B relative to star A along the x, y, and z axes of star A's body coordinate system, derived from the inter-satellite relative motion model combined with observational data;

[0080] 3) Calculate the orbital length of star A (semi-major axis) into 6 elements. eccentricity ,inclination Right ascension of ascending node Angular distance from perigee And true near point angle () is converted into position and velocity vectors.

[0081] 4) Position vector of satellite B in the geocentric inertial coordinate system and velocity vector It can be represented as:

[0082] ;

[0083] ;

[0084] in , These are the position vector and velocity vector of satellite A in the geocentric inertial coordinate system, respectively, obtained from the conversion of the orbital elements of satellite A in step three.

[0085] 5) Convert the position and velocity vectors of satellite B into orbital root numbers.

[0086] The angular momentum vector of star B for:

[0087] ;

[0088] orbital inclination of star B Right ascension of ascending node for:

[0089] ;

[0090] in , , These are the angular momentum vectors. Components of the x, y, and z axes in a geocentric inertial coordinate system;

[0091] ;

[0092] if If it is negative, then ,like If non-negative, then ; The corrected right ascension of the ascending node.

[0093] The semi-major axis of B-star and eccentricity for:

[0094] ;

[0095] in Let μ be the orbital mechanical energy of star B, and μ be the standard gravitational constant of the central body (for Earth, μ = 3.986 × 10⁻⁶). 14 m³ / s²);

[0096] ;

[0097] ;

[0098] ;

[0099] in, This is the semi-aperture of star B.

[0100] Perigee angular distance of satellite B Peace Angle for:

[0101] ;

[0102] in The corrected right ascension of the ascending node; , , : Position components of satellite B along the x, y, and z axes in the geocentric inertial coordinate system; , : Velocity components of satellite B along the x and y axes in the geocentric inertial coordinate system;

[0103] if If it is negative, then ;

[0104] ;

[0105] ;

[0106] in The mean anomaly of satellite B represents the change in the satellite's orbital position over time.

[0107] It is the true anomaly angle of star B, which can be calculated using the following formula:

[0108] .

[0109] Let be the magnitude of the position vector of satellite B in the geocentric inertial coordinate system.

[0110] Step 7: Process the observation information using the extended Kalman filter algorithm, and perform time and measurement updates. Repeat steps 3 to 6, and perform step 7 to update the time after each loop to obtain the orbit determination result of the target satellite at any time.

[0111] After obtaining the approximate position coordinates and orbit of satellite B based on the state equation, the coordinates of satellite A at time t1 are determined using an iterative Kalman filter. Coordinates of satellite B at time t2 The distance between star A and star B , direction vector The solution is obtained iteratively through the following steps. (C is the iteration stopping error threshold, such as position error ≤ 10) -3 m and n are the number of iterations). It is the state vector of star B, which includes the position and velocity components in the geocentric inertial coordinate system.

[0112] 1) When At that time, the coordinates of star B were the initial coordinates. (Acquired by ground-based telemetry and control), then ;

[0113] 2) When hour, Then update the coordinates of star B based on the results of the previous iteration; The Kalman gain matrix at time n, used to balance prediction and observation errors; The observation vector at time n, containing microwave ranging values. and optical imaging orientation information .

[0114] 3) According to Whether the condition is true or false serves as a criterion for determining whether the iteration should stop.

[0115] The recursive equation for the extended Kalman filter is established as follows:

[0116] Updated in time: , ;

[0117] Measurement Update: , , ;

[0118] in, The optimal estimate of the state at time k; The optimal estimate of the state at time k+1; It is the predicted state quantity at time k+1 derived from the state model; These are discretized system equations; Let k be the system disturbance quantity at time k; This is the gain matrix at time k+1; The observation at time k+1; Let k→k+1 be the prediction error variance matrix; This is the observation matrix at time k+1; Let K+1 be the observation noise variance matrix. Let K be the variance matrix of the filtering error at time k; z is the identity matrix; z is the observation. Let be the covariance matrix of the process noise.

[0119] The extended Kalman filter algorithm is used to process the observation information and perform time and measurement updates to obtain the autonomous orbit determination results of the target satellite at any time.

[0120] Initial state estimation The initial orbit data is determined based on the satellite's initial orbit after launch and insertion into orbit, including the initial position coordinates of satellite B (acquired by ground control and telemetry) and its initial velocity (calculated from launch and orbit insertion parameters). For example, satellite B's initial position in the geocentric inertial coordinate system is... The initial velocity is ,but .

[0121] Initial covariance matrix The accuracy and initial measurement error settings are determined based on the track. For example, position errors ±Δx, ±Δy, and ±Δz correspond to the diagonal elements. , , Speed ​​error ±Δ ±Δ ±Δ Corresponding diagonal elements , , When the error is independent, the off-diagonal elements are 0.

[0122] Based on the observation data from ground stations, combined with the station's position in a certain coordinate system... ,speed Information is used to determine the orbit of satellite A using the extended Kalman filter method, thereby obtaining the position and velocity information of satellite A;

[0123] Given the orbital data of satellite A, including its position in a certain coordinate system. ,speed The information, combined with optical imaging and microwave measurement data, yielded the position and velocity information of satellite B relative to satellite A. This information was then transformed into a geocentric inertial coordinate system or other required coordinate system through coordinate transformation. Finally, the orbit of satellite B was determined using the extended Kalman filter method, thus fully determining the position and velocity of satellite B.

[0124] Example:

[0125] The binary constellation consists of a reference satellite (A) and a target satellite (B). Satellite A's microwave ranging equipment uses the Ka-band, achieving a ranging accuracy better than 1.5m@100km. Its optical imaging equipment boasts high resolution and a wide field of view, enabling clear imaging of Satellite B's optical beacon. The optical beacon possesses distinctive features such as high brightness, unique color, and shape, facilitating identification and positioning by Satellite A's imaging equipment. During the flight of the A and B binary constellations, their maximum relative distance must be maintained within a certain range to ensure that Satellite A's optical imaging payload can image and identify Satellite B.

[0126] The microwave ranging device on satellite A transmits microwave signals to satellite B every 10 seconds, and the imaging device takes pictures of satellite B at the same time intervals. This frequency setting takes into account both the timeliness of data and avoids the burden on the processing system caused by excessive data volume, while meeting the monitoring needs of satellite relative motion in lunar orbit.

[0127] The images of the optical beacon of satellite B taken by satellite A are preprocessed, including noise reduction and contrast enhancement. For example, median filtering is used to remove salt-and-pepper noise from the image, and histogram equalization is used to enhance the image contrast, making the optical beacon more clearly visible in the image.

[0128] Features of the optical beacon are extracted from the preprocessed image, such as shape features (the beacon's outline is extracted using edge detection algorithms to determine its geometric shape), color features (the beacon's color distribution is determined through color space analysis), and brightness features (the beacon's brightness distribution is analyzed through grayscale value statistics). This feature information helps to accurately identify the location of the optical beacon in the image.

[0129] Based on the imaging device's intrinsic parameters (such as focal length and pixel size) and extrinsic parameters (such as the camera's mounting attitude and position on the satellite), and combined with the pixel coordinates of the optical beacon's characteristics in the image, the direction vector of satellite B in the coordinate system of satellite A is calculated through geometric relationships such as perspective transformation. For example, if the intrinsic parameter matrix of the imaging device is known... extrinsic parameter matrix and (rotation and translation matrices), and the pixel coordinates of the optical beacon in the image. It can be done through formula (Q is the imaging device intrinsic parameter matrix, including focal length, pixel size, etc.) This yields the direction vector of star B in the camera coordinate system. Then, through coordinate transformation, it is converted to the coordinate system of satellite A, and then used for satellite formation positioning.

[0130] The time synchronization scheme for the two satellites employs a "measurement + calibration" approach. Because the measurement results contain random and dynamic errors, direct correction would lead to ranging time jitter. Therefore, the time calibration process requires preprocessing the measured synchronization error values ​​by filtering. The measurement and calibration process for the synchronization error between the two satellites is iterative.

[0131] The time difference elimination process consists of three steps, controlled by a state machine. Satellite A's communication unit is the clock reference source, and Satellite B's clock is aligned with that of Satellite A.

[0132] 1) Fast synchronization process:

[0133] The timeframe is the initial power-on phase.

[0134] Initial state: Satellite A's communication unit powers on, receives the time broadcast transmitted by the integrated power supply, and sends ranging frames to Satellite B's Ka-band communication payload at the designated integer seconds. Satellite B's Ka-band communication payload powers on, receives the radio frequency signals transmitted by Satellite A's communication unit, and first sends radio frequency signals to Satellite A's communication unit at the power-on time, taking the integer seconds as the starting time. After locking, it parses the Satellite A timestamp contained in the ranging frame and regards that time as the starting time of the integer seconds for this device.

[0135] Required time: The duration of this process is the synchronization time of the Ka-band payload frame of Satellite B, which is approximately 10~15 seconds.

[0136] Result: Time synchronization between satellites A and B was achieved to the distance / speed of light (0.6ms based on a distance of 200km).

[0137] 2) Clock calibration and tracking process:

[0138] Initial state: The time error between satellite A and satellite B has been reduced to within 1ms.

[0139] Operating status: Satellites A and B begin exchanging and receiving radio frequency signals during this phase, using a bidirectional comparison method to measure time difference. They receive 10-second time difference measurement data, smooth and filter it, and then use the smoothed time difference data to correct the time scale of Satellite B's Ka-band conduction payload.

[0140] Time synchronization result: Since the time error in this stage is approximately on the order of milliseconds or less, the clock correction error in this stage mainly originates from the error caused by relative motion within this time. Based on a Doppler range of 300kHz and a carrier frequency of 25GHz, the corresponding relative motion speed is calculated to be (300×10³×3×10⁸) / (25×10⁹) = 3600m / s. Therefore, after the fast synchronization in the previous stage, the transmission path difference caused by relative motion within a maximum time error of 1ms is 3.6m, and the calculated time error is approximately 3.6 / (3×10⁸) = 12ns. Thus, at the end of this state, the time synchronization error has been reduced to 10~15ns.

[0141] 3) Distance maintenance process:

[0142] Initial state: The time error between satellites A and B has been reduced to within 15ns.

[0143] Operating status: When the time error is within 15ns, the timing effect caused by relative motion can be ignored. At this time, the time error almost entirely originates from the clock difference and clock drift between the two satellites. Satellites A and B measure each other to obtain two pseudoranges. By solving the equations, the true distance and time difference between the two satellites can be obtained. To achieve synchronization, this time difference needs to be eliminated.

[0144] During this phase, both Satellite A and Satellite B employed a two-way comparison method for time difference measurement. Because the data may contain dynamic errors, measurements were taken every second, accumulating 10 seconds of time difference data each time. A filter was then used to smooth the 10-second data, and the smoothed time difference data was used to correct the time scale of Satellite B's Ka-band conduction payload. The correction was performed using loop control, adjusting the resolution to 10M / 256 = 0.14ns.

[0145] Required time: The initial stage of this maintenance process involves the transition between clock correction and ranging maintenance, requiring 2-3 minutes of convergence time. Afterwards, the loop is dynamically adjusted, and the error can be maintained within 3ns.

[0146] Time synchronization results: Since the strategy adopted is real-time adjustment, the synchronization error gradually decreases during the convergence process, and the final synchronization error is within 3ns.

[0147] 4) Power-on / off restart process: After each power-on / off, the above process needs to be repeated, that is, the ranging value needs a 2-3 minute re-stabilization process after each power-on / off.

[0148] 5) Self-timekeeping design:

[0149] During ranging and time synchronization, satellite A initially receives a 1pps time signal from the integrated power supply and uses that moment as the starting time for 1 second, subsequently maintaining self-time at a constant 10MHz. Satellite B, on the other hand, uses loop-filtered time difference data to correct its local clock. This clock is used only for self-time synchronization within seconds for satellite B, with an adjustment accuracy of 10MHz / 256 = 0.14ns. This adjustment method ensures that the relative error in self-time synchronization between satellites B and A is constantly being adjusted, achieving high accuracy without causing link interruptions and demonstrating feasibility. This process is repeated after each power-on to maintain self-time synchronization.

[0150] Given the orbital data of star A, including its position in a certain coordinate system. ,speed And orbital elements (such as semi-major axis, eccentricity, etc.).

[0151] Based on the position and velocity information of satellite B relative to satellite A obtained through Kalman filtering, it is transformed to a geocentric inertial coordinate system or other required coordinate system through coordinate transformation. For example, if the position vector of satellite A at a certain moment is known... Position vector of star B relative to star A (Estimated by Kalman filtering), then the position vector of star B in this coordinate system .

[0152] To derive the orbital root numbers of satellite B based on the six orbital root numbers of satellite A and the azimuth, distance, and velocity vectors of satellite B relative to satellite A, this invention requires a series of coordinate transformations and calculations. The detailed derivation process and mathematical expressions are as follows:

[0153] Convert the orbital root numbers of satellite A into position and velocity vectors. The orbital root numbers of satellite A include: semi-major axis. eccentricity ,inclination Right ascension of ascending node Angular distance from perigee And true near point angle .

[0154] First, calculate the semi-circular diameter of star A. :

[0155] ;

[0156] Then, calculate the true anomaly angle of star A. :

[0157] ;

[0158] ;

[0159] Next, calculate the position vector of star A. and velocity vector Representation in the orbital plane coordinate system:

[0160] ;

[0161] ;

[0162] ;

[0163] in, It is the standard gravitational constant.

[0164] Then, the position and velocity vectors of satellite A are transformed from the orbital plane coordinate system to the geocentric inertial coordinate system (ECI):

[0165] ;

[0166] ;

[0167] in , Let A be the position and velocity vector of satellite A in the geocentric inertial coordinate system.

[0168] The next step is to calculate the position and velocity vectors of satellite B relative to satellite A.

[0169] The next step is to calculate the position and velocity vector of satellite B in the geocentric inertial coordinate system.

[0170] This section details the process and mathematical expression for deriving the orbital element numbers of satellite B based on the six orbital elements of satellite A and the azimuth, distance, and velocity vectors of satellite B relative to satellite A. It should be noted that these calculations involve complex coordinate transformations and trigonometric function operations, and in practical applications, numerical methods may be necessary to solve them.

[0171] In iterative updates using Kalman filtering, the initial state is estimated. The initial orbital data is determined based on the satellite's initial orbital data after launch and insertion into orbit, including the initial position coordinates of satellite B (obtainable through ground-based telemetry and control) and its initial velocity (calculated based on the orbital parameters). For example, if the initial position of satellite B in a geocentric inertial coordinate system is known to be... The initial velocity is ,but .

[0172] Initial covariance matrix The initial covariance of the position coordinates is determined based on the accuracy of the satellite orbit determination and the initial measurement error. For example, the diagonal element values ​​can be set based on the position error range of ground control and telemetry; the covariance of the velocity components can be set based on the measurement error range of the injection velocity. Assuming the initial measurement error of the position coordinates is within... , , Within the range, the speed measurement error is within , , Within the range, then The diagonal elements can be set to , , , , , Off-diagonal elements can be set according to the correlation between errors (if the errors are independent, the off-diagonal elements are 0).

[0173] At some point Distance measured by microwave ranging equipment Imaging equipment captures images and analyzes them to obtain direction vectors. .

[0174] Calculate the predicted state based on the state equation. and predicted covariance For example, for the simple two-body motion assumption, if The matrix is ​​determined based on Newton's law of universal gravitation and the kinematic relationship of satellites. , ,in Determined based on the characteristics of the process noise (such as noise power spectral density).

[0175] Calculate Kalman gain and update the state estimate. Covariance .For example, , , ,here Determined based on the relationship between the observation equation and the state vector Determined based on the measurement noise characteristics of microwave ranging and optical imaging.

[0176] As time progresses, repeat the above steps to continuously update the state estimate of satellite B.

[0177] Kalman filtering is a recursive algorithm used to estimate the state of linear dynamic systems. For nonlinear systems, extended Kalman filtering (EKF) or unscented Kalman filtering (UKF) can be used.

[0178] Through these steps, this invention can calculate the number of orbital elements for satellite B and use Kalman filtering for iterative updates to the state estimate. This process requires combining a specific system model with observational data.

[0179] This invention also provides a device for determining the orbit of a binary satellite formation based on inter-satellite relative optical imaging and microwave measurement, comprising the following modules:

[0180] The state equation establishment module establishes satellite motion equations based on orbital dynamics, thereby obtaining the state equations;

[0181] The observation equation establishment module establishes observation equations based on inter-satellite relative optical imaging and microwave measurements.

[0182] The measurement module determines the precise orbit of the reference satellite A.

[0183] The imaging module uses reference star A to perform optical imaging of target star B, obtaining the orientation information of star B relative to star A.

[0184] The ranging module performs microwave ranging between the reference satellite A and the target satellite B to obtain relative distance information.

[0185] The calculation module calculates the orbital parameters of satellite B. Based on the observed coordinates of satellite B relative to satellite A, it obtains the position coordinates and satellite orbit of satellite B according to the state equation.

[0186] The orbit determination module uses the extended Kalman filter algorithm to process observation information and performs time and measurement updates to obtain the orbit determination result of the target satellite at any time.

[0187] This invention also provides an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, it implements the steps of the above-described method for determining the orbit of a binary satellite formation based on inter-satellite relative optical imaging and microwave measurement.

[0188] This invention also provides a non-transitory computer-readable storage medium storing a computer program thereon, which, when executed by a processor, implements the steps of the above-described method for determining the orbit of a dual-satellite formation based on inter-satellite relative optical imaging and microwave measurement.

[0189] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code. The solutions in the embodiments of the present invention can be implemented using various computer languages, such as the object-oriented programming language Java and the interpreted scripting language JavaScript.

[0190] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart illustrations and / or block diagrams. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0191] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0192] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0193] Although preferred embodiments of the invention have been described, those skilled in the art, upon learning the basic inventive concept, can make other changes and modifications to these embodiments. Therefore, the appended claims are intended to be interpreted as including both the preferred embodiments and all changes and modifications falling within the scope of the invention.

[0194] Obviously, those skilled in the art can make various modifications and variations to this invention without departing from its spirit and scope. Therefore, if these modifications and variations fall within the scope of the claims of this invention and their equivalents, this invention also intends to include these modifications and variations.

Claims

1. A method for determining the orbit of a dual-satellite formation based on inter-satellite relative optical imaging and microwave measurement, characterized in that, Comprising the following steps: Step one: establishing satellite motion equation based on orbit dynamics, and then obtaining state equation; Step two: establishing observation equation based on inter-satellite relative optical imaging and microwave measurement; Step three: measuring accurate orbit of reference star A; Step four: reference star A imaging target star B to obtain azimuth information of B relative to A; Step five: reference star A and target star B microwave ranging to obtain relative distance information; Step six: calculating orbit parameters of B, obtaining position coordinates of B and satellite orbit according to observation coordinate values of B relative to A and state equation; Step seven: processing observation information by extended Kalman filter algorithm, and performing time update and measurement update, repeating steps three to six, and performing time update and measurement update by step seven after each cycle to obtain orbit determination result of target satellite at any time. 2.The method of determining the orbit of a dual-satellite formation satellite based on inter-satellite relative optical imaging and microwave measurement according to claim 1, characterized in that, The step one comprises: Let the state vector where is the position coordinate of the satellite in a certain coordinate system, is the corresponding velocity component, t denotes time, and the superscript T denotes the transpose of a matrix; Assuming that the motion of the satellite is consistent with the dynamic model, the state equation can be expressed as where is the state transition function, is the process noise, which is usually assumed to be Gaussian white noise, and the covariance matrix is , Δt is the time interval between adjacent state updates, and for a simple two-body problem, can be derived according to Newton's law of motion and the law of universal gravitation. 3.The method of determining the orbit of a dual-satellite formation satellite based on inter-satellite relative optical imaging and microwave measurement according to claim 1, characterized in that, The step two comprises: Establishing A observation equation: Observation vector The observation equation is where is the observation matrix, is the observation noise, also assumed to be Gaussian white noise with covariance matrix The range equation is expressed as where is the position coordinates of the ground station, is the ranging noise; Establishing B observation equation: Observation vector The observation equation is where is the observation matrix, is the B-star observation noise, also assumed to be Gaussian white noise with covariance matrix ; the range equation can also be expressed as where is the position coordinates of the A-star, is the range noise. 4.The method of determining the orbit of a dual-satellite formation satellite based on inter-satellite relative optical imaging and microwave measurement according to claim 1, characterized in that, The step three comprises: The ground station receives the A-star observation data and calculates the distance between the A-star and the station according to the signal propagation time and other orbit information, wherein t represents time; the orbit of the A-star is obtained, including semi-major axis , eccentricity , inclination , ascending node right ascension , near point angle distance and true near point angle . 5.The method of determining the orbit of a dual-satellite formation satellite based on inter-satellite relative optical imaging and microwave measurement according to claim 1, characterized in that, The step four comprises: The imaging device of the A star takes pictures of the optical beacon of the B star in the same time interval, processes the pictures, identifies the pixel position and pose of the optical beacon in the picture, combines the pose information of the A star itself and the internal and external parameters of the imaging device, and calculates the directional vector of the B star relative to the A star . 6.The method of determining the orbit of a dual-satellite formation satellite based on inter-satellite relative optical imaging and microwave measurement according to claim 1, characterized in that, The step five comprises: A star microwave ranging equipment to B star at certain time interval to send microwave signal, and receives the signal of B star emission, through two-way contrast ranging method, according to the propagation time of microwave signal calculates the distance value between A star and B star , c is the speed of light, It is the two-way propagation time, that is, the total time interval of signal from A star emission to receiving B star retransmission signal (remove internal electronic response timing). 7.The method of determining the orbit of a dual-satellite formation satellite based on inter-satellite relative optical imaging and microwave measurement according to claim 1, characterized in that, The step six comprises the following steps: 1) determining a position vector of satellite B relative to satellite A and a velocity vector ; 2) Position vector of B star relative to A star is represented as: ； where d is the relative distance of the A, B stars obtained in step five , is the azimuth angle of the B star relative to the A star, and β is the elevation angle of the B star relative to the A star. The velocity vector of B star relative to A star is represented as: ； wherein , , are the components of the velocity of B star relative to A star in the x, y, z axes of the body coordinate system of A star, derived from the interstellar relative motion model combined with observation data; 3) converting the orbital elements of the Astar, 6 in number, semi-major axis , eccentricity , inclination , longitude of the ascending node , argument of perigee and true anomaly , into position and velocity vectors; 4) Position vector of B-star in the geocentric inertial coordinate system and velocity vector is expressed as: ； ； wherein , are the position vector and velocity vector of A-star in the Earth-centered inertial coordinate system, respectively, obtained from the conversion of the 6 orbital elements of A-star in step three 5) converting position and velocity vector of B into 6 orbit elements: The momentum vector of B-star is: ； The orbit inclination of B star and the longitude of the ascending node is: ； wherein , , are the components of the angular momentum vector in the x, y, z axes of the geocentric inertial coordinate system. ； If is negative, then , if is non-negative, then ; is the corrected right ascension of the ascending node; The semi-major axis of B-star and eccentricity is: ； wherein is the orbital mechanical energy of the B-star, μ is the standard gravitational constant of the central body, μearth= 3.986 x 10 14 m³ / s²; ； ； wherein is the semi-major axis of the B-star orbit, describing the shape of the orbit; ； The argument of perigee of B-star And the argument of perihelion Is: ; , , : Position components of B-star in the x, y, z axes of the geocentric inertial coordinate system; , : Velocity components of B-star in the x, y axes of the geocentric inertial coordinate system; If is negative, then ; ； ； wherein is the mean anomaly of the satellite B, used to characterize the position of the satellite on the orbit as a function of time; is the true anomaly of the B-star, calculated by the following equation: ； is the magnitude of the position vector of B-star in the geocentric inertial coordinate system. 8.The method of determining the orbit of a dual-satellite formation satellite based on inter-satellite relative optical imaging and microwave measurement according to claim 1, characterized in that, The step seven comprises: After obtaining the approximate position coordinates and orbit of satellite B based on the state equation, the coordinates of satellite A at time t1 are determined using an iterative Kalman filter. Coordinates of satellite B at time t2 The distance between star A and star B , direction vector The solution is obtained iteratively through the following steps. , Let be the state vector of satellite B, containing position and velocity components in the geocentric inertial coordinate system; where C is the set iteration stopping error threshold, and n is the number of iterations. 1) when B-star coordinates are initial coordinates then ; 2) when , , update the B-star coordinate according to the result of the previous iteration; : Kalman gain matrix at time n, used to balance the prediction error and observation error; : observation vector at time n, containing microwave ranging value and optical imaging azimuth information ; 3) According to whether or not it is true, as a condition for judging whether or not iteration is stopped; Establishing recursive equation of extended Kalman filter, and the extended Kalman filter is as follows: Time update: , ; Measurement update: , ， ； wherein, is the optimal state estimate at time k; is the optimal state estimate at time k+1; is the predicted state quantity at time k+1 calculated from the state model; is the discretized system equation; is the system disturbance quantity at time k; is the gain matrix at time k+1; is the observation quantity at time k+1; is the prediction error variance matrix at time k which is transferred to time k+1; is the observation matrix at time k+1; is the observation noise variance matrix at time k+1; is the filter estimation error variance matrix at time k; is the identity matrix; is the covariance matrix of the process noise; Processing observation information by extended Kalman filter algorithm, and performing time update and measurement update to obtain autonomous orbit determination result of target satellite at any time.

9. A device for determining the orbit of a dual-satellite formation satellite based on inter-satellite relative optical imaging and microwave measurement, characterized in that, Comprising the following modules: State equation establishing module, establishing satellite motion equation based on orbit dynamics, and then obtaining state equation; Observation equation establishing module, establishing observation equation based on inter-satellite relative optical imaging and microwave measurement; Measuring module, measuring accurate orbit of reference star A; Imaging module, reference star A imaging target star B to obtain azimuth information of B relative to A; Ranging module, reference star A and target star B microwave ranging to obtain relative distance information; Calculating module, calculating orbit parameters of B, obtaining position coordinates of B and satellite orbit according to observation coordinate values of B relative to A and state equation; Orbit determination module, processing observation information by extended Kalman filter algorithm, and performing time update and measurement update to obtain orbit determination result of target satellite at any time.

10. An electronic device comprising a memory, a processor, and a computer program stored on the memory and executable on the processor, characterized in that, The processor implements the steps of the double-satellite formation satellite orbit determination method based on inter-satellite relative optical imaging and microwave measurement according to any one of claims 1 to 8 when executing the program. 11.A non-transitory computer-readable storage medium having stored thereon a computer program. The computer program implements the steps of the double-satellite formation satellite orbit determination method based on inter-satellite relative optical imaging and microwave measurement according to any one of claims 1 to 8 when executed by the processor.

Images