CN105371851B

CN105371851B - A kind of satellite attitude model construction method based on frequency-domain analysis

Info

Publication number: CN105371851B
Application number: CN201510798994.0A
Authority: CN
Inventors: 谢俊峰; 唐新明; 莫凡; 付兴科; 范大昭; 高小明; 唐洪钊; 祝小勇; 窦显辉; 申二华
Original assignee: SATELLITE SURVEYING AND MAPPING APPLICATION CENTER NASG
Current assignee: Ministry Of Natural Resources Land Satellite Remote Sensing Application Center
Priority date: 2015-11-20
Filing date: 2015-11-20
Publication date: 2018-08-07
Anticipated expiration: 2035-11-20
Also published as: CN105371851A

Abstract

The invention discloses one kind taking attitude of satellite variational regularity and periodic satellite attitude model construction method into account, and this method carries out reasonably building satellite attitude model by way of fully taking the intrinsic posture vibration frequency of satellite into account.The attitude mode of structure includes two parts：General polynomial part and trigonometric polynomial part, cover the smooth in satellite attitude signal in general polynomial part, trigonometric polynomial part includes periodicity and regular part in satellite attitude signal, finally realizes the foundation for optimizing satellite attitude model with joint mapping by the above two-part independent fitting.This method is under the premise of fully taking the peculiar rule of satellite itself into account, has the characteristics that high-precision attitude fitting, the interpolation of the attitude of satellite can be carried out using the model, and can effectively inhibit the exceptional value in attitude data, and the white noise to mixing in posture can also effectively inhibit.

Description

Satellite attitude model construction method based on frequency domain analysis

Technical Field

The invention relates to the technical field of satellite attitude post-processing, in particular to a satellite attitude model construction method based on frequency domain analysis, which is applied to ground positioning of satellite photogrammetry.

Background

At present, with the continuous improvement of the ground resolution of a high-resolution remote sensing satellite, the satellite attitude precision becomes one of the main bottlenecks which restrict the development of satellite mapping towards large-scale mapping. Due to the reasons of energy distribution of loads on the satellite, limitation of hardware equipment of an attitude measurement system, bandwidth of a data downloading channel, capacity limitation of ground storage equipment and the like, attitude data cannot be measured and recorded line by line according to images, and the attitude data is acquired by adopting a fixed frequency sampling method. The problem of acquisition and transmission and the like can be well solved through equidistant attitude sampling, for discrete attitude data obtained through rarefaction, discrete numerical values need to adopt an interpolation method to obtain data among discrete points during positioning mapping calculation, so that distortion of the attitude data is inevitably caused, namely errors are brought when ground geometric positioning calculation is carried out by using the interpolated attitude data, and therefore the true mapping attitude of the satellite needs to be approximated through construction of a proper attitude model.

The conventional gesture obtaining method generally includes two types: lagrange interpolation [ Gonghui ] research on high-resolution satellite remote sensing image positioning theory and method based on quaternion [ D ]. Zheng ]: university of information engineering, 2011] [ CLIVEFRASER, JULIANG SHAO. external organization Determination of MOMS-02 Three-LineImagey: experiences with Australian Testfield Data [ J ]. International Archivesof Photogrammetry and Remote Sensing, 1996, 31 (3): 207 + 214 ] and general polynomial fitting [ Yuan Xiaxian, Cao Jinshan. high-resolution satellite remote sensing accurate target location theory and method [ M ], Beijing: scientific press, 2012: 98-110.] [ Dongsook Shin, John K. Pollard, Jan-Peter Muller. "[ Accurategeometric correction of ATSR images,". IEEE TRANSACTIONS 0N GEOSCIENCE ANDREMOTE SENSING, VOL.35, N0.4, July 1997 ]. There are also more mathematically applied interpolation methods such as: a continuous fractional interpolation method, a unitary whole-interval interpolation, a unitary three-point interpolation, an aitin step-by-step difference value, a smooth interpolation and the like [ van city, wangmi, zhuingying and the like. 15-21, 30.].

In the process of satellite in-orbit flight, due to the comprehensive influences of the rotation of the on-satellite solar wing sailboard, the inherent frequency of the momentum wheel and the like, the satellite attitude shows high-frequency periodic change, the attitude change quantity is very small, and the change quantity cannot be reflected by a general polynomial fitting method when the overall fitting of the attitude value is carried out. For the lagrange interpolation method, although the attitude values around the interpolation point are considered in the high-order interpolation, the calculated amount is large, the high-order interpolation is easily influenced by the abnormal value, and the periodic term in the attitude cannot be reflected at the same time; the advantages and disadvantages of the current pose interpolation method are shown in table 1.

