Method for calculating visible window of optical satellite to ground target

Abstract

The invention relates to a method for calculating a visible window of an optical satellite to a ground target, which comprises the following steps of: establishing a visible judgment function of the optical satellite to the ground target by considering the orbit attitude motion of the optical satellite and the installation and view field parameters of an on-satellite optical camera, fitting the judgment function by adopting a piecewise cubic polynomial function, and calculating the visible window of the optical satellite to the ground target; calculating a visible time window of the optical satellite to the ground point target; calculating a visible time window of the optical satellite to the ground line target by adopting a discretization method; considering the motion law of the optical satellite, and calculating a visible time window of the optical satellite to the ground triangular area target; a polygon triangulation method is adopted to calculate a visible time window of the optical satellite to any polygon area target on the ground.

Description

Technical Field

[0001] The present invention relates to the field of optical satellite technology, and in particular to a method for calculating a visible time window of an optical satellite to a ground target. Background Art

[0002] With the development of aerospace technology, optical satellite remote sensing observation technology for the Earth has been widely used in many fields, including meteorological observation, environmental monitoring, disaster warning, and urban planning. Rapidly and accurately calculating the time window during which optical satellites are visible to ground targets is fundamental to arranging optical remote sensing satellites for ground observation missions. Ground targets can be point targets such as airports and ports, linear targets such as rivers and borders, or arbitrary polygonal ground areas. Furthermore, while the optical camera is powered on, optical satellites may be scheduled to perform attitude maneuvers to enable complex optical remote sensing observations of the Earth. Meeting these practical application requirements requires rapidly and accurately calculating the time window during which optical satellites are visible to ground targets, a relatively complex problem.

[0003] Currently, the traditional method for calculating the visibility time window of an optical satellite to a ground target is the tracking propagation method, also known as the "brute force method." This method uses a set step size and calculates the visibility of the optical satellite to the ground target at each moment, thereby obtaining the visibility time window of the optical satellite to the ground target. The tracking propagation method is intuitive and simple to operate, but its calculation accuracy and efficiency are highly dependent on the selected step size, resulting in low computational efficiency. Summary of the Invention

[0004] In order to solve the technical problems existing in the above-mentioned prior art, the purpose of the present invention is to provide a method for calculating the visible window of an optical satellite to a ground target, which can realize the rapid and accurate calculation of the optical satellite to the ground point, line and area targets.

[0005] To achieve the above-mentioned object of the invention, the present invention provides a method for calculating the visibility window of an optical satellite to a ground target, comprising the following steps:

[0006] Step S1: Considering the orbital attitude motion of the optical satellite and the installation and field of view parameters of the optical camera on the optical satellite, a visibility judgment function of the optical satellite to the ground point target is established, and a piecewise cubic polynomial function is used to fit the judgment function to calculate the time window in which the optical satellite is visible to the ground point target;

[0007] Step S2: Based on the time window in which the optical satellite is visible to the ground point target, a discretization method is used to calculate the time window in which the optical satellite is visible to the ground line target;

[0008] Step S3, based on the time window in which the optical satellite is visible to the ground line target and taking into account the motion law of the optical satellite, calculating the time window in which the optical satellite is visible to the ground triangular area target;

[0009] Step S4: Based on the visible time window of the optical satellite to the ground triangular area target, a polygon triangulation method is used to calculate the visible time window of the optical satellite to the ground arbitrary polygon area target.

[0010] According to a technical solution of the present invention, in step S1, it specifically includes:

[0011] Step S11: Calculating the Earth-fixed position vector and the Earth-fixed velocity vector of the optical satellite at any time within the time interval of the optical satellite Earth-fixed ephemeris based on the optical satellite Earth-fixed ephemeris;

[0012] Step S12: Calculating the attitude data of the optical satellite at any time within the time interval of the optical satellite Earth-fixed system ephemeris based on the optical satellite attitude sequence;

[0013] Step S13, calculating the conversion matrix from the ground-fixed system to the camera installation coordinate system of the optical satellite at any time according to the optical satellite attitude data and the installation parameters of the optical camera of the optical satellite;

[0014] Step S14: Considering the optical satellite attitude motion and the optical satellite optical camera installation and field of view parameters, a determination function is established as to whether the optical satellite is visible to the ground point target;

[0015] Step S15: establishing a first-order derivative function and a second-order derivative function of the determination function based on the central difference quotient;

[0016] Step S16: Adopting an adaptive step size strategy, constructing a piecewise cubic Hermite interpolation polynomial function to fit the determination function;

[0017] Step S17, finding roots of the piecewise cubic Hermite interpolation polynomial function, and calculating a time window set that satisfies the determination function;

[0018] Step S18: performing a set operation on the set that satisfies the judgment function to obtain a time window in which the optical satellite is visible to the ground point target.

[0019] According to a technical solution of the present invention, in step S11, the optical satellite earth-fixed system ephemeris is represented as t i represents the ephemeris time node in the optical satellite ground-fixed ephemeris, R i and represents the ephemeris time node t in the optical satellite fixed system ephemeris iThe corresponding Earth-fixed position vector and Earth-fixed velocity vector of the optical satellite;

[0020] For any time t∈[T s , T e ], the ground-fixed position vector and ground-fixed velocity vector of the optical satellite are expressed as R(t)=(x(t), y(t), z(t)) and The ephemeris time node in the optical satellite ground fixed system corresponding to time t is obtained by segmented 5-times Lagrangian interpolation. i+k |k=0,1…5,t i ≤t≤t i+5}, k represents the number of Lagrange interpolation, then

[0021]

[0022] Among them, l k (t) is the kth Lagrange interpolation basis function, expressed as

[0023]

[0024] In step S12, the optical satellite attitude sequence is represented as Att={(t i , Q i )|i=0,1…N,t i ∈[T s , T e ]}, where Q i =(qa i ,qb i ,qc i ,qd i ) is t i The optical satellite body coordinate system at the moment is relative to the unit quaternion of the VVLH coordinate system. The first three are the vector part of the quaternion, and the fourth is the scalar part of the quaternion.

