A prediction method for imaging windows of Earth observation satellites

Abstract

This invention belongs to the field of satellite imaging algorithms, specifically to a method for predicting the imaging window of an imaging satellite for Earth observation. This method is used to image an imaging target and includes the following steps: Step 1: Rough Search; Step 2: Image Estimation; and Step 3: Image Correction. Compared to existing technologies, this method addresses the issues of imaging satellite mission planning, which often lack consideration of constraints, resulting in large computational errors, complex algorithms, and poor imaging performance. This solution considers perturbations and errors caused by factors such as Earth's oblateness, orbital eccentricity, and atmospheric drag, and further utilizes a simple and fast iterative algorithm to solve the satellite imaging window, thereby improving the efficiency and accuracy of autonomous missions.

Description

Technical Field

[0001] The present invention belongs to the field of satellite imaging algorithms, and in particular relates to a method for predicting an imaging window for an imaging satellite to observe the earth. Background Art

[0002] Among the many operational satellites, Earth observation satellites are particularly common. Earth observation satellites (EOS) are a type of satellite used to monitor and study Earth's surface conditions and atmospheric environment. These satellites collect diverse data across the visible spectrum, infrared spectrum, and electromagnetic wavebands using a variety of onboard sensors. Earth observation satellites currently play a vital role in scientific research, military operations, and reconnaissance. In scientific research, these satellites help scientists study Earth's climate change, environmental changes, ecosystems, and natural disasters. In military reconnaissance, they provide high-resolution surface imagery and other data to aid in monitoring and analyzing military operations. In agriculture, Earth observation satellites monitor crop growth, land use, and water resource management. Due to the complex variety and large number of these satellites, they play a key role in various fields. Different types of sensors onboard provide different types of data, meeting diverse mission requirements. Some Earth observation imaging satellites are used for disaster monitoring, such as detecting the frequency of forest fires. Forest fires currently account for one-fifth of global carbon dioxide emissions, making the use of Earth observation satellites to monitor forest fires a crucial task.

[0003] Satellite Earth observation activities require mission planning. Constraints such as resource allocation, satellite orbit information, and ground station resources are factored in. Within these constraints, mission planning is solved for the satellite and the optimal mission execution plan is developed. This process is a complex combinatorial optimization problem, involving multiple constraints and resource scheduling, as well as multiple tasks, multiple time windows, and multiple optimization objectives. Throughout the history of satellite mission planning, ground stations have played a crucial role. Most satellites must immediately transmit collected information to ground stations. Upon receiving the information, the ground stations combine the satellite coordinates, orbital information, and platform payload with the parsed mission planning constraints. The ground stations then upload the information, along with pre-defined mission instructions, to the space-based satellite platform. The satellite platform then completes the mission step by step according to the uploaded instructions and then transmits the collected mission information back to the ground station. This is the traditional operational model for satellite mission planning. The continuous information transmission back to the ground station represents a large loop between space (satellite) and ground (ground station). It's easy to see that traditional approaches consume both time and human resources, leading to inefficient satellite missions. This inability to gather sufficient information within a limited timeframe leaves little room for emergency response. Furthermore, traditional approaches rely on a passive task scheduling process, lacking real-time access to satellite information and resource status, and preventing the ability to adjust observation plans on the fly, ultimately leading to a waste of satellite resources.

[0004] In reality, humanity's pursuit of satellite efficiency and resource utilization has never ceased. For increasingly complex space missions, most satellite system command and control instructions, such as attitude maneuvers, data downlinks, orbital maneuvers, and resource scheduling, will no longer rely solely on ground-based operations. However, for previous Earth observation satellites, fully automated autonomous mission planning was beyond the technological capabilities of the time. Furthermore, early satellites were large in size and carried a vast number of sensors, necessitating significant trial-and-error costs. Therefore, engineers have long employed improved solutions based on traditional scheduling schemes that relied solely on ground stations. These improvements increased the integration of mission instructions, optimized mission planning algorithms, and reduced the number of calculations required for equations related to satellite orbits to achieve improved efficiency.

[0005] In recent years, with the rapid development of microsatellites, the number of satellites in orbit has increased. These satellites are often required to form constellations for space missions. Furthermore, more and more users are placing higher demands on the complexity and diversity of satellite observation missions. Under such circumstances, if satellites still rely entirely on ground-based base stations for observation missions, the ground-based base stations will become overwhelmed. At the same time, the increasing variety of remote sensor types and the growing demand for user information services have brought tremendous prospects and opportunities for satellite observation missions, but also brought unprecedented challenges. These challenges are mainly reflected in the following aspects:

[0006] 1) Observation task requests are sudden, and it is necessary to rationally utilize the relationship between existing planning schemes and randomly arriving observation task requests.

[0007] 2) The weather conditions at the observation point determine whether the satellite can observe the target, and this type of information inconsistency needs to be processed.

[0008] 3) Communication interference and delay affect the reliability and robustness of the control system. Summary of the Invention

[0009] The purpose of the present invention is to provide a method for predicting the imaging window of an imaging satellite for Earth observation in order to solve at least one of the above problems, so as to address the problems in the existing technology that imaging satellite mission planning takes less consideration of constraints, has large calculation errors, complex algorithms, and poor imaging effects. This solution takes into account the perturbations and errors caused by factors such as the Earth's oblateness, orbital eccentricity, and atmospheric drag, and further uses a simple and fast iterative algorithm to solve the satellite imaging window to improve the efficiency and accuracy of autonomous missions.

[0010] The purpose of the present invention is achieved through the following technical solutions:

[0011] A method for predicting an imaging window for an imaging satellite to observe the Earth, for imaging an imaging target located on the Earth's surface, comprises the following steps:

[0012] Step 1: A rough search is performed using an analytical method, treating the Earth as a perfect sphere and the satellite orbit as a perfect circle, to search for a transiting satellite that meets the imaging window conditions.

[0013] Step 2: Image estimation, using analytical methods, considering the Earth as an ideal sphere and the satellite orbit as an ideal circle, and estimating the satellite's image point based on the results of a rough search;

[0014] Step 3: Imaging correction, using numerical algorithms to compensate for imaging estimation errors through orbit extrapolation models;

[0015] When the satellite does not have the ability to maneuver or the imaging target is beyond the satellite's attitude maneuvering range: Steps 1 and 2: Use the satellite's camera optical axis pointing as the initial position. Based on the satellite camera's inclination and FOV, obtain the angular deviation through the great circle, then obtain the longitude deviation on the latitude TLL of the imaging target, and finally obtain the satellite imaging window estimate. Step 3: Calculate the precise coordinates of the satellite using an orbital extrapolation model. This orbital extrapolation model considers secular perturbations, long-period perturbations, and short-period perturbations caused by the Earth's oblateness, orbital eccentricity, and atmospheric drag. Substitute the estimated imaging window value into the model and use the iterative equation to solve the imaging window.

[0016] The conditions of the imaging window are: L ·V L =0, where r L V is the position vector of the line connecting the imaging target and the projection point of the camera optical axis, L is the velocity vector of the projection point of the camera optical axis;

[0017] When the satellite has attitude maneuverability: Step 1: Determine the limit imaging angle of the satellite camera and search for a transiting satellite that meets the imaging window conditions; Step 2: The great circle containing the imaging target and the satellite's projection point on the Earth is perpendicular to the satellite orbit plane. The parameters of the great circle projected on the ground by the imaging target relative to the camera's optical axis are obtained, the number of satellite imaging windows is obtained, and the approximate imaging point is estimated; Step 3: Apply the orbit extrapolation model to obtain the imaging window;

[0018] The imaging window condition is: the moment when the camera optical axis coincides with the target tilt vector, and the target tilt vector is the position vector of the line connecting the satellite and the imaging target.

[0019] Preferably, when the satellite does not have attitude maneuvering capability, a two-body motion model is used:

