Method for predicting earth observation imaging window of imaging satellite

Abstract

The invention belongs to the field of satellite imaging algorithms, and particularly relates to a method for predicting an earth observation imaging window of an imaging satellite, which is used for imaging an imaging target and comprises the following steps of: 1, roughly searching; step 2, performing imaging estimation; and step 3, performing imaging correction. Compared with the prior art, the imaging satellite task planning method solves the problems that in the prior art, imaging satellite task planning takes little consideration to constraint conditions, calculation errors are large, algorithms are complex, and the imaging effect is poor, and perturbation and errors caused by the factors such as the earth oblateness, the orbit eccentricity rate and the atmospheric resistance are considered in the scheme. And a simple and rapid iterative algorithm is further applied to solve a satellite imaging window, so that the efficiency and the accuracy of an autonomous task are improved.

Description

Technical Field

[0001] The 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 working satellites, earth observation satellites are particularly common. Earth Observation Satellite (EOS) is a type of satellite used to monitor and study the surface conditions and atmospheric environment of the earth. This type of satellite can collect different types of data through various sensors on board, including visible spectrum, infrared spectrum, electromagnetic wave band, etc. Nowadays, various earth observation satellites have played a very important role in scientific research, military, reconnaissance and other fields. In scientific research, these satellites can help scientists study the earth's climate change, environmental change, ecosystem and natural disasters. In military reconnaissance, they can provide high-resolution surface images and other data to help monitor and analyze military operations. In agriculture, earth observation satellites can monitor crop growth, land use and water resource management. Due to the complex types and large number of these satellites, they play a key role in various fields. Different types of sensors on the satellite can provide different types of data to meet various mission requirements. Some earth observation imaging satellites have the function of disaster monitoring, such as detecting the frequency of forest fires. Carbon dioxide emissions generated by global forest fires now account for one-fifth of global carbon dioxide emissions, so using earth observation satellites to monitor forest fires is one of the most important tasks.

[0003] To carry out satellite earth observation activities, it is necessary to plan the mission of the satellite, in which resource allocation, satellite orbit information, ground station resources and other conditions are constraints. Under these constraints, the satellite mission planning is solved, and then the best mission execution plan is formulated. The above process is actually a complex combinatorial optimization problem, involving the scheduling of multiple constraints and resources, as well as multiple tasks, multiple time windows and multiple optimization goals. Looking at the development history of satellite mission planning, ground base stations have played a very important role. After most satellites collect information, they need to immediately transmit the information to the ground station. After receiving the satellite information, the ground station will upload it to the space base, that is, the satellite platform, based on the satellite coordinates, satellite orbit information, satellite platform payload and other information, plus the already parsed mission planning constraints, together with the pre-planned satellite mission instructions. The satellite platform needs to complete the task step by step according to the uploaded instructions, and then transmit the collected mission information back to the ground base station again. This is the traditional satellite mission planning operation mode. From the perspective of continuous information transmission, this is a large loop mode of sky (satellite)-ground (ground station). It is not difficult to see that the traditional approach consumes both time and human resources, which results in low efficiency of satellite missions, inability to collect sufficient information within a limited time, and almost no ability to respond to emergencies. In addition, the traditional solution is also a passive task scheduling process, which cannot obtain satellite information and resource status in real time, nor can it temporarily adjust the observation plan, which ultimately leads to a waste of satellite resources.

[0004] In fact, mankind has never stopped pursuing satellite efficiency and satellite resource utilization. For increasingly complex space missions, most of the command and control instructions of satellite systems, such as attitude maneuvers, data downlink, orbital maneuvers, resource scheduling, etc., will no longer rely entirely on ground issuance. However, for previous earth observation satellites, enabling satellites to have fully automatic autonomous mission planning was beyond the technological level at the time, and early satellites had a large size and a huge number of sensors, and the cost of satellite trial and error was huge. Therefore, engineers have been using improved solutions based on the traditional scheduling solution that relies entirely on ground stations for a long time, improving the integration of mission instructions and optimizing the mission planning algorithm, reducing the number of calculations of related equations such as satellite orbits to achieve the purpose of improving efficiency.