TABLE 1 common pose interpolation algorithm comparison

Aiming at the characteristics of the traditional attitude interpolation method, the invention provides a novel attitude acquisition method, overcomes the defects of poor stability, low accuracy of general polynomial and poor interpolation adaptability of a mathematical method in a Lagrange interpolation method, fully considers the characteristics of a satellite platform and constructs a function model of the satellite attitude.

Disclosure of Invention

To this end, the present invention is directed to a method for constructing a satellite attitude model based on frequency domain analysis that substantially obviates one or more of the problems due to limitations and disadvantages of the related art.

Additional advantages, objects, and features of the invention will be set forth in part in the description which follows and in part will become apparent to those having ordinary skill in the art upon examination of the following or may be learned from practice of the invention. The objects and advantages of the invention may be realized and attained by the structure particularly pointed out in the written description and claims hereof as well as the appended drawings.

The invention provides a satellite attitude model construction method based on frequency domain analysis, wherein the method comprises the following steps:

step 1, acquiring satellite attitude data expressed in a normalized quaternion form;

step 2, performing frequency domain analysis on the satellite attitude data obtained in the step 1, wherein the step 2 specifically comprises the following substeps:

step 2.1, converting the expression form of the satellite attitude data from the form of normalized quaternion into the form of Euler angle according to the following formula;

wherein q is₁、q₂And q is₃Representing the imaginary part of a normalized quaternion, q₄Representing the real part of a normalized quaternion, q₄Is calculated by the formulaRoll, Pitch and Yaw are respectively Roll, Pitch and Yaw angles expressed in the form of Euler angles, A_eularRepresenting the transformed three-axis euler angle matrix.

Step 2.2, the satellite attitude data in the Euler angle form obtained in the step 2.1 is converted from an inertial coordinate system to an orbit coordinate system to obtain an attitude value angle matrix A under the orbit coordinate system_{body_orbit}＝[Roll_oriPitch_oriYaw_ori]^T(ii) a Wherein, Roll_ori、Pitch_oriAnd Yaw_oriRespectively the roll angle, pitch angle and yaw angle after the processing of steps 2.1 and 2.2.

Step 2.3, fitting the satellite attitude data by using a general polynomial to obtain attitude data containing frequency and corresponding amplitude; step 2.3 specifically comprises the following steps:

step 2.3.1, determining a fitting formula of a general polynomial:

f(x)＝a₀+a₁x+L+a_nxⁿ(1)

wherein f (x) represents the fitting result, a₀～a_nFor the fitting coefficient, n is the fitting order, and preferably, n is 8:

step 2.3.2, solving polynomial coefficients by using least square estimation method Andwherein, andrespectively representing coefficients of general polynomials corresponding to three angles of Roll, Pitch and Yaw;

step 2.3.3, polynomial coefficients obtained in step 2.3.2 Andsubstituting equation (1) can obtain fitting equation (2):

wherein t represents a time variable, f_Roll(t)、f_Pitch(t) and f_Yaw(t) each represents Roll_ort、 Pitch_ortAnd Yaw_oriFitting the attitude value;

step 2.3.4, on f obtained by step 2.3.3_Roll(t)、f_Pitch(t) and f_Yaw(t) performing discrete processing to obtain a fitted discrete attitude value:

wherein,representing the fitted discrete pose values,andrespectively representing fitted discrete attitude values of a roll angle, a pitch angle and a yaw angle;

step 2.3.5, obtaining the attitude value A through the step 2.2_{body_orbit}And the fitted discrete attitude value obtained by the step 2.3.4A difference calculation was performed as shown below:

in the formula,a difference value representing satellite attitude data, the difference value being attitude data comprising a frequency and a corresponding amplitude;

step 2.4, adopting fast Fourier transform pairPerforming frequency domain analysis to detect a frequency region of the satellite attitude data;

step 3, constructing a satellite attitude model, wherein the step 3 specifically comprises the following substeps:

step 3.1, determining a frequency reservation region;

step 3.2, constructing a triangular polynomial in the satellite attitude model according to the frequency reservation region determined in the step 3.1, as shown in the following formula:

in the formula, g_Roll(t)、g_Pitch(t) and g_Yaw(t) is a trigonometric polynomial relating to the time variable t constructed for Roll, Pitch and Yaw respectively,andthe dc variables are all dc variables corresponding to Roll, Pitch and Yaw, K is a frequency reserving region, K is a frequency corresponding to a K frequency region (the first frequency region K is 1 corresponding to 27-term trigonometric function, the second frequency region K is 2 corresponding to 38-term trigonometric function),andrespectively represents the Fourier coefficients, psi, of Roll, Pitch and Yaw in the K frequency region^Roll、ψ^RollAnd psi^RollRespectively, the initial phases of Roll, Pitch and Yaw in the K frequency region are shown.

Step 3.3, constructing a satellite attitude model mixed by a general polynomial and a triangular polynomial, wherein the constructed satellite attitude model is shown as the following formula:

wherein H_Roll(t)、H_pitch(t)、H_Yaw(t) pose models constructed for Roll, Pitch, and Yaw, respectively, f_Roll(t)、f_pitch(t)、f_Yaw(t) general polynomial model established for step 2.3, g_Roll(t)、g_Pitch(t)、g_Yaw(t) is the trigonometric polynomial model established in step 3.2.

Preferably, the step 1 specifically includes: randomly extracting data downloaded from an orbit satellite to obtain a group of discrete satellite attitude data expressed in a normalized quaternion form.

Preferably, the step 2.2 specifically includes:

step 2.2.1, coordinate conversion is carried out according to the formula (3) to obtain a rotation matrix

Wherein,is A_eularThe rotation matrix of (a) is,is a rotation matrix from an inertial coordinate system to an orbital coordinate system,a rotation matrix of the body coordinate system relative to the track coordinate system;

wherein, X_s、Y_s、Z_sThe position of the satellite centroids in the spatially-fixed inertial reference frame is determined, the velocity of the satellite centroid in the spatially-fixed inertial reference frame is determined;

(X₂)_X、(X₂)_Yand (X)₂)_ZRepresenting the vector soughtThe components in the X, Y and Z directions respectively,

(Y₂)_X、(Y₂)_Yand (Y)₂)_ZRepresenting the vector soughtThe components in the X, Y and Z directions respectively,

(Z₂)_X、(Z₂)_Yand (Z)₂)_ZRepresenting the vector soughtThe components in the X, Y and Z directions, respectively.

Step 2.2.2, rotating the matrix according to equation (4)Converted into corresponding angle matrix A_{body_orbit}：

Preferably, the step 2.3.2 specifically includes:

the formula of the least squares adjustment is:

AX＝B (5)

from equation (5):

X＝(A^TA)^-1A^TB (6)

the values of X and B are respectively substituted into formula (6) to obtain the polynomial coefficient And

preferably, in step 2.3.4, the pair f is sampled_Roll(t)、f_Pitch(t) and f_Yaw(t) discrete processing is performed and the sampling time interval is the same as the sampling time interval inherent to the satellite.

Preferably, in step 3.1, the intersection of the frequency region detected in step 2.4 with the frequency region inherent to the satellite is retained, and the discarded other frequencies are truncated and removed.

Preferably, the real and imaginary parts of the complex numbers of the discarded other frequency parts are zeroed out.

Preferably, in step 2.4, the frequency region of the detected satellite attitude data is a frequency interval, and the frequency reservation region determined in step 3.1 is also the corresponding frequency interval. Preferably, in step 3.2, the different frequency regions K correspond to frequencies K relative to the region.

The invention provides a satellite attitude model construction method considering the regularity and periodicity of satellite attitude change. The constructed attitude model comprises two parts: the general polynomial part covers a smooth part in the satellite attitude signal, the triangular polynomial part comprises a periodic part and a regular part in the satellite attitude signal, and the establishment of the optimized satellite attitude model is finally realized through the independent fitting and the joint construction of the two parts. The method has the characteristic of high-precision attitude fitting on the premise of fully considering the special regularity of the satellite, can utilize the model to perform interpolation of the attitude of the satellite, can effectively inhibit abnormal values in attitude data, and can also effectively inhibit white noise mixed in the attitude. Based on a large amount of experimental verification, the method fully considers more than 95% of effective periodic attitude jitter in attitude signals, so that the triaxial attitude interpolation precision is better than 0.3 arc second, and the constructed attitude model meets the precision requirement of general satellite surveying and mapping.

Drawings

FIG. 1 is a flowchart of a method for constructing a satellite attitude model based on frequency domain analysis according to an embodiment of the present invention;

FIG. 2A is the initial attitude Roll angle Roll of the body system relative to the rail system after pre-processing in Euler angle form;

