A method for short-arc orbit determination of spacecraft based on visual measurement

Abstract

A method for short arc orbit determination of spacecraft based on visual measurement. The present invention relates to a method for short arc orbit determination of spacecraft based on visual measurement. The purpose of the present invention is to solve the problem that the existing short arc orbit determination method is greatly affected by the length of observation time. The process is: 1: obtain the orbit parameters of the observation spacecraft and the relevant parameters of the camera; 2: determine the timestamp interval for capturing the target spacecraft image; 3: obtain the pixels where the target spacecraft image is located for processing to obtain the processed image; 4: obtain the image after grayscale superposition; 5: solve the direction vector of the target spacecraft relative to the observation spacecraft in the camera coordinate system; 6: solve the direction vector of the target spacecraft relative to the observation spacecraft in the geocentric inertial coordinate system; 7: use the position vector sequence of the observation spacecraft and the direction vector sequence of the target spacecraft in the geocentric inertial coordinate system as input to optimize the six orbital numbers of the target spacecraft. The present invention relates to the field of spacecraft orbit parameter estimation.

Description

Technical Field

[0001] The present invention relates to the field related to spacecraft orbit parameter estimation, and in particular to a spacecraft short arc orbit determination method based on visual measurement. Background Art

[0002] Short-arc orbit determination, also known as initial orbit calculation or orbit calculation, refers to the calculation of the orbit in the sense of two bodies based on a small amount of observation data from a single station and a single loop without prior information. The orbit calculation of angle measurement data is a classic problem in celestial mechanics. The research on orbit calculation methods based on measured data has been going on for more than two hundred years. With the emergence of new celestial bodies and new observation methods, new needs and difficulties have also arisen, which need to be further explored.

[0003] At present, the problem of initial orbit calculation has never been satisfactorily solved. It is fundamentally caused by the intrinsic ill-conditioning of the equations to be solved, that is, the short-arc orbit determination problem is ill-posed. Faced with this difficulty, in order to obtain better results, we can only rely on collecting arc segments as long as possible. For optical observations, arc segments of 5-6 minutes are sufficient to obtain better results. In practice, it is generally required to collect arc segments of at least 3 minutes continuously. For new sky survey observation methods, its observation arc length is much less than 3 minutes. In this case, the traditional orbit determination method is no longer applicable. The extremely short arc orbit determination becomes particularly important, and its problems are more prominent. There is currently no strict definition of the specific arc length of the extremely short arc. Some literature uses 2 minutes as the standard, some use 1 minute, and some are even shorter. In most cases, the observation arc segment that cannot obtain reasonable orbit determination results using classical methods is called an extremely short arc, to distinguish it from the short arc orbit determination in the traditional sense. Summary of the invention

[0004] The purpose of the present invention is to solve the problem that the existing short-arc orbit determination method is greatly affected by the length of observation time, and to propose a spacecraft short-arc orbit determination method based on visual measurement.

[0005] A spacecraft short arc orbit determination method based on visual measurement. The specific process is as follows:

[0006] Step 1: Obtain the orbital parameters of the observation spacecraft and relevant parameter information of the camera;

[0007] When the target spacecraft appears in the field of view of the observation spacecraft, the grayscale image of the target spacecraft is continuously captured by the onboard camera of the observation spacecraft;

[0008] Step 2: Determine the time stamp interval for the onboard camera on the observation spacecraft to select and capture the target spacecraft image;

[0009] The position coordinates of the observed spacecraft corresponding to the timestamp in the geocentric inertial coordinate system are calculated according to the six orbital elements of the observed spacecraft and the timestamp interval of the grayscale image captured by the target spacecraft;

[0010] Step 3: Obtain the pixels where the target spacecraft image captured by the observation spacecraft is located, and modify the grayscale values ​​of the pixels other than the pixels where the target spacecraft is located to 0 to obtain a processed image;

[0011] Step 4: grayscale overlay the processed image obtained in step 3 according to the time stamp interval in step 2 to obtain a grayscale overlay image;

[0012] Step 5: Calculate the direction vector of the target spacecraft relative to the observation spacecraft in the camera coordinate system based on the grayscale superimposed image;

[0013] Step 6: Calculate the direction vector of the target spacecraft relative to the observation spacecraft in the geocentric inertial coordinate system;

[0014] Step 7: Take the position vector sequence of the observation spacecraft and the direction vector sequence of the target spacecraft in the geocentric inertial coordinate system as input, and iteratively optimize the six orbital elements of the target spacecraft through the particle swarm algorithm.

[0015] Preferably, the orbital parameters of the observed spacecraft in step 1 are six orbital parameters, including the orbital semi-major axis, eccentricity, right ascension of the ascending node, orbital inclination, true anomaly angle, and perigee amplitude.

[0016] Preferably, the relevant parameter information of the camera in step 1 is focal length, resolution, field of view, and camera orientation.