[0005] In recent years, with the rapid development of micro-satellites, the number of satellites in orbit has increased. These satellites usually need to form constellations to carry out space mission activities, and more and more users have put forward higher requirements for the complexity and diversification of satellite observation tasks. Under this situation, if satellites still rely entirely on ground base stations for observation tasks, the ground base stations will become overwhelmed. At the same time, the types of remote sensors have become more diverse, and the growing demand for user information services has brought great prospects and opportunities for satellite observation missions, but also brought unprecedented challenges. It is mainly reflected in the following points:

[0006] 1) Observation task requests are sudden, and it is necessary to reasonably use the relationship between the existing planning scheme and the 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 prediction method for an imaging satellite's earth observation imaging window in order to solve at least one of the above problems, so as to solve the problems in the prior art that imaging satellite mission planning has little consideration of constraints, large calculation errors, complex algorithms, and poor imaging effects. This scheme takes into account the perturbations and errors caused by factors such as the earth's flattening, 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 of an imaging satellite for earth observation, used for imaging an imaging target located on the surface of the earth, comprises the following steps:

[0012] Step 1: Rough search, using analytical methods, considering the earth as an ideal sphere and the satellite orbit as an ideal circle, searching for transit satellites that meet the imaging window conditions;

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

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

[0015] When the satellite does not have the attitude maneuvering capability or the imaging target is beyond the satellite attitude maneuvering range: Step 1 and Step 2, take the satellite's camera optical axis pointing as the initial position, obtain the angle deviation through the great circle according to the satellite camera inclination and FOV, and 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 through the orbit extrapolation model, the orbit extrapolation model considers the secular perturbation, long-period perturbation and short-period perturbation caused by the earth's oblateness, orbital eccentricity and atmospheric drag, substitute the imaging window estimate 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 the ability of attitude maneuvering: Step 1, determine the limit imaging angle of the satellite camera, and search for a transit satellite with 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, obtain the parameters of the great circle of the imaging target relative to the camera's optical axis projected on the ground, obtain the number of satellite imaging windows and estimate the approximate imaging point; Step 3, apply the orbit extrapolation model to obtain the imaging window;

[0018] The condition of the imaging window 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 the ability to maneuver, a two-body motion model is used:

[0020] Step 1: When the longitude offset angles of the imaging target and the satellite trajectory at TLL satisfy: Δv1≥dv≥Δv2, the satellite has the imaging window condition; where Δv1 and Δv2 are the longitude offset angles at 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 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 and 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 meets:

[0026]

[0027] In the formula, is the Earth's rotation rate, n is the angular velocity of the mean anomaly, 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 θ1 and θ2 respectively; φ 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, Γ is the rotation angle between TD and CD in the spherical triangle TCD; sign(θ) function is to take the sign of θ.

[0036] Preferably, when calculating the limiting imaging angle, a margin of 0.5-1 degree 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 deviation angle Δv is taken as the actual longitude difference angle dv;

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

[0041] The rotation angle Σ between CD and CB in the 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π, and 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 circular 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] In the formula, α 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 successively 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 that transforms ECEF coordinates to LT coordinates is as follows:

[0078]

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

[0080] The unit direction vector of the camera optical axis is calculated based on the information of the satellite camera, and then transformed into the LO, ECI, ECEF and LT coordinate systems in turn to obtain r L and V L ;

[0081] Among them, V L Calculate 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 angles of the imaging target and the satellite trajectory at TLL satisfy: Δv1≥dv≥Δv2, the satellite has the imaging window condition; where Δv1 and Δv2 are the longitude offset angles at 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 to one side of 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 at the TLL, and the satellite 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 for iteration 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 to the side of 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, s is the tilt vector from the satellite to the imaging target in the earth-fixed non-inertial system ECEF, and 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 then 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, 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 meets:

[0100]

[0101] In the formula, is the Earth's rotation rate, n is the angular velocity of the mean anomaly, 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 task planning solution 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, the constraints such as the camera inclination angle, FOV, yaw and roll attitude maneuvers are considered, and the analytical modeling of spatial geometry and the numerical results of the iterative algorithm are used to simplify the task 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, the perturbations and errors caused by factors such as the earth's flattening, orbital eccentricity, and atmospheric drag are considered and compensated, and the prediction results are optimized by relying on a more refined earth shape and accurate orbital extrapolation equations.

