Postoperative precise attitude determination method, system and equipment for optical remote sensing satellite and storage medium

Abstract

The invention discloses an optical remote sensing satellite postmortem precision attitude determination method, system and device and a storage medium, and belongs to the technical field of aerospace. The frequencies and timestamps of the auxiliary data stored by the star sensor, the auxiliary data stored by the gyroscope and the auxiliary data stored by the flywheel are respectively aligned; respectively converting the aligned auxiliary data stored by the gyroscope and the auxiliary data stored by the flywheel to obtain a gyroscope angular velocity sequence and a flywheel angular momentum sequence; traversing auxiliary data stored by all star sensors according to a time sequence, and determining a whole-satellite attitude quaternion sequence by using star sensor optical axis double-vector attitude determination; respectively carrying out forward filtering processing on the whole satellite attitude quaternion sequence, the gyro angular velocity sequence and the disturbance torque based on the flywheel angular momentum sequence; performing reverse filtering by taking the forward filtering result as an initial value; and smoothing state variable sequences obtained by forward filtering and reverse filtering to finish precise attitude determination.

Description

Technical Field

[0001] The present invention relates to the field of aerospace technology, and particularly to a method, system, device, and storage medium for post-precision attitude determination of an optical remote sensing satellite. Background Art

[0002] With the continuous development of low-orbit satellite technology, low-orbit optical remote sensing satellites have been widely used in the fields of earth observation, map mapping, environmental monitoring, urban planning, etc. Precision attitude determination has always been a key research content in the field of sun-synchronous orbit optical remote sensing satellites. For a sun-synchronous orbit remote sensing satellite with an orbital altitude of 500 km, when there is no side-sway angle push-broom imaging of the earth, an attitude error of every 0.01° will cause the imaging center point to move 87 m on the ground. Therefore, high-precision attitude determination during push-broom imaging is very necessary, which helps to improve the ground mapping efficiency and the geometric accuracy of remote sensing images.

[0003] Limited by the computing power of the on-board computer and the communication rate of the on-board hardware system, relatively simple attitude recursion and filtering algorithms are usually used for attitude determination during the on-orbit operation of the satellite to ensure the reliability and real-time performance of attitude determination. In the prior art, Chinese Patent CN113701755A discloses "a method for attitude determination of an optical remote sensing satellite without a high-precision gyro", which obtains a linearized state equation and measurement equation by designing the kinematic and dynamic equations of the deviation quantity, and then applies the Kalman filter to obtain an estimation method for high-precision attitude quaternion and angular velocity. The simulation and test results prove the effectiveness of this method, and a low-quality, low-power, low-cost, and high-precision optical remote sensing satellite attitude determination system can be realized. However, it mainly focuses on how to achieve relatively high-precision satellite attitude determination in the case of lacking a high-precision gyro (i.e., high-precision angular velocity measurement input). For the case where the satellite is equipped with a high-precision gyro, the above algorithm cannot effectively utilize all measurement information (star sensor, gyro, flywheel), and the accuracy of indirectly estimating the angular velocity through the flywheel is not as high as that directly estimated from the gyro measurement value.

[0004] In summary, the existing methods cannot effectively utilize all measurement information and have low accuracy in estimating measurement values. Summary of the Invention

[0005] The present invention solves the problems that the existing methods cannot effectively utilize all measurement information and have low accuracy in estimating measurement values.

[0006] A method for post-precision attitude determination of an optical remote sensing satellite according to the present invention includes the following steps:

[0007] Step S1: Based on the attitude filtering period, align the frequencies and timestamps of the auxiliary data stored in all star sensors, gyroscopes, and flywheels respectively.

[0008] Step S2: Convert the aligned auxiliary data stored in the gyroscope and flywheel respectively to obtain the gyro angular velocity sequence and the flywheel angular momentum sequence.

[0009] Step S3: Traverse all the auxiliary data stored in the star sensors in chronological order, and use the star sensor double-vector attitude determination method to determine the whole-star attitude quaternion sequence.

[0010] Step S4: Based on the UKF algorithm, perform forward filtering on the interference torque for the whole-star attitude quaternion sequence, gyro angular velocity sequence, and flywheel angular momentum sequence respectively.

[0011] Step S5: Use the final value of the forward filtering as the initial value, and perform backward filtering on this initial value based on the UKF algorithm.

[0012] Step S6: Smooth the state variable sequences obtained from the forward filtering and backward filtering, and then the precise attitude determination is completed.

[0013] Further, in an embodiment of the present invention, in step S1, the frequencies of all star sensors are the same as the attitude filtering frequency, and the auxiliary data stored in all star sensors does not need to be processed.

[0014] Further, in an embodiment of the present invention, in step S1, downsample the auxiliary data stored in the gyroscope, specifically:

[0015]