[0020] Step 1: When the longitude offset angle between the imaging target and the satellite trajectory at the TLL satisfies: Δv1 ≥ dv ≥ Δv2, the satellite has the imaging window condition; where Δv1 and Δv2 are the longitude offset angles at the TLL corresponding to the extreme imaging angles of the satellite camera, and dv is the longitude difference angle between the satellite trajectory and the imaging target after the satellite orbits the Earth N times;

[0021] Step 2: Estimate the imaging point based on the satellite's phase angle α at the TLL, the initial longitude offset angle Δv between the satellite and the imaging target at the TLL at time t0, and the satellite's orbit inclination angle i;

[0022] Step 3: Obtain the precise position coordinates of the satellite through the orbit extrapolation model, calculate the time error Δt from the current satellite position to the position required by the imaging window conditions, add Δt to the estimated time, and substitute it into the orbit extrapolation model again for iteration until Δt is within the allowable tolerance range. The iteration is completed.

[0023] Preferably, in step 1:

[0024] On the surface, the imaging target is set as T, the intersection of the satellite trajectory and the TLL is set as C, and the intersection of the satellite trajectory and the equatorial plane is set as D;

[0025] DV satisfies:

[0026]

[0027] Where, is the Earth's rotation rate, n is the mean anomaly angular velocity, v s is the longitude of the satellite trajectory, v T is the longitude of the imaging target, and N is the number of satellite revolutions since time t0;

[0028] The limiting imaging angles θ1 and θ2 satisfy:

[0029]

[0030] Where τ is the inclination angle of the satellite camera, γ is the FOV of the satellite camera, R is the radius of the earth, and h is the satellite altitude;

[0031] Convert the limiting imaging angle to the longitude offset angle Δv at the TLL:

[0032]

[0033] in,

[0034]

[0035] where θ is taken as θ1 and θ2 respectively; φ is the latitude of the imaging target, i is the satellite orbit inclination angle, p, η and l are the great circle angles of TC, TD and CD respectively, Γ is the rotation angle between TD and CD in the spherical triangle TCD; sign(θ) function is used to take the sign of θ.

[0036] Preferably, when calculating the limiting imaging angle, a margin of 0.5-1 degrees is added to the FOV of the satellite camera.

[0037] Preferably, in step 2:

[0038] On the surface, the following settings are made: the imaging target is T, the intersection of the satellite trajectory and the TLL is C, the intersection of the satellite trajectory and the equatorial plane is D, the imaging point is S, the intersection of the meridian of point C and the equatorial plane is B, and the Earth's pole is N;

[0039] The initial longitude offset angle Δv is used as the actual longitude difference angle dv;

[0040] cosp=sin 2 φ+cos 2 φ·cosdv;

[0041] The rotation angle Σ between CD and CB in spherical triangle BCD is:

[0042]

[0043] The rotation angle Λ between CN and CT in the spherical triangle NCT is:

[0044]

[0045] The rotation angle π between CT and CS in the spherical triangle TCS is:

[0046] Π=π-Λ-Σ;

[0047] The value m of the large fillet CS is:

[0048] 1) When 0≤α S <π / 2 or π≤α S <3π / 2, and dv≥0, and 0<i<π / 2:

[0049] m=tan -1 [tanp·cos(Λ+Σ)];

[0050] 2) When 0≤α S <π / 2 or π≤α S <3π / 2, and dv<0, and 0<i<π / 2:

[0051] m=tan -1 [tanp·cos(Λ-Σ)];

[0052] 3) When π / 2≤α S <π or 3π / 2≤α S <2π, and dv≥0, and 0<i<π / 2:

[0053] m=-tan -1 [tanp·cos(Λ-Σ)];

[0054] 4) When π / 2≤α S <π or 3π / 2≤α S<2π, and dv<0, and 0<i<π / 2:

[0055] m=-tan -1 [tanp·cos(Λ+Σ)];

[0056] 5) When 0≤α S <π / 2 or π≤α S <3π / 2, and dv≥0, and π / 2<i<π:

[0057] m=tan -1 [tanp·cos(Λ-Σ)];

[0058] 6) When 0≤α S <π / 2 or π≤α S When <3π / 2, and dv<0, and π / 2<i<π:

[0059] m=tan -1 [tanp·cos(Λ+Σ)];

[0060] 7) When π / 2≤α S <π or 3π / 2≤α S <2π, and dv≥0, and π / 2<i<π:

[0061] m=-tan -1 [tanp·cos(Λ+Σ)];

[0062] 8) When π / 2≤α S <π or 3π / 2≤α S <2π, dv<0, and π / 2<i<π:

[0063] m=-tan -1 [tanp·cos(Λ-Σ)];

[0064] Where φ is the latitude of the imaging target, i is the satellite orbit inclination angle, p and l are the large circle angles of TC and CD respectively, and α S is the phase angle of the satellite at point S;

[0065] The imaging point S satisfies:

[0066] α S =α C +m;

[0067] Where, α C are the phase angles of the satellite at point C respectively.

[0068] Preferably, in step 3:

[0069] The orbital extrapolation model is:

[0070]

[0071] Where a is the semi-major axis of the orbit, ρ is the radius variation parameter, A is the amplitude of the orbit radius variation, α is the phase angle, and α p is the true anomaly, κ is the secular offset constant, θ is the orbital precession constant, χ is the long-period perturbation constant, Δ is the short-period perturbation constant, B is the atmospheric drag factor, Ω is the right ascension of the ascending node of the satellite orbit, i is the orbital inclination angle, λ is the satellite latitude parameter, r is the radial distance of the satellite in the orbital coordinate system, and n is the angular velocity;

[0072] The satellite coordinates obtained by the orbit extrapolation model are time functions based on the orbital coordinate system LO. The obtained satellite coordinates are sequentially transformed into the Earth-centered inertial system ECI, the Earth-fixed non-inertial system ECEF, and the local tangent coordinate system LT of the imaging target.

[0073] Specifically, the formula for converting LO coordinates to ECEF is as follows:

[0074] ξ=rcosλ;η=rsinλ;

[0075]

[0076] Where, v is the ephemeris time;

[0077] The coordinate transformation matrix for converting ECEF coordinates to LT coordinates is as follows:

[0078]

[0079] Where λ is the target longitude and φ is the target latitude;

[0080] Calculate the unit direction vector of the camera optical axis according to the information of the satellite camera, and then transform it into the LO, ECI, ECEF and LT coordinate systems in turn to obtain r L and V L ;

[0081] Among them, V L Calculated according to the following formula:

[0082]

[0083] Where Δt' is a given small time interval, for example, 1 second or 2 seconds;

[0084] Calculate the time error Δt from the current satellite position to the position required by the imaging window conditions:

[0085]

[0086] The time error Δt is added to the estimated time and substituted into the orbit extrapolation model again for iteration.

[0087] Preferably, when the satellite has attitude maneuverability, the imaging area of ​​the satellite camera is a large circular belt:

[0088] Step 1: When the longitude offset angle between the imaging target and the satellite trajectory at the TLL satisfies: Δv1 ≥ dv ≥ Δv2, the satellite has the imaging window condition; where Δv1 and Δv2 are the longitude offset angles at the TLL corresponding to the extreme imaging angles of the satellite camera, and dv is the longitude difference angle between the satellite trajectory and the imaging target after the satellite orbits the Earth N times;

[0089] According to the relationship between the position of the imaging target and the large circular ring, it is divided into 0, 1 and 2 imaging windows;

[0090] When the imaging window is 0 times, the imaging target cannot be imaged;

[0091] When the imaging window is 1 time,

[0092] Step 2: Obtain the initial yaw angle Ψ0 so that the satellite points the camera toward the imaging target. Estimate the imaging point based on the satellite's phase angle α at the TLL, the initial longitude offset angle Δv between the satellite and the imaging target at time t0, and the satellite's orbit inclination angle i.

[0093] Step 3: Obtain the precise position coordinates of the satellite through the orbit extrapolation model, calculate the time error Δt from the current satellite position to the position required by the imaging window conditions, add Δt to the estimated time and substitute it into the orbit extrapolation model again to iterate until Δt is within the allowable tolerance range, and the iteration is completed;

[0094] When the imaging window is 2 times,