[0132] Many traditional studies on imaging satellite mission planning problems pay little attention to constraints. When designing the algorithm, the present invention considers constraints such as the inclination angle of the camera installed on the satellite, the availability of satellite attitude maneuvers, and the number of imaging times of a single satellite transit 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 the analytical method of solid geometry to solve the initial values, and then uses the iterative method 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] The algorithm was evaluated by computational experiments based on actual satellite imaging mission scenarios. The results show that the results of the satellite imaging algorithm can achieve the expected goals and can obtain the optimal solution of the imaging window in each working mode. This algorithm has relatively simple calculation steps while ensuring accuracy. It can be attempted to perform local operations on a microsatellite platform to further verify its performance. 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 It is an auxiliary schematic diagram of spherical triangulation solution in the 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) in the embodiment on the LH plane as it rotates around its zenith axis.

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

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

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

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

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

[0146] Fig.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] Fig.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] Fig.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 is described in detail below with reference to the accompanying drawings and specific embodiments. All matters not covered in the present invention are well-known technologies.

[0150] Example

[0151] Definition 1: This scheme generally defines the satellite camera imaging window as two cases: the first case is when the optical axis of the satellite camera can accurately point to the specified 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 maneuvers are available;

[0154] Premise 2: In the second case, the satellite has no available attitude maneuvers or the target's position is outside the satellite's attitude maneuvering tolerance, for example, if 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 when the optical axis of the camera on the satellite coincides with the target tilt vector. The target tilt vector is defined as the vector connecting the satellite and the ground target in this scheme. A more intuitive statement is that the optical axis of the satellite camera can directly face the target, and the deviation of its projection point on the ground from the target is zero.

[0156] Definition 3: In the second case, the solution of the imaging window becomes complicated. When the satellite flies along its 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 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, so 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 given longitude and latitude on the surface of the earth, the x-axis points to the east, the y-axis points to the north, and the axis is given by the right-hand system, pointing vertically to the direction of 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 has no attitude maneuvers:

[0161] Assumption 1: Most satellite orbit solutions can be approximated using two-body motion. In the rough search process, this scheme 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 through the same point exactly in an inertial coordinate system (i.e., a coordinate system that does not follow the rotation of the Earth), such as Figure 2As the Earth rotates along its axis, the target will move along the target latitude line (TLL) mentioned above, which is defined as a small circle of constant latitude that passes through the target location. Therefore, in inertial space, each time the satellite re-passes the TLL, the Earth's rotation will bring the target closer to the longitude of the satellite on 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 length 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 over the 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 angle 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. As the satellite flies along its orbit, 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 inclination angle, γ is the camera FOV, R is the radius of the earth, and h is the satellite altitude.

[0170] Assumption 2: Since R and r (the radial distance of the satellite in the orbital coordinate system) change as the satellite moves along its orbit, these two angles are not constant with respect to time. However, for low-Earth orbit (LEO) Earth observation satellites, their orbits are usually near-circular polar orbits. For simplicity, the orbit semi-major axis and the Earth's equatorial radius are used in the previous equations to obtain a constant deviation angle. In this way, 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 to 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 fillet 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, targeting the cases of θ>0 and φ>0, and 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 currently obtained can already infer the feasibility of satellite imaging, when the longitude offset angle of 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, which means that this particular satellite qualifies as a potential imaging window. Therefore, the approximate time (or phase angle) when the satellite crosses the TLL within the relevant longitude range can be quickly determined without involving the satellite’s orbital extrapolation formula.

[0185] Step 1.5: Due to the assumption that the earth is a perfect sphere and the orbit is a perfect circle, the rough search algorithm may ignore some transit satellites that would provide imaging windows. According to some studies, this only occurs when the target is at the edge of the camera FOV. On the other hand, when using accurate earth models and orbit information, there are still some transit satellites that meet the rough search process and potentially provide imaging windows, but the targets are actually outside the FOV of these satellite cameras. To compensate for these shortcomings, a small margin (0.5-1 degree) is added to the FOV angle in formula (3) when calculating the imaging strip boundary angle through the equation. Then the error distance between the target and the camera optical axis projection point in the LH plane is 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 transits are determined in the form of TLL transit times. 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 Near point S in , where the great circle connecting point S and the target (T) should be perpendicular to the orbital plane CD. Again assuming that the earth is a perfect sphere, we can calculate the great circle angle CS from spherical trigonometry and get the imaging window, that 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), and it is noted that the longitude deviation angle Δv is regarded as 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, the spherical triangle CTS can be used to 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] In the formula, α 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) change slightly. Table 1 lists the algorithm for calculating m for all possible cases.