[0025] For any time t∈[T s , T e ], the attitude data of the optical satellite is expressed as Q(t) = (qa(t), qb(t), qc(t), qd(t)), the attitude data of the optical satellite is obtained by piecewise spherical linear interpolation according to the optical satellite attitude sequence, and the time node in the optical satellite attitude sequence corresponding to time t is t i , t i+1 (t i ≤t≤t i+1 ), then

[0026]

[0027] Where k0, k1 and θ are the starting weight, ending weight and the angle between the two quaternions, respectively, expressed as

[0028]

[0029] θ=acos(qa i ×qa i+1 +qb i ×qb i+1 +qc i ×qc i+1 +qd i ×qd i+1 );

[0030] In step S13, the conversion matrix from the ground-fixed system of the optical satellite to the camera installation coordinate system is expressed as

[0031] M(t)=M sb M bo (t)M of (t)

[0032] Among them, M sb is the camera installation matrix, which is a constant matrix; M bo (t) is the transformation matrix from the VVLH coordinate system to the optical satellite body coordinate system, M of (t) is the transformation matrix from the ground-fixed system to the VVLH coordinate system of the optical satellite, M b。 (t) is calculated from the attitude data Q(t) of the optical satellite at time t and is expressed as

[0033]

[0034]

[0035] M of (t) is the ground-fixed position vector R(t) and ground-fixed velocity vector of the optical satellite at time t Calculated, expressed as

[0036]

[0037] Among them, I x (t), I y (t), I z (t) and V(t) represent the coordinates of the unit vector of the VVLH coordinate axis in the ground-fixed system, w=(0, 0, 7.29211581560071e-5) is the Earth's rotation angular velocity vector.

[0038] According to a technical solution of the present invention, in step S14, the following steps are specifically included:

[0039] Step S141: Obtain the geodetic coordinates P(L, B, H) of the ground point target, and convert the geodetic coordinates of the ground point target into rectangular coordinates to obtain the ground fixed vector R of the ground point target. p =(x, y, z), expressed as

[0040]

[0041] Wherein, L is longitude, B is latitude, H is altitude; x, y and z represent the rectangular coordinates of the ground point target fixed vector respectively; N is the radius of curvature of the circle, a and e are the semi-major axis and eccentricity of the Earth's ellipsoid, respectively;

[0042] Step S142: transform the ground fixed system vector of the ground target point into the camera installation coordinate system, and construct a ground point target determination function within the optical satellite field of view, expressed as

[0043] r(t)=M(t)(R(t)-R p )

[0044]

[0045] Wherein, α0 and β0 represent the vertical half angle and horizontal half angle of the field of view of the optical satellite respectively; r y (t), r x (t), r z (t) respectively represent the coordinate components of the vector from the optical satellite to the ground point target in the camera installation coordinate system; F1(t) and F2(t) are the ground point target determination functions within the field of view of the optical satellite;

[0046] Step S143: construct a judgment function for the minimum elevation angle of the optical satellite to the ground point target, expressed as

[0047]

[0048] Among them, n p Represents the geoid normal vector of the ground point target, n p =(cosBcosL, cosBsinL, sinB); θ0 is the minimum elevation angle of the optical satellite to the ground point target;

[0049] Step S144: construct a ground solar altitude angle determination function, expressed as

[0050]

[0051] Where θ1 is the minimum solar altitude angle on the ground; R sun(t) is the sun position vector of the Earth-fixed system at time t;

[0052] Step S145: simultaneously satisfy the judgment function F i1 When (t+Δt)≤0, i1=1, 2, 3, 4, the optical satellite is visible to the ground point target.

[0053] According to a technical solution of the present invention, in step S15, the first-order derivative function and the second-order derivative function of the optical satellite to ground point target visibility determination function are established based on the central difference quotient, which are expressed as follows:

[0054]

[0055] Wherein, Δt is the difference quotient step size, which is fixed at 0.1 seconds.

[0056] According to a technical solution of the present invention, in step S16, constructing a piecewise cubic Hermite interpolation polynomial function specifically includes:

[0057] Step S1601: Set the maximum interpolation step length mdt to one quarter of the orbital period;

[0058] Step S1602: Set the current interpolation node t c =T s , calculate the judgment function at t c The function value and derivative value F(t c ), F′(t c )、F″(t c ), save the current interpolation node data t c 、F(t c ), F′(t c )、F″(t c );

[0059] Step S1603: Calculate the interpolation step length δt=min(mdt, T e -t c );

[0060] Step S1604: Set the subsequent interpolation node t n =t c +δt, calculate the judgment function at t n The function value and derivative value F(t n ), F′(t n )、F″(t n );

[0061] Step S1605: In the interpolation interval [t c , t n ], construct the two-point cubic Hermite interpolation polynomial I h(t), calculate I h (t) The second-order derivative value I″ at the end point of the interpolation interval h (t c ), I″ h (t n );

[0062] Step S1606: Calculate the fitting error D = |I″ h (t c )-F″(t c )|+|I″ h (t n )-F″(t n )|, if D≥D error Then shorten the interpolation step δt=0.5δt and go to step S1604;

[0063] Step S1607: Calculate I′ h (t)=0 in the open interval (t k , t k+1 ) in the root, if the equation has roots t k1 , t k2 And t k1 ≤t k2 , then adjust the interpolation step δt=t k1 -t c , go to step S1604;

[0064] Step S1608: Save subsequent interpolation node data t n 、F(t n ), F′(t n )、F″(t n );

[0065] Step S1609: Update the current interpolation node data t c =t n 、F(t c )=F(t n ), F′(t c )=F′(t n )、F″(t c )=F″(t n );

[0066] Step S1610: If t c <T e , then go to step S1603, otherwise end the calculation process and obtain the difference node set T s =t0<t1…t n-1 <t n =T e , where n represents the number of interpolation nodes.

[0067] According to a technical solution of the present invention, in step S17, the following steps are specifically included:

[0068] Step S171, set the time window flag isInWin=false;

[0069] Step S172: The difference node set T s =t0<t1…t n-1 <t n =T e The first interpolation node data in is used as the current interpolation node data t c 、F(t c ), F(t c )<0, set the window start time t s =t c , set isInWin = true;