[0016] where N = floor(0.5 × f 2 / f 1 ), f 1 represents the frequency of the auxiliary data stored in all star sensors, f 2 represents the frequency of the auxiliary data stored in the gyroscope, floor represents rounding down. For each star sensor timestamp t k , find t′ k in the auxiliary data stored in the gyroscope such that |t′ k - t k | ≤ 0.51 / f 2 . Search N numbers forward and backward respectively, and perform average sampling with 2N + 1 data. For the initial and end times, search N numbers backward and forward respectively for average sampling. The sampling result at time t′ k is denoted as [ω ib Ml (t′​k )), the subscript Ml of the angular velocity sequence [ω ib represents the angular velocity sequence after reducing the data frequency, and the subscript Mh represents the original high-frequency sampling sequence. Ml Specifically, in an embodiment of the present invention, in step S1, downsampling the auxiliary data stored in the flywheel is performed as follows:

[0017] Specifically, in an embodiment of the present invention, in step S1, downsampling the auxiliary data stored in the flywheel is performed as follows:

[0018]

[0019] where, v h = [v 1 v 2 v 3 v 4 , T where, v 1 , v 2 , v 3 , v 4 are all the rotational speeds of the four flywheels, N = floor(0.5 × f 3 / f 1 ), f 1 represents the frequency of all the auxiliary data stored in the star sensors, f 3 represents the frequency of the auxiliary data stored in the flywheel. For the sampling result at time t′ k , it is denoted as v l (t′ k ). The subscripts l and h of the flywheel rotational speed v represent the flywheel rotational speed sequence after reducing the data frequency and the original flywheel auxiliary data sequence, respectively.

[0020] Specifically, in an embodiment of the present invention, in step S2, the gyro angular velocity sequence and the flywheel angular momentum sequence are as follows:

[0021]

[0022] where, ω ib represents the gyro angular velocity sequence, h represents the flywheel angular momentum sequence, is the coordinate transformation matrix from the gyro measurement coordinate system to the body coordinate system, is the coordinate transformation matrix from the flywheel to the body coordinate system, J is the moment of inertia of the reaction flywheel, Ml represents the angular velocity sequence after reducing the data frequency, and the subscript l of the flywheel rotational speed v represents the flywheel rotational speed sequence after reducing the data frequency.

[0023] Specifically, in an embodiment of the present invention, in step S4, the forward filtering process is as follows:

[0024]

[0025] Among them, Z k represents the observed variable at time k, represents the one-step predicted value of the state variable from time k - 1 to time k, K k represents the Kalman filter gain matrix, represents the one-step predicted value of the observed variable from time k - 1 to time k obtained based on the sigma points and the observation equation, P k / k-1 represents the one-step predicted mean square error of the state variable from time k - 1 to time k, P (ZZ)k / k-1 represents the one-step predicted mean square error of the observed variable from time k - 1 to time k, represents the transpose of the Kalman filter gain matrix. The state variable estimation in each attitude cycle of the forward filtering process is The covariance matrix is P f,k ;

[0026] In the step S5 described above, the reverse filtering process is specifically as follows:

[0027]

[0028] Among them, represents the one-step predicted value of the state variable from time k + 1 to time k, represents the one-step predicted value of the observed variable from time k + 1 to time k obtained based on the sigma points and the observation equation, P k / k+1 represents the one-step predicted mean square error of the state variable from time k + 1 to time k, P (ZZ)k / k+1 represents the one-step predicted mean square error of the observed variable from time k + 1 to time k, represents the transpose of the Kalman filter gain matrix. The state variable estimation in each attitude cycle of the reverse filtering process is The covariance matrix is P b,k .

[0029] Furthermore, in an embodiment of the present invention, in the step S6, the smoothing process of the state variable sequences obtained by the forward filtering and the reverse filtering is specifically as follows:

[0030]

[0031]

[0032] Among them, P s,k is the covariance matrix after the smoothing process, are all the state variables after the smoothing process.

[0033] An optical remote sensing satellite post-precision attitude determination system described in the present invention includes the following modules:

[0034] Module S1 aligns the frequencies and timestamps of the auxiliary data stored in all star sensors, the auxiliary data stored in gyroscopes, and the auxiliary data stored in flywheels respectively based on the attitude filtering period.

[0035] Module S2 converts the aligned auxiliary data stored in gyroscopes and the aligned auxiliary data stored in flywheels respectively to obtain a gyro angular velocity sequence and a flywheel angular momentum sequence.

[0036] Module S3 traverses all the auxiliary data stored in star sensors in chronological order and uses the star sensor double-vector attitude determination to determine the whole-satellite attitude quaternion sequence.

[0037] Module S4 performs forward filtering on the interference torque based on the UKF algorithm for the whole-satellite attitude quaternion sequence, the gyro angular velocity sequence, and the flywheel angular momentum sequence respectively.

[0038] Module S5 uses the final value of the forward filtering as the initial value and performs backward filtering on this initial value based on the UKF algorithm.

[0039] Module S6 smooths the state variable sequences obtained from the forward filtering and the backward filtering, thus completing the precise attitude determination.

[0040] An electronic device according to the present invention includes a processor, a communication interface, a memory, and a communication bus. Among them, the processor, the communication interface, and the memory complete communication with each other through the communication bus.

[0041] The memory is used to store computer programs.

[0042] When the processor is used to execute the program stored in the memory, it implements the method steps described in any of the above methods.

[0043] A computer-readable storage medium according to the present invention stores a computer program therein, and when the computer program is executed by a processor, it implements the method steps described in any of the above methods.

[0044] The present invention solves the problem that the existing methods cannot effectively utilize all measurement information and the accuracy of the measured value estimation is relatively low. The specific beneficial effects include:

[0045] An optical remote sensing satellite post-flight precise attitude determination method according to the present invention can, based on the auxiliary data stored by star sensors, gyroscopes, and flywheels during the mission, use the UKF method on the ground to achieve high-precision estimation of the satellite attitude quaternion, angular velocity, and interference torque, thus solving the problem that the existing methods cannot effectively utilize all measurement information and the accuracy of the measured value estimation is relatively low.

[0046] An optical remote sensing satellite post - mission precise attitude determination method according to the present invention performs data fusion through the unscented Kalman filter (UKF) algorithm based on the data stored in relevant sensors and control components of the attitude control system during the satellite mission, and performs post - processing on the ground to obtain high - precision satellite attitude information. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] The above - mentioned and / or additional aspects and advantages of the present invention will become obvious and easy to understand from the following description of the embodiments in conjunction with the drawings, where:

[0048] Figure 1 is the diagram of the overall satellite attitude quaternion sequence described in Embodiment 1;

[0049] Figure 2 is the diagram of the X - axis attitude filtering result described in Embodiment 1;

[0050] Figure 3 is the diagram of the Y - axis attitude filtering result described in Embodiment 1;

[0051] Figure 4 is the diagram of the Z - axis attitude filtering result described in Embodiment 1;

[0052] Figure 5 is the diagram of the X - axis angular velocity filtering result described in Embodiment 1;

[0053] Figure 6 is the diagram of the Y - axis angular velocity filtering result described in Embodiment 1;

[0054] Figure 7 is the diagram of the Z - axis angular velocity filtering result described in Embodiment 1;

[0055] Figure 8 is the diagram of the X - axis disturbance torque estimation described in Embodiment 1;

[0056] Figure 9 is the diagram of the Y - axis disturbance torque estimation described in Embodiment 1;

[0057] Figure 10 is the diagram of the Z - axis disturbance torque estimation described in Embodiment 1. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0058] The following will clearly and completely describe various embodiments of the present invention in conjunction with the drawings. The embodiments described by referring to the drawings are exemplary and are intended to explain the present invention, and should not be construed as a limitation of the present invention.