[0017] Preferably, in step 3, the pixels where the target spacecraft image captured by the observation spacecraft is located are obtained, and the gray values ​​of the pixels other than the pixels where the target spacecraft is located are modified to 0 to obtain a processed image; the specific process is:

[0018] The pixels where the target spacecraft image captured by the observation spacecraft is located are obtained through the target segmentation algorithm or Photoshop image processing software, and the gray values ​​of the pixels except the pixels where the target spacecraft is located are modified to 0 to obtain the processed image.

[0019] Preferably, in step 5, the direction vector of the target spacecraft relative to the observation spacecraft in the camera coordinate system is calculated based on the grayscale superimposed image; the specific process is:

[0020] The initial frame target spacecraft mass center position and the last frame target spacecraft mass center position are manually marked on the grayscale overlay image, and the pixel coordinates (u i′ ,v i′ ),i′=1,2,…,N;

[0021] Where N represents the number of timestamps;

[0022] The projection relationship is used to obtain the direction vector corresponding to each pixel in the two-dimensional image, that is, the unit direction vector pointing from the coordinate origin to the center of mass of the target spacecraft in the camera coordinate system.

[0023] Preferably, the unit direction vector calculation process is as follows:

[0024] According to the projection relationship of camera imaging, the following formula is obtained:

[0025]

[0026]

[0027] Among them, w x and w y are the horizontal and vertical field of view angles, respectively; X and Y are the width and height of the camera target surface, respectively; f′ is the focal length of the camera;

[0028] Further, according to the position of the pixel point, we can get

[0029]

[0030]

[0031] Among them, (x, y) is the coordinate of the center of mass of the target spacecraft in the image physical coordinate system, (u, v) is the coordinate of the center of mass of the target spacecraft in the pixel coordinate system, and p x and p y are the number of pixels in each row and column of the image respectively;

[0032] Because (w x ,w y ), (u,v) and p x 、p y Known, then combining the above formulas we can get:

[0033]

[0034]

[0035] Then there is

[0036]

[0037] Finally, the vector Normalize to get the unit direction vector corresponding to a frame of image;

[0038] Calculate each frame of image to obtain the unit direction vector corresponding to each frame of image.

[0039] Preferably, in step 6, the direction vector of the target spacecraft relative to the observation spacecraft in the geocentric inertial coordinate system is calculated; the specific process is:

[0040] Define the origin of the observation spacecraft orbital coordinate system as the center of the earth, the x-axis points from the center of the earth to the spacecraft, the z-axis is the normal direction of the orbital plane, and the y-axis direction is determined according to the right-hand rule;

[0041] (1) Obtain the transformation matrix from the camera coordinate system to the observation spacecraft orbit coordinate system:

[0042] Calculate the transfer matrix from the camera coordinate system to the observation spacecraft orbit coordinate system, denoted as C o1 ;

[0043] (2) Obtain the transformation matrix from the observation spacecraft orbital coordinate system to the geocentric inertial coordinate system:

[0044] True anomaly of the observing spacecraft Argument of perigee Calculate the latitude argument of the observing spacecraft

[0045]

[0046] Then the transformation matrix from the observation spacecraft orbital coordinate system to the geocentric inertial coordinate system is

[0047]

[0048] in, is the right ascension of the ascending node of the observing spacecraft, To observe the orbital inclination of the spacecraft, is the latitude argument of the observing spacecraft;

[0049] The following equations represent the rotation angle α around the x-axis, y-axis, and z-axis of the coordinate system respectively:

[0050]

[0051] R x (α) is the rotation angle θ around the x-axis of the coordinate system, R y (α) is the rotation angle α around the y-axis of the coordinate system, R z (α) is the rotation angle α around the z-axis of the coordinate system, α is the rotation angle, R x (-90°) is the rotation angle of -90° around the x-axis of the coordinate system, R z (90°) is the rotation angle of 90° around the z-axis of the coordinate system. is the rotation angle around the z-axis of the coordinate system Rotation around the x-axis of the coordinate system Rotation around the z-axis of the coordinate system

[0052] (3) Obtain the transformation matrix from the camera coordinate system to the geocentric inertial coordinate system:

[0053] Assuming that the unit direction vector corresponding to each frame of the target spacecraft is v in the camera coordinate system and v' in the geocentric inertial coordinate system, we can get

[0054]

[0055] Where T is the transpose.

[0056] Preferably, in step 7, the position vector sequence of the observation spacecraft and the direction vector sequence of the target spacecraft in the geocentric inertial coordinate system are used as input, and the six orbital elements of the target spacecraft are iteratively optimized by a particle swarm algorithm; the specific process is:

[0057] 1) First layer optimization

[0058] From the target spacecraft's (a, e, M0), we get the M of the jth frame j ;

[0059] M j =M0+nΔt j

[0060] Where Δt j represents the time difference between the jth frame and the initial frame, current n represents the angular velocity of the target spacecraft in the current state; a represents the semi-major axis of the target spacecraft; e represents the eccentricity of the target spacecraft; M0 represents the true anomaly angle of the target spacecraft at the initial moment;

