A method for determining payload coverage strip considering satellite attitude maneuvers

Abstract

The present invention discloses a method for determining a payload coverage strip taking into account satellite attitude maneuvers, comprising the following steps: discretizing the field boundary of an imaging payload on a satellite into a number of payload view elements, and calculating the position of the projection point of each payload view element on the earth's spherical surface; calculating the speed of each payload view element projection point based on the speed change caused by the satellite's circular motion around the earth and the speed change caused by the satellite's attitude maneuvers; calculating an index function for each payload view element projection point based on the position and speed of each payload view element projection point, and screening out projection points that can serve as coverage area envelope points based on the index function, thereby obtaining a payload coverage strip. The present invention relates to the field of space-based observation technology, has good effectiveness and efficiency, and can be effectively applied to the payload coverage range calculation requirements of ultra-agile earth observation satellites.

Description

Technical Field

[0001] The present invention relates to the field of space-based observation technology, and in particular to a method for determining a payload coverage strip taking satellite attitude maneuvers into consideration. Background Art

[0002] Large-scale low-orbit satellite constellations, characterized by wide coverage, high throughput, and low latency, have attracted widespread attention and are rapidly developing in applications such as communications, navigation, and remote sensing. The rapid growth of constellations poses significant challenges to the rapid and accurate analysis of payload coverage. Ground coverage refers to the satellite's coverage of ground targets over a specific period of time or at a specific moment. The actual coverage of the payload must be considered during satellite operations. In the field of Earth observation, advances in satellite attitude control capabilities and optical sensor technology have led to the emergence of a new generation of agile satellites capable of performing push-broom imaging during maneuvers. These ultra-agile satellites, such as the Gaofen Multi-Mode Integrated Imaging Satellite, offer in-motion imaging capabilities. Traditional remote sensing satellites maintain a fixed attitude during image acquisition, relying solely on their orbital motion relative to the Earth to achieve swath coverage aligned with the satellite's flight direction. However, ultra-agile satellites can image in a direction other than the satellite's flight direction. Because payload operation is subject to both orbital motion and attitude changes, payload coverage calculation becomes more complex, significantly increasing the challenges of payload coverage analysis.

[0003] Satellite payload coverage capability analysis is the basis for resource planning and task scheduling, and has attracted widespread attention from scholars at home and abroad. Current research methods can be mainly divided into geometric analysis and numerical root-finding methods. Geometric analysis provides the coverage of satellite payloads by establishing geometric relationships based on the relative positions of the satellite and the ground. Mai et al. proposed a two-stage search algorithm for low-orbit sun-synchronous orbits to determine the visibility window of satellite payloads for ground targets. Wu et al. improved this method into a more refined three-phase search method. However, the above method is only applicable to near-circular orbits and has low calculation accuracy. Ulybyshev proposed the concept of a two-dimensional coverage mapping plane, mapping the satellite trajectory into line segments on a plane, and then performing visibility analysis. This method is only applicable to circular orbit spacecraft, circular payload fields of view, and spherical Earth models. Hu Yasi et al. defined the imaging payload by the half-angle of the field of view and determined the coverage range of the field of view on the spherical surface through coordinate transformation. The method has good universality, but the calculation accuracy will decrease in high-latitude areas. Li Dong et al. proposed a great circle approximation algorithm for satellite observation area targets for circular and rectangular payload fields of view. The calculation accuracy of this method is low and it is only applicable to two-body orbits. Under the assumption that the earth is a sphere, Dai et al. proposed a spherical subdivision method to calculate the continuous coverage range of the satellite. Zuo et al. calculated the coverage boundary of the satellite conical payload based on the envelope curve theory, without considering the satellite's attitude maneuver. Subsequently, Han et al. proposed a field of view mapping method, which can map the field of view geometric parameters to a plane composed of the longitude and latitude angles of the ascending node, and then determine the visible window of the target.

[0004] Numerical methods for coverage analysis focus on reducing computational overhead by designing iterative strategies. Alfano et al. used a hybrid parabola method to approximate the visibility curve of conical payloads. Song et al. proposed instantaneous visibility conditions based on a grid point method to calculate satellite coverage performance. Shen Xin et al. equated the satellite coverage area to a characteristic cone to determine whether the imaging conditions were met. This method is only applicable to payloads with a conical field of view. Gu et al. used Kriging interpolation to achieve rapid satellite-to-ground coverage analysis and obtain the satellite's visibility window for ground targets. Han et al. combined adaptive interpolation technology to propose a visibility calculation method that reduces computational time. E Zhibo et al. proposed a basis for the visibility of circular and rectangular fields of view for regional targets and determined the visibility window based on discrete step lengths. Wang et al. proposed a coverage analysis method for circular field of view payloads under satellite attitude maneuvers, which can determine the payload's coverage of regional targets over a period of time. Gu et al. proposed the concept of ground elevation view units and, combining geometric analysis with numerical search methods, solved the problem of calculating satellite payload ground coverage windows under ground elevation constraints.

[0005] However, current approaches still have the following limitations. First, most studies have not examined how to determine payload coverage when satellite orbital motion and attitude maneuvers are simultaneously present. Second, satellite-borne payloads have diverse fields of view, and most research focuses solely on the coverage characteristics of a single field of view type, making it difficult to meet the coverage analysis requirements for multiple payload types in large-scale constellations. Summary of the Invention

[0006] In response to the above-mentioned deficiencies in the prior art, the present invention provides a method for determining the payload coverage strip taking into account satellite attitude maneuvers, which has good effectiveness and efficiency and can be effectively applied to the payload coverage range calculation requirements of ultra-agile earth observation satellites.