[0095] Step 2: Obtain the initial yaw angle Ψ0 so that the satellite points the camera toward the imaging target. Estimate the imaging point based on the distance d1 from the projection point of the camera optical axis to the satellite trajectory and the distance d from the imaging target to the satellite trajectory.

[0096] Step 3: Define the imaging condition, -r·s=|r||s|cosτ, where r is the position vector of the satellite in the Earth-fixed non-inertial system ECEF, and s is the tilt vector from the satellite to the imaging target in the Earth-fixed non-inertial system ECEF. Both r and s are functions of time. Calculate the distance e1 from the camera optical axis projection point to the imaging target and the distance e2 from the camera optical axis projection point to the satellite trajectory in the LH plane. Then calculate the distance error Δd in the LH plane and obtain the time error Δt. Substitute the time error Δt into the imaging condition and iterate until the tilt vector s is close to the camera optical axis. Then convert the tilt vector s to the Earth-centered inertial system ECI and the orbital coordinate system LO to obtain the yaw angle, and finally obtain the total yaw angle Ψ.

[0097] Preferably, in step 1:

[0098] On the surface, the imaging target is set as T, the intersection of the satellite trajectory and the TLL is set as C, and the intersection of the satellite trajectory and the equatorial plane is set as D;

[0099] DV satisfies:

[0100]

[0101] Where, is the Earth's rotation rate, n is the mean anomaly angular velocity, v s is the longitude of the satellite trajectory, v T is the longitude of the imaging target, and N is the number of satellite revolutions since time t0;

[0102] The limiting imaging angles θ1 and θ2 satisfy:

[0103]

[0104] Where τ is the inclination angle of the satellite camera, γ is the FOV of the satellite camera, R is the radius of the earth, and h is the satellite altitude;

[0105] Convert the limiting imaging angle to the longitude offset angles Δv1 and Δv2 at the TLL:

[0106]

[0107] in,

[0108]

[0109] where φ is the latitude of the imaging target, i is the satellite orbit inclination angle, p, η, and l are the large circle angles of TC, TD, and CD, respectively, and Γ is the rotation angle between TD and CD in the spherical triangle TCD.

[0110] Preferably, in step 2:

[0111] Get the position vector r of the line connecting the imaging target and the subsatellite point S and the subsatellite velocity vector V S And the judgment condition z is calculated by vector cross product:

[0112] z=|r S ×V S |=r x V y -r y V x ;

[0113] According to the calculated sign of z and the satellite camera inclination angle τ, the initial yaw angle Ψ0 is determined:

[0114] When τ and z have the same sign, then Ψ0 = 0;

[0115] When τ and z have opposite signs, then Ψ0 = π.

[0116] Preferably, in step 3:

[0117] When the imaging window is 2 times,

[0118] Calculate the angle ξ between r and s, and then calculate the angle difference ∈ based on s and the direction of the satellite camera:

[0119]

[0120] Where, σ is the elevation angle from the imaging target to the satellite;

[0121] Based on the angle ν between the distance from the ideal camera optical axis projection point to the imaging target and the satellite velocity, Δd is calculated:

[0122] When ν<π / 2:

[0123]

[0124] When ν>π / 2:

[0125]

[0126] Then calculate Δt:

[0127]

[0128] Total yaw angle Ψ:

[0129] Ψ=Ψ'+Ψ0.

[0130] Compared with the prior art, the present invention has the following beneficial effects:

[0131] In the above scheme, the present invention can realize the solution of mission planning for imaging satellite earth observation, and proposes a fast prediction algorithm for predicting the optimal imaging window of the imaging satellite and giving the attitude angle when possible. When establishing the mathematical model, constraints such as camera tilt angle, FOV, and yaw and roll attitude maneuvers are taken into account. The numerical results of analytical modeling of spatial geometry and iterative algorithms are used to simplify the mission planning process. The imaging prediction algorithm is mainly divided into three steps: rough search, imaging estimation, and imaging correction. In the first two steps, the satellite orbit is assumed to be an ideal circle and the earth is assumed to be an ideal sphere. In the third step, perturbations and errors caused by factors such as the earth's oblateness, orbital eccentricity, and atmospheric drag are taken into account and compensated. The prediction results are optimized by relying on a more detailed earth shape and accurate orbit extrapolation equations.

[0132] Many traditional studies on imaging satellite mission planning problems have paid little attention to constraints. When designing the algorithm, the present invention considers constraints such as the tilt angle of the camera installed on the satellite, the availability of satellite attitude maneuvers, and the number of imaging times of a single satellite pass to establish a model and obtain the imaging window.

[0133] When studying satellite imaging mission planning, the present invention first assumes the earth and satellite orbits to be ideal spheres and ideal circles, uses analytical methods of solid geometry to solve initial values, and then uses iterative methods to solve the exact numerical solution, which greatly speeds up the calculation speed, reduces the demand for computing resources, and improves the efficiency and accuracy of autonomous tasks.

[0134] Computational experiments were conducted to evaluate the algorithm in real-world satellite imaging scenarios. The results demonstrate that the algorithm achieves the desired results and can find the optimal solution for the imaging window in each operating mode. This algorithm maintains accuracy while simplifying the computational steps. Further validation of its performance is possible through local computation on microsatellite platforms. BRIEF DESCRIPTION OF THE DRAWINGS

[0135] Figure 1 Schematic diagram of imaging conditions in the LH plane in the embodiment.

[0136] Figure 2 is the change of the satellite trajectory on the TLL in the embodiment.

[0137] Figure 3 Schematic diagram of the limiting imaging angle of the imaging strip in the embodiment.

[0138] Figure 4 2 is an auxiliary schematic diagram for solving spherical triangulation in an embodiment.

[0139] Figure 5 LT coordinate system schematic diagram in the embodiment.

[0140] Figure 6It is a circular strip obtained by projecting the optical axis of the camera of the yaw maneuverable satellite (with attitude maneuverability) around its zenith axis onto the LH plane in the embodiment.

[0141] Figure 7 (a) is a spatial schematic diagram in the LT coordinate system in the embodiment, and (b) is a spatial schematic diagram in the LH plane in the embodiment.

[0142] Figure 8 is the comparison of the right ascension Ω of the ascending node of the satellite orbit in the test case.

[0143] Figure 9 is the comparison of the satellite orbit inclination angle i in the test case.

[0144] Figure 10 is the comparison of satellite dimensionality parameter (equivalent phase angle) λ in the test case.

[0145] Figure 11 is the change in the radial distance r of the satellite in the orbital coordinate system in the experimental example.

[0146] Figure 12 It is the residual change during the iteration of the imaging window at 2000-09-08 12:21:40.265 in the test case.

[0147] Figure 13 It is the residual change during the iteration of the imaging window at 2000-09-10 12:04:38.710 in the test case.

[0148] Figure 14 It is the residual change during the iteration of the imaging window at 2000-09-10 12:05:48.870 in the test case. DETAILED DESCRIPTION

[0149] The present invention will be described in detail below with reference to the accompanying drawings and specific embodiments.

[0150] Example

[0151] Definition 1: This scheme generally defines the satellite camera imaging window as two cases: the first case is when the satellite camera's optical axis can accurately point to the designated target (imaging target, located on the surface, hereinafter referred to as the target); the second case is when the target is within the camera's field of view (FOV) and the projection distance between the target and the camera's optical axis on the ground reaches the minimum time point.

[0152] Both of these situations can be considered as satellite imaging windows, and they meet different prerequisites:

[0153] Premise 1: In the first case, the satellite's attitude maneuvering is available;

[0154] Premise 2: In the second case, the satellite has no available attitude maneuver or the target position exceeds the satellite's attitude maneuvering tolerance, for example, when the satellite is specified to have a maximum roll angle due to image distortion considerations.

[0155] Definition 2: In the first case, the imaging window is defined as the point where the satellite's camera optical axis coincides with the target's tilt vector. In this scheme, the target tilt vector is defined as the vector connecting the satellite and the ground target. A more intuitive statement is that the satellite's camera optical axis is directly aligned with the target, and its projection on the ground has zero deviation from the target.

