A high-precision lunar angle prediction method for space remote sensing instruments

Abstract

A high-precision lunar angle prediction method for a space remote sensing instrument includes the following steps: according to the forecast start time broadcast by a satellite platform and the required forecast time interval, the lunar geocentric apparent right ascension and apparent declination at the current forecast time are calculated; according to the satellite orbit instant root and perturbation motion equation at the forecast start time broadcast by the satellite platform, the satellite orbit instant root at the current forecast time is calculated; according to the satellite orbit instant root at the current forecast time, the transformation matrix from the geocentric inertial coordinate system to the satellite orbit coordinate system is calculated; according to the satellite attitude angle at the forecast start time broadcast by the satellite platform, the satellite attitude correction matrix at the current forecast time is calculated; according to the above results and the calibrated satellite body coordinate system to the instrument body coordinate system conversion matrix, the lunar angle in the instrument body coordinate system at the current forecast time is calculated. The present invention simplifies the forecast process of the lunar angle, is simple to calculate, and can perform long-term and high-precision lunar angle forecasts.

Description

Technical Field

[0001] The invention relates to the technical field of space remote sensing instrument control, and in particular to a high-precision moon angle prediction method for space remote sensing instruments. Background Art

[0002] Space remote sensing instruments are widely used in various fields, and their accurately measured remote sensing data can provide scientific basis for many fields. Since space remote sensing instruments face long-term attenuation problems on orbit, in order to ensure the accuracy of remote sensing data, space remote sensing instruments will design calibration methods for calibration operations. The results of on-orbit calibration can be used to monitor the long-term performance changes of space remote sensing instruments and directly reflect the working status of space remote sensing instruments on orbit. At the same time, the calibration results on-orbit are also the basis for correcting remote sensing data. According to the different calibration light sources, on-orbit calibration can be divided into on-orbit internal calibration system and on-orbit external calibration system. The light source of the on-orbit internal calibration system is the built-in light source of the space remote sensing instrument; while the light source of the on-orbit external calibration system is the sun, the moon, stars, etc. Among them, the reflection characteristics of the moon are stable and the brightness is low. When the space remote sensing instrument observes the moon, the incident light avoids the influence of scattering and absorption of the atmosphere. Therefore, the moon is suitable as a low-brightness radiation calibration reference source.

[0003] Since the satellite is affected by various perturbations during its in-orbit operation, the instantaneous root of the satellite orbit will change over time. At the same time, due to the movement of the moon itself, its incident angle will also change with the seasons. As a commonly used calibration method for space payloads, lunar calibration usually requires the satellite's attitude to be adjusted to ensure that the lunar angle meets the requirements during calibration. Therefore, in order to ensure the smooth progress of the calibration and the accuracy of the results, it is necessary to master the timing and sustainable duration of the lunar calibration in advance, that is, to predict the lunar angle. At present, space remote sensing instruments usually use semi-analytical theory to calculate the lunar angle. The theory contains a set of calculation methods and corresponding iterative periodic terms. However, due to the large number of these iterative periodic terms, the semi-analytical theory calculation process becomes complicated. At the same time, some iterative periodic terms are very small values, and the gain to the result is limited. In addition, due to the limited computing resources of the processor that meets the radiation resistance index, it is impossible to perform high-precision and complex calculations. Therefore, space remote sensing instruments can only predict the lunar angle with low accuracy in a short time. Summary of the invention

[0004] The purpose of the present invention is to provide a high-precision lunar angle prediction method for space remote sensing instruments to solve the problems of complex calculations and only short-term and low-precision lunar angle prediction in the existing semi-analytical theory for lunar angle prediction.

[0005] The technical solution adopted by the present invention to solve the technical problem is as follows:

[0006] A high-precision lunar angle prediction method for a space remote sensing instrument of the present invention comprises the following steps:

[0007] Step S1: Calculate the lunar geocentric apparent right ascension and apparent declination at the current forecast time according to the forecast start time broadcast by the satellite platform and the required forecast time interval;

[0008] Step S2: Calculate the instantaneous root of the satellite orbit at the current forecast time according to the instantaneous root of the satellite orbit at the forecast start time broadcast by the satellite platform and the pre-established perturbation dynamic motion equation;

[0009] Step S3: Calculate the transformation matrix from the geocentric inertial coordinate system to the satellite orbit coordinate system at the current forecast time according to the instantaneous root of the satellite orbit at the current forecast time;

[0010] Step S4: Calculate the satellite attitude correction matrix at the current forecast time according to the satellite attitude angle at the forecast start time broadcast by the satellite platform;

[0011] Step S5: Calculate the moon angle in the instrument body coordinate system at the current forecast time based on the above results and the calibrated conversion matrix from the satellite body coordinate system to the instrument body coordinate system.

[0012] Further, in step S1, according to the forecast start time t0 broadcast by the satellite platform and the required forecast time interval Δt, the lunar geocentric apparent right ascension and apparent declination at the current forecast time t0+m*Δt are calculated, where m is the number of forecast time intervals Δt, i.e., compared with the forecast start time t0, the current forecast time t0+m*Δt has increased m forecast time intervals Δt; m starts from 0, and after the lunar angle forecast at the current moment is completed, m increases by 1 and is used for the lunar angle forecast at the next moment.