[0007] To achieve the above object, the present invention provides a method for determining a payload coverage strip taking into account satellite attitude maneuvers, comprising the following steps:

[0008] Step 1: discretize the field boundary of the imaging payload on the satellite into a number of payload view elements, and calculate the position of the projection point of each payload view element on the earth's spherical surface;

[0009] Step 2: Calculate the velocity of each payload viewpoint projection point based on the velocity change caused by the satellite's circular motion around the earth and the velocity change caused by the satellite's attitude maneuver;

[0010] Step 3: Based on the position and velocity of each load viewer projection point, calculate the index function for each load viewer projection point respectively, and screen out the projection points that can serve as the envelope points of the coverage area based on the index function, thereby obtaining the load coverage strip.

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

[0012] 1. This invention addresses the problem of calculating the instantaneous coverage of a satellite payload's field of view on the Earth's surface and proposes a coverage range determination method based on satellite payload view elements. By calculating the position and velocity of the projection points of the payload view elements on the Earth's spherical surface, the instantaneous coverage area of ​​the payload's field of view on the Earth's surface can be obtained based on the projection point results of different payload view elements, effectively ensuring the accuracy and efficiency of the calculation.

[0013] 2. In order to determine the payload coverage range within a period of time, the present invention proposes a coverage strip determination method based on envelope points, and in the preferred scheme, proposes envelope point criteria for conical loads and rectangular loads respectively. It can connect multiple instantaneous coverage ranges into payload coverage strips, and realize the rapid calculation of the payload coverage strip on the ground under the combined motion of satellite attitude and orbit. BRIEF DESCRIPTION OF THE DRAWINGS

[0014] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on the structures shown in these drawings without paying any creative work.

[0015] Figure 1 Flowchart of a method for determining a payload coverage strip considering satellite attitude maneuvers according to an embodiment of the present invention;

[0016] Figure 2 Schematic diagram of the view element description of the satellite payload field of view in an embodiment of the present invention;

[0017] Figure 3 Schematic diagram of the conversion between the orbital coordinate system and the Earth-centered inertial system in an embodiment of the present invention;

[0018] Figure 4 Schematic diagram of instantaneous coverage of satellite payload field of view in an embodiment of the present invention;

[0019] Figure 5 Schematic diagram of satellite payload projection point velocity calculation in an embodiment of the present invention;

[0020] Figure 6 Schematic diagram of the continuous coverage area of ​​a satellite payload according to an embodiment of the present invention, wherein: (a) is a schematic diagram of a satellite without maneuvering, and (b) is a schematic diagram of a satellite with maneuvering;

[0021] Figure 7 Schematic diagram of the cone load field of view envelope in an embodiment of the present invention;

[0022] Figure 8 Schematic diagram of the envelope points of the rectangular load field of view in an embodiment of the present invention;

[0023] Figure 9 Schematic diagram of the envelope point criterion for rectangular load coverage in an embodiment of the present invention;

[0024] Figure 10 Schematic diagram of the instantaneous coverage of a cone load at different posture angles in an embodiment of the present invention;

[0025] Figure 11 Schematic diagram of the instantaneous coverage of a rectangular load at different posture angles in an embodiment of the present invention;

[0026] Figure 12 This is a comparison diagram of the instantaneous coverage results of rectangular loads in an embodiment of the present invention;

[0027] Figure 13 Schematic diagram of the conical load coverage strip result during the non-attitude maneuver in an embodiment of the present invention;

[0028] Figure 14 Schematic diagram of the conical load coverage strip result during the non-attitude maneuver in an embodiment of the present invention;

[0029] Figure 15 Schematic diagram of the conical load coverage strip results during attitude maneuvering in an embodiment of the present invention;

[0030] Figure 16 This is a schematic diagram of the rectangular load covering strip results during attitude maneuvering in an embodiment of the present invention.

[0031] The purpose, features and advantages of the present invention will be further described with reference to the accompanying drawings and in conjunction with the embodiments. DETAILED DESCRIPTION

[0032] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. All other embodiments obtained by ordinary technicians in this field based on the embodiments of the present invention without making any creative efforts shall fall within the scope of protection of the present invention.

[0033] In addition, the technical solutions between the various embodiments of the present invention can be combined with each other, but it must be based on the fact that ordinary technicians in this field can implement it. When the combination of technical solutions is mutually contradictory or cannot be implemented, it should be deemed that such a combination of technical solutions does not exist and is not within the scope of protection required by the present invention.

[0034] This embodiment discloses a method for determining payload coverage strips taking into account satellite attitude maneuvers, which is mainly used to achieve rapid calculation of ground coverage strips for multiple types of satellite payloads while taking into account satellite orbital motion and attitude maneuvers. Figure 1 In this embodiment, the method for determining the payload coverage strip considering satellite attitude maneuvers specifically includes the following steps:

[0035] Step 1: discretize the field boundary of the imaging payload on the satellite into a number of payload view elements, and calculate the position of the projection point of each payload view element on the earth's spherical surface;

[0036] Step 2: Based on the velocity change caused by the satellite's circular motion around the earth and the velocity change caused by the satellite's attitude maneuver, the velocity of each payload viewpoint projection point is calculated;

[0037] Step 3: Based on the position and velocity of each load viewer projection point, the index function is calculated for each load viewer projection point, and the projection points that can be used as the envelope points of the coverage area are screened out based on the index function to obtain the load coverage strip.