[0059] Embodiment 1. An optical remote sensing satellite post - mission precise attitude determination method described in this embodiment includes the following steps:

[0060] Step S1: Based on the attitude filtering period, align the frequencies and timestamps of the auxiliary data stored in all star sensors, gyroscopes, and flywheels respectively.

[0061] Step S2: Convert the aligned auxiliary data stored in the gyroscope and flywheel respectively to obtain the gyro angular velocity sequence and the flywheel angular momentum sequence.

[0062] Step S3: Traverse all the auxiliary data stored in the star sensors in chronological order, and use the star sensor double-vector attitude determination to determine the whole satellite attitude quaternion sequence.

[0063] Step S4: Based on the UKF algorithm, perform forward filtering on the whole satellite attitude quaternion sequence, the gyro angular velocity sequence, and the disturbance torque based on the flywheel angular momentum sequence respectively.

[0064] Step S5: Use the final value of the forward filtering as the initial value, and perform backward filtering on this initial value based on the UKF algorithm.

[0065] Step S6: Smooth the state variable sequences obtained from the forward filtering and backward filtering, and then the precise attitude determination is completed.

[0066] In this embodiment, in step S1, the frequencies of all star sensors are consistent with the attitude filtering frequency, and the auxiliary data stored in all star sensors do not need to be processed.

[0067] In this embodiment, in step S1, downsample the auxiliary data stored in the gyroscope. Specifically:

[0068]