[0156] Definition 3: In the second case, the solution of the imaging window becomes complicated. When the satellite flies along its own orbit, the projection point of the camera optical axis on the ground can be considered as a straight line in a short period, such as Figure 1 As shown. The projection point of the camera optical axis will move along this straight line at a speed of V L Move, and then in the imaging window, the distance between the projection point and the target should be the minimum; because the perpendicular line from the point to the straight line is the shortest, the position vector r of the projection point at this time L Should be perpendicular to its velocity vector V L Therefore, the conditional equation of the imaging window in this case is:

[0157] r L ·V L =0 (1).

[0158] Based on this:

[0159] Theory 1: Since the LT coordinate system is widely used in this scheme, its elements and characteristics are briefly described below. LT represents the local tangent coordinate system, whose origin is at the target of a given longitude and latitude on the surface of the earth, with the x-axis pointing due east and the y-axis pointing due north. The axes are given by the right-hand system and point vertically to the sky, such as Figure 5 In addition, its xy plane is the so-called LH plane.

[0160] 1. For the case where the satellite does not have attitude maneuvers:

[0161] Assumption 1: Most satellite orbit solutions can be approximated using two-body motion. During the rough search process, this solution also builds a model based on two-body motion.

[0162] The basic principle of coarse search is as follows:

[0163] Theory 2: For a given satellite, after each orbital period T, it will re-pass the same point in the inertial coordinate system (that is, the coordinate system that does not follow the rotation of the earth), such as Figure 2As the Earth rotates along its axis, the target moves along the target latitude line (TLL) mentioned above, which is defined as a small circle of constant latitude passing through the target location. Therefore, in inertial space, each time the satellite re-passes the TLL, the Earth's rotation will cause the target to move closer to the longitude of the satellite at the TLL.

[0164] Step 1.1: Due to the principle of relative motion, on the other hand, the target is or The speed of the approach to its own orbit, is the Earth's rotation rate, and n is the mean anomaly angular velocity. Therefore, in the Earth-fixed coordinate system (ECEF), for the target, the longitude of the satellite when it passes through the TLL has a step size of This relationship can be described by the following formula:

[0165]

[0166] Where: Δv is the initial longitude difference (longitude offset angle) between the satellite and the specific ground target when it first passes through the TLL at time t0, v s is the longitude of the satellite's trajectory on Earth, v T is the longitude of the target, N is the number of times the satellite has orbited the Earth since t0, and dv is the longitude difference between the satellite trajectory and the target after the satellite has orbited N times.

[0167] Step 1.2: Next, consider the processing of satellite camera information. The satellite flies along its orbit, and the camera's field of view (FOV) covers a long strip on the Earth's surface, with a cross section like Figure 3 As shown. In this way, the two strip boundaries passing through the great circle, that is, the limit imaging angles θ1 and θ2 can be expressed by the following formula:

[0168]

[0169] Where: τ is the camera tilt angle, γ is the camera FOV, R is the radius of the earth, and h is the satellite altitude.

[0170] Assumption 2: Because R and r (the radial distance of the satellite in the orbital coordinate system) vary as the satellite moves along its orbit, these two angles are not constant with respect to time. However, for low-Earth Earth observation satellites (LEO), their orbits are typically near-circular polar orbits. For simplicity, the orbital semi-major axis and the Earth's equatorial radius are used in the previous equations to obtain a constant deviation angle. Thus, as described above, the Earth is assumed to be a perfect sphere, and the satellite's orbit is assumed to be a perfect circle.

[0171] Step 1.3: In order to connect the imaging boundary angle (limit imaging angle) with the satellite longitude position on the TLL, the boundary angle is further converted into the longitude offset angle along the TLL, as Figure 4 As shown (forming a spherical triangle TCD). The following equation is used to convert the imaging boundary angles θ1 and θ2 into offset angles:

[0172]

[0173] cosp=coslcosη+sinlsinηcosΓ (5);

[0174]

[0175] Where: φ is the target latitude; i is the satellite orbit inclination, p, η and l are the corresponding large circle angles along TC, TD and CD respectively, and Γ is the rotation angle between TD and CD in the spherical triangle TCD, as shown in Figure 4 shown.

[0176] In particular, the above formula is a preliminary result derived from geometric relationships, which is aimed at the cases of θ>0 and φ>0. Other cases will be considered next.

[0177] For the case of θ<0, formula (4) should be modified as follows:

[0178]

[0179] For the case of φ<0, formulas (6) and (8) should be modified as follows:

[0180]

[0181] Formulas (5) and (7) remain unchanged.

[0182] Step 1.4: The information obtained so far can be used to infer the feasibility of satellite imaging. When the longitude offset angle between the target and the satellite trajectory on the TLL satisfies the relationship:

[0183] Δv1≥dv≥Δv2 (12);

[0184] The target will then be within the camera's FOV, meaning that this particular satellite qualifies as a potential imaging window. Therefore, the approximate time (or phase angle) when the satellite will cross the TLL within the relevant longitude range can be quickly determined without having to use the satellite's orbital extrapolation formula.

[0185] Step 1.5: Due to the assumption that the Earth is a perfect sphere and its orbit is a perfect circle, the rough search algorithm may ignore some passing satellites that would provide imaging windows. According to some studies, this only occurs when the target is at the edge of the camera's FOV. On the other hand, when using accurate Earth models and orbit information, there are still some passing satellites that meet the rough search process and potentially provide imaging windows, but the target is actually outside the FOV of these satellite cameras. To compensate for these shortcomings, when calculating the imaging strip boundary angle through the equation, a small margin (0.5-1 degree) is added to the FOV angle in formula (3). The error distance between the target and the camera optical axis projection point in the LH plane is then calculated in the final result, and those error distances that exceed the actual FOV range of the camera should be discarded.

[0186] Theory 3: During the coarse search, the relevant satellite transit time is determined in the form of TLL transit time. The default position of the camera optical axis is in a great circle perpendicular to the orbital plane and coinciding with the nadir direction. The imaging point should be located at Figure 4 The image is taken from point S in the image above, where the great circle connecting point S and the target (T) should be perpendicular to the orbital plane CD. Assuming the Earth is a perfect sphere, we can calculate the great circle angle CS from spherical trigonometry and obtain the imaging window, which is the satellite coordinates that meet the imaging conditions.

[0187] Step 2.1: First, the large circle angle TC can be determined from formula (9), noting that the longitude offset angle Δv is considered to be the actual declination angle dv:

[0188] cosp=sin 2 φ+cos 2 φ·cosdv(13);

[0189] Step 2.2: From spherical triangle BCD, the rotation angle Σ (between great circle CD and great circle CB) is:

[0190]

[0191] Step 2.3: Then in the spherical triangle NTC, the rotation angle Λ is:

[0192]

[0193] Step 2.4: Finally, using the spherical triangle CTS, we can determine the value m of the rotation angle π and the large fillet angle CS:

[0194] Π=π-Λ-Σ (16); m=tan -1 [tanp·cos(Λ+Σ)] (17);

[0195] Step 2.5: Then, the estimated imaging point S is determined as:

[0196] α S =α C +m (18);

[0197] Where, α S and α C They represent the phase angles of the satellite at points S and C respectively.

[0198] Step 2.6: The above algorithm is for 0≤α S <π / 2, Δv≥0 and i<90°, such as Figure 4 As shown. When the satellite phase angle is in other ranges, Δv < 0, and i > 90°, the above geometric relationship will change and the formula needs to be fine-tuned. Formulas (13-15) and (18) remain unchanged, while formulas (16) and (17) undergo slight changes. Table 1 lists the algorithm for calculating m for all possible cases.

[0199] Table 1m: General categories of calculation methods

[0200]

[0201]

[0202] Step 3.1: To optimize the estimation results, accurate satellite orbit extrapolation equations are required to obtain accurate satellite position coordinates and velocity information at any time. The following extrapolation equations use four redundant coordinates to describe the satellite orbit extrapolation, which is the basic model that the subsequent correction step modeling relies on:

[0203]

[0204] Where Ω is the right ascension of the ascending node of the satellite orbit, i is the orbit inclination angle, λ is the satellite latitude parameter, r is the radial distance of the satellite in the orbital coordinate system, a is the semi-major axis of the orbit, A is the amplitude of the orbital radius change, α is the phase angle, and α p is the true anomaly, n is the angular velocity, and the other parameters are the influencing factors of various factors. Their types and causes are given in Table 2.

[0205] Table 2 Constants in orbital extrapolation equations

[0206]

[0207]

[0208] For LEO satellites with orbital eccentricity not exceeding 0.005, the above orbit extrapolation equations give direct solutions to satellite positions. Satellite coordinates are single-valued functions of phase angle or time, and their accuracy is better than the results obtained using the SGP4 method.

[0209] Step 3.2: The satellite coordinates are given as a function of time in the orbital coordinate system (LO) according to the orbit extrapolation equations mentioned above, and then transformed into the Earth-centered inertial system (ECI), the Earth-fixed non-inertial system (ECEF), and the local tangent coordinate system (LT).

[0210] Specifically, the formula for converting LO coordinates to ECEF is as follows:

[0211]

[0212] Where, v is the ephemeris time;

[0213] The coordinate transformation matrix for converting ECEF coordinates to LT coordinates is as follows:

[0214]

[0215] Where λ is the target longitude and φ is the target latitude.

[0216] Step 3.3: Next, based on the given satellite camera information (camera tilt, etc.), calculate the value of the unit direction vector of the camera optical axis in the satellite body axis system, and then transform it to the LO, ECI, ECEF and LT coordinate systems in sequence. Therefore, the tilt vector from the target to the satellite is the satellite position vector in the LT coordinate system. From the azimuth vector of the camera optical axis, the position vector r of the camera optical axis projection point on the LH plane can be easily determined. L ,like Figure 1 shown.

[0217] Step 3.4: Velocity vector V of the camera optical axis projection point on the LH plane L It is also one of the evaluation conditions of the imaging condition of formula (1), and can be calculated by differentiating the orbit extrapolation equation and performing similar coordinate transformation.

[0218] In this scheme, a simpler method is adopted for calculation. First, a very small time interval δt is specified, and then two consecutive projection points of the camera optical axis on the LH plane before and after δt are found. Since the time difference is very small, the distance between the two projection points divided by δt can be used to obtain an average velocity vector. When the time step is set small enough, it can be considered as instantaneous velocity.

[0219] Step 3.5: The time error Δt from the current satellite position to the position that meets the imaging condition (1) can be calculated as:

[0220]

[0221] Step 3.6: Add this time error Δt to the estimated time and enter the orbit extrapolation formula to evolve the new satellite coordinates, and then generate the new satellite camera optical axis projection point in the LH plane. Apply the iterative process based on this algorithm until the time error Δt is reduced to the allowable tolerance range, and the iterative process is completed:

[0222]

[0223] Where Δt' is a given small time interval, such as 1 second or 2 seconds.

[0224] 2. For the case where the satellite has yaw attitude maneuverability:

[0225] Theory 4: When the satellite performs yaw maneuvers, the imaging coverage area of ​​the camera with an inclined angle is greatly increased. Due to the existence of FOV, the projection of the camera's complete imaging area on the LH plane becomes a large ring belt, such as Figure 6 Using simple geometric principles, we can consider that a straight line and a circle may have 1, 2, or 0 intersections, so the satellite may have one or two imaging opportunities in this operating mode.

[0226] Step 1: For a rough search, the limiting imaging angle of the Earth's surface imaging area is modified. Formula (3) should be adjusted accordingly and replaced by Formula (22). Here, since the satellite can rotate around its nadir axis, the formula presents a symmetrical form. After determining the limiting imaging angle, the subsequent search algorithm remains the same as in the case without attitude maneuvers:

[0227]

[0228] For the imaging estimation step, the algorithm discussed in the previous algorithm can be applied to this model. First, determine the imaging approximate point when the great circle connecting the satellite trajectory and the target is perpendicular to the orbital plane.

[0229] Since satellite imaging has more complex imaging conditions in the second working mode, the following three steps are further used to estimate the approximate imaging point.

[0230] Step 2.1: Point the camera to the target side. Since the yaw maneuver allows the satellite to rotate around its nadir axis, and thus the imaging angle of the camera changes continuously, the initial camera pointing can point to the target side or the opposite side of the target. Next, an algorithm needs to be constructed to derive the relative position relationship between the camera and the target. Here, an initial yaw angle Ψ0 (0 or π) is defined so that the satellite points the camera to the target side. First, the position vector r of the satellite trajectory in the LH plane is calculated using the method discussed aboveS =(r x ,r y ) and velocity vector V S =(V x ,V y ). Then give the judgment condition and use vector cross product as follows:

[0231] z=|r S ×V S |=r x V y -r y V x (twenty three);

[0232] According to the calculated sign of z and the camera tilt angle, the initial yaw angle Ψ0 can be determined according to Table 3.

[0233] Table 3 Determination of Ψ0 value

[0234]

[0235] Step 2.2: Determine how many imaging windows a given transiting satellite has. For a transiting satellite with imaging conditions, the number of imaging windows has three cases: 0, 1, or 2. Figure 6 As shown. The second step is the most critical point for determining the number of imaging times. Project the optical axis of the camera and its outer half of the FOV (with an inclination angle of τ+γ / 2) onto the LH plane to Figure 6 Then calculate the distances from the target T, the camera optical axis point C, and the external boundary limit point B to the satellite trajectory S, which are d, d1, and d2, respectively, as follows: Figure 6 As shown in Table 4, for the case where the target is within the circle drawn by the camera's optical axis, the satellite has two potential imaging points. In summary, the relationship between these three distance values ​​determines the number of potential imaging windows during this satellite transit, as listed in Table 4.

[0236] Table 4 Single imaging times of transit satellites

[0237]

[0238] Step 2.3: Estimate the approximate imaging point. The different situations given in the above table are discussed in a classified manner.

[0239] In Case 1, the target position is within the outer half of the camera's FOV. The camera's optical axis cannot accurately point toward the target, indicating imaging deviation. This situation is similar to Case 1 (where the satellite does not perform attitude maneuvers). As mentioned above, in Case 1, the camera's placement is already at the extreme imaging point, so yaw maneuvers are no longer necessary. The estimated point S can be directly incorporated into the orbit extrapolation formula for iterative correction.

[0240] For case 2, the target position is within the range of the camera's optical axis maneuvering, which corresponds to the situation where the circle and the line have two intersection points, so there are two imaging opportunities, and further yaw maneuvers are required to accurately align the camera's optical axis with the target. As an approximate estimate, it is assumed that the camera's optical axis projection on the LH plane during yaw maneuvers is a circle with its center at the satellite trajectory projection point and along its track line, as shown in Figure 6 Then, the approximate distance in the LH plane will be calculated by the following formula:

[0241]

[0242] Among them, V S is the velocity of the satellite trajectory on the LH plane, as shown in formula (23).

[0243] The satellite coordinates of the imaging window are determined by the time difference Δt, distance difference Δd and the relative position of point S in formula (25).

[0244] For case 3, the satellite was unable to image the target during this transit.

[0245] Step 3.1: Apply the orbit extrapolation equation to obtain the accurate imaging time. For the first two imaging opportunities listed in Table 4, the correction process is different.

[0246] For Case 1, the imaging conditions remain the same as in Equation (1), so the correction process and algorithm are identical to those described above. The position and velocity vector of the camera's optical axis projected onto the LH plane are used in the update step of the iterative process, using Equation (20). The required yaw angle in this case is the initial yaw angle in Table 3.

[0247] For Case 2, the situation becomes more complicated. The imaging condition requires that the camera's optical axis accurately point to the target. Since the camera tilt angle in this operating mode is fixed, the angle between the satellite vertical direction and the tilt vector from the satellite to the target should be equal to the camera's tilt angle, which can be described by the following equation:

[0248] -r·s=|r||s|cosτ (26);

[0249] Where r is the position vector of the satellite in ECEF, and s is the tilt vector in ECEF from the satellite to the target.

[0250] Theory 5: This formula is a newly defined imaging condition for this operating mode, not the imaging condition in formula (1). Since both r and s are functions of time, formula (26) is a nonlinear equation with respect to time. This equation has no analytical solution, or even if it does, it is very complex, so a numerical iterative algorithm is used here to solve the equation.

[0251] Step 3.2: From Figure 7 As can be seen in (a), the angle between r and s is first calculated as ξ, and then the angle between s and the camera direction ( Figure 7 The angle difference ∈ of SC in a):

[0252]

[0253] Step 3.3: The distance e1 from the camera optical axis projection point to the target and the distance e2 to the satellite trajectory in the LH plane can then be calculated:

[0254]

[0255] Where: σ is the elevation angle from the target to the satellite, which can be calculated by converting s from ECEF to LT and calculating its elevation angle relative to the LH plane.

[0256] Step 3.4: However, e1 is not the distance difference along the satellite orbit. Figure 7 As shown in (b), there is an angle ν between the ideal range difference direction and the satellite velocity direction, and this angle ν is not 90 degrees in this yaw maneuver mode. From the geometric relationship in the LH plane, the range error Δd parallel to the satellite motion direction can be calculated as:

[0257]

[0258] Step 3.5: The iterative time error Δt is obtained;

[0259]

[0260] Step 3.6: The above algorithm is used in the imaging correction process until the time error is reduced to the specified tolerance. Computational tests show that the above iterative algorithm converges very quickly, and only two or three iterations are needed to achieve millisecond-level accuracy. For the case where ∈ and e1 are negative, the above algorithm remains unchanged. However, for the case where ν>π / 2, that is, the second imaging window point, formula (30) needs to be modified to:

[0261]

[0262] Step 3.7: Through iteration, the tilt vector s gradually converges as Δt decreases until it approaches the camera optical axis. The final yaw pose required for image capture is the azimuth of the tilt vector s in the LO coordinate system. Therefore, by converting the tilt vector s from ECEF to ECI and LO, the second yaw angle can be calculated. Finally, this angle is added to the initial yaw angle listed in Table 3 to obtain the final yaw angle:

[0263] Ψ=Ψ'+Ψ0 (33).

[0264] Test example

[0265] Because satellite imaging algorithms are developed based on practical mission requirements, validating them using actual satellite imaging missions is essential. This test case uses the imaging mission of the Tsinghua-1 satellite, which entered orbit on June 28, 2000, for validation. Imaging window prediction is performed under the first through third operating modes defined in this paper. Specific validation steps are presented separately for the three algorithmic steps: coarse search, image estimation, and image correction.

[0266] The present invention uses programming to implement the algorithm and uses the MATLAB software package under the Windows 11 system for program package development to realize the imaging window prediction technology of the "Tsinghua-1" satellite, and finally realizes the function of outputting the appropriate satellite imaging window and satellite attitude angle.

[0267] This program package reads the NORAD orbit data file 'sat000026385.csv' of the Tsinghua-1 satellite. The applicant of the present invention applied to Dr. TSKelso of the CelesTrak platform to access the satellite orbit parameter TLE (Two-Line Element) file and csv format file for a total of 10 days from September 5, 2000 to September 15, 2000. Part of the content of the csv format file is given in Table 5:

[0268] Table 5 Partial content of the NORAD orbit data file 'sat000026385.csv' of the Tsinghua-1 satellite

[0269]

[0270] The camera data of the Tsinghua-1 satellite are summarized in Table 6:

[0271] Table 6 “Tsinghua-1” satellite camera information

[0272]

[0273] Calculation results demonstrate that the algorithm developed in this paper has excellent estimation accuracy and fast convergence. After a rough search and estimation algorithm, the imaging time estimated by the imaging method is within one minute of the actual time, and the imaging time after correction is within a few seconds of the actual time. The iterative algorithm used in the imaging correction step achieves time residuals at the millisecond level after only one or two iterations. The calculation results for each algorithm step are presented below.

[0274] (1) Rough search steps

[0275] The present invention uses the epoch time: 11:42:58:306176 on September 5, 2000 as the orbit input time of the algorithm, and uses 21:27:13.444704 on September 12, 2000 as the end time, that is, the 1007th to 1115th orbits of the "Tsinghua-1" satellite around the earth since it was put into orbit as the reference range for the rough search.

[0276] In the case of no attitude control of the satellite, two satellite passes with imaging potential are searched. Since the time when the satellite is at the imaging point is very close to the time when the satellite passes the target latitude (TLL), the present invention defines the time when the satellite passes the TLL as the result of the rough search process, which is given in Table 7.

[0277] Table 7 The time from the satellite passing the target latitude (TLL) when the satellite is at the imaging point in the first case

[0278]

[0279] For the case where the satellite is capable of yaw maneuvering, six satellite passes with imaging possibilities were found, and the time when the satellite passed the TLL is given in Table 8.

[0280] Table 8 The time from the satellite passing the target latitude (TLL) when the satellite is at the imaging point in the second case

[0281]

[0282] It can be seen that, under the premise that the satellite has a tilted camera, yaw maneuvers and roll maneuvers will provide the satellite with more imaging windows.

[0283] (2) Imaging estimation steps

[0284] Next, the present invention performs image estimation based on the roughly searched imaging window. In this step, the Earth is assumed to be a perfect sphere, and the satellite orbit is assumed to be a perfect circle. The algorithm relies entirely on spherical geometry. The imaging window obtained in the final image estimation step is determined by analytical methods and serves as the initial value for the imaging correction process.

[0285] For the first working mode, that is, the two imaging windows obtained when the satellite has no available attitude control, the algorithm gives the estimated value of the imaging window, which is given in Table 9.

[0286] Table 9 Satellite imaging window estimates for the first case

[0287]

[0288] For the second working mode, that is, the satellite has 6 imaging windows obtained when the yaw maneuver is performed, the algorithm gives the estimated value of the imaging window, which is given in Table 10.

[0289] Table 10 Satellite imaging window estimates for the second case

[0290]

[0291]

[0292] (3) Imaging correction steps

[0293] The core of the imaging correction algorithm is the ability to determine the specific coordinate position of the satellite at a specific moment. The orbit extrapolation equation described above is the key to solving this problem. To verify the feasibility and accuracy of the algorithm developed based on this orbit extrapolation equation, the epoch time (11:42:58:306176, September 5, 2000) was used as the algorithm's orbit input time, and the end time (21:51:53.535456, September 15, 2000) was used as the algorithm's orbit input time. This refers to the 1007th to 1159th orbit of the Tsinghua-1 satellite since its entry into orbit, spanning approximately 10 days. The present invention uses a total of 17 orbit information recording times as the actual satellite coordinate parameters as a control group. The satellite orbit parameters at the epoch time (11:42:58:306176, September 5, 2000) are substituted into the orbit extrapolation equation algorithm as the initial values. The satellite coordinate parameters at the 17 recording points are obtained as the calculation results, and the two sets of data are compared. The reference satellite coordinate parameters include the right ascension Ω of the orbital ascending node, the orbital inclination angle i, the satellite latitude parameter (equivalent phase angle) λ, and the radial distance r of the satellite in the orbital coordinate system is also given.

[0294] For the comparison of the right ascension Ω of the ascending node of the satellite orbit, Figure 8 The blue line in the figure is the calculated value, and the red dotted line is the actual value. It can be seen that the prediction of the satellite's orbital ascending node right ascension Ω is quite accurate. The error reaches the maximum at the 17th recording point (about the tenth day), which is 0.079°.

[0295] For the comparison of satellite orbit inclination angle i, Figure 9The blue line in the figure is the calculated value, and the red dotted line is the actual value. It can be seen that the error of the satellite's orbital inclination at most recording points is kept at 10 -4 The accuracy is on the order of deg. This is because most of the recording points are recorded near the ascending node of the satellite, which gives the highest calculation accuracy. At the 10th and 11th recording points, since both of these recording points are in the mid-latitudes of the Northern Hemisphere, the accuracy will decrease. The maximum error occurs at the 10th recording point, which is 8.2×10 -3 deg.