FIG. 2B is the original attitude Pitch angle Pitch of the pre-processed body system relative to the rail system in Euler angle form;

FIG. 2C is the preprocessed original attitude Yaw angle Yaw of the body system relative to the orbital system in Euler angle form;

FIG. 3A is a diagram of the attitude residuals after the Roll angle Roll is subjected to 8 th order polynomial fitting;

FIG. 3B is a diagram of the attitude residuals after the Pitch angle Pitch is subjected to 8-order polynomial fitting;

FIG. 3C is a diagram of the attitude residuals after the Yaw angle Yaw is subjected to 8-order polynomial fitting;

FIG. 4A is a spectrogram obtained by fast Fourier transforming Roll attitude residuals;

fig. 4B is a frequency spectrum diagram obtained by fast fourier transform of the Pitch angle Pitch attitude residual error;

FIG. 4C is a spectrogram obtained by fast Fourier transforming the Yaw angle Yaw attitude residual error;

FIG. 5A is a diagram of the effect of attitude fitting after cross Roll angle Roll frequency domain screening;

FIG. 5B is a diagram of the attitude fitting effect after the Pitch angle Pitch frequency domain screening;

FIG. 5C is a diagram of the attitude fitting effect after frequency domain screening of the Yaw angle Yaw;

FIG. 6A is a residual error plot of Roll angle Roll build model calculated values versus original attitude values;

FIG. 6B is a residual error plot of the Pitch angle Pitch model calculation value versus the original attitude value;

FIG. 6C is a residual plot of Yaw angle Yaw model-constructed calculations versus the original attitude values.

Detailed Description

The present invention now will be described more fully hereinafter with reference to the accompanying drawings, in which exemplary embodiments of the invention are shown.

As shown in fig. 1, the method for constructing a satellite attitude model based on frequency domain analysis provided by the present invention specifically includes the following steps:

step 1, acquiring satellite attitude data expressed in a normalized quaternion form; the step 1 specifically comprises the following steps: randomly extracting data downloaded from an orbit satellite to obtain a group of discrete satellite attitude data expressed in a normalized quaternion form.

According to the embodiment of the invention, satellite attitude data can be randomly extracted; preferably, the embodiment of the present invention takes the attitude data downloaded from the resource satellite three as an example, but is not limited thereto. Randomly extracting the attitude data of 00381 th orbit, which is the attitude data of the orbit at the initial transmitting stage, from the satellite, wherein the acquisition time is 2, 3 and 2 months in 2012, and the data comprises 2188 discrete attitude data values expressed in the form of quaternions.

2.1, converting the expression form of the satellite attitude data from the normalized quaternion form into an Euler angle form;

the satellite attitude is generally expressed in the form of an Euler angle or a quaternion, and as the quaternion has the advantages of nonsingularity, smooth interpolation, angular displacement inversion and quick connection, most remote sensing satellites describe attitude data in the form of quaternion, the attitude described in the form of quaternion has the defect of poor physical property and intuition, the attitude described in the form of Euler angle directly reflects the variation of the satellite three-axis attitude, and the physical significance is obvious; for ease of analysis, the present invention first unifies the gesture description form into the euler angles.

The 00381 th orbit attitude data downloaded from the resource No. three satellite is expressed in the form of normalized quaternion, which is specifically shown in the table 1:

TABLE 1 Star upload and download attitude data

The sequence number columns in table 1 represent the order of the attitude data, which corresponds to the time of acquisition of the attitude of the satellite, with a time interval of 0.25s between two sequence numbers, i.e. acquisition of the attitude data of the satellite at a time interval of 0.25s, q₁、q₂And q is₃The columns represent satellite-up-and-down recorded attitude values in the form of normalized quaternion, where a set of q₁、q₂And q is₃The value describes an attitude, q, of the satellite₁、q₂And q is₃The imaginary part, the real part q, of the attitude values in quaternion form₄Since it can pass throughAnd obtaining, in order to avoid embellishment residue, a real part is not downloaded.

Converting the expression form of the attitude data from the form of the normalized quaternion into the form of an Euler angle according to a conversion formula (10):

wherein q is₁、q₂And q is₃Denotes the imaginary part of the normalized quaternion, which is also the corresponding q in Table 1₁、q₂And q is₃Column data, q₄The real part of the normalized quaternion is expressed and can be obtained by imaginary part calculation, and the calculation formula isRoll, Pitch and Yaw are respectively Roll, Pitch and Yaw angles expressed in the form of Euler angles, A_eularRepresenting the transformed three-axis euler angle matrix.