[0070] Step S173: Obtain subsequent interpolation node data t n 、F(t n );

[0071] Step S174: If F(t c )×F(t n )≤0, calculate I h (t) in [t c , t n ] root t r Otherwise, go to Step 5; if isInWin is true, set the window end time t e =t r And save the window data; if isInWin is false, set the window start time t s =t r ; Modify the time window flag isInWin = !isInWin;

[0072] Step S175: Update the current interpolation node data t c =t n 、F(t c )=F(t n );

[0073] Step S176: If t c <T e , then go to step S173, otherwise go to step S177;

[0074] Step S177: If isInWin is true, set the window end time t e =t c And save the window data.

[0075] According to a technical solution of the present invention, in step S2, the optical satellite to ground line target visibility window algorithm includes:

[0076] The two endpoints of the ground line target are respectively P a (L a , B a ), P b (L b , B b ), the equidistant division method is used along the geodesic line to discretize the ground line target into a set of ground point targets {P i2 (L i2 ,B i2 ,0)|i2=1,2,…,N,N≥2}, apply the optical satellite to ground point target visibility window algorithm to obtain the visible time window of the optical satellite to each ground point target in the ground point target set, and perform union calculation to obtain the visible time window of the optical satellite to the ground line target, which is expressed as:

[0077]

[0078] According to a technical solution of the present invention, in step S3, the optical satellite visible window algorithm for the ground triangle area specifically includes:

[0079] The optical satellite-to-ground line target visibility window algorithm is applied to the three sides of the ground triangle area target, respectively, to obtain three time window sets, which are expressed as:

[0080]

[0081] First, the three time window sets are unioned, expressed as:

[0082]

[0083] Then traverse W tri For each time window in the adjacent two time windows [ts i ,te i ]、[ts i+1 ,te i+1 ], if te i+1 -ts i If the time window is less than half an orbital period, the two adjacent time windows are merged into one time window [ts i ,te i+1 ];

[0084] Repeat the above process until W tri There is no longer any te in the two adjacent time windows i+1 -ts iWhen the time window is less than half an orbital period, the final time window W tri It is the time window during which the optical satellite is visible to the ground triangle area target.

[0085] According to a technical solution of the present invention, in step S4, it specifically includes:

[0086] Let G = {P j1 |j1=1,2,…,M,M≥3}, where P j1 =(L j1 ,B j1 ) are the geodetic coordinates of the polygon vertices;

[0087] The longitude and latitude of the vertices of the polygonal area are used as the horizontal and vertical coordinates respectively, and the plane polygon triangulation method is used to divide the polygonal area into M-2 triangular areas, which are expressed as

[0088] T={T j2 ={P j2a ,P j2b ,P j2c |j2=1,2,…,M-2}}

[0089] Among them, P j2a =(L j2a ,B j2a ), P j2b =(L j2b ,B j2b ), P j2c =(L j2c ,B j2c ) are the geodetic coordinates of the vertices of the triangular area;

[0090] By applying the optical satellite visible window algorithm for the ground triangular area, a time window set of each triangular area visible to the optical satellite can be obtained, and the time window of the optical satellite visible to the arbitrary polygonal area target on the ground can be obtained by performing a union process, which is expressed as

[0091]

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

[0093] The present invention provides a method for calculating the visibility window of an optical satellite for ground targets. First, the method considers the optical satellite's orbital attitude motion and the optical satellite's optical camera installation and field of view parameters to establish a visibility judgment function for the optical satellite to the ground point target. A piecewise cubic polynomial function is used to fit the judgment function to calculate the visibility time window of the optical satellite to the ground point target. Secondly, based on the visibility window algorithm for ground point targets, a discretization method is used to calculate the visibility time window of the optical satellite to the ground line target. Then, based on the visibility window algorithm for ground line targets, the optical satellite's motion law is considered to calculate the visibility time window of the optical satellite to the ground triangular area target. Finally, based on the visibility window algorithm for the ground triangular area target, a polygon triangulation method is used to calculate the visibility time window of the optical satellite to the ground arbitrary polygonal area target. The present invention considers the optical satellite's attitude maneuvering and payload installation and pointing constraints, and adopts a unified technical framework to achieve fast and accurate calculation of the optical satellite for ground point, line, and area targets. This solves the problem that the calculation accuracy and efficiency of the traditional tracking propagation method are heavily dependent on the selected step size, resulting in low calculation efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0094] 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 the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be derived from these drawings without inventive effort.

[0095] Figure 1 A flowchart schematically illustrates a method for calculating a visible window of an optical satellite to a ground target provided in an embodiment of the present invention;

[0096] Figure 2 Schematically shows a specific flow chart of calculating the visibility time window of an optical satellite to a ground point target provided in an embodiment of the present invention;

[0097] Figure 3 A schematic diagram illustrating the principle of determining visibility of a ground point target by an optical satellite according to an embodiment of the present invention;

[0098] Figure 4 Schematically showing a flow chart of constructing a piecewise cubic interpolation polynomial provided in an embodiment of the present invention;

[0099] Figure 5 Schematically shows a flow chart of generating a time window in which an optical satellite is visible to a ground point target according to an embodiment of the present invention;

[0100] Figure 6 Schematically showing a schematic diagram of discretizing a ground line target into a set of point targets provided in an embodiment of the present invention;

[0101] Figure 7 The diagram schematically shows the division of any polygon on the ground into triangles according to an embodiment of the present invention. DETAILED DESCRIPTION

[0102] The description of the embodiments in this specification should be combined with the corresponding drawings, which should be considered a complete part of this specification. In the drawings, the shapes and thicknesses of the embodiments may be exaggerated and indicated for simplicity or convenience. Furthermore, the various structural components in the drawings will be described separately. It is worth noting that components not shown in the drawings or not described in words are known to those of ordinary skill in the art.

[0103] The description of the embodiments herein and any references to directions and orientations are for ease of description only and are not to be construed as limiting the scope of the present invention. The following description of the preferred embodiments may involve combinations of features, which may exist independently or in combination. The present invention is not specifically limited to the preferred embodiments. The scope of the present invention is defined by the claims.