[0013] Furthermore, the specific operation process of step S1 is as follows:

[0014] Step S1.1: Calculate the Julian day J according to the current forecast time t0+m*Δt D , then the Julian century number T is:

[0015] T=(J D -2451545) / 36525 (1)

[0016] Step S1.2: Calculate the lunar geocentric ecliptic longitude period term and ecliptic latitude period term at the current forecast time t0+m*Δt;

[0017] Step S1.2.1: The expressions of the lunar geocentric ecliptic longitude and latitude in the ecliptic coordinate system are:

[0018] λ=A×sin(θ) (2)

[0019] β=A×sin(θ) (3)

[0020] Where, λ is the lunar geocentric ecliptic longitude, β is the lunar geocentric ecliptic latitude, A is the amplitude, and θ is the argument;

[0021] The calculation formula of the argument θ is:

[0022] θ=a×D+b×M+c×M'+d×F (4)

[0023] Where a, b, c, d are all argument coefficients, all in degrees; D is the sun-moon distance angle, M is the sun's mean near-earth angle, M' is the moon's mean near-earth angle, and F is the moon's ascending node's mean angular distance;

[0024] Step S1.2.2: Correct the lunar geocentric ecliptic period term;

[0025] If the angle θ contains M or -M, each lunar geocentric ecliptic period term is multiplied by the coefficient E1 for correction; if the angle θ contains 2M or -2M, each lunar geocentric ecliptic period term is multiplied by the square of the coefficient E1 for correction;

[0026] The calculation formula of the coefficient E1 is:

[0027] E1=1-0.002516×T-0.0000074×T 2 (9)

[0028] The sum of the lunar geocentric ecliptic period terms ∑ I The calculation formula is:

[0029] ∑I=I1+I2+...+I 60 (10)

[0030] In the formula, I1, I2, ... I 60 They are the first to the sixtieth geocentric ecliptic cycle items of the moon respectively;

[0031] Step S1.2.3: Correct the lunar geocentric apparent ecliptic latitude period term;

[0032] The correction method of the lunar geocentric ecliptic latitude period term is the same as step S1.2.2; the calculation method of the sum of lunar geocentric ecliptic latitude period terms ∑b is the same as the calculation method of the sum of lunar geocentric ecliptic longitude period terms ∑I;

[0033] Step S1.3: Venus perturbation correction;

[0034] Step S1.4: nutation correction;

[0035] Step S1.5: According to the above results, the lunar geocentric ecliptic longitude and ecliptic latitude at the current forecast time t0+m*Δt are obtained:

[0036] λ=L'+∑I / 1000000+Δψ (21)

[0037] β=∑b / 1000000 (22)

[0038] Where, λ is the lunar geocentric ecliptic longitude, L' is the lunar mean ecliptic longitude, ∑I is the sum of the lunar geocentric ecliptic longitude periodic terms, Δψ is the nutation of the ecliptic longitude, β is the lunar geocentric ecliptic latitude, ∑b is the sum of the lunar geocentric ecliptic latitude periodic terms;

[0039] Step S1.6: Calculate the lunar geocentric apparent right ascension and declination at the current forecast time t0+m*Δt;

[0040] The calculation formulas for the lunar geocentric apparent right ascension and declination are:

[0041]

[0042] Where, α is the lunar geocentric right ascension, δ is the lunar geocentric declination, λ is the lunar geocentric ecliptic longitude, β is the lunar geocentric ecliptic latitude, ε is the true ecliptic obliquity, which is equal to the sum of the mean ecliptic obliquity ε0 and the nutation Δε.

[0043] The calculation formula of the mean obliquity angle ε0 is:

[0044]

[0045] According to the above results, the calculation formula of the lunar vector M1 in the geocentric inertial coordinate system is:

[0046]

[0047] Furthermore, in step S1.2.1, the calculation formulas for the sun-moon distance angle D, the sun's mean near-earth angle M, the moon's mean near-earth angle M', and the moon's ascending node's mean angular distance F are respectively:

[0048]

[0049]

[0050]

[0051]

[0052] Furthermore, the specific operation process of step S1.3 is as follows:

[0053] ΔI=3958×sin(A1)+1962×sin(L'-F)+318×sin(A2) (11)

[0054]

[0055] Wherein, ΔI is the correction term for the sum of the lunar geocentric ecliptic longitude period terms. When making corrections, the sum of the lunar geocentric ecliptic longitude period terms ∑I is added to ΔI; Δb is the correction term for the sum of the lunar geocentric ecliptic latitude period terms. When making corrections, the sum of the lunar geocentric ecliptic latitude period terms ∑b is added to Δb; L' is the lunar mean ecliptic longitude, and its calculation formula is formula (13); A1, A2 and A3 all represent the perturbation angle corrections, and their calculation formulas are (14), (15) and (16) respectively:

[0056]

[0057] A1=119.75+131.849×T (14)

[0058] A2=53.09+479264.29×T (15)

[0059] A3=313.45+481266.484×T (16).

[0060] Furthermore, the specific operation process of step S1.4 is as follows:

[0061] Step S1.4.1: Calculation of ecliptic nutation;

[0062] Each ecliptic nutation period term S includes the calculation of the sine coefficient of the ecliptic nutation and the linear combination coefficient of the five basic angular distances for calculating the nutation amplitude angle; the calculation formula for the nutation amplitude angle θ1 is:

[0063] θ1=a1×D+b1×M+c1×M’+d1×F+e1×Ω1 (17)

[0064] Where a1, b1, c1, d1, and e1 are all nutation angular coefficients, and Ω1 is the ecliptic longitude of the ascending node of the ecliptic and the lunar orbit. The calculation formula is:

[0065]

[0066] The calculation formula of the ecliptic nutation period term S is:

[0067] S=(S1+S2×T)×sin(θ1) (19)

[0068] Where θ1 is the nutation amplitude angle, S1 and S2 are the sinusoidal coefficients of the ecliptic nutation;

[0069] The ecliptic nutation period term S is accumulated to obtain the ecliptic nutation Δψ, whose unit is 0.0001';

[0070] Step S1.4.2: Calculation of the nutation angle;

[0071] Each cross angle nutation period term C includes the cosine coefficient of the cross angle nutation and the linear combination coefficient of the five basic angular distances of the nutation amplitude; the calculation formula of the cross angle nutation period term C is:

[0072] C=(C1+C2×T)×cos(θ1) (20)

[0073] Where θ1 is the nutation amplitude angle, and C1 and C2 are the cosine coefficients of the nutation angle.

[0074] Furthermore, the specific operation process of step S2 is as follows:

[0075] Introducing the auxiliary quantities eccentric anomaly angle and flat anomaly angle, the relationship between the eccentric anomaly angle, flat anomaly angle and true anomaly angle is:

[0076]

[0077] M3=Ee×sinE (27)

[0078]

[0079] Where, f is the true anomaly, e is the eccentricity, E is the eccentric anomaly, and M3 is the mean anomaly;

[0080] Assuming that the perturbation function corresponding to the satellite in orbit is R, the Lagrangian perturbation equation is:

[0081]

[0082] Where a2 is the semi-major axis, e is the eccentricity, M3 is the mean anomaly, i is the inclination, Ω is the right ascension of the ascending node, w is the argument of perigee, and t is the time;

[0083] The perturbation motion equation of the satellite is obtained from the Lagrange perturbation equation:

[0084]

[0085] Where, J2 is the band harmonic coefficient; R e is the average radius of the earth's equator; t is time; n is the average angular velocity of the satellite, and the calculation formula is:

[0086]

[0087] Where G is the gravitational constant, M is e is the mass of the Earth;

[0088] According to formula (26) to formula (31), the instantaneous root of the satellite orbit at the forecast start time t0 is used to calculate the instantaneous root of the satellite orbit at the current forecast time t0+m*Δt.

[0089] Furthermore, the specific operation process of step S3 is as follows:

[0090] According to the forecast result of step S2, the transformation matrix T1 from the geocentric inertial coordinate system to the satellite orbit coordinate system at the current forecast time t0+m*Δt is obtained:

[0091]

[0092] In the formula, i is the inclination, Ω is the right ascension of the ascending node, and u is the orbital angle. The calculation formula is:

[0093] u=f+w (33)

[0094] Where f is the true anomaly angle and w is the argument of perigee angle.

[0095] Furthermore, the specific operation process of step S4 is as follows:

[0096] The satellite attitude angle is the angle of the satellite relative to its own motion position, which is based on the satellite's center of mass as the coordinate origin. It includes the roll angle, pitch angle, and yaw angle. Based on the satellite attitude angle at the forecast start time t0 broadcast by the satellite platform, the satellite attitude correction matrix T2 at the current forecast time t0+m*Δt is calculated.

[0097]

[0098] In the formula, is the satellite roll angle, ψ is the satellite yaw angle, and θ2 is the satellite pitch angle.

[0099] Furthermore, the specific operation process of step S5 is as follows:

[0100] According to the lunar vector M1 in the geocentric inertial coordinate system, the transformation matrix T1 from the geocentric inertial coordinate system to the satellite orbit coordinate system, the satellite attitude correction matrix T2 and the conversion matrix T3 from the calibrated satellite body coordinate system to the instrument body coordinate system, the lunar vector M2 in the instrument body coordinate system is obtained by matrix calculation:

[0101] M2=T3T2T1M1 (35).

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

[0103] The invention discloses a method for predicting the high-precision lunar angle of a space remote sensing instrument. The method uses the time broadcast by the satellite platform to calculate the geocentric apparent right ascension and declination of the moon. The instantaneous root of the satellite orbit broadcast by the satellite platform and the simplified perturbation motion equation are then used to predict the instantaneous root of the satellite orbit, thereby obtaining a transformation matrix from the geocentric inertial coordinate system to the satellite orbit coordinate system. The satellite attitude correction matrix is ​​calculated by the satellite attitude angle, and the transformation matrix from the calibrated satellite body coordinate system to the instrument body coordinate system is combined to finally obtain the lunar angle in the instrument body coordinate system, thereby completing the high-precision lunar angle prediction of the space remote sensing instrument. Compared with the prior art, the invention has the following advantages:

[0104] 1. Simple calculation: The present invention optimizes the periodic terms of the lunar geocentric ecliptic longitude and ecliptic latitude, simplifies the Lagrangian perturbation equation, and uses the satellite attitude angle at the start of the forecast to calculate the satellite attitude correction matrix during the entire forecast duration, thereby simplifying the calculation process and reducing the complexity of the calculation.

[0105] 2. High prediction accuracy: During the calculation process, the present invention performs Venus perturbation correction on the sum of the lunar geocentric ecliptic longitude period terms and the sum of the lunar geocentric ecliptic latitude period terms. After completing the Venus perturbation correction, nutation correction is also performed. The correction operation can improve the accuracy of the lunar angle prediction and make the prediction result more accurate.

[0106] 3. Long-term prediction: The present invention simplifies the prediction process of the lunar angle and can perform long-term and high-precision lunar angle prediction, which is more suitable for space remote sensing instrument applications. BRIEF DESCRIPTION OF THE DRAWINGS

[0107] Figure 1 The present invention is a general flow chart of a high-precision lunar angle prediction method for a space remote sensing instrument.

[0108] Figure 2 This is a flow chart for calculating the lunar geocentric right ascension and lunar geocentric declination in step S1.

[0109] Figure 3 This is the relationship diagram between the eccentric anomaly angle E and the true anomaly angle f. DETAILED DESCRIPTION

[0110] The present invention is further described in detail below in conjunction with the accompanying drawings.

[0111] like Figure 1 As shown, a high-precision lunar angle prediction method for a space remote sensing instrument of the present invention mainly comprises the following steps:

[0112] Step S1: Calculate the lunar geocentric right ascension and lunar geocentric declination at the current forecast time t0+m*Δt according to the forecast start time t0 broadcast by the satellite platform and the required time interval Δt; wherein m is the number of forecast time intervals Δt, i.e., compared with the forecast start time t0, the current forecast time t0+m*Δt has increased m forecast time intervals Δt; m starts from 0, and after the lunar angle forecast at the current time is completed, m increases by 1 and is used for the lunar angle forecast at the next time, and so on.

[0113] The present invention uses a simplified ELP-2000 / 82 theory to calculate the lunar geocentric ecliptic longitude and lunar geocentric ecliptic latitude, and on this basis, improves the accuracy of the results through Venus perturbation correction and nutation correction, and finally converts them into lunar geocentric right ascension and lunar geocentric declination. The specific calculation process is as follows Figure 2 shown.

[0114] Step S1.1: First, the Julian day J is calculated based on the current forecast time t0+m*Δt D , according to the Julian date D Calculate the Julian century number T:

[0115] T=(J D -2451545) / 36525 (1)

[0116] Step S1.2: Calculate the lunar geocentric ecliptic longitude period term and lunar geocentric ecliptic latitude period term at the current forecast time t0+m*Δt. The specific calculation process is as follows:

[0117] Step S1.2.1: The expressions of the lunar geocentric ecliptic longitude and lunar geocentric ecliptic latitude in the ecliptic coordinate system are:

[0118] λ=A×sin(θ) (2)

[0119] β=A×sin(θ) (3)

[0120] Where, λ is the lunar geocentric ecliptic longitude, β is the lunar geocentric ecliptic latitude, A is the amplitude, and θ is the argument;

[0121] The calculation formula of the angle θ is:

[0122] θ=a×D+b×M+c×M'+d×F (4)

[0123] Where a, b, c, d are all argument coefficients, all in degrees, D is the sun-moon distance angle, M is the sun's mean near-earth angle, M' is the moon's mean near-earth angle, and F is the lunar ascending node's horizontal angular distance; the four argument coefficients a, b, c, d are given by the lunar geocentric ecliptic longitude cycle term or lunar geocentric ecliptic latitude cycle term coefficient table, and the calculation formulas for the remaining parameters D, M, M', and F are:

[0124] D is the distance angle between the sun and the moon:

[0125]

[0126] M is the sun's mean angle of approach to the earth:

[0127]

[0128] M' is the moon's angle of approach to the Earth:

[0129]

[0130] F is the horizontal angular distance of the lunar ascending node:

[0131]

[0132] Step S1.2.2: Correct the lunar geocentric ecliptic period term. The specific correction method is as follows:

[0133] If the argument θ contains M or -M (M is the sun's mean angle near the earth), each lunar geocentric ecliptic longitude cycle term needs to be multiplied by the coefficient E1 for correction; if the argument θ contains 2M or -2M (M is the sun's mean angle near the earth), each lunar geocentric ecliptic longitude cycle term needs to be multiplied by the square of the coefficient E1 for correction; the calculation formula of the coefficient E1 is:

[0134] E1=1-0.002516×T-0.0000074×T 2 (9)

[0135] The calculation formula for the sum of the lunar geocentric ecliptic period terms ∑I is:

[0136] ∑I=I1+I2+...+I 60 (10)

[0137] In the formula, I1, I2, ... I 60 They are the first to sixtieth geocentric ecliptic longitude cycle items of the moon respectively.

[0138] Step S1.2.3: Correct the lunar geocentric apparent ecliptic latitude period term. The specific correction method is as follows:

[0139] The calculation method for the sum of the lunar geocentric ecliptic latitude periodic terms ∑b is the same as the calculation method for the sum of the lunar geocentric ecliptic longitude periodic terms ∑I. Each lunar geocentric ecliptic latitude periodic term also contains an amplitude A and four angular coefficients a, b, c, and d. When the angular coefficient b of the sun's mean near-earth angle M is not 0, it also needs to be multiplied by the coefficient E1 or the square of the coefficient E1 for correction. Therefore, the specific correction method for the lunar geocentric ecliptic latitude periodic term is the same as step S1.2.2.

[0140] Step S1.3: Venus perturbation correction;

[0141] Due to the perturbation of Venus and Jupiter on the moon, it is necessary to make Venus perturbation correction for the calculated sum of the lunar geocentric ecliptic longitude period term ∑I and the sum of the lunar geocentric ecliptic latitude period term ∑b. The specific correction method is as follows:

[0142] ΔI=3958×sin(A1)+1962×sin(L'-F)+318×sin(A2) (11)

[0143]

[0144] Wherein, ΔI is the correction term for the sum of the lunar geocentric ecliptic longitude period terms. When making corrections, the sum of the lunar geocentric ecliptic longitude period terms ∑I and ΔI can be added. Δb is the correction term for the sum of the lunar geocentric ecliptic latitude period terms. When making corrections, the sum of the lunar geocentric ecliptic latitude period terms ∑b and Δb can be added. F is the lunar ascending node horizontal angle, M' is the lunar mean perigee angle, and L' is the lunar mean ecliptic longitude. The calculation formula is shown in formula (13). A1, A2 and A3 all represent the perturbation angle corrections. The calculation formulas are shown in formula (14), formula (15) and formula (16) respectively.

[0145]

[0146] A1=119.75+131.849×T (14)

[0147] A2=53.09+479264.29×T (15)

[0148] A3=313.45+481266.484×T (16)

[0149] Step S1.4: nutation correction;

[0150] After completing the perturbation correction of Venus, nutation correction is also required. Nutation can be decomposed into ecliptic nutation Δψ and angle nutation Δε. The present invention uses the optimized IAU1980 nutation theory to calculate nutation, reducing the cost of computing resources and improving the computing speed. The specific correction method is as follows:

[0151] Step S1.4.1: Calculation of ecliptic nutation;

[0152] Each ecliptic nutation period term S includes the calculation of the sine coefficient of the ecliptic nutation and the linear combination coefficient of the five basic angular distances for calculating the nutation amplitude angle θ1; the calculation formula of the nutation amplitude angle θ1 is:

[0153] θ1=a1×D+b1×M+c1×M’+d1×F+e1×Ω1 (17)

[0154] Where a1, b1, c1, d1, e1 are all nutation angular coefficients, D is the sun-moon distance angle, M is the sun's mean near-earth angle, M' is the moon's mean near-earth angle, F is the moon's ascending node's horizontal angular distance, Ω1 is the ecliptic and the moon's mean orbit ascending node's ecliptic longitude, the calculation formula for the ecliptic and the moon's mean orbit ascending node's ecliptic longitude Ω1 is:

[0155]

[0156] According to formula (19), the nutation period term S of the ecliptic longitude is calculated according to the nutation amplitude angle θ1:

[0157] S=(S1+S2×T)×sin(θ1) (19)

[0158] Where θ1 is the nutation amplitude angle, S1 and S2 are the sinusoidal coefficients of the ecliptic nutation.

[0159] By accumulating each ecliptic nutation period term S, we can get the ecliptic nutation Δψ, whose unit is 0.0001'.

[0160] Step S1.4.2: Calculation of the nutation angle;

[0161] Each cross angle nutation period term C includes the cosine coefficient of the cross angle nutation and the linear combination coefficient of the five basic angular distances of the nutation amplitude; the calculation formula of the cross angle nutation period term C is:

[0162] C=(C1+C2×T)×cos(θ1) (20)

[0163] Where θ1 is the nutation amplitude angle, and C1 and C2 are the cosine coefficients of the nutation angle.

[0164] By accumulating each angular nutation period term C, the angular nutation Δε can be obtained, and its unit is 0.0001'.

[0165] The invention carries out Venus correction and nutation correction on the predicted moon angle, thereby improving the accuracy of the prediction.

[0166] Step S1.5: After all corrections are completed, the lunar geocentric ecliptic longitude λ and lunar geocentric ecliptic latitude β at the current forecast time t0+m*Δt are finally obtained according to the calculation results of the above steps 1.1 to 1.4:

[0167] λ=L'+∑I / 1000000+Δψ (21)

[0168] β=∑b / 1000000 (22)

[0169] Where L' is the lunar mean ecliptic longitude, Δψ is the nutation of the ecliptic longitude, ∑b is the lunar geocentric apparent ecliptic latitude periodic term, and ∑I is the sum of the lunar geocentric apparent ecliptic longitude periodic terms.

[0170] Step S1.6: Calculate the lunar geocentric right ascension and lunar geocentric declination at the current forecast time t0+m*Δt; the calculation formulas for the lunar geocentric right ascension α and lunar geocentric declination δ are:

[0171]

[0172] Wherein, α is the lunar geocentric right ascension, δ is the lunar geocentric declination, λ is the lunar geocentric ecliptic longitude, β is the lunar geocentric ecliptic latitude, ε is the true ecliptic obliquity, which is equal to the sum of the mean ecliptic obliquity ε0 and the nutation Δε; the calculation formula of the mean ecliptic obliquity ε0 is:

[0173]

[0174] According to the calculation results of formula (1) to formula (24), the lunar geocentric apparent right ascension and declination in the geocentric inertial coordinate system at the current forecast time t0+m*Δt can be obtained, and then the lunar vector M1 in the geocentric inertial coordinate system can be obtained as follows:

[0175]

[0176] Where α is the lunar geocentric right ascension, and δ is the lunar geocentric declination.

[0177] Step S2: Calculate the instantaneous root of the satellite orbit at the current forecast time t0+m*Δt based on the instantaneous root of the satellite orbit at the forecast start time t0 broadcast by the satellite platform and the pre-established perturbation dynamic motion equation.

[0178] The instantaneous root of satellite orbit is a set of parameters used to describe the properties of satellite orbit shape, position and motion, including semi-major axis a2, eccentricity e, inclination i, right ascension of ascending node Ω, argument of perigee w and true anomaly f. Since the on-orbit computing resources of space remote sensing instruments are limited, high-precision complex calculations cannot be performed. Therefore, under the premise of taking into account the simplicity and accuracy of calculation, the present invention uses a simple high-precision satellite orbit instantaneous root prediction method. To facilitate the calculation, two auxiliary quantities are introduced, the eccentric anomaly E and the mean anomaly M3. The relationship between the eccentric anomaly E and the true anomaly f is as follows: Figure 3 As shown. Figure 3 ( Figure 3 In the figure, O is the center of the instantaneous elliptical orbit, O E is the center of the earth, b2 is the semi-minor axis), we can see that the relationship between the eccentric anomaly E, the mean anomaly M3 and the true anomaly f is:

[0179]

[0180] M3=Ee×sinE (27)

[0181]

[0182] Where f is the true anomaly, e is the eccentricity, E is the eccentric anomaly, and M3 is the mean anomaly.

[0183] Assuming that the perturbation function corresponding to the satellite in orbit is R, the Lagrangian perturbation equation is:

[0184]

[0185] Where a2 is the semi-major axis, e is the eccentricity, M3 is the mean anomaly, i is the inclination, Ω is the right ascension of the ascending node, w is the argument of perigee, and t is the time.

[0186] To simplify the calculation, only the influence of the earth's flattening perturbation J2 is considered. From the Lagrangian perturbation equation (29), the satellite's perturbation motion equation when only the earth's flattening perturbation J2 is considered is:

[0187]

[0188] Where, J2 is the band harmonic coefficient, constant; R e is the average radius of the earth's equator, a constant; t is time; n is the average angular velocity of the satellite, and the calculation formula for n is:

[0189]

[0190] Where G is the gravitational constant, M is e is the mass of the Earth, both are constants.

[0191] According to formula (26) to formula (31), the instantaneous root of the satellite orbit at the forecast start time t0 can be calculated.

[0192] Step S3: According to the forecast result of step S2, the transformation matrix T1 from the geocentric inertial coordinate system to the satellite orbit coordinate system at the current forecast time t0+m*Δt is obtained:

[0193]

[0194] Where u is the orbital angle, and its calculation formula is:

[0195] u=f+w (33)

[0196] Where f is the true anomaly angle and w is the argument of perigee angle.

[0197] Step S4: Calculate the satellite attitude correction matrix T2 at the current forecast time t0+m*Δt based on the satellite attitude angle at the forecast start time t0 broadcast by the satellite platform; the satellite attitude angle is based on the satellite center of mass as the coordinate origin, and describes the angle of the satellite relative to its own motion position, which are the roll angle Pitch angle θ2 and yaw angle ψ. When the satellite is in orbit, the attitude control method is used to keep the satellite attitude in a stable state. Therefore, it is considered that the satellite attitude angle is unchanged within the forecast duration. Therefore, according to the satellite attitude angle at the forecast start time t0 broadcast by the satellite platform, the satellite attitude correction matrix T2 within the entire forecast duration can be calculated:

[0198]

[0199] In the formula, is the satellite roll angle, ψ is the satellite yaw angle, and θ2 is the satellite pitch angle.

[0200] Step S5: Based on the calculation results of steps S1 to S4 and the conversion matrix from the calibrated satellite body coordinate system to the instrument body coordinate system, calculate the moon angle in the instrument body coordinate system at the current forecast time t0+m*Δt, and stop the calculation until the forecast end time t1.

[0201] Based on the calculated lunar vector M1 in the geocentric inertial coordinate system, the transformation matrix T1 from the geocentric inertial coordinate system to the satellite orbit coordinate system, the satellite attitude correction matrix T2, and the calibrated transformation matrix T3 from the satellite body coordinate system to the instrument body coordinate system, the lunar vector M2 in the instrument body coordinate system can be obtained by matrix calculation:

[0202] M2=T3T2T1M1 (35)

[0203] The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principle of the present invention. These improvements and modifications should also be regarded as the scope of protection of the present invention.

Claims

1. A high-precision lunar angle prediction method for space remote sensing instruments, characterized in that: The following steps are involved: Step S1: Calculate the lunar geocentric apparent right ascension and apparent declination at the current forecast time according to the forecast start time broadcast by the satellite platform and the required forecast time interval; Step S2: Calculate the instantaneous root of the satellite orbit at the current forecast time according to the instantaneous root of the satellite orbit at the forecast start time broadcast by the satellite platform and the pre-established perturbation dynamic motion equation; Step S3: Calculate the transformation matrix from the geocentric inertial coordinate system to the satellite orbit coordinate system at the current forecast time according to the instantaneous root of the satellite orbit at the current forecast time; Step S4: Calculate the satellite attitude correction matrix at the current forecast time according to the satellite attitude angle at the forecast start time broadcast by the satellite platform; Step S5: Calculate the moon angle in the instrument body coordinate system at the current forecast time based on the above results and the calibrated conversion matrix from the satellite body coordinate system to the instrument body coordinate system.

2. A high-precision lunar angle prediction method for space remote sensing instruments according to claim 1, characterized in that: In step S1, the lunar geocentric apparent right ascension and declination at the current forecast time t0+m*Δt are calculated based on the forecast start time t0 broadcast by the satellite platform and the required forecast time interval Δt, where m is the number of forecast time intervals Δt, i.e., compared with the forecast start time t0, the current forecast time t0+m*Δt has increased m forecast time intervals Δt; m starts from 0, and after the lunar angle forecast at the current moment is completed, m increases by 1 and is used for the lunar angle forecast at the next moment.

3. A high-precision lunar angle prediction method for space remote sensing instruments according to claim 2, characterized in that: The specific operation process of step S1 is as follows: Step S1.1: Calculate the Julian day J according to the current forecast time t0+m*Δt D , then the Julian century number T is: T=(J D -2451545) / 36525 (1) Step S1.2: Calculate the lunar geocentric ecliptic longitude period term and ecliptic latitude period term at the current forecast time t0+m*Δt; Step S1.2.1: The expressions of the lunar geocentric ecliptic longitude and latitude in the ecliptic coordinate system are: λ=A×sin(θ) (2) β=A×sin(θ) (3) Where, λ is the lunar geocentric ecliptic longitude, β is the lunar geocentric ecliptic latitude, A is the amplitude, and θ is the argument; The calculation formula of the argument θ is: θ=a×D+b×M+c×M'+d×F (4) Where a, b, c, d are all argument coefficients, all in degrees; D is the sun-moon distance angle, M is the sun's mean near-earth angle, M' is the moon's mean near-earth angle, and F is the moon's ascending node's mean angular distance; Step S1.2.2: Correct the lunar geocentric ecliptic period term; If the angle θ contains M or -M, each lunar geocentric ecliptic period term is multiplied by the coefficient E1 for correction; if the angle θ contains 2M or -2M, each lunar geocentric ecliptic period term is multiplied by the square of the coefficient E1 for correction; The calculation formula of the coefficient E1 is: E1=1-0.002516×T-0.0000074×T 2 (9) The calculation formula of the sum of the lunar geocentric ecliptic period terms ∑I is: ∑I=I1+I2+...+I 60 (10) In the formula, I1, I2, ... I 60 They are the first to the sixtieth geocentric ecliptic cycle items of the moon respectively; Step S1.2.3: Correct the lunar geocentric apparent ecliptic latitude period term; The correction method of the lunar geocentric ecliptic latitude period term is the same as step S1.2.2; the calculation method of the sum of lunar geocentric ecliptic latitude period terms ∑b is the same as the calculation method of the sum of lunar geocentric ecliptic longitude period terms ∑I; Step S1.3: Venus perturbation correction; Step S1.4: nutation correction; Step S1.5: According to the above results, the lunar geocentric ecliptic longitude and ecliptic latitude at the current forecast time t0+m*Δt are obtained: λ=L'+∑I / 1000000+Δψ (21) β＝∑b / 1000000 (22) Where, λ is the lunar geocentric ecliptic longitude, L' is the lunar mean ecliptic longitude, ∑I is the sum of the lunar geocentric ecliptic longitude periodic terms, Δψ is the nutation of the ecliptic longitude, β is the lunar geocentric ecliptic latitude, ∑b is the sum of the lunar geocentric ecliptic latitude periodic terms; Step S1.6: Calculate the lunar geocentric apparent right ascension and declination at the current forecast time t0+m*Δt; The calculation formulas for the lunar geocentric apparent right ascension and declination are: Where, α is the lunar geocentric right ascension, δ is the lunar geocentric declination, λ is the lunar geocentric ecliptic longitude, β is the lunar geocentric ecliptic latitude, ε is the true ecliptic obliquity, which is equal to the sum of the mean ecliptic obliquity ε0 and the nutation Δε. The calculation formula of the mean obliquity angle ε0 is: According to the above results, the calculation formula of the lunar vector M1 in the geocentric inertial coordinate system is:

4. A high-precision lunar angle prediction method for space remote sensing instruments according to claim 3, characterized in that: In step S1.2.1, the calculation formulas for the sun-moon distance angle D, the sun's mean near-earth angle M, the moon's mean near-earth angle M', and the moon's ascending node's mean angular distance F are respectively:

5. A high-precision lunar angle prediction method for space remote sensing instruments according to claim 4, characterized in that: The specific operation process of step S1.3 is as follows: ΔI=3958×sin(A1)+1962×sin(L'-F)+318×sin(A2) (11) Where ΔI is the correction term for the sum of the lunar geocentric ecliptic period terms. When making the correction, the sum of the lunar geocentric ecliptic period terms ∑I is added to ΔI. Δb is the correction term for the sum of the lunar geocentric apparent ecliptic latitude period terms. During the correction, the sum of the lunar geocentric apparent ecliptic latitude period terms ∑b is added to Δb. L' is the lunar mean ecliptic longitude, and its calculation formula is formula (13). A1, A2 and A3 all represent the perturbation angle correction, and their calculation formulas are (14), (15) and (16) respectively: A1=119.75+131.849×T (14) A2=53.09+479264.29×T (15) A3=313.45+481266.484×T (16).

6. A high-precision lunar angle prediction method for space remote sensing instruments according to claim 5, characterized in that: The specific operation process of step S1.4 is as follows: Step S1.4.1: Calculation of ecliptic nutation; Each ecliptic nutation period term S includes the calculation of the sine coefficient of the ecliptic nutation and the linear combination coefficient of the five basic angular distances for calculating the nutation amplitude angle; the calculation formula for the nutation amplitude angle θ1 is: θ1=a1×D+b1×M+c1×M'+d1×F+e1×Ω1 (17) Where a1, b1, c1, d1, and e1 are all nutation angular coefficients, and Ω1 is the ecliptic longitude of the ascending node of the ecliptic and the lunar orbit. The calculation formula is: The calculation formula of the ecliptic nutation period term S is: S=(S1+S2×T)×sin(θ1) (19) Where θ1 is the nutation amplitude angle, S1 and S2 are the sinusoidal coefficients of the ecliptic nutation; The ecliptic nutation period term S is accumulated to obtain the ecliptic nutation Δψ, whose unit is 0.0001'; Step S1.4.2: Calculation of the nutation angle; Each cross angle nutation period term C includes the cosine coefficient of the cross angle nutation and the linear combination coefficient of the five basic angular distances of the nutation amplitude; the calculation formula of the cross angle nutation period term C is: C=(C1+C2×T)×cos(θ1) (20) Where θ1 is the nutation amplitude angle, and C1 and C2 are the cosine coefficients of the nutation angle.

7. A high-precision lunar angle prediction method for space remote sensing instruments according to claim 6, characterized in that: The specific operation process of step S2 is as follows: Introducing the auxiliary quantities eccentric anomaly angle and flat anomaly angle, the relationship between the eccentric anomaly angle, flat anomaly angle and true anomaly angle is: M3=Ee×sinE (27) Where, f is the true anomaly, e is the eccentricity, E is the eccentric anomaly, and M3 is the mean anomaly; Assuming that the perturbation function corresponding to the satellite in orbit is R, the Lagrangian perturbation equation is: Where a2 is the semi-major axis, e is the eccentricity, M3 is the mean anomaly, i is the inclination, Ω is the right ascension of the ascending node, w is the argument of perigee, and t is the time; The perturbation motion equation of the satellite is obtained from the Lagrange perturbation equation: Where, J2 is the band harmonic coefficient; R e is the average radius of the earth's equator; t is time; n is the average angular velocity of the satellite, and the calculation formula is: Where G is the gravitational constant, M is e is the mass of the Earth; According to formula (26) to formula (31), the instantaneous root of the satellite orbit at the forecast start time t0 is used to calculate the instantaneous root of the satellite orbit at the current forecast time t0+m*Δt.

8. A high-precision lunar angle prediction method for space remote sensing instruments according to claim 7, characterized in that: The specific operation process of step S3 is as follows: According to the forecast result of step S2, the transformation matrix T1 from the geocentric inertial coordinate system to the satellite orbit coordinate system at the current forecast time t0+m*Δt is obtained: In the formula, i is the inclination, Ω is the right ascension of the ascending node, and u is the orbital angle. The calculation formula is: u=f+w(33) Where f is the true anomaly angle and w is the argument of perigee angle.

9. A high-precision lunar angle prediction method for space remote sensing instruments according to claim 8, characterized in that: The specific operation process of step S4 is as follows: The satellite attitude angle is the angle of the satellite relative to its own motion position, which is based on the satellite's center of mass as the coordinate origin. It includes the roll angle, pitch angle, and yaw angle. Based on the satellite attitude angle at the forecast start time t0 broadcast by the satellite platform, the satellite attitude correction matrix T2 at the current forecast time t0+m*Δt is calculated. In the formula, is the satellite roll angle, ψ is the satellite yaw angle, and θ2 is the satellite pitch angle.

10. A high-precision lunar angle prediction method for space remote sensing instruments according to claim 9, characterized in that: The specific operation process of step S5 is as follows: According to the lunar vector M1 in the geocentric inertial coordinate system, the transformation matrix T1 from the geocentric inertial coordinate system to the satellite orbit coordinate system, the satellite attitude correction matrix T2 and the conversion matrix T3 from the calibrated satellite body coordinate system to the instrument body coordinate system, the lunar vector M2 in the instrument body coordinate system is obtained by matrix calculation: M2=T3T2T1M1 (35).

Images