[0061] Further calculation is performed to obtain the true near point angle f corresponding to the jth frame j and the eccentric anomaly angle E j , there is the following relationship;

[0062] M j =E j -esinE j

[0063]

[0064] Since it is impossible to pass M j Direct calculation of E j , so an iterative calculation method is used to repeatedly iterate the following formula:

[0065]

[0066] Until the following is satisfied

[0067] |M j +esinE j -E j |<ε

[0068] Among them, ε is the threshold, which can be set to 10 -9 ;

[0069] Further calculation:

[0070] r j =a(1-ecosE j )

[0071]

[0072]

[0073]

[0074] in, represents the geocentric position vector of the target spacecraft in the jth frame in the geocentric inertial coordinate system, r j represents the geocentric position vector modulus of the target spacecraft in the jth frame in the geocentric inertial coordinate system, ρ j represents the slant distance of the target spacecraft in the jth frame in the geocentric inertial coordinate system, represents the geocentric position vector of the observing spacecraft in the geocentric inertial coordinate system, R j represents the geocentric position vector modulus of the observing spacecraft in the geocentric inertial coordinate system, represents the unit direction vector of the target spacecraft relative to the observing spacecraft in the geocentric inertial coordinate system; β j express and The angle between , · represents the vector multiplication;

[0075] For the jth frame and the kth frame, k>j, we get (r k ,r j ) and (f k ,f j );

[0076] Therefore, the fitness function is defined as follows:

[0077]

[0078] in, represents the geocentric position vector of the target spacecraft in the geocentric inertial coordinate system in the kth frame, r k represents the geocentric position modulus of the target spacecraft in the kth frame in the geocentric inertial coordinate system, f kIndicates the true anomaly angle corresponding to the kth frame in the geocentric inertial coordinate system, and N indicates the number of timestamps;

[0079] The fitness function value is taken as the minimum corresponding to the semi-major axis a of the target spacecraft, the eccentricity e of the target spacecraft, and the true anomaly M of the target spacecraft at the initial moment as the prediction value 1;

[0080] 2) Second layer optimization

[0081] According to the predicted value of (a, e, M), the true near point angle f corresponding to the jth frame is obtained j and the eccentric anomaly angle E j

[0082] There are the following relationships;

[0083] M j =E j -esinE j

[0084]

[0085] Since it is impossible to pass M j Direct calculation of E j , so an iterative calculation method is used to repeatedly iterate the following formula:

[0086]

[0087] Until the following is satisfied

[0088] |M j +esinE j -E j |<ε

[0089] Among them, ε is the threshold, which can be set to 10 -9 ;

[0090] Further calculation:

[0091] k j =a(1-e 2 ) / (1+ecosf j )

[0092]

[0093]

[0094]

[0095] Among them, k jrepresents the position vector modulus of the target spacecraft in the jth frame in the geocentric inertial coordinate system, Ω represents the right ascension of the ascending node of the target spacecraft, ω represents the perigee argument of the target spacecraft, i represents the orbital inclination of the target spacecraft, and f j Indicates the true anomaly angle corresponding to the jth frame, represents the geocentric position vector of the observing spacecraft in the geocentric inertial coordinate system, represents the position vector of the target spacecraft in the jth frame in the geocentric inertial coordinate system, · represents the vector multiplication, θ j It represents the angular difference between the estimated value and the true value of the direction vector of the target spacecraft with respect to the observation spacecraft in the jth frame in the geocentric inertial coordinate system;

[0096] The fitness function is defined as follows:

[0097]

[0098] The fitness function value is taken as the minimum corresponding to the ascending node right ascension Ω, the perigee argument ω, and the orbital inclination i of the target spacecraft as the prediction value 2;

[0099] The six orbital numbers of the target spacecraft are obtained based on the predicted values ​​1 and 2.

[0100] The beneficial effects of the present invention are:

[0101] The invention discloses a spacecraft short arc orbit determination method based on a particle swarm algorithm, and belongs to the technical field of spacecraft orbit measurement.

[0102] 1. The present invention proposes a method for short-arc orbit determination of spacecraft based on particle swarm algorithm. Compared with other existing methods, this method can achieve orbit determination through observation information in a very short period of time (seconds).

[0103] 2. This method is less affected by the short arc observation time, and the trade-off between accuracy and real-time performance can be achieved by adjusting the particle swarm algorithm parameters. This method is more dependent on the accuracy of the target spacecraft center of mass position at the initial moment and the target spacecraft center of mass position at the end moment. When the marking is accurate, the observation time has almost no effect. If the marking is inaccurate, the longer the observation time, the smaller the impact. BRIEF DESCRIPTION OF THE DRAWINGS

[0104] Figure 1 It is a flow chart of the present invention;