Taking a group of data with the sequence number of 0 as an example to show the process of converting the expression form of the attitude data from the form of the normalized quaternion into the form of the Euler angle, q₁＝-0.76259744，q₂0.58472162 and q₃0.25788319, thenSubstituting the formula (10) to obtain the Euler angle matrix A corresponding to the q values_eular＝[0.9625 0.8794 -1.1201]^T。

2.2, converting the satellite attitude data in the Euler angle form obtained in the step 2.1 from an inertial coordinate system to an orbit coordinate system;

the coordinate system of the attitude data downloaded from the resource number three satellite is a J2000 inertial coordinate system, and in order to better reflect the change of the attitude of the satellite platform, the coordinate system needs to be converted into an orbit coordinate system (the orbit coordinate system only reflects the angle of the satellite platform, and does not relate to other angle interference).

The step 2.2 specifically comprises:

step 2.2.1, coordinate conversion is carried out according to the formula (11) to obtain a rotation matrix

Step 2.2.2, rotating the matrix according to equation (13)Converted into corresponding angle matrix A_{body_orbit}。

The rotation matrix is obtained from equation (11)A_{eular_orbit}Is a rotation matrixAngle matrix of (1), wherein angle matrix A and rotation matrixCan be represented by the following relationship:

if angle matrixThen a corresponds to the rotation matrix as shown in equation (12):

rotation matrix obtained by equation (11)Using the angle matrix A and the rotation matrixThe rotation matrix can be obtained by converting the relationship betweenConverted into corresponding angle matrix A_{body_orbit}The specific calculation is shown in formula (13):

in equation (13), Roll, Pitch, and Yaw are Roll angle, Pitch angle, and Yaw angle expressed in the form of euler angles, respectively.

Equation (13) is a rotation matrix obtained by equation (11)And then, the rotation matrix and the angle matrix are converted according to the relation of the formula (12).

After the attitude data is preprocessed (the expression form of the attitude data is converted into an Euler angle form from a normalized quaternion form, and the attitude data is converted from a body coordinate system relative to an inertial coordinate system to the body coordinate system relative to an orbit coordinate system), an attitude value angle matrix A in the orbit coordinate system described by the Euler angle can be obtained_{body_orbit}＝[Roll_oriPitch_oriYaw_ori]^T. Similarly, taking the data with serial number 0 as an example, the process of transformation of coordinate system is explained, and step 2.1 obtains A_eular＝[0.9625 0.8794 -1.1201]^TAnd calculating the rotation matrix by combining the rotation matrix from the inertial system to the orbital system given by the international geodetic society through the calculation of the formula (11) and the formula (13) to obtain A_{eular_orbit}＝[-0.00001245 -0.00000509 0.03356850]^T) Calculating all data to obtain A_{eular_orbit}The attitude data is shown in table 2, and the attitude overall change is shown in fig. 2A, 2B, and 2C, respectively.

TABLE 2 post-preprocessing attitude data

In Table 2, the sequence number is given as the order of the preprocessed pose data, Roll_ori、Pitch_oriAnd Yaw_oriThe columns are the values (in radians) of roll angle, pitch angle and yaw angle, respectively, after preprocessing, corresponding to fig. 2A, 2B and 2C, respectively.

And 2.3, fitting the satellite attitude data by using a general polynomial to obtain attitude data containing frequency and corresponding amplitude. Step 2.3 specifically comprises the following steps:

step 2.3.1, determining a fitting formula of a general polynomial:

f(x)＝a₀+a₁x+L+a_nxⁿ(14)

wherein f (x) represents the fitting result, a₀～a_nFor fitting coefficient, n is fitting order, and 8-order general polynomials are preferably adopted for Roll according to experience_ori、Pitch_oriAnd Yaw_oriThe fitting is performed, i.e. in an embodiment of the invention, preferably n-8. The above formula is a theoretical formula for fitting the satellite attitude data, and according to the formula, the parameters of a general polynomial are firstly required to be solved when fitting is carried out.

Step 2.3.2, solving polynomial coefficients by using least square estimation method Andwherein,andthe coefficients of general polynomials for three angles, Roll, Pitch and Yaw, are shown. Step 2.3.2 specifically includes:

the formula of the least squares adjustment is:

AX＝B (15)

from the least squares principle, we can derive from equation (15):

X＝(A^TA)^-1A^TB (16)

the values of X and B are respectively substituted into the formula (16) to obtain the polynomial coefficient And

