Gyroscope error parameter rapid calibration method based on intelligent optimization maneuvering

Abstract

The invention discloses a gyroscope error parameter rapid calibration method based on intelligent optimization maneuvering, and belongs to the field of spacecraft attitude determination. The method specifically comprises the following steps: formulating a quantum genetic coding mode for an attitude maneuver sequence of a satellite attitude determination system, and coding and decoding a current chromosome; secondly, establishing a gyro error parameter calibration model based on star sensor measurement assistance and a satellite angular velocity model, so as to establish a fitness function of quantum genetic optimization; selecting the current-generation optimal chromosome according to the fitness function, judging whether the fitness difference between the current-generation optimal chromosome and the parent-generation optimal chromosome is smaller than a threshold value or not or whether the fitness difference reaches the number of iterations or not, and if yes, extracting the gene sequence of the chromosome as a final optimization result; otherwise, updating the chromosome by using the quantum revolving door, and returning to execute the decoding of the next-generation chromosome. And finally, executing the optimized attitude maneuver sequence by the satellite, and estimating a final calibration value of the gyroscope error parameter. According to the invention, the precision of gyroscope error parameter calibration is improved.

Description

Technical Field

[0001] The present invention relates to the field of spacecraft attitude determination, and specifically to a method for rapidly calibrating gyro error parameters based on intelligent optimization maneuvers. Background Art

[0002] High-precision attitude information is the key for a satellite to complete complex observation and flight missions. The main factor for a satellite to achieve high-precision attitude determination depends on the accuracy level of attitude sensors.

[0003] Currently, the commonly used satellite attitude determination system mainly consists of on-board gyroscopes and star sensors. Among them, on-board gyroscopes can output the angular velocity information of the satellite during on-orbit operation in real time, and then obtain the satellite attitude through calculation, with good autonomous reliability and a high data update rate. Star sensors obtain high-precision attitude information through calculation by measuring the starlight vectors of stars and combining reference vectors obtained from prior information such as star catalogs. Due to the influence of complex space environments during the on-orbit operation of the satellite and the deficiencies of the sensors themselves, relatively large measurement errors are generated. In order to achieve high-precision attitude determination of the satellite, it is necessary to effectively calibrate and compensate for the measurement errors of attitude sensors.

[0004] The gyro measurement error terms mainly include constant drift, installation error, scale factor error, and random measurement noise, etc. By calibrating and compensating the key error parameters of the gyro, the attitude accuracy of its calculation can be effectively improved. Although the star sensor has a relatively low data update rate and is vulnerable to dynamic environments, its accuracy is relatively high under long-term operating conditions and the measurement error does not accumulate over time. Therefore, using the measurement information of the star sensor to assist in calibrating and compensating the main error parameters of the gyro is an effective means to improve the satellite attitude determination accuracy.

[0005] For a satellite with a basically constant or slowly changing angular velocity in the normal operating state, the error parameters of the on-board gyroscope are difficult to be effectively estimated due to insufficient observability. If the satellite performs an attitude maneuver with a variable angular velocity, at this time, error terms such as the installation and scale factor of the gyro will have different forms of influence on the measurement output, thereby improving the observability of the gyro error parameters. When improving the observability of gyro error parameters through satellite attitude maneuvers, not only the system's requirement for the rapidity of the calibration process should be met, but also various constraint conditions such as maneuver amplitude limitations, easy implementation of maneuver forms, and effective measurement by star sensors should be considered.

[0006] At present, a maneuvering method in which satellites sequentially perform several rounds of fixed maneuvers around each axis of this system is usually adopted. However, this method does not consider the excitation effect of the satellite's own running angular velocity on the error parameters, but treats each axis equally, which causes a certain degree of waste of on-board resources and burdens the attitude control system. Therefore, how to design a satellite attitude maneuvering scheme with low cost and easy implementation to effectively stimulate the observability of various gyro error parameters and then improve the calibration effect is a problem that needs to be solved.

[0007] Bio-inspired optimization algorithms represented by genetic algorithms (GA) are effective means to solve the optimal design problem of maneuvering sequences. The genetic algorithm simulates the survival-of-the-fittest rule in the biological evolution process and searches for the optimal individual through three basic operations: chromosome selection, crossover, and mutation. This algorithm uses the objective function to perform global search under the guidance of probability, is not restricted by factors such as the nature of the problem and the form of the optimization criterion, and has good robustness and wide applicability. However, the conventional genetic algorithm is very sensitive to the ways of chromosome selection, crossover, and mutation, and is extremely prone to problems such as a large number of iterative evolution times and slow convergence speed. In addition, this heuristic optimization algorithm is also prone to falling into local optimal values.

[0008] Aiming at the problems existing in the conventional GA, the quantum genetic algorithm (QGA) combines the rapidity of quantum computing and the optimization ability of the genetic algorithm, adopts quantum bit encoding, and uses quantum gate updates to achieve search evolution. It has the characteristics of a small population size, fast convergence speed, and strong global search ability, making it more advantageous in dealing with complex optimization problems and has been widely used in planning problems in fields such as robots and unmanned aerial vehicles.

[0009] In summary, using the fast-converging iterative evolution algorithm QGA to optimize the design of the satellite attitude maneuvering sequence can, while fully stimulating the observability of the gyro error parameters, use limited star sensor measurement information for assistance to achieve high-precision on-orbit rapid calibration of the on-board gyro error parameters. Summary of the Invention

[0010] Aiming at the on-orbit rapid high-precision calibration requirements of on-board gyro error parameters, the present invention provides a fast calibration method for gyro error parameters based on intelligent optimization maneuvers. By substituting the measured data of the attitude sensor collected into the gyro error parameter calibration model, high-precision error parameter solution is realized, and the gyro measurement error is compensated to improve the on-orbit measurement accuracy of the angular velocity information.

[0011] The specific steps of the fast calibration method for gyro error parameters based on intelligent optimization maneuvers are as follows:

[0012] Step 1: For the satellite's attitude determination system, determine the basic maneuver forms and parameters that make up the attitude maneuver sequence.

[0013] The basic maneuver forms include six ways of changing the angular rate magnitude in the form of triangle, trapezoid, sine, and their combinations, which are used as the basic building blocks of the attitude maneuver sequence.

[0014] The parameters reflecting the attitude maneuver sequence include maneuver form, maximum maneuver angular rate amplitude, single maneuver period, single-direction maneuver duration, total number of maneuver rounds, and total maneuver time, etc.

[0015] Step 2: Develop a quantum genetic coding method to encode and decode the contemporary chromosomes of the quantum genetic algorithm for subsequent iterative optimization of the attitude maneuver sequence.

[0016] The present invention uses the decimal integer method for encoding. Each chromosome represents the encoding of an attitude maneuver sequence. The attitude maneuver parameters participating in the encoding and optimization include maneuver type, the axis where the current maneuver is implemented, maneuver angular rate amplitude, etc.

[0017] The number of genes in each chromosome is determined by the number of maneuvers that make up the attitude maneuver sequence. For a chromosome with a gene length of m, the encoding method of its k-th gene "abcd" is as follows

[0018] The thousands digit a is the encoding of the basic maneuver form, 1 ≤ a ≤ n m , n m is the number of candidate maneuver forms.

[0019] The hundreds digit b is the encoding of the maximum maneuver angular rate amplitude η max , and the specific calculation formula is as follows:

[0020]

[0021] where η max and are respectively the lower and upper limits of the value range of the maximum maneuver angular rate amplitude η max .

[0022] In the decoding stage, the maximum maneuver angular rate amplitude η max is decoded according to the following formula:

[0023]

[0024] The tens digit c represents the maneuver axis: 1 corresponds to the x-axis, i.e., the pitch axis, 2 corresponds to the y-axis, i.e., the yaw axis, 3 corresponds to the z-axis, i.e., the roll axis, 1 ≤ c ≤ 3.

[0025] The units digit d is the flag encoding for whether a maneuver is executed in the current period. When d is odd, it indicates an attitude maneuver, and the maneuver flag is set to 1; when d is even, it indicates no attitude maneuver, and the maneuver flag is set to 0.

[0026] Convert the above attitude maneuver sequence in decimal integer encoding form into binary encoding, then introduce the state vector expression of the quantum, and use quantum bits to encode the genes.

[0027] Step 3: Establish a calibration model for gyro error parameters assisted by star sensor measurements, and use the least squares method to estimate the gyro error parameters.

[0028] The gyro error parameters include the constant drift b, installation error Δ, scale factor error δs, and random measurement noise w g ;

[0029] The specific steps are as follows:

[0030] Step 301: Establish the measurement model of the three-axis gyro assembly as follows:

[0031]

[0032] where ω g is the angular velocity measured and output by the gyro assembly; ω b is the true angular velocity of the satellite rotating relative to the inertial system in this coordinate system; is the installation matrix of the gyro assembly; b is the constant drift vector; Λ is a diagonal matrix with the components of the scale factor error δs as the diagonal elements; is the installation error matrix of the gyro assembly; w g is the random measurement noise vector; I is the identity matrix.