[0199] Table 1m Calculation methods

[0200]

[0201]

[0202] Step 3.1: In order 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 mainly relied on in the following correction steps:

[0203]

[0204] In the formula, Ω is the right ascension of the satellite orbit ascending node, i is the orbit inclination angle, λ is the satellite latitude parameter, r is the radial distance of the satellite in the orbit coordinate system, a is the semi-major axis of the orbit, A is the amplitude of the orbit radius change, α is the phase angle, and α p is the true anomaly angle, 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 orbital extrapolation equations give a direct solution to the satellite position. The 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 to the Earth-centered inertial system (ECI), the Earth-fixed non-inertial system (ECEF) and the local tangent coordinate system (LT) in sequence.

[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 that transforms 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 inclination, 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 turn. Therefore, the inclination 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: The 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 used 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 get 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 bring it into the orbit extrapolation formula to evolve the new satellite coordinates, and then generate a new satellite camera optical axis projection point in the LH plane. Apply an 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] Wherein, Δ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 maneuvers:

[0225] Theory 4: When the satellite has a yaw maneuver, 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, 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 imaging area on the Earth's surface is modified, and 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 mode. First, the imaging approximate point when the great circle connecting the satellite trajectory and the target is perpendicular to the orbital plane is determined.

[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 points.

[0230] Step 2.1: Point the camera to the side of the target. Since the yaw maneuver allows the satellite to rotate around its nadir axis, and the imaging angle of the camera changes continuously, the initial camera can be pointed to the target side or the opposite side of the target. Next, an algorithm needs to be built 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 the 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 inclination, the initial yaw angle Ψ0 can be determined according to Table 3.

[0233] Table 3 Determination of Ψ0 value

[0234]

[0235] Step 2.2: Find out how many imaging windows a given transit satellite has. For a transit satellite with imaging conditions, the number of imaging windows has three cases, 0, 1 or 2, such as Figure 6 As shown in the figure, 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. 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 optical axis, the satellite has two potential imaging points. In summary, the relationship between the three distance values ​​determines how many potential imaging windows there are 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. Different situations are discussed in the above table.

[0239] For case 1, the target position is within the outer half of the camera's FOV. The camera's optical axis cannot accurately point to the target, which means there is an imaging deviation. This situation is similar to the first case (the case where the satellite has no attitude maneuvers). As mentioned above, the camera's placement in the first case is already the extreme imaging point, so yaw maneuvers are no longer required. The estimated point S can be directly brought 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 intersections, 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 accompanied by yaw maneuvers is a circle with its center at the satellite trajectory projection point and along its track line, as 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 will be 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 equations to obtain accurate imaging times. For the first two imaging opportunities listed in Table 4, the correction process is different.

[0246] For case 1, the imaging condition is still formula (1), so the correction process and algorithm are the same as the previous algorithm. The position and velocity vector of the camera optical axis projection on the LH plane are used in the update step of the iterative process, and the algorithm is based on formula (20). The yaw angle required 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 optical axis of the camera is accurately pointed at the target. Since the camera tilt angle in this working 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 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 in this working mode, not the imaging condition in formula (1). Since r and s are both functions of time, formula (26) is a nonlinear equation about time. There is no analytical solution for this equation, or even if there is, it is very complicated, 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: Then 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 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 distance difference direction and the satellite velocity direction, and this angle ν is not 90 degrees in this yaw maneuvering mode. From the geometric relationship in the LH plane, the distance error Δd parallel to the satellite motion direction can be calculated as:

[0257]

[0258] Step 3.5: Then the iteration 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 allowable range. Computational tests show that the above iterative algorithm converges very quickly and only two or three iterations are required to achieve millisecond-level accuracy. For the case where ∈ and e1 are negative, the above algorithm remains unchanged, while for the case where ν>π / 2, that is, the second imaging window point, formula (30) needs to be modified as follows:

[0261]

[0262] Step 3.7: Through the iteration results, the tilt vector s will gradually converge as Δt becomes smaller until it is close to the camera optical axis. The final yaw attitude 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 value of 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] Since the satellite imaging algorithm is proposed based on actual mission requirements, it is very necessary to use the actual satellite imaging mission to verify the algorithm. This test example will use the imaging mission of the "Tsinghua No. 1" (TZINGHUA1) satellite that entered orbit on June 28, 2000 for verification, and perform imaging window prediction on it in the first to third working modes defined in the present invention. The specific verification content is presented according to the three algorithm steps of rough search, imaging estimation and imaging correction.