step 2.3.3, polynomial coefficients obtained in step 2.3.2 Andinstead of formula (14), fitting formula (17) can be obtained:

in the formula, t represents a time variationAmount f_Roll(t)、f_Pitch(t) and f_Yaw(t) each represents Roll_ori、Pitch_oriAnd Yaw_oriAnd (5) fitting the attitude value. The values obtained by the fitting are shown in table 3:

TABLE 3

In Table 3, the coefficient columns represent polynomial coefficients, column f_Roll(t)_coefficient、f_Pitch(t)_coefficient、 f_Yaw(t) _ coefficient represents a polynomial coefficient value corresponding to the polynomial fit value.

wherein,representing the fitted discrete pose values,andand respectively representing fitted discrete attitude values of a roll angle, a pitch angle and a yaw angle.

In particular, in step 2.3.4, preferably, f is sampled_Roll(t)、f_Pitch(t) and f_Yaw(t) discrete processing is performed and the sampling time interval is the same as the sampling time interval inherent to the satellite.

That is, will obtainGeneral polynomial f_Roll(t)、f_Pitch(t)、f_Yaw(t) performing data discrete processing at the time sampling interval of the original data, namely substituting each fixed time into a formula to calculate an attitude value, wherein the discrete attitude value is usedExpressed, the discretization process can be expressed as:

step 2.3.5, obtaining the attitude value A through the step 2.2_{body_orbit}And the fitted discrete attitude value obtained by the step 2.3.4The difference calculation is performed as shown in the following formula (19):

in the formula,a difference value representing the attitude data of the satellite, the difference value being the attitude data including a frequency and a corresponding amplitude.

The linear part in the original value can be fitted through polynomial fitting, and subtracting the fitted value from the original value is equivalent to offsetting the linear part, so that a model of the linear part is constructed, and a posture part containing certain frequency and corresponding amplitude is left, thereby facilitating further analysis.

As shown in fig. 3A, 3B and 3C, compared to fig. 2A, 2B and 2C, the linear part in the original data is substantially cancelled,the change of the posture along the horizontal axis can be visually seen.

using fast Fourier transform pairsThe three axes respectively carry out frequency domain analysis, the theoretical basis of the fast Fourier transform is Fourier series, the theory of the Fourier series is visual, and the formula for carrying out infinite approximation on data is as follows:

wherein, a_kAnd b_kAre Fourier coefficients, corresponding to both the real and imaginary parts of the Fourier transform, a₀For the dc variable, k is the number of the fourier series (frequency dependent), f (x) is the result of the fit of the fourier series, and x is the data sequence.

Because the theoretical basis of the fast Fourier transform is Fourier series, the fast Fourier transform algorithm is a fast algorithm of discrete Fourier transform, the writing of a formula is complicated, and the fast Fourier transform algorithm belongs to theoretical knowledge and engineering algorithms disclosed in the field of computers.

To pairPerforming three times of fast Fourier transform at different angles respectively to obtainThe Roll part in (1) is an example, namely data shown in fig. 3A, and is subjected to fast fourier transform:

in the formula, m and n respectively represent a real part and an imaginary part obtained after fast Fourier transform, j is a serial number of a complex number group,to representIn the Roll section, the FFT is a fast fourier transform.

By the same method, forThe Pitch and Yaw sections of (1) are fast fourier transformed.

A series of complex number groups [ m + ni ] related to frequency can be obtained through fast Fourier transform]From the properties of the fast fourier transform, it can be seen that the sequence number j has the following correspondence with the frequency: f is j/NT, where f is frequency, j is 1, 2, L, N, T is sampling time interval, N is the number of sampling points, and the complex number and the amplitude have the following correspondence:where a is the amplitude. The conversion from the time domain to the frequency domain is completed through a formula (21), a theoretical basis is provided for the next processing, namely, a frequency region and an amplitude value thereof contained in the satellite attitude data are obtained through the calculation in the step, and in order to more intuitively display and facilitate the processing in the step 3.1, a spectrogram is obtained by using the amplitude and the frequency, as shown in fig. 4A; similarly, the spectrograms of Pitch and Yaw, namely FIG. 4B and FIG. 4C, can be obtained according to the above method, and signals of 0.0-0.05Hz and 0.64-0.71Hz exist in the satellite attitude data can be obtained from FIG. 4A, FIG. 4B and FIG. 4C, namely, the frequency regions of the satellite attitude data are 0.0-0.05Hz and 0.64-0.71Hz, and the natural frequency in the satellite attitude of resource No. three is covered by the coverIncluding both regions.