[0104] like Figure 1 As shown, the present invention provides a method for calculating the visible window of an optical satellite to a ground target, comprising the following steps:

[0105] Step S1: Considering the orbital attitude motion of the optical satellite and the installation and field of view parameters of the optical camera on the optical satellite, a visibility judgment function of the optical satellite to the ground point target is established, and a piecewise cubic polynomial function is used to fit the judgment function to calculate the time window in which the optical satellite is visible to the ground point target;

[0106] In step S1, it specifically includes:

[0107] Step S11, calculating the ground-fixed position vector and the ground-fixed velocity vector of the optical satellite at any time within the time interval of the optical satellite ground-fixed ephemeris based on the optical satellite ground-fixed ephemeris;

[0108] In this embodiment, there is no constraint on the satellite orbit type, and no calculation of the satellite orbit is involved. Instead, the satellite Earth-fixed system ephemeris is directly used to describe the satellite orbit motion state.

[0109] In step S11, the satellite ephemeris is in the time interval [T s , T e ] is a set of discrete satellite position and velocity vectors, and the optical satellite fixed system ephemeris is expressed as t i Represents the ephemeris time node in the optical satellite ground-fixed ephemeris, R i and Represents the ephemeris time node t in the optical satellite ground-fixed ephemeris iThe corresponding Earth-fixed position vector and Earth-fixed velocity vector of the optical satellite;

[0110] For any time t∈[T s , T e ], the ground-fixed position vector and ground-fixed velocity vector of the optical satellite are expressed as R(t) = (x(t), y(t), z(t)) and It is obtained by segmented 5-times Lagrangian interpolation; the ephemeris time node in the optical satellite ground fixed system corresponding to time t is {t i+k |k=0,1…5,t i ≤t≤t i+5}, k represents the number of Lagrange interpolation, then

[0111]

[0112] Among them, l k (t) is the kth Lagrange interpolation basis function, expressed as

[0113]

[0114] Step S12: Calculating the attitude data of the optical satellite at any time within the time interval of the optical satellite Earth-fixed system ephemeris based on the optical satellite attitude sequence;

[0115] This embodiment has no constraints on the satellite attitude mode and does not involve the calculation of the satellite attitude. Instead, the satellite attitude motion state is directly described using the attitude sequence of the satellite itself relative to the VVLH coordinate system.

[0116] In step S12, the satellite attitude sequence is in the time interval [T s , T e ] is a series of discrete satellite attitude quaternion sets, and the optical satellite attitude sequence is expressed as Att = {(t i , Q i )|i=0,1…N ， t i ∈[T s , T e ]}, where Q i =(qa i ,qb i ,qc i ,qd i ) is t i The optical satellite body coordinate system at the moment is relative to the unit quaternion of the VVLH coordinate system. The first three are the vector part of the quaternion, and the fourth is the scalar part of the quaternion.

[0117] For any time t∈[T s , T e], the attitude data of the optical satellite is expressed as Q(t) = (qa(t), qb(t), qc(t), qd(t)), and the attitude data of the optical satellite is obtained by piecewise spherical linear interpolation according to the optical satellite attitude sequence. The time node in the optical satellite attitude sequence corresponding to time t is t i , t i+1 (t i ≤t≤t i+1 ), then

[0118]

[0119] Where k0, k1 and θ are the starting weight, ending weight and the angle between the two quaternions, respectively, expressed as

[0120]

[0121] θ=acos(qa i ×qa i+1 +qb i ×qb t+1 +qc i ×qc i+1 +qd i ×qd i+1 );

[0122] Step S13: Calculate the conversion matrix from the ground-fixed coordinate system of the optical satellite to the camera installation coordinate system at any time based on the optical satellite attitude data and the optical camera installation parameters of the optical satellite;

[0123] In step S13, the conversion matrix from the ground-fixed system of the optical satellite to the camera installation coordinate system is expressed as

[0124] M(t)=M sb M bo (t)M of (t)

[0125] Among them, M sb is the camera installation matrix, which is a constant matrix; M bo (t) is the transformation matrix from the VVLH coordinate system to the optical satellite body coordinate system, M of (t) is the transformation matrix from the ground-fixed system to the VVLH coordinate system of the optical satellite, M bo (t) is calculated from the attitude data Q(t) of the optical satellite at time t and is expressed as

[0126]

[0127] M of (t) is calculated from the ground-fixed position vector B(t) and ground-fixed velocity vector R(t) of the optical satellite at time t, and is expressed as

[0128]

[0129] I x (t) = I y (t)×I z (t)

[0130] Among them, I x (t), I y (t) and I z (t) represents the coordinates of the unit vectors of the V, V, L, and H axes in the Earth-fixed frame, V(t) represents the sum of the satellite Earth-fixed velocity vector and the Earth's rotation tangential velocity vector. ω=(0, 0, 7.29211581560071e-5) is the Earth's rotation angular velocity vector.

[0131] Step S14: Considering the optical satellite's attitude motion and the optical satellite's optical camera installation and field of view parameters, a judgment function is established to determine whether the optical satellite is visible to the ground point target;

[0132] In step S14, it specifically includes:

[0133] Step S141: Obtain the geodetic coordinates P (L, B, H) of the ground point target, and convert the geodetic coordinates of the ground point target into rectangular coordinates to obtain the ground fixed system vector R of the ground point target. p =(x, y, z), expressed as

[0134]

[0135] Where L is longitude, B is latitude, and H is altitude; x, y, and z represent the rectangular coordinates of the ground point target’s fixed vector; N is the radius of curvature of the y-axis. a and e are the semi-major axis and eccentricity of the Earth's ellipsoid, respectively;

[0136] Step S142: transform the ground-fixed vector of the ground target point into the camera installation coordinate system, and construct a judgment function for the ground point target within the optical satellite field of view, expressed as

[0137] r(t)=M(t)(R(t)-R p )

[0138]

[0139] Where α0 and β0 represent the vertical half angle and horizontal half angle of the field of view of the optical satellite, respectively; r y (t), r x (t), r z(t) respectively represent the coordinate components of the vector from the optical satellite to the ground point target in the camera installation coordinate system; F1(t) and F2(t) are the judgment functions of the ground point target in the field of view of the optical satellite;