[0266] The present invention implements the algorithm by programming, uses the MATLAB software package under the Windows 11 system to develop the program package, realizes 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] The program package reads the NORAD orbit data file 'sat000026385.csv' of the Tsinghua-1 satellite. The applicant of the present invention applies to Dr.TSKelso of the CelesTrak platform to call the satellite orbit parameter TLE (Two-Line Element) file and csv format file for a total of 10 days from 2000-09-05 to 2000-09-15, wherein part of the csv format file is given in Table 5:

[0268] Table 5 Partial content of 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] The calculation results show that the algorithm developed by the present invention has good estimation accuracy and fast convergence. After the rough search and estimation algorithm, the imaging time estimated by the imaging is within 1 minute of the actual time, and the imaging time after the imaging correction is within a few seconds of the actual time. In the iterative algorithm used in the imaging correction step, the time residual reaches the millisecond level after only 1-2 rounds of iterations. The calculation results for each algorithm step are presented below.

[0274] (I) 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 where the satellite has no attitude control, two satellite passes with imaging possibilities 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] In the case where the satellite is capable of yaw maneuvers, 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] (II) Imaging estimation steps

[0284] Next, the present invention performs imaging estimation based on the roughly searched imaging window. In this step, the earth is considered as an ideal sphere, the satellite orbit is considered as an ideal circle, and the establishment of the algorithm is completely based on spherical geometry. The imaging window obtained in the final imaging estimation step is given by an analytical method as the initial value of 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] (III) Imaging correction steps

[0293] For the algorithm of the imaging correction process, the core content is to be able to obtain the specific coordinate position of the satellite at a certain moment, and the orbit extrapolation equation described above is the key to solving this problem. In order to verify the feasibility and accuracy of the algorithm developed based on this orbit extrapolation equation, the epoch time (Epoch time): 11:42:58:306176 on September 5, 2000 is used as the orbit input time of the algorithm, and 21:51:53.535456 on September 15, 2000 is used as the termination time, that is, the 1007th to 1159th laps of the "Tsinghua No. 1" satellite orbiting the earth since entering orbit, with a span of about 10 days. The present invention uses a total of 17 orbit information recording times as actual satellite coordinate parameters as a control group; the satellite orbit parameters of the epoch time 11:42:58:306176 on September 5, 2000 are substituted into the orbit extrapolation equation algorithm as initial values, and the satellite coordinate parameters at the 17 recording points are obtained as the calculation results, and the two groups of data are compared. The reference satellite coordinate parameters include the orbit ascending node right ascension Ω, the orbit 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 satellite orbit ascending node, 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, and 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, Fig. 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 of the recording points is maintained at 10 -4 The error is in the deg order of magnitude. This is because most of the recording points are recorded near the satellite ascending node, which results in the highest calculation accuracy. At the 10th and 11th recording points, since both 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) λ, Fig.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, and the error reaches the maximum at the 17th recording point (about the tenth day), which is 0.7966°.

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

[0298] Under the premise that the algorithm can ensure the accuracy of satellite coordinates, the 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 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 Fig.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 2000-09-08 12:21:37.265 Imaging window correction results