[0105] Figure 2 To obtain the position relationship diagram of the observation spacecraft and the target spacecraft, represents the position vector of the target spacecraft in the geocentric inertial coordinate system, ρ represents the slant range of the target spacecraft in the geocentric inertial coordinate system, represents the position vector of the observing spacecraft in the geocentric inertial coordinate system, Represents the direction vector of the target spacecraft relative to the observing spacecraft in the geocentric inertial coordinate system;

[0106] Figure 3 This is a schematic diagram of image overlay;

[0107] Figure 4 It is a projection diagram;

[0108] Figure 5 Schematic diagram of the two-dimensional projection principle. DETAILED DESCRIPTION

[0109] Specific implementation method 1: This implementation method is a method for determining the short arc orbit of a spacecraft based on visual measurement. The specific process is as follows:

[0110] Step 1: Obtain the orbital parameters of the observation spacecraft and the relevant parameter information of the camera (assuming that the orbital parameters of the observation spacecraft (six orbital numbers) and the relevant parameter information of the camera are known);

[0111] When the target spacecraft appears in the field of view of the observation spacecraft, the grayscale image of the target spacecraft is continuously captured by the onboard camera of the observation spacecraft;

[0112] Step 2: Determine the timestamp interval for the onboard camera on the observation spacecraft to capture the target spacecraft image, such as 2s;

[0113] If the total time of capturing the grayscale image of the target spacecraft is long (total time is greater than or equal to 60s), the timestamp interval can be appropriately increased to 2-5 times the original value (reasonable selection can be made according to the total time); if the total time of capturing the grayscale image of the target spacecraft is short (total time is less than 3s), only the image and timestamp data corresponding to the initial frame and the last frame of the grayscale image of the target spacecraft are selected;

[0114] The position coordinates of the observed spacecraft corresponding to the timestamp in the geocentric inertial coordinate system are calculated according to the orbital parameters of the observed spacecraft and the timestamp interval of the grayscale image captured by the target spacecraft;

[0115] Step 3: Obtain the pixels where the target spacecraft image captured by the observation spacecraft is located, and modify the grayscale values ​​of the pixels other than the pixels where the target spacecraft is located to 0 to obtain a processed image;

[0116] Step 4: grayscale overlay the processed image (including the target spacecraft and the background image) obtained in step 3 according to the time stamp interval in step 2 to obtain a grayscale overlay image;

[0117] Step 5: Calculate the direction vector of the target spacecraft relative to the observation spacecraft in the camera coordinate system based on the grayscale superimposed image;

[0118] Step 6: Calculate the direction vector of the target spacecraft relative to the observation spacecraft in the geocentric inertial coordinate system;

[0119] Step 7: Take the position vector sequence of the observation spacecraft (obtained in step 2) and the direction vector sequence of the target spacecraft in the geocentric inertial coordinate system (obtained in step 6) as input, and iteratively optimize the six orbital elements of the target spacecraft through the particle swarm algorithm.

[0120] Specific implementation method 2: The difference between this implementation method and specific implementation method 1 is that the orbital parameters of the observed spacecraft in step 1 are six orbital parameters, including the semi-major axis of the orbit, eccentricity, right ascension of the ascending node, orbital inclination, true anomaly angle, and perigee amplitude.

[0121] The other steps and parameters are the same as those in the first embodiment.

[0122] Specific implementation method three: This implementation method is different from specific implementation methods one or two in that the relevant parameter information of the camera in step 1 is focal length, resolution, field of view, and camera orientation.

[0123] The other steps and parameters are the same as those in the first or second embodiment.

[0124] Specific implementation method 4: This implementation method is different from any one of specific implementation methods 1 to 3 in that: in step 3, the pixels where the target spacecraft image captured by the observation spacecraft is located are obtained, and the gray values ​​of the pixels other than the pixels where the target spacecraft is located are modified to 0 to obtain the processed image; the specific process is:

[0125] The pixels where the target spacecraft image captured by the observation spacecraft is located are obtained through the target segmentation algorithm or Photoshop image processing software, and the gray values ​​of the pixels except the pixels where the target spacecraft is located are modified to 0 to obtain the processed image.

[0126] The other steps and parameters are the same as those in Specific Embodiments 1 to 3.

[0127] Specific implementation method 5: This implementation method is different from any one of the specific implementation methods 1 to 4 in that: in step 5, the direction vector of the target spacecraft relative to the observation spacecraft in the camera coordinate system is solved based on the grayscale superimposed image; the specific process is:

[0128] Since the distance between the target spacecraft and the observation spacecraft is extremely small relative to their respective semi-major axes, it can be approximately considered that the target spacecraft is moving in a uniform straight line within the image sequence captured by the camera.

[0129] The initial frame target spacecraft mass center position and the last frame target spacecraft mass center position are manually marked on the grayscale overlay image, and the pixel coordinates (ui′ ,v i′ ),i′=1,2,…,N;

[0130] Where N represents the number of timestamps (images);

[0131] According to the camera's resolution, field of view and other information, the projection relationship is used to obtain the direction vector corresponding to each pixel in the two-dimensional image, that is, the unit direction vector pointing from the coordinate origin to the center of mass of the target spacecraft in the camera coordinate system, such as Figure 4 shown.