[0140] Step S143: construct a judgment function for the minimum elevation angle of the optical satellite to the ground point target, expressed as

[0141]

[0142] Among them, n p The geoid normal vector representing the ground point target, n p =(cosBcosL, cosBsinL, sinB); θ0 is the minimum elevation angle of the optical satellite to the ground point target;

[0143] Step S145: Construct a ground solar altitude angle determination function, expressed as

[0144]

[0145] Where θ1 is the minimum solar altitude angle on the ground; R sun (t) is the sun position vector of the Earth-fixed system at time t;

[0146] Step S145: simultaneously satisfy the judgment function F i1 When (t+Δt)≤0 and i1=1, 2, 3, 4, the optical satellite is visible to the ground point target.

[0147] Step S15: establishing a first-order derivative function and a second-order derivative function of the determination function based on the central difference quotient;

[0148] In step S15, the first-order derivative function and the second-order derivative function of the optical satellite visibility judgment function for the ground point target are established based on the central difference quotient, which are expressed as

[0149]

[0150] Wherein, Δt is the difference quotient step size, which is fixed at 0.1 seconds.

[0151] Step S16: Adopting an adaptive step-size strategy to construct a piecewise cubic Hermite interpolation polynomial function fitting judgment function;

[0152] The judgment function F(t) is complex in form and it is difficult to determine its complete analytical expression. s , T e ] To find the root of the decision equation F(t)=0, we can only use the mechanical search method to solve it, which is computationally intensive and difficult to ensure that no roots are lost.

[0153] This embodiment uses two-point cubic Hermite interpolation to construct a piecewise cubic polynomial Ih (t) Function fitting judgment function. Construct interpolation node set T s =t0<t1…t n-1 <t n =T e , so that I h (t) Satisfy the following conditions:

[0154] aI h (t)∈C′[T s , T e ]

[0155] ikB h (t k )=F(t k ), I′ h (t k )=F′(t k )

[0156] cI h (t) = α3t 3 +α2t 2 +α1t+α0,t∈[t k , t k+1 ]

[0157] Therefore, the root-finding problem of the equation F(t)=0 is transformed into the root-finding problem of the interpolation polynomial, which can effectively reduce the difficulty of solving the problem and improve the efficiency of solving the problem.

[0158] In the time interval [T s , T e ] select a reasonable number of interpolation points so that the interpolation function I h (t) can accurately approximate the decision function F(t) and improve the computational efficiency as much as possible. This embodiment adopts the following principles to construct the interpolation nodes:

[0159] 1) Use the second-order derivative difference at the interpolation interval endpoints to control the interpolation step size

[0160] In the interval [t k , t k+1 ] is defined on the fitting error, expressed as

[0161] D(k)=D k+1 +D k

[0162]

[0163] Given error limit D error , if D(k)≥D error , then shorten the interpolation interval [t k , t k+1 ].

[0164] 2) Ensure that I h (t) is a monotonic function

[0165] Given an interpolation interval [t k , t k+1 ], construct the interpolation polynomial, expressed as

[0166] I h (t) = α3t 3 +α2t 2 +α1t+α0,t k ≤t≤t k+1

[0167] to I h (t) Take the derivative and we get

[0168] I′ h (t) = 3α3t 2 +2α2t+α1,t k ≤t≤t k+1

[0169] Find the quadratic polynomial equation I′ h (t)=0 in the open interval (t k , t k+1 ) in the root. If the equation has no roots, then the interpolating polynomial is in the interval [t k , t k+1 ] is monotonic within the range, accept the interpolation interval; if the equation has roots t k1 , t k2 (t k1 ≤t k2 ), then [t k , t k1 ] as the interpolation interval.

[0170] Following the above principles, such as Figure 4 As shown, in step S16, constructing a piecewise cubic Hermite interpolation polynomial function specifically includes:

[0171] Step S1601: Set the maximum interpolation step length mdt to one quarter of the orbital period;

[0172] Step S1602: Set the current interpolation node t c =T s , calculate the judgment function at t c The function value and derivative value F(t c ), F′(t c )、F″(t c ), save the current interpolation node data t c 、F(t c ), F′(tc )、F″(t c );

[0173] Step S1603: Calculate the interpolation step length δt=min(mdt, T e -t c );

[0174] Step S1604: Set the subsequent interpolation node t n =t c +δt, calculate the judgment function at t n The function value and derivative value F(t n ), F′(t n )、F″(t n );

[0175] Step S1605: In the interpolation interval [t c , t n ], construct the two-point cubic Hermite interpolation polynomial I h (t), calculate I h (t) The second-order derivative value I″ at the end point of the interpolation interval h (t c ), I″ h (t n );

[0176] Step S1606: Calculate the fitting error D = |I″ h (t c )-F″(t c )|+|I″ h (t n )-F″(t n )|, if D≥D error Then shorten the interpolation step δt=0.5δt and go to step S1604;

[0177] Step S1607: Calculate I′ h (t)=0 in the open interval (t k , t k+1 ) in the root, if the equation has roots t k1 , t k2 And t k1 ≤t k2 , then adjust the interpolation step δt=t k1 -t c , go to step S1604;

[0178] Step S1608: Save subsequent interpolation node data t n 、F(t n ), F′(t n )、F″(t n );

[0179] Step S1609: Update the current interpolation node data t c =t n 、F(t c )=F(t n ), F′(t c )=F′(t n )、F″(t c )=F″(t n );

[0180] Step S1610: If t c <T e , then go to step S1603, otherwise end the calculation process and obtain the difference node set T s =t0<t1…t n-1 <t n =T e , where n represents the number of interpolation nodes.

[0181] Step S17: Find the roots of the piecewise cubic Hermite interpolation polynomial function and calculate the time window set that satisfies the judgment function;

[0182] In each interpolation interval [t k , t k+1 ], if I h (t k )×I h (t k+1 )≤0, then there is only one root in this interval, which can be found by the cubic polynomial root formula. By traversing all interpolation nodes, we can calculate I h (t) in [T s , T e ] to generate a set of time windows that satisfy the judgment function.