[0038] refer to Figure 2 , the origin S represents the center of mass of the spacecraft, and this coordinate system is the payload installation coordinate system. Vector n s The direction of the ray starting from the center of mass of the spacecraft is called the payload view element. The satellite payload field of view boundary can be determined by a set of payload view elements, and then different satellite payload fields of view can be characterized by discrete payload view elements. Specifically, Figure 2 The load view vector n in the coordinate system shown s To define and characterize, we can define n s is a unit vector, that is, only its vector direction is taken. The angle α in the figure represents n s The angle between the projection on the Sxy plane and the x-axis is recorded as the azimuth; the angle β represents n s The angle with the plane Sxy is recorded as the pitch angle. When N discrete payload view element vectors are used to describe the satellite payload field of view boundary, the angle definition of the payload field of view can be recorded as:

[0039] FOV N ={(α1,β1),(α2,β2),...,(α N ,β N )} (1)

[0040] Medium FOV N is the angle of the payload field of view, α1, ···, α N is the azimuth of the 1st,...,Nth load view vector, β1,...,β N is the pitch angle of the 1st,...,Nth payload view element vector.

[0041] Similar to the projection parameters, the field of view of any simply connected domain can be defined by the field of view angle parameter. Therefore, for the load field of view defined based on the angle method, there is a certain relationship between the pitch angle β and the azimuth angle α, for example:

[0042] β=β(α) (2)

[0043] Formula (2) can be called the angular characteristic function of the payload field of view.

[0044] In order to facilitate unified representation and subsequent derivation, the payload installation coordinate system is assumed to coincide with the satellite body coordinate system. The payload visual element n s The coordinates in the satellite system are expressed as:

[0045]

[0046] Among them, {n s} b is the load view element n s Representation in the satellite system;

[0047] When the satellite undergoes attitude maneuvers, the coordinate transformation matrix between its body coordinate system and the orbital coordinate system (VVLH) is R bo , which is related to the attitude maneuver angle and sequence of the satellite. Assume that the attitude transformation angles of the satellite along the body coordinate axis z, y, and x are ψ, θ, Then the transformation matrix R bo The expression is:

[0048]

[0049] Among them, R x 、R y 、R z Respectively represent the coordinate transformation matrices of the corresponding angles around the x, y, and z axes;

[0050] Assume that the position and velocity vector of the satellite in the Earth-centered inertial coordinate system (ECI) are r s and v s , are expressed as column vectors, such as Figure 3 As shown, where Sx o y o z o is the orbital coordinate system, Sx i y i z i For the ECI coordinate system, we know that:

[0051]

[0052] In formula (6), vector x o 、y o and z o is the direction vector of each axis of the VVLH coordinate system, then the transformation matrix R from the ECI coordinate system to the VVLH coordinate system oi for:

[0053] h s =r s ×v s (7)

[0054]

[0055] Among them, h s is a symbolic representation of the cross product of the position vector and the velocity vector;

[0056] Based on the coordinate transformation matrix R bo and the coordinate transformation matrix R oi , we can get the direction vector of the load view element in the geocentric inertial coordinate system, which is:

[0057] r l ={n s}i =R o ' i R b ' o {n s} b (9)

[0058] Among them, {n s} i 、r l Both are direction vectors of the load view element in the geocentric inertial coordinate system.

[0059] At a specific moment, the position of the satellite in space is fixed. Figure 4 As shown in the figure, the satellite body covers the surface of the earth in a certain attitude. The projection of its payload on the spherical surface of the earth can be called the instantaneous coverage range, that is, the field of view projection area of ​​the satellite payload on the surface of the earth. This area is determined by the payload field of view information, the orbital height and attitude angle of the satellite.

[0060] In order to determine the instantaneous coverage of the satellite payload, based on the payload view element description, we first need to determine the projection point of the payload view element vector on the sphere. Let r s is the satellite's position vector, the satellite's position vector r s and the direction vector r of the load view element l In the Cartesian coordinate system, it can be expressed as:

[0061]

[0062] Among them, x s 、y s 、z s are the position vector r s Projection on the x, y, and z axes, x l 、y l 、z l are direction vectors r l Projection on the x, y, and z axes;

[0063] The parametric equation is used to characterize the ray where the load view element is located, which is:

[0064]

[0065] Among them, x L 、y L 、z L are the projections of the ray where the payload view element is located on the x, y, and z axes respectively. d represents the distance from the satellite payload center to any point on the straight line, and d ≥ 0, which means that the payload view element starting from the payload center is a ray.

[0066] Considering the Earth as a standard ellipsoid model under the WGS84 framework, the calculation accuracy is higher because the effect of the Earth's flattening is taken into account compared to the method of considering the Earth as an ideal sphere. The Earth's equatorial radius is denoted as R E , the Earth's polar radius is denoted as R P , then the analytical expression of the ellipsoid in three-dimensional space is:

[0067]

[0068] Since the Earth-centered Earth-fixed coordinate system (ECEF) is fixed to the ellipsoid, the equations of the ellipsoid remain unchanged in the Earth-centered Earth-fixed coordinate system. It should be noted that the plane Ox i y i It coincides with the Earth’s equatorial plane, and the origin of the ECI coordinate system coincides with the center of the Earth. The Earth only rotates around the z-axis relative to the ECI coordinate system. Therefore, the expression of the ellipsoid equation of the WGS84 frame in the ECI coordinate system is still as shown in Equation (12).

[0069] By expressing both Equations (11) and (12) in the ECI coordinate system, constructing the projection point solution equations together, and solving them, we can obtain the intersection of the ray where the load view element is located and the ellipsoid in three-dimensional space, that is, the projection point of the load view element, and take this projection point as P. From Equations (11) and (12), we can get:

[0070]

[0071] Formula (13) can be transformed into a quadratic equation, which is:

[0072] Ad 2 +Bd+C=0 (14)

[0073]

[0074] Among them, A, B, and C are coefficient parameters;

[0075] By calculating the discriminant Δ=B 2 -4AC, we get two solutions d1 and d2 for the unknown number d, which are:

[0076]

[0077] According to the value of the discriminant, the following situations can be analyzed:

[0078] (1) When Δ>0, there are two roots in the equation solved at the projection point. At this time, the ray where the load view element is located intersects the earth ellipsoid. According to the geometric meaning, only the case where d takes a small positive value is considered at this time;

[0079] (2) When Δ = 0, there is one root in the equation for the projection point. At this time, the ray where the load view element is located is tangent to the earth ellipsoid, and the solution with d ≥ 0 is taken;

[0080] (3) When Δ<0, there is no solution to the equation for the projection point. This embodiment does not consider this situation.

[0081] Therefore, when Δ≥0, a smaller positive value d can be obtained as the final solution, that is, the smaller positive value of d1 and d2 is selected as the distance d between the projection point of the payload view element and the satellite payload center. t , d t ∈{d1,d2}, and then the coordinates of the load view element projection point in the geocentric inertial coordinate system can be obtained:

[0082]

[0083] Among them, x L 、y L 、z L are the coordinate values ​​of the load view element projection point on the x, y, and z axes of the geocentric inertial coordinate system respectively;

[0084] Finally, based on the current moment, the projection point of the payload viewpoint in the Earth-centered inertial coordinate system is converted to the Earth-centered Earth-fixed coordinate system to obtain the latitude and longitude of the projection point of the payload viewpoint on the Earth's spherical surface. After calculating all discrete payload viewpoints of the satellite imaging payload, the set of all projection points can be obtained. Then, by connecting the lines in the order of the payload viewpoints, the instantaneous mapping domain of the payload field of view on the Earth's surface can be obtained.

[0085] After obtaining the instantaneous coverage of the satellite payload field of view, the velocity of each projection point can be further analyzed and calculated to better analyze the coverage of the payload field of view on the spherical surface. Figure 5 As shown in the figure, in inertial space, there are two types of motions of satellite S: 1) The satellite body moves in a circular motion along the orbit, which can be described in the ECI coordinate system, and its speed is v s ; 2) The attitude rotation of the satellite body can be expressed in the VVLH coordinate system, and its angular velocity is The movement of the satellite will cause the position of the payload viewpoint projection point P to change. Figure 5 The symbol d represents the distance from the satellite to point P at this time, r p is the position vector of point P.

[0086] Assume that the attitude conversion angular velocities of the satellite body coordinate axes z, y, and x are and The calculation formula of the satellite's rotation angular velocity in the satellite's body coordinate system is:

[0087]

[0088] in, is the expression of the satellite's rotational angular velocity in the satellite's body coordinate system;

[0089] The coordinate transformation matrix R between the satellite's own system and the VVLH coordinate system bo , we can get:

[0090]

[0091] in, is the satellite rotation angular velocity expressed in the VVLH coordinate system;

[0092] Correspondingly, the velocity of point P in the ECI coordinate system consists of two parts: one is the velocity caused by the circular motion of satellite S around the earth, and the other is the velocity change caused by the satellite attitude maneuver. Let P be the moving point, ECI coordinate system as the fixed coordinate system, and VVLH coordinate system as the moving coordinate system, then the absolute velocity v of point P in the ECI coordinate system is p It can be expressed as:

[0093]

[0094] in, is the velocity of the instantaneous coincidence point of the load view projection point P in the VVLH coordinate system, is the velocity of the load view element projection point P relative to the VVLH coordinate system;

[0095] The motion type of the instantaneous coincidence point is consistent with that of the satellite S, both of which are circular motions around the earth, and the calculation formula of the satellite speed is:

[0096]

[0097] in, is the equivalent earth angular velocity;

[0098] The velocity of the instantaneous coincidence point is:

[0099]

[0100] The velocity of the load view element projection point P relative to the VVLH coordinate system is:

[0101]

[0102] According to the aforementioned transformation matrix R bo and the transformation matrix R oi , we can get the velocity in ECI coordinate system and relative speed The absolute velocity of the projection point {vP} i .

[0103] Since the load viewpoint projection point P is a projection point on the earth's spherical surface, its motion will also be restricted to the spherical surface, so its velocity direction is only parallel to the tangent plane of the point, and the velocity perpendicular to the tangent plane will be restricted to 0. When assuming that the earth is a sphere, r p is the normal direction, then the velocity perpendicular to the tangent plane is for:

[0104]

[0105] Based on the absolute speed v p Velocity perpendicular to the tangent plane The true velocity of the load viewpoint projection point can be calculated for:

[0106]

[0107] Performing the above calculations for each payload view element of the satellite yields the velocity set for all view element projection points. Velocity calculations are crucial for analyzing the movement and changes in an instantaneous coverage area, serving as the link between two instantaneous coverage areas.

[0108] Satellite coverage is the area on the ground that is covered by the field of view of the satellite's onboard payload during normal operation in orbit. Traditional Earth observation satellites operate in push-broom mode. The satellite's attitude does not change during imaging, and it relies on the satellite's orbital motion to continuously cover the ground area. At this time, the attitude of the Earth-pointing satellite is fixed, and its coverage range is as follows: Figure 6 As shown in (a). With the advancement of satellite technology, more and more agile satellites have the ability to perform imaging in motion, such as Figure 6 As shown in (b), the satellite can control the area swept by its payload through various attitude maneuvers, thereby observing different areas on the Earth's surface.