[0296] For the comparison of satellite latitude parameters (equivalent phase angle) λ, Figure 10 The blue line in the figure is the calculated value, and the red dotted line is the actual value. It can be seen that the prediction of the satellite latitude parameter (equivalent phase angle) λ is quite accurate, with the error reaching the maximum of 0.7966° at the 17th recording point (approximately the tenth day).

[0297] For the change of the radial distance r of the satellite in the orbital system, Figure 11 given.

[0298] Under the premise that the algorithm can ensure the accuracy of satellite coordinates, an iterative algorithm is used to give the satellite imaging window. The present invention will take a typical imaging window for each of the three working modes for verification.

[0299] For the first working mode, the present invention corrects the imaging window of 2000-09-08 12:21:37.265 in the imaging estimation step, and the final accurate value of the imaging window is 2000-09-08 12:21:40.298, as shown in Table 11. The residual in the iterative process is given by Figure 12 As can be seen from the figure, after two rounds of iterations, the time residual has reached 0.036 seconds, and the convergence is good.

[0300] Table 11 Correction results of the imaging window at 2000-09-08 12:21:37.265

[0301]

[0302] For the second working mode, the present invention corrects the imaging window of 2000-09-10 12:05:10.789 in the imaging estimation step. 2000-09-10 12:05:10.789 is the estimated time when the satellite orbit is perpendicular to the great circle connecting the sub-satellite point and the target. According to the algorithm analysis described above, d = 7.9×10 4 m, d1=2.6×10 5m, so d < d1. In this case, the satellite has two imaging windows. The estimated values ​​of the imaging windows are given as 2000-09-10 12:04:36.650 and 2000-09-10 12:05:51.929. The final accurate values ​​of the imaging windows are 2000-09-10 12:04:38.710 and 2000-09-10 12:05:48.870, as shown in Table 12. The residual in the iterative process is given by Figure 13 、 Figure 14 Present.

[0303] Table 12 Correction results of the imaging window at 2000-09-10 12:05:10.789

[0304]

[0305] Depend on Figure 13 、 14 As can be seen from the figure, after two iterations, the time residuals have reached 0.249 seconds and -0.34 seconds, indicating good convergence. The algorithm also gives the satellite's yaw angles in the imaging window as 74.59° and -74.43°, respectively, as shown in Table 13.

[0306] Table 13 Yaw angle of satellite imaging window

[0307]

[0308] The above description of the embodiments is intended to facilitate understanding and use of the invention by those skilled in the art. It will be apparent that those skilled in the art can readily make various modifications to these embodiments and apply the general principles described herein to other embodiments without requiring inventive effort. Therefore, the present invention is not limited to the above-described embodiments. Improvements and modifications made by those skilled in the art based on the disclosure of the present invention, without departing from the scope of the present invention, should be within the scope of protection of the present invention.

Claims

1. A method for predicting an imaging window of an imaging satellite for earth observation, used for imaging an imaging target located on the earth's surface, characterized in that: The steps include: Step 1: A rough search is performed using an analytical method, treating the Earth as a perfect sphere and the satellite orbit as a perfect circle, to search for a transiting satellite that meets the imaging window conditions. Step 2: Image estimation, using analytical methods, considering the Earth as an ideal sphere and the satellite orbit as an ideal circle, and estimating the satellite's image point based on the results of a rough search; Step 3: Imaging correction, using numerical algorithms to compensate for imaging estimation errors through orbit extrapolation models; When the satellite does not have the ability to maneuver or the imaging target is beyond the satellite's attitude maneuvering range: Steps 1 and 2: Use the satellite's camera optical axis pointing as the initial position. Based on the satellite camera's inclination and FOV, obtain the angular deviation through the great circle, then obtain the longitude deviation on the latitude TLL of the imaging target, and finally obtain the satellite imaging window estimate. Step 3: Calculate the precise coordinates of the satellite using an orbital extrapolation model. This orbital extrapolation model considers secular perturbations, long-period perturbations, and short-period perturbations caused by the Earth's oblateness, orbital eccentricity, and atmospheric drag. Substitute the estimated imaging window value into the model and use the iterative equation to solve the imaging window. The imaging window conditions are: L ·V L =0, where r L V is the position vector of the line connecting the imaging target and the projection point of the camera optical axis, L is the velocity vector of the projection point of the camera optical axis; When the satellite has attitude maneuverability: Step 1: Determine the limit imaging angle of the satellite camera and search for a transiting satellite that meets the imaging window conditions; Step 2: The great circle containing the imaging target and the satellite's projection point on the Earth is perpendicular to the satellite orbit plane. The parameters of the great circle projected on the ground by the imaging target relative to the camera's optical axis are obtained, the number of satellite imaging windows is obtained, and the approximate imaging point is estimated; Step 3: Apply the orbit extrapolation model to obtain the imaging window; The imaging window condition is: the moment when the camera optical axis coincides with the target tilt vector, and the target tilt vector is the position vector of the line connecting the satellite and the imaging target.

2. The method for predicting an imaging window for an imaging satellite to observe the earth according to claim 1, wherein: When the satellite does not have attitude maneuvering capability, it is modeled based on two-body motion: Step 1: When the longitude offset angle between the imaging target and the satellite trajectory at the TLL satisfies: Δv1 ≥ dv ≥ Δv2, the satellite has the imaging window condition; where Δv1 and Δv2 are the longitude offset angles at the TLL corresponding to the extreme imaging angles of the satellite camera, and dv is the longitude difference angle between the satellite trajectory and the imaging target after the satellite orbits the Earth N times; Step 2: Estimate the imaging point based on the satellite's phase angle α at the TLL, the initial longitude offset angle Δv between the satellite and the imaging target at the TLL at time t0, and the satellite's orbit inclination angle i; Step 3: Obtain the precise position coordinates of the satellite through the orbit extrapolation model, calculate the time error Δt from the current satellite position to the position required by the imaging window conditions, add Δt to the estimated time, and substitute it into the orbit extrapolation model again for iteration until Δt is within the allowable tolerance range. The iteration is completed.

3. The method for predicting an imaging window for an imaging satellite to observe the earth according to claim 2, wherein: In step 1: On the surface, the imaging target is set as T, the intersection of the satellite trajectory and the TLL is set as C, and the intersection of the satellite trajectory and the equatorial plane is set as D; DV satisfies: Where, is the Earth's rotation rate, n is the mean anomaly angular velocity, v s is the longitude of the satellite trajectory, v T is the longitude of the imaging target, and N is the number of satellite revolutions since time t0; The limiting imaging angles θ1 and θ2 satisfy: Where τ is the inclination angle of the satellite camera, γ is the FOV of the satellite camera, R is the radius of the earth, and h is the satellite altitude; Convert the limiting imaging angle to the longitude offset angles Δv1 and Δv2 at the TLL: in, cosp=coslcosη+sinlsinηcosΓ; where φ is the latitude of the imaging target, i is the satellite orbit inclination angle, p, η, and l are the large circle angles of TC, TD, and CD, respectively, and Γ is the rotation angle between TD and CD in the spherical triangle TCD.

4. The method for predicting an imaging window for earth observation by an imaging satellite according to claim 3, wherein: When calculating the limiting imaging angle, a margin of 0.5-1 degrees is added to the FOV of the satellite camera.