[0183] In step S17, it specifically includes:

[0184] Step S171, set the time window flag isInWin=false;

[0185] Step S172: The difference node set T s =t0<t1…t n-1 <t n =T e The first interpolation node data in is used as the current interpolation node data t c 、F(t c ), F(t c )<0, set the window start time t s =t c , set isInWin = true;

[0186] Step S173: Obtain subsequent interpolation node data t n 、F(t n );

[0187] Step S174: If F(t c )×F(t n )≤0, calculate I h (t) in [t c , t n ] root t r Otherwise, go to Step 5; if isInWin is true, set the window end time t e =t r And save the window data; if isInWin is false, set the window start time t s =t r ; Modify the time window flag isInWin = !isInWin;

[0188] Step S175: Update the current interpolation node data t c =t n 、F(t c )=F(t n );

[0189] Step S176: If t c <T e , then go to step S173, otherwise go to step S177;

[0190] Step S177: If isInWin is true, set the window end time t e =t c And save the window data.

[0191] Step S18: Perform a set operation on the set that satisfies the judgment function to obtain the time window in which the optical satellite is visible to the ground point target.

[0192] The four judgment functions of the visibility of the optical satellite to the ground point target are fitted by piecewise cubic polynomial functions, and four time window sets can be obtained.

[0193] W F1 ={[ts i ,te i ]|i=0,1…n1}

[0194] W F2 ={[ts i ,te i ]|i=0,1…n2}

[0195] W F3={[ts i ,te i ]|i=0,1…n3}

[0196] W F4 ={[ts i ,te i ]|i=0,1…n4}

[0197] Performing intersection operation on these four time window sets can obtain the visible time window of the optical satellite to the ground point target

[0198] W site =W F1 ∩W F2 ∩W F3 ∩W F4

[0199] Step S2: Based on the time window in which the optical satellite is visible to the ground point target, a discretization method is used to calculate the time window in which the optical satellite is visible to the ground line target;

[0200] In step S2, the optical satellite to ground line target visibility window algorithm includes:

[0201] The two endpoints of the ground line target are P a (L a , B a ), P b (L b , B b ), the equidistant division method is used along the geodesic line to discretize the ground line target into a set of ground point targets {P i2 (L i2 , B i2 , 0)|i2=1,2,...,N,N≥2}, apply the optical satellite to ground point target visibility window algorithm to obtain the optical satellite to each ground point target in the ground point target set visibility time window, and perform union calculation to obtain the optical satellite to the ground line target visibility time window, expressed as:

[0202]

[0203] Step S3: Based on the time window in which the optical satellite is visible to the ground line target and taking into account the motion law of the optical satellite, the time window in which the optical satellite is visible to the ground triangular area target is calculated;

[0204] In step S3, the optical satellite visible window algorithm for the ground triangle area specifically includes:

[0205] For the three sides of the ground triangle area target, the optical satellite to ground line target visible window algorithm is applied respectively, and three time window sets can be obtained, which are expressed as:

[0206]

[0207] First, the three time window sets are unioned, expressed as:

[0208]

[0209] Then traverse W tri For each time window in the adjacent two time windows [ts i ,te i ]、[ts i+1 ,te i+1 ], if te i+1 -ts i If the time window is less than half an orbital period, the two adjacent time windows are merged into one time window [ts i ,te i+1 ];

[0210] Repeat the above process until W tri There is no longer any te in the two adjacent time windows i+1 -ts i When the time window is less than half an orbital period, the final time window W tri It is the time window in which the optical satellite is visible to targets in the triangular area on the ground.

[0211] Step S4: Based on the visible time window of the optical satellite to the ground triangular area target, a polygon triangulation method is used to calculate the visible time window of the optical satellite to the ground arbitrary polygon area target.

[0212] In step S4, it specifically includes:

[0213] Let G = {P j1 |j1=1,2,…,M,M≥3}, where P j1 =(L j1 ,B j1 ) are the geodetic coordinates of the polygon vertices;

[0214] The longitude and latitude of the vertices of the polygonal area are used as the horizontal and vertical coordinates respectively, and the plane polygon triangulation method is used to divide the polygonal area into M-2 triangular areas, which can be expressed as

[0215] T={T j2 ={P j2a ,P j2b ,P j2c |j2=1,2,…,M-2}}

[0216] Among them, P j2a =(L j2a ,Bj2a ), P j2b =(L j2b ,B j2b ), P j2c =(L j2c ,B j2c ) are the geodetic coordinates of the vertices of the triangular area;

[0217] Applying the visible window algorithm of optical satellite to the ground triangular area, we can obtain the time window set of optical satellite visible to each triangular area, and perform union processing to obtain the time window of optical satellite visible to the target in any polygonal area on the ground, which is expressed as

[0218]

[0219] It should be noted that, in this document, the terms "comprises," "includes," or any other variations thereof are intended to encompass non-exclusive inclusion, such that a process, method, article, or terminal device that includes a series of elements includes not only those elements but also other elements not explicitly listed, or also includes elements inherent to such process, method, article, or terminal device. In the absence of further limitations, an element defined by the phrase "comprises a..." does not exclude the presence of other identical elements in the process, method, article, or terminal device that includes the element.

[0220] Finally, it should be noted that the above is a preferred embodiment of the present invention. It should be noted that although the preferred embodiment of the present invention has been described, it is clear that those skilled in the art, once they understand the basic inventive concept of the present invention, can make various improvements and modifications without departing from the principles of the present invention. Such improvements and modifications should also be considered as within the scope of protection of the present invention. Therefore, the appended claims are intended to be interpreted as including the preferred embodiment and all changes and modifications that fall within the scope of the embodiments of the present invention.

Claims