[0069] where N = floor(0.5 × f 2 / f 1 ), f 1 represents the frequency of the auxiliary data stored in all star sensors, f 2 represents the frequency of the auxiliary data stored in the gyroscope, floor represents rounding down. For each star sensor timestamp t k , find t′ k in the auxiliary data stored in the gyroscope such that |t′ k - t k | ≤ 0.51 / f 2 . Search N numbers forward and backward respectively, and perform average sampling with 2N + 1 data. For the initial moment and the end moment, search N numbers backward and forward respectively for average sampling. The sampling result at the t′ k moment is denoted as [ω ib Ml (t′ k ), and the angular velocity sequence [ω​ib Ml The subscript Ml of [] represents the angular velocity sequence after reducing the data frequency, and the subscript Mh represents the original high-frequency sampling sequence.

[0070] In this embodiment, in step S1, downsampling the auxiliary data stored in the flywheel is specifically as follows:

[0071]

[0072] where v h =[v 1 v 2 v 3 v 4 T , v 1 , v 2 , v 3 , v 4 are all the rotational speeds of the four flywheels. N = floor(0.5 × f 3 / f 1 ), f 1 represents the frequency of all the auxiliary data stored in the star sensors, f 3 represents the frequency of the auxiliary data stored in the flywheel. For the sampling result at time t′ k , it is denoted as v l (t′ k ). The subscripts l and h of the flywheel rotational speed v represent the flywheel rotational speed sequence after reducing the data frequency and the original flywheel auxiliary data sequence respectively.

[0073] In this embodiment, in step S2, the gyro angular velocity sequence and the flywheel angular momentum sequence are specifically as follows:

[0074]

[0075] where ω ib represents the gyro angular velocity sequence, h represents the flywheel angular momentum sequence, is the coordinate transformation matrix from the gyro measurement coordinate system to the body coordinate system, is the coordinate transformation matrix from the flywheel to the body coordinate system, J is the moment of inertia of the reaction flywheel, Ml represents the angular velocity sequence after reducing the data frequency, and the subscript l of the flywheel rotational speed v represents the flywheel rotational speed sequence after reducing the data frequency.

[0076] In this embodiment, in step S4, the forward filtering process is specifically as follows:

[0077]

[0078] where Z k represents the observation variable at time k,​​ Denotes the one-step predicted value of the state variable from time k-1 to time k, K k Denotes the Kalman filter gain matrix, Denotes the one-step predicted value of the observation variable from time k-1 to time k obtained based on the σ points and the observation equation, P k / k-1 Denotes the one-step predicted mean square error of the state variable from time k-1 to time k, P (ZZ)k / k-1 Denotes the one-step predicted mean square error of the observation variable from time k-1 to time k, Denotes the transpose of the Kalman filter gain matrix. The state variable estimation in each attitude cycle of the forward filtering process is The covariance matrix is P f,k ;

[0079] In the step S5, the reverse filtering process is specifically as follows:

[0080]

[0081] Wherein, Denotes the one-step predicted value of the state variable from time k+1 to time k, Denotes the one-step predicted value of the observation variable from time k+1 to time k obtained based on the σ points and the observation equation, P k / k+1 Denotes the one-step predicted mean square error of the state variable from time k+1 to time k, P (ZZ)k / k+1 Denotes the one-step predicted mean square error of the observation variable from time k+1 to time k, Denotes the transpose of the Kalman filter gain matrix. The state variable estimation in each attitude cycle of the reverse filtering process is The covariance matrix is P b,k 。

[0082] In this embodiment, in the step S6, the smoothing process of the state variable sequences obtained by the forward filtering and the reverse filtering is specifically as follows:

[0083]

[0084]

[0085] Wherein, P s,k Is the covariance matrix after the smoothing process, Is the state variable after the smoothing process.

[0086] Existing methods mainly target the situation where the satellite carries a high-precision gyroscope, and cannot effectively utilize all measurement information (star sensors, gyroscopes, flywheels), and the accuracy of the obtained measurement estimates is relatively low.

[0087] To solve the above technical problems, that is, to further improve the attitude determination accuracy of the satellite, a feasible concept is proposed in this embodiment: using the sensors carried by the satellite to collect and store data in orbit, and then downloading the data during the mission to the ground for post-processing. Since the ground no longer has the above hardware conditions and real-time limitations, more complex and higher-precision attitude determination algorithms can be selected for high-precision post-attitude determination during post-processing.

[0088] Therefore, based on the unscented Kalman filter, this embodiment proposes a method for post-precision attitude determination estimation of an optical remote sensing satellite. In this method, the attitude quaternion, angular velocity, and disturbance torque are used as state variables to build a UKF filter. The method includes the following steps:

[0089] Step S1, set the attitude filtering period T, and align the auxiliary data alignment frequency and timestamp of different sensors.

[0090] Generally, the frequencies at which different sensors on the satellite can support data collection are different. Denote the frequencies of the star sensor collecting the attitude quaternion, the gyroscope collecting the angular velocity, and the reaction wheel collecting the rotational speed as f 1 、f 2 and f 3 . The attitude filtering frequency can be selected to be consistent with the sensor with the lowest sampling frequency. Generally, among the above sensors of the remote sensing satellite, the star sensor has the lowest sampling frequency. Therefore, the attitude filtering period can be set as T = 1 / f 1 .

[0091] The auxiliary data stored by the star sensor is the attitude quaternion for attitude determination and the corresponding UTC time. Since the attitude filtering frequency is consistent with the star sensor frequency, the auxiliary data stored by the star sensor does not require additional processing.

[0092] The auxiliary data stored by the gyroscope is the three-axis angular velocity [ω ib Mh in the gyro measurement coordinate system and the corresponding UTC time. To make the frequency of the angular velocity data consistent with the frequency of the attitude quaternion solved by the star sensor, downsampling is performed on the gyro auxiliary data, and average sampling is performed according to their frequencies. The method is as follows:

[0093]

[0094] where N = floor(0.5 × f 2 / f 1 ), and floor means rounding down. For each star sensor timestamp t k , find t k in the gyro auxiliary data such that |t′ k - t k | ≤ 0.51 / f 2 ​, Search forward and backward for N numbers respectively according to the method of formula (1), and perform average sampling with 2N + 1 data. For the initial moment and the end moment, search backward and forward for N numbers respectively for average sampling. For the sampling result at time t′ k is denoted as [ω ib Ml (t′ k ).

[0095] The auxiliary data stored in the flywheel is the flywheel speed v h and the corresponding UTC time. For the four obliquely mounted flywheels, v h = [v 1 v 2 v 3 v 4 T , where v 1 , v 2 , v 3 , v 4 are the speeds of the four flywheels; the downsampling method for the flywheel auxiliary data is similar to that of the gyroscope:

[0096]

[0097] where N = floor(0.5 × f 3 / f 1 ), and the sampling result at time t′ k is denoted as v l (t′ k ).

[0098] Step S2, Calculate the body angular velocity sequence and the flywheel angular momentum sequence according to the downsampled gyroscope and flywheel data.

[0099] The downsampled gyroscope data and flywheel data still cannot be directly used and need to be transformed to obtain the body angular velocity sequence ω ib and h required for attitude dynamics. Calculate according to the following method:

[0100]

[0101] where is the coordinate transformation matrix from the fiber optic gyro measurement coordinate system to the body coordinate system, is the coordinate transformation matrix from the four obliquely mounted flywheels to the body coordinate system (this matrix is a 3×4 matrix), and J is the moment of inertia of the reaction flywheel.

[0102] Step S3, Traverse all star sensor attitude quaternions in chronological order, and use the star sensor optical axis double vector attitude determination to determine the whole - star attitude quaternion sequence Q ib .

[0103] ​​Select the two star sensors with the highest and second highest priority as the basis for attitude determination. For the moments when the attitude quaternion cannot be solved, recursion is performed based on the attitude quaternion of the adjacent successful solution moments. For example, when star sensor A fails to successfully solve the quaternion at time k, and quaternion Q has been solved at time k-1 A (k-1), then use the angular velocity information ω obtained in step S2 ib (k-1), estimate the attitude quaternion at time k:

[0104]

[0105] in, symbol represents quaternion multiplication operation, Represents the coordinate transformation matrix from the star body coordinate system to the star sensor A measurement coordinate system.

[0106] Through the above method, the attitude quaternions generated by the two star sensors with the highest priority can be made continuous and valid. Here, let’s assume that the corresponding attitude quaternion sequences are Q A and Q B , and then we can traverse from front to back in time, and use the double vector attitude determination based on the optical axis of star sensor A and star sensor B at each moment to uniquely determine the quaternion sequence Q of the whole star attitude ib Since the time of star sensitivity A and star sensitivity B is not necessarily strictly consistent, the angular velocity sequence ω is also introduced here ib In the process of dual-vector attitude determination, the attitude quaternion is compensated, and the compensation principle is the same as that of formula (5).

[0107] Step 4: Forward UKF filtering.

[0108] Through the above three steps, the auxiliary data of the star sensor, gyroscope, and flywheel have been converted into attitude quaternion sequences Q with the same length and aligned timestamps. ib , angular velocity sequence ω ib And the flywheel angular momentum sequence h, assuming the sequence length is N.

[0109] Taking the quaternion of the star's attitude, angular velocity, and disturbance torque as state variables,

[0110] X=[Q ib ω ib T d ] T ; (6)

[0111] The dynamic equation of the system can be expressed as:

[0112]

[0113] where the subscript k represents the current k-th moment, and u k is the control variable input at the current moment, and w k is the process noise vector. For the differential equation of the state variable, the following expansions are carried out for the three parts respectively:

[0114]

[0115] where J s is the whole-satellite inertia matrix. Combining equations (8)-(10), the differential equation corresponding to equation (7) can be determined by the state variable X k at the k-th moment, the flywheel angular momentum h, and the control torque T c . The sequence of flywheel angular momentum has been obtained in equation (4). Since the control torque of the satellite is actually provided by the flywheel, T c (k) = -(h k+1 -h k ) / T. Therefore, the control u k at the k-th moment can be regarded as a combination of the flywheel angular momentum h and the control torque T c . The observables of the system are selected as the attitude quaternion and the angular velocity:

[0116] Z k = h(X k ) + V k ; (11)

[0117] That is, h(X k ) = [Q ib ω ib T , and V k is the measurement noise vector.

[0118] Then, the UKF algorithm can be carried out in the following order:

[0119] Step S401, initialization:

[0120] X 0 is the value of the state variable X at the initial moment.

[0121] Step S402, generate σ points:

[0122]

[0123] where is the i-th column of the lower triangular Cholesky decomposition of the matrix (n + λ)P k-1 , n is the dimension of the vector X, and in this embodiment, n = 10. λ = 3α 2 - n, and α is a very small positive number satisfying 10 -4 ≤ α ≤ 1.​

[0124] Step S403, determine the weight:

[0125]

[0126] For a normal distribution, β = 2.

[0127] Step S404, perform a one-step prediction:

[0128]

[0129] where Q k-1 is the process noise matrix at time k-1.

[0130] Step S405, calculate the one-step prediction sample point at time k:

[0131]

[0132] Step S406, calculate P (XZ)k / k-1 and P (ZZ)k / k-1 :

[0133]

[0134] where R k is the measurement noise variance matrix.

[0135] Step S407, calculate the filter gain matrix:

[0136]

[0137] Step S408, calculate the filtering result:

[0138]

[0139] where Z k is the observed variable at time k.

[0140] Through the above steps, filtering processing for the above attitude quaternion and angular velocity sequence can be achieved based on UKF, and the disturbance torque can be estimated. Record the state variable estimation and covariance matrix P k in each attitude cycle of the forward filtering process, and denote them as and P f,k .

[0141] Step S5, inverse UKF filtering.

[0142] In order to make more effective use of the auxiliary data, a backward filtering link is cascaded after the forward filtering link is completed. The initial value of the backward filter is the final value of the forward filter. Mark the backward filtering link with subscript b, then:

[0143]

[0144] Subsequently, recursive filtering is performed from back to front in the opposite direction. At this time, the one-step prediction in Equation (14) becomes:

[0145]

[0146] The remaining calculation methods are the same as those of the forward filter. Record the state variable estimates and the covariance matrix P k in each attitude cycle of the backward filtering process, and record them as and P b,k .

[0147] Step S6, perform optimal smoothing on the forward and backward filtering results.

[0148] The errors in the front and back parts of the measurement sequence are uncorrelated. Therefore, the optimal smoothing estimation result at time k can be obtained by using the information fusion formula as follows:

[0149]

[0150] Therefore, in this embodiment, for the problem of post-precision attitude determination of on-orbit data of an optical remote sensing satellite, a high-precision post-attitude determination method based on the unscented Kalman filter algorithm and combined with "forward-backward filtering" for optimal smoothing is proposed.

[0151] To better illustrate the post-precision attitude determination method of an optical remote sensing satellite described in this embodiment, it is described in detail through the following embodiments:

[0152] Numerical simulation was carried out based on the UKF algorithm using the auxiliary data accumulated during the on-orbit imaging mission of a certain type of optical remote sensing satellite. The orbital altitude of this satellite is 520 km, the orbital inclination is 97.54°, the attitude update frequency is consistent with the sampling frequency of the star sensor auxiliary data, and the sampling frequencies of the gyroscope and flywheel are higher than those of the star sensor. The specific implementation process is as follows.

[0153] Step S1, preprocessing of sensor data:

[0154] Preprocess the auxiliary data accumulated during the on-orbit imaging mission of this type of optical remote sensing satellite. Use the frequency of the auxiliary data stored in the two star sensors carried on it as the attitude filtering frequency, downsample the gyroscope and flywheel auxiliary data, and generate an angular velocity sequence and a momentum moment sequence.

[0155] Step S2, generate the whole-satellite attitude quaternion sequence:

[0156] Traverse the auxiliary data of two star sensors, extract the attitude determination quaternion, and then generate the attitude quaternion sequence according to the star sensor optical axis double-vector attitude determination. Perform quaternion compensation calculation on the star maps with failed solutions based on adjacent valid star maps and angular velocity. Finally, obtain the whole-satellite quaternion attitude sequence during the mission as Figure 1 shown.

[0157] Step S3, forward-backward filtering and optimal smoothing:

[0158] The parameters of the embodiment are shown in Table 1. In the initial value setting of the forward filter, the attitude quaternion and angular velocity directly use the sequence values at the initial moment, and the disturbance torque at the initial moment is set to 0. After the filter is started, it can converge quickly. For the convenience of comparison, the deviation of the attitude quaternion is characterized by the Euler angle deviation on the X, Y, and Z axes, as Figures 2 - 4 shown; the estimation results of the three-axis angular velocity and the estimation results of the three-axis disturbance torque are respectively as Figures 5 - 7 shown, Figures 8 - 10 shown.

[0159] Table 1

[0160]

[0161] In summary, the results of this embodiment show that the method provided by this embodiment can achieve high-precision estimation of the attitude quaternion and angular velocity, and can also estimate the disturbance torque.

[0162] Embodiment 2. An optical remote sensing satellite post-precision attitude determination system described in this embodiment includes the following modules:

[0163] Module S1, based on the attitude filtering period, align the frequencies and timestamps of the auxiliary data stored in all star sensors, the auxiliary data stored in the gyroscope, and the auxiliary data stored in the flywheel respectively;

[0164] Module S2, respectively convert the aligned auxiliary data stored in the gyroscope and the auxiliary data stored in the flywheel to obtain the gyro angular velocity sequence and the flywheel angular momentum sequence respectively;

[0165] Module S3, traverse all the auxiliary data stored in the star sensors in chronological order, and use the star sensor optical axis double-vector attitude determination to determine the whole-satellite attitude quaternion sequence;

[0166] Module S4, based on the UKF algorithm, perform forward filtering on the whole-satellite attitude quaternion sequence, the gyro angular velocity sequence, and the disturbance torque based on the flywheel angular momentum sequence respectively;

[0167] Module S5, use the final value of the forward filtering as the initial value, and perform backward filtering on the initial value based on the UKF algorithm;

[0168] Module S6 smooths the final value of the forward filtering and the final value of the backward filtering, thus completing precise attitude determination.

[0169] Embodiment 3. An electronic device according to this embodiment includes a processor, a communication interface, a memory, and a communication bus. Among them, the processor, the communication interface, and the memory complete communication with each other through the communication bus;

[0170] The memory is used to store a computer program;

[0171] The processor is used to implement the method steps described in Embodiment 1 when executing the program stored on the memory.

[0172] Embodiment 4. A computer-readable storage medium according to this embodiment stores a computer program in the computer-readable storage medium. When the computer program is executed by a processor, the method steps described in Embodiment 1 are implemented.

[0173] The above has introduced in detail a method, a system, a device, and a storage medium for post-event precise attitude determination of an optical remote sensing satellite proposed by the present invention. Specific examples are used in this article to elaborate on the principle and implementation manner of the present invention. The description of the above embodiments is only used to help understand the method and its core idea of the present invention; at the same time, for those of ordinary skill in the art, according to the idea of the present invention, there will be changes in the specific implementation manner and application scope. In summary, the content of this specification should not be construed as a limitation to the present invention.

Claims

1. A method for post-precision attitude determination of an optical remote sensing satellite, characterized in that: The following steps are involved: Step S1, based on the attitude filtering period, align the frequencies and timestamps of all auxiliary data stored by the star sensors, the auxiliary data stored by the gyroscopes, and the auxiliary data stored by the flywheels; Step S2, converting the auxiliary data stored in the aligned gyro and the auxiliary data stored in the flywheel to obtain a gyro angular velocity sequence and a flywheel angular momentum sequence respectively; Step S3, traversing all auxiliary data stored by the star sensors in chronological order, and using the star sensor optical axis double vector attitude determination to determine the quaternion sequence of the entire star attitude; Step S4, based on the UKF algorithm, forward filtering is performed on the whole satellite attitude quaternion sequence, the gyro angular velocity sequence, and the interference torque based on the flywheel angular momentum sequence; Step S5, using the final value of the forward filtering as the initial value, and performing reverse filtering on the initial value based on the UKF algorithm; Step S6, smoothing the state variable sequence obtained by forward filtering and reverse filtering, thereby completing precise attitude determination.

2. The method for post-precision attitude determination of an optical remote sensing satellite according to claim 1, characterized in that: In the step S1, all star sensor frequencies are consistent with the attitude filtering frequency, and the auxiliary data stored in all star sensors do not need to be processed.

3. The method for post-precision attitude determination of an optical remote sensing satellite according to claim 1, characterized in that: In the step S1, downsampling the auxiliary data stored in the gyroscope is performed, specifically: Where N = floor (0.5 × f2 / f1), f1 represents the frequency of auxiliary data stored by all star sensors, f2 represents the frequency of auxiliary data stored by the gyroscope, and floor represents rounding down. For each star sensor timestamp t k , find t′ in the auxiliary data stored by the gyro k So that |t′ k -t k |≤0.51 / f2, search N numbers forward and backward respectively, and take 2N+1 data for average sampling. For the initial time and the end time, search N numbers backward and forward respectively for average sampling. For t′ k The sampling result at the time is recorded as [ω ib ] Ml (t′ k ), angular velocity sequence [ω ib ] Ml The subscript M1 represents the angular velocity sequence after the data frequency is reduced, and the subscript Mh represents the original high-frequency sampling sequence.

4. The method for post-precision attitude determination of an optical remote sensing satellite according to claim 1, characterized in that: In the step S1, downsampling the auxiliary data stored in the flywheel is performed as follows: Among them, v h =[v1 v2 v3 v4] T , v1, v2, v3, v4 are the rotation speeds of the four flywheels, N = floor (0.5 × f3 / f1), f1 represents the frequency of the auxiliary data stored by all star sensors, f3 represents the frequency of the auxiliary data stored by the flywheel, for t′ k The sampling result at time v l (t′ k ), the subscripts l and h of the flywheel speed v represent the flywheel speed sequence after the data frequency is reduced and the original flywheel auxiliary data sequence, respectively.

5. The method for post-precision attitude determination of an optical remote sensing satellite according to claim 1, characterized in that: In step S2, the gyro angular velocity sequence and the flywheel angular momentum sequence are specifically: Among them, ω ib represents the gyro angular velocity sequence, h represents the flywheel angular momentum sequence, is the coordinate transformation matrix from the gyro measurement coordinate system to the stellar coordinate system, is the coordinate transformation matrix from the flywheel to the planetary coordinate system, J is the moment of inertia of the reaction flywheel, Ml represents the angular velocity sequence after the data frequency is reduced, and the subscript l of the flywheel speed v represents the flywheel speed sequence after the data frequency is reduced.

6. The method for post-precision attitude determination of an optical remote sensing satellite according to claim 1, characterized in that: In the step S4, the forward filtering process is specifically as follows: Among them, Z k represents the observed variable at time k, represents the one-step prediction value of the state variable from time k-1 to time k, K k represents the Kalman filter gain matrix, represents the one-step prediction value of the observed variable from time k-1 to time k based on the σ point and the observation equation, P k / k-1 represents the one-step prediction mean square error of the state variable from time k-1 to time k, P (ZZ)k / k-1 It represents the one-step prediction mean square error of the observed variable from time k-1 to time k. represents the transpose of the Kalman filter gain matrix. The state variable in each attitude cycle of the forward filtering process is estimated as The covariance matrix is ​​P f,k ; In step S5, the reverse filtering process is specifically as follows: in, represents the one-step predicted value of the state variable from time k+1 to time k, represents the one-step prediction value of the observed variable from time k+1 to time k based on the σ point and the observation equation, P k / k+1 represents the one-step prediction mean square error of the state variable from time k+1 to time k, P (ZZ)k / k+1 It represents the one-step prediction mean square error of the observed variable from time k+1 to time k. represents the transpose of the Kalman filter gain matrix. The state variable in each attitude cycle of the reverse filtering process is estimated as The covariance matrix is ​​P b,k .

7. The method for post-precision attitude determination of an optical remote sensing satellite according to claim 1, characterized in that: In the step S6, the state variable sequence obtained by forward filtering and reverse filtering is smoothed, specifically: Among them, P s,k is the covariance matrix after smoothing, All are state variables after smoothing.

8. An optical remote sensing satellite post-precision attitude determination system, characterized in that: Includes the following modules: Module S1, based on the attitude filtering cycle, aligns the frequencies and timestamps of all auxiliary data stored by the star sensor, the auxiliary data stored by the gyro, and the auxiliary data stored by the flywheel; Module S2, converting the auxiliary data stored in the aligned gyro and the auxiliary data stored in the flywheel to obtain the gyro angular velocity sequence and the flywheel angular momentum sequence respectively; Module S3, traverses all auxiliary data stored by star sensors in chronological order, and uses the star sensor optical axis double vector attitude determination to determine the quaternion sequence of the entire star attitude; Module S4, based on the UKF algorithm, performs forward filtering on the whole satellite attitude quaternion sequence, the gyro angular velocity sequence, and the interference torque based on the flywheel angular momentum sequence; Module S5, using the final value of the forward filtering as the initial value, and performing reverse filtering on the initial value based on the UKF algorithm; Module S6 performs smoothing on the state variable sequence obtained by forward filtering and reverse filtering to complete precise attitude determination.

9. An electronic device, characterized in that: It includes a processor, a communication interface, a memory and a communication bus, wherein the processor, the communication interface and the memory communicate with each other through the communication bus; Memory, used to store computer programs; A processor, for implementing the method steps described in any one of claims 1 to 7 when executing a program stored in a memory.

10. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the method steps described in any one of claims 1 to 7 are implemented.

Images