step 3.1, determining a frequency reservation region;

through a large number of experimental verifications and documents, the inherent attitude jitter frequency of the satellite in-orbit flight comprises five frequency regions of 0-0.22Hz, 0.28Hz, 0.58Hz, 0.6-0.75Hz and 1.12Hz, the five frequency regions are the jitter frequency existing in the satellite platform, the frequency region detected in the step 2.4 is reserved if the frequency region intersects with the inherent frequency region of the satellite (the frequency regions are reserved at 0.0-0.05Hz and 0.64-0.71Hz), because the intersection part is the effective jitter of the satellite, other parts are considered to be random noise caused by unknown reasons, and other frequencies are truncated and removed (the real part and the imaginary part of the complex number of the frequency part to be discarded are cleared).

Step 3.2, constructing a triangular polynomial in the satellite attitude model according to the frequency reservation regions (0.0-0.05Hz and 0.64-0.71Hz) determined in the step 3.1, as shown in the following formula (22):

in the formula, g_Roll(t)、g_Pitch(t) and g_Yaw(t) is a trigonometric polynomial relating to the time variable t constructed for Roll, Pitch and Yaw respectively,andthe direct current variables correspond to Roll, Pitch and Yaw, K is a frequency reserving region, when K is 1, the frequency reserving region is 0.0-0.05Hz, when K is 2, the frequency reserving region is 0.64-0.71Hz, and K is the frequency corresponding to the K frequency region (the first frequency region K is 1 corresponds to 27-term trigonometric function formula, and the second frequency region K is 2 corresponds to 38-term trigonometric function formula), andrespectively represents the Fourier coefficients, psi, of Roll, Pitch and Yaw in the K frequency region^Roll、ψ^RollAnd psi^RollRespectively, the initial phases of Roll, Pitch and Yaw in the K frequency region are shown.

The coefficients of the trigonometric polynomial (in Roll for example) are shown in table 4:

TABLE 4 trigonometric polynomial coefficients

In the table, the sequence number column indicates the sequence number of the trigonometric function,the columns represent the amplitude of the trigonometric function, k corresponds to the frequency,. phi^RollIndicating the phase.

3.3, constructing a satellite attitude model mixed by a general polynomial and a triangular polynomial;

constructing a trigonometric polynomial based on the relationship between the fourier series and the trigonometric function according to the frequency region provided in the step 3.1 and the amplitude and the phase of the frequency region provided in the step 2.4, namely constructing a trigonometric function formula by using the corresponding amplitudes and phases of the frequency domains of 0.0-0.05Hz and 0.64-0.71 Hz; using a general polynomial f_Roll(t)、f_Pitch(t)、f_Yaw(t) and trigonometric polynomial g_Roll(t)、g_Pitch(t)、g_Yaw(t) the hybrid linear model constructed as shown in equation (23):

in the formula, H_Roll(t)、H_Pitch(t)、H_Yaw(t) pose models constructed for Roll, Pitch, and Yaw, respectively, f_Roll(t)、f_Pitch(t)、f_Yaw(t) general polynomial model established for step 2.3, g_Roll(t)、 g_Pitch(t)、g_Yaw(t) is the trigonometric polynomial model established in step 3.2, andthe coefficients of the polynomial expression are represented by,the dc variables of Roll, Pitch and Yaw, K is a frequency reserving region, K is 1 corresponding to the frequency domain of 0.0-0.05Hz, K is 2 corresponding to the frequency domain of 0.64-0.71Hz, K is the frequency corresponding to the K frequency region (the first frequency region K is 1 corresponding to 27-term trigonometric function formula, the second frequency region K is 2 corresponding to 38-term trigonometric function formula),andrespectively represents the Fourier coefficients, psi, of Roll, Pitch and Yaw in the K frequency region^Roll、ψ^RollAnd psi^RollRespectively, the initial phases of Roll, Pitch and Yaw in the K frequency region are shown.

The constructed models are shown in fig. 5A, 5B, and 5C, respectively, and the relative residuals of the models (the difference between the values calculated using the constructed models and the preprocessed values) are shown in fig. 6A, 6B, and 6C, respectively.

The attitude model is constructed by combining the advantages of the traditional general polynomial for constructing the attitude, and the triangular polynomial is introduced to utilize regular parts which are difficult to express by the general polynomial, so that the precision of the traditional general polynomial for constructing the satellite attitude is improved, and meanwhile, the constructed model has a certain extrapolation effect and can meet the mapping requirement with higher precision.