1. A method for calculating the visibility window of an optical satellite to a ground target, characterized in that: The following steps are involved: Step S1: Considering the orbital attitude motion of the optical satellite and the installation and field of view parameters of the optical camera on the optical satellite, a visibility judgment function of the optical satellite to the ground point target is established, and a piecewise cubic polynomial function is used to fit the judgment function to calculate the time window in which the optical satellite is visible to the ground point target; Step S2: Based on the time window in which the optical satellite is visible to the ground point target, a discretization method is used to calculate the time window in which the optical satellite is visible to the ground line target; Step S3, based on the time window in which the optical satellite is visible to the ground line target and taking into account the motion law of the optical satellite, calculating the time window in which the optical satellite is visible to the ground triangular area target; Step S4: Based on the visible time window of the optical satellite to the ground triangular area target, a polygon triangulation method is used to calculate the visible time window of the optical satellite to the ground arbitrary polygon area target.

2. The method for calculating the visible window of an optical satellite to a ground target according to claim 1, characterized in that: In the step S1, it specifically includes: Step S11: Calculating the Earth-fixed position vector and the Earth-fixed velocity vector of the optical satellite at any time within the time interval of the optical satellite Earth-fixed ephemeris based on the optical satellite Earth-fixed ephemeris; Step S12: Calculating the attitude data of the optical satellite at any time within the time interval of the optical satellite Earth-fixed system ephemeris based on the optical satellite attitude sequence; Step S13, calculating the conversion matrix from the ground-fixed system to the camera installation coordinate system of the optical satellite at any time according to the optical satellite attitude data and the installation parameters of the optical camera of the optical satellite; Step S14: Considering the optical satellite attitude motion and the optical satellite optical camera installation and field of view parameters, a determination function is established as to whether the optical satellite is visible to the ground point target; Step S15: establishing a first-order derivative function and a second-order derivative function of the determination function based on the central difference quotient; Step S16: Adopting an adaptive step size strategy, constructing a piecewise cubic Hermite interpolation polynomial function to fit the determination function; Step S17, finding roots of the piecewise cubic Hermite interpolation polynomial function, and calculating a time window set that satisfies the determination function; Step S18: performing a set operation on the set that satisfies the judgment function to obtain a time window in which the optical satellite is visible to the ground point target.

3. The method for calculating the visible window of an optical satellite to a ground target according to claim 2, characterized in that: In step S11, the optical satellite earth-fixed system ephemeris is expressed as t i represents the ephemeris time node in the optical satellite ground-fixed ephemeris, R i and represents the ephemeris time node t in the optical satellite fixed system ephemeris i The corresponding Earth-fixed position vector and Earth-fixed velocity vector of the optical satellite; For any time t∈[T s , T e ], the ground-fixed position vector and ground-fixed velocity vector of the optical satellite are expressed as R(t)=(x(t), y(t), z(t)) and The ephemeris time node in the optical satellite ground fixed system corresponding to time t is obtained by segmented 5-times Lagrangian interpolation. i+k |k=0,1…5,t i ≤t≤t i+5 }, k represents the number of Lagrange interpolation, then Among them, l k (t) is the kth Lagrange interpolation basis function, expressed as In step S12, the optical satellite attitude sequence is represented as Att={(t i , Q i )|i=0,1…N,t i ∈[T s , T e ]}, where Q i =(qa i ,qb i ,qc i ,qd i ) is t i The optical satellite body coordinate system at the moment is relative to the unit quaternion of the VVLH coordinate system. The first three are the vector part of the quaternion, and the fourth is the scalar part of the quaternion. For any time t∈[T s , T e ], the attitude data of the optical satellite is expressed as Q(t) = (qa(t), qb(t), qc(t), qd(t)), the attitude data of the optical satellite is obtained by piecewise spherical linear interpolation according to the optical satellite attitude sequence, and the time node in the optical satellite attitude sequence corresponding to time t is t i , t i+1 (t i ≤t≤t i+1 ), then Where k0, k1 and θ are the starting weight, ending weight and the angle between the two quaternions, respectively, expressed as θ=acos(qa i ×qa i+1 +qb i ×qb i+1 +qc i ×qc i+1 +qd i ×qd i+1 ); In step S13, the conversion matrix from the ground-fixed system of the optical satellite to the camera installation coordinate system is expressed as M(t)=M sb M bo (t)M of (t) Among them, M sb is the camera installation matrix, which is a constant matrix; M bo (t) is the transformation matrix from the VVLH coordinate system to the optical satellite body coordinate system, M of (t) is the transformation matrix from the ground-fixed system to the VVLH coordinate system of the optical satellite, M bo (t) is calculated from the attitude data Q(t) of the optical satellite at time t and is expressed as M of (t) is the ground-fixed position vector R(t) and ground-fixed velocity vector of the optical satellite at time t Calculated, expressed as I x (t)=I y (t)×I z (t) Among them, I x (t), I y (t), I z (t) represents the coordinates of the unit vectors of the V, V, L, and H axes in the Earth-fixed frame, V(t) represents the sum of the satellite Earth-fixed velocity vector and the Earth's rotation tangential velocity vector. ω=(0, 0, 7.29211581560071e-5) is the Earth's rotation angular velocity vector.

4. The method for calculating the visible window of an optical satellite to a ground target according to claim 3, characterized in that: In the step S14, it specifically includes: Step S141: Obtain the geodetic coordinates P(L, B, H) of the ground point target, and convert the geodetic coordinates of the ground point target into rectangular coordinates to obtain the ground fixed vector R of the ground point target. p =(x, y, z), expressed as Wherein, L is longitude, B is latitude, H is altitude; x, y and z represent the rectangular coordinates of the ground point target fixed vector respectively; N is the radius of curvature of the circle, a and e are the semi-major axis and eccentricity of the Earth's ellipsoid, respectively; Step S142: transform the ground fixed system vector of the ground target point into the camera installation coordinate system, and construct a ground point target determination function within the optical satellite field of view, expressed as r(t)=M(t)(R(t)-R p ) Wherein, α0 and β0 represent the vertical half angle and horizontal half angle of the field of view of the optical satellite respectively; r y (t), r x (t), r z (t) respectively represent the coordinate components of the vector from the optical satellite to the ground point target in the camera installation coordinate system; F1(t) and F2(t) are the ground point target determination functions within the field of view of the optical satellite; Step S143: construct a judgment function for the minimum elevation angle of the optical satellite to the ground point target, expressed as Among them, n p represents the geoid normal vector of the ground point target, n p =(cosB cos L, cos B sin L, sin B); θ0 is the minimum elevation angle of the optical satellite to the ground point target; Step S144: construct a ground solar altitude angle determination function, expressed as Where θ1 is the minimum solar altitude angle on the ground; R sun (t) is the sun position vector of the Earth-fixed system at time t; Step S145: simultaneously satisfy the judgment function F i1 When (t+Δt)≤0 and i1=1, 2, 3, 4, the optical satellite is visible to the ground point target.