[0109] Normally, points inside the area cannot affect the coverage range, so the coverage area can be completely determined by the boundary points. For two consecutive coverage areas, there is a one-to-one correspondence between their boundary discrete points, that is, the point at time t+Δt is obtained by the movement of the point at time t. However, between two consecutive coverage areas, it is necessary to determine their outer envelope, such as Figure 6 (a) Figure 6The outer envelope connects the two coverage areas with a smooth curve. By segmenting the individual instantaneous coverage areas with the envelope line, a continuous load coverage strip is obtained. The key to determining the envelope line is to identify the envelope points within a single instantaneous field of view. Once the envelope points are determined, the coverage strip can be quickly obtained by connecting the discrete envelope points.

[0110] In this embodiment, when the field of view of the imaging payload on the satellite is a conical payload field of view, the specific implementation process of step 3 is:

[0111] The cone load field of view has a continuous and smooth characteristic, and the angular characteristic function of the load field of view is continuously differentiable, that is, there exists a differential function D in equation (2) that satisfies:

[0112]

[0113] Furthermore, the direction vector D of the coverage area at each load viewpoint projection point can be determined based on the specific shape of the conical load field of view. t ,like Figure 7 As shown, it can be seen that for the conical field of view, the direction vector D at point P1 t exists and is related to the numbering order of the load view element. After obtaining the coordinates of the projection point of each load view element, the direction vector D of the smooth field of view at each point can be calculated by the numbering order of the projection point. t ;

[0114] from Figure 7 It can be seen that when the load coverage area moves from the position defined by P1, A1, and B1 to the position defined by P2, A2, and B2, the speed of point P1 moving to point P2 can be expressed as vector It is not difficult to find that point P1 is not an envelope point, because its velocity vector is included in the coverage area and is not tangent to the boundary of the adjacent mapping area. Figure 6 Points A1 and B1 in the figure are envelope points. During the movement of the coverage area, the movement of the two points is tangent to the coverage area, that is, the coverage area is enclosed within the envelope. When the adjacent time interval Δt is small enough, there exists:

[0115]

[0116] Among them, v A is the velocity vector from point A1 to point A2, |A1A2| is the distance between points A1 and A2, v B is the velocity vector from point B1 to point B2, |B1B2| is the distance between point B1 and point B2;

[0117] Therefore, the index function e of the load coverage area envelope can be defined as:

[0118] e=D t ×v t (30)

[0119] Among them, v t is the velocity vector of the load view element projection point;

[0120] The direction of motion of the envelope point is collinear with the direction vector of the boundary of the load field of view. Based on this feature, the envelope point determination condition can be obtained as follows:

[0121] e=0 (31)

[0122] Based on the position and velocity of the cone load viewpoint on the Earth's surface at each moment, the index function e can be calculated for each projection point, thereby determining the envelope point of the coverage area. By connecting the corresponding envelope points on the instantaneous coverage area at different times on the same side, the load band over a period of time can be obtained.

[0123] In the specific application process, the load field of view needs to be a continuous and smooth field of view. At this time, the direction vector D of the boundary of the load view element coverage area is t This is true for a conical load field of view. However, there will be problems for a rectangular load field of view commonly used in applications. The elevation angle of the rectangular field of view at the vertex is not differentiable with respect to the azimuth angle, resulting in its projection area being continuous but not smooth at the vertex view element, and thus it is impossible to calculate the direction vector D along the tangent line. t Therefore, the envelope point determination condition given in the form of formula (31) cannot be used for the load coverage area under the rectangular field of view.

[0124] The coverage of the rectangular load field of view over a period of time is as follows: Figure 8 As shown in the figure, the outer bounding box formed by the thick solid line is the load coverage strip, and the dots are the envelope points of the coverage area during this period. It is not difficult to see that the coverage strips are all inside the envelope points. From time t to time t+Δt, the four vertices A1, B1, C1, and D1 of the coverage range move to positions A2, B2, C2, and D2, respectively. For example, at point B1, although the coverage boundary is continuous, there is a discontinuity between A1, B1 and B1, C1 at point B1, which means that their direction vector cannot be calculated.

[0125] For non-smooth continuous fields of view similar to rectangular load fields of view, new envelope point determination conditions need to be proposed. Assume that there are very small guide angles at the four vertices of the rectangular field of view, thereby achieving the differentiability of the elevation angle to the azimuth angle at the guide angle. Since the differential function at the guide angle changes dramatically, it can be approximately regarded as the existence of two sets of tangent vectors. Figure 9 As shown, along the order of the four vertices A1, B1, C1, and D1, there are direction vectors at point B1. and There are direction vectors at point C1 and At time t, the velocities of B1 and C1 are and Depend on Figure 8 It can be seen that point C1 is the envelope point.

[0126] According to the judgment condition of formula (31), it is determined whether the vertex of the field of view is an envelope point, which is equivalently converted into two sets of tangent vectors D at the guide angle. t1 、D t2 With the vertex velocity vector v t The relationship is determined when:

[0127] (D t1 ×v t )·(D t2 ×v t )≤0 (32)

[0128] According to the continuity theorem, there must be a t1 With D t2 The vector between

[0129]

[0130] Where λ is the proportionality coefficient;

[0131] satisfy:

[0132]

[0133] Therefore, the small guide angle assumption is adopted to make the vertex of the field of view continuous and differentiable, and the envelope point determination method under the continuous smooth field of view is modified. The envelope determination basis of the rectangular load field of view is given by formula (32), and the method for determining the continuous coverage strip of the load is improved. Figure 8 Calculate the vertices B1 and C1 in the image separately, and you can get:

[0134]

[0135] According to the calculation results of formula (35), C1 is the vertex envelope point and B1 is the non-envelope point, which is consistent with Figure 8 The results shown in are consistent.

[0136] Therefore, in this embodiment, when the field of view of the imaging payload on the satellite is a conical payload field of view, the specific implementation process of step 3 is:

[0137] Based on the specific shape of the conical load field of view, determine the direction vector D of the load view element projection point at the corner of the coverage area t1 、D t2 ,

[0138] Calculate the index function e corresponding to the projection point of each load view element at the corner point of the coverage area, which is e=(D t1 ×v t )·(D t2 ×v t );

[0139] The load viewer projection points with the index function e≤0 are screened out, which are the envelope points of the coverage area. The envelope points at different times on the same side are connected to obtain the load coverage strip within a period of time.

[0140] The method for determining the payload coverage strip considering satellite attitude maneuvers in this embodiment is further described below with reference to specific examples.

[0141] The start time of the simulation scenario is set to 0:00 on July 1, 2024, and the simulation duration is 12 hours. The orbital parameters of the satellite at the initial moment are shown in Table 1. The orbit prediction model is set to J2, considering the satellite's attitude maneuverability from -30° to +30°.

[0142] Table 1 Satellite orbit parameters

[0143]

[0144] The satellite payload's field of view was set to either conical or rectangular. The conical payload's field of view half-angle was set to 10°, while the rectangular payload's horizontal and vertical half-angles were set to 8°. The experiment involved analyzing the satellite payload's regional coverage. A rectangular region R was added, located in Africa, with vertex latitudes and longitudes at (5.33, 30.26), (5.33, 9.32), (-5.17, 9.32), and (-5.17, 30.26), respectively. In the absence of attitude maneuvers, the satellite payload's visibility window for this region was calculated. The conical payload's visibility of the target was [37696.867s, 37900.902s], while the rectangular payload's visibility was [37698.518s, 37899.342s]. The visibility windows for the two payloads were similar, but there were differences. Subsequent experiments were conducted based on this scenario.

[0145] To verify the accuracy of the instantaneous coverage of the satellite payload, this example verifies the conical and rectangular fields of view respectively, calculates the payload coverage mapping domain results at different attitude angles, and compares them with the instantaneous mapping domain of the payload in STK. Taking the instantaneous time as 37696.867s, the payload coverage of the satellite at different attitude maneuvering angles is as follows: Figure 10 and Figure 11 As shown, Figure 10 and Figure 11The coverage of conical and rectangular payloads in the attitude angle range from -30° to +30° is shown respectively.

[0146] like Figure 10 As shown in the figure, different line types represent different attitude angles. It can be found that the payload coverage range varies in different attitudes. For example, the coverage area is significantly larger at -30° and +30°. This is because the area distortion caused by the satellite's side swing. The smaller the attitude maneuvering angle, the smaller the distortion of the coverage range. At the same time, it can be seen that when the attitude maneuvering angle interval is 20°, the coverage area is tangent, as shown in the figure. Figure 10 Point A is the tangent point of the instantaneous coverage range under -30° and -10° attitude maneuvers. This is because the semi-cone angle of the load is set to 10°. Figure 10 Point A is located at the rightmost end of the thick solid line area on the left and at the leftmost end of the third thin solid line area from the left. At this time, the satellite is at a side swing angle of -20°.

[0147] The instantaneous coverage results of the rectangular load at different posture angles are as follows: Figure 11 As shown in the figure, it is not difficult to find that when there is attitude maneuvering, the coverage range at the vertex of the rectangular payload field of view will be greatly distorted, and the distortion will be aggravated as the side swing angle increases. Figure 11 The coverage distortion is most serious when the angle is between -30° and +30°, and the coverage area is much larger than when the side swing angle is 0°. Figure 10 and Figure 11 The time consumption of the 7 groups of instantaneous coverage areas is less than 0.01s, so the calculation of the instantaneous coverage area is very efficient, which verifies the efficiency of the method of this embodiment.