The above description is only a preferred embodiment of the present invention, and for those skilled in the art, the present invention should not be limited by the description of the present invention, which should be interpreted as a limitation.

Claims

1. A satellite attitude model construction method based on frequency domain analysis is characterized by comprising the following steps:

wherein q is₁、q₂And q is₃Representing the imaginary part of a normalized quaternion, q₄Representing the real part of a normalized quaternion, q₄Is calculated by the formulaRoll, Pitch and Yaw are respectively Roll, Pitch and Yaw angles expressed in the form of Euler angles, A_eularRepresenting the converted triaxial Euler angle matrix;

step 2.2, the satellite attitude data in the Euler angle form obtained in the step 2.1 is converted from an inertial coordinate system to an orbit coordinate system to obtain an attitude value angle matrix A under the orbit coordinate system_{body_orbit}＝[Roll_oriPitch_oriYaw_ori]^T(ii) a Wherein, Roll_ori、Pitch_oriAnd Yaw_oriRespectively the roll angle, the pitch angle and the yaw angle processed in the steps 2.1 and 2.2;

step 2.3.1, determining a fitting formula of a general polynomial:

f(x)＝a₀+a₁x+…+a_nxⁿ(1)

wherein f (x) represents the fitting result, a₀～a_nN is a fitting order, wherein n is 8;

step 2.3.2, solving polynomial coefficients by using least square estimation methodAndwherein,andrespectively representing coefficients of general polynomials corresponding to three angles of Roll, Pitch and Yaw;

step 2.3.3, polynomial coefficients obtained in step 2.3.2Andsubstituting equation (1) can obtain fitting equation (2):

wherein t represents a time variable, f_Roll(t)、f_Pitch(t) and f_Yaw(t) each represents Roll_ori、Pitch_oriAnd Yaw_oriFitting the attitude value;

step 3.1, determining a frequency reservation region;

in the formula, g_Roll(t)、g_Pitch(t) and g_Yaw(t) is a trigonometric polynomial relating to the time variable t constructed for Roll, Pitch and Yaw, respectively, a^Roll、a^PitchAnd a^YawAre direct current variables corresponding to Roll, Pitch and Yaw respectively, K is a frequency reservation region, K is a frequency corresponding to a K frequency region,andrespectively represents the Fourier coefficients, psi, of Roll, Pitch and Yaw in the K frequency region^Roll、ψ^RollAnd psi^RollRespectively representing initial phases of Roll, Pitch and Yaw corresponding to a K frequency region;

step 3.3, constructing a satellite attitude model mixed by general polynomial and trigonometric polynomial,

the constructed satellite attitude model is shown as the following formula:

2. The method according to claim 1, wherein step 1 specifically comprises: randomly extracting data downloaded from an orbit satellite to obtain a group of discrete satellite attitude data expressed in a normalized quaternion form.

3. The method according to claim 1, characterized in that said step 2.2 comprises in particular:

in the formula (3), the reaction mixture is,

wherein (X)₂)_X、(X₂)_YAnd (X)₂)_ZRepresenting the vector soughtComponents in the three directions of X, Y and Z, respectively, (Y)₂)_X、(Y₂)_YAnd (Y)₂)_ZRepresenting the vector soughtComponents in the three directions X, Y and Z, respectively, (Z)₂)_X、(Z₂)_YAnd (Z)₂)_ZRepresenting the vector soughtThe components in the X, Y and Z directions, respectively;

wherein,

wherein, X_s、Y_s、Z_sThe position of the satellite centroids in the spatially-fixed inertial reference frame is determined,the velocity of the satellite centroid in the spatially-fixed inertial reference frame is determined;

4. Method according to claim 1, characterized in that in step 2.3.4, f is sampled by_Roll(t)、f_Pitch(t) and f_Yaw(t) discrete processing is performed and the sampling time interval is the same as the sampling time interval inherent to the satellite.

5. Method according to claim 1, characterized in that in step 3.1, the intersection of the frequency region detected in step 2.4 with the frequency region inherent to the satellite is retained, while the other frequencies discarded are truncated.

6. A method as claimed in claim 5, characterized in that the real and imaginary parts of the complex numbers of the discarded other frequency parts are zeroed out.

7. A method according to claim 1, characterized in that in step 2.4 the frequency region of the detected satellite attitude data is a frequency interval and the frequency reservation region determined in step 3.1 is also the corresponding frequency interval.

8. The method according to any of claims 1-7, characterized in that in step 3.2, different frequency regions K correspond to frequencies K relative to this region.