5. The method for calculating the visible window of an optical satellite to a ground target according to claim 4, characterized in that: In step S15, the first-order derivative function and the second-order derivative function of the optical satellite to ground point target visibility determination function are established based on the central difference quotient, which are expressed as Wherein, Δt is the difference quotient step size, which is fixed at 0.1 seconds.

6. The method for calculating the visible window of an optical satellite to a ground target according to claim 5, characterized in that: In step S16, constructing a piecewise cubic Hermite interpolation polynomial function specifically includes: Step S1601: Set the maximum interpolation step length mdt to one quarter of the orbital period; Step S1602: Set the current interpolation node t c =T s , calculate the judgment function at t c The function value and derivative value F(t c ), F′(t c )、F″(t c ), save the current interpolation node data t c 、F(t c ), F′(t c )、F″(t c ); Step S1603: Calculate the interpolation step length δt=min(mdt, T e -t c ); Step S1604: Set the subsequent interpolation node t n =t c +δt, calculate the judgment function at t n The function value and derivative value F(t n ), F′(t n )、F″(t n ); Step S1605: In the interpolation interval [t c , t n ], construct the two-point cubic Hermite interpolation polynomial I h (t), calculate I h (t) The second-order derivative value I″ at the end point of the interpolation interval h (t c ), I″ h (t n ); Step S1606: Calculate the fitting error D = |I″ h (t c )-F″(t c )|+|I″ h (t n )-F″(t n )|, if D≥D error Then shorten the interpolation step δt=0.5δt and go to step S1604; Step S1607: Calculate I′ h (t)=0 in the open interval (t k , t k+1 ) in the root, if the equation has roots t k1 , t k2 And t k1 ≤t k2 , then adjust the interpolation step δt=t k1 -t c , go to step S1604; Step S1608: Save subsequent interpolation node data t n 、F(t n ), F′(t n )、F″(t n ); Step S1609: Update the current interpolation node data t c =t n 、F(t c )=F(t n ), F′(t c )=F′(t n ), F(t″ c )=F″(t n ); Step S1610: If t c <T e , then go to step S1603, otherwise end the calculation process and obtain the difference node set T s =t0<t1…t n-1 <t n =T e , where n represents the number of interpolation nodes.

7. The method for calculating the visible window of an optical satellite to a ground target according to claim 6, characterized in that: In the step S17, it specifically includes: Step S171, set the time window flag isInWin=false; Step S172: The difference node set T s =t0<t1…t n-1 <t n =T e The first interpolation node data in is used as the current interpolation node data t c 、F(t c ), F(t c )<0, set the window start time t s =t c , set isInWin = true; Step S173: Obtain subsequent interpolation node data t n 、F(t n ); Step S174: If F(t c )×F(t n )≤0, calculate I h (t) in [t c , t n ] root t r Otherwise, go to Step 5; if isInWin is true, set the window end time t e =t r And save the window data; if isInWin is false, set the window start time t s =t r ; Modify the time window flag isInWin = !isInWin; Step S175: Update the current interpolation node data t c =t n 、F(t c )=F(t n ); Step S176: If t c <T e , then go to step S173, otherwise go to step S177; Step S177: If isInWin is true, set the window end time t e =t c And save the window data.

8. The method for calculating the visible window of an optical satellite to a ground target according to claim 1, characterized in that: In step S2, the optical satellite to ground line target visibility window algorithm includes: The two endpoints of the ground line target are respectively P a (L a , B a ), P b (L b , B b ), the equidistant division method is used along the geodesic line to discretize the ground line target into a set of ground point targets {P i2 (L i2 ,B i2 ,0)|i2=1,2,…,N,N≥2}, apply the optical satellite to ground point target visibility window algorithm to obtain the visible time window of the optical satellite to each ground point target in the ground point target set, and perform union calculation to obtain the visible time window of the optical satellite to the ground line target, which is expressed as:

9. The method for calculating the visible window of an optical satellite to a ground target according to claim 8, characterized in that: In step S3, the optical satellite visible window algorithm for the ground triangle area specifically includes: The optical satellite-to-ground line target visibility window algorithm is applied to the three sides of the ground triangle area target, respectively, to obtain three time window sets, which are expressed as: First, the three time window sets are unioned, expressed as: Then traverse W tri For each time window in the adjacent two time windows [ts i ,te i ]、[ts i+1 ,te i+1 ], if te i+1 -ts i If the time window is less than half an orbital period, the two adjacent time windows are merged into one time window [ts i ,te i+1 ]; Repeat the above process until W tri There is no longer any te in the two adjacent time windows i+1 -ts i When the time window is less than half an orbital period, the final time window W tri It is the time window during which the optical satellite is visible to the ground triangle area target.

10. The method for calculating the visible window of an optical satellite to a ground target according to claim 9, characterized in that: In the step S4, it specifically includes: Let G = {P j1 |j1=1,2,…,M,M≥3}, where P j1 =(L j1 ,B j1 ) are the geodetic coordinates of the polygon vertices; The longitude and latitude of the vertices of the polygonal area are used as the horizontal and vertical coordinates respectively, and the plane polygon triangulation method is used to divide the polygonal area into M-2 triangular areas, which are expressed as T={T j2 ={P j2a ,P j2b ,P j2c |j2=1,2,…,M-2}} Among them, P j2a =(L j2a ,B j2a ), P j2b =(L j2b ,B j2b ), P j2c =(L j2c ,B j2c ) are the geodetic coordinates of the vertices of the triangular area; By applying the optical satellite visible window algorithm for the ground triangular area, a time window set of each triangular area visible to the optical satellite can be obtained, and the time window of the optical satellite visible to the arbitrary polygonal area target on the ground can be obtained by performing a union process, which is expressed as

Images