[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. Finally, the precise values ​​of the imaging windows are 2000-09-1012: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 Fig.13 , Fig.14 Presentation.

[0303] Table 12 2000-09-10 12:05:10.789 Imaging window correction results

[0304]

[0305] Depend on Fig.13 , 14 It can be seen that after two rounds of iterations, the time residual has reached 0.249 seconds and -0.34 seconds, with good convergence. The algorithm also gives the satellite's yaw angle in the imaging window, which is 74.59° and -74.43°, respectively, as shown in Table 13.

[0306] Table 13 Yaw angle during satellite imaging window

[0307]

[0308] The above description of the embodiments is to facilitate the understanding and use of the invention by those skilled in the art. It is obvious that those skilled in the art can easily make various modifications to these embodiments and apply the general principles described herein to other embodiments without creative work. Therefore, the present invention is not limited to the above embodiments, and 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 surface of the earth, characterized in that: The steps include: Step 1: Rough search, using analytical methods, considering the earth as an ideal sphere and the satellite orbit as an ideal circle, searching for transit satellites that meet the imaging window conditions; Step 2: Imaging estimation, using analytical methods, considering the earth as an ideal sphere and the satellite orbit as an ideal circle, and estimating the satellite's imaging point based on the results of a rough search; Step 3: Imaging correction, using numerical algorithms to compensate for the imaging estimation errors through the orbit extrapolation model; When the satellite does not have the attitude maneuvering capability or the imaging target is beyond the satellite attitude maneuvering range: Step 1 and Step 2, take the satellite's camera optical axis pointing as the initial position, obtain the angle deviation through the great circle according to the satellite camera inclination and FOV, and 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 through the orbit extrapolation model, the orbit extrapolation model considers the secular perturbation, long-period perturbation and short-period perturbation caused by the earth's oblateness, orbital eccentricity and atmospheric drag, substitute the imaging window estimate and use the iterative equation to solve the imaging window; 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; When the satellite has the ability of attitude maneuvering: Step 1, determine the limit imaging angle of the satellite camera, and search for a transit satellite with 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, obtain the parameters of the great circle of the imaging target relative to the camera's optical axis projected on the ground, obtain the number of satellite imaging windows and estimate the approximate imaging point; Step 3, apply the orbit extrapolation model to obtain the imaging window; The condition of the imaging window 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 of an imaging satellite for earth observation according to claim 1, characterized in that: When the satellite does not have the ability to maneuver, it is modeled based on two-body motion: Step 1: When the longitude offset angles of the imaging target and the satellite trajectory at TLL satisfy: Δv1≥dv≥Δv2, the satellite has the imaging window condition; where Δv1 and Δv2 are the longitude offset angles at 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 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 and the iteration is completed.

3. The method for predicting an imaging window of an imaging satellite for earth observation according to claim 2, characterized in that: 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 meets: In the formula, is the Earth's rotation rate, n is the angular velocity of the mean anomaly, 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 of an imaging satellite for earth observation according to claim 3, characterized in that: When calculating the limiting imaging angle, a margin of 0.5-1 degree is added to the FOV of the satellite camera.

5. The method for predicting an imaging window of an imaging satellite for earth observation according to claim 2, characterized in that: 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 deviation angle Δv is taken as the actual longitude difference angle dv; cosp=sin 2 φ+cos 2 φ·cosdv; The rotation angle Σ between CD and CB in the 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 circular 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； In the formula, α C are the phase angles of the satellite at point C respectively.

6. The method for predicting an imaging window of an imaging satellite for earth observation according to claim 2, characterized in that: 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 successively 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; The unit direction vector of the camera optical axis is calculated based on the information of the satellite camera, and then transformed into the LO, ECI, ECEF and LT coordinate systems in turn to obtain r L and V L ; Among them, V L Calculate 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 of an imaging satellite for earth observation according to claim 1, characterized in that: 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 angles of the imaging target and the satellite trajectory at TLL satisfy: Δv1≥dv≥Δv2, the satellite has the imaging window condition; where Δv1 and Δv2 are the longitude offset angles at 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 to one side of 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 at the TLL, and the satellite 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, 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 to the side of 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, s is the tilt vector from the satellite to the imaging target in the earth-fixed non-inertial system ECEF, and 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 then 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, obtain the yaw angle, and finally obtain the total yaw angle Ψ.

8. The method for predicting an imaging window of an imaging satellite for earth observation according to claim 7, characterized in that: 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 meets: In the formula, is the Earth's rotation rate, n is the angular velocity of the mean anomaly, 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 of an imaging satellite for earth observation according to claim 7, characterized in that: 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 of an imaging satellite for earth observation according to claim 7, characterized in that: 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