[0132] The other steps and parameters are the same as those in Specific Embodiments 1 to 4.

[0133] Specific implementation method 6: This implementation method is different from the specific implementation methods 1 to 5 in that: the unit direction vector calculation process is as follows:

[0134] according to Figure 5 The projection relationship of the camera imaging shown is as follows:

[0135]

[0136]

[0137] Among them, w x and w y are the horizontal and vertical field of view angles, respectively; X and Y are the width and height of the camera target surface, respectively; f′ is the focal length of the camera;

[0138] Further, according to the position of the pixel point, we can get

[0139]

[0140]

[0141] Among them, (x, y) is the coordinate of the center of mass of the target spacecraft in the image physical coordinate system, (u, v) is the coordinate of the center of mass of the target spacecraft in the pixel coordinate system, and p x and p y are the number of pixels in each row and column of the image respectively;

[0142] Because (w x ,w y ), (u,v) and p x 、p y Known, then combining the above formulas we can get:

[0143]

[0144]

[0145] Then there is

[0146]

[0147] Finally, the vector Normalize to get the unit direction vector corresponding to a frame of image;

[0148] Calculate each frame of image to obtain the unit direction vector corresponding to each frame of image.

[0149] The other steps and parameters are the same as those in Specific Implementation Methods 1 to 5-1.

[0150] Specific implementation method 7: This implementation method is different from any one of the specific implementation methods 1 to 6 in that: in step 6, the direction vector of the target spacecraft relative to the observation spacecraft in the geocentric inertial coordinate system is calculated; the specific process is:

[0151] Calculate the unit direction vector corresponding to each frame timestamp, and calculate the transfer matrix from the camera coordinate system to the geocentric inertial coordinate system, and transform the direction vector into the geocentric inertial coordinate system;

[0152] The direction vector sequence obtained above is based on the camera coordinate system and needs to be converted to the geocentric inertial coordinate system. In general, first, the transformation matrix from the camera coordinate system to the observation spacecraft orbital coordinate system is obtained based on the camera placement and other information, and then the transformation matrix from the spacecraft orbital coordinate system to the geocentric inertial coordinate system is calculated based on the six orbital numbers of the observation spacecraft. Using the above transformation matrix, the direction vector is converted to the geocentric inertial coordinate system as the reference system.

[0153] The specific process is as follows:

[0154] Define the origin of the observation spacecraft orbital coordinate system as the center of the earth, the x-axis points from the center of the earth to the spacecraft, the z-axis is the normal direction of the orbital plane, and the y-axis direction is determined according to the right-hand rule;

[0155] (1) Obtain the transformation matrix from the camera coordinate system to the observation spacecraft orbit coordinate system:

[0156] The transfer matrix from the camera coordinate system to the observation spacecraft orbit coordinate system is calculated using the camera's swing attitude and other information, denoted as C o1 ;

[0157] (2) Obtain the transformation matrix from the observation spacecraft orbital coordinate system to the geocentric inertial coordinate system:

[0158] True anomaly of the observing spacecraft Argument of perigee Calculate the latitude and angle of the observing spacecraft

[0159]

[0160] Then the transformation matrix from the observation spacecraft orbital coordinate system to the geocentric inertial coordinate system is

[0161]

[0162] in, is the right ascension of the ascending node of the observing spacecraft, To observe the orbital inclination of the spacecraft, is the latitude argument of the observing spacecraft;

[0163] The following equations represent the rotation angle α around the x-axis, y-axis, and z-axis of the coordinate system respectively:

[0164]

[0165] R x (α) is the rotation angle θ around the x-axis of the coordinate system, R y (α) is the rotation angle α around the y-axis of the coordinate system, R z (α) is the rotation angle α around the z-axis of the coordinate system, α is the rotation angle, R x (-90°) is the rotation angle of -90° around the x-axis of the coordinate system, R z (90°) is the rotation angle of 90° around the z-axis of the coordinate system. is the rotation angle around the z-axis of the coordinate system Rotation around the x-axis of the coordinate system Rotation around the z-axis of the coordinate system

[0166] (3) Obtain the transformation matrix from the camera coordinate system to the geocentric inertial coordinate system:

[0167] Assuming that the unit direction vector corresponding to each frame of the target spacecraft is v in the camera coordinate system and v' in the geocentric inertial coordinate system, we can get

[0168]

[0169] Where T is the transpose.

[0170] The other steps and parameters are the same as those in Specific Implementations 1 to 6-1.

[0171] Specific implementation eight: This implementation differs from any one of specific implementations one to seven in that: in step 7, the position vector sequence of the observation spacecraft (obtained in step 2) and the direction vector sequence of the target spacecraft in the geocentric inertial coordinate system (obtained in step 6) are used as inputs, and the six orbital elements of the target spacecraft are iteratively optimized by a particle swarm algorithm; the specific process is:

[0172] For orbit calculation, optimizing 6 parameters at the same time makes the dimension of solution space too high, which is not conducive to solving. For this reason, 3 orbital six numbers (a, e, M) are selected as the optimization variables, so that the optimization result is directly the orbital quantity, which achieves the decoupling of calculation results and observation data while keeping the problem dimension low. After determining (a, e, M), these three values ​​are used to continue iterative calculation (i, ω, Ω) using particle swarm algorithm. The details are as follows:

[0173] 1) First layer optimization

[0174] From the target spacecraft's (a, e, M0), we get the M of the jth frame j ;

[0175] M j =M0+nΔt j

[0176] Where Δt j represents the time difference between the jth frame and the initial frame, current n represents the angular velocity of the target spacecraft in the current state; a represents the semi-major axis of the target spacecraft; e represents the eccentricity of the target spacecraft; M0 represents the true anomaly angle of the target spacecraft at the initial moment;

[0177] Further calculation is performed to obtain the true near point angle f corresponding to the jth frame j and the eccentric anomaly angle E j , there is the following relationship;

[0178] M j =E j -esinE j

[0179]

[0180] M j Given, find E j , according to E j Find f j ;

[0181] Since it is impossible to pass M j Direct calculation of E j , so an iterative calculation method is used to repeatedly iterate the following formula:

[0182]

[0183] Until the following is satisfied

[0184] |M j +esinE j -E j |<ε

[0185] Among them, ε is the threshold, which can be set to 10 -9 ;

[0186] Further calculation:

[0187] r j =a(1-ecosE j )

[0188]

[0189]

[0190]

[0191] in, represents the geocentric position vector of the target spacecraft in the jth frame in the geocentric inertial coordinate system, r j represents the geocentric position vector modulus of the target spacecraft in the jth frame in the geocentric inertial coordinate system, ρ j represents the slant distance of the target spacecraft in the jth frame in the geocentric inertial coordinate system, represents the geocentric position vector of the observing spacecraft in the geocentric inertial coordinate system, R j represents the geocentric position vector modulus of the observing spacecraft in the geocentric inertial coordinate system, represents the unit direction vector of the target spacecraft relative to the observing spacecraft in the geocentric inertial coordinate system; β j express and The angle between , · represents the vector multiplication;

[0192] For the jth frame and the kth frame, k>j, we get (r k ,r j ) and (f k ,f j );

[0193] Therefore, the fitness function is defined as follows, where the smaller the fitness function, the better the parameters:

[0194]

[0195] in, represents the geocentric position vector of the target spacecraft in the geocentric inertial coordinate system in the kth frame, r k represents the geocentric position modulus of the target spacecraft in the kth frame in the geocentric inertial coordinate system, f k It represents the true anomaly angle corresponding to the kth frame in the geocentric inertial coordinate system, and N represents the number of timestamps (images);

[0196] The fitness function value is taken as the minimum corresponding to the semi-major axis a of the target spacecraft, the eccentricity e of the target spacecraft, and the true anomaly M of the target spacecraft at the initial moment as the prediction value 1;

[0197] 2) Second layer optimization

[0198] According to the predicted value of (a, e, M), the true near point angle f corresponding to the jth frame is obtained j and the eccentric anomaly angle E j

[0199] There are the following relationships;

[0200] M j =E j -esinE j

[0201]

[0202] M j Given, find E j , according to E j Find f j ;

[0203] Since it is impossible to pass M j Direct calculation of E j , so an iterative calculation method is used to repeatedly iterate the following formula:

[0204]

[0205] Until the following is satisfied

[0206] |M j +esinE j -E j |<ε

[0207] Among them, ε is the threshold, which can be set to 10 -9 ;

[0208] Further calculation:

[0209] k j =a(1-e 2 ) / (1+ecosf j )

[0210]

[0211]

[0212]

[0213] Among them, k j represents the position vector modulus of the target spacecraft in the jth frame in the geocentric inertial coordinate system, Ω represents the right ascension of the ascending node of the target spacecraft, ω represents the perigee argument of the target spacecraft, i represents the orbital inclination of the target spacecraft, and f jIndicates the true anomaly angle corresponding to the jth frame, represents the geocentric position vector of the observing spacecraft in the geocentric inertial coordinate system, represents the position vector of the target spacecraft in the jth frame in the geocentric inertial coordinate system, · represents the vector multiplication, θ j It represents the angular difference between the estimated value and the true value of the direction vector of the target spacecraft with respect to the observation spacecraft in the jth frame in the geocentric inertial coordinate system;

[0214] The fitness function is defined as follows, where the smaller the fitness function, the better the parameters:

[0215]

[0216] The fitness function value is taken as the minimum corresponding to the ascending node right ascension Ω, the perigee argument ω, and the orbital inclination i of the target spacecraft as the prediction value 2;

[0217] The six orbital numbers of the target spacecraft are obtained based on the predicted values ​​1 and 2.

[0218] The other steps and parameters are the same as those in Specific Embodiments 1 to 7.