[0148] In order to verify the accuracy of the instantaneous coverage area calculated by the method of this embodiment, the rectangular payload coverage result when the satellite side swing angle is +30° is compared with the effect of STK, as shown in the following figure: Figure 12 The area enclosed by the dotted line is the coverage calculated by the method of this embodiment (corresponding to Figure 11 The area inside the dotted line is the actual coverage of the satellite payload in STK. Figure 12 The lower left corner shows the coverage of the three-dimensional scene in STK. It can be found that the results calculated by the method of this embodiment are completely consistent with the load coverage range displayed by STK, which verifies the accuracy of the method of this embodiment.

[0149] In addition, this example also verifies the accuracy of payload coverage strip calculation through simulation tests. The coverage strips of the conical and rectangular payloads carried by the satellite within the visible range of the target area R are calculated respectively. The calculation interval of adjacent instantaneous coverage areas is set to 5s. The test results are shown as follows: Figure 13 and Figure 14 shown.

[0150] The results of the cone load covering strip are as follows Figure 13 As shown in the figure, the area enclosed by the thin solid line is the target area R, and the area enclosed by the thick solid line is the strip calculated by the method of this embodiment. By importing the strip results into STK, the following comparison can be performed. The three groups of sub-figures on the right show the instantaneous load coverage at different times, including the start, end, and midpoint of the interval [37696.867s, 37900.902s]. Based on the load coverage results, it is not difficult to see that the instantaneous coverage area is tangent to the strip, verifying the accuracy of the coverage strip calculation.

[0151] The rectangular load covering strip results are as follows Figure 14 As shown in the figure, it can be found that at the end points of the interval [37698.518s, 37899.342s], the two sides of the rectangular load field of view completely coincide with the coverage strip area, and the two vertices of the instantaneous coverage area at the midpoint are located on the coverage strip, which can verify the accuracy of the rectangular coverage strip calculation. Figure 13 and Figure 14 The time taken for the load to cover the strip is 0.053s and 0.025s respectively, reflecting the high efficiency of the calculation method of this embodiment.

[0152] In order to further verify the calculation results of the method of this embodiment numerically, Table 2 counts the projection latitude and longitude coordinates of the four vertices of the rectangular load field of view on the spherical surface at different times and compares them with the STK results. By numerically solving the equation, the accurate projection results of each view element can be directly obtained. Table 2 retains 4 decimal places; in STK, the number of decimal places to be retained needs to be selected, and three decimal places are selected in the simulation. It can be found from Table 2 that the calculation results of the projection coordinates of the vertices of the rectangular load field of view are completely consistent with the STK results, which is consistent with the Figure 14 The results shown are consistent, verifying the accuracy of the coverage range calculation method in this embodiment.

[0153] Table 2 Results of longitude and latitude mapping of the vertex of the rectangular load field of view

[0154]

[0155] To verify the payload coverage strip results when the satellite is in combined attitude and orbit motion, the satellite motion time interval is set to [0s, 300s], and the satellite attitude angle parameters at different times are set as shown in Table 3. The satellite attitude angle is linearly varied between adjacent times, and the interval for calculating the satellite instantaneous coverage is set to 5s.

[0156] Table 3 Satellite attitude maneuvering angles

[0157]

[0158] The payload coverage strip results when the satellite attitude angle changes are as follows: Figure 15 and Figure 16 As shown in the figure, it can be found that the coverage strip is not parallel to the satellite sub-satellite point trajectory, but forms a zigzag change as the satellite attitude maneuvers.

[0159] Figure 15 The area enclosed in the figure is the coverage strip calculated by the method of this embodiment. The corresponding data is imported into the STK software to obtain the area in the figure. The dotted line range is the coverage strip of the satellite payload in the interval [0s, 300s] displayed in the STK. Figure 15 The circular colored area in the figure represents the instantaneous coverage of the conical payload. The three small figures on the right show the relationship between the instantaneous coverage area of ​​the payload and the calculated strip at 50s, 200s, and 250s, respectively. It is not difficult to see that the instantaneous coverage area is tangent to the strip, verifying the accuracy of the method for determining the coverage strip of the satellite attitude-orbit composite motion payload. Furthermore, the results at 150s show that the edge of the instantaneous coverage range of the payload is tangent to the red line, while the dotted line lies within the instantaneous coverage area. This indicates that there are errors in the strip calculation within STK. However, the method in this embodiment has good adaptability in determining the strip at the inflection point of the attitude maneuver.

[0160] Figure 16 The coverage strip results of the rectangular payload under the satellite's attitude maneuver are shown. The rectangular color block area is the instantaneous coverage range of the rectangular payload. Figure 16 The degree of agreement between the solid and dashed lines verifies the accuracy of the coverage strip calculation. It can also be seen that at 50s, 100s, and 200s, both vertices of the rectangular field of view are located on the strip, and at 50s, one edge of the payload field of view is tangent to the coverage strip, verifying the accuracy of the rectangular payload field of view coverage strip calculation method. It can also be seen that the envelope points of the coverage area are all composed of the mapping points of the vertex view elements. This is because during the satellite's motion, as the rectangular payload field of view's instantaneous ground coverage area changes over time, it always first contacts or finally leaves the instantaneous coverage area at the vertex of the instantaneous coverage domain.

[0161] Calculated Figure 15 and Figure 16 The time required to cover the strip within 300 seconds is related to the time interval for calculating the instantaneous coverage area. The calculation time statistics based on the simulation results are shown in Table 4.

[0162] Table 4 Calculation time of coverage strip results

[0163]

[0164] It can be seen that as the instantaneous coverage area calculation interval shortens, the calculation time of the coverage strip increases linearly. After comparison with STK, the results are very accurate when the instantaneous coverage area calculation interval is 5 seconds. When the interval is 1 second, the calculation time does not exceed 0.7 seconds. This verifies the accuracy and efficiency of the method of this embodiment. In actual applications, different calculation intervals can be selected according to the needs of calculation accuracy and efficiency.

[0165] It can also be found that the calculation time for the rectangular payload coverage strip is shorter. This is because in the envelope point determination phase, the rectangular payload field of view only needs to determine the vertex, while the conical payload requires calculation and judgment of all discrete view elements, which will increase the calculation time. It should be pointed out that in the simulation experiment, only the satellite's roll angle change is set, which facilitates the comparison and display of the results. However, as shown in formula (20), in the specific application process, the method of this embodiment can also be applied to the calculation of coverage strips under arbitrary satellite attitude maneuvers.

[0166] The above description is only a preferred embodiment of the present invention and does not limit the patent scope of the present invention. All equivalent structural transformations made by using the contents of the present invention description and drawings under the inventive concept of the present invention, or direct / indirect application in other related technical fields are included in the patent protection scope of the present invention.

Claims

1. A method for determining payload coverage strips considering satellite attitude maneuvers, characterized in that: The steps include: Step 1: discretize the field boundary of the conical or rectangular imaging payload on the satellite into a number of payload view elements, and calculate the position of the projection point of each payload view element on the earth's spherical surface; Step 2: Calculate the velocity of each payload viewpoint projection point based on the velocity change caused by the satellite's circular motion around the earth and the velocity change caused by the satellite's attitude maneuver; Step 3: Based on the position and velocity of each load viewme projection point, calculate the index function for each load viewme projection point, and screen out the projection points that can serve as the envelope points of the coverage area based on the index function, thereby obtaining the load coverage strip; When the field of view of the imaging payload on the satellite is a conical payload field of view, step 3 specifically includes: Based on the specific shape of the conical load field of view, determine the direction vector of the coverage area at each load viewpoint projection point D t ; Calculate the index function corresponding to each load viewpoint projection point e ,for ,in, is the velocity vector of the load view element projection point; Filter out indicator functions The load view element projection point is the envelope point of the coverage area. Connecting the envelope points at different times on the same side will give the load coverage strip within a period of time. When the field of view of the imaging payload on the satellite is a rectangular payload field of view, step 3 specifically includes: Based on the specific shape of the conical load field of view, determine the direction vector of the load view element projection point at the corner of the coverage area D t1 、 D t2 ; Calculate the index function corresponding to the projection point of each load view element at the corner point of the coverage area e ,for ,in, is the velocity vector of the load view element projection point; Filter out indicator functions The load view element projection point is the envelope point of the coverage area. By connecting the envelope points at different times on the same side, the load coverage strip within a period of time is obtained.

2. The method for determining the payload coverage strip considering satellite attitude maneuvers according to claim 1, characterized in that: The position of the projection point of the load view element on the earth's spherical surface specifically includes: Get the satellite's position vector in the geocentric inertial coordinate system and the direction vector of the load view element ; The position vector and direction vector They are represented in the Cartesian coordinate system as follows: in, 、 、 are position vectors Projection on the x, y, and z axes, 、 、 are direction vectors Projection on the x, y, and z axes; The parametric equation is used to characterize the ray where the load view element is located, which is: in, x L 、 y L 、 z L are the projections of the ray where the load view element is located on the x, y, and z axes, respectively. d ≥ 0 indicates the distance from the satellite payload center to any point on the straight line; Based on the standard ellipsoid model of the earth and the ray where the load view element is located, the projection point solution equation is constructed as follows: in, R E is the equatorial radius of the Earth, R P is the Earth's polar radius; Solve the projection point equation to get the distance between the projection point of the payload view element and the satellite payload center. d t , and the coordinates of the load view element projection point in the geocentric inertial coordinate system are obtained as follows: in, x L 、 y L 、 z L are the coordinate values ​​of the load view element projection point on the x, y, and z axes of the geocentric inertial coordinate system respectively; According to the current time, the projection point of the load view element in the geocentric inertial coordinate system is converted to the geocentric earth-fixed coordinate system to obtain the latitude and longitude of the projection point of the load view element on the earth's spherical surface.

3. The method for determining the payload coverage strip considering satellite attitude maneuvers according to claim 2, characterized in that: The direction vector of the load view element is obtained in the geocentric inertial coordinate system r l , specifically: Assume that the payload installation coordinate system coincides with the satellite body coordinate system, and obtain the payload view element n s The expression in the satellite system is: in,{ n s } b Load view element n s In the satellite system, α Load view element n s The azimuth of β Load view element n s Pitch angle; Get the coordinate transformation matrix between the satellite body coordinate system and the orbit coordinate system R bo ,for: in, ψ 、 θ 、 φ are the attitude conversion angles of the satellite along the z, y, and x axes of the satellite coordinate system, respectively. R x 、 R y 、 R z Respectively represent the coordinate transformation matrices of the corresponding angles around the x, y, and z axes; Get the coordinate transformation matrix between the orbital coordinate system and the geocentric inertial coordinate system R oi ,for: in, v s is the satellite's velocity vector, h s is a symbolic representation of the cross product of the position vector and the velocity vector; Based on the coordinate transformation matrix R bo and the coordinate transformation matrix R oi , we can get the direction vector of the load view element in the geocentric inertial coordinate system: in,{ n s } i 、 r l Both are direction vectors of the load view element in the geocentric inertial coordinate system.

4. The method for determining the payload coverage strip considering satellite attitude maneuvers according to claim 2, characterized in that: The distance between the projection point of the payload view unit and the satellite payload center d t The specific calculation process is: The projection point solution equation is organized into a quadratic equation form, which is: in, A 、 B 、 C is the coefficient parameter; By calculating the discriminant , get the unknown number d Two solutions of d 1. d 2, for: choose d 1. d The smaller positive value of the two is used as the distance between the projection point of the payload view element and the satellite payload center. d t .

5. The method for determining a payload coverage strip considering satellite attitude maneuvers according to any one of claims 1 to 4, characterized in that: Step 2 is as follows: Taking the geocentric inertial coordinate system as the fixed coordinate system and the orbital coordinate system as the moving coordinate system, the absolute velocity of the load view element projection point in the geocentric inertial coordinate system is obtained. v p ,for: in, is the velocity of the instantaneous coincidence point of the projection point of the load visual element in the orbital coordinate system, is the velocity of the load view element projection point relative to the orbit coordinate system; Calculate the velocity of the load view element projection point on the vertical tangent plane in the geocentric inertial coordinate system ,for: in, r p is the position vector of the load view element projection point; Based on absolute speed v p Velocity perpendicular to the tangent plane , the true velocity of the load view element projection point is calculated as: in, is the true velocity of the load view element projection point.

Images