[0033] Step 302: Establish the measurement model of the star sensor as follows:

[0034]

[0035] where Q st is the attitude quaternion measured and output by the star sensor; is the true attitude quaternion of the satellite, that is, the attitude representation of the satellite's body coordinate system relative to the inertial system; is the quaternion corresponding to the installation matrix of the star sensor; Q ξ is the random measurement noise ξ represented in quaternion form.

[0036] Step 303: Use the measurement models of the three-axis gyro assembly and the star sensor to construct a calibration model for gyro error parameters;

[0037] First, use the measurement information (Q st ) obtained by the star sensor at times t and t + nΔtt and (Q st ) t+nΔt ; Calculate the attitude change of the star sensor during [t, t + nΔt).

[0038] where

[0039] Then, according to the angular velocity sequence {ω b} measured by the three-axis gyroscope during [t, t + nΔt), use the Picard integral method to recursively calculate the attitude quaternion change G({ω g )Q 0 , where Q 0 is the attitude quaternion at the initial moment; G({ω g ) represents the continuous multiplication of the angular velocity exponential matrices.

[0040] Finally, use the attitude change of the star sensor and the attitude quaternion change of the gyroscope measurement error to construct the gyroscope error parameter calibration model as follows:

[0041]

[0042] where Δt is the time interval of measurement sampling; Δ is the gyroscope installation error; δs is the gyroscope scale factor error; In the formula, n is the total number of samplings, and are the coefficients corresponding to the Picard solution method; is the angular velocity measured by the gyroscope at the moment t + jΔt; represents the formed skew-symmetric matrix; is the diagonal matrix composed of the elements of ; V t represents the model equivalent noise caused by gyroscope random measurement noise, star sensor random measurement noise, model linearization error, etc.

[0043] Step 304. Simplify the gyroscope error parameter calibration model to obtain:

[0044]

[0045] where X = [b T Δ T δs T T is the gyroscope error parameter vector to be calibrated; is the vector part on the left side of the gyroscope error parameter calibration model; is the vector part of the equivalent noise term in the gyroscope error parameter calibration model on the right side,​ is the measurement matrix.

[0046] Step 305: Collect the measurement data of the star sensor and the gyroscope, and establish an error calibration equation for the entire attitude maneuver period interval;

[0047] The star sensor has made (N + 1) effective measurements at times t k (k = 1, 2, …, N + 1); during two adjacent star sensor measurements [t k , t k+1 , the gyroscope measures n k times, and at this time t k+1 = t k + n k Δt.

[0048] Based on this, multiple calibration equations are established and combined to obtain the error calibration equation within the entire maneuver period as follows:

[0049] Z = HX + v (7) where,

[0050] Step 306: Use the least - square criterion combined with the error calibration equation to find the parameter values that minimize the residual of the gyro error parameter calibration model, and then the estimation of the gyro error parameters to be calibrated can be obtained;

[0051] The formula is as follows:

[0052]

[0053] where the estimated value From this, the estimated values of various error parameters can be obtained, including the constant bias estimation the installation error estimation and the scale factor error estimation

[0054] Step Four: Establish a satellite angular velocity model in the presence of attitude maneuvers;

[0055] The specific formula is as follows:

[0056]

[0057] where, is the projection of the angular velocity of the orbital system relative to the inertial Earth system in the inertial Earth system, that is, the angular velocity during the normal operation of the satellite; is the transformation matrix from the inertial Earth system to the orbital coordinate system; is the transformation matrix from the orbital coordinate system to the satellite body coordinate system; is the projection of the rotational angular velocity of the satellite's body system relative to the orbital system in the body system, representing the angular velocity when the satellite executes an additional attitude maneuver sequence around the body axis in addition to the attitude angle motion during normal satellite operation; is the projection of the rotational angular velocity of the satellite's body system relative to the inertial Earth system in the satellite's body system.

[0058] Step Five: Establish a fitness function optimized by quantum genetics using the gyro error parameter calibration model and the satellite angular velocity model;

[0059] The fitness function is as follows:

[0060]

[0061] where trace(·) represents the trace operation of the matrix; t m is the total time consumed during the satellite maneuver; θ m is the total rotation angle during the satellite maneuver; η m is the maximum maneuver angular rate amplitude during the satellite maneuver; is the time when the satellite operates at the maximum maneuver angular rate amplitude; w o 、w t 、w θ and w η are the weight coefficients reflecting the contribution degrees of the observable level, maneuver time consumption, maneuver rotation angle, and maneuver angular rate amplitude in the fitness function.

[0062] Step Six: Select the chromosome with the highest fitness function value from the current generation of chromosomes, and determine whether the fitness difference between the current best chromosome and the best chromosome of the parent generation is less than a preset threshold, or whether the number of iterative evolution times has reached the preset maximum value. If so, end the evolution process offline and transfer to Step Eight; otherwise, execute Step Seven to achieve the global search for the optimal maneuver sequence.

[0063] Before iterative optimization, the chromosome population is artificially initialized, and the fitness values of each chromosome in the initial population are calculated to determine the initial value of the best chromosome of the parent generation.

[0064] Step Seven: Update the chromosome using the quantum rotation gate, and return to Step Two to continue iteratively executing the decoding of the next generation of chromosomes.

[0065] The specific update process is as follows:

[0066] Taking the best individual in the parent chromosome as the evolution target, use the quantum rotation gate to update the population, is the rotation gate angle of the i-th chromosome individual in the n ga -th generation.

[0067] The process of updating qubits using quantum rotation gates is expressed as:

[0068]

[0069] Among them, and are respectively the two probability amplitudes corresponding to the j-th gene position of the i-th chromosome in the n-th ga generation of evolution.

[0070] Taking into account both the convergence speed and the convergence accuracy, the following strategy is adopted to adaptively adjust the rotation gate angle:

[0071]

[0072] Among them, Δθ max and Δθ min are respectively the maximum and minimum rotation angles, N max is the maximum number of iterations, is the cosine similarity between the i-th chromosome in the n-th ga generation and the best chromosome of the parent generation, and w c is the weight coefficient of the cosine similarity.

[0073] The quantum NOT gate mutation process is expressed as:

[0074]

[0075] Among them, and are respectively the qubit positions before and after the quantum NOT gate mutation; U not is the quantum NOT gate,

[0076] Step 8: The gene sequence of the best chromosome in the current generation is the final optimized result of the offline satellite attitude maneuver sequence. Decode according to the conversion relationship between qubit encoding, binary encoding, and decimal encoding to extract the optimized attitude maneuver sequence.

[0077] Step 9: In the on-orbit application stage, the satellite executes the optimized attitude maneuver sequence, collects the measurement data of the gyroscope and star sensor, and processes the measurement data using the established gyroscope error parameter calibration model to estimate the final calibration value of the gyroscope error parameter.

[0078] The advantages of the present invention are:

[0079] 1. A rapid calibration method for gyro error parameters based on intelligent optimization maneuver of the present invention has good rapidity. By using the quantum genetic algorithm to iteratively optimize the satellite attitude maneuver sequence, the maneuver times in each axial direction can be reasonably allocated, and the observability of various gyro error parameters can be maximally improved under various established constraint conditions of the attitude determination and control system, thereby improving the calibration accuracy of the gyro error parameters.

[0080] 2. A rapid calibration method for gyro error parameters based on intelligent optimization maneuver of the present invention uses a variety of reasonably designed basic maneuver forms for the design of the attitude maneuver sequence, which can meet the excitation requirements for different types of gyro error parameters and is easy to implement in engineering.

[0081] 3. A rapid calibration method for gyro error parameters based on intelligent optimization maneuver of the present invention adopts additional evolutionary strategies such as adaptive adjustment of the rotation angle and quantum catastrophe on the basic quantum genetic algorithm, effectively avoiding the premature convergence of the population, greatly improving the convergence speed of the maneuver parameter optimization process, and providing the possibility for on-orbit implementation.

[0082] 4. A rapid calibration method for gyro error parameters based on intelligent optimization maneuver of the present invention establishes a fitness function based on the multiplicative error quaternion between the gyro and star sensor measurement outputs to optimize the maneuver sequence, which has certain engineering practicability. The present invention is easy to implement and is not only applicable to the error calibration of satellite attitude sensors, but also can be extended to other fields. BRIEF DESCRIPTION OF THE DRAWINGS

[0083] Figure 1 is a flow chart of a rapid calibration method for gyro error parameters based on intelligent optimization maneuver of the present invention;

[0084] Figure 2 is a schematic diagram of the satellite attitude maneuver form adopted by the present invention;

[0085] Figure 3 is a schematic diagram of the chromosome using integer coding adopted by the present invention;

[0086] Figure 4 is a schematic diagram of the satellite body coordinate system and the orbital coordinate system adopted by the present invention;

[0087] Figure 5 is a flow chart of the quantum genetic algorithm adopted by the present invention;

[0088] Figure 6 is a schematic diagram of the three-axis angular velocity of the satellite before performing the maneuver in the embodiment of the present invention;

[0089] Figure 7 is a schematic diagram of the optimization result of the maneuver sequence when the maneuver duration is 10 min in the embodiment of the present invention;

[0090] Figure 8 Schematic diagram of the optimized maneuver sequence when the maneuver duration is 4 minutes in the embodiment of the present invention;

[0091] Figure 9 Schematic diagram of the fitness convergence curve based on QGA in the embodiment of the present invention;

[0092] Figure 10 Schematic diagram of the fitness convergence curve based on IQGA in the embodiment of the present invention;

[0093] Figure 11 Schematic diagram of the conventional maneuver sequence before optimization in the embodiment of the present invention

[0094] Figure 12 Schematic diagram of the optimized maneuver sequence under multiple optimization objectives in the embodiment of the present invention. Detailed implementation manners

[0095] To facilitate the understanding and implementation of the present invention by those of ordinary skill in the art, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Obviously, the described embodiments are only part of the embodiments of the present invention, rather than all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0096] The present invention provides a method for quickly calibrating gyro error parameters based on intelligent optimized maneuvers. For a satellite attitude determination system composed of a gyro and a star sensor, under the conditions of meeting various constraints such as the dynamic range limit of attitude maneuvers, easy implementation, and the satellite returning to the normal working state after maneuvers, the basic maneuver form and the maneuver sequence coding method are designed; with the basic goal of fully stimulating the observability of each error term parameter of the on-board gyro, the quantum genetic algorithm is used to iteratively optimize the maneuver sequences in each axis; during the iterative evolution of the maneuver sequences, an adaptive adjustment strategy for the rotation angle is introduced into the quantum rotation gate update process, and additional optimization strategies such as quantum NOT gate mutation and quantum catastrophe are adopted to accelerate the convergence speed and avoid falling into local optimal solutions; considering the requirements of rapidity and accuracy for on-board gyro error calibration, a least squares calibration model for gyro error parameters is established based on the measurement models of the gyro and the star sensor; in order to effectively determine the evolution direction of the maneuver sequences, a fitness function based on multiplicative error quaternions is established using the measurement outputs of the gyro and the star sensor to assist in realizing the global search for the optimal maneuver sequences, which has certain engineering applicability; the optimized satellite maneuver sequence is extracted using the gene coding of the best individual in the final generation of evolution; finally, the satellite executes the optimized maneuver sequence in orbit to achieve the rapid calibration of gyro error parameters.

[0097] As Figure 1 shown, the specific steps are as follows:

[0098] Step 1: For a satellite with an attitude determination system composed of a gyroscope and a star sensor, determine the basic maneuver forms and parameters that make up the attitude maneuver sequence.

[0099] When the satellite performs attitude maneuvers to excite the observability of gyro error parameters, the following constraints are considered: parameter amplitude limits such as the maximum angular velocity and angular acceleration of the satellite attitude maneuver; the maneuver form should be simple and easy to implement; the sampling interval for the star sensor to perform effective measurements should be reserved during the maneuver; the satellite should reset to the normal working state after the attitude maneuver ends.

[0100] Based on the above conditions, satellite maneuver strategies are usually designed using maneuver angular rate change forms such as triangular, trapezoidal, and sinusoidal. The basic maneuver forms include six ways of changing the magnitude of the angular rate of triangular, trapezoidal, sinusoidal and their combinations, which are used as the basic components of the attitude maneuver sequence.

[0101] In order to enrich the maneuver form to fully stimulate the role of gyro error parameters in the measurement output under the condition of meeting the maneuver constraints, the present invention uses a method of compounding multiple maneuver forms to design the maneuver sequence, and optimizes the parameters of the maneuver sequence to maximize the observability of various error parameters of the gyro.

[0102] The specific attitude maneuver form is as Figure 2 shown. In the figure, η represents the maneuver angular rate of the satellite. The star sensor needs to perform at least one effective sampling in the maneuver pause interval after the forward maneuver and the reverse maneuver end. Assume that the durations of the forward maneuver and the reverse maneuver of the satellite around a certain axis are and respectively, and the time for pausing sampling between the forward and reverse attitude maneuvers is t s , then the complete single-round maneuver time period is

[0103] The parameters reflecting the attitude maneuver sequence include maneuver form, maximum maneuver angular rate amplitude, single maneuver period, single-direction maneuver duration, total number of maneuver rounds, and total maneuver time, etc.

[0104] Step 2: Develop a quantum genetic coding method to encode and decode the contemporary chromosomes of the quantum genetic algorithm for subsequent iterative optimization of the attitude maneuver sequences of each axis.

[0105] To represent the maneuver sequence more intuitively and simplify the code complexity of algorithm implementation, an integer coding method is adopted to represent the satellite's attitude maneuver sequence. Each chromosome represents the coding of a feasible satellite attitude maneuver sequence, and each gene of the chromosome is obtained by coding the maneuver parameters corresponding to the corresponding time period in the sequence. To represent the maneuver sequence intuitively, the present invention uses the decimal integer method for coding. The attitude maneuver parameters participating in coding and optimization include maneuver type, the axis where the current maneuver is implemented, the magnitude of the maneuver angular rate, etc.

[0106] As the types of maneuver parameters to be optimized increase, the gene length increases. The total length of chromosome coding will increase with the increase of gene length and maneuver sequence length, and the optimization convergence speed will also slow down accordingly.

[0107] Figure 3 The chromosome with gene length m is shown. The arrangement order of genes in the chromosome represents the sequence of satellite execution of each round of attitude maneuvers. The number of genes in each chromosome is determined by the number of maneuvers constituting the attitude maneuver sequence. For the chromosome with gene length m, the coding method of its k-th gene "abcd" is as follows

[0108] The thousands digit a is the coding of the basic maneuver form, 1 ≤ a ≤ n m , n m is the number of candidate maneuver forms. When using the Figure 2 shown maneuver form as the candidate maneuver form for optimization, n m = 6.

[0109] The hundreds digit b is the coding of the maximum maneuver angular rate magnitude η max , and the specific calculation formula is as follows:

[0110]

[0111] where η max and are respectively the lower limit and upper limit of the value range of the maximum maneuver angular rate magnitude η max .

[0112] In the decoding stage, the maximum maneuver angular rate magnitude η max is decoded according to the following formula:

[0113]

[0114] The tens digit c represents the maneuver axis: 1 corresponds to the x-axis, i.e., the pitch axis, 2 corresponds to the y-axis, i.e., the yaw axis, 3 corresponds to the z-axis, i.e., the roll axis, 1 ≤ c ≤ 3.

[0115] The units digit d is the flag encoding for whether to perform a maneuver in the current period. When d is odd, it indicates that there is an attitude maneuver. The maneuver flag is set to 1, and the relevant parameter decoding for the current maneuver is continued and the maneuver is executed. When d is even, it indicates that there is no attitude maneuver. The maneuver flag is set to 0, that is, no maneuver is executed at the current moment. On this basis, quantum encoding and collapse are performed, and further quantum genetic optimization is realized.

[0116] Convert the attitude maneuver sequence in the above decimal integer encoding form into binary encoding, then introduce the state vector expression of the quantum, and use quantum bits to encode the genes.

[0117] Introduce the state vector expression of the quantum into the genetic encoding, adopt the quantum bit encoding method, and use quantum bits to store and express a gene. This enables each chromosome to simultaneously express the superposition of multiple states, increasing the diversity.

[0118] The representation of the quantum bit is as follows:

[0119]

[0120] Among them, |0> and |1> respectively represent the spin-down state and the spin-up state. The quantum bit can simultaneously be in the superposition state of two quantum states. α and β are the probability amplitudes corresponding to the states |0> and |1> respectively, and satisfy the normalization condition |α| 2 +|β| 2 =1.

[0121] Thus, the chromosome structure is represented as:

[0122]

[0123] Among them, represents the chromosome of the j-th individual in the t-th generation, corresponding to a maneuver sequence; m is the number of genes of the chromosome, which corresponds to the number of cycles of a single-round maneuver in the executed maneuver sequence in the present invention; k is the number of quantum bits included in each gene. When initializing the quantum encoding, the probability amplitudes (α i , β i ) of all genes are set to That is, the value probabilities of all possible solutions are the same.

[0124] The gene encoded by the quantum bit enables each chromosome to simultaneously express the superposition of multiple states, increasing the diversity. During the decoding process, the quantum bit encoding is converted into binary encoding after quantum collapse, then the binary encoding is converted into decimal encoding, and finally the relevant information of the satellite attitude maneuver sequence is extracted from the decimal encoding according to the encoding protocol agreed upon above.

[0125] Quantum collapse is a probability measurement, that is, the quantum bit is taken with |α| 2The probability collapses to the classical bit 0 state with |β| 2 The probability collapses to the classical bit 1 state, expressed as follows:

[0126]

[0127] where m r = rand, and rand represents a random number between [0, 1].

[0128] Step 3: Establish a calibration model for gyro error parameters based on star sensor measurement assistance, and use the least squares method to estimate the gyro error parameters.

[0129] Based on the establishment of the measurement models of the gyro and star sensor and the attitude determination model, a calibration model for gyro error parameters adapted to the rapid attitude maneuver situation is constructed. Then, using the measurement data of the gyro and star sensor obtained during the rapid attitude maneuver of the satellite, the estimation of the gyro error parameters is realized.

[0130] By simplifying the gyro measurement model to different degrees, calibration models with different precisions between the measurement data of the gyro and star sensor and the gyro error parameters to be estimated can be established. Here, the attitude change amount measured by the high-precision star sensor over a period of time is used as a reference. Using the difference between it and the attitude change amount calculated from the gyro measurement angular velocity affected by various error sources, a calibration model for gyro error parameters is established under the attitude maneuver condition, and the estimation of the gyro error parameters is realized by solving this model. When solving the calibration model, optimization algorithms based on the minimum variance criterion such as least squares and Gauss-Newton, as well as intelligent optimization algorithms such as sparrow and ant colony, can be adopted.

[0131] From the perspective of improving the rapidity of gyro error calibration and ensuring the calibration accuracy at the same time, the least squares method is taken as an example to elaborate the method. Below, based on the establishment of the measurement models of two types of attitude sensors, the gyro and star sensor, a calibration model for gyro error parameters is established based on the measurement data of the gyro and star sensor, and the least squares method is used for estimation.

[0132] Define the satellite body coordinate system (b system) and the orbital coordinate system (o system), as Figure 4 shown. The x b y b z b axis of the satellite body coordinate system (ox b y b z b ) corresponds to the pitch axis, the y b axis corresponds to the yaw axis, and the z b axis corresponds to the roll axis. The coordinate plane of the orbital coordinate system (ox o y o z o ) is the satellite orbital plane, and the zo The axis points from the satellite's centroid to the Earth's center, and the x o axis points in the direction of the satellite's velocity within the orbital plane and is perpendicular to the z o axis, and the y o axis is parallel to the normal of the orbital plane and is right-handed orthogonal to the other two axes.

[0133] The unit vectors of the axes of the orbital coordinate system can be expressed in terms of the geocentric radius vector r and the velocity vector v of the satellite's motion as

[0134]

[0135] Then, the transformation matrix from the geocentric inertial coordinate system (i-system) to the orbital coordinate system is

[0136] The gyro error parameters include the constant drift b, the installation error Δ, the scale factor error δs, and the random measurement noise w g ; its measurement model can be represented by the combination of the true angular velocity and the above-mentioned error terms

[0137] The specific steps are as follows:

[0138] Step 301: Establish the measurement model of the three-axis gyro assembly as follows:

[0139]

[0140] where ω g is the angular velocity measured and output by the gyro assembly; ω b is the true angular velocity of the satellite's rotation relative to the inertial system in this system, is the installation matrix of the gyro assembly; b is the constant drift vector; Λ is a diagonal matrix with the components of the scale factor error δs as the diagonal elements, i.e., the error matrix; is the installation error matrix formed by the installation error angles of the gyro assembly; when the installation error angle Δ is a small angle, 9Δ×] represents the skew-symmetric matrix formed by the components of Δ; w g is the random measurement noise vector; I is the identity matrix.

[0141] Assume that the scale factor errors of the gyro assembly on the x, y, and z axes are δs x , δs y , and δs z , respectively. Then, the diagonal matrix formed by the scale factor error vector δs = [δs x δs y δs z T is Λ; assume that the installation error angles of the gyro around the x, y, and z axes are and are all small angles. From the formed anti-symmetric matrix is [Δ×], then the installation error matrix can be approximated as Therefore, the relationship between the true angular velocity and the gyro-measured angular velocity can be established as:

[0142]

[0143] Without loss of generality, assume that the sensitive axes of the gyro components are all installed along the body axes, then the installation matrix is the identity matrix. Substituting it into the above formula, after arrangement and neglecting higher-order small quantities, the true rotational angular velocity of the satellite can be approximated as:

[0144] ω b ≈(I - [Δ×] + Λ)(ω g - b - w g ) (9)

[0145] Step 302: Considering the random measurement noise as the main error source of the star sensor, establish the measurement model of the star sensor;

[0146] The star sensor is the most accurate sensor in the current satellite attitude determination system. According to needs, it can output satellite attitude information in various forms such as Euler angles, quaternions, and starlight vectors, and these forms can be converted to each other. In this example, the satellite measurement information is output in the form of quaternions. Due to factors such as star trailing in the dynamic situation, which will affect the accuracy of star map processing, therefore, in order to ensure the attitude measurement accuracy of the star sensor, it should avoid operating in a large attitude maneuver state.

[0147] In the gyro error calibration stage, since the calibration time of the gyro error parameters is much shorter than the low-frequency error period of the star sensor, the low-frequency error and other long-period error sources in the star sensor measurement have little impact on the measurement accuracy, and the installation error of the star sensor has been calibrated in advance. Therefore, mainly considering the random measurement noise as the main error source of the star sensor, the measurement model is established as follows:

[0148]

[0149] where Q st is the attitude quaternion measured and output by the star sensor; is the true attitude quaternion of the satellite, that is, the attitude representation of the satellite body coordinate system relative to the inertial system; is the quaternion corresponding to the star sensor installation matrix; Q ξ is the random measurement noise ξ represented in the form of quaternions.

[0150] Step 303: Construct a gyro error parameter calibration model by using the measurement models of the three-axis gyro assembly and the star sensor;

[0151] According to the satellite attitude kinematic equation, the relationship between the satellite angular velocity ω b and the attitude quaternion can be expressed by the following quaternion differential equation

[0152]

[0153] where Ξ(ν) represents the matrix formed by the vector ν = [ν 1 ν 2 ν 3 T Similarly, Ξ(ω ) is the matrix formed by ω b b .

[0154] Since the satellite attitude change during normal operation is relatively stable, especially for a three-axis stabilized satellite used for earth-pointing observation or communication, it can be considered that the attitude angular velocity of the satellite remains basically unchanged within a short time interval Δt. The Picard approximation method is used to solve the above satellite attitude kinematic equation. According to the attitude quaternion at time t, using the satellite angular velocity sequence values within (t, t + nΔt), the attitude quaternion at time t + nΔt can be approximately calculated as

[0155]

[0156] Let {ω b} represent the angular velocity sequence Define

[0157]

[0158] The attitude change quaternion within the time period (t, t + nΔt) is

[0159]

[0160] From (12) and (14) above, the relationship between and {ω b} can be established as follows

[0161]

[0162] where Q 0 is the attitude quaternion at the initial time.

[0163] ​​According to the Picard solution method, the angular velocity exponential matrix can be approximately expressed by the angular increment as follows:

[0164]

[0165] Among them, and are the coefficients corresponding to the Picard solution method, which are determined by the angular increment within Δt, where Thus, we can obtain For different orders of the Picard solution algorithm, different and can be approximately obtained. When using the third-order Picard solution algorithm,

[0166] Substituting Equation (16) into Equation (13) and expanding, we get:

[0167]

[0168] Among them,

[0169] According to Equation (9), ignoring the high-order small terms of the second order and above, the angular velocity of the satellite at time t can be obtained as:

[0170]

[0171] Among them, is the equivalent measurement noise. It can be seen from Equation (16) that the values of the Picard solution coefficients and are related to . Since cannot be obtained and is determined here, determined by is approximately used to replace for calculating the coefficients:

[0172]

[0173] A calibration model for gyro error parameters is established, specifically:

[0174] Using the measurement information (Q st ) t and (Q st ) t+nΔt obtained by the star sensor at times t and t + nΔt; calculate the attitude change amount

[0175] of the star sensor during [t, t + nΔt). Among them

[0176] Then, according to the angular velocity sequence {ω b} measured by the three-axis gyroscope during the period [t, t + nΔt), the attitude quaternion variation G({ω g})Q 0 affected by the gyroscope measurement error is recursively calculated using the Picard integral method, where Q 0 is the attitude quaternion at the initial moment; G({ω g}) represents the consecutive multiplication of the angular velocity exponential matrices.

[0177] Finally, using the attitude variation of the star sensor and the attitude quaternion variation of the gyroscope measurement error, the gyroscope error parameter calibration model is constructed as follows:

[0178]

[0179] where Δt is the time interval of measurement sampling; Δ is the gyroscope installation error; δs is the gyroscope scale factor error; In the formula, n is the total number of samplings, and are the coefficients corresponding to the Picard solution method; is the angular velocity measured by the gyroscope at the moment t + jΔt; represents the formed skew-symmetric matrix; is the diagonal matrix formed by the elements of ; V t represents the model equivalent noise caused by gyroscope random measurement noise, star sensor random measurement noise, model linearization error, etc.

[0180] In each term on the right side of the above formula, the first row of M 0 is a zero vector, so the first element of the left column vector is 0, which is not affected by the gyroscope error parameters. The role of the gyroscope error parameters is reflected by the vector part of the error quaternion. Therefore, only the vector part is used to establish the calibration model.

[0181] Step 304: Simplify the gyroscope error parameter calibration model by using the measurement information of the gyroscope and the star sensor obtained under the condition of continuous excitation of the attitude maneuver sequence;

[0182] Substitute the attitude measurement values successively output by the star sensor at times t k and (t k + n k Δt) and the sequence {ω k} composed of n g gyroscope measurements obtained within this time interval into Equation (20), and take the vector part to establish the calibration equation of the error parameters as follows:

[0183]

[0184] where X = [b T Δ T δs T T is the vector of gyro error parameters to be calibrated; is the vector part on the left side of Equation (20); is the vector part of the equivalent noise term on the right side of Equation (20), is the measurement matrix:

[0185]

[0186] Step 305: Collect the measurement data of the star sensor and the gyro, and establish an error calibration equation for the entire attitude maneuver period interval;

[0187] If during the maneuver pause period of the entire attitude maneuver interval, the star sensor makes (N + 1) effective measurements at times t k (k = 1, 2,..., N + 1); during the period between two adjacent star sensor measurements [t k , t k+1 , the gyro measures n k times, and at this time t k+1 = t k + n k Δt.

[0188] Based on this, multiple calibration equations are established and combined to obtain the error calibration equation for the entire maneuver period as follows:

[0189] Z = HX + v (22) where,

[0190] Step 306: In the case where the noise statistical characteristics are difficult to accurately obtain, use the least squares criterion in combination with the error calibration equation to find the parameter values that minimize the residual of the gyro error parameter calibration model, and thus obtain the estimation of the gyro error parameters to be calibrated, that is, the estimated value of the parameter vector

[0191] The formula is as follows:

[0192]

[0193] where the estimated values of each error parameter can be obtained from the estimated value, including the estimated value of the constant bias the estimated value of the installation error and the estimated value of the scale factor error

[0194] ​Step 4: Establish a satellite angular velocity model in the presence of attitude maneuvers; and use it to simulate and generate satellite angular velocity information under different maneuver sequences during the offline optimization process.

[0195] During the offline optimization process, calculate the actual angular velocity of the satellite under additional attitude maneuvers based on the maneuver sequence extracted from a certain contemporary individual as follows:

[0196]

[0197] where, is the projection of the angular velocity of the orbital system relative to the inertial Earth system in the inertial Earth system, i.e., the angular velocity when the satellite is operating normally; is the transformation matrix from the inertial Earth system to the orbital coordinate system; is the transformation matrix from the orbital coordinate system to the satellite body system; is the projection of the angular velocity of the satellite body system relative to the orbital system in the body system, representing the angular velocity when the satellite executes additional attitude maneuver sequences around the body axes in addition to the attitude angle motion during normal satellite operation; is the projection of the angular velocity of the satellite body system relative to the inertial Earth system in the satellite body system.

[0198] During the offline optimization maneuver phase, historical satellite operation data over a period of time or desired on-orbit motion data can be used as During the on-orbit error calibration phase, the angular velocity of the satellite during normal operation is

[0199] During the offline optimization process, the angular velocity of the satellite with attitude maneuvers is obtained using Equation (24) Then, according to the established gyro measurement model, simulate and generate measurement values as follows:

[0200]

[0201] where, and are respectively a set of empirical reference values for the gyro scale factor error diagonal matrix, installation error angle vector, and constant bias vector.

[0202] Taking into account various factors such as the accuracy of prior information, the cognitive level of engineering personnel, and the complexity of the actual satellite operating environment, multiple sets of gyro error reference values can be randomly generated and substituted into the above formula to generate multiple sets of gyro measurement data. On this basis, a new fitness value is generated by weighted fusion of the fitness values obtained based on each set of measurement data, and on this basis, the maneuver sequence is iteratively optimized.

[0203] Theoretically speaking, the fusion fitness value calculated based on the measurement data of multiple groups of error reference values can better reflect the best direction of evolution. However, the computational complexity will increase significantly and the convergence time will be long. Therefore, in order to reduce the computational complexity and accelerate the convergence speed, this embodiment takes a single group of error reference values as an example to carry out the method.

[0204] Step Five: Establish a fitness function optimized by quantum genetics using the gyro error parameter calibration model and the satellite angular velocity model; and calculate the fitness function values corresponding to different attitude maneuver sequences during the gyro error parameter estimation accordingly during the iterative evolution process.

[0205] In practical engineering applications, satellite attitude maneuvers should minimize the requirements for the dynamic performance of the attitude control system. That is, the maneuver angular velocity, maximum maneuver amplitude, maneuver time, etc. of the satellite should meet various limiting conditions. At the same time, it should also ensure sufficient excitation of various sensor errors to improve the observability of error parameters.

[0206] Considering comprehensive indicators such as calibration accuracy, maneuver cost, and rapidity, the fitness function is as follows:

[0207]

[0208] Among them, trace(·) represents the trace operation of a matrix; t m is the total time consumed during the satellite maneuver, with the unit of s; θ m is the total rotation angle during the satellite maneuver, with the unit of degree; η m is the maximum maneuver angular rate amplitude during the satellite maneuver, with the unit of degree / s; is the time for the satellite to run at the maximum maneuver angular rate amplitude, with the unit of s; w o 、w t 、w θ and w η are weight coefficients reflecting the contribution degrees of the observable level, maneuver time consumption, maneuver rotation angle, and maneuver angular rate amplitude in the fitness function.

[0209] Considering various costs in practical engineering applications, the fitness function described in Equation (26) can be adaptively adjusted.

[0210] On the premise of meeting the requirements of engineering applications, the values of some maneuver parameters can be determined in advance based on prior knowledge. This can reduce the number of parameters to be optimized, simplify the optimization problem, and thus shorten the time for optimization solution. For example, the maximum maneuver angular velocity that the satellite can execute can be used as the angular rate amplitude for a single maneuver, which can increase the excitation intensity and meet the requirement of quickly changing the attitude angle; according to the requirement for rapidity during the on-orbit calibration process, the total duration of the maneuver and the time interval for a single maneuver can also be reasonably limited in advance.

[0211] With the total maneuvering time and the maximum maneuvering angular rate amplitude fixed in advance, during the optimization process, the observability index trace((H T H) -1 ) does not change significantly. At this time, an objective function reflecting the calibration accuracy is designed based on the multiplicative error quaternion to replace the first term on the right side of equation (26). According to equation (8), the angular velocity represented by the calibrated error parameters can be obtained as follows:

[0212]

[0213] Using the angular velocity compensated by the calibrated error parameters within the time period (t k , t k + n k Δt) and the satellite attitude quaternion at time t k According to equation (12), the attitude estimation at time t = t k+1 = t k + n k Δt is:

[0214]

[0215] Using the star sensor in the attitude maneuvering interval t 1 , …, t N Measured The multiplicative error quaternion at the corresponding time can be obtained Based on this, the fitness function based on the multiplicative error quaternion is constructed as follows:

[0216]

[0217] Among them, is The vector part of, w q is the transformation factor determined according to the accuracy level of the attitude sensor. In the subsequent embodiments, take w q = 10 6 .

[0218] Step Six: Select the chromosome with the highest fitness function value from the contemporary chromosomes, and determine whether the fitness difference between the contemporary best chromosome and the parental best chromosome is less than the preset threshold, or whether the number of iterative evolutions has reached the preset maximum value. If so, end the evolutionary process offline and transfer to Step Eight; otherwise, execute Step Seven to achieve the global search for the optimal maneuvering sequence.

[0219] Before the iterative optimization, the chromosome population is initialized manually. By calculating the fitness values of each chromosome in the initial population, the initial value of the parental best chromosome is determined.

[0220] Step 7: Update the chromosome using the quantum rotation gate, and return to Step 2 to continue iterating and performing the decoding of the next-generation chromosome.

[0221] By adopting the quantum genetic algorithm and replacing the selection, crossover, and mutation operators in the traditional genetic algorithm with the quantum rotation gate operator, only the optimal individual needs to be retained, which greatly improves the search and optimization efficiency and reduces the probability of premature convergence of the population.

[0222] Combine the calculation of the adaptive rotation angle to update the chromosomes of the contemporary population using the quantum rotation gate, then perform quantum NOT gate mutation with a set probability, and then judge that the algorithm enters premature convergence and take timely measures according to the fact that the fitness of the best individual has not changed for a long time or is lower than the preset quantum catastrophe threshold for a certain period of time. At the same time, introduce the quantum catastrophe strategy to force the algorithm to jump out of the local optimal solution. After completing the above evolutionary operations, return to execute the decoding process of the next-generation chromosome in Step 2.

[0223] The specific update process is as follows:

[0224] Take the best individual in the parental chromosome as the evolutionary goal and use the quantum rotation gate to update the population, where ga is the rotation gate angle of the i-th chromosome individual in the n-th

[0225] The process of updating the quantum bit using the quantum rotation gate is expressed as:

[0226]

[0227] where and are respectively the 2 probability amplitudes corresponding to the j-th gene position of the i-th chromosome in the n-th ga generation of evolution.

[0228] The magnitude of the rotation angle is a key factor affecting the convergence accuracy and convergence speed of the optimization algorithm. An overly small rotation angle will prolong the convergence time, while an overly large rotation angle is likely to cause the algorithm to fall into the local optimal solution, thereby affecting the convergence accuracy. Therefore, considering both the convergence speed and convergence accuracy, the following strategy is adopted to adaptively adjust the rotation gate angle:

[0229]

[0230] where Δθ max and Δθ min are respectively the maximum and minimum rotation angles, N max is the maximum number of iterations, is the cosine similarity between the i-th chromosome in the n-th ga generation and the best parental chromosome, and w cis the weight coefficient of cosine similarity.

[0231] In the early stage of evolution, a larger rotation angle can be adopted to achieve a large-range and fast search. In the later stage of evolution, a smaller rotation angle is adopted to achieve a more accurate search in a small range. By introducing cosine similarity as an index to measure the similarity between individuals, the individuals evolve towards the best direction with more reasonable self-performance based on the quantum coding method.

[0232] Denote the rotation direction of the quantum bit as Determine the rotation direction according to the conditions shown in Table 1.

[0233] Table 1

[0234]

[0235] In the table, represents the j-th gene position of the i-th individual in the n-th ga generation of evolution, represents the j-th gene position of the best individual in the (n ga -1)-th generation, and respectively represent the fitness values corresponding to the i-th individual in the n-th ga generation and the best individual in the (n ga -1)-th generation.

[0236] On the basis of updating the chromosome by the quantum rotation gate, quantum mutation is adopted to make some individuals deviate from the current evolution direction to avoid falling into the local optimal solution. Specifically, traverse all the individuals in the current population, and use the quantum NOT gate to exchange the probability amplitudes of the quantum bits according to a predetermined probability to reverse the evolution direction of the individuals. The quantum NOT gate mutation process can be expressed as:

[0237]

[0238] Among them, and are the quantum bit positions before and after the quantum NOT gate mutation respectively; U not is the quantum NOT gate,

[0239] To avoid premature convergence of the optimization process, a quantum catastrophe strategy is further introduced.

[0240] The adopted catastrophe strategy is that when the fitness value of the best individual in the evolution process remains unchanged for N e generations, first sort all the individuals in the current population from high to low according to the fitness value, and then randomly re-initialize the last one-third of the inferior individuals to force the algorithm to jump out of the local optimal solution. N max can be determined according to N e , Ne = λ e N max , according to experience, the coefficient λ e is taken as 0.1. After evolutionary operations such as updating the chromosome through the quantum rotation gate, continue to iteratively execute the decoding operation of the next-generation chromosome in Step 2.

[0241] To distinguish it from the conventional QGA, the quantum genetic optimization algorithm improved by the above adaptive rotation angle, quantum NOT gate mutation, and quantum catastrophe strategy is denoted as IQGA (Improved Quantum Genetic Algorithm). The above process of quantum genetic optimization for the satellite attitude maneuver sequence is as Figure 5 shown.

[0242] Step 8: The gene sequence of the current best chromosome is the final optimized result of the offline satellite attitude maneuver sequence. Decode it according to the conversion relationship between qubit encoding, binary encoding, and decimal encoding to extract the optimized attitude maneuver sequence.

[0243] In the offline optimization maneuver phase, use a period of satellite operation historical data or expected on-orbit motion data as The satellite angular velocity with attitude maneuvers is obtained according to Equation (24) Then, using the gyro measurement model Equation (7), the empirical reference values of the gyro scale factor error, installation error angle, and constant bias can be simulated and generated and under the

[0244]

[0245] Considering various factors such as the accuracy of prior information, the cognitive level of engineering personnel, and the complexity of the actual satellite operating environment, multiple groups of gyro error reference values can be randomly generated. Substitute them into Equation (7) to generate gyro measurement data under multiple groups of gyro error parameters. On this basis, the maneuver sequence is iteratively optimized according to the fitness value calculated by Equation (26).

[0246] Based on the convergence of the fitness values of the current population and the number of iteration steps, a discrimination condition is established for terminating the evolutionary iteration: If the current iteration generation number n GA < N max , N max is the maximum number of iterations, and the difference between the fitness values of the best individuals in two adjacent generations and is greater than the threshold ζ, that is then continue the iterative optimization; otherwise, the iterative evolution ends, and the gene sequence of the current best individual is the final optimized result of the offline satellite attitude maneuver sequence.

[0247] After the iterative evolution ends, convert the quantum encoding of the best individual in the last generation of the population into binary encoding, and then convert the binary encoding into decimal encoding to determine the optimized attitude maneuver sequence.

[0248] Step Nine: In the on-orbit application stage, the satellite executes the optimized attitude maneuver sequence, collects the measurement data of the gyroscope and star sensor, and processes the measurement data using the established gyroscope error parameter calibration model to estimate the final calibration value of the gyroscope error parameter.

[0249] The satellite executes the optimized attitude maneuver sequence on orbit, simultaneously collects the measurement data of the gyroscope and star sensor, then executes the gyroscope error parameter calibration modeling method and error parameter estimation method, and thus can estimate the gyroscope error parameter. Use the estimated parameter to compensate for the measurement error of the gyroscope to improve the on-orbit measurement accuracy of the gyroscope.

[0250] Embodiment

[0251] In summary, the present invention is a method for quickly calibrating gyroscope error parameters based on intelligent optimized maneuvers, and the specific steps are as follows:

[0252] Step One: Determine the design constraints such as the basic maneuver forms, parameters, effective parameter ranges, longest maneuver execution time, and maneuverable axes that make up the attitude maneuver sequence;

[0253] Step Two: Develop a quantum genetic encoding method for the attitude maneuver sequence to be optimized, and encode the chromosomes of the quantum genetic algorithm accordingly.

[0254] According to the decimal encoding method of the attitude maneuver sequence, establish the conversion relationship between decimal encoding, binary encoding, and quantum bit encoding. During the iterative evolution process, use this conversion relationship to realize the decoding from quantum genetic encoding to the attitude maneuver sequence.

[0255] Step Three: Establish a gyroscope error parameter calibration model assisted by star sensor measurement, and use the model to solve the gyroscope error parameter by methods such as least squares.

[0256] On the basis of establishing the measurement models of the gyroscope and star sensor and the attitude determination model, a gyroscope error parameter calibration model adapted to the rapid attitude maneuver situation is constructed, and then the measurement data of the gyroscope and star sensor during the rapid attitude maneuver of the satellite is used to realize the estimation of the gyroscope error parameter.

[0257] Step Four: Establish a satellite angular velocity model in the presence of attitude maneuvers, and simulate and generate the satellite angular velocity information under different maneuver sequences during the offline optimization process.

[0258] Step 5: Establish a fitness function optimized by quantum genetics, and calculate the fitness function values for gyro error parameter estimation corresponding to different attitude maneuver sequences during the iterative evolution process. On this basis, determine the next task of the optimization process according to the judgment of the iterative evolution termination condition. When the number of iterative evolution generations reaches the preset maximum value, or the fitness difference between the current best individual and the best individual of the parent generation is less than the preset threshold, end the evolution process and proceed to Step 7; otherwise, continue to execute Step 6.

[0259] Step 6: Determine the chromosome update strategy and update the chromosomes of the current population accordingly. After completing the above evolutionary operations, perform the decoding process in Step 2.

[0260] Update the chromosomes of the current population using a quantum rotation gate combined with an adaptive rotation angle calculation strategy, then perform quantum NOT gate mutation with a certain probability, and judge that the algorithm enters premature convergence and take timely measures according to the fact that the fitness of the best individual has not changed for a long time or is lower than the preset quantum catastrophe threshold for a certain period of time. At the same time, introduce a quantum catastrophe strategy to force the algorithm to jump out of the local optimal solution.

[0261] Step 7: Determine the optimized attitude maneuver sequence. Take the gene coding of the best individual in the final generation population of the iterative evolution, and extract the optimized maneuver sequence using the conversion relationship between quantum bit coding, binary coding, and decimal coding.

[0262] Step 8: Execute the optimized fast attitude maneuver in orbit to calibrate and compensate for the gyro error parameters.

[0263] The satellite executes the optimized attitude maneuver sequence, collects the measurement data of the gyro and star sensor, and processes the measurement data using the established gyro error parameter calibration model to estimate the final calibration value of the gyro error parameters.

[0264] Taking a Geostationary Earth Orbit (GEO) satellite as an example, conduct simulation experiments for verification and analysis. The specific orbital parameters of the selected GEO satellite are shown in Table 2. Use a satellite orbit simulation toolkit to achieve high-precision simulation of the satellite orbit operation state, and generate the position, velocity, and attitude motion data of the satellite during on-orbit operation.

[0265] Table 2

[0266] Semi-major axis of orbit / km Eccentricity Orbit inclination / ° Argument of perigee / ° Right ascension of ascending node / ° True anomaly / ° 42371.00 0.00 0.00 91.45 0.00 360.00

[0267] The upper limit of the attitude maneuver angular velocity of each axis of the satellite is 2° / s; the forward maneuver time and the reverse maneuver time Both are 10 s. A pause for additional maneuvers is made between forward and reverse maneuvers, and the time for the star sensor to make measurements is ts 2 s. Then, the time T for a single-axis to perform a single complete maneuver is 24 s. The simulation experiment parameters of the spaceborne gyro and the star sensor are shown in Table 3.

[0268] Table 3

[0269]

[0270] Under the above conditions, using the Monte Carlo experiment method, multiple groups of simulation experiments are carried out under different noise sequences. In each group of experiments, the sensor measurement data generated by the optimized attitude maneuver sequence are used to calibrate the gyro error parameters, and finally, the statistical analysis of the calibration results of each error parameter is carried out.

[0271] Let the number of Monte Carlo simulation experiments be N mc , and denote the calibration result and the true value of the i-th error parameter in the n-th experiment as and X i , respectively. Then, the relative error of calibrating the i-th error parameter in this experiment is:

[0272] J n,i =|ΔX n,i / X i |×100% (i = 1, …, 9, n ≤ N mc )

[0273] Among them, Then, according to the N mc experimental results, the mean value of the relative error of the i-th error parameter calibrated is:

[0274]

[0275] The mean value of the relative estimation errors of all error parameters is:

[0276]

[0277] Here, the number of Monte Carlo simulation experiments N mc = 10.

[0278] According to the satellite orbit parameters listed in Table 2, the three-axis angular velocities during satellite operation are simulated as Figure 6 shown. It can be seen from Figure 6 that when this satellite is in orbit, in order to maintain the earth-pointing orientation state, it rotates around the x-axis, i.e., the pitch axis, at a certain angular velocity, while there is no obvious angular motion information on the y-axis and the z-axis, i.e., the yaw axis and the roll axis, and the angular velocity magnitude is only of the order of 10 -3 times that of the x-axis angular velocity.

[0279] To simplify the optimization problem, only the maneuver form and the applied axial direction are taken as the parameters to be optimized first. When the total duration of satellite attitude maneuvers is 10 minutes and 4 minutes respectively, the maneuver sequences optimized by the quantum genetic algorithm are as shown in Figure 7 and Figure 8 respectively. In the figure, η x , η y , η z are the maneuver angular velocities of the x, y, and z axes respectively; f s is the star sensor sampling flag, and f s =1 indicates that there is an effective measurement output from the star sensor at this moment.

[0280] As can be seen from the maneuver optimization results shown in Figure 7 and Figure 8 , the number of maneuver rounds in the x - axial direction is the least, and the number of maneuver rounds in the y - axial direction is the most. This is because when the GEO satellite operates normally in orbit in the three - axis stabilized earth - pointing mode, its rotational angular velocity basically remains unchanged or changes slowly. In the case of an orbital inclination of 0°, the satellite attitude change is manifested as a large change in the pitch angle, while the yaw angle and roll angle almost remain unchanged. According to the satellite angular velocity signal before maneuver shown in Figure 6 , since the satellite itself has a certain operating angular velocity around the x - axis, only a few rounds of attitude maneuvers are required to meet the observability requirements for calibrating the gyro error parameters in this axial direction. When performing more rounds of attitude maneuvers on the y - axis and z - axis, the scale factor error and installation error terms of the gyro in the corresponding axial directions can be effectively excited to produce different forms of influence in the measurement output; at the same time, the number of effective measurement samples of the star sensor increases before and after the maneuver in this axial direction, effectively improving the observability of the gyro installation error and scale factor error in this axial direction.

[0281] Considering the rapidity requirement for error calibration during satellite on - orbit operation, subsequent simulation experiments are all carried out with the total duration of satellite attitude maneuvers being 4 minutes as an example. During the process of intelligently optimizing the maneuver sequence parameters, the fitness convergence curves based on the QGA and IQGA algorithms are as shown in Figure 9 and Figure 10 respectively. It can be seen that the QGA algorithm and the IQGA algorithm achieve the convergence of the fitness curve at about 150 and 30 generations respectively, thus verifying that the adoption of the adaptive rotation angle strategy plays a significant role in accelerating convergence and, to a certain extent, avoiding the phenomenon of premature convergence of the population.

[0282] To verify the advantage of the proposed intelligent optimization maneuver sequence determination method in improving the estimation accuracy of gyro error parameters, the maneuver method that forms a conventional maneuver sequence based on a single maneuver form (the satellite performs trapezoidal maneuvers in the order of x-y-z around the three main axes of the carrier system) will be used as the control group, and will be compared and analyzed with the quantum genetic intelligent optimization maneuver sequence method based on the composite maneuver form.

[0283] The relative errors of the gyro error parameters calibrated by the two methods are listed in Table 4. Here, the maximum maneuver angular velocity, maneuver duration, total number of maneuver rounds, etc. are kept consistent under the two maneuver sequence acquisition methods. It can be seen from the data in the table that compared with the conventional single-form maneuver sequence, when the optimized maneuver sequence based on the composite maneuver form is adopted, the calibration of the gyro error parameters has a certain degree of improvement in the calibration accuracy of each error term and its error components in each axial direction. Especially for the calibration accuracy of the three-axis constant bias error, y-axis installation error, scale factor error, and z-axis scale factor error, the relative error of the error parameter estimation drops significantly, which verifies that the intelligent optimization maneuver of the method proposed in the present invention can more fully stimulate the observability of various gyro error parameters while ensuring the rapidity of on-orbit estimation, and thus effectively improve the gyro calibration error accuracy.

[0284] Table 4

[0285]

[0286] When using the composite objective function model for iterative evolution, the error calibration scheme based on maneuver intelligent optimization is compared with the error calibration scheme based on the conventional maneuver method to verify the effectiveness of the proposed method in the case of multi-objective parameter optimization. The maximum maneuver time is set to 10 min, and the maneuver form, maneuver axis, maneuver angular rate amplitude, and maneuver time are used as the parameters to be optimized. In the control group, the satellite performs triangular maneuvers around the three main axes of the carrier system in the order of x-y-z, and the specific form of the maneuver sequence is as Figure 11 shown. Based on experience, the weight coefficients of the evaluation indexes of the fitness function are set. Here, w o =5, w t =0.1, w θ =0.7, w η =1. The obtained maneuver sequence after the optimization is as Figure 12 shown. It can be seen that the number of maneuver rounds of the optimized maneuver sequence is significantly reduced compared with the conventional maneuver sequence as Figure 11 shown, and the amplitude of the maximum maneuver angular velocity is also slightly reduced.

[0287] Table 5 statistics the calibration effects of various gyro error parameters under different maneuver schemes. Among them, Scheme A represents the intelligent optimization maneuver scheme, and Scheme B representsFigure 11 The conventional maneuvering scheme shown, where Scheme C represents the use of conventional maneuvering methods. The total number of maneuvering rounds and the maneuvering time are the same as the parameters of Scheme A, but the maneuvering amplitude uses the maximum angular rate amplitude.

[0288] Analyzing the data in the table, it can be seen that since Scheme B has the most maneuvering rounds and the longest maneuvering time, the observability excitation degree of its error is the largest, and the calibration accuracy of each gyro error parameter is the highest. Scheme C uses the optimized maneuvering time and the number of maneuvering rounds, and the maneuvering amplitude uses the maximum angular rate amplitude, making its error calibration parameters in some axes slightly better than those of Scheme A. However, from the overall calibration effect, the calibration effect of Scheme A in multiple axes is better than that of Scheme C, which verifies the necessity of the maneuvering form.

[0289] Table 5

[0290]

[0291] Table 6 statistics the three evaluation indicators of the total maneuvering time, the total maneuvering angle, and the mean relative error under two maneuvering schemes. Analyzing the data in the table, it can be seen that although the relative calibration accuracy of the gyro error slightly decreases due to the reduction in the number of maneuvering rounds in the optimized maneuvering sequence, the duration and the corner cost of the maneuvering process are improved by about 60% and 45% respectively, greatly reducing the cost of satellite maneuvering and avoiding unnecessary resource waste. Since the fitness function used is a comprehensive evaluation index composed of the observability level of error parameters, maneuvering time, amplitude, and corner cost, the data in the table shows that the obtained result is an optimized result of comprehensively weighing each evaluation index.

[0292] Table 6

[0293]

Claims

1. A method for rapid calibration of gyro error parameters based on intelligent optimization maneuvers, characterized in that: The specific steps are as follows: Step 1: For the satellite's attitude determination system, initialize the basic maneuver forms and parameters that constitute the attitude maneuver sequence; and formulate a quantum genetic coding method to encode and decode the contemporary chromosomes of the quantum genetic algorithm; Decimal integer encoding is used to convert the decimal integer encoding form of the posture maneuver sequence into binary encoding, and then the quantum state vector expression is introduced, and quantum bits are used to encode the gene; Step 2: Establish a gyro error parameter calibration model based on star sensor measurement assistance, and use the least squares method to estimate the gyro error parameters; The specific steps are as follows: Step 201: Establish a measurement model of the three-axis gyro assembly as follows: Among them, ω g The angular velocity measured by the gyro component; ω b is the true angular velocity of the satellite relative to the inertial system in this system; is the installation matrix of the gyro assembly; b is the constant drift vector; Λ is a diagonal matrix with the scale factor error δs component as the diagonal element; is the installation error matrix of the gyro assembly; w g is the random measurement noise vector; I is the unit matrix; Step 202: Establish a measurement model of the star sensor as follows: Among them, Q st The attitude quaternion output by the star sensor measurement; is the true attitude quaternion of the satellite, that is, the attitude representation of the satellite system relative to the inertial system; The quaternion corresponding to the star sensor installation matrix; Q ξ is the random measurement noise ξ expressed in quaternion form; Step 203: construct a gyro error parameter calibration model using the measurement model of the three-axis gyro assembly and the measurement model of the star sensor; First, the measurement information (Q st ) t and(Q st ) t+nΔt ; Calculate the attitude change of the star sensor during [t, t+nΔt) in Then, according to the angular velocity sequence {ω b }, the attitude quaternion change G({ω g })Q0, where Q0 is the attitude quaternion at the initial moment; G({ω g }) represents the continuous multiplication of the angular velocity exponent matrix; Finally, the gyro error parameter calibration model is constructed using the attitude change of the star sensor and the attitude quaternion change of the gyro measurement error as follows: Wherein, Δt is the time interval of measurement sampling; Δ is the gyro installation error; δs is the gyro scale factor error; Where n is the total number of sampling times, and is the coefficient corresponding to the Picard solution; is the angular velocity measured by the gyroscope at time t+jΔt; Indicated by The antisymmetric matrix formed; Reason The diagonal matrix of the elements of t represents the model equivalent noise caused by gyro random measurement noise, star sensor random measurement noise, and model linearization error; Step 204: Simplify the gyro error parameter calibration model to obtain: Where X = [b T Δ T δs T ] T is the gyro error parameter vector to be calibrated; The vector part on the left side of the model is used to calibrate the gyro error parameters; is the vector part of the equivalent noise term on the right side of the gyro error parameter calibration model, is the measurement matrix; Step 205: collect the measurement data of the star sensor and the gyroscope, and establish the error calibration equation for the entire attitude maneuvering period; The star sensor is at t k (k=1,2,…,N+1) valid measurements are performed at time (k=1,2,…,N+1); during the period between two consecutive star sensor measurements [t k ,t k+1 ], gyro measurement n k times, at this time t k+1 =t k +n k Δt; Based on this, multiple calibration equations are established and combined to obtain the error calibration equation for the entire maneuvering period as follows: Z=HX+v in, Step 206, using the least squares criterion combined with the error calibration equation, find the parameter value that minimizes the residual of the gyro error parameter calibration model, and obtain an estimate of the gyro error parameter to be calibrated; The formula is as follows: Among them, the estimated value From this, we can get the estimated values ​​of various error parameters, including the constant deviation estimate Installation error estimation and scale factor error estimates Step 3: Establish the satellite angular velocity model in the presence of attitude maneuvers; The specific formula is as follows: in, It is the projection of the angular velocity of the orbital system relative to the Earth's inertial system in the Earth's inertial system, that is, the angular velocity of the satellite during normal operation; is the transformation matrix from the Earth's inertial system to the orbital coordinate system; is the transformation matrix from the orbital coordinate system to the satellite's own system; It is the projection of the angular velocity of the satellite's own system relative to the orbital system in the own system, indicating the angular velocity of the satellite when performing an additional attitude maneuver sequence around its own axis in addition to the attitude angular motion of the satellite during normal operation; is the projection of the angular velocity of the satellite system relative to the Earth's inertial system in the satellite system; Step 4: Use the gyro error parameter calibration model and satellite angular velocity model to establish the fitness function of quantum genetic optimization; Step 5: Select the chromosome with the highest fitness function value from the current chromosome as the best chromosome of the current generation, and judge whether the fitness difference between the best chromosome of the current generation and the best chromosome of the parent generation is less than the preset threshold, or whether the number of iterative evolutions reaches the preset maximum value. If so, end the evolution process offline and proceed to step 7; otherwise, execute step 6 to achieve global search for the optimal maneuver sequence; Step 6: Use the quantum rotating gate to update the chromosome, and return to step 2 to continue iteratively decoding the next generation of chromosomes; Step 7: The gene sequence of the contemporary best chromosome is the final optimization result of the offline satellite attitude maneuver sequence. It is decoded according to the conversion relationship between quantum bit coding, binary coding, and decimal coding to extract the optimized attitude maneuver sequence; Step 8: During the on-orbit application phase, the satellite executes the optimized attitude maneuver sequence, collects the measurement data of the gyroscope and star sensor, processes the measurement data using the established gyroscope error parameter calibration model, and estimates the final calibration value of the gyroscope error parameter.

2. A method for rapid calibration of gyro error parameters based on intelligent optimization maneuvers as claimed in claim 1, characterized in that: In the step 1, the basic maneuvering forms include six angular velocity magnitude change modes of triangle, trapezoid, sine and their combination, which are used as the basic component units of the attitude maneuvering sequence; The parameters reflecting the attitude maneuver sequence include maneuver form, maximum maneuver angular rate amplitude, single maneuver cycle, single-direction maneuver duration, total maneuver rounds and total maneuver time.

3. A method for rapid calibration of gyro error parameters based on intelligent optimization maneuvers as claimed in claim 1, characterized in that: In the step 1, the number of genes in each chromosome is determined by the number of maneuvers that constitute the attitude maneuver sequence; each chromosome represents the encoding of an attitude maneuver sequence; the attitude maneuver parameters involved in the encoding and optimization include the maneuver type, the axis of the current maneuver, and the maneuver angular velocity amplitude.

4. A method for rapid calibration of gyro error parameters based on intelligent optimization maneuvers as claimed in claim 1, characterized in that: In step 3, for a chromosome with a gene length of m, the encoding method of the kth gene "abcd" is as follows The thousands digit a is the code of the basic maneuver form, 1≤a≤n m , n m is the number of candidate maneuvers; The hundreds digit b is the maximum maneuvering angular velocity amplitude η max The specific calculation formula is as follows: in, η max and are the maximum maneuvering angular velocity amplitude η max The lower and upper limits of the value range; In the decoding stage, the maximum maneuvering angular rate amplitude η is obtained by decoding according to the following formula max : The tens digit c represents the maneuvering axis: 1 corresponds to the x-axis or pitch axis, 2 corresponds to the y-axis or yaw axis, 3 corresponds to the z-axis or roll axis, 1≤c≤3; The last digit d is a flag code indicating whether a maneuver is to be performed in the current period. If d is an odd number, it indicates that an attitude maneuver is performed and the maneuver flag is set to 1; if d is an even number, it indicates that there is no attitude maneuver and the maneuver flag is set to 0.

5. A method for rapid calibration of gyro error parameters based on intelligent optimization maneuvers as claimed in claim 1, characterized in that: In step 2, the gyro error parameters include constant drift b, installation error Δ, scale factor error δs and random measurement noise w g .

6. A method for rapid calibration of gyro error parameters based on intelligent optimization maneuvers as claimed in claim 1, characterized in that: In step 4, the fitness function is as follows: Where, trace(·) represents the matrix trace operation; H is the measurement matrix; t m is the total time of the satellite maneuver process; θ m is the total rotation angle of the satellite during maneuver; η m is the maximum angular velocity amplitude of the satellite during maneuvering; is the time for the satellite to operate at the maximum maneuvering angular velocity amplitude; w o 、w t 、w θ and w η It is the weight coefficient that reflects the contribution of observability level, maneuvering time, maneuvering angle and maneuvering angular rate amplitude in the fitness function.

7. A method for rapid calibration of gyro error parameters based on intelligent optimization maneuvers as claimed in claim 1, characterized in that: In the step 5, after the chromosome population before iterative optimization is manually initialized, the fitness value of each chromosome in the initial population is calculated to determine the initial value of the best chromosome of the parent generation.

8. A method for rapid calibration of gyro error parameters based on intelligent optimization maneuvers as claimed in claim 1, characterized in that: The specific updating process of step 6 is as follows: Taking the best individual in the parent chromosome as the evolutionary goal, using the quantum revolving door Update the population, For nth ga The revolving door angle of the i-th chromosome individual in the generation; The process of updating a quantum bit using a quantum rotating gate is expressed as: in, and The nth ga The two probability amplitudes corresponding to the j-th gene position of the i-th chromosome in the generation evolution; Taking into account the convergence speed and convergence accuracy, the following strategy is used to adaptively adjust the revolving door angle: Among them, Δθ max and Δθ min The maximum and minimum rotation angles, N max is the maximum number of iterations, For nth ga The cosine similarity between the i-th chromosome of the generation and the best chromosome of the parent generation, w c is the weight coefficient of cosine similarity; The quantum NOT gate mutation process is expressed as: in, and are the quantum bits before and after the quantum NOT gate mutation; U not is a quantum NOT gate,

Images