5. The method for predicting an imaging window for earth observation by an imaging satellite according to claim 2, wherein: In step 2: On the surface, the following settings are made: the imaging target is T, the intersection of the satellite trajectory and the TLL is C, the intersection of the satellite trajectory and the equatorial plane is D, the imaging point is S, the intersection of the meridian of point C and the equatorial plane is B, and the Earth's pole is N; The initial longitude offset angle Δv is used as the actual longitude difference angle dv; cosp=sin 2 φ+cos 2 φ·cosdv; The rotation angle Σ between CD and CB in spherical triangle BCD is: The rotation angle Λ between CN and CT in the spherical triangle NCT is: The rotation angle π between CT and CS in the spherical triangle TCS is: Π=π-Λ-Σ; The value m of the large fillet CS is: 1) When 0≤α S <π / 2 or π≤α S <3π / 2, and dv≥0, and 0<i<π / 2: m=tan -1 [tanp·cos(Λ+Σ)]; 2) When 0≤α S <π / 2 or π≤α S <3π / 2, and dv<0, and 0<i<π / 2: m=tan -1 [tanp·cos(Λ-Σ)]; 3) When π / 2 ≤ α S < π or 3π / 2 ≤ α S < 2π, and dv ≥ 0, and 0 < i < π / 2: m=-tan -1 [tanp·cos(Λ-Σ)]; 4) When π / 2≤α S <π or 3π / 2≤α S <2π, and dv<0, and 0<i<π / 2: m=-tan -1 [tanp·cos(Λ+Σ)]; 5) When 0 ≤ α S < π / 2 or π ≤ α S < 3π / 2, and dv ≥ 0, and π / 2 < i < π: m=tan -1 [tanp·cos(Λ-Σ)]; 6) When 0≤α S <π / 2 or π≤α S When <3π / 2, and dv<0, and π / 2<i<π: m=tan -1 [tanp·cos(Λ+Σ)]; 7) When π / 2 ≤ α S < π or 3π / 2 ≤ α S < 2π, and dv ≥ 0, and π / 2 < i < π: m=-tan -1 [tanp·cos(Λ+Σ)]; 8) When π / 2 ≤ α S < π or 3π / 2 ≤ α S < 2π, and dv < 0, and π / 2 < i < π: m=-tan -1 [tanp·cos(Λ-Σ)]; Where φ is the latitude of the imaging target, i is the satellite orbit inclination angle, p and l are the large circle angles of TC and CD respectively, and α S is the phase angle of the satellite at point S; The imaging point S satisfies: α S ＝α C +m； Where, α C are the phase angles of the satellite at point C respectively.

6. The method for predicting an imaging window for an imaging satellite to observe the earth according to claim 2, wherein: In step 3: The orbital extrapolation model is: Where a is the semi-major axis of the orbit, ρ is the radius variation parameter, A is the amplitude of the orbit radius variation, α is the phase angle, and α p is the true anomaly, κ is the secular offset constant, θ is the orbital precession constant, χ is the long-period perturbation constant, Δ is the short-period perturbation constant, B is the atmospheric drag factor, Ω is the right ascension of the ascending node of the satellite orbit, i is the orbital inclination angle, λ is the satellite latitude parameter, r is the radial distance of the satellite in the orbital coordinate system, and n is the angular velocity; The satellite coordinates obtained by the orbit extrapolation model are time functions based on the orbital coordinate system LO. The obtained satellite coordinates are sequentially transformed into the Earth-centered inertial system ECI, the Earth-fixed non-inertial system ECEF, and the local tangent coordinate system LT of the imaging target. Calculate the unit direction vector of the camera optical axis according to the information of the satellite camera, and then transform it into the LO, ECI, ECEF and LT coordinate systems in turn to obtain r L and V L ; Among them, V L Calculated according to the following formula: Where Δt' is a given small time interval; Calculate the time error Δt from the current satellite position to the position required by the imaging window conditions: The time error Δt is added to the estimated time and substituted into the orbit extrapolation model again for iteration.

7. The method for predicting an imaging window for earth observation by an imaging satellite according to claim 1, wherein: When the satellite has the ability to maneuver, the imaging area of ​​the satellite camera is a large circular belt: Step 1: When the longitude offset angle between the imaging target and the satellite trajectory at the TLL satisfies: Δv1 ≥ dv ≥ Δv2, the satellite has the imaging window condition; where Δv1 and Δv2 are the longitude offset angles at the TLL corresponding to the extreme imaging angles of the satellite camera, and dv is the longitude difference angle between the satellite trajectory and the imaging target after the satellite orbits the Earth N times; According to the relationship between the position of the imaging target and the large circular ring, it is divided into 0, 1 and 2 imaging windows; When the imaging window is 0 times, the imaging target cannot be imaged; When the imaging window is 1 time, Step 2: Obtain the initial yaw angle Ψ0 so that the satellite points the camera toward the imaging target. Estimate the imaging point based on the satellite's phase angle α at the TLL, the initial longitude offset angle Δv between the satellite and the imaging target at time t0, and the satellite's orbit inclination angle i. Step 3: Obtain the precise position coordinates of the satellite through the orbit extrapolation model, calculate the time error Δt from the current satellite position to the position required by the imaging window conditions, add Δt to the estimated time and substitute it into the orbit extrapolation model again to iterate until Δt is within the allowable tolerance range, and the iteration is completed; When the imaging window is 2 times, Step 2: Obtain the initial yaw angle Ψ0 so that the satellite points the camera toward the imaging target. Estimate the imaging point based on the distance d1 from the projection point of the camera optical axis to the satellite trajectory and the distance d from the imaging target to the satellite trajectory. Step 3: Define the imaging condition, -r·s=|r||s|cosτ, where r is the position vector of the satellite in the Earth-fixed non-inertial system ECEF, and s is the tilt vector from the satellite to the imaging target in the Earth-fixed non-inertial system ECEF. Both r and s are functions of time. Calculate the distance e1 from the camera optical axis projection point to the imaging target and the distance e2 from the camera optical axis projection point to the satellite trajectory in the LH plane. Then calculate the distance error Δd in the LH plane and obtain the time error Δt. Substitute the time error Δt into the imaging condition and iterate until the tilt vector s is close to the camera optical axis. Then convert the tilt vector s to the Earth-centered inertial system ECI and the orbital coordinate system LO to obtain the yaw angle, and finally obtain the total yaw angle Ψ.

8. The method for predicting an imaging window for earth observation by an imaging satellite according to claim 7, wherein: In step 1: Set the imaging target as T, the intersection of the satellite trajectory and TLL as C, and the intersection of the satellite trajectory and the equatorial plane as D; DV satisfies: Where, is the Earth's rotation rate, n is the mean anomaly angular velocity, v s is the longitude of the satellite trajectory, v T is the longitude of the imaging target, and N is the number of satellite revolutions since time t0; The limiting imaging angles θ1 and θ2 satisfy: Where τ is the inclination angle of the satellite camera, γ is the FOV of the satellite camera, R is the radius of the earth, and h is the satellite altitude; Convert the limiting imaging angle to the longitude offset angles Δv1 and Δv2 at the TLL: in, cosp=coslcosη+sinlsinηcosΓ; where φ is the latitude of the imaging target, i is the satellite orbit inclination angle, p, η, and l are the large circle angles of TC, TD, and CD, respectively, and Γ is the rotation angle between TD and CD in the spherical triangle TCD.

9. The method for predicting an imaging window for earth observation by an imaging satellite according to claim 7, wherein: In step 2: Get the position vector r of the line connecting the imaging target and the subsatellite point S and the velocity vector V of the subsatellite point S And the judgment condition z is calculated by vector cross product: z=|r S ×V S |=r x V y -r y V x ; According to the calculated sign of z and the satellite camera inclination angle τ, the initial yaw angle Ψ0 is determined: When τ and z have the same sign, then Ψ0 = 0; When τ and z have opposite signs, then Ψ0 = π.

10. The method for predicting an imaging window for an imaging satellite to observe the earth according to claim 7, wherein: In step 3: When the imaging window is 2 times, Calculate the angle ξ between r and s, and then calculate the angle difference ∈ based on s and the direction of the satellite camera: Where, σ is the elevation angle from the imaging target to the satellite; Based on the angle ν between the distance from the ideal camera optical axis projection point to the imaging target and the satellite velocity, Δd is calculated: When ν<π / 2: When ν>π / 2: Then calculate Δt: Total yaw angle Ψ: Ψ=Ψ'+Ψ0.

Images