[0219] The present invention may also have many other embodiments. Without departing from the spirit and essence of the present invention, those skilled in the art may make various corresponding changes and modifications based on the present invention, but these corresponding changes and modifications should all fall within the scope of protection of the claims attached to the present invention.

Claims

1. A method for short-arc orbit determination of spacecraft based on visual measurement, characterized in that: The specific process of the method is: Step 1: Obtain the orbital parameters of the observation spacecraft and relevant parameter information of the camera; When the target spacecraft appears in the field of view of the observation spacecraft, the grayscale image of the target spacecraft is continuously captured by the onboard camera of the observation spacecraft; Step 2: Determine the time stamp interval for the onboard camera on the observation spacecraft to select and capture the target spacecraft image; The position coordinates of the observed spacecraft corresponding to the timestamp in the geocentric inertial coordinate system are calculated according to the orbital parameters of the observed spacecraft and the timestamp interval of the grayscale image captured by the target spacecraft; Step 3: Obtain the pixels where the target spacecraft image captured by the observation spacecraft is located, and modify the grayscale values ​​of the pixels other than the pixels where the target spacecraft is located to 0 to obtain a processed image; Step 4: grayscale overlay the processed image obtained in step 3 according to the time stamp interval in step 2 to obtain a grayscale overlay image; Step 5: Calculate the direction vector of the target spacecraft relative to the observation spacecraft in the camera coordinate system based on the grayscale superimposed image; Step 6: Calculate the direction vector of the target spacecraft relative to the observation spacecraft in the geocentric inertial coordinate system; Step 7: Take the position vector sequence of the observation spacecraft obtained in step 2 and the direction vector sequence of the target spacecraft in the geocentric inertial coordinate system obtained in step 6 as input, and iteratively optimize the six orbital elements of the target spacecraft through the particle swarm algorithm; The orbital parameters of the observed spacecraft in step 1 are six orbital parameters, including orbital semi-major axis, eccentricity, right ascension of ascending node, orbital inclination, true anomaly, and argument of perigee; The relevant parameter information of the camera in step 1 is focal length, resolution, field of view, and camera orientation; In step 3, the pixels where the target spacecraft image captured by the observation spacecraft is located are obtained, and the gray values ​​of the pixels other than the pixels where the target spacecraft is located are modified to 0 to obtain the processed image; the specific process is: The pixels where the target spacecraft image captured by the observation spacecraft is located are obtained by using a target segmentation algorithm or Photoshop image processing software, and the gray values ​​of the pixels other than the pixels where the target spacecraft is located are modified to 0 to obtain a processed image; In step 5, the direction vector of the target spacecraft relative to the observation spacecraft in the camera coordinate system is solved based on the grayscale superimposed image; the specific process is: The initial frame target spacecraft mass center position and the last frame target spacecraft mass center position are manually marked on the grayscale overlay image, and the pixel coordinates (u i′ ,v i′ ),i′=1,2,…,N; Where N represents the number of timestamps; The projection relationship is used to obtain the direction vector corresponding to each pixel in the two-dimensional image, that is, the unit direction vector pointing from the coordinate origin to the center of mass of the target spacecraft in the camera coordinate system; The unit direction vector calculation process is as follows: According to the projection relationship of camera imaging, the following formula is obtained: Among them, w x and w y are the horizontal and vertical field of view angles, respectively; X and Y are the width and height of the camera target surface, respectively; f′ is the focal length of the camera; Further, according to the position of the pixel point, we can get Among them, (x, y) is the coordinate of the center of mass of the target spacecraft in the image physical coordinate system, (u, v) is the coordinate of the center of mass of the target spacecraft in the pixel coordinate system, and p x and p y are the number of pixels in each row and column of the image respectively; Because (w x ,w y ), (u,v) and p x 、p y Known, then combining the above formulas we can get: Then there is Finally, the vector Normalize to get the unit direction vector corresponding to a frame of image; Calculate each frame of image to obtain the unit direction vector corresponding to each frame of image; In step 6, the direction vector of the target spacecraft relative to the observation spacecraft in the geocentric inertial coordinate system is calculated; the specific process is: Define the origin of the observation spacecraft orbital coordinate system as the center of the earth, the x-axis points from the center of the earth to the spacecraft, the z-axis is the normal direction of the orbital plane, and the y-axis direction is determined according to the right-hand rule; (1) Obtain the transformation matrix from the camera coordinate system to the observation spacecraft orbit coordinate system: Calculate the transfer matrix from the camera coordinate system to the observation spacecraft orbit coordinate system, denoted as C o1 ; (2) Obtain the transformation matrix from the observation spacecraft orbital coordinate system to the geocentric inertial coordinate system: True anomaly of the observing spacecraft Argument of perigee Calculate the latitude argument of the observing spacecraft Then the transformation matrix from the observation spacecraft orbital coordinate system to the geocentric inertial coordinate system is in, is the right ascension of the ascending node of the observing spacecraft, To observe the orbital inclination of the spacecraft, is the latitude argument of the observing spacecraft; The following equations represent the rotation angle α around the x-axis, y-axis, and z-axis of the coordinate system respectively: R x (α) is the rotation angle α around the x-axis of the coordinate system, R y (α) is the rotation angle α around the y-axis of the coordinate system, R z (α) is the rotation angle α around the z-axis of the coordinate system, α is the rotation angle, R x (-90°) is the rotation angle of -90° around the x-axis of the coordinate system, R z (90°) is the rotation angle of 90° around the z-axis of the coordinate system. is the rotation angle around the z-axis of the coordinate system Rotation around the x-axis of the coordinate system Rotation around the z-axis of the coordinate system (3) Obtain the transformation matrix from the camera coordinate system to the geocentric inertial coordinate system: Assuming that the unit direction vector corresponding to each frame of the target spacecraft is v in the camera coordinate system and v' in the geocentric inertial coordinate system, we can get Where T is the transpose; In step 7, the position vector sequence of the observation spacecraft and the direction vector sequence of the target spacecraft in the geocentric inertial coordinate system are used as inputs, and the six orbital elements of the target spacecraft are iteratively optimized by the particle swarm algorithm; the specific process is: 1) First layer optimization From the target spacecraft's (a, e, M0), we get the M of the jth frame j ; M j =M0+nΔt j Where Δt j represents the time difference between the jth frame and the initial frame, current n represents the angular velocity of the target spacecraft in the current state; a represents the semi-major axis of the target spacecraft; e represents the eccentricity of the target spacecraft; M0 represents the true anomaly angle of the target spacecraft at the initial moment; Further calculation is performed to obtain the true near point angle f corresponding to the jth frame j and the eccentric anomaly angle E j , there is the following relationship; I j =Yes j -esinE j Since it is impossible to pass M j Direct calculation of E j , so an iterative calculation method is used to repeatedly iterate the following formula: Until the following is satisfied |M j +esinE j -ITS j |<ε Among them, ε is the threshold, which can be set to 10 -9 ; Further calculation: r j =a(1-ecosE j ) in, represents the geocentric position vector of the target spacecraft in the jth frame in the geocentric inertial coordinate system, r j represents the geocentric position vector modulus of the target spacecraft in the jth frame in the geocentric inertial coordinate system, ρ j represents the slant distance of the target spacecraft in the jth frame in the geocentric inertial coordinate system, represents the geocentric position vector of the observing spacecraft in the geocentric inertial coordinate system, R j represents the geocentric position vector modulus of the observing spacecraft in the geocentric inertial coordinate system, represents the unit direction vector of the target spacecraft relative to the observing spacecraft in the geocentric inertial coordinate system; β j express and The angle between , · represents the vector multiplication; For the jth frame and the kth frame, k>j, we get (r k ,r j ) and (f k ,f j ); Therefore, the fitness function is defined as follows: in, represents the geocentric position vector of the target spacecraft in the geocentric inertial coordinate system in the kth frame, r k represents the geocentric position modulus of the target spacecraft in the kth frame in the geocentric inertial coordinate system, f k Indicates the true anomaly angle corresponding to the kth frame in the geocentric inertial coordinate system, and N indicates the number of timestamps; The fitness function value is taken as the minimum corresponding to the semi-major axis a of the target spacecraft, the eccentricity e of the target spacecraft, and the true anomaly M0 of the target spacecraft at the initial moment as the prediction value 1; 2) Second layer optimization According to the predicted value of (a, e, M0), the true near point angle f corresponding to the jth frame is obtained j and the eccentric anomaly angle E j There are the following relationships; I j =Yes j -esinE j Since it is impossible to pass M j Direct calculation of E j , so an iterative calculation method is used to repeatedly iterate the following formula: Until the following is satisfied |M j +sinE j -ITS j |<ε Among them, ε is the threshold, which can be set to 10 -9 ; Further calculation: k j =a(1-e 2 ) / (1+ecosf j ) Among them, k j represents the position vector modulus of the target spacecraft in the jth frame in the geocentric inertial coordinate system, Ω represents the right ascension of the ascending node of the target spacecraft, ω represents the perigee argument of the target spacecraft, i represents the orbital inclination of the target spacecraft, and f j Indicates the true anomaly angle corresponding to the jth frame, represents the geocentric position vector of the observing spacecraft in the geocentric inertial coordinate system, represents the position vector of the target spacecraft in the jth frame in the geocentric inertial coordinate system, · represents the vector multiplication, θ j It represents the angular difference between the estimated value and the true value of the direction vector of the target spacecraft with respect to the observation spacecraft in the jth frame in the geocentric inertial coordinate system; The fitness function is defined as follows: The fitness function value is taken as the minimum corresponding to the ascending node right ascension Ω, the perigee argument ω, and the orbital inclination i of the target spacecraft as the prediction value 2; The six orbital numbers of the target spacecraft are obtained based on the predicted values ​​1 and 2.

Images